# Revealing a charge-density-wave gap in the predicted weak topological
insulator HoSbTe

J. L. Liu  Center for Advanced Quantum Studies and Department of Physics,
Beijing Normal University, Beijing 100875, China  R. Liu  Center for Advanced
Quantum Studies and Department of Physics, Beijing Normal University, Beijing
100875, China  M. Yang  Beijing National Laboratory for Condensed Matter
Physics and Institute of Physics, Chinese Academy of Sciences, Beijing 100190,
China  Center of Materials Science and Optoelectronics Engineering, University
of Chinese Academy of Sciences, Beijing 100049, China  School of Physical
Sciences, University of Chinese Academy of Sciences, Beijing 100190, China  L.
Y. Cao  Center for Advanced Quantum Studies and Department of Physics, Beijing
Normal University, Beijing 100875, China  B. X. Gao  Center for Advanced
Quantum Studies and Department of Physics, Beijing Normal University, Beijing
100875, China  L. Wang  Center for Advanced Quantum Studies and Department of
Physics, Beijing Normal University, Beijing 100875, China  A. F. Fang
Department of Physics, Beijing Normal University, Beijing 100875, China  Y. G.
Shi  Beijing National Laboratory for Condensed Matter Physics and Institute of
Physics, Chinese Academy of Sciences, Beijing 100190, China  Center of
Materials Science and Optoelectronics Engineering, University of Chinese
Academy of Sciences, Beijing 100049, China  School of Physical Sciences,
University of Chinese Academy of Sciences, Beijing 100190, China  Z. P. Yin
Center for Advanced Quantum Studies and Department of Physics, Beijing Normal
University, Beijing 100875, China  R. Y. Chen  Center for Advanced Quantum
Studies and Department of Physics, Beijing Normal University, Beijing 100875,
China

###### Abstract

HoSbTe was predicted to be a weak topological insulator, whose spin-orbit
coupling (SOC) gaps are reported to be as large as hundreds of meV. Utilizing
infrared spectroscopy, we find that the compound is of metallic nature from
350 K down to 10 K. Particularly, both of its itinerant carrier density and
scattering rate are demonstrated to decrease with temperature cooling, which
is responsible for the appearance of a broad hump feature in the temperature
dependent resistivity around 200 K. More importantly, we reveal the appearance
of a charge-density-wave (CDW) gap in addition to the SOC related gap. The
energy scale of the CDW gap is identified to be 364 meV at 10 K, which shift
to 252 meV at 350 K. The coexistence of CDW and SOC gaps in the same compound
paves a new avenue to explore more intriguing physics.

##  I introduction

Topological materials with nontrivial electronic bands have drawn tremendous
interests in condensed matter physics due to their exotic physical properties
and their potential applicationsQi and Zhang (2011); Hasan and Kane (2010);
Nayak  _et al._ (2008). Huge progress has been made in the field and a large
number of topological materials have been recognized. However topological
materials with strong electronic correlation are much less explored and might
host more interesting physicsPesin and Balents (2010); Deng  _et al._ (2017);
Schoop  _et al._ (2018); Peters  _et al._ (2018); Fu and Kane (2008); Wang and
Zhang (2013). The $WHM$ ($W$=Zr, Hf, or Lanthanides, $H$=Si, Ge, Sn, or Sb,
and $M$=O, S, Se, or Te) family, initially noticed by theorist in the search
for two-dimensional topological insulatorsXu  _et al._ (2015), can serve as a
good platform to study the interactions between topology, magnetism and other
emergent instabilities.

Most of the $WHM$ compounds are reported to be Dirac nodal line semimetalsGuo
_et al._ (2019); Liu  _et al._ ; Wang  _et al._ (2021); Yen  _et al._ (2019);
Nakamura  _et al._ (2019); Hu  _et al._ (2016); Fu  _et al._ (2019) or weak
topological insulatorsXu  _et al._ (2015); Lou  _et al._ (2016). However when
$W$ are Lanthanides with $f$ electrons, magnetic orderings and Kondo effect
could play a role in the system. For example, LaSbTe was suggested to be a
topological insulatorSingha  _et al._ (2017) or nodal-line semimetalWang  _et
al._ (2021) by different experiments. Its counterpart CeSbTe with $f$
electrons from Ce is found to be a low-carrier-density Kondo semimetal with a
Kondo temperature of 10 KLv  _et al._ (2019), and its electronic structure
could be severely modified by long range magnetic ordersSchoop  _et al._
(2018). Meanwhile GdSbTe is demonstrated to be a nodal-line semimetal with a
robust Dirac-like band structure across the AFM transition at 13 KHosen  _et
al._ (2018).

Among these $WHM$ materials HoSbTe was theoretically predicted to be a weak
topological insulatorYang  _et al._ (2020), which is corroborated afterwards
by angle-resolved photoemission spectroscopy (ARPES) experiments exhibiting
energy gaps much larger than 200 meV along certain momentum directionsYue  _et
al._ (2020). However the temperature dependent resistivity of HoSbTe exhibits
a broad hump-like feature at around 200 K, indicating a bad-metal-like state
at low temperatures, the underlying physics of which is yet to be revealedYang
_et al._ (2020). The 4$f$ electrons of Ho atoms are believed to be generally
localized and far away from the Fermi level, but an extremely large Sommerfeld
coefficient $\gamma\sim$ 382.2 mJ/mol-1/K2 is obtained by a specific heat
measurementYang  _et al._ (2020), which is much larger than that of CeSbTe
($\gamma$ = 41 mJ/mol-1/K2)Lv  _et al._ (2019), and comparable to many typical
heavy fermion materials like CeCoIn5Petrovic  _et al._ (2001). This seems to
infer a very large Kondo hybridization strength and thus the itinerancy of the
$f$ electrons in HoSbTe.

In order to resolve the existing enigma and further investigate the role
played by 4$f$ electrons, we performed infrared spectroscopy measurements on
single crystalline HoSbTe samples. Our results suggest that the hump feature
in resistivity could be well explained by the competing effect of carrier
density and scattering rate. In addition to the expected SOC gap we observe
for the first time another temperature dependent gap whose value varies from
364 meV at 10K to 252 meV at 250K. We propose that this newly observed gap is
caused by CDW ordering similar to those in CeSbTeLi  _et al._ (2021).

##  II experimental

Single crystalline HoSbTe samples were synthesized by Sb-flux methodYang  _et
al._ (2020). The temperature dependent resistivity $\rho(T)$ and magnetization
$\chi(T)$ up to 350 K were measured in a Quantum Design physical property
measurement system.

Infrared spectroscopic studies were performed with a Bruker IFS 80 V in the
frequency range from 30 to 50 000 cm-1 , on the as-grown shiny surfaces of
HoSbTe which is the $ab$ plane of this quasi-two dimensional compound. In the
measurement of frequency dependent reflectivity $R(\omega)$, either gold or
aluminum coating techniques are adopted in order to eliminate the impact of
microscopic surface textures of the single crystalline compound. The real part
of the optical conductivity $\sigma_{1}(\omega)$ are derived through Kramers-
Kronig transformation of the reflectivity R($\omega$), which is extrapolated
by a Hagen-Rubens relation to zero at the low frequency end, and by
$\omega^{-1.5}$ from 50 000 cm-1 to 800 000 cm-1 , and then by $\omega^{-4}$
for frequency higher than 800 000 cm-1 .

Theoretical calculations were carried out using the linearized-augmented
plane-wave method implemented in WIEN2KBlaha  _et al._ . The Perdew-Burke-
Ernzerhof version of generalized-gradient-approximation (GGA) to the exchange
correlation functionalPerdew  _et al._ (1996) was used. To move the $f$
orbitals of Ho atoms away from Fermi level, GGA+U method was adopted in our
calculations, and we applied an effective Hubbard $U_{eff}=U-J_{H}$ = 7 eV to
the Ho f orbitals. Because the AFM transition temperature of HoSbTe ($T_{N}$
$\sim$ 4 K) is very low, we ignore the magnetic order and force the
calculations to non-magnetic solution. In self-consistent calculation we used
27$\times$27$\times$12 k-mesh grid, and in order to clearly show Fermi
surfaces sections, 500$\times$500 k-mesh grid was used in each specified $kz$
plane. Spin-orbit coupling (SOC) was taken into account in our calculation.

##  III results and discussion

Figure 1: The in-plane (a) resistivity and (b) magnetic susceptibility as a
function of temperature. The magnetic susceptibility is measured with H = 1000
Oe. The first order derivatives of resistivity and magnetic susceptibility are
displayed in the inset of (a) and (b), respectively.

