# Emergence of unstable avoided crossing in the collective excitations of
spin-1 spin-orbit coupled Bose-Einstein condensates

Sanu Kumar Gangwar  Department of Physics, Indian Institute of Technology,
Guwahati 781039, Assam, India  Rajamanickam Ravisankar  Corresponding author:
<EMAIL_ADDRESS>Department of Physics, Zhejiang Normal University,
Jinhua 321004, China  Zhejiang Institute of Photoelectronics & Zhejiang
Institute for Advanced Light Source, Zhejiang Normal University, Jinhua,
Zhejiang 321004, China  Henrique Fabrelli  Centro Brasileiro de Pesquisas
Físicas, 22290-180 Rio de Janeiro, RJ, Brazil  Paulsamy Muruganandam
Department of Physics, Bharathidasan University, Tiruchirappalli 620024,
Tamilnadu, India  Pankaj Kumar Mishra  Corresponding author:
<EMAIL_ADDRESS>Department of Physics, Indian Institute of
Technology, Guwahati 781039, Assam, India

###### Abstract

We present the analytical and numerical results on the collective excitation
spectrum of quasi-one-dimensional spin-orbit (SO) coupled spin-1 spinor
ferromagnetic Bose-Einstein condensates. The collective excitation spectrum,
using Bogoliubov-de-Gennes theory, reveals the existence of a diverse range of
phases in the SO and Rabi ($k_{L}-\Omega$) coupling plane. Based on the nature
of the eigenvalue of the excitation spectrum, we categorize the $k_{L}-\Omega$
plane into three distinct regions, namely I, II, and III. In region I, a
stable mode with phonon-like excitations is observed. In region IIa, single
and multi-band instabilities are noted with a gapped mode, while multi-band
instability accompanied by a gapless mode between low-lying and first excited
states is realized in region IIb, which also provides evidence of unstable
avoided crossing between low-lying and first excited modes, responsible for
the $I_{o}$ type of oscillatory non-equilibrium dynamical pattern formation.
The gap between low-lying and first-excited states increases upon increasing
the Rabi coupling while decreases upon increase of SO coupling. Using
eigenvector analysis, we confirm the presence of the spin-dipole mode in the
spin-like modes in Region II. We corroborate the nature of the collective
excitation through real-time dynamical evolution of the ground state perturbed
with the quench of the trap using the mean-field Gross-Pitaevskii model. This
analysis suggests the presence of dynamical instability leading to the
disappearance of the $0$-th component of the condensate. In Region III, mainly
encompassing $\Omega\sim 0$ and finite $k_{L}$, we observe phonon-like
excitations in both the first excited and the low-lying state. The
eigenvectors in this region reveal alternative in- and out-of-phase behaviours
of the spin components. Numerical analysis reveals the presence of a super
stripe phase for small Rabi coupling in this region, wherein the eigenvector
indicates the presence of more complicated spin-like-density mixed modes.

##  I Introduction

Ultracold atoms offer a versatile but precise platform for exploring quantum
matter in the presence of diverse synthetic fields Zhang  _et al._ (2016). The
recent experimental realization of spin-orbit (SO) coupling in spinor
ultracold atoms Jacob  _et al._ (2012); Campbell  _et al._ (2016) has sparked
great interest among researchers to explore the SO coupling-related physics in
large-spin systems, which is challenging to achieve in the electronic
materials system of condensed matter physics. Over the last few decades,
research on spinor Bose-Einstein condensates (BECs) has captured the
community’s attention due to its magnanimity in showing the interesting
physical phenomena such as vortices Kasamatsu  _et al._ (2005), turbulence
state triggered by the modulational instabilities as a result of the
counterflow motion between the binary component Fujimoto and Tsubota (2012),
quantum phase transitions Yu (2016); Martone  _et al._ (2016), modulational
instability Li  _et al._ (2017), solitons Malomed (2018), magnetized vector
solitons Peng  _et al._ (2019), and supersolid-like behaviour in quasi-2D
Adhikari (2021a).

In the field of ultracold gases, numerical simulations have become an
indispensable tool for exploring numerous complex phases and the underlying
detailed mechanisms that appear in those systems as low-lying excitation
phases. To this end, the mean-field model Gross-Pitaevskii equations-based
models have been widely used to investigate various aspects of SO-coupled
spinor BECs, such as phase separation Gautam and Adhikari (2014), phase
separation of vector solitons Adhikari (2019), stability, and dynamics of
solitons Gautam and Adhikari (2015); Adhikari (2020); Gautam and Adhikari
(2021), spin precession and the separation between the spin component owing to
the anomalous spin-dependent velocities Mardonov  _et al._ (2015, 2018, 2019),
etc. Some studies have explored the dynamics of modulation instability Mithun
and Kasamatsu (2019), preparation of stripe states Cabedo  _et al._ (2021),
and condensate flow past an obstacle Zhu  _et al._ (2020) in SO-coupled BECs.

One of the intriguing features of the SO-coupled BECs is the presence of a
modulated ground state that appears as stripe phases for high values of SO
couplings Lin  _et al._ (2011); Ravisankar  _et al._ (2021a). In recent years,
there have been several studies indicating an intimate connection between the
stripe phase and the supersolid behaviour of the BECs. For instance, Adhikari
Adhikari (2021a) reported supersolid-like states in quasi-two-dimensional
trapped SO-coupled spinor BECs. In multicomponent condensates, Goldstein et
al. Goldstein and Meystre (1997) have initiated the study of the Hartree-
Bogoliubov theory to determine the resulting quasi-particle frequency spectrum
and observed the interferences resulting from cross-coupling among the
condensate leading to the reversal of the sign of two-body interaction which
triggers the onset of spatial instabilities. There are some studies that
suggest domain formation through the dynamical stability analysis Zhang  _et
al._ (2005).

The investigation of collective excitations, which are low-lying excitations
in BECs, has played a pivotal role in comprehending the fundamental
characteristics of these quantum degenerate gases. This includes aspects such
as the stability of various ground state phases, superconductivity, and
superfluidity Pethick and Smith (2008); Pitaevskii and Stringari (2016). In
1941, Landau coined the term “excitations” to elucidate the emergence of
superfluids as a combined effect of a weakly interacting mixture of phonons
and rotons quasi-particles. Subsequently, Bogoliubov employed Landau’s
excitation spectrum tools to explicate the superfluidity in BECs Bogoliubov
(1947). Experimentally, the realization of low-lying excitations has been
achieved in dilute gases of rubidium Jin  _et al._ (1996) and sodium Mewes
_et al._ (1996) atoms. For trapless spinor condensates, the study of
collective excitations involved considering equal Rashba and Dresselhaus
couplings. This investigation revealed the presence of maxon-roton excitations
Martone  _et al._ (2012); Li  _et al._ (2013). Utilizing Bragg spectroscopy,
Khamehchi et al. provided experimental evidence of collective excitations by
realizing roton-maxon modes in SO coupled spinor systems Khamehchi  _et al._
(2014).

Following the experimental realization of the collective excitation mode for
SO coupled BECs in the laboratory, a series of theoretical and numerical works
have been undertaken in the recent past that reveal the more complex nature of
these excitations. For instance, Chen et al. Chen  _et al._ (2017)
demonstrated the collective excitation spectrum of Raman-induced SO coupled
spinor BECs confined in a quasi-1D harmonic trap. By tuning the Raman coupling
strength, they projected the presence of three distinct phases: stripe (ST),
plane-wave (PW), and zero-momentum (ZM). Additionally, they identified the
Spin dipole and breathing modes of collective excitations, indicating clear
phase boundaries. Numerical confirmation of these features was obtained
through quench dynamics. Few studies indicate the presence of the PW, ZM, and
ST phases of the condensate as quadratic Zeeman coupling is introduced in the
SO coupled spin-1 BECs Chen  _et al._ (2022); He and Lin (2023). Ozawa et al.
Ozawa  _et al._ (2013) numerically analyzed the dynamical and energetic
instabilities of quasi-one-dimensional SO-coupled BECs. Ravisankar et al.
Ravisankar  _et al._ (2021b) analytically and numerically analyzed the
stability of the spin-1/2 binary SO-coupled BECs in quasi-2D using collective
excitation spectrum and found the presence of phonon-maxon-roton mode in the
spectrum for finite SO and Rabi couplings. Katsimiga et al. Katsimiga  _et
al._ (2021) analyzed the non-linear solitary wave excitation that appear in
the form of different combinations of the dark and bright solitons for the
trapped spin-1 condensate in the harmonic trap. Subsequently different sorts
of vortex-bright type excitations in the 2D harmonically confined spin-1 BECs
were also explored Katsimiga  _et al._ (2023). Rajat et al. Rajat  _et al._
(2022) theoretically and numerically analyzed the collective excitation of SO
coupled spin-1 BECs in a cigar-shaped trap both at zero and finite
temperature, and demonstrated the presence of density and spin excitation,
which exhibits qualitatively different features at finite temperature than
those at zero temperature.

So far, the studies of collective excitation in the spinor BECs are mainly
restricted to the specified range of coupling parameters. However, a
comprehensive picture of the instabilities arising as a result of the
collective excitation has not been understood clearly. Moreover, from the
dynamics point of view, it has been demonstrated that in the presence of
finite detuning between the two spin states unstable avoided crossing between
the low-lying and first excited states appears upon perturbation. Such
physical mechanism is responsible for constructing a class of $I_{o}$
(oscillatory patterns) dynamically oscillatory non-equilibrium pattern
formation in the quantum systems Cross and Hohenberg (1993); Bernier  _et al._
(2014) which remains to be unexplored in the spin-1 coupled BECs along with
spin-orbit and linear Rabi coupling. In this paper, we attempt to address
these issues considering the ferromagnetic spin-1 BECs by scanning a wide
range of the coupling parameters regime using analytical computation of the
eigenspectrum and eigenvectors. Based on the detailed nature of the
eigenspectrum, we divide $k_{L}-\Omega$ plane into three regimes: (i) a stable
region characterized by a phonon-like mode, (ii) unstable regimes broadly
divided into two parts characterized by the presence of a gap (IIa) and
gapless (IIb) behaviour between the low-lying (ll) and first-excited (fe)
spectra, and also yields the evidence of unstable avoided crossing modes, and
(iii) for coupling values near $\Omega\sim 0$, we find that the system
exhibits a phonon-like symmetric maxon-roton mode, which is also a gapless
mode between the low-lying and first-excited spectra. In the presence of a
relatively small imaginary frequency, we observe both a weak instability and
the emergence of a superstripe phase. In all the regions, we corroborate the
analytical observation of the eigenspectrum with the real-time dynamics of the
ground state using the GP equations.

The paper is organized as follows. In Section II, we present our theoretical
model for investigating the collective excitations and instabilities of SO-
coupled spin-1 ferromagnetic BECs. Section III illustrates the single-particle
spectrum, followed by the analysis of the collective excitation spectrum using
the Bogoliubov-de-Gennes theory in Section IV. In Section V, we analytically
and numerically report the collective excitation spectrum, identifying various
stability regions. Each subsection within Section V provides the excitation
spectrum and corresponding dynamical observations using the GP equations.
Finally, Section VI summarizes the findings.

##  II Mean-Field Model

We consider a quasi-1D spin-$1$ spinor SO coupled BECs realized by tight
confinement in the transverse direction Salasnich  _et al._ (2002). The non-
dimensional dynamical equations for quasi-1D SO coupled spin-$1$ BECs are
given by Gautam and Adhikari (2015); Katsimiga  _et al._ (2021); Adhikari
(2021b); Kawaguchi and Ueda (2012),

$\displaystyle\mathrm{i}\frac{\partial\psi_{\pm 1}}{\partial t}$
$\displaystyle=\bigg{[}-\frac{1}{2}\frac{\partial^{2}}{\partial
x^{2}}+V+c_{0}\rho\bigg{]}\psi_{\pm
1}\mp\frac{k_{L}}{\sqrt{2}}\frac{\partial\psi_{0}}{\partial x}$
$\displaystyle+c_{2}^{\pm}\psi_{\pm 1}+\psi_{0}^{2}\psi_{\mp
1}^{*}+\frac{\Omega}{\sqrt{2}}\psi_{0},$ (1a)
$\displaystyle\mathrm{i}\frac{\partial\psi_{0}}{\partial t}$
$\displaystyle=\bigg{[}-\frac{1}{2}\frac{\partial^{2}}{\partial
x^{2}}+V+c_{0}\rho\bigg{]}\psi_{0}+\frac{k_{L}}{\sqrt{2}}\bigg{[}\frac{\partial\psi_{+1}}{\partial
x}-\frac{\partial\psi_{-1}}{\partial x}\bigg{]}$
$\displaystyle+c_{2}^{0}\psi_{0}+2\psi_{0}^{*}\psi_{+1}\psi_{-1}+\frac{\Omega}{\sqrt{2}}(\psi_{1}+\psi_{-1})$
(1b)

