# Dynamics of quantum soliton in Lee-Huang-Yang spin-orbit coupled Bose-
Einstein condensates

Sonali Gangwar Department of Physics, Indian Institute of Technology, Guwahati
781039, Assam, India Rajamanickam Ravisankar Department of Physics, Indian
Institute of Technology, Guwahati 781039, Assam, India Institute of Atomic and
Molecular Sciences, Academia Sinica, Taipei 10617, Taiwan Paulsamy
Muruganandam Department of Physics, Bharathidasan University, Tiruchirappalli
620024, Tamilnadu, India Pankaj Kumar Mishra Department of Physics, Indian
Institute of Technology, Guwahati 781039, Assam, India

###### Abstract

We present the numerical results of the structure and dynamics of the self-
bound ground state arising solely because of the presence of beyond mean field
quantum fluctuation in spin-orbit coupled binary Bose-Einstein condensates in
one dimension. Depending upon spin-obit (SO) and Rabi couplings, we observe
that the ground state exhibits either quantum-bright (plane) or quantum-stripe
soliton nature. We find an analytical soliton solution for non-zero SO
coupling that matches quite well with the numerical results. Further, we
investigate the dynamical stability of these solitons by adopting three
protocols, such as (i) adding initial velocity to each component, (ii)
quenching the SO and Rabi coupling parameters at initial and finite time, and
(iii) allowing collision between the two spin-components by giving equal and
opposite direction velocity to them. Many interesting dynamical features of
the solitons, like, multi-fragmented, repelling, and breathing in time and
space-time, are observed. For given Rabi coupling frequency, the breathing
frequency of the soliton increases upon the increase in SO coupling, attaining
a maximum at the critical SO coupling where the phase transition from the
bright to stripe soliton occurs. We observe that the maximum breathing
frequency exhibits power law dependence on the Rabi coupling frequency with an
exponent $\sim 0.16$. The quenching of SO and Rabi coupling frequency leads to
dynamical phase transition from bright to stripe soliton and from breathing
stripe to a bright soliton. In the absence of the SO and Rabi couplings,
depending upon the velocity of the up-spin and down-spin components, the
collision between them is either elastic or inelastic, which is consistent
with the earlier numerical and experimental observation. However, in the
presence of coupling parameters, the collision appears to be inelastic and
quasi-elastic in nature.

##  I Introduction

An exciting aspect of Bose-Einstein condensates (BECs) is the formation of the
quantum droplet, which appears in the form of small clusters of atoms bound
together owing to the balance between the attractive and repulsive forces. In
a binary mixture of the condensate, the quantum droplet manifests in a
spherical shape due to the competition between the effective short-range
attractive interaction between the atoms and the repulsive interaction arising
solely due to quantum fluctuation, which is also responsible for its
stabilization [1]. Theoretically, one can model the contribution of the
quantum fluctuations by considering the beyond mean-field (BMF) term,
popularly known as the Lee-Huang-Yang (LHY) correction in the mean-field
Hamiltonian of the condensate [2]. BMF terms were corroborated first time in
experiments with homogeneous and single component Bose gas of 85Rb [3] and 7Li
[4], and recently in binary 39K [5]. In general, the LHY correction term is
attractive in quasi-one dimension, while it is repulsive in two and three
dimension [1, 6].

The advancement in the state-of-art technology and simulation tools in
ultracold BECs led to an upsurge in the research in quantum droplets during
the past few decades. The quantum droplet was first reported experimentally in
the dipolar BECs, for which a single component-dipolar BECs made of dysprosium
(164Dy) produced elongated quantum droplets in one direction [7, 8, 9]. Chomaz
et al. observed similar features with the condensate of erbium atoms (166Er)
[10]. Subsequently, several groups realized quantum droplets in binary
mixtures of isotropic BECs [11, 12, 13, 14], Bose-Fermi mixture [15] and also
in binary magnetic gases [16]. These quantum droplets are commonly observed in
three dimensions, while it is possible to increase the lifetime of the
droplets in the lower dimensions. The formation of one-dimensional quantum
droplets is due to the balance between the repulsive mean-field (MF)
contribution to the energy per particle, which is linear in the density ($n$)
of the gas, and the attractive BMF correction, proportional to $-n^{1/2}$ [6,
17].

Since the proposition of the non-trivial attractive nature of the LHY term in
quasi-one dimensional quantum droplet [17], it has caught the great attention
of the scientific community. A significant number of experimental [18], as
well as theoretical and numerical works, have been performed in recent years
using the effective one component [19], binary [20] and Spin-orbit (SO) and
Rabi coupled binary BECs [21, 22] that explored the structure and dynamics of
the Quantum droplet. Tononi et al., in SO-coupled BECs, demonstrated that the
self-bound states are solitonic in nature for vanishing mean-field
contribution [21]. Depending upon the Rabi and SO coupling parameters range,
these soliton-like states are either of single peak (bright soliton) or
multiple peaks (stripe soliton) nature [23, 21]. Some of these works
demonstrate the existence of soliton and droplet nature of the self-bound
state for the effective-one component binary BECs in one dimension [24] as
well as in the binary mixture with SO coupling in two dimension [25, 22]. The
transition between the quantum soliton and quantum droplet regime depends on
several parameters, like, the atom number ($N$), interaction strength, the
strength of the confining potential, etc. [24]. For large $N$, a highly dense
droplet solution exists, while low dense bright soliton occurs for small $N$
[19]. A recent experiment shows that the attractive mixture of BECs confined
in an optical waveguide could exhibit both quantum soliton and quantum
droplets. At large $N$, the droplet exhibits similar characteristics as those
of the classical droplets [26].

Though lots of emphasis on exploring the structure and stability of quantum
droplets in the recent past, only a limited number of works are available in
the literature that focuses on the dynamical aspect of quantum droplets.
Ferioli et al. analyzed the dynamics of quantum droplets in a binary mixture
by allowing the droplets to collide with each other [14]. In classical droplet
collisions, there are two possibilities, either they merge into a single one,
or the colliding droplets separate two or more ones after collision [27].
Ferioli et al. have investigated the collision between the quantum droplets by
experimental and numerical means for the binary mixture of hyperfine states in
39K and found the merging and separation between the droplets depending on
their velocity. Critical velocity for the transition between the elastic and
inelastic collision strongly depends upon the atom number $N$ [14]. A few
theoretical studies are available that establish the robustness of the quantum
soliton or quantum droplet during the collision. For instance, Astrakharchik
and Malomed numerically investigated the static and dynamical properties of
quantum droplets using the mean-field theory and studied the collision
properties of two counter-propagating droplets [19]. They showed that
collisions between the tiny droplets are quasi-elastic, indicating solitonic
behaviour. While depending on the velocity, the large drops may undergo
merging or fragmentation, which shows the collision as inelastic. Recently,
Young and Adhikari numerically studied the collision properties of bright
solitons in two-dimensional dipolar BECs [28].

In recent years, one may witness numerous theoretical and numerical works
focusing on the structure and dynamics of the droplets in binary BECs.
However, only a few studies are available on the droplets in the spinor BECs,
especially in spin-orbit coupled binary BECs [21, 29]. In particular, the role
of SO and Rabi coupling parameters in dictating the shape and dynamics of the
quantum droplet is not well understood. In this paper, we present a detailed
numerical investigation to understand the effect of the SO and Rabi coupling
parameters on the stability and shape of the quantum soliton. Although Tononi
et al. [21] demonstrates the presence of the self-bound quantum soliton due to
the LHY term with vanishing mean-field contributions, several aspects, like
dynamical robustness of the ground state for different perturbations, have not
been explored yet. In this paper, we have performed a systematic analysis of
the effect of initial velocity, quenching of the SO and Rabi coupling
parameters, and allowing the collision between the components by initially
perturbing the components with equal and opposite speeds. All these protocols
facilitate us to obtain a variety of dynamical phases that include breathing
(both in space and time), repelling, multi-fragmented solitons, etc.

The structure of our paper is as follows. In Sec. II, we present governing
equations and numerical simulation details and outline a possible scenario to
connect our numerical parameters with the experiment. We illustrate an
analytical solution of the ground state of the quantum soliton in the presence
of SO coupling in Sec. III. Following this in Sec. IV, we present a detailed
analysis of the ground state structure followed by its dynamics which are set
up in the system by different procedures. First, we discuss the different
sorts of dynamics that arise due to the initial velocity given to the soliton,
followed by the dynamics due to the quenching of coupling parameters for both
zero and finite speeds. Further, we highlight some of the pronounced dynamical
behaviour shown by the solitons in the presence of collisions. Finally, we
conclude our work in Sec. V.

##  II Beyond Mean-field model for SO coupled BECs

We consider a pseudo-spin-$1/2$ Bose-Einstein condensate confined in strong
transverse confinement modelled using quasi-one-dimensional spin-orbit and
Rabi coupled condensates with LHY term. The corresponding coupled Gross-
Pitaevskii (GP) equations in dimensionless form are given by [21]:

$\displaystyle\mathrm{i}\partial_{t}\psi_{\uparrow}=\bigg{[}$
$\displaystyle-\frac{1}{2}\partial_{x}^{2}-\mathrm{i}k_{L}\partial_{x}+g\lvert\psi_{\uparrow}\rvert^{2}+g_{\uparrow\downarrow}\lvert\psi_{\downarrow}\rvert^{2}$
$\displaystyle-\dfrac{g_{\text{LHY}}^{3/2}}{\pi}\sqrt{\lvert\psi_{\uparrow}\rvert^{2}+\lvert\psi_{\downarrow}\rvert^{2}}\bigg{]}\psi_{\uparrow}+\Omega\psi_{\downarrow},$
(1a) $\displaystyle\mathrm{i}\partial_{t}\psi_{\downarrow}=\bigg{[}$
$\displaystyle-\frac{1}{2}\partial_{x}^{2}+\mathrm{i}k_{L}\partial_{x}+g_{\downarrow\uparrow}\lvert\psi_{\uparrow}\rvert^{2}+g\lvert\psi_{\downarrow}\rvert^{2}$
$\displaystyle-\dfrac{g_{\text{LHY}}^{3/2}}{\pi}\sqrt{\lvert\psi_{\uparrow}\rvert^{2}+\lvert\psi_{\downarrow}\rvert^{2}}\bigg{]}\psi_{\downarrow}+\Omega\psi_{\uparrow},$
(1b)

