Magnetic and nematic phases in a Weyl type spin-orbit-coupled spin-1 Bose gas

We present a variational study of the spin-1 Bose gases in a harmonic trap with three-dimensional spin-orbit coupling of Weyl type. For weak spin-orbit coupling, we treat the single-particle ground states as the form of perturbational harmonic oscillator states in the lowest total angular momentum manifold with $j=1, m_j=1,0,-1$. When the two-body interaction is considered, we set the trail order parameter as the superposition of three degenerate single-particle ground-states and the weight coefficients are determined by minimizing the energy functional. Two ground state phases, namely the magnetic and the nematic phases, are identified depending on the spin-independent and the spin-dependent interactions. Unlike the non-spin-orbit-coupled spin-1 Bose-Einstein condensate for which the phase boundary between the magnetic and the nematic phase lies exactly at zero spin-dependent interaction, the boundary is modified by the spin-orbit-coupling. We find the magnetic phase is featured with phase-separated density distributions, 3D skyrmion-like spin textures and competing magnetic and biaxial nematic orders, while the nematic phase is featured with miscible density distributions, zero magnetization and spatially modulated uniaxial nematic order. The emergence of higher spin order creates new opportunities for exploring spin-tensor-related physics in spin-orbit coupled superfluid.

For a trapped Weyl coupled spin-1 BEC, one expects the emergence of topological objects with even higher spin order, i.e. the nematic order. An immediate question is that, what will the phase diagram look like, and in what a way will the 3D skyrmion and (or) the nematic order manifest themselves in the individual phases. Here we consider a spin-1 bosonic system subject to a weak Weyl SO coupling of s·p type in a harmonic trap. We implement a standard variational approach [3,4,24,25,52] to give a phase diagram of the system in different interaction regimes. A magnetic phase and a nematic phase are predicted in the ground state, and the latter is entirely new and has no analogue in the 3D SO coupled pseudo-spin-1/2 bosonic gases [44,46,48,52].
The paper is organized as follows. In Sec. II the energy functional for the 3D SO coupled spin-1 condensate is introduced in harmonic oscillator units. In the weak coupling limit we construct our variational order parameter in Sec. III and calculate the energy functional by means of the irreducible tensor method. The ground state phase diagram is determined by numerically minimizing the energy functional with respect to the variational parameters, and Sec. IV is devoted to an explicit illustration and detailed discussion of the densities and spin orders in the two ground state phases. We summarize our main results in Sec. V.

II. MODEL
We start from the mean-field Gross-Pitaevskii (GP) energy functional of spin-1 bosons with a Weyl type 3D SO coupling in the presence of a harmonic trap where the single particle part is with m the atomic mass and ω the trap frequency. Ψ = (ψ 1 , ψ 0 , ψ −1 ) T denotes the spinor order parameters for bosons with hyperfine components 1, 0, −1 respectively, and s = (s x , s y , s z ) are spin-1 matrices. The Weyl SO coupling s · p with strength λ is a 3D analogy of the Rashba configuration [55]. Without the trapping potential, the momentum p and the helicity s · p/ |p| are constants of motion, and the single-particle ground states are highly degenerate, i.e. the lower-energy helical branch achieves a minimum along a sphere of radius p SO = mλ, called SO sphere [44]. The presence of a trapping potential will partially lift this degeneracy, but at least a two-fold degeneracy related to the time-reversal symmetry will still remain [2]. The single particle Hamitonian is also invariant under simultaneous rotation of spin and coordinate space SO(3) R+S that leaves the total (instead of spin) angular momentum a good quantum number [46]. This breaking of rotational symmetry in spin space leads to spin-textured ground states with magnetic and nematic orders, as previously studied in Ref. [19,21,24,25,52]. The interaction energy functional is formulated in the standard form [57][58][59] where n (r) = n 1 (r) + n 0 (r) + n −1 (r) is the total particle density and n 1,0,−1 (r) = |ψ 1,0,−1 (r)| 2 are densities for the three components, respectively, and S = Ψ † sΨ is spin density. c 0,2 are spin-independent and spindependent interaction strengthes respectively, which are related to the two-body scattering lengths in the total spin-0 and spin-2 channels as c 0 = 4π(a 0 + 2a 2 )/3m and c 2 = 4π(a 2 − a 0 )/3m. The interaction is time-reversal (TR) symmetric under the operation T = e −iπsy K with K the complex conjugate. Besides, this two-body interaction is also SU(2) spin-rotation symmetric, which is different from the spin-1/2 bosons [44,46,48,52]. In the latter case the two body interaction is SU(2) symmetric only under the condition g ↑↑ = g ↓↓ = g ↑↓ , i.e. c = g ↑↓ /g ↑↑ = 1. Later we will see that this highly symmetric interaction will lead to highly degenerate ground states. In harmonic oscillator units, the system has length scale l T = /mω, energy scale ω, and SO coupling strength is in unit of ω/m. If we normalize the order parameter to unity, i.e., Ψ → N/l 3 T Ψ with N the total particle number in the condensate, the energy functional per particle is obtained as

