Bright-Dark Soliton Complexes in Spinor Bose-Einstein Condensates

# Bright-Dark Soliton Complexes in Spinor Bose-Einstein Condensates

H. E. Nistazakis Department of Physics, University of Athens, Panepistimiopolis, Zografos, Athens 15784, Greece    D.J. Frantzeskakis Department of Physics, University of Athens, Panepistimiopolis, Zografos, Athens 15784, Greece    P.G. Kevrekidis Department of Mathematics and Statistics, University of Massachusetts, Amherst MA 01003-4515, USA    B.A. Malomed Department of Interdisciplinary Studies, Faculty of Engineering, Tel Aviv University, Tel Aviv 69978, Israel    R. Carretero-González Nonlinear Dynamical Systems Group, Department of Mathematics and Statistics, and Computational Science Research Center, San Diego State University, San Diego CA, 92182-7720, USA
Submitted to Phys. Rev. A, May 2007
###### Abstract

We present bright-dark vector solitons in quasi-one-dimensional spinor () Bose-Einstein condensates. Using a multiscale expansion technique, we reduce the corresponding nonintegrable system of three coupled Gross-Pitaevskii equations (GPEs) to a completely integrable Yajima-Oikawa system. In this way, we obtain approximate solutions for small-amplitude vector solitons of dark-dark-bright and bright-bright-dark types, in terms of the spinor components, respectively. By means of numerical simulations of the full GPE system, we demonstrate that these states indeed feature soliton properties, i.e., they propagate undistorted and undergo quasi-elastic collisions. It is also shown that, in the presence of a parabolic trap of strength , the bright component(s) is (are) guided by the dark one(s), and, as a result, the small-amplitude vector soliton as a whole performs harmonic oscillations of frequency in the shallow soliton limit. We investigate numerically deviations from this prediction, as the depth of the solitons is increased, as well as when the strength of the spin-dependent interaction is modified.

## I Introduction

The development of far-off-resonant optical techniques for trapping of ultracold atomic gases has opened new directions in the studies of Bose-Einstein condensates (BECs), allowing one to confine atoms regardless of their spin (hyperfine) state, see, e.g., Ref. ket0 (). One of major achievements in this direction was the experimental creation of spinor BECs ket1 (); cahn (), in which the spin degree of freedom (frozen in magnetic traps) comes into play. This gave rise to the observation of various phenomena that are not present in single-component BECs, including formation of spin domains spindomain () and spin textures spintext ().

A spinor condensate formed by atoms with spin is described by a -component macroscopic wave function. Accordingly, a number of theoretical works have been dealing with multi-component (vector) solitons in spinor BECs. Bright wad1 (); boris (); zh () and dark wad2 () solitons, as well as gap solitons ofyspin (), have been predicted in this context (the latter type requires the presence of an optical lattice). However, mixed vector soliton solutions (in particular, ones composed of bright and dark components) for the respective system of coupled Gross-Pitaevskii equations (GPEs) have not been reported yet, to the best of our knowledge. Actually, compound solitons of the mixed type may be of particular interest, as they would provide for the possibility of all-matter-wave waveguiding, with the dark soliton component forming an effective guide for the bright component, similar to the all-optical waveguiding studied in detail (chiefly, theoretically) in nonlinear optics kiag (). Such waveguides could be useful for applications, such as quantum switches and splitters emulating their optical counterparts bld ().

On the other hand, mixed solitons of the dark-bright type were considered in a model of a two-component condensate, described by two coupled GPEs BA (). Actually, the model also assumed that the two components represented different spin states of the same atomic species, with equal scattering lengths of the intra-component and inter-component atomic collisions (i.e., the matrix of the nonlinear coefficients in the coupled GPEs was of the Manakov’s type, that makes the system integrable in the absence of the external potentials).

In this work we consider a quasi-one-dimensional (quasi-1D) spinor condensate, described by a system of three coupled GPEs. In the physically relevant case of Rb and Na atoms with , which are known to form spinor condensates of ferromagnetic and polar types, respectively (the definitions are given below), the system includes a naturally occurring small parameter, , namely, the ratio of the strengths of the spin-dependent and spin-independent interatomic interactions kemp (); greene (). In the case of and without the external potential, the system of the three coupled GPEs becomes the so-called Manakov system manakov (), which is completely integrable intman (). Exploiting the smallness of , we will develop a multiscale-expansion method to asymptotically reduce the nonintegrable GPE system to another completely integrable one, viz., the Yajima-Oikawa (YO) system. The latter was originally derived to describe the interaction of Langmuir and sound waves in plasmas yo () and has been used in studies of vector solitons in the context of optics ofyol () and binary BECs maximo (). The asymptotic reduction is valid for homogeneous polar spinor BECs (such as Na), that are not subject to the modulational instability ofy2 (); boris (). Borrowing exact soliton solutions of the YO system, we predict two types of vector-soliton complexes in the spinor condensate, viz., dark-dark-bright (DDB) and bright-bright-dark (BBD) ones for the spin components. Numerical simulations of the underlying (full) GPE system show that these solitary pulses (including ones with moderate, rather than small amplitudes) emulate genuine solitons (solitary waves in integrable systems) quite well, propagating undistorted for long times and undergoing quasi-elastic collisions (quasi-elastic collisions of solitons in the nonintegrable two-component system were mentioned earlier in Ref. BA ()).