where $\psi_{\uparrow}$ and $\psi_{\downarrow}$ are the wavefunctions of the
spin-up and spin-down components respectively, $k_{L}$ is the spin-orbit
coupling strength, $\Omega$ is the Rabi coupling frequency, $g$ is the
intraspecies interaction, and $g_{\uparrow\downarrow}$ are the interspecies
interaction strengths. We consider the interaction term due to the LHY
correction as $g_{\text{LHY}}=g$. The wave functions are subjected to the
following normalization condition:

$\displaystyle\int\limits_{-\infty}^{\infty}\left(\lvert\psi_{\uparrow}\rvert^{2}+\lvert\psi_{\downarrow}\rvert^{2}\right)\,dx=1,$
(2)

and also that the total density
$n=n_{\uparrow}+n_{\downarrow}\equiv\lvert\psi_{\uparrow}\rvert^{2}+\lvert\psi_{\downarrow}\rvert^{2}$
is conserved.

Equations (1) are non-dimensionalized using transverse harmonic oscillator
length $a_{0}=\sqrt{\hbar/(m\omega_{\perp})}$ as a characteristic length scale
(where, $\omega_{\perp}$ is the trap frequency in the transverse direction),
$\omega_{\perp}^{-1}$ as a time scale, and $\hbar\omega_{\perp}$ as an energy
scale. The other interaction parameters are considered as $g$ =
$2Na_{\uparrow\uparrow}/a_{0}$ and
$g_{\uparrow\downarrow}=2Na_{\uparrow\downarrow}/a_{0}$. Here,
$a_{\uparrow\uparrow}$ and $a_{\uparrow\downarrow}$ represent the scattering
length corresponding to intra and inter components, respectively. The SO and
the Rabi coupling parameters have been rescaled as $k_{L}\to k_{L}a_{0}$ and
$\Omega\to{\Omega}/\omega_{\perp}$, respectively, while the wave function is
rescaled as
$\psi_{\uparrow,\downarrow}=\psi_{\uparrow,\downarrow}\sqrt{a_{0}}$.

To make the numerical simulation experimentally viable, we choose the
parameters same as considered in the recent realization of the quantum droplet
in the binary hyperfine states mixture of 39K condensates [11, 12, 14].
Following this we consider $N\sim 10^{4}$ atoms confined in the harmonic trap
potential with frequencies $\omega_{x}=2\pi\times 50$ Hz,
$\omega_{\perp}=2\pi\times 800$ Hz in perpendicular and axial directions,
respectively. Using this the characteristic length scale can be obtained as
$a_{\perp}\sim 0.6\,\mu\mbox{m}$. Generally, in the experiment two internal
hyperfine states $\lvert F=1,m_{F}=-1\rangle$ and $\lvert F=1,m_{F}=0\rangle$
are considered which can be attributed, respectively, to the pseudo-spin up
$\lvert\uparrow\rangle$, and pseudo-spin down $\lvert\downarrow\rangle$ states
of our model. These two spin states have an equal number of atoms, and intra-
and interspecies interaction strengths can be controlled by tuning $s$-wave
scattering lengths through Feshbach resonance and by varying the magnetic
field. Following the experiment we set
$a_{\uparrow\uparrow}=a_{\downarrow\downarrow}=a_{\uparrow\downarrow}=0.2686a_{0}$
($a_{0}$ is the Bohr radius) which gives the dimensionless interaction
strengths as $g=g_{\text{LHY}}\approx 0.5$. Note that here we have taken all
the scattering lengths to be the same to make the contribution of the mean-
field term on the structure and dynamics of the soliton vanishingly small.
However, the recent experiment on the quantum droplet with ${}^{39}K$
indicates the presence of soliton nature of the condensate is only possible
when all the scattering lengths are not equal [12]. We expect our results will
not alter much upon considering the unequal scattering lengths, as observed in
the experiment. Another parameter in this system is the Rabi coupling
frequency ($\Omega$), commonly used for coupling the spin states by tuning the
frequency of Raman lasers ranging from
$\Omega=2\pi\times\\{0.080-40\\}\mbox{kHz}$, which implies that the
dimensionless Rabi coupling frequency range as $\Omega=\\{0.1,50\\}$. It is
possible to vary the SO coupling strength ($k_{L}$) with the laser wavelength
and their geometry. In the present work we have considered the dimensionless
SO coupling range $k_{L}=\\{0.1,8\\}$, where the laser wavelength varies as
$\lambda_{L}=17.8~{}\mu\mbox{m}-223.2~{}\mbox{nm}$. Note that the LHY
interaction is not equal to the MF interaction because the LHY nonlinear
interaction term appears with an extra factor $\pi$. We have considered
$g_{\text{LHY}}=g=0.5$ and $1.0$, for which the quantum fluctuation term
contributes approximately of the order of $11\%$ and $30\%$ respectively.
Interestingly albeit of small contribution, we have obtained a large variety
of dynamical features for our system without any confinement.

We employ the imaginary time propagation method with the aid of a split-step
Crank-Nicolson scheme [30, 31, 32] to numerically solve the coupled GP
equations (1a) and (1b). The box size $[-153.6:153.6]$ with spatial resolution
as $dx=0.025$ is chosen for all the simulation runs, and we use the Gaussian
initial condition with anti-symmetric profiles on the components, i.e.,
${\psi}_{\uparrow}(x)=-{\psi}_{\downarrow}(-x)$. The time step is fixed at
$dt=10^{-5}$. We also consider various interaction strength ($g=0.5,1,2$) for
our studies.

##  III Analytical soliton solution

Here, we present the analytical solution of the ground state for our model.
For brevity, we wish to obtain a solution for zero Rabi coupling frequency
($\Omega=0$) and finite $k_{L}$. For this purpose, we consider the following
transformation to eliminate the SO coupling term from Eq. (1),

$\displaystyle\psi_{\uparrow}(x,t)$
$\displaystyle=\widetilde{\psi}_{\uparrow}(x,t)\exp{\left[\frac{\mathrm{i}k_{L}}{2}\left(k_{L}t-2x\right)\right]}$
(3a) $\displaystyle\psi_{\downarrow}(x,t)$
$\displaystyle=\widetilde{\psi}_{\downarrow}(x,t)\exp{\left[\frac{\mathrm{i}k_{L}}{2}\left(k_{L}t+2x\right)\right]}$
(3b)

Using the single field approximation we get
$\widetilde{\psi}_{\uparrow}=\widetilde{\psi}_{\downarrow}=\phi(x)\exp{\left(\mathrm{i}\mu
t\right)}$. Based on the system parameters $g_{\uparrow\downarrow}=-g$, we
find that the mean-field interaction term cancels each other, which yields the
stationary state solution of the form

$\displaystyle-\mu\phi=-\frac{1}{2}\partial_{x}^{2}\phi-\frac{\sqrt{2}}{\pi}g^{3/2}\phi^{2}$
(4)

which gives

$\displaystyle\phi^{\prime 2}=\alpha\phi^{2}-\beta\phi^{3}$ (5)

where prime (′) represents the spatial derivative with respect to $x$ and
$\alpha=2\mu,\beta=2\sqrt{2}g^{3/2}/3\pi$. Now the solution of the Eq. (5)
becomes

$\displaystyle\phi(x)=\phi(0)~{}\text{sech}^{2}(\sqrt{\mu/2}~{}x)$ (6)

where, $\phi(0)=\alpha/\beta$. Here, the stationary state solution depends on
the chemical potential. Using the normalization condition
$\int\lvert\psi\rvert^{2}dx=1$ we obtain
$2\lvert\mu\rvert=(g^{2}/3^{2/3}\pi^{4/3})$. Thus the final solution of the
Eq. (1) have the form as,

$\displaystyle\psi_{\uparrow}=\left(\frac{\alpha}{\beta}\right)\text{sech}^{2}\left(\sqrt{\frac{\mu}{2}}~{}x\right)\exp{\left[\frac{\mathrm{i}k_{L}}{2}(k_{L}t-2x)+\mathrm{i}\mu
t\right]}$
$\displaystyle\psi_{\downarrow}=\left(\frac{\alpha}{\beta}\right)\text{sech}^{2}\left(\sqrt{\frac{\mu}{2}}~{}x\right)\exp{\left[\frac{\mathrm{i}k_{L}}{2}(k_{L}t+2x)+\mathrm{i}\mu
t\right]}$ (7)

In Fig. 1, we plot the analytical (red line) solution [Eqs. (7)] as well
numerical simulation (black dots) results of the stationary ground state of
the Eq. (1) for $g_{LHY}=g=0.5$, with $\Omega=0$ and $k_{L}=1$.

Figure 1: A comparison between the numerically obtained soliton profiles
(dotted-black line) and analytical (solid-red line) solution (7) for
interaction strengths $g=0.5$, and coupling parameters $k_{L}=1$ and
$\Omega=0$. Real parts of the spin components: (a)
$\mathrm{Re}(\psi_{\uparrow})$ and (b) $\mathrm{Re}(\psi_{\downarrow})$,
imaginary parts of the spin components: (c) $\mathrm{Im}(\psi_{\uparrow})$ and
(d) $\mathrm{Im}(\psi_{\downarrow})$, and the densities of the (e) up
$\lvert\psi_{\uparrow}\rvert^{2}$ and (f) down
$\lvert\psi_{\downarrow}\rvert^{2}$ components.

Both analytical and numerical simulation results match reasonably well for the
real and imaginary parts and the total density of components.

