# The Sun as a target for axion dark matter detection

Elisa Todarello<EMAIL_ADDRESS>Dipartimento di Fisica,
Università di Torino, Via P. Giuria 1, 10125 Torino, Italy  Istituto Nazionale
di Fisica Nucleare, Sezione di Torino, Via P. Giuria 1, 10125 Torino, Italy
Marco Regis<EMAIL_ADDRESS>Dipartimento di Fisica, Università di
Torino, Via P. Giuria 1, 10125 Torino, Italy  Istituto Nazionale di Fisica
Nucleare, Sezione di Torino, Via P. Giuria 1, 10125 Torino, Italy  Marco Taoso
<EMAIL_ADDRESS>Istituto Nazionale di Fisica Nucleare, Sezione di Torino, Via
P. Giuria 1, 10125 Torino, Italy  Maurizio Giannotti<EMAIL_ADDRESS>Centro de Astropartículas y Física de Altas Energías (CAPA), Universidad de
Zaragoza, Zaragoza, 50009, Spain<EMAIL_ADDRESS>Jaime Ruz<EMAIL_ADDRESS>Centro de Astropartículas y Física de Altas Energías (CAPA), Universidad de
Zaragoza, Zaragoza, 50009, Spain  Julia K. Vogel<EMAIL_ADDRESS>Centro de
Astropartículas y Física de Altas Energías (CAPA), Universidad de Zaragoza,
Zaragoza, 50009, Spain

###### Abstract

The exploration of the parameter space of axion and axion-like particle dark
matter is a major aim of the future program of astroparticle physics
investigations. In this context, we present a possible strategy that focuses
on detecting radio emissions arising from the conversion of dark matter axions
in the Sun’s magnetic field, including conversion in sunspots. We demonstrate
that near-future low-frequency radio telescopes, such as the SKA Low, may
access regions of unexplored parameter space for masses $m_{a}\lesssim
10^{-6}$ eV.

QCD axions are hypothetical particles introduced to solve the strong CP
problem, namely the non-observation of CP violation from QCD interactions [1,
2, 3, 4], see e.g. refs. [5, 6] for recent reviews. They are pseudo-scalar
Nambu-Goldstone bosons, associated with the spontaneous breaking of the so-
called Peccei-Quinn symmetry. Additionally, QCD axions are also a viable dark
matter (DM) candidate, possibly accounting for the entirety of cold DM in the
Universe. They can be produced through different thermal and non-thermal
mechanisms in the early Universe, before or after inflation [7, 8, 9]. The
allowed range of the Peccei-Quinn scale, $f_{a}$, to obtain the observed DM
relic density is $10^{10}-10^{12}$ GeV, corresponding to masses $m_{a}$ in the
rage $10^{-6}-10^{-3}$ eV. The QCD DM axion stands as a key focus in the
ongoing and future program of astroparticle physics investigations. Numerous
experimental searches, both in the laboratory and in the sky have been
designed to look for this particle [10, 11, 12].

Furthermore, several beyond-the-standard-model theories, such as string theory
models [13], predict in many cases one or more axion-like particles (ALPs),
with properties very similar to those of axions but not necessarily related to
the strong CP-problem [14]. ALPs can have masses as low as $10^{-22}$ eV and
are typically very weakly coupled to the SM. Moreover, they are also good DM
candidates in some portions of the parameter space [15].

Experimental searches for both the QCD axion and ALPs often rely on their
coupling to photons [16]. In both cases, the coupling is described by the
effective Lagrangian term $\mathcal{L}\propto
g_{a\gamma}aF_{\mu\nu}\tilde{F}^{\mu\nu}$.

One of the most promising ALP 111From now on, we will use the terms axion and
ALP interchangeably. searches with astrophysical probes involve the radio
signal associated with ALP-photon conversion in magnetized astrophysical
plasma. This signature is based on the above-mentioned ALP-photon coupling,
and it would appear as a nearly monochromatic spectral line. In order to have
a significant probability of conversion, large magnetic fields and dense
plasma are required, so to resonantly amplify the process. In the past few
years, the search has been focused especially on the magnetosphere of neutron
stars [18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34].

In this paper, we instead highlight that the Sun and its sunspots are
promising targets. Sunspots are magnetic structures located within active
regions of the Sun. They appear darker on the solar surface than the normal
solar photosphere (quiet Sun), and are composed of an inner, darker part,
called the umbra, and an outer, less dark part, called the penumbra. It is
well known that the magnetic field can be very intense above the photosphere
of sunspots, reaching a few thousand Gauss in the umbra [35]. The magnetic
field decreases with height in the solar atmosphere, but it could still be
rather strong in the chromosphere and corona. For example, Ref. [36] reported
a magnetic field of about 4000 G at the base of the corona above a sunspot,
with 1000 G at a height of $10^{4}$ km. Sunspots might be environments hosting
copious axion-photon conversions in their magnetic field.

The idea of considering the Sun was proposed in [37] for a similar process,
i.e. investigating dark photon-photon mixing. A search for the associated
monochromatic radio signal in the solar observation data of the LOFAR
telescope was then conducted in [38], which discusses also the application to
the axion-photon conversion.

Besides non-relativistic DM axions, one could also consider the conversion of
relativistic axions thermally produced in the Sun’s interior and then
converted into photons in the Sun’s atmosphere [39]. In this case, they would
lead to X-ray signatures which we plan to investigate in a forthcoming
analysis [40].

This paper aims to provide a simple and general description of the radio
emission arising from the conversion of ambient DM axions in the Sun’s
magnetic field. The goal is to present robust and broad arguments, with
computations that can be easily reproduced by the reader. These arguments
should provide a framework that can later be applied to specific examples in
dedicated studies. As we shall see, the analysis of the radio signal from
axion conversion in the solar magnetic field may provide an opportunity to
access regions of unexplored parameter space with the next generation of radio
telescopes.

For a non-relativistic axion, a weakly magnetized plasma, and assuming the
magnetic field to be static, the conversion probability can be approximated as
[33]

$P_{a\rightarrow\gamma}\simeq\frac{\pi}{2}\,\frac{g_{a\gamma}^{2}\,B_{\perp}^{2}}{v_{a}\,\omega^{\prime}_{q|res}}\;,$
(1)

where $\omega_{q}$ is the plasma frequency and $\omega^{\prime}_{q|res}$ is
its gradient $d\omega_{q}/dr$ evaluated at the resonance location, i.e., at
the radius $r_{c}$ where $\omega_{q}\simeq m_{a}$. In eq. (1) and in the
following we assume spherical symmetry and, for simplicity, radial emission.
The radius $r_{c}$ at which the conversion occurs, for different ALP masses,
can be estimated from the resonance condition and the relation
$\omega_{q}(r)=1.17\,\mu{\rm eV}\,\sqrt{n_{e}(r)/(10^{9}\,{\rm cm^{-3}})}$,
where $n_{e}$ is the electron plasma density.

For a uniform conversion over a surface $\Delta A$, and assuming the ALP
distribution to be isotropic, the ALP-induced photons flux per unit frequency
on Earth can be estimated as

$\displaystyle S$ $\displaystyle=$
$\displaystyle\int\,\frac{d\Omega}{4\pi\,\Delta\nu}\,\rho_{a}\,v_{a}\,P_{a\rightarrow\gamma}\,e^{-\tau}$
(2) $\displaystyle\simeq$ $\displaystyle\frac{\Delta
A}{4\pi\,\,\Delta\nu\,d^{2}}\,\rho_{a}\,v_{a}\,P_{a\rightarrow\gamma}\,e^{-\tau}\;,$

where $d$ is the distance of the conversion surface from us, $\rho_{a}$ is the
DM density at the conversion surface, $\tau$ is the photon optical depth, and
$\Delta\nu$ is the bandwidth of the signal.

To gain some intuition about the expected signal, let us consider a sunspot,
i.e. a structure with strong magnetic fields appearing on the solar surface.
Considering a radius of the sunspot $\ell_{s}$, the conversion area is
computed as $\Delta A=\pi\,\ell_{s}^{2}.$ For an analytical estimate, we can
approximate the plasma profile as a power-law $\omega_{p}\propto h^{\alpha},$
with $\alpha\simeq 0.5$ in the range of conversion radii of interest. Then we
have $\omega^{\prime}_{q|res}=\alpha\,\omega_{p}/h_{c},$ where $h_{c}$ is the
distance of the conversion surface from the photosphere. Plugging in typical
values for a large sunspot [35], we find

$\displaystyle S$ $\displaystyle=$ $\displaystyle\frac{\Delta
A}{8\,\,\Delta\nu\,d^{2}}\,\rho_{a}\,\frac{g_{a\gamma}^{2}\,B^{2}}{\omega^{\prime}_{q|res}}\,e^{-\tau}=0.7\,{\rm
mJy}\,\left(\frac{10^{-6}}{\Delta\nu/\nu}\right)$ (3) $\displaystyle\times$
$\displaystyle\left(\frac{\ell_{s}}{4\times 10^{4}{\rm
km}}\right)^{2}\,\left(\frac{\rho_{a}}{1.0\,{\rm
GeV/cm^{3}}}\right)\left(\frac{g_{a\gamma}}{10^{-12}{\rm
GeV^{-1}}}\right)^{2}$ $\displaystyle\times$
$\displaystyle\left(\frac{B_{\perp}}{10\,{\rm
G}}\right)^{2}\,\left(\frac{\mu{\rm
eV}}{m_{a}}\right)^{2}\,\left(\frac{0.5}{\alpha}\right)\,\left(\frac{h_{c}}{3\times
10^{3}{\rm km}}\right)\,e^{-\tau}\;.$

Figure 1: In the red region the electron cyclotron resonance causes a strong
absorption of the photons produced by ALPs conversion during their propagation
in the solar atmosphere. Assuming that the magnetic field decreases
monotonically with distance from the photosphere, this effect does not occur
when the magnetic field at the ALP conversion location lies in the white
region. The optical depth associated with thermal bremsstrahlung, the other
main absorption mechanism, is shown in blue (right-axis). The top axis shows
the electron density at which the plasma frequency matches the ALP mass
(bottom axis) and the resonant conversion can occur. Vertical lines mark
different regions of the solar atmosphere. The orange horizontal line marks
the typical value of 1 G of the magnetic field in the corona for the quiet
Sun.

