# Effect of Rashba spin-orbit and Rabi couplings on the excitation spectrum of
binary Bose-Einstein condensates

Rajamanickam Ravisankar  Department of Physics, Indian Institute of
Technology, Guwahati 781039, Assam, India  Department of Physics,
Bharathidasan University, Tiruchirappalli 620024, Tamilnadu, India  Henrique
Fabrelli  Instituto de Física, Universidade de São Paulo, 05508-090 São Paulo,
Brazil  Arnaldo Gammal  Instituto de Física, Universidade de São Paulo,
05508-090 São Paulo, Brazil  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 collective excitation spectrum analysis of binary Bose-Einstein
condensates (BECs) with spin-orbit (SO) and Rabi couplings in a quasi-two-
dimensional system. In particular, we investigate the role of SO and Rabi
coupling strengths in determining the dynamical stability of the coupled BECs
using Bogoliubov-de Gennes (BdG) theory. Using the eigenergy of BdG spectrum,
we confirm the existence of phonon, roton, and maxon modes with weak repulsive
intra- and inter-species contact interactions. The depth of the minimum
corresponding to the roton mode depends strongly on the coupling strength. We
find that the increase of the SO coupling leads to instability, while the
increase in the Rabi coupling stabilizes the system. Also the eigenvectors of
BdG spectrum indicates the presence of density like mode in the stable regime
and spin like modes in unstable regimes. A phase diagram demonstrating the
stability regime in the plane of SO and Rabi coupling strengths is obtained.
Finally, we complement the observation of the excitation spectrum with the
direct numerical simulation results of coupled Gross-Pitaevskii equations.

##  I Introduction

The experimental realization of Bose-Einstein condensates (BECs) in dilute
atomic gases has triggered immense interest in the physics of ultracold matter
Anderson  _et al._ (1995); Davis  _et al._ (1995); Bradley  _et al._ (1995).
Since then, BECs has become an excellent system for manipulating many of the
macroscopic phenomena through the controlled environment in the quantum
regime. Last few decades have seen an upsurge in the research of BEC in
particular understanding the fundamental and dynamical aspects of solitons,
its behaviour under the optical lattice and disordered potentials trap,
superfluid-Mott insulator phase transition, presence of localization, dipolar
and spin-orbit (SO) coupled BECs, etc. Gerton  _et al._ (2000); Greiner  _et
al._ (2002); Morsch and Oberthaler (2006); Lewenstein  _et al._ (2007); Roati
_et al._ (2008); Chin  _et al._ (2010); Griesmaier  _et al._ (2005); Lu  _et
al._ (2011); Aikawa  _et al._ (2012); Lin  _et al._ (2011); Galitski and
Spielman (2013). In this paper we investigate the effect of the Rabi and
Rashba spin-orbit coupling on the dynamical stability of the binary BEC
system.

In BECs, the internal atomic states can be manipulated to produce quite novel
systems like binary and multi-component condensates. These binary or multi-
component BECs consist of different isotopes or hyperfine states of the same
or the different atomic species, which are coherently coupled by the external
fields Ho (1998); Papp  _et al._ (2008). This facilitates the testing ground
of plethora of macroscopic quantum many-body phenomena such as quantum
turbulence Takeuchi  _et al._ (2010), quantum phase transitions Sabbatini  _et
al._ (2011), quantized vortices, matter wave solitons Law  _et al._ (2010),
vortex-antivortex Wen  _et al._ (2013) etc. BECs of bosons in two different
hyperfine states, designated as spin-$1/2$ bosons, have opened up a new way
for synthetic SO coupling, a key ingredient for many important condensed
matter phenomena. In this connection, realization of SO coupled Bose and Fermi
gases at ultra-low temperature has paved the way for many important physical
phenomena of current interest such as measurement of spin Hall effect Galitski
and Spielman (2013); Aidelsburger  _et al._ (2013), topological insulators
Goldman  _et al._ (2010), topological superfluids Wu  _et al._ (2016),
atomtronics (or spintronics) Seaman  _et al._ (2007), and quantum computing
Andrianov and Moiseev (2014).

Numerical simulations have played an important role in unravelling different
interesting phases in coupled BECs. The stability and dynamics of matter-wave
bright and dark solitons in one-dimensional SO coupled BECs have been
investigated quite extensively using coupled Gross-Pitaevskii equations (GPEs)
Ravisankar2020sol ; Achilleos2013-bs ; Achilleos  _et al._ (2013). In two
dimension Rashba SO coupled BECs with weak harmonic trap exhibits plane and
stripe wave phases while under strong harmonic trap it displays the presence
of vortex pairs, honeycomb-lattice and half-quantum vortices upon varying the
Rabi (zero momentum) coupling strength Jin  _et al._ (2014). The SO coupled
BEC’s confined in the optical lattice reveal different ground state structures
as the form of vortex-antivortex pair Li and Sakaguchi (2013). Using the
numerical and variational analysis Cheng et al. demonstrated that SO coupled
BECs trapped under bichromatic optical lattices show Anderson Localization
Cheng  _et al._ (2014). It was found that the Rabi coupling stabilizes the
superfluid phase in coupled BECs in optical lattices He2021 . A variety of
collective modes, namely, Nambu-Goldstone, slosh, bifurcation modes were
studied in the trapped two-component quasi-two-dimensional (quasi-2D) BECs
Pal2017 ; Pal2018 .

The stability of different ground state phases of SO coupled BECs could be
well understood by analyzing the spectrum of elementary excitation. For
instance, the Bogoliubov-de Gennes (BdG) spectrum, much related to macroscopic
quantum phenomena, such as, superfluidity and superconductivity, provide the
fundamental information about the condensate dynamics. The excitation spectrum
of BECs with Rashba-Dresselhaus SO coupling found to exhibit roton-maxon
structures Martone  _et al._ (2012); Zheng  _et al._ (2013); Khamehchi  _et
al._ (2014); Ji  _et al._ (2015). Using Green’s function technique various
interesting features, like, multi-criticality, metastability, and the roton
were found to exist in three dimensional Rashba SO coupled BEC Liao  _et al._
(2015). The metastability could be understood as a result of the absence of
imaginary frequencies in the BdG spectrum. Ozawa et al. numerically
investigated the dynamical and energetic instabilities in quasi-1D SO and
Raman coupled BECs Ozawa2013 . Spin-dipole and breathing modes of the
collective excitation spectrum give clear picture of phase boundaries, which
was also confirmed by the quenching dynamics numerically Chen2017 . The effect
of SO, Rabi couplings and nonlinear interactions provide collective
oscillations, which has a transition from harmonicity to anharmonicity Yu2018
. In a recent work Geier et al. found the signature of Goldstone modes in
harmonically trapped SO coupled BECs Geier2021 .

The application of BdG spectrum was studied in the context of exploring the
superfluid phase in SO coupled BECs Zhu2012 . Further Yu et al. used this idea
in obtaining the ground state phase diagram particularly, excitations of zero
momentum phases in quasi-2D SO coupled BECs Yu (2013). In the liquid phase it
was demonstrated that the quasi-2D BECs exhibits different excitation modes,
like, phonon, roton and double roton modes with different SO couplings as the
interaction strengths are varied Sahu2020 .

It may be noted that most of the studies on the collective excitations are
mainly focused on quasi-1D spin-orbit coupled BECs Martone  _et al._ (2012);
Zheng  _et al._ (2013); Khamehchi  _et al._ (2014); Ji  _et al._ (2015) and
limited explorations are available on the stability of the superfluid phases
in two dimensions Wang  _et al._ (2010). For obvious reasons, it would be more
appropriate to study collective excitations in higher spatial dimensions. In
this paper we present a detailed study on the stability analysis of SO coupled
BECs in two dimensions from the excitation spectrum of Bogoliubov-de Gennes
(BdG) equations. In particular, we carry out a systematic analysis on
collective excitations of SO coupled BECs in quasi-two dimensions using the
dispersion relations obtained by the application of Bogoliubov theory to the
coupled Gross-Pitaevskii equations (GPEs). These dispersion relations are then
used for the stability analysis of plane waves, phonon-maxon-roton excitations
and the interplay between SO and Rabi couplings. The roton mode is the
precursor of the stripe phase with periodic fringes Ji  _et al._ (2015) and
instability is a fundamental ingredient of the existence of matter-wave
solitons.

The paper is organized as follows. In Sec. II, we introduce mean-field
theoretical model used for the study. In Sec. III, we analytically derive the
single particle dispersion relation. Following this, using Bogoliubov-de-
Gennes method, we investigate the collective excitation of Rashba SO coupled
BECs analytically and numerically by computing the eigenspectrum and
eigenvectors in Sec. IV. The numerical simulation of the stability analysis of
the ground state obtained by solving the coupled GPEs of quasi-two dimensional
Rashba SO coupled BEC is presented in Sec. V. Finally in Sec. VI we conclude
our observation of the SO coupled BEC.

##  II Mean-field model of coupled BECs

We consider a pseudospin-$1/2$ Bose-Einstein condensate with Rashba spin-orbit
and Rabi coupling which Hamiltonian is given by Jin  _et al._ (2014)

$\displaystyle\mathcal{H}=\mathcal{H}_{0}+\mathcal{H}_{I},$ (1)

with

$\displaystyle\mathcal{H}_{0}=$
$\displaystyle\int\Psi^{\dagger}\left[\frac{\boldsymbol{p}^{2}}{2m}+{V}+{k_{L}^{\prime}}\boldsymbol{p}\cdot\boldsymbol{\sigma}+{\Omega^{\prime}}\sigma_{x}\hbar\right]\Psi
d\mathbf{r},$ (2) $\displaystyle\mathcal{H}_{I}=$
$\displaystyle\int\left[\frac{{g}_{\uparrow\uparrow}}{2}|\Psi_{\uparrow}|^{2}+\frac{{g}_{\downarrow\downarrow}}{2}|\Psi_{\downarrow}|^{2}+{g}_{\uparrow\downarrow}|\Psi_{\uparrow}||\Psi_{\downarrow}|\right]d\mathbf{r}.$
(3)

Here $\Psi=\left(\Psi_{\uparrow}\;\Psi_{\downarrow}\right)^{T}$ is the two
component spinor normalized wave functions that satisfy the condition
$\int\left(|\Psi_{\uparrow}|^{2}+|\Psi_{\downarrow}|^{2}\right)dr=N$, with $N$
being the total number of particles, $m$ is the atomic mass, and
$\boldsymbol{p}=-\mathrm{i}\hbar\boldsymbol{\nabla}$ is the momentum operator,
$\boldsymbol{\sigma}=(\sigma_{x},\sigma_{y})$ are the $2\times 2$ Pauli
matrices, ${g}_{\uparrow\uparrow}$ and ${g}_{\downarrow\downarrow}=4\pi
a\hbar^{2}/m$, $(a=a_{\uparrow\uparrow}=a_{\downarrow\downarrow})$,
${g}_{\uparrow\downarrow}=4\pi a_{\uparrow\downarrow}\hbar^{2}/m$, are the
intra- and inter-species contact interaction strengths with $a$ and
$a_{\uparrow\downarrow}$ being the intra- and inter-species $s$-wave
scattering lengths, respectively. We consider the condensates confined in the
harmonic trap potential with form
$V=m[\lambda^{2}x^{2}+\kappa^{2}y^{2}+\eta^{2}z^{2}]/2$, where the trap aspect
ratios are $\lambda=\omega_{x}/\omega_{\perp}$,
$\kappa=\omega_{y}/\omega_{\perp}$ and $\eta=\omega_{z}/\omega_{\perp}\gg 1$.

In order to study the ground state and dynamical properties of Rashba SO
coupled BECs with strong axial traps, we consider the two-dimensional coupled
GPEs in dimensionless form as

$\displaystyle\mathrm{i}\frac{\partial\psi_{\uparrow}}{\partial t}=$
$\displaystyle\left[-\frac{1}{2}\nabla^{2}+V_{2D}(x,y)+\alpha|\psi_{\uparrow}|^{2}+\beta|\psi_{\downarrow}|^{2}\right]\psi_{\uparrow}$
$\displaystyle-\Lambda_{+}^{SO}\psi_{\downarrow},$ (4a)
$\displaystyle\mathrm{i}\frac{\partial\psi_{\downarrow}}{\partial t}=$
$\displaystyle\left[-\frac{1}{2}\nabla^{2}+V_{2D}(x,y)+\beta|\psi_{\uparrow}|^{2}+\alpha|\psi_{\downarrow}|^{2}\right]\psi_{\downarrow}$
$\displaystyle-\Lambda_{-}^{SO}\psi_{\uparrow},$ (4b)

where $\nabla^{2}=\partial_{x}^{2}+\partial_{y}^{2}$,
$V_{2D}(x,y)=\left(\lambda^{2}x^{2}+\kappa^{2}y^{2}\right)/2$ is the harmonic
trap potential and
$\Lambda_{\pm}^{SO}=\left(k_{L}\left(\mathrm{i}\partial_{x}\pm\partial_{y}\right)+|\Omega|\right)$.
In the above equations (4), length is measured in units of harmonic oscillator
length $a_{0}=\sqrt{\hbar/(m\omega_{\perp})}$, time in the units of
$\omega^{-1}_{\perp}$, and energy in the units of $\hbar\omega_{\perp}$. The
parameters $\alpha$ = $\sqrt{8\pi\eta}Na/a_{0}$ and
$\beta=\sqrt{8\pi\eta}Na_{\uparrow\downarrow}/a_{0}$ represent intra- and
inter-species contact interaction strengths, respectively. The Rashba SO
coupling and the Rabi coupling parameters have been rescaled as
$k_{L}=k_{L}^{\prime}/a_{0}\omega_{\perp}$ and
$\Omega={\Omega^{\prime}}/\omega_{\perp}$, respectively, while the wave
function is rescaled as
$\psi_{\uparrow,\downarrow}=\Psi_{\uparrow,\downarrow}a_{0}^{3/2}/\sqrt{N}$.
We consider the Rabi coupling as $\Omega=|\Omega|e^{\mathrm{i}\theta}$ that
minimizes the energy when $\Omega=-|\Omega|$ for $\theta=\pi$ Abad and Recati
(2013). The wave functions are subjected to the following normalization
condition,