In brief, we found the approximate soliton solution for the system of coupled
Eqs. (1) by eliminating the SO coupling term using the transformation with
zero Rabi coupling frequency. In the process, we have used single field
approximation, which provides $\text{sech}^{2}(x)$ soliton solution due to
quadratic nonlinearity [21].

##  IV Numerical results

In this section, first, we present our numerical results for different ground
state solitons obtained using the imaginary-time propagation scheme. Our main
emphasis is to ascertain the role of mean-field and beyond mean-field LHY term
on the overall shape and structure of the self-bound quantum soliton state.
Following this, we focus on the soliton dynamics by employing the real-time
propagation of the governing dynamical GP equations. Finally, we present a
detailed analysis of the dynamics of quantum solitons by giving initial
velocity, quenching the coupling parameters and allowing the collision between
the components by initially imparting them with an equal speed in the opposite
direction.

###  IV.1 Stationary ground states of quantum solitons

To better understand the role of the beyond mean-field term ($g_{\text{LHY}}$)
along with the SO and Rabi coupling frequency on the shape and size of the
ground state, we now proceed to get the different ground state of the solitons
by assuming $g=-g_{\uparrow\downarrow}=g_{\text{LHY}}$. Here, we use MF and
LHY ground states to represent ground states in the absence and presence of
LHY correction, respectively. In general, depending upon the coupling
parameters range, the ground states are either bright soliton (BS) or stripe
soliton (SS) in nature [21, 33].

We begin the analysis of the structure of BS in absence of $g_{\text{LHY}}$
for the fixed parameters $k_{L}=\Omega=1$ and $g=-g_{\uparrow\downarrow}=0.5$.
In Fig. 2(a), the solid black line represents the BS state. In the presence of
$g_{\text{LHY}}$ magnitude of the BS gets doubled, and it tends to become more
localized even though the width of the soliton appears to get reduced, as
represented by the dashed-red line in Fig. 2(a).

Figure 2:  Ground state density profiles of the spin-up component in the
presence and absence of LHY term. Bright soliton for $\Omega=1$ and $k_{L}=1$
with (a) $g=0.5$ and (b) $g=1.0$, and stripe soliton for $\Omega=1$ and
$k_{L}=4$ with (a) $g=0.5$ and (d) $g=1.0$. The spin density gets more
confined when the LHY term (dashed-red line) is present. The insets in (c) and
(d) show the magnified views of the stripe patterns in the presence and
absence of the LHY term.  Figure 3: Temporal evolution of mean-field (MF)
bright soliton (BS) at different initial velocities: (a) $v=0$, (b) $v=0.2$,
(c) $v=0.5$, and (d) $v=1$ for $g=-g_{\uparrow\downarrow}=0.5$ and
$\Omega=k_{L}=1$. Upon increasing the velocity from $v=0$ to $v=0.2$, a
transition from the expanding soliton to an oscillating soliton takes place.
For large velocity multi-soliton (at $v=0.5$) and bifurcated solitons (at
$v=1$) are observed. Figure 4: Temporal evolution of MF stripe soliton for
different velocities: (a) $v=0$ and (b) $v=0.4$ with
$g=-g_{\uparrow\downarrow}=0.5$, $\Omega=1$, and $k_{L}=4$. For $v=0$ shape of
the soliton remains unchanged until $t\sim 200$. An expansion occurs beyond
this time. For finite velocity ($v=0.4$), the stripe soliton shows propagation
along the initial velocity direction.

Further upon increasing the interaction strengths to
$g=-g_{\uparrow\downarrow}=1$, we find a negligible change in the amplitude
and width of the mean-field Bright soliton (MF-BS) [see Fig. 2(b) (solid-black
line)]. However, in the presence of LHY interaction ($g_{\text{LHY}}=1$) a
well localized quantum bright soliton (QBS) width gets reduced by half, while,
its amplitude is increased by a factor $\approx 2.66$ [see Fig. 2(b) (dashed-
red line)].

Figure 5: Plots of the total density,
$\sum_{j=\uparrow,\downarrow}\lvert\psi_{j}\rvert^{2}$, showing dynamics of
quantum bright soliton (QBS) with LHY correction at different velocities: (a)
soliton propagates with no shape change for $v=0$, (b) breathers for finite
velocity $v=0.5$, and (c) oscillation and bifurcation for $v=1.0$. Increase in
the velocity leads transitions from bright soliton to breather which further
gets transformed into multi-solitons.  Figure 6: Plots of the total density,
$\sum_{j=\uparrow,\downarrow}\lvert\psi_{j}\rvert^{2}$, showing the dynamics
of quantum stripe soliton (QSS) with LHY correction at different velocities:
(a) $v=0$, (b) $v=0.2$, and (c) $v=1.0$ for $g=-g_{\uparrow\downarrow}=0.5$
and $\Omega=1,k_{L}=4$. For $v=0$ stable stripe soliton is observed. However,
for finite velocity breathers ($v=0.2$) and bifurcated ($v=1.0$) solitons are
observed.

Next, we consider the shape and structure of SS both in presence and absence
of LHY correction term which was obtained for $\Omega=1$ and $k_{L}=4$ at
$g=0.5,1$. A phase transition from BS to SS states occurs as we change the SO
coupling strength from $k_{L}=1$ to $k_{L}=4$. The fringes that appear in the
density profiles characterize the stripe wave pattern. One may note that
several local maxima appear in the density profile for both the MF and quantum
solitons. In the presence of the LHY term with interaction strength $g=0.5$,
the number of stripes, as well as the width, gets reduced, and amplitude gets
doubled compared to those without LHY, indicating the localization of the
soliton [see Fig. 2(c)]. On increasing the interaction strength to $g=1$, we
find that the mean-field stripe soliton (MF-SS) exhibits similar nature as
shown with the black line in Fig. 2(d)]. However, the LHY-SS or quantum SS
(QSS) phase width gets reduced by half, and an increase in the amplitude by a
factor $\approx 4$, accompanied by a loss in the number of stripes and an
increment in the localization [shown with the red line in Fig. 2(d)] is
observed. Overall, we noticed that the condensate shape hardly changes upon
increasing the MF interaction strengths. However, we note a significant change
in the shape and amplitude when we consider the LHY correction
($g_{\text{LHY}}\neq 0$). With the LHY correction on increasing the nonlinear
interaction ($g$), the size of the soliton gets reduced while the amplitude
increases.

###  IV.2 Dynamics of different phases of quantum soliton

In this section, we present the dynamics of the ground state of the MF and
quantum soliton by solving the governing equation [cf. Eqs. (1)] with the help
of real-time propagation. The main aim here is to investigate the dynamics of
the MF and quantum solitons by giving some initial velocity to the condensate
by making a uniform change in the phase of the soliton wavefunction [34, 21].
Further, we demonstrate breathing solitons, fragmented soliton, free expansion
of soliton, etc., by manipulating the magnitude of the initial velocity.

Figure 7: (a) Variation of breathing frequency ($\omega_{bf}$) with the SO
coupling strength ($k_{L}$) for different $\Omega$ with fixed velocity
$v=0.1$. (b) The Variation of $\omega_{bf}$ with $k_{L}$ for different Rabi
coupling ($\Omega=1,4$) and velocities $v=0.1,0.2$ with $g=0.5$.
$\omega^{c}_{bf}$ at the critical $k_{L}$ where PW-SS phase transition happens
and remains independent of the velocity. (c) Variation of $\omega_{bf}$ with
$k_{L}$ for different interaction strengths $g=0.5,1.0,2.0$ at $\Omega=4$ and
$v=0.1$. To bring the graph to the same scale, we multiply the data for
$g=0.5$ and $g=1.0$ by a factor of two.

Fig. 3 depicts the temporal evolution of MF-BS for different velocities,
$v=0,0.2,0.5$ and $1$ with $g=-g_{\uparrow\downarrow}=0.5$ and
$\Omega=k_{L}=1$. For $v=0$, the soliton propagates without distortion until
$t\sim 150$, and after this time a slight expansion is observed [see Fig.
3(a)]. This feature is quite evident as the mean-field contribution vanishes,
and the system does not have any trap. In this case, the soliton maintains its
shape over a definite extent of time, beyond which the interactions could not
stabilise the shape, resulting in the expansion indicating the appearance of
dynamical instability [35]. For a small but finite initial velocity (for
example, $v=0.2$), we notice setting up a time-dependent oscillation in the
soliton that manifests as an undulating motion of the density in time and
space. Beyond $t\sim 30$, the oscillation vanishes and stable soliton is
observed [cf. Fig. 3(b)]. This oscillation can be associated with the
instability that soliton displays upon making a sudden change in the phase by
a factor $\exp{(\pm\mathrm{i}vx)}$, where $v$ is the resultant acquired
velocity. Similar to the case $v=0$, we observe an expansion in the soliton
upon waiting for a longer period. For $v=0.5$, soliton appears to display
breathers-like oscillation in which it exhibits expansion followed by
compression at a periodic time interval [see Fig. 3(c)]. The spin component
exhibits interesting dynamical behaviour. We noticed the absence of the
breather on one side while it is present on the other side. However, the
middle lobe exhibits oscillation for a shorter time and expansion in a long
duration indicates the presence of an oscillation death-like behaviour in the
soliton [35]. In addition, we also observe the presence of multi-soliton-like
features. At $v=1$, fragmented solitons with breather-like modes get generated
and propagate in the respective directions for the spin components,
accompanied by the appearance of bifurcation-like solitons [cf. Fig. 3(d)].

We observe similar nature of dynamical expansion in the SS as those observed
for BS in the absence of the LHY term. In Fig. 4 we show the temporal
evolution of SS in absence of LHY term for $\Omega=1$ and $k_{L}=4$ with
$g=0.5$ for (a) $v=0$ and (b)$v=0.4$. For $v=0$, soliton does not have a
stable structure and exhibits expansion for a long time. However, at higher
velocity ($v=0.4$), bifurcation of the stripe soliton into two propagating
solitons occurs, which exhibits expansion with time.