The effect of the harmonic trapping potential on the solitons is also studied, analytically and numerically (the potential also breaks the exact integrability of the coupled GPE equations, even with ). It is shown that the small-amplitude vector solitons of the mixed types perform harmonic oscillations in the presence of the trap of strength , at frequency , which coincides with the well-known frequency of oscillations of a dark soliton in the single-component BEC oscfreq (). In particular, the bright component(s) of the vector soliton oscillates at the same frequency, following their dark counterpart(s) (ordinary bright solitons in single-component BECs oscillate at frequency , rather than st (), according to the Kohn’s theorem kohn ()). As a matter of fact, it is a manifestation of the guidance of the bright component by the dark one. The oscillations of moderate- and large-amplitude solitons are studied as well and we find that the oscillation frequency is down-shifted from its value as the soliton depth is increased. In the particular case of large-amplitude solitons (that perform small-amplitude oscillations around the trap center), the frequency down-shift can be very well approximated by the prediction of Ref. BA ().

We also investigate the effect of a larger normalized spin-dependent interaction strength on the stability of the predicted vector soliton solutions. In fact, we consider a value of an order of magnitude larger than the physically relevant one of the polar sodium spinor condensate ( rather than ) and investigate numerically the respective coupled GPEs. The result is that, generally, under such a strong perturbation the solitons emit stronger radiation and are eventually destroyed. However, even in such a case of a large , small- and moderate-amplitude DDB solitons are found to persist up to relatively large times, of order of more than 300 ms in physical units. Given that for the physically relevant small value of pertaining to the sodium spinor BEC, the respective lifetime was four times as large, we believe that the vector solitons derived in this work have a good chance to be observed experimentally.

The paper is organized as follows: In Sec. II we present the model, expound our analytical approach for the homogeneous system, and derive solutions for the bright-dark soliton complexes. Section III is devoted to the presentation of the numerical and analytical results for the dynamics of the solitons in both the homogeneous and inhomogeneous (harmonically confined) system. Finally, Sec. IV concludes the paper.

## Ii The model and its analytical consideration

### ii.1 The model

At sufficiently low temperatures (finite temperature effects have been considered recently in Refs. Guilleumas ()), and in the framework of the mean-field approach, the spinor BEC with is described by the vector order parameter, , with the components corresponding to the three values of the vertical spin projection, . Assuming that the condensate is confined in a highly anisotropic trap with frequencies , we may assume the wave functions approximately separable, , where the transverse wave function is the ground state of the respective harmonic oscillator. Then, averaging of the underlying system of the coupled three-dimensional (3D) GPEs in transverse plane gpe1d () leads to the following system of coupled 1D equations for the longitudinal components of the wave functions (see also Refs. wad1 (); wad2 (); boris (); zh (); ofyspin ()):

 iℏ∂tψ±1 = ^Hsiψ±1+c(1D)2(|ψ±1|2+|ψ0|2−|ψ∓1|2)ψ±1 (1) + c(1D)2ψ20ψ∗∓1, iℏ∂tψ0 = ^Hsiψ0+c(1D)2(|ψ−1|2+|ψ+1|2)ψ0 (2) + 2c(1D)2ψ−1ψ∗0ψ+1,

where the asterisk denotes the complex conjugate, and is the spin-independent part of the Hamiltonian, with being the total density ( is the atomic mass). The nonlinearity coefficients have an effectively 1D form, namely and , where is the transverse harmonic oscillator length, which defines the size of the transverse ground state. Finally, coupling constants and which account, respectively, for spin-independent and spin-dependent collisions between identical spin- bosons, are given by (in the mean-field approximation)

 c0=4πℏ2(a0+2a2)3m,c2=4πℏ2(a2−a0)3m, (3)

where and are the -wave scattering lengths in the symmetric channels with total spin of the colliding atoms and , respectively. Note that the spinor condensate may be either ferromagnetic (such as the Rb), characterized by , or polar (such as the Na), with ho (); ohmi ().

Measuring time, length and density in units of , and , respectively (where is the peak density), we cast Eqs. (1)–(2) in the dimensionless form,

 i∂tψ±1 = Hsiψ±1+δ(|ψ±1|2+|ψ0|2−|ψ∓1|2)ψ±1 (4) + δψ20ψ∗∓1, i∂tψ0 = Hsiψ0+δ(|ψ−1|2+|ψ+1|2)ψ0 (5) + 2δψ−1ψ∗0ψ+1,

where , the normalized harmonic trap strength is given by

 Ωtr=ωxω⊥32(a0+2a2)n0, (6)

and we define

 δ≡c(1D)2c(1D)0=a2−a0a0+2a2. (7)

According to what was said above, and correspond, respectively, to ferromagnetic and polar spinor BECs. In the relevant cases of Rb and Na atoms with , kemp () and greene (), respectively, i.e., in either case, plays the role of a small parameter of Eqs. (4)–(5).

Generally, Eqs. (4)–(5) give rise to spin-mixing states sm (). However, there also exist non-spin-mixing, or spin-polarized states, which are stable stationary solutions of Eqs. (4) and (5) ho (); yi (). Here, we first consider the spatially homogeneous system (), and focus on such solutions having at least one component equal to zero, the remaining ones being continuous waves (CWs). The corresponding exact stationary solutions are

 ψ−1=ψ+1=√μ2exp(−iμt),ψ0=0, (8)