III. VARIATIONAL APPROACH
Here we first introduce the single-particle ground states in the form of perturbational 3D harmonic oscillator in the weak SO coupling limit, the superposition of which constitutes the variational order parameter. We then calculate the energy functional of Eq. (4) using the irreducible tensor method.

A. Variational Order Parameter
The variational order parameter is just a spin-1 version of that presented in our recent work on the 3D SO coupled pseudo-spin-1/2 system [52]. We include some necessary contents here for completeness. It is expected that in the case of weak SO coupling the single particle energy is dominated by the 3D harmonic oscillator part, which has been verified in Ref. [44,50,52]. The solution of the 3D harmonic oscillator is well-known with the energy eigenvalues ǫ nrl = 2n r + l + 3 2 and eigenfunctions φ nr lm l (r, θ, ϕ) = R nr l (r) Y lm l (θ, ϕ). Here R nrl is the radial wave function and Y lm l is the spherical harmonics with n r the radial quantum number, l the orbital angular momentum quantum number, and m l its magnetic quantum number. In order to take into account the SO coupling term s · p later, it is convenient to choose the coupled representation of angular momentum for spin-1 particles, i.e. the complete set of commutative operators includes l 2 , s 2 , j 2 , j z where j = l + s and j z denote the total angular momentum and its z-component, respectively. The eigenfunction φ nr lm l (r, θ, ϕ) should be combined with the spin wave function χ ms (m s = 1, 0, −1) in the coupled representation as where Y l jmj (Ω) = m l ,ms C jmj lm l 1ms Y lm l χ ms is the vector spherical harmonics [60] with j = l + 1, l, |l − 1| (if l = 0, j = 1 only) and C jmj lm l 1ms the Clebsch-Gordan coefficients. In the coupled representation, the ground state wave function has n r = l = 0. This gives a total angular momentum j = 1 with m j = 1, 0, −1 and the three ground states are and respectively. The lowest few levels of the 3D harmonic oscillator is shown in Fig. 1. The single particle spectrum may be understood as a weak perturbation (slight level mixing) of the harmonic-oscillator energy levels in the case of weak SO coupling, which conserves the total angular momentum j and j z but l is no longer a good quantum number. This means SO coupling term would couple s, p and even d waves into the ground state with the same j and j z as illustrated in Fig. 1. In the case of strong SO coupling, the energy spectrum is weakly dispersive or nearly flat [32,44,50] which will not be considered in this work. The ground state wave functions are thus the superposition of the lowest s, p and d wave states with j = 1. The state with m j = 1 is Here A l (l = 0, 1, 2) are the weight coefficients with the normalization constraint 2 l=0 |A l | 2 = 1, and i in front of A 1 comes from the pure imaginary matrix elements of SO coupling between φ 0011 (φ 0211 ) and φ 0111 [44,50,52]. Explicitly this state is a spinor Moreover, in the single particle level, the energies are irrelevant to the magnetic quantum number of j z . Thus the other two states with m j = 0 and −1 are respectively. Note that the state Φ j=1,mj =−1 is the time reversal of Φ j=1,mj =1 , and T Φ j=1,mj=0 = −Φ j=1,mj =0 . We have neglected contribution from the excited 1s wave state φ 1011 shown in Fig. 1 (gray dashed), which can be absorbed into the lowest 0s state φ 0011 owing to the same angular-spin wave function. It has been verified that inclusion of this excited 1s state in the calculation will not alter our main conclusion for weak two-body interaction.
In the language of Raman-induced SO coupling [1,8,14], the single particle ground state has three minima corresponding to three states Φ j=1,mj =0,±1 . At the manybody level, the interaction will determine which minimum or minima the Bose gases will condensate to by minimizing the GP energy functional [18]. We therefore set the variational wavefunction as with the normalization constraint c 2 1+ +c 2 10 +c 2 1− = 1 that ensures the conservation of the particle numbers. The coefficients c 1+ , c 10 , c 1− , to be determined by the interaction, are generally complex. For the sake of simplicity, we restrict them to be real here. Later, we will discuss the consequences of such restriction. This variational wave function ansatz is extensively encountered in SO coupled cold atom gases. [3,4,18,19,24,25,27,30,53]