In the previous equation we have considered a local DM density
$\rho_{\infty}=0.3\,\rm{GeVcm}^{-3}$ far away from the Sun and accounted for
the gravitational focusing of the Sun following [41]. This leads to a DM
density at the solar surface of $\rho_{a}\simeq 1\,\rm{GeVcm}^{-3}.$ Moreover,
taking a local DM velocity dispersion $v\simeq 10^{-3}$ leads to a intrinsic
bandwidth $\Delta\nu\simeq 10^{-6}.$ In the following, in our numerical
analysis performed to derive the sensitivity forecasts, the plasma frequency
and its gradient are computed from the solar electron density distribution of
Ref. [42]. We find that adopting a different model, e.g. [43], leads to a
similar sensitivity on the axion signal.

Let us now analyze the absorption term. There are two main contributions to
the opacity in the solar atmosphere at these wavelengths: the gyro-resonance
absorption and the thermal bremsstrahlung. The first effect occurs when the
resonant condition $m_{a}\simeq n\,\Omega_{B}(r)$ is met, where
$\Omega_{B}=e\,B/m_{e}\simeq 0.012\,\mu{\rm eV}\,B/{\rm G}$ is the electron
cyclotron frequency and the integer $n$ corresponds to the different energy
levels. To compute the location where the cyclotron resonance is realized, we
assume that the magnetic field decays with height following a power-law
behavior after the ALP conversion point. Then, we extract the electron number
density and the temperature of the plasma at the cyclotron resonance point.
Finally, with these ingredients and following Ref. [44], we obtain the
cyclotron optical depth. We find that for any reasonable choice of the
magnetic field profile in the solar atmosphere, the optical depth $\tau_{g}$
is very large for the first three transitions $n=1,2,3,$ while
$\tau_{g}\lesssim 1$ for $n=4$ and negligible for higher levels. Therefore,
assuming a monotonically decreasing magnetic field and plasma density, we
conservatively assume that photons from axion conversion can reach the Earth
only if the conversion happens after the surface of gyro-resonance absorption
associated with the $n=4$ level. This condition translates into an upper limit
of the magnetic field at the ALP conversion point, which is represented in
fig. 1, with the _no-signal region_ shown in red.

Concerning the thermal bremsstrahlung, the absorption coefficient reads [45]:

$\alpha_{b}\simeq 0.018\,\frac{T}{{\rm K}}\,\frac{n_{e}}{10^{6}\,{\rm
cm^{-3}}}\frac{Z^{2}_{eff}n_{i}}{10^{6}\,{\rm cm^{-3}}}\left(\frac{\nu}{{\rm
MHz}}\right)^{-2}\,G_{F}\;,$ (4)

where $n_{i}$ is the number density of ions with effective charge $Z_{eff}$,
$T$ is the temperature of the plasma, all taken from [42], and $G_{F}$ is the
Gaunt factor, for which we referred to [46]. In the corona, the plasma is
essentially fully ionized, and $Z_{eff}\approx 1$, while this is not the case
in the chromosphere and transition region. The bremsstrahlung optical depth
$\tau_{b}=\int_{h_{c}}\,d\ell\,\alpha_{b}$ is shown as a blue line in fig. 1.
In the chromosphere, we have $\tau_{b}\gg 1$ therefore the ALP signal is
completely absorbed. On the other hand, in the corona
$\tau_{b}\lesssim\mathcal{O}(1)$.

A crucial ingredient is the magnetic field in the corona, which determines the
conversion probability as described in eq. (1). The magnetic field of the
quiet Sun has a strength between $1-4$ G, on average [47]. We will take $B=1$
G when estimating the ALP signal from across the whole Sun. As mentioned
before, the magnetic field of sunspots can be very intense, even above the
photosphere. Sunspot observations accumulated over the years show great
variability of their properties, with sizes ranging from $10$ to $10^{5}$ km,
and magnetic fields from few G up to few thousand G [35].

Here, in the spirit of providing a general treatment, we will assume sunspots
to populate the white triangle in fig. 1, delimited by the red region and
orange line. Namely, we assume that there exist sunspots with associated
magnetic fields in the corona and transition region, where the resonance ALP
conversion occurs, larger than the magnetic field of the quiet Sun, and up to
the maximal value allowed by the cyclotron absorption, see fig. 1. For the
sake of concreteness, in the following, we will consider two cases: one where
the magnetic field maximizes the signal, i.e. lies at the border between the
red and white regions in fig. 1, and one where we simply assume $B=4$ G at the
ALP conversion radius above the sunspot.

To examine observational prospects, we derive forecasts for the first phase of
the SKA observatory, which is expected to be completed in its final
configuration in 2029. We follow [48], in particular their Table 6 provides
the line sensitivity for sources with an angular size within a certain range.
For larger extensions, we assume a linear degradation of the sensitivity with
the angular size of the source, as suggested by their Table 4 for our range of
interest.

The angular radius of a sunspot is typically smaller than a few arcmin, e.g.
$\ell_{s}=4\times 10^{4}$ km, which we will use for our estimates, corresponds
to $0.9^{\prime}$. On the other hand, radio photons are affected by scattering
processes during their propagation in the solar plasma, which leads to an
angular smearing of the signal. We estimate the broadening of a point-like
source as described in [49]. In the case of the entire Sun, this broadening is
instead marginal. The same conclusion has also been reported in [50], finding
that the observed radio size of the Sun at 100-240 MHz is only 25–30% larger
in area than the intrinsic one.

Finally, the line sensitivity reported in [48] refers to observations of
regions of the sky without strong radio emissions. Therefore, to include the
large radio background produced by the Sun, we correct the quoted sensitivity
adding the flux density of the Sun to the system-equivalent flux density of an
SKA station. Concretely, we follow the estimates provided in [51].

Figure 2: Projected sensitivity on the ALP-photon coupling $g_{a\gamma}$
versus ALP mass $m_{a}$, assuming 100 hours of observations with SKA-1 Low, of
a sunspot having the maximum magnetic field allowed by the gyro-resonance
absorption condition at the conversion location (see fig. 1) (solid), a
sunspot with a magnetic field of 4 G at the conversion location (dashed), and
the whole Sun (dotted). The radius of the conversion surface above a sunspot
is taken to be $\ell_{s}=4\times 10^{4}$ km. In red, we show the laboratory
bounds from ADMX [52] and CAST [53], while in green we show astrophysical
bounds: the R2-parameter in globular cluster [54], and the most conservative
limits from polarization of magnetic dwarf [55] and ALP production in pulsar
polar caps [56]. Other relevant constraints in this mass range (not shown in
the figure for clarity) come from the R-parameter in globular clusters [57],
X-ray observations of magnetic white dwarfs [58] and from $\gamma$-ray
observation of the blazar Markarian 421 (Mrk 421) [59, 60].

Following the procedure explained above, we obtain the results presented in
fig. 2. The projected sensitivity is at 95% C.L., namely, they correspond to
couplings providing a conversion signal equal to twice the SKA sensitivity.
For ALP masses in the range $10^{-6}-10^{-7}$ eV, the proposed search strategy
can improve existing laboratory bounds, and it is competitive with existing
astrophysical constraints. We find that the forecasted sensitivities for
observations of the entire Sun and a sunspot with a large magnetic field are
comparable. A technical advantage in the latter case is offered by the small
angular size of the emission region, in contrast to the case of the entire
Sun, where the large angular extension requires care in the data manipulation.
On the other hand, the properties of a sunspot needed to compute the
prediction of the expected flux might not be straightforward to be derived
from the observations.

The curves in fig. 2 are computed assuming an observing time of $\Delta t=100$
hr. Considering that solar physics is one of the key scientific goals of the
SKAO [61], hundreds of hours of observations of the Sun are foreseen. In the
case of sunspots, the $100$ hr value is guaranteed if intended as a
combination of different sunspots, whilst it is probably difficult to achieve
for a single one. Note that, in any case, the sensitivity to the ALP-photon
coupling scales very moderately with observing time, as $\Delta t^{1/4}$.

In fig. 2, we assume a bandwidth of $\Delta\nu/\nu=10^{-6}$, matching the one
induced by the dispersion of the DM velocity. This is achievable by the SKA
Low. On the other hand, this value might be somewhat optimistic, since
neglects possible smearing effects 222We postpone a full ray-tracing
computation to simulate photon propagation in the Sun to future work. Here we
accounted for absorption and spatial broadening as described in the text. The
bandwidth smearing is instead neglected for simplicity.. Note however that,
again, the sensitivities scale as $(\Delta\nu/\nu)^{1/4}$, so even
significantly larger bandwidths can still lead to interesting prospects.

To conclude, in this work, we highlighted that the Sun can be a promising
target to look for radio waves from ALP conversion. A discussion in this
direction was recently reported also in [38]. Here we presented a simple and
general derivation of the expected signal, with forecasts for the SKA, and
added the case of sunspots. The Sun can be thus considered as
complementary/alternative to the case of neutron stars, which are the target
that has been extensively investigated in connection to the ALP conversion
[18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34]. As a
simple comparison let us consider the case of a nearby isolated neutron star,
keeping in mind however that neutron stars could also be located in
environments hosting DM densities much larger than the local one, such as the
galactic center. Focusing on the main ingredients determining the radio flux,
we have

$\displaystyle S$ $\displaystyle\propto$
$\displaystyle\left(\frac{h_{c,s}\,\ell_{s}^{2}}{5\times 10^{12}\,{\rm
km^{3}}}\right)\left(\frac{B_{s}}{5\,{\rm G}}\right)^{2}\left(\frac{1.5\times
10^{8}\,{\rm km}}{d_{s}}\right)^{2}$ (5) $\displaystyle\simeq$
$\displaystyle\left(\frac{r_{c,{NS}}}{200\,{\rm
km}}\right)^{3}\left(\frac{B_{NS}}{10^{12}\,{\rm
G}}\right)^{2}\left(\frac{1\,{\rm kpc}}{d_{NS}}\right)^{2}\;,$

where $s$ denotes a large sunspot and $NS$ a neutron star. Therefore, naively,
we see that the expected signals are roughly comparable in order of magnitude.
However, there are important differences between the two systems, in terms of
systematic uncertainties and values of the relevant parameters. Hence, the two
approaches present complementary ways to explore the parameter space
corresponding to ALP dark matter. Since the properties of the Sun are much
better known and certainly more easily accessible to future measurements, what
is presented here could be a more robust and accessible method to probe the
ALP parameter space in the sub $\mu$eV region. Unfortunately, absorption
effects limit the range of masses that can be probed to $m_{a}\lesssim
10^{-6}$ eV, as discussed in the text. This also implies that low-frequency
radio observations will be the key to test the hypothesis discussed here.

