Quantum corrections to a spin-orbit coupled Bose-Einstein Condensate

TL;DR

This paper investigates how quantum fluctuations, amplified by spin-orbit coupling, alter the properties and phase boundaries of a Bose-Einstein condensate, providing detailed calculations of damping rates and quantum corrections.

Contribution

It provides a systematic analysis of quantum corrections in spin-orbit coupled BECs, highlighting shifts in phase boundaries and decay processes not previously detailed.

Findings

Quantum fluctuations modify superfluid density, spin polarizability, and sound velocity.

The phase boundary shifts to smaller transverse fields due to quantum effects.

Landau damping dominates quasiparticle decay even at low temperatures.

Abstract

We study systematically the quantum corrections to a weakly interacting Bose-Einstein condensate with spin-orbit coupling. We show that quantum fluctuations, enhanced by the spin-orbit coupling, modify quantitatively the mean-field properties such as the superfluid density, spin polarizability, and sound velocity. We find that the phase boundary between the plane wave and zero momentum phases is shifted to a smaller transverse field. We also calculate the Beliaev and Landau damping rates and find that the Landau process dominates the quasiparticle decay even at low temperature.

Equations

$$h_0=\frac{(p_x-k_0{\sigma }_z{)}^2+p_y^2+p_z^2}{2}+\frac{\Omega }{2}{\sigma }_x,$$

$$H_{\mathrm{int}}=\frac{1}{2}\sum _{\sigma {\sigma }^{\prime }}\int {\mathrm{d}}^3\mathbf{r}{g}_{\sigma {\sigma }^{\prime }}{n}_{\sigma }(\mathbf{r})n_{{\sigma }^{\prime }}(\mathbf{r}),$$

$$\displaystyle\psi=\frac{e^{i\phi}}{\sqrt{2}}\left[\begin{array}[]{c}\sqrt{\rho+\zeta_{1}}e^{i\varphi}\\ -\sqrt{\rho+\zeta_{2}}e^{-i\varphi}\end{array}\right],$$

$$\mathcal{L}$$

$${\mathcal{L}}_{\mathrm{mf}}$$

$$\displaystyle\mathcal{G}^{-1}_{0}(i\omega_{n},\mathbf{q})=\left[\begin{array}[]{cccc}\mathcal{A}(\mathbf{q})&-\omega_{n}&0&-ik_{0}q_{x}\\ \omega_{n}&\mathcal{B}(\mathbf{q})&ik_{0}q_{x}&0\\ 0&-ik_{0}q_{x}&\mathcal{C}(\mathbf{q})&-\omega_{n}\\ ik_{0}q_{x}&0&\omega_{n}&\mathcal{D}(\mathbf{q})\end{array}\right],$$

$${\epsilon }_{\mathrm{ph}}=c_0\sqrt{q^2-\frac{2k_0^2{q}_x^2}{\Omega +2G_2}}=c_{0}q\sqrt{1-\frac{2k_0^2{\cos }^2\theta }{\Omega +2G_2}}\equiv c_{\theta }q,$$

$${\chi }_M={\mathcal{G}}_{0,{\zeta }_{-}{\zeta }_{-}}(q_x\to 0)/\rho =\frac{2}{\Omega -{\Omega }_{c,\mathrm{mf}}},$$

$${\mathcal{L}}_{\mathrm{eff}}=(K{\omega }_n^2+{\rho }_i{q}_i^2)|\varphi {|}^2,$$

$${\mathcal{L}}_{\mathrm{fluct}}$$

$$\delta \rho =\rho \left(⟨{\varphi }^2⟩+⟨{\phi }^2⟩\right)+\frac{⟨{\zeta }_{+}^2⟩+⟨{\zeta }_{-}^2⟩}{4\rho }.$$

$${\rho }_x$$

$$\delta ⟨{\sigma }_x⟩=2⟨{\phi }^2⟩+\frac{⟨{\zeta }_{-}^2⟩}{2{\rho }^2}.$$

$${\Omega }_c+2G_2-2\rho {\Sigma }_{{\zeta }_{-}{\zeta }_{-}}(0)=2k_0^2,$$

$${\chi }_M=\frac{2}{\Omega -{\Omega }_{c,\mathrm{mf}}-2\rho {\Sigma }_{{\zeta }_{-}{\zeta }_{-}}(0)},$$

$${\gamma }_{\mathrm{B}}$$

$${\gamma }_{\mathrm{L}}=\frac{3{\pi }^{3}q{T}^4}{40\rho c_{\theta }^4}{\left[1-\frac{2\Omega k_0^2{\cos }^2\theta }{(\Omega +2G_2{)}^2}\right]}^2\sqrt{1+\frac{2k_0^2{\sin }^2\theta }{\Omega -{\Omega }_{c,\mathrm{mf}}}}.$$

$${\epsilon }_{\mathrm{ph}}^2(\mathbf{q})=\frac{a(\mathbf{q})-\sqrt{a^2(\mathbf{q})-b(\mathbf{q})}}{2},$$

$${\epsilon }_{\mathrm{sp}}^2(\mathbf{q})=\frac{a(\mathbf{q})+\sqrt{a^2(\mathbf{q})-b(\mathbf{q})}}{2},$$

$${\epsilon }_{\mathrm{sp}}={\Delta }_0+\frac{q^2}{2m_{\mathrm{sp}}},$$

$${\epsilon }_{\mathrm{ph}}=c_{\theta }q+d_{\theta }{q}^3,$$

$${\Delta }_0=\sqrt{\Omega (\Omega +2G_2)},$$

$$m_{\mathrm{sp}}^{-1}=\frac{(G_2+\Omega )(2G_2+\Omega )+2k_0^2(\Omega +G_1+2G_2){\cos }^2\theta }{\sqrt{\Omega }(2G_2+\Omega {)}^{3/2}},$$

$$c_{\theta }=c_0\sqrt{1-\frac{2k_0^2{\cos }^2\theta }{\Omega +2G_2}},$$

$$d_{\theta }=\frac{1}{8c_{\theta }}\left(1-\frac{4k_0^2[(\Omega +2G_2)({\Omega }^2+(G_1+3G_2)\Omega +2(G_1+G_2{)}^2){\cos }^2\theta -k_0^2(\Omega +2(G_1+G_2){)}^2{\cos }^4\theta ]}{\Omega (\Omega +2G_2{)}^3}\right),$$

