Phase Diagram, Stability and Magnetic Properties of Nonlinear Excitations in Spinor Bose-Einstein Condensates

We present the phase diagram, the underlying stability and magnetic properties as well as the dynamics of nonlinear solitary wave excitations arising in the distinct phases of a harmonically confined spinor $F=1$ Bose-Einstein condensate. Particularly, it is found that nonlinear excitations in the form of dark-dark-bright solitons exist in the antiferomagnetic and in the easy-axis phase of a spinor gas, being generally unstable in the former while possessing stability intervals in the latter phase. Dark-bright-bright solitons can be realized in the polar and the easy-plane phases as unstable 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 means of 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 others the coherent evolution of magnetic dark-bright, nematic dark-bright-bright and dark-dark-dark solitons. Moreover, for the unstable cases unmagnetized or magnetic droplet-like configurations and spin-waves consisting of regular and magnetic solitons are seen to dynamically emerge remaining thereafter robust while propagating for extremely large evolution times. Our investigations pave the wave for a systematic production and analysis involving spin transfer processes of such waveforms which have been recently realized in ultracold experiments.


I. INTRODUCTION
Ultracold atoms constitute ideal platforms for investigating the nonlinear behavior of quantum many-body systems due to their high degree of controllability and isolation from the environment [1][2][3]. A principal example has been the exploration of dark and bright solitons and their dynamical manifestations, as well as their multi-dimensional and multi-component extensions in Bose-Einstein condensates (BECs) [4,5]. Indeed, a variety and admixtures of these types of excitations are nowadays known to exist in scalar [6][7][8][9][10][11], pseudospinor [12][13][14][15][16][17] and spinor [18][19][20][21][22] BECs. Importantly, in recent years, many of these works have featured experimental realizations of such excitations. However, up to now the majority of both theoretical and experimental endeavors has been mainly focused on studying solitons in single and pseudo-spinor BEC systems and also within the so-called Manakov limit [23]. The latter assumes the intra-and the inter-species coupling to be on equal footing. As such the physics of nonlinear excitations outside this limit is less explored although there is an ongoing theoretical effort in this direction over the past few years [24][25][26][27][28].
Arguably, even less explored appears to be the connection between regular e.g. vector solitons and magnetic solitons or higher spin objects such as F = 1 spinors and spin-waves. Namely, nonlinear structures for which the magnetic interactions between the species are a crucial component. For instance, a magnetic soliton typically residing in a spin balanced density back-ground [29] is characterized by a localized spin magnetization, measured as the difference between the population of the participating components. Such nonlinear polarization waves have also been studied earlier [30] in binary BECs for parametric variations lying outside the Manakov limit both in the absence and in the presence of a Rabi coupling between the ensuing components [31]. Case examples of magnetic solitons are dark and dark-bright (DB) matter waves that differ from their regular or standard counterparts in a two-fold manner: (i) they exist for unequal intra and interspin couplings and (ii) their width scales according to the spin-healing length [29]. They can also have the form of dark-antidark solitons, with the latter being density bumps on top of the BEC background, that have been very recently experimentally monitored [32][33][34][35]. F = 1 spinor BECs offer the possibility for studying not only regular solitons but also magnetic ones and admixtures thereof. In particular, owing to the far richer phase diagram exhibited by such gases [36] (see, also, [37] for a recent discussion and [38] for the impact of many-body effects) already several works have been devoted in studying a variety of nonlinear excitations that arise in them [39][40][41][42][43][44]. These include for instance spin domains [45,46], spin textures [47,48], the very recently experimentally observed dark-dark-bright (DDB) and dark-bright-bright (DBB) solitons [49] (and variants [50], as well as interactions [22] thereof) and even twisted magnetic solitons [51].
In the present work we exploit this momentum and by following recent experiments [49] we first map out arXiv:2008.00475v1 [nlin.PS] 2 Aug 2020 the phase diagram of nonlinear excitations that arise in a one-dimensional (1D) harmonically confined spinor F = 1 BEC. More precisely, DDB, DBB and dark-darkdark (DDD) solitons constitute the principal ingredients of this phase diagram. DDB solutions exist in the antiferromagnetic (AF) and the easy-axis (EA) phase, DBB solitons 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 properties of the solitonic solutions involved in each phase and read out the magnetic properties of the participating solitons. Finally, we consider transitions between the distinct phases of the spinor system, by means of varying the quadratic Zeeman energy shift, and study deformations of the above-identified matter waves. Dynamical evolution of stable and unstable configurations reveals among others: the coherent evolution of magnetic DB solitons and changes in the magnetic properties of the evolved entities including the formation of composite spin objects. The latter are composed of regular solitons and spin-waves. Additionally we observe metastable states evolving into periodically recurring unmagnetized Thomas-Fermi (TF)-droplet configurations and also magnetized entities with droplets occupying the symmetric spin sublevels -with a domain wall (DW) separating them-and a localized wavefunction hosted in the remaining spin-component. The latter nearly periodic structures closely resemble magnon drops [52,53], while in both cases domain-walls (DW) are imprinted in the local magnetization.
Our work is structured as follows. In section II the relevant mean-field theoretical framework is introduced. The ground state phase diagram of a harmonically trapped 1D spin-1 BEC is initially discussed in Section III and we then proceed to the presentation and systematic exploration of the relevant phase diagram of nonlinear excitations in the form of DDD, DDB and DBB solitons. Section IV addresses the existence, the stability properties, by means of Bogoliubov de-Gennes (BdG) linearization analysis, and subsequently the dynamics of the different solitonic waveforms that arise in the distinct phases of the spinor system. Finally, in section V we summarize our findings and also provide future perspectives.