## Acknowledgements

We would like to thank Gianni Bernardi for his help in determining the SKA
sensitivity, and Sam Witte for useful comments.

MT acknowledges the research grant “Addressing systematic uncertainties in
searches for dark matter No. 2022F2843L” funded by MIUR. ET, MR, and MT
acknowledge the project “Theoretical Astroparticle Physics (TAsP)” funded by
Istituto Nazionale di Fisica Nucleare (INFN). MR and ET acknowledge support
from the research grant ‘From Darklight to DM: understanding the galaxy/matter
connection to measure the Universe’ No. 20179P3PKJ funded by miur, and from
the “Grant for Internationalization” of the University of Torino. ET thanks
the University of Zaragoza for its hospitality during the initial phases of
this work. This article/publication is based upon work from COST Action COSMIC
WISPers CA21106, supported by COST (European Cooperation in Science and
Technology)”.

## References

  * Peccei and Quinn [1977a] R. D. Peccei and H. R. Quinn, CP Conservation in the Presence of Instantons, Phys. Rev. Lett. 38, 1440 (1977a), [,328(1977)].
  * Peccei and Quinn [1977b] R. D. Peccei and H. R. Quinn, Constraints Imposed by CP Conservation in the Presence of Instantons, Phys. Rev. D16, 1791 (1977b).
  * Weinberg [1978] S. Weinberg, A New Light Boson?, Phys. Rev. Lett. 40, 223 (1978).
  * Wilczek [1978] F. Wilczek, Problem of Strong p and t Invariance in the Presence of Instantons, Phys. Rev. Lett. 40, 279 (1978).
  * Kim and Carosi [2010] J. E. Kim and G. Carosi, Axions and the Strong CP Problem, Rev. Mod. Phys. 82, 557 (2010), [Erratum: Rev.Mod.Phys. 91, 049902 (2019)], arXiv:0807.3125 [hep-ph] .
  * Di Luzio  _et al._ [2020] L. Di Luzio, M. Giannotti, E. Nardi, and L. Visinelli, The landscape of QCD axion models, Phys. Rept. 870, 1 (2020), arXiv:2003.01100 [hep-ph] .
  * Preskill  _et al._ [1983] J. Preskill, M. B. Wise, and F. Wilczek, Cosmology of the Invisible Axion, Phys. Lett. B120, 127 (1983), [,URL(1982)].
  * Abbott and Sikivie [1983] L. F. Abbott and P. Sikivie, A Cosmological Bound on the Invisible Axion, Phys. Lett. B120, 133 (1983), [,URL(1982)].
  * Dine and Fischler [1983] M. Dine and W. Fischler, The Not So Harmless Axion, Phys. Lett. B120, 137 (1983), [,URL(1982)].
  * Chou  _et al._ [2022] A. S. Chou _et al._ , Snowmass Cosmic Frontier Report, in  _Snowmass 2021_ (2022) arXiv:2211.09978 [hep-ex] .
  * Adams  _et al._ [2022] C. B. Adams _et al._ , Axion Dark Matter, in  _Snowmass 2021_ (2022) arXiv:2203.14923 [hep-ex] .
  * Antel  _et al._ [2023] C. Antel _et al._ , Feebly Interacting Particles: FIPs 2022 workshop report, in  _Workshop on Feebly-Interacting Particles_ (2023) arXiv:2305.01715 [hep-ph] .
  * Arvanitaki  _et al._ [2010] A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper, and J. March-Russell, String Axiverse, Phys. Rev. D81, 123530 (2010), arXiv:0905.4720 [hep-th] .
  * Jaeckel and Ringwald [2010] J. Jaeckel and A. Ringwald, The Low-Energy Frontier of Particle Physics, Ann. Rev. Nucl. Part. Sci. 60, 405 (2010), arXiv:1002.0329 [hep-ph] .
  * Arias  _et al._ [2012] P. Arias, D. Cadamuro, M. Goodsell, J. Jaeckel, J. Redondo, and A. Ringwald, WISPy Cold Dark Matter, JCAP 06, 013, arXiv:1201.5902 [hep-ph] .
  * Sikivie [1983] P. Sikivie, Experimental Tests of the Invisible Axion,  _Particle physics and cosmology: Dark matter_ , Phys. Rev. Lett. 51, 1415 (1983), [,321(1983)].
  * Note [1] From now on, we will use the terms axion and ALP interchangeably.
  * Pshirkov and Popov [2009] M. S. Pshirkov and S. B. Popov, Conversion of Dark matter axions to photons in magnetospheres of neutron stars, J. Exp. Theor. Phys. 108, 384 (2009), arXiv:0711.1264 [astro-ph] .
  * Hook  _et al._ [2018] A. Hook, Y. Kahn, B. R. Safdi, and Z. Sun, Radio Signals from Axion Dark Matter Conversion in Neutron Star Magnetospheres, Phys. Rev. Lett. 121, 241102 (2018), arXiv:1804.03145 [hep-ph] .
  * Safdi  _et al._ [2019] B. R. Safdi, Z. Sun, and A. Y. Chen, Detecting Axion Dark Matter with Radio Lines from Neutron Star Populations, Phys. Rev. D 99, 123021 (2019), arXiv:1811.01020 [astro-ph.CO] .
  * Huang  _et al._ [2018] F. P. Huang, K. Kadota, T. Sekiguchi, and H. Tashiro, Radio telescope search for the resonant conversion of cold dark matter axions from the magnetized astrophysical sources, Phys. Rev. D97, 123001 (2018), arXiv:1803.08230 [hep-ph] .
  * Battye  _et al._ [2020] R. A. Battye, B. Garbrecht, J. I. McDonald, F. Pace, and S. Srinivasan, Dark matter axion detection in the radio/mm-waveband, Phys. Rev. D 102, 023504 (2020), arXiv:1910.11907 [astro-ph.CO] .
  * Leroy  _et al._ [2020] M. Leroy, M. Chianese, T. D. P. Edwards, and C. Weniger, Radio Signal of Axion-Photon Conversion in Neutron Stars: A Ray Tracing Analysis, Phys. Rev. D 101, 123003 (2020), arXiv:1912.08815 [hep-ph] .
  * Foster  _et al._ [2020] J. W. Foster, Y. Kahn, O. Macias, Z. Sun, R. P. Eatough, V. I. Kondratiev, W. M. Peters, C. Weniger, and B. R. Safdi, Green Bank and Effelsberg Radio Telescope Searches for Axion Dark Matter Conversion in Neutron Star Magnetospheres, Phys. Rev. Lett. 125, 171301 (2020), arXiv:2004.00011 [astro-ph.CO] .
  * Prabhu and Rapidis [2020] A. Prabhu and N. M. Rapidis, Resonant Conversion of Dark Matter Oscillons in Pulsar Magnetospheres, JCAP 10, 054, arXiv:2005.03700 [astro-ph.CO] .
  * Witte  _et al._ [2021] S. J. Witte, D. Noordhuis, T. D. P. Edwards, and C. Weniger, Axion-photon conversion in neutron star magnetospheres: The role of the plasma in the Goldreich-Julian model, Phys. Rev. D 104, 103030 (2021), arXiv:2104.07670 [hep-ph] .
  * Battye  _et al._ [2021] R. A. Battye, B. Garbrecht, J. I. McDonald, and S. Srinivasan, Radio line properties of axion dark matter conversion in neutron stars, JHEP 09, 105, arXiv:2104.08290 [hep-ph] .
  * Battye  _et al._ [2022] R. A. Battye, J. Darling, J. I. McDonald, and S. Srinivasan, Towards robust constraints on axion dark matter using PSR J1745-2900, Phys. Rev. D 105, L021305 (2022), arXiv:2107.01225 [astro-ph.CO] .
  * Wang  _et al._ [2021] J.-W. Wang, X.-J. Bi, and P.-F. Yin, Detecting axion dark matter through the radio signal from Omega Centauri, Phys. Rev. D 104, 103015 (2021), arXiv:2109.00877 [astro-ph.HE] .
  * Foster  _et al._ [2022] J. W. Foster, S. J. Witte, M. Lawson, T. Linden, V. Gajjar, C. Weniger, and B. R. Safdi, Extraterrestrial Axion Search with the Breakthrough Listen Galactic Center Survey, Phys. Rev. Lett. 129, 251102 (2022), arXiv:2202.08274 [astro-ph.CO] .
  * Zhou  _et al._ [2022] Y.-F. Zhou _et al._ , Searching for axion dark matter with the MeerKAT radio telescope, Phys. Rev. D 106, 083006 (2022), arXiv:2209.09695 [astro-ph.CO] .
  * McDonald and Witte [2023] J. I. McDonald and S. J. Witte, Generalized ray tracing for axions in astrophysical plasmas, Phys. Rev. D 108, 103021 (2023), arXiv:2309.08655 [hep-ph] .
  * McDonald  _et al._ [2023] J. I. McDonald, B. Garbrecht, and P. Millington, Axion-photon conversion in 3D media and astrophysical plasmas, JCAP 12, 031, arXiv:2307.11812 [hep-ph] .
  * Tjemsland  _et al._ [2023] J. Tjemsland, J. McDonald, and S. J. Witte, Adiabatic Axion-Photon Mixing Near Neutron Stars, (2023), arXiv:2310.18403 [hep-ph] .
  * Solanki [2003] S. K. Solanki, Sunspots: An overview, A&A Rev. 11, 153 (2003).
  * Anfinogentov  _et al._ [2019] S. A. Anfinogentov, A. G. Stupishin, I. I. Mysh’yakov, and G. D. Fleishman, Record-breaking Coronal Magnetic Field in Solar Active Region 12673, ApJL 880, L29 (2019), arXiv:1907.06398 [astro-ph.SR] .
  * An  _et al._ [2021] H. An, F. P. Huang, J. Liu, and W. Xue, Radio-frequency Dark Photon Dark Matter across the Sun, Phys. Rev. Lett. 126, 181102 (2021), arXiv:2010.15836 [hep-ph] .
  * An  _et al._ [2023] H. An, X. Chen, S. Ge, J. Liu, and Y. Luo, Searching for Ultralight Dark Matter Conversion in Solar Corona using LOFAR Data, arXiv e-prints , arXiv:2301.03622 (2023), arXiv:2301.03622 [hep-ph] .
  * Carlson and Tseng [1996] E. D. Carlson and L.-S. Tseng, Pseudoscalar conversion and x-rays from the sun, Phys. Lett. B 365, 193 (1996), arXiv:hep-ph/9507345 .
  * [40] E. Todarello _et al._ , Paper in preparation, .
  * Alenazi and Gondolo [2006] M. S. Alenazi and P. Gondolo, Phase-space distribution of unbound dark matter near the Sun, Phys. Rev. D 74, 083518 (2006), arXiv:astro-ph/0608390 .
  * Selhorst  _et al._ [2005] C. L. Selhorst, A. V. R. Silva, and J. E. R. Costa, Solar atmospheric model with spicules applied to radio observation, A&A 433, 365 (2005).
  * De La Luz  _et al._ [2008] V. De La Luz, A. Lara, E. Mendoza, and M. Shimojo, 3D Simulations of the Quiet Sun Radio Emission at Millimeter and Submillimeter Wavelengths, Geofisica Internacional 47, 197 (2008).
  * Dulk [1985] G. A. Dulk, Radio emission from the sun and stars., ARA&A 23, 169 (1985).
  * Rybicki and Lightman [1986] G. B. Rybicki and A. P. Lightman,  _Radiative Processes in Astrophysics_ (1986).
  * Chluba  _et al._ [2020] J. Chluba, A. Ravenni, and B. Bolliet, Improved calculations of electron-ion bremsstrahlung Gaunt factors for astrophysical applications, MNRAS 492, 177 (2020), arXiv:1911.08861 [astro-ph.CO] .
  * Yang  _et al._ [2020] Z. Yang, C. Bethge, H. Tian, S. Tomczyk, R. Morton, G. Del Zanna, S. W. McIntosh, B. B. Karak, S. Gibson, T. Samanta, J. He, Y. Chen, and L. Wang, Global maps of the magnetic field in the solar corona, Science 369, 694 (2020), arXiv:2008.03136 [astro-ph.SR] .
  * Braun  _et al._ [2019] R. Braun, A. Bonaldi, T. Bourke, E. Keane, and J. Wagg, Anticipated Performance of the Square Kilometre Array – Phase 1 (SKA1), arXiv e-prints , arXiv:1912.12699 (2019), arXiv:1912.12699 [astro-ph.IM] .
  * Kontar  _et al._ [2019] E. P. Kontar, X. Chen, N. Chrysaphi, N. L. S. Jeffrey, A. G. Emslie, V. Krupar, M. Maksimovic, M. Gordovskyy, and P. K. Browning, Anisotropic Radio-wave Scattering and the Interpretation of Solar Radio Emission Observations, ApJ 884, 122 (2019), arXiv:1909.00340 [astro-ph.SR] .
  * Sharma and Oberoi [2020] R. Sharma and D. Oberoi, Propagation Effects in Quiet Sun Observations at Meter Wavelengths, ApJ 903, 126 (2020), arXiv:2009.10604 [astro-ph.SR] .
  * Macario  _et al._ [2022] G. Macario, G. Pupillo, G. Bernardi, P. Bolli, P. Di Ninni, G. Comoretto, A. Mattana, J. Monari, F. Perini, M. Schiaffino, M. Sokolowski, R. Wayth, J. Broderick, M. Waterson, M. Grazia Labate, R. Chiello, A. Magro, T. Booler, A. Mcphail, D. Minchin, and R. Bhushan, Characterization of the SKA1-Low prototype station Aperture Array Verification System 2, Journal of Astronomical Telescopes, Instruments, and Systems 8, 011014 (2022), arXiv:2109.11983 [astro-ph.IM] .
  * Asztalos  _et al._ [2010] S. J. Asztalos, G. Carosi, C. Hagmann, D. Kinion, K. van Bibber, M. Hotz, L. J. Rosenberg, G. Rybka, J. Hoskins, J. Hwang, P. Sikivie, D. B. Tanner, R. Bradley, and J. Clarke, Squid-based microwave cavity search for dark-matter axions, Physical Review Letters 104, 10.1103/physrevlett.104.041301 (2010).
  * CAS [2017] New cast limit on the axion–photon interaction, Nature Physics 13, 584–590 (2017).
  * Dolan  _et al._ [2022] M. J. Dolan, F. J. Hiskens, and R. R. Volkas, Advancing globular cluster constraints on the axion-photon coupling, Journal of Cosmology and Astroparticle Physics 2022 (10), 096.
  * Dessert  _et al._ [2022a] C. Dessert, D. Dunsky, and B. R. Safdi, Upper limit on the axion-photon coupling from magnetic white dwarf polarization, Physical Review D 105, 10.1103/physrevd.105.103034 (2022a).
  * Noordhuis  _et al._ [2023] D. Noordhuis, A. Prabhu, S. J. Witte, A. Y. Chen, F. Cruz, and C. Weniger, Novel constraints on axions produced in pulsar polar-cap cascades (2023), arXiv:2209.09917 [hep-ph] .
  * Ayala  _et al._ [2014] A. Ayala, I. Domínguez, M. Giannotti, A. Mirizzi, and O. Straniero, Revisiting the bound on axion-photon coupling from globular clusters, Physical Review Letters 113, 10.1103/physrevlett.113.191302 (2014).
  * Dessert  _et al._ [2022b] C. Dessert, A. J. Long, and B. R. Safdi, No evidence for axions from chandra observation of the magnetic white dwarf re j0317-853, Physical Review Letters 128, 10.1103/physrevlett.128.071102 (2022b).
  * Li  _et al._ [2021] H.-J. Li, J.-G. Guo, X.-J. Bi, S.-J. Lin, and P.-F. Yin, Limits on axionlike particles from mrk 421 with 4.5-year period observations by argo-ybj and fermi-lat, Physical Review D 103, 10.1103/physrevd.103.083003 (2021).
  * Li  _et al._ [2022] H.-J. Li, X.-J. Bi, and P.-F. Yin, Searching for axion-like particles with the blazar observations of magic and fermi-lat *, Chinese Physics C 46, 085105 (2022).
  * Nindos  _et al._ [2019] A. Nindos, E. P. Kontar, and D. Oberoi, Solar physics with the Square Kilometre Array, Advances in Space Research 63, 1404 (2019), arXiv:1810.04951 [astro-ph.SR] .
  * Note [2] We postpone a full ray-tracing computation to simulate photon propagation in the Sun to future work. Here we accounted for absorption and spatial broadening as described in the text. The bandwidth smearing is instead neglected for simplicity.


