{'fraction_non_alphanumeric': 0.04672699370270413, 'fraction_numerical': 0.017674763189924326, 'mean_word_length': 4.414191576735746, 'pattern_counts': {'":': 0, '<': 24, '<?xml version=': 0, '>': 23, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Stability and dynamics across magnetic phases of vortex-bright type excitations in spinor Bose-Einstein condensates

The static properties, i.e., existence and stability, as well as the quench-induced dynamics of vortex-bright type excitations in two-dimensional harmonically confined spin-1 Bose-Einstein condensates are investigated. Linearly stable vortex-bright-vortex and bright-vortex-bright solutions arise in both antiferromagnetic and ferromagnetic spinor gases upon quadratic Zeeman energy shift variations. Their deformations across the relevant transitions are exposed and discussed in detail evincing also that emergent instabilities can lead to pattern formation. Spatial elongations, precessional motion and spiraling of the nonlinear excitations when exposed to finite temperatures and upon crossing the distinct phase boundaries, via quenching of the quadratic Zeeman coefficient, are unveiled. Spin-mixing processes triggered by the quench lead, among others, to changes in the waveform of the ensuing configurations. Our findings reveal an interplay between pattern formation and spin-mixing processes being accessible in contemporary cold atom experiments.


I. INTRODUCTION
It is nowdays possible to controllably create Bose-Einstein condensates (BECs) possessing internal degreesof-freedom [1][2][3][4][5]. These multi-component systems, due to the Zeeman splitting of the involved magnetic sublevels are known as spinor condensates and have been discussed in dedicated reviews [6,7] and books [8][9][10]. Among spinors with hyperfine spin F = 1 or 2, spin-1 BECs represent arguably the most studied class. The two-body interaction of spin-1 bosons features density (or interparticle) and spin-interactions. By engineering the internal states using optical and magnetic fields, various magnetic ground states and the related to them first and second order phase transitions are now accessible [6]. For instance a 23 Na spinor gas experiences antiferromagnetic (AF) interactions [1,2] whilst 87 Rb [3,4] and 7 Li [5,11] feature weak and strong ferromagnetic (FM) ones.
Given the enhanced theoretical and experimental [19,20,22,24] recent interest in spinor BECs and the different excitations that can form in their distinct magnetic phases, we hereby consider harmonically trapped quasi-2D, spin-1 BECs featuring either AF or FM spininteractions. Concerning the static properties of the two setups under consideration, we tackle spinorial stationary states that bear at least one vortex component being filled by bright solitons. In comparison to earlier studies [39], a central feature of our work is that we consider vortical states of the same charge and zero net magnetization. Also, a key property of the structures of interest herein is the filling of vortices with bright components when the parameters of the system permit it (see details below). The understanding of the stability properties of such configurations, being addressed herein via a gener-alized Bogoliubov-de Gennes (BdG) theory [10,51,52], is still far from complete. Only partial results of this kind exist, as e.g. in the recent study of [49] where the maximal growth rate of the so-called nematic spin vortex state is provided. Here, we build on earlier findings based on simpler one-dimensional settings [23,53], in order to obtain the phase diagram of the identified vortical states, coupled with their corresponding potential instabilities. Triggering the latter can also be valuable, as it is strongly suggested by recent single-component BEC experiments [54], for designing certain topological states in the different spinor phases examining thereafter their dynamical response and spin-mixing processes [54,55].
Our findings indicate that vortex-bright-vortex (VBV) and bright-vortex-bright (BVB) excitations exist as stable configurations for either AF or FM spin-dependent interactions [see Fig. 1(a)-(b)]. These excitations experience structural deformations upon quadratic Zeeman (QZ) energy shift variations and importantly they feature narrow QZ intervals where oscillatory instabilities occur [51,56].
Dynamical evolution of perturbed VBV (BVB) entities entails, among others, their irregular (regular) precessional motion, nucleation of cross-shaped spinor patterns, and potential spiraling of the ensuing waveforms. These are findings evincing that spinor BECs provide a fruitful platform for probing instability-related spontaneous pattern formation [57,58]. Further, quenchinduced spin-mixing processes are unveiled under QZ energy shift variations at finite temperatures. The inclusion of thermal effects is inspired by their relevance in recent experiments [59,60]. Specifically, population transfer mechanisms are shown to be enhanced for larger values of the QZ coefficient and higher temperatures. Finally, the nonequilibrium dynamics of the vortical spinor configurations reveals the generic activation of their precessional motion, but also deformations where spinors simultaneously exhibit characteristic spatially anisotropic elongations.
The workflow of the present effort is as follows. Section II sets up the model mean-field equations of motion and the linearization method utilized herein. Section III contains our main findings regarding the existence, stability and dynamics of AF and FM spin-1 BECs. Their quench dynamics at finite temperatures is discussed in Sec. IV. In Sec. V we provide a summary of our results and a list of interesting perspectives for future investigations. Appendix B, elaborates on the impact of highercharge vorticity generalizing earlier instability findings occurring in single-component settings [61], demonstrating also dynamical triangular pattern formation.

A. Mean-Field equations
We consider a spin-1 BEC of N = 10 4 87 Rb [22,62] or 23 Na [19] atoms of mass M . [63] A uniform magnetic field B is applied along the transversal z-direction, and the system is confined in a quasi-2D harmonic trap. The quasi-2D trap is of the form V (x, y, z) = M ω 2 (x 2 + y 2 )/2 + M ω 2 z z 2 /2, obeying the condition ω z >> ω. Here ω z denotes the out-of-plane oscillator frequency, i.e., the one along the z-direction, and ω refers to the frequency in the (x − y)−plane (alias in-plane oscillator frequency). The corresponding three-component wave function, Ψ(r; t) = (Ψ 1 (r; t), Ψ 0 (r; t), Ψ −1 (r; t)) with r ≡ {x, y, z}, represents the distinct spin-components, m F = ±1, 0, of a spin-1 BEC. Additionally, throughout this work we choose as characteristic length and energy scales the in-plane oscillator length l osc = /M ω and ω respectively. Accordingly, space and time coordinates are rescaled as x ′ = x/l osc , y ′ = y/l osc , z ′ = z/l osc and t ′ = ωt respectively and the wave function as Ψ mF (x ′ , y ′ , z ′ ) = (l 3 osc /N )Ψ mF (x, y, z). However, due to the quasi-2D geometry of the potential considered herein (i.e. ω z >> ω) the aforementioned threedimensional wave function can be factorized as follows is the normalized GS wave function in the z-direction, and Ψ mF (x ′ , y ′ , t) is the quasi-2D wave function. The latter, with the above choices and rescaling (and dropping the primes for convenience) is described within the mean-field framework by the following dimensionless system of three coupled Gross-Pitaevskii equations (GPE) [7,18,22] In the above equations, H = − 1 2 ∂ 2 x + ∂ 2 y +V (x, y) is the single particle Hamiltonian with V (x, y) = (x 2 + y 2 )/2 denoting the 2D harmonic potential. Moreover, c 0 and c 1 are the so-called spin-independent and spin-dependent interaction coefficients given by c 0 = 2N √ 2πκ(a0+2a2) 3losc and respectively, in the units adopted herein. κ = ω z /ω is the anisotropy parameter, while the scattering lengths a 0 and a 2 account for collisions between two atoms belonging to the scattering channels with total spin F = 0 and F = 2 respectively. Additionally, c 0 > 0 (c 0 < 0) accounts for repulsive (attractive) interatomic interactions, while c 1 > 0 and c 1 < 0 designate AF and FM spin-interactions, respectively. Furthermore, the QZ energy shift, q, can be determined via the relation q = µ 2 B B 2 /(4 ωE hfs ), where µ B denotes the Bohr magneton and E hfs is the hyperfine splitting. Notably, q can be tuned experimentally either by adjusting the external magnetic field B [64] or the hyperfine splitting E hfs by utilizing a microwave dressing field [65,66].
Moreover, the total number of particles, 1 ≡ mF dx dy |Ψ mF (x, y, t)| 2 , is preserved with the population fraction of each spin component being defined as and satisfying 0 ≤ n mF ≤ 1. Throughout this work we prescribe that the (similarly conserved quantity of the) net magnetization along the z-direction i.e., M z = dx dy |Ψ +1 | 2 − |Ψ −1 | 2 , remains zero. This, in turn, implies that there is no population imbalance between the symmetric m F = ±1 components.
Below, the in-plane trapping frequency is set to ω = 2π × 20 Hz and the transverse one to ω z = 2π × 400 Hz. This leads to an anisotropy parameter κ = 20 inspired by recent 2D BEC experiments, see, e.g., Ref. [57]. Additionally, for AF interactions, a BEC of 23 Na atoms is considered having mass M = 23amu, s−wave scattering lengths a 0 = 2.52862nm, a 2 = 2.77196nm and therefore, c 0 ≈ 0.013N and c 1 ≈ 0.00039N [6,7]. For FM interactions, a BEC of 87 Rb atoms is employed with mass M = 87amu, a 0 = 5.387nm, a 2 = 5.313nm and thus c 0 ≈ 0.05N and c 1 ≈ −0.00023N . The QZ coefficient, q, is typically varied within the interval [−3, 3]. The latter, has been identified to be a representative interval of the principal phenomenology of interest. Unless stated otherwise, the total particle number and the vortex charge are fixed to N = 10 4 and S = 1 respectively.

B. Vortex-bright spinor ansatz and BdG approach
Initially [Sec. III], we focus on obtaining stationary solutions of the spinor system of Eqs. (1)-(3) in the form of vortex-bright (VB) solitons [52,[67][68][69] that can occupy all three hyperfine components by utilizing a Newton-Krylov iterative scheme [70]. Specifically, in order to introduce a vortex (V) of charge S and a bright (B) soliton in the desired m F component, the following ansatz is applied to the relevant wave functions In Eq. (5), H m (x) = (−1) m e x 2 d m dx m e −x 2 and H n (y) = (−1) n e y 2 d n dy n e −y 2 are the mth-and nth-order Hermite polynomials respectively. A singly quantized vortex can be obtained by employing as an initial guess the (m, n) = (1, 0) polynomial namely the first excited state for the real part of the relevant wave function, and the (m, n) = (0, 1) for the imaginary part, respectively. In a similar vein, e.g. a doulby quantized vortex (S = 2) is realized by a suitable combination of (m, n) i.e. by using (m, n) = (2, 0) − (0, 2) for the real part, while (m, n) = 2(1, 1) for the imaginary one. Subsequently, in sections III and IV, the stability properties and the quench-induced dynamics of the previously identified equilibrium states are investigated. Notice that we restrict our investigations to the case where the components contain vortices of the same charge S. However, it would be worthwhile to consider in the future also cases in which e.g. the symmetric spin states include oppositely charged vortices in order to unravel the creation of patterns analogous to the monopoles appearing in threedimensions.
For studying the stability of the VBV and BVB configurations found herein, a spectral BdG analysis suitably generalized for 2D spinorial BECs is performed [10,23,51,52]. In delineating the latter, we note that it consists of perturbing the iteratively identified stationary states, Ψ 0 mF (x, y), of each phase via the ansatz Here, ǫ is a small amplitude perturbation parameter and µ mF with m F = 0, ±1 is the chemical potential of each spin-component. Ω and (a mF , b ⋆ mF ) T denote, respectively, the eigenfrequencies and eigenfunctions of the resulting eigenvalue problem that one obtains upon substituting Eqs. (7) into the system of Eqs. (1)-(3) and keeping terms of order O(ǫ) [10,51,52]. Namely, In the above expression λ ≡ −iΩ and M j (with j = 1, . . . , 9) are 2 × 2 matrices whose explicit form is provided in Appendix A. The resulting eigenvalue problem of Eq. (8) is subsequently solved numerically. Note that in 2D spinor condensates BdG analysis of vortical configurations bearing also a bright soliton component is still elusive and only partial results to that effect are available, to the best of our knowledge. On the dynamical side, in order to study alterations of the stationary states existing in a specific phase when crossing a phase boundary [71], a quench of the QZ energy shift is applied. The quench is performed from an initial (pre-quench) q ≡ q i to a final (post-quench) value q ≡ q f in a way that assures penetration to a different phase. To seed population transfer in the quench dynamics, the commutator of the total spin operator with the Hamiltonian has to be nonzero and we achieve this by including dissipation into the system. Such dissipation, can naturally arise in BEC experiments when a non-negligible thermal gas component is present in the system. Furthermore, in the large particle limit that we operate it is expected, in line with recent spin-1 BEC experiments [19,20,22,24], that quantum fluctuations are suppressed. For the dynamical evolution of the spinorial system a fourth-order (in time) Runge-Kutta method is used with temporal and spatial discretization dt = 10 −4 and dx = dy = 0.05 respectively, while a (2nd order) finite difference scheme is utilized for the spatial derivatives.

III. STATIC PROPERTIES OF VBV AND BVB SPINOR EXCITATIONS
A. Antiferromagnetic vortex-bright type configurations To tackle the nonlinear excitations of the VB form that arise in the distinct phases of 2D harmonically confined spin-1 BECs, an initial guess provided by Eqs. (5)-(6) is introduced to the time-independent version of the system of Eqs. (1)-(3). Specifically, for AF interactions (c 1 > 0), it is well-known that two distinct phases exist depending on the value of the QZ energy shift [6,7]. Namely, for q < 0 the AF phase is realized while for q > 0 the system resides in the PO phase. In the former phase and at the GS level, only the symmetric m F = ±1 spin-components are populated.
Thus, a natural choice for accessing the corresponding excited states is to consider an initial guess where vortices (bright solitons) are embedded in the m F = ±1 hyperfine states and a bright soliton (vortex) occupies the m F = 0 component. It turns out that among these two, i.e. VBV and BVB, configurations only VBV excitations exist in the AF phase. Representative density profiles are illustrated as insets in Fig. 2(a 1 ) − (a 3 ). We note in passing that for all vortex entities to be presented throughout we have verified that they are accompanied by the expected 2πS phase winding (with S denoting the vortex charge). Recall that the polarization, is a measure of population transfer phenomena. It obeys −1 ≤ P ≤ 1, when all three m F components (3C) are populated, but P = 1 (P = −1) if only the m F = 0 (m F = ±1) state(s) is (are) populated yielding a single (two) component, 1C (2C), configuration. The aforementioned 3C stationary states exhibit polarizations −1 < P < 1 (see orange line in Fig. 2) and their interval of existence is provided in the last column of Table I in Fig. 1(b). VBV excitations are further found to deform upon a q variation into highly localized vortices occupying the symmetric spincomponents as q decreases [ Fig. 2(a 4 ), (a 5 )]. These 2C vortices are indeed characterized by P = −1 and they exist for all values of q < −2.5 that we have checked, see also second column of Table I in Fig. 1(b). Yet another deformation occurs for the VBV configurations but upon increasing q. In this case, each vortex core gradually becomes wider in order to effectively trap [68] the accompanying wider bright soliton of the m F = 0 spincomponent. This alteration holds until the 1C GS of the PO phase is reached that is, in turn, characterized by P = 1 [first column of Table I in Fig. 1(b)]. Notice the abrupt jump of P from P ≈ −0.48, q = −0.3 to P = 1,  Having examined the existence of VBV excitations along with their relevant structural deformations we next explore the stability properties of such configurations. In contrast to earlier predictions mostly focused on energy based considerations [39,[47][48][49] below we utilize a generalized BdG theory to microscopically determine the involved internal modes. As stated earlier, to perform the BdG analysis the ansatz of Eq. (7) is used for this specific stationary solution. The relevant BdG spectra, obtained upon solving the eigenvalue problem of Eq. (8) associated to the VBV solutions, are depicted in Fig. 3(a 1 ) − (b 1 ). Note, that there exist in the spectrum three different pairs of modes lying at the Re(Ω) axis around the origin of the Re(Ω) − Im(Ω)-plane, i.e., at Re(Ω) = Im(Ω) = 0. These zero eigenfrequencies, not visible in the scales shown, are generated by continuous symmetries. The spinor system under study preserves the total particle number (phase invariance of the equations of motion), the magnetization and further has rotational symmetry, thus explaining the existence of these three pairs. Besides the aforementioned modes, two additional negative energy ones appear among the remaining modes of the discrete spectra that are denoted by light blue circles. The two distinct trajectories, obtained with respect to q, of these so-called anomalous modes (AMs) can be discerned in Fig. 3(a 1 ). Each of these modes is known to correspond to the precession of each of the two vortices within the parabolic trap [67,68]. Additionally, these AMs are quantified through their negative energy or negative Krein signature [51] which for the 2D spinor system reads It should be marked here that the existence of these modes is an immediate byproduct of the fact that the stationary states found herein are excited states of the spinor system. Namely, such modes would be absent in the case of the system's GS. Moreover, as long as these eigenfrequencies maintain their real nature, then their negative Krein signature further indicates that while a stationary solution is dynamically stable, it is simultaneously unstable thermodynamically [10]. The latter, in turn, implies that given a channel of energy dissipation, as in the case of the dissipative spinor system that will be discussed below, these eigendirections will be activated leading to an instability of the ensuing configuration. Notice that upon increasing q so as to reach the phase transition point (q = 0), in the vicinity of the latter, the aforementioned negative energy modes decrease in frequency, with both crossing the zero frequency axis around q ≈ −0.2. At the same time also a decreasing in frequency positive energy mode crosses Ω = 0 and leads to the appearance of the BdG spectra of 3C (a1)-(b1) VBV and (a2)-(b2) BVB stationary states upon a q variation for c1 > 0. In both cases the anomalous modes (AMs) are depicted by light blue circles while the background ones by black circles. Notice the two AM present for VBV structures when compared to the single pair occurring for BVB configurations. finite imaginary part, Im(Ω) = 0, shown in Fig. 3 The destabilization of the deformed VBV configuration is followed by a change in the Krein signature of the two (previously) negative energy modes from negative (light blue circles) to positive (black circles).
In addition to the above stability analysis results, there exist narrow intervals of q where oscillatory instabilities [56] take place for the VBV solution. In general, this type of instability stems from collision events involving pairs of positive and negative Krein signature modes resulting in eigenfrequency quartets and also possessing a finite imaginary component Im(Ω) = 0 [23,56]. We must emphasize here, that this is yet another key feature related to the theory of AMs: namely, their role in the manifestation of instabilities even in the absence of finite temperatures. Three such collision events can be readily seen in the BdG spectrum of Fig. 3(a 1 ) appearing e.g. at q = −1.5, q = −0.6 and q = −0.4. The first two are associated with the higher-lying anomalous mode whose absence for these values of q is transparent while the last one entails the collision and disappearance of both negative energy modes.
Two case examples are considered below for q = −0.6, demonstrating the activation of e.g. the lower-lying anomalous mode (AM 1 ) along with exploring the oscillatory instability present for this value of q. Particularly, Fig. 3(c 1 ) − (g 1 ) illustrate 2D contours quantifying the density difference between a perturbed and an equilibrium solution ∆Ψ mF (x, y) ≡ |Ψ mF (x, y)| 2 −|Ψ mF (x, y)| 2 . The perturbation here, consists of adding to the VBV stationary state the eigenvector associated either with AM 1 or with the eigenfrequency quartet identified for q = −0.6. Notice the two-lobe structure imprinted in ∆Ψ mF (x, y) resembling a 2p orbital-like configuration. The lobes are centered around the origin of the (x − y)− plane being parallel to the y = 0 axis. They are further found to be asymmetric with respect to x = 0 with ∆Ψ −1 (x > 0, y) > 0, Fig. 3 Fig. 3(d 1 )]. Moreover, the m F = +1 component (not shown) has the same effect with that of m F = −1 when the VBV is perturbed via AM 1 but ∆Ψ +1 (x, y) is complementary to ∆Ψ −1 (x, y) when the VBV is perturbed via the AM 2 mode. However, this is not the case when considering the quartet scenario [ Fig. 3(e 1 ) − (g 1 )]. The predominant effect of this mode is the asymmetric distribution of ∆Ψ ±1 (x, y) with respect to y = 0 being Both components are azimuthally deformed exhibiting a counterclockwise rotation. The m F = 0 one practically remains unaffected, with ∆Ψ 0 (x, y) ∼ 10 −4 featuring an asymmetry along x = −y. Finally, it is worth commenting here, that dynamical evolution of the excited, with AM 1 , VBV entity leads to its precessional motion where the entire VBV rotates around the trap. Whilst, exciting the configuration with AM 2 results in a rotating cross-shaped pattern in which the vortex components perform an anti-phase oscillation among each other and the m F = 0 bright component remains unaltered. This anti-phase vibration leads, in turn, to an overall breathing of the BEC background.
For AF interactions but for q > 0, namely within the PO phase, the preferable configuration consists of a solely occupied m F = 0 spin-component. Since this component, according to the GS of the system [7], is expected  to become the majority one, in our search for nonlinear excitations arising in this phase we choose to imprint a vortex on it. Consequently, bright solitons are plugged in the remaining m F = ±1 spin-components. With such an initial guess, indeed, BVB stationary solutions are captured for 0.3 ≤ q < 1.95 [see also the relevant third column of Table I in Fig. 1(b)]. Characteristic density contours of such a BVB structure are presented as insets in Fig. 2 . Notice that similarly to the VBV configurations, the BVB stationary states are characterized by −1 < P < 1 (see the purple line in Fig. 2) and they also experience two deformations with respect to q. One deformation is rather gradual as captured by the slope of the polarization as q increases, leading to a single highly localized vortex occupying the m F = 0 hyperfine state [see Fig. 2(b 6 ) and the relevant first column of Table I in Fig. 1(b)]. On the contrary, as q decreases towards the first order transition boundary (q = 0) separating the PO and the AF phase, an abrupt deformation of the BVB configuration to the 2C state, reminiscent of the GS of the AF phase, occurs [see Fig. 2(b 4 ), (b 5 ) and the second column of Table I in Fig. 1(b)] around q = 0.3. Notice, that both the VBV and the BVB configurations feature smooth deformations towards the 2C and the 1C vortex state respectively. In the opposite q direction a sharp transition takes place when the relevant phase boundary is approached to 1C and 2C zero vortex states respectively. This behavior of the polarization is in direct contrast to the corresponding sharp transition occurring on the GS level, i.e. in the absence of nonlinear excitations (results not shown here for brevity) [6,7].
BVB excitations turn out to be linearly stable configurations for all values of q ∈ (0.4, 1.95), with a relevant example shown in the BdG spectrum of Fig. 3(b 2 ) e.g. for q = 0.5. Due to the single vortex contained in this configuration, only a single pair of negative energy modes is present in this spectrum. According to our discussion above, when activated, i.e., upon adding the associated to it eigenvector to the BVB solution, this mode leads to the precessional motion of the BVB structure. It is only for significantly deformed BVB configurations, namely for states where the bright soliton dominates the configuration corresponding to q ≤ 0.4, that oscillatory instabilities [like the one depicted in Fig. 3(a 2 )] appear. In order to appreciate the effect of the emergent eigenfrequency quartet on the BVB solution, we have added to the latter the corresponding quartet eigenvector. A close inspection of the associated density difference ∆Ψ mF (x, y) illustrated in Fig. 3(c 2 ) − (e 2 ), reveals that such an addition leads to an asymmetric across the anti-diagonal (x = −y) BVB structure having ∆Ψ ±1 (x > 0, y) > 0 and ∆Ψ 0 (x > 0, y) < 0. In all cases a counterclockwise rotation takes place that is in turn related to the precessional motion of the entire BVB entity observed in the dynamics. Finally, the anomalous mode ceases to exist for q < 0.3 signaling the transition to the GS of the AF phase. Moreover, we emphasize at this point that the robustness of stable VBV and BVB stationary states  has been also dynamically confirmed by monitoring their spatiotemporal evolution for times up to t = 2.0s.

B. Ferromagnetic VBV and BVB spinors
Turning to FM spin-interactions (c 1 < 0) three phases can be realized as q is varied, supporting GS with an occupancy ranging from 1C to 3C [6]. In particular, the socalled 1C fully magnetized along the +z (−z)-direction easy-axis (EA) phase exists for q < 0. Since we operate in the regime where the harmonic oscillator length is smaller than the spin-healing length, phase separation is absent in our setting. The case where the relevant inequality is reversed, while interesting in its own right, is outside the scope of the present work and hence deferred to future studies. The 3C easy-plane (EP) phase occurs for 0 < q < q T and the 1C PO phase is characterized by q ≥ q T [6,12,23]. In the latter two inequalities q T = 2c 1 n (which equals 0.05, for our chosen parameters) designates the threshold between the involved phases with n being the peak density at the trap center. In this FM spinor setting, VBV stationary states are identified for −2.74 ≤ q < −0.14 [third row of Table  I in Fig. 1(b)]. These states possess zero net magnetization and −1 < P < 1 as shown in Fig. 4. They also have density profiles, |Ψ mF (x, y)| 2 , similar to their AF siblings [ Fig. 4(a 1 ) − (a 3 )]. Strikingly, FM VBV waves are more persistent configurations when compared to their AF counterparts. They are seen to penetrate deeper into the EA phase before deforming into a 2C vortex [ Fig. 4(a 4 ), (a 5 )] structure for smaller q values [third row of Table I in Fig. 1(b)]. They further transform slower to the PO GS [ Fig. 4(a 6 )] following an increment of q towards the phase transition boundary (q = 0). As such the corresponding polarization curve is found to be rightshifted thus being closer to the origin when compared to the relevant AF one.
For 0 < q < q T , i.e. within the EP phase, the existence of BVB stationary states is also unveiled and presented in Fig. 4. It is noteworthy that FM BVB structures also feature larger q intervals of existence in comparison to their AF analogues [fourth row of Table I in Fig. 1(b)]. These structures penetrate the PO regime with the underlying 3C densities as depicted in the insets of Fig. 4 . Recall that at the GS level the PO phase exists for q ≥ q T . Eventually, the 3C BVB structure deforms into the 1C vortex configuration illustrated, e.g., for q = 2.5 in Fig. 4(b 6 ). The existence of these states is (parametrically) prolonged also following a decrease of q until a 2C Thomas-Fermi state is reached within the EA phase [ Fig. 4(b 4 ), (b 5 )]. This has as a result, a left-shifted polarization curve that is closer to the origin when compared to the relevant AF one.
Investigating the stability of both configurations we find that, as their AF counterparts, VBV and BVB stationary states experience stable intervals of existence. This result can be verified by inspecting the BdG spectra shown in Fig. 5(a 1 ) for the VBV solution and in Fig. 5(a 2 ) and (b 2 ) for the BVB one. Notice that in both cases and for the parametric intervals shown, all eigenfrequencies maintain their real nature, i.e., Im(Ω) = 0. However, these structures further feature narrow q intervals where oscillatory instabilities occur. One such example is presented regarding the VBV entity for q = −1.2 in the BdG spectrum of Fig. 5(b 1 ). Similarly to the AF cases discussed above, also here the emergence of an eigenfrequency quartet is observed, that owes its presence to the collision of the higher-lying negative energy mode, AM 3 , with a positive energy one. Importantly though, and also in sharp contrast to the AF VBV solutions, three instead of two AMs appear in the spectrum of this configuration. As stated earlier, since two vortices participate in this configuration two anomalous mode pairs are to be expected for this stationary state. Thus, we initially investigate further the presence of the lowest-lying mode, namely AM 1 . This mode appears remarkably close to the zero eigenfrequency axis and remains near the latter as q is varied till its destabilization slightly below the threshold separating the EA and the EP, i.e. at q = −0.2. ∆Ψ ±1 (x, y) has a four lobe spatial distribution closely resembling a 3d xy orbital configuration which is further found to be slightly rotated counterclockwise with respect to the x = 0 axis. A similar outcome is evidenced for S = 2 VBV solutions as discussed in Appendix B and visualized e.g in Fig. 10 (h 1 ). On the other hand, AM 1 has a vanishing impact on the relevant bright component, with ∆Ψ 0 (x, y) ∼ 10 −8 . Particularly, AM 1 leads dynamically to an anisotropic spatial elongation of the two vortices that perform a precessional type of motion but with the vortices in the m F = ±1 components rotating with a π phase difference among each other and a bright soliton that remains put throughout the evolution. As such, this is a mode involving inter-component dynamics, rather than the intra-component ones, associated with the vorticity of the VBV structure. Next, we appreciate the effect that the remaining two AMs have on VBV solutions while we note that their destabilization takes place at q = −0.05. Considering the eigenvector related to AM 2 results in an asymmetric 2p orbital-like distribution of ∆Ψ mF (x, y), with the two lobes oriented along the anti-diagonal x = −y as showcased in Fig. 5(c 1 ), (d 1 ). It also holds that ∆Ψ −1 (x > 0, y) < 0 and ∆Ψ 0 (x > 0, y) > 0. Note that a similar 2p orbital configuration is also obtained for FM S = 2 VBV spinors (see the relevant discussion around AM 3 and AM 4 in Appendix B). This mode leads upon activation to the normal or regular precession of the VBV structure. Namely, the two vortices are on the same side and oscillate around the trap center with the bright soliton following their motion. A much more drastic deformation is evidenced when the solution is perturbed through the eigenvector of AM 3 leading to an asymmetric azimuthally rotated ∆Ψ mF (x, y) for the symmetric m F = ±1 vortex components analogous to the one found for AF VBV equilibrium states [see Fig. 3(e 1 ), (g 1 )]. Also here, ∆Ψ 0 (x, y) ∼ 10 −6 has a vanishing effect for the bright soliton component. As we shall show in the dynamics below, once excited, the mode AM 3 leads to a different form of precession of the VBV solution.
Here, the precession of the VBV consists of two vortices hosted in the m F = ±1 being anti-diametrically located with respect to the center and performing oscillations that have a π phase difference with respect to one another, while the m F = 0 bright soliton component remains intact. However, this motion becomes responsible for an instability when AM 3 collides with a positive Krein background mode. Recall that whenever such a collision takes place an eigenfrequency quartet occurs in the BdG spectrum instead of the ensuing AM pair. Indeed, notice that e.g. AM 3 is absent in Fig. 5(b 1 ) giving rise to the observed quartet. In this latter case as it is shown in Fig. 5(e 1 )-(g 1 ), a spiral is imprinted in the density difference ∆Ψ ±1 (x, y) being of a complementing nature among these two hyperfine components, yet minuscule for the m F = 0 one [ Fig. 5(f 1 )]. This leads in turn dynamically, to a spiraling of the 2D VBV entity, an outcome caused by the oscillatory instability.
As an example for the BVB solution, we choose the one of a significantly deformed, i.e., close to threshold, BVB excitation [ Fig. 5(c 2 )-(e 2 )]. It turns out that, the bright soliton hosted in the m F = ±1 spin-components dominates the configuration for q = 0.2. This bright dominated entity is additionally found to be significantly  broadened. Its width becomes comparable to the size of the background cloud, suggesting that the BVB character of this solution is lost. Perturbing this state with the eigenvector associated with the single -in this case-AM pair, leads to a two lobe asymmetric density difference resembling a 2p orbital for all three hyperfine states. The two lobes are oriented along the diagonal but experience an asymmetry, with ∆Ψ ±1 (x, y) < 0 for x = −y and ∆Ψ 0 (x, y) > 0, along the anti-diagonal. Featuring in this way, a similar yet inverted behavior to the one found for FM VBVs but also to FM S = 2 BVB spinors (Appendix B) when perturbed by Snapshots during the spatiotemporal evolution of this perturbed entity are provided in Fig. 6(a 1 ) − (c 5 ). As expected, the precessional motion of the entire BVB structure is observed from the initial stages of the dynamics, with the bright soliton m F = ±1 components remaining trapped in the course of the evolution around the vortex core, see Fig. 6(a 1 )−(c 5 ). For comparison, in the bottom panels of Fig. 6(d 1 )− (f 5 ), a perturbed VBV excitation via the eigenvector of AM 3 is presented for q = −1.0. Two key findings are worth commenting here. The one concerns the fact that even though the amplitude of the perturbation for both structures is the same, the precession of the VBV excitation is not as pronounced as the one observed for the deformed BVB solution. However, and even more importantly irregular precession is featured by the VBV structure with the two vortices being out-of-phase throughout their motion. This is an outcome that has a drastic effect also on the bright soliton which, contrary to the BVB state, now remains unaffected. Finally, in order to emulate the presence of a finite thermal fraction being usually present in cold atom experiments we introduced the following ansatz Ψ pert mF = Ψ 0 mF (x, y) [1 + εδ(x, y)] to the m F component wave function [72]. In this expression, ǫ accounts for the thermal fraction and δ(x, y) denotes a normally distributed perturbation with zero mean and variance unity [57]. Generically, this ansatz allows for the activation of the respective AM in the course of the evolution. Additionally, it should be noted that the AMs are converted to unstable eigendirections in the presence of a thermal fraction, correspondingly dominating the BEC dynamics, similarly to what is known, e.g., for two-component condensates [73]. This way, the destabilization mechanisms found above would be evident in a corresponding experimental realization.

IV. QUENCH DYNAMICS ACROSS MAGNETIC PHASES
Having explicated the static properties of VBV and BVB nonlinear excitations, in the following we aim at addressing alterations of the ensuing waveforms being subjected to quenches of the q parameter in order to cross the distinct magnetic phase boundaries (see also Fig. 2 and Fig. 4). To monitor the quench-induced dy- . Spin-mixing is triggered in all cases being more suppressed for VBV excitations as compared to the BVB ones and for both types of entities when the relevant transition entails quenches within the same phase. In all cases, the initial state configuration, having a pre-quench value q ≡ qi, refers to the underlying in each phase 3C stationary state transitioning either from 3C→2C or from 3C→1C states being characterized by different post-quench QZ coefficients q ≡ q f (see legends). All quenches are considered for fixed γ that is either γ = 0.0023 or γ = 0.01 (see legends namical evolution of the spinor gases at hand in an experimentally relevant fashion [5], we expose them to finite temperatures. Note that quenches are routinely utilized in spin-1 ultracold atom experiments to probe transition boundaries [59], spin-turbulence and the related to it half-quantum vortex generation [60] but also to study matter-wave jet formation [11]. Contrary to the above, here we use quenches at finite temperatures i) to activate the internal motion of the identified vortical spinors, ii) facilitate population transfer among the components and iii) study structural deformations of both BVB and VBV configurations across the distinct magnetic phases.
Representative examples among the extensive investigations performed herein, are presented in Fig. 7(a 1 )-(c 4 ) and Fig. 7(d 1 )-(f 4 ) regarding the density evolution for AF and FM spin-interactions respectively with γ = 0.0023. In the former case, we monitor the dynamics of an AF BVB excitation once quenched from the PO phase having q i = 0.3 towards the AF phase with postquench QZ coefficient q f = −1.0. It becomes apparent that population transfer from the m F = 0 to the m F = ±1 states takes place [see also Fig. 8(b 1 )] from the initial stages of the quench-induced dynamics triggering the precession of an initially stationary spinorial BVB structure. This motion is accompanied by a prominent elongation along with the instantaneous rotation of all three spin constituents. Moreover, the vortex experiences a structural deformation reminiscent of a doughnut-like pattern: an outcome that is further captured by the two mode motion of the relevant temporal evolution of the populations of the individual components illustrated in Fig. 8(b 1 ). This two mode motion is characterized by rapid oscillations of the populations and a long-time transfer (not shown in the presented timescales) where exchange of the populations between the m F = 0 and m F = ±1 takes place. Notice that the bright soliton m F = ±1 components remain trapped around the vor-tex core, following its composite motion throughout the evolution. Turning to FM interactions and upon considering a quench from q i = 0.15 (EP phase) to q f = 4.5 (PO phase) it is observed that the precessional motion constitutes the dominant dynamical mode, entailing an arguably faster spin-mixing process when compared to the aforementioned AF scenario.
In order to shed light onto the underlying spin-mixing processes triggered by the quench, a close inspection of the temporal evolution of the population of the individual components, n mF (t), is performed. Specifically, Fig. 8(a 1 ) − (c 2 ) and Fig. 8(d 1 ) − (f 2 ) capture the essence of our findings for a wide selection of pre-and postquench QZ energies and for distinct γ values. AF (c 1 > 0) and FM (c 1 < 0) condensates are treated on equal footing. For both spinor settings, transitions across the distinct magnetic phases are initiated from the relevant in each phase 3C VBV and BVB stationary states towards the corresponding 2C or 1C configuration.
Particularly, our key observations are the following. Irrespectively of the spinorial BEC system, spin-mixing processes are activated from the initial stages of the quench-induced dynamics. We find that population transfer occurs faster for larger post-quench values q f accessing this way states that are deeper in the relevant magnetic phase [ Fig. 8(a 1 ), (b 2 ) and Fig. 8(d 1 ), (e 2 )]. However, it is found to be more suppressed for VBV excitations as compared to BVB ones. This suppression occurs also for both types of entities when the relevant transition entails quenches within the same phase when compared to transitions that cross distinct phase boundaries. Additionally, spin-mixing is accelerated for a larger dissipation parameter γ being in turn related to higher temperatures, see for instance Fig. 8(c 1 ), (c 2 ) and Fig. 8(f 1 ), (f 2 ). We also remark that slightly enhanced intercomponent population transfer arises for AF rather than FM interactions as can be inferred by comparing Fig. 8(a 2 ) and Fig. 8(d 2 ) due to the larger spin-spin interaction in the former case. Finally, it is important to note here, that similar to the aforementioned findings occur during the nonequilibrium dynamics of higher charge excitations. However, in this case, the spin-mixing processes discussed above, are found to be relatively accelerated.

V. CONCLUSIONS AND FUTURE PERSPECTIVES
In the present work the existence, stability as well as the quench-induced dynamics of VB-type nonlinear excitations arising in 2D harmonically trapped spin-1 antiferromagnetic and ferromagnetic BECs have been explored. Our investigation has been focusing on variations of the quadratic Zeeman energy shift so as to access and subsequently cross the distinct magnetic phases of such settings. A systematic Bogoliubov de-Gennes linearization analysis has been utilized for the extraction of the sta-bility properties of the considered nonlinear excitations.
In particular, the existence of VBV and BVB stationary states has been exemplified, with the former being present in the antiferromagnetic and the easyplane phases for antiferromagnetic and ferromagnetic spin-interactions respectively. On the contrary, BVB solutions appear in the polar phase of either antiferromagnetic or ferromagnetic spinors. In this latter scenario, stable BVB structures are also found within the easyplane phase. In both settings deformations of the ensuing waveforms as the associated transition boundary is approached are explicated complementing this way the phase diagram of this type of nonlinear excitations in the (c 1 , q)−plane.
It turns out that independently of their flavor and also of their charge, the aforementioned entities exhibit stable intervals of existence that can be interrupted by narrow windows where oscillatory instabilities take place. Indeed, we have elaborated on the number of anomalous mode eigendirections that the structures bear and thus the number of potential instabilities, as well as illustrated when these instabilities may materialize as a result of collision of these anomalous modes with positive energy ones. We have also monitored the dynamical outcome of excitation of the different anomalous modes. The robustness or unstable dynamics of the above-described entities are confirmed accordingly, demonstrating for instance the precessional motion of VBV and BVB spinors and their structural deformation towards -among otherstriangular-shaped patterns.
We have further investigated the quench-induced dynamical evolution of the aforementioned threecomponent spinors at finite temperatures so as to appreciate the system's dynamical response. Here, it is found that spin-mixing processes occur faster for larger postquench quadratic Zeeman energy shifts and an increasing dissipation parameter. Also, population transfer is slightly enhanced when considering antiferromagnetic instead of ferromagnetic spin-dependent interactions. Monitoring the nonequilibrium dynamics reveals, among others, the activation of the precessional motion along with a spatial elongation of the spinorial nonlinear excitations, irrespectively of their specific nature and spin-interactions. The above processes are accelerated when higher charge vortices are contained in the spinorial configuration. The latter also bear a significantly larger number of anomalous modes and, thus, potentially unstable eigendirections.
There exist several extensions of the present work worth pursuing in future endeavors. A straightforward generalization would be to study the quench dynamics in a 7 Li spin-1 BEC where the strong ferromagnetic spininteraction would certainly enhance the spin-mixing processes which might be possibly associated with a richer pattern formation. A detailed investigation of e.g. S = 3 vortical spinors, that we barely touched upon herein, in symmetry broken settings would facilitate the engineering of exotic pattern formation with atomic orbital-like signatures. Additionally, exploring the interaction effects of vortex lattices as well as their stability and dynamics in spinor setups is of direct relevance, due to the potential of inclusion of external rotation [8,9]. Indeed, it is already of significant recent interest to explore the interaction of vortical patterns, as has been done in two-component settings, e.g., in [75,76] (see also references therein). Moreover, in the current setup the inclusion of three-body recombination processes as a dissipative mechanism in selective spin-channels constitutes a situation that accounts for possible experimental imperfections [5]. Yet another fruitful perspective is to consider domain-walls formed by two out of the three spin-components with the remaining one being a nonlinear excitation of different flavor, e.g. a vortex [77]. This setting will enable one to devise particular spin-mixing channels and consequently study dynamical pattern formation. The corresponding matrix elements A ij read Recall, that Ψ 0 mF (x, y) denotes the relevant for each magnetic phase equilibrium solution. Substituting Eqs. (A1) and Eqs. (A2) in the eigenvalue problem of Eq. (8) leads, upon numerical evaluation, to the BdG spectra given in the main text.
Appendix B: Impact of larger system sizes and higher-charge vorticity Here, we aim to generalize our findings presented in the main text by considering different system sizes and vortex charges. In particular, in the former case we systematically vary the total number of particles within the range N ∈ [1 × 10 3 , 2 × 10 4 ] while in the latter situation vortices of S = 2, 3 are explored. Experimentally higher-charge vortices can be realized using the topological phase-imprinting technique [78]. Remarkably enough, by monitoring the polarization of the FM spinor system under (q, N ) variations reveals that it remains insensitive under such parametric changes independently of the stationary configuration (not shown for brevity). Sizable deviations are only present when higher charge vortices are contained either in a VBV or a BVB equilibrium solution. Indeed, as presented in Fig. 9(a)-(b), P experiences drastic changes under a (q, S) variation. Particularly, while S increases an overall shift of P towards more positive (negative) q values is observed for BVB (VBV) solutions altering in this way the distinct magnetic phase transition boundaries. Since S = 2, 3 vortices are structures having significantly wider cores, see Fig. 9(c 1 )-(c 3 ) and Fig. 9(d 1 )-(d 3 ) respectively, when compared to the S = 1 configurations (see the insets in Fig. 4), the abovementioned shift can be explained as follows. Initially, we should recall that bright solitons can only be sustained in repulsive environments, via their effective trapping by nonlinear excitations such as the vortices studied herein [68]. Thus, higher charge vortices can effectively trap in a more efficient manner the bright soliton component leading in turn to persistent over wider parametric intervals 3C entities.
Even though it is known that multiply-quantized vortices are prone to decay into singly quantized vortex pairs in scalar [61,[78][79][80][81] and two-component BECs [82], the fate of such higher charge entities in spinorial BEC systems remains still elusive [35]. As such, below we further investigate the stability properties of these configurations. Specifically, we focus on the simplest case scenario, namely the one involving spinors in which the vortices have charge S = 2. Our stability analysis reveals that doubly quantized FM VBV and BVB are, in principle, linearly stable configurations for values of q ∈ [−4.0, −0.5) and q ∈ [0.2, 4.0] respectively that we have checked and for the particular particle number chosen. Narrow windows where oscillatory instabilities are identified, giving rise to a finite imaginary contribution of the order of Im(Ω) ∼ 10 −3 − 10 −2 , occur for the VBV configuration e.g. for q ∈ [−0.75, 0.65] and q = −0.9. Remarkably, seven negative energy modes, AM i (i = 1, 2, . . . , 7), are found in the BdG spectrum of this structure as can be seen for instance in Fig. 10(a 1 ) for q = −1.0.
Among these modes the lowest-lying one, AM 1 , resides close to the zero frequency axis, as in the FM S = 1 scenario. In order to visualize the effect that the perturbation has on the VBV excitation, we invoke, as in the main text, the density difference, ∆Ψ mF (x, y) ≡ |Ψ mF (x, y)| 2 − |Ψ mF (x, y)| 2 . It turns out that contrary to the S = 1 case here ∆Ψ ±1 (x, y) develops an eight lobe dumbbell-shaped structure centered around the origin of the (x − y)−plane (x = y = 0) and being a reminiscent of a 5g z 3 x orbital. This density difference is further found to acquire its maximum/minimum value in an alternating fashion among the distinct lobes. Importantly though, also for higher charges, the number of negative Krein modes is greater than the one anticipated for an S = 2 VBV solution. Indeed, it is known [82] that since the two vortices are doubly quantized in this case one can assign two anomalous mode pairs to each of the two participating vortices. These yield in turn four anomalous mode pairs for such a state rather than the seven identified herein.
Thus, in what follows ∆Ψ mF (x, y) is evaluated and shown in Fig. 10(b 1 ) − (j 1 ) for three out of the seven modes that VBV solutions possess. Notice that in all three cases the bright soliton of the m F = 0 spincomponent is not altered as captured by ∆Ψ 0 (x, y) ∼ 10 −6 . This is in contrast to the vortices of the m F = ±1 spin-components that complement one another. Evidently, perturbing the VBV solution with the eigenvector related to AM 4 results in an asymmetric two lobe ∆Ψ ±1 (x, y) configuration resembling a 2p orbital. The latter, is oriented along the anti-diagonal x = −y but is slightly shifted from it counterclockwise. ∆Ψ −1 (x > 0, y) > 0 while ∆Ψ +1 (x > 0, y) < 0. AM 5 leads to a centered around the origin 4f xz 2 orbital-like configuration. Namely, a six dumbbell-shaped lobe structure [ Fig. 10(e 1 ) and (g 1 )]. Notice that the density difference maximizes and minimizes in an alternating manner as we go from one lobe to the other. Here, dynamical activation of AM 4 unveils the formation of anti-phase triangular patterns in the vortex m F = ±1 components which along with an intact bright soliton m F = 0 component precess around the trap center. Contrary to the above dynamics, perturbing the VBV entity with AM 5 leads to the formation and robust propagation of a deformed structure. The two vortices perform an irregular out-of-phase precession leaving in this way the bright soliton in the m F = 0 component intact, but instead of forming triangles, they feature dipolarly elongated density distributions being inverted between the m F = +1 and m F = −1 components. However, addition of the eigenvector associated with AM 7 entails a completely different deformation. ∆Ψ ±1 (x, y) develops a 3d xz orbitallike pattern [ Fig. 10(h 1 ) and (j 1 )]. That is, a four lobe cloverleaf distribution with the symmetric hyperfine components complementing one another.
Dynamical activation of AM 7 leads to a breathing core VBV structure that performs an irregular (out-ofphase) precession having spatially anisotropic and oppositely elongated with respect to each other symmetric spin components. The remaining eigenvectors associated with AM 2 , AM 3 and AM 6 result respectively in a ∆Ψ mF (x, y) that has a 4d xz orbital structure in all three hyperfine components, having ∆Ψ 0 (x, y) ∼ 10 −8 and being centered at the origin of the (x − y)−plane. AM 3 leads to a 2p orbital distribution like the one found for the AM 4 mode but with the two lobes being slightly shifted with respect to each other while residing anti-diametrically along the diagonal x = y. Here, ∆Ψ −1 (x, y) = ∆Ψ +1 (x, y) and both are complementary to the m F = 0 bright soliton component. Additionally, the effect of AM 6 closely resembles that found for AM 7 but with the symmetric vortex components having now exactly the same structure while being complementary to ∆Ψ 0 (x, y) which is now finite. Finally, we note that AM 1 and AM 2 perform an eigenfrequency zero crossing at q = −0.4 but are not responsible for an instability (Im(Ω = 0)). The rest of the AM , i.e. AM i with i = 3, . . . , 7, decrease in frequency but only around q = −0.05 cross the zero frequency axis signaling the termination of this nonlinear excitation.
On the other hand, S = 2 BVB solutions destabilize via two eigenfrequency zero crossings of the two principal AMs present in the BdG spectrum of this configuration. Namely, AM 3 which is the higher-lying negative energy mode and AM 2 being the lowest-lying one. These destabilizations take place at q = 0, i.e., at the threshold (q = 0) separating the EP and the EA phases, and q = 0.15. However, among the two only the second destabilization produces a sizable imaginary component being of the order of Im(Ω) ∼ 10 −2 . Also an oscillatory instability is identified for the S = 2 BVB entity appearing at around q = 0.7. This is an instability that owes its existence to the collision of AM 3 with a positive Krein mode giving rise to an eigenfrequency quartet similar to those identified for the S = 1 structures. There exists also a third anomalous mode for this BVB configuration. Namely AM 1 , that stems from a change in sign of a background mode from positive to negative. This mode appears in the BdG spectrum for q = 0.6 and remains present as q is further lowered towards the phase transition point.
The above-discussed modes are illustrated in Fig. 10(a 2 ) while their activation leads to deformations of the stationary S = 2 BVB state, an effect that is measured via ∆Ψ mF (x, y) shown in Fig. 10(b 2 ) − (j 2 ). Notice that ∆Ψ mF (x, y) is finite irrespectively of which mode, i.e. AM 1 , AM 2 and AM 3 , is activated. Particularly, for the first mode at hand, ∆Ψ mF (x, y) acquires a 4f xz 2 orbital-like distribution as the one found for the perturbed via AM 5 VBV entity. Here though, ∆Ψ −1 (x, y) = ∆Ψ +1 (x, y) while both are complementary to the m F = 0 vortex component. Likewise, the density difference assumes a 2p orbital-like structure once AM 2 is taken into account, a result similar to the one found for the VBV solutions when AM 3 was triggered.
Finally, the eigenvector related to AM 3 is responsible for a 3d xz deformation imprinted in ∆Ψ mF (x, y) like the one found for the VBV structure when perturbed with the eigenvector associated with AM 7 [see here Fig. 10(h 2 )-(j 2 ). However here, ∆Ψ −1 (x, y) = ∆Ψ +1 (x, y) while both are complementary to ∆Ψ 0 (x, y) that is also finite in this case.
A case example showcasing the dynamical evolution of a perturbed S = 2 configuration is provided in Fig. 10(k 1 ) − (m 3 ) for q = 0.6. Notice the structural deformation of the ensuing BVB structure caused by the addition of the eigenvector related to AM 1 . Evidently, already at t = 0s a triangular pattern [54,57,58], breaking the radial symmetry of the trap along the azimuthal direction, is seen in Fig. 10(k 1 ) − (k 3 ) whose precessional motion is then followed for times up to t = 1.0s [ Fig. 10(m 1 ) − (m 3 )]. An outcome verifying that indeed, this deformation is caused by the above-identified azimuthal mode with triangular symmetry (i.e., an e 3iθ perturbation mode). It is also worthwhile to mention that similar findings are also present for AF spinor BECs (not shown).