# Quantum phases and spectrum of collective modes in a spin-1 BEC with spin-
orbital-angular-momentum coupling

Paramjeet Banger<EMAIL_ADDRESS>Department of Physics, Indian
Institute of Technology Ropar, Rupnagar-140001, Punjab, India  Rajat
<EMAIL_ADDRESS>Department of Physics, Indian Institute of
Technology Ropar, Rupnagar-140001, Punjab, India  Arko Roy
<EMAIL_ADDRESS>School of Physical Sciences, Indian Institute of
Technology Mandi, Mandi-175075 (H.P.), India  Sandeep Gautam
<EMAIL_ADDRESS>Department of Physics, Indian Institute of Technology
Ropar, Rupnagar-140001, Punjab, India

###### Abstract

Motivated by the recent experiments [Chen et al., Phys. Rev. Lett 121, 113204
(2018), Chen et al., Phys. Rev. Lett. 121, 250401 (2018)], we investigate the
low-lying excitation spectrum of the ground-state phases of spin-orbital-
angular-momentum-coupled (SOAM-coupled) spin-1 condensates. At vanishing
detuning, a ferromagnetic SOAM-coupled spin-1 BEC can have two ground-state
phases, namely coreless and polar-core vortex states, whereas an
antiferromagnetic BEC supports only polar-core vortex solution. The angular
momentum per particle, longitudinal magnetization, and excitation frequencies
display discontinuities across the phase boundary between the two phases. The
low-lying excitation spectrum evaluated by solving the Bogoliubov-de-Gennes
equations is marked by avoided crossings and hence the hybridization of the
spin and density channels. The spectrum is further confirmed by the dynamical
evolution of the ground state subjected to a perturbation suitable to excite a
density or a spin mode and a variational analysis for the density-breathing
mode.

##  I Introduction

The experimental realization of spin-orbit (SO) coupling marked an important
milestone in the field of quantum degenerate Bose gases Lin  _et al._ (2011);
Campbell  _et al._ (2016); Luo  _et al._ (2016). The SO coupling in these
experiments, coupling the spin and linear momentum of electrically neutral
bosons, is created by controlling the interaction between atoms and light Lin
_et al._ (2011); Campbell  _et al._ (2016); Luo  _et al._ (2016); Goldman  _et
al._ (2014). The rich ground-state phase diagram of SO-coupled spin-1 BECs
besides having plane-wave, standing-wave, and zero-momentum phases Wang  _et
al._ (2010); Martone  _et al._ (2016) also admit half-quantum vortex
Ramachandhran  _et al._ (2012), vortex-lattice states Ruokokoski  _et al._
(2012), etc. Apart from the ground-state phase diagram, collective excitations
in trapped SO-coupled BECs, another aspect of fundamental interest Pethick and
Smith (2008); Pitaevskii and Stringari (2016), have been studied
experimentally Khamehchi  _et al._ (2014); Ji  _et al._ (2015) as well as
theoretically Chen  _et al._ (2017); Geier  _et al._ (2021) in harmonically-
trapped SO-coupled pseudospinor BECs. Recently, we studied the collective
excitations in a quasi-one-dimensional SO-coupled spin-1 BEC with
antiferromagnetic interactions at zero and finite temperatures Rajat  _et al._
(2022).

For the last few years, there has been a growing interest in coupling the
orbital angular momentum of atoms’ center-of-mass with their internal spin
states using a pair of copropagating Laguerre-Gaussian laser beams with
opposite winding numbers. Commonly known as the spin-orbital-angular-momentum
(SOAM) coupling, this feature has been independently demonstrated by two
experimental groups by coupling two Zhang  _et al._ (2019) or three magnetic
sub-levels of $F=1$ manifold of 87Rb atoms Chen  _et al._ (2018a, b), thus
affirming the validity of earlier theoretical proposals Hu  _et al._ (2015);
DeMarco and Pu (2015); Sun  _et al._ (2015); Qu  _et al._ (2015); Chen  _et
al._ (2016); Peng  _et al._ (2022). In the context of SOAM-coupled
pseudospin-1/2 BECs, polarized and zero-momentum phases have been observed
experimentally. Besides these, a stripe, an annular-stripe, a two-vortex
molecule, and a vortex-antivortex molecule phases have also been studied
theoretically Qu  _et al._ (2015); Duan  _et al._ (2020); Chen  _et al._
(2020a); Chiu  _et al._ (2020). Theoretical studies on the effects of the
ring-trapping potential on the annular-stripe phase in SOAM-coupled
pseudospin-1/2 condensate have also been carried out Sun  _et al._ (2015);
Bidasyuk  _et al._ (2022).

Along with studies on equilibrium ground-state phase diagrams, spectroscopic
studies have been carried out on SOAM-coupled pseudospin-1/2 BEC Chen  _et
al._ (2020a); Vasić and Balaž (2016); Chen  _et al._ (2020b). In particular,
the low-lying excitation spectrum, including breathing and dipole modes, have
been studied for the half-skyrmion and vortex-antivortex phases Vasić and
Balaž (2016); Chen  _et al._ (2020a). Additionally, the ground-state phases
and excitation spectrum have been studied for pseudospin-1/2 BEC with higher
order SOAM coupling Chen  _et al._ (2020b).

However, the detailed phase diagrams and excitation spectrums of SOAM-coupled
spin-1 BECs with polar and ferromagnetic spin-exchange interactions have not
been yet theoretically studied. This sets the stage for the present work. With
inspiration drawn from the experimental research reported in Chen  _et al._
(2018b) and an aim to bridge the research gap, our objective is to study the
excitation spectrum of the ground-state phases observed in SOAM-coupled spin-1
BECs. The excitation spectrum, calculated by solving the Bogoliubov-de Gennes
(BdG) equations, is supported by the time evolution of the expectation of the
physical observables with an aptly chosen perturbation being added to the
Hamiltonian at time $t=0$. For the sake of comprehensiveness, we additionally
employ the variational method to analytically calculate the frequency of the
density-breathing mode.

The paper is organized as follows. In Sec. II, we present the Hamiltonian
describing an SOAM-coupled spin-1 BEC in cylindrical coordinates and the
reduction to a quasi-two-dimensional (quasi-2D) formulation through a set of
coupled Gross-Pitaevskii equations (GPEs). In Sec. III, we discuss the ground-
state phases of SOAM-coupled ferromagnetic and polar BECs in the limit of
vanishing detuning. In Sec. IVA, we discuss the spectrum of the noninteracting
SOAM-coupled spin-1 BEC, and follow it with the collective excitations of the
interacting SOAM-coupled spin-1 BECs in IVB. In Sec. IVC, we study real-time
dynamics of the perturbed ground state to illustrate the ensuing dynamics in
the density and spin channels. In Sec. IVD, the variational method to study a
few low-lying modes is discussed, which is followed by the summary of key
results in Sec. V.

##  II Model

In this work, we consider SOAM-coupled spin-1 BECs in which the orbital
angular momentum of the center of the mass of the atoms is synthetically
coupled to their internal spin states Chen  _et al._ (2018a); Zhang  _et al._
(2019). In the cylindrical coordinate system, the non-interacting (single-
particle) part of the Hamiltonian for the spinor BEC is Chen  _et al._ (2018a,
b)

$\displaystyle H_{\rm s}=$
$\displaystyle\Bigg{[}-\frac{\hbar^{2}}{2M}\frac{\partial}{r\partial
r}\left(r\frac{\partial}{\partial
r}\right)+\frac{L_{z}^{2}}{2Mr^{2}}-\frac{\hbar^{2}}{2M}\frac{\partial^{2}}{\partial
z^{2}}$ (1) $\displaystyle+V({\bf
r})\Bigg{]}\mathds{I}+{\Omega(r)\cos(\phi)S_{x}-\Omega(r)\sin(\phi)S_{y}+\delta
S_{z}},$

where $\mathds{I}$ is a 3$\times$3 identity matrix, $V({\bf
r})=M\omega_{0}^{2}r^{2}/2+M\omega_{z}^{2}z^{2}/2$ constitutes the external
harmonic potential to trap the atoms of mass $M$,
$L_{z}=-i\hbar\partial/\partial\phi$ is the angular momentum operator,
$\Omega(r)=\Omega_{0}\sqrt{e}(r/r_{0})e^{-r^{2}/2r_{0}^{2}}$ is the Raman-
coupling strength with $\Omega_{0}$ and $r_{0}$ as the Rabi frequency and the
radius of the maximum-intensity (cylindrical) surface, respectively, $\delta$
is the Raman detuning, and $S_{x},S_{y}$ and $S_{z}$ are irreducible
representations of the spin-1 angular momentum operators Chen  _et al._
(2018a, b). Under mean-field approximation, the interacting part of the
Hamiltonian $H_{\rm int}$ is given by Kawaguchi and Ueda (2012)

$H_{\rm int}=\frac{c_{0}}{2}\rho+\frac{c_{1}}{2}{\bf F}.{\bf S}$ (2)

with $c_{0}$ and $c_{1}$ as the mean-field interaction parameters. The total
density of the system is given by $\rho$, ${\bf F}=(F_{x},F_{y},F_{z})$ is the
spin-density vector, and ${\bf S}=(S_{x},S_{y},S_{z})$. Since the SOAM
coupling is restricted to the radial plane, and we consider
$\omega_{z}\gg\omega_{0}$, the dominant dynamics is constrained to the same
plane with frozen axial degrees of freedom. We can then integrate out the $z$
degree of freedom from the condensate wave function and describe the system as
quasi-2D on the radial $r$-$\phi$ plane. Starting from the Hamiltonian
$H=H_{\rm s}+H_{\rm int}$, in polar coordinates, we obtain the following
coupled quasi-2D GPEs in dimensionless form

$\displaystyle i\frac{\partial\psi_{\pm 1}}{\partial t}$ $\displaystyle=$
$\displaystyle\mathcal{H}\psi_{\pm 1}+{c_{1}}(\rho_{0}\pm\rho_{-})\psi_{\pm
1}+{c_{1}}\psi_{\mp 1}^{*}\psi_{0}^{2}$ (3a) $\displaystyle\pm\delta\psi_{\pm
1}+\frac{\Omega(r)}{\sqrt{2}}e^{\pm i\phi}\psi_{0},$ $\displaystyle
i\frac{\partial\psi_{0}}{\partial t}$ $\displaystyle=$
$\displaystyle\mathcal{H}\psi_{0}+{c_{1}}\rho_{+}\psi_{0}+2{c_{1}}\psi_{+1}\psi_{-1}\psi_{0}^{*}$
(3b)
$\displaystyle+\frac{\Omega(r)}{\sqrt{2}}(e^{-i\phi}\psi_{+1}+e^{i\phi}\psi_{-1}),$
where  $\displaystyle\mathcal{H}$
$\displaystyle=-\frac{1}{2}\frac{\partial}{r\partial
r}\left(r\frac{\partial}{\partial r}\right)+\frac{L_{\rm
z}^{2}}{2r^{2}}+\frac{r^{2}}{2}+{c_{0}}\rho,$ $\displaystyle\rho$
$\displaystyle=\sum_{j=\pm
1,0}\rho_{j},~{}\rho_{j}=|\psi_{j}|^{2},~{}\rho_{\pm}=\rho_{+1}\pm\rho_{-1}.$

Under geometric renormalization, in terms of $s$-wave scattering lengths
$a_{0}$ and $a_{2}$ in the total spin $0$ and $2$ channels, respectively,
$c_{0}$ and $c_{1}$ take the form

$c_{0}=\sqrt{8\pi\alpha}\frac{N(a_{0}+2a_{2})}{3{a_{\rm osc}}},\quad
c_{1}=\sqrt{8\pi\alpha}\frac{N(a_{2}-a_{0})}{3{a_{\rm osc}}}$ (5)

denoting the spin-independent and spin-dependent interactions, respectively.
The anisotropy parameter $\alpha=\omega_{z}/\omega_{0}$ is defined to be the
trapping frequency ratio along the axial to the radial direction, and $N$ is
the total number of atoms. The units of length, time, energy, and energy
eigenfunctions are considered to be $a_{\rm osc}=\sqrt{\hbar/(M\omega_{0})}$,
$\omega_{0}^{-1}$, $\hbar\omega_{0}$, and $a_{\rm osc}^{-1}$, respectively,
and $\int r\rho(r)drd\phi=1$.

##  III Ground state quantum phases of SOAM coupled spinor BEC