# Reward Once, Penalize Once: Rectifying Time Series Anomaly Detection

Keval Doshi  Department of Electrical Engineering)
University of South Florida)
Tampa, USA
<EMAIL_ADDRESS>Shatha Abudalou  Department of Electrical Engineering)
University of South Florida)
Tampa, USA
<EMAIL_ADDRESS>Yasin Yilmaz  Department of Electrical Engineering)
University of South Florida)
Tampa, USA
<EMAIL_ADDRESS>
###### Abstract

While anomaly detection in time series has been an active area of research for
several years, most recent approaches employ an inadequate evaluation
criterion leading to an inflated F1 score. We show that a rudimentary Random
Guess method can outperform state-of-the-art detectors in terms of this
popular but faulty evaluation criterion. In this work, we propose a proper
evaluation metric that measures the timeliness and precision of detecting
sequential anomalies. Moreover, most existing approaches are unable to capture
temporal features from long sequences. Self-attention based approaches, such
as transformers, have been demonstrated to be particularly efficient in
capturing long-range dependencies while being computationally efficient during
training and inference. We also propose an efficient transformer approach for
anomaly detection in time series and extensively evaluate our proposed
approach on several popular benchmark datasets.

###### Index Terms:

Time series, anomaly detection, sequential anomalies, self-attention,
transformer

##  I Introduction

Time series analysis is used to perform important tasks such as predicting the
future values of a variable (e.g., stock market price) and detecting anomalies
in sequential data. Time series anomaly detection methods aim to identify
abnormal data patterns in temporal data. For instance, in health care, ECG
signals are analyzed to determine if the patient suffers from a heart disease
[1]. Similarly, in cybersecurity, the data traffic over time in a computer
network is monitored to detect cyber-attacks [2]. Time series anomaly
detection methods are used in various domains by companies, e.g., Extensible
Generic Anomaly Detection System (EGADS) by Yahoo [3] and SR/CNN developed by
Microsoft [4].

Time series anomalies are typically classified into two main categories, point
anomalies and sequential anomalies. A point anomaly, also known as outlier, is
a single data instance with unexpected value with respect to the nominal
baseline. In many applications, anomalies continue for a duration with
successive anomalous data instances, which is called a sequential anomaly. [5]
presents a detailed survey about time series anomaly detection by evaluating
twenty different methods based on statistical and deep learning approaches on
univariate time series datasets.