$${\epsilon }_{\mathrm{ph}}=\sqrt{G_1{q}_{\parallel }^2+\frac{1}{4}{q}_{\parallel }^4+\frac{G_2{q}_x^4}{4G_2+2\Omega }-\frac{(G_1+G_2)[2(G_1+G_2)+\Omega ]}{2\Omega (2G_2+\Omega )}{q}_{\parallel }^2{q}_x^2},$$

$$q_x\ll \sqrt{\Omega -2k_0^2}.$$

$$q_y,q_z\ll \sqrt{G_1}.$$

$$T\ll \Omega -2k_0^2.$$

$${\chi }_M={\mathcal{G}}_{0,{\zeta }_{-}{\zeta }_{-}}(q_x\to 0)/\rho =\frac{2}{\Omega +2G_2-2k_0^2}.$$

$$S_d(\mathbf{q})=\int \mathrm{d}\omega {\mathcal{G}}_{0,{\zeta }_{+}{\zeta }_{+}}(\omega ,\mathbf{q})/\rho ,S_s(\mathbf{q})=\int \mathrm{d}\omega {\mathcal{G}}_{0,{\zeta }_{-}{\zeta }_{-}}(\omega ,\mathbf{q})/\rho .$$

$${\mathcal{L}}_{\mathrm{eff}}$$

$${\mathcal{G}}_{0,\mathrm{eff}}^{-1}$$

Full text

I Introduction

The spin-orbit coupling, arising due to the interaction of a particle’s spin with its motion in an electric field plays a crucial role in various branches of physics, including topological insulators Hasan and Kane (2010); Qi and Zhang (2011), topological semimetals Yan and Felser (2017); Armitage et al. (2018), and Majorana fermions Elliott and Franz (2015). In bosonic systems, the interplay of the interparticle interaction and spin-orbit coupling gives rise to exotic Bose-Einstein condensates which have been investigated in a rich variety of systems, including magnons Rüegg et al. (2003); Sirker et al. (2004); Demokritov et al. (2006), excitons Hakioğlu and Şahin (2007); Can and Hakioğlu (2009); High et al. (2012, 2013), exciton-polaritons Carusotto and Ciuti (2013); Byrnes et al. (2014); Sala et al. (2015); Whittaker et al. (2018); Klembt et al. (2017); Zezyulin et al. (2018), and ultracold atoms Lin et al. (2011); Zhang et al. (2012); Ji et al. (2015); Wu et al. (2016).

For a weakly interacting Bose-Einstein condensate, the mean-field theory provides a reliable description of various physical properties Dalfovo et al. (1999). To reveal beyond mean-field effects, one method is to reach the strongly interacting regime, which can be achieved in exciton-polaritons because of the strong coupling between the exciton and photon Carusotto and Ciuti (2013); Kasprzak et al. (2006); Balili et al. (2007); Rodriguez et al. (2017); Fink et al. (2018), and for ultracold atoms strong interactions are accessible by means of Feshbach resonances Papp et al. (2008); Pollack et al. (2009); Navon et al. (2011). However, strong interactions reduce the lifetime of Bose-Einstein condensates significantly. Another method is to fine tune the interaction parameters such that the mean-field interactions almost cancel out Petrov (2015); Li et al. (2017); Cabrera et al. (2018); Semeghini et al. (2018); Jørgensen et al. (2018), making the quantum fluctuations unmasked. The spin-orbit coupling provides an alternative way to enhance interaction effects due to the increased density of states Zhai (2015). However, only a handful of theoretical studies have addressed the beyond mean-field effects Ozawa and Baym (2012); Cui and Zhou (2013); Zheng et al. (2013); Kawasaki and Holzmann (2017); Wu and Liang (2018), and a thorough analysis of the quantum fluctuations in spin-orbit coupled bosonic systems is still lacking.

In this paper, we systematically investigate the quantum corrections to a spin-orbit coupled Bose-Einstein condensate. We study a model system that is simple and general, potentially realizable in various platforms and already implemented with ultracold atom experiments Lin et al. (2011); Zhang et al. (2012); Ji et al. (2015). The model shows three novel condensation phases Lin et al. (2011); Li et al. (2012), namely the stripe, plane wave, and zero momentum phases. To demonstrate the interplay between interaction and spin-orbit coupling, we focus on the zero momentum phase, which is the simplest case capturing the essential physics of spin-orbit coupling and interactions.

We calculate quantum corrections to a number of physical properties, including the superfluid density, spin polarizability, and sound velocity. The superfluid density at the mean-field phase transition point between the plane wave and zero momentum phases becomes nonzero due to quantum fluctuations, and as a result, the phase transition point is shifted towards a smaller transverse field. The spin polarizability diverges at the corrected phase transition point but remains finite at the mean-field phase boundary, which seems to be consistent with a recent experiment Zhang et al. (2012). The sound velocity also acquires quantitative corrections, which may be detected in current ultracold atom experiments and provides a way to explore the beyond mean-field effects. Finally, we obtain an analytical result for the Landau decay rate of phonons at low temperature. Unlike the Beliaev decay predicted in Wu and Liang (2018), the Landau damping is not suppressed in the direction of spin-orbit coupling, making it the dominant mechanism for the quasiparticle decay.

II The model system

We consider a generic model of a spin-1/2 Bose gas with spin-orbit coupling, described by the single particle Hamiltonian (we set $\hbar=m=1$)

[TABLE]

where $\sigma_{i}$ with $i=x,y,z$ are the $2\times 2$ Pauli matrices. The one dimensional spin-orbit coupling, characterized by $k_{0}$, appears in many realistic systems, including ultracold gases Lin et al. (2011); Zhang et al. (2012); Ji et al. (2015); Wang et al. (2012); Cheuk et al. (2012) and semiconducting nanowires Quay et al. (2010); Mourik et al. (2012); Das et al. (2012). The model applies to several systems but to compare with experiments, we consider the cold atom setup where $k_{0}$ is given by the momentum transfer from the two Raman laser beams and $\Omega$ is the Rabi frequency of the Raman beams. The interaction between the particles can be written as

[TABLE]