II. SPINOR SETUP AND MAGNETIZATION MEASURES
A spin-1 BEC composed of the magnetic sublevels m F = 0, ±1 of the hyperfine state F = 1, either of a 87 Rb [49] or a 23 Na [54] atom gas being confined in a 1D harmonic trap is considered. A cigar-shaped geometry is employed that has been very recently realized experimentally [49] utilizing a highly anisotropic trap with the longitudinal and transverse trapping frequencies obeying ω x ω ⊥ . In the mean-field framework the dynamics of such a spinor system can be described by the following coupled dimensionless Gross-Pitaevskii equations (GPEs) of motion [40,42,55,56] for the m F = 0 magnetic sublevel, while the symmetric m F = ±1 spin-components obey In Eqs. (1)-(2), Ψ m F (x, t) denotes the wavefunction of the |F = 1, m F = 0 and |F = 1, m F = ±1 spincomponents respectively. The single particle Hamilto- 2 Ω 2 x 2 being the 1D harmonic potential. Here, Ω ≡ ω x /ω ⊥ plays the role of the longitudinal over the transverse trapping frequency and is typically a small parameter i.e., Ω 1 [1,2]. Additionally, q denotes the quadratic Zeeman energy shift parameter that leads to an effective detuning of the m F = ±1 spin-components with respect to the m F = 0 one. It is quadratically proportional to an external magnetic field applied along the spin-z direction [36,57] and can be experimentally tuned by either adjusting the applied magnetic field [58] or by using a microwave dressing field [59,60].
Moreover, c 0 and c 1 are the so-called spinindependent and spin-dependent interaction coefficients. The former accounts for attractive (repulsive) interatomic interactions upon taking negative (positive) values and the latter is positive (c 1 > 0) for antiferromagnetic and negative (c 1 < 0) for ferromagnetic interactions. Both c 0 and c 1 are expressed in terms of the s-wave scattering lengths a 0 and a 2 , accounting for two atoms in the scattering channels with total spin F = 0 and F = 2 respectively, via the relations c 0 = (a0+2a2) 3a ⊥ and c 1 = (a2−a0) 3a ⊥ [21,40]. Here, a ⊥ = /M ω ⊥ is the transverse harmonic oscillator length with M denoting the mass, e.g., of a 87 Rb atom. Eqs. (1)-(2) have been made dimensionless by measuring length, energy and time in units of /(M ω ⊥ ), ω ⊥ and ω −1 ⊥ respectively. Consequently, the corresponding interaction strengths are expressed in terms of 3 ω ⊥ /M . In the adopted units and for a ferromagnetic, i.e., c 1 < 0, spinor BEC of 87 Rb atoms [21,37], the experimentally measured spin-dependent and spin-independent couplings also used herein are c 1 ≈ −5 × 10 −3 3 ω ⊥ /M and c 0 = 1 3 ω ⊥ /M . Additionally, the population of each spin-component is defined as Here, N = m F dx|Ψ m F | 2 denotes the total number of particles that is a conserved quantity for the spinorial system of Eqs. (1)- (2). Evidently, 0 ≤ n m F ≤ 1 is satisfied. Furthermore, in order to quantify first-and second-order transitions between the distinct phases of the spin-1 BEC system as well as to monitor the magnetic properties of the emergent nonlinear excitations during evolution we utilize the magnetization along the spin-z-axis that reads M z essentially measures the population imbalance between the symmetric m F = ±1 components and −1 ≤ M z ≤ 1. For instance, a fully magnetized state along the +z or −z spin direction corresponds to M z = +1 or M z = −1 respectively. To encounter also possible population transfer between the m F = 0 and the m F = ±1 spin states we invoke the polarization of the spinorial setting defined as follows [36] It can be easily deduced that −1 ≤ P ≤ 1. As we will unveil later on, P also accounts for alterations in the magnetic properties of the spinor system and allows us to distinguish among fully magnetized and unmagnetized spin configurations as we cross, by means of varying the quadratic Zeeman energy shift q, a phase transition boundary. Finally, for the numerical investigations that follow, the trapping frequency is fixed to Ω = 0.1, but we note that the results presented herein are not altered even for trapping frequencies of the order of Ω = 0.01 that are used in recent spin-1 BEC experiments [49]. Only slight deviations of the corresponding transition boundaries are observed. For instance, for the EP to PO transition, while q th ≡ 2n|c 1 | = 0.02 (with n denoting the peak density) for Ω = 0.01 it is q th ≈ 0.017 for Ω = 0.1. Additionally, c 0 = 1, c 1 = ±5×10 −3 and we choose the chemical potentials of the different components µ 0,±1 = 2. It is also important to mention that we have checked that the results to be presented below are robust also for c 0 = 1 and c 1 = 3.6 × 10 −2 , namely for the experimentally relevant interaction coefficient parameter ratio corresponding to 23 Na gas and also for larger chemical potentials, i.e. µ 0,±1 = 3 and µ 0,±1 = 5. To access the distinct phases of the spinor system, we typically vary q within the intervals [−1.5, 0.5] and [−0.5, 1.5]. Moreover, in order to identify the existence of stationary states a fixed-point numerical iteration scheme, based on Newton's method, is employed [61]. To simulate the dynamical evolution of the distinct DDD, DDB and DBB solitons governed by Eqs. (1)-(2), a fourth-order Runge-Kutta integrator is utilized while a second-order finite differences method is used for the spatial derivatives. The spatial and time discretization are dx = 0.05 and dt = 0.001 respectively. Our numerical computations are restricted to a finite region by employing hardwall boundary conditions. Particularly, in the dimensionless units adopted herein, the hard-walls are located at x ± = ±80 and we do not observe any appreciable density for |x| > 20.