Figure 1: The commonly used adjusted instance-based evaluation method. A threshold of $0.5$ is applied to an example sequence of anomaly scores produced by a detection algorithm. The traditional instance-based evaluation compares the anomaly/no-anomaly decision for each instance with the ground truth to determine true/false positive/negative decisions. In the recently proposed and widely used adjusted instance-based evaluation, while errors are penalized once as in the traditional evaluation approach, true detections are greatly amplified by considering all instances in an anomalous sequence as multiple true positives if an alarm is raised during the anomalous sequence. This amplification of true positives causes an artificially inflated F1 score (see Table I). Dataset | SMD | MSL | SMAP | SWaT | PSM
---|---|---|---|---|---
Metric | P | R | F1 | P | R | F1 | P | R | F1 | P | R | F1 | P | R | F1
OCSVM [6] | 44.34 | 76.72 | 56.19 | 59.78 | 86.87 | 70.82 | 53.85 | 59.07 | 56.34 | 45.39 | 49.22 | 47.23 | 62.75 | 80.89 | 70.67
IsolationForest [7] | 42.31 | 73.29 | 53.64 | 53.94 | 86.54 | 66.45 | 52.39 | 59.07 | 55.53 | 49.29 | 44.95 | 47.02 | 76.09 | 92.45 | 83.48
LOF [8] | 56.34 | 39.86 | 46.68 | 47.72 | 85.25 | 61.18 | 58.93 | 56.33 | 57.60 | 72.15 | 65.43 | 68.62 | 57.89 | 90.49 | 70.61
Deep-SVDD [9] | 78.54 | 79.67 | 79.10 | 91.92 | 76.63 | 83.58 | 89.93 | 56.02 | 69.04 | 80.42 | 84.45 | 82.39 | 95.41 | 86.49 | 90.73
DAGMM [10] | 67.30 | 49.89 | 57.30 | 89.60 | 63.93 | 74.62 | 86.45 | 56.73 | 68.51 | 89.92 | 57.84 | 70.40 | 93.49 | 70.03 | 80.08
LSTM-VAE [11] | 75.76 | 90.08 | 82.30 | 85.49 | 79.94 | 82.62 | 92.20 | 67.75 | 78.10 | 76.00 | 89.50 | 82.20 | 73.62 | 89.92 | 80.96
BeatGAN [12] | 72.90 | 84.09 | 78.10 | 89.75 | 85.42 | 87.53 | 92.38 | 55.85 | 69.61 | 64.01 | 87.46 | 73.92 | 90.30 | 93.84 | 92.04
OmniAnomaly [13] | 83.68 | 86.82 | 85.22 | 89.02 | 86.37 | 87.67 | 92.49 | 81.99 | 86.92 | 81.42 | 84.30 | 82.83 | 88.39 | 74.46 | 80.83
InterFusion [14] | 87.02 | 85.43 | 86.22 | 81.28 | 92.70 | 86.62 | 89.77 | 88.52 | 89.14 | 80.59 | 85.58 | 83.01 | 83.61 | 83.45 | 83.52
THOC [15] | 79.76 | 90.95 | 84.99 | 88.45 | 90.97 | 89.69 | 92.06 | 89.34 | 90.68 | 83.94 | 86.36 | 85.13 | 88.14 | 90.99 | 89.54
Anomaly Transformer [16] | 89.40 | 95.45 | 92.33 | 92.09 | 95.15 | 93.59 | 94.13 | 99.40 | 96.69 | 91.55 | 96.73 | 94.07 | 96.91 | 98.90 | 97.89
_Random Guess_ | 90.54 | 96.43 | 93.39 | 94.95 | 98.13 | 96.51 | 95.43 | 97.74 | 96.48 | 93.21 | 99.47 | 96.24 | 98.83 | 98.14 | 98.48
TABLE I: Performance of _Random Guess_ with probability of alarm $p=0.01$
dominates the state-of-the-art methods in terms of the commonly used adjusted
instance-based evaluation, demonstrating its inherent flaw.

Predictive models provide an intuitive way to detect anomalies [17]. The
motivation is that a predictive model trained on nominal data should give
statistically similar prediction error for nominal test data, whereas the
prediction error is expected to be larger when encountered with anomalous
data. Recurrent Neural Networks (RNN) replaced the classical statistical
methods such as Autoregressive (AR) and Autoregressive Moving Average (ARMA)
models in many applications. While the early RNN structures could not utilize
long-term dependencies in time series data due to the diminishing gradient
problem, Long Short-Term Memory (LSTM) network overcome this problem by
introducing a more complex memory unit [18]. Recently, the attention mechanism
was proposed to improve the predictive performance of LSTM [19]. However, more
recently, a completely new deep neural network structure called self-
attention, also known as transformer, significantly outperformed the
LSTM+attention model to become the state-of-the-art predictive model in many
applications. While both the attention and self-attention mechanisms are
originally proposed for Natural Language Processing (NLP), they were shown to
be effective in various other time series data domains [20, 21].

Evaluating performance for detecting sequential anomalies has been
traditionally done in the same way as point anomalies using instance-based
detection metrics such as AUC and F1 score. However, consecutive anomalous
instances are typically caused by the same event, and in real-world
applications raising an alarm for such anomalous event is what matters. For
instance, after successfully raising alarm for an anomalous event, it is not
important in practice to label every anomalous instance. Hence, recent state-
of-the-art methods in time series anomaly detection [15, 14, 13, 22, 23, 24,
16] used a different performance evaluation system instead of the traditional
instance-based evaluation, which is suitable for point anomalies. In this new
evaluation system, all instances in an anomalous segment are considered as
true positives if a single alarm is raised in the entire segment, as shown in
Fig. 1. While this evaluation system rightfully focuses on the detection of
anomalous events/sequences, it fails to provide a meaningful metric for
sequential anomalies since it still uses the instance-based F1 score for
performance evaluation. In this _adjusted instance-based evaluation_ , since
errors are penalized once, but detection is rewarded generously, it leads to
an inflated F1 metric.

To illustrate the serious flaw of the adjusted instance-based F1 metric, let
us consider using a _Random Guess_ method which randomly (e.g., from a uniform
distribution) raises an alarm for each instance with probability $p$. Note
that this _Random Guess_ method makes arbitrary decisions independent of the
dataset. The probability of raising a true alarm increases with the duration
of anomalous sequence. The exact probabilities are given in Section III. Since
its each random detection is amplified by the anomalous sequence duration
while its errors are penalized once, such a rudimentary approach is able to
achieve very high adjusted instance-based F1 score on the popular benchmark
datasets and even outperform the state-of-the-art models, as shown in Table I.

Motivated by this limitation, we propose a proper performance metric for
sequential anomalies, which also evaluates the timeliness of alarm. Leveraging
the high potential of transformer architectures in capturing the long-term
dependencies in time series data, we also propose a novel transformer-based
time series anomaly detector. Our contributions in this paper can be
summarized as follows:

  * •

A thorough analysis of the inherent flaw in the existing evaluation metric and
a proper metric that measures the timeliness and precision of raised alarms.

  * •

A novel end-to-end trained transformer algorithm for time series anomaly
detection with asymptotic false alarm rate analysis and closed-form expression
for detection threshold.

After discussing related works in Section II and limitations of the existing
evaluation system in Section III, we present the proposed performance metric
and the detector in Section IV and the experimental results in Section V. The
paper is concluded in Section VI.

##  II Related Work

The detection of anomalies in time series has been extensively investigated
for a long time and remains an active subject due to the great need for more
robust methods in complex real-world scenarios [25, 26]. Popular approaches
for detecting anomalies in time series data include CNN models [27, 28, 29],
RNN models [30, 31, 32, 33], and spectral residuals [4, 34]. The most recent
approaches use the attention-based mechanisms for time series forecasting [35,
36, 37].

Self-attention, which is also known as transformer, is an attention-based
neural network architecture without the sequential structure. Current state-
of-the results in the NLP domain are obtained by transformer models, e.g.,
GPT-3 [38]. The transformer encoder in [39] is used for unsupervised
representation learning of multivariate time series data. It outperforms the
state-of-the-art time series classification and regression methods.

The Temporal Fusion Transformer (TFT) presented in [40] is an attention-based
model used for multi-horizon forecasting. TFT employs recurrent layers for
local processing and interpretable self-attention layers for long-term
dependency to learn temporal correlations at various scales. In addition, TFT
curbs any pointless components by using specialized elements that pick the
critical characteristics and a succession of gating layers to get significant
performance in various applications.

Based on Generative Adversarial Network (GAN), in [41], the authors develop
the Adversarial Sparse Transformer, which acts as a generator for learning
sparse attention mappings for specific time steps to enhance time series
forecasting. Using graph learning with the transformer-based network [42]
proposes learning the graph structure for IoT systems to learn sensor
dependencies in multivariate time series datasets automatically. Informer is
another time series forecasting model based on self-attention, which can be
used in anomaly detection [43, 44].

##  III Problem Formulation

Figure 2: Proposed TiSAT architecture

Sequential Anomaly Detection: Consider a time series
$\\{X_{1},\ldots,X_{t},\ldots\\}$, which may include anomalous sequences
starting and ending at unknown times. Denote the unknown starting time of the
$i$th anomalous sequence with $\tau_{i}$. Since in real-world applications,
such as cybersecurity, surveillance, etc. such anomalous sequences are caused
by potentially hazardous anomalous events, it is critical to detect such
sequences in a timely manner. Controlling the number of false alarms is also
crucial to ensure the reliability of the detection system. Instead of the
traditional instance-based evaluation of performance, which is commonly used
for point anomalies and other standard machine learning tasks (e.g.,
classification, regression), the performance of sequential anomaly detection
methods should be evaluated in terms of true/false detection of anomalous
sequences.

Flaw of Adjusted Instance-based Evaluation: The adjusted instance-based
evaluation method has been extensively used in the recent literature (e.g.,
[15, 14, 13, 22, 23, 24, 16]) to compare the state-of-the-art deep learning
methods. However, this evaluation method is severely flawed, as illustrated by
the high performance of the _Random Guess_ algorithm (Table I). Assuming
_Random Guess_ raises an alarm with probability $p$, the expected number of
false alarms is equal to $Np$, where $N$ is the number of nominal instances.
The probability of raising a true alarm within the $i$th anomaly sequence of
length $M_{i}$ is given by $1-Binom(M_{i},0,p)$, where
$Binom(M_{i},0,p)=(1-p)^{M_{i}}$ is the binomial probability mass function for
zero success with $M_{i}$ trials and $p$ success probability. With the
adjusted labeling of the entire sequence as true positive in case of a true
alarm, the expected number of true positives is equal to
$\sum_{i=1}M_{i}[1-(1-p)^{M_{i}}]$. Consequently, the expected number of false
negatives is given by $\sum_{i=1}M_{i}(1-p)^{M_{i}}$. Hence, the expected
precision and recall are

$\displaystyle\begin{split}\text{Precision}_{\text{RandomGuess}}&=\frac{\sum_{i=1}M_{i}[1-(1-p)^{M_{i}}]}{\sum_{i=1}M_{i}[1-(1-p)^{M_{i}}]+Np}\\\
\text{Recall}_{\text{RandomGuess}}&=\frac{\sum_{i=1}M_{i}[1-(1-p)^{M_{i}}]}{\sum_{i=1}M_{i}}.\end{split}$
(1)

As the duration $M_{i}$ of anomalous sequences increase, the expected number
of true positives increases significantly, making both precision and recall
approach to $1$. Note the insignificant effect of false positives ($Np$) since
they are penalized once. In the popular benchmark datasets, there are long
anomalous sequences lasting thousands of instances, and thus we observe the
high precision, recall, and F1 scores in Table I with $p=0.01$.

##  IV Proposed Approach

In this section, we first present a proper performance metric for sequential
anomalies and then our transformer based anomaly detection approach.

###  IV-A Performance Evaluation