After discussing the dynamics of the soliton in the presence of MF term, we
now analyse the dynamical evolution when we consider the LHY term for both QBS
and QSS. In Fig. 5, we illustrate the dynamics of QBS with the LHY correction
term for different initial velocities $v=0,0.5$ and $v=1.0$ at $g=-g_{ab}=0.5$
and $\Omega=k_{L}=1$. For zero velocity soliton exhibits stable nature with no
change in its shape and size with time as depicted in Fig. 5(a). The MF
dynamical instability gets under control due to the presence of
$g_{\text{LHY}}$ interaction, stabilsing the expanding condensate without any
trap. On increasing the velocity to $v=0.5$, the soliton displays breathers-
like expansion and compression behaviour at periodic intervals, as shown in
Fig. 5(b). The frequency of these breathing modes depends on the coupling
parameters, which we shall discuss in detail in the later part of the paper.
Note that Tononi et al. reported similar breathers like soliton [21], which
attained by a slight perturbation in the width of the ground state. As we
increase the initial velocity to $v=1$, the soliton develops periodic breather
oscillations and exhibits bifurcation-like features in which the centre of the
soliton moves in the $\pm x$ directions. Note that the up-spin component moves
in the $-x$ direction while the other component moves in the $+x$ direction.
This feature is noticeable for a finite time $t\lesssim 10$. Beyond this time,
the breathers expand and span the entire space while the centred soliton
exhibits expansion. Apart from this, we also notice that it exhibits bright-
dark soliton for a shorter time, which finally gets converted into multi-
soliton at a longer period [19] as displayed in Fig. 5(c).

Further, we focus on the dynamics of QSS with LHY correction at different
velocities $v=0,0.2,$ and $1.0$ for the parameters $g=-g_{ab}=0.5$,
$\Omega=1$, and $k_{L}=4$. At $v=0$, QSS remains stable and maintains its
shape and size for a longer duration [see Fig. 6(a)]. However with finite
velocity ($v=0.2$) QSS propagates for a while in the $\pm x$ directions,
thereafter, for $t\gtrsim 50$ they start displaying decaying (dilatation)
feature until $t\lesssim 300$. Following this the dilatation slowly gets
diminished and for $t\gtrsim 500$ soliton gets revived [see Fig. 6(b)]. The
whole revival and expansion of the soliton with time is similar to those
obtained for the QBS [see Fig. 5(b)]. For $v>0.5$, the QSS evolves with time
in the respective directions, and after separation, they do not combine again,
as shown for $v=1.0$ in Fig. 6(c).

Figure 8: (a) Variation of breathing frequency ($\omega^{c}_{{}_{bf}}$) at
critical $k_{L}$ with $\Omega$ for fixed velocity $v=0.1$. The other
parameters are the same as those in Fig. 7. The critical breathing frequency
increases and exhibits power law behaviour with $\Omega$ as
$\omega^{c}_{{}_{bf}}\sim\Omega^{0.16}$. The dashed line is drawn as a guide
to the eye to show the power law nature of the $\omega^{c}_{{}_{bf}}$ with
$\Omega$. (b) Variation of breathing frequency $\omega_{bf}$ with initial
velocity $v$ for different $k_{L}$ keeping Rabi frequency fixed to $\Omega=2$.
The $\omega_{bf}$ is small for higher $k_{L}$ for all the velocity.

Next, we analyze the dependence of the breathing frequency for the quantum
soliton as observed in Figs. 5(b) and 6(b) on the SO and Rabi coupling
parameters. We compute the breathing frequency ($\omega_{bf}$) using the time
series of the density oscillation of the soliton. In Fig. 7(a) we show
dependence of $\omega_{bf}$ on the SO couplings ($k_{L}$) for various $\Omega$
when the soliton was assigned the initial velocity as $v=0.1$. Following Ref.
[21], here we consider the initial wavefunction as Gaussian, which is
different than that of the exact solution of the ground state. The breathing
frequency is around $0.025$ for $k_{L}=0$ for all values of $\Omega$. At
$\Omega=0$ the $\omega_{bf}$ remains unchanged upon increasing $k_{L}$. For
finite Rabi coupling frequency ($\Omega=1,2,4$), breather starts appearing at
periodic intervals. For fixed $\Omega$ as $k_{L}^{2}<\Omega$ the breathing
frequency increases upon increasing in $k_{L}$. It attains the maximum at
critical SO coupling ($k^{2}_{L}=\Omega$) where the QBS to QSS phase
transition occurs. For $k_{L}^{2}>\Omega$, $\omega_{bf}$ decreases slowly in
the regime of quantum stripe soliton. At higher $k_{L}$, $\omega_{bf}$ have
the same value as those at $k_{L}\sim 0$ for all $\Omega$. The critical
breathing frequency ($\omega^{c}_{bf}$) increases upon an increase in the Rabi
frequency. Fig. 7(b) illustrates the variation of $\omega_{bf}$ with $k_{L}$
for different velocities, $v=0.1$ and $0.2$, at $g=0.5$. We have considered
two sets of Rabi couplings $\Omega=1$ and $4$ for the analysis. The overall
variation in $\omega_{bf}$ with $k_{L}$ appears to be the same for both the
velocities ($v=0.1$ and $0.2$), except that at lower velocity ($v=0.1$) it has
higher value compared to those at large velocity ($v=0.2$). Interestingly,
$\omega^{c}_{bf}$ remains the same for all the velocities ($v=0.1$ and $0.2$)
as shown in Fig. 7 (b), indicating that the breathing frequency at the
critical point remains unchanged upon increasing the velocity. To better
understand the dependence of the $\omega_{bf}$ on the interaction strengths in
Fig. 7(c) we show its variation with $k_{L}$ for different $g$
($g=0.5,1.0,2.0$) at $\Omega=4$ and $v=0.1$. We observe that the $\omega_{bf}$
for a given $k_{L}$ shows increasing behaviour on increasing the interaction
strength. The critical SO coupling ($k^{c}_{L}$) and corresponding critical
breathing frequency ($\omega^{c}_{bf}$) also increase upon the increase in
$g$.

In Fig. 8(a), we plot the variation of $\omega^{c}_{bf}$ with $\Omega$ for
$v=0.1$. We find that $\omega^{c}_{bf}$ shows power law dependence on $\Omega$
with an exponent, $0.16$. In Fig. 8(b), we plot the variation of $\omega_{bf}$
with initial velocity for fixed Rabi coupling ($\Omega=2$) and for two sets of
SO couplings $k_{L}=1.2$ (solid red line) and $k_{L}=1.8$ (dashed blue line).
Here, we notice that $\omega_{bf}$ decreases upon increasing the initial
velocity for both the $k_{L}$. The fall in the breathing frequency with the
velocity becomes sharper upon the increase in $k_{L}$. However, for a given
initial velocity, the $\omega_{bf}$ is small for $k_{L}=1.8$, which happens to
be in QSS than those for $k_{L}=1.2$ at which we have QBS ground state. We
note that the chemical potential in the QBS regime does not change with
$k_{L}$ until $k_{L}<k^{c}_{L}$, which indicates that the soliton is
energetically stable in the BS regime. On the other hand, when we increase the
Rabi frequency, $\Omega$ to $2$, the critical $k_{L}$ for the transition from
the bright to stripe soliton also increases. Overall we find that the ground
state energy of the quantum soliton gets lower upon increasing $\Omega$ for a
given $k_{L}$, suggesting a better way of obtaining an energetically stable
soliton.

Another notable feature is the transition of the breathing soliton phase into
the moving soliton upon the increase in the velocity as depicted in Figs. 5
and 6. We can understand the instability that leads to the transition from
breathers to the moving soliton by looking at the SO, Rabi and the kinetic
part of the chemical potential. While the kinetic term varies as $v^{2}$, the
SO and Rabi part will have linear and square dependence on the velocity
respectively. Due to this, one witnesses competition between the kinetic
energy, which is positive and the SO and Rabi parts will have negative values
upon changing the velocity. Overall, we expect the domination of the breathing
soliton region by the SO and Rabi coupling contribution, while the moving
soliton region kinetic energy dominates. The above also indicates that the
chemical potential attains minimum at the critical velocity, which is
consistent with our numerical observation.

After analyzing the dynamics of the solitons by imparting an initial velocity,
in the following section, we present different dynamical features arising due
to making an instantaneous quench in the coupling parameters ($\Omega$ and
$k_{L}$).

###  IV.3 Quench dynamics of quantum soliton

To better understand the effect of the SO and Rabi couplings on the overall
dynamics (without giving any initial velocity) of the quantum solitons in this
section, we employ the instantaneous quench of the coupling parameters and
analyze the resultant dynamics. With this protocol, we have obtained a variety
of notable features like the generation of secondary and repulsive solitons,
the dynamical phase transition from the QBS to the QSS and vice versa, etc.

Figure 9: Dynamics of the quantum soliton appearing due to different quenching
protocols at $g=-g_{\uparrow\downarrow}=0.5$ as it was initially prepared for
$\Omega=0,k_{L}=2$. (a)-(c): When the Rabi frequency is quenched as
$\Omega=0\to 1$ at $t=20$, the initial QBS phase changes to QSS, and (d)-(f):
Quenching of Rabi frequency as $\Omega=0\to 5$ at $t=20$ transforms the QBS
phase into repulsive solitons. Figure 10: Dynamics of the quantum soliton
appearing due to different quenching protocols at $t=20$ with
$g=-g_{\uparrow\downarrow}=0.5$ as it was initially prepared for $\Omega=0$,
$k_{L}=2$. (a)-(c): When both Rabi frequency and SO coupling are quenched as
$\Omega=0\to 50,k_{L}=2\to 4$ at $t=20$, which results the transformation of
non-moving QBS soliton into moving multi soliton with large repulsion angle.
(d)-(f): For quenching $\Omega=0\to 50,k_{L}=2\to 8$ at $t=20$, QBS undergoes
bifurcation into moving multi soliton state with small repulsive angle.

