Half-quantum vortex state in a spin-orbit coupled Bose-Einstein condensate

2.8 MB·

by  B. Ramachandhran

Half-quantum vortex state in a spin-orbit coupled Bose-Einstein condensate
B. Ramachandhran1, Bogdan Opanchuk2, Xia-Ji Liu2, Han Pu1, Peter D. Drummond2, and Hui Hu2
1Department of Physics and Astronomy, and Rice Quantum Institute, Rice University, Houston, TX 77251, USA
2ARC Centres of Excellence for Quantum-Atom Optics and Centre for Atom Optics and Ultrafast Spectroscopy,
Swinburne University of Technology, Melbourne 3122, Australia
(Dated: January 9, 2012)
We investigate theoretically the condensate state and collective excitations of a two-component
Bosegasintwo-dimensionalharmonictrapssubjecttoisotropicRashbaspin-orbitcoupling. Inthe
weaklyinteractingregimewhentheinter-speciesinteractionislargerthantheintra-speciesinterac-
tion(g >g),wefindthatthecondensategroundstatehasahalf-quantum-angular-momentumvor-
↑↓
2 texconfigurationwithspatialrotationalsymmetryandskyrmion-typespintexture. Uponincreasing
1 the interatomic interaction beyond a threshold gc, the ground state starts to involve higher-order
0 angular momentum components and thus breaks the rotational symmetry. In the case of g < g,
↑↓
2 thecondensatebecomesunstabletowardsthesuperpositionoftwodegeneratehalf-quantumvortex
states. Both instabilities (at g > g and g < g) can be determined by solving the Bogoliubov
n c ↑↓
equations for collective density oscillations of the half-quantum vortex state, and by analyzing the
a
J softening of mode frequencies. We present the phase diagram as functions of the interatomic inter-
actions and the spin-orbit coupling. In addition, we directly simulate the time-dependent Gross-
6
Pitaevskii equation to examine the dynamical properties of the system. Finally, we investigate the