Sequence Detection Delay: Given $\tau_{i}$ as the starting time of an
anomalous event $i$ and $T_{i}\geq\tau_{i}$ as the alarm time, we can
empirically formulate the average detection delay as

$\text{ADD}=\frac{1}{S}\sum_{i=1}^{S}(T_{i}-\tau_{i}),$ (2)

where $S$ denotes the number of anomalous events. Since most anomalies
indicate critical incidents, it might be essential to detect an anomalous
event within a certain time period. Hence, if no alarm is raised within the
duration $[\tau_{i},\tau_{i}+\delta_{\text{max}}]$ after anomalous activity
$i$ happens, we set the delay to the maximum tolerable delay
$\delta_{\text{max}}$. Here, it is important to note that minimizing the
detection delay is analogous to the the more commonly used objective of
maximizing the true positive rate, except it assigns a more specific cost of
detection delay $\delta_{i}=T_{i}-\tau_{i}$.

Sequence Alarm Precision: Our second objective emphasizes on maximizing the
number of anomalous events being detected with respect to the total number of
alarms, similar to the well-known precision metric. However, in contrast to
the instance-based precision metric, our metric focuses on the detection of
true anomalous _sequences_ , and hence only focuses on detecting the anomalous
event onset accurately. If an alarm is raised before an alarm even begins,
i.e., $T^{j}\leq\tau_{i}$, then it is considered as a false alarm.

Empirically, the alarm precision is computed as

$P=\frac{1}{\hat{S}}\sum_{j=1}^{\hat{S}}\mathbbm{1}_{\\{=\\}}{T_{j}\in\cup[\tau_{i},\tau_{i}+\delta_{\text{max}}]},$
(3)

where $\mathbbm{1}_{\\{\cdot\\}}$ denotes the indicator function,
$\hat{S}=|\\{T_{j}\\}|$ is the number of all alarms, and $|\cdot|$ denotes the
cardinality of a set.

Sequence Precision Delay: Finally, we present a new metric called _Sequence
Precision Delay_ that combines the sequence based average detection delay with
sequence alarm precision in order to achieve a single metric for easily
comparing time series algorithms. The SPD statistic quantifies the area under
the Precision vs. normalised ADD (NADD) curve, much like the common AUC metric
does for TPR and FPR. To map ADD into $[0,1]$, we normalize it by the maximum
delay, i.e., $\text{NADD}=\text{ADD}/\delta_{\text{max}}$. Mathematically, SPD
is given by

$\text{SPD}=\int_{0}^{1}P(\alpha)~{}\text{d}\alpha,$ (4)

where $\alpha$ denotes NADD, and $P$ denotes the precision. A highly
successful algorithm with an SPD value close to $1$ must have high precision
and low delay in its alarms.

Most existing approaches leverage an RNN based model for time series
forecasting, and compute the _residuals_ , i.e. the prediction or
reconstruction error to determine if an observation is anomalous or not.
However, it is shown in [43] that RNN based approaches suffer in long sequence
time series forecasting. To this end, we propose a novel transformer based
approach called _Time Series Anomaly Transformer (TiSAT)_ , which is superior
in capturing long range temporal dependencies. We next discuss our proposed
approach in detail.

###  IV-B Time Series Anomaly Transformer (TiSAT)

The overall structure of TiSAT is shown in Fig. 2. The proposed approach
utilizes the the ProbSparse mechanism discussed in the Informer architecture
[43] as compared to the self-attention mechanism proposed in the Vanilla
Transformer [45] for reducing computational complexity. Traditionally, the
self attention mechanism for a $d$-dimensional input is defined as

$A(Q,K,V)=Softmax\left(\frac{QK^{T}}{\sqrt{d}}\right)V$ (5)

where $Q$, $K$ and $V$ represent the query, key and value respectively.
However, it was recently observed that only a few key and value pairs
contribute to the attention score, rendering a majority of the computed dot
products as worthless. Hence, we propose using a probabilistic attention
mechanism, since it reduces the computational complexity from
$\mathcal{O}(L^{2})$ to $\mathcal{O}(LlogL)$. Particularly, since only a
subset of the query/value tensors require costly operations, ProbSparse
attention allows each key to focus on the most important queries rather than
all of them. The ProbSparse attention is given by

$A(Q,K,V)=Softmax\left(\frac{\bar{Q}K^{T}}{\sqrt{d}}\right)V$ (6)

where $\bar{Q}$ is a sparse matrix consisting of the top queries. As shown in
[43], the measurement of sparsity between a query $q_{i}$ and all keys $K$ is
given by

$M(q_{i},K)=\log\left(\sum_{j=1}^{L_{K}}e^{\frac{q_{i}k_{j}^{T}}{\sqrt{d}}}\right)-\frac{1}{L_{K}}\sum_{j=1}^{L_{K}}\frac{q_{i}k_{j}^{T}}{\sqrt{d}}$
(7)

Encoder: The input to the proposed architecture at a time instance $t$ is a
matrix representation of the $d$-dimensional time series over a sequence of
$L$ instances given by $\mathcal{Z}^{t}$. Following the dilated convolution
approach proposed in [46, 47], we employ a distilling procedure to extract a
focused self-attention feature map. The distilling function from the $j^{th}$
layer towards the $(j+1)^{th}$ layer is given by

$\mathcal{Z}_{j+1}^{t}=\text{MaxPool ELU}(Conv1d([Z_{j}]))$ (8)

where $\left[\cdot\right]$ represents the attention block, which is followed
by a 1-D convolution filter with a kernel width of 3. The number of self-
attention blocks are progressively decreased and then finally concatenated to
form the final representation of the encoder.

Decoder: We leverage the canonical decoder structure proposed by [45] and
consists of two similar self-attention modules. The input to the decoder
architecture is given by concatenating the start token ($Z_{token})$ and a
placeholder for the target sequence ($Z_{target}$) as follows

$Z^{t}_{de}=Concat(Z_{token},Z_{target})$ (9)

This is then passed to a dense fully connected layer. The network is trained
using the mean squared error loss by propagating it through the decoder and
encoder.

###  IV-C Anomaly Detection Framework

We propose an online and non-parametric detection approach for detecting
persistent and abrupt anomalies using the transformer output. Due to the
sequential (persistent) nature of time series anomalies, we need an approach
which accumulates the evidence over time and then makes a decision instead of
a hard threshold on the anomaly score for each instant. To this end, we
propose using a nonparametric sequential algorithm based on $k$ nearest
neighbors ($k$NN). First, the algorithm trains on a set of nominal historic
observations in an offline fashion and then tests the incoming observations
until it detects a change in the observations with respect to the nominal
baseline. In the training phase, assuming a training set $\mathcal{X}_{N}$
consisting of $N$ nominal data instances, it randomly partitions
$\mathcal{X}_{N}$ into two sets $\mathcal{X}_{N_{1}}$ and
$\mathcal{X}_{N_{2}}$, where $N_{1}+N_{2}=N$. Then, for each point in
$\mathcal{X}_{N_{1}}$ it finds the Euclidean distance to each point in
$\mathcal{X}_{N_{2}}$ as $\\{d_{1},d_{2},\dots,d_{N_{1}}\\}$. The
(1-$\alpha$)th percentile $d_{\alpha}$ is used as a baseline statistic during
the testing phase, where $\alpha$ is the statistical significance level (e.g.,
$0.05$). During testing, we compute the anomaly evidence for each instance as

$D_{t}=d_{t}^{m}-d_{\alpha}^{m},$ (10)

where $d_{t}$ is the $k$NN distance between the test point at time instance
$t$ and $\mathcal{X}_{N_{2}}$ and $m$ is the dimensionality of the transformer
output.

The anomaly evidences are accumulated over time

$s_{t}=\max\\{s_{t}+D_{t},0\\},s_{0}=0.$ (11)

and an alarm is raised when the anomaly statistic $s_{t}$ exceeds a threshold
$h$, i.e., at time

$T=\min\\{t:s_{t}^{k}\leq h\\}.$ (12)

Here, $T$ is the minimum time required for the accumulated evidence
$s_{t}^{k}$ to be sufficiently high to raise an alarm based on a detection
threshold $h$. The detection threshold $h$ manifests a trade-off between
minimizing the detection delay and minimizing the false alarm.

For an anomaly detection algorithm to be implemented in a practical setting, a
clear procedure is necessary for selecting the decision threshold such that it
satisfies a desired false alarm rate. The reliability of an algorithm in terms
of false alarm rate is crucial for minimizing human involvement. To provide
such a performance guarantee for the false alarm rate, we derive an asymptotic
upper bound on the average false alarm rate of the proposed algorithm.

###### Theorem 1

The false alarm rate of the proposed algorithm is asymptotically (as
$N_{2}\rightarrow\infty$) upper bounded by

$FAR\leq e^{-\omega_{0}h},$ (13)

where $h$ is the decision threshold, and $\omega_{0}>0$ is given by

$\displaystyle\omega_{0}$
$\displaystyle=v_{m}-\theta-\frac{1}{\phi}\mathcal{W}\left(-\phi\theta
e^{-\phi\theta}\right),$ (14) $\displaystyle\theta$
$\displaystyle=\frac{v_{m}}{e^{v_{m}d_{\alpha}^{m}}}.$

In (14), $\mathcal{W}(\cdot)$ is the Lambert-W function,
$v_{m}=\frac{\pi^{m/2}}{\Gamma(m/2+1)}$ is the constant for the
$m$-dimensional Lebesgue measure (i.e., $v_{m}d_{\alpha}^{m}$ is the
$m$-dimensional volume of the hyperball with radius $d_{\alpha}$), and $\phi$
is the upper bound for $D_{t}$.

Proof. In [48][page 177], for CUSUM-like algorithms with independent
increments, such as TiSAT with independent $D_{t}$, a lower bound on the
average false alarm period is given as follows

$E_{\infty}[T]\geq e^{\omega_{0}h},$

where $h$ is the detection threshold, and $\omega_{0}\geq 0$ is the solution
to $E[e^{\omega_{0}D_{t}}]=1$.

To analyze the false alarm period, we need to consider the nominal case. The
anomaly evidence in the nominal case does not necessarily depend on the
previous selections due to lack of an anomaly which could correlate the
evidences. Hence, in the nominal case, it is safe to assume that $D_{t}$ is
independent over time.

We firstly derive the asymptotic distribution of the instance-level anomaly
evidence $D_{t}$ in the absence of anomalies. Its cumulative distribution
function is given by