where, $c_{2}^{\pm}=\left(\rho_{\pm 1}+\rho_{0}-\rho_{\mp 1}\right)$,
$c_{2}^{0}=\left(\rho_{+1}+\rho_{-1}\right)$, $\psi_{j}$; $j=+1,0,-1$, are the
spinor condensate wavefunctions that satisfy the normalization condition
$\int_{-\infty}^{\infty}\rho dx=1$, where $\rho_{j}=|\psi_{j}|^{2}$, and
$\rho=|\psi_{+1}|^{2}+|\psi_{0}|^{2}+|\psi_{-1}|^{2}$. The dimensionless
equations (1a)-(1b) are obtained by non-dimensionalized the time, length, and
energy with respect to $\omega^{-1}$, $l_{0}=\sqrt{\hbar/m\omega}$, and
$\hslash\omega$ respectively, where $\omega=\omega_{x}$ is the trap frequency
along the $x$-axis. The resulting condensate wavefunction takes the form as
$\psi_{j}=\sqrt{\frac{l_{0}}{N}}\tilde{\psi}_{j}$, and $N$ is the total number
of atoms. $V(x)=x^{2}/2$ is the trap potential,
$c_{0}=2Nl_{0}(a_{0}+2a_{2})/3l_{\perp}^{2}$ is the density-density
interaction strength, and $c_{2}=2Nl_{0}(a_{2}-a_{0})/3l_{\perp}^{2}$, spin-
exchange interaction strength with $a_{0}$ and $a_{2}$ are the s-wave
scattering length in total spin channels $0$ and $2$, respectively . The
nature of interaction strength depends on the sign of $c_{2}$. $c_{2}<0$
represents the ferromagnetic condensate, while $c_{2}>0$ denotes the anti-
ferromagnetic condensate Kawaguchi and Ueda (2012); Stamper-Kurn and Ueda
(2013). Here, $l_{\perp}=\sqrt{\hbar/(m\omega_{\perp})}$ is the harmonic
oscillator length in the transverse direction with
$\omega_{\perp}=\sqrt{\omega_{y}\omega_{z}}$. The SO and Rabi coupling
strengths are given by $k_{L}=\tilde{k_{L}}/\omega_{x}l_{0}$,
$\Omega=\tilde{\Omega}/(\hbar\omega_{x})$, respectively. In the above
description, the quantities with tilde represent dimension-full quantities.

One important entity that characterizes the miscibility of different spin
components is the magnetization of spin-1 spinor condensates defined by,

$\displaystyle\mathcal{M}=\int_{-\infty}^{\infty}\bigg{\\{}{\rho}_{+1}({x})-\rho_{-1}({x})\bigg{\\}}d{x}$
(2)

The energy functional corresponding to the coupled GP Eqs.(1a)-(1b) is given
by Gautam and Adhikari (2014); Ravisankar  _et al._ (2021c),

$\displaystyle E=$
$\displaystyle\frac{1}{2}\int_{-\infty}^{\infty}dx\bigg{\\{}\sum_{j}\lvert\partial_{x}\psi_{j}\rvert^{2}+2V(x)\rho+c_{0}\rho^{2}$
$\displaystyle+c_{2}[\rho_{+1}^{2}+\rho_{-1}^{2}+2(\rho_{+1}\rho_{0}+\rho_{-1}\rho_{0}-\rho_{+1}\rho_{-1}$
$\displaystyle+\psi_{-1}^{*}\psi_{0}^{2}\psi_{+1}^{*}+\psi_{-1}\psi_{0}^{*2}\psi_{+1})]+\sqrt{2}\Omega[(\psi_{+1}^{*}+\psi_{-1}^{*})\psi_{0}$
$\displaystyle+\psi_{0}^{*}(\psi_{+1}+\psi_{-1})]+\sqrt{2}k_{L}[(\psi_{-1}^{*}-\psi_{+1}^{*})\partial_{x}\psi_{0}$
$\displaystyle+\psi_{0}^{*}(\partial_{x}\psi_{+1}-\partial_{x}\psi_{-1})]\bigg{\\}}$
(3)

In order to make the dimensionless set of parameters used for our simulations
viable for an experiment, here, we outline the dimensionful experimental
parameters. For the ferromagnetic system, we consider 87Rb BECs with $N\sim
2\times 10^{4}$ number of atoms. The axial trap frequency is
$\omega_{x}=2\pi\times 50$ Hz while the transverse trap frequencies are
considered to be $\omega_{y}=\omega_{z}=2\pi\times 500$ Hz. The spin-dependent
and spin-independent interactions can be achieved by controlling the $s$-wave
scattering lengths through Feshbach resonance; by varying the magnetic field,
we can tune the $s$-wave scattering lengths Inouye  _et al._ (1998); Marte
_et al._ (2002); Chin  _et al._ (2010). We can achieve the SO coupling
strengths in the range of $k_{L}=\\{0.1-5\\}$ by changing the laser
wavelengths from $\\{68.86\mu$m - $1377.22\text{nm}\\}$. Further, one can also
control the Rabi coupling strength in a similar range by changing the
frequency of the laser, $2\pi\hbar\times\\{5-250\\}$ Hz.

##  III Single-particle spectrum

In this section, we present the single-particle spectrum Wen  _et al._ (2012);
Peng  _et al._ (2019); Rajaswathi  _et al._ (2023) of the non-interacting
spinor condensate for trapless SO and Rabi-coupled spin-1 spinor BECs.
Following that, we introduce the collective excitation spectrum for the
interacting systems.

For $V(x)=0$, and $c_{0}=c_{2}=0$ as we substitute $\psi_{0,\pm 1}=\phi_{0,\pm
1}e^{\mathrm{i}({q_{x}x-\omega t})}$ in the Eqs.(1a)-(1b)) we obtain

$\displaystyle\mathcal{L}_{sp}\begin{pmatrix}\phi_{1}\\\ \phi_{0}\\\
\phi_{-1}\end{pmatrix}$ $\displaystyle=\omega\begin{pmatrix}\phi_{1}\\\
\phi_{0}\\\ \phi_{-1}\end{pmatrix}$ (4)

with

$\displaystyle\mathcal{L}_{sp}=\frac{1}{2}\begin{pmatrix}q_{x}^{2}&L&0\\\
S&q_{x}^{2}&L\\\ 0&S&q_{x}^{2}\\\ \end{pmatrix},$ (5)

where
$L=-\sqrt{2}(\mathrm{i}k_{L}q_{x}-\Omega),S=\sqrt{2}(\mathrm{i}k_{L}q_{x}+\Omega)$.
Upon diagonalizing the Eq.(4), we obtain the single-particle spectrum as

$\displaystyle\omega_{0}=$ $\displaystyle\frac{q_{x}^{2}}{2}$ (6a)
$\displaystyle\omega_{\pm}=$ $\displaystyle\frac{1}{2}(q_{x}^{2}\pm
2\sqrt{\Omega^{2}+k_{L}^{2}q_{x}^{2}})$ (6b)

From Eq. 6(a) and Eq. 6(b), it is quite evident that the single-particle
spectrum exhibits three branches, namely, $\omega_{0}$ and $\omega_{\pm}$.
Among them, $\omega_{0}$ is independent of SO and Rabi couplings, a constant
branch of the single-particle spectrum. However, the detailed nature of the
other two branches $\omega_{\pm}$ depends on the SO and Rabi couplings, which
have been designated as positive and negative branches of the spectrum,
respectively, throughout the paper.

Figure 1: Single-particle energy spectrum in momentum space for the different
set of SO and Rabi coupling strengths $(k_{L},\Omega)$: (a)-(f): $(0,0)$,
$(0.7,0)$, $(0,0.5)$, $(0.7,0.5)$, $(0.7,0.7)$, $(1.0,0.5)$. The thick-solid-
red, dash-dotted green line and thin-solid-blue line represent the energy
eigenspectrum for $\\{-1,0,+1\\}$ components of the spin, respectively.
Varying the coupling parameters lead to the changes in the eigenspectrum of
different components, where $\omega_{-}$, only having the lowest minimum for
$\Omega>k_{L}^{2}$, whereas $\Omega<k_{L}^{2}$, the $\omega_{-}$ exhibits
double minima appearance of plane-wave phase, indicating the presence of
stripe phase.

To start with, we consider both the SO ($k_{L}$) and the Rabi ($\Omega$)
coupling strengths as zero. We obtain the non-degenerate parabolic spectrum of
the SO-coupled spin-1 system as couplings are set to zero [see Fig. 1(a)]. As
$\omega_{0}$ of the spectrum is independent of both the SO and Rabi couplings,
it does not show any change upon the variation of coupling strengths as shown
in Figs. 1(a)- 1(e) with the dashed green line. To understand the behaviour of
the positive branch ($\omega_{+}$) and negative branch ($\omega_{-}$) of the
single-particle spectrum, we introduce the finite value of SO coupling in the
absence of Rabi coupling and vice-versa. First, we introduce a finite value to
SO coupling strength, $k_{L}=0.7$, in the absence of Rabi coupling strength,
$\Omega=0.0$. For this parameter, the negative branch of the spectrum gets
transformed from the parabolic to double minima, which appears at $q_{x}=\pm
k_{L}$ as shown in Fig. 1(b). Here the energy minima of $\omega_{-}=-0.245$
occurs at $q_{x}=\pm 0.7$. As Rabi coupling is increased to a finite value,
$\Omega=0.5$ while keeping $k_{L}=0.0$, we find a change from double minima to
a single minimum ($\Omega>k_{L}^{2}$) and also leading to opening up a gap
between the $\omega_{+}-\omega_{0}$ branches. The magnitude of the gap between
the $\omega_{+}-\omega_{0}$ branches, and $\omega_{0}-\omega_{-}$ branches, is
of the order of $\Omega$. It appears to be of $2\Omega$ between the
$\omega_{+}-\omega_{-}$ branches as shown in Fig. 1(c) which exhibits a single
minimum with $\omega_{-}=-0.5$ at $q_{x}=0$. Further, we analyze the spectrum
by introducing both couplings strengths,i.e., $k_{L}=0.7$, $\Omega=0.5$ for
which we obtain the global minima in the spectrum, also satisfies the
condition, $\Omega\approx k_{L}^{2}$. Upon increasing the Rabi coupling value
from $\Omega=0.5\to\Omega=0.7$, we obtain a similar kind of spectrum, but the
gap between the branches increases with an increase in Rabi coupling as shown
in Fig. 1(d) and Fig. 1(e), with energies $\omega_{-}=-0.5$, $\omega_{-}=-0.7$
at $q_{x}\approx 0.0$. On the other hand, upon increasing the SO coupling from
$k_{L}=0.7$ to $k_{L}=1.0$, we obtain the transition of single minimum nature
of $\omega_{-}$ into a double minimum as shown in Fig. 1(f). At a later stage,
we will see that this feature (occurrence of double minima) in the single
particle spectrum is responsible for the appearance of the stripe phase for
the situation when $\Omega<k_{L}^{2}$. In Fig. 1(f), the lowest energy minima
is $\omega_{-}=-0.625$ at $q_{x}\approx\pm 0.9$. From the observation we find
that, when $\Omega>k_{L}^{2}$, the energy depends on only the Rabi coupling,
while if $\Omega<k_{L}^{2}$ energy depends on SO coupling strength. Overall,
we find that, upon variation of SO coupling, the positive branch exhibits a
single minimum upon fixing the Rabi coupling at a finite value. However, the
negative branch exhibits a transition from the single minimum to a double
minimum upon varying the SO coupling, keeping the Rabi coupling fixed. This
typical feature of the negative branch of the spectrum appears to indicate the
change in the phases of the condensate upon varying the coupling parameters,
and therefore, for later parts of our discussion, we will focus on analyzing
the typical features of the $\omega_{-}$ branch.

In the next section, we discuss the effect of both coupling terms (SO and
Rabi) on the collective excitation spectrum for the interacting SO and Rabi-
coupled spin-1 BECs.

##  IV Collective Excitation spectrum: Bogoliubov-de-Gennes Analysis

In the previous section, we have discussed the single-particle spectrum of SO-
coupled spin-1 ferromagnetic BECs. In this section, we shall present the
collective excitation spectrum of the condensate using the Bogoliubov-de-
Gennes (BdG) analysis. The excitation wave function in terms of perturbation
term $\delta\phi_{j}$ and ground state wave function $\phi_{j}$ can be
represented as Zhu  _et al._ (2020); Goldstein and Meystre (1997),

$\displaystyle\psi_{j}(x,t)=e^{-i\mu_{j}t}[\phi_{j}+\delta\phi_{j}(x,t)]$ (7)

where,

$\displaystyle\delta\phi_{j}(x,t)=u_{j}e^{i(qx-\omega
t)}+v_{j}^{*}e^{-i(qx-\omega^{*}t)}$ (8)