where $n_{\sigma}$ is the density of particles with spin $\sigma=\uparrow,\downarrow$, and $g_{\sigma\sigma^{\prime}}=4\pi a_{\sigma\sigma^{\prime}}$ are the interaction strengths in different spin channels, with $a_{\sigma\sigma^{\prime}}$ being the corresponding $s$-wave scattering lengths in case of ultracold quantum gases. In the following we assume $g_{\uparrow\uparrow}=g_{\downarrow\downarrow}\equiv g$ and $g_{\uparrow\downarrow}\equiv g^{\prime}$, and correspondingly, $a_{\uparrow\uparrow}=a_{\downarrow\downarrow}\equiv a$ and $a_{\uparrow\downarrow}\equiv a^{\prime}$. It is convenient to define interaction parameters $G_{1}=g_{+}\rho$ and $G_{2}=g_{-}\rho$ with $g_{\pm}=(g\pm g^{\prime})/2$ and $\rho=N/V$ being the total particle density. In recent experiments Lin et al. (2011); Zhang et al. (2012); Ji et al. (2015), 87Rb atoms are employed and the interaction is almost $SU(2)$ invariant, with $G_{2}/G_{1}\approx 10^{-3}$. The typical interaction parameter is $G_{1}\approx 0.24k^{2}_{0}$ with the peak density $\rho\approx 0.57k^{3}_{0}$ Zhang et al. (2012). The dimensionless parameter $\sqrt{a^{3}\rho}\approx 0.004$ is small, ensuring that the condensate is in the weakly interacting regime and the perturbation calculations are controlled.

The mean-field phase diagram of this model has been extensively investigated, for a review see Zhai (2015). For small Rabi frequency, the condensate wave function is a superposition of two plane waves with different momenta, characterizing the stripe phase with density modulations in the ground state. In this phase, both the translational and $U(1)$ symmetries are broken, and therefore there are two branches of gapless excitations. Increasing the Rabi frequency $\Omega$, the system enters the plane wave phase, in which the bosons condense in a single plane wave state. There is only one branch of gapless excitations in this phase and the energy dispersion contains a roton minimum at finite momentum. Further increasing the Rabi frequency such that $\Omega>\Omega_{c,\mathrm{mf}}\equiv 2k^{2}_{0}-2G_{2}$, the system enters the zero momentum phase, where the roton minimum disappears and the phonon excitation spectrum resembles that of a Bose-Einstein condensate without spin-orbit coupling. To reveal the essential effect of interactions, we focus on the simplest zero momentum phase to reduce the effect of nontrivial mean-field energy dispersions in the plane wave and stripe phases.

The ground state wave function for the zero momentum phase is described by a spinor $\psi=\sqrt{\rho/2}[1,-1]^{T}$. To characterize excitations on top of the condensate, we introduce phase and number fluctuations, and write the spinor field as

[TABLE]

where $\phi$ is the total and $\varphi$ is the relative phase fluctuations of the condensate, and $\zeta_{1}$ and $\zeta_{2}$ are the density fluctuations for spin up and spin down particles, respectively.

We use the imaginary time path integral formalism. The Lagrangian density is obtained through the Hamiltonian as

[TABLE]

with $\mu$ being the chemical potential. It is convenient to introduce the density and spin fluctuations, $\zeta_{+}=(\zeta_{1}+\zeta_{2})/2$ and $\zeta_{-}=(\zeta_{1}-\zeta_{2})/2$, which are conjugate to $\phi$ and $\varphi$, respectively. We then expand $\mathcal{L}$ in terms of the new variables, and up to the second order, we get the mean-field Lagrangian density

[TABLE]

where the mean-field Green’s function in the momentum and frequency representation is

[TABLE]

here $\omega_{n}=2\pi nT$ is the Matsubara frequency (we set $k_{B}=1$), $\mathcal{A}(\mathbf{q})=\rho\mathbf{q}^{2}$, $\mathcal{B}(\mathbf{q})=\mathbf{q}^{2}/(4\rho)+g_{+}$, $\mathcal{C}(\mathbf{q})=\rho\mathbf{q}^{2}+2\Omega\rho$, and $\mathcal{D}(\mathbf{q})=\mathbf{q}^{2}/(4\rho)+g_{-}+\Omega/(2\rho)$. The chemical potential is determined by requiring $\langle\zeta_{+}\rangle=0$, and at the mean-field level we find $\mu=G_{1}+(k^{2}_{0}-\Omega)/2$, so the first order term of $\zeta_{+}$ vanishes and the mean-field Lagrangian density is quadratic. The diagonal elements of $\mathcal{G}_{0}$ are represented by Feynman diagrams shown in Fig. 1.

III Mean-field results

Before studying the beyond mean-field corrections, we first present the mean-field predictions of the physical properties we are interested in. These results are readily obtained from the mean-field Green’s function.

The mean-field excitation energy is determined by $\det\mathcal{G}^{-1}_{0}=0$. In the low momentum limit, we find the gapless phonon dispersion to be

[TABLE]

where $c_{0}=\sqrt{G_{1}}$ and $c_{\theta}$ is the sound velocity which depends on the angle $\theta$ between the directions of the momentum $\mathbf{q}$ and the $x$ axis. The mean-field sound velocities $c_{y}$ and $c_{z}$ are the same as the usual Bogoliubov sound velocity $c_{0}$. An intriguing feature is that the mean-field sound velocity in the $x$ direction, $c_{x}$, vanishes at the phase transition point between the plane wave and zero momentum phases. Besides the gapless phononic mode, there also exists a gapped mode which is dominated by spin excitation, with the mean-field gap given by $\Delta_{0}=\sqrt{\Omega(\Omega+2G_{2})}$.

The density and spin response functions are given by the Green’s functions $\mathcal{G}_{\zeta_{+}\zeta_{+}}$ and $\mathcal{G}_{\zeta_{-}\zeta_{-}}$, respectively. From the spin response function, the spin polarizability Martone et al. (2012); Li et al. (2012) can be obtained, and at the mean-field level, we get

[TABLE]

which diverges at the mean-field phase transition point.

An important quantity characterizing superfluidity is the superfluid density, which governs the total phase fluctuations. To get the superfluid density, we integrate out the $\varphi$, $\zeta_{-}$ and $\zeta_{+}$ fields and obtain an effective theory of $\phi$ (see Appendix A). In the low energy and long wave length limit, we find