Fig. 1(a) and (b) display the temperature dependent resistivity $\rho(T)$ and
magnetization $\chi(T)$, respectively. Although we have measured $\rho(T)$ and
$\chi(T)$ up to a higher temperature, there is not much difference with the
previous reportYang  _et al._ (2020): $\rho(T)$ shows a broad hump around 200
K and $\chi(T)$ follows the Curie-Wiess law down to at least 100 K; an
antiferromagnetic phase transition occurs at around 4 K. In order to check if
there are any other phase transitions that may lead to the ``insulator to
metal'' change inferred by resistivity, the first-order derivatives
$d\rho(T)/dT$ and $d\chi(T)/dT$ are calculated and displayed in the inset of
Fig. 1(a) and (b), respectively. The smooth derivative curves at relatively
high temperatures obviously do not support existence of such transitions.

Figure 2: The frequency dependent (a) reflectivity R($\omega$) and (b) optical
conductivity $\sigma_{1}(\omega)$ at several temperatures up to 10 000 cm-1 .
The two insets show R($\omega$) and $\sigma_{1}(\omega)$ in an expanded range
up to 50 000 cm-1 , respectively. (c) The renormalized optical conductivity as
a function of frequency up to 5000 cm-1 . (d ) The spectral weight as a
function of frequency. The inset shows renormalized SW at different
temperatures.

The frequency dependent reflectivity $R(\omega)$ at several temperatures is
shown in Fig. 2 (a). In the far-infrared (IR) region, metallic responses in
the whole temperature range are clearly evidenced by the increasing of
$R(\omega)$ with decreasing energy and its approaching to unity at zero
frequency. Although HoSbTe was predicted to be a weak topological insulator,
it could still be a metal since the Fermi level may not lie in the topological
gaps and may pass through other trivial energy bands. A diamond-shaped Fermi
surface (FS) around $\Gamma$ is indeed observed by ARPES experiments on
HoSbTeYue  _et al._ (2020). However this is somewhat inconsistent with the
insulating behavior of $\rho(T)$ above 200 K, which will be discussed in the
following.

In the mid-IR region, a platform-like feature can be clearly seen in
$R(\omega)$ at 350 K, which is gradually suppressed at lower temperatures and
ultimately replaced by two dip structures as indicated by the red arrows in
Fig. 2(a). Such kind of low energy absorptions in $R(\omega)$ are usually
linked to the loss of density of state at the Fermi level. The growing
spectral weights of the two dips with temperature cooling further imply that
they might be related to some kind of gap openings. For even higher
frequencies, $R(\omega)$ decreases with increasing frequency up to 18 000 cm-1
and exhibits only mild temperature dependence. A few weak humps could be
resolved in this frequency range.

The evolution of the electronic structure could be more straightforwardly
reflected by the optical conductivity $\sigma_{1}(\omega)$, which is directly
related to the joint density of states for electron transitions from occupied
to unoccupied states. As displayed in Fig. 2(b), the direct current value of
$\sigma_{1}(\omega\rightarrow 0)$ first reduces slightly with decreasing
temperature above 250K, then increases rapidly down to 10 K, roughly matching
with the result of $\rho(T)$. The presence of the Drude components centered at
zero frequency is another evidence of the metallic nature of HoSbTe. In
correspondence with the dip structures in $R(\omega)$, two extremely
pronounced peaks could be identified in the optical conductivity
$\sigma(\omega)$ at 10 K, where a crossover can be observed at around 2400
cm-1 between the two peaks. As temperature increases, the two individual peaks
tend to get closer to each other and seem to turn into a very broad single
peak above room temperature. Considering the two-dip structure in $R(\omega)$
can persist up to 300 K, we believe the broad peak in $\sigma_{1}(\omega)$ is
actually composed by two independent peaks.

To better analyze the evolvement of the two peaks, we plot the renormalized
optical conductivity spectra in Fig.2(c). It is clearly shown that the energy
scale of both peaks shift to higher frequency as temperature decreasing.
Supposing they are originated from SOC gaps as revealed by the ARPES
experiments or trivial interband transitions, the spectral weight of them
should undergo a subtle decrease upon temperature cooling, due to the
reduction of thermal excitations. While the lower energy peak is slightly
suppressed as expected, the higher energy one gets much more pronounced when
temperature decreases, suggesting a totally different origination. We further
present the spectral weight SW=$\int_{0}^{\omega}\sigma_{1}(\omega)d\omega$ as
a function of frequency and temperature in Fig.2(d). At the low frequencies
end, the SW increases as temperature cooling due to the narrowing of the Drude
component. Then the higher temperature SW increases and surpasses the lower
temperature value roughly above 700 cm-1 . At even higher energies, SW at
different temperatures almost merge together above 5000 cm-1 , complying with
the conservation law. Additionally, the renormalized SW reaches a minimum at
around 1200 cm-1 at 10 K, as can be seen in the inset of Fig.2(d), indicating
the SW transfer from frequencies lower than 1200 cm-1 to higher energy region.
Notably, these behaviors are similar with some CDW materials, such as CeTe3Hu
_et al._ (2011).

In order to unravel the individual mechanism of the two peaks, we performed
first-principals calculations on the optical conductivity
$\sigma_{1-cal}(\omega)$ and band structure of HoSbTe, as shown in Fig.3(b)
and (d) respectively. Surprisingly, there are some distinct differences
between the calculated and experimental optical conductivity: the overall
calculated spectra are much higher than the experimental data, partially due
to the high-energy extrapolation of $R(\omega)$; $\sigma_{1-cal}(\omega)$ only
shows a single peak in the mid-IR region at around 1450 cm-1 , which
corresponds to interband transition across an SOC associated gap, as indicated
by the light green arrow in Fig.3(d). These disagreements clearly indicate a
yet unknown reconstruction of the electronic structure. Nevertheless, as the
lowest energy interband transition, the SOC gap in $\sigma_{1-cal}(\omega)$
locates very close to the lower energy peak in the experimental data, pointing
to a similar origination.

Figure 3: (a) Crystal structure of HoSbTe. A unit cell is indicated by black
lines. (b) The calculated conductivity $\sigma_{1-cal}$ and experimental data
at 10 K. (c) Brillouin zone of HoSbTe. (d) The calculated bulk band structure
with orbital characters for HoSbTe. The light green arrow marks the optical
interband transition that give rise to the mid-IR peak in the optical
conductivity. Te 5$p$, Sb 5$p$ and Ho 4$f$ orbital characters are presented by
red, blue, and green circle, respectively.

We further analyze the optical conductivity quantitatively by decomposing it
into several different terms according to the Drude-Lorentz model:

$\sigma_{1}(\omega)=\frac{\omega_{p}^{2}}{4\pi}\frac{\gamma_{D}}{\omega^{2}+\gamma_{D}^{2}}+\sum_{j}\frac{S_{j}^{2}}{4\pi}\frac{\gamma_{j}\omega^{2}}{(\omega_{j}^{2}-\omega^{2})^{2}+\omega^{2}\gamma_{j}^{2}},$
(1)

where $\omega_{p}$ and $\gamma_{D}$ are the plasma frequency and scattering
rate of the itinerant carriers, while $\omega_{j}$, $\gamma_{j}$ and $S_{j}$
are the resonant frequency, damping and strength of the $j$th Lorentz
oscillators, respectively. The first term at the right side of equation (1) is
the Drude component modeling intraband transitions of free electrons while the
second Lorentz term describes interband transitions across energy gaps. A
presentative example of the fitting results of $\sigma_{1}(\omega)$ at 10 K
are plotted in Fig.4, where one Durde and seven Lorentz terms are employed.
Part of the fitting parameters are listed in Table 1.

Figure 4: The Drude-Lorentz model fitting result of optical conductivity $\sigma_{1}(\omega)$ at 10 K, where one Drude term and seven Lorentz terms are used. T | $\omega_{P}$ | $\gamma_{D}$ | $\omega_{1}$ | $\gamma_{1}$ | $S_{1}$ | $\omega_{2}$ | $\gamma_{2}$ | $S_{2}$
---|---|---|---|---|---|---|---|---
10K | 11311(2) | 269 | 1647 | 1276(2) | 19455(16) | 2934 | 2173(1) | 24967
100K | 12417(2) | 305 | 1492 | 1224(1) | 19830(14) | 2700 | 2183(1) | 24163(14)
200K | 12983(2) | 455 | 1405 | 1504(1) | 21539(15) | 2523 | 2177(1) | 22211(15)
250K | 14698(7) | 752 | 1363 | 1465(2) | 18848(23) | 2432 | 2223 | 23560(16)
300K | 19175 | 1246 | 1396 | 1730 | 18581 | 2335 | 2171 | 20089
| (46) | (6) | (1) | (6) | (83) | (1) | (1) | (39)
350K | 18066 | 1077 | 1089 | 1938 | 18925 | 2032 | 2247 | 21267
| (1855) | (223) | (5) | (53) | (1900) | (2) | (5) | (135)
Table 1: The fitting parameters at different temperatures. $\omega_{P}$ is the
plasma frequency and $\gamma_{D}$=1/$\tau_{D}$ is the scattering rate of the
Drude term. The center frequency $\omega_{j}$, width $\gamma_{j}$=1/$\tau_{j}$
and square root of the oscillator strength $S_{j}$ of the first two Lorentz
components are listed. The error bars of all fitting parameters are included.