represents the ground state wave function with
$\phi_{j}=(1/2,-1/\sqrt{2},1/2)^{T}$. Here, $\mu_{j}$ denotes the chemical
potential, $u_{j}$, $v_{j}$ are the Bogoliubov amplitudes, with $j=$ $+1$, $0$
$-1$ corresponding to the three component of spinor density of condensates. As
we substitute Eq. (7) in the dynamical equations (Eqs. (1a)-(1b)) we obtain,

$\displaystyle\mathcal{L}\begin{pmatrix}u_{+1}&v_{+1}&u_{0}&v_{0}&u_{-1}&v_{-1}\\\
\end{pmatrix}^{T}$
$\displaystyle=\omega\begin{pmatrix}u_{+1}&v_{+1}&u_{0}&v_{0}&u_{-1}&v_{-1}\\\
\end{pmatrix}^{T}$ (9)

where $T$ represents the transpose of the matrix and $\mathcal{L}$ is $6\times
6$ matrix given by,

$\displaystyle\mathcal{L}=\begin{pmatrix}H_{+}-\mu_{+}&\mathcal{L}_{12}&\mathcal{L}_{13}&\mathcal{L}_{14}&\mathcal{L}_{15}&\mathcal{L}_{16}\\\
\mathcal{L}_{21}&-H_{+}+\mu_{+}&\mathcal{L}_{23}&\mathcal{L}_{24}&\mathcal{L}_{25}&\mathcal{L}_{26}\\\
\mathcal{L}_{31}&\mathcal{L}_{32}&H_{0}-\mu_{0}&\mathcal{L}_{34}&\mathcal{L}_{35}&\mathcal{L}_{36}\\\
\mathcal{L}_{41}&\mathcal{L}_{42}&\mathcal{L}_{43}&-H_{0}+\mu_{0}&\mathcal{L}_{45}&\mathcal{L}_{46}\\\
\mathcal{L}_{51}&\mathcal{L}_{52}&\mathcal{L}_{53}&\mathcal{L}_{54}&H_{-}-\mu_{-}&\mathcal{L}_{56}\\\
\mathcal{L}_{61}&\mathcal{L}_{62}&\mathcal{L}_{63}&\mathcal{L}_{64}&\mathcal{L}_{65}&-H_{-}+\mu_{-}\end{pmatrix}$
(10)

The matrix elements of $\mathcal{L}$ are given in appendix A. Bogoliubov
coefficients follows the normalization condition,

$\displaystyle\int(\lvert u_{j}\rvert^{2}-\lvert v_{j}^{*}\rvert^{2})dx=1$
(11)

The simplified form of the BdG excitation spectrum is obtained by calculating
the determinant of the matrix $\mathcal{L}$ and equate it with zero, i.e.,
$det~{}\mathcal{L}=0$. The characteristic equation can be written as,

$\displaystyle\omega^{6}+b\omega^{4}+c\omega^{2}+d=0$ (12)

where the coefficient of $b$, $c$, and $d$ are given in the appendix A.

Figure 2: Stability phase diagram in the $k_{L}$\- $\Omega$ plane for the
interaction parameter $c_{0}=0.5$, $c_{2}=-0.1$ which is being considered for
all the simulation runs. Based on the different characteristics of the
eigenspectrum and ground state, the phase diagram is divided into regions I,
IIa, IIb, and III. While Region I is stable, regions IIa, IIb, and III are
unstable. The dotted white line along with green dots represent the
$\Omega=k_{L}^{2}$ curv that separates region I from IIa. Blue dots with a
white solid line representing $\Omega\approx 0.1365k_{L}^{2}-0.0686$ separates
the regions IIa and IIb. The horizontal line with a filled half square for
$\Omega\sim 0$ denotes the region III after the cutoff, $k_{L}^{c}=0.68$,
which is a tricritical point for regions IIa-IIb-III, indicated as a red
triangle.

##  V Effect of Rabi and SO coupling on the excitation spectrum

To better understand the stability of the ground state phases of the system,
we plot the stability diagram of the system in $k_{L}$ \- $\Omega$ plane with
interaction parameters $c_{0}=0.5$, and $c_{2}=-0.1$ Wang  _et al._ (2010). We
obtain the stability phase plot by solving the collective excitation spectrum.
The real frequencies of the BdG matrix imply the dynamically stable phase
while the imaginary frequency indicates a dynamically unstable phase Kawaguchi
and Ueda (2012); Goldstein and Meystre (1997); Ravisankar  _et al._ (2021b);
Ozawa  _et al._ (2013); Zhu  _et al._ (2012). In Fig. 2, we show the stability
diagram in the $k_{L}-\Omega$ plane. Based upon the eigenvalue, the plane is
divided mainly into three regions, viz., stable region I, unstable regions II,
and III. White dots with the dotted line separate the stable and unstable
regions, which hold the phase transition boundary, i.e., $\Omega=k_{L}^{2}$
Ravisankar  _et al._ (2021b). Further, we divide the unstable region into two
separate parts: regions IIa and IIb. Blue dots with white solid line separate
the regions IIa and IIb. The phase transition critical points for region-IIa
to region-IIb can be obtained by a polynomial fit given by $\Omega\approx
0.1365k_{L}^{2}-0.0686$.

Figure 3: Excitation spectrum (a) and eigenvectors (b) and (c) for the
coupling parameters $k_{L}=0.5$, $\Omega=2.0$. (a) solid magenta line
represents $\mathrm{Re}(\omega_{-})$, and dashed-dotted green line represents
the $|{\mathrm{Im}(\omega_{-})}|$, here solid and dash-dotted lines are the
analytical results of BdG Eq. 12 and open circles are numerical results
obtained by solving Eq. 9. In the eigenvectors, $|u_{+1}|^{2}$(red open
circles), $|u_{-1}|^{2}$ (black diamonds), $|u_{0}|^{2}$ (magenta open
hexagons), $|v_{0}|^{2}$ (cyan open stars ), $|v_{+1}|^{2}$ (green open
squares), and $|v_{-1}|^{2}$ (magenta triangles), (b-c) the eigenvectors
corresponding to eigenspectrum has shown, which we obtained by solving Eq. 9
numerically. The first column of the figure shows the presence of phonon mode
in the eigenspectrum. Also, it has no negative or complex eigenfrequency thus,
it is energetically and dynamically stable.

In region IIa, we obtain a clear gap between the low-lying (negative branch)
and first-excited (positive branch) branches of the spectrum. Also, for a
fixed value of Rabi coupling strength ($\Omega$), upon increasing the SO
coupling strength ($k_{L}$), the number of instability bands changes from one
to two or vice versa. In region IIb, we obtain overlap among the low-lying and
first-excited branches of the eigenspectrum. Thus, the line delineates regions
IIa and IIb, encompassing points where instability is present by means of
multiple instability bands and gap openings between low and first excited
states. Furthermore, in region IIb, we observe a spectrum characterized by
multiple instability bands. In this region for a fixed value, Rabi strength
($\Omega$), upon increasing the SO coupling strength $k_{L}$, we achieve
gapless mode between the low-lying and first-excited branches of the spectrum
for some range of quasi-momentum. By increasing the Rabi coupling strength
with fixed SO coupling, we can again achieve the gapped mode between the low-
lying and first excited branches.

The transition from the region-IIa to IIb can be obtained after a certain
critical coupling point. We obtain the gapless behaviour of eigenspectrum
either increasing SO coupling ($k_{L}>k_{L}^{c}$) for a fixed Rabi coupling
($\Omega>\Omega^{c}$) or decreasing Rabi coupling ($\Omega^{c}<\Omega$) for a
fixed SO coupling ($k_{L}>k_{L}^{c}$), with the critical coupling parameters
as $k_{L}^{c}=0.68$, and $\Omega^{c}=0.001$. This point is the origin of the
line which separates regions IIa and IIb. The horizontal line at which Rabi
coupling strength is approximately zero, ($\Omega\approx 0.0$), we realize the
presence of phonon mode in the low-lying as well as in the first-excited
branch of the spectrum, along with the presence of instability band. These
features present in IIa and IIb are completely different compared to those
obtained for the region I. The cutoff value of SO coupling strength to obtain
this behaviour along the line is $k_{L}^{c}=0.68$. We analyze the behaviour of
the system along the horizontal line separately under Region III. Here,
$k_{L}^{c}$ is the tricritical point for regions IIa, IIb, and III, which is
indicated by the red triangle in the phase plot [see Fig. 2]. Based on the
above observations, we divide the phase plot into three different regions:
regions I, IIa, IIb, and III. In Appendix C, we show a similar kind of phase
diagram for two other different sets of interaction strengths $c_{0}=5.0$,
$c_{2}=-0.1$, and $c_{0}=885.72$, $c_{2}=4.09$. With these interactions, we
also obtain similar kinds of phase plots [see Appendix. B] as those in Fig. 2.

As we divide the phase diagram of the $k_{L}-\Omega$ plane into three regions,
we discuss their behaviour in each region individually. The detailed structure
of the presentation in different regions is as follows. We perform collective
excitation spectrum calculations to understand the stability of the system,
both dynamically and energetically. As the real but negative eigenfrequencies
of the BdG matrix suggests energetically unstable nature of the condensate
while complex eigenfrequencies indicates dynamically unstable nature of the
condensate Kawaguchi and Ueda (2012); Ozawa  _et al._ (2013); Ravisankar  _et
al._ (2021c). To validate these claims, we perform numerical simulations using
the GP equations. We provide the eigenvalue spectrum, eigenvectors, and
numerical simulation in every subsection.

Following this, we corroborate the analytical results for the excitation
spectrum by numerically solving the BdG equations (10) from which we also
obtain the eigenvectors as a function of $q_{x}$. First, we consider a
$[-1000:1000]$ grid in real space with step size $h_{x}=0.2$. Then, we use the
Fourier collocation method, where we numerically compute the Fourier transform
of the BdG equations and obtain a truncated reduced BdG matrix, which is
subsequently diagonalized using the LAPACK package Anderson  _et al._ (1999).
In momentum space, we consider $[-700:700]$ modes in the $q_{x}$ direction,
with a grid step size of $h_{q_{x}}=0.0157$.

Further we provide the real-time dynamics analysis corresponding to the
collective excitation spectrum by solving coupled GP equations (Eqs. (1a) -
(1b)) numerically. We use the imaginary-time propagation method to obtain the
ground state. Then, we evolve the ground state in time using the real-time
propagation method. For this purpose, we adopt the split-step Crank-Nicholson
scheme outlined in Refs. Muruganandam and Adhikari (2009); Young-S  _et al._
(2016); Ravisankar  _et al._ (2021c). We have considered a grid of $1280$
space points with space step $dx=0.05$, time step $dt=0.00025$ for imaginary-
time-propagation, and $dt=0.0005$ for real-time-propagation. Initially, we
obtain the ground states according to the regions in the phase plot [see Fig.
2] by using imaginary-time propagation. After calculating the ground state, we
evolve it using real-time propagation by quenching the trap strength.

Figure 4: Time evolution of ground state density profiles illustrating
dynamical stability of (a) $+1$, (b) $0$, and (c) $-1$ components of the
condensate, where the parameters are similar to Fig. 3. All three components
$|\psi_{+1}|^{2}$, $|\psi_{0}|^{2}$, and $|\psi_{-1}|^{2}$ depict stable
breather-like dynamics, which validate the dynamical stability of the
condensate.

###  V.1 Excitation spectrum in region I

#### Excitation spectrum:

We present the collective excitation spectrum for the stable case and
corresponding eigenvectors in Fig. 3. Here, we consider the coupling
strengths, $k_{L}=0.5$ and $\Omega=2.0$. Fig. 3(a) shows only the presence of
$\mathrm{Re}(\omega)$ in the eigenspectrum. The above suggests the presence of
phonon mode Pethick and Smith (2008). The absence of negative and imaginary
eigenfrequency in the spectrum indicates that the phase is energetically and
dynamically stable throughout the region-I.

We present the eigenvectors corresponding to the eigenspectrum in Fig.3 (b),
(c). We have three sets of eigenvectors: $\lvert u_{+1}\rvert^{2}$, $\lvert
u_{-1}\rvert^{2}$, $\lvert u_{0}\rvert^{2}$, $\lvert v_{0}\rvert^{2}$, $\lvert
v_{+1}\rvert^{2}$, and $\lvert v_{-1}\rvert^{2}$. The eigenvectors exhibit in-
phase behaviour in all components for all wavenumber, which means that $\lvert
u_{+1}(q_{x})\rvert^{2}=\lvert u_{-1}(q_{x})\rvert^{2}$ and $\lvert
v_{+1}(q_{x})\rvert^{2}=\lvert v_{-1}(q_{x})\rvert^{2}$, indicating the
presence of a density-like (in-phase) excitations. [see Figs. 3(b) and 3(c)].
We demonstrate that both $u$ and $v$ eigenvectors are in-phase. For
$q_{x}\approx 0$, the eigenvector components approach each other, and at
$q_{x}=0$, they have equal values, signifying the presence of a phonon mode,
as mentioned earlier. As anticipated, the eigenspectrum and eigenvectors
exhibit the same behaviour for other sets of parameters ($k_{L}$, $\Omega$) in
the same region (region I of Fig. 2) of the $k_{L}$-$\Omega$ plane. The
critical change in the eigenspectrum is that with increasing coupling
strengths, the real eigenspectrum widens similarly to how the eigenvectors
also widen after the junction at $q_{x}=0$.