[TABLE]

where $K$ is the zero momentum static density response function with its mean-field value being $1/g_{+}$, and $\rho_{y}=\rho_{z}=\rho$, $\rho_{x}=\rho[1-2k^{2}_{0}/(\Omega+2G_{2})]$ are the mean-field superfluid densities. From the effective Lagrangian, we see that the sound velocity is related to the superfluid density through $c_{i}=\sqrt{K^{-1}\rho_{i}}$.

Note that the superfluid density in the $x$ direction vanishes when $\Omega=2k^{2}_{0}-2G_{2}$. Formally, for smaller $\Omega$, the superfluid density becomes negative, which means that a state with nonzero phase gradient, i.e. the plane wave phase, is energetically more favorable. In other words, a vanishing superfluid density indicates a second order phase transition from the zero momentum phase to the plane wave phase.

In Zhang et al. (2016), the superfluid density $\rho_{x}$ is calculated from the current-current correlation function, which can be written in terms of the transverse spin polarization $\langle\sigma_{x}\rangle$ and the excitation gap $\Delta$ as $\rho_{x}=\rho\left(1+2k^{2}_{0}\Omega\langle\sigma_{x}\rangle/\Delta^{2}\right)$. Substituting the mean-field values $\Delta=\Delta_{0}$ and $\langle\sigma_{x}\rangle=-1$, we obtain from this the same result as given by the effective theory method above. Note that if the gap becomes larger or $\sigma_{x}$ is not fully polarized, the superfluid density $\rho_{x}$ will increase.

IV Beyond mean-field corrections

To study the lowest order (one-loop) beyond mean-field corrections, we expand the Lagrangian density up to the fourth order of the fields,

[TABLE]

The Feynman diagrams corresponding to the vertices are given in Fig. 2. Without the spin-orbit coupling, the one-loop corrections can be calculated analytically, and the results are given in Appendix B. In the main text we focus on the more interesting situation with nonzero spin-orbit coupling and calculate the one-loop corrections numerically. Since the parameter $G_{2}$ is small, we take it to be zero unless otherwise mentioned.

IV.1 Quantum depletion

Due to the quantum fluctuations, the condensate is depleted by a fraction of the total density. Up to the lowest order, the quantum depletion is given by (see Appendix B)

[TABLE]

Fig. 3 shows the quantum depletion as a function of the interaction strength and spin-orbit coupling for different transverse fields. The quantum depletion increases with the interaction strength. We find that it also increases with the spin-orbit coupling strength, which is consistent with previous results Ozawa and Baym (2012); Cui and Zhou (2013). As Fig. 3 shows, the quantum depletion increases with decreasing $\Omega$, which means that the quantum fluctuations are enhanced as the system approaches the phase transition point.

IV.2 Lee-Huang-Yang correction and chemical potential shift

We study the correction to the mean-field energy density, which is known as the Lee-Huang-Yang (LHY) correction Lee et al. (1957) $\mathcal{E}_{\mathrm{LHY}}$, and can be viewed as the zero point energy of the excitations Andersen (2004). With increasing $k_{0}$, the phonon mode softens, and therefore the zero point energy decreases. Fig. 4 (a) shows this behavior clearly. Remarkably, we find that $\mathcal{E}_{\mathrm{LHY}}$ becomes negative for large enough spin-orbit coupling. This leads to a non-monotonic dependence of $\mathcal{E}_{\mathrm{LHY}}$ on $G_{1}$: If we fix $k_{0}$ and increase $G_{1}$ from zero, then for small $G_{1}$ (large $k^{2}_{0}/G_{1}$), the LHY correction decreases from zero to negative; increasing $G_{1}$ further, the LHY correction will increase since it becomes positive for small $k^{2}_{0}/G_{1}$. The non-monotonic behavior of $\mathcal{E}_{\mathrm{LHY}}$ is most clearly seen at the phase transition point, see Fig. 4 (b).

We then calculate the correction to the chemical potential, which is given by the tadpole diagrams shown in Fig. 5. The numerical results of $\delta\mu$ are shown in Figs. 4 (c) and (d). As the LHY correction, the chemical potential shift decreases with increasing of $k_{0}$ and depends non-monotonically on $G_{1}$. This is expected, because the chemical potential shift can also be obtained as the first order derivative of the LHY energy with respect to the density.

IV.3 Superfluid density, phase boundary shift, and spin polarizability

To get the correction to the superfluid density, we first calculate the one-loop self-energy and then integrate out the massive fields $\varphi$, $\zeta_{-}$ and $\zeta_{+}$ to get the effective Lagrangian of the total phase fluctuations. The superfluid density in the $x$ direction is found to be

[TABLE]

where $\Sigma_{\zeta_{-}\zeta_{-}}(0)$ is the self-energy at zero frequency and momentum. There is no correction to $\rho_{y}$ and $\rho_{z}$ at zero temperature, consistent with the general result of superfluid density in Galilean invariant superfluids Leggett (1998).

Our numerical calculations show that $\Sigma_{\zeta_{-}\zeta_{-}}(0)$ is nonzero at the mean-field transition point. Consequently, the superfluid density also becomes nonzero at $\Omega=\Omega_{c,\mathrm{mf}}$, see Fig. 6 (a). Physically, this can be explained by the decrease of the transverse polarization $\langle\sigma_{x}\rangle$ and the increase of the spin gap $\Delta$. Because of the spin-orbit coupling, the spin of excited particles is not perfectly along the $x$ direction, and therefore the magnitude of the transverse spin polarization is reduced. Up to the lowest order, the deviation of spin polarization is (see Appendix B)

[TABLE]

We plot the numerical result of $\delta\langle\sigma_{x}\rangle$ in Fig. 6 (b). Another quantity that determines $\rho_{x}$ is the excitation gap. We obtain from the one-loop self-energy the correction to the mean-field gap and find it is positive, see Fig. 6 (c). Combining the behavior of $\delta\langle\sigma_{x}\rangle$ and $\Delta$, the non-monotonic dependence of $\rho_{x}$ on $k^{2}_{0}$ can be explained: The superfluid density increases with increasing $\delta\langle\sigma_{x}\rangle$ and $\Delta$, and with increasing $k_{0}$, $\delta\langle\sigma_{x}\rangle$ increases but $\Delta$ decreases. As a result, the superfluid density first increases and then decreases with increasing the spin-orbit coupling strength.