Among the fitting parameters, the plasma frequency $\omega_{p}$ is related to
the itinerant carrier density $n$ and effective mass $m^{*}$ by
$\omega_{p}^{2}=4\pi ne^{2}/m^{*}$. The variation of $\omega_{p}$ with
temperature is shown in Fig.5 (a), where $\omega_{p}$ decreases monotonically
from 19 175 cm-1 at 300 K to 11 311 cm-1 at 10 K, indicating either a
continuous losing of the free carrier density or a constant enhancement of the
effective mass. It is worth to remark that the data of 350 K seems to exhibit
a deviation from the above mentioned trend. This is because the Drude pattern
of $\sigma_{1}(\omega)$ at 350 K is very broad and highly overlap with the
first Lorentz term. Therefore, an extremely huge error bar is generated in the
fitting procedure. The temperature dependent scattering rate $\gamma_{D}$,
which is the width at half maximum of the Drude spectra, is plotted in Fig.5
(a) as well. Notably, $\gamma_{D}$ evolves in a very similar fashion with
$\omega_{p}$. It is well known that the electrical conductivity of a material
is proportional to its density of free carriers $n$ and relaxation time
$\tau=1/\gamma_{D}$. Here for HoSbTe, $\tau$ obviously increases with lowering
temperature, in favor of a diminishing resistivity. As a consequence, how the
free carrier density $n$ evolve with temperature is crucial to explain the
hump feature in $\rho(T)$, which will be further elaborated later.

Figure 5: (a) The temperature dependent plasma frequency $\omega_{p}$ and
scattering rate $\gamma_{D}$ . (b) The energy scales of the two gaps as a
function of temperature.

The two prominent peaks of $\sigma_{1}(\omega)$ in the mid-IR region could be
well fitted by two Lorentz oscillators, corresponding to transitions across
two different gaps. The obtained resonant frequencies, as displayed in Fig.5
(b), can be considered as the excitation gap energies, which are identified to
be $2\Delta_{1}\simeq$ 204 meV and $2\Delta_{2}\simeq$ 364 meV at 10 K,
respectively. In a previous reportYue  _et al._ (2020), ARPES measurements
have revealed SOC gaps as large as hundreds of meV, which seems to be
consistent with both $2\Delta_{1}$ and $2\Delta_{2}$. As discussed above,
although this SOC scenario or trivial interband transitions can work for the
lower energy gap, they can not explain the drastic decrease of $2\Delta_{2}$
to 252 meV at 350 K and the increase of its SW at lower temperatures. Taking
the closeness between the calculated SOC gap and $2\Delta_{1}$ into
consideration, we propose this lower energy gap to be an SOC gap, while the
higher energy gap requires an alternative mechanism which will be discussed in
the following.

In consideration of the exceptional large Sommerfeld coefficient $\gamma$Yang
_et al._ (2020), there is a possibility that the higher energy gap arises from
the Kondo effect, in which case the interaction between located 4$f$ electrons
and conduction bands would generate a mid-IR peak. The Drude spectral weight
as well as the plasma frequency $\omega_{p}$ will be severely suppressed owing
to the enhancement of effective mass $m^{*}$. This is in excellent agrement
with our spectroscopic results. Even the drop of $\rho(T)$ below 200 K could
be well justified because the initial localized $f$ electrons will be enclosed
by the Fermi surface and contribute carrier density as a result of Kondo
screening, which reduces the scattering rate of itinerant carriers at the same
time. However, even if we adopt a constant carrier density, the enhancement of
$m^{*}$ from 350 K to 10 K is estimated to be only about
$m^{*}_{10K}/m^{*}_{350K}=\omega^{2}_{p(350K)}/\omega^{2}_{p(10K)}\simeq
2.55$, which disagrees violently with the huge Sommerfeld coefficient $\gamma$
= 382.2 mJ/mol-1/K2Yang  _et al._ (2020). Besides, the sustaining of
$2\Delta_{2}$ up to 350 K would indicate an even higher Kondo temperature
$T_{K}$, below which the magnetic susceptibly is supposed to deviate from the
Curie-Weiss law along with a downward bending resistivity. Such characters is
obviously absent in our measurements. Hence we have ruled out the existence of
Kondo effect in this system .

Another option for $2\Delta_{2}$ is a density-wave gap, since the transfer of
spectral weight from the free electron Drude response to higher energy peaks
agrees perfectly with the electrodynamic prediction of density-wave orderings
with a case-I coherent factor. Notably, a charge density wave gap about 0.3 eV
is observed in the Dirac semimetal CeSbTeLi  _et al._ (2021), driven by
electron-phonon coupling due to FS nesting. The FSs of CeSbTe are constituted
by a small enclosed hole pocket at $\Gamma$ and two 2D rhombus-shaped sheets,
the latter of which contain largely parallel segments that could be connected
by a single wave vector, in favor of FS nesting and thus the development of
density wave orderings. For comparison, we plot the FS section of HoSbTe at
several specified $kz$ planes, as shown in Fig.6. It also contains two
rhombus-shape sheets at $kz$=0, but the two sheets are closer to each other
and disconnected along the $\Gamma-M$ direction due to SOC gaps. When $kz$
increases, they gradually shrinks and eventually disappears as the SOC gaps
shift closer to the Fermi level. As a result, the two sheets are
indistinguishable within the energy resolution limitation of ARPES
measurementsYue  _et al._ (2020). In reality, however, it is most likely the
nearly perfect nesting of these two sheets induces a density wave instability,
and the corresponding single particle energy gap is captured by our infrared
spectroscopy. As there are no additional signatures in support of spin-
density-wave instabilities, like AFM orderings above 10 K, we speculate that
$2\Delta_{2}$ is actually a charge density wave gap. Since the two rhombus-
shape are so close to each other, an extremely small nesting wave vector is
expected.

Figure 6: Fermi surface sections in (a) $kz$ = 0, (b) $kz$ = $\pi$/3c, (c)
$kz$ = 2$\pi$/3c and (d) $kz$ = $\pi$/c plane of HoSbTe. Red, green and blue
curves represent Fermi surfaces formed by different bands which are divided by
the order of energy. The black lines mark the boundaries of Brillouin zone.

Assuming that $m^{*}$ is temperature independent and equal to the electron
rest mass $m_{e}$, $n$ is estimated to be about 3.62 $\times$ 1021 $cm^{-3}$
at 350 K, which reduce to 1.42 $\times$ 1021 $cm^{-3}$ at 10 K, meaning a
60.8$\%$ loss of the FS due to the gap enlargement. The residual FS at the
lowest temperature manifests that the CDW gaps are partially opened. As
addressed above, the scattering rate $\gamma_{D}$ also decreases with
temperature lowering. Consequently, the broad hump feature appeared in
$\rho(T)$ is probably caused by the competition of the simultaneous reduction
of $n$ and $\gamma_{D}$, which contribute oppositely to the resistivity. This
phenomena is also observed in the three dimensional Dirac semimtal ZrTe5Chen
_et al._ (2015).

Although we can principally ascertain the occurrence of a CDW ordering in
HoSbTe, we have not been able to determine its phase transition temperature
and ordering wave vector from our current results. The large gap energies at
350 K suggest that the CDW transition takes place at a much higher
temperature, which is inaccessible for our current measurement. In the weak
coupling mean field theory framework, the gap energy 2$\Delta$ =
3.52$k_{B}T_{C}$ for strict 1D systems, where $T_{C}$ is the phase transition
temperature. By this relation the CDW $T_{C}$ of HoSbTe would be over 2000 K
with gap value 2$\Delta_{2}=364$ meV. However $T_{C}$ will be severely
suppressed when the effect of fluctuations is taken into account, especially
in strongly anisotropic materials. Therefore the 2$\Delta$/$k_{B}T_{C}$ ratio
is typically much larger than 3.52. For example, another Ho-based material
HoTe3 undergoes a CDW transition at 288 K with an energy gap of 380 meVPfuner
_et al._ (2010), very close to the energy scale of $2\Delta_{2}$. We therefore
conjecture that the transition temperature of HoSbTe is probably of the same
order of magnitude with HoTe3.