Figure 5: Energy variation in the time evolution for the stable region I. The
parameter values considered in this figure are similar to Fig. 3. During time
evolution, the energy of the condensate decreases in the beginning and further
shows stable behaviour, which signifies that the condensate is energetically
stable due to the lack of negative eigenfrequency.

#### Dynamical stability:

We have considered the same coupling parameters, $(k_{L},\Omega)=(0.5,2)$, to
examine the dynamical behaviour of region-I, and the interaction strengths are
$c_{0}=0.5$, and $c_{2}=-0.1$ by solving the GP equations (Eqs. (1a) - (1b)).
Here, we obtain the plane-wave phase as a ground state of the condensate. To
understand the stability of the ground state, quench the trap strength to half
and further evolve the condensate. The density of the condensate shows stable
breather-like motion throughout the dynamical evolution [See Fig. 4(a)-(c)].
In addition, the system density profile is symmetric in nature and doesn’t
show any oscillations, as shown in Figs. 4(a)-(c). Thus, it is clear that the
system is not showing polarization behaviour as we expected from the in-phase
nature of eigenvectors shown in Figs. 3(b)-(c).

In Fig. 5, we plot the total energy of the condensate. In real-time evolution,
we quench the trap strength to one-half its original value. Hence, in the
beginning, the energy of the condensate starts decreasing. After a while, the
energy gets stabilized and constant throughout the dynamical evolution.
Therefore, we find that the ground state of the condensate is energetically
and dynamically stable, as shown in the stability phase plot Fig. 2
corresponding to the excitation spectrum. We also investigate the dynamical
spin texture of this region, which does not change its behaviour during the
dynamical evolution.

Figure 6:  Stability excitation spectrum eigenvalues and eigenvectors: The SO
and Rabi coupling strength are $(k_{L},\Omega)$ $=(2.0,2.0)$ (top row), and
$(2.35,4.0)$ (bottom row), respectively. (a1), (a2) shows eigenenergy, (b1),
(c1), and (b2), (c2 ) corresponding eigenvectors, figure labels are as like
Fig. 3. The eigenenergy spectrum shows the presence of complex
eigenfrequencies in terms of two (one) bands in the top (bottom) row, which
indicates the system is dynamically unstable. Whereas the eigenvectors display
the spin-dipole along with spin-like mode along the $q_{x}$ momentum
direction, also symmetric along the axis.

###  V.2 Excitation spectrum in region IIa

#### Excitation spectrum:

To investigate region II, we divide it into two sub-parts, region IIa and IIb.
In the first sub-part, we consider two points with fixed SO and Rabi coupling
strengths: $(k_{L},\Omega)=(2,2)$, and $(k_{L},\Omega)=(2.35,4)$,
respectively. We plot the dispersion relation of collective excitations and
corresponding eigenvectors at these values of the coupling parameters in Fig.
6. Fig. 6(a1) shows the presence of the imaginary part of the eigenspectrum at
larger values of quasi-momentum $q_{x}$ having two imaginary bands. As we
increase the value of Rabi coupling to larger extents than SO coupling
strength, the double imaginary eigen-bands change into a single band [see Fig.
6(a2)]. It is also shown in Fig. 2, where the second point falls on the
boundary of regions I and IIa. Therefore, we can state that the second point
is approaching towards the stable regime.

The eigenvector corresponding to the eigenenergy is presented in the
subsequent column of Fig. 6. Unlike the behaviour observed in Region I, we
find complex behaviour in the eigenvector components here. The eigenvectors
shows in momentum direction the transition from density-like (in-phase)
excitations, where $\lvert u_{+1}(q_{x})\rvert^{2}=\lvert
u_{-1}(q_{x})\rvert^{2}$ and $\lvert v_{+1}(q_{x})\rvert^{2}=\lvert
v_{-1}(q_{x})\rvert^{2}$, to the spin-like (out-of-phase) excitations, where
$\lvert u_{+1}(q_{x})\rvert^{2}=\lvert v_{-1}(q_{x})\rvert^{2}$ and $\lvert
v_{+1}(q_{x})\rvert^{2}=\lvert u_{-1}(q_{x})\rvert^{2}$. In Fig. 6(a1), the
eigenspectrum has double band instability; corresponding to this, dual
transitions from density-like mode to the spin-like mode for $\pm 1$
components of eigenvectors. Still, the $0$-th component exists only in
density-like mode. Also, there is an amplitude difference among all three
components of the eigenvector; the $\pm 1$ components have half of the
amplitude value than the $0$-th component [see Figs. 6(b1) and 6(c1)].

Figure 7: Dynamics: (a-c) and (d-f) depict temporal evolution of ground state
density profile of $|\psi_{+1}|^{2}$, $|\psi_{0}|^{2}$, and $|\psi_{-1}|^{2}$
for $(k_{L},\Omega)$: top row $(2.0,2.0)$ and bottom row $(2.35,4.0)$,
respectively. The density profile holds stripe wave behavior in both cases for
$t<150$, and $t<180$ for the top and bottom row, respectively. Further, it
gets fragmented into several small domains. This behavior signifies the change
in shape and amplitude of the density profile, where the bottom row $0$-th
component gets diminished, whereas the $\pm 1$ components gained the density.

For the second point in the region [see Fig. 6(a2)], the eigenspectrum has
only single-band instability. And there is a constant transition from the
density-like mode to the spin-like mode of $\pm 1$ components of eigenvectors
[see Fig. 6(b2)]. The amplitude of $\pm 1$ and $0$ appears the same as for the
previous point. It is observed from the above discussion whenever the
$\mathrm{Im}(\omega_{-})$ appears, the density-like mode (where we have
$\mathrm{Re}(\omega_{-})=0$) changes into the spin-like mode. From the
transition point of view, when comparing these two points in this region-IIa,
one can notice that the multi-band instability has mixed density-like and
spin-like mode behaviour in the eigenvectors. In detail, by comparing Figs.
6(b1) and 6(b2), the second only has the density-like mode in the major
portion of the momentum space. But first one has a mixed type of mode in the
same range of momentum space, where imaginary eigenvalues appeared at
$0<q_{x}<2$, it has spin-like mode when $2<q_{x}<3$ eigenvalues are real, then
it exhibits density-like mode, for $3<q_{x}<3.5$ again spin-like mode,
$q_{x}>3.5$ density-like mode (behaviour is symmetric about the momentum
axis). Such complex-mixed-mode has not been explored in the spinor SO coupled
system. We highlight that this study will be helpful for the experimental
researcher to achieve the stable ferromagnetic spinor SO coupled BEC
experimentally.

#### Dynamical stability:

We divide region II into two subregions: IIa and IIb. We focus on the system’s
dynamics in region IIa, specifically at two different points within the region
IIa. For each point, we analyzed the dynamic behaviour of the system
separately. We discuss the dynamic behaviour of the system for each point
separately. We choose $(k_{L},\Omega)=(2,2)$. We obtain the ground state using
imaginary-time propagation, which is the stripe-wave and evolve the ground
state using real-time propagation. The dynamical evolution of the state shows
the stripe wave phase holds up to $t=150$. Furthermore, the density profile of
all three components fragmented into non-periodic domains [see Figs. 7(a)-(c)]
Mithun and Kasamatsu (2019). If the system has a nonvanishing $0$-th component
in the dynamical evolution, then the system tends to be less subjected to
instability. Conversely, if the $0$-th component vanishes, the system becomes
more prone to instability Tasgal and Band (2015). In this case, the appearance
of the $0$-th component in the dynamical evolution confirms that the chosen
coupling point has a lower risk of instability. We plot the total energy of
the condensate in Fig. 8(a). The energy of the condensate starts with a value
of $-2.471$. It changes sign from negative to positive at $t=205$, due to the
quench similar kind of agitation happens, and further, it shows stable
behaviour near $t=365$. For this point, as shown in Fig. 6(a1), the BdG
excitation spectrum has complex eigenfrequencies, and the eigenvector shows
the out-of-phase behaviour for the certain range of quasi-momentum. The above
indicates that the system is dynamically unstable and generates non-linear
wave patterns, which matches well with numerical simulation results.

Figure 8: (a) and (b) represent the total energy of the condensate for the
density evolution in Fig. 7. In both cases, the energy of the condensate
increases at the beginning of dynamical evolution, but finally shows well-
settled stable behavior, such kind of sudden change happened due to the
external perturbation. This behavior validates the energetic stability of the
condensate, however, the system is dynamically unstable.

We choose another point in the region IIa at which we consider the value of SO
coupling strength $k_{L}=2.35$ and $\Omega=4.0$. We achieve the ground state,
which is stripe-wave-phase. The excitation spectrum has complex
eigenfrequencies, as shown in Fig. 6(a2). The presence of complex
eigenfrequencies indicates that the system is dynamically unstable. To verify
this claim numerically, we evolve the ground state of the system using the
real-time-propagation method as illustrated in Figs. 7(d)-(f). The stripe wave
behaviour of the condensate holds its shape for a while. Then, it gets
fragmented into several small domains Mithun and Kasamatsu (2019). In
addition, we observe that with time, the zeroth component of density starts
diminishing, and finally, it disappears [see Fig. 7(e)]. Other components of
the condensate density ($\lvert\psi_{\pm 1}\rvert^{2}$) boosted the magnitude
after the zeroth component of the density disappears [see Figs. 7(d), and
(f)]. As one considers the quadratic Zeeman effect, then such kind of
dynamical disappearance does not take place Tasgal and Band (2015). However,
in this work, we see such sort of dynamical instability in the excitation of
spin-$1$ BECs in the presence of SO coupling alone. So far we find that the
density pattern of the system changes its shape and amplitude accompanied with
disappearance of the $0$-th component, which is a signature of dynamical
instability. Further, we show the total energy of the condensate in Fig. 8(b).
It starts from $-4.182$ and changes its sign from negative to positive at
$t=1080$, which is higher than the previous point, where the $0$-th component
vanishes, and then the particles are shared by the spin-up and down
components, there, the energy gets settled. Such kind of $0$-th component
vanishing behaviour is not observed for the previous case due to the
appearance of spin and density-mixed modes. Upon investigating the dynamical
spin texture of this region, we notice that the initial pattern does not
maintain its shape as time progresses and changes its spin texture during the
dynamical evolution [see Appendix. C].

Figure 9: Excitation spectrum (a), and (b) for the coupling strengths,
$(k_{L},\Omega)=(3.1,1.17)$ and $(5.0,2.5)$ depict the low-lying and first-
excited branches of the spectrum. The solid magenta line represents the
$\mathrm{Re}(\omega_{-})$ and the dashed green line represents the
$|\mathrm{Im}(\omega_{-})|$ which is low-lying branch of the spectrum, while
the blue dash-dotted line represents $\mathrm{Re}(\omega_{+})$ and black
dotted line represents $|\mathrm{Im}(\omega_{+})|$ corresponding to the first-
excited branch of the spectrum, which obtained from Eq. 12. The eigenspectrum
in (a, b) shows multiple instability bands as we obtained in 6(a1), however,
(a) and (b) show one and two gapless modes. Insets are for the better
visualization of gapless mode, the secondary $y$-axis scale is for
representing the imaginary eigenfrequencies magnitudes. Figure 10: The
eigenvectors top and bottom rows corresponding to Fig. 9(a), and Fig. 9(b),
respectively. The eigenvectors represents as $|u_{+1}|^{2}_{ll}$ (red dot),
$|u_{-1}|^{2}_{ll}$ (black hexagons), $|u_{0}|^{2}_{ll}$ (magenta hexagons),
$|v_{0}|^{2}_{ll}$ (cyan open stars), $|v_{+1}|^{2}_{ll}$ (green open
squares), and $|v_{-1}|^{2}_{ll}$ (magenta triangles). Whereas,
$|u_{+1}|^{2}_{fe}$ (blue open pentagon), $|u_{-1}|^{2}_{fe}$ (orange stars),
$|u_{0}|^{2}_{fe}$ (olive open plus), $|v_{0}|^{2}_{fe}$ (red open x with the
dotted line), $|v_{+1}|^{2}_{fe}$ (cyan open diamond), and $|v_{-1}|^{2}_{fe}$
(maroon dots). We depict the $+q_{x}$, however, the eigenfrequencies and
eigenvectors are symmetric about the axis. Here we found that the eigenvectors
show the spin-like mode when the complex-frequency exhibits, the density-like
mode when it is only real.

###  V.3 Excitation spectrum in region IIb

#### Excitation spectrum:

In this subpart of region II, we choose two different points at which we
consider the coupling strengths as $k_{L}=3.1$, $\Omega=1.17$, and
$k_{L}=5.0$, $\Omega=2.5$, respectively. The collective excitation spectrum of
the condensate for these sets of parameters is given in Fig. 9. We obtain
complex eigen frequencies for both points [see Figs. 9(a)- 9(b)]. Still, the
number of instability bands along the quasimomentum direction and its
amplitude differ for both points. In Fig. 9(a), the number of low-lying
instability bands is three, and the positions and amplitudes are
$\\{q_{x},\omega\\}$ = $\\{0.27,0.29\\}$, $\\{1.48,0.11\\}$, and
$\\{5.75,0.11\\}$, respectively. Further, we show the first excited
($\omega_{fe}$) and low-lying ($\omega_{ll}$) branches of the eigenspectrum,
which overlap with each other, where we also find that the complex
eigenfrequency is not similar to the eigenspectrum of the previous case
region-IIa. This results in gapless modes between the first- and low-lying
excited spectrum. However, we find three complex eigenfrequencies in which the
overlapping point shows complex eigenfrequencies for both spectrums. The
phenomenon that has been observed in the spinor BECs is related to the
unstable avoided crossing between spin and charge modes by considering
$\theta=\pi$, the phase difference between the perturbed wave fields, which
also provides the class $I_{o}$ oscillatory patterns as a result of the
dynamical instability Cross and Hohenberg (1993); Bernier  _et al._ (2014).
However, in this work, we have realized the same phenomenon only for the
symmetric case $\theta=0$, even though we consider only SO and Rabi couplings
without Raman detuning. Due to the $I_{o}$ dynamical instability, we find out-
of-equilibrium dynamics in the density patterns, and it is the route to the
nonlinear patterns formation, which we shall discuss later in this subsection.
In Fig. 9(b), the number of instability bands is four, and the positions and
amplitudes are $\\{q_{x},\omega\\}$ = $\\{0.27,0.33\\}$, $\\{1.39,0.18\\}$,
$\\{3.48,0.16\\}$, and $\\{9.45,0.11\\}$ (not shown in figure), respectively.
The appearance of instability bands for both points is symmetric about the
momentum ($q_{x}$) axis. In addition to this, we show the first excited state
in the same plot. Here, the first excited and low-lying branches overlap
twice, where we observed that double-unstable avoided crossing in SO-coupled
spinor ferromagnetic BECs. For this case, the first excited spectrum also has
two complex eigenfrequencies, while the unstable-crossing avoided points.

Figure 11: Dynamics: (a)-(c) and (d)-(f) represent the dynamics of the ground
state density profile of $|\psi_{+1}|^{2}$, $|\psi_{0}|^{2}$, and
$|\psi_{-1}|^{2}$ components for $(k_{L},\Omega)$ top row $(3.1,1.17)$, and
bottom row $(5.0,2.5)$. The density profile holds stripe wave behavior in both
cases for $t<100$ and $t<50$ units for the top and bottom row, respectively.
Top row: $t>100$, it shows non-periodic oscillation and non-uniform density in
space-time, whereas the $0$-th component gets disappeared. Bottom row: $t>50$,
it further gets fragmented into several small domains. This behaviour
signifies the change in shape and amplitude of the density profile, along with
the disappearance of the $0$-th component, which confirms the dynamical
instability.

The eigenvector corresponding to the eigenspectrum is shown in the subsequent
columns of Fig. 10. For the choice of parameters $(k_{L},\Omega)=(3.1,1.17)$,
we observed no overlap between the low-lying and the first excited spectrum
eigenvectors. However, the gap between them closes, and the overlap happens
for certain values of quasimomentum $q_{x}$. The interval $q_{x}=(-0.4,0.4)$,
where the low-lying spectrum exhibits imaginary eigenvalues, for this case,
the eigenvectors correspond to the spin-like mode, as we see from the out-of-
phase behaviour between the low-lying eigenvectors. Further, when $q_{x}>\pm
0.4$ up to the gapless point, we notice that the low-lying eigenvectors
display density-like (in-phase) mode. However, the first-excited state only
have a density-like mode up to this gapless point. Nevertheless, in the region
of the spectra where the low-lying and the first excited branches overlap, we
see that the eigenvectors of both the low-lying and first-excited states are
out-of-phase among them. When the gapless mode transforms as gapped mode, all
the eigenvectors are in phase [see Fig. 10(a1)], where the $0$-th component
shows only density-like mode [see Fig. 10(a2)].

For the second point of this region, all eigenvector components have the same
amplitude value. Choosing the coupling parameters $(k_{L},\Omega)=(5,2.5)$
reveals a substantial difference. With these parameters, we observed double
overlaps between the low-lying and the first excited spectrum, closing the gap
between them. These overlaps occur for certain values of quasimomentum
$q_{x}$. The first two eigenvectors exhibit behaviour similar to that observed
at the previous point.

As the gapless mode transforms into a gapped mode, all the eigenvectors become
in phase. We also observe similar feature for the next gapless mode and
further momentum directions. Additionally, we find that when the real
eigenvalues fall to zero and then rise, the eigenvectors are flipped but
remain in phase as shown in Fig. 10(a1). The obvious change in the spectrum
from the gapped to the gapless mode of the eigenvectors demonstrates the
transition from a constant density-like mode to a hybridized density-spin-like
mode Abad and Recati (2013) [see Figs. 10(b1, b2)].

#### Dynamical stability:

We consider the points $(k_{L},\Omega)=(3.1,1.17)$, with interaction
parameters, $c_{0}=0.5$ and $c_{2}=-0.1$. The collective excitation spectrum
corresponding to this point is shown in Fig. 9(a), which has complex
eigenfrequencies in terms of multiple bands. The presence of complex
eigenfrequencies makes the system dynamically unstable. To verify that the
system is dynamically unstable, we perform numerical simulation. We obtain the
ground state, which is the stripe-wave phase. We evolve the ground state by
quenching the trap, and during dynamical evolution, the density profile holds
stripe wave behaviour initially. After a while, the $\pm 1$ components are
fragmented into small domains and show non-periodic oscillation in space-time
evolution, shown in the density profiles in Figs. 11(a) and 11(c)]. On the
other hand, the zeroth component of density starts diminishing as time
progresses [see Fig. 11(b)]. Thus, we state that the density pattern of the
condensate changes its shape and amplitude during dynamical evolution, which
confirms the dynamical instability as the $0$-th component disappears due to
the instability of the system Tasgal and Band (2015).

Figure 12: (a) and (b) depict the total energy of the condensate for the
dynamical evolution for the Fig. 11. From the energy saturation, we understood
that the system is energetically stable but dynamically unstable.

We investigate the second point by choosing $(k_{L},\Omega)=(5,2.5)$. For
this, we obtain the ground state as the stripe wave and evolve the ground
state numerically. The dynamical evolution of this phase shows the stripe wave
phase holds for a while. Further, the density profile of all three components
fragmented into non-periodic domains [see Figs. 11(a) and 11(c)] Mithun and
Kasamatsu (2019). Also, the density of the zeroth component starts
diminishing, explicitly showing that the system holds dynamical instability.
We also observe that the density pattern of the system changes its shape and
amplitude during dynamical evolution.

We plot the total energy of the condensate in the Fig. 12. The ground state
energy of the condensate starts with the values (a) $E_{0}=-4.845$, (b)
$E_{0}=-12.593$. It changes sign from negative to positive, then stable, after
the sudden change in the parameter during time evolution. We observed that the
second point is more energetic than the first point as well as in the ground
state and excited state, i.e., (a) $E_{ex}\approx 0.432$, (b) $E_{ex}\approx
1.708$. Further, we investigated the dynamical spin texture, where the initial
pattern is not observed, and at the final time, they changed their spin-
textures in the dynamical evolution [see Appendix. C].

###  V.4 Excitation spectrum in region-III

#### Excitation spectrum:

In region-III, we choose points along the horizontal axis to understand the
system’s behaviour in the absence of Rabi coupling strength ($\Omega=0$), and
we consider the SO coupling strength as $k_{L}=5.0$. The collective excitation
spectrum for this is demonstrated in Fig. 13(a, b), where differently from the
previous cases, in region III we observe the presence of roton-maxon-like
modes. As demonstrated in previous studies Chen  _et al._ (2022); Khamehchi
_et al._ (2014); Lyu and Zhang (2020) when $\Omega\rightarrow 0$ there is the
appearance of roton-maxon modes and the subsequent softening of the roton mode
when $\Omega$ approaches to zero. The softening of the roton mode favors the
appearance of superstripe structures. The low-lying excitation has complex
eigenfrequencies alongside the phonon and symmetric roton-like mode,
indicating that the system is dynamically unstable due to the presence of the
imaginary eigenfrequencies.

Figure 13: Excitation spectrum (a) low-lying and (b) first-excited branches,
for the coupling strengths $k_{L}=5.0$, $\Omega=0.0$. The label representation
is similar to the Fig. 9. Insets are maximized for the visualization of the
complex frequency mode. Whereas the gapless mode also appeared along with the
symmetric roton mode.

Further, the first-excited state also shows similar behaviour, i.e., the
presence of the phonon mode along with the complex eigenfrequencies but no
symmetric roton-like structure. Although, the unstable avoided-crossing is
noticed at $q_{x}\approx 4.8499$, between the first- and low-lying states.
After that point, the real eigenfrequency decreases and reaches zero, where it
exhibits the imaginary eigenfrequency. Beyond this point, the real frequency
gains the growth. Based on this observation, one has to understand such
behaviour using the eigenvectors. Additionally, the complex eigenfrequency
presents multiple instability bands along the quasimomentum axis. The
positions of the low-lying excitation spectrum instability bands are
$q_{x}=0.06,4.91,9.97$, and the corresponding complex eigenfrequencies are
$\omega_{-}=0.1037,0.1404,0.0993$, respectively. The first excited state has
two instability bands, and its position and complex eigenfrequency are the
same as those of the low-lying state. Notably, the instability bands are
symmetric about the $x$-axis.

Fig. 14 represents the eigenvector components, in this region, the low-lying
excitation eigenvectors exhibits in-phase mode between $|u_{+1}|_{ll}^{2}$ and
$|u_{-1}|_{ll}^{2}$ as well as between $|v_{+1}|_{ll}^{2}$ and
$|v_{-1}|_{ll}^{2}$. This feature continues up to the point of real
eigenvalues, when the imaginary eigenvalues exhibit out-of-phase between them
($|u_{+1}|_{ll}^{2}-|u_{-1}|_{ll}^{2}\neq 0$). Although at $q_{x}=0$, the
eigenvalues become zero $\omega=0$, a phonon-mode is observed, which is
explicitly shown in the form of sharp spike-like-eigenvectors appeared in
$|u_{+1}|_{ll}^{2}$ and $|v_{+1}|_{ll}^{2}$. However, in this complex
eigenvalues regime, we have both phonon-mode as well as complicated spin-like
mode. In addition to this, the $0$-th component eigenvectors show a mixed-mode
behaviour, and as a result of this, we observed different dynamical phases
from the ground state. Also, the first excited state exhibits $\omega=0$ for
$q_{x}=0$, where all the eigenvectors met at a single point as before. The
first-excited state behaves like a spin-like mode with an asymmetric nature,
indicating the existence of the superstripe phase. Further, the density mode
appears while the real eigenfrequency is present, also evident from the $0$-th
component eigenvector behaviour.

Figure 14: Eigenvectors top and bottom rows are corresponding to the Fig.
13(a) and, Fig. 13(b) respectively. The label representation is similar to the
Fig. 10. In this region, both low-lying and the first excited state assume the
value $\omega=0$ for $q_{x}=0$ with eigenvectors are in-phase with the low-
lying and first-excited states, which indicates the existence of the phonon
mode for these branches.

#### Dynamical stability:

We consider a case with SO coupling strength $k_{L}=5.0$ and zero Rabi
coupling strength $\Omega=0.0$. The corresponding BdG excitation spectrum in
Figs. 13(a) and (b) reveal the presence of phonon-gapless complex
eigenfrequencies indicating the dynamical unstable nature of the condensate.
To investigate this numerically, we obtained the ground state of the
condensate, which falls within the unpolarized stripe-wave phase. We then
evolved this state in time using a real-time propagation method. The density
profile initially maintains its unpolarized stripe-wave behaviour however
subsequently breaks into two parts at $t=85$. Following this, the density of
$\pm 1$ components oscillates around the trap centre, exhibiting accumulated
breather behaviour for a period. Ultimately, the fragmented oscillating waves
merge and form several small domains [see Figs. 15(a) and (c)] Mithun and
Kasamatsu (2019), while the zeroth component of the density disappears [see
Fig. 15(b)]. This behaviour signifies that the dynamical instability also
changes its shape and size Tasgal and Band (2015). We examined the energy
during its dynamic evolution, which eventually stabilizes, indicating
energetic stability. This finding is consistent with the analysis of the
excitation spectrum, suggesting that the numerical simulation accurately
reflects the theoretical predictions.