$P(D_{t}\leq y)=P(d_{t}^{m}\leq d_{\alpha}^{m}+y).$

It is sufficient to find the probability distribution of $d_{t}^{m}$, the
$m$th power of the $k$NN distance at time $t$. As discussed above, we have
independent $m$-dimensional vectors $\\{X_{t}\\}$ over time, which form a
Poisson point process. The nearest neighbor ($k=1$) distribution for a Poisson
point process is given by

$P(d_{t}\\}\leq r)=1-\exp(-\Lambda(b(X_{t},r)))$

where $\Lambda(b(X_{t},r))$ is the arrival intensity (i.e., Poisson rate
measure) in the $m$-dimensional hypersphere $b(X_{t},r)$ centered at $X_{t}$
with radius $r$ [49]. Asymptotically, for a large number of training instances
as $N_{2}\to\infty$, under the null (nominal) hypothesis, $d_{t}$ of $X_{t}$
takes small values, defining an infinitesimal hyperball with homogeneous
intensity $\lambda=1$ around $X_{t}$. Since for a homogeneous Poisson process
the intensity is written as $\Lambda(b(X_{t},r))=\lambda|b(X_{t},r)|$ [49],
where $|b(X_{t},r)|=\frac{\pi^{m/2}}{\Gamma(m/2+1)}r^{m}=v_{m}r^{m}$ is the
Lebesgue measure (i.e., $m$-dimensional volume) of the hyperball $b(X_{t},r)$,
we rewrite the nearest neighbor distribution as

$\displaystyle P(d_{t}\leq r)$ $\displaystyle=1-\exp\left(-v_{m}r^{m}\right),$

where $v_{m}=\frac{\pi^{m/2}}{\Gamma(m/2+1)}$ is the constant for the
$m$-dimensional Lebesgue measure. Now, applying a change of variables we can
write the probability density of $d_{t}^{m}$ and $D_{t}$ as

$\displaystyle f_{d_{t}^{m}}(y)$ $\displaystyle=\frac{\partial}{\partial
y}\left[1-\exp\left(-v_{m}y\right)\right],$ (15)
$\displaystyle=v_{m}\exp(-v_{m}y),$ (16) $\displaystyle f_{D_{t}}(y)$
$\displaystyle=v_{m}\exp(-v_{m}d_{\alpha}^{m})\exp(-v_{m}y)$ (17)

Using the probability density derived in (17), $E[e^{\omega_{0}D_{t}}]=1$ can
be written as

$\displaystyle 1$
$\displaystyle=\int_{-d_{\alpha}^{m}}^{\phi}e^{\omega_{0}y}v_{m}e^{-v_{m}d_{\alpha}^{m}}e^{-v_{m}y}dy,$
(18) $\displaystyle\frac{e^{v_{m}d_{\alpha}^{m}}}{v_{m}}$
$\displaystyle=\int_{-d_{\alpha}^{m}}^{\phi}e^{(\omega_{0}-v_{m})y}dy,$ (19)
$\displaystyle=\frac{e^{(\omega_{0}-v_{m})y}}{\omega_{0}-v_{m}}\Biggr{|}_{-d_{\alpha}^{m}}^{\phi},$
(20)
$\displaystyle=\frac{e^{(\omega_{0}-v_{m})\phi}-e^{(\omega_{0}-v_{m})(-d_{\alpha}^{m})}}{\omega_{0}-v_{m}},$
(21)

where $-d_{\alpha}^{m}$ and $\phi$ are the lower and upper bounds for
$D_{t}=d_{t}^{m}-d_{\alpha}^{m}$. The upper bound $\phi$ is obtained from the
training set.

As $N_{2}\to\infty$, since the $m$th power of $(1-\alpha)$th percentile of
nearest neighbor distances in training set goes to zero, i.e.,
$d_{\alpha}^{m}\to 0$, we have

$\displaystyle e^{(\omega_{0}-v_{m})\phi}$
$\displaystyle=\frac{e^{v_{m}d_{\alpha}^{m}}}{v_{m}}(\omega_{0}-v_{m})+1.$
(22)

We next rearrange the terms to obtain the form of $e^{\phi x}=a_{0}(x+\theta)$
where $x=\omega_{0}-v_{m}$, $a_{0}=\frac{e^{v_{m}d_{\alpha}^{m}}}{v_{m}}$, and
$\theta=\frac{v_{m}}{e^{v_{m}d_{\alpha}^{m}}}$. The solution for $x$ is given
by the Lambert-W function [50] as $x=-\theta-\frac{1}{\phi}\mathcal{W}(-\phi
e^{-\phi\theta}/a_{0})$, hence

$\displaystyle\omega_{0}=v_{m}-\theta-\frac{1}{\phi}\mathcal{W}\left(-\phi\theta
e^{-\phi\theta}\right).$ (23)

Finally, since the false alarm rate (i.e., frequency) is the inverse of false
alarm period $E_{\infty}[T]$, we have

$FAR\leq e^{-\omega_{0}h},$

where $h$ is the detection threshold, and $\omega_{0}$ is given above.

Figure 3: Actual false alarm periods vs. derived lower bounds for the PSM,
SMAP and SWAT datasets respectively.

Specifically, $v_{m}$ is directly computed using the dimensionality $m$,
$d_{\alpha}$ comes from the training phase, $\phi$ is also found in training,
and finally there is a built-in Lambert-W function in popular programming
languages such as Python and Matlab. Hence, given the training data,
$\omega_{0}$ can be easily computed, and based on Theorem 1, the threshold $h$
can be chosen to asymptotically achieve the desired false alarm period as
follows

$h=\frac{-\log(FAR)}{\omega_{0}}$ (24)

We finally present the comparison between the bound for false alarm rate
derived in Theorem 1 and the empirical false alarm period in Fig. 3. The
figure depicts the logarithm of false alarm period, which is the inverse of
false alarm rate, for clarity. Hence, the upper bound on the false alarm rate
becomes the lower bound on the false alarm period in this scenario.

##  V Experiments

Datasets: Most of the existing works evaluate their performance on the
following datasets:

  * •

SMD (Server Machine Dataset): The SMD dataset is collected from a large
internet company and consists of data collected over 5 weeks with 38
dimensions.

  * •

PSM (Pooled Server Metrics): The PSM dataset is proposed by eBay and consists
of 26 dimensional data captured internally from application server nodes.

  * •

MSL (Mars Science Laboratory rover) and SMAP (Soil Moisture Active Passive
satellite): The MSL and SMAP datasets are provided by NASA and consists of
telemetry data and anomalies featuring 55 and 25 dimensions respectively.
Since most of the dimensions are categorical, we only focus on the telemetry
values.

  * •

SWaT: The SWaT dataset is collected in an industrial setting and features data
collected from a sewage water treatment facility. The dataset is collected
over an entire week and consists of 51 dimensions, where the anomalies are
caused due to cyberattacks.

Each dataset includes training, validation and testing subsets. Anomalies are
only labeled in the testing subset.

Implementation Details: Following the sliding window approach used in existing
works, the input to the proposed TiSAT model is a sub-series with a window
size of 100. The TiSAT model encoder consists of a 3-layer stack followed by a
1-layer stack and the decoder consists of a 2-layer stack. We train the model
using Adam optimizer, and the learning rate is set as 1$e^{-4}$. We train the
model for 4 epochs with a batch size of 64. We normalize all the datasets
between [0,1] using minmax normalization.

Results: We extensively evaluate the performance of the proposed approach on
the five publicly available benchmark datasets using the proposed SPD metric.
We compare our approaches with classical algorithms such as OC-SVM,
IsolationForest and LOF, as well as time series forecasting based approaches
such as ARIMA and LSTM. As shown in Table II, the proposed TiSAT model
significantly outperforms all other approaches. While there have been several
neural network based approaches [15, 14, 13, 22, 23, 24, 16] proposed
recently, all of them present their results using the faulty adjusted
instance-based metric shown in Fig. 1 (See Table I for the inherent flaw of
this metric). Since most of them do not have their codes available, it is not
feasible for us to compare with them.

Dataset | SMD | MSL | SMAP | SWaT | PSM
---|---|---|---|---|---
Metric (SPD) |  |  |  |  |
OCSVM [6] | 21.08 | 20.15 | 12.82 | 17.21 | 19.44
IsolationForest [7] | 19.38 | 16.50 | 9.26 | 12.83 | 18.98
LOF [8] | 17.25 | 14.29 | 18.21 | 13.68 | 20.14
ARIMA | 27.64 | 24.13 | 22.87 | 29.35 | 26.47
LSTM [33] | 28.31 | 23.85 | 27.54 | 32.97 | 28.41
Ours | 52.70 | 33.40 | 29.35 | 46.23 | 38.58
TABLE II: Comparison with existing approaches using the proposed SPD metric.

##  VI Conclusion

In this work we identified a crucial shortcoming in the existing evaluation
criterion used by most recent approaches for time series anomaly detection. To
rectify the evaluation method, we presented a novel performance metric which
measures the timeliness and precision of detection methods. Moreover, we
proposed a novel transformer based approach called TiSAT for unsupervised time
series anomaly detection and provided an asymptotic false alarm rate analysis
for TiSAT. This analysis leads to a closed-form expression for the detection
threshold, which was emppirically corroborated on benchmark datasets. We
comprehensively evaluated the proposed approach and showed that TiSAT is able
to achieve state-of-the-art performance on benchmark datasets.