As we have explained before (see also Appendix C), the phase transition between the zero momentum and plane wave phases is characterized by the vanishing superfluid density, so Eq. (18) means that the phase transition point is shifted by quantum fluctuations. The new phase boundary is determined through

[TABLE]

where $\Sigma_{\zeta_{-}\zeta_{-}}(0)$ should be evaluated at $\Omega_{c}$. The solid lines in Fig. 6 (d) show the relative phase transition shift as a function of $k^{2}_{0}/G_{1}$ for different $G_{2}/G_{1}$. The shift becomes larger with decreasing $k^{2}_{0}/G_{1}$ and reaches its maximum at a critical spin-orbit coupling strength $k_{0,c}$, below which the plane wave phase is preempted by the stripe phase Li et al. (2012). We plot the phase boundary shift for $k_{0}$ larger than the mean-field critical value $k_{0,c}=\sqrt{2G_{2}(1+G_{2}/G_{1})}$ Li et al. (2012). It is possible that the mean-field critical spin-orbit coupling strength is shifted by quantum fluctuations, but this is beyond the scope of this paper and we expect that it does not change the results presented in Fig. 6 (d) qualitatively. We also calculate the phase boundary by minimizing the ground state energy $\mathcal{E}_{\mathrm{mf}}+\mathcal{E}_{\mathrm{LHY}}$. The technical details are given in Appendix C, and the phase boundary shifts obtained in this way are presented by the dots in Fig. 6 (d). As can be seen, the two methods predict the same results.

The self-energy $\Sigma_{\zeta_{-}\zeta_{-}}(0)$ also gives a correction to the spin polarizability,

[TABLE]

which diverges at the corrected phase boundary but becomes finite at the mean-field phase transition point. We have checked numerically that around $\Omega_{c}$, the dependence of the self-energy $\Sigma_{\zeta_{-}\zeta_{-}}(0)$ on $\Omega$ is weak, and therefore $\chi_{M}$ diverges as $1/(\Omega-\Omega_{c})$ close to the phase boundary, as predicted by the mean-field theory. The spin polarizability has been measured Zhang et al. (2012), and it seems that our one-loop result agrees better with the experimental data than the mean-field theory, see Fig. 7. However, the current experimental data cannot lead to a decisive conclusion and future experiments are required to verify our prediction.

IV.4 Sound velocity and damping rate

Using the one-loop results for the static density response $K^{-1}$ and the superfluid density $\rho_{x}$, we obtain the quantum corrected sound velocity in the $x$ direction, $c_{x}=\sqrt{K^{-1}\rho_{x}}$. At the corrected phase transition point, the sound velocity $c_{x}$ vanishes because of the vanishing superfluid density $\rho_{x}$. This is different from the result in Chen et al. (2017), where a nonzero sound velocity at the phase boundary has been predicted within the Hartree-Fock-Bogoliubov-Popov approximation.

Since the sound velocity goes to zero slower than the superfluid density, it is easier to detect the beyond mean-field effects through the measurement of the sound velocity. In Fig. 8 we plot the $c_{x}$ against $\sqrt{\delta\Omega}/k_{0}$, with $\delta\Omega=\Omega-\Omega_{c,\mathrm{mf}}$. For typical experimental parameters Lin et al. (2011); Zhang et al. (2012); Ji et al. (2015), the one-loop prediction deviates clearly from the mean-field behavior when $\sqrt{\delta\Omega}/k_{0}<0.1$. The sound velocity has been measured Ji et al. (2015), but the parameters are not close enough to the phase transition point. However, our prediction should be observable with current experimental methods.

Finally, we calculate the damping rate of phonons, for details see Appendix D. At zero temperature, the damping is due to the Beliaev process Beliaev (1958), i.e., an excitation decays into two with lower energy. In the small momentum limit ($q_{y},q_{z}\ll\sqrt{G_{1}}$ and $q_{x}\ll\sqrt{\Omega-2k^{2}_{0}}$), we find

[TABLE]

which coincides with the result obtained in Wu and Liang (2018). The Beliaev damping is strongly suppressed along the direction of the spin-orbit coupling.

At finite temperature, the Landau damping Hohenberg and Martin (1965) arises because the phonon couples to thermal excitations. The Landau damping is experimentally more relevant since it is responsible for damping in trapped Bose gases Pitaevskii and Stringari (1997); Liu (1997); Fedichev et al. (1998). In the low temperature and small momentum limit ($c_{\theta}q\ll T\ll\Omega-2k^{2}_{0}$), we obtain

[TABLE]

Because of the extra $c_{\theta}$ dependence, the Landau damping rate, unlike the Beliaev decay, is not suppressed in the direction of spin-orbit coupling, which means that the Landau process is the dominant damping mechanism even for uniform systems at very low temperature.

V Conclusions

We calculate systematically the one-loop corrections to a spin-orbit coupled Bose-Einstein condensate. We find that quantum fluctuations cause quantitative modifications to the superfluid density, spin polarizability, sound velocity, and damping rate. The quantum depletion increases while the LHY energy decreases with the transverse field in the zero momentum phase. The phase boundary between the plane wave and zero momentum phases is shifted to a smaller transverse field. The superfluid density vanishes and the spin polarizability diverges at the one-loop phase transition point. But at the mean-field phase boundary, the spin polarizability remains finite, consistent with an experimental measurement Zhang et al. (2012). We also point out that the beyond mean-field corrections may be detected through the measurement of the sound velocity, and give the parameter regime in which the deviation from the mean-field behavior is visible. We calculate the Beliaev and Landau damping rates and identify the Landau damping as the dominant mechanism of quasiparticle decay. Our results show that the spin-orbit coupling leads to, even for moderate interactions, quantum fluctuations strong enough to make detectable modifications to the properties of a macroscopic quantum state such as a Bose-Einstein condensate. The results can be readily tested in ultracold quantum gases, and in the future, in spin-orbit coupled Bose-Einstein condensates realized in other systems.

VI Acknowledgements

This work was supported by theAcademy of Finland under Projects No. 303351, No. 307419, No. 318987, and by the European Research Council (ERC-2013-AdG-340748-CODE). L.L. would like to acknowledge the Aalto Centre for Quantum Engineering for support.