Figure 15: Dynamics: (a)-(c) represent the time evolution of super-stripe
ground state density profiles of $|\psi_{+1}|^{2}$, $|\psi_{0}|^{2}$, and
$|\psi_{-1}|^{2}$ components of the condensate. The coupling strengths are
$k_{L}=5.0$, $\Omega=0.0$ with interaction parameters $c_{0}=0.5$ and
$c_{2}=-0.1$. During dynamics, the unpolarized super-stripe wave holds its
behavior for $t<30$. Further, it breaks into two domains and oscillated either
sides of the condensates of $\pm 1$ components, and the $0$-th component
disappears, where $\pm 1$ components densities are obtaining the gain, which
shows the change in the density profile confirms the dynamical instability.

###  V.5 Effect of coupling strength on the Band gap

So far, we have discussed the nature of collective excitations, eigenvectors
and further the dynamical behaviour of the condensate in different regions of
the $k_{L}-\Omega$ plane, next we consider the effect of SO and Rabi coupling
strengths on the band gap $\Delta_{g}$ between the positive and negative
branches of the spectrum, where we define the band gap as $\Delta_{g}$. In
Fig. 16(a1), we show the band gap variation upon varying the Rabi coupling
strength ($\Omega$) for a fixed $k_{L}$. In what follows we discuss it for two
distinct SO strengths, for instance, $k_{L}=1.14$ and $3.10$. The gap closes
at a certain value of $\Omega$ for these SO couplings. For the point
$k_{L}=1.14$, $\Delta_{g}=0$ for $\Omega=0.107$ and similarly for
$k_{L}=3.10$, it closes at $\Omega=1.218$ for the quasi-momentum values
$q=0.521$, $1.486$, respectively. Starting from the point $\Delta_{g}=0$, upon
increasing the $\Omega$, we observe that the gap becomes non-zero even for a
relatively small increment in the Rabi coupling strength, and the magnitude of
the band gap increases with an increase in $\Omega$. Also, we provide the
logarithmic scale plot corresponding to Fig. 16(a1) wherein we show the
fitting curve $\Delta_{g}=a\Omega^{b}$ with the maroon and orange dotted line
for $k_{L}=1.14$ and $k_{L}=3.10$, respectively. The coefficient for the
fitting function is $a=0.981$, $b=0.997$ and $a=1.083$, $b=1.026$ for the SO
coupling strength $k_{L}=1.14$ and $k_{L}=3.10$, respectively. Here, we wish
to highlight that the band gap gets increased upon increasing $\Omega$ for
fixed $k_{L}$. For higher $k_{L}$, the gap attains sudden maximum even for the
lesser $\Omega$ while comparing small $k_{L}$, evident from the red dashed
line crosses the green dashed-dotted line in Fig. 16(a1). Moreover, for higher
$k_{L}$ the gap closes at higher $\Omega$, compare starting points in Fig.
16(a1).

Figure 16: Variation of band gap $\Delta_{g}$ between negative and positive
branches of the spectrum. (a1) For the fixed value SO coupling strength upon
variation of Rabi coupling strength, vice versa in (b1). In (a1), the dash-
dotted green and the dashed red lines represent gap calculation for
$k_{L}=1.14$ and $k_{L}=3.10$, respectively. In (b1), the dash-dotted blue and
the dashed black lines represent gap calculation for $\Omega=1.3$ and
$\Omega=3.1$, respectively. The interaction parameter strength are
$c_{0}=0.5$, and $c_{2}=-0.1$. (a2) and (b2) are the logarithmic scale plots
corresponding to (a1) and (b1), respectively. Further, in (a2) and (b2) the
maroon (orange) dotted lines represent the fitting curve
$\Delta_{g}=a\Omega^{b}$ with the parameters $a=0.9811,b=0.9974$,
($a=1.0825,b=1.0258$), and $a=126.2992,b=-4.9568$, ($a=46.4438,b=-2.2146$),
respectively.

Further, we discuss the variation of band gap upon varying the SO coupling for
fixed Rabi coupling strength. We choose two different values of Rabi coupling
strength $\Omega=1.3,3.1$. Since Rabi coupling strength is responsible for
opening the gap among the branches of the collective excitation spectrum (as
discussed above). As result of this the band gap $\Delta_{g}$ value is higher
for a lesser value of $k_{L}$ at fixed $\Omega$ [see Fig. 16(b1)]. Upon
increasing $k_{L}$, the $\Delta_{g}$ starts decreasing, and further, it
reaches zero at the certain quasi-momentum value $q_{x}=1.5325,2.3381$ with
$k_{L}=3.20,4.86$ for the Rabi coupling strengths $\Omega=1.3,3.1$,
respectively. We find that the band gap $\Delta_{g}$ decreases upon increasing
the $k_{L}$ at fixed $\Omega$. In Fig. 16(b2), we show the fitting function
$\Delta_{g}=c\Omega^{d}$ with the maroon and orange dotted line for
$\Omega=1.3$ and $\omega=3.1$, respectively. The coefficient of fitting
function is $c=126.30,d=-4.96$ and $c=46.44,d=-2.22$ for the SO coupling
strength $\Omega=1.3$ and $k_{L}=3.1$, respectively.

##  VI Summary and Conclusions

In this paper, we have investigated the collective excitation spectrum of SO
coupled spin-1 ferromagnetic BECs with Rabi coupling in quasi-one dimension
using Bogoliubov-de Gennes theory and numerical simulation. Firstly we have
analyzed the single-particle spectrum for different sets of SO and Rabi
coupling strengths. Upon increasing the strength of coupling parameters, we
noticed separation among different negative ($\omega_{-}$) and positive
($\omega_{+}$) branches of the eigenspectrum, with $\omega_{-}$ having the
lowest minima. However for $\Omega<k_{L}^{2}$, the negative branch of the
spectrum shows double minima indicating the presence of a stripe phases.

Based upon collective excitation calculation, we have obtained the phase
diagram in the $k_{L}-\Omega$ plane which have been broadly divided into three
regions: stable phase (region I) and unstable (region II, III). We have
divided the unstable region II into two subparts according the nature of
eigenspectrum, which is based upon the appearance of the the gapped and
gapless (unstable-avoided crossing) mode. In region I, the excitation spectrum
exhibits solely real positive eigenvalues with the presence of a phonon mode.
The eigenvectors corresponding to these eigenvalues substantiate the existence
of the phonon mode by exhibiting the density like mode in the wavenumber
space. The absence of the negative and complex eigenfrequencies in the
spectrum confirms dynamically and energetically stable condensate in this
region.

As the coupling strengths increase, we find a transition from Region I to the
unstable region II, attributed to the relation $\Omega<k_{L}^{2}$. In this
region, complex eigenfrequencies emerge throughout the spectrum, accompanied
by gapped and gapless modes between the low-lying $(\omega_{-})$ and first-
excited $(\omega_{+})$ branches. Based on the presence of gapped and gapless
modes region II is divided into two subregions, IIa and IIb, respectively. The
eigenspectrum in region IIa exhibits gapped behaviour across the entire
region. Notably, the presence of complex eigenfrequencies is reflected in the
corresponding eigenvectors, indicating a transition from spin-mode to density-
mode nature of the eigenvectors as the eigenfrequency shifts from complex to
real. This phenomenon is evident in the $\pm 1$ components of the
eigenspectrum, while the zeroth component exclusively displays a density-like
mode. The combined analysis of eigenvectors and dynamic behaviour reveals a
mixed-spin-density mode.

Region IIb is separated from Region IIa by a line, the latter of which
exhibits a multi-instability band across its entirety. Upon increasing the SO
coupling for a point along the separation line, we observe gapless behaviour
in the eigenspectrum for a particular range of quasi-momentum among the low-
lying and first-excited branches of the spectrum. Similarly, decreasing the
Rabi coupling strength for any point on the separation line will also yield
gapless behaviour. Similar to the previous region, we also identify the mixed-
spin-density eigenvector mode in this region. For the gapless region, the
eigenvector of the low-lying and first-excited branches exhibits a spin-like
mode for all its components. We also characterize the gapless mode as an
unstable avoided crossing between the low-lying and first-excited spectrum.
This instability manifests as the emergence of complex eigenfrequencies when
the two branches come into contact. The corresponding dynamical evolution also
exhibits an unstable stripe state, which we have examined numerically. We
observed that the $0$-th component is non-zero in region IIa but disappears in
region IIb. This suggests that a non-zero $0$-th component tends to be less
susceptible to dynamical instability, while a vanishing $0$-th component tends
to be higher than the previous one.

In addition to these two regions, we have identified the region III which
corresponds to the $\Omega=0$ line. In this region, both the low-lying and the
first excited states assume the value $\omega=0$ for $q_{x}=0$ with
eigenvectors in phase with the low-lying state. This indicates the existence
of a phonon mode. The low-lying and first-excited spectra also exhibit
gapless, unstable-avoided crossing behaviour between them. Further analysis of
the eigenvectors reveals not only the phonon character but also the presence
of a spin-like mode in both states. Numerical simulations for this case reveal
the presence of super stripe ground state, and their respective dynamics
exhibit unstable dynamical-nonlinear patterns. Our observations suggest that
this work would be a valuable resource for experimental researchers seeking to
achieve stable and unstable SO coupled spinor ferromagnetic BECs. Here, we
have provided a comprehensive numerical simulation of regions with diverse
dynamical behaviours.

In this work we find the absence of the PW phase in the excitation spectrum.
In the recent studies it has been shown that in presence of the Zeeman
coupling the double degeneracy in the excitation spectrum gets lifted off
hence the PW phase appears Chen  _et al._ (2022); He and Lin (2023). We would
be interested to extend the present work in presence of the Zeeman field and
explore the effect of the field on the nature of overall eigenspectrum and
eigenvectors. One of the another interesting features that have emerged in the
SO coupled BECs is the appearance of self-bound quantum droplet in the binary
BECs Cabrera  _et al._ (2018); Petrov (2015) in which the repulsive nature of
the mean-field interaction is balanced by the beyond mean-field quantum
fluctuation in one-dimension Petrov (2015); Gangwar  _et al._ (2022). It would
be interesting to extend the present work of spin-1 considering the effect of
quantum fluctuation on different phases of the collective excitation modes in
the similar line as presented for spin-1/2 SO coupled BECs Gangwar  _et al._
(2023).

###### Acknowledgements.

We gratefully acknowledge our super-computing facility Param-Ishan (IITG),
where all the simulation runs were performed. R.R. acknowledges the
postdoctoral fellowship supported by Zhejiang Normal University, China, under
Grants No. YS304023964. H. F. acknowledges the financial support from Programa
de Capacitação Institucional (PCI-CBPF) [No. 31.7202/2023-5]. P.M.
acknowledges the financial support from MoE RUSA 2.0 (Bharathidasan University
- Physical Sciences).

##  Appendix A Relevant terms of the BdG matrix of collective excitations

In this appendix, we provide an explicit form of the matrix elements of the
BdG matrix Eq. (10) used in Section V. The matrix elements of Eq. (10) read as

