Bose-Einstein condensates in the presence of Weyl spin-orbit coupling

We consider two-component Bose-Einstein condensates subject to Weyl spin-orbit coupling. We obtain mean-field ground state phase diagram by variational method. In the regime where interspecies coupling is larger than intraspecies coupling, the system is found to be fully polarized and condensed at a finite momentum lying along the quantization axis. We characterize this phase by studying the excitation spectrum, the sound velocity, the quantum depletion of condensates, the shift of ground state energy, and the static structure factor. We find that spin-orbit coupling and interspecies coupling generally leads to competing effects.


I. INTRODUCTION
The creation of synthetic gauge fields in ultracold atomic gases provides fascinating opportunities for exploring quantum many-body physics [1]. Of particular interest is the realization of non-Abelian spin-orbit coupling (SOC) [2][3][4]. Spin-orbit coupling is crucial for realizing intriguing phenomena such as the quantum spin Hall effect [5], new materials classes such as topological insulators and superconductors [6][7][8]. In bosonic systems, the presence of SOC may lead to novel ground states that have no known analogs in conventional solidstate materials [9][10][11]. In cold atomic gases, spin-orbit coupling can be implemented by Raman dressing of atomic hyperfine states [12,13]. The tunability of the Raman coupling parameters promises a highly flexible experimental platform to explore interesting physics resulting from spin-orbit coupling [14]. Recently, twodimensional SOC has been experimentally realized in cold atomic gases [15,16].
In anticipation of immediate experimental relevance, intense theoretical attention has been paid to the physics of ultracold atomic gases in the presence of SOC [3,4]. In the absence of interparticle interactions, the lowlying density of states is two-dimensional for Rashbatype SOC [9]. In particular, the single-particle energy minimum featured a Rashba-ring, which has important consequences on the ground state and finite-temperature properties of SOC Bose gases [17][18][19][20][21][22][23][24][25][26], as the role of quantum fluctuations gets enhanced due to huge degeneracies at the lowest-lying states. The three-dimensional analog of Rashba-type SOC is interesting because it is expected to stabilize a long-sought skyrmion mode in the ground state of trapped Bose-Einstein condensates (BECS) [19,27,28]. This Weyl-type SOC can be implemented following the proposals [29][30][31] by using powerful quantum technology. Although there is currently no evidence for Weyl fermions to exist as fundamental particles in our universe, Weyl-like quasiparticles have been detected recently in condensed-matter systems [32,33]. In light of these discoveries, the study of Weyl SOC in ultracold atom systems becomes particularly relevant, since the ability of manipulate the Weyl-SOC strength creates interesting opportunities for the exploration of effects not predicted in the realm of particle physics. In addition, the study of the effects of SOC may reveal some interesting physics unexplored in conventional binary Bose condensates [34,35]. In this work, we shall examine the physics of two-component Bose gases subject to Weyltype SOC. Firstly, we will introduce the model and determine the mean-field ground state by variation approach. Secondly, we will set out to study a particular realization of ground state where quantum fluctuation plays an essential role. Specially, we will investigate the interplay of spin-orbit coupling and interspecies interaction upon the ground state properties of the system. Finally, we will come to a summary.