## References

  * [1] S. Chauhan and L. Vig, “Anomaly detection in ecg time signals via deep long short-term memory networks,” in 2015 IEEE International Conference on Data Science and Advanced Analytics (DSAA), pp. 1–7, IEEE, 2015.
  * [2] K. Doshi, Y. Yilmaz, and S. Uludag, “Timely detection and mitigation of stealthy ddos attacks via iot networks,” IEEE Transactions on Dependable and Secure Computing, 2021.
  * [3] N. Laptev, S. Amizadeh, and I. Flint, “Generic and scalable framework for automated time-series anomaly detection,” in Proceedings of the 21th ACM SIGKDD international conference on knowledge discovery and data mining, pp. 1939–1947, 2015.
  * [4] H. Ren, B. Xu, Y. Wang, C. Yi, C. Huang, X. Kou, T. Xing, M. Yang, J. Tong, and Q. Zhang, “Time-series anomaly detection service at microsoft,” in Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pp. 3009–3017, 2019.
  * [5] M. Braei and S. Wagner, “Anomaly detection in univariate time-series: A survey on the state-of-the-art,” arXiv preprint arXiv:2004.00433, 2020.
  * [6] B. Schölkopf, J. C. Platt, J. Shawe-Taylor, A. J. Smola, and R. C. Williamson, “Estimating the support of a high-dimensional distribution,” Neural computation, vol. 13, no. 7, pp. 1443–1471, 2001.
  * [7] F. T. Liu, K. M. Ting, and Z.-H. Zhou, “Isolation forest,” in 2008 eighth ieee international conference on data mining, pp. 413–422, IEEE, 2008\.
  * [8] M. M. Breunig, H.-P. Kriegel, R. T. Ng, and J. Sander, “Lof: identifying density-based local outliers,” in Proceedings of the 2000 ACM SIGMOD international conference on Management of data, pp. 93–104, 2000.
  * [9] L. Ruff, R. Vandermeulen, N. Goernitz, L. Deecke, S. A. Siddiqui, A. Binder, E. Müller, and M. Kloft, “Deep one-class classification,” in International conference on machine learning, pp. 4393–4402, PMLR, 2018.
  * [10] B. Zong, Q. Song, M. R. Min, W. Cheng, C. Lumezanu, D. Cho, and H. Chen, “Deep autoencoding gaussian mixture model for unsupervised anomaly detection,” in International conference on learning representations, 2018.
  * [11] D. Park, Y. Hoshi, and C. C. Kemp, “A multimodal anomaly detector for robot-assisted feeding using an lstm-based variational autoencoder,” IEEE Robotics and Automation Letters, vol. 3, no. 3, pp. 1544–1551, 2018.
  * [12] B. Zhou, S. Liu, B. Hooi, X. Cheng, and J. Ye, “Beatgan: Anomalous rhythm detection using adversarially generated time series.,” in IJCAI, pp. 4433–4439, 2019.
  * [13] S. Li, X. Jin, Y. Xuan, X. Zhou, W. Chen, Y.-X. Wang, and X. Yan, “Enhancing the locality and breaking the memory bottleneck of transformer on time series forecasting,” Advances in Neural Information Processing Systems, vol. 32, 2019.
  * [14] Z. Li, Y. Zhao, J. Han, Y. Su, R. Jiao, X. Wen, and D. Pei, “Multivariate time series anomaly detection and interpretation using hierarchical inter-metric and temporal embedding,” in Proceedings of the 27th ACM SIGKDD Conference on Knowledge Discovery & Data Mining, pp. 3220–3230, 2021.
  * [15] L. Shen, Z. Li, and J. Kwok, “Timeseries anomaly detection using temporal hierarchical one-class network,” Advances in Neural Information Processing Systems, vol. 33, pp. 13016–13026, 2020.
  * [16] J. Xu, H. Wu, J. Wang, and M. Long, “Anomaly transformer: Time series anomaly detection with association discrepancy,” arXiv preprint arXiv:2110.02642, 2021.
  * [17] K. Hundman, V. Constantinou, C. Laporte, I. Colwell, and T. Soderstrom, “Detecting spacecraft anomalies using lstms and nonparametric dynamic thresholding,” in Proceedings of the 24th ACM SIGKDD international conference on knowledge discovery & data mining, pp. 387–395, 2018.
  * [18] S. Hochreiter and J. Schmidhuber, “Long short-term memory,” Neural computation, vol. 9, no. 8, pp. 1735–1780, 1997.
  * [19] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, and I. Polosukhin, “Attention is all you need,” in Advances in neural information processing systems, pp. 5998–6008, 2017.
  * [20] S. Chaudhari, V. Mithal, G. Polatkan, and R. Ramanath, “An attentive survey of attention models,” arXiv preprint arXiv:1904.02874, 2019.
  * [21] H. Song, D. Rajan, J. J. Thiagarajan, and A. Spanias, “Attend and diagnose: Clinical time series analysis using attention models,” in Thirty-second AAAI conference on artificial intelligence, 2018.
  * [22] H. Ren, B. Xu, Y. Wang, C. Yi, C. Huang, X. Kou, T. Xing, M. Yang, J. Tong, and Q. Zhang, “Time-series anomaly detection service at microsoft,” in Proceedings of the 25th ACM SIGKDD international conference on knowledge discovery & data mining, pp. 3009–3017, 2019.
  * [23] Y. Su, Y. Zhao, C. Niu, R. Liu, W. Sun, and D. Pei, “Robust anomaly detection for multivariate time series through stochastic recurrent neural network,” in Proceedings of the 25th ACM SIGKDD international conference on knowledge discovery & data mining, pp. 2828–2837, 2019.
  * [24] H. Xu, W. Chen, N. Zhao, Z. Li, J. Bu, Z. Li, Y. Liu, Y. Zhao, D. Pei, Y. Feng, et al., “Unsupervised anomaly detection via variational auto-encoder for seasonal kpis in web applications,” in Proceedings of the 2018 world wide web conference, pp. 187–196, 2018.
  * [25] V. Hodge and J. Austin, “A survey of outlier detection methodologies,” Artificial intelligence review, vol. 22, no. 2, pp. 85–126, 2004.
  * [26] V. Chandola, A. Banerjee, and V. Kumar, “Anomaly detection: A survey,” ACM computing surveys (CSUR), vol. 41, no. 3, pp. 1–58, 2009.
  * [27] T. Wen and R. Keyes, “Time series anomaly detection using convolutional neural networks and transfer learning,” arXiv preprint arXiv:1905.13628, 2019\.
  * [28] M. Canizo, I. Triguero, A. Conde, and E. Onieva, “Multi-head cnn–rnn for multi-time series anomaly detection: An industrial case study,” Neurocomputing, vol. 363, pp. 246–260, 2019.
  * [29] M. Munir, S. A. Siddiqui, A. Dengel, and S. Ahmed, “Deepant: A deep learning approach for unsupervised anomaly detection in time series,” Ieee Access, vol. 7, pp. 1991–2005, 2018.
  * [30] L. Bontemps, V. L. Cao, J. McDermott, and N.-A. Le-Khac, “Collective anomaly detection based on long short-term memory recurrent neural networks,” in International conference on future data and security engineering, pp. 141–152, Springer, 2016.
  * [31] L. Zhu and N. Laptev, “Deep and confident prediction for time series at uber,” in 2017 IEEE International Conference on Data Mining Workshops (ICDMW), pp. 103–110, IEEE, 2017.
  * [32] N. Laptev, J. Yosinski, L. E. Li, and S. Smyl, “Time-series extreme event forecasting with neural networks at uber,” in International conference on machine learning, vol. 34, pp. 1–5, 2017.
  * [33] P. Malhotra, A. Ramakrishnan, G. Anand, L. Vig, P. Agarwal, and G. Shroff, “Lstm-based encoder-decoder for multi-sensor anomaly detection,” arXiv preprint arXiv:1607.00148, 2016.
  * [34] S.-V. Oprea, A. Bâra, F. C. Puican, and I. C. Radu, “Anomaly detection with machine learning algorithms and big data in electricity consumption,” Sustainability, vol. 13, no. 19, p. 10963, 2021.
  * [35] Y. Liang, S. Ke, J. Zhang, X. Yi, and Y. Zheng, “Geoman: Multi-level attention networks for geo-sensory time series prediction.,” in IJCAI, vol. 2018, pp. 3428–3434, 2018.
  * [36] Y. Qin, D. Song, H. Chen, W. Cheng, G. Jiang, and G. Cottrell, “A dual-stage attention-based recurrent neural network for time series prediction,” arXiv preprint arXiv:1704.02971, 2017.
  * [37] K. Zhou, W. Wang, T. Hu, and K. Deng, “Time series forecasting and classification models based on recurrent with attention mechanism and generative adversarial networks,” Sensors, vol. 20, no. 24, p. 7211, 2020\.
  * [38] T. B. Brown, B. Mann, N. Ryder, M. Subbiah, J. Kaplan, P. Dhariwal, A. Neelakantan, P. Shyam, G. Sastry, A. Askell, et al., “Language models are few-shot learners,” arXiv preprint arXiv:2005.14165, 2020.
  * [39] G. Zerveas, S. Jayaraman, D. Patel, A. Bhamidipaty, and C. Eickhoff, “A transformer-based framework for multivariate time series representation learning,” in Proceedings of the 27th ACM SIGKDD Conference on Knowledge Discovery & Data Mining, pp. 2114–2124, 2021.
  * [40] B. Lim, S. Ö. Arık, N. Loeff, and T. Pfister, “Temporal fusion transformers for interpretable multi-horizon time series forecasting,” International Journal of Forecasting, 2021.
  * [41] S. Wu, X. Xiao, Q. Ding, P. Zhao, W. Ying, and J. Huang, “Adversarial sparse transformer for time series forecasting,” 2020.
  * [42] Z. Chen, D. Chen, X. Zhang, Z. Yuan, and X. Cheng, “Learning graph structures with transformer for multivariate time series anomaly detection in iot,” IEEE Internet of Things Journal, 2021.
  * [43] H. Zhou, S. Zhang, J. Peng, S. Zhang, J. Li, H. Xiong, and W. Zhang, “Informer: Beyond efficient transformer for long sequence time-series forecasting,” in Proceedings of AAAI, 2021.
  * [44] L. Guo, R. Li, and B. Jiang, “A data-driven long time-series electrical line trip fault prediction method using an improved stacked-informer network,” Sensors, vol. 21, no. 13, p. 4466, 2021.
  * [45] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, and I. Polosukhin, “Attention is all you need,” in Advances in neural information processing systems, pp. 5998–6008, 2017.
  * [46] F. Yu, V. Koltun, and T. Funkhouser, “Dilated residual networks,” in Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 472–480, 2017.
  * [47] A. Gupta and A. M. Rush, “Dilated convolutions for modeling long-distance genomic dependencies,” arXiv preprint arXiv:1710.01278, 2017.
  * [48] M. Basseville and I. Nikiforov, Detection of abrupt changes: theory and application, vol. 104.  Prentice Hall, Englewood Cliffs, 1993.
  * [49] S. N. Chiu, D. Stoyan, W. S. Kendall, and J. Mecke, Stochastic geometry and its applications.  John Wiley & Sons, 2013.
  * [50] T. C. Scott, G. Fee, and J. Grotendorst, “Asymptotic series of generalized lambert w function,” ACM Communications in Computer Algebra, vol. 47, no. 3/4, pp. 75–83, 2014.