Finally we would like to compare HoSbTe with LaSbTe and CeSbTe. Although all
of the three compounds bear strikingly similar FSsYue  _et al._ (2020); Wang
_et al._ (2021); Li  _et al._ (2021), their physical properties are much
different from each other. Our results on HoSbTe and previous reports on
CeSbTe demonstrate that these two materials have CDW instabilities, while
LaSbTe is reported to be a nodal line semimetalWang  _et al._ (2021) with no
CDW observed by so far. For CeSbTe the CDW gap is estimated to be about 0.3 eV
with a transition temperature above 338 K and a weak Kondo effect is also
identifiedLi  _et al._ (2021). We observe a CDW gap of about 364 meV in HoSbTe
with no Kondo effect. This rich phase diagram supplies an ideal platform to
investigate the influence of chemical pressure on the CDW ordering and its
interplay with Kondo effect by substituting $Ln$SbTe with different lanthanide
elements.

##  IV conclusion

In summary, we have investigated the charge carrier dynamics of a predicted
topological insulator compound HoSbTe. The temperature dependent reflectivity
and optical conductivity both exhibit metallic behaviors from 350 K to 10 K,
with a decreasing carrier density and scattering rate, the competing
contribution of which gives rise to the hump feature in resistivity. We
identify an expected SOC gap of around 204 meV, and an additional temperature
dependent gap whose value shifts from 364 meV at 10 K to 252 meV at 350 K.
This enormous change along with its sizable increase of spectral weight upon
cooling are clearly incompatible with SOC band gaps or trivial interband
transitions. As Kondo hybridization is also excluded as the origin of this
gap, due to a minor enhancement of the effective mass, we ultimately propose
it to stem from a CDW ordering. This complex system together with its LnSbTe
cousins provide a promising arena to explore the entanglement of CDW, Kondo
effect and topological orders.

ACKNOWLEDGMENTS

This work was supported by the National Science Foundation of China (No.
12074042), the Young Scientists Fund of the National Natural Science
Foundation of China (Grant No. 11704033), the National Key Research and
Development Program of China (No. 2016YFA0302300), and the Fundamental
Research Funds for the Central Universities.

## References

  * Qi and Zhang (2011) X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83, 1057 (2011).
  * Hasan and Kane (2010) M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
  * Nayak  _et al._ (2008) C. Nayak, S. H. Simon, A. Stern, M. Freedman, and S. Das Sarma, Rev. Mod. Phys. 80, 1083 (2008).
  * Pesin and Balents (2010) D. Pesin and L. Balents, Nature Physics 6, 376 (2010).
  * Deng  _et al._ (2017) M.-X. Deng, W. Luo, R.-Q. Wang, L. Sheng, and D. Y. Xing, Phys. Rev. B 96, 155141 (2017).
  * Schoop  _et al._ (2018) L. M. Schoop, A. Topp, J. Lippmann, F. Orlandi, L. Müchler, M. G. Vergniory, Y. Sun, A. W. Rost, V. Duppel, M. Krivenkov, S. Sheoran, P. Manuel, A. Varykhalov, B. Yan, R. K. Kremer, C. R. Ast, and B. V. Lotsch, Science Advances 4, eaar2317 (2018).
  * Peters  _et al._ (2018) R. Peters, T. Yoshida, and N. Kawakami, Phys. Rev. B 98, 075104 (2018).
  * Fu and Kane (2008) L. Fu and C. L. Kane, Phys. Rev. Lett. 100, 096407 (2008).
  * Wang and Zhang (2013) Z. Wang and S.-c. Zhang, Physical Review B 87, 161107(R) (2013).
  * Xu  _et al._ (2015) Q. Xu, Z. Song, S. Nie, H. Weng, Z. Fang, and X. Dai, Physical Review B 92, 205310 (2015).
  * Guo  _et al._ (2019) L. Guo, T.-W. Chen, C. Chen, L. Chen, Y. Zhang, G.-Y. Gao, J. Yang, X.-G. Li, W.-Y. Zhao, S. Dong, and R.-K. Zheng, ACS Applied Electronic Materials 1, 869 (2019).
  * Liu  _et al._ (0) X. Liu, C. Yue, S. V. Erohin, Y. Zhu, A. Joshy, J. Liu, A. M. Sanchez, D. Graf, P. B. Sorokin, Z. Mao, J. Hu, and J. Wei, Nano Letters 0, null (0).
  * Wang  _et al._ (2021) Y. Wang, Y. Qian, M. Yang, H. Chen, C. Li, Z. Tan, Y. Cai, W. Zhao, S. Gao, Y. Feng, S. Kumar, E. F. Schwier, L. Zhao, H. Weng, Y. Shi, G. Wang, Y. Song, Y. Huang, K. Shimada, Z. Xu, X. J. Zhou, and G. Liu, Physical Review B 103, 125131 (2021).
  * Yen  _et al._ (2019) Y. Yen, C.-L. Chiu, P.-H. Lin, R. Sankar, F. Chou, T.-M. Chuang, and G.-Y. Guo, ``Dirac nodal line and rashba splitting surface states in nonsymmorphic zrgete,'' (2019), arXiv:1912.07002 .
  * Nakamura  _et al._ (2019) T. Nakamura, S. Souma, Z. Wang, K. Yamauchi, D. Takane, H. Oinuma, K. Nakayama, K. Horiba, H. Kumigashira, T. Oguchi, T. Takahashi, Y. Ando, and T. Sato, Phys. Rev. B 99, 245105 (2019).
  * Hu  _et al._ (2016) J. Hu, Z. Tang, J. Liu, X. Liu, Y. Zhu, D. Graf, K. Myhro, S. Tran, C. N. Lau, J. Wei, and Z. Mao, Phys. Rev. Lett. 117, 016602 (2016).
  * Fu  _et al._ (2019) B.-B. Fu, C.-J. Yi, T.-T. Zhang, M. Caputo, J.-Z. Ma, X. Gao, B. Q. Lv, L.-Y. Kong, Y.-B. Huang, P. Richard, M. Shi, V. N. Strocov, C. Fang, H.-M. Weng, Y.-G. Shi, T. Qian, and H. Ding, Science Advances 5, eaau6459 (2019).
  * Lou  _et al._ (2016) R. Lou, J.-Z. Ma, Q.-N. Xu, B.-B. Fu, L.-Y. Kong, Y.-G. Shi, P. Richard, H.-M. Weng, Z. Fang, S.-S. Sun, Q. Wang, H.-C. Lei, T. Qian, H. Ding, and S.-C. Wang, Phys. Rev. B 93, 241104(R) (2016).
  * Singha  _et al._ (2017) R. Singha, A. Pariari, B. Satpati, and P. Mandal, Phys. Rev. B 96, 245138 (2017).
  * Lv  _et al._ (2019) B. Lv, J. Chen, L. Qiao, J. Ma, and X. Yang, J.phys.:condens.matter 31, 355601 (2019).
  * Hosen  _et al._ (2018) M. M. Hosen, G. Dhakal, K. Dimitri, P. Maldonado, A. Aperis, F. Kabir, C. Sims, P. Riseborough, P. M. Oppeneer, D. Kaczorowski, T. Durakiewicz, and M. Neupane, Scientific Reports 8, 13283 (2018).
  * Yang  _et al._ (2020) M. Yang, Y. Qian, D. Yan, Y. Li, Y. Song, Z. Wang, C. Yi, H. L. Feng, H. Weng, and Y. Shi, Physical Review Materials 4, 094203 (2020).
  * Yue  _et al._ (2020) S. Yue, Y. Qian, M. Yang, D. Geng, C. Yi, S. Kumar, K. Shimada, P. Cheng, L. Chen, Z. Wang, H. Weng, Y. Shi, K. Wu, and B. Feng, Physical Review B 102, 155109 (2020).
  * Petrovic  _et al._ (2001) C. Petrovic, P. G. Pagliuso, M. F. Hundley, R. Movshovich, J. L. Sarrao, J. D. Thompson, Z. Fisk, and P. Monthoux, Journal of Physics: Condensed Matter 13, L337 (2001).
  * Li  _et al._ (2021) P. Li, B. Lv, Y. Fang, W. Guo, Z. Wu, Y. Wu, and D. Shen, Science China Physics, Mechanics & Astronomy 64, 237412 (2021).
  * (26) P. Blaha, K. Schwarz, G. K. Madsen, D. Kvasnicka, and J. Luitz, WIEN2K, An Augmented Plane Wave+ Local Orbitals Program for Calculating Crystal Properties, (Karlheinz Schwarz, Vienna University of Technology, Austria, 2001) .
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Hu  _et al._ (2011) B. F. Hu, P. Zheng, R. H. Yuan, T. Dong, B. Cheng, Z. G. Chen, and N. L. Wang, Physical Review B 83, 155113 (2011).
  * Chen  _et al._ (2015) R. Y. Chen, S. J. Zhang, J. A. Schneeloch, C. Zhang, Q. Li, G. D. Gu, and N. L. Wang, Physical Review B 92, 075107 (2015).
  * Pfuner  _et al._ (2010) F. Pfuner, P. Lerch, J.-H. Chu, H.-H. Kuo, I. R. Fisher, and L. Degiorgi, Phys. Rev. B 81, 195110 (2010).