$\displaystyle\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\left(|\psi_{\uparrow}|^{2}+|\psi_{\downarrow}|^{2}\right)\,dx\,dy=1,$
(5)

The stationary state solution of the wave function is given by

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

where, $j\in\\{\uparrow,\downarrow\\}$, $\psi_{jR}$ and $\psi_{jI}$ are the
real and imaginary part of the stationary wave function respectively,
$\mu_{\uparrow,\downarrow}$ are the chemical potential of spin-up and down
components respectively. Now using Eqs. (4a), (4b) and (6) the ground state
energy of the Rashba SO coupled BECs can be obtained as

$\displaystyle
E_{num}=\sum_{j=\uparrow,\downarrow}\frac{\iint\left(E_{j}^{2C}+E_{j}^{SO}\right)\,dx\,dy}{\iint\psi_{jR}^{2}\,dx\,dy}$
(7)

The detailed form of the $E_{j}^{2C}$ and $E_{j}^{SO}$ are given in the
Appendix A.

At first, we provide a detailed analysis on single-particle spectrum both in
$k_{x}$ and $k_{y}$ momentum directions. Next, we present the analytical and
numerical studies of the collective excitation spectrum to investigate the
stability of the ground states obtained from the coupled GP equations (4a) and
(4b) as the small fluctuation is added in the ground state. It is followed by
the detailed numerical analysis of the dynamics of the ground states.

##  III Single-particle spectrum

In this section, first we present the calculation of two-component coupled GP
equations with Rashba SO coupling for non-interacting BEC without any trapping
potential, which gives the “single particle spectrum”. Following this we
include the intra- and inter-species nonlinear contact interactions with small
perturbation, which gives rise to the “excitation spectrum” of the coupled
system. This also includes the interesting features that appear due to
interplay of Rashba SO ($k_{L}$) and Rabi couplings ($\Omega$).

Let’s consider the Eqs. (4a) and (4b) in the absence of trap and contact
interaction strengths, (i.e, $V=\alpha=\beta=0$) and use the plane wave
solution
$\psi_{\uparrow,\downarrow}=\phi_{\uparrow,\downarrow}\textrm{e}^{\mathrm{i}(k_{x}x+k_{y}y-\nu
t)}$ in the computational basis, we get

$\displaystyle\nu=\begin{pmatrix}\frac{1}{2}\left(k^{2}_{x}+k^{2}_{y}\right)&k_{L}\left(k_{x}-\mathrm{i}k_{y}\right)-|\Omega|\\\
k_{L}\left(k_{x}+\mathrm{i}k_{y}\right)-|\Omega|&\frac{1}{2}\left(k^{2}_{x}+k^{2}_{y}\right)\end{pmatrix},$
(8)

which gives the single particle energy spectrum $\nu$ from equation (8) as,

$\displaystyle\nu(\mathbf{k})_{\pm}=\frac{1}{2}\left(k^{2}_{x}+k^{2}_{y}\right)\pm\left(\sqrt{\left(k_{L}k_{x}-|\Omega|\right)^{2}+k_{L}^{2}k_{y}^{2}}\right).$
(9)

The single particle spectrum has two branches. First is the positive branch
($\nu_{+}$) which always have the single minimum and second is the negative
branch ($\nu_{-}$) that makes transition from the single minimum to the double
minima as the SO coupling strength ($k_{L}$) is increased for a fixed Rabi
coupling. In what follows we focus our study on analyzing the negative branch
as it exhibits transition between the single minimum to the double minima. We
analyze the spectrum with respect to $k_{x}$ and $k_{y}$ momentum directions
which will be useful to analyse the phase transition clearly.

Figure 1: Single particle energy spectrum in the $\\{k_{x},k_{y}\\}$ momentum
space for different set of coupling parameters indicated in the figure. Left
column indicates $\nu_{\pm}(k_{x},0)$ with $k_{y}=0$ and right column shows
$\nu_{\pm}(0,k_{y})$ for fixed $k_{x}=0$. Solid red line shows
$\nu_{-}(k_{x},k_{y})$ and green dashed line is for $\nu_{+}(k_{x},k_{y})$.

In $k_{x}$-direction: for zero Rabi and SO couplings ($\Omega=k_{L}=0$) the
spectrum is a non-degenerate parabolic single-particle dispersion spectrum
(Fig. 1(a)) which can also be seen from the Eq. 9. However, the spectrum
exhibits double minima characteristics for finite values of the SO coupling
strength ($k_{L}$). The minima are located at the position $k_{x}=\pm k_{L}$
(Fig. 1(c)). The value of energy minima increases with the SO coupling
strength as $\pm k_{L}^{2}/2$. For nonzero Rabi coupling ($\Omega\neq 0$) and
$k_{L}=0$ the spectrum exhibits a single minimum and develops an energy gap
between $\nu_{\pm}$ which is equal to $2\Omega$ (Fig. 1(e)). However, for non-
zero $k_{L}$ spectrum develops asymmetric double minima as shown in Fig.
1(g)). As the Rabi coupling is increased, the system attains the minimized
energy state. For example, with $k_{L}=0.7$, and $\Omega=0.5$, the spectrum
have a global minimum with $\nu_{-}=-0.745$ (See Fig. 1(g)). Further increase
in Rabi coupling from $\Omega=0.5$ to $0.7$ results in lowering the energy
state to $\nu_{-}=-0.945$ (in Fig. 1(i)). As $\Omega$ is increased further,
the spectrum develops a global minimum and the corresponding energy varies as
$k_{L}^{2}/2+\Omega$.

For the fixed Rabi coupling $\Omega=0.5$, as we change the SO coupling
strength from $k_{L}=0.7$ to $k_{L}=1.0$ we find that the energies change from
single minimum, i.e., $\nu_{-}=-0.745$ to double minima $\nu_{-}=\\{-1.0,0\\}$
respectively. It suggests that in general the atoms will get condensate in the
lowest energy state, i.e, global minimum (See Fig. 1(i)). While increasing
$k_{L}$ we found two minima with opposite wavevectors which will form stripe
like patterns (See Fig. 1(k)).

In $k_{y}$-direction: The spectrum exhibits two minima in the absence of Rabi
coupling with a finite $k_{L}$ (See Fig. 1(d)). The spectrum typically known
as Rashba ring in 2D momentum space Wang  _et al._ (2010); Jin  _et al._
(2014). As the Rabi coupling is increased for fixed $k_{L}$ a transition from
the double-minima state to the single minimum state is observed at critical
$\Omega(\approx k_{L}^{2})$ while symmetry of the system remains unchanged. We
find that the energy gap between these two energy spectrum is $2\Omega$ as
shown in the Figs. 1(h, j). The energy of the system gets lowered as Rabi
coupling is increased. As for an example, as the coupling parameter is changed
from $\Omega=0.5$ to $\Omega=0.7$ for a fixed $k_{L}=0.7$, the spectrum
energies change from $\nu_{-}=-0.5$ and $\nu_{-}=-0.7$ (See Figs. 1(h) and
1(j)). The energy strongly depends on Rabi coupling strengths for the cases
$\Omega>k_{L}^{2}$. However, for $\Omega<k_{L}^{2}$ the energy depends on SO
coupling strengths.

As we perform the energy comparison between $k_{x}$ and $k_{y}$ direction we
find that in presence of the Rabi coupling the energy have the lowest value in
the $k_{x}$ direction indicating the breaking of the rotational symmetry of
the single-particle energy spectrum. This feature suggests that plane wave
spectrum will get sustained with finite momentum for $\Omega>k_{L}^{2}$.
However, for $\Omega<k_{L}^{2}$ stripe wave will exist. We do not observe any
zero momentum (ZM) phase for $k_{L}\neq 0$ which was realized in one
dimensional SO coupled BEC Li2012 ; Ravisankar2020  and in two-dimensional
BECs Jin  _et al._ (2014); Bhuvaneswari2018  with nonlinear contact
interactions.

In the following section, we investigate the effect of couplings on the
excitation spectrum of the SO coupled BECs.

##  IV Analysis of the excitation spectrum of Rashba SO coupled BECs

In this section we present our analytical and numerical investigation of the
excitation spectrum of the coupled BEC with SO and Rabi couplings. In 1941
Landau initiated the concept of elementary excitation to explain the
Superfluid behaviour in 4He. The mathematical derivation of excitation
spectrum for the Bose gases was given by Bogoliubov in 1947 Bogolyubov (1947).
Collective excitation spectrum of BECs gives the basic information about the
dynamics of the quantum systems such as superfluid helium Zilsel (1950),
superconductors Rickayzen (1959), phonon-like excitation observed in 87Rb atom
in optical trap Jin  _et al._ (1996), and magnetic trap Mewes  _et al._
(1996), etc. Here our aim is to investigate the excitation spectrum
analytically which at the later part of the paper will be complemented with
the numerical simulation. At first, we analytically transform the pure ground
state wavefunction ($\psi_{\uparrow,\downarrow}$) by adding a small
perturbation ($\delta\psi_{\uparrow,\downarrow}$). After a direct algebraic
manipulations, the system has four different excitation branches, namely, two
positive and two negative branches. We analyze the characteristic of those
energy branches by varying the SO and Rabi coupling parameters and further
investigate the stability of the different modes.

###  IV.1 Analytical and numerical description of excitation spectrum

In order to understand the stability of our system, we calculate the
excitation spectrum of plane wave solutions using Bogoliubov theory. Let’s
assume that the total density of the system is $n=n_{\uparrow}+n_{\downarrow}$
and the chemical potential is $\mu$. Therefore the stationary state evolution
can be written as Goldstein and Meystre (1997); Abad and Recati (2013):

$\displaystyle\Psi_{j}$ $\displaystyle=\mathrm{e}^{-\mathrm{i}\mu
t}\left[\psi_{j}+\delta\psi_{j}\right],$ (10) $\displaystyle\delta\psi_{j}$
$\displaystyle=u_{j}\mathrm{e}^{\mathrm{i}(k_{x}x+k_{y}y-t\omega)}+v^{*}_{j}\mathrm{e}^{-\mathrm{i}(k_{x}x+k_{y}y-t\omega^{*})},$
(11)

where $\psi_{j}=\sqrt{n_{j}}\mathrm{e}^{\mathrm{i}\varphi_{j}}$,
$j=(\uparrow,\downarrow)$ is the ground state wave functions, $u_{j}$ and
$v_{j}$ are the amplitudes, $n_{j}$ and $\varphi_{j}$ are the density and
phase respectively. The Bogoliubov coefficients $u$’s and $v$’s could be
obtained by substituting Eq. (10) in Eqs. (4a) and (4b). Therefore, we have

$\displaystyle\omega\begin{pmatrix}u_{\uparrow}&v_{\uparrow}&u_{\downarrow}&v_{\downarrow}\end{pmatrix}^{T}=\mathcal{L}\begin{pmatrix}u_{\uparrow}&v_{\uparrow}&u_{\downarrow}&v_{\downarrow}\end{pmatrix}^{T}$
(12)

where the superscript $T$ denotes transpose of matrix and

$\displaystyle\mathcal{L}=\begin{pmatrix}f(n_{\uparrow},n_{\downarrow})&\alpha
n_{\uparrow}&L_{13}&\beta\sqrt{n_{\uparrow}n_{\downarrow}}\\\ -\alpha
n_{\uparrow}&-f(n_{\uparrow},n_{\downarrow})&-\beta\sqrt{n_{\uparrow}n_{\downarrow}}&-L_{24}\\\
L_{31}&\beta\sqrt{n_{\uparrow}n_{\downarrow}}&g(n_{\downarrow},n_{\uparrow})&\alpha
n_{\uparrow}\\\ -\beta\sqrt{n_{\uparrow}n_{\downarrow}}&-L_{42}&-\alpha
n_{\downarrow}&-g(n_{\downarrow},n_{\uparrow})\end{pmatrix},$ (13)

and other coefficients are given in the Appendix B. The normalization
condition yields

$\displaystyle\iint(|u_{j}|^{2}-|v^{*}_{j}|^{2})\,dx\,dy=1.$ (14)

The simplified form of Bogoliubov-de-Gennes (BdG) equation under the condition
$\text{det}\,\mathcal{L}=0$ will have the form for the interacting case with
$n_{\uparrow}=n_{\downarrow}=1/2$ as

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

By direct mathematical manipulation of the equation (15), we obtained four
dispersion relations. The complicated expressions of the coefficients
($b,c,d$) are given in the Appendix C.

Following this we corroborate the analytical results for the excitation
spectrum by numerically solving the BdG equations from which we also obtain
the eigenvectors as a function of $k_{x}$ and $k_{y}$. First we consider a
$[-1000:1000][-1000:1000]$ grid in real space with step size
$h_{x}=h_{y}=0.05$. Then we use the Fourier collocation method where we
numerically perform the Fourier transformation of BdG equations and obtained a
truncated reduced BdG matrix, which was subsequently diagonalized using the
LAPACK package Anderson1999 . In momentum space we consider $[-50:50][-50:50]$
modes in $k_{x}$, $k_{y}$ directions with a grid step size of
$h_{k_{x}}=h_{k_{y}}=0.0628$.

###  IV.2 Comparison of excitation spectrum with and without interactions