B. Nematic Order
Prior to the calculation of the energy functional, we first elucidate the meaning of our variational order parameter ansatz. Since our single-particle hamiltonian respects the SO(3) R+S symmetry, it is convenient to introduce the polarization operator [60] to describe the spin order. The polarization operators for spin-s system T (l) m l (s) with l = 0, 1, . . . , 2s and m l = −l, −l + 1, . . . l, are a set of (2s + 1) 2 operators which act on the spin wave functions and transform under the coordinate system rotation according to the irreducible representation of SO(3) group. In this sense, it is also an irreducible tensor of rank l.
For spin-1 objects, the nine polarization operators T (l) m l (s) with l = 0, 1, 2 constitute a complete set of square 3 × 3 matrices and are generators of the unitary group U (3) in rotationally covariant Racah form [61]. The rank-0 operator is the unit 3 × 3 matrix I The rank-1 operators T (1) m l (s) are proportional to the irreducible rank-1 spin tensor with components where the spherical components of the irreducible rank-1 spin tensor s (1) m l are related to the cartesian ones as The rank-2 polarization operators T m l (s) are where A (m) B (n) (k) defines the rank-k tensor product of rank-m irreducible tensor A (m) and rank-n irreducible tensor B (n) . T m l (s) are equivalent to a symmetric traceless cartesian tensor of rank-2, i.e. the spin nematic tensor or quadrupole tensor N through the relation [60] The unit matrix I is a spin-rotation invariant scalar which contains important information regarding the charge (density) order, three matrices s x , s y , s z form a vector which represents the local spin (magnetic) order, and five matrices N ij form a symmetric traceless tensor which represents the local spin fluctuations or nematic order [62,63]. The spin-1 systems therefore support spin nematic order in addition to the charge and magnetic ones. The magnetic and the nematic orders compete with each other as increasing one of them requires reducing the other [62].
With the polarization operators, our variational wavefunction Eq. (12) can be written as where ζ = (c 1+ , c 10 , c 1− ) T is a normalized spinor and the position-dependent transformation on the spinor ζ leads to the spin-textured ground states in the SO coupled cold atom gases. Here the dot defines the scalar product of the modified spherical harmonics C (l) and T (l) . This local modulation operator is expanded in series of C (l) · T (l) with the highest order l = 2s. The hamiltonian is invariant under the simultaneous rotation in spin and coordinate space SO(3) R+S , therefore trapped spinor condensate with SO coupling inevitably carries angular momentum by twisting its spinor order parameter [52,62]. It is also intuitive to understand the modulation of spin-1/2 objects discussed in detail in Ref. [52], in which case the four polarization operators T where the density and spin modulations are apparently separated owing to the absence of nematic order. The spin-1 system thus enables to explore spin-tensor-related physics in the SO coupling superfluid, which has fundamentally different rotation properties as in spin-1/2 system.