# Anomaly Detection of Underwater Gliders
Verified by Deployment Data ††thanks:  The research work is supported by ONR
grants N00014-19-1-2556 and N00014-19-1-2266; AFOSR grant FA9550-19-1-0283;
NSF grants GCR-1934836, CNS-2016582 and ITE-2137798; and NOAA grant
NA16NOS0120028.

Ruochu Yang  School of Electrical and Computer Engineering
Georgia Institute of Technology
Atlanta, United States  Mengxue Hou  College of Engineering
Purdue university
West Lafayette, United States  Chad Lembke  College of Marine Science
University of South Florida
St.Petersburg, United States  Catherine Edwards  Skidaway Institute of
Oceanography
University of Georgia
Savannah, United States  Fumin Zhang  School of Electrical and Computer
Engineering
Georgia Institute of Technology
Atlanta, United States

###### Abstract

This paper utilizes an anomaly detection algorithm to check if underwater
gliders are operating normally in the unknown ocean environment. Glider pilots
can be warned of the detected glider anomaly in real time, thus taking over
the glider appropriately and avoiding further damage to the glider. The
adopted algorithm is validated by two valuable sets of data in real glider
deployments, the University of South Florida (USF) glider Stella and the
Skidaway Institute of Oceanography (SkIO) glider Angus.

###### Index Terms:

anomaly detection, glider navigation

##  I Introduction

Autonomous underwater profiling gliders, hereafter referred to as gliders,
have been widely used in oceanography for a range of applications, including
harmful algal bloom monitoring, improving hurricane prediction, passive and
active acoustics, ocean observing systems, and targeted science questions [1,
2, 3, 4, 5]. Similar progress has been made to improve the performance and
value of the data collected through coordinated control of vehicles over long
durations [6, 7, 8]. Despite these advances, unpredictable events like shark
or vessel strikes [9], wing loss, or interaction with echneids (remoras) can
lead to abnormal flight behavior, prematurely end a glider’s mission, or even
result in the loss of the vehicle. Human pilots cannot always directly track
gliders’ abnormal behavior from the subsets of measured data that are
telemetered to shore, especially when gliders confront external disturbances
[10, 11], and sending more data or adding monitoring devices to detect
abnormal behavior is not always feasible or helpful. Therefore, anomaly
detection algorithms can be developed to aid human pilots, and their design
must rely on a limited subset of glider data.

Anomaly detection has been studied extensively in the robotics community [12,
13, 14]. Specifically in the domain of marine robots, some anomaly detection
algorithms focus on individual components that frequently degrade like
propellers, thrusters, and rotors [15, 16, 17, 18]. Other anomaly detection
algorithms utilize robot motion like monitoring pitch angle and depth angle to
detect possible deviation from expected motion [19]. However, research to date
has not yet resolved issues such as inaccessible hardware dynamics,
incapability of performing online detection, and lack of experimental
verification. Detected anomalies must also be checked against the case of
false alarm [20]. One common scenario occurs when flow speed exceeds the
maximum robot speed; the performance of marine robots is degraded in terms of
speed over the ground [6], but the anomaly occurs due to environmental
conditions rather than flight performance of the robot itself. Therefore, it
is significant to compare marine robot speed with ocean flow speed; this
comparison is not trivial since flow speed is not directly measured by the
glider but is instead estimated through dead reckoning.

Controlled Lagrangian particle tracking (CLPT) is a theoretical framework that
can simultaneously estimate ocean flow speed and robot speed [21]. Based on
these two estimates, it is feasible to detect abnormal motion while providing
robust false alarm avoidance by excluding unexpected ocean flow. The anomaly
detection algorithm applied in this paper is adopted from the work [22], which
uses trajectory data and heading angle data that are available and generic in
all types of gliders. This similar technique of leveraging glider navigational
data can also be seen in previous work [23, 24, 25]. Given a trajectory and
heading angles, the adopted algorithm generates both the estimated flow speed
and the estimated glider speed in real time. The estimated glider speed can be
compared with the normal speed range to check for anomalies, and the estimated
flow speed is compared with the glider-estimated flow speed to check if the
detected anomaly is a false alarm.

In general, anomaly detection is difficult to validate in field experiments.
In abnormal conditions, gliders are likely to encounter issues of GPS search,
thus imposing risk and cost to human rescue. Anomalies occur unexpectedly,
there is lack of “ground truth” to determine the major factor causing the
anomaly. Further, anomalies can accumulate over time, and it can be difficult
for glider pilots to notice in real time even with close monitoring during
field experiments. Given only a subset of data, small changes in flight
performance can be too subtle for pilots to identify. The main contribution of
this paper is to verify the anomaly detection algorithm in [22] using real-
life glider deployment data for the first time. The algorithm can infer
whether or not the glider may have had a shark hit, wing loss, or remora
attachment simply from a subset of flight data rather than large amounts of
data. The University of South Florida (USF) glider Stella and the Skidaway
Institute of Oceanography (SkIO) glider Angus provide two valuable sets of
experimental data, in which anomalies can be attributed to marine bio-hazards
and strong evidence is available to confirm their cause. Both data sets show
promising detection results to validate the anomaly detection algorithm in
[22].

This paper is organized as follows. Section II describes the working
principles of Slocum glider and the experiment setup of two coastal glider
deployments in the Gulf of Mexico and south east Atlantic Ocean. Section III
illustrates the framework of the anomaly detection algorithm. Section IV
verifies the algorithm by showing results of detected anomalies. Section V
provides conclusions and future work.

##  II Glider and Experiment Setup

This section introduces the Slocum glider, its method of navigation, and the
details of deployments of two gliders, Stella and Angus, that were subjected
to anomalous conditions while deployed.

###  II-A Slocum Glider

The Slocum glider is an autonomous underwater vehicle (AUV) that moves by
changing its buoyancy and center of gravity [26]. Capable of collecting data
in the ocean for long duration missions of up to 6 months without a recharge,
gliders typically travel at low speed (approximately $0.25–0.35m\cdot
s^{-1}$). During the deployment, the glider surfaces at pre-defined intervals
to transmit flight and science data, as well as receive updates and other
instructions via satellite connection. In order to minimize time at the
surface, where the glider is vulnerable to ship traffic and other dangers,
only a small portion of heavily subsetted flight data are telemetered to shore
in small binary data (SBD) files. For example, in The trajectory and heading
information of sbd file data is processed in real time by the anomaly
detection algorithm presented in Section III. After the mission is over, all
the full data sets are downloaded off the glider into Dinkum binary files (DBD
and EBD for flight and science data, respectively). The anomaly detection
algorithm is validated in Section IV by comparing the anomaly detected from
the real-time SBD file date with the anomaly directly shown from the full DBD
file data.

While in mission, the glider estimates its horizontal position using pitch and
depth data through dead reckoning. From one surfacing to the next, the
deviation the glider makes from the dead reckoned estimate is attributed to
the influence of flow. The difference between these two positions are used to
estimate the average flow speed along the glider trajectory. GPS positions at
the beginning and end of the surfacing are used to adjust the effect of
surface wind drift on estimated flow velocity. This glider-estimated flow is
checked against the flow estimate given by the anomaly detection algorithm in
Section III.

###  II-B Field Deployments