In the following we analyze the effect of the couplings on the stability of
the negative branch of the energy spectrum. Note that the negative eigenenergy
of the excitation spectrum implies that the system is energetically unstable,
while, the imaginary or complex eigenenergies indicate the dynamical
instability Ozawa2013 . Bogoliubov-de-Gennes excitation spectrum obtained from
equation (15), without SO coupling ($k_{L}=0$) is similar as obtained in Ref.
Abad and Recati (2013); Goldstein and Meystre (1997). Apart from this the
dispersion relation has two dimensional Rashba SO and Rabi couplings. Fig. 2
shows system of collective BdG excitation spectrum. The comparison between the
noninteracting (solid red lines) and interacting ($\alpha=\beta=1$) cases
(shown with the green and blue lines) indicate that the former case does not
have any imaginary part while later case have the imaginary part. This
particular features suggests that the metastability of the plane wave phase of
SO coupled BECs is destroyed in the $k_{y}$-direction which makes the system
dynamically unstable indicating the system lacks any superfluid behaviour
Zhu2012 .

Figure 2: Collective excited dispersion of $\Re(\omega_{-})$ (solid red line)
for non-interacting and of $\Re(\omega_{-})$ (solid green line) and
$\Im(\omega_{-})$ (solid blue line) for interacting ($\alpha=\beta=1$) case.
The parameters are $\Omega=0.5$, $k_{L}=1$. (a) Variation of negative spectrum
along the $k_{x}$ direction for $k_{y}=0$ and (b) Variation of negative
spectrum along the $k_{y}$ direction for $k_{x}=0$. Solid lines represent the
analytical results obtained from BdG Eq. (15) and open circles denote the
numerical solution of Eq. (12). Figure 3: The eigenvectors obtained by solving
Eq. (12): Red dots for $|u_{\uparrow}|$ component, Black open diamond for
$|u_{\downarrow}|$, green open square for $|v_{\uparrow}|$ and magenta
triangle for $|v_{\downarrow}|$. Top row is for the noninteracting case and
bottom row is for the interacting case. Left column shows the variation of
eigenvectors along the $k_{x}$ directions and right column those variation
along the $k_{y}$ directions. All the other parameters are same as in Fig. 2.
For non-interacting case plane wave, density mode dominates for all range of
the wave number, while for the interacting case a transition from the density
to spin mode occurs where the negative spectrum becomes complex in $k_{y}$
momentum direction.

In Fig. 3, we plot the eigenvectors corresponding to the eigenenergy spectrum
as shown in the Fig. 2. For non-interacting case the in-phase of the
eigenvector components for all wave number indicate the presence of only
density like modes (in-phase) in both momentum directions (See Fig. 3(a,b)).
However, for interacting case ($\alpha=\beta=1$) presence of some complex
patterns are observed. There is transition from the density like modes (in-
phase) to the spin-like mode (out-phase) happens in $k_{y}$ direction as shown
in the Fig. 3(c,d). We find that both eigen components of the eigenvectors
$u$’s and $v$’s are in-phase in the $k_{x}$ direction. For $k_{x}\approx 0$
eigenvectors approach towards each other and have equal values at $k_{x}=0$,
which indicates the presence of phonon mode. Beyond this there is sudden
increase in the value of the eigenvectors. At finite $k_{x}$ again there is a
change in the curve that corresponds to the maxon point which is followed by
returning to the density-like mode. In $k_{y}$ direction, we found two types
of behaviour. First up to the $\Re(\omega_{-})$ the $u$’s and $v$’s are in-
phase, while the presence of $\Im(\omega_{-})$ changes $u$’s and $v$’s in out-
phase which indicate the presence of spin-like mode Abad and Recati (2013).

###  IV.3 Effect of Rabi coupling on the excitation spectrum

Theoretically it was noticed that excitation spectrum of quasi-one-dimensional
Raman SO coupled BECs consists of roton like minimum for the case of finite
Rabi coupling Martone  _et al._ (2012); Zheng  _et al._ (2013) which was also
confirmed experimentally in Khamehchi  _et al._ (2014); Ji  _et al._ (2015).
In this section we present the effect of Rabi coupling on the excitation
spectrum in two dimension coupled BECs.

Figure 4: Excitation spectrum in $k_{x}$ (first row) and $k_{y}$ (second row)
momentum directions for different $\Omega$. $\Omega$ varies in column (a)-(d)
= ($0,0.3,0.5,0.6$) and other fixed parameters are $k_{L}=0.75$,
$\alpha=\beta=1$. First row: Solid red line and solid green line indicate
$\omega_{-}(k_{x},0)$ and $\omega_{+}(k_{x},0)$ respectively; Second row: Red
solid line and solid blue line depicts that $\Re(\omega_{-}(0,k_{y}))$ and
$|\Im(\omega_{-}(0,k_{y}))|$ respectively. Solid lines are analytical results
from BdG Eq. (15) and open circles are obtained numerically solving Eq. (12).
First row shows roton softening upon increase in $\Omega$, which reflects as a
disappearance of the imaginary modes for high $\Omega$.

As Rabi coupling is decreased we find that beyond the critical value of
$\Omega$ roton-like minimum starts softening, it develops a negative
frequency, which indicates the appearance of instability in the system. We fix
the SO coupling parameter to $k_{L}=0.75$ with $\alpha=\beta=1$, $k_{y}=0$,
and vary the Rabi strength $\Omega$. We consider $\Omega=0,0.3,0.5,$ and
$0.6$. In Fig. 4 we show the dispersion behaviour for different $\Omega$ by
keeping the SO coupling parameter fixed. We notice the presence of minimum
attributed to the phonon, maxon and roton. At zero Rabi coupling ($\Omega=0$)
strength the excitation spectrum possesses negative frequency in $k_{x}$
direction indicating that the system is energetically unstable. However, in
$k_{y}$ direction the system exhibits dynamical instability due to the
presence of the complex excitation frequency. Owing to this feature at zero
Rabi coupling system does not exhibit any superfluid behaviour Zhu2012 . As
the Rabi coupling strength is increased to $\Omega=0.3$ the eigenfrequency
suggests the presence of minimum related to phonon-maxon-roton in $k_{x}$
direction while in the transverse direction ($k_{y}$) it suggests the presence
of the complex eigen-frequencies (See Fig. 4(b1) and Fig. 4(b2)). This
behaviour indicates the dynamically unstable state for this SO coupling
strength. For $\Omega=0.5$, the system remains dynamically unstable, however,
decrease in the roton minimum and amplitude of the imaginary frequency is
observed (See Fig. 4(c1) and Fig. 4(c2)).

As the Rabi coupling strength is increased beyond a threshold value
$\Omega\geq 0.56$ the system does not show any complex frequency in $k_{y}$
direction indicating the stable behaviour. Only the presence of real frequency
in the $k_{y}$ direction indicates metastable state, also the axial symmetry
is preserved (See Fig. 4(d1) and Fig. 4(d2)). Overall we find that the
increase in the Rabi coupling strength for a fixed SO coupling leads the
stabilization of the system.

Figure 5: The eigenvectors corresponding to the eigenspectrum of Fig.4. The
representation of the eigenvector components is same as in Fig 3. Top panel
represents the eigenvectors corresponding to $\omega_{+}(k_{x},0)$ (solid
green line), middle panel for $\omega_{-}(k_{x},0)$ (solid red line) and
bottom panel for $\omega_{-}(0,k_{y})$. Maxon mode in the top panel, phonon-
maxon mode in the middle panel and spin to density mode transition in the
bottom panel are noticeable.

From the excitation spectrum we analyzed the phonon-maxon-roton modes. Now in
order to get more insight about the stability of these modes we move to the
analysis of the characteristics of the eigenvectors in the momentum space.
Figs. 5 shows the eigenvectors corresponding the eigen spectrum given in the
Fig. 4. The first row displays the eigenvectors corresponding to
$\omega_{+}(k_{x},0)$ while second and third row demonstrate the the
eigenvectors of $\omega_{-}(k_{x},0)$, $\omega_{-}(0,k_{y})$, respectively.
For $\Omega=0$ with $k_{L}=0.75$ both eigenvectors in $k_{x}$ momentum
direction are in-phase except near $k_{x}=0$. At $k_{x}=0$ both $u$’s and
$v$’s becomes equal indicating the presence of phonon mode, and absence of
maxon-mode. However, for $\omega_{-}$ we observe a complicated spin-flipping
behaviour in both $k_{x}$ and $k_{y}$ directions of eigenvectors, which is
with respect to the excitation spectrum’s negative and imaginary eigen
frequency.

Similar to the earlier studies Goldstein and Meystre (1997); Abad and Recati
(2013); Tommasini2003 ; Recati2019  we observe the two branches in the
spectrum. First one represents to the gapless density mode (in-phase mode)
corresponding to the Goldstone mode with U(1) symmetry and second one denotes
the gapped spin mode (out-of-phase mode). As we carefully analyze their
corresponding eigenmodes we notice some important differences. At one hand for
the density mode we obtain
$|u_{1}|^{2}-|u_{2}|^{2}=|v_{1}|^{2}-|v_{2}|^{2}=0$, which in terms of spin
language is an unpolarized mode, while at the other hand for the spin mode we
find breaking of the $\mathbb{Z}_{2}$ symmetry, meaning that
$|u_{1}|^{2}-|u_{2}|^{2}\neq 0$ and $|v_{1}|^{2}-|v_{2}|^{2}\neq 0$ that
implies a polarized mode in the spin language Recati2019 .

In Fig. 4, we display the eigenmodes corresponding to the excitation spectra
for different values of $\Omega$. For $\Omega<0.6$ (Figs. 4(a)-4(c)) the
excitation spectra exhibit complex eigenvalues for small values of $k_{y}$
which indicates that the system is dynamically unstable for such parameter
values. As we look at the eigenmodes of such excitation spectra we find that
those instabilities can be attributed to the generation of the spin modes.
Increase in value of $\Omega$ to $\Omega=0.6$ results complete suppression of
such unstable spin mode and emergence of the density mode in the system. The
same scenario appear to be present for the parameters of Fig. 7, where
increase of $k_{L}$ generates dynamical instabilities in the system and that
manifest as the emergence of the spin modes.

Figure 6: Variation of different quantities with $\Omega$ for $\alpha=1$,
$\beta=1$ and $k_{L}=0.75$. (a) Variation of Maxon (red open inverted
triangles) and Roton minima (blue open diamond), (b) Decay in roton depth
($\Delta_{RD}$), (c) Position of maxon (red open circles) and roton minima
(black dotted lines) in $k_{x}$ momentum direction, (d) Amplitude loss in
instability (cyan open circles) and corresponding bandwidth (magenta circles)
in the $k_{y}$ momentum space.

We also analytically compute the effect of $\Omega$ with fixed parameters
$k_{L}=0.75,\alpha=\beta=1$ on the maxon and roton minimum from the collective
excitation spectrum. In Fig. 6(a) solid red line with triangle, we show the
maxon and blue line with diamond we indicate the roton minimum. The points of
maxon and roton are well fitted with the straight line beyond the critical
value of $\Omega$ both maxon-roton points are merged as shown in the Fig.
6(a). As we understand that the difference between roton minimum and maxon
gives the roton gap/depth ($\Delta_{RD}$) as displayed in Fig. 6(b). As the
Rabi coupling is increased, the roton minima disappears. Beyond the critical
value only maxon in $k_{x}$ direction is observed. We find that the outcome of
increase in Rabi coupling strength is the decay of roton gap/depth.

In Fig. 6(c) we show the positions of maxon and roton minimum positions along
$k_{x}$ momentum space. We find that the roton minima are always present at
$k_{x}=0.75$, which appears to be similar as we we observed for fixed Rashba
SO coupling strength $k_{L}=0.75$. However, upon looking the behaviour
carefully we find that the position of maxon increases upon increase of the
Rabi coupling strength and finally at $\Omega\approx 0.8$ the maxon and the
roton gets merged. As a consequence of this above this Rabi coupling strength
only maxon remains present in the system. In the similar way, we notice the
loss in the amplitude of instability upon increasing of $\Omega$ denoted as
$A_{L}$ (amplitude loss). In Fig. 6(d) we plot the variation of $A_{L}$ and
instability bandwidth ($B_{W}$) with the Rabi coupling which clearly
complement our results attributed to stabilization of the system beyond
$\Omega\approx 0.8$.

Figure 7: Variation of excitation spectrum in $k_{x}$ (first row) and $k_{y}$
(second row) momentum directions for different SO coupling strengths $k_{L}$,
varies in column (a-e) = $(0,0.5,1.0,1.5,1.75)$. Other parameters are
$\Omega=0.75,\alpha=\beta=1$. First row: Green and solid red lines indicate
$\omega_{+}(k_{x},0)$ and $\omega_{-}(k_{x},0)$ respectively; Second row:
solid red and blue lines show $\Re(\omega_{-}(0,k_{y}))$ and
$|\Im(\omega_{-}(0,k_{y}))|$. Solid lines are analytical results from Eq. (15)
and open circles are numerically obtained from Eq. (12). The increase in
$k_{L}$ leads to generation of imaginary energy frequency modes in the $k_{y}$
direction.

###  IV.4 Role of Rashba spin-orbit coupling ($k_{L}$)

In this section, we investigate the effect of SO coupling strength on the
excitation spectrum by fixing the other parameters, like, $\Omega=0.75$ and
$\alpha=\beta=1$. In absence of SO coupling ($k_{L}=0$) we find that the
spectrum only have the real frequencies which are symmetric in both $k_{x}$
and $k_{y}$ directions, also exploring the zero momentum phase (see Fig. 7(a1)
and Fig. 7(a2)). As the SO coupling is raised to $k_{L}=0.1$ the maxon mode
appears in the system. On further increase in $k_{L}$ leads to the generation
of roton minimum in the spectrum as shown in Fig. 7(b1). For $k_{L}>0.86$, we
find transition from stable to the unstable state which is quite evident by
the presence of complex frequencies in $k_{y}$ directions of
$\omega_{-}(0,k_{y})$ spectrum as shown in the Fig. 7(c2). Further increase in
$k_{L}$ leads to increase in the magnitude of the instability (see Fig.
7(d2)). For $k_{L}>1.5$, we notice the negative frequency with two minima in
$k_{x}$ direction, which gives two stability bands in both side of $k_{y}$
direction (See Fig. 7(e1) and Fig. 7(e2)). We find that the role of SO
coupling here is to gradually increase the phonon-maxon and roton minimum up
to the critical value of $k_{L}$ in $k_{x}$ momentum space. Increasing $k_{L}$
beyond a critical value leads to the loss of phonon-maxon and roton minimum
apart of having the negative energy that shows the system becomes
energetically unstable and thus lacks any superfludity Zhu2012 .