Appendix A Mean-field results

In this section we present the mean-field results of the excitation energy, density and spin response function, and superfluid density with some detailed derivations.

A.1 Excitation energy

The excitation energy is determined by $\det\mathcal{G}^{-1}_{0}=0$, which gives

[TABLE]

where $a(\mathbf{q})=\mathcal{A}(\mathbf{q})\mathcal{B}(\mathbf{q})+\mathcal{C}(\mathbf{q})\mathcal{D}(\mathbf{q})+2k^{2}_{0}q^{2}_{x}$, $b(\mathbf{q})=4[\mathcal{A}(\mathbf{q})\mathcal{D}(\mathbf{q})-k^{2}_{0}q^{2}_{x}][\mathcal{B}(\mathbf{q})\mathcal{C}(\mathbf{q})-k^{2}_{0}q^{2}_{x}]$, $\varepsilon_{\mathrm{ph}}$ is the gapless phonon mode, and $\varepsilon_{\mathrm{sp}}$ is the gapped mode which is dominated by spin excitations. In the small momentum limit,

[TABLE]

where

[TABLE]

with $c_{0}=\sqrt{G_{1}}$ being the usual Bogoliubov sound velocity for a weakly interacting single component Bose-Einstein condensate. In the absence of spin-orbit coupling, the sound velocity is the same as $c_{0}$. In the presence of spin-orbit coupling, it depends on $\theta$, which is the angle between the momentum $\mathbf{q}$ and direction of the spin-orbit coupling. When $\Omega=2k^{2}_{0}-2G_{2}$, the sound velocity along the $x$ direction becomes zero, and the phonon dispersion along the $x$ direction becomes quadratic,

[TABLE]

with $q^{2}_{\parallel}=q^{2}_{y}+q^{2}_{z}$.

Knowing the low energy dispersion relation of the phonons, we can define the momentum region in which the dispersion is linear. When the momentum is along the $x$ direction, by requiring $c_{x}q_{x}\gg d_{x}q^{3}_{x}$, we find the condition

[TABLE]

When the momentum is along the $y$ or $z$ direction, the condition is

[TABLE]

At finite temperature, the linear dispersion region also requires that the dispersion of the thermal excitations is linear, and this leads to the condition

[TABLE]

These conditions are used in deriving the analytical expressions for Beliaev and Landau damping rates.

A.2 Density and spin response functions

In the modulus-phase representation, the density and spin response functions are given by the Green’s functions $\mathcal{G}_{\zeta_{+}\zeta_{+}}$ and $\mathcal{G}_{\zeta_{-}\zeta_{-}}$, respectively. So the spin polarizability defined in Martone et al. (2012); Li et al. (2012) is simply given by $\mathcal{G}_{\zeta_{-}\zeta_{-}}(q_{x}\to 0)/\rho$, and at the mean-field level,

[TABLE]

The mean-field density and spin static structure factors are given by

[TABLE]

We show the mean-field static structure factors for different spin-orbit coupling strength in Fig. 9. As comparison, the contribution of the phonon branch are also shown. Without spin-orbit coupling, the density and spin excitations are decoupled and the phonon branch does not contribute to the spin structure factor. In the presence of spin-orbit coupling, a density perturbation along the $x$ direction also induces a spin response and vice versa, so the density and spin structure factors are carried by both the phonon and gapped excitations. In the large momentum limit, the total static structure factors approach to 1 and the phonon branch contributes to one half. Remarkably, we find a peak in the total spin static structure factor. When the parameter approaches to the phase transition point, the peak becomes higher and its location moves to the zero momentum. By contrast, the peak is not observed in the total density structure factor, although there is peak in the contribution of the phonon branch.

A.3 Superfluid density

To get the superfluid density, we integrate out the $\varphi$ and $\zeta_{-}$ fields and obtain an effective theory of $\phi$ and $\zeta_{+}$

[TABLE]

where

[TABLE]

which in the low energy limit is

[TABLE]

Integrating out the $\zeta_{+}$ field, we arrive at an effective Lagrangian of the phase fluctuation, and in the low energy and long wave length limit,

[TABLE]

where $K$ is the zero momentum static density response function whose mean-field value is $1/g_{+}$, and the mean-field superfluid densities are

[TABLE]

Appendix B Analytical results of one-loop corrections in the absence of spin-orbit coupling

Without the spin-orbit coupling, we can calculate the one-loop corrections analytically. It is useful to calculate the following integral,

[TABLE]

where ${}_{2}F_{1}(a,b;c;z)$ is the hypergeometric function. To get the above result we have used dimensional regularization.

The condensate fraction is

[TABLE]

so the quantum depletion is

[TABLE]

with $x=\Omega/(g_{-}\rho)$ and $E(z)$ and $K(z)$ are the complete elliptic integral of the second and first kind, respectively. The quantum depletion increases with increasing $g_{+}$ and $g_{-}$, but decreases with increasing $\Omega$.

The transverse spin polarization is

[TABLE]

so

[TABLE]

The Lee-Huang-Yang correction Lee et al. (1957) can be obtained as the zero point energy of the system Andersen (2004), and we find

[TABLE]

The first term in Eq. (60) is the same as the result for a weakly-interacting spinless Bose gas Lee et al. (1957). The second term comes from the spin excitation. The function ${}_{2}F_{1}\left(-\frac{1}{2},\frac{3}{2};3;\frac{2}{2+x}\right)$ depends weakly on $x$, with ${}_{2}F_{1}\left(-\frac{1}{2},\frac{3}{2};3;1\right)=32/(15\pi)$ and ${}_{2}F_{1}\left(-\frac{1}{2},\frac{3}{2};3;0\right)=1$, so the second term increases with increasing $g_{-}$ and $\Omega$.

The chemical potential shift is given by the tadpole diagrams shown in Fig. 5. Evaluating the integrals, we find

[TABLE]

which increases with $g_{+}$, $g_{-}$, and $\Omega$. Another way to calculate the chemical potential shift is to take derivative of the LHY energy density $\mathcal{E}_{\mathrm{LHY}}$ with respect to $\rho$, $\delta\mu=\partial_{\rho}\mathcal{E}_{\mathrm{LHY}}$, and the result is the same as Eq. (61).