C. Energy Functional
We now have six variational parameters A 0 , A 1 , A 2 , c 1+ , c 10 , c 1− in the trail variational order parameter. The energy functional of Eq. (4) can be calculated analytically using the proposed order parameter Eq. (12).
The single particle part of the energy functional consists of the kinetic energy, the trapping potential of the 3D harmonic oscillator, and the SO coupling term. As in the spin-1/2 case [52], the matrix elements for the kinetic energy and trapping potential are non-vanishing for states with the same parity, while the SO coupling term s · p will mix states with opposite parities. The result is The contribution from the d-wave states is collected in ∆ 0 (see Supplemental Material), and we have used φ 0011 |s · p| φ 0111 = −i, φ 0111 |s · p| φ 0211 = i 5 6 .
We refer to Ref. [52] for details of the integral calculation in which the method of irreducible tensor algebra is employed [60].
The calculation of the interaction is tedious but straightforward which yields and . Here ∆ n and ∆ s again denote the contributions from the d-wave states (see Supplemental Material). The energy functional per particle is simply the summation of Eqs. (21), (22) and (23).

IV. GROUND STATE PHASE DIAGRAM
The ground state phase diagram for a given coupling strength λ is obtained by minimizing the variational energy with respect to the parameters A l and x under two constraints 2 l=0 |A l | 2 = 1 and c 2 1+ + c 2 10 + c 2 1− = 1. The latter further restricts x to the region [0, 1]. Should the spin-independent interaction vanish c 0 = 0, we see from Eq. (23) that the optimized parameter x either takes value of 0 for c 2 < 0, or 1 for c 2 > 0 for negligibly small contribution ∆ s from the d-wave states. The analysis below will show that this assumption holds generally except for extremely strong interaction c 0 . The variational ansatz characterizes two quantum phases: (I) the magnetic phase with a ferromagnetic manifold ζ for c 2 < 0, as x = 0 means that (c 1+ + c 1− ) 2 = 1; (II) the polar or nematic phase with a polar manifold for c 2 > 0, as x = 1 means (c 1+ + c 1− ) 2 = 0 or 2. This is consistent with the conventional spin-1 BEC [57][58][59] where for c 2 < 0 (c 2 > 0) a ferromagnetic (polar) spinor is needed to minimize the mean-field energies. If c 0 > 0, the boundary between magnetic and nematic phases will drift a little to the positive c 2 side because the spin-independent interaction (22) which prefers x = 0, prevails in the c 2 > 0 regime over the spin-dependent one (23) which prefers x = 1. For the values of optimized parameters A l , generally, the weights of p and d waves becomes more and more important with increasing interaction c 0 , which effectively diminishes the partition of s wave. In the range of c 0 = 0 ∼ 20 we considered here, |A 0 | 2 decreases from 0.884 to 0.723, |A 1 | 2 increases from 0.110 to 0.260, while |A 2 | 2 from 0.003 to 0.014 which is always negligibly small. A typical phase diagram is shown in Fig. 2 where the phase boundary is calculated in three successive approximations: "sp" -only the lowest 0s and 0p states are considered in the approximation; "spd" -the 0d energy level is added; "spds" -the 1s energy level is added further. We can find that the boundary does not alter significantly when we include the excited s wave states φ 101mj in the variational order parameter. The boundary between the magnetic phase and the nematic phase is determined in successive approximations in which we include 0s, 0p, 0d and 1s orbits step-by-step. Density distributions and spin textures of the two phases are shown in Fig. 3 to Fig. 5 respectively.
It has been pointed out that the SO coupling manifests itself in a way that the modulation of the ferromagnetic and polar spin textures in the pseduspin space could be transferred to patterned structures in orbit space even in the ground states [18]. The reason of this modulation lies in that, in the presence of SO coupling [64] or dipolar interaction [65], the spin-dependent interaction would inevitably influence the spatial motion, which leads to rich density pattern. We discuss this in the following for the magnetic and polar phases explicitly.