Figure 8: Eigenvectors associated with Fig. 7, which are obtained numerically
from Eq. (15). The representation of the eigenvector components is same as in
Fig 3. First row: $\omega_{+}(k_{x},0)$; Second row: $\omega_{-}(k_{x},0)$;
Third row: $\omega_{-}(0,k_{y})$. First row shows the evolution of maxon
modes, second row displays phonon-maxon to complicated spin-flipping modes,
third row clearly shows the transition from density to spin-mode.

In Fig. 8 we plot the eigenvectors corresponding to eigen spectrum as shown in
the Fig. 7. For $k_{L}\neq 0$ cases, in $\omega_{+}(k_{x},0)$ spectrum’s $u$’s
($v$’s) exhibits decreasing (increasing) trend to reach the point of maxon-
mode which also shows the coupling point between local minimum of
$\omega_{+}(k_{x},0)$ and maxon of $\omega_{-}(k_{x},0)$. In the $k_{y}$
momentum direction we initially notice the presence of density modes that get
transmuted to spin-density modes as $k_{L}$ is increased for fixed $\Omega$
(Fig. 8 third row). Through the eigenvector analysis we establish the
existence of phonon-maxon modes.

Figure 9: Stability phase diagram illustrating the stable and unstable phases
in $k_{L}-\Omega$ plane with fixed weak repulsive nonlinear contact
interaction strengths, $\alpha=\beta=1$. Green dots represent the phase
transition points of unstable and stable boundary obtained from eigen spectrum
using Eq. (15).

###  IV.5 Stability phase diagram

In order to provide a detailed picture of the stability of ground state phases
we identify a phase transition from stable to unstable phase in $k_{L}-\Omega$
plane using excitation spectrum and illustrate this in the phase-diagram (see
Fig. 9). By simultaneously varying $\Omega$ and $k_{L}$ real and complex
frequencies where present in the $k_{y}$ direction those critical points are
noted for a fixed $k_{x}=0$. We show the critical points in the $k_{L}-\Omega$
plane with the green dots. We characterize the phase as stable one if the
frequencies are real otherwise phase is denoted as unstable for the complex
frequencies Goldstein and Meystre (1997); Ozawa2013 ; Zhu2012 .

After having the fair understanding about the different phases of excitation
spectrum of the coupled BEC using BdG equations now in the following section
we complement our observation by directly solving the set of SO coupled GPE.

##  V Numerical simulation results

In this section, we present the numerical results which is obtained by solving
the dynamical equations (Eqs. 4) of the coupled BEC. We use the imaginary-time
propagation (ITP) method to obtain the ground state of the system.
Subsequently, we evolve the ground state wave-function by using the
conventional real-time propagation (RTP), where in both the methods, we
adopted the split-step Crank-Nicholson scheme Muruganandam and Adhikari
(2009); Young2016 ; Ravisankar2021 ; Muruganandam2021 . We have considered the
grid sizes $400\times 400$ with space steps $dx=dy=0.1$, and time step
$dt=0.005$ is used. Initially, we obtained the ground state wave-functions
using ITP method with respect to stable ($k_{L}=0,\Omega=0.75$) and unstable
($k_{L}=1.5,\Omega=0.75$) regime of the stability phase diagram of Fig. 9.
Once the ground state wave-function is obtained, we then evolve it employing
the RTP method, in which, we quench the system by reducing the trap strength
as $\lambda=\kappa=1.0\to 0.5$ at $t=0$.

Figure 10: Time evolution of the ground state density for $\Omega=0.75$,
$k_{L}=0$ and $\alpha=\beta=1$. First column (a) is the actual ground state
wavefunctions obtained from imaginary-time propagation method at $t=0$. (b-e)
are their dynamical density patterns at different time instant
$t=50,100,200,300$ (column wise). Top row indicates spin-up
($|\psi_{\uparrow}|^{2}$) density, spin-down ($|\psi_{\downarrow}|^{2}$) in
the middle, and total density ($|\psi_{\uparrow}|^{2}$ \+
$|\psi_{\downarrow}|^{2}$) in the bottom row. Figure 11: One-dimensional
density patterns corresponding to Fig. 10. First row represents the
$x$-direction densities where $y=0$, and second row for $y$-direction
densities for $x=0$, $|\psi_{\uparrow}|^{2}$ (red line),
$|\psi_{\downarrow}|^{2}$ (blue line), $|\psi_{\uparrow}|^{2}$ \+
$|\psi_{\downarrow}|^{2}$ (black line).

To begin with, we consider the parameters $\Omega=0.75$, $k_{L}=0$ and
$\alpha=\beta=1$ with $\lambda=\kappa=1$ that lie in the stable region. For
these parameters we find that the ground state wave-function is of plane wave
(PW) nature (See Fig. 10(a)). In order to analyze the stability of this state
we evolved the ground state wave function using RTP. In Figs. 10(b-h) we show
the temporal evolution of the ground state upto $t=300$. As we look at the
evolution carefully we find that as time progresses condensate experiences the
expansion due to the presence of repulsive interactions. So to stabilize
against the expansion trap is added Ravisankar2020sol . At $t=0$ owing to
quenching the trap strength by half the condensate expands for a time interval
$0<t<50$ (See Fig. 7(b)). For $t>50$ the size of the condensate does not
change and system attains the stable steady state (Figs. 10(c-e)). This is
evident from the evolution of the density of each spin component as shown in
the top and middle rows of Figs. 10. Also the similar feature can be seen for
the total density (bottom row of Figs. 10). Moreover we find that the density
profiles of both the components are symmetric in nature and do not exhibit any
oscillations (Fig. 11). Thus they lack any polarization which is also quite
evident from the behaviour of the eigenvectors that clearly show the presence
of the Goldstone mode (in-phase-mode) Abad and Recati (2013); Recati2019  (See
first column of Fig 8).

Figure 12: (a) Plot depicts time evolution of root-mean-squared (rms) size in
$\langle x^{2}\rangle$ (solid red line), $\langle y^{2}\rangle$ (solid green
line), and $\langle r^{2}\rangle$ (solid blue line) directions, (b) illustrate
the maximum of density $\textrm{max}(|\psi_{\uparrow,\downarrow}|^{2})$ (solid
red and green lines),
$\textrm{max}(|\psi_{\uparrow}|^{2}+|\psi_{\downarrow}|^{2})$ (solid blue
line) and (c) total system energy, the parameters are same as in Fig. 10.

In order to verify the dynamics of the condensate next we calculate the root-
mean-square (rms) size of the condensate. Initially, we find that the
condensate diffuses after the trap is quenched. The rms size of $\langle
x^{2}\rangle$, $\langle y^{2}\rangle$, and $\langle r^{2}\rangle$ are
illustrated in Fig. 12(a). The peak densities of both components
$\textrm{max}(|\psi_{\uparrow,\downarrow}|^{2})$ and
$\textrm{max}(|\psi_{\uparrow}|^{2}+|\psi_{\downarrow}|^{2})$ decrease for a
while, thereafter, they attain a steady state which is shown in the Fig.
12(b). However, the condensate size in both $x$ and $y$ directions is the same
as well the maximum density of both components have the similar feature. In
Fig. 12(c) we plot the temporal evolution of the energy which confirms that
the system attains the steady state after $t>50$ consistent with the evolution
of the density profile. As we analyze the condensate evolution more carefully,
as shown in the Figs. 10 and 12, we find that the ground state is dynamically
and energetically stable. It also coincides with the stability regime in the
phase diagram of Fig. 9, which we get from the BdG excitation spectrum.

Figure 13: Dynamics of ground state density in $x-y$ plane, same as in the
Fig. 10 for parameters $k_{L}=1.5$, $\Omega=0.75$ with $\alpha=\beta=1$.
Figure 14: One-dimensional density patterns as in Fig. 11. The parameters
correspond to Fig. 13.

In order to get more insight about the behaviour of the condensate we also
analyze the dynamics of instability regime which is investigated by fixing the
parameters as $\Omega=0.75$, $k_{L}=1.5$ and $\alpha=\beta=1$ with trap
strength $\lambda=\kappa=1$. We obtain the ground state as elongated plane
wave phase as illustrated in the Fig. 13(a). However, the BdG excitation
spectrum have negative energy minimum as well as complex eigen frequencies
(see Fig. 7(d)) in $k_{x}$ and $k_{y}$ directions respectively. This feature
indicates that the system is energetically (See Fig. 15(c)) and dynamically
(See Fig. 13(b-e)) unstable. Apart from this, we also capture the temporal
evolution of the the density corresponding the each spin component along with
the total density as shown in Fig. 14. Note that these wave functions are also
related to the spin like mode Abad and Recati (2013); Recati2019  as quite
evident from the nature of eigenvectors in Fig. 8(see fourth column). In Fig.
15 we show the dynamical evolution of different physical quantities like,
root-mean-square size of the condensate, maximum of the density and total
energy of the system. As we look at the evolution of the condensate rms size
in the $x$\- and $y$-direction as illustrated in Fig. 15(a), we find the
condensate size decreases in the $x$-direction while it increases in the
$y$-direction with time. However, the spin-components are polarized that leads
to the emergence of interference pattern as time progresses. This suggests the
presence of instability in the system which manifestation can also be seen in
the temporal evolution of the total density profile as shown in the bottom row
of Fig. 13. Here the density profile that initially was quite symmetric in
both directions becomes elongated along the $y$-direction at $t=300$. The
other way through which we can characterize the unstable nature of the
condensate is by looking at the temporal evolution of the maximum of the the
density corresponding to the spin component. Fig. 15(b) depicts the evolution
of the maximum of the spin component as well as maximum corresponding to the
total density. After $t\gtrsim 150$ we notice a gradual decrease in their
values signifying the unstable behaviour. Finally we show the temporal
evolution of the total condensate energy in Fig. 15(c). We find a sharp
increase in the total energy beyond $t\gtrsim 20$, that clearly indicates
towards dynamically unstable nature of the condensate.