and

 ψ−1=ψ+1=0,ψ0=√μexp(−iμt). (9)

As we demonstrate below, weakly nonlinear perturbations around these solutions take the form of three-component dark-bright soliton complexes. In particular, perturbations around solutions (8) and (9) will lead to solitons of the DDB and BBD types, respectively, for components and . Since the analytical approach and the derivation of the soliton solutions for the two cases are quite similar, we focus below on the DDB solitons, and discuss the BBD ones in a brief form.

It is relevant to note that, for solutions with , two equations (4) coalesce into a single one for wave function . Then, the transformation , casts the system of two equations (4) and (5) into the form of two coupled GPEs, with nonlinear cross-coupling coefficients , which was introduced in Ref. BA () (also with the objective to study bound complexes of dark-bright solitons, carried by fields and , respectively). In fact, the analysis in Ref. BA () was limited to the case of , while in this work we focus on effects generated by small . It should also be stressed that, although stationary equations presented below may indeed be found from the system of two, rather than three, coupled GPEs, the stability tests for the solutions are performed against general perturbations (see below), which include those with (i.e., the full system of the three equations was employed in the stability simulations).

### ii.2 Linear analysis

Aiming to find solutions of Eqs. (4)–(5) close to the CW solution given by Eq. (8), we start the analysis by adopting the following ansatz,

 ψ−1=ψ+1 = √n(x,t)exp[−iμt+iϕ(x,t)], ψ0 = Φ0(x,t)exp(−iμt), (10)

where and are real density and phase of fields , while is, generally, a complex function. Substituting Eq. (10) into Eqs. (4)–(5), we derive the following a system,

 −n[12(∂xϕ)2−12√n∂2x√n+δΦ20e−2iϕ]=0, (11)
 i∂tΦ0 + 12∂2xΦ0−(2n−μ+|Φ0|2)Φ0 (12) − 2δn(Φ0+Φ∗0e−2iϕ)=0.

The CW state (8) corresponds to an obvious solution of Eqs. (11) and (12) with , , . Next, we linearize the equations around this state, looking for a solution as , and , where is a formal small parameter. At order , the linearization leads to the following system:

 i(∂t~n+μ2∂2x~ϕ)−μ(∂t~ϕ+2~n−μ4∂2x~n)=0, (13) i∂t~Φ0+12∂2x~Φ0−δ⋅μ(~Φ0+~Φ∗0)=0. (14)

Combining real and imaginary parts of Eq. (13), we arrive at a dispersive wave equation,

 ∂2t~n−μ∂2x~n+(μ2/8)∂4x~n=0, (15)

which gives rise to a stable dispersion relation between wavenumber and frequency (the absence of complex roots for at real implies the modulational stability of the underlying CW state):

 ω2=μk2(1+μk2/8). (16)

It follows from Eq. (16) that, in the long-wave limit (), small-amplitude waves can propagate on top of CW solution (8) with the speed of sound,

 c=√μ. (17)

A similar analysis for Eq. (14), which is decoupled from Eq. (13), leads to dispersion relation

 ω2=k2(δμ+k2/4). (18)

It is clear from here that, for (which corresponds to the polar state), Eq. (18) has no complex roots for , hence the trivial solution to Eq. (12), , is modulationally stable. However, (corresponding to the ferromagnetic state) gives rise to modulational instability of the solution against the perturbations whose wavenumbers belong to the instability band, . Note that these results comply with those reported in Ref. ofy2 (). Below, we focus on the modulationally stable case, which pertains to the polar state with .

### ii.3 Asymptotic soliton solutions

To consider solutions for weakly nonlinear deviations from the CW state, we recall that is a small parameter of the system (4)–(5), which suggests to define the stretched variables,

 X≡δ1/2(x−√μt),T≡δt. (19)

Then, we seek for solutions of Eqs. (11)–(12) as

 n=(μ/2)+δ⋅ρ,ϕ=δ1/2α,Φ0=δ3/4q, (20)
 q≡q1cos(Kx−Ωt)+iq2sin(Kx−Ωt), (21)

where , , , while and are unknown wavenumber and frequency. Substituting Eqs. (20) in Eq. (11), at the leading order in , which is , we derive a relation between density and phase ,

 √μ∂Xα=2ρ. (22)

At the next order , the resulting equation is complex:

 −(iμ/4)(2∂Xρ−√μ∂2Xα)+∂Tα+|q|2=0. (23)

The imaginary part of the expression on the left-side of Eq. (23) vanishes due to the validity of Eq. (22), while the real part leads to equation . The compatibility condition of the latter with Eq. (22) leads to

 ∂Tρ=−(√μ/2)∂X(|q|2). (24)

We now proceed to Eq. (12), which, to leading order in , i.e., at , yields the following system:

 Ωq1−(K2/2)q2 = 0, −[(K2/2)+2μ]q1+Ωq2 = 0. (25)

Nontrivial solutions to Eqs. (25) are possible when the following dispersion relation for and holds:

 Ω2=K2(μ+K2/4). (26)

Next, to order , Eq. (12) leads to system

 −√μ∂Xq1+K∂Xq2 = 0, −K∂Xq1+√μ∂Xq1 = 0,