To understand the intercomponent phase relationship imposed by various
competing terms in the Hamiltonian, we consider a generic circularly symmetric
ansatz, $\psi_{j}=f_{j}(r)e^{i(w_{j}\phi+\beta_{j})}$, for the component
wavefunctions, where $w_{j}$ and $\beta_{j}$ are, respectively, the phase-
winding number and constant phase associated with the radially-symmetric real
function $f_{j}$. The phase-dependent part of the interaction energy is
minimized, provided Isoshima and Machida (2002)

$\displaystyle w_{+1}-2w_{0}+w_{-1}$ $\displaystyle=0,$ (6a)
$\displaystyle\beta_{+1}-2\beta_{0}+\beta_{-1}$
$\displaystyle=\begin{cases}2n\pi\quad{\rm for}~{}c_{1}<0,\\\
(2n^{\prime}+1)\pi\quad{\rm for}~{}c_{1}>0,\end{cases}$ (6b)

where $n$ and $n^{\prime}$ are integers. Similarly, SOAM-part of the energy is
minimized if

$\displaystyle w_{+1}-w_{0}$ $\displaystyle=1,\quad w_{0}-w_{-1}=1,$ (7a)
$\displaystyle\beta_{+1}-\beta_{0}$
$\displaystyle=(2p+1)\pi,\quad\beta_{0}-\beta_{-1}=(2p^{\prime}+1)\pi,$ (7b)

where $p$ and $p^{\prime}$ are again integers. If the conditions on the
winding numbers in Eq. (7a) are satisfied, the condition in Eq. (6a) is
satisfied too. On the other hand, conditions between the constant phase
factors in Eqs. (6b) and (7b) can be simultaneously satisfied for $c_{1}<0$
only.

To further substantiate the intercomponent phase relationships imposed by
SOAM-coupling, we extract ${\cal S}=S_{x}\cos\phi-S_{y}\sin\phi$ Duan  _et
al._ (2020) from the single-particle Hamiltonian $H_{\rm s}$. In the limit
when $\Omega_{0}$ is large, $c_{1}$ dependent part of the Hamiltonian can be
neglected, and the phase structure of the emergent ground-state solution is
mainly determined by ${\cal S}$ via its minimum energy eigen spinor. The
normalized eigen spinor of ${\cal S}$ with minimum eigen energy $-1$ can be
written as $(e^{i(m+1)\phi},~{}-\sqrt{2}e^{im\phi},~{}e^{i(m-1)\phi})^{T}/2$
with $m$ being any integer. The phase structure of this eigenspinor is
consistent with phase relations in Eqs. (7a)-(7b). With an increase in $m$,
there is an energy cost from the phase-dependent part of the kinetic energy
suggesting that only small values of phase-winding numbers may emerge.
Numerical results confirm this, where we obtain a solution corresponding to
$m=0$ in large $\Omega_{0}$ limit irrespective of the nature of spin-exchange
interactions. The spinor part of the ground state in this limit tends to
approach the aforementioned eigenstate of ${\cal S}$ with $m=0$.

Various numerical techniques have been employed in the literature to study
spinor BECs in quasi-one-dimensional, quasi-two-dimensional, and three-
dimensional settings Roy  _et al._ (2020); Kaur  _et al._ (2021); Ravisankar
_et al._ (2021); Banger  _et al._ (2021, 2022). In practice, we choose the
finite-difference method and choose different initial guess solutions as an
input to Eqs. (3a)-(3b) to arrive at ground-state solutions. As an example, we
take initial states $\Psi\sim
e^{-r^{2}/2}\times(e^{i(m+1)\phi},~{}-\sqrt{2}e^{im\phi},~{}e^{i(m-1)\phi})^{T}/2$,
with different values of $m$. Besides these initial states, we consider a
random initial guess where $\psi_{j}(r)$ are complex Gaussian random numbers.

At the outset, motivated by the experimental realization of the SOAM-coupled
BECs Chen  _et al._ (2018a, b) using spin-1 87Rb atoms, we validate our
numerical simulations to study and emulate the observed ground-state quantum
phases of the ferromagnetic system in the absence of detuning $\delta=0$.
Similar to the experiment, we consider the 87Rb atoms confined in an
anisotropic harmonic trap with $\omega_{0}=2\pi\times 140$ Hz and $r_{0}=15$
$\mu$m Chen  _et al._ (2018b). However, we take $\omega_{z}=2\pi\times 2400$
Hz enabling us to perform quasi-2D simulations. Here $a_{0}=101.8a_{B}$ and
$a_{2}=101.4a_{B}$ with $a_{B}$ as the Bohr radius Van Kempen  _et al._
(2002). The ground-state densities and phase distributions, obtained
numerically by solving the coupled GPEs (3a)-(3b) with imaginary-time
propagation, for given $\Omega_{0}$ and $N$, are in qualitative agreement with
the experimental results. The ground-state densities calculated for a pair of
$\Omega_{0}$ values with $N=5000$ are shown in Figs. 1(a)-(b).

Figure 1: (Color online) Ground-state densities of the SOAM-coupled 87Rb
spin-1 BEC with $c_{0}=121.28$ and $c_{1}=-0.56$ corresponding to $N=5000$ for
(a) $\Omega_{0}=0.25$ and (b) $\Omega_{0}=1$. The $j=+1,0,$ and $-1$ spin
components carry phase winding numbers of $+2,+1,$ and $0$, respectively, in
(a) and $+1,0,$ and $-1$, respectively, in (b). The radial coordinate and
densities are in the units of $a_{\rm osc}$ and $a_{\rm osc}^{-2}$,
respectively, where $a_{\rm osc}=0.91$ $\mu$m. Figure 2: (Color online) (a)
$\left<\cal S\right>$ as a function of SOAM-coupling strength $\Omega_{0}$ for
87Rb with $c_{0}=121.28$ and $c_{1}=-0.56$ and 23Na with $c_{0}=121.35$ and
$c_{1}=3.80$. Inset in (a): $\langle L_{z}\rangle$ for 87Rb as a function of
$\Omega_{0}$. (b) $|f_{z}|$ and $f$ for 87Rb and 23Na as a function of SO
coupling strength $\Omega_{0}$. The $c_{0}$ and $c_{1}$ for 87Rb and 23Na are
the same as those in (a). Figure 3: (Color online) The ground-state phase
diagrams in (a) $c_{1}/c_{0}$-$\Omega_{0}$ and (b) $N$-$\Omega_{0}$ planes. In
(a) $c_{0}$ was kept fixed at $121.28$ while varying $c_{1}$. In (b)
$c_{1}/c_{0}=-0.0046$ corresponding to 87Rb.

For $\Omega_{0}=0.25$, the solutions with $(+2,+1,0)$ and $(0,-1,-2)$ phase-
winding numbers are two degenerate ground states, and with $\Omega_{0}=1$,
$(+1,0,-1)$ state is obtained as the ground state solution. As we vary
$\Omega_{0}$ from 0 to 20, at small $\Omega_{0}$, due to the co-action of
spin-dependent interaction term and SOAM coupling, (+2,+1,0)-type solution
appears as the ground state. After a critical value of coupling strength (say
$\Omega_{0}^{c}$), $\Omega{\cal S}$ primarily dictates the nature of the
solution to result in $(+1,0,-1)$-type phase. The condition $\langle{\cal
S}\rangle\approx-1$ is satisfied in this latter phase for sufficiently large
$\Omega_{0}$ as shown in Fig. 2(a), which indicates that no further phase can
be expected with higher $\Omega_{0}$. We term these two phases I and II. In
contrast to 87Rb, $(+1,0,-1)$-type is the single ground state phase for 23Na
with $c_{1}>0$. In this case too, $\langle{\cal S}\rangle\approx-1$ at large
$\Omega_{0}$ as shown in Fig. 2(a).

Longitudinal magnetization per particle $f_{z}=\int{F_{z}}d{\bf r}$, spin
expectation per particle $f=\int{|{\bf F}|}d{\bf r}$ where $|{\bf
F}|={\sqrt{F_{x}^{2}+F_{y}^{2}+F_{z}^{2}}}$, and angular momentum per particle
$\langle L_{z}\rangle$ can be used to characterize these ground-state phases.
In the ferromagnetic domain with $c_{0}=121.28$ and $c_{1}=-0.56$, for
$\Omega_{0}<=\Omega_{0}^{c}=0.3$, i.e. in phase I, $\langle L_{z}\rangle\neq
0$ and increases continuously as shown in the inset of Fig. 2(a), whereas
$|f_{z}|\approx 1$ and $f=1$ as shown in Fig. 2(b). For
$\Omega_{0}>\Omega_{0}^{c}$, the transition to phase II is accompanied by
discontinuities in $\langle L_{z}\rangle,|f_{z}|$, and $f$, whereas the former
two reduce to zero, the latter becomes less than one. In the antiferromagnetic
domain, e.g. with $c_{0}=121.35$ and $c_{1}=3.8$, there is no phase transition
with an increase in $\Omega_{0}$ resulting in smooth behavior of the same
quantities. Here $f$ asymptotically approaches one, whereas $|f_{z}|$ and
$\langle L_{z}\rangle$, expectantly, remain zero.

Figure 4: (Color online) (a) and (b) show the spin-texture for 5000 atoms of
87Rb system with the coupling strength $\Omega_{0}=0.25$ and $\Omega_{0}=1$,
respectively.

Furthermore, we calculate the ground state phase diagrams in
$c_{1}/c_{0}$-$\Omega_{0}$ plane, where we fix $c_{0}=121.28$ and vary
$c_{1}$, and $N$-$\Omega_{0}$ plane for fixed $c_{1}/c_{0}=-0.0046$ which
corresponds to 87Rb. The ratio $c_{1}/c_{0}$ may be manipulated experimentally
by tuning one of the scattering lengths by optical Feshbach resonance Theis
_et al._ (2004). These two are respectively shown in Figs. 3(a)-(b), thus
again illustrating that an antiferromagnetic BEC has one ground-state phase in
contrast to the ferromagnetic one. It can be seen that with a decrease in
$c_{1}$ (keeping $c_{0}$ fixed) in the ferromagnetic phase, the domain of
phase I increases, whereas with an increase in the number of atoms (keeping
$c_{1}/c_{0}$ fixed), it decreases. Phase I and II also have distinctive
topological spin textures  ${\bf F}=(F_{x},F_{y},F_{z})$ . For the solutions
in Figs. 1(a) and (b) spin-textures are shown in Figs. 4(a) and (b),
respectively. This allows the identification of phase I and II with the
coreless vortex and polar-core vortex states, respectively Chen  _et al._
(2018b).

##  IV Collective Excitation Spectrum

To study the excitation spectrum, we exploit the innate circular symmetry of
the Hamiltonian. To this end, we perform a local spin rotation about $\hat{z}$
by the azimuthal angle -$\phi$ to remove the $\phi$-dependence from the
Hamiltonian. As a result, the order parameter
$\Psi=(\psi_{+1},\psi_{0},\psi_{-1})^{T}$ is transformed to
$e^{-iS_{z}\phi}{\Psi}=(e^{-i\phi}\psi_{+1},\psi_{0},e^{i\phi}\psi_{-1})^{T}$,
and the transformed Hamiltonian takes form

$\displaystyle{H}=$ $\displaystyle\left[-\frac{1}{2}\frac{\partial}{r\partial
r}\left(r\frac{\partial}{\partial
r}\right)+\frac{(L_{z}+S_{z})^{2}}{2r^{2}}+V(r)\right]\mathds{I}$
$\displaystyle+\Omega(r)S_{x}+H_{\rm int},$ (8)

where ${H}_{\rm int}=c_{0}\rho/2+c_{1}{\bf{F}}.{\bf S}/2$. The Hamiltonian in
Eq. (IV) is circularly symmetric, and one can seek the simultaneous eigen
functions of $H$ and $L_{z}$ with fixed angular momentum $l_{z}=0,1,\ldots$.
For example, the solutions presented in Figs. 1(a) and (b) can now be seen as
corresponding to $l_{z}=1$ and $0$, respectively. We use the Bogoliubov
approach to study the excitation spectrum.  In which we consider the
fluctuations to the ground state by writing the perturbed order parameter as

${\Psi}(r,\phi,t)=e^{-i\mu t+i(l_{z}+S_{z})\phi}[\Psi_{\rm
eq}(r)+\delta{\Psi}(r,t)e^{il_{q}\phi}],$ (9)

where $\Psi_{\rm eq}(r)=(R_{+1}(r),R_{0}(r),R_{-1}(r))^{T}$ is the radial part
of the order parameter with $R_{j}$ as the radial wavefunction corresponding
to the $j^{\rm th}$ spin component, $\mu$ is the chemical potential, and
$l_{q}=0,\pm 1,\pm 2,\ldots$ is the magnetic quantum number associated with
the angular momentum of the quasiparticle excitations. The details of the
Bogoliubov-de Gennes (BdG) analysis are given in the Appendix.