In Fig. 9, we show the dynamical evolution of the soliton which generally
arises due to the quench of Rabi coupling frequency by fixing the interaction
strength to $g=-g_{ab}=0.5$. The ground state soliton was initially prepared
for $\Omega=0,k_{L}=2$. As we quench the Rabi frequency from $\Omega=0\to 1$
at finite time ($t\sim 20$) we notice that the QBS gets transformed to the QSS
[see Figs. 9(a-c)]. However, upon employing $\Omega=0\to 5$, the QBS phase of
the soliton gets transformed into a phase that displays expansion and certain
characteristics of the repulsive soliton as shown in Fig. 9(b-d).

Figure 11: Dynamics of the quantum soliton prepared with
$g=-g_{\uparrow\downarrow}=0.5$, $\Omega=0$, and $k_{L}=2$ appearing due to
double quenching protocols. (a)-(c): first quenching is performed at $t=0$
with change in SO coupling ($k_{L}=2\to 0$) and second quenching is done at
$t=20$ when both coupling parameters are changed as $\Omega=0\to 1$ and
$k_{L}=0\to 2$ which results as transforming the moving soliton to the
stationary soliton. (d)-(f): At $t=0$ first quenching is same as those in
(a)-(c), while, at $t=10$ second quenching is done as $\Omega=0\to 5$ that
result transforming the moving soliton into soliton with secondary waves.
Figure 12: Dynamics of the quantum PW soliton prepared with
$g=-g_{\uparrow\downarrow}=0.5$ and $\Omega=2$, $k_{L}=0$ as quenching of both
the coupling parameters are performed. (a)-(c): Quenching is done as
$k_{L}=0\to 2$, $\Omega=2\to 0$ that results the transition of bright soliton
into moving soliton. (d)-(f): With quenching $k_{L}=0\to 2$, $\Omega=2\to 4$
breather soliton gets bifurcated into filament like soliton beyond $t\sim
250$.

In Fig. 10, we show the dynamical evolution of the solitons for two kinds of
quenching protocols: (a-c) when $\Omega=0\to 50$ and $k_{L}=2\to 4$ and (d-f)
when $\Omega=0\to 50,k_{L}=2\to 8$. In the first case, we find the appearance
of multisolitons with interference-like patterns. Following the quenching,
solitons get bifurcated into four that propagate away from each other and
behave like repelling solitons. In the density profile, the outermost solitons
show breathing-like features, and the innermost one inherits pure, stable
soliton-like characteristics, i.e., they do not exhibit any change in shape
and size with time as shown in Fig. 10(a-c). For the second quenching protocol
($\Omega=0\to 50$ and $k_{L}=2\to 8$), again, multisoliton behaviour appears,
which has some differences compared to the previous case. For this case, we
notice that the outermost soliton-breather gets transformed into solitons with
a relatively larger repulsive angle, while the innermost solitons transform
into the stripe solitons and show attraction towards each other as depicted in
Fig. 10(d)-(f). Thus, quenching of the coupling parameters (either one of Rabi
and SO or both simultaneously) leads to secondary solitons.

Figure 13: Dynamics of the quantum stripe soliton prepared with
$g=-g_{\uparrow\downarrow}=0.5$ and $\Omega=1,k_{L}=2$ and the quenching is
performed at finite time. (a)-(c): Quenching is performed as $k_{L}=2\to
0.2,\Omega=1\to 0.5$ at $t=50$ that results the transition from the stripe
soliton to the space-time breather bright soliton. (d)-(f): Quenching is
performed as $\Omega=1\to 2,k_{L}=2$ at $t=20$ resulting the transition from
the stripe soliton to the breathing stripe soliton.

Next we discuss the effect of double quench protocols on the dynamics of the
soliton. In Fig. 11 we show the dynamics of the quantum solitons as we
implement double quenching protocols [38] of the coupling parameters for the
initial ground state prepared at $\Omega=0,k_{L}=2$ for
$g=-g_{\uparrow\downarrow}=0.5$. In Fig. 11(a)-(c), we show the evolution when
first quench is performed by instantaneously switching off the SO coupling,
that is, making the change $k_{L}=2\to 0$ at $t=0$. Subsequently soliton is
allowed to evolve until $t=20$ when another quench was executed by changing
Rabi frequency as $\Omega=0\to 1$ and switching on the SO coupling as
$k_{L}=0\to 2$. The first quenching in which the SO coupling is switched off
initially allows the solitons corresponding to the up and down components to
display repulsive behaviour and move in the space until the second quenching
is done. After the second quench ($\Omega=0\to 1,k_{L}=0\to 2$) at $t=20$, we
notice a switching off the soliton position from $x=0$ to some finite value
which further remains invariant with time. Overall we see that the double
quench protocol provides more stability to the solitons. This particular
feature is one of the promising potential application for quantum information
and quantum computing, which based on mixing and demixing of qubits [36, 37].
As we perform the first quenching same as previous case while in second quench
only change in Rabi frequency is executed with $\Omega=0\to 5$ at $t=10$ we
find the appearance of moving BS with generation of the second harmonics [see
Fig. 11(d)-(f)].

Next we consider $\Omega=2,k_{L}=0$ with $v=0$ which shows only the stable BS.
As we quench $k_{L}=0\to 2$ and $\Omega=2\to 0$ the soliton repel each other
[see Fig. 12(a)-(c)]. Upon quenching the Rabi coupling frequency from
$\Omega=2\to 4$ we find that the soliton gets transformed into filaments after
a finite time ($t\gtrsim 150$), which also appears to repel each other at
($t\gtrsim 250$) as depicted in Fig. 12 (d)-(f).

Figure 14: Dynamics of quantum stripe soliton prepared with
$g=-g_{\uparrow\downarrow}=1.0$ and $\Omega=k_{L}=2$ as the quenching is
performed on the Rabi frequency at $t=10$. (a) Soliton gets bifurcated into
repulsive and attractive solitons as Rabi coupling is quenched as (a)
$\Omega=2\to 1$ and (b) $\Omega=2\to 3$.

After discussing different dynamical behaviour arising in the quantum bright
soliton due to the quench of the coupling parameter in what follows, we
present the related dynamics for the quantum stripe soliton. In Fig. 13, we
depict the dynamics that arise due to the quenching of both coupling
parameters when the ground state is quantum stripe soliton nature.

In Fig. 13(a)-(c) we show the dynamics of soliton when quenching of both
coupling parameters are performed at $t=50$. The quenching protocol for SO
coupling is implemented as $k_{L}=2\to 0.2$ while for Rabi coupling frequency
it is $\Omega=1\to 0.5$. We noticed that as a result of the quenching, stripes
solitons get bifurcated into two parts wherein one shows weak breather-
soliton, while the other exhibits stronger behaviour. Further decrease of Rabi
coupling frequency leads to the disappearance of breathers and the generation
of repulsive solitons. In Fig. 13(d)-(f) we depict the dynamics of soliton
when quenching of Rabi coupling frequency is executed from $\Omega=1\to 2$ at
$t=20$. The quenching generates the stripe soliton breathers. Similar kinds of
breathing solitons have been observed upon quenching the interaction strength
so that the mean field contribution becomes quite small in the binary mixture
of the condensates [20].

Figure 15: (a) The total density plots depicting the inelastic collision
dynamics between the quantum soliton for $g=-g_{\uparrow\downarrow}=0.5$ and
$\Omega=k_{L}=0$ as the velocity given to the up and down component is in the
direction of $-x$ and $x$ direction respectively, with magnitude $v=0.2$. (b)
Evolution of the different energies, $E_{ki}$, $E_{lhy}$, $E_{mf}$ and total
energy ($E_{T}$) with time for collision as shown in (a). There is a
significant increase in the $E_{lhy}$ and total energy ($E_{T}$) after the
collision, indicating the inelastic nature of the same. (c) plots of the total
density depicting the elastic collision dynamics of quantum soliton for
$v=0.9$. (d) Evolution of the different energies with time for the collision
as shown in (c). Before and after the collision, all the energy remains the
same, indicating the elastic collision nature.

In Fig. 14 we show the dynamics that solely arises due to the quenching of the
Rabi coupling for the ground states initially prepared for
$g=-g_{\uparrow\downarrow}=1.0$ and $\Omega=k_{L}=2$. We quench the parameters
as (a) $\Omega=2\to 1$ and in (b) $\Omega=2\to 3$ at $t=10$. Here we notice
that the stripe soliton exhibits expansion upon the quenching of the
parameters. We find that the breathing oscillation frequency increases on
increasing the interaction strength, as shown in Fig. 7(c). Further we
consider another stripe soliton with $\Omega=k_{L}=2$ in which the Rabi
coupling is quenched as $\Omega=2\to 1$. For this, we find that the frequency
of the striped soliton breather becomes $\omega_{bf}\sim 0.08$ upon increasing
the interaction strength. When we change $\Omega=2\to 3$, we obtain that the
initial profile, which has two maximum peaks and two sidelobes, gets
transformed into the multiple side lobes after the quench. Apart from this,
the innermost soliton appears to repel each other after a finite time as they
approach each other. A similar kind of dynamical behaviour gets repeated over
time.