which has nontrivial solutions if . In combination with Eq. (26), the latter relation selects the frequency, . Finally, at order , Eq. (12) leads to equation

 i∂Tq+12∂2Xq−2ρq=0. (27)

Equations (24) and (27), which are the basic result of our analysis, constitute the Yajima-Oikawa (YO) system. It describes the interaction of low-frequency and high-frequency waves, and was originally derived in the context of plasma physics, where it applies to Langmuir (high-frequency) waves, forming a packet (soliton) moving at velocities close to the speed of sound, and thus strongly coupled to the ion-acoustic (low-frequency) waves yo (). As shown in Ref. yo (), the YO system is integrable by means of the inverse-scattering-transform, and gives rise to soliton solutions. The solitons have the shape for field , and shape for , which correspond to a density dip for components and a bright soliton for , as per Eqs. (20). According to Eq. (22), the phase profile of the components, in the form of , is associated to the density dip, hence the patterns in these components, generated by the exact solution of the YO system, are genuine dark solitons. The full form of the approximate (asymptotic) DDB soliton solution to Eqs. (4)–(5), into which the YO soliton is mapped by Eqs. (10), (19), and (20), is

 ×exp[−iμt−2iη√δ/μtanh(2√δηZ)], (28) ψ0(x,t)=23/2δ3/4ημ−1/4√ξsech(2√δηZ) ×exp[−iμt+i√μx−2i√δξZ+2iδ(η2−ξ2)t], (29)

where , while and are arbitrary parameters of order .

A similar analysis can be performed to derive asymptotic soliton solutions of the BBD type. In that case, starting from CW solution (9), we seek for solutions of Eqs. (4)–(5) in the form of

 ψ−1 = ψ+1=Φ0(x,t)exp(−iμt), ψ0 = √n(x,t)exp[−iμt+iϕ(x,t)]. (30)

Next, following the same analytical approach which has led above to the DDB soliton, we again end up with the YO system, in a form similar to Eqs. (24) and (27), namely,

 ∂Tρ=−2√μ∂X|q|2, i∂Tq+12∂2Xq−ρq=0. (31)

Eventually, the approximate BBD soliton solution to Eqs. (4)–(5), generated by the YO soliton, is

 ψ±1(x,t)=2δ3/4ημ−1/2√ξsech(2√δηZ) ×exp[−iμt+i√μx−2i√δξZ+2iδ(η2−ξ2)t], (32) ψ0(x,t)=√(μ/2)−4δη2sech2(2√δηZ) ×exp[−iμt−2iη√δ/μtanh(2√δηZ)]. (33)

As the latter solution is quite similar to the DDB one, given by Eqs. (28)–(29), below we only deal with the dynamics of the DDB solitons. It is worthwhile to note in passing that both types of solutions are genuinely traveling ones, i.e., they do not exist with zero speed.

## Iii Dynamics of DDB spinor solitons

### iii.1 Numerical Results

In order to test the prediction of the existence of the DDB solitons in the underlying spinor BEC model, we turn to numerical integration of the original GPEs (4)–(5). In particular, we fix (the value corresponding to Na) and use the following initial conditions for the densities,

 |ψ±1(x,t=0)|2 = 12[μ−ν sech2(√νx)], (34) |ψ0(x,t=0)|2 = ν3/2ξη√μsech2(√νx), (35)

while initial phase profiles are similar to those in Eqs. (28)–(29), the parameter that determines the initial width of the soliton being

 ν≡4η2δ. (36)

Other parameters are chosen as , , and (for the homogeneous condensate), or for the trapped condensate. In physical terms, this choice corresponds to the spinor condensate of sodium atoms with the peak 1D density m, which contains atoms confined in the trap with frequencies Hz; in this case, the time and space units are, respectively, ms and m.

Choosing value of the arbitrary parameter introduced above, and substituting , , we have checked that these values indeed provide for very good agreement of the analytical predictions with numerical results. However, we will display numerical results obtained for an essentially larger value of , viz., [which, as seen from Eq. (36) corresponds to and hence, from Eq. (34), to a soliton complex with deeper and narrower dark components, and similarly taller and narrower bright components]. In this way, we intend to showcase the really wide range of validity of the analytical approach, and the robustness of the obtained solitary wave solutions.

More specifically, we first check if the spinor DDB soliton complexes indeed behave as genuine solitons, in the small-amplitude limit. To this end, we take initial conditions in the form of the superposition of two different pulses,

 ψ±1(x,0)=√μ2−ν2[sech2(x+)+sech2(x−)] ×exp(−i√νμtanh(x+)+i√νμtanh(x−)), (37)
 ψ0(x,0) = ν3/4√ξη√μ[sech(x+)e+i√μνx+−i(ξ/η)x+ (38) +sech(x−)e−i√μνx−+i(ξ/η)x−],

where and are positions of centers of the two pulses. As seen in Eqs. (37)–(38), the soliton components are lent opposite initial momenta and, as a result, they propagate in opposite directions, as shown in the top panel of Fig. 1. We stress that even though a small amount of radiation is emitted (see the four bottom panels of Fig. 1), the two dark solitons in the fields, coupled to their bright counterparts in the component, propagate practically undistorted, and around they undergo a quasi-elastic collision; moreover, it is clearly observed that the solitons remain unscathed after the collision. This result is consistent with our asymptotic calculations performed above, indicating that the small-amplitude limit (for sufficiently small ) of the nonintegrable system of Eqs. (4)–(5) behaves like the integrable YO system.