###  IV.1 Non-interacting system

To understand the effect of coupling strength, we first study the single-
particle excitation spectrum. The ground-state solution has phase-winding
numbers $(\pm 1,0)$ in $j=\pm 1,0$ spin states, respectively.

Figure 5: (Color online) Single particle excitation spectrum for spin-1 BEC as
a function of SOAM coupling strength $\Omega_{0}$.

The excitation spectrum is shown in Fig. 5. For $\Omega_{0}=0$, the $n^{\rm
th}$ energy level is $3(n+1)$-fold degenerate, as the single-particle
Hamiltonian is identical with a system of three decoupled isotropic two-
dimensional harmonic oscillators. For example, excitations with energies $0$
and $1$ are three- and six-fold degenerate, respectively. The SOAM-coupling
lifts the degeneracies partially. For example, for $\Omega_{0}\neq 0$, there
is only one zero energy excitation; similarly, the red lines in the spectrum
in Fig. 5 correspond to non-degenerate excitations, whereas the black ones to
two-fold degenerate modes. The non-degenerate modes have the magnetic quantum
number of the excitation $l_{q}=0$, whereas modes with two-fold degeneracy
have $l_{q}\neq 0$.

###  IV.2 Interacting spin-1 BEC

Here we study the excitation spectrum (a) as a function of $\Omega_{0}$ for
fixed $c_{0}$ and $c_{1}$ and (b) as a function of $N$ for fixed $\Omega_{0}$
and $c_{0}/c_{1}$ ratio. Both $\Omega_{0}$ and $N$ can be varied in an
experiment Chen  _et al._ (2018a, b). As was discussed in Sec. III, for
$c_{1}<0$, both phase I and II can appear as the ground-state phases with a
variation of either $\Omega_{0}$ or $N$. We primarily consider 87Rb BEC in the
following discussion.

Figure 6: (Color online) Low-lying excitation spectrum of 87Rb SOAM-coupled
spin-1 BEC with $c_{0}=121.18$ and $c_{1}=-0.56$ as a function of coupling
strength $\Omega_{0}$ of phase I with $l_{q}=0,\pm 1$ in (a), $l_{q}=\pm 2,\pm
3,\pm 4,\ldots$ in (b), and of phase II in (c). Among the named modes,
$l_{q}=0$ for density- and spin-breathing, $l_{q}=+1$ for density-dipole,
$l_{q}=-1$ for spin-dipole, $l_{q}=+2$ for density-quadrupole, and $l_{q}=-2$
for spin-quadrupole modes. In (a) and (c), the dashed (magenta-colored) line
is the variational estimate for the density-breathing mode.

Phase I: Here we consider $c_{0}=121.18$ and $c_{1}=-0.56$ and vary
$\Omega_{0}$. The excitation spectrum for phase I is shown in Fig. 6(a) for
$l_{q}=0,\pm 1$ and (b) for $|l_{q}|\geq 2$. The modes with frequencies $1$
and $2$ are, respectively, dipole and density-breathing modes in Fig. 6(a) .
This identification of a mode is based on the real-time evolution of the
expectation of a suitably chosen observable, as will be discussed in the next
subsection. The presence of ferromagnetic interactions further aids the
lifting of the degeneracy, in this case between the modes with magnetic
quantum numbers $\pm l_{q}$, which are degenerate at the single-particle
level. We have confirmed this, for example, by examining the excitation
spectrum of a system with $c_{0}=121.18$ and $c_{1}=-0.6c_{0}\ll-0.56$ (not
shown here), where the non-degenerate nature of the spectrum is clearly seen.
In phase I, there are two zero-energy Goldstone modes corresponding to two
broken continuous symmetries, namely gauge and rotational symmetry. The latter
corresponds to the symmetry transformation generated by $L_{z}$.

Phase II: As already mentioned in Sec. III, the transition from phase I to II
occurs at $\Omega_{0}>0.3$ for $c_{0}=121.18$ and $c_{1}=-0.56$. The
transition is accompanied by the discontinuities in the excitation spectrum.
The excitation spectrum for phase II is shown in Fig. 6(c). Here among the
low-lying modes are dipole and breathing modes corresponding to both density
and spin channels. Both density- and spin-dipole modes are doubly degenerate
corresponding to magnetic quantum number $l_{q}=\pm 1$. On the other hand,
both density- and spin-breathing modes are non-degenerate with $l_{q}=0$. At
small values of $\Omega_{0}$, the energies of the spin modes are less than
their density-mode analogues. There is a single zero-energy mode due to the
broken gauge symmetry in this phase. Besides these modes, the density- and
spin-quadrupole modes are also marked in the excitation spectrum in Figs.
6(a)-(c).

Additionally, the variation in SOAM-coupling strength leads to avoided
crossings between the pairs of excitations, a few of which are identified by
the blue circles in Fig. 6(c). We observe that the avoided crossing occur
between the density and spin oscillations associated with the same magnetic
quantum number $l_{q}$. In the vicinity of the avoided crossing, the roles of
the density and spin modes are interchanged as shown in Fig. 6(c). We study
this mode-mixing by examining the density $(\delta\rho)$ and spin fluctuations
$(\delta F_{x},\delta F_{y},\delta F_{z})$ yielded by the perturbed order
parameter and defined as

$\displaystyle\delta\rho=$ $\displaystyle 2{\rm
Re}\sum_{j}\psi_{j}\delta\psi_{j}^{*},$ (10a) $\displaystyle\delta F_{x}=$
$\displaystyle\sqrt{2}{\rm
Re}(\psi_{+1}\delta\psi_{0}^{*}+\psi_{0}\delta\psi_{+1}^{*}+\psi_{-1}\delta\psi_{0}^{*}+$
$\displaystyle\psi_{0}\delta\psi_{-1}^{*}),$ (10b) $\displaystyle\delta
F_{y}=$ $\displaystyle-\sqrt{2}{\rm
Im}(-\psi_{+1}\delta\psi_{0}^{*}+\psi_{0}\delta\psi_{+1}^{*}+\psi_{-1}\delta\psi_{0}^{*}-$
$\displaystyle\psi_{0}\delta\psi_{-1}^{*}),$ (10c) $\displaystyle\delta
F_{z}=$ $\displaystyle 2{\rm
Re}(\psi_{+1}\delta\psi_{+1}^{*}-\psi_{-1}\delta\psi_{-1}^{*}).$ (10d)

The order-parameter fluctuation $\delta\Psi(r,\phi,t)$, and hence density and
spin fluctuations, can be constructed with the Bogoliubov quasiparticle
amplitudes $u$ and $v$ corresponding to the frequency $\omega$ of the mode  as
$\delta\psi_{j}(r,\phi,t)\propto
e^{i(l_{z}+S_{z}+l_{q})\phi}\left[u_{j}(r)e^{-i\omega
t}-v_{j}^{*}(r)e^{i\omega t}\right]$. In the excitation spectrum in Fig. 6(c)
at $\Omega_{0}=1$, the density- and spin-dipole modes’ frequencies are
$\omega_{\rm D}=1$ and $\omega_{\rm SD}=0.08$, respectively, and the density-
and spin-breathing modes’ frequencies are $\omega_{\rm B}=1.97$ and
$\omega_{\rm SB}=0.37$, respectively. One can see that the density-dipole,
density-breathing, and spin-dipole modes encounter avoided crossings, whereas
the spin-breathing mode does not. This observation agrees with the density and
spin-density fluctuations evaluated along the $\phi=0$ line and shown in Figs.
7(A)-(F). For the density-dipole mode with $\omega_{\rm D}=1$, both density
and spin channels are excited as is seen from $\delta\rho(r,\phi=0,t)$ and
$\delta F_{\nu}(r,\phi=0,t)$ in Fig. 7(A), where $\nu=x,y,z$. Similarly,
number density, longitudinal, and transverse magnetization densities oscillate
in time, corresponding to the spin-dipole mode in Fig. 7(B), and density-
breathing mode ends up exciting both the number and transverse magnetization
densities in Fig. 7(C). On the other hand, the spin-breathing mode excites the
spin channel alone in Fig. 7(D). The density- and spin-quadrupole modes too
excite both the density and spin fluctuations which are not shown. The
nomenclature of the modes in Figs. 6(a)-(c) is consistent with the density,
$\delta\rho(x,y,t)$, and longitudinal magnetization density, $\delta
F_{z}(x,y,t)$, fluctuations corresponding to density, breathing, and
quadrupole modes in Fig. 8 shown at $t=0,T/4,T/2,3T/4$, and $T$ instants,
where $T$ is the period of the collective excitation.

Figure 7: (Color online) (A) shows the density fluctuations,
$\delta\rho(r,\phi=0,t)$, and spin-density fluctuations, $\delta
F_{\nu}(r,\phi=0,t)$, with $\nu=x,y,z$ corresponding to $\omega_{\rm D}=1$.
(B)-(F) present the same for $\omega_{\rm SD}=0.08$, $\omega_{\rm B}=1.97$,
$\omega_{\rm SB}=0.37$, respectively. The radial and time extents in each
subfigure are $4a_{\rm osc}$ and $5T$, respectively, where $T=2\pi/\omega$ is
the time period of the corresponding mode with $\omega$ frequency. Figure 8:
(Color online) Density and longitudinal magnetization density fluctuations at
$t=0,T/4,T/2,3T/4,$ and $T$ where $T=2\pi/\omega$ with $\omega$ as the mode-
frequency: (A) $\delta\rho(x,y,t)$ for the density-dipole mode with
$\omega_{\rm D}=1$, (B) $\delta F_{z}(x,y,t)$ for the spin-dipole mode with
$\omega_{\rm SD}=0.08$, (C) $\delta\rho(x,y,t)$ for the density-breathing mode
with $\omega_{\rm B}=1.97$, (D) $\delta F_{z}(x,y,t)$ for the spin-breathing
modes with $\omega_{\rm SB}=0.37$, (E) $\delta\rho(x,y,t)$ for the density-
quadrupole mode with $\omega_{\rm Q}=1.46$, and (F) $\delta F_{z}(x,y,t)$ for
the spin-quadrupole mode with $\omega_{\rm SQ}=0.46$. The box size in each
subfigure is $6.4a_{\rm osc}\times 6.4a_{\rm aosc}$. Figure 9: (Color online)
Low-lying excitation spectrum for 87Rb spin-1 BEC with $c_{1}/c_{0}=-0.0046$
as a function of the number of atoms $N$: (a)-(b) for $\Omega_{0}=0.3$ with a
phase transition from phase I to II at $N$ = 5700 and (c) $\Omega_{0}=3$. (a)
corresponds to the spectrum of phase I, whereas (b) and (c) correspond to the
spectrum of phase II. The different colors in (a) signify non-degenerate modes
with different $l_{q}$, while in (b) and (c), red, black, green, blue, and
brown colors correspond, respectively, to the modes with $l_{q}=0,\pm 1,\pm
2$, $\pm 3$ and $\pm 4$.

Next, we study the excitation spectrum as a function of $N$ for
$c_{1}/c_{0}=-0.0046$. Here first, we fix $\Omega_{0}$ to $0.3$, where a phase
transition from phase I to II occurs at $N=5700$. The excitation spectrum, in
this case, for phase I and II are shown in Fig. 9(a) and (b). The same for
$\Omega_{0}=3$ is shown in Fig. 9(c), where phase II is the ground state phase
with no phase transition. The modes in phase II are, again, either non-
degenerate or with two-fold degeneracy. For SOAM-coupled 23Na BEC with
$c_{0}=121.35$ and $c_{1}=3.8$ the excitation spectrum, which is not shown
here, is similar to the spectrum in Fig. 9(c) with some quantitative
differences attributable to different $c_{1}$ values.

###  IV.3 Dynamics

We examine the nature of low-lying collective excitations through the time
evolution of the expectation of physical observables, which also serves to
validate our calculation of the excitation spectrum from BdG equations. Here,
we consider the Hamiltonian with an appropriately chosen time-independent
perturbation, say $H_{\rm s}^{\prime}$ added to its single-particle part
$H_{\rm s}$. This modifies the coupled GP Eqs. (3a)-(3b) with an added term
corresponding to $H_{\rm s}^{\prime}\Psi(r,\phi,t)$ in each equation. We then
solve these resultant GPEs over a finite period of time by considering
previously obtained ground-state solutions as the initial solutions at $t=0$.
Numerically, one needs to consider a two-dimensional spatial grid over here,
for which we choose the Cartesian $x$-$y$ grid.