Figure 15: Time evolution of the root-mean-square (rms) size in $\langle
x^{2}\rangle$ (solid red line), $\langle y^{2}\rangle$ (solid green line) and
$\langle r^{2}\rangle$ (solid blue line) directions, (b) corresponds to the
maximum of density $\text{max}(|\psi_{\uparrow,\downarrow}|^{2}$ (solid red
and green lines),
$\textrm{max}(|\psi_{\uparrow}|^{2}+|\psi_{\downarrow}|^{2})$ (solid blue
line) and (c) total system energy, the parameters are same as in Fig. 13.

##  VI Summary and Conclusions

In this paper we have investigated the collective excitation spectrum of
Rashba SO coupled Bose-Einstein condensate, with Rabi mixing in two-dimensions
using the Bogoliubov-de-Gennes theory. First we have analyzed the dispersion
of single-particle system in the momentum space, in which as consequences of
Rabi coupling the system has broken rotational symmetry, as $k_{x}$ momentum
space has lower ground state energy compared with $k_{y}$ momentum direction.
Following this, we have presented the excitation spectrum with and without
nonlinear contact interactions. In the case of noninteracting system,
excitations have real negative frequencies, implying that the system is
energetically unstable. For interacting case, increases in the Rabi strength
for fixed $k_{L}$ leads the system to make a transition from the unstable
phase to the stable phase. We find that the phonon-maxon-roton modes disappear
upon increase of the Rabi strengths that leads to the loss of the roton
gap/depth. The position of roton minimum only present at $k_{x}=k_{L}$, and
the maxon position approximately equal to the Rabi strengths. The increase of
Rabi coupling strength resulted in the loss of instability and bandwidth.
However, in contrast to the Rabi coupling, the SO coupling has a destabilizing
effect. Further, the Rashba SO coupling turns a symmetric system to an
asymmetric one. By increasing the SO strengths, the phonon-maxon-roton modes
and instability are revealed in $k_{x}$ and $k_{y}$ momentum directions. After
a critical value of SO strength, we have negative frequency and additional
increase in $k_{L}$ shows double minimum in the $k_{x}$ and two instability
bands in the $k_{y}$ momentum direction.

We confirm the observation of the BdG spectrum by directly solving the coupled
GPEs. We find the presence of plane wave and elongated plane wave. Their size
and maximum of density and energies are studied in the time evolution, which
clearly explain about the dynamic and energetic stability and instability
phases, also confirms the phenomena from excitation spectrum.

###### Acknowledgements.

R.R. acknowledges DST-SERB (Department of Science & Technology - Science and
Engineering Research Board) for the financial support through Project No.
ECR/2017/002639 and UGC (University Grants Commission) for financial support
in the form of UGC-BSR-RFSMS Research Fellowship scheme (2015-2020). A.G.
acknowledges FAPESP 2016/17612-7 and CNPq (Conselho Nacional de
Desenvolvimento Científico e Tecnológico) No. 306920/2018-2. The work of P.M.
is supported by CSIR (Council of Scientific and Industrial Research) under
Grant No. 03(1422)/18/EMR-II, DST-SERB under Grant No. CRG/2019/004059, FIST
(Department of Physics), DST-PURSE and MHRD RUSA 2.0 (Physical Sciences)
Programmes, and FAPESP (Fundação de Amparo à pesquisa do Estado de São Paulo)
Grant No. 2016/00269-8. P.K.M. acknowledges DST-SERB for the financial support
through Project No. ECR/2017/002639.

##  Appendix A A calculation of energy of SO coupled BECs

In this Appendix we provide the detailed steps to obtain the total energy of
the SO coupled BECs. We substitute the stationary state form of the wave
function given in Eq. (6) in the Eq. (4) and further separate the real and
imaginary parts, which yields

$\displaystyle\mu_{\uparrow}\psi_{\uparrow R}=$
$\displaystyle\left[-\frac{1}{2}\frac{\partial^{2}}{\partial
x^{2}}-\frac{1}{2}\frac{\partial^{2}}{\partial
y^{2}}+V_{2D}(x,y)+\alpha\lvert\psi_{\uparrow}\rvert^{2}+\beta\lvert\psi_{\downarrow}\rvert^{2}\right]\psi_{\uparrow
R}+k_{L}\left(\frac{\partial\psi_{\downarrow I}}{\partial
x}-\frac{\partial\psi_{\downarrow R}}{\partial
y}\right)+\Omega\psi_{\downarrow R},$ (16a)
$\displaystyle\mu_{\downarrow}\psi_{\downarrow R}=$
$\displaystyle\left[-\frac{1}{2}\frac{\partial^{2}}{\partial
x^{2}}-\frac{1}{2}\frac{\partial^{2}}{\partial
y^{2}}+V_{2D}(x,y)+\beta\lvert\psi_{\uparrow}\rvert^{2}+\alpha\lvert\psi_{\downarrow}\rvert^{2}\right]\psi_{\downarrow
R}+k_{L}\left(\frac{\partial\psi_{\uparrow I}}{\partial
x}+\frac{\partial\psi_{\uparrow R}}{\partial y}\right)+\Omega\psi_{\uparrow
R},$ (16b) and for the imaginary parts we have
$\displaystyle\mu_{\uparrow}\psi_{\uparrow I}=$
$\displaystyle\left[-\frac{1}{2}\frac{\partial^{2}}{\partial
x^{2}}-\frac{1}{2}\frac{\partial^{2}}{\partial
y^{2}}+V_{2D}(x,y)+\alpha\lvert\psi_{\uparrow}\rvert^{2}+\beta\lvert\psi_{\downarrow}\rvert^{2}\right]\psi_{\uparrow
I}-k_{L}\left(\frac{\partial\psi_{\downarrow R}}{\partial
x}+\frac{\partial\psi_{\downarrow I}}{\partial
y}\right)+\Omega\psi_{\downarrow I},$ (16c)
$\displaystyle\mu_{\downarrow}\psi_{\downarrow I}=$
$\displaystyle\left[-\frac{1}{2}\frac{\partial^{2}}{\partial
x^{2}}-\frac{1}{2}\frac{\partial^{2}}{\partial
y^{2}}+V_{2D}(x,y)+\beta\lvert\psi_{\uparrow}\rvert^{2}+\alpha\lvert\psi_{\downarrow}\rvert^{2}\right]\psi_{\downarrow
I}-k_{L}\left(\frac{\partial\psi_{\uparrow R}}{\partial
x}-\frac{\partial\psi_{\uparrow I}}{\partial y}\right)+\Omega\psi_{\uparrow
I},$ (16d)

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. (16a) with $\psi_{\uparrow
R}$ and Eq. (16b) with $\psi_{\downarrow R}$, and integrating we get

$\displaystyle\mu_{\uparrow}\int dx\,dy\psi_{\uparrow R}^{2}=$
$\displaystyle\int dx\,dy\psi_{\uparrow
R}\left\\{\left[-\frac{1}{2}\frac{\partial^{2}}{\partial
x^{2}}-\frac{1}{2}\frac{\partial^{2}}{\partial
y^{2}}+V_{2D}(x,y)+\alpha\lvert\psi_{\uparrow}\rvert^{2}+\beta\lvert\psi_{\downarrow}\rvert^{2}\right]\psi_{\uparrow
R}+k_{L}\left(\frac{\partial\psi_{\downarrow I}}{\partial
x}-\frac{\partial\psi_{\downarrow R}}{\partial
y}\right)+\Omega\psi_{\downarrow R}\right\\},$ (17a)
$\displaystyle\mu_{\downarrow}\int dx\,dy\psi_{\downarrow R}^{2}=$
$\displaystyle\int dx\,dy\psi_{\downarrow
R}\left\\{\left[-\frac{1}{2}\frac{\partial^{2}}{\partial
x^{2}}-\frac{1}{2}\frac{\partial^{2}}{\partial
y^{2}}+V_{2D}(x,y)+\beta\lvert\psi_{\downarrow}\rvert^{2}+\alpha\lvert\psi_{\downarrow}\rvert^{2}\right]\psi_{\downarrow
R}-k_{L}\left(\frac{\partial\psi_{\uparrow I}}{\partial
x}+\frac{\partial\psi_{\uparrow R}}{\partial y}\right)+\Omega\psi_{\uparrow
R}\right\\},$ (17b)

Rearranging the above equations yield

$\displaystyle\mu_{\uparrow}=$ $\displaystyle\frac{1}{\int\psi_{\uparrow
R}^{2}dx\,dy}\int\left[\frac{1}{2}\left(\frac{\partial\psi_{\uparrow
R}}{\partial x}\right)^{2}+\frac{1}{2}\left(\frac{\partial\psi_{\uparrow
R}}{\partial
y}\right)^{2}+\left(V_{2D}(x,y)+\alpha\lvert\psi_{\uparrow}\rvert^{2}+\beta\lvert\psi_{\downarrow}\rvert^{2}\right)\psi_{\uparrow
R}^{2}\right]dx\,dy$ $\displaystyle+\frac{1}{\int\psi_{\uparrow
R}^{2}dx\,dy}\int\left[k_{L}\left(\frac{\partial\psi_{\downarrow I}}{\partial
x}-\frac{\partial\psi_{\downarrow R}}{\partial
y}\right)+\Omega\psi_{\downarrow R}\right]\psi_{\uparrow R}dx\,dy,$ (18a)
$\displaystyle\mu_{\downarrow}=$ $\displaystyle\frac{1}{\int\psi_{\downarrow
R}^{2}dx\,dy}\int\left[\frac{1}{2}\left(\frac{\partial\psi_{\downarrow
R}}{\partial x}\right)^{2}+\frac{1}{2}\left(\frac{\partial\psi_{\downarrow
R}}{\partial
y}\right)^{2}+\left(V_{2D}(x,y)+\beta\lvert\psi_{\uparrow}\rvert^{2}+\alpha\lvert\psi_{\uparrow}\rvert^{2}\right)\psi_{\downarrow
R}^{2}\right]dx\,dy$ $\displaystyle+\frac{1}{\int\psi_{\downarrow
R}^{2}dx\,dy}\int\left[k_{L}\left(\frac{\partial\psi_{\uparrow I}}{\partial
x}+\frac{\partial\psi_{\uparrow R}}{\partial y}\right)+\Omega\psi_{\uparrow
R}\right]\psi_{\downarrow R}dx\,dy,$ (18b)

From this we get the total energy is given by
$E_{num}=\sum_{j=\uparrow,\downarrow}\frac{\iint\left(E_{j}^{2C}+E_{j}^{SO}\right)\,dx\,dy}{\iint\psi_{jR}^{2}\,dx\,dy}$
can be extracted from the above chemical potential equations as:

$\displaystyle E_{\uparrow}^{2C}$
$\displaystyle=\left(\partial_{x}^{2}\psi_{\uparrow
R}\right)^{2}+\left(\partial_{y}^{2}\psi_{\uparrow
R}\right)^{2}+\left(V_{2D}+\alpha|\psi_{\uparrow}|^{2}/2+\beta|\psi_{\downarrow}|^{2}/2\right)\psi_{\uparrow
R}^{2}$ (19a) $\displaystyle E_{\downarrow}^{2C}$
$\displaystyle=\left(\partial_{x}^{2}\psi_{\downarrow
R}\right)^{2}+\left(\partial_{y}^{2}\psi_{\downarrow
R}\right)^{2}+\left(V_{2D}+\beta|\psi_{\uparrow}|^{2}/2+\alpha|\psi_{\downarrow}|^{2}/2\right)\psi_{\downarrow
R}^{2}$ (19b) $\displaystyle E_{\uparrow}^{SO}$
$\displaystyle=k_{L}\left(\psi_{\uparrow R}\partial_{x}\psi_{\downarrow
I}-\psi_{\uparrow R}\partial_{y}\psi_{\downarrow R}\right)+\psi_{\uparrow
R}\Omega\psi_{\downarrow R}$ (19c) $\displaystyle E_{\downarrow}^{SO}$
$\displaystyle=k_{L}\left(\psi_{\downarrow R}\partial_{x}\psi_{\uparrow
I}+\psi_{\downarrow R}\partial_{y}\psi_{\uparrow R}\right)+\psi_{\downarrow
R}\Omega\psi_{\uparrow R}$ (19d)

##  Appendix B Elements of the BdG matrix

$\displaystyle f(n_{\uparrow},n_{\downarrow})$
$\displaystyle=\frac{\left(k^{2}_{x}+k^{2}_{y}\right)}{2}+2\alpha
n_{\uparrow}+\beta
n_{\downarrow}-\frac{1}{2}\left[\alpha(n_{\uparrow}+n_{\downarrow})+\beta
n-\frac{n\Omega}{\sqrt{n_{\uparrow}n_{\downarrow}}}\right],$ (20a)
$\displaystyle g(n_{\uparrow},n_{\downarrow})$
$\displaystyle=\frac{\left(k^{2}_{x}+k^{2}_{y}\right)}{2}+\alpha
n_{\uparrow}+2\beta
n_{\downarrow}-\frac{1}{2}\left[\beta(n_{\uparrow}+n_{\downarrow})+\alpha
n-\frac{n\Omega}{\sqrt{n_{\uparrow}n_{\downarrow}}}\right],$ (20b)
$\displaystyle L_{13(24)}$
$\displaystyle=\beta\sqrt{n_{\uparrow}n_{\downarrow}}\pm
k_{L}\left(k_{x}-\mathrm{i}k_{y}\right)-\Omega,$ (20c) $\displaystyle
L_{31(42)}$ $\displaystyle=\beta\sqrt{n_{\uparrow}n_{\downarrow}}\pm
k_{L}\left(k_{x}+\mathrm{i}k_{y}\right)-\Omega,$ (20d)

##  Appendix C Coefficients of the BdG excitation spectrum

$\displaystyle b=$
$\displaystyle-\left(k_{x}^{2}+k_{y}^{2}\right)\left[2\left(k_{L}^{2}+\Omega\right)\alpha+\frac{1}{2}\left(k_{x}^{2}+k_{y}^{2}\right)\right]-2\Omega\left[(\alpha-\beta)+2\Omega\right],$
(21a) $\displaystyle c=$ $\displaystyle
2k_{L}k_{x}\Big{[}2\Omega\left[(\alpha-\beta)+2\Omega\right]-\left(k_{x}^{2}+k_{y}^{2}\right)(\beta-2\Omega)\Big{]},$
(21b) $\displaystyle d=$
$\displaystyle\frac{1}{16}\bigg{[}\left(k_{x}^{8}+k_{y}^{8}\right)+4P\left(k_{x}^{6}+k_{y}^{6}\right)+2Q\left(k_{x}^{4}+k_{y}^{4}\right)+k_{x}^{2}k_{y}^{2}\left[\left(k_{x}^{4}+k_{y}^{4}\right)+12P\left(k_{x}^{2}+k_{y}^{2}\right)+6k_{x}^{2}k_{y}^{2}+4Q\right]\bigg{]}+\left(Rk_{x}^{2}+Sk_{y}^{2}\right)$
(21c)

where

$\displaystyle P$ $\displaystyle=\alpha-2k_{L}^{2}+2\Omega$ (22a)
$\displaystyle Q$
$\displaystyle=2\left[\alpha^{2}-\beta^{2}+\alpha\left(6\Omega-4k_{L}^{2}\right)+2\beta\Omega+4\left(k_{L}^{2}-\Omega\right)^{2}\right]$
(22b) $\displaystyle R$
$\displaystyle=\Omega\left(\alpha+\beta-2k_{L}^{2}\right)(\alpha-\beta+2\Omega)$
(22c) $\displaystyle S$
$\displaystyle=\Omega(\alpha+\beta)\left(\alpha-\beta-2k_{L}^{2}+2\Omega\right)$
(22d)

## References

  * Anderson  _et al._ (1995) M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. Wieman, and E. A. Cornell, Science 269, 198 (1995).
  * Davis  _et al._ (1995) K. B. Davis, M. O. Mewes, M. R. Andrews, N. J. van Druten, D. S. Durfee, D. M. Kurn, and W. Ketterle, Phys. Rev. Lett. 75, 3969 (1995).
  * Bradley  _et al._ (1995) C. C. Bradley, C. A. Sackett, J. J. Tollett, and R. G. Hulet, Phys. Rev. Lett. 75, 1687 (1995).
  * Gerton  _et al._ (2000) J. M. Gerton, D. Strekalov, I. Prodan, and R. G. Hulet, Nature 408, 692 (2000).
  * Greiner  _et al._ (2002) M. Greiner, O. Mandel, T. Esslinger, T. W. Hänsch, and I. Bloch, Nature 415, 39 (2002).
  * Morsch and Oberthaler (2006) O. Morsch and M. Oberthaler, Rev. Mod. Phys. 78, 179 (2006).
  * Lewenstein  _et al._ (2007) M. Lewenstein, A. Sanpera, V. Ahufinger, B. Damski, A. Sen(De), and U. Sen, Adv. Phys. 56, 243 (2007).
  * Roati  _et al._ (2008) G. Roati, C. D’Errico, L. Fallani, M. Fattori, C. Fort, M. Zaccanti, G. Modugno, M. Modugno, and M. Inguscio, Nature 453, 895 (2008).
  * Chin  _et al._ (2010) C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, Rev. Mod. Phys. 82, 1225 (2010).
  * Griesmaier  _et al._ (2005) A. Griesmaier, J. Werner, S. Hensler, J. Stuhler, and T. Pfau, Phys. Rev. Lett. 94 (2005).
  * Lu  _et al._ (2011) M. Lu, N. Q. Burdick, S. H. Youn, and B. L. Lev, Phys. Rev. Lett. 107 (2011).
  * Aikawa  _et al._ (2012) K. Aikawa, A. Frisch, M. Mark, S. Baier, A. Rietzler, R. Grimm, and F. Ferlaino, Phys. Rev. Lett. 108, 210401 (2012).
  * Lin  _et al._ (2011) Y.-J. Lin, K. Jiménez-García, and I. B. Spielman, Nature 471, 83 (2011).
  * Galitski and Spielman (2013) V. Galitski and I. B. Spielman, Nature 494, 49 (2013).
  * Ho (1998) T.-L. Ho, Phys. Rev. Lett. 81, 742 (1998).
  * Papp  _et al._ (2008) S. B. Papp, J. M. Pino, and C. E. Wieman, Phys. Rev. Lett. 101, 040402 (2008).
  * Takeuchi  _et al._ (2010) H. Takeuchi, S. Ishino, and M. Tsubota, Phys. Rev. Lett. 105, 205301 (2010).
  * Sabbatini  _et al._ (2011) J. Sabbatini, W. H. Zurek, and M. J. Davis, Phys. Rev. Lett. 107, 230402 (2011).
  * Law  _et al._ (2010) K. J. H. Law, P. G. Kevrekidis, and L. S. Tuckerman, Phys. Rev. Lett. 105, 160405 (2010).
  * Wen  _et al._ (2013) L. Wen, Y. Qiao, Y. Xu, and L. Mao, Phys. Rev. A 87, 033604 (2013).
  * Aidelsburger  _et al._ (2013) M. Aidelsburger, M. Atala, M. Lohse, J. T. Barreiro, B. Paredes, and I. Bloch, Phys. Rev. Lett. 111, 185301 (2013).
  * Goldman  _et al._ (2010) N. Goldman, I. Satija, P. Nikolic, A. Bermudez, M. A. Martin-Delgado, M. Lewenstein, and I. B. Spielman, Phys. Rev. Lett. 105, 255302 (2010).
  * Wu  _et al._ (2016) Z. Wu, L. Zhang, W. Sun, X.-T. Xu, B.-Z. Wang, S.-C. Ji, Y. Deng, S. Chen, X.-J. Liu, and J.-W. Pan, Science 354, 83 (2016).
  * Seaman  _et al._ (2007) B. T. Seaman, M. Krämer, D. Z. Anderson, and M. J. Holland, Phys. Rev. A 75, 023615 (2007).
  * Andrianov and Moiseev (2014) S. N. Andrianov and S. A. Moiseev, Phys. Rev. A 90, 042303 (2014).
  * (26) R. Ravisankar, T. Sriraman, L. Salasnich, and P. Muruganandam, J. Phys. B 53, 195301 (2020).
  * (27) V. Achilleos, D. J. Frantzeskakis, P. G. Kevrekidis, and D. E. Pelinovsky, Phys. Rev. Lett. 110, 264101 (2013).
  * Achilleos  _et al._ (2013) V. Achilleos, J. Stockhofe, P. G. Kevrekidis, D. J. Frantzeskakis, and P. Schmelcher, Europhys. Lett. 103, 20002 (2013).
  * Jin  _et al._ (2014) J. Jin, S. Zhang, and W. Han, J. Phys. B 47, 115302 (2014).
  * Li and Sakaguchi (2013) B. Li and H. Sakaguchi, J. Low. Temp. Phys. 175, 243 (2013).
  * Cheng  _et al._ (2014) Y. Cheng, G. Tang, and S. K. Adhikari, Phys. Rev. A 89, 063602 (2014).
  * (32) H. He and Y. Zhang, Phys. Rev. A 103, 053322 (2021).
  * (33) S. Pal, A. Roy, and D. Angom, J. Phys. B: At. Mol. Opt. Phys. 50, 195301 (2017).
  * (34) S. Pal, A. Roy, and D. Angom, J. Phys. B: At. Mol. Opt. Phys. 51, 085302 (2018).
  * Martone  _et al._ (2012) G. I. Martone, Y. Li, L. P. Pitaevskii, and S. Stringari, Phys. Rev. A 86, 063621 (2012).
  * Zheng  _et al._ (2013) W. Zheng, Z.-Q. Yu, X. Cui, and H. Zhai, J. Phys. B 46, 134007 (2013).
  * Khamehchi  _et al._ (2014) M. A. Khamehchi, Y. Zhang, C. Hamner, T. Busch, and P. Engels, Phys. Rev. A 90, 063624 (2014).
  * Ji  _et al._ (2015) S.-C. Ji, L. Zhang, X.-T. Xu, Z. Wu, Y. Deng, S. Chen, and J.-W. Pan, Phys. Rev. Lett. 114, 105301 (2015).
  * Liao  _et al._ (2015) R. Liao, O. Fialko, J. Brand, and U. Zülicke, Phys. Rev. A 92, 043633 (2015).
  * (40) T. Ozawa, L. P. Pitaevskii, and S. Stringari, Phys. Rev. A 87, 063610 (2013).
  * (41) L. Chen, H. Pu, Z-Q. Yu, and Y. Zhang, Phys. Rev. A 95, 033616 (2017).
  * (42) Z-F. Yu and J-K. Xue, Europhys. Lett. 121, 20003 (2018).
  * (43) K. T. Geier, G. I. Martone, P. Hauke, and S. Stringari, arXiv:2102.02221v1 (2021).
  * (44) Q. Zhu1, C. Zhang, and B. Wu, Europhys. Lett. 100, 50003 (2012).
  * Pu and Bigelow (1998) H. Pu and N. P. Bigelow, Phys. Rev. Lett. 80, 1130 (1998).
  * Yu (2013) Z.-Q. Yu, Phys. Rev. A 87, 051606 (2013).
  * (47) S. Sahu and D. Majumder, J. Phys.B: At. Mol. Opt. Phys. 53, 095301 (2020).
  * Wang  _et al._ (2010) C. Wang, C. Gao, C.-M. Jian, and H. Zhai, Phys. Rev. Lett. 105, 160403 (2010).
  * Goldstein and Meystre (1997) E. V. Goldstein and P. Meystre, Phys. Rev. A 55, 2935 (1997).
  * Abad and Recati (2013) M. Abad and A. Recati, Eur. Phys. J. D 67, 148 (2013).
  * (51) P. Tommasini, E. J. V. de Passos, A. F. R. de Toledo Piza, M. S.Hussein, and E. Timmermans Phys. Rev. A 67, 023606 (2003).
  * (52) A. Recati and F. Piazza, Phys. Rev. B 99, 064505 (2019).
  * (53) Y. Li, L. P. Pitaevskii, and S. Stringari, Phys. Rev. Lett. 108, 225301 (2012).
  * (54) R. Ravisankar, T. Sriraman, and P. Muruganandam, AIP Conf. Proc. 2265, 030022 (2020).
  * (55) S. Bhuvaneswari, K. Nithyanandan, and P. Muruganandam, J. Phys. Commun. 2, 025008 (2018).
  * Bogolyubov (1947) N. N. Bogolyubov, J. Phys. (USSR) 11, 23 (1947).
  * Zilsel (1950) P. R. Zilsel, Phys. Rev. 79, 309 (1950).
  * Rickayzen (1959) G. Rickayzen, Phys. Rev. 115, 795 (1959).
  * 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).
  * (61) E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra,J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, D. Sorensen, LAPACK Users’ Guide, Society for Industrial and Applied Mathematics, Third Ed (1999).
  * Muruganandam and Adhikari (2009) P. Muruganandam and S. K. Adhikari, Comput. Phys. Commun. 180, 1888 (2009).
  * (63) L. E. Young-S., D. Vudragović, P. Muruganandam, S. K. Adhikari, and A. Balaž, Comput. Phys. Commun. 204, 209 (2016).
  * (64) R Ravisankar, D Vudragović, P Muruganandam, A Balaž, SK Adhikari, Comput. Phys. Commun. 259, 107657 (2021).
  * (65) P. Muruganandam, S. K. Adhikari, and A. Balaž, Comput. Phys. Commun. 264, 107926 (2021).