Next, we consider the confined system, with . In this case, strictly speaking, the asymptotic reduction of the primary system of Eqs. (4)–(5) to the YO system [Eqs. (24) and (27)] is not valid. However, even in the presence of the potential term, the solutions obtained with may be used as an initial configuration set near the bottom of the trap, to generate DDB- (or BBD-) like solutions of the inhomogeneous system. To that end, we first integrate Eqs. (4)–(5) in imaginary time, finding a ground state of the Thomas-Fermi (TF) type for fields , which is approximated by the well-known analytical density profile pit (), . Then, at , the initial conditions for the components are taken as the numerically found TF profiles multiplied by the dark soliton as in Eq. (34), while the initial configuration of the field is taken as the bright soliton in Eq. (35).

In such a case, and given that the spinor DDB solitons were found above to be robust objects behaving similar to genuine solitons of an integrable system, one would expect that the solitons would perform harmonic oscillations in the presence of the (sufficiently weak) parabolic trap. Figure 2 shows that this is the case indeed: The DDB soliton, which was initially placed at the trap’s center, oscillates as a whole without significant deformations of its components up to large times [while the figure extends to (which is seconds in physical units), a similar behavior continues at still larger times]. This is a clear indication to the fact that the predicted DDB complexes have a good chance to be observed in the experiment. A noteworthy feature of the numerical data is that the bright-soliton component is guided by the dark ones, the entire soliton complex oscillating at a single frequency. The value of the frequency is estimated below.

### iii.2 The soliton’s oscillation frequency

The effect of the harmonic trap on the spinor-soliton dynamics can be studied analytically, using the asymptotic multiscale expansion, similar to how it was done in Refs. huang (); tonks1 (); tonks2 (). In those works, asymptotic reductions of various scalar GPE-based models, which included the trapping potential, led to Korteweg-de Vries equations with variable coefficients [instead of constant ones that would be the case for the homogeneous (untrapped) system]. In the present situation, we may expect, accordingly, that the inclusion of the harmonic potential may lead to a YO system with variable coefficients.

An important result that can be produced by such an analysis is the oscillation frequency of the solitons in the presence of the harmonic trap. As shown in Refs. huang (); tonks1 (); tonks2 () by means of the adiabatic perturbation theory for solitons RMP (), the oscillation frequency can be obtained from the equation (which is valid to leading order in the small parameter characterizing the inhomogeneity of the system),

 d~Xdt=~c(~X), (39)

where is a properly chosen slow spatial variable, and is the speed of sound, which is now position-dependent due to the background density of the trapped condensate. In fact, Eq. (39) is a straightforward generalization of the result obtained for the homogeneous system, where it has been found that the speed of sound is [see Eq. (17)], while the soliton’s velocity is , implying that (for sufficiently small ) .

To adopt this approach to the present case, we need to find the local speed of sound when the harmonic potential term, , is included in the spin-independent part of the Hamiltonian, . Taking into regard that the potential varies slowly on the soliton’s spatial scale, (see, e.g., Fig. 2), we define the above-mentioned slow spatial variable as , where [recall that given in Eq. (6) is of order ], and is an auxiliary scale parameter. This way, the trapping potential takes the form of , i.e., it depends only on slow variable .

Then, the local speed of sound can easily be derived upon considering the linearization of Eqs. (11) and (12), which are modified by the inclusion of the term in the left-hand side of Eq. (11). The ground state of this system can easily be found by setting and . Then, since Eq. (11) implies that is time-independent in the ground state, we assume that and, to the leading order in , we obtain

 n0(~X)=(1/2)[μ−V(~X)], (40)

in the region where , and outside. Equation (40), which is the TF approximation for the density profile, also implies that, for the axial size of the trapped condensate is . Similarly to the analysis presented above in Sec. II.2, we now consider the linearization around the ground state and seek for respective solutions to Eqs. (11)–(12) as , , and , with . This way, we obtain the following dispersion relation for the inhomogeneous system,

 ω2=2n0(~X)k2+(1/4)k4, (41)

and, accordingly, the local speed of sound:

 c(~X)=√2n0(~X), (42)

which bears resemblance to the sound propagation in weakly nonuniform media landau (); in the homogeneous case, Eq. (42) is reduced to Eq. (17).

Next, we substitute Eq. (42) in Eq. (39) and, taking into regard the density profile given by Eq. (40), we integrate the resulting first-order differential equation. The result is

 ~X=Lsin[(~Ωtr/√2)t], (43)

which demonstrates that a sufficiently shallow spinor dark soliton oscillates with frequency

 ωosc=Ωtr√2, (44)

similarly to the known result for the oscillations of dark solitons in a single-component BEC oscfreq () [note that, for the sake of simplicity, we dropped the tilde in Eq. (44)].