Figure 16: Collisional dynamics of the quantum solitons prepared with
$g=-g_{\uparrow\downarrow}=0.5$ and $\Omega=1$, $k_{L}=0$ and at $t=0$ with
velocity given to the up and down component with opposite value and magnitude
as (a)-(c): $v=0$ and (d)-(f): $v=1$. For $v=0$, the soliton displays spatio-
temporal breathers with propagation along a line. However, for $v=1$ soliton
exhibits space-time breathers with quasi-elastic collision around $t\sim 20$.
Figure 17: Collisional dynamics of the quantum solitons prepared with
$g=-g_{\uparrow\downarrow}=0.5$ $\Omega=0$ and $k_{L}=0.5$ and the individual
components are given equal and opposite velocity ($v=0.5$) at $t=0$. Also
quenching on SO coupling parameter is performed as (a)-(c): $k_{L}=0.5\to 0.1$
and (d)-(f): $k_{L}=0.5\to 2$. Upon quenching the SO coupling, the solitons
undergo inelastic collision in (a)-(c), while solitons repel each other for
(d)-(f). Figure 18: Collisional dynamics of the quantum solitons prepared with
the same parameters and given same velocity at $t=0$ as those for Fig. 17. At
$t=0$ the Rabi coupling is quenched as (a)-(c): $\Omega=0\to 0.1$ (d)-(f):
$\Omega=0\to 0.8$. Upon quenching, Rabi coupling quantum soliton undergoes
inelastic collision and displays an interference pattern (a)-(c). However, the
soliton displays quasi-elastic collision with secondary solitons generation in
(d)-(f).

###  IV.4 Collisional dynamics of quantum soliton

After studying the quenching dynamics of the quantum soliton by perturbing it
with an initial velocity, in this section, we investigate how the solitons
undergo different kinds of collisions that affect its overall stability
depending upon the magnitude of the initial velocity. In particular, we find
elastic, inelastic, repulsive and space-time breather soliton against
collisions.

In absence of the trap, the quantum solitons are set in motion by giving an
initial velocity $\mp v$ to the stationary ground state $(\psi_{1},\psi_{2})$
with a multiplier $\exp{(\mp ivx)}$ respectively [19]. Here we start for the
case when the solitons are positioned at $\pm 50$ for $g=-g_{ab}=0.5$ and
$\Omega=k_{L}=0$. For a weak velocity ($v=0.2$) the solitons undergo inelastic
collision as depicted in Fig. 15(a). We find that the densities before the
collision and after the collision are not the same. This indicates that after
the collision, the solitons are unable to retain their shape and size same as
those before the collision. As the initial positions of the soliton are at
$\pm 50$ the expected time for collision is $t=x/v=250$ which can be clearly
seen in Fig. 15(a). We find that after the collision ($t\sim 400$), the
solitons start exhibiting the expansion, which makes them quite different in
shape and size than they were before the collision. In order to understand the
inelastic collision in better way in Fig. 15(b) we plot the time evolution of
the different energy contribution, like, kinetic energy ($E_{ki}$), mean-field
energy ($E_{mf}$), LHY energy ($E_{lhy}$) and the total energy ($E_{T}$). The
explicit form of energies is provided in the Eq. (13) of Appendix A. We find
that the mean-field energy, kinetic, LHY and total energy after the collision
does not remain same as they are before the collisions which is consistent
with the inelastic nature of the collision observed in the classical case.

For large velocity ($v=0.9$), we observe the elastic collision with no change
in the shape of the soliton after the collision. In Fig. 15 (c) we show the
total density corresponding to the soliton when each component has been given
an initial equal velocity ($v=0.9$) but opposite direction with their initial
positions at $\pm 50$ for $g=-g_{ab}=0.5$ and $\Omega=k_{L}=0$. For this, the
expected collision time will be at $t=55.5$. We find that the densities before
and after the collision are almost similar. In Fig. 15 (d), we plot the time
evolution of the different energy contributions to confirm the elastic nature
of the collision. It is observed that the mean-field energy, kinetic, LHY and
total energy remain the same after the collisions what they were before the
collision, which are consistent with the elastic nature of the collision. This
behaviour is consistent with the experimental observations of the collision
between the quantum droplet in three dimension [14].

So far, we have analyzed the collisional dynamics of the binary BECs with LHY
correction and in the absence of Rabi and SO coupling parameters. However, as
we consider it for finite SO and Rabi couplings, we find that the solitons get
degenerated and thus, it could not provide the relevant collision dynamics.
Considering what follows, we present the collision dynamics by keeping either
SO or Rabi coupling frequency to be finite. Here as we obtain the spatially
separated soliton, we employ the quenching of velocity and the coupling
parameter, which was zero while preparing the states. To study the collisional
dynamics with finite coupling parameters, we prepare the initial ground state
with $\Omega=1$ and $k_{L}=0$. In Fig. 16(a)-(c) we show the dynamics with
zero velocity. For this case, we find each spin component degenerates in the
$\pm x$ direction and further exhibits the space-time breathers, which
experience expansion after the finite time ($t\sim 40$). As we pass on the
finite velocity, the solitons start interacting with each other. As for an
example, for $v=1$ the soliton starts moving towards each other that results
in the elastic collision between them at $t\sim 20$ [see Fig. 16(d)-(f)]. Note
that here, after the collision, there is a small change in the density, which
is different from the collision behaviour observed in Fig. 15(c). Therefore,
we refer to this as the quasi-elastic collision.

Next, we present the collisional dynamics by preparing the ground states with
finite SO coupling and zero Rabi coupling frequency. We consider the system
with $\Omega=0$, $k_{L}=0.5$ with $g=0.5$ and finite velocity $v=0.5$. Once we
prepare the ground state, we analyse the collisional dynamics by quenching the
SO coupling for two different situations. When $v>k_{L}/2$ non-degenerated
soliton with decaying or expansion feature undergo inelastic collision as
depicted in Fig. 17(a)-(c). However, for $v<k_{L}/2$ we observe that the
soliton appears to repel each other upon progression of time [see Fig.
17(d)-(f)].

Further, we analyse the collision dynamics by quenching the Rabi coupling.
Like in the previous case, we study the dynamics of the solitons for two
different situations. When $\Omega<k_{L}^{2}$ non-degenerated soliton undergo
inelastic collision. They appear to show the stripe fringes as depicted in
Fig. 18(a)-(c). For $\Omega>k_{L}^{2}$, we find a collision of bright soliton
is quasi-elastic with the generation of interference stripe fringes. We also
witness the generation of secondary solitons, as shown in Fig. 18(d)-(f).

##  V Summary and Conclusions

Using the mean-field model along with the LHY correction, we have numerically
investigated the structure and dynamics of the ground state of the self-bound
state in 1D spin-orbit coupled BECs. Depending upon the nature of Rabi and SO
couplings, the ground state is either quantum bright or stripe soliton. We
deduced an analytical solution of the quantum bright soliton without Rabi
coupling. We further studied the dynamics of the soliton using the three
protocols, namely, by giving the initial velocity to the soliton, quenching
the SO and Rabi coupling parameters and allowing the component to collide by
attributing an equal and opposite initial velocity. We found that the velocity
perturbations generate breathing-like soliton. The breathing frequency
increases upon the increase in the SO coupling for a given Rabi coupling
frequency attains a maximum at the SO coupling, where the phase transition
from the quantum bright to stripe soliton takes place. The magnitude of
breathing frequency at critical SO coupling remains independent of the initial
velocity and exhibits a power law dependence on the Rabi coupling frequency
with an exponent $\sim 0.16$. We found that the critical breathing frequency
increases upon increasing the interaction strengths. Using the contribution of
the kinetic, SO, and Rabi energy terms, we find that while the domination of
the breathing soliton phase by the SO and Rabi coupling, the moving soliton
phase is the result of the kinetic energy part. This particular feature leads
the attainment of the minima of the chemical potential at the critical
velocity where change in the phase takes place. We have realized the presence
of several dynamical phase transitions, like, QBS to QSS and the multisoliton
behaviour depending upon the quenching protocol of the Rabi coupling
parameters. By quenching both Rabi and the SO couplings, we can control the
direction/angle of the inner- and outer soliton. The quenching of Rabi and SO
couplings facilitate several interesting dynamical phases like a repulsive
soliton, space-time breathers, filamentation, etc.

We have also analyzed the collision dynamics of solitons. Depending upon the
velocity of the soliton, we observed the presence of elastic and inelastic
collisions. An inelastic collision occurs for low velocities, while it
exhibits elastic collision at high velocities consistent with earlier
experimental observations [14]. We have also complemented the collision
dynamics by analyzing the nature of different energy terms like kinetic
energy, mean-field energy, and LHY. Also, we analyzed the collision dynamics
of the solitons by quenching either the velocity or the SO/Rabi coupling. For
the former, we found the presence of quasi-elastic collision, while for the
latter, we realize inelastic collision.

In this paper, we have considered the dynamical evolution of the quantum
solitons for the situation when the mean-field contribution is negligible. It
would be interesting to explore the dynamics in the similar line presented in
this work in the presence of significant mean-field contribution along with
the LHY correction in the quantum soliton.

###### Acknowledgements.

We thank Cesar Cabrera for his comments on the initial draft of our
manuscript. S.G. would like to acknowledge the financial support from the
University Grants Commission - Council of Scientific and Industrial Research
(UGC-CSIR), India. R.R. acknowledge support from the Ministry of Science and
Technology (MOST), Taiwan, under Grant No. MOST-111-2119-M-001-002. P.K.M.
acknowledges the Department of Science and Technology - Science and
Engineering Research Board (DST-SERB) India for the financial support through
Project No. ECR/2017/002639. The work of P.M. is supported by DST-SERB under
Grant No. CRG/2019/004059, DST-FIST under Grant No. SR/FST/PSI-204/2015(C),
and MoE RUSA 2.0 (Physical Sciences).

##  Appendix A Calculation of energy in the SO coupled BECs with LHY
correction

In this appendix, we provide the detailed steps to obtain the total energy of
the SO coupled BECs with LHY correction. The stationary state solution is
given by

$\displaystyle\psi_{j}(x,y)=\big{(}\psi_{jR}+\mathrm{i}\psi_{jI}\big{)}e^{-\mathrm{i}\mu_{j}t}$
(8)