Research Report Fall 2020
Indian Institute of Technology Bombay
Department of Computer Science and Engineering

Deep Neural Matching Models for
Graph Retrieval

_Author:_

Yash Jain
<EMAIL_ADDRESS>Chitrank Gupta
<EMAIL_ADDRESS>
_Supervisors:_
Prof. Soumen Chakrabarti
<EMAIL_ADDRESS>
Prof. Abir De
<EMAIL_ADDRESS>
August 27, 2024

###### Abstract

Graph Retrieval has witnessed continued interest and progress in the past few
years. In this report, we focus on neural network based approaches for Graph
matching and retrieving similar graphs from a corpus of graphs. We explore
methods which can soft predict the similarity between two graphs. Later, we
gauge the power of a particular baseline (Shortest Path Kernel) and try to
model it in our product graph random walks setting while making it more
generalised.

###### Contents

  1. 1 Introduction
    1. 1.1 Objectives
    2. 1.2 Contributions
  2. 2 Deep Neural Matching Models for Graph Retrieval
    1. 2.1 Introduction
    2. 2.2 Notation
    3. 2.3 Method
    4. 2.4 Results
  3. 3 Datasets
    1. 3.1 SQuAD [6]
    2. 3.2 COX2
    3. 3.3 PTC
    4. 3.4 MUTAG
    5. 3.5 Clevr VQA
    6. 3.6 Pascal VOC Keypoints
  4. 4 Baselines
    1. 4.1 Non-trainable Heuristics based
      1. 4.1.1 Random Prediction Baseline
      2. 4.1.2 Graph Kernels
      3. 4.1.3 Shortest Path Kernel [1]
      4. 4.1.4 Reweighted Random Walk kernel [3]
    2. 4.2 Trainable Methods
      1. 4.2.1 GMN (Graph Matching Networks) [5]
  5. 5 Biased/ Trainable walks in GxNet
    1. 5.1 Introduction
    2. 5.2 Method
    3. 5.3 Results
  6. 6 Neural matching with Shortest path product graph
    1. 6.1 Introduction
    2. 6.2 Method
    3. 6.3 Results
    4. 6.4 Conclusion

###### List of Figures

  1. 4.1 Left: Siamese network using graph embedding method, Right: GMN model

###### List of Tables

  1. 2.1 Mean Average Precision Values across various baselines.
  2. 5.1 Mean Average Precision values of various approaches. The best performing approaches are emphasised, suggesting superior performance of GxNET. GMN (Graph Matching Networks), RRW (Reweighted Random walk kernel) SPK (shortest path kernel) $\phi$ refers to the trainable parameters of learning transition probabilities in the random walks
  3. 6.1 Mean Average Precision values. The best performing approaches are emphasised, suggesting the superior performance of GMN. Note that GMN (Graph Matching Networks), RW (Random walk kernel) SPK (shortest path kernel), SPRW (Shortest Path Random Walk) has three settings (a) unbiased random walks (b) biased random walks with $\alpha$ = 1 and (c) biased random walks with $\alpha$ = $inf$

##  Chapter 1 Introduction

###  1.1 Objectives

The objective of this study was to:

  * •

Improve the performance of the SIGIR paper "Deep Neural Matching Models for
Graph Retrieval" [4] by incorporating trainable random walks in GxNet

  * •

Implement new graph kernel baselines apart from the GxNet paper [4]

  * •

Survey new graph datasets for robust testing

  * •

Develop novel shortest path product graph approach to generalise Shortest path
kernel [2] in random walk setting

###  1.2 Contributions

During this study, we developed a thorough understanding of graph based
information retrieval systems. We learnt the pitfalls in various graphs based
IR systems and tried to improve the performance of GxNet by incorporating
trainable random walks in its architecture. Further, we tried to generalise
the Shortest Path Kernel approach by integrating its ideas into product graph
regime to match its performance with Shortest Path Kernel in random walk
framework on several datasets.

##  Chapter 2 Deep Neural Matching Models for Graph Retrieval

###  2.1 Introduction

The graph retrieval problem is to return ‘relevant’ or ‘good’ response graphs
from a corpus, given a query graph. The solution to this problem requires
scoring (and thereby rank) the corpus graphs with respect to a query graph
using some technique. Though in each dataset, the concept of relevance is
different, i.e. either semantic or structural and hence, a single heuristic
will not work for both types of such datasets. Most prevalent graph kernels
quantify the structural similarity between graphs, providing high similarity
scores to isomorphic graphs. That’s why this check of similarity between
graphs is also called isomorphism check. On the other hand, in many other
applications like drug active-sites identification or question answering from
documents, we may want a corpus graph to be prompted relevant if a subgraph of
the corpus graph is isomorphic to the query graph. Thus, GxNet proposes an E2E
supervised approach which implicitly learns such isomorphic subgraphs.

###  2.2 Notation

Given a set of query graphs
$\mathcal{Q}=\left\\{G_{q}=(V_{q},E_{q},F_{q},T_{q})\right\\}$ and a set of
corpus graphs $\mathcal{C}=\\{G_{c}=(V_{c},E_{c},F_{c},T_{c})\\}$ . Each query
graph $G_{q}$ is associated with a set of _good_ (relevant) corpus graphs

$\displaystyle\mathcal{C}_{q+}$
$\displaystyle=\\{G_{c}\in\mathcal{C}|\text{Relevance}(\mathcal{G}_{q},\mathcal{G}_{c})=1\\}$
(2.1) and a set of _bad_ (irrelevant) corpus graphs
$\displaystyle\mathcal{C}_{q-}$
$\displaystyle=\\{G_{c}\in\mathcal{C}|\text{Relevance}(\mathcal{G}_{q},\mathcal{G}_{c})=0\\}.$
(2.2)