The correction to $K^{-1}$ is given by $\Sigma_{\zeta_{+}\zeta_{+}}(0)$, and in the absence of spin-orbit coupling,

[TABLE]

Note that $\delta K^{-1}$ can be related to $\delta\mu$ through $\delta K^{-1}=\partial_{\rho}\delta\mu$.

To calculate the correction to the mean-field excitation gap $\Delta_{0}$, we need to compute the self-energies $\Sigma_{\varphi\varphi}(\Delta_{0})$, $\Sigma_{\zeta_{-}\zeta_{-}}(\Delta_{0})$, and $\Sigma_{\varphi\zeta_{-}}(\Delta_{0})$, which can also be done analytically in the absence of the spin-orbit coupling, and we find that the one-loop correction to the gap is zero. In the presence of spin-orbit coupling, we calculate the self-energies numerically, and find the one-loop correction increases the gap slightly, see the main text.

Appendix C Phase boundary between the plane wave and zero momentum phases: the effect of the LHY energy

In this section we study the phase boundary between the zero momentum and plane wave phases by minimizing the ground state energy. As we will show, this also provides another way to calculate the superfluid density.

We only consider the plane wave and zero momentum phases, and in general the field operator can be written as

[TABLE]

where we have introduced the phase fluctuations $\phi$ and $\varphi$, and the density fluctuations $\zeta_{1}$ and $\zeta_{2}$, which are simply set to be zero in the mean-field approximation. The parameters $k_{1}$ and $\alpha$ should be determined by minimizing the ground state energy, and $k_{1}\neq 0$ characterizes the plane wave phase while $k_{1}=0$ gives the zero momentum phase.

Substituting Eq. (65) to the Lagrangian density

[TABLE]

and up to the quadratic order of the fluctuations, we find

[TABLE]

where $\zeta_{+}=\zeta_{1}\cos^{2}{\alpha}+\zeta_{2}\sin^{2}{\alpha}$, $\zeta_{-}=\zeta_{1}\cos^{2}{\alpha}-\zeta_{2}\sin^{2}{\alpha}$, and $\mathcal{G}^{-1}_{0}$ in the momentum and frequency representation reads

[TABLE]

We choose the renormalization condition $\langle\zeta_{+}\rangle=\langle\zeta_{-}\rangle=0$, which gives two conditions at the mean-field level

[TABLE]

The first condition Eq. (73) determines the mean-field chemical potential and the second condition Eq. (74) gives a relation between $\alpha$ and $k_{1}$. Note that for small $k_{1}$, we have $\cos{2\alpha}\propto k_{1}$.

The mean-field energy density is given by the first line in Eq. (67),

[TABLE]

Note that Eq. (74) can also be obtained by minimizing the energy with respect to $\alpha$. In Li et al. (2012), a relation between $k_{1}$ and $\alpha$ is obtained by minimizing the energy with respect to $k_{1}$, which leads to $\alpha=\arccos{(k_{1}/k_{0})}/2$. This relation and Eq. (74) determine the mean-field value of $\alpha$ and $k_{1}$ and therefore the mean-field phase boundary, which are the same as the results in Li et al. (2012). However, $\alpha=\arccos{(k_{1}/k_{0})}/2$ no longer holds when the LHY energy is taken into account because in this case there will be extra contribution to the energy density depending on $k_{1}$. In contrast, Eq. (74) is still valid up to at least one-loop since there is no one-loop correction proportional to $\zeta_{-}$ and therefore $\langle\zeta_{-}\rangle=0$ leads to the same condition. Therefore, to include the effects of the LHY energy, we should utilize the condition Eq. (74) instead of the form used in Li et al. (2012).

Using Eq. (74), we can rewrite $\mathcal{E}_{\mathrm{mf}}$ in terms of $k_{1}$, and then we can view the resultant expression as a Landau functional in terms of the ‘order parameter’ $k_{1}$. The disordered phase corresponds to the zero momentum phase while the ordered phase is the plane wave phase. Technically, it is simpler to use $x\equiv\cos{2\alpha}$ as the order parameter (because $x\propto k_{1}$ for small $k_{1}$) and we have

[TABLE]

By minimizing the above express with respect to $x$, we can determine the mean-field phase diagram.

Expanding Eq. (76) around $x=0$ and rewriting the result in terms of $k_{1}$, we get

[TABLE]

It is then clear that the mean-field phase transition point is determined by $1-\frac{2k^{2}_{0}}{\Omega+2G_{2}}=0$. In the zero momentum phase, the coefficient before $k_{1}$ measures the energy cost of the phase fluctuations, and therefore it is by definition the superfluid density $\rho_{x}/2$. From the point view of the Landau theory of phase transitions, the superfuid density is the coefficient of the quadratic term in the order parameter expansion. A negative superfluid density simply means that the zero momentum phase is unstable, and $k_{1}$ will acquire a nonzero expectation value such that the system enters the plane wave phase. In the plane wave phase, the superfluid density becomes positive again.

To calculate the correction to the mean-field phase boundary, we include the LHY contribution to the ground state energy density and minimize $\mathcal{E}_{\mathrm{mf}}+\mathcal{E}_{\mathrm{LHY}}$ as a function of $k_{1}$. The LHY energy is obtained through the excitation energy determined by $\det{\mathcal{G}^{-1}_{0}}=0$ with $\mathcal{G}^{-1}_{0}$ given by Eq. (72). The minimization can be done in the following way: We first calculate numerically the LHY energy for small $k_{1}$, and then extract the coefficient of the $k^{2}_{1}$ term in $\mathcal{E}_{\mathrm{LHY}}$. This coefficient gives a correction to the coefficient of $k^{2}_{1}$ in Eq. (77), and the new phase boundary is determined by requiring the corrected coefficient to be zero. As shown in Fig. 6 (d), the phase boundary determined in this way agrees perfectly with the one determined through the one-loop result of $\rho_{x}$.

Before closing this section, we mention that the same method can be used to get the superfluid density in the plane wave phase. Assuming $\mathcal{E}_{\mathrm{mf}}$ reaches its minimal at $k_{1,c}$, then the superfluid density is obtained by expanding the mean-field energy Eq. (76) around $k_{1,c}$,

[TABLE]