III. PHASE DIAGRAM OF NONLINEAR EXCITATIONS
Before delving into the details of the phase diagram of nonlinear excitations in the form of DDD, DDB and DBB solitons that arise in spin-1 BECs, we first briefly revisit the relevant ground state phase diagram of a harmonically confined 1D spin-1 BEC [37]. This description will enable us to qualitatively expose the effect of embedding nonlinear structures into the different magnetic phases.

A. Ground state phase diagram
A schematic representation of the ground state phase diagram is illustrated in Fig. 1(a). As it has been recently demonstrated [36][37][38] different phases can be realized for such a confined spin-1 system. They stem from the interplay between the sign of the spin-dependent interaction coefficient c 1 and the strength of the quadratic Zeeman term q. Specifically, for c 1 > 0, q < 0 the system is in the AF phase with equally populated m F = ±1 spin-components thus having an unmagnetized ground state [see (Eq. 4)]. The latter is indeed characterized by M z = 0 and P = −1. A first order phase transition [62,63] separates this phase from the PO one that can be reached upon increasing q. The transition point appears at q = 0 and the PO phase is characterized again by an unmagnetized ground state but with all atoms populating the m F = 0 spin-component. Therefore M z = 0 and P = 1. On the other hand, for c 1 < 0, q < 0 the system resides in the EA phase. Its ground state is fully magnetized either along the +z or the −z spin-direction, i.e. either the m F = +1 or m F = −1 spin state is populated. As a result M z = +1 or M z = −1 respectively and P = −1. Upon increasing q a second-order phase transition occurs at q = 0 and for 0 < q < q th ≈ 0.017 the system enters the EP phase with its ground state having all three m F components populated. Particularly here, M z = 0 reflecting the fact that the m F = ±1 spin states are equally populated while P ∈ (−1, 1). Finally, for q > q th ≈ 0.017 yet another second order phase transition takes place which leads 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. Note here, that in order to obtain the above-discussed ground state phase diagram TF profiles are employed as initial guesses, within our fixed point algortithm [61], for the distinct m F = 0, ±1 states having the form When this expression is used here and below, it is implied to be valid when the quantity under the radical is non-negative (and the relevant wavefunction is padded with zeros outside that region). In Eq. (6), µ m F denotes the chemical potential of each spin-component while a stationary state satisfies the phase matching condition µ 0 = (µ +1 + µ −1 ) /2 [37,64].

B. Phase diagram of solitonic excitations
In order to unravel the phase diagram of nonlinear excitations depicted in Fig. 1(b), the spin-1 system is initialized in each of the above-identified phases embedding dark and bright solitons as wavefunctions for each of the m F = 0, ±1 states. In particular, the standard, stationary solitonic waveforms used read [49] In the above expressions Ψ D (x) and Ψ B (x) denote the wavefunctions utilized for a dark and a bright soliton configuration respectively. In Eq. (7) the quantity under the square root denotes the customary used TF background needed for dark solitons to be embedded on. Moreover, D and η refer, respectively, to the common inverse width considered for each spinorial soliton component and the amplitude of the bright soliton configuration (see our detailed discussion in Sec. IV). It is found that DDB solitons, being unmagnetized configurations, exist within the AF phase for all values of the quadratic Zeeman energy shift lying within the interval q ∈ (−1.5, 0), with the dark solitons effectively trapping the bright one appearing in the m F = 0 spincomponent. This trapping mechanism becomes progressively less effective. Namely, as q increases towards the phase transition point (q = 0) the bright soliton gradually becomes the dominant configuration before morphing into a TF one. We remark that the existence of a DDB soliton in the AF phase already presents fundamental deviations from its ground state properties. Indeed, in the latter case the m F = 0 magnetic sublevel is unpopulated (of course also the m F = ±1 states do not feature a dark soliton in the relevant ground state). On the contrary, DBB solitons, being again unmagnetized configurations, are identified in the PO phase, namely for q ∈ [0, 0.994) which deform towards a single dark soliton occupying the m F = 0 component for q > 0.994. Remarkably these states persist even upon decreasing q so as to enter the AF phase until a critical value of the quadratic Zeeman energy shift, i.e. q 3 ≈ −0.005, is reached [see dashed green box in Fig. 1(b)]. Note that such DBB configurations also constitute excited states within the PO phase since for the ground state only the m F = 0 state is occupied. Turning to c 1 < 0, stationary solutions of the DDB type are realized within the ferromagnetic EA phase existing within the parametric  region q ∈ (−0.515, 0.007). These DDB states deform as q decreases further into fully magnetized, i.e. M z = +1 (M z = −1), DB solitons that occupy the m F = 0 and m F = +1 (m F = 0 and m F = −1) components. Once again this is far from the ground state of the EA featuring only m F = +1 (or m F = −1) populations. Moving to q > 0, namely entering the EP phase, two types of solitonic solutions are found to exist for the spin-1 system. These excitations can have the form of unmagnetized spinor DBB or DDD solitons, a result that is permitted by the relevant ground state where all three magnetic components are occupied. The former solitonic entities appear to be significantly broader when compared to the more localized DDD configurations and become highly localized as we enter the PO phase. Recall that the PO ground state supports population only in the m F = 0 magnetic sublevel. The transition point for the DBB configuration appears at q 4 ≈ 0.994 while it occurs significantly earlier, q 4 ≈ 0.016, for the DDD state. Decreasing the quadratic Zeeman term, q, in order to enter the EA phase reveals that the DBB configuration deforms fast, around q ≈ −0.007, to a metastable state with two TF wavefunctions occupying the m F = ±1 components. Contrary to this deformation, DDD solitons continue to exist within the EA phase for values up to q 2 ≈ −0.016 before their transitioning towards two darks that occupy the symmetric m F = ±1 components [see dashed purple box in Fig. 1(b)].