where $j\in\\{\uparrow,\downarrow\\}$, $\psi_{jR}$ and $\psi_{jI}$ are the
real and imaginary parts of the stationary wavefunction respectively, and
$\mu_{\uparrow,\downarrow}$ are the chemical potential of the spin-up and down
components respectively [39]. As we insert this solution in the Eqs. (1) we
get

$\displaystyle\mu_{\uparrow}\psi_{\uparrow R}=$
$\displaystyle\left[-\frac{1}{2}\frac{\partial^{2}}{\partial
x^{2}}+g\lvert\psi_{\uparrow}\rvert^{2}+g_{\uparrow\downarrow}\lvert\psi_{\downarrow}\rvert^{2}-\frac{g^{3/2}}{\pi}\sqrt{\lvert\psi_{\uparrow}\rvert^{2}+\lvert\psi_{\downarrow}\rvert^{2}}\right]\psi_{\uparrow
R}+k_{L}\left(\frac{\partial\psi_{\uparrow I}}{\partial
x}\right)+\Omega\psi_{\downarrow R},$ (9a)
$\displaystyle\mu_{\downarrow}\psi_{\downarrow R}=$
$\displaystyle\left[-\frac{1}{2}\frac{\partial^{2}}{\partial
x^{2}}+g\lvert\psi_{\downarrow}\rvert^{2}+g_{\downarrow\uparrow}\lvert\psi_{\uparrow}\rvert^{2}-\frac{g^{3/2}}{\pi}\sqrt{\lvert\psi_{\uparrow}\rvert^{2}+\lvert\psi_{\downarrow}\rvert^{2}}\right]\psi_{\downarrow
R}-k_{L}\left(\frac{\partial\psi_{\downarrow I}}{\partial
x}\right)+\Omega\psi_{\uparrow R},$ (9b) and for the imaginary parts we have
$\displaystyle\mu_{\uparrow}\psi_{\uparrow I}=$
$\displaystyle\left[-\frac{1}{2}\frac{\partial^{2}}{\partial
x^{2}}+g\lvert\psi_{\uparrow}\rvert^{2}+g_{\uparrow\downarrow}\lvert\psi_{\downarrow}\rvert^{2}-\frac{g^{3/2}}{\pi}\sqrt{\lvert\psi_{\uparrow}\rvert^{2}+\lvert\psi_{\downarrow}\rvert^{2}}\right]\psi_{\uparrow
I}-k_{L}\left(\frac{\partial\psi_{\uparrow R}}{\partial
x}\right)+\Omega\psi_{\downarrow I},$ (9c)
$\displaystyle\mu_{\downarrow}\psi_{\downarrow I}=$
$\displaystyle\left[-\frac{1}{2}\frac{\partial^{2}}{\partial
x^{2}}+g\lvert\psi_{\downarrow}\rvert^{2}+g_{\downarrow\uparrow}\lvert\psi_{\downarrow}\rvert^{2}-\frac{g^{3/2}}{\pi}\sqrt{\lvert\psi_{\uparrow}\rvert^{2}+\lvert\psi_{\downarrow}\rvert^{2}}\right]\psi_{\uparrow
I}+k_{L}\left(\frac{\partial\psi_{\downarrow R}}{\partial
x}\right)+\Omega\psi_{\uparrow I},$ (9d)

where $\lvert\psi_{\uparrow}\rvert^{2}=\psi_{\uparrow R}^{2}+\psi_{\uparrow
I}^{2}$ and $\lvert\psi_{\downarrow}\rvert^{2}=\psi_{\downarrow
R}^{2}+\psi_{\downarrow I}^{2}$. Multiplying Eq. (9a) with $\psi_{\uparrow R}$
and Eq. (9b) with $\psi_{\downarrow R}$, and integrating we get

$\displaystyle\mu_{\uparrow}\int dx\,\psi_{\uparrow R}^{2}=$
$\displaystyle\int dx\,\psi_{\uparrow
R}\left\\{\left[-\frac{1}{2}\frac{\partial^{2}}{\partial
x^{2}}+g\lvert\psi_{\uparrow}\rvert^{2}+g_{\uparrow\downarrow}\lvert\psi_{\downarrow}\rvert^{2}-\frac{g^{3/2}}{\pi}\sqrt{\lvert\psi_{\uparrow}\rvert^{2}+\lvert\psi_{\downarrow}\rvert^{2}}\right]\psi_{\uparrow
R}+k_{L}\left(\frac{\partial\psi_{\uparrow I}}{\partial
x}\right)+\Omega\psi_{\downarrow R}\right\\},$ (10a)
$\displaystyle\mu_{\downarrow}\int dx\,\psi_{\downarrow R}^{2}=$
$\displaystyle\int dx\,\psi_{\downarrow
R}\left\\{\left[-\frac{1}{2}\frac{\partial^{2}}{\partial
x^{2}}+g\lvert\psi_{\downarrow}\rvert^{2}+g_{\downarrow\uparrow}\lvert\psi_{\uparrow}\rvert^{2}-\frac{g^{3/2}}{\pi}\sqrt{\lvert\psi_{\uparrow}\rvert^{2}+\lvert\psi_{\downarrow}\rvert^{2}}\right]\psi_{\downarrow
R}-k_{L}\left(\frac{\partial\psi_{\downarrow I}}{\partial
x}\right)+\Omega\psi_{\uparrow R}\right\\},$ (10b)

Similarly, from the imaginary part equations (Eqs. 9c and 9d) we obtain

$\displaystyle\mu_{\uparrow}\int dx\,\psi_{\uparrow I}^{2}=$
$\displaystyle\int dx\,\psi_{\uparrow
I}\left\\{\left[-\frac{1}{2}\frac{\partial^{2}}{\partial
x^{2}}+g\lvert\psi_{\uparrow}\rvert^{2}+g_{\uparrow\downarrow}\lvert\psi_{\downarrow}\rvert^{2}-\frac{g^{3/2}}{\pi}\sqrt{\lvert\psi_{\uparrow}\rvert^{2}+\lvert\psi_{\downarrow}\rvert^{2}}\right]\psi_{\uparrow
I}-k_{L}\left(\frac{\partial\psi_{\uparrow R}}{\partial
x}\right)+\Omega\psi_{\downarrow I}\right\\},$ (11a)
$\displaystyle\mu_{\downarrow}\int dx\,\psi_{\downarrow I}^{2}=$
$\displaystyle\int dx\,\psi_{\downarrow
I}\left\\{\left[-\frac{1}{2}\frac{\partial^{2}}{\partial
x^{2}}+g\lvert\psi_{\downarrow}\rvert^{2}+g_{\downarrow\uparrow}\lvert\psi_{\uparrow}\rvert^{2}-\frac{g^{3/2}}{\pi}\sqrt{\lvert\psi_{\uparrow}\rvert^{2}+\lvert\psi_{\downarrow}\rvert^{2}}\right]\psi_{\downarrow
I}+k_{L}\left(\frac{\partial\psi_{\downarrow R}}{\partial
x}\right)+\Omega\psi_{\uparrow I}\right\\},$ (11b)

Upon rearranging the above equations (Eqs. 10b-11b) we obtain

$\displaystyle\mu_{\uparrow}=$
$\displaystyle\frac{1}{\int\lvert\psi_{\uparrow}\rvert^{2}dx}\int\left[\frac{1}{2}\left|\frac{\partial\psi_{\uparrow}}{\partial
x}\right|^{2}+\left\\{g\lvert\psi_{\uparrow}\rvert^{2}+g_{\uparrow\downarrow}\lvert\psi_{\downarrow}\rvert^{2}-\frac{g^{3/2}}{\pi}\sqrt{\lvert\psi_{\uparrow}\rvert^{2}+\lvert\psi_{\downarrow}\rvert^{2}}\right\\}\lvert\psi_{\uparrow}\rvert^{2}\right]dx$
$\displaystyle+\frac{1}{\int|\psi_{\uparrow}\rvert^{2}dx}\int\left[k_{L}\left(\frac{\partial\psi_{\uparrow
I}}{\partial x}\right)+\Omega\psi_{\downarrow R}\right]\psi_{\uparrow
R}dx+\frac{1}{\int|\psi_{\uparrow}\rvert^{2}dx}\int\left[-k_{L}\left(\frac{\partial\psi_{\uparrow
R}}{\partial x}\right)+\Omega\psi_{\downarrow I}\right]\psi_{\uparrow I}dx,$
(12a) $\displaystyle\mu_{\downarrow}=$
$\displaystyle\frac{1}{\int\lvert\psi_{\downarrow}\rvert^{2}dx}\int\left[\frac{1}{2}\left|\frac{\partial\psi_{\downarrow}}{\partial
x}\right|^{2}+\left\\{g\lvert\psi_{\downarrow}\rvert^{2}+g_{\downarrow\uparrow}\lvert\psi_{\downarrow}\rvert^{2}-\frac{g^{3/2}}{\pi}\sqrt{\lvert\psi_{\uparrow}\rvert^{2}+\lvert\psi_{\downarrow}\rvert^{2}}\right\\}\lvert\psi_{\downarrow}\rvert^{2}\right]dx$
$\displaystyle+\frac{1}{\int|\psi_{\downarrow}\rvert^{2}dx}\int\left[-k_{L}\left(\frac{\partial\psi_{\downarrow
I}}{\partial x}\right)+\Omega\psi_{\uparrow R}\right]\psi_{\downarrow
R}dx+\frac{1}{\int\lvert\psi_{\downarrow}\rvert^{2}dx}\int\left[k_{L}\left(\frac{\partial\psi_{\downarrow
R}}{\partial x}\right)+\Omega\psi_{\uparrow I}\right]\psi_{\downarrow I}dx,$
(12b)

From the above equations we obtain the different contribution in the energy as