For each query or corpus graph $G_{\bullet}$ we also observe a set of node
features $F_{\bullet}=\\{f_{u}\\}_{u\in V_{\bullet}}$ and a set of edge types
$T_{\bullet}=\\{t_{e}\\}_{e\in E_{\bullet}}$, where $\bullet$ may be one of
$q$ and $c$. The meaning of the underlying node features and edge types vary
across datasets.

##### Product graph

For each pair of query $G_{q}$ and corpus graph $G_{c}$, we construct the
tensor product graph ${\mathcal{G}}^{qc}=G_{q}\times G_{c}$ :

$\displaystyle{\mathcal{G}}^{qc}=G_{q}\star
G_{c}=({\mathcal{V}}^{qc},{\mathcal{E}}^{qc},{\mathcal{F}}^{qc},{\mathcal{T}}^{qc}),\;\text{where}$
(2.3)
$\displaystyle{\mathcal{V}}^{qc}\\!=\\!V_{q}\\!\times\\!V_{c},\;{\mathcal{E}}^{qc}\\!=\\!E_{q}\\!\times\\!E_{c},\;{\mathcal{F}}^{qc}=F_{q}\times
F_{c},\;{\mathcal{T}}^{qc}$ $\displaystyle=T_{q}\times T_{c}.$

Edges are included in the product graph $\mathcal{G}^{qc}$ according to the
rule

$\displaystyle((u,u^{\prime}),(v,v^{\prime}))\in{\mathcal{E}}^{qc}\\!\iff\\!(u,v)\in
E_{q}\ \text{and}\ (u^{\prime},v^{\prime})\in{E_{c}};$ (2.4)

the node feature for a node $(u,u^{\prime})$ in $\mathcal{G}^{qc}$ are set as
follows:

$\displaystyle f_{(u,u^{\prime})}=(f_{u},f_{u^{\prime}})\ \forall u\in
G_{q},u^{\prime}\in G_{c};$ (2.5)

and the edge labels in $\mathcal{G}^{qc}$ are set as follows:

$\displaystyle
t_{(u,u^{\prime}),(v,v^{\prime})}=(t_{u,v},t_{u^{\prime},v^{\prime}})\
\text{if }(u,v)\in E_{q}\ \text{and}\ (u^{\prime},v^{\prime})\in E_{c}$ (2.6)

###  2.3 Method

GxNet computes an embedding $\bm{z}_{u,u^{\prime}}$ for each node
$(u,u^{\prime})$ in the product graph $\mathcal{G}^{qc}$. In addition, GxNet
also computes an embedding $\bm{r}_{ee^{\prime}}$ for each edge in the product
graph based on the edge features $\bm{t}_{e},\bm{t}_{e^{\prime}}$. More
specifically, we have:

$\displaystyle\bm{z}_{u,u^{\prime}}$
$\displaystyle=\bm{g}_{\text{node}}(\bm{f}_{u},\bm{f}_{u^{\prime}})$ (2.7)
$\displaystyle\bm{r}_{e,e^{\prime}}$
$\displaystyle=\bm{g}_{\text{edge}}(\bm{t}_{e},\bm{t}_{e^{\prime}})$ (2.8)

In the next step, GxNet constructs the embedding vector for each (source node,
edge, target node) triple in the product graph $\mathcal{G}^{qc}$,

$\displaystyle\text{denoted}\quad\tau$
$\displaystyle=\big{(}(u,u^{\prime})\overset{ee^{\prime}}{\longrightarrow}(v,v^{\prime})\big{)},$
(2.9) $\displaystyle\text{as}\quad\bm{h}_{\tau}$
$\displaystyle=\bm{\rho}(\bm{z}_{u,u^{\prime}},\bm{r}_{e,e^{\prime}},\bm{z}_{v,v^{\prime}}),$
(2.10)

which in turn gives us a score or scalar weight of $\tau$, given by

$\displaystyle s_{\tau}=\sigma(\bm{W}_{s,h}^{\top}\bm{h}_{\tau}),$ (2.11)

where $s_{\tau}$ indicates an affinity score for the edge $\tau$. Here
$\bm{g}_{\text{node}}$, $\bm{g}_{\text{edge}}$, $\bm{\rho}$, $\sigma$ are
suitable standard networks with appropriate non-linearities. The parameter
$\theta$ stands for the set of all the trainable parameters in these networks.
The final score $y_{\theta}(G_{q},G_{c})$ is calculated by aggregating the
scores $s_{\tau}$ of a Subgraph $S$ in the product graph $\mathcal{G}_{qc}$.
The idea is that if $G_{q}$ and $G_{c}$ are positively related then there will
be a subgraph in both of them (and subsequently $\mathcal{G}_{qc}$) which
highly ‘matches’ with each other. The neural networks hopes to learn this
latent substructure.

Since finding such a $S$ is tedious and almost impractical GxNet employs
Random walks of length $K$ (a hyperparameter) to sample a part of the product
graph and get an aggregated score with the hope of sampling the similar
subgraph in $G_{q}$ and $G_{c}$. Formally, for each $k$ $\epsilon$ [$K$], it
performs $n$ random walks with restarts in the product graph. The start
probability of all nodes is kept equal. At any node, transition probabilities
for all edges are equal. After sampling a random walk, it samples the subgraph
as the induced subgraph of a random walk. After sampling subgraphs
$S_{i},i\leq n$ computes the approximate maximum of the aggregated affinity
scores in the following way.

$\displaystyle\text{score}_{k}(G_{q},G_{c})=\sum_{i\in[n]}\dfrac{\sum_{\tau\in\mathcal{S}_{i}||\mathcal{S}_{i}|=k}s_{\tau}\exp(\sum_{\tau\in\mathcal{S}_{i}||\mathcal{S}_{i}|=k}s_{\tau})}{\sum_{i\leq
n}\exp(\sum_{\tau\in\mathcal{S}_{i}||\mathcal{S}_{i}|=k}s_{\tau})}$ (2.12)

Finally, GxNet compute the average affinity over the different lengths of the
walk, penalized by the length of path $k$. Hence, the final confidence score
of the relevance of $G_{q}$ with $G_{c}$ is

$\displaystyle
y_{\theta}(G_{q},G_{c})=\frac{1}{K}\sum_{k\in[K]}\frac{\text{score}_{k}(G_{q},G_{c})}{k},$
(2.13)

To train the parameter $\theta$, GxNet uses pairwise ranking loss function

$\displaystyle\underset{\theta}{\text{argmin}}\sum_{q\in\mathcal{Q}}\sum_{\begin{subarray}{c}G_{c}\in\mathcal{C}_{q-},\\\
G_{c^{\prime}}\in\mathcal{C}_{q+}\end{subarray}}\\!\\!\\!\text{ReLU}\Bigl{(}\Delta+y_{\theta}(G_{q},G_{c})-y_{\theta}(G_{q},G_{c^{\prime}})\Bigr{)}$
(2.14)

where $\Delta$ is a tuned margin and ReLU is $\max\\{\cdot,0\\}$.

###  2.4 Results

GxNet was evaluated on two datasets, i.e. VQA and SQuAD, and was compared
against the baselines GMN and RRWM. Please refer to the table 2.1 for detailed
results.

Dataset | GxNet | GMN | RRW |
---|---|---|---|---
| VQA (Clevr)
---
0.98 | 0.38 | 0.65 |
SQuAD | 0.34 | 0.27 | 0.37 |
COX2 | 0.31 | 0.35 | 0.37 |
Table 2.1: Mean Average Precision Values across various baselines.

Preliminary experiments suggested that GxNet showed great promise in graph
neural matching and was highly resilient to noise in the data compared to its
counterparts. We will delve in further about GxNet in future chapters.

##  Chapter 3 Datasets

###  3.1 SQuAD [6]

Stanford Question Answering Dataset (SQuAD) is a reading comprehension
dataset, consisting of questions posed by crowdworkers on a set of Wikipedia
articles, where the answer to every question is a segment of text, or span,
from the corresponding reading passage, or the question might be unanswerable.
To make it amiable for the graph retrieval task, each natural language
sentence was converted into its dependency parse using spacy dependency
parser. The corpus graph was made up of the rest of the spans from the reading
passage.

###  3.2 COX2

COX2 is a molecular dataset where every graph is classified into one of two
classes.

###  3.3 PTC

PTC is a molecular dataset where molecules are classified according to their
carcinogenicity on rodents. This dataset is classified into 4 datasets
PTC_FM,PTC_FR,PTC_MM and PTC_MR where M/F stands for male/female and M/R
stands for mice/rodent.

###  3.4 MUTAG

This dataset is a molecular dataset where molecules are classified according
to their mutagenicity on bacteria.

###  3.5 Clevr VQA

The Clevr VQA benchmark consists of an image corpus where each image is
associated with a scene graph. The nodes of a scene graph have attributes
concerning shape, color, texture and size. Edges of a scene graph represent
spatial relations like ‘left of’ and ‘behind’. Clevr queries are originally in
natural language and usually parsed into a graph. The query graphs were
generated by placing five objects randomly on the scene. Object attributes
(shape, color, size, material) are sampled uniformly randomly from predefined
choices. Next, spatial relations (left, right, front, behind) are computed
between object pairs. The label of an edge is a 1-hot vector of length 4,
corresponding to the above relations. Positive corpus graphs are generated by
adding more objects to the scene whose attributes are decided uniformly
randomly, and edges are added by computing relationships between object pairs
as before. Negative corpus graphs are sampled from the positive samples of
other queries since the probability of a corpus graph matching two query
graphs well is very low. To test the resilience of various algorithms, noise
is added to the color (RGB ∈ R 3 ) and size (∈ R) node attributes. Each
component is in [0, 1], to which an independent Gaussian noise with varying
standard deviation was added.

###  3.6 Pascal VOC Keypoints

Pascal Visual Object Classification challenge consist of 20 image classes,
where each image was parsed into an image graph using keypoints as nodes. The
nodes signify the VGG16 features of the bounding box of the part of image
local to the keypoint. The dataset was later rejected for the experiments as
the node features within the classes does not make sense semantically, i.e.,
nodes of an image-graph of aeroplane cannot be compared to the nodes of a car
but only be compared to an another instance of aeroplance image-graph. This
was unsuitable for our proposed task of deep graph matching with a query and a
set of corpus graphs.

##  Chapter 4 Baselines

###  4.1 Non-trainable Heuristics based

####  4.1.1 Random Prediction Baseline

Here we arbitrarily predict the similarity between the graphs.

####  4.1.2 Graph Kernels

A kernel is a function that takes two vectors and returns the similarity
between these two vectors as the dot product between some high dimensional
mapping of the two vectors. The two conditions equivalent to a kernel function
are symmetricity and positive semi-definite property. In context of graph
kernels, the higher dimensional feature space most often corresponds to the
decomposition of the graph into some specific substructures like walks, paths,
sub-graphs which may or may not be further refined. Here we consider two
different graph kernels- reweighted random walk kernel and shortest path
kernel.

####  4.1.3 Shortest Path Kernel [1]

This kernel essentially calculates the similarity between the graphs by
comparing the two graphs on the basis of the no of shortest paths of some
specific length and connecting two nodes of some specific node labels.
Formally let $\phi_{(}n,l_{u},l_{v})(G)$ denote the number of shortest paths
in graph G of length $n$ that exist between two nodes with node labels $l_{u}$
and $l_{v}$. Then we have,

$\begin{split}\phi_{(n,l_{u},l_{v})}(G)=\lvert(u,v):\mathtt{label}(u)&=l_{u},\mathtt{label}(v)=l_{v}\\\
&\mathtt{shortest\\_path}(u,v)=n\;\forall u,v\in V_{G}\rvert\end{split}$ (4.1)

We do this for all such possible $k$,$l_{u}$ and $l_{v}$ which can be done in
polynomial time thanks to algorithms like Floyd-warshall. Once done this, then
to compute the similarity of some graph in the corpus $G_{c}$ with that a
query graph $G_{q}$ one can do the following.

$\begin{split}k_{shortest\\_path}(G_{q},G_{c})=\sum_{(n,l_{u},l_{v})}\phi_{(n,l_{u},l_{v})}(G_{q}).\phi_{(n,l_{u},l_{v})}(G_{c})\end{split}$
(4.2)

One can rewrite this formula as the dot product between two vectors, each
representing the high dimensional vector mapping of graphs.

$\begin{split}k_{shortest\\_path}(G_{q},G_{c})=\Phi(G_{q}).\Phi(G_{c})\end{split}$
(4.3)

Normalised shortest path kernel[7]\- One can normalise the vectors $\Phi(G)$
and obtain the normalised shortest path kernel values.

$\begin{split}k_{Norm\\_shortest\\_path}(G_{q},G_{c})=\frac{\Phi(G_{q}).\Phi(G_{c})}{\left\lvert\Phi(G_{q})\right\rvert_{2}\left\lvert\Phi(G_{c})\right\rvert_{2}}\end{split}$
(4.4)

####  4.1.4 Reweighted Random Walk kernel [3]

The main idea in a random walk kernel for graph similarity is to count the
number of matching walks between the two graphs where the nodes are matched on
the basis of their node labels. But one can obtain this kernel value in the by
using the following formulation.

$\begin{split}k_{random\\_walk}(G_{q},C_{c})=\sum\limits_{i,j=1}^{\lvert
V_{x}\rvert}[\sum\limits_{n=0}^{N}\lambda_{n}A_{x}^{n}]_{ij}\end{split}$ (4.5)

But they are not effective in realf-life datasets when there are a lot of
noisy and outlier nodes present in the association graph of two similar
graphs. To solve this, a reweighted random walk kernel iteratively updates the
confidence in the candidate correspondences between nodes of two graphs,
hopefully becoming robust against noisy and outlier nodes.

###  4.2 Trainable Methods

####  4.2.1 GMN (Graph Matching Networks) [5]

This method outputs the similarity score by using a cross graph attention
mechanism. This is in contrast to siamese networks where the two graphs are
independently mapped to two vectors and then compared. In GMN, the information
of how similar or dis-similar a node in one graph is with every other node in
the other graph is also used while calculating the node embeddings.