where $\delta k=k_{1}-k_{1,c}$. To find $k_{1,c}$ we minimize Eq. (76) with respect to $x$ and find the position $x_{c}$ at which the energy takes minimum. Then using Eq. (74), we find $k_{1,c}=k_{0}\sqrt{1-\Omega^{2}/\Omega^{2}_{c,\mathrm{mf}}}$. Expanding Eq. (76) around $x_{c}$ and change the variable from $x-x_{c}$ to $k_{1}-k_{1,c}$, we obtain the mean-field superfluid density in the plane wave phase

[TABLE]

which is the same as the result in Zhang et al. (2016). By taking into account the LHY contribution, we can also obtain the correction to the mean-field superfluid density in the plane wave phase.

Appendix D The damping rate at zero and finite temperature

The damping rate $\gamma$, i.e., the imaginary part of the phonon excitation energy, is determined by

[TABLE]

where $\Sigma(\varepsilon_{\mathrm{ph}}+i0^{+},\mathbf{q})$ is the one-loop self-energy evaluated at the phonon frequency.

To solve Eq. (80), we first integrate out the $\varphi$ and $\zeta_{-}$ fields and obtain an effective theory for the low energy mode [c.f. Eq. (38)]

[TABLE]

where $\mathcal{G}^{-1}_{\mathrm{eff}}$ can be written as

[TABLE]

And then from Eq. (86), the damping rate is obtained

[TABLE]

We focus on the linear dispersion regime defined through Eqs. (33)-(35). By analyzing the low energy and momentum behavior of all the one-loop self-energies, we find that it is enough to consider the Feynman diagrams constructed from only two vertices Figs. 2 (a) and (d), and the momentum dependence of vertex Fig. 2 (d) can be neglected. Therefore the relevant parts of the effective self-energy matrix is

[TABLE]

As an example, we calculate $\Sigma_{\phi\phi}(i\omega_{n},\mathbf{q})$ explicitly. The Feynman diagrams are shown in Fig. 10.

[TABLE]

We write the noninteracting Green’s function explicitly

[TABLE]

Since we are studying the damping rate in the linear regime, the gapped branch can be neglected, and it is enough to know the low momentum behavior of $B_{11}$, $B_{22}$, and $B_{12}$,

[TABLE]

Evaluating the Matsubara frequency summation, $\Sigma_{\phi\phi}(i\omega_{n},\mathbf{q})$ can be written as

[TABLE]

with

[TABLE]

which is nonzero even if the temperature is zero and is relevant to the Beliaev damping rate, and

[TABLE]

which is nonzero only at finite temperature and is relevant to the Landau damping rate.

We calculate the imaginary part of $\Sigma_{\phi\phi,1}(\varepsilon_{\mathrm{ph}}(\mathbf{q})+i0^{+},\mathbf{q})$ at zero temperature,

[TABLE]

To calculate the above integral, we need to solve the internal $\mathbf{k}$ allowed by the energy and momentum conservation. We can scale the momentum as $c_{x}k_{x}\equiv c_{0}k^{\prime}_{x}$ and $k_{y/z}=k^{\prime}_{y/z}$, and then the phonon dispersion can be written as

[TABLE]

The momentum and energy conservation can be solved in terms of the new variables in the small $\mathbf{q}$ limit ($\theta^{\prime}$ is the angle between $\mathbf{k}^{\prime}$ and $\mathbf{q}^{\prime}$),

[TABLE]

with the restrition $k^{\prime}<q^{\prime}$. This means that $\mathbf{k}^{\prime}$ and $\mathbf{q}^{\prime}$ are along the same direction and $k^{\prime}<q^{\prime}$ and therefore $\mathbf{k}$ and $\mathbf{q}$ are also along the same direction and $k<q$. Under this condition,

[TABLE]

so

[TABLE]

We now calculate $\Im\Sigma_{\phi\phi,2}(\varepsilon_{\mathrm{ph}}(\mathbf{q})+i0^{+},\mathbf{q})$ at finite temperature,

[TABLE]

where

[TABLE]

To get Eq. (116) we have assumed $c_{\theta}q/T\ll 1$ and expand $n(\varepsilon_{\mathrm{ph}}(\mathbf{k+q}))-n(\varepsilon_{\mathrm{ph}}(\mathbf{k}))$ to the lowest order. In general it is difficult to solve the energy and momentum conserving condition $\delta(\varepsilon_{\mathrm{ph}}(\mathbf{q})+\varepsilon_{\mathrm{ph}}(\mathbf{k})-\varepsilon_{\mathrm{ph}}(\mathbf{k+q}))$ even if $\mathbf{q}$ is small, because $\mathbf{k}$ is not necessarily small and for general $\mathbf{k}$, the phonon dispersion is very complicated. However, if we focus on the low temperature region such that the corresponding phonon dispersion is linear, then we can replace $\varepsilon_{\mathrm{ph}}(\mathbf{k})$ by the linear dispersion because $\frac{\partial n(\varepsilon_{\mathrm{ph}}(\mathbf{k}))}{\partial\varepsilon_{\mathrm{ph}}(\mathbf{k})}$ decays rapidly when $\varepsilon_{\mathrm{ph}}(\mathbf{k})>T$. In this region the momentum and energy conservation is easily solved: $\mathbf{k}$ and $\mathbf{q}$ are along the same direction and the length of $\mathbf{k}$ is unrestricted. Under this condition $g(\mathbf{q},\mathbf{k})$ also takes a simple form

[TABLE]

and

[TABLE]

To get Eq. (121) from Eq. (120), we have used the condition $c_{\theta}q/T\ll 1$.

We can calculate other self-energies in the similar way, and here we just summarize the final results,

[TABLE]

From the above results we get the Beliaev damping rate at zero temperature

[TABLE]

and the Landau damping rate at finite temperature

[TABLE]

The Beliaev damping rate takes the same form as the result in Wu and Liang (2018), where a different method was used. The analytical expression for the Landau damping rate is obtained here for the first time.

If $G_{2}=g_{-}\rho=0$, the damping rates can be further simplified as

[TABLE]

Since $c_{\theta}=c_{0}\sqrt{1-\frac{2k^{2}_{0}\cos^{2}{\theta}}{\Omega+2G_{2}}}$, the Beliaev damping is strongly suppressed when the momentum is along the direction of the spin-orbit coupling. However, the Landau damping is not suppressed.