II. MODEL AND FORMALISM
We consider a 3D homogeneous interacting twocomponent Bose gas subject to Weyl-type spin-orbit coupling, described by the Hamiltonian H = H 0 + H I , with Here Ψ(r) = (ψ ↑ , ψ ↓ ) T is a two-component spinor field, σ =xσ x +ŷσ y +ẑσ z ,p is the momentum operator, n σ = ψ † σ ψ σ is the density for component σ ∈ {↑, ↓}, λ is the strength of the spin-orbit coupling, and the strength for the intraspecies interaction and interspecies interaction is g and g ↑↓ , respectively. For brevity, we set = 2m = 1 from now on.
Diagonalization of H 0 yields the two-branch singleparticle energy spectrum E ± (p) = p 2 ± λp, and the corresponding eigenfunctions are given by where V is the volume of the system. The lowest-energy state for a given propagating direction parameterized by θ p and ϕ p is from the "-" branch and occurs at momentum p = λ 2 (sin θ p cos ϕ p , sin θ p sin ϕ p , cos θ p ). To determine the ground state of an interacting system, as routinely done in the literature [10,23,[36][37][38], we assume that the system has condensed into a coherent superposition of two plane-wave states with opposite momenta with magnitude p = λ/2. Thus the condensate wave function adopts the form Φ 0 = C + Φ − (p) + C − Φ − (−p), where C + and C − are two complex numbers to be determined and subject to normalization condition |C + | 2 + |C − | 2 = n 0 . Without loss of generality, the normalization condition suggests the parametrization |C + | 2 = n 0 cos 2 (α/2) and |C − | 2 = n 0 sin 2 (α/2), with α ∈ [0, π]. Upon substitution into E G =< Φ 0 |H|Φ 0 >, the variational ground state energy per particle is evaluated as where f (θ p , α) = sin 2 θ p +sin 2 α−3 sin 2 θ p sin 2 α/2. Minimization of the ground state energy with respect to θ p and α, one obtains the ground state phase diagram, summarized in Fig.1. When g ↑↓ − g > 0, the system is found to be in the phase of PW-Polar, which is a fully polarized phase with condensation momentum lying along the quantization axis. When g ↑↓ − g < 0, at mean-field level, there are two degenerate phases: one is unpolarized PW-Axial phase, which is condensed at one plane-wave with momentum lying in the x-y plane; the other one is the SP-Polar phase, which is striped phase mixing of two opposite momentum along the z-axis. There exists a critical point when g ↑↓ − g = 0. In this case the system enjoys a SU(2) pseudo-spin rotation symmetry. To determine which phase the system prefers requires calculation going beyond mean field, and in principle it is believed to lead to a unique ground state via the mechanism of "order from disorder" [18,39]. Within the imaginary-time field integral, the partition function of the system may be cast as [40] , where β = 1/T is the inverse temperature and µ is the chemical potential introduced to fix the total particle number. Here, for simplicity, we restrict ourself to studying the PW-Polar phase. Without loss of generality, we further assume that the condensation occurs at momentum κ = (0, 0, −λ/2), then the ground state wave function is determined as It is a fully polarized phase with condensation momentum aligning antiparallel with the quantization axis. We split the Bose field into the mean-field part φ 0σ and the fluctuating part φ qσ as ψ qσ = φ 0σ δ q κ + φ qσ . After substitution, the action can be formally written as S = S 0 + S f , where is the mean-field contribution and S f denotes a contribution from the fluc- Panel (a): for g ↑↓ − g > 0, the ground state is in PW-Polar phase, where it is a plane wave with the condensation momentum being parallel to the z-axis. For g ↑↓ − g < 0, the system may be in the phase of either PW-Axial or SP-Polar. Here PW-Axial phase stands for one plane wave with the condensation momentum lying in the x-y plane, and SP-Polar phase stands for the condensation at two opposite momenta along the z-axis. Panel (b): schematic representation of the PW-Polar and PW-Axial phases. For one plane-wave condensation at either the north pole or the south pole is called the PW-Polar phase, it is a fully polarized phase as only one component is allowed. For one plane-wave condensation in the x-y plane is called PW-Axial, it is an unpolarized phase. tuating fields. The chemical potential may be determined via saddle point condition ∂S 0 /∂n 0 = 0, yielding µ = − λ 2 4 + 2gn 0 . At this point, the action is exact. However, it contains terms of cubic and quartic orders in fluctuating fields. To proceed, we resort to the celebrated Bogoliubov approximation, where only terms of zeroth and quadratic orders in the fluctuating fields are retained. By defining a four-dimensional column vec- , we can bring the fluctuating part of the action into the compact form q,σ ǫ qσ ,where w n = 2πn/β is the bosonic Matsubar frequencies, and the inverse Green's function G −1 (q, iw n ) is defined as where ǫ q↑ = q 2 +2gn 0 , ǫ q↓ = q 2 −2λq z +λ 2 +2(g ↑↓ −g)n 0 , and R q = λ(q x − iq y ). Throughout our calculation, we will choose gn 0 as a basic energy scale and √ gn 0 as the corresponding momentum scale. To characterize the strength of interspecies coupling, we define a dimensionless parameter η = g ↑↓ /g.

III. CALCULATION AND RESULTS
The excitation spectrum of the system can be found by examining the poles of the Green's function G(q, iw n ). To achieve this , one proceeds by evaluating the determinant of G −1 (q, iw n ), where ω 10 = q q 2 + 4gn 0 and ω 20 = λ 2 + q 2 + 2(g ↑↓ − g)n 0 , and q ⊥ = q 2 x + q 2 y . By solving the secular equation DetG −1 (q, ω qs ) = 0, one finds two branches of excitation spectrum ω q± . As seen from Eq. (5), the excitation spectrum enjoys the azimuthal symmetry. Therefore we only plot the spectrum along two typical directions in Fig. 2. Along the z-axis, the lower branch show the features of roton-maxon structure, indication of the tendency toward crystallization [36]. Such roton-maxon spectrum has been detected in recent Bragg spectroscopy experiments [41][42][43], and the spectrum is asymmetrical with respect to reversing the direction. In the x-y plane, the two branches are well separated as the upper branch is gapped while the lower branch becomes gapless as it approaches the origin q = (0, 0, 0).
Aside from the roton mode discussed above, the lower branch of the excitation spectrum also contains important information about the photon mode. Along z direction where q ⊥ = 0, it is straight forward to ana- lytically derive two branches of solutions from Eq. (5): ω − = ω 10 and ω + = −2λq z + ω 20 . The sound velocity along this direction is v s z = 2 √ gn 0 . In the x-y plane, lowenergy expansion around the gapless point (0, 0, 0) yields . Numerically we compute the sound velocity via v s (q) = lim q→0 ω − (q)/q. We find that the sound velocity varies with the polar angle θ q , as shown in Fig. 3. The sound velocity enjoys a symmetry of v s (θ q ) = v s (π − θ q ), with the maximum sound velocity achieved along z-axis and the minimum one in the x-y plane. Away from the critical point where η = 1, the spin-orbit coupling suppresses the sound velocity along any polar direction except for θ q = 0 and π, as indicated in Fig. 3(a). Interestingly, as seen in Fig. 3(b), suppression of sound velocity due to spin-orbit coupling could be mitigated by increasing the interspecies coupling, an indication of competing effects of spin-orbit coupling and interspecies coupling.
Being an intrinsic property of a BEC, the quantum depletion of the condensates provides vital information concerning the robustness of the superfluid state. The number density of exited particles can be evaluated by employing the quasi-particle's Green's function We show the density of the excited particles out of the condensates due to quantum fluctuation in Fig. 4. At a fixed interspecies coupling η, the quantum depletion is monotonically enhanced by spin-orbit coupling, and it reduces to the case of spinless Bose gases with n ex = (gn 0 ) 3/2 /(3π 2 ) in the absence of spin-orbit coupling [44], as seen in Fig. 4(a). At a fixed spin-orbit coupling strength, the interspecies coupling actually suppresses quantum depletion, signifying the competing effects of spin-orbit coupling and interspecies coupling upon quantum depletion. When the spin-orbit coupling is small, the effect of interspecies coupling decreases as well, as indicated in Fig. 4(b). This is quite remarkable, because there is only one species of condensation. In the absence of spin-orbit coupling, we do not expect that the the interspecies coupling plays any role in quantum depletion. We attribute this behavior to stemming from quantum fluctuation enhanced by spin-orbit coupling. The thermodynamic potential of this system is given by Ω = − ln Z/β = Ω 0 + Ω f , where the mean-field part is Ω 0 = −V gn 2 0 and the fluctuating part is Ω f = 1 2β T rlnG −1 − 1 2 qσ ǫ qσ . The thermodynamic potential Ω possesses an ultraviolet divergence, an artifact of zero range interaction, which can be removed either by replacing the bare interaction g with a T matrix [45] or by subtracting counter-terms [46]. At zero temperature, the ground-state energy becomes E G = Ω + µN , renormalized as Here E MF = V (gn 2 0 − λ 2 4 ) is the mean-field energy. We show the shift of ground state energy due to quantum fluctuation ∆E G = E G − E MF in Fig. 5. As seen in panel (a), at a fixed interspecies coupling η, the shift of the ground state energy ∆E G decreases monotonically with the strength of spin-orbit coupling λ. In the absence of the spin-orbit coupling and interspecies interaction, we have checked that the ground state energy E G recovers the well-known Lee-Huang-Yang result [47] for spinless and weakly-interacting Bose gases with E G /V = µn 2 (1 + a finite spin-orbit coupling, the shift of the ground state energy increases with interspecies coupling η, evidently shown in panel (b). The static structure factor S(q) probes density fluctuations of a system. It provides information on both the spectrum of collective excitations, which could be investigated at low momentum transfer, and the momentum distribution, which characterizes the behavior of the system at high momentum transfer, where the response is dominated by single-particle effects. At the Bogoliubov level, it can be evaluated as where A(q, iw n ) = (iw n + 2λq z ) 2 − ω 2 20 + λ 2 ω 20 sin 2 θ q . It is quite clear that the static structure factor possesses the cylindrical symmetry S(q) = S(q, θ q ). At q ⊥ = 0, the static structure factor adopts a close form as follows In this case, it recovers the Feynman relation [48,49], which connects the static structure factor to the excitations spectrum of a Bose system with time-reversal symmetry. We show the behavior of the static structure in Fig. 6. In the upper panel, we show the in-plane static structure factor S(q, θ q = π/2) in terms of in-plane momentum q ⊥ . It decreases as the spin-orbit coupling strength is increased, but increases as the interspecies coupling is increased. Such reversing trend signifies that spin-orbit coupling and interspecies coupling act with reversal role in the density response of the system. In the lower panel, we show angular dependence of the static structure factor at q = √ gn 0 . It is interesting to notice that S(q) is also symmetrical with reflection about the x-y plane, namely S(q, θ q ) = S(q, π − θ q ). The static structure factor develops its minimum along θ q = π/2. The spin-orbit coupling suppresses the density response greatly in the x-y plane, as seen in panel (c). In turn, the interspecies coupling enhances the density response greatly along the direction of θ q = π/2.

IV. SUMMARY AND CONCLUSIONS
To sum up, we have studied two-component Bose gases in the presence of Weyl-type SOC. We obtain the phase diagram via a variational approach. We find competing effects between spin-orbit coupling and interspecies coupling strength upon various properties of the PW-Polar phase. There is one crucial difference between them: spin-orbit coupling allows the process of pseudospin flipping process, while interspecies interaction does not permit that. This has far-reaching consequence in the quantum depletion of the condensates. In addition to cylindrical symmetry endorsed by the ground state where the condensation momentum lying along the quantization axis, the sound velocity and the static structure factor also enjoy a reflection symmetry with respect to x-y plane. We hope that our work will contribute to a deeper understanding of SOC BECs and the role of quantum fluctuations.