$\displaystyle E_{ki}=$
$\displaystyle\int\left[\frac{1}{2\int\lvert\psi_{\uparrow}\rvert^{2}dx}\left|\frac{\partial\psi_{\uparrow}}{\partial
x}\right|^{2}+\frac{1}{2\int\lvert\psi_{\downarrow}\rvert^{2}dx}\left|\frac{\partial\psi_{\downarrow}}{\partial
x}\right|^{2}\right]dx$ (13a) $\displaystyle E_{mf}=$
$\displaystyle\int\left[\frac{1}{\int\lvert\psi_{\uparrow}\rvert^{2}dx}\left[\frac{g}{2}\lvert\psi_{\uparrow}\rvert^{2}+\frac{g_{\uparrow\downarrow}}{2}\lvert\psi_{\downarrow}\rvert^{2}\right]\lvert\psi_{\uparrow}\rvert^{2}+\frac{1}{\int\lvert\psi_{\downarrow}\rvert^{2}dx}\left[\frac{g_{\downarrow\uparrow}}{2}\lvert\psi_{\uparrow}\rvert^{2}+\frac{g}{2}\lvert\psi_{\downarrow}\rvert^{2}\right]|\psi_{\downarrow}\rvert^{2}\right]dx$
(13b) $\displaystyle E_{lhy}=$
$\displaystyle\int\left[-\frac{1}{\int\lvert\psi_{\uparrow}\rvert^{2}dx}\left[\frac{2}{3\pi}g^{3/2}\sqrt{\lvert\psi_{\uparrow}\rvert^{2}+\lvert\psi_{\downarrow}\rvert^{2}}\right]\lvert\psi_{\uparrow}\rvert^{2}-\frac{1}{\int\lvert\psi_{\downarrow}\rvert^{2}dx}\left[\frac{2}{3\pi}g^{3/2}\sqrt{\lvert\psi_{\uparrow}\rvert^{2}+\lvert\psi_{\downarrow}\rvert^{2}}\right]\lvert\psi_{\downarrow}\rvert^{2}\right]dx$
(13c) $\displaystyle E_{SO}=$
$\displaystyle\int\left\\{\frac{1}{\int\lvert\psi_{\uparrow}\rvert^{2}dx}\left[\left(k_{L}\frac{\partial\psi_{\uparrow
I}}{\partial_{x}}\right)\psi_{\uparrow
R}+\left(-k_{L}\frac{\partial\psi_{\uparrow
R}}{\partial_{x}}\right)\psi_{\uparrow
I}\right]+\frac{1}{\int\lvert\psi_{\downarrow}\rvert^{2}dx}\left[\left(-k_{L}\frac{\partial\psi_{\downarrow
I}}{\partial_{x}}\right)\psi_{\downarrow
R}+\left(k_{L}\frac{\partial\psi_{\downarrow
R}}{\partial_{x}}\right)\psi_{\downarrow I}\right]\right\\}dx$ (13d)
$\displaystyle E_{Rabi}=$
$\displaystyle\int\left\\{\frac{1}{\int\lvert\psi_{\uparrow}\rvert^{2}dx}\left[\left(\Omega\psi_{\downarrow
R}\right)\psi_{\uparrow R}+\left(\Omega\psi_{\downarrow
I}\right)\psi_{\uparrow
I}\right]+\frac{1}{\int\lvert\psi_{\downarrow}\rvert^{2}dx}\left[\left(\Omega\psi_{\uparrow
R}\right)\psi_{\downarrow R}+\left(\Omega\psi_{\uparrow
I}\right)\psi_{\downarrow I}\right]\right\\}dx$ (13e)

where $E_{ki}$ represents the kinetic energy, $E_{mf}$ mean-field
contribution, $E_{lhy}$ is the energy due to the LHY correction, and $E_{SO}$
is the contribution due to the SO coupling. The total energy
$E_{T}=E_{ki}+E_{mf}+E_{lhy}+E_{SO}$.

## References

  * [1] Z-H. Luo, W. Pang, B. Liu, Y-Y. Li, and B. A. Malomed, Front. Phys. 16(3), 32201 (2021).
  * [2] T. D. Lee, K. Huang, and C. N. Yang, Phys. Rev. 106, (1957).
  * [3] S. B. Papp, J. M. Pino, R. J. Wild, S. Ronen, C. E. Wieman, D. S. Jin, and E. A. Cornell, Phys. Rev. Lett. 101, 135301 (2008).
  * [4] N. Navon, S. Piatecki, K. Günter, B. Rem, T. C. Nguyen, F. Chevy, W. Krauth, and C. Salomon, Phys. Rev. Lett. 107, 135301 (2011).
  * [5] T. G. Skov, M. G. Skou, N. B. Jørgensen, and J. J. Arlt, Phys. Rev. Lett. 126 230404 (2021).
  * [6] D. S. Petrov, Phys. Rev. Lett. 115, 155302 (2015).
  * [7] I. Ferrier-Barbut, H. H. Kadau, M. Schmitt, M. Wenzel, and T. Pfau, Phys. Rev. Lett. 116, 215301 (2016).
  * [8] M. Schmitt, M. Wenzel, F. Böttcher, I. Ferrier-Barbut, and T. Pfau, Nature 539, 259 (2016).
  * [9] I. Ferrier-Barbut, M. Schmitt, M.Wenzel, H. Kadau, and T. Pfau, J. Phys. B 49, 214004 (2016).
  * [10] L. Chomaz, S. Baier, D. Petter, M. J. Mark, F. Wächtler, L. Santos, and F. Ferlaino, Phys. Rev. X 6, 041039 (2016).
  * [11] C. R. Cabrera, L. Tanzi, J. Sanz, B. Naylor, P. Thomas, P. Cheiney, and L. Tarruell, Science 359, 301 (2018).
  * [12] P. Cheiney, C. R. Cqabrera, J. Sanz, B. Naylor, L. Tanzi and L. Tarruell, Phys. Rev. Lett. 120, 135301 (2018).
  * [13] G. Semeghini, G. Ferioli, L. Masi, C. Mazzinghi, L. Wolswijk, F. Minardi, M. Modugno, G. Modugno, M. Inguscio, and M. Fattori, Phys. Rev. Lett. 120, 235301 (2018).
  * [14] G. Ferioli, G. Semeghini, L. Masi, G. Giusti, G. Modugno, M. Inguscio, A. Gallemi, A. Recati, and M. Fattori, Phys. Rev. Lett. 122, 090401 (2019).
  * [15] D. Rakshit, T. Karpiuk, M. Brewczyk, and M. Gajda, SciPost Phys. 6, 079 (2019).
  * [16] J. C. Smith, D. Baillie, and P. B. Blakie, Phys. Rev. Lett. 126, 025302 (2021).
  * [17] D. S. Petrov and G. E. Astrakharchik, Phys. Rev. Lett. 117, 100401 (2016).
  * [18] A. Frolian, C. S. Chisholm, E. Neri, C. R. Cabrera, R. Ramos, A. Celi, L. Tarruell, Nature 608, 293 (2022).
  * [19] G. E. Astrakharchik and B. A. Malomed, Phys. Rev. A 98, 013631 (2018).
  * [20] S. I. Mistakidis, T. Mithun, P. G. Kevrekidis, H. R. Sadeghpour, and P. Schmelcher, Phys. Rev. Res. 3, 043128 (2021).
  * [21] A. Tononi, Y. Wang, and L. Salasnich, Phys. Rev. A 99, 063618 (2019).
  * [22] S. Sahu and D. Majumdar, J. Phys. B: At. Mol. Opt. Phys. 53, 095301 (2020).
  * [23] E. Chiquillo, Phys. Rev. A 97, 013614 (2018).
  * [24] M. R. Pathak and A. Nath, Sci. Rep. 12, 6904 (2022).
  * [25] Y. Li, Z. Luo, Y. Liu, Z. Chen, C. Huang, S. Fu, H. Tan, and B. A. Malomed, New J. Phys. 19, 113043 (2017).
  * [26] P. Cheiney, C. R. Cabrera, J. Sanz, B. Naylor, L. Tanzi, and L. Tarruell, Phys. Rev. Lett. 120, 135301 (2018).
  * [27] K.-L. Pan, C. K. Law, K. Chung, and B. Zhou, J. Appl. Phys. 103, 064901 (2008).
  * [28] L. E. Young-S and S. K. Adhikari, Commn. Non. Sci. and Num. Simu. 106, 106094 (2022).
  * [29] Y. J. Lin, K. Jimenez-Garcia, and I. Spielman, Nature 471, 83 (2011).
  * [30] P. Muruganandam, and S. K. Adhikari, Comput. Phys. Commun. 180, 1888 (2021).
  * [31] L. E. Young-S., D. Vudragović, P. Muruganandam, S. K. Adhikari, and A. Balaž, Comput. Phys. Commun. 204, 209 (2016).
  * [32] R Ravisankar, D Vudragović, P Muruganandam, A Balaž, S. K. Adhikari, Comput. Phys. Commun. 259, 107657 (2021).
  * [33] R. Ravisankar, T. Sriraman, and P. Muruganandam, AIP Conf. Proc. 2265, 030022 (2020).
  * [34] J. Denschlag, J. E. Simsarian, D. L. Feder, C. W. Clark, L. A. Collins, J. Cubizolles, L. Deng et al., Science 287, 5450 (2000).
  * [35] R. Ravisankar, T. Sriraman, L. Salasnich, and P. Muruganandam, J. Phys. B: At. Mol. Opt. Phys. 53, 195301 (2020).
  * [36] S. Mardonov, M. Modugno, and E. Y, Sherman, J. Phys. B: At. Mol. Phys. 48, 115302 (2015).
  * [37] R. Ravisankar, T. Sriraman, R. K. Kumar, P. Muruganandam and P. K. Mishra, J. Phys. B: At. Mol. Opt. Phys. 54, 225301 (2021).
  * [38] C-H. Li, et al, Nat. Commun. 10, 375 (2019).
  * [39] R. Ravisankar, H. Fabaralli, A. Gammal, P. Muruganandam, and P. K. Mishra, Phys. Rev. A 104, 053315 (2021).