stability of the half-quantum vortex state against both the trap anisotropy and anisotropy in the
]
s spin-orbit coupling term.
a
g
PACSnumbers: 05.30.Jp,03.75.Mn,67.85.Fg,67.85.Jk
-
t
n
a I. INTRODUCTION Bose gas into a half-quantum vortex configuration was
u
first suggested by Congjun Wu and co-workers in 2008
q
. Owing to the unprecedented control in interatomic in- and its existence was discussed under the condition that
at teraction, geometry and purity, atomic quantum gases the interaction is SU(2) symmetric, i.e., g = g↑↓ [13].
m have proven to be an ideal many-body platform for Here, we explore systematically the parameter space for
- exploring fundamental quantum states, such as Bose- the half-quantum vortex state and analyze its stability.
d Einsteincondensates(BEC)[1],stronglyinteractinguni- Wepresentaphasediagramforthehalf-quantumvortex
n tary Fermi superfluids [2, 3] and Mott-insulating states stateasfunctionsoftheSOcouplingandtheinteratomic
o [4]. One of the latest achievement concerns the spin- interaction strengths. We also investigate the dynamical
c
[ orbit (SO) coupling in an ultracold spinor Bose gas of properties of the half-quantum vortex state by directly
87Rb atoms [5], induced by the so-called “synthetic non- simulating the time-dependent GPE. Finally, the stabil-
1 Abelian gauge fields”. Novel quantum states may be an- ityofthehalf-quantumvortexstateagainstboththetrap
v ticipated in the presence of SO coupling [6–18]. Indeed, anisotropyandanisotropyinthespin-orbitcouplingterm
1
for a homogeneous SO coupled spin-1/2 Bose gas with is examined.
7
4 intra- and inter-species interactions (g and g ), a sin- Our main results are summarized in Fig. 1. The
↑↓
1 gle plane-wave or a density-stripe condensate state has half-quantum vortex state (the phase I) is the ground
. been predicted [8], depending on whether g is smaller or state if the intra-species interaction is smaller than the
1
0 larger than g↑↓. Interesting density patterns have been inter-species interaction (g < g↑↓) and if the interaction
2 observed in the theoretical simulations for an SO cou- strength is below a threshold (g < g ). Otherwise, it
c
1 pled spinor condensate, in the absence [8, 11, 12, 17, 18] becomes energetically unstable towards a superposition
: or presence [14–16] of rotation. The phenomenon of self- state of two degenerate half-quantum vortex states (the
v
i trapped BECs has also been proposed, in particular, in phaseIIA),orastateinvolvinghigher-orderangularmo-
X one-dimensional (1D) geometry [7]. mentum components (the phase IIB). With decreasing
ar In this work, we show that in a Rashba SO coupled, the dimensionless SO coupling strength λSO, the thresh-
weakly interacting spin-1/2 Bose gas in two-dimensional old g becomes exponentially large, leading to a large
c
(2D) harmonic traps, all bosons may condense into a parameter space for the half-quantum vortex state (see
non-trivial half-integer angular momentum state (or a Fig. 14). Itisthereforefeasibletobeobservedinthecur-
half-quantum vortex state) with a skyrmion-type spin rent experiments with ultracold SO coupled spinor Bose
texture. We solve the mean-field Gross-Pitaevskii equa- gases of 87Rb atoms.
tion(GPE)foritsdensitydistributionsandspintextures, The rest of the paper is organized as follows. In the
and obtain its collective excitation spectrum by solving next section, we outline the model Hamiltonian and dis-
theBogoliubovequationand bydirectlysimulatingreal- cuss briefly the existence of half-quantum vortex state
time propagation of the GPE ground state under per- in the non-interacting limit. In Sec. III, we present the
turbation. The condensation of an SO coupled spin-1/2 numerical procedure of solving the GPE and Bogoliubov
2

remarks.
2.0
I
1.5 IIB II. THEORETICAL FRAMEWORK
g
/

g1.0 We consider a two-component Bose gas confined in
a 2D isotropic harmonic trap V(ρ) = Mω2(x2 +
⊥
y2)/2 = Mω2ρ2/2 with a Rashba SO coupling V =
0.5 IIA ⊥ SO
−iλ (σˆ ∂ − σˆ ∂ ), where λ is the Rashba SO cou-
(a) =1 R x y y x R
SO pling strength and σˆx, σˆy, and σˆz are´the 2×2 Pauli
0.0
0 10 20 30 40 50 matrices. ThemodelHamiltonianH= dr[H0+Hint]is
given by,
2.0
1.5 I IIB H0 = Ψ†(cid:20)−(cid:126)22M∇2 +V (ρ)+VSO−µ(cid:21)Ψ, (1)
g
/ H = (g/2) (cid:88) Ψ†Ψ†Ψ Ψ +g Ψ†Ψ Ψ†Ψ , (2)
g 1.0 int σ σ σ σ ↑↓ ↑ ↑ ↓ ↓
σ=↑,↓

0.5 IIA where r = (x,y) and Ψ = [Ψ↑(r),Ψ↓(r)]T denotes the
(b) =4 spinor Bose field operators in a collective way, and the
SO
0.0 chemical potential µ is to´be determined by the total
0.0 0.2 0.4 0.6 0.8 1.0 number of bosons N, i.e., drΨ†Ψ=N. For simplicity,
2
g(N-1)a /( ) wehaveassumedequalintra-speciesinteractionstrength
g = g = g. In experiments, the two-dimensionality
Figure 1: (color online). Phase diagram at two dimension- ↑↑ ↓↓
can be readily realized by imposing a strong harmonic
less SO coupling strengths, λ = 1 (a) and λ = 4 (b).
SO SO
The half-quantum vortex state (the phase I) becomes un- potential V(z) = Mωz2z2/2 along axial direction, in
stable when the intra-species interaction is larger than the such a way that µ,kBT (cid:28) (cid:126)ωz [19]. For the realis-
inter-species interaction (g > g , the phase IIA) or when tic case of 87Rb atoms, the interaction strengths can
↑↓
the interatomic interactions are sufficient strong (g >g , the be calculated from the two s-wave scattering lengths
c √
phaseIIB).Theinsetsshowsthedensitypatternsofthespin- a (cid:39) 100a and a , using g = 8π((cid:126)2/M)(a/a ) and
√B ↑↓ z
up and spin-down bosons in the phases I and IIA. We note g = 8π((cid:126)2/M)(a /a ), respectively. Here a =
that,thecriticalinteractionstrengthgcincreasesrapidlywith (cid:112)↑↓(cid:126)/(Mω ) is the cha↑r↓actzeristic oscillator length inz z-
decreasing the SO coupling strength λ . z
SO
direction.
For a weakly interacting Bose gas at zero tempera-
equations and discuss the typical density distributions ture, we assume that all the bosons condense into a
andcollectivemodebehaviorsofthehalf-quantumvortex single quantum state Φ(r)=[Φ↑(r),Φ↓(r)]T. Following
state. The collective excitation spectrum obtained from the standard mean-field theory [20], we separate the
the Bogoliubov equation is compared to a direct simula- field operator into a condensate and a fluctuation part,
tion of the time-dependent GPE. In Sec. IV, we analyze Ψσ(r)=Φσ(r)+Ψ˜σ(r). Keeping up to the´ quadratic
the stability of the half-quantum vortex state by mon- terms inΨ˜ (r), thisseparation leadsto H= dr[H +
σ GP
itoring the softening of collective mode frequencies and H ], where the condensate part is given by,
T
by comparing the energy with that of some competing
states. The phase diagram is then constructed as func- H = Φ†[H +V −µ]Φ
GP osc SO
tions of interatomic interactions and SO coupling. The g (cid:16) (cid:17)
+ |Φ |4+|Φ |4 +g |Φ Φ |2, (3)
stabilityagainsttheanisotropyintrappingpotentialand 2 ↑ ↓ ↑↓ ↑ ↓
inspin-orbitcouplingtermisalsocarefullyexamined. Fi-
nally, we summarize in Sec. V and give some concluding and the fluctuation part H =Ψ˜†H Ψ˜ with
T Bog
 H +g|Φ |2 V +g Φ Φ∗ gΦ2 g Φ Φ 
s↑ ↑ so ↑↓ ↑ ↓ ↑ ↑↓ ↑ ↓
V† +g Φ∗Φ H +g|Φ |2 g Φ Φ gΦ2 
 so ↑↓ ↑ ↓ s↓ ↓ ↑↓ ↑ ↓ ↓ 
HBog = g(cid:16)Φ∗(cid:17)2 g Φ∗Φ∗ H +g|Φ |2 −V† +g Φ∗Φ . (4)
 ↑ ↑↓ ↑ ↓ s↑ ↑ so ↑↓ ↑ ↓ 
 (cid:16) (cid:17)2 
g Φ∗Φ∗ g Φ∗ −V +g Φ Φ∗ H +g|Φ |2
↑↓ ↑ ↓ ↓ so ↑↓ ↑ ↓ s↓ ↓
3
HereH ≡−(cid:126)2∇2/(2M)+V (ρ),H ≡H +g|Φ |2+ In the absence of interatomic interactions, the single-
osc s↑ osc ↑
Vg↑↓|Φ≡↓|2−−iλµ (a∂nd+Hsi↓∂≡) Hanodsc +V†g↑↓≡|Φ↑−|2iλ+g(∂|Φ↓−|2i−∂µ),, pgiavretniclbey,wave-function [φ↑(r),φ↓(r)]T with energy (cid:15) is
so R y x so R y x
and we have introduced a 4 × 4 Nambu spinor Ψ˜ = (cid:20) (cid:21)(cid:20) (cid:21) (cid:20) (cid:21)
H −iλ (∂ +i∂ ) φ φ
[Ψ˜↑(r),Ψ˜↓(r),Ψ˜†↑(r),Ψ˜†↓(r)]T. −iλR(∂oysc−i∂x) RHoysc x φ↑↓ =(cid:15) φ↑↓ .
The condensate wave-function can be obtained from
(7)
the GP equations δH /δΦ(r)=0 [20], or explicitly,
GP In polar coordinates (ρ,ϕ), we have −i(∂ ± i∂ ) =
y x
(cid:20) (cid:21)(cid:20) (cid:21) e∓iϕ[±∂/∂ρ−(i/ρ)∂/∂ϕ]. Because of the isotropic har-
H −iλ (∂ +i∂ ) Φ (r)
s↑ R y x ↑ =0. monic potential V (ρ), the single-particle wave-function
−iλ (∂ −i∂ ) H Φ (r)
R y x s↓ ↓ may have a well-defined azimuthal angular momentum
(5) l =m and may take the form,
z
Atzerotemperature,weassumeasinglecondensatestate
with zero quantum depletion, ´so that the condensate (cid:20) φ (ρ) (cid:21)eimϕ
wave-function is normalized by dr[|Φ↑|2+|Φ↓|2] = N, φm(r)= φ↓(↑ρ)eiϕ √2π. (8)
where N is the total number of bosons. The equation
becomes simplified if we write Φ = N1/2φ and Φ = This state also has a well-defined total angular momen-
↑ ↑ ↓
N1/2φ↓ and use accordingly the interaction strengths tumjz =lz+sz =m+1/2. Ingeneral,wemaydenotethe
g(N −1)´and g↑↓(N −1). The normalization condition energyspectrumas(cid:15)nm,wheren=(0,1,2...)isthequan-
becomes dr[|φ |2+|φ |2]=1. tum number for the transverse (radial) direction. There
↑ ↓
The quasi-particle wave-functions with energy (cid:126)ω sat- is an interesting two-fold degeneracy of the energy spec-
isfy the Bogoliubov equations [20], trum: any eigenstate φ(r) = [φ↑(r),φ↓(r)]T is degener-
ate with its time-reversal partner Tφ(r)≡(iσ C)φ(r)=
y
u↑(r) +u↑(r) [φ∗↓(r),−φ∗↑(r)]T. Here C is the complex conjugate oper-
u (r) +u (r) ation. This Kramer doublet is the direct consequence of
HBog v↓(r) =(cid:126)ω −v↓(r) , (6) thetime-reversalsymmetrysatisfiedbythemodelHamil-
↑ ↑
v (r) −v (r) tonian. Itpreservesaswellinthepresenceofinteratomic
↓ ↓
interactions. As a result, we may restrict the quantum
´
andisnormalizedby dr[|u |2+|u |2−|v |2−|v |2]=1. number m to be non-negative integers, as a negative m
↑ ↓ ↑ ↓
can always be regarded as the time-reversal partner for
These Bogoliubov quasi-particles correspond to the dif-
a state with m≥0.
ferentcollectivedensityoscillationmodesaroundthecon-
To solve numerically the single-particle spectrum, we
densate with the frequency ω [21]. It is easy to see that
the wave-function [v∗(r),v∗(r),u∗(r),u∗(r)]T is also a adoptabasis-expansionmethod. Tothisend,weexpand
↑ ↓ ↑ ↓ first,
solution of Eq. (6), but with energy −(cid:126)ω. This is antic-
ipated for the usual Bogoliubov transformation. Phys- (cid:88)
φ (ρ) = A R (ρ), (9)
ically, we should restrict to a non-negative mode fre- ↑ k km
quency, ω ≥0. k
(cid:88)
In harmonic traps, it is natural to use the trap units, φ (ρ) = B R (ρ), (10)
↓ k km+1
i.e. to take (cid:126)ω⊥ as the unit for energy and the harmonic k
oscillatorlengtha =(cid:112)(cid:126)/(Mω )astheunitforlength.
⊥ ⊥
This is equivalent to set (cid:126)=k =M =ω =1. For the where
B ⊥
Sc(cid:126)o2Ou/(pcMloinuλgpRlis)ntgrae,nndwgetchoinnλtsSreoOqduu=ecnetalay⊥n/daSeλfiOn=ecoa(cid:112)upd(liMimng/en(cid:126)les3in)oλgnRtlhe/s√asλωS⊥=O. Rkm = a1⊥(cid:115)(k+2k|m! |)!(cid:18)aρ⊥(cid:19)|m|e−2ρa22⊥Lk|m|(aρ2⊥2) (11)
InanSOcoupledspin-1/2BECof87Rbatomsasrealized
is the radial wave-function of a 2D harmonic oscillator
recently by the NIST group [5], λ is about 10. In
SO H withenergy(2k+|m|+1)(cid:126)ω ,andL|m| istheasso-
the typical experiment for 2D spin-1/2 87Rb BECs [19], osc ⊥ k
ciatedLegendrepolynomial. Then,wehavethefollowing
the interatomic interaction strengths are about g(N −
1) ≈ g (N −1) = 102 ∼ 103(cid:126)ω /a2. These coupling secular matrix,
↑↓ ⊥ ⊥
strengths, however, can be precisely tuned by properly (cid:20) (cid:21)(cid:20) (cid:21) (cid:20) (cid:21)
H MT A A
choosing the parameters of the laser fields that lead to osc↑ k =(cid:15) k , (12)
M H B B
the harmonic confinement and the SO coupling. osc↓ k k
where the matrix elements are given by (for m≥0)
A. Single-particle solutions Hosc↑,kk(cid:48) = (cid:126)ω⊥[2k+m+1]δkk(cid:48),
Hosc↓,kk(cid:48) = (cid:126)ω⊥[2k+(m+1)+1]δkk(cid:48),
The appearance of the half-quantum vortex state may (cid:104)√ √ (cid:105)
be easily understood in the non-interacting limit [13]. Mkk(cid:48) = (cid:126)ω⊥λSO k(cid:48)+m+1δkk(cid:48) + k(cid:48)δkk(cid:48)−1 .
4
Diagonalization of the secular matrix Eq. (12) leads Here α and β are two arbitrary complex numbers satis-
to the single-particle spectrum and single-particle wave- fying |α|2+|β|2 =1.
functions. In numerical calculations, it is necessary to In the presence of very weak interatomic interactions
impose a cut-off kmax for the radial quantum number k such that g(N −1)a2⊥,g↑↓(N −1)a2⊥ (cid:28)∆(cid:15), where ∆(cid:15) is
ofthe2Dharmonicoscillator. Forλ ≤20,wefindthat the energy difference between the single-particle ground
SO
kmax = 256 is already sufficiently large to have an accu- state φ0(r) and the first excited state φ1(r), we may de-
rate energy spectrum. With this cut-off, the dimension termine the superposition coefficients α´and β by min-
of the secular matrix in Eq. (12) is 2kmax =512. imizing the GP energy, EGP[φs(r)] = drHGP[φs(r)].
After a simple algebra, we find that,

8
∆E = E [φ (r)]−E [φ (r)], (14)
) GP s GP 0
6 = (g −g)(N −1)|αβ|2W[φ (r)], (15)
↑↓ 0
( 4
/m
n where the W-function is given by,
2
(a) =1 ˆ
SO
0-10 -8 -6 -4 -2 0 2 4 6 8 W[φ(r)]= dr[(|φ↑|2−|φ↓|2)2−2φ2↑φ2↓]. (16)
m
e
ofil 2 (b) SO=1, m=0 0.2 (c) m=0 Therefore, a half-quantum vortex state is preferable if
ensity pr 1 nntot W0.1 (t|αigo↑|n↓=−o|fβgt|)wW=o 1d>/e√g0e.2neOwrtialhtleebrhewaitslhfe-e,quagnaronuetqunumdaslv-twaotreetige.xhAtsstsaustpheoesrwpwnoistiihn-
D n
Fig. 2c,theW-functionforφ (r)ispositiveforarbitrary
0
0 0.0 SOcoupling. Wethusconcludethatahalf-quantumvor-
0 1 2 3 0 1 2 3 4
/a tex state should appear at weak interatomic interactions
SO
provided that the inter-species interaction is larger than
Figure 2: (color online). (a) Single-particle energy spectrum
the intra-species interaction (g >g).
at λ = 1. (b) The density profiles for the single-particle ↑↓
SO
statewithm=0atλ =1. (c)TheW-functionforthem=
SO
0 single-particle state as a function of SO coupling strength.
It is always positive at arbitrary SO coupling strength. III. DENSITY DISTRIBUTIONS AND
COLLECTIVE EXCITATIONS
InFig. 2a,weshowthesingle-particleenergyspectrum
at λSO = 1. For arbitrary SO interaction strength, we Let us now consider finite interatomic interactions,
find numerically that the doublet single-particle ground by solving the GPE for density distributions and spin-
state always occurs at m = 0 (or m = −1 for its time- textures, and the Bogoliubov equation for the collective
reversal partner state). density excitations.
B. Appearance of the half-quantum vortex state
A. GPE solutions of the half-quantum vortex state
The single-particle state with m = 0, φ (r) =
√ 0
For the half-quantum vortex condensate state with
[φ (ρ),φ (ρ)eiϕ]T/ 2π, has a half-quantum vortex con-
↑ ↓ m = 0, the GP equation becomes L [φ (ρ),φ (ρ)] =
figuration [13, 22], as the spin-up component stays in GP ↑ ↓
0, where
the s-state while the spin-down component in the p-
state and the resulting spin texture is of skyrmion type
(cid:20)H +g¯φ2+g¯ φ2 λ (∂ +1/ρ) (cid:21)
(see Fig. 2b for density distributions and Sec. IIIB for L = s,0 ↑ ↑↓ ↓ R ρ ,
more discussions on spin-texture). In the absence of in- GP λR(−∂ρ) Hs,1+g¯↑↓φ2↑+g¯φ2↓
teractions, however, there is a degenerate time-reversal (17)
√
state,Tφ0(r)=[φ↓(ρ)e−iϕ,−φ↑(ρ)]T/ 2π,whichisalso g¯ ≡ g(N − 1)/(2π) and g¯↑↓ ≡ g↑↓(N − 1)/(2π), and
a half-quantum vortex state. Therefore, in general, the Hs,m ≡−[(cid:126)2/(2M)][∂2/∂ρ2+(1/ρ)∂ρ−m2/ρ2]+V(ρ)−µ.
ground single-particle state is a superposition of two de- The numerical procedure for solving GPE is very similar
generatehalf-quantumvortexstatesofφ (r)andTφ (r), to that for single-particle states in Eq. (12). We expand
0 0 (cid:80) (cid:80)
which takes the form φs(r) = αφ0(r)+βTφ0(r), or ex- φ↑(ρ) = kAkRk0(ρ) and φ↓(ρ) = kBkRk1(ρ), and
plicitly, obtain the secular matrix (with m=0),
1 (cid:20)αφ (ρ)+βφ (ρ)e−iϕ (cid:21) (cid:20)H +I MT (cid:21)(cid:20) A (cid:21) (cid:20) A (cid:21)
φ (r)= √ ↑ ↓ . (13) osc↑ ↑ k =µ k , (18)
s 2π αφ↓(ρ)eiϕ−βφ↑(ρ) M Hosc↓+I↓ Bk Bk
5
where that the density distributions are flattened significantly
ˆ by interatomic interactions, as anticipated. The 2D con-
∞
I↑,kk(cid:48) = ρdρRk0(ρ)(cid:0)g¯φ2↑+g¯↑↓φ2↓(cid:1)Rk(cid:48)0(ρ),(19) tour plot of the spin-up and spin-down density patterns
of the half-quantum vortex state is shown in the inset of
ˆ0
∞ Fig. 1 (in the phase I).
I↓,kk(cid:48) = ρdρRk1(ρ)(cid:0)g¯↑↓φ2↑+g¯φ2↓(cid:1)Rk(cid:48)1(ρ).(20) Togainmoreinsightsofthehalf-quantumvortexstate,
0
it is useful to calculate the spin vector
The chemical potential is given by the lowest eigenvalue
of the secular matrix. Due to the non-linear terms of 1Φ†σΦ
S(r)= (21)
I↑,kk(cid:48) andI↑,kk(cid:48),wehavetoupdatethecondensatewave- 2 |Φ|2
functionsanddensitiesiteratively. Toovercomethelarge
non-linearity,weuseasimplemixingschemebysettinga and the skyrmion density
smallparameter0<γ <1andreplacethepreviousden-
sityφ2 by(1−γ)φ2 +γφ2,whereφ2 isthedensity 8
σ,old σ,old σ σ n (r)= S·[∂ S×∂ S]. (22)
calculated in the current step [23]. The choice of γ de- skyrmion 4π x y
pends on the interaction strengths. It becomes smaller
The skyrmion density is a measure of the winding of the
for larger g¯ and g¯ . We run the iteration until conver-
↑↓ spin profile. If it integrates to 1 or −1, a topological
genceisachievedwithinasettolerance. Wehavechecked
stable knot exists in the spin texture.
that this procedure of solving GPE is stable for interac-
tion strengths up to g(N−1),g (N−1)<103(cid:126)ω /a2.
↑↓ ⊥ ⊥
For even larger non-linearity, it seems to be impractical 4 4 4 -0.5
to expand the condensate wave-function using the 2D (a) <Sx> (b) <Sy> (c) <Sz>
2 2 2
harmonic oscillator basis. Therefore for large interac-
a 0
y/0 0 0
tion strengths, we use a time-splitting spectral method
(TSSP)techniquetosolvethecoupledGPequationsand -2 -2 -2
0.5
obtain the ground state by imaginary-time propagation
-4 -4 -4
-4 -2 0 2 4 -4 -2 0 2 4 -4 -2 0 2 4
[24,25]. Forsmallinteractionstrengths,resultsobtained x/a x/a x/a
fromTSSPareidenticaltothoseobtainedfromthebasis-
Figure 4: (color online). Contour plots of the three compo-
expansion method. nents of spin vector S(r) at λ =1, g(N −1)=40(cid:126)ω /a2
SO ⊥ ⊥
and g /g=1.1.
↑↓
B. Density distributions and spin textures

4 4 -0.05
(a) spin texture (b) skyrmion density
0.6 6 2 2
0.5 (a) SO=1 5 (b) SO=4 y/a0 y/a0
ofile0.4 ntot ofile4 ntot -2 -2 0.1
pr n pr n
y y -4 -4
ensit0.3 n ensit 3 n -4 -2 x/0 a 2 4 -4 -2 x/0a 2 4
D D
0.2 2 Figure 5: (color online) (a) Two-dimensional vector plot
of the transverse spin vector (S S ) at λ = 1, g(N −
x, y SO
0.1 1 1) = 40(cid:126)ω⊥/a2⊥ and g↑↓/g = 1.1. The color and length
of arrows give respectively the orientation and the magni-
tude of (S S ). (b) The corresponding skyrmion density
0.0 0 x, y
0 1 2 3 4 0.0 0.5 1.0 1.5 2.0 n (r).
skyrmion
/ a /a
Figure3: (coloronline). DensitydistributionsatλSO =1and In Fig. 4, we report the three components of the
g(cid:126)(ωN⊥/−a2⊥1)(b=).4H0(cid:126)eωre⊥,/tah2⊥e r(aat)ioangd↑↓/agt λ=S1O.1=. 4 and g(N −1) = gsp↑↓in/gve=cto1r.1.atTλhSeOtr=an1sv,egrs(eNsp−in1)te=xtu4r0e(cid:126)iωs⊥s/hao2⊥wnanind
Fig. 5a by arrows, with color and length representing
InFig.3,wepresenttheradialdensitydistributionsof the orientation and the magnitude of the transverse spin
thehalf-quantumvortexcondensatestateattwoSOcou- vector (S S ), respectively. It is readily seen that the
x, y
pling strengths: λ = 1 and λ = 4. The increase of spin vector spirals in space and form a skyrmion-type
SO SO
theSOcouplingleadstomoreoscillationsintheradialdi- texture. Quantitatively,thisismostclearlyillustratedin
rection. BycomparingFig.3(a)withFig.2(b),onefinds Fig. 5b, where we plot the skyrmion density.
6
C. Solutions of Bogoliubov equations We should take the positive branch. We note also that
the Bogoliubov quasi-particles at a negative azimuthal
Given the wave-function of the half-quantum vor- quantum number m may be obtained from the negative
√
tex state, [φ (ρ),φ (ρ)eiϕ]T/ 2π, we now turn to branchofthesolutionwithm>0,becauseoftheduality.
↑ ↓
consider its collective excitations, as described by
the coupled Bogoliubov equations (6). As a re-
1. Breathing modes
sult of rotational symmetry, it is easy to see
that, the Bogoliubov wave-functions have a good az-
imuthal quantum number m and can be√written as, In the case of the breathing mode (m=0), where
[u↑(ρ),u↓(ρ)eiϕ,v↑(ρ),v↓(ρ)e−iϕ]Teimϕ/ 2π. There- (cid:20)L +U U (cid:21)
fore, we have H = GP , (31)
Bog U L +U
GP
u (ρ) +u (ρ)
↑ ↑ we may have an alternative way to solve the Bogoliubov
u (ρ) +u (ρ)
H  ↓ =(cid:126)ω ↓ , (23) equation, following Hutchinson, Zaremba, and Griffin
Bog v (ρ)   −v (ρ) 
↑ ↑ (HZG) [26]. By denoting collectively u = [u (ρ),u (ρ)]
v (ρ) −v (ρ) ↑ ↓
↓ ↓ and v =[v (ρ),v (ρ)], we have,
↑ ↓
where (L +2U)(u+v) = (cid:126)ω(u−v), (32)
GP
(cid:20)L +U U (cid:21) L (u−v) = (cid:126)ω(u+v). (33)
H = m , (24) GP
Bog U L +U
−m
Let us now expand the wave-functions u ± v in terms
with of the eigenfunctions ψα of LGP with energy (cid:15)α (i.e.,
L ψ =(cid:15) ψ ),
(cid:20)H +g¯φ2+g¯ φ2 λ [∂ +(m+1)/ρ] (cid:21) GP α α α
Lm = λs,Rm(−∂ρ↑+m/↑ρ↓) ↓ HsR,m+ρ1+g¯↑↓φ2↑+g¯φ2↓(25,) u−v = α(cid:88)(cid:54)=0(cid:15)c1αα/2ψα, (34)
and
U =(cid:20)g¯ g¯φφ2↑φ g¯↑↓g¯φφ↑2φ↓ (cid:21). (26) u+v = α(cid:88)(cid:54)=0(cid:15)1α(cid:126)/2ωcαψα. (35)
↑↓ ↑ ↓ ↓
Here, the lowest eigenstate of L with zero energy
GP
TosolvetheBogoliubovequation,asbeforeweexpand should be removed, as it corresponds exactly to the con-
the wave-functions using 2D harmonic oscillator basis, densatemode. Itiseasytoseethat(L +2U)L (u−
GP GP
v)=((cid:126)ω)2(u−v)andL (L +2U)(u+v)=((cid:126)ω)2(u+
(cid:88) GP GP
u↑(ρ) = akRkm(ρ), (27) v). Inserting the expansion of u−v or u+v, one finds
k the secular equation,
(cid:88)
u↓(ρ) = bkRkm+1(ρ), (28) (cid:88)(cid:110)(cid:15)2δ +2(cid:15)1/2U (cid:15)1/2(cid:111)c =((cid:126)ω)2c , (36)
k α αβ α αβ β β α
(cid:88) β
v (ρ) = c R (ρ), (29)
↑ k km
where
k
ˆ
(cid:88)
v (ρ) = d R (ρ). (30) ∞
↓ k km−1 U = ρdρψ† (ρ)Uψ (ρ). (37)
αβ α β
k
0
This leads to a secular matrix of H , whose ele- By diagonalizing the secular matrix, one obtains the
Bog
ments can be calculated directly using the 2D har- mode frequency ω and the coefficients c . The latter
α
monic oscillator basis. We note that, to obtain the Bo- should be normalized as (cid:80) c2 =(cid:126)ω, in accord with the
α α
goliubov quasiparticles we cannot diagonalize directly normalization condition for u and v.
the secular matrix, because of the minus sign before We have checked numerically that the HZG so-
v↑(ρ) and v↓(ρ) at the right-hand side of Eq. (23). lution leads to exactly the same result as the di-
Instead, we should diagonalize a non-symmetric ma- rect diagonalization of the non-symmetric matrix
trixDiag{+1,+1,−1,−1}HBog andno´rmalizethequasi- Diag{+1,+1,−1,−1}HBog, if we discard the zero-
particle wave-functions according to ∞ρdρ[u2 + u2 − frequency condensate mode in the latter method.
0 ↑ ↓
v2 − v2] = 1. The number of resulting eigenvalues is
↑ ↓
two times the number that we want. There are two
branches of eigenvalues, one is positive and the other D. Collective excitations
negative, as a result of the duality between the so-
lution [u (r),u (r),v (r),v (r)]T (with energy +(cid:126)ω) In Fig. 6, we report the breathing (m = 0) and the
↑ ↓ ↑ ↓
and [v∗(r),v∗(r),u∗(r),u∗(r)]T (with energy −(cid:126)ω). dipole mode (m = ±1) frequencies as a function of the
↑ ↓ ↑ ↓
7

1.0
(a) m=-1 (b) m=0 (c) m=1
4 (A) n = 0 u( ) (C) n = 2
0.5 u( )
3 ns
o0.0
/ cti
2 un-0.5 v( )
f
-
ve v( )
1 wa1.0
e
d0.5
0 o
0 1g0(N-1)2a0 2/( 3)0 0 1g0(N-1)2a02 /( 3)0 0 10g(N-1)2a 02 /( 3)0 40 M0.0
Figure 6: (color online). The mode frequency of breathing
-0.5
(m = 0) and dipole (m = ±1) modes as a function of in-
teraction strength at a fixed SO coupling λ = 1 and at (B) n = 1 (D) n = 3
SO -1.0
g =1.1g. 0 1 2 3 4 0 1 2 3 4 5
↑↓
/a /a
Figure 8: (color online). Bogoliubov wave-functions of the
lowest four breathing modes at λ = 1, g(N − 1) =
SO
40(cid:126)ω /a2 and g = 1.1g. The mode frequencies are indi-
⊥ ⊥ ↑↓
4 (a) m=-1 (b) m=0 (c) m=1 cated in Fig. 7b by solid symbols.
3 D
E. Dynamical Calculations
/ C
2
B To investigate the dynamical properties of the sys-
1 tem, we also perform direct simulations of the system
by real-time propagation of the ground state under per-
0 A turbation. To do this, firstly we obtain the ground state
0.0 0.2 0.4 0.6 0.8 0.0 0.2 0.4 0.6 0.8 0.0 0.2 0.4 0.6 0.8 1.0
S O SO SO by solving the coupled GP equations in Eqn. (3) using
the TSSP technique. The half-quantum vortex ground
Figure 7: (color online). The mode frequency of breathing
state is perturbed in various ways. We observe that the
(m = 0) and dipole (m = ±1) modes as a function of SO
couplingatafixedinteractionstrengthg(N−1)=40(cid:126)ω /a2 modefrequenciesobtainedbydynamicalsimulationagree
⊥ ⊥ wellwiththoseobtainedbysolvingBogoliubovequations
and at g =1.1g.
↑↓
(shown in Fig. 6).
Breathing mode analysis, m = 0: We excite the
monopole mode by weak relaxation of the trapping fre-
quencyattimet=0,andlettingthesystempropagatein
interaction strength. With increasing interaction, the real-time. As the breathing mode excitation is isotropic
mode frequency decreases and seems to saturate at suffi- in x-y space, it is sufficient to observe the dynamic re-
ciently large interactions. This may be anticipated from sponse of the collective coordinate along one axis, say,
the point of view of two-fluid hydrodynamic behavior in the x-axis. Here, we pick the mean square of the center-
the Thomas-Fermi regime. In Fig. 7, we report the de- of-mass coordinate as the quantity of interest:
pendence of the mode frequencies on SO coupling. In ´
the absence of SO coupling, the breathing mode with (cid:104)x2(cid:105) = ´|φσ|2x2dxdy ,
ω = 2ω⊥ and the dipole mode with ω = ω⊥ are the σ |φσ|2dxdy
exact solutions of quantum many-body systems in har-
monic traps. At a finite SO coupling, however, we find where σ = ↑,↓-spin components. In Fig. 9 (a),(b), we
that these two solutions are no longer exact. The rela- plot the time response of (cid:104)x2(t)(cid:105)σ for a typical param-
tive deviations of the breathing mode and dipole mode eter set. In Fig. 9 (c),(d), we show the correspond-
at λ =1 are about 10% and 30%, respectively. ing frequency response by plotting the single-sided am-
SO
plitude spectrum |(cid:104)x2(ω)(cid:105)| , which are just the Fourier
σ
In Fig. 8, we plot the Bogoliubov wave-functions of transforms of (cid:104)x2(t)(cid:105) . We observe frequency peaks at
σ
the lowest four breathing modes at λ =1, g(N−1)= ω/ω (cid:39)0.46,1.8,2.18 and at 3.40 (not shown). We note
SO ⊥
40(cid:126)ω /a2 and g = 1.1g. We find that the density re- that these values exactly match with the mode frequen-
⊥ ⊥ ↑↓
sponseismainlycarriedbyu (ρ)andu (ρ)components. ciesobtainedforthisparametersetbysolvingBogoliubov
↑ ↓
With increasing mode frequency, more and more nodes equations, shown in Fig. 6(b).
appear in u (ρ) and u (ρ). In contrast, the response in Dynamical calculations also reveal the coupling be-
↑ ↓
v (ρ) and v (ρ) is relatively weak and the curve shape is tweenthecenter-of-massmotionandtheinternalspinde-
↑ ↓
nearly unchanged as the mode frequency increases. greesoffreedom,atrademarksignatureofspin-orbitcou-
8
0.03(a) <x2(t)>↑ ( c )|<x2(ω)>|↑ center-of-mass coordinate in x-direction:
.01
´
0
|φ |2xdxdy
(cid:104)x(cid:105) = ´ σ .
σ |φ |2dxdy
-0.03 σ
0
(b) <x2(t)>↓ ( d )|<x2(ω)>|
↓
0.05
In Fig. 11 (a),(b), we plot the time response of
.02
0 this collective coordinate in x-direction of ↑- and ↓-
spin components for a typical parameter set. In
-0.05
Fig. 11 (c),(d), we show the corresponding frequency
response by plotting the single-sided amplitude spec-
0 2ω t 4 0 ω1 /ω 2
┴ trum |(cid:104)x(ω)(cid:105)| . We observe frequency peaks at ω/ω (cid:39)
┴ σ ⊥
0.05,0.43,0.70,1.25,1.34, (shown) and at 2.5, 2.64, 2.76
Figure 9: (color online). (a),(b): Dynamic response of the (not shown). We note that these values exactly agree
meansquareofthecenter-of-masscoordinateinx-directionof the mode frequencies obtained for this parameter set by
↑- and ↓- spin components respectively. We have shifted the solving Bogoliubov equations, shown in Fig. 6(a),(c).
curves by subtracting the time-averaged (cid:104)x2(t)(cid:105) . Without
σ
thisshift,theFourierspectrumasshownin(c)and(d)willbe In the inset of Fig. 11(a), we show the dynamics of
dominated by a large peak at ω =0. (c),(d): Corresponding the center-of-mass coordinate. It is important to note
single-sided amplitude spectrum of the collective coordinate. that even though the trap is displaced only in the x-
Parameters used: λSO =1.0,g(N −1)=40(cid:126)ω⊥/a2⊥,g↑↓/g= direction, we also observe a similar dynamic response in
1.1. y-directionofbothspincomponents(only↑-spincompo-
nent shown). This behavior occurs due to the vorticity
induced by the spin-orbit coupling — the vortex state
pledsystems. Weshallnowdiscuss´thedynamicresponse experiences a Magnus force that is perpendicular to its
ofthepopulationdifference∆n= dr(|φ↑|2−|φ↓|2). In motion. Hence a displacement in the x-direction induces
Fig. 10(a), we plot the time response of ∆n(t) for the amotionalongthey-direction. Furthermore,thetraceof
same parameter set mentioned in Fig. 9. In Fig. 10(b), the center-of-mass and its magnitude are affected by the
we show the corresponding frequency response by plot- strength of the inter-particle interactions and the spin-
ting the single-sided amplitude spectrum |∆n(ω)|. We orbit coupling induced population transfer, as observed
observe frequency peaks at ω/ω⊥ (cid:39)0.46,1.8,2.18 and at inthecaseofthebreathingmodeexcitation,betweenthe
3.40 (not shown), exactly matching with the mode fre- ↑- and ↓- spin components.
quencies obtained in Fig. 9. This analysis clearly shows
that the population transfer between the two spin com-
ponents shares a similar dynamic response with the col-
lective motional coordinate. In this aspect, response of (a) <x(t)>↑ ( c )|<x(ω)>↑|
0
∆n in a spin-orbit coupled spinor BEC (shown here) is .02
similar to the effects observed in the presence of internal 0.05
Josephsoncouplinginmulti-componentcondensates[27]. -0.05 y<(t)>↑ 0
0.15 -0.05-0.<1x(t)>0↑ 0
.04
0.08 (a)Δn(t) (b)Δn(ω) (d)|<x(ω)>↓|
0
0 .02 .02
-0.15 (b) <x(t)>↓
-0.08 0
0 5 10 0.5 1 1.5
0 2ω t 4 0 1ω/ω 2 ω┴t ω/ω┴
┴ ┴
Figure11: (coloronline). Parametersused: λ =1.0,g(N−
SO
Figure 10: (color online). (a) Dynamic response and (b) 1) = 40(cid:126)ω⊥/a2⊥,g↑↓ = 1.1g. (a),(b): Dynamic response of
single-sided amplitude spectrum of population difference ∆n the center-of-mass coordinate in x-direction of ↑- and ↓- spin
for the same parameter set used in Fig. 9. componentsrespectively. Theinsetin(a)showsthedynamics
ofthecenter-of-masscoordinateover12trapperiodsandthe
filled (red) marker denotes the initial position. (c),(d): Cor-
Dipole mode analysis, m = ±1: We excite the dipole responding single-sided amplitude spectrum of the collective
modes by displacing the trap in x-direction by a small coordinate.
amount at time t = 0, and letting the system propagate
in real-time. We observe the dynamic response of the
9
IV. INSTABILITY ANALYSIS AND PHASE interactions. It always appears to be positive, though
DIAGRAM theinteractionstendtodecreaseitsabsolutemagnitude.
Hence, there must be a quantum phase transition occur-
We are now ready to analyze the parameter space for ring at the isotropic point g = g↑↓. Once g > g↑↓, a
the existence of half-quantum vortex state. It will be- superposition state with density pattern,
come unstable with respect to increasing the interaction
(cid:34) (cid:35)
strength or decreasing the ratio g↑↓/g. The instability n = 1 φ2↑+φ2↓ ±φ φ cos(ϕ−ϕ ) , (40)
could be indicated from some energy considerations and ↑,↓ 2π 2 ↑ ↓ 0
from the softening of collective density modes.
becomes preferable. The 2D contour plot of this density
pattern with ϕ = 0 is schematically shown in the inset
0
A. Superposition instability of Fig. 1 (in the phase IIA).
In general, in passing the quantum phase transition
As we mentioned earlier, for any half-quantum vor- point, we would observe softening of a particular mode
√
tex state, φ(r) = [φ (ρ),φ (ρ)eiϕ]T/ 2π, we always frequency. As the superposition state involves a time-
↑ ↓
have a degenerate time-reversal partner state, Tφ(r) = reversal state with angular momentum m = −1, the
√
[φ (ρ)e−iϕ,−φ (ρ)]T/ 2π. There is an instability for lowest dipole mode with m = −1 may become unsta-
↓ ↑
half-quantum vortex state with respect to a superposi- ble. In Fig. 12(c), we plot the lowest dipole mode fre-
tion state, which with equal weight takes the form, quency ωm=−1 as a function of g↑↓/g at λSO = 1 and
g(N −1) = 20(cid:126)ω /a2. Indeed, with decreasing g /g,
φ (r)= √1 (cid:20)φ↑(ρ)+φ↓(ρ)e−i(ϕ−ϕ0) (cid:21). (38) the mode frequenc⊥y ω⊥m=−1 decreases and approach↑e↓s to
s 4π φ↓(ρ)ei(ϕ−ϕ0)−φ↑(ρ) zero exactly at the phase transition point.
Here ϕ is an arbitrary azimuthal angle. The energy
0
difference between the superposition state and the half- B. Instability to high-order angular momentum
quantum vortex state is given by, components
(g −g)(N −1)
∆E = ↑↓ W[φ(r)]. (39) There is another instability for the half-quantum vor-
GP 4 tex state, occurring with increasing the interatomic in-
teractions. With sufficiently large interactions, we an-
Therefore, if W[φ(r)]>0, the half-quantum vortex state
ticipate that the state with high-order azimuthal angu-
is stable only when g <g .
↑↓ lar momentum will energetically become favorable. For
example, let us consider a condensate state with an az-

0.2 0.04 imuthal angular momentum m = 1 (the 3/2-quantum
2
(a) g(N-1)= /a (b) SO=1 vortex state), which has the form,
W W (cid:20) (cid:21)
0.1 0.02 φ (r)= √1 φ↑(ρ)eiϕ . (41)
m=1 2π φ↓(ρ)ei2ϕ
0.0 0.00 The GP energy of this state can be obtained by solv-
0 1 2SO 3 4 0 1g0(N-1)2a 02/( 30) 40 ing the GPE equation as before, except that we need
0.08 to take R (ρ) and R (ρ) as the expansion functions
2 k1 k2
0.06 (c) SO=1, g(N-1)=20( /a ) for φ↑(ρ) and φ↓(ρ), respectively. Its degenerate time-
/1 reversal partner state has an azimuthal angular momen-
m=-0.04 superposition instability tum m=−2.
0.02 It is easy to see from Fig. 13(a) that beyond a critical
interaction strength the condensate state with m = 1,
0.00
0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 φ (r), is lower in energy than the half-quantum vor-
m=1
g /g tex state, φ (r). We note, however, that the critical
m=0
interaction strength determined in this way is not accu-
Figure12: (coloronline). (a)TheW-functionasafunctionof
SO coupling at g(N −1)=(cid:126)ω⊥/a2⊥ and g↑↓ =1.1g. (b) The rate, as a superposition state of φm=0(r) and φm=1(r)
W-function as a function of interaction strength at λ = may already become energetically more preferable than
SO
1 and g = 1.1g. (c) The instability of the lowest dipole φ (r) at a smaller interaction strength.
↑↓ m=1
mode frequency ωm=−1 with decreasing g↑↓/g at λSO = 1 An accurate determination of the threshold could be
and g(N −1)=20(cid:126)ω⊥/a2⊥. obtained by monitoring the instability in a particular
collective mode. As the condensate state may preserve
InFigs.12(a)and(b), wechecktheW-functionofthe a well-defined parity, we find that the instability oc-
half-quantum vortex state in the presence of interatomic curs in the lowest quadrupole mode with m = −2. In
10
starts to involve high-order angular momentum compo-
2.0
nents. The critical interaction strength g can be simply
2 c
2 / (a) determined from the softening of the mode frequency,
)
O
S ω (g =g )=0.
( 1.5 m=−2 c
+ In Fig. 14, we present critical interacting strength as
)
a function of SO coupling at g =g and g =2g. The
↑↓ ↑↓
(N 1.0 m=0 solid line at the isotropic point g↑↓/g has been recently
E/ calculated by Xiang-Fa Zhou and Congjun Wu by using
m=1 (or -2)
animaginarytimeevolutionmethod[13,28]. Ourresults
0.5 are in excellent agreement with theirs. We find that at
0 2 4 6 8 10 12 14 16 18 20
2 smaller SO coupling the critical interaction strength de-
g(N-1)a /( ) creases rapidly with increasing g /g.
0.4 ↑↓
(b) Re( )
0.3
Im( )
C. Instability against anisotropy in SO coupling
/m=-20.2 quadrupole instability strength
0.1
So far we have focused our attention on the half-
0.0 quantum vortex state supported by an isotropic 2D har-
monic trap subject to an isotropic Rashba SO coupling.
-0.1
0 1 2 3 4 5 6 7 8 9 10 Here we discuss the effect of the anisotropy in SO cou-
2
g(N-1)a /( ) plingstrengthλR onthestabilityofhalf-quantumvortex
state. The effect of the trap anisotropy will be discussed
Figure 13: (color online). (a) GP energy of the 3/2-quantum
in the next subsection. In the context of ultracold gases,
vortex state φ (r) and of the half-quantum vortex state
m=1 anisotropicRashbaspin-orbitcoupledwasfirstdiscussed
φ (r) as a function of interaction strength at λ = 2
m=0 SO
in Ref. [6] and the coupled GP equations were solved for
and g /g = 1.1. Beyond a critical interaction strength as
↑↓
indicated by an arrow, φ (r) becomes energetically favor- a many-body system in the absence of the trap and in
m=1
able. (b) The corresponding lowest quadrupole mode fre- the restricted scenario when g↑↓ = g. Here, we move
quency ω . It becomes unstable beyond a threshold g . beyond these restrictions and discuss the ground state
m=−2 c
of the system. We write the SO coupling term in the
form VSO = −i(λyσˆx∂y −λxσˆy∂x), where λx,λy are SO
20 coupling strengths in the two perpendicular directions.
By including this SO coupling term and solving the cou-
g =g
pled GP equations under the Hamiltonian as given in
) 15 g =2g Eq. (3) using the TSSP technique, we obtain the ground
state wavefunction at various values of anisotropy in SO
(
2 / coupling represented by λx/λy. In Fig. 15, we plot the
)a 10 IIB correspondinggroundstatedensityprofilesof↓-spincom-
N-1 ponent for an SO coupling strength of λx =4.0, and for
g( various values of λx/λy.
5 We see from Fig. 15(a) that the half-quantum vortex
I
state is indeed the ground state (already mentioned in
Fig.1(b))fortheparameterset: g(N−1)=0.1(cid:126)ω /a2,
⊥ ⊥
0 g /g = 1.1, λ = 4.0 and λ /λ = 1.0. We shall now
0 1 2 3 4 ↑↓ x x y
analyze the pattern in which the density profile changes
SO with anisotropy in SO coupling strength as shown in
Figure 14: (color online). Phase diagram at g = g and Fig. 15(b)-(d). It is evident from the density distribu-
↑↓
g↑↓ =2g. Thecriticalinteractionstrengthhasbeenshownas tions in Fig. 15, that the half-quantum vortex state is
a function of SO coupling. unstable even against small anisotropy in SO coupling
strength. Adopting a similar method as presented in
Ref. [29], we analyze this systematically by expanding
Fig. 13(b), we report the lowest quadrupole mode fre- the wavefunction of ↓-component in an orthogonal ba-
quency ω as a function of the interaction strength. sis set of the form: Φ (ρ) = Σ f (ρ)ei(2n+1)ϕ, where
m=−2 ↓ n n
As the interaction increases, the real part of mode fre- n measures the vorticity, and f (ρ) absorbs the nth
n
quency decreases down to zero and then, the imaginary mode’s contribution in radial direction. We quantify the
part becomes positive, indicating clearly that this mode weightsof´thewavefunctioninthenthmodebycomput-
willexponentiallygrowifthecondensateisinitiallyinthe ing a = dρ|f (ρ)|2. In Fig. 16, we plot the weights
n n
half-quantumvortexconfiguration. Thecondensatethen a relativetoa computedforhalf-quantumvortexstate
n 0