Phase Diagram, Stability and Magnetic Properties of Nonlinear Excitations in Spinor Bose-Einstein Condensates
G. C. Katsimiga, S. I. Mistakidis, P. Schmelcher, P. G. Kevrekidis
PPhase Diagram, Stability and Magnetic Properties of Nonlinear Excitationsin Spinor Bose-Einstein Condensates
G. C. Katsimiga, S. I. Mistakidis, P. Schmelcher,
1, 2 and P. G. Kevrekidis Center for Optical Quantum Technologies, Department of Physics,University of Hamburg, Luruper Chaussee 149, 22761 Hamburg, Germany The Hamburg Centre for Ultrafast Imaging, University of Hamburg,Luruper Chaussee 149, 22761 Hamburg, Germany Department of Mathematics and Statistics, University of Massachusetts Amherst, Amherst, MA 01003-4515, USA (Dated: August 4, 2020)We present the phase diagram, the underlying stability and magnetic properties as well as thedynamics of nonlinear solitary wave excitations arising in the distinct phases of a harmonically con-fined spinor F = 1 Bose-Einstein condensate. Particularly, it is found that nonlinear excitationsin the form of dark-dark-bright solitons exist in the antiferomagnetic and in the easy-axis phaseof a spinor gas, being generally unstable in the former while possessing stability intervals in thelatter phase. Dark-bright-bright solitons can be realized in the polar and the easy-plane phases asunstable and stable configurations respectively; the latter phase can also feature stable dark-dark-dark solitons. Importantly, the persistence of these types of states upon transitioning, by meansof tuning the quadratic Zeeman coefficient from one phase to the other is unraveled. Additionally,the spin-mixing dynamics of stable and unstable matter waves is analyzed, revealing among oth-ers the coherent evolution of magnetic dark-bright, nematic dark-bright-bright and dark-dark-darksolitons. Moreover, for the unstable cases unmagnetized or magnetic droplet-like configurations andspin-waves consisting of regular and magnetic solitons are seen to dynamically emerge remainingthereafter robust while propagating for extremely large evolution times. Our investigations pave thewave for a systematic production and analysis involving spin transfer processes of such waveformswhich have been recently realized in ultracold experiments. I. INTRODUCTION
Ultracold atoms constitute ideal platforms for inves-tigating the nonlinear behavior of quantum many- bodysystems due to their high degree of controllability andisolation from the environment [1–3]. A principal ex-ample has been the exploration of dark and bright soli-tons and their dynamical manifestations, as well as theirmulti-dimensional and multi-component extensions inBose-Einstein condensates (BECs) [4, 5]. Indeed, avariety and admixtures of these types of excitationsare nowadays known to exist in scalar [6–11], pseudo-spinor [12–17] and spinor [18–22] BECs. Importantly,in recent years, many of these works have featured ex-perimental realizations of such excitations. However,up to now the majority of both theoretical and experi-mental endeavors has been mainly focused on studyingsolitons in single and pseudo-spinor BEC systems andalso within the so-called Manakov limit [23]. The latterassumes the intra- and the inter-species coupling to beon equal footing. As such the physics of nonlinear exci-tations outside this limit is less explored although thereis an ongoing theoretical effort in this direction over thepast few years [24–28].Arguably, even less explored appears to be the con-nection between regular e.g. vector solitons and mag-netic solitons or higher spin objects such as F = 1spinors and spin-waves. Namely, nonlinear structuresfor which the magnetic interactions between the speciesare a crucial component. For instance, a magnetic soli-ton typically residing in a spin balanced density back- ground [29] is characterized by a localized spin magne-tization, measured as the difference between the pop-ulation of the participating components. Such non-linear polarization waves have also been studied ear-lier [30] in binary BECs for parametric variations lyingoutside the Manakov limit both in the absence and inthe presence of a Rabi coupling between the ensuingcomponents [31]. Case examples of magnetic solitonsare dark and dark-bright (DB) matter waves that differfrom their regular or standard counterparts in a two-foldmanner: (i) they exist for unequal intra and interspincouplings and (ii) their width scales according to thespin-healing length [29]. They can also have the formof dark-antidark solitons, with the latter being densitybumps on top of the BEC background, that have beenvery recently experimentally monitored [32–35]. F = 1 spinor BECs offer the possibility for study-ing not only regular solitons but also magnetic onesand admixtures thereof. In particular, owing to the farricher phase diagram exhibited by such gases [36] (see,also, [37] for a recent discussion and [38] for the impactof many-body effects) already several works have beendevoted in studying a variety of nonlinear excitationsthat arise in them [39–44]. These include for instancespin domains [45, 46], spin textures [47, 48], the very re-cently experimentally observed dark-dark-bright (DDB)and dark-bright-bright (DBB) solitons [49] (and vari-ants [50], as well as interactions [22] thereof) and eventwisted magnetic solitons [51].In the present work we exploit this momentum andby following recent experiments [49] we first map out a r X i v : . [ n li n . PS ] A ug the phase diagram of nonlinear excitations that arisein a one-dimensional (1D) harmonically confined spinor F = 1 BEC. More precisely, DDB, DBB and dark-dark-dark (DDD) solitons constitute the principal ingredientsof this phase diagram. DDB solutions exist in the anti-ferromagnetic (AF) and the easy-axis (EA) phase, DBBsolitons arise in the polar (PO) and the easy-plane (EP)phases and DDD waves are realized in the EP phase too.Subsequently, we exhaustively study the stability prop-erties of the solitonic solutions involved in each phaseand read out the magnetic properties of the participat-ing solitons. Finally, we consider transitions betweenthe distinct phases of the spinor system, by means ofvarying the quadratic Zeeman energy shift, and studydeformations of the above-identified matter waves. Dy-namical evolution of stable and unstable configurationsreveals among others: the coherent evolution of mag-netic DB solitons and changes in the magnetic proper-ties of the evolved entities including the formation ofcomposite spin objects. The latter are composed ofregular solitons and spin-waves. Additionally we ob-serve metastable states evolving into periodically recur-ring unmagnetized Thomas-Fermi (TF)-droplet config-urations and also magnetized entities with droplets oc-cupying the symmetric spin sublevels –with a domainwall (DW) separating them– and a localized wavefunc-tion hosted in the remaining spin-component. The lat-ter nearly periodic structures closely resemble magnondrops [52, 53], while in both cases domain-walls (DW)are imprinted in the local magnetization.Our work is structured as follows. In section IIthe relevant mean-field theoretical framework is intro-duced. The ground state phase diagram of a harmon-ically trapped 1D spin-1 BEC is initially discussed inSection III and we then proceed to the presentation andsystematic exploration of the relevant phase diagram ofnonlinear excitations in the form of DDD, DDB andDBB solitons. Section IV addresses the existence, thestability properties, by means of Bogoliubov de-Gennes(BdG) linearization analysis, and subsequently the dy-namics of the different solitonic waveforms that arisein the distinct phases of the spinor system. Finally, insection V we summarize our findings and also providefuture perspectives. II. SPINOR SETUP AND MAGNETIZATIONMEASURES
A spin-1 BEC composed of the magnetic sublevels m F = 0 , ± F = 1, either of a Rb [49] or a Na [54] atom gas being confined in a1D harmonic trap is considered. A cigar-shaped geom-etry is employed that has been very recently realizedexperimentally [49] utilizing a highly anisotropic trapwith the longitudinal and transverse trapping frequen-cies obeying ω x (cid:28) ω ⊥ . In the mean-field frameworkthe dynamics of such a spinor system can be described by the following coupled dimensionless Gross-Pitaevskiiequations (GPEs) of motion [40, 42, 55, 56] i∂ t Ψ = H Ψ + c (cid:0) | Ψ +1 | + | Ψ | + | Ψ − | (cid:1) Ψ + c (cid:0) | Ψ +1 | + | Ψ − | (cid:1) Ψ + 2 c Ψ +1 Ψ ∗ Ψ − , (1)for the m F = 0 magnetic sublevel, while the symmetric m F = ± i∂ t Ψ ± = H Ψ ± + c (cid:0) | Ψ +1 | + | Ψ | + | Ψ − | (cid:1) Ψ ± + c (cid:0) | Ψ ± | + | Ψ | − | Ψ ∓ | (cid:1) Ψ ± + q Ψ ± + c Ψ ∗∓ Ψ . (2)In Eqs. (1)-(2), Ψ m F ( x, t ) denotes the wavefunctionof the | F = 1 , m F = 0 (cid:105) and | F = 1 , m F = ± (cid:105) spin-components respectively. The single particle Hamilto-nian term is H ≡ − ∂ x + V ( x ), with V ( x ) = Ω x being the 1D harmonic potential. Here, Ω ≡ ω x /ω ⊥ plays the role of the longitudinal over the transversetrapping frequency and is typically a small parameteri.e., Ω (cid:28) q denotes the quadraticZeeman energy shift parameter that leads to an effectivedetuning of the m F = ± m F = 0 one. It is quadratically proportional toan external magnetic field applied along the spin- z direc-tion [36, 57] and can be experimentally tuned by eitheradjusting the applied magnetic field [58] or by using amicrowave dressing field [59, 60].Moreover, c and c are the so-called spin-independent and spin-dependent interaction coeffi-cients. The former accounts for attractive (repulsive)interatomic interactions upon taking negative (positive)values and the latter is positive ( c >
0) for antiferro-magnetic and negative ( c <
0) for ferromagnetic inter-actions. Both c and c are expressed in terms of the s -wave scattering lengths a and a , accounting for twoatoms in the scattering channels with total spin F = 0and F = 2 respectively, via the relations c = ( a +2 a )3 a ⊥ and c = ( a − a )3 a ⊥ [21, 40]. Here, a ⊥ = (cid:112) (cid:126) /M ω ⊥ is thetransverse harmonic oscillator length with M denotingthe mass, e.g., of a Rb atom. Eqs. (1)-(2) have beenmade dimensionless by measuring length, energy andtime in units of (cid:112) (cid:126) / ( M ω ⊥ ), (cid:126) ω ⊥ and ω − ⊥ respectively.Consequently, the corresponding interaction strengthsare expressed in terms of (cid:112) (cid:126) ω ⊥ /M . In the adoptedunits and for a ferromagnetic, i.e., c <
0, spinorBEC of Rb atoms [21, 37], the experimentally mea-sured spin-dependent and spin-independent couplingsalso used herein are c ≈ − × − (cid:112) (cid:126) ω ⊥ /M and c = 1 (cid:112) (cid:126) ω ⊥ /M .Additionally, the population of each spin-componentis defined as n m F = 1 N (cid:90) dx | Ψ m F | , m F = 0 , ± . (3)Here, N = (cid:80) m F (cid:82) dx | Ψ m F | denotes the total numberof particles that is a conserved quantity for the spino-rial system of Eqs. (1)-(2). Evidently, 0 ≤ n m F ≤ | Ψ +1 | | Ψ − | | Ψ | c or q th DDD DBB q Polaror qq th Easy−Planeor orDDB c DBB q q q Easy−AxisAntiferromagneticDDB q (a) PolarEasy−Axis Easy−Plane DBAntiferromagnetic (b) D FIG. 1. Schematic illustration of (a) the ground state phase diagram of a harmonically trapped spin-1 BEC in the ( c , q )plane. (b) The corresponding phase diagram of nonlinear excitations having the form of DDD, DDB and DBB solitons. Inboth cases characteristic density profiles, | Ψ m F | with m F = 0 , ±
1, of each of the phases that can be realized in such asystem are provided (see legend). For c > q < q > c < q <
0, in the EP phase if 0 < q < q th ≡ n | c | ≈ .
017 and in the PO oneif q ≥ q th . In both phase diagrams solid black lines mark the individual phase transition boundaries. Vertical dashed blacklines designate the regions from where and on the distinct configurations deform, i.e., the DDB into DB ( q ≈ − . q ≈ . q ≈ . q ≈ − .
016 ( q ≈ − . µ , ± = 2. satisfied. Furthermore, in order to quantify first- andsecond-order transitions between the distinct phases ofthe spin-1 BEC system as well as to monitor the mag-netic properties of the emergent nonlinear excitationsduring evolution we utilize the magnetization along thespin- z -axis that reads M z = 1 N (cid:90) dx (cid:0) | Ψ +1 | − | Ψ − | (cid:1) . (4) M z essentially measures the population imbalance be-tween the symmetric m F = ± − ≤ M z ≤
1. For instance, a fully magnetized state alongthe + z or − z spin direction corresponds to M z = +1 or M z = − m F = 0 and the m F = ± P = 1 N (cid:90) dx (cid:2) | Ψ | − (cid:0) | Ψ +1 | + | Ψ − | (cid:1)(cid:3) . (5)It can be easily deduced that − ≤ P ≤
1. As wewill unveil later on, P also accounts for alterations inthe magnetic properties of the spinor system and al-lows us to distinguish among fully magnetized and un-magnetized spin configurations as we cross, by means of varying the quadratic Zeeman energy shift q , a phasetransition boundary.Finally, for the numerical investigations that follow,the trapping frequency is fixed to Ω = 0 .
1, but we notethat the results presented herein are not altered even fortrapping frequencies of the order of Ω = 0 .
01 that areused in recent spin-1 BEC experiments [49]. Only slightdeviations of the corresponding transition boundariesare observed. For instance, for the EP to PO transition,while q th ≡ n | c | = 0 .
02 (with n denoting the peakdensity) for Ω = 0 .
01 it is q th ≈ .
017 for Ω = 0 .
1. Addi-tionally, c = 1, c = ± × − and we choose the chem-ical potentials of the different components µ , ± = 2.It is also important to mention that we have checkedthat the results to be presented below are robust alsofor c = 1 and c = 3 . × − , namely for the experi-mentally relevant interaction coefficient parameter ratiocorresponding to Na gas and also for larger chemicalpotentials, i.e. µ , ± = 3 and µ , ± = 5. To access thedistinct phases of the spinor system, we typically vary q within the intervals [ − . , .
5] and [ − . , . dx = 0 . dt = 0 .
001 respectively. Our numerical computa-tions are restricted to a finite region by employing hard-wall boundary conditions. Particularly, in the dimen-sionless units adopted herein, the hard-walls are locatedat x ± = ±
80 and we do not observe any appreciabledensity for | x | > III. PHASE DIAGRAM OF NONLINEAREXCITATIONS
Before delving into the details of the phase diagramof nonlinear excitations in the form of DDD, DDB andDBB solitons that arise in spin-1 BECs, we first brieflyrevisit the relevant ground state phase diagram of a har-monically confined 1D spin-1 BEC [37]. This descriptionwill enable us to qualitatively expose the effect of em-bedding nonlinear structures into the different magneticphases.
A. Ground state phase diagram
A schematic representation of the ground state phasediagram is illustrated in Fig. 1(a). As it has been re-cently demonstrated [36–38] different phases can be real-ized for such a confined spin-1 system. They stem fromthe interplay between the sign of the spin-dependent in-teraction coefficient c and the strength of the quadraticZeeman term q . Specifically, for c > q < m F = ± M z = 0 and P = −
1. A first order phase tran-sition [62, 63] separates this phase from the PO onethat can be reached upon increasing q . The transi-tion point appears at q = 0 and the PO phase is char-acterized again by an unmagnetized ground state butwith all atoms populating the m F = 0 spin-component.Therefore M z = 0 and P = 1. On the other hand,for c < q < z or the − z spin-direction, i.e. either the m F = +1 or m F = − M z = +1or M z = − P = −
1. Upon increasing q a second-order phase transition occurs at q = 0 andfor 0 < q < q th ≈ .
017 the system enters the EP phasewith its ground state having all three m F componentspopulated. Particularly here, M z = 0 reflecting the factthat the m F = ± P ∈ ( − , q > q th ≈ .
017 yet an-other second order phase transition takes place whichleads to an unmagnetized ground state having only the m F = 0 spin component populated, i.e. M z = 0 and P = 1. In this case, once more the PO phase is reached. −1.5 −1 −0.5 0 0.5−1−0.8−0.6−0.4−0.200.20.40.60.81 q P M z −1.5 −1 −0.5 0 0.500.51 q n n +1 n − −1.5−1−0.500.500.51 w w +1 w − −20 0 2000.511.522.5 x −20 0 2000.511.52 x −20 0 2000.511.52 x | Ψ | | Ψ +1 | | Ψ − | | Ψ tot | −1.5 −1 −0.5 0 0.50102030 q (a) (b)(c) FIG. 2. (a) Polarization P , and magnetization, M z , of thespin-1 system for a DDB state existing within the AF phaseupon increasing the quadratic Zeeman coefficient q in orderto enter the PO phase. Vertical dashed-dotted (green) linemarks this transition. Insets from bottom left to top rightillustrate characteristic soliton profiles for q = − . q = − .
015 and q = 0 .
1. (b) Population, n m F ( q ), and (c) solitonwidth, w m F ( q ), with m F = 0 , ±
1, of the DDB solitons asa function of q (see legends). Other parameters used areΩ = 0 . µ , ± = 2, c = 5 × − , and c = 1. Note here, that in order to obtain the above-discussedground state phase diagram TF profiles are employed asinitial guesses, within our fixed point algortithm [61], forthe distinct m F = 0 , ± m F ( x, t = 0) = (cid:113) c − (cid:0) µ m F − V ( x ) (cid:1) . (6)When this expression is used here and below, it is im-plied to be valid when the quantity under the radical isnon-negative (and the relevant wavefunction is paddedwith zeros outside that region). In Eq. (6), µ m F denotesthe chemical potential of each spin-component while astationary state satisfies the phase matching condition µ = ( µ +1 + µ − ) / B. Phase diagram of solitonic excitations
In order to unravel the phase diagram of nonlinearexcitations depicted in Fig. 1(b), the spin-1 system isinitialized in each of the above-identified phases embed-ding dark and bright solitons as wavefunctions for eachof the m F = 0 , ± D ( x, t = 0) = (cid:113) c − [ µ m F − V ( x )] tanh ( Dx ) , (7)Ψ B ( x, t = 0) = η sech ( Dx ) . (8)In the above expressions Ψ D ( x ) and Ψ B ( x ) denote thewavefunctions utilized for a dark and a bright solitonconfiguration respectively. In Eq. (7) the quantity underthe square root denotes the customary used TF back-ground needed for dark solitons to be embedded on. −1 −0.5 000.20.40.60.81 q R e ( ω ) −1 −0.5 000.020.040.060.080.1 q I m ( ω ) (a) (b) FIG. 3. BdG spectrum of DDB stationary states identifiedin the AF phase as q increases towards the PO phase. (a)Real part, Re( ω ), of the corresponding eigenfrequencies, ω ,as a function of the quadratic Zeeman coefficient q . Thetrajectories of the two anomalous modes (AMs) present inthis spectrum are indicated by the light blue dots (see text).(b) Imaginary part, Im( ω ), of the respective eigenfrequen-cies. The eigenfrequency zero crossings occur at q = 0.The remaining system parameters correspond to Ω = 0 . µ , ± = 2, c = 5 × − , and c = 1. Moreover, D and η refer, respectively, to the commoninverse width considered for each spinorial soliton com-ponent and the amplitude of the bright soliton configu-ration (see our detailed discussion in Sec. IV).It is found that DDB solitons, being unmagnetizedconfigurations, exist within the AF phase for all valuesof the quadratic Zeeman energy shift lying within theinterval q ∈ ( − . , m F = 0 spin-component. This trapping mechanism becomes progres-sively less effective. Namely, as q increases towards thephase transition point ( q = 0) the bright soliton grad-ually becomes the dominant configuration before mor-phing into a TF one. We remark that the existence ofa DDB soliton in the AF phase already presents funda-mental deviations from its ground state properties. In-deed, in the latter case the m F = 0 magnetic sublevel isunpopulated (of course also the m F = ± q ∈ [0 , . m F = 0 component for q > . q so as to enter the AF phase until a critical value ofthe quadratic Zeeman energy shift, i.e. q ≈ − . m F = 0 state is occupied. Turning to c <
0, station-ary solutions of the DDB type are realized within theferromagnetic EA phase existing within the parametric x −20020 x −20020 t x x −20020 x −20020 t x t n n +1 n − t n n +1 n − | Ψ +1 ( x, t ) + u AM | q = 0 | Ψ ( x, t ) + u AM | | Ψ − ( x, t ) + u AM | (c)(b)(a) (d)(e)(f) q = − . | Ψ − ( x, t ) + u AM | | Ψ +1 ( x, t ) + u AM | | Ψ ( x, t ) + u AM | (h)(g) FIG. 4. (a)-(f) Dynamical evolution of the density, | Ψ m F + u AM | , with m F = 0 , ± u AM , associated with thelowest lying anomalous mode (AM) appearing in the BdGspectrum of Fig. 3. Evolution of the perturbed DDB statefor (a)-(c) q = − .
01 and (d)-(f) q = 0. (g), (h) Temporalevolution of the populations, n m F ( t ), for the above selectionof q ’s. Other parameters used are Ω = 0 . µ , ± = 2, c =5 × − , and c = 1. region q ∈ ( − . , . q decreases further into fully magnetized, i.e. M z = +1( M z = − m F = 0 and m F = +1 ( m F = 0 and m F = −
1) components. Onceagain this is far from the ground state of the EA featur-ing only m F = +1 (or m F = −
1) populations. Movingto q >
0, namely entering the EP phase, two types ofsolitonic solutions are found to exist for the spin-1 sys-tem. These excitations can have the form of unmagne-tized spinor DBB or DDD solitons, a result that is per-mitted by the relevant ground state where all three mag-netic components are occupied. The former solitonic en-tities appear to be significantly broader when comparedto the more localized DDD configurations and becomehighly localized as we enter the PO phase. Recall thatthe PO ground state supports population only in the m F = 0 magnetic sublevel. The transition point for theDBB configuration appears at q ≈ .
994 while it oc-curs significantly earlier, q ≈ . q , in order toenter the EA phase reveals that the DBB configurationdeforms fast, around q ≈ − . m F = ± q ≈ − .
016 before their transitioning towards twodarks that occupy the symmetric m F = ± −0.5 0 0.5 1 1.5−1−0.8−0.6−0.4−0.200.20.40.60.81 q P M z −0.5 0 0.5 1 1.500.51 q n n +1 n − −20 0 2000.511.52 x | Ψ | | Ψ +1 | | Ψ − | | Ψ tot | −20 0 2000.511.52 x −20 0 2000.511.52 x −0.500.511.500.51 q w w +1 w − −0.5 0 0.5 1 1.50102030 q (b)(a) (c) −20 0 20012 x FIG. 5. (a) Polarization P , and magnetization, M z , for DBBsolutions existing within the PO phase upon decreasing thequadratic Zeeman coefficient q towards the AF phase. Insetsfrom bottom left to top right illustrate characteristic solitonprofiles for q = − . q = − . q = 0 . q = 1.(b) Populations, n m F ( q ) and (c) soliton width, w m F , forvarying q with m F = 0 , ± q ≈ − .
005 and q ≈ .
994 mark theboundaries of deformation of the DBB for different q ’s. In allcases the remaining parameters refer to Ω = 0 . µ , ± = 2, c = 5 × − , and c = 1. IV. STABILITY ANALYSIS AND DYNAMICSOF SPINOR SOLITONS
Our aim in what follows is not only to illustrate theexistence of stationary spinor solitons of the DDD, DDBand DBB type existing in a 1D harmonically confinedspin-1 BEC composed e.g. of Rb atoms and obeyingEqs. (1)-(2), but also to systematically investigate theirstability properties. We remark that for ferromagnetic(AF) BECs we consider c = − × − ( c = 5 × − ) as representative example and vary the quadraticZeeman energy shift to access the underlying magneticphases. A. Antiferromagnetic DDB matter waves
For instance, in order to infer about the existence ofDDB solitons within the AF phase shown in the phasediagram of Fig. 1(b), the matter wave dark solitons ofEq. (7) are embedded as initial guesses for the m F = ± m F = 0 spin state. Employing the aboveansatz, and using the iterative scheme discussed above,DDB stationary states are found within the AF phase,i.e. for c = 5 × − and for values of q ∈ ( − . , | Ψ , ± | , are pre-sented as insets in Fig. 2(a). However, upon increasing q towards the transition point ( q = 0) above which the POphase is realized, the DDB solitons deform into stateswhere the bright structure in the m F = 0 componentoverfills/dominates the dark wells. Also the total den- q R e ( ω ) (a) q I m ( ω ) (b) FIG. 6. BdG spectrum of DBB stationary states in the POphase upon decreasing q towards the AF phase. (a) Realpart, Re( ω ), of the corresponding eigenfrequencies as a func-tion of q . The trajectory of the AM appearing in this spec-trum is indicated by light blue dots (see text). (b) Imaginarypart, Im( ω ), of the respective eigenfrequencies. The desta-bilization of the DBB state occurs at q cr ≈ − .
004 while for q ≥ m F = 0 component.Other parameters used are Ω = 0 . µ , ± = 2, c = 5 × − ,and c = 1. sity, | Ψ tot | , of the spinor system exhibits a TF profileinstead of the dark-shaped density appearing deep inthe AF phase. This altered nature of the DDB con-figuration, which remains unmagnetized ( M z = 0) forall values of q , is naturally accompanied by a change inthe polarization of this configuration. The DDB soli-tons possess P = − q ≤ − . m F = ± − < P ≤ q ≈ − .
02 all threecomponents are equally populated having significantlywider [65] stationary states [Fig. 2(c)] as compared tothe ones for larger negative q values. This broadeningsuggests that the DDB character of the relevant states islost. Importantly, at q = 0 an abrupt population trans-fer to the m F = 0 component [Fig. 2(b)] associatedwith the drastic deformation of this latter configurationto the ground state of the PO phase manifests itself; seee.g. the right uppermost inset of Fig. 2(a).In order to extract the stability properties of theaforementioned DDB stationary states (as well as forthe DBB and DDD solitons to be presented below), alinear stability or BdG analysis is performed. The latterconsists of perturbing the iteratively identified in eachphase stationary solutions Ψ m F ( x ) (with m F = 0 , ± m F ( x, t ) = Ψ m F ( x ) + (cid:15) (cid:16) a m F ( x ) e − iωt + b ∗ m F ( x ) e iω ∗ t (cid:17) . (9)By inserting this ansatz into the system of Eqs. (1)-(2)and linearizing with respect to the small amplitude pa- q R e ( ω ) q I m ( ω ) (a) (b) FIG. 7. BdG spectrum of PO dark solitons for varying q .(a) Real part, Re( ω ), of the underlying eigenfrequencies asa function of q . The trajectory of the AM is illustrated bylight blue dots (see text). (b) Imaginary part, Im( ω ), ofthe respective eigenfrequencies. Notice that the dark wavedestabilizes for q ∈ [0 . ,
1] coexisting in a subcritical pitch-fork bifurcation with the DBB solitons identified also withinthe PO phase. Dashed green line in (b) is used as a guideto the eye for the ensuing bifurcation. The remaining sys-tem parameters are Ω = 0 . µ , ± = 2, c = 5 × − , and c = 1. rameter (cid:15) leads to an eigenvalue problem for the eigen-frequencies ω , or equivalently eigenvalues λ ≡ iω , andeigenfunctions ( a , b , a +1 , b +1 , a − , b − ) T that is solvednumerically. For further details on the BdG analysis werefer the reader to Refs. [5, 66, 67]. Due to the gen-erally complex nature of the ensuing eigenfrequencies,it becomes apparent that the following possibilities canarise: if modes with purely real eigenvalues or equiva-lently imaginary eigenfrequencies or complex eigenval-ues/eigenfrequencies are identified, these are responsiblefor the existence of an instability [66]. Moreover, dueto the Hamiltonian structure of the system investigatedherein, quartets of such eigenfrequencies can occur [67].Namely if ω is an eigenfrequency so are − ω and ± ω ∗ .As such, if Im( ω ) (cid:54) = 0, then there will always exist amode leading to the growth and eventual deformationof the examined in each phase solitonic configuration.The BdG analysis outcome for the DDB soliton so-lutions is shown in Fig. 3(a), (b). DDB solitons con-stitute excited states of the spin-1 system, exactly liketheir two-component dark-bright analogue [26], a fea-ture that is reflected in their linearization spectra viathe emergence of the so-called anomalous modes (AMs).These eigenstates are quantified via the negative energyor negative Krein signature [67] defined for the spinorsystem as K = Ω (cid:90) (cid:16) | a | − | b | + | a +1 | − | b +1 | + | a − | − | b − | (cid:17) dx. (10)The existence of these modes is central to our stabil- −505 −505 t t t x x −20020 x −200200 2000 400000.51 t t t n n +1 n − n n +1 n − n n +1 n − q = 0 . | Ψ ( x, t ) + u AM | (l)(j) q = 0 . q = − . (k) | Ψ − ( x, t ) + u AM | | Ψ +1 ( x, t ) + u AM | | Ψ ( x, t ) + u AM | | Ψ +1 ( x, t ) | | Ψ − ( x, t ) | | Ψ ( x, t ) | | Ψ − ( x, t ) + u AM | | Ψ +1 ( x, t ) + u AM | (a) (d) (g)(b) (e) (h)(i)(f)(c) FIG. 8. Density evolution of (a)-(c) | Ψ m F ( x, t ) | referring toa deformed DBB soliton in the AF phase and (d)-(i) | Ψ m F + u AM | of a DBB soliton within the PO phase but beingexcited by the eigenvector associated with the AM appearingin the BdG spectrum of Fig. 6. The different columns depictthe evolution of the DBB state for q = − . q = 0 . q = 0 . n m F ( t ), for the aforementionedvalues of q . In all cases Ω = 0 . µ , ± = 2, c = 5 × − and c = 1. ity analysis since their potential collision with posi-tive Krein signature modes can give rise to stability-changing events in the form of oscillatory instabilitiesor Hamiltonian-Hopf bifurcations [67]. In the presentsetting the AMs are denoted by light blue dots. TheDDB solution possesses, due to the presence of two darksolitons [68], two such modes [Fig. 3(a)] that cross theorigin of the spectral plane at q cr = 0 signalling thedestabilization of the DDB wave [Fig. 3(b)]. Interest-ingly, and also for all values of q ∈ ( − . , q closer to the critical point,defining the AF to PO transition boundary, the desta-bilization of both modes leads, irrespectively of whichmode we excite (namely the first lower-lying one or thesecond), to a breathing motion of the DDB configura-tion. Its response is visualized in the spatio-temporalevolution of the densities, | Ψ m F ( x, t ) + u AM | , pre-sented in Fig. 4(a)-(c) entailing both the particle-likeoscillation of the DDB soliton but predominantly the −1.5 −1 −0.5 0 0.5−1−0.8−0.6−0.4−0.200.20.40.60.81 q P M z −1.5 −1 −0.5 0 0.500.51 q n n +1 n − −1.5−1−0.500.500.51 w w +1 w − −0.05 0 0.0500.10.2 q −20 0 20012 x −20 0 2000.511.52 x | Ψ | | Ψ +1 | | Ψ − | | Ψ tot | −1.5 −1 −0.5 0 0.50102030 q (c) (b)(a) −20 0 200123 x FIG. 9. (a) Polarization, P , and magnetization, M z , of thespinor system for a DDB state existing within the EA phaseupon varying the quadratic Zeeman coefficient q so as to en-ter the EP and PO phases. Vertical dashed-dotted (green)lines mark the boundaries of deformation of the DDB as q isvaried (see text). Insets from bottom left to top right illus-trate characteristic wave profiles for q = − . q = − . q = 0 .
1. (b) Populations, n m F ( q ) and (c) soliton width w m F ( q ) as a function of q with m F = 0 , ± . µ , ± = 2, c = − × − , and c = 1. overall breathing of the state. Although the first AMis excited in this case, the dynamics is dominated bythe breathing mode. The nature of this composite mo-tion is also reflected in the irregular oscillation of thepopulation, n m F ( t ), of each m F component illustratedin Fig. 4(g). At the transition/destabilization point theprevailing feature of the perturbed DDB configurationis its breathing as can be seen by monitoring the evolu-tion of the densities illustrated in Fig. 4(d)-(f) togetherwith the coherent oscillation of the relevant populations[Fig. 4(h)]. In both of the aforementioned cases the os-cillatory character of n m F ( t ) implies a weak amplitudespin-mixing dynamics. B. Polar DBB solitons
Next we turn to the PO phase which is character-ized by c = 5 × − and q > q ∈ [0 , . m F = 0 spin-component, while bright solitons given by Eq. (8) areconsidered for the remaining symmetric m F = ± − < P < −0.6 −0.4 −0.2 000.20.40.6 q R e ( ω ) −0.4 −0.2 000.0020.0040.0060.0080.01 q I m ( ω ) (b)(a) FIG. 10. BdG spectrum of a DDB stationary state withinthe EA phase as q increases towards the EP and PO phases.(a) Real part, Re( ω ) and (b) imaginary part, Im( ω ), of thecorresponding eigenfrequencies with respect to q . In (a) thetrajectories of the two emergent AMs are indicated by lightblue dots. In (b) the destabilization windows (i.e. Im( ω ) (cid:54) =0) occur for q ∈ : [ − . , − . − . , − . − . , − . − . , − . , . . µ , ± = 2, c = − × − , and c = 1. In particular, P = 1 for q >
1, i.e. deep in the POphase, and it gradually decreases as q → + all the wayto P = − q <
0. This latter behavior of P reveals inturn that despite the fact that the ground state configu-ration does not support all three m F components to bepopulated this is not the case for the respective nonlin-ear excitations [Fig. 5(b)]. Selected DBB soliton profilesare depicted as insets in Fig. 5(a). From these profilesit can be deduced that these unmagnetized DBB wavesexist not only within the above-provided q interval butalso at ( q = 0) and below ( q ∈ [0 , − . q → + the DBB states deform towards wider configu-rations [Fig. 5(c)] featuring a pronounced bright solitoncomponent that dominates. This dominant bright com-ponent results in turn to a | Ψ tot | that has a TF profileinstead of a tanh-shaped one occurring for values of q well inside the PO phase. For q < − .
005 an abrupttransition leads to a metastable configuration in whichthe m F = ± P = − q > .
994 yet anotherbut gradual this time deformation of the DBB matterwaves towards a dark soliton with maximal polarization( P = +1) occupying the m F = 0 spin state occurs [topright inset of Fig. 5(a)].By investigating the stability properties of the abovesolitonic entities, it is found that two destabilizationpoints exist for the DBB configuration one residing inthe AF phase and one deep in the PO phase. Specif- x −10010 x −10010 t x x −20020 x −20020 t x t | Ψ − ( x, t ) + u AM | | Ψ − ( x, t ) + u AM | (b) | Ψ +1 ( x, t ) + u AM | | Ψ +1 ( x, t ) + u AM | | Ψ ( x, t ) + u AM | | Ψ ( x, t ) + u AM | q = − . (e)(f)(c) (d)(a) | Ψ ( x, t ) + u AM | | Ψ +1 ( x, t ) + u AM | | Ψ − ( x, t ) + u AM | (i)(h)(g) q = − . FIG. 11. (a)-(f) Spatio-temporal evolution of | Ψ m F + u AM i | with m F = 0 , ± i = 1 , q are q = − .
002 and q = − .
08, i.e.lying in [ − . , − . − . , − . . µ , ± = 2, c = − × − , and c = 1. ically for q < q cr ≈ − . q andthereafter. A result that is further supported by the fi-nite growth rate, Im( ω ) (cid:54) = 0, shown for these negativequadratic Zeeman energies in Fig. 6(b). This destabi-lization is related to a composite motion of the DBBstructure in the parabolic trap that we will soon tracein the dynamics. Contrary to the above destabilizationyet another critical point occurs for the DBB solution forpositive values of q . The latter appears at q cr ≈ . q ∈ [0 , . q ∈ [0 . ,
1] and leading to the loop bi-furcation shown in Fig. 7(b). Within this q intervalalso the Krein signature changes sign from negative be-fore the lower bound to positive after the upper bound.This, in turn, means that the dark soliton destabilizesslightly below unity and restabilizes for q >
1. It isin this interval that indeed the above identified DBBsolitons coexist with the single dark ones in a subcrit- t x t t t n n +1 n − −20 0 2000.511.5 x | Ψ ( x, t ) | | Ψ +1 ( x, t ) | | Ψ − ( x, t ) | −20 0 20−0.0200.02 x M z ( x, t ) (a) t (d) | Ψ ( x, t ) | t = 3690 | Ψ +1 ( x, t ) | t (b) | Ψ − ( x, t ) | t (c) t = 3690 (e) (f) FIG. 12. (a)-(c) Spatio-temporal evolution of the densities, | Ψ m F ( x, t ) | (see legends), of a stationary DDB soliton for q = 0 . , . n m F ( t ). In all cases m F = 0 , ± . µ , ± = 2, c = − × − , and c = 1. ical pitchfork bifurcation. Namely, dark solitons existas stable configurations, for q < . q ∈ [0 . , q >
1, the relevant real eigenvalue pair returns to theimaginary axis, restabilizing the relevant dark state.Direct evolution of the above-identified configurationsslightly below and above the phase transition thresh-old at q = 0 reveals that the DBB solitons undergo inboth cases an overall breathing motion [Fig. 8(a)-(c) and(d)-(f) respectively]. Notice the multi-frequency evolu-tion of the DBB stationary state [Fig. 8(d)-(f)] excitedalong its most unstable eigendirection when comparedto the absence of such an excitation [Fig. 8(a)-(c)]. Onecan correspondingly observe the irregular [Fig. 8(k)]versus coherent population transfer between the spin-components [Fig. 8(j)]. In sharp contrast to the abovedynamics for well-defined DBB solitons, i.e. away fromthe transition and the critical point, a well-defined, in-trap oscillation of the perturbed DBB wave is observed[Fig. 8(g)-(h)] for evolution times up to t = 5 × withthe respective populations, n m F ( t ), remaining constantfor all times. C. Easy-Axis symmetry broken DDB solitons
Moving on to the EA phase, i.e., for c = − × − and q <
0, again DDB stationary states are success-fully identified. However, contrary to the DDB solu-tions found in the AF phase here the DDB waves ex-hibit unequally populated m F = ± m F = − q values within theregion of existence, i.e. q ∈ ( − . , . −0.5 0 0.5 1−1−0.8−0.6−0.4−0.200.20.40.60.81 q P M z −20 0 20012 x −20 0 20012 x −20 0 20012 x −0.5 0 0.5 1 1.500.51 q n n +1 n − −0.5 0 0.5 100.51 q w w +1 w − −0.5 0 0.5 1 1.50102030 q (a) (c)(b) | Ψ | | Ψ +1 | | Ψ − | | Ψ tot | FIG. 13. (a) Polarization P , and magnetization, M z , of aDBB wave within the EP phase upon varying q towards theEA and PO phases. Insets from bottom left to top rightillustrate characteristic soliton profiles for q = − . q = − .
02 and q = 0 .
1. (b) Populations, n m F ( q ) and (c) solitonwidth, w m F ( q ), for different values of q , with m F = 0 , ± . µ , ± = 2, c = − × − , and c = 1. configuration. Notably, such waves preserve the sym-metry (i.e., equal population) of the m F = ± q ∈ [ − . , . q = 0) that separates the EA andthe EP phases. For the remaining quadratic Zeemanenergies lying in the aforementioned q interval the sym-metry is partially preserved, i.e. the m F = − < M z < − < P < P = 1 ( M z = 0) for q > . q < P = − M z = +1) for q < −
1, i.e.deep in the EA phase [Fig. 9(a)]. These symmetric DDBsolitons are fundamentally different (structurally) thanthe relevant ground state configuration in this paramet-ric regime. The latter, according to the phase diagramof Fig. 1(a), favors symmetry broken states that arefully magnetized along the + z - or − z -spin direction,i.e. configurations that have either the m F = +1 orthe m F = − q ∈ ( − . , − .
2] the dark soliton of the m F = − w − ( q ), shown in Fig. 9(c). In particular,around q ≈ − .
515 the configuration is deformed to asymmetry broken almost fully magnetized ( M z ≈ +1)DB soliton that exists for q ∈ [ − . , − . m F = +1 and m F = 0 spin components (similarly,of course, there is a state occupying the the m F = − m F = 0 states). Additionally, as the transitionpoint is approached from below, q → − , also the pop- −0.5 0 0.5 1−1−0.8−0.6−0.4−0.200.20.40.60.81 q P M z −0.5 0 0.5 100.51 q n n +1 n − −20 0 2000.511.52 x | Ψ | | Ψ +1 | | Ψ − | | Ψ tot | −20 0 2000.511.52 x −20 0 2000.511.52 x −20 0 2000.511.52 x −0.500.5100.51 q w w +1 w − −0.5 0 0.5 100.511.5 q (a) (c)(b) FIG. 14. (a) Polarization P , and magnetization, M z , ofthe spinor system for a DDD state existing within the EPphase upon varying q to enter the EA and PO phases. Insetsfrom bottom left to top right illustrate characteristic solitonprofiles for q = − . q = − . q = 0 .
003 and q = 0 . q ≥ . m F = 0 spin state, while for q = − .
017 within the EA phaseonly two darks equally populate the m F = ± n m F ( q ), and (c) the width, w m F ( q ), of thesolitons for varying q and with m F = 0 , ± . µ , ± = 2, c = − × − and c = 1. ulation, n ( q ), and width, w ( q ), of the bright compo-nent increases and the DDB solitons deform even fur-ther. Specifically, for q ∈ ( − . , . m F components to be simultaneouslyoccupied, but already at q ≈ .
008 the unmagnetizedground state of the PO phase is reached.The BdG analysis of the above-discussed soliton so-lutions illustrated in Fig. 10(a), (b) reveals that DDBsolitons possess potentially unstable eigendirections (al-though they also possess stability intervals). This re-sult can be inferred by the finite imaginary eigenfre-quencies (or instability growth rates), Im( ω ), shownin Fig. 10(b). The relevant instability intervals forthe DDB wave are [ − . , − . − . , − . − . , − . − . , − . − . ,
0] re-spectively. Notice that the last bifurcation possessesalso the larger instability growth rate that stems froman eigenfrequency zero crossing of both the higher-and the lower-lying AMs appearing at q cr = 0. For − . < q < − .
384 the state remains linearly stablehaving, however, a minuscule m F = − q , the fully magnetized linearlystable DB solitons are present in the spin-1 system. Toconfirm the above stability analysis findings we havemonitored the dynamical evolution of the DDB solu-tions in all of the above-identified instability intervals1and our results can be summarized as follows. Amongthe two modes that appear in the BdG spectrum ofFig. 10(a) the lower one is related to the weak amplitudein trap oscillation of the DDB wave. The higher modeis responsible for the larger in amplitude anti-phase os-cillation of the involved dark solitons. Additionally, itturns out that even when these states are found to be dy-namically unstable they have remarkably long lifetimesthat suggest their potential experimental observation inexisting spinor settings [49]. A case example showcas-ing the particle-like oscillations that a perturbed DDBstationary state undergoes is presented for q = − . t = 5 × . More specifically, by perturbingthe DDB soliton with the eigenvector associated withthe first AM leads to an oscillation of the wave withinthe trap [Fig. 11(a)-(c)]. On the other hand, the sec-ond mode results in the formation of two atomic blobsto which the dark solitons split the entire condensate[Fig. 11(d)-(f)]. These blobs execute an anti-phase oscil-lation alternating across the two dark components andacross the two sides (left and right) of the dark solitarywave in each component [Fig. 11(d)-(f)]. This latteranti-phase oscillation becomes even more pronouncedespecially for the m F = − q decreasesfurther towards the formation of DB solitons that oc-cupy the m F = 0 and m F = +1 magnetic sublevels[Fig. 11(g)-(i)]. Evidently, as q decreases further ande.g. for q = − .
08 illustrated in Fig. 11(g)-(i), the pop-ulation of the m F = − q = 0). For instance here, by moni-toring | Ψ m F ( x, t ) | for quadratic Zeeman energy shiftsthat lie within the last bifurcation [Fig. 10(b)] revealsthat these transient states for evolution times of theorder of t ≈ m F = 0 component. These localized density blobs arenot of permanent character as is evident in Fig. 12(a)but they revive in an almost periodic manner. In ev-ery recurrence event, the corresponding symmetric spin-components bear droplet-like configurations that appearin an alternating fashion either in the m F = +1 or inthe m F = − m F = 0 and m F = ± M z ( x ), e.g. at t = 3690 where this dynami-cally formed state emerges for the first time [Fig. 12(f)], q R e ( ω ) q I m ( ω ) (b)(a) FIG. 15. BdG spectrum of DBB waves existing in the EPphase as we vary q towards both the EA and the PO phases.(a) Real part, Re( ω ), of the corresponding eigenfrequenciesas a function of the quadratic Zeeman coefficient q . Thetrajectory of the AM present in this spectrum is shown bylight blue dots (see text). (b) Imaginary part, Im( ω ), ofthe respective eigenfrequencies. The eigenvalue zero crossingoccurs at q cr = 0 leading to Im( ω ) (cid:54) = 0. Solid (black) anddashed (green) lines are used as a guide to the eye. Otherparameters used are Ω = 0 . µ , ± = 2, c = − × − and c = 1. reveals that such a configuration bears indeed a DWcharacter across which M z ( x ) changes sign [64]. Such amagnetic entity holds close similarities to the so-calledmagnon drop, namely a soliton-like object that has thedirection of magnetization in each core opposite to itssurroundings [52, 53]. D. Nematic DBB and DDD solitons
Subsequently we study the properties of nonlinearstructures in the EP phase. The latter as per the phasediagram of Fig. 1(b) corresponds to c = − × − and0 < q < q th supporting both DBB and DDD stationarystates. These distinct nonlinear excitations illustratedrespectively in the insets of Fig. 13(a) and 14(a), ap-pear to be unmagnetized since M z = 0 in both caseswhile having a nontrivial polarization as q is varied.Moreover, the DBB entities are found to be significantlybroader around q = 0 when compared to the highly lo-calized DDD solitons (see top left insets in Fig. 13(a)and 14(a)). Interestingly, DBB solitons deform rapidly,i.e. soon after the transition point separating the EPto EA phases is crossed and for q ≈ − . q ≈ − . m F = ± m F = 0component has disappeared and only a Thomas-Fermi2 −0.1 0 0.1 0.200.20.40.60.81 R e ( ω ) q −0.1 0 0.1 0.200.010.020.030.040.05 I m ( ω ) q (a) (b) FIG. 16. BdG spectrum of DDD solitons existing in the EPphase but upon varying the quadratic Zeeman coefficient, q ,so as to enter both the EA and the PO phases. (a) Realpart, Re( ω ), of the corresponding eigenfrequencies ω . Thetrajectories of all three AMs present in this spectrum areindicated by light blue dots. (b) Imaginary part, Im( ω ), ofthe respective eigenfrequencies. Other parameters used areΩ = 0 . µ , ± = 2, c = − × − , and c = 1. type profile remains in the m F = ± q is increased so as to approachthe critical point q = q th that separates EP and the POphases, a rather sharp transitioning takes place for DDDsolitons when compared to the significantly smootherone exhibited by the DBB stationary states. This sharpversus smooth transition can be inferred by inspectingthe relevant slopes of the polarization for 0 < q < q th .Specifically, it is found that DDD solitons morph faster,i.e. for q = 0 . m F = 0 component. This observation is in agree-ment with the prediction from the ground state analysisthreshold value of the quadratic energy term, which inturn suggests that for q = q th , the PO phase shouldbe reached [top right inset in Fig. 14(a) and Fig. 14(b),(c)]. Contrary to this deformation, it is only around q = 1 that the polarization measured for DBB solitonsasymptotes to P = 1 [Fig. 13(a)]. The latter togetherwith the relevant negligible populations, n m F ( q ) [Fig. 13(b)], and widths, w m F ( q ) [Fig. 13(c)], of the bright mat-ter waves hosted in the symmetric m F = ± m F = ± m F = 0 component.Our BdG results reveal that both DBB and DDDsolitons are stable configurations within the EP phaseas can be inferred by the zero imaginary part shownin Figs. 15(b) and 16(b) respectively. Also stable are −0.1 0 0.100.050.10.150.2 q R e ( ω ) −0.1 0 0.100.0050.010.0150.02 q I m ( ω ) (b)(a) FIG. 17. BdG spectrum of the EA metastable states forvarying the quadratic Zeeman coefficient, q , in order to crossboth the EP and the PO phases. (a) Real part, Re( ω ), and(b) imaginary part, Im( ω ), of the eigenfrequencies ω . Acascade of bifurcations in (b) marks the destabilization ofthe metastable state leading to a change in the Krein signdesignated by the trajectories of six anomalous modes shownin (a) by light blue dots. The remaining parameters areΩ = 0 . µ , ± = 2, c = − × − and c = 1. the single dark solitons (into which the above DBBsand DDDs morph) in the PO phase, whose stabilityanalysis simply leads to the standard oscillatory mo-tion, with oscillation frequency ω osc = Ω / √ q decreases soas to enter the EA phase our stability analysis showsthat an eigenfrequency zero crossing occurs right at thetransition point ( q = 0) suggesting the destabilizationof the DBB wave. Below this point and specifically for q < − .
007 different types of stationary states exist forthe spin-1 system. These new metastable states con-sist of an unpopulated m F = 0 component and twonearly TF density profiles occurring in the other twoequally populated symmetric magnetic sublevels. Thefinite growth rate, Im( ω ), depicted in Fig. 15(b) un-veils the emergence of these new unstable configura-tions. For completeness, we examine the stability ofthese effective two-component metastable states, by ini-tializing them in our fixed point algorithm. The rele-vant BdG spectrum is shown in Fig.17(a)-(b). A cas-cade of bifurcations in addition to a persistent imag-inary eigenfrequency is observed for these metastablesymmetric states rendering them unstable for all valuesof q ∈ [ − . , .
1] that we have considered.Turning to the DDD solitons, for these negative3 t x −20020 x −20020 t x x −20 0 2000.51 x t n n +1 n − | Ψ ( x, t ) | | Ψ +1 ( x, t ) | | Ψ − ( x, t ) | t (l) | Ψ − ( x, t ) | | Ψ +1 ( x, t ) | | Ψ ( x, t ) | | Ψ +1 ( x, t ) + u AM | | Ψ ( x, t ) + u AM | | Ψ − ( x, t ) + u AM | | Ψ − ( x, t ) + u AM | | Ψ +1 ( x, t ) + u AM | (k)(j) q = − . q = 0 . (a)(b)(c) (e) (h)(i)(g)(f)(d) | Ψ ( x, t ) + u AM | t = 3030 t = 3030 q = 0 . FIG. 18. Spatio-temporal evolution | Ψ m F ( x, t ) | of (a)-(c) ametastable state in the EA phase and of a stationary DBBsoliton, | Ψ m F + u AM | , in (d)-(f) the EP and (g)-(i) the POphases upon perturbing it with the eigenvector associatedwith the AM appearing in the BdG spectrum of Fig. 15 (seelegends). (j)-(k) Profile snapshots of the densities of themetastable states shown in (a)-(c) at t = 3030 to illustratethe droplet formation and (l) the corresponding temporalevolution of the respective populations, n m F ( t ). In all cases m F = 0 , ± q = − . q = 0 .
05 and q = 0 . . µ , ± = 2, c = − × − and c = 1. quadratic Zeeman energies, we can easily deduce thatalso these waves gradually deform. Their destabiliza-tion as detected by the finite growth rate observed inFig. 16(b) occurs at q ≈ − .
009 rendering also theseDDD solitons unstable for q ∈ [ − . , − . ω ) (cid:54) = 0 even deeper in the EA phasethis further implies that also the DD solitons that areformed for q < − .
017 exist as unstable configurationsfor this value of q onward within the EA phase. In-teresting the instability of these states is caused by animaginary eigenfrequency reflecting the co-existence ofthese two components, analogously to Fig. 17(b).Confirmation of the above-obtained stability analysisresults is provided in Fig. 18(a)-(i) for the DBB solu-tions and in Fig. 19(a)-(i) and Fig. 20(a)-(f) for theDDD waves. Notice the coherent particle-like oscilla-tions observed for the stable DBB [Fig. 18(d)-(i)] andDDD [Fig. 19(a)-(i)] solitons for q > q <
0. Thespatio-temporal evolution of the metastable states de-picted in Fig. 18(a)-(c) is apparently rather similar tothe one found for the relevant states upon crossing theEA to EP phase boundary [see Fig. 12(a)-(c)]. Here,however, the two symmetric nonzero m F componentslose atoms towards the m F = 0 state in a nearly periodicfashion as a result of the instability. This dynamicalevolution leads to states featuring a flat-top, droplet-likeprofile [Fig. 18(j)] and nearly TF wavefunctions occupy- t t
100 200 t x x −505 x −505 00.5100.10.20.3 (a)(b) (h)(e)(d) (g) | Ψ − ( x, t ) + u AM | | Ψ − ( x, t ) + u AM | | Ψ − ( x, t ) + u AM | | Ψ +1 ( x, t ) + u AM | | Ψ +1 ( x, t ) + u AM | | Ψ +1 ( x, t ) + u AM | | Ψ ( x, t ) + u AM | | Ψ ( x, t ) + u AM | | Ψ ( x, t ) + u AM | (c) (f) (i) q = 0 . q = 0 . FIG. 19. (a)-(i) Spatio-temporal evolution of | Ψ m F ( x, t ) + u AM i | with m F = 0 , ± i = 1 , ,
3, for a stationary DDDsoliton within the EP phase but upon perturbing it with theeigenvector associated with each of the three AMs that ap-pear in the BdG spectrum of Fig. 16 (see legends). Columns(a)-(f), and (g)-(i), depict the evolution of the perturbedDDD state for q = 0 .
005 and q = 0 . . µ , ± = 2, c = − × − ,and c = 1. ing, respectively, the m F = 0 and m F = ± t > × (when this apparentperiodicity is modified).Next we monitor the relevant unstable evolution ofthe perturbed DD waves. Strikingly, their dynamicsentails completely new features as shown in Fig. 20(a)-(d). In this case, on top of the perturbed DD solitons,localized states having widths significantly larger thanthe healing length, which is the characteristic lengthscale of regular solitons, develop. These localized matterwaves consist of density humps followed by density dipsbuilding on top of the BEC background. They furtheremerge in an alternating fashion not only within butalso between the symmetric spin-components [see thedensity profiles at t = 2500 shown in Fig. 20(e)]. Thesestructures are reminiscent of phase separated states thathave been widely considered in multi-component con-densates [2, 5]. This is also reflected in the antisym-metric extended spatial profile of M z as can be seen inFig. 20(f). This quantity reflects the distinct spin do-mains formed among the m F = ± m F = ± x −20020 t x x −20 −10 0 10 20012 x | Ψ +1 ( x, u AM | | Ψ − ( x, u AM | t M z ( x, (e) (f) | Ψ − ( x, t ) + u AM | | Ψ +1 ( x, t ) + u AM | (d) | Ψ − ( x, t ) + u AM | (b) | Ψ +1 ( x, t ) + u AM | (a) (c) FIG. 20. (a)-(d) Spatio-temporal evolution of the densities, | Ψ m F ( x, t ) + u AM i | with m F = ± i = 1 ,
3, of an excitedDD soliton in the EA phase for q = − . m F = ± M z ( x ),at t = 2500 when the emergent spin-wave is spontaneouslynucleated (see legends). In all cases Ω = 0 . µ , ± = 2, c = − × − and c = 1. V. CONCLUSIONS
The phase diagram of solitonic nonlinear excitationsthat arise in 1D spin-1 harmonically trapped BECs hasbeen explored in detail. Particularly, we tackled theexistence of spinor matter waves in the form of DDD,DDB and DBB solitons and their stability in the ( c , q )-plane, motivated at least in part by recent experimentalrealizations of some of these states [22, 49]. Stationarystates of the DDB type are identified within the AF andthe EA phases, DBB solitons exist in the PO and in theEP phases and DDD solitons are found in the EP phase.Among the above solitonic entities DDB solitons ap-pear to be generally unstable configurations within theAF phase and to possess windows of stability withinthe EA phase. They were found to morph into widerstructures and into metastable states, respectively, asthe quadratic Zeeman energy increases towards the rele-vant for each phase transition boundary. In both phasesthe unstable dynamics of these configurations far fromthe associated transition point entails predominantlyparticle-like translational or breathing oscillations of theDDB waves. Such states feature remarkably long life-times which, in turn, is in line with their observabil-ity in existing experimental settings [49]. For the AFDDB solitons coherent spin-mixing dynamics occurs, asthe transition boundary separating the AF and the POphase is reached. For the ferromagnetic EA DDB waves,two distinct features take place. A stability interval oc-curs for quadratic Zeeman energies lying within the EAphase. Here, the solitary waves initially exist in con-figurations populating all the magnetic sublevels beforetheir deformation into symmetry broken DDBs and ulti- mately (deep within the EA regime) linearly stable DBsolitons. Secondly, in the vicinity of the transition point q → − , separating the EA and the EP phases, yet an-other deformation of the DDB configuration towards ametastable state appears. Remarkably, the unstable ir-regular spin-mixing dynamics of these latter states leadsto the formation of localized spin configurations consist-ing of Gaussian- and a droplet-like states occupying inan alternating fashion the distinct magnetic sublevels.Particularly here, the Gaussian serves as a matter wavebarrier pushing outwards, with respect to the trap cen-ter, the symmetric spin-components which develop be-tween them a DW (reminiscent of although differentfrom those of Ref. [64]). This DW bears a clear imprinton the local magnetization of the ensuing spin configura-tion which changes sign across it. Nearly periodic recur-rences of the aforementioned magnetic wave, that sharessimilarities to the so-called magnon drops [52, 53], areobserved while its persistence for extremely large prop-agation times suggests its potential observation.Turning to the DBB soliton solutions, it is found thatthese states penetrate the underlying first and secondorder phase transition boundaries. Specifically, unmag-netized DBB solitons enter the AF phase ( c >
0) up toa critical point below which they abruptly deform intoground state of the AF phase. The latter is character-ized by equally populated symmetric m F = ± c < m F = 0 spin-component takesplace upon increasing this time the quadratic coefficient.Interestingly, the DBB solitons exist as unstable con-figurations for antiferromagnetic interactions, predomi-nantly in the PO phase ( c >
0) but they are linearlystable within the EP phase which is characterized by c < m F = 0spin-component. This morphing, while analogous tothe one observed for the DBB entities, here it appearsfor smaller quadratic Zeeman energies. Additionally,it is found that these dark waves, being rather persis-tent configurations, further enter the EA phase. Thispenetration holds up to a critical point, in terms of de-creasing the quadratic Zeeman coefficient, below whichthe DDD waves gradually deform into two dark solitonsequally populating the symmetric spin sublevels. Theseunstable DD solitons dynamically evolve into spin-waveswhich nevertheless retain their shared dark soliton atthe center. The aforementioned spin-waves exhibit afinite spatial magnetization having a spatial extent pro-portional to the spin healing length.5Several extensions of the present work can be putforth. As a first step one can consider quenches betweenthe first and second order phase transition boundariesto investigate how the relevant in each phase spin-1 soli-ton solutions migrate or transform from one phase to theother. Yet another interesting perspective would be tostudy interactions [22, 43, 50] between the spinor soli-tons identified within each phase of the above-obtainedphase diagram or even unravel lattices consisting of mul-tiple DBB, DDB and DDD spinorial solitons in analogyto the two-component settings e.g. of Refs. [21, 25]. Insome cases where different solutions co-exist (e.g. theDBB and the DDD in the EP phase), one could evenconsider collisions between different types of entities.Another aspect that the present work motivates furtherconcerns the study of domain wall configurations in suit-able regimes of the relevant phase diagram (e.g., within the EA state). Furthermore, generalizing the phase dia-gram of nonlinear excitations extracted herein in higherdimensions where vortex-bright states [69, 70] are ex-pected to form would be also a fruitful future direction. VI. ACKNOWLEDGEMENTS
P.G.K. is grateful to the Leverhulme Trust and to theAlexander von Humboldt Foundation for support and tothe Mathematical Institute of the University of Oxfordfor its hospitality. G.C.K and P.S. gratefully acknowl-edge financial support by the Deutsche Forschungsge-meinschaft (DFG) in the framework of the SFB 925“Light induced dynamics and control of correlated quan-tum systems”. S.I.M. gratefully acknowledges financialsupport in the framework of the Lenz-Ising Award ofthe University of Hamburg. [1] C.J. Pethick and H. Smith,
Bose-Einstein condensa-tion in dilute gases , Cambridge University Press (Cam-bridge, 2002).[2] L.P. Pitaevskii and S. Stringari,
Bose-Einstein Conden-sation , Oxford University Press (Oxford, 2003).[3] I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. , 885 (2008).[4] F.Kh. Abdullaev, A. Gammal, A.M. Kamchatnov, L.Tomio, Int. J. Mod. Phys. B , 3415 (2005).[5] P. G. Kevrekidis, D. J. Frantzeskakis, and R. Carretero-Gonz´alez, The Defocusing Nonlinear Schr¨odinger Equa-tion , SIAM (Philadelphia, 2015).[6] S. Burger, K. Bongs, S. Dettmer, W. Ertmer, K. Seng-stock, A. Sanpera, G. V. Shlyapnikov, and M. Lewen-stein, Phys. Rev. Lett. , 5198 (1999).[7] J. Denschlag et al., Science , 97 (2000).[8] L. Khaykovich, F. Schreck, G. Ferrari, T. Bourdel, J.Cubizolles, L. D. Carr, Y. Castin, and C. Salomon, Sci-ence , 1290 (2002).[9] A. Weller, J. P. Ronzheimer, C. Gross, D. J.Frantzeskakis, G. Theocharis, P. G. Kevrekidis, J. Es-teve, and M. K. Oberthaler, Phys. Rev. Lett. ,130401 (2008).[10] D. J. Frantzeskakis, J. Phys. A: Mathematical and The-oretical, (21), 213001 (2010).[11] G. Lamporesi, S. Donadello, S. Serafini, F. Dolfovo, andG. Ferrari, Nat. Phys. , 656 (2013).[12] T. Busch, and J. R. Anglin, Phys. Rev. Lett. , 010401(2001).[13] C. Becker, S. Stellmer, P. Soltan-Panahi, S. D¨orscher,M. Baumert, E. M. Richter, J. Kronj¨ager, K. Bongs,and K. Sengstock, Nat. Phys. , 496 (2008).[14] C. Hamner, J. J. Chang, P. Engels, and M. A. Hoefer,Phys. Rev. Lett. , 065302 (2011).[15] D. Yan, J. J. Chang, C. Hamner, M. Hoefer, P. G.Kevrekidis, P. Engels, V. Achilleos, D. J. Frantzeskakis,and J. Cuevas, J. Phys. B: At. Mol. Opt. Phys. ,115301 (2012). [16] D. Garrett, T. Klotz, B. Prinari, and F. Vitale, Applic.Anal. , 379 (2013).[17] B. Prinari, F. Vitale, and G. Biondini, J. Math. Phys. , 071505 (2015).[18] J. Ieda, T. Miyakawa, and M. Wadati, Phys. Rev. Lett. , 194102 (2004).[19] J. Ieda, T. Miyakawa, and M. Wadati, J. Phys. Soc.Jpn. , 2996 (2004).[20] M. Uchiyama, J. Ieda, and M. Wadati, J. Phys. Soc.Jpn. , 064002 (2006).[21] A. Romero-Ros, G. C. Katsimiga, P. G. Kevrekidis, andP. Schmelcher, Phys. Rev. A , 013626 (2019).[22] S. Lannig, C. M. Schmied, M. Pr¨ufer, P. Kunkel, R.Strohmaier, H. Strobel, T. Gasenzer, P. G. Kevrekidis,M. K. Oberthaler, arXiv:2005.13278 (2020).[23] S. V. Manakov, , Sov. Phys. JETP , 248 (1974).[24] E. T. Karamatskos, J. Stockhofe, P. G. Kevrekidis, andP. Schmelcher, Phys. Rev. A , 043637 (2015).[25] D. Yan, F. Tsitoura, P. G. Kevrekidis, and D. J.Frantzeskakis, Phys. Rev. A , 023619 (2015).[26] G. C. Katsimiga, J. Stockhofe, P. G. Kevrekidis, and P.Schmelcher, Phys. Rev. A , 013621 (2017).[27] G. C. Katsimiga, P. G. Kevrekidis, B. Prinari, G. Bion-dini, and P. Schmelcher, Phys. Rev. A , 043623(2018).[28] H. Kiehn, S. I. Mistakidis, G. C. Katsimiga, P.Schmelcher Phys. Rev. A , 023613 (2019).[29] C. Qu, L. P. Pitaevskii, and S. Stringari, Phys. Rev.Lett. , 160402 (2016).[30] A. M. Kamchatnov, Y. V. Kartashov, P. -´E. Larr´e, andN. Pavloff Phys. Rev. A , 033618 (2014).[31] C. Qu, M. Tylutki, S. Stringari, and L. P. Pitaevskii,Phys. Rev. A , 033614 (2017).[32] I. Danaila, M. A. Khamehchi, V. Gokhroo, P. Engels,and P. G. Kevrekidis, Phys. Rev. A , 053617 (2016).[33] A. Farolfi, D. Trypogeorgos, C. Mordini, G. Lamporesi,and G. Ferrari, Phys. Rev. Lett. , 030401 (2020).[34] X. Chai, D. Lao, K. Fujimoto, R. Hamazaki, M. Ueda,and C. Raman, Phys. Rev. Lett. , 030402 (2020). [35] G. C. Katsimiga, S. I. Mistakidis, T. M. Bersano, M. K.H. Ome, S. M. Mossman, K. Mukherjee, P. Schmelcher,P. Engels, and P. G. Kevrekidis, arXiv:2003.00259 (2020).[36] Y. Kawaguchi, and M. Ueda, Phys. Rep. , 253(2012).[37] C.-M. Schmied, T. Gasenzer, M. K. Oberthaler, andP. G. Kevrekidis, Comm. Nonlin. Sc. Num. Sim. ,105050 (2020).[38] K. M. Mittal, S. I. Mistakidis, P. G. Kevrekidis, and P.Schmelcher, Phys. Rev. A , 013302 (2020).[39] L. Li, Z. Li, B. A. Malomed, D. Mihalache, and W. M.Liu, Phys. Rev. A , 033611 (2005).[40] B. J. Dabrowska-W¨uster, E. A. Ostrovskaya, T. J.Alexander, and Yu. S. Kivshar, Phys. Rev. A , 023617(2007).[41] W. Zhang, ¨O. E. M¨ustecaplıo˘glu, and L. You, Phys.Rev. A , 043601 (2007).[42] H. E. Nistazakis, D. J. Frantzeskakis, P. G. Kevrekidis,B. A. Malomed, and R. Carretero-Gonz´alez, Phys. Rev.A , 033612 (2008).[43] P. Szankowski, M. Trippenbach, and E. Infeld, Eur.Phys. J. D , 49 (2011).[44] D. M. Stamper-Kurn, and M. Ueda, Rev. Mod. Phys. , 1191 (2013).[45] H.-J. Miesner, D. M. Stamper-Kurn, J. Stenger, S. In-ouye, A. P. Chikkatur, and W. Ketterle, Phys. Rev.Lett. , 2228 (1999).[46] T. ´Swis(cid:32)locki, and M. Matuszewski, Phys. Rev. A, ,023601 (2012).[47] T. Ohmi, and K. Machida, J. Phys. Soc. Jpn , 1822(1998).[48] S. W. Song, L. Wen, C. F. Liu, S. C. Gou, and W. M.Liu, Frontiers of Physics , 302 (2013).[49] T. M. Bersano, V. Gokhroo, M. A. Khamehchi, J.D’Ambroise, D. J. Frantzeskakis, P. Engels, and P. G.Kevrekidis, Phys. Rev. Lett. , 063202 (2018).[50] L. Z. Meng, Y. H. Qin, and L. C. Zhao, arXiv:1912.00182 (2019).[51] K. Fujimoto, R. Hamazaki, and M. Ueda, Phys. Rev.Lett. , 173001 (2019).[52] F. Maci`a, D. Backes, and A. D. Kent, Nat. Nanotech. , 992 (2014).[53] B. Divinskiy, S. Urazhdin, V. E. Demidov, A. Kozhanov,A. P. Nosov, A. B. Rinkevich, and S. O. Demokritov,Phys. Rev. B , 224419 (2017).[54] J. Stenger, S. Inouye, D. Stamper-Kurn, H. Miesner, A.Chikkatur, and W. Ketterle, Nat. , 345 (1998).[55] J.-P. Martikainen, A. Collin, and K-A. Suominen, Phys.Rev. A , 053604 (2002).[56] Bo Xiong, and Jiangbin Gong, Phys. Rev. A , 033618(2010).[57] W. Zhang, S. Yi, and L. You, New J. Phys. , 77 (2003).[58] L. Santos, M. Fattori, J. Stuhler, and T. Pfau, Phys.Rev. A , 053606 (2007).[59] S. R. Leslie, J. Guzman, M. Vengalattore, J. D. Sau, M.L. Cohen, and D. M. Stamper-Kurn, Phys. Rev. A ,043631 (2009).[60] E. M. Bookjans, A. Vinit, and C. Raman, Phys. Rev.Lett. , 195306 (2011).[61] C. T. Kelley. Solving Nonlinear Equations with Newton’sMethod , Society for Industrial and Applied Mathemat-ics, Philadelphia, 1995. [62] L. Carr,
Understanding quantum phase transitions (CR-Cpress, 2010).[63] S. Sachdev, Physics World , 33 (1999).[64] H. E. Nistazakis, D. J. Frantzeskakis, P. G. Kevrekidis,B. A. Malomed, R. Carretero-Gonz´alez, and A. R.Bishop, Phys. Rev. A , 063603 (2007).[65] Note that in order to estimate the solitonic widths wemeasure the full-width-at-half-maximum (FWHM) ofthe ensuing matter waves for varying q .[66] P. G. Kevrekidis, and D. J. Frantzeskakis, Reviews inPhysics , 140 (2016).[67] Dmitry V. Skryabin, Phys. Rev. A , 013602 (2000).[68] T. Kapitula, and P. G. Kevrekidis, Chaos , 037114(2005).[69] K. J. H. Law, P. G. Kevrekidis, and L.S. TuckermanPhys. Rev. Lett. , 160405 (2010).[70] K. Mukherjee, S. I. Mistakidis, P. G. Kevrekidis, and P.Schmelcher, J. Phys. B: At. Mol. Opt. Phys.53