IV. STABILITY ANALYSIS AND DYNAMICS OF SPINOR SOLITONS
Our aim in what follows is not only to illustrate the existence of stationary spinor solitons of the DDD, DDB and DBB type existing in a 1D harmonically confined spin-1 BEC composed e.g. of 87 Rb atoms and obeying Eqs. (1)- (2), but also to systematically investigate their stability properties. We remark that for ferromagnetic (AF) BECs we consider c 1 = −5 × 10 −3 (c 1 = 5 × 10 −3 ) as representative example and vary the quadratic Zeeman energy shift to access the underlying magnetic phases.

A. Antiferromagnetic DDB matter waves
For instance, in order to infer about the existence of DDB solitons within the AF phase shown in the phase diagram of Fig. 1(b), the matter wave dark solitons of Eq. (7) are embedded as initial guesses for the m F = ±1 spin-components and the bright soliton of Eq. (8) is utilized for the m F = 0 spin state. Employing the above ansatz, and using the iterative scheme discussed above, DDB stationary states are found within the AF phase, i.e. for c 1 = 5 × 10 −3 and for values of q ∈ (−1.5, 0). Characteristic DDB density profiles, |Ψ 0,±1 | 2 , are presented as insets in Fig. 2(a). However, upon increasing q towards the transition point (q = 0) above which the PO phase is realized, the DDB solitons deform into states where the bright structure in the m F = 0 component overfills/dominates the dark wells. Also the total den- sity, |Ψ tot | 2 , of the spinor system exhibits a TF profile instead of the dark-shaped density appearing deep in the AF phase. This altered nature of the DDB configuration, which remains unmagnetized (M z = 0) for all values of q, is naturally accompanied by a change in the polarization of this configuration. The DDB solitons possess P = −1 for q ≤ −1.5 reflecting the fact that deep in the AF phase only the m F = ±1 components bearing dark solitons are populated [ Fig. 2(b)], while the polarization takes values −1 < P ≤ 0 as we approach the transition point. At q ≈ −0.02 all three components are equally populated having significantly wider [65] stationary states [ Fig. 2(c)] as compared to the ones for larger negative q values. This broadening suggests that the DDB character of the relevant states is lost. Importantly, at q = 0 an abrupt population transfer to the m F = 0 component [ Fig. 2(b)] associated with the drastic deformation of this latter configuration to the ground state of the PO phase manifests itself; see e.g. the right uppermost inset of Fig. 2(a).
In order to extract the stability properties of the aforementioned DDB stationary states (as well as for the DBB and DDD solitons to be presented below), a linear stability or BdG analysis is performed. The latter consists of perturbing the iteratively identified in each phase stationary solutions Ψ 0 By inserting this ansatz into the system of Eqs. (1)- (2) and linearizing with respect to the small amplitude pa- rameter leads to an eigenvalue problem for the eigenfrequencies ω, or equivalently eigenvalues λ ≡ iω, and For further details on the BdG analysis we refer the reader to Refs. [5,66,67]. Due to the generally complex nature of the ensuing eigenfrequencies, it becomes apparent that the following possibilities can arise: if modes with purely real eigenvalues or equivalently imaginary eigenfrequencies or complex eigenvalues/eigenfrequencies are identified, these are responsible for the existence of an instability [66]. Moreover, due to the Hamiltonian structure of the system investigated herein, quartets of such eigenfrequencies can occur [67]. Namely if ω is an eigenfrequency so are −ω and ±ω * . As such, if Im(ω) = 0, then there will always exist a mode leading to the growth and eventual deformation of the examined in each phase solitonic configuration. The BdG analysis outcome for the DDB soliton solutions is shown in Fig. 3(a), (b). DDB solitons constitute excited states of the spin-1 system, exactly like their two-component dark-bright analogue [26], a feature that is reflected in their linearization spectra via the emergence of the so-called anomalous modes (AMs). These eigenstates are quantified via the negative energy or negative Krein signature [67] defined for the spinor system as The existence of these modes is central to our stabil- FIG. 8. Density evolution of (a)-(c) |Ψm F (x, t)| 2 referring to a deformed DBB soliton in the AF phase and (d)-(i) |Ψm F + uAM 1 | 2 of a DBB soliton within the PO phase but being excited by the eigenvector associated with the AM appearing in the BdG spectrum of Fig. 6. The different columns depict the evolution of the DBB state for q = −0.002, q = 0.024 and q = 0.5 respectively. Panels (j)-(l) show the temporal evolution of the populations, nm F (t), for the aforementioned values of q. In all cases Ω = 0.1, µ0,±1 = 2, c1 = 5 × 10 −3 and c0 = 1.
ity analysis since their potential collision with positive Krein signature modes can give rise to stabilitychanging events in the form of oscillatory instabilities or Hamiltonian-Hopf bifurcations [67]. In the present setting the AMs are denoted by light blue dots. The DDB solution possesses, due to the presence of two dark solitons [68], two such modes [ Fig. 3(a)] that cross the origin of the spectral plane at q cr = 0 signalling the destabilization of the DDB wave [ Fig. 3(b)]. Interestingly, and also for all values of q ∈ (−1.5, 0), it is found that the eigenvector associated with the lowest lying AM causes an overall shift when added to the stationary DDB solution. This in turn implies that a perturbed, with this eigenvector, DDB soliton will perform an oscillatory motion within the parabolic trap. On the contrary, the eigenvector corresponding to the higher lying AM, besides a weak displacement, further leads to an asymmetric DDB configuration. This asymmetry, as we shall show later on, is responsible for the breathing motion of the DDB entity and its effect is dominant with respect to the aforementioned shift. It is this higher lying AM that is responsible for the generic instability shown in Fig. 3(b).
However, for values of q closer to the critical point, defining the AF to PO transition boundary, the destabilization of both modes leads, irrespectively of which mode we excite (namely the first lower-lying one or the second), to a breathing motion of the DDB configuration. Its response is visualized in the spatio-temporal evolution of the densities, |Ψ m F (x, t) + u AM1 | 2 , presented in Fig. 4(a)-(c) entailing both the particle-like oscillation of the DDB soliton but predominantly the overall breathing of the state. Although the first AM is excited in this case, the dynamics is dominated by the breathing mode. The nature of this composite motion is also reflected in the irregular oscillation of the population, n m F (t), of each m F component illustrated in Fig. 4(g). At the transition/destabilization point the prevailing feature of the perturbed DDB configuration is its breathing as can be seen by monitoring the evolution of the densities illustrated in Fig. 4(d)-(f) together with the coherent oscillation of the relevant populations [ Fig. 4(h)]. In both of the aforementioned cases the oscillatory character of n m F (t) implies a weak amplitude spin-mixing dynamics.