Figure 10: (Color online) (a) shows the center of mass oscillations, i.e.
$x_{\rm cm}(t)$ as a function of time and (b) corresponding Fourier transform
with a dominant peak at $\omega=1$ for 87Rb spin-1 BEC with $c_{0}=121.28$,
$c_{1}=-0.56$, and $\Omega_{0}=1$. (c) shows the oscillations in the mean
square size of the system $r^{2}_{\rm ms}(t)$ and (d) the corresponding
Fourier transform with a dominant peak at $\omega=1.99$ for the same
interaction and coupling strengths.  Figure 11: (Color online) (a) shows
$d_{x}(t)$ as a function of time and (b) corresponding Fourier transform with
a dominant peak at $\omega=0.1$ for 87Rb spin-1 BEC with $c_{0}=121.28$,
$c_{1}=-0.56$, and $\Omega_{0}=1$. Similarly, (c) and (d) show the
$d^{2}_{x}(t)$ and its Fourier transform with a dominant peak at $\omega=0.37$
for the same interaction and coupling strengths.

We consider $c_{0}=121.28$, $c_{1}=-0.56$, and $\Omega_{0}=1$, which yielded
the ground-state phase in Fig. 1(b), as an example set of parameters to study
the dynamics. To excite the density-dipole mode, we take the perturbation
$H_{\rm s}^{\prime}=\lambda x$, where $\lambda\ll 1$. We then examine the
dynamics of the center of mass of the BEC via $x_{\rm cm}(t)=\langle
x\rangle=\sum_{j=\pm 1,0}\int x\rho_{j}(x,y,t)dxdy$ which is plotted in Fig.
10(a). We also compute its Fourier transform $\widehat{x}_{\rm cm}(\omega)$ to
demonstrate that the dominant frequency resonates at $\omega=1$ as can be seen
in Fig. 10(b) and matches with $\omega_{\rm D}=1$ in the BdG spectrum in Fig.
6(c). We could have chosen $H_{s}^{\prime}=\lambda y$ and then calculated
$y_{\rm cm}(t)$ giving us the same excitation frequency. This is a consequence
of the two-fold degeneracy in the density-dipole mode. We have checked that
this mode can also be excited by shifting the minima of the external trapping
potential. Similarly, to examine the excitation of the density-breathing mode
with $H_{\rm s}^{\prime}=\lambda(x^{2}+y^{2})$, where the relevant observable
is $r^{2}=x^{2}+y^{2}$, we calculate mean square radius $r_{\rm
ms}^{2}(t)=\langle r^{2}\rangle$ as a function of time, which is plotted in
Fig. 10(c). The Fourier transform $\widehat{r_{\rm ms}^{2}}(\omega)$ of
$r_{\rm ms}^{2}(t)$ reveals a dominant peak at $\omega=1.99$ which is close to
BdG result of $\Omega_{\rm B}=1.97$. This mode, again, can be excited by
perturbing the trap strength. Similarly, the spin-dipole mode can be excited
by adding a perturbation $H^{\prime}=\lambda xS_{z}$ or $\lambda yS_{z}$ with
$xS_{z}$ or $yS_{z}$ as the pertinent observable corresponding to the spin-
dipole mode. The two possible observables again reflect the two-fold
degeneracy of spin-dipole modes. The time-variation of $d_{x}(t)=\langle
xS_{z}\rangle=\sum_{j=+1,-1}\int x\rho_{j}(x,y,t)dxdy$ is shown in Fig. 11(a)
and its Fourier transform in Fig. 11(b) has a dominant peak at $\omega=0.1$,
which corresponds to the spin-dipole mode labeled in Fig. 6(c) with
$\omega_{\rm SD}=0.08$. Similarly, the spin-breathing mode corresponds to
observable $r^{2}S_{z}$. In Figs. 11(c) and (d), we show the dynamics of
$d_{r}^{2}(t)=\langle r^{2}S_{z}\rangle$, i.e. the relative difference in the
mean-square radii of the $j=\pm 1$ components and the associated Fourier
transform, respectively, with a dominant peak at $\omega=0.37$, in agreement
with $\omega_{\rm SB}$ in Fig. 6(c). Finally, the density- and spin-quadrupole
modes’ frequencies calculated from the time evolution of $\langle xy\rangle$
and $\langle xyS_{z}\rangle$ are in agreement with the numbers in Fig. 6(c).

###  IV.4 Variational analysis

For an SOAM-coupled spin-1 system, a few low-lying modes can be studied using
a time-dependent variational method Pérez-García  _et al._ (1996). For
example, to calculate the density-breathing mode, we consider the following
variational ansatz

$\Psi=\frac{r}{2\sqrt{\pi}\sigma(t)^{2}}\exp\left[-\frac{r^{2}}{2\sigma(t)^{2}}+i\alpha(t)r^{2}\right]\times\begin{pmatrix}e^{i(m+1)\phi}\\\
-\sqrt{2}e^{im\phi}\\\ e^{i(m-1)\phi},\end{pmatrix}$ (11)

where $\sigma(t)$ and $\alpha(t)$ are time-dependent variational parameters
used to denote the width of condensate and chirp of Gaussian pulse,
respectively, and $m=\pm 1$ for phase I or $0$ for phase II. The Lagrangian of
the system is given by

$L=\sum_{j}\int
drd\phi\frac{i}{2}\left(\psi_{j}^{*}\frac{\partial\psi_{j}}{\partial
t}-\psi_{j}\frac{\partial\psi_{j}^{*}}{\partial t}\right)-E,$ (12)

where energy $E$ is defined as

$\displaystyle E=$
$\displaystyle\int_{0}^{\infty}\int_{0}^{2\pi}\left[\sum_{j}\psi_{j}^{*}\left(-\frac{1}{2}\frac{\partial}{r\partial
r}\left(r\frac{\partial}{\partial r}\right)+\frac{L_{\rm
z}^{2}}{2r^{2}}+\frac{r^{2}}{2}\right)\psi_{j}\right.$
$\displaystyle\left.+\frac{c_{0}}{2}\rho^{2}+\frac{c_{1}}{2}(\rho_{1}+\rho_{0}-\rho_{-1})\rho_{-1}+\frac{c_{1}}{2}(\rho_{1}+\rho_{-1})\rho_{0}+\right.$
$\displaystyle\left.c_{1}(\psi_{-1}^{*}\psi_{0}^{2}\psi_{1}^{*}+\psi_{-1}{\psi_{0}^{2}}^{*}\psi_{1})+\frac{\Omega(r)}{\sqrt{2}}(\psi_{1}^{*}e^{i\phi}\psi_{0}+\right.$
$\displaystyle\left.\psi_{1}e^{-i\phi}\psi_{0}^{*}+\psi_{0}^{*}e^{i\phi}\psi_{-1}+\psi_{-1}^{*}e^{i\phi}\psi_{0})\right]drd\phi.$
(13)

For $m=\pm 1$, the (coupled) Euler-Lagrange equations are

$\displaystyle\ddot{\sigma(t)}=$
$\displaystyle\frac{\sigma}{2}\left(\frac{6\sqrt{2\pi}\sqrt{e}\Omega_{0}\sqrt{\frac{1}{r_{0}^{2}}+\frac{2}{\sigma^{2}}}\left(r_{0}^{7}-2r_{0}^{5}\sigma^{2}\right)}{\left(2r_{0}^{2}+\sigma^{2}\right)^{4}}-2\right)$
$\displaystyle+\frac{c_{0}+c_{1}+10\pi}{8\pi\sigma^{3}},$ (14a)
$\displaystyle\alpha=$ $\displaystyle\frac{\dot{\sigma}}{2\sigma},$ (14b)

where $\dot{}$ denotes the time derivative. The equilibrium width $\sigma_{0}$
of the condensate satisfies

$\frac{c_{0}+c_{1}+10\pi}{4\pi\sigma^{4}_{0}}+\frac{6\sqrt{2\pi}\sqrt{e}\Omega_{0}\sqrt{\frac{1}{r_{0}^{2}}+\frac{2}{\sigma^{2}_{0}}}\left(r_{0}^{7}-2r_{0}^{5}\sigma^{2}_{0}\right)}{\left(2r_{0}^{2}+\sigma^{2}_{0}\right)^{4}}=2.$

The frequency of the oscillation in width calculated by linearizing Eq. (14a)
about equilibrium width $\sigma_{0}$ is

$\displaystyle\omega_{\rm B}^{\rm I}=$
$\displaystyle\left[\frac{15\sqrt{2\pi}r_{0}^{4}\sqrt{e}\sigma_{0}\Omega_{0}(3r_{0}^{2}-2\sigma^{2}_{0})\sqrt{{2r_{0}^{2}}+{\sigma^{2}_{0}}}+1}{(2r_{0}^{2}+\sigma^{2}_{0})^{5}}+\right.$
$\displaystyle\left.\frac{3(c_{0}+c_{1}+10\pi)}{8\pi\sigma^{4}_{0}}\right]^{1/2}.$
(15)

Similarly, for $m=0$ in Eq. (11), the density breathing mode is

$\displaystyle\omega_{\rm B}^{\rm II}=$
$\displaystyle\left[\frac{15\sqrt{2\pi}r_{0}^{4}\sqrt{e}\sigma_{0}\Omega_{0}(3r_{0}^{2}-2\sigma^{2}_{0})\sqrt{2{r_{0}^{2}}+{\sigma^{2}_{0}}}+1}{(2r_{0}^{2}+\sigma^{2}_{0})^{5}}+\right.$
$\displaystyle\left.\frac{3(c_{0}+c_{1}+6\pi)}{8\pi\sigma^{4}_{0}}\right]^{1/2}.$
(16)

The variationally calculated density-breathing mode’s frequency agrees with
the values in the BdG spectrum as demonstrated in Figs. 6(a) and (c) for
phases I and II, respectively.

##  V Summary and conclusions

We have investigated the low-lying collective excitations of the coreless
vortex and the polar-core vortex phases supported by the spin-1 BECs with SOAM
coupling. The existence of the two phases is seen in the full phase diagrams
in the ratio of interaction strengths versus coupling strength and also the
number of atoms versus coupling strength planes. We have studied the
excitation spectrum as a function of two experimentally controllable
parameters, namely coupling strength and the number of atoms. The excitation
spectrums are characterized by the discontinuities across the phase boundary
between the two phases and within a phase by avoided crossings between the
modes with the same magnetic quantum number of excitations. The avoided
crossings signal the hybridization of the density and spin channels; the
nature of spin and density fluctuations has indeed confirmed this. Among the
low-lying modes, we identify dipole, breathing, and quadrupole modes for
density and spin channels. The frequencies of these named modes are further
validated from the time evolution of the expectations of the physical
observables when an apt time-independent perturbation is added to the system’s
Hamiltonian. An analytic estimate for the density-breathing modes has also
been obtained using the variational analysis.

## Acknowledgments

AR acknowledges the support of the Science and Engineering Research Board
(SERB), Department of Science and Technology, Government of India under the
project SRG/2022/000057 and IIT Mandi seed-grant funds under the project
IITM/SG/AR/87. AR acknowledges National Supercomputing Mission (NSM) for
providing computing resources of PARAM Himalaya at IIT Mandi, which is
implemented by C-DAC and supported by the Ministry of Electronics and
Information Technology (MeitY) and Department of Science and Technology (DST),
Government of India. S.G. acknowledges support from the Science and
Engineering Research Board, Department of Science and Technology, Government
of India through Project No. CRG/2021/002597.

## Appendix

Bogoliubov-de Gennes (BdG) analysis: The fluctuation $\delta{\Psi(r,t)}$ to
the equilibrium order parameter in Eqn. (9) is
$\delta{\Psi}(r,t)=u(r)e^{-i\omega t}-v^{*}(r)e^{i\omega t}$, where $u(r)$ and
$v(r)$ are Bogoliubov amplitudes and $\omega$ is the excitation frequency.
Linearization of the three coupled Gross-Pitaevskii Eqs. (3a)-(3b) and the
conjugate set of equations using perturbed order parameter in Eq. (9) yields
following six-coupled BdG equations:

$\displaystyle\omega u_{+1}$ $\displaystyle=$
$\displaystyle\left[-\frac{\nabla^{2}_{r}}{2}+\frac{r^{2}}{2}+\delta+\frac{(l_{q}+l_{z}+1)^{2}}{2}-\mu+c_{0}R_{+1}^{2}+c_{1}(2R_{+1}^{2}+R_{0}^{2}-R_{-1}^{2})\right]u_{+1}$
(17a)
$\displaystyle+\left[\frac{\Omega(r)}{\sqrt{2}}+R_{+1}R_{0}(c_{0}+c_{1})+2c_{1}R_{0}R_{-1}\right]u_{0}+R_{+1}^{2}(c_{0}+c_{1})v_{+1}+R_{+1}R_{0}(c_{0}+c_{1})v_{0}$
$\displaystyle+R_{+1}R_{-1}(c_{0}-c_{1})u_{-1}+(R_{+1}R_{-1}(c_{0}-c_{1})+2c_{1}R_{0}^{2})v_{-1},$
$\displaystyle-\omega v_{+1}$ $\displaystyle=$
$\displaystyle\left[-\frac{\nabla^{2}_{r}}{2}+\frac{r^{2}}{2}+\delta+\frac{(l_{q}+l_{z}+1)^{2}}{2}-\mu+c_{0}R_{+1}^{2}-c_{1}(2R_{+1}^{2}+R_{0}^{2}-R_{-1}^{2})\right]v_{+1}$
(17b)
$\displaystyle+\left[\frac{\Omega(r)}{\sqrt{2}}+R_{+1}R_{0}(c_{0}+c_{1})+2c_{1}R_{0}R_{-1}\right]v_{0}+R_{+1}^{2}(c_{0}+c_{1})u_{+1}+R_{+1}R_{0}(c_{0}+c_{1})u_{0}$
$\displaystyle+R_{+1}R_{-1}(c_{0}-c_{1})v_{-1}+\big{[}R_{+1}R_{-1}(c_{0}-c_{1})+2c_{1}R_{0}^{2}\big{]}u_{-1},$
$\displaystyle\omega u_{0}$ $\displaystyle=$
$\displaystyle\left[-\frac{\nabla^{2}_{r}}{2}+\frac{r^{2}}{2}+\frac{(l_{q}+l_{z})^{2}}{2}-\mu+c_{0}R_{0}^{2}-c_{1}(R_{+1}^{2}+R_{-1}^{2})\right]u_{0}+\left[\frac{\Omega(r)}{\sqrt{2}}-R_{+1}^{2}(c_{0}+c_{1})\right]u_{+1}$
(17c)
$\displaystyle+R_{+1}R_{0}(c_{0}+c_{1})v_{+1}+(c_{0}R_{0}^{2}+2c_{1}R_{+1}R_{-1})v_{0}+\left[\frac{\Omega(r)}{\sqrt{2}}+R_{0}R_{-1}(c_{0}+c_{1})-2c_{2}R_{+1}R_{0}\right]u_{-1}$
$\displaystyle+R_{+1}R_{-1}(c_{0}+c_{1})v_{-1}$ $\displaystyle-\omega v_{0}$
$\displaystyle=$
$\displaystyle\left[-\frac{\nabla^{2}_{r}}{2}+\frac{r^{2}}{2}+\frac{(l_{q}+l_{z})^{2}}{2}-\mu+c_{0}R_{0}^{2}-c_{1}(R_{+1}^{2}+R_{-1}^{2})\right]v_{0}+\left[\frac{\Omega(r)}{\sqrt{2}}-R_{+1}^{2}(c_{0}+c_{1})\right]v_{+1}$
(17d)
$\displaystyle+R_{+1}R_{0}(c_{0}+c_{1})u_{+1}+(c_{0}R_{0}^{2}+2c_{1}R_{+1}R_{-1})u_{0}+\left[\frac{\Omega(r)}{\sqrt{2}}+R_{0}R_{-1}(c_{0}+c_{1})-2c_{1}R_{+1}R_{0}\right]v_{-1}$
$\displaystyle+R_{+1}R_{-1}(c_{0}+c_{1})u_{-1}$ $\displaystyle\omega u_{-1}$
$\displaystyle=$
$\displaystyle\left[-\frac{\nabla^{2}_{r}}{2}+\frac{r^{2}}{2}-\delta+\frac{(l_{q}+l_{z}-1)^{2}}{2}-\mu+c_{0}R_{-1}^{2}+c_{1}(2R_{-1}^{2}+R_{0}^{2}-R_{+1}^{2})\right]u_{-1}$
(17e)
$\displaystyle+\left[\frac{\Omega(r)}{\sqrt{2}}+R_{0}R_{-1}(c_{0}+c_{1})+2c_{1}R_{+1}R_{0}\right]u_{0}+(c_{0}-c_{1})R_{+1}R_{-1}u_{+1}$
$\displaystyle+(R_{+1}R_{-1}(c_{0}-c_{1})+2c_{1}R_{0}^{2})v_{+1}+R_{-1}^{2}(c_{0}+c_{1})v_{-1}+R_{+1}R_{0}(c_{0}+c_{1})v_{0}$
$\displaystyle-\omega v_{-1}$ $\displaystyle=$
$\displaystyle\big{[}-\frac{\nabla^{2}_{r}}{2}+\frac{r^{2}}{2}-\delta+\frac{(l_{q}+l_{z}-1)^{2}}{2}-\mu+c_{0}R_{-1}^{2}+c_{1}(2R_{-1}^{2}+R_{0}^{2}-R_{+1}^{2})\big{]}v_{-1}$
(17f)
$\displaystyle+\big{[}\frac{\Omega(r)}{\sqrt{2}}+R_{0}R_{-1}(c_{0}+c_{1})+2c_{1}R_{+1}R_{0}\big{]}v_{0}+(c_{0}-c_{1})R_{+1}R_{-1}v_{+1}$
$\displaystyle+(R_{+1}R_{-1}(c_{0}-c_{1})+2c_{1}R_{0}^{2})u_{+1}+R_{-1}^{2}(c_{0}+c_{1})u_{-1}+R_{+1}R_{0}(c_{0}+c_{1})u_{0}$

where $\nabla^{2}_{r}=-\partial^{2}/(2\partial r^{2})-{\partial}/(2r\partial
r)$, $l_{z}=1$ for phase I and $0$ for phase II. To solve coupled Eqs.
(17a)-(17f), we use finite-difference method to discretize these equations
over the spatial radial grid Gao and Cai (2020), thus transforming the BdG
equations to a matrix eigenvalue equation which can be solved using standard
matrix diagonalization subroutines.

## References

  * Lin  _et al._ (2011) Y.-J. Lin, K. Jiménez-García, and I. B. Spielman, Nature 471, 83 (2011).
  * Campbell  _et al._ (2016) D. Campbell, R. Price, A. Putra, A. Valdés-Curiel, D. Trypogeorgos, and I. Spielman, Nat. Commun. 7, 10897 (2016).
  * Luo  _et al._ (2016) X. Luo, L. Wu, J. Chen, Q. Guan, K. Gao, Z.-F. Xu, L. You, and R. Wang, Scientific reports 6, 18983 (2016).
  * Goldman  _et al._ (2014) N. Goldman, G. Juzeliūnas, P. Öhberg, and I. B. Spielman, Reports on Progress in Physics 77, 126401 (2014).
  * Wang  _et al._ (2010) C. Wang, C. Gao, C.-M. Jian, and H. Zhai, Phys. Rev. Lett. 105, 160403 (2010).
  * Martone  _et al._ (2016) G. I. Martone, F. V. Pepe, P. Facchi, S. Pascazio, and S. Stringari, Phys. Rev. Lett. 117, 125301 (2016).
  * Ramachandhran  _et al._ (2012) B. Ramachandhran, B. Opanchuk, X.-J. Liu, H. Pu, P. D. Drummond, and H. Hu, Phys. Rev. A 85, 023606 (2012).
  * Ruokokoski  _et al._ (2012) E. Ruokokoski, J. A. M. Huhtamäki, and M. Möttönen, Phys. Rev. A 86, 051607 (2012).
  * Pethick and Smith (2008) C. J. Pethick and H. Smith,  _Bose–Einstein Condensation in Dilute Gases_ (Cambridge University Press, 2008).
  * Pitaevskii and Stringari (2016) L. Pitaevskii and S. Stringari,  _Bose-Einstein condensation and superfluidity_ , Vol. 164 (Oxford University Press, 2016).
  * 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).
  * Chen  _et al._ (2017) L. Chen, H. Pu, Z.-Q. Yu, and Y. Zhang, Phys. Rev. A 95, 033616 (2017).
  * Geier  _et al._ (2021) K. T. Geier, G. I. Martone, P. Hauke, and S. Stringari, Phys. Rev. Lett. 127, 115301 (2021).
  * Rajat  _et al._ (2022) Rajat, A. Roy, and S. Gautam, Phys. Rev. A 106, 013304 (2022).
  * Zhang  _et al._ (2019) D. Zhang, T. Gao, P. Zou, L. Kong, R. Li, X. Shen, X.-L. Chen, S.-G. Peng, M. Zhan, H. Pu, and K. Jiang, Phys. Rev. Lett. 122, 110402 (2019).
  * Chen  _et al._ (2018a) H.-R. Chen, K.-Y. Lin, P.-K. Chen, N.-C. Chiu, J.-B. Wang, C.-A. Chen, P. Huang, S.-K. Yip, Y. Kawaguchi, and Y.-J. Lin, Phys. Rev. Lett. 121, 113204 (2018a).
  * Chen  _et al._ (2018b) P.-K. Chen, L.-R. Liu, M.-J. Tsai, N.-C. Chiu, Y. Kawaguchi, S.-K. Yip, M.-S. Chang, and Y.-J. Lin, Phys. Rev. Lett. 121, 250401 (2018b).
  * Hu  _et al._ (2015) Y.-X. Hu, C. Miniatura, and B. Grémaud, Phys. Rev. A 92, 033615 (2015).
  * DeMarco and Pu (2015) M. DeMarco and H. Pu, Phys. Rev. A 91, 033630 (2015).
  * Sun  _et al._ (2015) K. Sun, C. Qu, and C. Zhang, Phys. Rev. A 91, 063627 (2015).
  * Qu  _et al._ (2015) C. Qu, K. Sun, and C. Zhang, Phys. Rev. A 91, 053630 (2015).
  * Chen  _et al._ (2016) L. Chen, H. Pu, and Y. Zhang, Phys. Rev. A 93, 013629 (2016).
  * Peng  _et al._ (2022) S.-G. Peng, K. Jiang, X.-L. Chen, K.-J. Chen, P. Zou, and L. He, AAPPS Bulletin 32, 36 (2022).
  * Duan  _et al._ (2020) Y. Duan, Y. Bidasyuk, and A. Surzhykov, Phys. Rev. A 102, 063328 (2020).
  * Chen  _et al._ (2020a) X.-L. Chen, S.-G. Peng, P. Zou, X.-J. Liu, and H. Hu, Phys. Rev. Res. 2, 033152 (2020a).
  * Chiu  _et al._ (2020) N. Chiu, Y. Kawaguchi, S. Yip, and Y. Lin, New J. Phys. 22, 093017 (2020).
  * Bidasyuk  _et al._ (2022) Y. Bidasyuk, K. Kovtunenko, and O. Prikhodko, Phys. Rev. A 105, 023320 (2022).
  * Vasić and Balaž (2016) I. Vasić and A. Balaž, Phys. Rev. A 94, 033627 (2016).
  * Chen  _et al._ (2020b) K.-J. Chen, F. Wu, J. Hu, and L. He, Phys. Rev. A 102, 013316 (2020b).
  * Kawaguchi and Ueda (2012) Y. Kawaguchi and M. Ueda, Phys. Rep. 520, 253 (2012).
  * Isoshima and Machida (2002) T. Isoshima and K. Machida, Phys. Rev. A 66, 023602 (2002).
  * Roy  _et al._ (2020) A. Roy, S. Pal, S. Gautam, D. Angom, and P. Muruganandam, Comput. Phys. Commun. 256, 107288 (2020).
  * Kaur  _et al._ (2021) P. Kaur, A. Roy, and S. Gautam, Comput. Phys. Commun. 259, 107671 (2021).
  * Ravisankar  _et al._ (2021) R. Ravisankar, D. Vudragović, P. Muruganandam, A. Balaž, and S. K. Adhikari, Comput. Phys. Commun. 259, 107657 (2021).
  * Banger  _et al._ (2021) P. Banger, P. Kaur, and S. Gautam, Int. J. Mod. Phys. C 33, 2250046 (2021).
  * Banger  _et al._ (2022) P. Banger, P. Kaur, A. Roy, and S. Gautam, Comput. Phys. Commun. 279, 108442 (2022).
  * Van Kempen  _et al._ (2002) E. Van Kempen, S. Kokkelmans, D. Heinzen, and B. Verhaar, Phys. Rev. Lett. 88, 093201 (2002).
  * Theis  _et al._ (2004) M. Theis, G. Thalhammer, K. Winkler, M. Hellwig, G. Ruff, R. Grimm, and J. H. Denschlag, Phys. Rev. Lett. 93, 123001 (2004).
  * Pérez-García  _et al._ (1996) V. M. Pérez-García, H. Michinel, J. I. Cirac, M. Lewenstein, and P. Zoller, Phys. Rev. Lett. 77, 5320 (1996).
  * Gao and Cai (2020) Y. Gao and Y. Cai, J. Comput. Phys. 403, 109058 (2020).