A. Magnetic Phase
This phase lies in the lower part of the parameter space with the states featured as c 1+ + c 1− = ±1. The corresponding spinor ζ in Eq. (18) denotes a ferromagnetic state with magnetization along any axis in the xz plane. This asymmetry between x, z and y directions is a consequence of our simplified treatment which restricts the coefficients c 1+ , c 10 , c 1− to be real, such that the spinor ζ is unable to describe the state with magnetization along y direction. If we relax this restriction to allow complex coefficients c 1+ , c 10 , c 1− , all states with spinors ζ magnetized along any spatial direction belong to this magnetic phase, which leads to an infinitely degenerate ground states. Two typical spinors are which are longitudinally and transversely magnetized ferromagnetic states respectively. This phase spontaneously breaks the timereversal symmetry which leads to spontaneous magnetization along the xz plane. The situation is just like spin half case we have studied before [52]. The difference is that the x and z ferromagnetic states are not degenerate for the spin-half case as the interaction is not SU(2) symmetric, leaving only two-fold Kramers degeneracy originating from the time-reversal symmetry there. For spin-1 gases here, the ground states are infinitely degenerate resulted from the SU(2) symmetric interaction.
We first consider the longitudinally magnetized state ζ z = (1, 0, 0) T with the time reversal state being ζ −z = (0, 0, 1) T . We plot the ground state density distribution for the three components in Fig. 3 (a), which show clearly cylindrical symmetry. While the spin component 1 is dominantly occupied in the center, which allows the condensate to develop a longitudinal magnetization, the components 0 and −1 form two toruses surrounding the central part. It can be seen that the outermost shell of spin −1 density is negligibly small which is entirely attributed to the involvement of the d wave fraction in the order parameter. To visualize the d wave nature we need to zoom out 100 times in the density plot. Generally only less populated spin component can develop d-wave characters in the ground states because these complex structure in the high density spin components would cost too much kinetic energy [66]. For its time-reversal degenerate state ζ −z the spin components 1 and −1 are inversely populated.
The spin texture of the longitudinally magnetized state is plotted in Fig. 3 (b) where we find a synchronous modulation between the particle density and the spin density owing to the breaking of spin rotation symmetry. In the trap center, the spins are aligned along the z axis due to the dominant occupation of spin-1 component. The successive population of the 0 and −1 components forms a local spin texture S(r)/n(r) where the spin density vectors deflect continuously to the xy plane away from the center. This is a magnetic skyrmion-like texture similar to the half quantum vortex configuration in the Weyl SO coupled psedo-spin-1/2 BEC [25,46,52]. The spin vector field lines develop into a bundle of fountain-like streamlines close to the z axis around which a torus is formed near the xy plane [52].
The density distribution for the transversely magne- with the time reversal state looks quite differently as shown in Fig. 4 (a). With half of the atoms filled in the spin-0 component, the components ±1 becomes equally populated and spatially separated, which leads to a significant transverse magnetization [66]. The density distributions of three components are separated in an alternative way, i.e., the −1 (+1) component lies mainly in the −y (+y) half-space and its peak density center is along the direction joining the III and VIII octants (I and VI octants) [52], while the 0 component is embedded between them with a density profile along the x axis. This leads to a magnetic skyrmion-like texture with the spins in the trap center transversely aligned along the x axis as shown in Fig. 4 (b), and the increasing population of ±1 components away from the trap center makes the spin density vector forms a torus near the yz plane. Owing to the non-commutative nature of position-dependent transformation (19) and the spin rotation, the spin texture of the longitudinally magnetized state is different from the π/2 spin rotation around y axis of the transversely magnetized state, though the spinor wave functions themselves ζ z and ζ x are related by such a rotation. All other states degenerate with these two magnetic states have similar properties, except that the magnetization axes may be in any direction determined by the values of parameters c's.
The time-reversal symmetry is preserved in this phase, so the ground state has no spontaneous magnetization as the non-magnetic phase in the Raman induced SO coupled two-component Bose gases [1,4,8,14]. This is a new phase and no analogy in spin half case we considered before [52]. All states belong to this phase are again infinitely degenerate. The density distribution of the state ζ p1 is plotted in Fig. 5 (a). A signature of the d-wave is seen in the 100 times zoomed density of 0 component in the xz plane, which comes from the d-wave contribution in the variational order parameter. The phase separation among three components, which is possible only in the case of time-reversal symmetry breaking [67], i.e. the magnetic phase, is not seen for this nematic phase which preserves the time-reversal symmetry and the ±1 components are miscible.
We define the nematicity density tensor N = Ψ † N Ψ to characterize the nematic order due to the absence of the magnetization. In the nematic phase ζ p1 or ζ p2 , the normalized nematicity density tensor N/n has eigenvalues { 1 3 , 1 3 , − 2 3 } everywhere, therefore describing a uniaxial nematic state. The eigenvector associated with eigenvalue − 2 3 defines the nematic director which is plotted in Fig. 5 (b) together with the tensor magnetization density N zz = Ψ † N zz Ψ. The nematic directors form a lantern-like structure with the principle axis along the x direction. The spatially modulation of the nematic directors, shown as headless vectors in the xy, yz and xz planes, reflects indirectly the modulation of nematicity density tensor themselves. In the trap center, where the spinor wavefunction ζ p1 describes a transverse polar state, the nematic director points along the x axis. Away from the trap center, the 0 component is gradually populated and the nematic director are continuously modulated into concentric circles in the yz plane. For the state ζ p2 , the 0 component in the trap center will allow a longitudinal polar state and we find an alternative modulation of the nematic directors. The tensor magnetization N zz [17,20] has been adopted to resolve the order of the phase transition in the Raman-induced SO coupled spin-1 condensate, and the spatial modulation of N zz along direction of the SO coupling has been noticed [19,21].
The magnetic phase is also featured with a nematic order and the competing orders in both phases are shown in Fig. 6. These two spin orders are competing with each other [62] to meet the requirement This prevents us from writing the transformation Eq. (19) into a local spin rotation as Eq. (20). Diagonalizing the nematicity density tensor yields three distinct eigenvalues which are spatially modulated as well as the nematicity density tensor themselves [21]. Thus the magnetic phase has both magnetic and biaxial nematic orders, while the nematic phase exhibits only a uniaxial nematic order. Finally, we remind that our variational order parameter based on the perturbation expansion may not be ap-plicable in the limit of strong SO coupling [28,43,44,48] where a skyrmion-lattice-like ground state may appear. Moreover, when two-body interaction is strong enough, the variational order parameter starts to involve higher angular momentum j states which will break the SO(3) rotational symmetry, and the Bose gases will condense into the plane wave or stripe phases instead [26,44,53].

V. SUMMARY
We establish the ground state phase diagram of the weakly 3D spin-orbit coupled spin-1 bosons theoretically. The ground state may be in a magnetic or a nematic phase determined by the competing between the spinindependent and the spin-dependent two-body interaction. The nematic phase is a new phase that is absent in a 3D spin-orbit coupled pseudo-spin-1/2 bosonic system. We discuss the density distribution and spin orders of the two phases in detail. The magnetic phase permits both a magnetic order and a biaxial nematic order, while the nematic phase is featured by a uniaxial nematic order. These novel phases are in current experimental reach benefiting from the rapid progress of cold gases with artificial gauge field.