$\displaystyle H_{+}=$
$\displaystyle\frac{q_{x}^{2}}{2}+c_{0}(2\phi_{+1}^{2}+\phi_{0}^{2}+\phi_{-1}^{2})+c_{2}(2\phi_{+1}^{2}+\phi_{0}^{2}-\phi_{-1}^{2})$
(13a) $\displaystyle H_{0}=$
$\displaystyle\frac{q_{x}^{2}}{2}+c_{0}(\phi_{+1}^{2}+2\phi_{0}^{2}+\phi_{-1}^{2})+c_{2}(\phi_{+1}^{2}+\phi_{-1}^{2})$
(13b) $\displaystyle H_{-}=$
$\displaystyle\frac{q_{x}^{2}}{2}+c_{0}(\phi_{+1}^{2}+\phi_{0}^{2}+2\phi_{-1}^{2})+c_{2}(2\phi_{-1}^{2}+\phi_{0}^{2}-\phi_{+1}^{2})$
(13c) $\displaystyle\mu_{+}\phi_{+1}=$ $\displaystyle
c_{0}(\phi_{+1}^{2}+\phi_{0}^{2}+\phi_{-1}^{2})\phi_{+1}+c_{2}(\phi_{+1}^{2}+\phi_{0}^{2}-\phi_{-1}^{2})\phi_{+1}+c_{2}\phi_{0}^{2}\phi_{-1}^{*}+\frac{\Omega}{\sqrt{2}}\phi_{0}$
(13d) $\displaystyle\mu_{0}\phi_{0}=$ $\displaystyle
c_{0}(\phi_{+1}^{2}+\phi_{0}^{2}+\phi_{-1}^{2})\phi_{0}+c_{2}(\phi_{+1}^{2}+\phi_{-1}^{2})\phi_{0}+2c_{2}\phi_{0}^{*}\phi_{+1}\phi_{-1}+\frac{\Omega}{\sqrt{2}}(\phi_{+1}+\phi_{-1})$
(13e) $\displaystyle\mu_{-}\phi_{-1}=$ $\displaystyle
c_{0}(\phi_{+1}^{2}+\phi_{0}^{2}+\phi_{-1}^{2})\phi_{-1}+c_{2}(\phi_{-1}^{2}+\phi_{0}^{2}-\phi_{+1}^{2})\phi_{-1}+c_{2}\phi_{0}^{2}\phi_{+1}^{*}+\frac{\Omega}{\sqrt{2}}\phi_{0}$
(13f) $\displaystyle\mathcal{L}_{12}$
$\displaystyle=C^{+}\phi_{+1}^{2};\;\;\mathcal{L}_{13}=C^{+}\phi_{0}^{*}\phi_{+1}-\frac{k_{L}}{\sqrt{2}}iq_{x}+2c_{2}\phi_{0}\phi_{-1}^{*}+\frac{\Omega}{\sqrt{2}};\;\;\mathcal{L}_{14}=C^{+}\phi_{0}\phi_{+1};\;\;\mathcal{L}_{15}=C^{-}\phi_{-1}^{*}\phi_{+1}$
$\displaystyle\mathcal{L}_{16}$
$\displaystyle=C^{-}\phi_{-1}\phi_{+1}+c_{2}\phi_{0}^{2};\;\;\mathcal{L}_{21}=-C^{+}\phi_{+1}^{*2};\;\;\mathcal{L}_{23}=-C^{+}\phi_{0}^{*}\phi_{+1}^{*};\;\;\mathcal{L}_{24}=-C^{+}\phi_{0}\phi_{+1}^{*}+\frac{k_{L}}{\sqrt{2}}iq_{x}-2c_{2}\phi_{0}^{*}\phi_{-1}-\frac{\Omega}{\sqrt{2}};$
$\displaystyle\mathcal{L}_{25}$
$\displaystyle=-C^{-}\phi_{-1}^{*}\phi_{+1}{*}-c_{2}\phi_{0}^{*2};;\;\
\mathcal{L}_{26}=-C^{-}\phi_{-1}\phi_{+1}^{*};\;\;\mathcal{L}_{31}=C^{+}\phi_{+1}^{*}\phi_{0}+2c_{2}\phi_{0}^{*}\phi_{-1}+\frac{k_{L}}{\sqrt{2}}iq_{x}+\frac{\Omega}{\sqrt{2}};\;\;\mathcal{L}_{32}=C^{+}\phi_{+1}\phi_{0};$
$\displaystyle\mathcal{L}_{34}$
$\displaystyle=c_{0}\phi_{0}^{2}+2c_{2}\phi_{+1}\phi_{-1};\;\;\mathcal{L}_{35}=C^{+}\phi_{-1}^{*}\phi_{0}+2c_{2}\phi_{0}^{*}\phi_{+1}-\frac{k_{L}}{\sqrt{2}}iq_{x}+\frac{\Omega}{\sqrt{2}};\;\;\mathcal{L}_{36}=C^{+}\phi_{-1}\phi_{0};\;\;\mathcal{L}_{41}=-C^{+}\phi_{+1}^{*}\phi_{0}^{*}$
$\displaystyle\mathcal{L}_{42}$
$\displaystyle=-C^{+}\phi_{+1}\phi_{0}^{*}-\frac{k_{L}}{\sqrt{2}}iq_{x}-2c_{2}\phi_{0}\phi_{-1}^{*}-\frac{\Omega}{\sqrt{2}};\;\;\mathcal{L}_{43}=-c_{0}\phi_{0}^{*2}-2c_{2}\phi_{+1}^{*}\phi_{-1}^{*};\;\;\mathcal{L}_{45}=-C^{+}\phi_{-1}^{*}\phi_{0}^{*};\;\;$
$\displaystyle\mathcal{L}_{46}$
$\displaystyle=-C^{+}\phi_{-1}\phi_{0}^{*}+\frac{k_{L}}{\sqrt{2}}iq_{x}-2c_{2}\phi_{0}\phi_{+1}^{*}-\frac{\Omega}{\sqrt{2}};\;\;\mathcal{L}_{51}=C^{-}\phi_{+1}^{*}\phi_{-1};\;\;\mathcal{L}_{52}=C^{-}\phi_{+1}\phi_{-1}+c_{2}\phi_{0}^{2};\;\;$
$\displaystyle\mathcal{L}_{53}$
$\displaystyle=C^{+}\phi_{0}^{*}\phi_{-1}+\frac{k_{L}}{\sqrt{2}}iq_{x}+2c_{2}\phi_{0}\phi_{+1}^{*}+\frac{\Omega}{\sqrt{2}};\;\;\mathcal{L}_{54}=C^{+}\phi_{0}\phi_{-1};\;\;\mathcal{L}_{56}=C^{+}\phi_{-1}^{2};\;\;\mathcal{L}_{61}=-C^{-}\phi_{+1}^{*}\phi_{-1}^{*}-c_{2}\phi_{0}^{*2};$
$\displaystyle\mathcal{L}_{62}$
$\displaystyle=-C^{-}\phi_{+1}\phi_{-1}^{*};\;\;\mathcal{L}_{63}=-C^{+}\phi_{0}^{*}\phi_{-1}^{*};\;\;\mathcal{L}_{64}=-C^{+}\phi_{0}\phi_{-1}^{*}-\frac{k_{L}}{\sqrt{2}}iq_{x}-2c_{2}\phi_{0}^{8}\phi_{+1}-\frac{\Omega}{\sqrt{2}};\;\;\mathcal{L}_{65}=-C^{+}\phi_{-1}^{*2}.$

Also,

$\displaystyle C^{+}\equiv c_{0}+c_{2},\;C^{-}\equiv c_{0}-c_{2}$

The coefficients for the BdG characteristic equation (12) are given as
follows:

$\displaystyle b=$
$\displaystyle-5\Omega^{2}-4c_{2}^{2}-(2k_{L}^{2}+3\Omega+c_{0})q_{x}^{2}-\frac{3}{4}q_{x}^{4}+c_{2}(8\Omega+q_{x}^{2}),$
(14) $\displaystyle c=$ $\displaystyle
4\Omega^{4}+\Omega\left[(2\Omega(k_{L}^{2}+3\Omega)-(k_{L}^{2}-5\Omega)c_{0}\right]q_{x}^{2}+4c_{2}^{3}q_{x}^{2}+\frac{1}{2}\left[2k_{L}^{4}+9\Omega^{2}+(k_{L}^{2}+6\Omega)c_{0}\right]q_{x}^{4}+\frac{1}{2}(3\Omega+c_{0})q_{x}^{6}$
$\displaystyle+\frac{3}{16}q_{x}^{8}+4c_{2}^{2}\left[\Omega^{2}+(k_{L}^{2}-\Omega+c_{0})q_{x}^{2}\right]-\frac{1}{2}c_{2}\left[16\Omega^{3}+2\Omega(7k_{L}^{2}+5\Omega+8c_{0})q_{x}^{2}+(k_{L}^{2}+6\Omega+4c_{0})q_{x}^{4}+q_{x}^{6}\right],$
(15) $\displaystyle d=$
$\displaystyle-\frac{1}{64}q_{x}^{2}\left[(4k_{L}^{2}-4\Omega-
q_{x}^{2})(2\Omega+q_{x}^{2})+c_{2}(-8k_{L}^{2}+8\Omega+4q_{x}^{2})\right]\bigg{[}-4c_{0}\bigg{(}8\Omega^{2}-2(k_{L}^{2}-3\Omega)q_{x}^{2}+q_{x}^{4}-4c_{2}(2\Omega+q_{x}^{2})\bigg{)}$
$\displaystyle+(2\Omega+q_{x}^{2})\bigg{(}-16\Omega
c_{2}+16c_{2}^{2}-q_{x}^{2}(-4k_{L}^{2}+4\Omega+q_{x}^{2})\bigg{)}\bigg{]}.$
(16)

##  Appendix B Phase diagram for different interaction strengths

Here we present a detailed phase diagram in the $k_{L}-\Omega$ plane for the
interaction strengths other than $c_{0}=0.5$, and $c_{2}=-0.1$ as considered
in this work. For $c_{0}=0.5$, and $c_{2}=-0.1$ we have mainly observed three
different regions, on the basis of the collective excitation spectrum. Based
on the above observations, here, we show such kind of phase diagrams for two
other different sets of interaction strengths $c_{0}=5.0$, $c_{2}=-0.1$, and
$c_{0}=885.72$, $c_{2}=-4.09$. In Fig. B.1(a) we show the stability phase
diagram in the $k_{L}$-$\Omega$ plane for the interaction strength $c_{0}=5$,
$c_{2}=-0.1$. The stability phase diagram has been obtained by solving Eq.
(10) and analyzing the collective excitation spectrum for these set of
parameters. The phase diagram shows the presence of stable region (region I)
and unstable region, region II, which is given using the relation
$k_{L}^{2}=\Omega$. Region II has been divided into two subparts: region IIa,
and region IIb. The transition line that separates region IIa and IIb is given
using the relation $\Omega=0.134k_{L}^{2}-0.0668$. Where we observed that only
the boundary can be changed otherwise three regions are found.

Figure B.1: Phase diagram in the $k_{L}$\- $\Omega$ plane for different
interaction parameters: (a) For $(c_{0}=5.0$ and $c_{2}=-0.1$, and (b)
$(c_{0}=885.72$ and $c_{2}=-4.09)$. Based on the different characteristics of
the eigenspectrum and ground state, the phase diagram is divided into regions
I, IIa, IIb, and III. Region I is stable, while regions IIa and IIb are
unstable. The line $\Omega=k_{L}^{2}$ has been drawn, separating region I from
II. Blue dots separate the regions IIa and IIb. The horizontal line with
$\Omega\sim 0$ denotes the region III. The behaviour of all the regions remain
same as those for $c_{0}=0.5$, and $c_{2}=-0.1$.

The gapless behavior emerges after the point $k_{L}^{c}=0.72$,
$\Omega^{c}=0.001$, which is also the origin of the line that separates
regions IIa, and IIb. The horizontal line at which $\Omega=0.0$, we obtain
phonon mode in the low-lying as well as in the first-excited branch of the
spectrum, along with the presence of instability band, which is different than
previous cases. The cutoff value of SO coupling strength to obtain this
behavior along the line is $k_{L}=0.96$.

Figure B.1(b) illustrates a similar stability phase diagram $k_{L}$-$\Omega$
plane, featuring a very high interaction parameter strength for the
interaction parameters $c_{0}=885.72$ and $c_{2}=-4.09$. Both real and complex
eigenfrequencies coexist with gapped and gapless modes for this set of
interaction strengths. Varying the strengths of the spin-orbit and Rabi
couplings allows us to observe the transition from the stable region (region
I) to the unstable region (region II) in the phase diagram. The relationship
$k_{L}^{2}=\Omega$ defines the boundary between region I and region II. The
transition line that separates regions IIa and IIb is given by the relation
$\Omega=0.1317k_{L}^{2}-2.8512$. The origin point of the line is
$k_{L}^{c}=4.65$, $\Omega^{c}=0.02$. The horizontal line along which Rabi
coupling strength is zero, we treat it as region III, where the cut-off value
of SO coupling to achieve this behaviour is $k_{L}=8.1$.

##  Appendix C Spin-magnetization density vectors

The spatial distribution (or orientation) of spin magnetization density
vectors can be referred to as spin texture Sadler  _et al._ (2006). Comparing
the spin textures at the initial and finite times we have characterized the
dynamical behaviour of the stable (region I) and unstable regions (IIa, IIb,
and III) in the systemic way. The relevant quantities for the
characterization, the spin magnetization density vectors, are defined as,

$\displaystyle m_{x}$
$\displaystyle=\mathrm{Re}\left[\sqrt{2}(\psi_{+1}+\psi_{-1})^{*}\psi_{0}\right],$
(17a) $\displaystyle m_{y}$
$\displaystyle=\mathrm{Im}\left[\sqrt{2}(\psi_{+1}-\psi_{-1})^{*}\psi_{0}\right],$
(17b) $\displaystyle m_{z}$ $\displaystyle=|\psi_{+1}|^{2}-|\psi_{-1}|^{2}.$
(17c)

Figure C.1: Spin magnetization density components $(m_{x},m_{y},m_{z})$ for
the ground state at initial time $t=0$ (solid line), and later time $t=1000$
(dashed line) in the dynamically unstable region IIa in the $k_{L}-\Omega$
plane for different set of coupling parameters of the condensate with
ferromagnetic interactions $c_{0}=0.5$, and $c_{2}=-0.1$. (a1-c1) for
$k_{L}=2.0$, $\Omega=2.0$ and (a2-c2) for $k_{L}=2.35$, $\Omega=4.0$. Owing to
the dynamical instability, the spin texture deviates from its initial state
during the time evolution.

Figure C.2: Spin magnetization density components $(m_{x},m_{y},m_{z})$
spatial profile for the ground state(solid line) at initial time $t=0$ and
later time $t=1000$ (dashed line) in the dynamically unstable region IIb in
the $k_{L}-\Omega$ plane for different set of coupling parameters of the
condensate with ferromagnetic interactions $c_{0}=0.5$, and $c_{2}=-0.1$.
(a1-c1) for $k_{L}=3.1$ and $\Omega=1.17$ and (a2-c2) for $k_{L}=5.0$ and
$\Omega=2.5$. Owing to the dynamical instability, spin-texture shows deviation
from the initial behavior during time evolution.

In Figs. C.1 we show the components of the magnetization vector computed using
Eqs. 17a\- 17c for two sets of coupling parameters of the dynamically unstable
region IIa of the $k_{L}-\Omega$ plane as depicted in the Fig. 2. We show the
magnetization profile at two instant of time t=0 (for the ground state) and at
later instant t=1000 for $k_{L}=2.0$, $\Omega=2.0$ (see panel (a1-c1)) and
same for other parameter $k_{L}=2.35$, $\Omega=4.0$ (see panel (a2-c2)). For
both cases, all the magnetic components profile exhibits significant changes
at a later time ($t=1000$) compared to those for the ground state confirming
the dynamical instability of the condensate for these sets of the coupling
parameters in the region IIa. This trend continues for other set of points in
the region IIa. Similar trend of the magnetization vector component profile
has been observed for the two sets of the coupling parameters, $k_{L}=3.1$ and
$\Omega=1.17$ and $k_{L}=5.0$ and $\Omega=2.5$ of region IIb as depicted in
the Figs. C.2 (a1-c1) and (a2-c2) respectively.

## References

  * Zhang  _et al._ (2016) Y. Zhang, M. E. Mossman, T. Busch, P. Engels, and C. Zhang, Front. Phys. 11, 1 (2016).
  * Jacob  _et al._ (2012) D. Jacob, L. Shao, V. Corre, T. Zibold, L. De Sarlo, E. Mimoun, J. Dalibard, and F. Gerbier, Phys. Rev. A 86, 061601 (2012).
  * Campbell  _et al._ (2016) D. Campbell, R. Price, A. Putra, A. Valdés-Curiel, D. Trypogeorgos, and I. Spielman, Nat. commun. 7, 1 (2016).
  * Kasamatsu  _et al._ (2005) K. Kasamatsu, M. Tsubota, and M. Ueda, Int. J. Mod. Phys. B 19, 1835 (2005).
  * Fujimoto and Tsubota (2012) K. Fujimoto and M. Tsubota, Phys. Rev. A 85, 033642 (2012).
  * Yu (2016) Z.-Q. Yu, Phys. Rev. A 93, 033648 (2016).
  * Martone  _et al._ (2016) G. I. Martone, F. V. Pepe, P. Facchi, S. Pascazio, and S. Stringari, Phys. Rev. Lett. 117, 125301 (2016).
  * Li  _et al._ (2017) G.-Q. Li, G.-D. Chen, P. Peng, Z. Li, and X.-D. Bai, J. Phys. B: At. Mol. Opt. Phys. 50, 235302 (2017).
  * Malomed (2018) B. A. Malomed, EPL 122, 36001 (2018).
  * Peng  _et al._ (2019) P. Peng, G.-Q. Li, L.-C. Zhao, W.-L. Yang, and Z.-Y. Yang, Phys. Lett. A 383, 2883 (2019).
  * Adhikari (2021a) S. Adhikari, J. Phys.: Condens. Matter 33, 265402 (2021a).
  * Gautam and Adhikari (2014) S. Gautam and S. K. Adhikari, Phys. Rev. A 90, 043619 (2014).
  * Adhikari (2019) S. K. Adhikari, Phys. Rev. A 100, 063618 (2019).
  * Gautam and Adhikari (2015) S. Gautam and S. K. Adhikari, Laser Phys. Lett. 12, 045501 (2015).
  * Adhikari (2020) S. K. Adhikari, Physica E Low Dimens. Syst. Nanostruct. 118, 113892 (2020).
  * Gautam and Adhikari (2021) S. Gautam and S. Adhikari, Braz. J. Phys. 51, 298 (2021).
  * Mardonov  _et al._ (2015) S. Mardonov, M. Modugno, and E. Y. Sherman, Phys. Rev. Lett. 115, 180402 (2015).
  * Mardonov  _et al._ (2018) S. Mardonov, V. V. Konotop, B. A. Malomed, M. Modugno, and E. Y. Sherman, Phys. Rev. A 98, 023604 (2018).
  * Mardonov  _et al._ (2019) S. Mardonov, M. Modugno, E. Y. Sherman, and B. A. Malomed, Phys. Rev. A 99, 013611 (2019).
  * Mithun and Kasamatsu (2019) T. Mithun and K. Kasamatsu, J. Phys. B: At. Mol. Opt. Phys. 52, 045301 (2019).
  * Cabedo  _et al._ (2021) J. Cabedo, J. Claramunt, and A. Celi, Phys. Rev. A 104, L031305 (2021).
  * Zhu  _et al._ (2020) Q.-L. Zhu, L. Pan, and J. An, Phys. Rev. A 102, 053320 (2020).
  * Lin  _et al._ (2011) Y.-J. Lin, K. Jiménez-García, and I. B. Spielman, Nature (London) 471, 83 (2011).
  * Ravisankar  _et al._ (2021a) R. Ravisankar, T. Sriraman, R. K. Kumar, P. Muruganandam, and P. K. Mishra, J. Phys. B: At. Mol. Opt. Phys. 54, 225301 (2021a).
  * Goldstein and Meystre (1997) E. V. Goldstein and P. Meystre, Phys. Rev. A 55, 2935 (1997).
  * Zhang  _et al._ (2005) W. Zhang, D. L. Zhou, M.-S. Chang, M. S. Chapman, and L. You, Phys. Rev. Lett. 95, 180403 (2005).
  * Pethick and Smith (2008) C. J. Pethick and H. Smith,  _Bose–Einstein condensation in dilute gases_ (Cambridge university press, 2008).
  * Pitaevskii and Stringari (2016) L. Pitaevskii and S. Stringari,  _Bose-Einstein condensation and superfluidity_ , Vol. 164 (Oxford University Press, 2016).
  * Bogoliubov (1947) N. Bogoliubov, J. Phys 11, 23 (1947).
  * Jin  _et al._ (1996) D. S. Jin, J. R. Ensher, M. R. Matthews, C. E. Wieman, and E. A. Cornell, Phys. Rev. Lett. 77, 420 (1996).
  * Mewes  _et al._ (1996) M.-O. Mewes, M. R. Andrews, N. J. van Druten, D. M. Kurn, D. S. Durfee, C. G. Townsend, and W. Ketterle, Phys. Rev. Lett. 77, 988 (1996).
  * Martone  _et al._ (2012) G. I. Martone, Y. Li, L. P. Pitaevskii, and S. Stringari, Phys. Rev. A 86, 063621 (2012).
  * Li  _et al._ (2013) Y. Li, G. I. Martone, L. P. Pitaevskii, and S. Stringari, Phys. Rev. Lett. 110, 235302 (2013).
  * Khamehchi  _et al._ (2014) M. A. Khamehchi, Y. Zhang, C. Hamner, T. Busch, and P. Engels, Phys. Rev. A 90, 063624 (2014).
  * Chen  _et al._ (2017) L. Chen, H. Pu, Z.-Q. Yu, and Y. Zhang, Phys. Rev. A 95, 033616 (2017).
  * Chen  _et al._ (2022) Y. Chen, H. Lyu, Y. Xu, and Y. Zhang, New J. Phys. 24, 073041 (2022).
  * He and Lin (2023) J. He and J. Lin, New J. Phys. 25, 093041 (2023).
  * Ozawa  _et al._ (2013) T. Ozawa, L. P. Pitaevskii, and S. Stringari, Phys. Rev. A 87, 063610 (2013).
  * Ravisankar  _et al._ (2021b) R. Ravisankar, H. Fabrelli, A. Gammal, P. Muruganandam, and P. K. Mishra, Phys. Rev. A 104, 053315 (2021b).
  * Katsimiga  _et al._ (2021) G. Katsimiga, S. Mistakidis, P. Schmelcher, and P. Kevrekidis, New J. Phys. 23, 013015 (2021).
  * Katsimiga  _et al._ (2023) G. C. Katsimiga, S. I. Mistakidis, K. Mukherjee, P. G. Kevrekidis, and P. Schmelcher, Phys. Rev. A 107, 013313 (2023).
  * Rajat  _et al._ (2022) Rajat, A. Roy, and S. Gautam, Phys. Rev. A 106, 013304 (2022).
  * Cross and Hohenberg (1993) M. C. Cross and P. C. Hohenberg, Rev. Mod. Phys. 65, 851 (1993).
  * Bernier  _et al._ (2014) N. R. Bernier, E. G. Dalla Torre, and E. Demler, Phys. Rev. Lett. 113, 065303 (2014).
  * Salasnich  _et al._ (2002) L. Salasnich, A. Parola, and L. Reatto, Phys. Rev. A 65, 043614 (2002).
  * Adhikari (2021b) S. K. Adhikari, Phys. Rev. E 104, 024207 (2021b).
  * Kawaguchi and Ueda (2012) Y. Kawaguchi and M. Ueda, Phys. Rep. 520, 253 (2012).
  * Stamper-Kurn and Ueda (2013) D. M. Stamper-Kurn and M. Ueda, Rev. Mod. Phys. 85, 1191 (2013).
  * Ravisankar  _et al._ (2021c) R. Ravisankar, D. Vudragović, P. Muruganandam, A. Balaž, and S. K. Adhikari, Comput. Phys. Commun. 259, 107657 (2021c).
  * Inouye  _et al._ (1998) S. Inouye, M. Andrews, J. Stenger, H.-J. Miesner, D. M. Stamper-Kurn, and W. Ketterle, Nature (London) 392, 151 (1998).
  * Marte  _et al._ (2002) A. Marte, T. Volz, J. Schuster, S. Dürr, G. Rempe, E. Van Kempen, and B. Verhaar, Phys. Rev. Lett. 89, 283202 (2002).
  * Chin  _et al._ (2010) C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, Rev. Mod. Phys. 82, 1225 (2010).
  * Wen  _et al._ (2012) L. Wen, Q. Sun, H. Q. Wang, A. C. Ji, and W. M. Liu, Phys. Rev. A 86, 043602 (2012).
  * Rajaswathi  _et al._ (2023) K. Rajaswathi, S. Bhuvaneswari, R. Radha, and P. Muruganandam, Phys. Rev. A 108, 033317 (2023).
  * Wang  _et al._ (2010) C. Wang, C. Gao, C.-M. Jian, and H. Zhai, Phys. Rev. Lett. 105, 160403 (2010).
  * Zhu  _et al._ (2012) Q. Zhu, C. Zhang, and B. Wu, EPL 100, 50003 (2012).
  * Anderson  _et al._ (1999) E. Anderson, Z. Bai, C. Bischof, L. S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, _et al._ ,  _LAPACK users’ guide_ (SIAM, 1999).
  * Muruganandam and Adhikari (2009) P. Muruganandam and S. K. Adhikari, Comput. Phys. Commun. 180, 1888 (2009).
  * Young-S  _et al._ (2016) L. E. Young-S, D. Vudragović, P. Muruganandam, S. K. Adhikari, and A. Balaž, Comput. Phys. Commun. 204, 209 (2016).
  * Tasgal and Band (2015) R. S. Tasgal and Y. B. Band, Phys. Rev. A 91, 013615 (2015).
  * Abad and Recati (2013) M. Abad and A. Recati, Eur. Phys. J. D 67, 148 (2013).
  * Lyu and Zhang (2020) H. Lyu and Y. Zhang, Phys. Rev. A 102, 023327 (2020).
  * Cabrera  _et al._ (2018) C. R. Cabrera, L. Tanzi, J. Sanz, B. Naylor, P. Thomas, P. Cheiney, and L. Tarruell, Science 359, 301 (2018).
  * Petrov (2015) D. S. Petrov, Phys. Rev. Lett. 115, 155302 (2015).
  * Gangwar  _et al._ (2022) S. Gangwar, R. Ravisankar, P. Muruganandam, and P. K. Mishra, Phys. Rev. A 106, 063315 (2022).
  * Gangwar  _et al._ (2023) S. Gangwar, R. Ravisankar, S. Mistakidis, P. Muruganandam, and P. K. Mishra, arXiv preprint arXiv:2307.16742 (2023).
  * Sadler  _et al._ (2006) L. Sadler, J. Higbie, S. Leslie, M. Vengalattore, and D. Stamper-Kurn, Nature 443, 312 (2006).