Figure 4.1: Left: Siamese network using graph embedding method, Right: GMN
model

##  Chapter 5 Biased/ Trainable walks in GxNet

###  5.1 Introduction

The SIGIR GxNet paper was based on Uniform Random walks to sample the subgraph
within the product graph $\mathcal{G}_{qc}$. In a scenario where the product
graph is very huge there might be a situation where the uniform random walk
completely misses the target subgraph resulting in a low score and thereby
poor ranking. To avoid such a situation, we thought of training the walk paths
instead of doing a random sampling. We changed the start probability and
transition probabilities to depend on node or edge features in simple
parametric ways. Consider the edge $(u,u′)\to(v,v′)$ in the product graph
connecting $(u,u′)$ to $(v,v′)$. The corresponding node and edge features are
$f_{(u,u^{\prime})}$, $t_{(u,u^{\prime}),(v,v^{\prime})}$ and
$f_{(v,v^{\prime})}$. For random walk over the product graph, we designate the
weight of the edge as
$g_{\phi}(f_{(u,u^{\prime})},t_{(u,u^{\prime}),(v,v^{\prime})},f_{(v,v^{\prime})})$.
where $g_{\phi}$ is a neural network denoting a new set of non-linear
parameters $\phi$. Edge weights are then suitably scaled per-node to make the
walk probability matrix stochastic.

###  5.2 Method

Following the notation of GxNet, in the trainable walks the weight of the edge
is trained over a neural network whose parameters are denoted by $\phi$. Thus,
in the trainable walks the pairwise ranking loss function would become:

$\displaystyle\mathcal{L}(G,C,C^{\prime};\theta,\phi)=\sum_{G_{q}}\sum_{G_{c},G_{c}^{\prime}}\text{ReLU}(\Delta+y_{\theta,\phi}(G_{q},G_{c^{\prime}})-y_{\theta,\phi}(G_{q},G_{c}))$
(5.1)

We break the down the training process in three possible ways:

  * •

Train both $\theta$ and $\phi$ simultaneously

  * •

Train $\theta$ and $\phi$ alternatively

  * •

Train $\phi$ first and then train $\theta$ while keeping the parameters of
$\phi$ fixed

Along with the training procedure, we also experimented with the score
aggregation function for a random walk across two formulations:

  * •

MAX formulation: Selecting the max score across the $n$ walks of
$k\epsilon[K]$ lengths

$\displaystyle y_{\theta,\phi}(G_{q},G_{c},k)$
$\displaystyle=\max_{\begin{subarray}{c}\mathcal{S}\sim
P_{\phi}(\mathcal{S}|\mathcal{G}^{qc})|\mathcal{S}|=k\end{subarray}}(y_{\theta}(G_{q},G_{c},S))$
(5.2)
$\displaystyle\approx\frac{1}{\alpha}\log(\sum_{\begin{subarray}{c}\mathcal{S}\sim
P_{\phi}(\mathcal{S}|\mathcal{G}^{qc})|\mathcal{S}|=k\end{subarray}}\exp(\alpha
y_{\theta}(G_{q},G_{c},S)))$ (5.3)
$\displaystyle=\frac{1}{\alpha}\log(N\times\sum_{S||S|=k}\exp(\alpha\times
y_{\theta}(G_{q},G_{c},S))\frac{P_{\phi}(\mathcal{S}|\mathcal{G}^{qc})}{P_{\text{BiasedWalk}}(\mathcal{S}|\mathcal{G}^{qc})}*P_{\text{BiasedWalk}}(\mathcal{S}|\mathcal{G}^{qc}))$
(5.4)
$\displaystyle=\frac{1}{\alpha}\log(\sum_{S\sim\text{BiasedWalk}|S|=k}\exp(\alpha\times
y_{\theta}(G_{q},G_{c},S))\frac{P_{\phi}(\mathcal{S}|\mathcal{G}^{qc})}{P_{\text{BiasedWalk}}(\mathcal{S}|\mathcal{G}^{qc})})$
(5.5)

where $\alpha$ is a hyper-parameter

  * •

AVERAGE formulation: Averaging the scores across $n$ walks of $k\epsilon[K]$
lengths.
For a fixed walk length k,

$\displaystyle y_{\theta,\phi}(G_{q},G_{c},k)$
$\displaystyle=\mathbb{E}_{\begin{subarray}{c}\mathcal{S}\sim
P_{\phi}(\mathcal{S}|\mathcal{G}^{qc})\\\
|\mathcal{S}|=k\end{subarray}}[y_{\theta}(G_{q},G_{c},S)]$ (5.6)
$\displaystyle=\sum_{S:|S|=k}y_{\theta}(G_{q},G_{c},S)P_{\phi}(\mathcal{S}|\mathcal{G}^{qc})$
(5.7)
$\displaystyle=\sum_{S:|S|=k}y_{\theta}(G_{q},G_{c},S)\frac{P_{\phi}(\mathcal{S}|\mathcal{G}^{qc})}{P_{\text{BiasedWalk}}(\mathcal{S}|\mathcal{G}^{qc})}*P_{\text{BiasedWalk}}(\mathcal{S}|\mathcal{G}^{qc})$
(5.8)
$\displaystyle=\frac{1}{N}\times\sum_{\begin{subarray}{c}S\sim\text{BiasedWalk}\\\
|S|=k\end{subarray}}y_{\theta}(G_{q},G_{c},S)\frac{P_{\phi}(\mathcal{S}|\mathcal{G}^{qc})}{P_{\text{BiasedWalk}}(\mathcal{S}|\mathcal{G}^{qc})}$
(5.9)

Finally aggregating over walk lengths,

$\displaystyle
y_{\theta,\phi}(G_{q},G_{c})=\frac{1}{K}\sum\limits_{k\epsilon[K]}\lambda^{k}y_{\theta,\phi}(G_{q},G_{c},k)$
(5.10)

where $\lambda$ is a hyper-parameter, we kept $\lambda$ = 1 unless specified.

###  5.3 Results

The trainable walks certainly beat the uniform walks in the performance.
Though, when compared to other baselines, the trainable walks fall a little
short. In the next section we tried to improve upon this gap by taking a new
approach of shortest path product graphs.

Dataset |  | GxNET
---
(Unbiased walks)
| GxNET
---
(Biased walks)
fixed $\phi$
| GxNET
---
(Biased walks)
trainable $\phi$
GMN | RRW | SPK
| VQA-Clevr
---
(0.10 noise)
0.98 | 0.99 | 0.99 | 0.38 | 0.65 | 0.64
SQuAD | 0.34 | 0.41 | 0.42 | 0.27 | 0.37 | 0.29
COX2 | 0.31 | 0.34 | 0.36 | 0.35 | 0.37 | 0.25
Table 5.1: Mean Average Precision values of various approaches. The best
performing approaches are emphasised, suggesting superior performance of
GxNET. GMN (Graph Matching Networks), RRW (Reweighted Random walk kernel) SPK
(shortest path kernel) $\phi$ refers to the trainable parameters of learning
transition probabilities in the random walks

##  Chapter 6 Neural matching with Shortest path product graph

###  6.1 Introduction

After observing the performance of shortest path kernel on SQuAD dataset, we
decided to include the idea of shortest path kernel in our product graph
setting. The motive behind this was that in a random walk in a product graph
of original graphs the maximum hop that the walk can reach is limited while in
a product graph of shortest graphs, all nodes are within reach and the
information of the distance is also retained. We call this approach Shortest
Path Random Walk (SPRW).

###  6.2 Method

Give a query graph $G_{q}$ and corpus graph $G_{c}$, we first convert it into
its shortest path graphs $G_{qs}$ and $G_{cs}$ respectively. Note that the
shortest graph is a complete graph with the edges being the distance of the
two nodes in the original graphs. Thus evidently, if two nodes are not
connected in the original graph then in the shortest graph they will be
connected with a edge distance of $inf$ (infinity). Further, the tensor
product of $G_{qs}$ and $G_{cs}$, $\mathcal{G}_{qcs}$ was constructed.

To test out the advantages and the potential of shortest product graphs we
experimented with the biased walks by carefully curating the transition
probability to first match with the shortest path kernel baseline and then
surpass it in performance.

##### Start probability

Instead of a uniform random selection of the start node, we enforced the start
nodes to be the matching nodes in the query and corpus graphs, i.e. Give a
start node $(u,u^{\prime})$ in the product graph $\mathcal{G}_{qcs}$, its
start probability $P_{(u,u^{\prime})}$ will be

$\displaystyle
P_{(u,u^{\prime})}=\frac{\mathbbm{1}(l(u)=l(u^{\prime}))+\epsilon\mathbbm{1}(l(u)\neq
l(u^{\prime}))}{\sum_{(u,u^{\prime})\epsilon\mathcal{G}_{qcs}}{\mathbbm{1}(l(u)=l(u^{\prime}))+\epsilon\mathbbm{1}(l(u)\neq
l(u^{\prime}))}}$ (6.1)

where $l(u)$ = $\mathtt{label}(u)$ and $l(u^{\prime})$ =
$\mathtt{label}(u^{\prime})$ in the original graphs $G_{qs}$ and $G_{cs}$

##### Transition probability

Similarly, transitions to the matching nodes with the same shortest distance
in the query and corpus graphs was favoured over other transitions. Moreover,
by tweaking the shortest distance limit we could vary our method from shortest
path kernel to a more general setting. Formally, for a transition
$(u,u^{\prime})\to(v,v^{\prime})$ the transition probability
$P_{(u,u^{\prime})\to(v,v^{\prime})}$ is

$\displaystyle
P_{(u,u^{\prime})\to(v,v^{\prime})}=\frac{\mathbbm{1}(l(u)=l(u^{\prime})).\mathbbm{1}(l(v)=l(v^{\prime})).\mathbbm{1}(\texttt{sp}(u,v)=\texttt{sp}(u^{\prime},v^{\prime})).\mathbbm{1}(\texttt{sp}(u,v)\leq\alpha)}{\sum_{(u,u^{\prime})\to(v,v^{\prime})\epsilon\mathcal{G}_{qcs}}\mathbbm{1}(l(u)=l(u^{\prime})).\mathbbm{1}(l(v)=l(v^{\prime})).\mathbbm{1}(\texttt{sp}(u,v)=\texttt{sp}(u^{\prime},v^{\prime})).\mathbbm{1}(\texttt{sp}(u,v)\leq\alpha)}$
(6.2)

where $\texttt{sp}(u,v)$ is the shortest distance between nodes $u$ and $v$ in
$G_{qs}$ and $\alpha$ is a hyper-parameter denoting the maximum distance of
shortest path to consider for hopping. Hence, if $\alpha$ = $inf$, the above
formulation reduces to the shortest path kernel.

Note that there is no training involved in SPRW.

###  6.3 Results

Dataset |  | Random
---
Prediction
SPK | GMN |  | SPRW
---
(unbiased)
| SPRW
---
(biased)
$\texttt{SP}_{max}=1$
| SPRW
---
(biased)
$\texttt{SP}_{max}$=inf)
MUTAG | 0.56 | 0.70 | 0.78 | 0.648 | 0.675 | 0.67
PTC_MM | 0.53 | 0.55 | 0.57 | 0.54 | 0.54 | 0.54
PTC_FM | 0.52 | 0.55 | 0.57 | 0.35 | — | —
Table 6.1: Mean Average Precision values. The best performing approaches are
emphasised, suggesting the superior performance of GMN. Note that GMN (Graph
Matching Networks), RW (Random walk kernel) SPK (shortest path kernel), SPRW
(Shortest Path Random Walk) has three settings (a) unbiased random walks (b)
biased random walks with $\alpha$ = 1 and (c) biased random walks with
$\alpha$ = $inf$

###  6.4 Conclusion

From our research we draw following conclusions:

  * •

GxNet (biased walks with trainable $\phi$ ) performs better than GxNet (biased
walks with fixed $\phi$ ) which in turn is better than GxNet based on uniform
random walks (the original SIGIR paper) [4]

  * •

In many cases normalised Shortest Path Kernel beats the Gxnet approach as was
found with SQuAD dataset. Un-normalised Shortest Path Kernel still remains
beaten by GxNet approach

  * •

On TUDatasets, Shortest Path Kernel performs similar to Random Prediction
denoting that it fails to understand the semantic relations in these dataset

## References

  * [1] Borgwardt, K. M., and Kriegel, H.-P. Shortest-path kernels on graphs.  In Proceedings of the Fifth IEEE International Conference on Data Mining (USA, 2005), ICDM ’05, IEEE Computer Society, p. 74–81.
  * [2] Borgwardt, K. M., and Kriegel, H. P. Shortest-path kernels on graphs.  In Fifth IEEE International Conference on Data Mining (ICDM’05) (2005), pp. 8 pp.–.
  * [3] Cho, M., Lee, J., and Lee, K. M. Reweighted random walks for graph matching.  In Computer Vision – ECCV 2010 (Berlin, Heidelberg, 2010), K. Daniilidis, P. Maragos, and N. Paragios, Eds., Springer Berlin Heidelberg, pp. 492–505.
  * [4] Goyal, K., Gupta, U., De, A., and Chakrabarti, S. Deep neural matching models for graph retrieval.  In Proceedings of the 43rd International ACM SIGIR Conference on Research and Development in Information Retrieval (New York, NY, USA, 2020), SIGIR ’20, Association for Computing Machinery, p. 1701–1704.
  * [5] Li, Y., Gu, C., Dullien, T., Vinyals, O., and Kohli, P. Graph matching networks for learning the similarity of graph structured objects, 2019.
  * [6] Rajpurkar, P., Jia, R., and Liang, P. Know what you don’t know: Unanswerable questions for squad, 2018.
  * [7] Siglidis, G., Nikolentzos, G., Limnios, S., Giatsidis, C., Skianis, K., and Vazirgiannis, M. Grakel: A graph kernel library in python, 2020.