We should clearly stress here that we expect the above result to be valid in the case of shallow solitons, with a relative depth [see Eq. (34)]. For example, as seen in the left panel of Fig. 3, in the case of such a shallow soliton with (corresponding to the above mentioned physically relevant choice of and a value of the chemical potential ) the analytical prediction is quite accurate: the numerically found oscillation frequency (for ) is , while the analytical prediction of Eq. (44) is (the percentage error is ). On the other hand, as seen in the same figure, the amplitude of the soliton oscillation is , while (the respective percentage error is .

Although the result of Eq. (43) [and (44)] is “compatible” with our analytical approach developed in the previous section, which is also valid for small-amplitude solitons, it is also of interest to numerically investigate the oscillations of solitons of moderate and large amplitudes. In this respect, in the right panel of Fig. 3 we show the percentage error of the numerically found soliton oscillation frequency [relative to the prediction of Eq. (44)] as a function of the relative dark soliton depth . It is clear that the region corresponds to shallow solitons, while the limiting case of corresponds to a “black” soliton, with an initially zero intensity at the trap’s center (the latter is slightly displaced from the trap’s center to be set into motion). As seen in the figure, the analytical prediction of Eq. (44) is very good as concerns the estimation of the oscillation frequency, since the error is below for every . Notice that we have also computed the error in the oscillation amplitude (not shown here), which we have found to be larger (i.e., up to in this regime); this is due to the fact that the increasingly deeper solitons are not reflected at the rims of the condensate, but rather inside the cloud, as can be seen e.g., in Fig. 2.

In the case of moderate (and large) amplitude solitons the analytical prediction becomes worse: For example, in the particular case shown in Fig. 2 (in this case ), comparing the numerically found soliton’s oscillation frequency, to the aforementioned prediction of (again for ), we find a relatively large discrepancy () between the two values. However, an important observation regarding Fig. 2, is that the bright-soliton component performs oscillations at the same frequency, . This is a clear indication that the bright component is guided (effectively trapped) by the dark component of the DDB complex. Note that, in the single-component BEC, bright solitons oscillate in the parabolic potential with frequency, st () (as a consequence of the Kohn’s theorem kohn ()), rather than .

Furthermore, we observe that, naturally, the discrepancy becomes even larger in the case of large-amplitude (nearly-black) solitons, which perform small-amplitude oscillations around the trap’s center. For example, for () the numerically found values of the oscillation frequency deviate from the value of by (), while in the limiting case of the respective percentage error is . It is clear that such large deviations are due to the fact that the above numerical results pertain to non-small amplitude solitons with large values of , contrary to what is the case for the the analytical approach which requires , as mentioned above.

In such a case of large amplitude solitons, it is interesting to compare the numerically found oscillation frequencies to a different analytical prediction presented in Ref. BA (). In this work, the oscillation frequency may result from a nonlinear equation of motion for the bright-dark soliton in a binary BEC mixture (see Eq. (5) of BA ()). In fact, for shallow solitons the oscillation frequency is identical to the one given by Eq. (44), while, in the opposite limit of very deep dark solitons, may be approximated (in our language) as,

 ωosc=Ωtr√2(1−NB4√μ+(NB/4)2)1/2, (45)

where is the number of atoms of the bright-soliton component, which, employing Eq. (35), is easily found to be . It is clear that Eq. (45) shows that the oscillation frequency is down-shifted (as compared to the value of ) [i.e., the dark-bright pair executes slower oscillations, as the bright component is enhanced], thus having a quantitative agreement with our numerical observations. Perhaps more importantly, it also provides better quantitative estimates for the values of the soliton’s oscillation frequencies. In particular, for normalized soliton depths , and , the predictions of Eq. (45) deviate from the respective numerically found oscillation frequencies with percentage errors , and (recall that a respective comparison with the prediction of Eq. (44) led to percentage errors , , and , respectively). The above results indicate that a more detailed description of the motion of the bright-dark soliton complexes in the trapped spinor condensate, in the lines of Ref. BA (), would be very interesting. However, this is beyond the scope of the present work.

### iii.3 Effects of larger spin-dependent interactions

In the previous subsections we had kept the value of the parameter fixed, i.e., (corresponding to the polar spin-1 sodium BEC). Such a small value of validates our perturbative approach, which allowed us to find approximate DDB soliton solutions of the YO-type, and study their oscillations in the trapped spinor condensate. It is interesting, however, to test the stability and the dynamics of the obtained DDB solitons in the more general case of a stronger perturbation, considering also non-small values of . In this respect, we will now present numerical results obtained by direct numerical integration of Eqs. (4) and (5) for , which is an order of magnitude greater than the previous value. Note that we will study the evolution of DDB solitons with the same amplitudes as in the case of small , so as to directly compare the results pertaining to weak or moderate spin-dependent interaction strength.

In Fig. 4, we show the evolution of a moderate amplitude DDB soliton, characterized by the parameters , , and ; for the same values of the trap strength and chemical potential as before ( and ), the initial condition is the same as the one shown in the second-row left panel of Fig. 2 (in this case, recall that the normalized amplitude of the dark solitons in the components is ). As observed in Fig. 4, although the stronger perturbation induces emission of stronger radiation (see bottom left panel), the losses are not significant up to relatively large times, of order of (or ms in physical units): the dark (bright) soliton densities are only () smaller than their initial values, and it can safely be concluded that even for such a strong value of the DDB vector soliton has a good chance to be observed in an experiment (provided, of course, that such a magnitude of the spin-dependent interaction is experimentally achieved). However, for later times the continuous perturbation-induced emission of radiation results in the eventual destruction of the DDB complex. In particular, at , the solitons almost decay: at (see bottom right panel of Fig. 4), the density of the dark (bright) soliton is () smaller than its initial value.

The cases of large- and small-amplitude solitons, with normalized amplitudes (of the dark soliton) being and , respectively, was examined too (results not shown here). It was found that, naturally, the large-amplitude DDB soliton starts to accelerate immediately due to the strong radiation emitted and decays fast, being completely destroyed at : the soliton densities are smaller than of their initial values. On the other hand, the small-amplitude DDB soliton was found to be slightly more robust, showing a behavior similar to the one of the moderate amplitude soliton shown in Fig. 4, but for smaller times: in fact, the vector soliton persists up to (as per the above criterion, the densities are smaller than of their initial values) and then decays.

To conclude this subsection, it is clear that the small- and moderate-amplitude YO-type DDB vector soliton persist in the spin-1 condensate up to experimentally relevant times even for large spin-dependent interaction strength, namely an order of magnitude larger than the physically relevant value pertaining to the polar sodium spinor condensate.

## Iv Conclusions

We have studied bright-dark soliton complexes in polar spinor Bose-Einstein condensates, both analytically and numerically. Our analytical approach is based on the small amplitude-asymptotic reduction of the nonintegrable vector (three-component) system of the coupled Gross-Pitaevskii equations to a completely integrable model, viz., the Yajima-Oikawa (YO) system. Borrowing soliton solutions of the YO system and inverting the reduction, we have obtained an analytical approximation for small-amplitude vector solitons of the dark-dark-bright and bright-bright-dark types, in terms of the components, respectively. The analytical predictions were confirmed by direct numerical simulations. The so constructed approximate soliton states were found to propagate undistorted and undergo quasi-elastic collisions, featuring properties of genuine solitons.

The effect of the harmonic trapping potential (which also contributes toward the nonintegrability of the underlying equations) on the solitons was also studied numerically and analytically. It was found that even vector solitons of moderate (non-small) amplitudes maintain their identity in the presence of the parabolic trap, and perform harmonic oscillations, at large times ( seconds or even more, in physical units). We have found that the soliton’s oscillation frequency takes (in the analytical approximation) the value of , the same as featured by the dark soliton in the single-component repulsive trapped condensate. The numerical results verify the analytical prediction, for sufficiently shallow dark solitons: It was found that for initial soliton depths below of the chemical potential, the deviation of the analytical prediction from the numerically found oscillation frequency was below (the error in the estimation of the amplitude of the soliton oscillation was below ). For moderate and large amplitude solitons, the discrepancy in the frequency was larger (of the order of more than ); however, in the case of very deep dark solitons we found that the respective prediction of Ref. BA () led to a significantly smaller error, less than . Thus, an elaborated description of the bright-dark soliton motion (for solitons of arbitrary amplitude) in the trapped spinor condensate is certainly a challenge for future work.

We also tested the robustness of the derived vector soliton solutions in the case of a large normalized spin-dependent interaction strength, namely an order of magnitude larger than the physically relevant one corresponding to the polar sodium spinor condensate. We found that although the solitons are generically destroyed under such a strong perturbation, the lifetimes of small- and moderate-amplitude DDB solitons was more than 300 ms in physical units. Thus, the vector solitons obtained in this work have a good chance to be observed in an experiment.

The bright-soliton component(s) was found to be guided by their dark counterpart(s), oscillating with the frequency imposed by dark components. This is an example of the all-matter-wave soliton guidance, with potential applications in the design of quantum switches and splitters.

The work of H.E.N. and D.J.F. was partially supported by the Special Research Account of the University of Athens. H.E.N. acknowledges partial support from EC grants PYTHAGORAS-I. The work of B.A.M. was supported, in a part by the Israel Science Foundation through the Center-of-Excellence grant No. 8006/03, and German-Israel Foundation (GIF) No. 149/2006.

## References

• (1) D. M. Stamper-Kurn and W. Ketterle, cond-mat/0005001.
• (2) D. M. Stamper-Kurn, M. R. Andrews, A. P. Chikkatur, S. Inouye, H.-J. Miesner, J. Stenger, and W. Ketterle, Phys. Rev. Lett. 80, 2027 (1998).
• (3) M.-S. Chang, C. D. Hamley, M. D. Barrett, J. A. Sauer, K. M. Fortier, W. Zhang, L. You, and M. S. Chapman, Phys. Rev. Lett. 92, 140403 (2004).
• (4) J. Stenger, S. Inouye, D. M. Stamper-Kurn, H.-J. Miesner, A. P. Chikkatur, and W. Ketterle, Nature (London) 396, 345 (1998).
• (5) A. E. Leanhardt, Y. Shin, D. Kielpinski, D. E. Pritchard, and W. Ketterle, Phys. Rev. Lett. 90, 140403 (2003).
• (6) J. Ieda, T. Miyakawa, and M. Wadati, Phys. Rev. Lett. 93, 194102 (2004); J. Phys. Soc. Jpn. 73, 2996 (2004).
• (7) L. Li, Z. Li, B. A. Malomed, D. Mihalache, and W. M. Liu, Phys. Rev. A 72, 033611 (2005).
• (8) W. Zhang, Ö. E. Müstecaplioglu, and L. You, Phys. Rev. A 75, 043601 (2007).
• (9) M. Uchiyama, J. Ieda, and M. Wadati, J. Phys. Soc. Jpn. 75, 064002 (2006).
• (10) B. J. Dabrowska-Wüster, E. A. Ostrovskaya, T. J. Alexander, and Y. S. Kivshar, Phys. Rev. A 75, 023617 (2007).
• (11) Yu. S. Kivshar and G. P. Agrawal, Optical Solitons: From Fibers to Photonic Crystals (Academic Press, 2003).
• (12) B. Luther-Davies and X. Yang, Opt. Lett. 17, 496 (1992).
• (13) Th. Busch and J. R. Anglin, Phys. Rev. Lett. 87, 010401 (2001).
• (14) E. G. M. van Kempen, S. J. J. M. F. Kokkelmans, D. J. Heinzen, and B. J. Verhaar, Phys. Rev. Lett. 88, 093201 (2002).
• (15) N. N. Klausen, J. L. Bohn, and C. H. Greene, Phys. Rev. A 64, 053602 (2001).
• (16) S. V. Manakov, Zh. Eksp. Teor. Fiz. 65, 505 (1973) [Sov. Phys. JETP 38, 248 (1974)].
• (17) V. E. Zakharov and S. V. Manakov, Zh. Eksp. Teor. Fiz. 71, 203 (1976) [Sov. Phys. JETP 42, 842 (1976)]; V. E. Zakharov and E. I. Schulman, Physica D 4, 270 (1982); V. G. Makhankov and O. K. Pashaev, Theor. Math. Phys. 53, 979 (1982).
• (18) N. Yajima and M. Oikawa, Progr. Theor. Phys. 56, 1719 (1976).
• (19) Y. S. Kivshar, Opt. Lett. 17, 1322 (1992).
• (20) M. Aguero, D. J. Frantzeskakis, and P. G. Kevrekidis, J. Phys. A: Math. Gen. 39, 7705 (2006).
• (21) N. P. Robins, W. Zhang, E. A. Ostrovskaya, and Y. S. Kivshar, Phys. Rev. A 64, 021601(R) (2001).
• (22) Th. Busch and J. R. Anglin, Phys. Rev. Lett. 84, 2298 (2000); D. J. Frantzeskakis, G. Theocharis, F. K. Diakonos, P. Schmelcher, and Yu. S. Kivshar, Phys. Rev. A 66, 053608 (2002); V. V. Konotop and L. Pitaevskii, Phys. Rev. Lett. 93, 240403 (2004); G. Theocharis, P. Schmelcher, M. K. Oberthaler, P. G. Kevrekidis, and D. J. Frantzeskakis, Phys. Rev. A 72, 023609 (2005); D. E. Pelinovsky, D. J. Frantzeskakis, and P. G. Kevrekidis, Phys. Rev. E 72, 016615 (2005).
• (23) U. Al Khawaja, H. T. C. Stoof, R. G. Hulet, K. E. Strecker, and G. B. Partridge, Phys. Rev. Lett. 89, 200404 (2002); P. G. Kevrekidis, D. J. Frantzeskakis, R. Carretero-González, B. A. Malomed, G. Herring, and A. R. Bishop, Phys. Rev. A 71, 023614 (2005).
• (24) W. Kohn, Phys. Rev. 123, 1242 (1961); J. F. Dobson, Phys. Rev. Lett. 73, 2244 (1994).
• (25) M. Moreno-Cardoner, J. Mur-Petit, M. Guilleumas, A. Polls, A. Sanpera, and M. Lewenstein, Phys. Rev. Lett. 99, 020404 (2007); J. Mur-Petit, M. Guilleumas, A. Polls, A. Sanpera, and M. Lewenstein, K. Bongs, and K. Sengstock, Phys. Rev. A. 73, 013629 (2006).
• (26) V. M. Pérez-García, H. Michinel, and H. Herrero, Phys. Rev. A 57, 3837 (1998).
• (27) T.-L. Ho, Phys. Rev. Lett. 81, 742 (1998).
• (28) T. Ohmi and K. Machida, J. Phys. Soc. Jpn. 67, 1822 (1998).
• (29) H. Pu, C. K. Law, S. Raghavan, J. H. Eberly, N. P. Bigelow, Phys. Rev. A 60, 1463 (1999); H. Pu, S. Raghavan, and N. P. Bigelow, Phys. Rev. A 61, 023602 (2000).
• (30) S. Yi, Ö. E. Müstecaplioglu, C. P. Sun, and L. You, Phys. Rev. A 66, 011601(R) (2002).
• (31) L. P. Pitaevskii and S. Stringari, Bose-Einstein Condensation (Oxford University Press, Oxford, 2003).
• (32) G. Huang, J. Szeftel, and S. Zhu, Phys. Rev. A 65, 053605 (2002).
• (33) D. J. Frantzeskakis, N. P. Proukakis, and P. G. Kevrekidis, Phys. Rev. A 70, 015601 (2004).
• (34) D. J. Frantzeskakis, P. G. Kevrekidis, and N. P. Proukakis, Phys. Lett. A 364, 129 (2007).
• (35) Y. S. Kivshar and B. A. Malomed, Rev. Mod. Phys. 61, 763 (1989).
• (36) L. D. Landau and E. M. Lifshitz, Fluid Mechanics (Pergamon, New York, 1959).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters