Vector solitons in a spin-orbit coupled spin-2 Bose-Einstein condensate Sandeep Gautam∗1 and S. K. Adhikari†1 1Instituto de Física Teórica, Universidade Estadual Paulista - UNESP, 01.140-070 São Paulo, São Paulo, Brazil (Dated: July 28, 2018) Five-component minimum-energy bound states and mobile vector solitons of a spin-orbit-coupled quasi-one-dimensional hyperfine-spin-2 Bose-Einstein condensate are studied using the numerical solution and variational approximation of a mean-field model. Two distinct types of solutions with single-peak and multi-peak density distribution of the components are identified in different domains of interaction parameters. From an analysis of Galilean invariance and time-reversal symmetry of the Hamiltonian, we establish that vector solitons with multi-peak density distribution preserve time- reversal symmetry, but cannot propagate maintaining the shape of individual components. However, those with single-peak density distribution violate time-reversal symmetry of the Hamiltonian, but can propagate with a constant velocity maintaining the shape of individual components. PACS numbers: 03.75.Mn, 03.75.Hh, 67.85.Bc, 67.85.Fg I. INTRODUCTION Bright solitons are localized wave packets that can tra- verse at a constant speed without changing their shape due to a balancing of nonlinear attraction and dispersive effects. Solitons have been studied in a variety of sys- tems which include water waves, non-linear optics, Bose- Einstein condensates (BECs), etc [1]. In case of BECs, the nonlinear attraction arises because of the net attrac- tive inter-atomic interactions, whereas the kinetic energy term in the Hamiltonian is the source of dispersion. In binary mixtures of scalar BECs, the presence of two com- ponents allows for the emergence bright vector solitons [2]. Solitons have also been studied in spinor BECs with hyperfine spin f = 1 [3] and f = 2 [4] in the absence of spin-orbit (SO) coupling. In a neutral atom, there is no coupling between the spin and the center-of-mass motion of the atom [5]. The SO coupling can be generated by controlling the atom-light interaction, and hence creating artificial Abelian and non-Abelian gauge potentials cou- pled to the atoms [6]. The first experimental realization of SO coupling [7], with equal Rashba [8] and Dressel- haus [9] strengths, in a neutral atomic BEC by dressing two atomic spin states with a pair of lasers paved the way for other experimental studies on SO-coupled spinor BECs [10]. Recently, the effect of SO coupling on the solitonic structures has been studied in pseudospin-1/2 [11–13] and spin-1 BECs [14] In this paper, we study an SO-coupled spin-2 BEC in a quasi-one-dimensional (quasi-1D) geometry character- ized by three interaction strengths c0 ∝ (4a2 + 3a4)/7, c1 ∝ (a4 − a2)/7, and c2 = (7a0 − 10a2 + 3a4)/7, where a0, a2, and a4 are s-wave scattering lengths in total spin ftot = 0, 2, and 4 channels, respectively [16, 17]. Our the- oretical analysis based on the SO-coupled single-particle ∗sandeepgautam24@gmail.com †adhikari44@yahoo.com, URL http://www.ift.unesp.br/users/adhikari Hamiltonian allows us to identify two subdomains with distinct physical properties separated by the c2 = 20c1 line. In the c2 < 20c1 subdomain, the component den- sities of the wave function show a multi-peak structure, whereas in the c2 > 20c1 subdomain, a single-peak struc- ture of the same is found. Based on the solutions of the SO-coupled single-particle Hamiltonian, we identify an appropriate variational ansatz to calculate the minimum- energy bright soliton solutions for an SO-coupled BEC. The single-particle SO-coupled Hamiltonian is found to break Galilean invariance. By solving the Galilean trans- formed single-particle Hamiltonian dynamics, we demon- strate that only single-peak solitons can propagate with- out changing their shape. In order to find the moving vector solitons of the SO- coupled spinor BEC, we transform the SO-coupled GP equation using a Galilean transformation. Two degener- ate solutions of the GP equation in the rest frame evolve into two velocity-dependent nondegenerate solutions in the moving frame. We show that the multi-peak soliton of Ref. [14] is a linear combination of two such degen- erate solutions in the rest frame. In the moving frame the degeneracy is removed, and the multi-peak soliton cannot be formed by the same linear combination. The multi-peak soliton of Ref. [14] is thus velocity depen- dent and cannot propagate maintaining its shape. No such linear combination of two degenerate solutions in the rest frame is needed in the formation of the single- peak soliton identified in this paper. These solitons can emerge as the ground states of the coupled GP equation in a moving frame, and hence can move with a constant velocity maintaining their shape, including the densities and phase profiles of individual components, unlike the moving multi-peak solitons of Refs. [12, 14] which show spin-mixing dynamics. In Sec. II, we describe the mean-field formalism used to study the SO-coupled spin-2 BECs. Here, taking into account the effect of interactions on the solutions of the non-interacting and trapless system, we identify two types of ground state solutions with single-peak and ar X iv :1 50 5. 06 19 0v 2 [ co nd -m at .q ua nt -g as ] 2 5 M ay 2 01 5 http://www.ift.unesp.br/users/adhikari 2 multi-peak densities, respectively. In case of ground-state solutions with single-peak densities, the reduction of the coupled GP equation to a single nonlinear partial differ- ential equation using single-mode approximation (SMA) is also discussed in this section. In SMA all the compo- nent wavefunctions are assumed to have the same spatial dependence but can have different normalization [15]. In Sec. III, we study the minimum-energy stationary soli- tons by variationally minimizing the energy. We also investigate the moving solitons from a consideration of Galilean invariance of the Hamiltonian. In Sec. IV, we present the numerical results for both types of solitons. We conclude by providing a summary of the study in Sec. V. II. SPIN-ORBIT COUPLED BEC IN A QUASI-1D TRAP Although, unlike in the case of charged electrons in an atom, an SO coupling interaction does not naturally arise in the case of neutral atoms in a spinor BEC, an SO coupling can be realized by an engineering of atom-light interaction. An SO-coupled spinor BEC in a quasi-1D trap can be realized by making the trapping frequencies along the two axes, say y and z (ωy, ωz), much larger than that along the x axis (ωx). This leads to a BEC in which the dynamics is restricted only along x axis. The single-particle Hamiltonian of this BEC with in a quasi-1D trap is [7, 18] H0 = p2x 2m + V (x) + γpxΣx, (1) where px = −i~∂/∂x is the momentum operator along x axis, V (x) = mω2 xx 2/2 is the harmonic trapping potential along x axis, and Σx is the irreducible representation of x component of the spin operator: Σx =  0 1 0 0 0 1 0 √ 3/2 0 0 0 √ 3/2 0 √ 3/2 0 0 0 √ 3/2 0 1 0 0 0 1 0  . (2) The SO coupling corresponds to equal strengths of Rashba (Σxpx + Σypy) and Dresselhaus (Σxpx − Σypy) couplings [18–20], where Σy is the irreducible representa- tion of y component of the spin operator. An equivalence of these couplings with the usual Rashba (Σxpy −Σypx) [8] and Dresselhaus (Σxpy + Σypx) [9] couplings can be shown by a rotation in the spin space [18, 19]. In our previous studies [21, 22], we employed a dif- ferent SO coupling (pxΣz) again with equal Rashba [8] and Dresselhaus [9] strengths. The present single-particle Hamiltonian (1) already couples the different spin compo- nents. The previous SO interaction [21, 22] did not cou- ple the different spin components at the single-particle level and the coupling was achieved after introducing the nonlinear spinor interactions in the BEC. Taking interactions in the Hartree approximation and using the single-particle Hamiltonian (1), a quasi-1D spin-2 BEC can be described by the following set of five coupled mean-field partial differential equations for the wave-function components ψj [20] i~ ∂ψ±2 ∂t = ( − ~2 2m ∂2 ∂x2 + V (x) + c0ρ ) ψ±2 − i~γ ∂ψ±1 ∂x + c1 ( F∓ψ±1 ± 2Fzψ±2 ) + ( c2/ √ 5 ) Θψ∗∓2, (3) i~ ∂ψ±1 ∂t = ( − ~2 2m ∂2 ∂x2 + V (x) + c0ρ ) ψ±1 − i~γ ∂ψ±2 ∂x − i~γ √ 6 2 ∂ψ0 ∂x + c1 (√ 3/2F∓ψ0 + F±ψ±2 ± Fzψ±1 ) − ( c2/ √ 5 ) Θψ∗∓1, (4) i~ ∂ψ0 ∂t = ( − ~2 2m ∂2 ∂x2 + V (x) + c0ρ ) ψ0 − i~γ √ 6 2 ( ∂ψ1 ∂x + ∂ψ−1 ∂x ) + √ 6 2 c1 ( F−ψ−1 + F+ψ1 ) + c2√ 5 Θψ∗0 , (5) where c0 = 2~2(4a2 + 3a4)/(7ml2yz), c1 = 2~2(a4 − a2)/(7ml2yz), c2 = 2~2(7a0 − 10a2 + 3a4)/(7ml2yz), a0, a2 and a4 are the s-wave scattering lengths in the to- tal spin ftot = 0, 2 and 4 channels, respectively, ρj(x) = |ψj |2 with j = 2, 1, 0,−1,−2 are the component densi- ties, ρ(x) = ∑2 j=−2 ρj is the total density, and lyz =√ ~/(mωyz) with ωyz = √ ωyωz is the oscillator length in the transverse y − z plane and F+ = F ∗− = Fx + iFy = 2(ψ∗2ψ1 + ψ∗−1ψ−2) + √ 6(ψ∗1ψ0 + ψ∗0ψ−1), Fz = 2(|ψ2|2 − |ψ−2|2) + |ψ1|2 − |ψ−1|2, Θ = 2ψ2ψ−2 − 2ψ1ψ−1 + ψ2 0√ 5 , where Fx, Fy, Fz are the three components of the spin- density vector F, and Θ is the spin-singlet pair amplitude [20]. Before proceeding further, let us transform the Eqs. (3)-(5) into dimensionless form using t̃ = ωxt, x̃ = x l0 , φj(x̃, t̃) = √ l0√ N ψj(x̃, t̃), (6) where l0 = √ ~/(mωx) is the oscillator length along x axis, and N is the total number of atoms. Then, Eqs. 3 (3)-(5) in dimensionless form become i ∂φ±2 ∂t̃ = ( −1 2 ∂2 ∂x̃2 + Ṽ (x̃) + c̃0ρ̃ ) φ±2 − iγ̃ ∂φ±1 ∂x̃ + c̃1 ( F̃∓φ±1 ± 2F̃z̃φ±2 ) + ( c̃2/ √ 5 ) Θ̃φ∗∓2, (7) i ∂φ±1 ∂t̃ = ( −1 2 ∂2 ∂x̃2 + Ṽ (x̃) + c̃0ρ̃ ) φ±1 − iγ̃ ∂φ±2 ∂x̃ − iγ̃ √ 6 2 ∂φ0 ∂x̃ + c̃1 (√ 3/2F̃∓φ0 + F̃±φ±2 ± F̃z̃φ±1 ) − ( c̃2/ √ 5 ) Θ̃φ∗∓1, (8) i ∂φ0 ∂t̃ = ( −1 2 ∂2 ∂x̃2 + Ṽ (x̃) + c̃0ρ̃ ) φ0 − iγ̃ √ 6 2 ( ∂φ1 ∂x̃ + ∂φ−1 ∂x̃ ) + √ 6 2 c̃1 ( F̃−φ−1 + F̃+φ1 ) + c̃2√ 5 Θ̃φ∗0, , (9) where Ṽ = x̃2/2, γ̃ = ~kr/(mωxl0), c̃0 = 2N(4a2 + 3a4)l0/(7l 2 yz), c̃1 = 2N(a4 − a2)l0/(7l 2 yz), c̃2 = 2N(7a0 − 10a2 + 3a4)l0/(7l 2 yz), ρ̃j(x̃) = |φj |2 with j = 2, 1, 0,−1,−2, and ρ̃(x̃) = ∑2 j=−2 |φj |2 and F̃+ = F̃ ∗− =2(φ∗2φ1 + φ∗−1φ−2) + √ 6(φ∗1φ0 + φ∗0φ−1), (10) F̃z̃ =2(|φ2|2 − |φ−2|2) + |φ1|2 − |φ−1|2, (11) Θ̃ = 2φ2φ−2 − 2φ1φ−1 + φ20√ 5 . (12) The total density is now normalized to unity, i.e.,∫∞ −∞ ρ̃(x̃)dx̃ = 1. We will represent the dimensionless variables without tildes in the rest of the manuscript for notational simplicity. For a non-interacting system in the absence of a trap- ping potential, there are five linearly independent solu- tions of Eqs. (7)-(9): Φ1 = eikx 4 ( 1, 2, √ 6, 2, 1 )T , (13) Φ2 = eikx 4 ( 1,−2, √ 6,−2, 1 )T , (14) Φ3 = eikx 2 (−1,−1, 0, 1, 1) T , (15) Φ4 = eikx 2 (−1, 1, 0,−1, 1) T , (16) Φ5 = eikx √ 3 8 ( 1, 0,− √ 2 3 , 0, 1 )T , (17) where T stands for transpose. Here Φ1 and Φ2 are two degenerate solutions corresponding to energy per particle E/N = k2/2 ± 2kγ, which has a minimum Emin/N = −2γ2 at k = ∓2γ, where the upper and lower signs are for Φ1 and Φ2, respectively. Similarly, the eigen functions -2-4 0 2 4 0 -1 1 2 3 -3 -2 k�Γ E �HN Γ 2 L 21 3 4 5 FIG. 1: (Color online) Energy dispersion E/(Nγ2) vs. k/γ for the eigen functions of the single-particle SO-coupled Hamil- tonian for the five wave functions Φi, where i = 1, ..., 5. Φ3 and Φ4 are also degenerate with energy E/N = k2/2± kγ, which has a minimum Emin/N = −γ2/2 at k = ∓γ with upper and lower signs for Φ3 and Φ4, respectively; whereas the non-degenerate eigen function Φ5 has energy E/N = k2/2, which has a minimum Emin/N = 0 at k = 0. The E/(Nγ2) vs. k/γ dispersion curves are shown in Fig. 1 and we will consider the lowest-energy ground state Φ1,2 with k = ∓2γ. Compared to the SO-coupled pseudospin-1/2 [23, 24] and spin-1 condensates [23], there are two sets of degen- erate eigen functions in the spin-2 case. Due to this, two distinct types of muti-peak stationary profiles can emerge in SO coupled spin-2 condensate. The most general solution of Eqs. (7)-(9) with mini- mum energy for a BEC with a uniform density n in the absence of interactions and trapping potential can be ob- tained by the linear superposition of √ nΦ1 and √ nΦ2 √ nΦ =  φ2 φ1 φ0 φ−1 φ−2  = √ n (α1Φ1 + α2Φ2) , (18) = √ n 4  α1e −i2γx + α2e i2γx 2α1e −i2γx − 2α2e i2γx √ 6α1e −i2γx + √ 6α2e i2γx 2α1e −i2γx − 2α2e i2γx α1e −i2γx + α2e i2γx  , (19) where |α1|2 + |α2|2 = 1, so that Φ is normalized to unity. The solution (19) implies that the magnetization M≡ ∫ Fzdz = 0. It implies that we can only obtain sta- ble bright solitons with zero magnetization for the SO- coupled spinor condensate considered in this paper. On the other hand, in the absence of SO coupling, one can have stable bright solitons with non-zero magnetization for spin-2 condensate [4]. The time-reversal-symmetric single-particle Hamilto- nian (1) violates parity. This will have interesting conse- quences on its eigen functions. The effect of time-reversal 4 symmetry operator T acting on a spin or orbital angular- momentum state |j,m〉 is [20] (T φ)m = (−1)mφ∗m. (20) Time-reversal symmetry of a quantum state Φ requires that it should be the same as its time-reversed state apart from a phase factor. The degenerate states (13) and (14) violate time-reversal symmetry and are con- nected by time-reversal symmetry operator: T |Φ1〉 = −|Φ2〉, T |Φ2〉 = −|Φ1〉. In general, for arbitrary α1 and α2 the states (19) do not have time-reversal symmetry. In the presence of interactions, the interaction energy per particle is [20] εint = c0 2 n+ c1 2n |F|2 + c2 2n |Θ|2, = c0 2 n+ 2n [ c1 + |α1|2|α2|2 ( c2 − 20c1 5 )] .(21) If c2 < 20c1, the minimum of εint corresponds to |α1| = |α2| = 1/ √ 2. This state is nondegenerate and has multi- peak density distribution for the wave-function compo- nents and is time-reversal invariant. The spin expec- tation and absolute value of spin-singlet pair amplitude per particle for this state are, respectively, |F|/n = 0 and |Θ|/n = 1/ √ 5. On the other hand, for c2 > 20c1, the εint can be minimized if |α1| = 1, |α2| = 0 or |α1| = 0, |α2| = 1. These two plane-wave states are de- generate and are connected by the time-reversal operator and hence break time-reversal symmetry of the Hamil- tonian. The spin expectation per particle for this state |F|/n = 2, and the absolute value of spin-singlet pair am- plitude per particle |Θ|/n = 0. These states have single- peak density distribution for the wave-function compo- nents and can also be studied using SMA. In SMA, the order parameter can be written as Φ(x, t) = 1 4 (1, 2, √ 6, 2, 1)TφSMA(x, t), (22) based on the minimum-energy solutions (13) of the single-particle Hamiltonian. For SMA to be valid, all the spin-dependent interactions should be much smaller than the spin-independent interactions [20]. Hence, us- ing Eq. (22) in Eqs. (7)-(9) and neglecting the c1- and c2-dependent terms, one can obtain the single nonlinear differential equation i ∂φSMA ∂t = [ −1 2 ∂2 ∂x2 + V (x) + c0ρ− 2iγ ∂ ∂x ] φSMA,(23) where the tildes have been dropped. Equation (23) will be solved numerically. III. BRIGHT SOLITONS A. Variational analysis We just found that the physical properties of the sys- tem are different for c2 > 20c1 and c2 < 20c1, and now we present a variational analysis of the bright soliton with minimum energy in these two domains. We variationally minimize the energy of the trapless BEC [20, 25], given by E = N ∫ ∞ −∞ { 1 2 2∑ j=−2 ∣∣∣∣dφjdx ∣∣∣∣2 − iγφ∗2 dφ1dx − iγφ∗−2 dφ−1dx − iγ √ 6 2 φ∗0 ( dφ1 dx + dφ−1 dx ) − iγφ∗1 dφ2 dx − iγφ∗−1 dφ−2 dx − iγ √ 6 2 dφ0 dx ( φ∗1 + φ∗−1 ) + c0ρ 2 + c1|F|2 + c2|Θ|2 2 } dx, (24) in the two domains which are separated by the c2 = 20c1 line. An appropriate variational ansatz Φvar is con- structed by taking the product of the superposition of the eigen functions corresponding to minimum energy in Eqs. (13) and (14) with a localized spatial soliton, i.e., Φvar = √ σ 2 Φsech(σx), (25) where σ is the variational parameter characterizing the width and strength of the soliton, and |α1|2 + |α2|2 = 1. As discussed in Sec. II, the exact values of α1 and α2 differ in the two domains. In order to search for the stable stationary bright solitons with higher energies, one can also consider the following variational ansatz: Φvar = √ σ 2 (α1Φ3 + α2Φ4)sech(σx), (26) Φvar = √ σ 2 Φ5sech(σx). (27) Ansatz (26) and (27) can lead to bright solitons with suc- cessively higher energies than the one described by Eq. (25). Nevertheless, in the present paper, our emphasis is to calculate the minimum-energy bright solitons consis- tent with Eq. (25). If c2 < 20c1, we choose |α1| = |α2| = 1/ √ 2 in Eq. (25). Then, the energy of the soliton is E = N 30 ( −60γ2 + 5σ2 + 5σc0 + σc2 ) (28) with a minima at σ = − 1 10 (5c0 + c2) , (29) provided that 5c0 + c2 < 0. The energy and hence the shape of the soliton are not the functions of interaction parameter c1 due the fact F = 0 for |α1| = |α2| = 1 √ 2 as has been discussed in Sec. II. This parameter domain cor- responds to the antiferromagnetic phase (F = 0, |Θ|/n = 1/ √ 5) [17]. The state Φvar, defined by Eqs. (25) and (19) with |α1| = |α2| = 1/ √ 2 is time-reversal symmetric. 5 If c2 > 20c1, we choose |α1| = 1, |α2| = 0 or vice versa in Eq. (25). The energy of the soliton, then, is E = N 6 ( −12γ2 + σ2 + σc0 + 4σc1 ) (30) with a minima at σ = −1 2 (c0 + 4c1) , (31) provided that c0 + 4c1 < 0. The shape of the soliton is independent of interaction parameter c2 due vanishing Θ. This parameter domain corresponds to the ferromag- netic phase (|F|/n = 2,Θ = 0) [17]. The Φvar with |α1| = 1, |α2| = 0 or vice versa leads to single-peak density distribution for wave-function components and breaks time-reversal symmetry. B. Moving bright solitons The GP equation governing the dynamics and statics of a scalar BEC is Galilean invariant. The implication of this invariance can be understood by considering the scalar 1D GP equation in dimensionless form i ∂φ(x, t) ∂t = −1 2 ∂2φ(x, t) ∂x2 − c|φ(x, t)|2φ(x, t), (32) where the nonlinearity is considered attractive (c > 0) for the formation of a bright soliton. This equation has the analytic solution φ(x, t) = √ (σ/c)sech [ (x− vt) √ σ ] eivx+i(σ−v 2)t/2, (33) where v is the velocity and σ represents the width and the strength of the soliton. It implies that the stationary soliton φ(x, 0) = √ (σ/c)sech [x √ σ] moves as the bright soliton φM (x, t) defined by [12] φM (x, t) = φ(x− vt, t)eivx+i(σ−v 2)t/2, (34) maintaining the width and strength (σ) fixed. In the case of a multi-component spinor BEC, Eq. (34) remains valid with a multi-component wave function Φ replacing the single-component φ while the multiplying exponential factor remains unchanged. Now, let us examine the Galilean invariance of the SO- coupled Hamiltonian. Equation (34) implies that for the Galilean transformation x′ = x + vt, t′ = t, where v is the relative velocity of unprimed coordinate system with respect to primed coordinate system, the wave function Φ of Eqs. (7)-(9) should transform to ΦM related by Φ(x, t) = ΦM (x′, t′)e−ivx ′−i(σ−v2)t′/2. (35) Substituting Eq. (35) in Eqs. (7)-(9) and using ∂/∂x = ∂/∂x′, ∂/∂t = ∂/∂t′ + v∂/∂x′, we get i ∂ΦM (x′, t′) ∂t′ = [ −1 2 ∂2 ∂x′2 − γΣx ( i ∂ ∂x′ + v )] ΦM (x′, t′), (36) where, for the sake of simplicity, we have suppressed the terms proportional to c0, c1, and c2, which should remain unchanged. We have also dropped a σ dependent addi- tive term in the Hamiltonian which does not contribute to the dynamics. The extra term −γΣxvΦM on the right hand side of the equation indicates that the SO-coupled Hamiltonian is no longer Galilean invariant. The SO-coupled equation (36) has the solutions Φ1 and Φ2 of Eqs. (13) and (14) with energies E = −2N(vγ + γ2) and E = 2N(vγ − γ2), respectively. Hence, the degeneracy between Φ1 and Φ2 is removed for v 6= 0. This in turn implies that the superposition of Φ1 and Φ2 in Eq. (18) in the rest frame, which leads to a multi-peak soliton, is not possible in the moving frame. In other words, the multi-peak soliton cannot propagate with a constant velocity maintaining its shape. On the other hand, for a single-peak soliton no mixing between Φ1 and Φ2 is allowed based on energetic considerations as discussed in Sec. II and Sec. III. Hence a single-peak soliton can traverse with a constant velocity maintaining its shape. C. Numerical solutions We numerically solve the coupled Eqs. (7)-(9) using split-time-step Crank-Nicolson method [26, 27] in imagi- nary and real times. The ground state is determined by solving Eqs. (7)-(9) in imaginary time. In order to fix both the magnetization and normalization, we use the approach discussed in Ref. [21]. Accordingly, after each iteration in imaginary time, we renormalize the compo- nent wave functions as φj(x, τ + δτ) = djφj(x, τ), (37) where dj ’s satisfy the following relations [21] d1d−1 = d20, (38) d2d−2 = d20, (39) d2d 2 −1 = d30, (40) and d81N2+d20d 6 1N1 + d40d 4 1N0 + d60d 2 1N−1 + d80N−2 = 1, (41) 2d81N2 + d20d 6 1N1 − d60d21N−1 − 2d80N−2 =M. (42) Hence, one needs to solve the coupled set of non-linear Eqs. (41)-(42) to determine d0 and d1, which can be back substituted in Eqs. (38)-(40) to calculate the remaining normalization constants. The spatial and time steps em- ployed to generate the numerical results in this paper are δx = 0.05 and δt = 0.0000625, respectively. IV. NUMERICAL RESULTS First, we consider an SO-coupled f = 2 spinor BEC of 10000 23Na atoms trapped in a harmonic trapping po- tential with ωx/(2π) = 20 Hz, ωy/(2π) = ωz/(2π) = 400 6 -8 -4 40 8 0 0.05 0.1 x Ρ jHx L HaL c0=242.97 c1=12.06 c2=-13.03 Γ = 1Ρ0 Ρ±1 Ρ±2 ´´´´´´´´´´´´´´´´´´´´´´´´´´´´´´´´´´´´´´´´´´´´´´´´øøøøøøøøøøøø ø ø ø ø ø ø ø øøøøøøøøøøø ø ø ø ø ø ø ø øøøøøøøøøøøòòòòòòòòòòò ò ò ò ò ò ò ò ò ò òòòòòòòòò ò ò ò ò ò ò ò ò ò òòòòòòòòòò -4-8 40 8 0 0.03 0.06 x Ρ jHx L c0=201.36 c1=-1.8 c2=24.15 Γ=1 HbL Ρ±2´ ò Ρ±1 Ρ0 ø SMA òǿ -8 -4 0 84 0 0.03 0.06 x Ρ jHx L HbL num. Ρ0 Ρ±1 Ρ±2 FIG. 2: (Color online) Ground state density profiles of an SO-coupled trapped spinor BEC for (a) c0 = 242.97, c1 = 12.06, c2 = −13.03 < 20c1, and (b) c0 = 201.36, c1 = −1.8, c2 = 24.15 > 20c1. Hz. The oscillator lengths with these parameters are l0 = 4.69 µm and lyz = 1.05 µm. This value of l0 has been used in all calculations to write the dimensionless GP equations (7)-(9). The scattering lengths of 23Na in total spin ftot = 0, 2 and 4 channels are a0 = 34.9aB , a2 = 45.8aB , a4 = 64.5aB [17, 20], respectively, resulting in c0 = 242.97, c1 = 12.06, c2 = −13.03 < 20c1. Here aB is the Bohr radius. The ground state density profile with M = 0 is shown in Fig. 2(a). The multi-peak nature of the solution, obtained as result of the superposition of two counter-propagating plane waves, is consistent with analytic results obtained in Sec. II. The solution is dy- namically stable and retains its density and phase pro- file if real-time propagation of the dynamics is performed upon small perturbation. Next we consider the trapped BEC in the c2 > 20c1 do- main. For this, we consider a0 = 52.35aB , a2 = 45.8, and a4 = 43aB , which results in c0 = 201.36, c1 = −1.8, c2 = 24.15 > 20c1. The necessary modification of the scatter- ing lengths can be achieved by the optical and magnetic Feshbach resonance techniques [28]. The trapping po- tential parameters and number of atoms are the same as those in Fig. 2(a). The ground state densities with M = 0 are shown in Fig. 2(b). The component den- sities obtained by using SMA, i.e., by solving Eq. (23), are in good agreement with the full numerical solution of ´´´´´´´´´´´´´´´´´´ ´ ´ ´ ´´´´´´´´´´´´´´´´´´´øøøøøøøøøøøøøøøøøø ø ø ø øø ø ø ø ø ø ø ø øø ø ø ø øøøøøøøøøøøøøøøøøøøòòòòòòòòòòòòò ò ò ò òò ò ò òò ò ò ò ò ò ò ò ò ò òò ò ò òò ò ò ò òòòòòòòòòòòòòò -2-4 40 4 0 0.1 0.2 0.3 x Ρ jHx L c2=0.48 c1=0.39 c0=-1.55 Γ=1 HaL Ρ±2´ ò Ρ±1 Ρ0 ø var. òǿ -2-4 20 4 0 0.1 0.2 0.3 x Ρ jHx L HaL num. Ρ Ρ0 Ρ±1 Ρ±2 ´ ´ ´ ´ ´ ´ ´ ´ ´ ´ ´ ´ ´ ´ ´ ´ ´ ´ ´ ´ ´ ´ ´ ´ ´ ´ ´ ´ ´ ´øøøøøøøøøøøø ø ø ø ø ø øøø ø ø ø ø ø øøøøøøøøøøøøøò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò òòò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò -2-4 40 4 0 0.1 0.2 0.3 x Ρ jHx L c2=-8.06 c1=-1.8 c0=4.94 Γ=1 HbL Ρ±2´ ò Ρ±1 Ρ0 ø var. òǿ -2-4 20 4 0 0.1 0.2 0.3 x Ρ jHx L HbL num. Ρ0 Ρ±1 Ρ±2 ´ ´ ´ ´ ´ ´ ´ ´ ´ ´ ´ ´ ´ ´ ´ ´ ´ ´ ´ ´ ´ ´ ´ ´ ´ ´ ´ ´ ´ ´øøøøøøøøøø ø ø ø ø ø ø ø ø ø øøøøøøøøøøøò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò ò -2 -1 20 1 0 0.1 0.2 0.3 0.4 x Ρ jHx L c0=-2.74 Γ=1 HcL Ρ±2´ ò Ρ±1 Ρ0 ø var. òǿ � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � �õõõõõõõõõõõõõõ õ õ õ õ õõõ õ õ õ õ õ õõõõõõõõõõõõõõã ã ã ã ã ã ã ã ã ã ã ã ã ã ã ã ã ã ã ã ã ã ã ã ã ã ã ã ã ã ã ã ã ã ã ã ã ã ã ã -2 -1 20 1 0 0.1 0.2 0.3 0.4 x Ρ jHx L c1=-0.16 c2=1.77 Γ=1 HcL Ρ±2� ã Ρ±1 Ρ0 õ SMA ãõ� -2 -1 20 1 0 0.1 0.2 0.3 0.4 x Ρ jHx L HcL num. Ρ0 Ρ±1 Ρ±2 FIG. 3: (Color online) Numerical (num.) and variational (var.) densities of an SO-coupled f = 2 BEC soliton for (a) c0 = −1.55, c1 = 0.39, c2 = 0.48, (b) c0 = 4.94, c1 = −1.80, c2 = −8.06, and (c) c0 = −2.74, c1 = −0.16, and c2 = 1.77 > 20c1. In (c) the densities obtained by SMA are also shown. Eqs. (7)-(9). The solutions in this case are dynamically stable plane wave solutions, resulting in a single-peak density distribution for the wave-function components, consistent with the analytical analysis in Sec. II. In order to obtain the bright solitons (V (x) = 0) in SO-coupled spinor BECs, we consider the two cases high- lighted in Sec. III A: (a) c2 < 20c1 and 5c0 + c2 < 0, and (b) c2 > 20c1 and c0 + 4c1 < 0. In case (a), we consider a0 = −1.5aB , a2 = −1.2aB , a4 = 0, and N = 5000, which results in c0 = −1.55, c1 = 0.39, c2 = 0.48. The numerical results for the component densities of the soli- ton are illustrated in Fig. 3(a) together with the corre- sponding variational results, defined by Eqs. (25) and (29), withM = 0 and |α1| = |α2| = 1/ √ 2. The solution in this case has multi-peak density distribution, is time- 7 -3-5 30 5 0 0.1 0.2 x Ρ jHx L HaL c0=-1.55 c1=0.39 c2=0.48 Γ = 1 num. Ρ±1 Ρ±2 -2-4 20 4 0 0.1 0.2 0.3 x Ρ jHx L HbL c0=-1.55 c1=0.39 c2=0.48 Γ = 0 num. Ρ-1 Ρ2 FIG. 4: (Color online) (a) Nonzero numerical densities of an SO-coupled f = 2 BEC soliton with energy larger than the bright soliton shown in Fig. 3(a) . (b) Nonzero numerical densities for an f = 2 BEC soliton in cyclic phase in the absence of SO coupling. The interaction parameters for both (a) and (b) are c0 = −1.55, c1 = 0.39, c2 = 0.48. reversal symmetric, and dynamically stable. Despite the modulation in densities of the individual components, the total density profile (ρ) of the bright soliton as shown in Fig. 3(a) does not have any modulation unlike the bright solitons discussed in Ref. [11]. In this sense the present multi-peak solitons are spin-2 analogues of the stripe phase discussed in Ref. [23, 24]. Moreover, in Ref. [11] the authors noted component symmetries like Real(Ψ↑) = − Real(Ψ↓), and Imag(Ψ↑) = Imag(Ψ↓), for the real and imaginary parts, which is not possible in the present model. This is because of different SO coupling used in the two studies, e.g., γpxΣx in this paper and γpxσz in Ref. [11]. Had we used the SO coupling γpxΣz the component symmetries of Ref. [11] can be obtained as noted in Table 1 of Ref. [21]. In order to obtain the stationary bright soliton with higher energy, we use imaginary time propagation with (Φ3 + Φ4) exp(−x2/2)/( √ 2 √ π) as the initial guess for the order parameter. The dynamically stable bright soli- ton thus obtained with the parameters of Fig. 3(a) is shown in Fig. 4(a). The period of density modu- lation in this case is twice of the period in Fig. 3(a). Hence, there are two distinct stable multi-peak solitons in SO-coupled spin-2 condensate as compared to SO- coupled spin-1 condensate which can have only one type of multi-peak soliton [14]. There also exists a stable single-peak soliton, not shown here, with order parameter (a) ρ±2(x) ρ±1(x) ρ0(x) -25 -20 -15 -10 -5 0x 0 4 8 12 16 t 0 0.1 0.2 0.3 ρ (x ) (b) ρ±2(x) ρ±1(x) ρ0(x) -30 -20 -10 0 10 20 30 x 0 5 10 15 20 25 30 35 t 0 0.1 0.2 0.3 0.4 0.5 ρ (x ) FIG. 5: (Color online) (a) Dynamics of a soliton of Fig. 3(a) placed at x = −20 and set into motion with a con- stant velocity by multiplying the wave function by a phase eix. Here c0 = −1.55, c1 = 0.39, c2 = 0.48, M = 0, and v = 1. (b) Elastic collision between the two moving solitons of Fig. 3(b) placed at x = ±20 and set into mo- tion by multiplying by phase factors e∓ix, respectively. Here c0 = 4.94, c1 = −1.80, c2 = −8.06, N = 5000, M = 0, and v = 1. √ 3ρ(x)/8(1, 0,− √ 2/3, 0, 1) and energy higher than the aforementioned two multi-peak solitons. For the same set of parameters, these three types of solitons − minimum- energy ground and higher-energy excited states − have the same total density ρ(x). In the absence of SO coupling the minimum-energy soliton corresponds to time-reversal symmetry breaking cyclic phase (|F| = |Θ| = 0) for c1 > 0 and c2 > 0, c2 < 20c1 [17, 20, 21]. This is shown in Fig. 4(b) for the interaction parameters corresponding to Fig. 3(a) and γ = 0. The presence of SO coupling results in the time-reversal symmetric antiferromagnetic phase for which |F| = 0 and |Θ| > 0. Next we consider an SO-coupled soliton with a0 = 3.4aB , a2 = 4.58aB , a4 = −1.0aB , and N = 5000, which is equivalent to c0 = 4.94, c1 = −1.80, c2 = −8.06, c2 > 20c1. The nu- merical results in this case are contrasted in Fig. 3(b) with the variational results, defined by Eqs. (25) and (31) with |α1| = 1, |α2| = 0 or vice versa. The soli- tonic solution here breaks the time-reversal symmetry of the Hamiltonian, as can be seen, for example, from Eqs. (25) and (13) with α2 = 0. In the absence of SO coupling the density distribution remains unchanged in this case for the bright soliton with zero magnetization, but the 8 component wavefunctions have constant phase instead of constant phase gradient. The solution thus retains the time-reversal symmetry of the Hamiltonian. The SMA can not be applied in this case to calculate the bright solitonic solution as c0 > 0. To study the applicabil- ity of SMA, we consider c0 = −2.74, c1 = −0.16, and c2 = 1.77 > 20c1. The bright soliton obtained by the nu- merical solution of Eqs. (7)- (9), their variational approx- imation, viz. Eqs. (25) and (31) with |α1| = 1, |α2| = 0, and SMA, viz. Eq. (23, are shown in Fig. 3(c). In order to set the soliton into motion with a constant velocity v, we multiply the static solution by exp (ivx) with v = 1. We find that bright solitons with standing wave density profile in the c2 < 20c1 and 5c0 + c2 < 0 parameter domain, exhibit spin-mixing dynamics as is shown in Fig. 5(a) for the soliton with same interaction parameters as in Fig. 3(a) and moving with velocity v = 1 in dimensionless units. Thus, this type of soliton does not strictly preserve its shape while in motion. On the other hand, moving solitons obtained by multiplying the plane wave solutions in the c0 + 4c1 < 0 and c2 > 20c1 domain, by exp(±ivx) strictly preserve their shape. This is evident from Fig. 5(b), where we numerically study the collision between two bright solitons, initially located at x = ±20, and moving in opposite directions with a speed v = 1 (in dimensionless units). The solitons collide at x = 0 and pass through each other without suffering a change in their shapes as is evident from Fig. 5(b). The interaction parameters in this case are the same as in Fig. 3(b). As discussed in the Sec. III B, these moving five- component vector solitons do not exhibit any spin-mixing dynamics. V. SUMMARY We study the generation and propagation of five- component vector solitons in an SO-coupled spinor BEC with hyperfine spin f = 2 using a mean-field GP equa- tion with three interaction strengths: c0, c1, and c2. Two types of solitons − single-peak and multi-peak − emerge in this case for c2 > 20c1 and c2 < 20c1, respectively. In the former case, the solutions violate time-reversal sym- metry, whereas in the latter case, the solutions are time- reversal symmetric. The GP equation for this system is demonstrated to violate Galelian invariance. Analyzing the Galelian invariance of this equation, we show that the single-peak solitons can move with a constant veloc- ity maintaining constant component densities. On the other hand, the multi-peak SO-coupled solitons change the component densities during motion. Acknowledgments This study is financed by the Fundação de Amparo à Pesquisa do Estado de São Paulo (Brazil) under con- tracts 2013/07213-0 and 2012/00451-0 and also by the Conselho Nacional de Desenvolvimento Científico e Tec- nológico (Brazil) under project 303280/2014-0. [1] Y. S. Kivshar and B. A. Malomed, Rev. Mod. Phys. 61, 763 (1989); F. K. Abdullaev, A. Gammal, A. M. Kam- chatnov, and L. Tomio, Int. J. Mod. Phys. B 19, 3415 (2005). [2] V. M. Pèrez-Garc̀ia and J. B. Beitia, Phys. Rev. A 72, 033620 (2005); S. K. Adhikari, Phys. Lett. A 346, 179 (2005); Phys. Rev. A 72, 053608 (2005); L. Salasnich and B. A. Malomed, Phys. Rev. A 74, 053610 (2006). [3] J. Ieda, T. Miyakawa, and M. Wadati, J. Phys. Soc. Jpn. 73, 2996 (2004); J. Ieda, T. Miyakawa, and M. Wadati, Phys. Rev. Lett. 93, 194102 (2004); L. Li, Z. Li, B. A. Malomed, D. Mihalache, and W. M. Liu, Phys. Rev. A 72, 033611 (2005); M. Uchiyama, J. Ieda, and M.Wadati, J. Phys. Soc. Jpn. 75, 064002 (2006); W. Zhang, Ö. E. Müstecaplioǧlu, and L. You, Phys. Rev. A 75, 043601 (2007); B. J. Dąbrowska-Wüster, E. A. Ostrovskaya, T. J. Alexander, and Y. S. Kivshar, Phys. Rev. A 75, 023617 (2007); E. V. Doktorov, J. Wang, and J. Yang, Phys. Rev. A 77, 043617 (2008); B. Xiong and J. Gong; Phys. Rev. A 81, 033618 (2010); P. Szankowski, M. Trippenbach, E. Infeld, and G. Rowlands, Phys. Rev. Lett. 105, 125302 (2010). [4] M. Uchiyama, J. Ieda, and M. Wadati, J. Phys. Soc. Jpn. 76, 074005 (2007). [5] Y. Li, Giovanni I. Martone, and S. Stringari, arXiv:1410.5526; V. Galitski and I. B. Spielman, Nature (London) 494, 49 (2013). [6] K. Osterloh, M. Baig, L. Santos, P. Zoller, and M. Lewen- stein, Phys. Rev. Lett. 95, 010403 (2005); J. Ruseckas, G. Juzeliūnas, P. Öhberg, and M. Fleischhauer, Phys. Rev. Lett. 95, 010404 (2005); J. Dalibard, F. Gerbier, G. Juzeliūnas, and P. Öhberg, G. Juzeliūnas, J. Ruseckas, and J. Dalibard, Phys. Rev. A 81, 053403 (2010); Z. Lan and P. Öhberg, Phys. Rev. A 89, 023630 (2014); Rev. Mod. Phys. 83, 1523 (2011). [7] Y.-J. Lin , K. Jiménez-García, and I. B. Spielman, Nature (London) 471, 83 (2011). [8] Y. A. Bychkov E. I. Rashba, J. Phys. C 17, 6039 (1984). [9] G. Dresselhaus, Phys. Rev. 100, 580 (1955). [10] M. Aidelsburger, M. Atala, S. Nascimbène, S. Trotzky, Y.-A. Chen, and I. Bloch, Phys. Rev. Lett. 107, 255301 (2011); Z. Fu, P. Wang, and S. Chai, L. Huang, and J. Zhang, Phys. Rev. A 84, 043609 (2011); J.-Y. Zhang, S.- C. Ji, Z. Chen, L. Zhang, Z.-D. Du, B. Yan, G.-S. Pan, B. Zhao, Y.-J. Deng, H. Zhai, S. Chen, and J.-W. Pan, Phys. Rev. Lett. 109, 115301 (2012); C. Qu, C. Hamner, M. Gong, C. Zhang, and P. Engels, Phys. Rev. A 88, 021604(R) (2013). [11] V. Achilleos, D. J. Frantzeskakis, P. G. Kevrekidis, and D. E. Pelinovsky, Phys. Rev. Lett. 110, 264101 (2013). [12] Y. Xu, Y. Zhang, and B. Wu, Phys. Rev. A 87, 013614 (2013). http://arxiv.org/abs/1410.5526 9 [13] L. Salasnich, W. B. Cardoso, and B. A. Malomed, Phys. Rev. A 90, 033629 (2014); V. Achilleos, D. J. Frantzeskakis, and P. G. Kevrekidis, Phys. Rev. A 89, 033636 (2014); S. Cao, C.-J. Shan, D.-W. Zhang, X. Qin, and J. Xu, J. Opt. Soc. Am. B 32, 201 (2015). [14] Y.-K. Liu and S.-J. Yang, Eur. Phys. Lett., 108, 30004 (2014). [15] C. K. Law, H. Pu, and N. P. Bigelow, Phys. Rev. Lett. 81, 5257 (1998); H. Pu, C. K. Law, S. Raghavan, J. H. Eberly, and N. P. Bigelow Phys. Rev. A 60, 1463 (1999); S. Yi, Ö. E. Müstecaplioḡlu, C. P. Sun, and L. You, Phys. Rev. A 66, 011601(R) (2002). [16] M. Koashi and M. Ueda, Phys. Rev. Lett. 84, 1066 (2000). [17] C. V. Ciobanu, S.-K. Yip, and T.-L. Ho, Phys. Rev. A 61, 033607 (2000). [18] C. Wang, C. Gao, C.-M. Jian, and H. Zhai, Phys. Rev. Lett. 105, 160403 (2010). [19] H. Zhai, Int. J. of Mod. Phys. B 26, 1230001 (2012). [20] Y. Kawaguchi and M. Ueda, Phys. Rep. 520, 253 (2012). [21] S. Gautam and S. K. Adhikari, Phys. Rev. A 91, 013624 (2015). [22] S. Gautam and S. K. Adhikari, Phys. Rev. A 90, 043619 (2014). [23] C. Wang, C. Gao, C.-M. Jian, and H. Zhai, Phys. Rev. Lett. 105, 160403 (2010). [24] T.-L. Ho and S. Zhang, Phys. Rev. Lett. 107, 150403 (2011); S. Sinha, R. Nath, and L. Santos, Phys. Rev. Lett. 107, 270401 (2011); Y. Li, L. P. Pitaevskii, and S. Stringari, Phys. Rev. Lett. 108, 225301 (2012). [25] T. Kawakami, T. Mizushima, and K. Machida, Phys. Rev. A 84, 011607(R) (2011); Z. F. Xu, R. Lü, and L. You, Phys. Rev. A 83, 053602 (2011); Z. F. Xu, Y. Kawaguchi, L. You, and M. Ueda, Phys. Rev. A 86, 033628 (2012). [26] H. Wang, J. Comput. Phys., 230, 6155 (2011); 274, 473 (2014). [27] P. Muruganandam and S. K. Adhikari, Comput. Phys. Commun. 180, 1888 (2009); D. Vudragovic, I. Vidanovic, A. Balaz, P. Muruganandam, and S. K. Adhikari, Com- put. Phys. Commun. 183, 2021 (2012). [28] S. Inouye, M. R. Andrews, J. Stenger, H.-J. Miesner, D. M. Stamper-Kurn, and W. Ketterle, Nature (London) 392, 151 (1998). I Introduction II Spin-orbit coupled BEC in a quasi-1D trap III Bright solitons A Variational analysis B Moving bright solitons C Numerical solutions IV Numerical Results V Summary Acknowledgments References