B. Polar DBB solitons
Next we turn to the PO phase which is characterized by c 1 = 5 × 10 −3 and q > 0 [see also Fig. 1(a)]. According to the phase diagram of Fig. 1(b), here one can identify stationary DBB soliton solutions for values of q ∈ [0, 0.994). Specifically and so as to capture the occurrence of these solitonic waveforms we utilize, within our fixed point iteration scheme, the dark soliton ansatz of Eq. (7) as an initial guess for the m F = 0 spincomponent, while bright solitons given by Eq. (8) are considered for the remaining symmetric m F = ±1 components. These states are characterized by zero magnetization, preserving this way the magnetic properties of the ground state within this phase, but they have a polarization that acquires values −1 < P < 1 [ Fig. 5(a)]. In particular, P = 1 for q > 1, i.e. deep in the PO phase, and it gradually decreases as q → 0 + all the way to P = −1 for q < 0. This latter behavior of P reveals in turn that despite the fact that the ground state configuration does not support all three m F components to be populated this is not the case for the respective nonlinear excitations [ Fig. 5(b)]. Selected DBB soliton profiles are depicted as insets in Fig. 5(a). From these profiles it can be deduced that these unmagnetized DBB waves exist not only within the above-provided q interval but also at (q = 0) and below (q ∈ [0, −0.005)) the transition point that separates the PO and the AF phases. However, as we approach the transition point from above q → 0 + the DBB states deform towards wider configurations [ Fig. 5(c)] featuring a pronounced bright soliton component that dominates. This dominant bright component results in turn to a |Ψ tot | 2 that has a TF profile instead of a tanh-shaped one occurring for values of q well inside the PO phase. For q < −0.005 an abrupt transition leads to a metastable configuration in which the m F = ±1 are equally populated having also minimal polarization (P = −1), see the bottom left inset of Fig. 5(a). On the contrary, for q > 0.994 yet another but gradual this time deformation of the DBB matter waves towards a dark soliton with maximal polarization (P = +1) occupying the m F = 0 spin state occurs [top right inset of Fig. 5(a)].
By investigating the stability properties of the above solitonic entities, it is found that two destabilization points exist for the DBB configuration one residing in the AF phase and one deep in the PO phase. Specif-  ically for q < 0 the single in this case negative energy mode appearing in the BdG of Fig. 6(a) decreases in frequency and crosses the spectral plane at q cr ≈ −0.004 rendering these entities unstable for this value of q and thereafter. A result that is further supported by the finite growth rate, Im(ω) = 0, shown for these negative quadratic Zeeman energies in Fig. 6(b). This destabilization is related to a composite motion of the DBB structure in the parabolic trap that we will soon trace in the dynamics. Contrary to the above destabilization yet another critical point occurs for the DBB solution for positive values of q. The latter appears at q cr ≈ 0.994, i.e., the end point of the loop bifurcation illustrated in Fig. 6(b) above which DBB solitons cease to exist giving their place to a single dark solitary wave occupying the zeroth spin sublevel. This observation along with the second destabilization of the AM in this PO regime suggest the presence of a pitchfork bifurcation. In order to infer about the existence of the latter we perform the corresponding stability analysis of the PO dark states illustrated in the top right inset of Fig. 5(a). The relevant BdG spectrum is depicted in Fig. 7(a), (b). Interestingly enough, even though dark solitons exist for all quadratic Zeeman energies in q ∈ [0, 1.5], a narrow instability interval takes place for these stationary states occurring for q ∈ [0.994, 1] and leading to the loop bifurcation shown in Fig. 7(b). Within this q interval also the Krein signature changes sign from negative before the lower bound to positive after the upper bound. This, in turn, means that the dark soliton destabilizes slightly below unity and restabilizes for q > 1. It is in this interval that indeed the above identified DBB solitons coexist with the single dark ones in a subcrit-  ical pitchfork bifurcation. Namely, dark solitons exist as stable configurations, for q < 0.994, while their DBB counterparts are unstable. The collision of the two (and associated disappearance of the bright component of the DBB's) destabilizes the darks for q ∈ [0.994, 1], while for q > 1, the relevant real eigenvalue pair returns to the imaginary axis, restabilizing the relevant dark state. Direct evolution of the above-identified configurations slightly below and above the phase transition threshold at q = 0 reveals that the DBB solitons undergo in both cases an overall breathing motion [ Fig. 8(a) Fig. 8(j)]. In sharp contrast to the above dynamics for well-defined DBB solitons, i.e. away from the transition and the critical point, a well-defined, intrap oscillation of the perturbed DBB wave is observed [ Fig. 8(g)-(h)] for evolution times up to t = 5 × 10 3 with the respective populations, n m F (t), remaining constant for all times.

C. Easy-Axis symmetry broken DDB solitons
Moving on to the EA phase, i.e., for c 1 = −5 × 10 −3 and q < 0, again DDB stationary states are successfully identified. However, contrary to the DDB solutions found in the AF phase here the DDB waves exhibit unequally populated m F = ±1 spin-components as can be seen in the insets of Fig. 9(a) and also in the relevant populations of Fig. 9(b). Specifically, for these states the dark soliton of e.g. the m F = −1 spin state is suppressed for most of the q values within the region of existence, i.e. q ∈ (−0.515, 0.007), of this configuration. Notably, such waves preserve the symmetry (i.e., equal population) of the m F = ±1 components for q ∈ [−0.009, 0.007), namely including also the transition point (q = 0) that separates the EA and the EP phases. For the remaining quadratic Zeeman energies lying in the aforementioned q interval the symmetry is partially preserved, i.e. the m F = −1 is still populated. This result is encoded in the magnetization and the polarization properties of the DDB solutions which assume values 0 < M z < 1 and −1 < P < 1 respectively reflecting the non-negligible population of all three spin-components. More precisely, starting with P = 1 (M z = 0) for q > 0.007, the relevant quantity decreases (increases) when moving towards q < 0 and approaches the value of P = −1 (M z = +1) for q < −1, i.e. deep in the EA phase [ Fig. 9(a)]. These symmetric DDB solitons are fundamentally different (structurally) than the relevant ground state configuration in this parametric regime. The latter, according to the phase diagram of Fig. 1(a), favors symmetry broken states that are fully magnetized along the +z-or −z-spin direction, i.e. configurations that have either the m F = +1 or the m F = −1 component solely populated. As such, for q ∈ (−0.515, −0.2] the dark soliton of the m F = −1 spin state becomes narrower, in an almost exponentially decaying manner, as can be seen from the behavior of its width, w −1 (q), shown in Fig. 9(c). In particular, around q ≈ −0.515 the configuration is deformed to a symmetry broken almost fully magnetized (M z ≈ +1) DB soliton that exists for q ∈ [−1.5, −0.515] occupying the m F = +1 and m F = 0 spin components (similarly, of course, there is a state occupying the the m F = −1 and m F = 0 states). Additionally, as the transition point is approached from below, q → 0 − , also the pop- ulation, n 0 (q), and width, w 0 (q), of the bright component increases and the DDB solitons deform even further. Specifically, for q ∈ (−0.004, 0.007) a structure with a bright component that overfills the dark wells while gradually morphing into a TF profile can be identified, as shown in the upper left inset of Fig. 9(a). This deformed DDB structure enters the EP phase, which favors all three m F components to be simultaneously occupied, but already at q ≈ 0.008 the unmagnetized ground state of the PO phase is reached.
The BdG analysis of the above-discussed soliton solutions illustrated in Fig. 10(a), (b) reveals that DDB solitons possess potentially unstable eigendirections (although they also possess stability intervals). This result can be inferred by the finite imaginary eigenfrequencies (or instability growth rates), Im(ω), shown in Fig. 10 Notice that the last bifurcation possesses also the larger instability growth rate that stems from an eigenfrequency zero crossing of both the higherand the lower-lying AMs appearing at q cr = 0. For −0.515 < q < −0.384 the state remains linearly stable having, however, a minuscule m F = −1 component. For more negative values of q, the fully magnetized linearly stable DB solitons are present in the spin-1 system. To confirm the above stability analysis findings we have monitored the dynamical evolution of the DDB solutions in all of the above-identified instability intervals and our results can be summarized as follows. Among the two modes that appear in the BdG spectrum of Fig. 10(a) the lower one is related to the weak amplitude in trap oscillation of the DDB wave. The higher mode is responsible for the larger in amplitude anti-phase oscillation of the involved dark solitons. Additionally, it turns out that even when these states are found to be dynamically unstable they have remarkably long lifetimes that suggest their potential experimental observation in existing spinor settings [49]. A case example showcasing the particle-like oscillations that a perturbed DDB stationary state undergoes is presented for q = −0.002 in Fig. 11(a)-(c) and Fig. 11(d)-(f) respectively. Notice that indeed the DDB wave remains intact for all times up to t = 5 × 10 3 . More specifically, by perturbing the DDB soliton with the eigenvector associated with the first AM leads to an oscillation of the wave within the trap [ Fig. 11(a)-(c)]. On the other hand, the second mode results in the formation of two atomic blobs to which the dark solitons split the entire condensate [ Fig. 11(d)-(f)]. These blobs execute an anti-phase oscillation alternating across the two dark components and across the two sides (left and right) of the dark solitary wave in each component [ Fig. 11(d)-(f)]. This latter anti-phase oscillation becomes even more pronounced especially for the m F = −1 component as q decreases further towards the formation of DB solitons that occupy the m F = 0 and m F = +1 magnetic sublevels [ Fig. 11(g)-(i)]. Evidently, as q decreases further and e.g. for q = −0.08 illustrated in Fig. 11(g)-(i), the population of the m F = −1 magnetic sublevel becomes significantly suppressed.
Contrary to the above-described dynamics, the picture is drastically altered when considering the deformed DDB stationary states that exist near the transition point (q = 0). For instance here, by monitoring |Ψ m F (x, t)| 2 for quadratic Zeeman energy shifts that lie within the last bifurcation [ Fig. 10(b)] reveals that these transient states for evolution times of the order of t ≈ 4500 destabilize towards states that consist of Gaussian-like (localized) structures hosted in the m F = 0 component. These localized density blobs are not of permanent character as is evident in Fig. 12(a) but they revive in an almost periodic manner. In every recurrence event, the corresponding symmetric spincomponents bear droplet-like configurations that appear in an alternating fashion either in the m F = +1 or in the m F = −1 component [ Fig. 12(b), (c)]. This behavior essentially reflects the continuous spin transfer between the m F = 0 and m F = ±1 taking place during evolution [ Fig. 12(d)]. Inspecting the density profiles of the evolved states [ Fig. 12(e)] unveils that the wavefunction of the zeroth magnetic sublevel acts as a repulsive barrier pushing outwards, with respect to the trap center, the symmetric spin-components that develop in between them a DW [64]. Measuring the local magnetization, M z (x), e.g. at t 1 = 3690 where this dynamically formed state emerges for the first time [ Fig. 12(f) reveals that such a configuration bears indeed a DW character across which M z (x) changes sign [64]. Such a magnetic entity holds close similarities to the so-called magnon drop, namely a soliton-like object that has the direction of magnetization in each core opposite to its surroundings [52,53].

D. Nematic DBB and DDD solitons
Subsequently we study the properties of nonlinear structures in the EP phase. The latter as per the phase diagram of Fig. 1(b) corresponds to c 1 = −5 × 10 −3 and 0 < q < q th supporting both DBB and DDD stationary states. These distinct nonlinear excitations illustrated respectively in the insets of Fig. 13(a) and 14(a), appear to be unmagnetized since M z = 0 in both cases while having a nontrivial polarization as q is varied. Moreover, the DBB entities are found to be significantly broader around q = 0 when compared to the highly localized 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 EP to EA phases is crossed and for q ≈ −0.007, into the metastable state of the EA phase that has equally populated symmetric components ( Fig. 13(b)) when compared to the slower, around q ≈ −0.016, deformation of the DDD solitons into two dark ones equally populating the m F = ±1 spin states (Fig. 14(b)). Notice that in the former DDB case, the dark soliton in the m F = 0 component has disappeared and only a Thomas-Fermi type profile remains in the m F = ±1 components. Importantly though, as q is increased so as to approach the critical point q = q th that separates EP and the PO phases, a rather sharp transitioning takes place for DDD solitons when compared to the significantly smoother one exhibited by the DBB stationary states. This sharp versus smooth transition can be inferred by inspecting the relevant slopes of the polarization for 0 < q < q th . Specifically, it is found that DDD solitons morph faster, i.e. for q = 0.016, into a single dark state occupying the m F = 0 component. This observation is in agreement with the prediction from the ground state analysis threshold value of the quadratic energy term, which in turn suggests that for q = q th , the PO phase should be reached [top right inset in Fig. 14 the single dark solitons (into which the above DBBs and DDDs morph) in the PO phase, whose stability analysis simply leads to the standard oscillatory motion, with oscillation frequency ω osc = Ω/ √ 2 in the TF regime, known for harmonically trapped dark solitons [10]. Furthermore the existence of a single and three AMs pertaining to the DBB and the DDD configuration respectively can also be seen in the relevant real part of the spectrum illustrated in Figs. 15(a) and 16(a). Once again, it appears that the number of components bearing a dark solitary wave determines the number of AMs within the state of interest. However, and as far as the DBB solutions are concerned, as q decreases so as to enter the EA phase our stability analysis shows that an eigenfrequency zero crossing occurs right at the transition point (q = 0) suggesting the destabilization of the DBB wave. Below this point and specifically for q < −0.007 different types of stationary states exist for the spin-1 system. These new metastable states consist of an unpopulated m F = 0 component and two nearly TF density profiles occurring in the other two equally populated symmetric magnetic sublevels. The finite growth rate, Im(ω), depicted in Fig. 15(b) unveils the emergence of these new unstable configurations. For completeness, we examine the stability of these effective two-component metastable states, by initializing them in our fixed point algorithm. The relevant BdG spectrum is shown in Fig.17(a)-(b). A cascade of bifurcations in addition to a persistent imaginary eigenfrequency is observed for these metastable symmetric states rendering them unstable for all values of q ∈ [−0.1, 0.1] that we have considered.
Turning to the DDD solitons, for these negative  quadratic Zeeman energies, we can easily deduce that also these waves gradually deform. Their destabilization as detected by the finite growth rate observed in Fig. 16(b) occurs at q ≈ −0.009 rendering also these DDD solitons unstable for q ∈ [−0.016, −0.009]. However, since Im(ω) = 0 even deeper in the EA phase this further implies that also the DD solitons that are formed for q < −0.017 exist as unstable configurations for this value of q onward within the EA phase. Interesting the instability of these states is caused by an imaginary eigenfrequency reflecting the co-existence of these two components, analogously to Fig. 17(b). Confirmation of the above-obtained stability analysis results is provided in Fig. 18(a)-(i) for the DBB solutions and in Fig. 19(a)-(i) and Fig. 20(a)-(f) for the DDD waves. Notice the coherent particle-like oscillations observed for the stable DBB [ Fig. 18(d)-(i)] and DDD [ Fig. 19(a)-(i)] solitons for q > 0 when compared to the unstable evolution of the densities for q < 0. The spatio-temporal evolution of the metastable states depicted in Fig. 18(a)-(c) is apparently rather similar to the one found for the relevant states upon crossing the EA to EP phase boundary [see Fig. 12 19. (a)-(i) Spatio-temporal evolution of |Ψm F (x, t) + uAM i | 2 with mF = 0, ±1, i = 1, 2, 3, for a stationary DDD soliton within the EP phase but upon perturbing it with the eigenvector associated with each of the three AMs that appear in the BdG spectrum of Fig. 16 (see legends). Columns (a)-(f), and (g)-(i), depict the evolution of the perturbed DDD state for q = 0.005 and q = 0.1 respectively. Other parameters used are Ω = 0.1, µ0,±1 = 2, c1 = −5 × 10 −3 , and c0 = 1.
ing, respectively, the m F = 0 and m F = ±1 components [ Fig. 18(k)]. Coherent population transfer accompanies the periodic revival of these states [ Fig. 18(l)] which remain nematic, i.e. having zero magnetization, during evolution. This way, they preserve the magnetic properties expected for an EP configuration. Furthermore, these unmagnetized structures appear to robustly reemerge up to times t > 8 × 10 3 (when this apparent periodicity is modified).
Next we monitor the relevant unstable evolution of the perturbed DD waves. Strikingly, their dynamics entails 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 than the healing length, which is the characteristic length scale of regular solitons, develop. These localized matter waves consist of density humps followed by density dips building on top of the BEC background. They further emerge in an alternating fashion not only within but also between the symmetric spin-components [see the density profiles at t = 2500 shown in Fig. 20(e)]. These structures are reminiscent of phase separated states that have been widely considered in multi-component condensates [2,5]. This is also reflected in the antisymmetric extended spatial profile of M z as can be seen in Fig. 20(f). This quantity reflects the distinct spin domains formed among the m F = ±1 components, while the dark soliton remains in the middle being shared by the two otherwise phase-separated m F = ±1 components. (e) (f)

V. CONCLUSIONS
The phase diagram of solitonic nonlinear excitations that arise in 1D spin-1 harmonically trapped BECs has been explored in detail. Particularly, we tackled the existence of spinor matter waves in the form of DDD, DDB and DBB solitons and their stability in the (c 1 , q)plane, motivated at least in part by recent experimental realizations of some of these states [22,49]. Stationary states of the DDB type are identified within the AF and the EA phases, DBB solitons exist in the PO and in the EP phases and DDD solitons are found in the EP phase.
Among the above solitonic entities DDB solitons appear to be generally unstable configurations within the AF phase and to possess windows of stability within the EA phase. They were found to morph into wider structures and into metastable states, respectively, as the quadratic Zeeman energy increases towards the relevant for each phase transition boundary. In both phases the unstable dynamics of these configurations far from the associated transition point entails predominantly particle-like translational or breathing oscillations of the DDB waves. Such states feature remarkably long lifetimes which, in turn, is in line with their observability in existing experimental settings [49]. For the AF DDB solitons coherent spin-mixing dynamics occurs, as the transition boundary separating the AF and the PO phase is reached. For the ferromagnetic EA DDB waves, two distinct features take place. A stability interval occurs for quadratic Zeeman energies lying within the EA phase. Here, the solitary waves initially exist in configurations populating all the magnetic sublevels before their deformation into symmetry broken DDBs and ulti-mately (deep within the EA regime) linearly stable DB solitons. Secondly, in the vicinity of the transition point q → 0 − , separating the EA and the EP phases, yet another deformation of the DDB configuration towards a metastable state appears. Remarkably, the unstable irregular spin-mixing dynamics of these latter states leads to the formation of localized spin configurations consisting of Gaussian-and a droplet-like states occupying in an alternating fashion the distinct magnetic sublevels. Particularly here, the Gaussian serves as a matter wave barrier pushing outwards, with respect to the trap center, the symmetric spin-components which develop between them a DW (reminiscent of although different from those of Ref. [64]). This DW bears a clear imprint on the local magnetization of the ensuing spin configuration which changes sign across it. Nearly periodic recurrences of the aforementioned magnetic wave, that shares similarities to the so-called magnon drops [52,53], are observed while its persistence for extremely large propagation times suggests its potential observation.
Turning to the DBB soliton solutions, it is found that these states penetrate the underlying first and second order phase transition boundaries. Specifically, unmagnetized DBB solitons enter the AF phase (c 1 > 0) up to a critical point below which they abruptly deform into ground state of the AF phase. The latter is characterized by equally populated symmetric m F = ±1 spincomponents. Moreover, metastable states are identified also for the nematic DBB waves (c 1 < 0), originating from the deformation of these solitons as the quadratic Zeeman is varied towards the EA phase. Yet another deformation of the DBB into unmagnetized dark solitons solely occupying the m F = 0 spin-component takes place upon increasing this time the quadratic coefficient. Interestingly, the DBB solitons exist as unstable configurations for antiferromagnetic interactions, predominantly in the PO phase (c 1 > 0) but they are linearly stable within the EP phase which is characterized by c 1 < 0.
In this EP phase also nematic DDD solitons are realized. These solitary waves being stable configurations within the EP phase deform, as the quadratic coefficient is increased towards the PO phase, into stable unmagnetized dark solitons residing in the m F = 0 spin-component. This morphing, while analogous to the one observed for the DBB entities, here it appears for smaller quadratic Zeeman energies. Additionally, it is found that these dark waves, being rather persistent configurations, further enter the EA phase. This penetration holds up to a critical point, in terms of decreasing the quadratic Zeeman coefficient, below which the DDD waves gradually deform into two dark solitons equally populating the symmetric spin sublevels. These unstable DD solitons dynamically evolve into spin-waves which nevertheless retain their shared dark soliton at the center. The aforementioned spin-waves exhibit a finite spatial magnetization having a spatial extent proportional to the spin healing length.
Several extensions of the present work can be put forth. As a first step one can consider quenches between the first and second order phase transition boundaries to investigate how the relevant in each phase spin-1 soliton solutions migrate or transform from one phase to the other. Yet another interesting perspective would be to study interactions [22,43,50] between the spinor solitons identified within each phase of the above-obtained phase diagram or even unravel lattices consisting of multiple DBB, DDB and DDD spinorial solitons in analogy to the two-component settings e.g. of Refs. [21,25]. In some cases where different solutions co-exist (e.g. the DBB and the DDD in the EP phase), one could even consider collisions between different types of entities. Another aspect that the present work motivates further concerns the study of domain wall configurations in suitable regimes of the relevant phase diagram (e.g., within the EA state). Furthermore, generalizing the phase diagram of nonlinear excitations extracted herein in higher dimensions where vortex-bright states [69,70] are expected to form would be also a fruitful future direction.