# MarkerPose: Robust Real-time Planar Target Tracking for Accurate Stereo Pose
Estimation

Jhacson Meza Lenny A. Romero Andres G. Marrugo
Universidad Tecnológica de Bolívar  Cartagena  Colombia
{jmeza, lromero<EMAIL_ADDRESS>
###### Abstract

Despite the attention marker-less pose estimation has attracted in recent
years, marker-based approaches still provide unbeatable accuracy under
controlled environmental conditions. Thus, they are used in many fields such
as robotics or biomedical applications but are primarily implemented through
classical approaches, which require lots of heuristics and parameter tuning
for reliable performance under different environments. In this work, we
propose MarkerPose, a robust, real-time pose estimation system based on a
planar target of three circles and a stereo vision system. MarkerPose is meant
for high-accuracy pose estimation applications. Our method consists of two
deep neural networks for marker point detection. A SuperPoint-like network for
pixel-level accuracy keypoint localization and classification, and we
introduce EllipSegNet, a lightweight ellipse segmentation network for sub-
pixel-level accuracy keypoint detection. The marker’s pose is estimated
through stereo triangulation. The target point detection is robust to low
lighting and motion blur conditions. We compared MarkerPose with a detection
method based on classical computer vision techniques using a robotic arm for
validation. The results show our method provides better accuracy than the
classical technique. Finally, we demonstrate the suitability of MarkerPose in
a 3D freehand ultrasound system, which is an application where highly accurate
pose estimation is required. Code is available in Python and C++ at
https://github.com/jhacsonmeza/MarkerPose.

##  1 Introduction

Planar markers or fiducial markers are useful in robotics, computer vision,
medical imaging, and related fields for 6-degree-of-freedom (DoF) pose
estimation, mapping, localization, and other pose-related tasks. Tracking an
object’s pose by attaching a marker to it and using a vision system is
attractive because it can be highly accurate while still being a low-cost
alternative when compared to other strategies [3]. For instance, marker-less
techniques are not as robust or fast or not used for precise pose estimation
[16]. Moreover, marker-based methods are useful in different environments.
Robustly detecting these targets under different lighting or motion conditions
is essential for ensuring reliable system operation and accurate tracking for
spatial registration of acquired data [11, 12].

There are several applications where highly accurate pose estimation is
crucial. For example, for robot path planning and navigation [13], Unmanned
Aerial Vehicle (UAV) refueling [22], or for surgical instrument guidance [2].
These applications require a robust marker detection method under different
environments and lighting conditions to register accurate poses. Furthermore,
motion blur robustness is highly desired for real-time operation.

In this paper, we propose MarkerPose: a robust, real-time pose estimation
method using a planar target with concentric circles. Our method is inspired
by Deep ChArUco [9] and consists of two main stages. In the first stage, with
a SuperPoint-like network, the circles’ centers are estimated at a pixel
level. In the second stage, the points are estimated at sub-pixel accuracy
with a proposed lightweight encoder-decoder network called EllipSegNet. The
network segments the marker ellipses for sub-pixel center estimation. The pair
of stereo images are processed with the proposed method, and through
triangulation, the marker’s pose is determined. The experiments show the
highly accurate results of MarkerPose compared with a classical sub-pixel
center detection method [19]. Also, we show the potential and suitability of
our proposed method for pose estimation under low lighting conditions and high
motion blur.

Our main contributions are:

  1. 1.

Real-time, highly accurate pose estimation based on triangulation using a
marker of circles and stereo images.

  2. 2.

EllipSegNet: a lightweight network for ellipse segmentation meant for real-
time applications.

  3. 3.

Sub-pixel marker points detection at full resolution through ellipse fitting,
avoiding detection errors by detecting in a lower resolution and rescaling to
the original.

##  2 Related Work

Deep Convolutional Neural Networks (CNN) have been extensively used for object
detection with, for example, Faster R-CNN [18] or YOLO [17]. Additionally,
there are different examples for object segmentation like Mask-RCNN [8] or
DeepLab [4]. Detection in challenging conditions has also been addressed. For
example, Rozumnyi et al. [20] tackle fast-moving object detection in
challenging high-movement conditions. However, these methods are not readily
applicable for 6DoF object pose estimation. Alternatively, other works have
addressed object detection for marker-less pose estimation. For example,
Nakako [15] proposed a stereo vision for pose estimation using a robot
manipulator. Xiang et al. [24] proposed to jointly segment and estimate
different objects’ poses. Wang et al. [23] addressed pose estimation with
depth information using an RGB-D camera. Nevertheless, these approaches are
limited to moderately accurate pose estimation applications.

Pose estimation can be approached via keypoint detection. There are multiple
examples of keypoint detection with CNNs for marker-less pose estimation.
Human body pose estimation [1], hand pose estimation [6] or head pose
estimation [7] are popular examples in the computer vision community where 2D
points are detected using a single image for marker-less pose estimation.
There are methods for interest point detection using CNN. An example is
SuperPoint [5], where 2D locations and descriptors are estimated for reliable
performance of pipelines like Structure from Motion (SfM) or Simultaneous
Localization and Mapping (SLAM), where stable and repeatable keypoints are
required. By estimating the descriptors with SuperPoint the matching between
2D points is more easily tackled. Recently, SuperGlue [21] has been introduced
for robust matching of interest points in a pair of images, using a Graph
Neural Network. The interest points can be estimated with a CNN like
SuperPoint, or other approaches.

Different planar marker methods for pose tracking provide highly accurate pose
estimation compared to previously discussed approaches. However, the target’s
keypoints detection is typically carried out using classical computer vision
techniques, which require lots of heuristics and parameter tuning for reliable
performance under different environments. Recently, Deep ChArUco [9] was
proposed for pose estimation based on a planar marker that is robust in low
lighting and high movement conditions. It estimates 2D corners in a sub-pixel
accuracy with two CNNs. While being an attractive approach, it is focused on
the pose estimation of a ChArUco board.

##  3 Proposed Approach

Here we describe MarkerPose, our proposed pose estimation approach. For that,
we used a three-coplanar circle marker and a calibrated stereo vision system.
The 3D coordinates of the centers define the 3D frame of the target. The
circles, which due to perspective projection look like ellipses in the image,
are labeled as shown in Fig. 1. These labels help us to solve the
correspondences between both views for triangulation. But for that, we need to
estimate these centers in a sub-pixel-level accuracy and classify them with
the corresponding label or ID.

Figure 1: Fiducial marker used for pose estimation and the respective labels
of each circle.

The proposed MarkerPose method is shown in Fig. 2, and consists of three
stages. For pose estimation, we need first to detect the centers of the
ellipses in a sub-pixel-level accuracy for both views. We perform center
detection in the first two stages of our method: pixel centers detection and
sub-pixel centers estimation. In the first stage, marker points in a pixel-
level accuracy, and their IDs are estimated with a two-headed network based on
SuperPoint [5]. This SuperPoint variant is the same used in Deep ChArUco [9],
where an ID-head replaces the descriptor head, and the same point localization
head of SuperPoint is used. These classified rough 2D locations are used in
the second stage to extract three square patches from each image around the
pixel points that contain each ellipse. Then, we segment the ellipses to use
their contour for sub-pixel centroid estimation. For that, we introduce
EllipSegNet, a lightweight encoder-decoder network for fast ellipse contour
estimation. Finally, in the last stage, with the sub-pixel matches of both
views, triangulation is applied for 3D pose estimation. The correspondences
are established with the points’ IDs. Our method is robust to low lighting
conditions and motion blur, essential for real-time applications and usage in
different environments.

Figure 2: MarkerPose: the proposed approach for marker-based pose estimation,
which consists of three stages. The first and second stage is for sub-pixel
centers estimation, and the last stage for triangulation and target pose
estimation.

###  3.1 Pixel Centers Detection

Point detection from the first and second stages of our proposed method is
shown in more detail in Fig. 3 for a single image. In this section, we focus
on the first stage of MarkerPose which is based on a SuperPoint-like network.
We use $1280\times 1024$ images for training and testing. With the SuperPoint-
like network, 2D pixel-level accuracy center detection is performed at a lower
resolution of $320\times 240$. For that, the input $1280\times 1024$ images
are resized to $320\times 240$. The estimated rough 2D locations are then
scaled-back to full resolution.

The original SuperPoint network [5] consists of a backbone and two heads: a
head for point localization and another for point descriptors. Instead of the
descriptors head, an ID-head is used for point classification as shown Fig. 3.
The backbone is a shared encoder for both branches, which reduces the
dimensionality of an image $W\times H$ to a factor of eight $W/8\times H/8$,
using eight convolution layers and three $2\times 2$ max-pooling layers after
every two convolutions. The point localization head produces a $W/8\times
H/8\times 65$ tensor, which is a pixel location distribution of each $8\times
8$ region of the input image and an additional channel for background or non-
point locations. The ID head produces a $W/8\times H/8\times 4$ tensor, which
is a classification distribution for each $8\times 8$ region of the input
image, where we have three labels for each center (see Fig. 1) and an
additional label for background or non-center pixels. In this way, we have the
following labels: 0 for $c_{0}$, 1 for $c_{1}$, 2 for $c_{2}$, and 3 for
background. With the estimated IDs, we sort the resulting array of 2D points
from 0 to 2 labels, where the first row corresponds to the $c_{0}$ point, and
so on. With the output from this network, we solve the correspondence problem
between both views.

Figure 3: Sub-pixel centers detection for a single image.

The loss for this stage is the sum of the center location loss
$\mathcal{L}_{loc}$, and the ID classification loss $\mathcal{L}_{id}$,

$\mathcal{L}(\hat{\mathbf{P}},\mathbf{P},\hat{\mathbf{C}},\mathbf{C})=\lambda_{1}\mathcal{L}_{loc}(\hat{\mathbf{P}},\mathbf{P})+\lambda_{2}\mathcal{L}_{id}(\hat{\mathbf{C}},\mathbf{C})\enspace,$
(1)

where $\hat{\mathbf{P}}$ is the predicted $W/8\times H/8\times 65$ location
tensor by the network, and $\mathbf{P}$ is the ground truth.
$\hat{\mathbf{C}}$ is the predicted $W/8\times H/8\times 4$ classification
tensor, and $\mathbf{C}$ the ground truth. For both $\mathcal{L}_{loc}$ and
$\mathcal{L}_{id}$ cross-entropy loss is used.

###  3.2 Sub-pixel Centers Estimation

Once the circle center locations have been detected at a pixel-level
resolution, in this second stage we refine these points to a sub-pixel-level
accuracy for triangulation and accurate pose estimation. With the scaled-back
pixel-level points from the SuperPoint-like network to the original
resolution, we extract $120\times 120$ patches around each detected center
from the original $1280\times 1024$ image. The ellipse is supposed within
those patches for sub-pixel center estimation.

Since circles project as ellipses onto the image plane, we propose a
lightweight encoder-decoder network called EllipSegNet for segmenting the
extracted ellipse patches for contour estimation and sub-pixel centroid
calculation as shown Fig. 3. EllipSegNet is meant for real-time applications,
and its weight is just 1.1 MB, and a forward pass takes 2.3 ms approximately
with an NVIDIA Tesla P100. The encoder of the network consists of eight
$3\times 3$ convolution layers with 16-16-32-32-64-64-64-64 channels,
respectively. From the second convolution layer, three $2\times 2$ max pool
layers are used between each pair of convolutional layers. The decoder part of
EllipSegNet consists of bilinear upsampling with a factor of 2 and skip
connections via concatenation followed by two convolution operations. There
are a total of six $3\times 3$ convolution layers with 64-32-32-16-16-16
channels, respectively. Finally, a $1\times 1$ convolution layer is used to
map the 16 channels to 1 for binary segmentation.

The binary cross-entropy loss is used in this stage for the training of the
network,

$\mathcal{L}(\mathbf{Y},\hat{\mathbf{Y}})=-\sum_{i}\sum_{j}y_{ij}\log
p(\hat{y}_{ij})+(1-y_{ij})\log[1-p(\hat{y}_{ij})]\enspace,$ (2)

where $\hat{\mathbf{Y}}$ is the predicted mask, and $\mathbf{Y}$ the ground
truth. Furthermore, $p$ is the sigmoid function,

$p(x)=\dfrac{1}{1+e^{-x}}\enspace.$ (3)

For sub-pixel point estimation, we can compute the centroid of the white
pixels of the mask with OpenCV’s moments function, or we can use the contour
of the ellipse to fit an ellipse with the OpenCV fitEllipse function. Either
way, the detected center coordinates are estimated at sub-pixel accuracy. We
adopt using fitEllipse, because we obtain more robust and accurate results
than the other approach.

For sub-pixel estimation, the ellipse can be partially inside of the
$120\times 120$ patch. When there is more than one ellipse within the patch,
the nearest to the image center is segmented. Although ellipse contour
estimation is performed pixel-wise (discrete way), the sub-pixel calculation
is not discrete, i.e., the estimation is not restricted to a discrete grid of
positions, which enables higher accuracy. Furthermore, this sub-pixel
estimation is carried out at full resolution, which leads to highly accurate
pose estimation results because if we perform sub-pixel detection in a lower
resolution, we need to scale-back the points to the original size. In this
process, the small detection error in the lower resolution can be also scaled,
and incur pose estimation errors. For example, a sub-pixel detection error
less than a pixel at $320\times 240$ can become larger than a pixel after
scaling-back to the original resolution of $1280\times 1024$. Therefore, an
error in low resolution can increase in the original resolution and this is
critical in final pose estimation for high-accuracy applications.

###  3.3 Marker Pose Estimation

In the two previous sections, we described the sub-pixel 2D location
estimation of the target’s centers for a single image. This section describes
the last stage of MarkerPose which consists of estimating the pose of the
marker. This work is focused on pose estimation with a stereo vision system.
With the estimated IDs, we solve the point correspondence problem. Hence, with
the calibrated stereo vision system and the correspondences
$\mathbf{x}_{1}\leftrightarrow\mathbf{x}_{2}$, we can recover the 3D position
of the points. Following the camera pinhole model, the projection of a 3D
point $\mathbf{X}$ in both image planes can be written as,

$\mathbf{x}_{1}=\mathbf{P}_{1}\mathbf{X}\enspace,\qquad\mathbf{x}_{2}=\mathbf{P}_{2}\mathbf{X}\enspace,$
(4)

where $\mathbf{P}_{1}$ and $\mathbf{P}_{2}$ are the projection matrices for
the cameras, estimated through calibration. To solve for $\mathbf{X}$, the
homogeneous solution is given by

$\begin{bmatrix}x_{1}\mathbf{p}_{1}^{3\mathsf{T}}-\mathbf{p}_{1}^{1\mathsf{T}}\\\
y_{1}\mathbf{p}_{1}^{3\mathsf{T}}-\mathbf{p}_{1}^{2\mathsf{T}}\\\
x_{2}\mathbf{p}_{2}^{3\mathsf{T}}-\mathbf{p}_{2}^{1\mathsf{T}}\\\
y_{2}\mathbf{p}_{2}^{3\mathsf{T}}-\mathbf{p}_{2}^{2\mathsf{T}}\\\
\end{bmatrix}\mathbf{X}=\mathbf{0}\enspace,$ (5)

where $\mathbf{p}_{1}^{i}$ and $\mathbf{p}_{2}^{i}$ represent the $i$-th row
of $\mathbf{P}_{1}$ and $\mathbf{P}_{2}$ respectively, as column vectors. We
solve this with the OpenCV’s triangulate function.

With the three 3D points of the target, we can define a 3D coordinate system.
For that, we define the $c_{0}$ center as the origin of the target frame. The
3D position of this point gives us the translation vector $\mathbf{t}$ of the
target frame with respect to the cameras. Furthermore, the circles $c_{1}$ and
$c_{2}$ represent the direction of the $x$-axis and $y$-axis of the target
coordinate system, respectively. The unit vector $\hat{\mathbf{x}}$ of the
marker frame can be estimated with $c_{0}$ and $c_{1}$ 3D centers. Similarly,
the unit vector $\hat{\mathbf{y}}$ with $c_{0}$ and $c_{2}$ 3D coordinates.
Using cross product, the unit vector of the $z$-axis, that is perpendicular to
the target, is calculated
$\hat{\mathbf{z}}=\hat{\mathbf{x}}\times\hat{\mathbf{y}}$. Finally, the
rotation matrix of the target frame relative to the stereo vision system is
$\mathbf{R}=[\hat{\mathbf{x}}\;\>\hat{\mathbf{y}}\;\>\hat{\mathbf{z}}]$. In
this way, we track the position $\mathbf{t}$ and orientation $\mathbf{R}$ of
the target coordinate system. For accurate triangulation results, we
compensate the 2D sub-pixel points for lens distortion before triangulation.

###  3.4 Detection Method with a Classical Approach

We also developed a center detection strategy using classical computer vision
techniques. This method was designed to aid the labeling of training data and
for comparison with the proposed approach using convolutional neural networks.
This classical approach consists of binarizing and determining the image
contours. As the outer black circle and internal white circles are concentric,
their centers must have approximately the same image coordinates. We exploit
this criterion to estimate the contour of the black circle and its center.
This approach provides highly accurate results, but it is not robust. We need
controlled environments for successful detection results. Furthermore, it
requires parameter tuning to make it work in different environments and
lighting conditions. Finally, as our stereo vision system is calibrated, we
exploit epipolar constraints for solving the matching problem of points for
pair of images, required for pose estimation.

##  4 Implementation Details

This section explains the details about the datasets, training of the
SuperPoint-like and EllipSegNet networks, and the experimental setup.

###  4.1 Training Dataset

We acquired a total of 5683 grayscale images with a $1280\times 1024$
resolution. There are multiple positions and orientations of the target in our
dataset, and the images were acquired at different lighting conditions.
Furthermore, we also used two different target dimensions. Two example images
are shown in Fig. 4. For training the first stage of MarkerPose with the
SuperPoint-like network, these images were resized to $320\times 240$. For
supervision, the center estimation was performed using the classical approach
introduced in Section 3.4. For data augmentation, we randomly applied the
following operations: vertical and horizontal flip, random affine
transformations, contrast modifications, motion blur, Gaussian blur, and
Gaussian noise.

Figure 4: Example of training images for the SuperPoint-like network.

For EllipSegNet training, we extract $120\times 120$ patches from some of the
images used for training the SuperPoin-like network, resulting in 11010
patches. Fig. 5 shows some patch examples, where there are images where a
single ellipse is in the patch and cases where we have more than one ellipse.
We applied the same data augmentation transformations as described above for
the keypoints detection network to ensure the robustness of the entire
pipeline.

Figure 5: Example of training images for EllipSegNet.

###  4.2 Models Training

We implemented the point detection and ellipse segmentation networks with
PyTorch. A version was implemented with Python and another with C++ using
LibTorch. For the SuperPoint-like network and the lightweight EllipSegNet, we
used the Adam optimizer with a learning rate of 0.0001.

###  4.3 Stereo Vision Setup

The experimental stereo vision setup consisted of two monochromatic Basler
cameras acA1300-200um (1280x1024, 203 fps), and two lenses Edmund Optics 58001
(12 mm fixed focal length). We calibrated our stereo vision system with OpenCV
and a pattern of asymmetric circles. A total of 21 different poses of the
pattern were acquired for calibration. We obtained a mean reprojection error
of 0.1230 px and 0.1284 px for the first and second views, respectively. These
errors are small and appropriate for accurate pose estimation.

##  5 Experiments

For validation, we compared MarkerPose with the classical approach discussed
in Section 3.4. For the experiments presented in this section, the distance
between $c_{0}$ and $c_{1}$ centers of the target is 25 mm, and the distance
between $c_{0}$ and $c_{2}$ centers is 40 mm.

###  5.1 Validation Experiments with a Robotic Arm

We evaluated the proposed pose estimation system with two experiments using an
EPSON C3 robot (0.02 mm of precision). The robot was used for assessing the
tracking accuracy. We attached the target to the arm’s end effector, and we
evaluated the displacement and rotation of the target moved by the robot.

####  5.1.1 Displacement Evaluation

For this experiment, the end effector was displaced a fixed distance of
$\Delta x=10$ mm in a total of 20 positions. We measured the distance between
consecutive frames, which give us 19 displacement estimations of each center
of the target, i.e., a total of 57 displacement measurements. Fig. 6 shows an
example of the first 10 frames used for this experiment, acquired with the
left camera (first view).

Figure 6: First 10 frames of the sequence for distance evaluation. A total of
20 positions were captured.

Fig. 7 shows the results of the estimated relative distance between the 19
pairs of consecutive frames with MarkerPose and the classical approach. The
measurements can be compared with the reference displacement line of 10 mm.
Our proposed approach achieves high-accuracy results compared to the classical
method.

Figure 7: Displacement results from the 19 pairs of frames of MarkerPose and
the classical method.

We also estimated the mean absolute error (MAE), the RMS error, and the
(maximum error) infinity norm $\ell_{\infty}$ of the errors of the calculated
57 distances of all centers, using both methods. The results are reported in
Table 1, where the results of MarkerPose with MAE and RMS metrics are
comparable to the classical method. The maximum error is slightly larger with
our approach, but a maximum error of 0.1164 mm is low enough for pose
estimation in high precision applications.

Method | MAE (mm) | RMS (mm) |  $\mathbf{\ell_{\infty}}$ (mm)
---|---|---|---
MarkerPose | 0.0461 | 0.0544 | 0.1164
Classical | 0.0446 | 0.0508 | 0.1086
Table 1: Displacement experiment: performance metrics for the 57 distance
measurements.

Finally, we report the same metrics for the 3D distances between the centers
estimated with the triangulated points. In this case, we have 20 measurements
from the 20 acquired frames. These were obtained as the distance
$d(c_{0},c_{1})$ between $c_{0}$ and $c_{1}$ centers and the distance
$d(c_{0},c_{2})$ between $c_{0}$ and $c_{2}$ points, where $d(.)$ is the
Euclidean distance between two 3D points. Table 2 shows the results.
MarkerPose provides an overall better performance than the classical approach
in this case.

| MarkerPose | Classical method
---|---|---
| $d(c_{0},c_{1})$ | $d(c_{0},c_{2})$ | $d(c_{0},c_{1})$ | $d(c_{0},c_{2})$
MAE (mm) | 0.1315 | 0.2976 | 0.1558 | 0.3435
RMS (mm) | 0.1326 | 0.2982 | 0.1569 | 0.3448
$\mathbf{\ell_{\infty}}$ (mm) | 0.1692 | 0.3257 | 0.1833 | 0.4107

Table 2: Displacement experiment: metrics of the distance between the centers
$c_{0}$ and $c_{1}$, and the centers $c_{0}$ and $c_{2}$, using the frames
acquired for displacement estimation. Here, $d(.)$ is the Euclidean distance
between two 3D points.

####  5.1.2 Rotation Evaluation

The reference was established by rotating the robot’s end effector at a fixed
angle of $\Delta\theta=5^{\circ}$ in a sequence of 7 frames. We measured the
angle between pairs of consecutive frames, which gives us 6 rotation
estimations. The calculation of the rotation angle between frames is performed
with the following expression,

$\theta=\arccos\left[\dfrac{\mathrm{Tr}(\mathbf{R}_{12})-1}{2}\right]\enspace,$
(6)

where $\mathbf{R}_{12}$ is the rotation matrix between the first and second
frame, and can be estimated as
$\mathbf{R}_{12}=\mathbf{R}_{1}^{\mathsf{T}}\mathbf{R}_{2}$. That is,
$\mathbf{R}_{12}$ is the rotation matrix of the target coordinate system
between frame 1 and frame 2. The entire sequence captured by the left camera
is shown in Fig. 8.

Figure 8: Entire sequence of 7 frames from the left camera for angle
estimation.

The estimated relative angle between pairs of frames is shown in Fig. 9 for
MarkerPose and the classical method. By comparing with the reference line of
5∘, the proposed approach achieves better results than the classical approach
with less dispersion. Furthermore, the metrics in Table 3 show the same
behavior, where the MAE, RMS error, and the $\ell_{\infty}$ norm are
significantly lower with MarkerPose. These results show the potential of our
approach for highly accurate pose estimation applications.

Figure 9: Rotation angle measurement results of MarkerPose and the classical method. A total of 6 pairs of consecutive frames were used. Method | MAE | RMS | $\mathbf{\ell_{\infty}}$
---|---|---|---
MarkerPose | 0.0322° | 0.0413° | 0.0687°
Classical | 0.1107° | 0.1290° | 0.2317°
Table 3: Rotation experiment: angle metrics evaluation of the six pairs of
frames. Figure 10: Results of our center point detection method for
synthetically applied motion blur and lighting changes to a test input image.

Similarly to the displacement experiment, the 3D Euclidean distance between
centers is reported in Table 4, where using the seven frames acquired, the
distance between $c_{0}$ and $c_{1}$ centers and the distance between $c_{0}$
and $c_{2}$ centers, are estimated and the same metrics are used for these
estimations. The metrics between the $c_{0}$ and $c_{2}$ centers are better
with our proposed method. For the metrics between $c_{0}$ and $c_{1}$, the
classical method is better except for maximum error.

| MarkerPose | Classical method
---|---|---
| $d(c_{0},c_{1})$ | $d(c_{0},c_{2})$ | $d(c_{0},c_{1})$ | $d(c_{0},c_{2})$
MAE (mm) | 0.1537 | 0.2639 | 0.1315 | 0.3001
RMS (mm) | 0.1548 | 0.2642 | 0.1416 | 0.3011
$\mathbf{\ell_{\infty}}$ (mm) | 0.1896 | 0.2867 | 0.2150 | 0.3387

Table 4: Rotation experiment: metrics of the distance between the centers
$c_{0}$ and $c_{1}$, and the centers $c_{0}$ and $c_{2}$, using the frames
acquired for angle estimation.

###  5.2 Evaluation with Synthetic Effects

We evaluated the proposed method under low lighting conditions and motion
blur. The transformations were applied synthetically to the image, and the
results are shown in Fig. 10. For motion blur, six different kernels sizes
were applied to the original image, and center detection was carried out.
Synthetic darkness was applied by scaling the original image with a factor
$\alpha$. These $\alpha$ values range from 0.25 to 0.039. Those results show
the robustness of our point detection method.

The distance estimation experiment was carried out again using the same
dataset for different lighting changes with $\alpha$ values ranging from 0.5
to 0.0625. An $\alpha$ value of 1 does not represent a lighting change. The
results are reported in Table 5. These results show the robustness of
MarkerPose to low lighting conditions, where even for the darkest case
($\alpha=0.0625$) we obtain highly accurate results.

$\alpha$ | MAE (mm) | RMS (mm) | $\ell_{\infty}$
---|---|---|---
1 | 0.0461 | 0.0544 | 0.1164
0.5 | 0.0451 | 0.0539 | 0.1352
0.25 | 0.0476 | 0.0557 | 0.1404
0.125 | 0.0528 | 0.0603 | 0.1313
0.0625 | 0.0568 | 0.0723 | 0.2130
Table 5: Performance metrics for the distance estimation results for different
lighting scaling values $\alpha$ from 0.5 to 0.0625. For $\alpha=1$ there are
no changes in the images.

###  5.3 Pose Estimation for a 3D Freehand Ultrasound System

To showcase the highly accurate pose estimation result from MarkerPose, we
describe an experiment with the 3D freehand ultrasound (US) technique [10]. It
consists of tracking the pose of a US probe while simultaneously acquiring US
images. This procedure allows for mapping the 2D US scans to a volume by using
a predefined transformation obtained via calibration [14].

For this experiment, we used the same stereo vision setup and the US machine
Biocare iS20. To evaluate our proposed pose estimation method, we
reconstructed and estimated a cylinder’s diameter with the US images. This
cylinder was submerged in water for ultrasound visualization into another
cylindrical object for easy acquisition. An example of stereo and US images
during acquisition is shown in Fig. 11, where the measured object is inside.

Figure 11: Example of 3D freehand ultrasound acquisition with the estimated
pose with the MarkerPose method. The acquired ultrasound image shows a ring-
like structure corresponding to the measured cylinder object.

For this experiment, 94 ultrasound images were acquired along the central axis
of the cylinder. We semi-automatically segmented the images for 3D
reconstruction. The reconstruction results with our pose estimation method are
shown in Fig. 12, which is in agreement with the reference shape. Moreover,
with this point cloud, we estimated the cylinder’s outer diameter through
least squares. For that, we use the outer points of the 3D coordinates. The
estimated diameter with the point cloud was 15.999 mm, whereas the cylinder’s
diameter measured with a caliper was 15.82 mm. This result corresponds to an
absolute error of 0.18 mm, which shows the method’s suitability for highly
accurate pose estimation applications.

Figure 12: 3D reconstruction results of the measured cylinder object with the
3D freehand ultrasound technique and the proposed pose estimation method. 94
2D ultrasound slices were used for 3D reconstruction.

###  5.4 Runtime Evaluation

We measured execution time for both of our point detection stages
independently and for the entire MarkerPose pipeline. We use an NVIDIA Tesla
P100 in a cloud computing platform in Python for runtime evaluation. Our first
stage with the SuperPoint-like network takes 2.7 ms approximately for a single
$320\times 240$ image. A forward pass in the second stage with the EllipSegNet
runs in 2.3 ms approximately for a single $120\times 120$ array. Our entire
pose estimation method, for two input $1280\times 1024$ images, including 2D
points undistortion in both views, triangulation, and finally rotation matrix
and translation vector estimation, runs at around 62 FPS. These results show
the suitability of our proposed approach for real-time applications.

The same experiments were carried out with a GPU NIVIDA GTX 960M on a laptop
using Python. A forward pass with the SuperPoint-like network runs at about
32.6 ms. A forward pass with the EllipSegNet takes 6.1 ms approximately. The
entire pose estimation pipeline runs on average at 12 FPS. These results show
that our method can achieve real-time performance even in an environment with
fewer computational resources.

##  6 Conclusions

The need for highly accurate pose estimation is linked to applications that
require reliable object tracking, like in the case of freehand ultrasound and
other computer vision-based navigation and guidance systems. In most cases,
the use of infrared or electromagnetic tracking devices leads to expensive
systems which are not routinely used. Here, we have proposed MarkerPose: a
real-time, highly accurate pose estimation method based on stereo
triangulation using a planar marker detected with a deep neural network
approach. Our multi-stage detection architecture leads to fast marker keypoint
ID identification and sub-pixel keypoint detection at full image resolution by
exploiting ellipse fitting using EllipSegNet, an introduced lightweight
network. We showed the potential of the proposed method for accurate pose
estimation and tracking under different lighting and motion blur conditions,
and with a biomedical application. Future work involves designing novel
markers that further exploit the capabilities of the proposed method.

## Acknowledgement

This work has been partly funded by Universidad Tecnologica de Bolivar project
C2018P005 and Minciencias contract 935-2019. J. Meza thanks Universidad
Tecnologica de Bolivar for a post-graduate scholarship and MinCiencias, and
MinSalud for a “Joven Talento” scholarship.

## References

  * [1] Mykhaylo Andriluka, Umar Iqbal, Eldar Insafutdinov, Leonid Pishchulin, Anton Milan, Juergen Gall, and Bernt Schiele.  Posetrack: A benchmark for human pose estimation and tracking.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 5167–5176, 2018.
  * [2] Ehsan Basafa, Pezhman Foroughi, Martin Hossbach, Jasmine Bhanushali, and Philipp Stolka.  Visual tracking for multi-modality computer-assisted image guidance.  In Medical Imaging 2017: Image-Guided Procedures, Robotic Interventions, and Modeling, volume 10135, page 101352S. International Society for Optics and Photonics, 2017.
  * [3] Alisa JV Brown, Ali Uneri, Tharindu S De Silva, Amir Manbachi, and Jeffrey H Siewerdsen.  Design and validation of an open-source library of dynamic reference frames for research and education in optical tracking.  Journal of Medical Imaging, 5(2):021215, 2018.
  * [4] Liang-Chieh Chen, George Papandreou, Iasonas Kokkinos, Kevin Murphy, and Alan L Yuille.  Deeplab: Semantic image segmentation with deep convolutional nets, atrous convolution, and fully connected crfs.  IEEE transactions on pattern analysis and machine intelligence, 40(4):834–848, 2017.
  * [5] Daniel DeTone, Tomasz Malisiewicz, and Andrew Rabinovich.  Superpoint: Self-supervised interest point detection and description.  In Proceedings of the IEEE conference on computer vision and pattern recognition workshops, pages 224–236, 2018.
  * [6] Liuhao Ge, Zhou Ren, Yuncheng Li, Zehao Xue, Yingying Wang, Jianfei Cai, and Junsong Yuan.  3d hand shape and pose estimation from a single rgb image.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 10833–10842, 2019.
  * [7] Aryaman Gupta, Kalpit Thakkar, Vineet Gandhi, and PJ Narayanan.  Nose, eyes and ears: Head pose estimation by locating facial keypoints.  In ICASSP 2019-2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 1977–1981. IEEE, 2019.
  * [8] Kaiming He, Georgia Gkioxari, Piotr Dollár, and Ross Girshick.  Mask r-cnn.  In Proceedings of the IEEE international conference on computer vision, pages 2961–2969, 2017.
  * [9] Danying Hu, Daniel DeTone, and Tomasz Malisiewicz.  Deep charuco: Dark charuco marker pose estimation.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 8436–8444, 2019.
  * [10] Qinghua Huang and Zhaozheng Zeng.  A Review on Real-Time 3D Ultrasound Imaging Technology.  BioMed research international, 2017:6027029, 2017.
  * [11] Ho Chuen Kam, Ying Kin Yu, and Kin Hong Wong.  An improvement on aruco marker for pose tracking using kalman filter.  In 2018 19th IEEE/ACIS International Conference on Software Engineering, Artificial Intelligence, Networking and Parallel/Distributed Computing (SNPD), pages 65–69. IEEE, 2018.
  * [12] Jisung Kim, Youngdo Jeong, Hyojin Lee, and Hongsik Yun.  Marker-based structural displacement measurement models with camera movement error correction using image matching and anomaly detection.  Sensors, 20(19):5676, 2020.
  * [13] Ji Yang Lee and Cheol-soo Lee.  Path planning for scara robot based on marker detection using feature extraction and, labelling.  International Journal of Computer Integrated Manufacturing, 31(8):769–776, 2018.
  * [14] Jhacson Meza, Pedro Simarra, Sara Contreras-Ojeda, Lenny A Romero, Sonia H Contreras-Ortiz, Fernando Arámbula Cosío, and Andrés G Marrugo.  A low-cost multi-modal medical imaging system with fringe projection profilometry and 3d freehand ultrasound.  Proc. SPIE, 11330:1133004, 2020.
  * [15] Yoshihiro Nakano.  Stereo vision based single-shot 6d object pose estimation for bin-picking by a robot manipulator.  arXiv preprint arXiv:2005.13759, 2020.
  * [16] Tanmay Nath, Alexander Mathis, An Chi Chen, Amir Patel, Matthias Bethge, and Mackenzie Weygandt Mathis.  Using deeplabcut for 3d markerless pose estimation across species and behaviors.  Nature protocols, 14(7):2152–2176, 2019.
  * [17] Joseph Redmon, Santosh Divvala, Ross Girshick, and Ali Farhadi.  You only look once: Unified, real-time object detection.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 779–788, 2016.
  * [18] Shaoqing Ren, Kaiming He, Ross Girshick, and Jian Sun.  Faster r-cnn: towards real-time object detection with region proposal networks.  IEEE transactions on pattern analysis and machine intelligence, 39(6):1137–1149, 2016.
  * [19] C Romero, C Naufal, J Meza, and Andrés G Marrugo.  A validation strategy for a target-based vision tracking system with an industrial robot.  J. Phys.: Conf. Series, 1547(1):012018–8, 2020.
  * [20] Denys Rozumnyi, Jiri Matas, Filip Sroubek, Marc Pollefeys, and Martin R Oswald.  Fmodetect: Robust detection and trajectory estimation of fast moving objects.  arXiv preprint arXiv:2012.08216, 2020.
  * [21] Paul-Edouard Sarlin, Daniel DeTone, Tomasz Malisiewicz, and Andrew Rabinovich.  Superglue: Learning feature matching with graph neural networks.  In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pages 4938–4947, 2020.
  * [22] Yongbin Sun, Yimin Deng, Haibin Duan, and Xiaobin Xu.  Bionic visual close-range navigation control system for the docking stage of probe-and-drogue autonomous aerial refueling.  Aerospace Science and Technology, 91:136–149, 2019.
  * [23] He Wang, Srinath Sridhar, Jingwei Huang, Julien Valentin, Shuran Song, and Leonidas J Guibas.  Normalized object coordinate space for category-level 6d object pose and size estimation.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 2642–2651, 2019.
  * [24] Yu Xiang, Tanner Schmidt, Venkatraman Narayanan, and Dieter Fox.  Posecnn: A convolutional neural network for 6d object pose estimation in cluttered scenes.  arXiv preprint arXiv:1711.00199, 2017.