Data from two gliders will be used for field validation of the detection
anomaly algorithm: the glider Stella, operated by the University of South
Florida, and the glider Angus, operated by Skidaway Institute of Oceanography,
part of the University of Georgia. Stella ($\\#$772) was deployed by the
University of South Florida in the Gulf of Mexico on the West Florida Shelf
near $(27.6970^{\circ}N,83.8778^{\circ}W)$ on July 01, 2021 UTC, transited
offshore and then back inshore, and ended at the location
$(27.7318^{\circ}N,83.4136^{\circ}W)$ on July 19, 2021 UTC. The Google Earth
trajectory of Stella is shown in Fig. 1.

Figure 1: Google Earth trajectory of the 2021 Stella deployment.

As part of an acoustic telemetry and navigation experiment, the glider Angus
was deployed near Gray’s Reef National Marine Sanctuary. Its navigation was
supported by the Glider Environment Networked Information system (GENIoS)
[27], which generated optimal waypoints for the glider to use based on real-
time glider-estimated currents and position at each surfacing. The mission
started at the location $(31.4633^{\circ}N,80.7757^{\circ}W)$ at April 04,
2022, 15:30 UTC, and ended at the location
$(31.5072^{\circ}N,80.7368^{\circ}W)$ at April 25, 2022, 18:51 UTC. The Google
Earth trajectory of Angus is shown in Fig. 2.

Figure 2: Google Earth trajectory of the 2022 Angus deployment.

##  III Anomaly Detection Algorithm

The adopted anomaly detection algorithm can simultaneously generate glider
speed estimate and flow speed estimate in real time. If the glider speed
estimate is within the normal speed range, it is assumed that no anomaly has
occurred. Otherwise, the glider may not be operating normally. The algorithm-
estimated flow speed is compared with glider-estimated flow speed to compute
flow estimation error. If the error is too large, the detected anomaly is
false alarm. This section illustrates the main framework of anomaly detection
algorithm. More theoretical details can be found in the work [22].

The glider motion model is shown in (1),

$\displaystyle\dot{x}=F_{R}(x,t)+V_{R}(t)\Psi_{c}(t)$
$\displaystyle\Psi_{c}(t)=[\cos\psi_{c}(t),\sin\psi_{c}(t)]^{T}$ (1)

where $F_{R}$ is the actual flow field, $x$ is the true glider position,
$V_{R}$ is the true glider speed, and $\psi_{c}$ is the glider heading angle.

In practice, the glider speed is limited due to the glider control input
constraint. Equation (1) can be modified as the combination of (2) and (3),

$\displaystyle\dot{x}=F_{R}(x,t)+u\Psi_{c}(t),$ (2)

where

$\displaystyle u=\begin{cases}V_{R},&\mbox{if }V_{R}\leq u_{0}\\\
u_{0},&\mbox{if }V_{R}>u_{0}\end{cases}$ (3)

The maximum glider speed $u_{0}$ is determined by the hardware configuration
of the glider, as well as the water depth of operation. As shown in (4), the
flow field can be represented by spatial and temporal basis functions
(specifically Gaussian radial and tidal functions) [28],

$\displaystyle F_{R}(x,t)=\theta\phi(x,t)$ (4)

where

$\displaystyle\theta=\begin{bmatrix}\theta_{1}\\\
\theta_{2}\end{bmatrix}=\begin{bmatrix}\theta_{1}^{1}&\cdots&\theta_{1}^{N}\\\
\theta_{2}^{1}&\cdots&\theta_{2}^{N}\end{bmatrix}\in\mathbb{R}^{2\times N}$
(5)

is the unknown parameter to be estimated,

$\displaystyle\phi=\begin{bmatrix}\phi^{1}(x,t)&\cdots&\phi^{N}(x,t)\end{bmatrix}^{T}$
(6)
$\displaystyle\phi^{i}(x,t)=\exp^{-\frac{||x-c_{i}||}{2\sigma_{i}}}\cos(\omega_{i}t+\upsilon_{i})$
(7)

is the combined basis function, $c_{i}$ is the center, $\sigma_{i}$ is the
width, $\omega_{i}$ is the tidal frequency, $\upsilon_{i}$ is the tidal phase,
and $N$ is the selected number of basis functions. The heading $\Psi_{c}(t)$
and the true trajectory $x(t)$ are known for all time $t$ from the sbd file
data operationally sent by the glider. Given the true trajectory $x(t)$, the
adopted algorithm can estimate the unknown parameter $\theta$ in (4) and the
glider speed $V_{R}$ in (1) in real time. Define

$\displaystyle\xi(t)=\begin{bmatrix}\xi_{1}(t)\\\
\xi_{2}(t)\end{bmatrix}=\begin{bmatrix}\xi_{1}^{1}(t)&\cdots&\xi_{1}^{N}(t)\\\
\xi_{2}^{1}(t)&\cdots&\xi_{2}^{N}(t)\end{bmatrix}\in\mathbb{R}^{2\times N}$
(8)

as the estimate of the flow parameter $\theta$. Define $V_{L}(t)$ as the
estimate of the glider speed $V_{R}$. The adopted algorithm can generate
converged estimates, i.e., $\lim_{t\to\infty}\xi(t)\to\theta$ and
$\lim_{t\to\infty}V_{L}(t)\to V_{R}$ to ensure that the maximum trajectory
estimation error (CLLE) converges to zero. The adopted algorithm also designs
three gains $K$, $\bar{\gamma}_{1}$ and $\bar{\gamma}_{2}$ to accelerate the
estimating process.

Next, the glider speed estimate $V_{L}(t)$ can be used to decide whether
anomalies happened or not, and the flow speed estimate $\xi(t)$ can be used to
avoid false alarm. The adopted algorithm assumes that the maximum glider speed
$V_{max}$ and the minimum glider speed $V_{min}$ are known beforehand. If the
estimated glider speed $V_{L}(t)\in[V_{min},V_{max}]$, it is assumed that no
anomaly happens to the glider. Otherwise, the glider may not be operating
normally. In order to avoid false alarm, the adopted algorithm checks the
estimated flow speed as well. Let us introduce $F_{L}(t)$ as the flow speed
estimated by the adopted algorithm, and $F_{M}(t)$ as the flow speed estimated
by the glider. In practice, prior information like ocean models or sensor
measurements can help generate the glider-estimated flow speed $F_{M}(t)$. The
algorithm-estimated flow speed $F_{L}(t)$ is compared with the glider-
estimated flow speed $F_{M}(t)$ to determine if the detected anomaly is false
alarm or not. Defining $||F_{M}(t)-F_{L}(t)||$ as the flow estimation error,
the criteria $p_{E}$ to evaluate the flow estimation error is given as (9),

$\displaystyle
p_{E}=\frac{||F_{M}(t)-F_{L}(t)||}{2max(\hat{F}_{Lmax},\hat{F}_{Mmax})}$ (9)

where $\hat{F}_{Lmax}=max(||F_{L}(\tau)||_{\tau\in[0,t]})$ is the maximum
algorithm-estimated flow speed until time $t$, and
$\hat{F}_{Mmax}=max(||F_{M}(\tau)||_{\tau\in[0,t]})$ is the maximum glider-
estimated flow speed until time $t$. If $p_{E}>\gamma_{f}$, the detected
anomaly should not be trusted, where $\gamma_{f}$ is a false alarm threshold
selected beforehand. The pseudo-code is shown in Algorithm 1. For convenience,
we introduce $f_{a}$ as anomaly detection flag shown in (10).

$\displaystyle f_{a}=\begin{cases}0&\mbox{no anomaly}\\\ 1&\mbox{anomaly
detected}\\\ 2&\mbox{false alarm}\end{cases}$ (10)

Algorithm 1 Anomaly Detection

Input: Algorithm-estimated flow speed $F_{L}(t)$, glider-estimated flow speed
$F_{M}(t)$, false alarm threshold $\gamma_{f}$, estimated glider speed
$V_{L}(t)$, maximum glider speed $V_{max}$, minimum glider speed $V_{min}$
Output: Anomaly detection flag $f_{a}$

1: $\hat{F}_{Lmax}=max(||F_{L}(\tau)||_{\tau\in[0,t]})$

2: $\hat{F}_{Mmax}=max(||F_{M}(\tau)||_{\tau\in[0,t]})$

3: $p_{E}=\frac{||F_{M}(t)-F_{L}(t)||}{2max(\hat{F}_{Lmax},\hat{F}_{Mmax})}$

4: if $p_{E}>\gamma_{f}$ then

5: $f_{a}=2$

6: else if $V_{L}(t)>V_{max}$ or $V_{L}(t)<V_{min}$ then

7: $f_{a}=1$

8: else

9: $f_{a}=0$

10: end if

##  IV Experimental Results

For the purpose of validation, the adopted anomaly detection algorithm is
applied to two field deployments with verified anomalies in the coastal ocean
of the Gulf of Mexico and south east Atlantic Ocean. The estimated glider
speed is compared with the maximum and minimum glider speeds to decide if
anomalies can be identified. The estimated flow speed is checked against the
glider-estimated flow speed to avoid false alarm. For verification, the
anomaly detected by the algorithm is compared with the anomaly seen from
glider dbd file data and pilot logs. For reference, the designed parameters of
the algorithm are shown in TABLE I.

TABLE I: Parameters of Experiments Parameters |  Stella Deployment Data |  Angus Deployment Data
---|---|---
number of basis functions N |  4 |  4
width $\sigma_{i}$ |  1.3e4 |  2e4
tidal phase $\upsilon_{i}$ |  0 |  0
tidal frequency $\omega_{i}$ |  $2\pi$e-6 |  $2\pi$e-6
gain K |  $\begin{bmatrix}0.003&0\\\ 0&0.003\end{bmatrix}$ |  $\begin{bmatrix}0.002&0\\\ 0&0.002\end{bmatrix}$
gain $\bar{\gamma}_{1}$ |  1e-7 |  2e-5
gain $\bar{\gamma}_{2}$ |  35e-10 |  2e-8
false alarm threshold $\gamma_{f}$ |  0.6 |  0.7

###  IV-A Results of Stella Deployment

Upon recovery July 19th, 2021, the glider Stella missing one of its two
detachable wings, and scratches on the glider hull and the shark tooth
embedded in the thermoplastic aft glider cowling (Figs. refshark scratch and
3(b)) suggest a serious shark strike. Glider pilot analysis of the full DBD
data post deployment suggests that the shark hit of the glider likely occurred
around July 07, 2021, 16:00 UTC when it was descending at about 55m depth.
This post-mission analysis shows that the glider went from a $-4deg$ roll to
an $8deg$ roll, and finally settled back at about a $0deg$ roll (Fig. 4(a) and
Fig. 4(b)). The sudden change of roll and the consistent change through the
remainder of the deployment suggests that the wing was lost during this
process.

(a) Hull scratches documented after recovery.

(b) Shark tooth found embedded in the hull after recovery.

Figure 3: Post-recovery evidence of shark strike.

(a) Glider-measured roll (degrees) from post-recovery DBD data.

(b) Glider-measured depth (m) from post-recovery DBD data.

Figure 4: ground truth (Stella deployment).

Applying the anomaly detection algorithm to the telemetered SBD data from the
glider, the adopted algorithm can guarantee that the estimated trajectory
converges to the true trajectory since the persistent excitation matrix $W(t)$
is not singular. The comparison of the estimated trajectory with the true
trajectory is shown in Fig. 5. It should be noted that the 4 green circles
cover the whole trajectory, which means that the 4 combined basis functions
cover the whole flow fields. This coverage is necessary for the parameter
estimation convergence. As shown in Fig. 6, the maximum trajectory estimation
error (CLLE) is only $1.7m$, which is small enough considering the glider
moved approximately 0.25m/s over O(100)km in its deployment. Therefore, it can
be concluded that the CLLE converges to zero.

Figure 5: Comparison of the estimated (blue) and reported “true” (red)
trajectory for the 2021 Stella deployment.

Figure 6: CLLE (m) for the 2021 Stella deployment.

Given the CLLE convergence, the glider speed estimate and the flow speed
estimate should also converge. In this case, the anomaly can be detected based
on these two estimates. For better comparison, the flow is divided into $u$
(West-East, W-E, or zonal) and $v$ (North-South, N-S, or meridional)
components. As shown in Fig. 7(a), the algorithm-estimated W-E flow is close
to the corresponding glider-estimated W-E flow, which means that the $u$ flow
estimation error is small. Similarly as the N-S flow shown in Fig. 7(b). This
comparison suggests that the anomaly detection flag can be trusted when it is
signaled.

(a) W-E flow component.

(b) N-S flow component.

Figure 7: Comparison of glider-estimated and algorithm-estimated W-E ($u$,
upper) and N-S ($v$, lower) velocities.

The estimated glider speed is compared with the reported “true” speed value
collected by the glider. If the estimated speed is not in the normal speed
range, the anomaly is assumed to have occurred. As shown in Fig. 8, the
estimated glider speed drops out of the normal speed range (green dot line) at
around July 07, 2021, 18:00 UTC, and keeps dropping afterwards. As shown in
Fig. 9, the flag value switches from 0 to 1 at around July 07, 2021, 18:00 UTC
as well. The flag value never changes to 2 with the false alarm threshold
$\gamma_{f}=0.6$. The timestamp around which the anomaly is detected by the
adopted algorithm corresponds to the timestamp from both the glider team’s
report and the full post-deployment DBD file data. Therefore, the adopted
algorithm is verified valid and accurate by successfully detecting the
anomaly.

Figure 8: Comparison of estimated (red) and reported “true” (blue) speed of
Stella in the 2021 deployment.

Figure 9: Anomaly detection flag as calculated by the adopted algorithm for
the 2021 Stella deployment.

###  IV-B Results of Angus Deployment

According to the mission report, pilots noted that the glider Angus
periodically flew slower than expected, and suffered from unusually poor
satellite data transmission, with multiple broken connections per surfacing,
consistent with remora strikes. The problem was first identified by pilots
around April 17, 2022, 00:00 UTC. Since remoras are negatively buoyant, the
glider’s Iridium antenna is not fully out of the water, and multiple animals
attached to the instrument can induce too much of a buoyancy loss for the
glider’s pump to overcome, preventing flight altogether. The remora attachment
scene is difficult to catch in real-time missions. For illustration, Fig. 10
shows an example image of two remoras attached to a glider (photo credit: Chad
Lembke, College of Marine Science, University of South Florida). Based on the
post-recovery DBD file data, the roll and depth ground truth is shown in Fig.
11(a) and Fig. 11(b), respectively. The timestamp April 17, 2022, 00:00 UTC
when depth and roll data becomes abnormal exactly matches the glider team’s
reported timestamp.

Figure 10: Example image of remora attachment (photo credit: Chad Lembke,
College of Marine Science, University of South Florida).

(a) Glider-measured roll (degrees) from post-recovery DBD data.

(b) Glider-measured depth (m) from post-recovery DBD data.

Figure 11: ground truth (Angus deployment).

Using the SBD file offloaded from the glider, the calculated persistent
excitation matrix $W(t)$ is not singular, which means that the adopted
algorithm can guarantee that the estimated trajectory converges to the true
trajectory. The comparison of the estimated trajectory with the true
trajectory is shown in Fig. 12. It should be noted that the four green circles
cover the whole trajectory, which means that the four combined basis functions
cover the whole flow fields. This coverage is necessary for the parameter
estimation convergence. As shown in Fig. 13, the maximum trajectory estimation
error (CLLE) only reaches around $6m$, which is sufficiently small given the
glider speed and range. Therefore, it can be concluded that the CLLE converges
to zero.

Figure 12: Comparison of the estimated (blue) and reported “true” (red)
trajectory for the 2022 Angus deployment.

Figure 13: CLLE (m) for the 2022 Angus deployment.

Given the CLLE convergence, the glider speed estimate and the flow speed
estimate should converge as well. In this case, the anomaly can be detected
based on these two estimates. For better comparison, the flow is divided into
$u$ (West-East, W-E, or zonal) and $v$ (North-South, N-S, or meridional)
components. As shown in Fig. 14(a) and 14(b), the algorithm-estimated W-E and
N-S flow is close to the corresponding glider-estimated values, so both the
$u$ and $v$ flow estimation errors are small. Therefore, any anomaly detection
flag can be trusted.

(a) W-E flow component.

(b) N-S flow component.

Figure 14: Comparison of glider-estimated and algorithm-estimated W-E ($u$,
upper) and N-S ($v$, lower) velocities.

As shown in Fig. 15, the estimated glider speed sharply drops out of the
normal speed range(green dot line) at around April 16, 2022, 22:00 UTC.
Afterwards, the glider tries to recover its speed but continues below the
normal range. As shown in Fig. 16, the flag value switches from 0 to 1 at
around April 16, 2022, 22:00 UTC as well. The flag value never changes to 2
with the false alarm threshold $\gamma_{f}=0.7$. The timestamp around which
the anomaly is detected by the adopted algorithm corresponds to the timestamp
from both the glider team’s report and the DBD file data. Therefore, the
adopted algorithm is verified valid and accurate, and the anomaly is
successfully detected.

Figure 15: Comparison of estimated speed (red) and normal speed range (green)
of Angus in the 2022 deployment.

Figure 16: Anomaly detection flag as calculated by the adopted algorithm for
the 2022 Angus deployment.

##  V Conclusion

This paper implements and validates an anomaly detection algorithm for two
sets of data in real missions of gliders, the University of South Florida
(USF) glider Stella and the Skidaway Institute of Oceanography (SkIO) glider
Angus, where the anomalies of shark attack and remora attachment can be
verified. Based on the glider trajectory and the heading angle data, the
algorithm generates glider speed estimate and flow speed estimate in real
time. The estimated glider speed is compared with the normal glider speed
range to decide whether anomalies happen or not. The estimated flow speed is
compared with the glider-estimated flow speed to determine if the detected
anomaly is a false alarm. From a theoretical point of view, the current
algorithm focuses on the real-time estimation process by predicting and
updating the trajectory, which is known as the filtering framework, as the
ongoing deployment gives the real glider data as feedback. Future work will
improve the estimation accuracy by integrating the post-deployment glider data
as a smoothing framework. From a practical point of view, future work will
implement algorithms that account for dependence of glider speed on water
depth, and consider the dependency of error on glider heading angle, which can
be used to diagnose poor ballasting of the instrument.

## References

  * [1] T. N. Miles, D. Zhang, G. R. Foltz, J. A. Zhang, C. Meinig, F. Bringas, J. Triñanes, M. Le Hénaff, M. F. Aristizabal Vargas, S. Coakley _et al._ , “Uncrewed ocean gliders and saildrones support hurricane forecasting and research,” _Oceanography_ , vol. 34, no. 4, pp. 78–81, 2021\.
  * [2] P. Testor, B. de Young, D. L. Rudnick, S. Glenn, D. Hayes, C. M. Lee, C. Pattiaratchi, K. Hill, E. Heslop, V. Turpin, P. Alenius, C. Barrera _et al._ , “Oceangliders: A component of the integrated goos,” _Frontiers in Marine Science_ , vol. 6, 2019.
  * [3] R. H. Weisberg, Y. Liu, C. Lembke, C. Hu, K. Hubbard, and M. Garrett, “The coastal ocean circulation influence on the 2018 west florida shelf k. brevis red tide bloom,” _Journal of Geophysical Research: Oceans_ , vol. 124, no. 4, pp. 2501–2512, 2019.
  * [4] P. Cauchy, K. J. Heywood, N. D. Merchant, B. Y. Queste, and P. Testor, “Wind speed measured from underwater gliders using passive acoustics,” _Journal of Atmospheric and Oceanic Technology_ , vol. 35, no. 12, pp. 2305–2321, 2018.
  * [5] M. F. Baumgartner, J. Bonnell, P. J. Corkeron, S. M. Van Parijs, C. Hotchkin, B. A. Hodges, J. Bort Thornton, B. L. Mensi, and S. M. Bruner, “Slocum gliders provide accurate near real-time estimates of baleen whale presence from human-reviewed passive acoustic detection information,” _Frontiers in Marine Science_ , vol. 7, 2020.
  * [6] F. Zhang, D. M. Fratantoni, D. A. Paley, J. M. Lund, and N. E. Leonard, “Control of coordinated patterns for ocean sampling,” _International Journal of Control_ , vol. 80, no. 7, pp. 1186–1199, 2007.
  * [7] D. A. Paley, F. Zhang, and N. E. Leonard, “Cooperative control for ocean sampling: The glider coordinated control system,” _IEEE Transactions on Control Systems Technology_ , vol. 16, no. 4, pp. 735–744, 2008.
  * [8] M. Hou, S. Cho, H. Zhou, C. R. Edwards, and F. Zhang, “Bounded cost path planning for underwater vehicles assisted by a time-invariant partitioned flow field model,” _Frontiers in Robotics and AI_ , vol. 8, 2021.
  * [9] M. Stanway, B. Kieft, T. Hoover, B. Hobson, D. Klimov, J. Erickson, B. Raanan, D. Ebert, and J. Bellingham, “White shark strike on a long-range auv in monterey bay,” 05 2015.
  * [10] J. J. Gertler, _Fault detection and diagnosis in engineering systems_. CRC press, 2017.
  * [11] J. Chen and R. J. Patton, _Robust model-based fault diagnosis for dynamic systems_. Springer Science & Business Media, 2012, vol. 3.
  * [12] B.-Y. Raanan, J. G. Bellingham, Y. Zhang, M. Kemp, B. Kieft, H. Singh, and Y. Girdhar, “Automatic fault diagnosis for autonomous underwater vehicles using online topic models,” in _OCEANS 2016 MTS/IEEE Monterey_. IEEE, 2016, pp. 1–6.
  * [13] D. Park, Z. Erickson, T. Bhattacharjee, and C. C. Kemp, “Multimodal execution monitoring for anomaly detection during robot manipulation,” in _2016 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2016, pp. 407–414.
  * [14] R. Isermann, “Model-based fault-detection and diagnosis–status and applications,” _Annual Reviews in control_ , vol. 29, no. 1, pp. 71–85, 2005\.
  * [15] G. Fagogenis, V. De Carolis, and D. M. Lane, “Online fault detection and model adaptation for underwater vehicles in the case of thruster failures,” in _2016 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2016, pp. 2625–2630.
  * [16] Y.-s. Sun, X.-r. Ran, Y.-m. Li, G.-c. Zhang, and Y.-h. Zhang, “Thruster fault diagnosis method based on gaussian particle filter for autonomous underwater vehicles,” _International Journal of Naval Architecture and Ocean Engineering_ , vol. 8, no. 3, pp. 243–251, 2016.
  * [17] A. Caiti, F. Di Corato, F. Fabiani, D. Fenucci, S. Grechi, and F. Pacini, “Enhancing autonomy: Fault detection, identification and optimal reaction for over—actuated auvs,” in _OCEANS 2015-Genova_. IEEE, 2015, pp. 1–6.
  * [18] M. Caccia, R. Bono, G. Bruzzone, G. Bruzzone, E. Spirandelli, and G. Veruggio, “Experiences on actuator fault detection, diagnosis and accomodation for rovs,” _International Symposiyum of Unmanned Untethered Sub-mersible Technol_ , pp. 3–21, 2001.
  * [19] B. Y. Raanan, J. Bellingham, Y. Zhang, B. Kieft, M. J. Stanway, R. McEwen, and B. Hobson, “A real-time vertical plane flight anomaly detection system for a long range autonomous underwater vehicle,” in _OCEANS 2015 - MTS/IEEE Washington_ , 2015, pp. 1–6.
  * [20] V. Chandola, A. Banerjee, and V. Kumar, “Anomaly detection: A survey,” _ACM computing surveys (CSUR)_ , vol. 41, no. 3, pp. 1–58, 2009.
  * [21] K. Szwaykowska and F. Zhang, “Controlled lagrangian particle tracking: error growth under feedback control,” _IEEE Transactions on Control Systems Technology_ , vol. 26, no. 3, pp. 874–889, 2017.
  * [22] S. Cho, F. Zhang, and C. R. Edwards, “Learning and detecting abnormal speed of marine robots,” _International Journal of Advanced Robotic Systems_ , vol. 18, no. 2, p. 1729881421999268, 2021.
  * [23] B. T. Morris and M. M. Trivedi, “A survey of vision-based trajectory learning and analysis for surveillance,” _IEEE transactions on circuits and systems for video technology_ , vol. 18, no. 8, pp. 1114–1127, 2008.
  * [24] P.-R. Lei, “A framework for anomaly detection in maritime trajectory behavior,” _Knowledge and Information Systems_ , vol. 47, no. 1, pp. 189–214, 2016.
  * [25] O. Rosen and A. Medvedev, “An on-line algorithm for anomaly detection in trajectory data,” in _2012 American Control Conference (ACC)_. IEEE, 2012, pp. 1117–1122.
  * [26] O. Schofield, J. Kohut, D. Aragon, L. Creed, J. Graver, C. Haldeman, J. Kerfoot, H. Roarty, C. Jones, D. Webb, and S. Glenn, “Slocum gliders: Robust and ready,” _Journal of Field Robotics_ , vol. 24, no. 6, pp. 473–485, 2007.
  * [27] D. Chang, F. Zhang, and C. R. Edwards, “Real-time guidance of underwater gliders assisted by predictive ocean models,” _Journal of Atmospheric and Oceanic Technology_ , vol. 32, no. 3, pp. 562–578, 2015.
  * [28] X. Liang, W. Wu, D. Chang, and F. Zhang, “Real-time modelling of tidal current for navigating underwater glider sensing networks,” _Procedia Computer Science_ , vol. 10, pp. 1121–1126, 2012.