text:
\section{Introduction}

Quite recently, "buckyballs" have been observed for the first time in space by
Astronomers of University of Western Ontario \cite{nasajully10,jcjbsepsem2010}
using NASA's Spitzer Space Telescope. They identified spectral signatures
of fullerenes in a cloud of cosmic dust surrounding a distant star 6500 light
years away. The existence of fullerenes has thus been
confirmed as anticipated after the first observation in laboratory by
Kroto {\it et al.} at Rice University \cite{hwkjrhscobrfcres85}. One may therefore
expect observation of molecular species trapped in fullerenes in interstellar medium.
In this work we investigate on the changes that are induced in the spectral signatures
of species trapped in fullerene in order to determine how these changes can be used
to characterize the environment where fullerenes are observed. From inert matrix
isolation spectroscopy, one knows that infrared (IR) spectra of molecular
species trapped in a site are simplified in that, as rotational structures
are absent in the spectra, a vibrational transition
is observed as one peak blue or red shifted or sometimes two or more in case of
splitting of degenerate levels \cite{cgal89_1,prdibvrjlthcallam99,prdalhcjmc2006}.

Trapping of atoms or ions and small diatomic molecules in fullerene have been studied
both experimentally and theoretically
\cite{dsbrdjjrsmsdvcsy93,hstwdkbjh93,mshajvrjcrjp93,hajvrjcmsrjp94,
mshajvrjcsmmlgdeg94,rsrjcms97,rjcakms2003,jc91,ciwmawlp93,jhrjbjmgl96_1,jhrjbjmgl96_2,
ehtoavdapesw96,rjc2008,mxfszbrlnjt2008_1,mxfszbrlnjt2008_2,zsfulasn2005}.
However, except for CO \cite{jhrjbjmgl96_2,ehtoavdapesw96,rjc2008}
and H$_2$ \cite{mxfszbrlnjt2008_1,mxfszbrlnjt2008_2} molecules, only a few
of the theoretical work has been done in the infrared region. In the case
of NH$_3$ trapped in fullerene, calculations have been performed
from {\it ab-initio} methods to study the trapping
mechanism \cite{zsfulasn2005,selt2003,mdgmmog2009} but no experimental
work has been published. Erko\c{c} and T\"{u}rker \cite{selt2003} have
reported that up to six NH$_3$ molecules can be trapped
in C$_{60}$, on one hand and Ganji {\it et al.} \cite{mdgmmog2009} suggest
that only one NH$_3$ molecule can form a stable
complex NH$_3$@C$_{60}$ on the other hand. While Slanina {\it et al.}
\cite{zsfulasn2005} have reported a binding energy of -1830 cm$^{-1}$ with
the ammonia molecule oriented towards a pair of parallel pentagons.

The aim of this paper is to determine how fullerene modifies the spectral signature
of ammonia molecule, in particular in the region of the vibration-inversion
mode (namely umbrella mode).

A theoretical study of an ammonia molecule trapped in fullerene C$_{60}$ is thus
performed in order to simulate its infrared spectrum in the vibration-inversion
frequency region.
The choice of NH$_3$ is motivated by the fact that it contains
a nitrogen atom. CO$_2$ could otherwise have been chosen for the simulation.

To simulate the infrared spectra of NH$_3$ in fullerene, we use the site
inclusion model successfully applied to analyze spectra of CO$_2$ isotopologues
\cite{prdalhcjmc2006} isolated in rare gas matrix. When the molecule is
trapped in the nano-cage of fullerene, the interaction with carbon atoms in the
short range scale modifies the translational and rotational degrees of freedom
because the molecule is no longer free to move. As a result, the umbrella mode
which is characterized by the inversion doubling through tunneling of the
nitrogen atom across the plane defined by the three H atoms is perturbed
as discussed below.

To calculate the bar-spectra of NH$_3$ trapped in C$_{60}$ nano-cage a
renormalization procedure is applied on the total Hamiltonian of the system. The
dynamical coupling between the molecular degrees of freedom and those of C$_{60}$
nano-cage is then small enough, at least to a first approximation, to justify
an \textit{initial chaos hypothesis}.
As a result, the density matrix operator of the system can be factorized into a
product consisting of two terms, one related to the NH$_3$ renormalized optical
states and the other to the C$_{60}$ bath states spanned by its vibrational modes.
The dynamical line-shifts and line-widths of the bar-spectra can be disregarded and
C$_{60}$ nano-cage considered as rigid.

The interaction potential model used in our calculations is
presented in section 2 and the renormalization method applied to separate
the different motions from each other described. Results of the calculation
to determine the frequency shifts of the vibrational modes and the orientational
level schemes of the molecule are then given. Sections 3 and 4 are devoted to
a presentation of the formalism used for the construction of the infrared
spectra of NH$_3$ in fullerene and to a discussion of the infrared bar-spectrum
in the frequency region of the $\nu_2$ mode.

\section{The interaction potential energy}

\subsection{Potential energy model}

The interaction potential energy $V_{\text{MC}}$ between the trapped ammonia molecule NH$_3$
and the rigid fullerene molecule C$_{60}$ is modelled as the sum of 12-6 Lennard-Jones (LJ)
pairwise atom-atom potentials characterizing the repulsion-dispersion contributions and the
induction part due to the interaction between the permanent electric multipole moments of the
molecule and their images created at the positions of the polarized fullerene carbon atoms.
It can be expressed as

\begin{eqnarray}
V_{\text{MC}} = \sum\limits_{j=1}^{60}\sum\limits_{i=1}^44\epsilon_{ij}\left\{ \left(
\frac{\sigma _{ij}}{\left| \mathbf{r}_{ij}\right| }\right) ^{12}
-\left( \frac{\sigma _{ij}}{\left| \mathbf{r}_{ij}\right| }\right) ^6 \right\}
-\frac 12\sum\limits_{j=1}^{60} \textbf{E}_{\text{M}}^j :\mathbf{\alpha }^j :\textbf{E}_{\text{M}}^j,
\label{eq1}
\end{eqnarray}

\noindent where \textit{i} and \textit{j} denote the \textit{i}th atom of the ammonia
molecule and the \textit{j}th carbon atom of the fullerene molecule, respectively;
$\epsilon_{ij}$ and $\sigma _{ij}$ are the mixed LJ potential parameters, obtained from the
usual Lorentz-Berthelot combination rules $\epsilon _{ij}=\sqrt{\epsilon _{ii}\epsilon _{jj}}$
and $2\sigma_{ij}=\sigma _{ii}+\sigma _{jj}$ and $\textbf{r}_{ij}$ is the distance vector
between the \textit{i}th atom of the molecule and the \textit{j}th carbon atom.

In the second term of Eq. (\ref{eq1}) $\textbf{E}_{\text{M}}^j$ is the field generated by
the molecular permanent electric multipole moments ($\mathbf{\mu }$, $\mathbf{\Theta }$,...)
of the NH$_3$ molecule at position \textbf{r}$_0$ on the \textit{j}th carbon atom of the
fullerene molecule at position \textbf{r}$_j$ with polarizability
tensor $\mathbf{\alpha }^j$. It is expressed as

\begin{eqnarray}
\textbf{E}_{\text{M}}^j = \mathbf{\nabla\nabla }\left( \frac 1{\left| \mathbf{r}_j-\mathbf{r}_0\right|
}\right) \cdotp\mathbf{\mu} + \frac 13 \mathbf{\nabla\nabla
\nabla }\left( \frac 1{\left| \mathbf{r}_j-\mathbf{r}_0\right| }\right)
\mathbf{:\Theta }+... \label{eq2}
\end{eqnarray}

Note that the multipole moments of NH$_3$ and the polarizability tensor of each carbon atom
of C$_{60}$ are usually defined and given with respect to local frames with
their origins at the centre of mass of NH$_3$ (G,\textbf{x},\textbf{y},\textbf{z}) or
each carbon atom (C$_j$,\textbf{x}$_j$,\textbf{y}$_j$,\textbf{z}$_j$) respectively.
In the calculation, it is then necessary to express these quantities in an absolute
frame (O,\textbf{X},\textbf{Y},\textbf{Z}) defined with respect to the fixed C$_{60}$
molecule. We use a transformation rotation matrix following Rose's
convention \cite{mer67} and given in Appendix B.

Figure 1 indicates how to describe the internal degrees of freedom of the
molecule with respect to its frame (G,\textbf{x},\textbf{y},\textbf{z})
and its external orientational and translational degrees of freedom with respect
to the absolute frame, and Table 1 gives the various fullerene
and ammonia characteristic values used in our calculations.

\begin{figure}[h!]
\begin{center}
\includegraphics[width=16cm]{Figure1.eps}
\end{center}
\caption{Geometrical characteristics of a molecule trapped in the
fullerene molecule. (O,\textbf{X},\textbf{Y},\textbf{Z})
and (G,\textbf{x},\textbf{y},\textbf{z}) represent the absolute
frame and the molecular frame, respectively. \textit{q}
and $\beta$ define the equilibrium internal coordinates of the trapped molecule.}
\label{FIG. 1}
\end{figure}

\begin{table}[h!]
\begin{center}
\caption{Numerical parameters for NH$_3$ molecule and
for C$_{60}$ molecule used in our calculations.}
\begin{tabular}{llll}
\hline\hline
NH$_3$ & &  C$_{60}$  & \\
\hline
$q$ (\AA ) & 1.0156 &  $R^{(a)}$ (\AA ) & 3.556 \\
$\beta$ (deg.) & 68. &  d(C$-$C)$^{(a)}$ (\AA ) & 1.455 \\
 &  & d(C$=$C)$^{(a)}$ (\AA ) & 1.398 \\
$\mu ^{\mathrm{e}}$ (D) & 1.476 &  $\alpha _{\bot }^{\mathrm{c}}$
(\AA $^3$) & 1.44 \\
$\Theta ^{\mathrm{e}}$ (D\AA ) & -2.930 &  $\alpha _{\Vert }^{\mathrm{c
}$ (\AA $^3$) & 0.41 \\
$B$ (cm$^{-1}$) & 9.941 & &  \\
$C$ (cm$^{-1}$) & 6.309 & &  \\
\hline
\hline

 & C - C$^{(b)}$ &  N - N$^{(c)}$ & H - H$^{(c)}$   \\
\hline
$\epsilon$ (cm$^{-1}$) & 19.4 & 26.5 & 17.2 \\
$\sigma$ (\AA ) & 3.40 & 3.38 & 2.53 \\

\hline

\raggedright $(a)$ From Ref. \cite{khlhdbcabhcdrdjmdv91}.
\raggedright $(b)$ From Ref. \cite{hysxwz2009}.
\raggedright $(c)$ From Ref. \cite{al2000}.

\end{tabular}
\end{center}
\label{Table 1}
\end{table}

The distance vector $\mathbf{r}_{ij}$ in Eq.(\ref{eq1}) can be expressed in
terms of the position vectors \textbf{r}$_0$ of the molecular centre of mass (c.m.)
and $\mathbf{r}_j$ of the $j$th carbon atom of the fullerene molecule with respect
to the absolute frame (O,\textbf{X},\textbf{Y},\textbf{Z}) and of the position
vector $\mathbf{\eta }_i$ of the $i$th atom of the molecule with respect to its
associated frame (G,\textbf{x},\textbf{y},\textbf{z}):

\begin{equation}
\mathbf{r}_{ij} = \mathbf{r}_j-\mathbf{r}_0-\mathbf{\eta }_i, \label{eq3}
\end{equation}

The IR spectra of the trapped NH$_3$ molecule are usually given in terms of
molecular vibrational \{$Q$\} degrees of freedom.

Thus, to determine the influence of the surroundings on the molecular
internal motions, the position vectors $\mathbf{\eta }_i$ of the atoms in the
molecule, the molecular dipole moment vector $\mathbf{\mu }$, and the
quadrupole moment tensor $\mathbf{\Theta }$ can all be expressed in a series
expansion with respect to the molecular frame in terms of the vibrational
normal coordinates \{$Q$\}.

Taking only the first order terms gives

\begin{eqnarray}
\mathbf{\eta }_i &=&\mathbf{\eta }_i^{\text{e}}+\sum\limits_\nu \mathbf{a
_i^\nu Q_\nu , \nonumber \\
\mathbf{\mu } &=&\mathbf{\mu }^{\text{e}}+\sum\limits_\nu \mathbf{b}^\nu
Q_\nu , \nonumber \\
\mathbf{\Theta } &=&\mathbf{\Theta }^{\text{e}}+\sum\limits_\nu \mathbf{c
^\nu Q_\nu . \label{eq4}
\end{eqnarray}

In these expressions $\mathbf{\eta }_i^{\text{e}}$ is the vector position of
the \textit{i}th atom in the rigid molecule, and $\mathbf{\mu }^{\text{e}}$
and $\mathbf{\Theta }^{\text{e}}$ are its permanent multipole moments,
while $\mathbf{a}_i^\nu $, $\mathbf{b}^\nu $, and $\mathbf{c}^\nu $ are the
first derivatives of $\mathbf{\eta }_i$, $\mathbf{\mu }$, and $\mathbf{\Theta }$
with respect to the normal coordinate $Q_\nu $ which describes the $\nu $th
molecular vibrational mode with frequency $\omega _\nu $.
Note that the expression of $\mathbf{\mu }$ will also be used in the near and
far infrared spectra calculations.

However, it must be noticed that because of the \textit{dual nature} of the
$\nu_2$ mode of the ammonia molecule, that is, the \textit{high frequency vibration-low
frequency inversion}, this mode will be treated in a particular way.

Using the Born-Oppenheimer approximation (adiabatic approximation) to separate the
high frequency vibrational modes \{$Q$\} of the molecule from its low frequency
external modes, that is, the orientation $\mathbf{\Omega}=(\varphi,\theta,\chi)$ and
c.m. translation \textbf{r}$_0$ motions, the interaction
potential energy $V_{\text{MC}}$ can be written as

\begin{equation}
V_{\text{MC}} = V_{\text{M}}(\mathbf{r}_0,\mathbf{\Omega}) +
W_{\text{M}}(\{Q\}) + \bigtriangleup V_{\text{M}}(\mathbf{r}_0,\mathbf{\Omega},
\{Q\}), \label{eq5}
\end{equation}

\noindent where $V_{\text{M}}(\mathbf{r}_0,\mathbf{\Omega})$ represents the low frequency
motions dependent part of the potential energy for the non vibrating molecule,
and $W_{\text{M}}(\{Q\})$ characterizes the vibrational dependent part for
the molecule at its equilibrium position and orientation, and is generally a small
perturbation which involves only vibrational level shifts and/or splittings.
It can be developed in a Taylor series expansion up to second order with respect to
the vibrational normal coordinates \{$Q$\} as

\begin{equation}
W_{\text{M}}(\{Q\}) = \sum\limits_\nu \frac{\partial
V_{\text{MC}}}{\partial Q_\nu }Q_\nu +\frac 12\sum\limits_{\nu \nu ^{\prime }}\frac{
\partial ^2V_{\text{MC}}}{\partial Q_\nu \partial Q_{\nu ^{\prime }}}Q_\nu Q_{\nu
^{\prime }}+... \label{eq6}
\end{equation}

Finally the third term in Eq.(\ref{eq5}) characterizes the coupling between the
external and the internal modes which can lead to the vibrational energy relaxation.

\subsection{Potential energy surfaces}

\subsubsection{Orientation-translation motions}

The first step of the numerical calculations consists in determining the equilibrium
configuration of the rigid ammonia molecule into the C$_{60}$ cage. The potential
energy $V_{\text{M}}$ (see Eq.(\ref{eq5}) is minimized with
respect to both the c.m. displacement vector \textbf{r}$_0$ and
the molecular orientional coordinates $\varphi$, $\theta$, and $\chi$.

We find that the equilibrium configuration has an energy minimum value of
- 1460 cm$^{-1}$ for the \textbf{z} molecular axis (threefold symetry axis C$_3$)
oriented along the direction connected to each two
symmetrical carbon atoms. The c.m. is then displaced from the site centre by
about 0.184 $\text{\AA}$ in the same direction.
An energy maximum of - 1445 cm$^{-1}$ is obtained when the molecule is oriented and
its c.m. displaced by 0.188 $\text{\AA}$ along the pentagon centre directions.
There is also an intermediate energy value of - 1455 cm$^{-1}$ that is calculated for
the molecule oriented and its c.m. displaced by 0.180 $\text{\AA}$ along the hexagon
centre directions.

Thus, the calculations show that \textit{i}) in all the molecular orientational
configurations, the c.m. displacements \textit{r}$_0$ have always opposite directions
to the molecular \textbf{z} axis, and \textit{ii}) the energy changes are negligibly
small $\lesssim$ 3 cm$^{-1}$ when \textit{r}$_0$ changes by about 0.01 $\text{\AA}$
along these directions. So, in the following we will take the value
\textit{r}$_{0}^{\text{e}}$ = 0.184 $\text{\AA}$ as a mean value for every orientation.
This means that the trajectory of the molecular c.m. can be approximated by the
surface of a sphere of radius \textit{r}$_{0}^{\text{e}}$.

Note, moreover, that the induction and repulsion-dispersion contributions
represent, respectively, 35$\%$ (- 498 cm$^{-1}$) and 65$\%$ (- 963 cm$^{-1}$) of
the total energy.

As a result, the potential energy surface
$V_{\text{M}}(\mathit{r}_0(\mathbf{\Omega}))$ is nearly flat since the barrier
height is always $\lesssim$ 16 cm$^{-1}$ and the c.m. displacement vectors \textbf{r}$_0$
from the cage centre remain always nearly anticollinear to the molecular \textbf{z}
axis \textit{i.e.} \textbf{z}$\cdotp$ \textbf{r}$_0$ $\simeq$ - \textit{r}$_0$.
This is not surprising because of the repulsive character of the hydrogen atoms
of the molecule and reflects a very strong coupling between its orientation
and the location of its centre of mass.

Thus, on one hand, the axis of the molecule undergoes only oscillatory motions
($\varphi,\theta$) of small amplitude around the direction of each displacement vector
\textbf{r}$_{0}^{\text{e}}$. The corresponding potential energy surface which is given in
Figure 2 is a nearly harmonic two dimension librational oscillator.
The associated frequencies are $\omega_{\varphi}$ $\simeq$ 155 cm$^{-1}$ and
$\omega_{\theta}$ $\simeq$ 151 cm$^{-1}$.

\begin{figure}[h!]
\begin{center}
\includegraphics[width=10cm]{Figure2.eps}
\end{center}
\caption{Potential energy surface \textit{versus} ($\varphi,\theta$) angular motions
around the molecular equilibrium configuration
($\mathbf{r}_{0}^{\text{e}}$,$\mathbf{\Omega}$).}
\label{FIG. 2}
\end{figure}

On the other hand, for a given molecular orientation, the radial potential
function $V_{\text{M}}(\textit{r}_0)$ of the translation motion of the molecular c.m.
around the equilibrium position $r_{0}^{\text{e}}$ = 0.184 $\text{\AA}$ is
plotted in Figure 3. The curve appears as that of a nearly harmonic one dimensional
oscillator with a force
constant value \textit{k} = 2.4$\times$10$^4$ cm$^{-1} \ldotp \text{\AA}^{-2}$
corresponding to a frequency $\omega_{\text{rad}}$ $\simeq$ 218 cm$^{-1}$.

\begin{figure}[h!]
\begin{center}
\includegraphics[width=10cm]{Figure3.eps}
\end{center}
\caption{Potential energy function \textit{versus}
molecular c.m. displacement $r_0$ around its equilibrium position
vector \textbf{r}$_{0}^{\text{e}}$.}
\label{FIG. 3}
\end{figure}

On the basis of these features, we can conclude that ammonia should exhibit a
nearly free rotation motion simultaneously around the centre of the cage by describing
a spherical curve with radius \textit{r}$_{0}^{\text{e}}$ = 0.184 $\text{\AA}$, and
around its centre of mass. Note however that the spinning motion $\chi$ remains around
the threefold symmetry axis of the molecule, only.

Finally, it is clear that this rotation-c.m. translation motion is slow with respect to the
librational and the radial motions and can be separated from the latter using the
adiabatic approximation.

\subsubsection{Inversion-translation motions}

It is also clear, from the results above, that when the inversion motion of the molecule
proceeds, its centre of mass moves simultaneously from one position
(\textbf{r}$_{0}^{\text{e}},\mathbf{\Omega}$) to the opposite one
(-\textbf{r}$_{0}^{\text{e}},\mathbf{\Omega}$) with respect to the cage
centre where the molecule lies in its planar configuration.

\subsubsection{Few remarks}

It must be noticed that the interaction potential energy experienced by the ammonia
molecule trapped in the fullerene molecule strongly depends on the geometrical
characteristics such as the radius of the C$_{60}$ cage and the equilibrium
internal coordinates (\textit{q},$\beta$ see Fig. 1) of the ammonia molecule, and
also on the dispersion-repulsion potential parameters especially the
Lennard-Jones $\sigma$ parameters.

For instance, in one hand, the minimum depth is enhanced from - 1460
to - 750 cm$^{-1}$ when the fullerene radius is decreased
by 1$\%$, from 3.556 to 3.519 $\text{\AA}$, and it diminishes
to - 2000 and - 2411 cm$^{-1}$ when the
radius value is increased by 1$\%$ and 2$\%$, to 3.592 and 3.627 $\text{\AA}$,
respectively. The induction energy contribution undergoing only weakly variations.

In the other hand, changes of the $\sigma$ parameters of + 1$\%$ and - 1$\%$
involve, respectively, minimum depths of - 952 and - 1976 cm$^{-1}$.

However, in spite of these changes in the well depth values, the potential energy surface
$V_{\text{M}}(\mathbf{r}_{0}^{\text{e}},\mathbf{\Omega})$ remains nearly flat since
its barrier height does not exceed 24 cm$^{-1}$. Furthermore, the potential energy surfaces
associated with the various motions remain nearly unchanged.

Thus, from the previous analysis of the potential energy experienced by the
ammonia molecule into the fullerene molecule, we will assume that the c.m. position
vector \textbf{r}$_{0}$ remains permanently anticollinear to the \textbf{z} molecular
axis and separate it into a static position vector
\textbf{r}$_{0}^{\text{e}}(\varphi,\theta,Q_2)$ which parametrically depends on
the orientation and inversion degrees of freedom, and a dynamical
coordinate \textbf{u}$_{0}$ characterizing the c.m. vibration
around \textbf{r}$_{0}^{\text{e}}(\varphi,\theta,Q_2)$.

Finally, it is clear that inversion-c.m. translation motions
(\textit{r}$_{0}^{\text{e}}$,$Q_2$), on one hand, and
orientation-c.m. translation motions (\textbf{r}$_{0}^{\text{e}}$,$\mathbf{\Omega}$),
on the other hand, are strongly coupled; their kinetic operator must be classically
treated, as explained in appendix A, before their quantum mechanical treatment as
described below. The interaction potential energy $V_{\text{MC}}$ can be rewritten
as

\begin{equation}
V_{\text{MC}} = V_{\text{M}}(\mathbf{r}_{0}^{\text{e}},\mathbf{\Omega}) +
V_{\text{M}}(\textit{r}_{0}) + V_{\text{M}}(\textit{r}_{0},Q_2) +
W_{\text{M}}(\{Q\}) + \bigtriangleup V_{\text{M}}(\mathbf{r}_0,\mathbf{\Omega},
\{Q\}), \label{eq7}
\end{equation}

\section{Quantum-classical model}

As we have mentioned above, the purpose of this work is to determine the infrared
spectra of the NH$_3$ molecule trapped in a rigid C$_{60}$ molecule. The Hamiltonian
of the optically active molecule considered as quantum mechanical system can
be expressed from the potential energy surfaces calculated above and the quantum
mechanical model for the orientational and vibrational kinetic energy operators.
The latter are obtained from the correspondence principle applied to the classical
formulation described in Appendix A. Disregarding the dynamical couplings between
the orientational and vibrational degrees of freedom the Hamiltonian can thus be
written as

\begin{equation}
H_{\text{a}}^{\text{eff}} = H_{\text{orient}}^{\text{eff}} + H_{\text{vib}}^{\text{eff}};
\label{eq8}
\end{equation}

\noindent where $H_{\text{orient}}^{\text{eff}}$ and $H_{\text{vib}}^{\text{eff}}$
are, respectively, the renormalized Hamiltonians associated with the
orientation and the vibration modes which account for the couplings
(\textbf{r}$_{0}^{\text{e}}$,$\mathbf{\Omega}$)
and (\textit{r}$_{0}^{\text{e}}$,$Q_2$).

\subsection{The orientation modes}

Following expressions derived in Appendix A, the
quantum mechanical renormalized Hamiltonian $H_{\text{orient}}^{\text{eff}}$ for
the non-vibrating molecule can be written as

\begin{equation}
H_{\text{orient}}^{\text{eff}} = T_{\text{rot}}^{\text{eff}} +
V_{\text{M}}(\mathbf{r}_{0}^{\text{e}},\mathbf{\Omega}); \label{eq9}
\end{equation}

\noindent where $T_{\text{rot}}^{\text{eff}}$ is
the effective rotational kinectic operator

\begin{eqnarray}
T_{\text{rot}}^{\text{eff}} = -B^{\text{eff}}\left\{ \frac{\partial ^2}
\partial \theta ^2} + \cot \theta \frac
\partial {\partial \theta } + \frac 1{\sin ^2\theta }\frac{\partial ^2}
\partial \varphi ^2} + \left( \cot ^2\theta +
\frac {C}{B^{\text{eff}}}\right) \frac{\partial ^
}{\partial \chi ^2} - 2\frac{\cot \theta }{\sin \theta }\frac{\partial ^2}
\partial \varphi \partial \chi }\right\}. \label{eq10}
\end{eqnarray}

In this equation $B^{\text{eff}}$ = $\hbar^2/2I_{\text{B}}^{\text{eff}}$ is
the effective rotational constant connected to the rotation-c.m. translation motion of the
C$_3$ symmetry axis for the trapped molecule and $C$ its rotational constant of
the spinning motion around this axis for the gas-phase molecule
$C$ = 6.309 cm$^{-1}$. The value of $B^{\text{eff}}$ is equal
to 7.419 cm$^{-1}$ which is less than the gas-phase value of 9.941 cm$^{-1}$.

Furthermore, since the barrier height associated with the orientation-c.m. translation
motion is always $\lesssim$ 16 cm$^{-1}$, the potential
$V_{\text{M}}(\mathbf{r}_{0}^{\text{e}},\mathbf{\Omega})$ can be neglected.
The renormalized Hamiltonian becomes
$H_{\text{orient}}^{\text{eff}} \simeq T_{\text{rot}}^{\text{eff}}$.

The eigensolutions of this Hamiltonian are then
$E_{JMK}=B^{\text{eff}}J(J+1)+(C-B^{\text{eff}})K^2$ for the eigenvalues
(energy levels) and $|JMK\rangle$ for the eigenvectors where $J$, $M$, and $K$
are the usual quantum numbers connected to the free rotation motion of the
ammonia molecule with the conditions $-J\leq M \leq +J$ and $-J\leq K \leq +J$.
Note that these levels have (2$J$+1)-fold degeneracy on the quantum
number $M$ and twofold degeneracy ($\pm$) on the $K$ one, except when $K=0$ \cite{chtals55}.
The corresponding line structure of the far infrared (FIR) bar-spectrum is then
identical to that of NH$_3$ gas-phase one.

\subsection{The $\nu_2$ vibration-inversion mode}

\subsubsection{Gas-phase NH$_3$}

In the vibration-inversion mode of the NH$_3$ molecule, the nitrogen atom can pass
from one side of the plane defined by the three hydrogen atoms to the opposite one.
This mode is generally described, to a good approximation, as one-dimension
tunneling motion of one particle of mass equal
to the reduced mass $\mu_2$ of NH$_3$ (2.563 g.mol$^{-1}$) moving in a symmetrical
double-well potential function $V_{\nu_2}(Q_2)$ (left (L) and right (R) wells) as
shown in Figure 4a, with a high but finite hindering barrier
of 2047 cm$^{-1}$ \cite{jdsjai62}, such as

\begin{equation}
V_{\nu_2}(Q_2) = \frac 12 k_2 Q_{2}^2 + A_2 (\text{e}^{-a_{2}Q_{2}^2} - 1); \label{eq11}
\end{equation}

\noindent where $k_2$ = 58306 cm$^{-1} \ldotp \text{\AA}^{-2}$, $A_2$ = 12469 cm$^{-1}$ and
$a_{2}$ = 4.8224 $\text{\AA}^{-2}$.

The solution of such a problem leads to vibrational energy levels split into
doublets, ($+$) and ($-$), due to inversion. The vibration-inversion wavefunctions
are then described by $|v_{2}^{(\alpha)}\rangle$ ($\alpha = +,-$), instead of
$|v_2\rangle$, and are expressed as linear combinations (symmetric and
antisymmetric) of the vibrational wavefunctions $|v_{2}^{(\text{L})}\rangle$ and
$|v_{2}^{(\text{R})}\rangle$ associated with the left and right wells, respectively.
Moreover, the only permitted transitions are $(+) \longleftrightarrow (-)$ as shown
in Fig. 4a.

Note that the potential parameters given above were adjusted to give the fundamental
frequency and the level splittings corresponding to the experimental
values, that is, 950 cm$^{-1}$ for the frequency
and 0.8 cm$^{-1}$ for $v_2$ = 0 and 35.8 cm$^{-1}$ for $v_2$ = 1 level splittings.

\subsubsection{Trapped NH$_3$}

When the NH$_3$ molecule is trapped in C$_{60}$ cage, the gas-phase double-well
potential function is modified. In the calculations, the inversion-c.m. translation
potential energy $V_{\text{M}}(r_0,Q_2)$ (see Eq. (\ref{eq7})) is numerically
added to the potential function $V_{\nu_2}(Q_2)$ of Eq. (\ref{eq11}).
The new double-well potential function is shown in Fig. 4b as a function of
the inversion angular coordinate. Its symmetry is preserved but the hindering
barrier is increased by 404 cm$^{-1}$.
Moreover, from expressions given in Appendix A, the effective reduced
mass is equal to $\mu_{2}^{\text{eff}}$ = 6.569 g.mol$^{-1}$ leading to
an increase of 156$\%$ with respect to the gas-phase value of 2.563 g.mol$^{-1}$.

The corresponding Schr\"{o}dinger equation was solved numerically using
a discrete variable representation method \cite{jcliphjvl85}. The computed level
scheme is given in Fig. 4b.

In gas-phase, the $\nu_2$ mode is measured at 931.7 cm$^{-1}$ and 968.3 cm$^{-1}$.
In C$_{60}$, it is shifted by 47.0 cm$^{-1}$ (as given below) and because the
double-well inversion potential is modified with an increase in the hindering
barrier height and the effective reduced mass as discussed above, this
mode is quasi-degenerate and is calculated at 751.6 cm$^{-1}$, that is,
(704.6+47.0) cm$^{-1}$ (see Fig. 4b). This result is to be related to an
increase in the tunneling
time for both the fundamental and the first excited levels in C$_{60}$, since
the values are equal to 20.85 ps and 0.47 ps in gas-phase compared
to 55594 ps and 333.56 ps in C$_{60}$, respectively. We can infer from these
results that in C$_{60}$, localization of the N atom occurs in one of the
wells of the potential function and that the $\nu_2$ mode is quasi-degenerate.

\begin{figure}[h!]
\begin{center}
\includegraphics[width=7cm]{Figure4a.eps}
\includegraphics[width=7cm]{Figure4b.eps}
\end{center}
\caption{Vibration-inversion double-well potential energy of NH$_3$ as a function of
the angle between the N$-$H bond and the plane of the H atoms. The obtained energy
levels and the allowed fundamental transitions are also shown (a) for the
gas-phase and (b) for the trapped one in the fullerene nano-cage.}
\label{FIG. 4}
\end{figure}

\subsection{The vibrational frequency shifts}

The second term of Eq. (\ref{eq8}) is the renormalized vibrational Hamiltonian
which accounts for the vibrational dependent part
$W_{\text{M}}(\{Q\})$ and inversion-c.m. translation part $V_{\text{M}}(r_0,Q_2)$ of the
interaction potential energy

\begin{equation}
H_{\text{vib}}^{\text{eff}} = H_{\text{vib}} + W_{\text{M}}(\{Q\}) + V_{\text{M}}(r_0,Q_2);
\label{eq12}
\end{equation}

\noindent where $H_\text{vib}$ characterizes the gas-phase molecular vibrational
Hamiltonian written as

\begin{equation}
H_{\text{vib}} = \sum\limits_\nu \frac{P_\nu ^2}{2\mu _\nu }+\frac 12 \sum\limits_\nu k_\nu
Q_\nu ^2+ \bigtriangleup V_{\text{vib}}^{\text{anh}}(\{Q\}). \label{eq13}
\end{equation}

In this expression $\mu _\nu$ and $k _\nu$ are, respectively, the reduced mass and
the harmonic force constant connected to the $\nu$th vibrational mode, with associated
normal coordinate $Q_\nu$ and conjugate momentum $P_\nu$, and
$\bigtriangleup V_{\text{vib}}^{\text{anh}}(\{Q\})$ represents the anharmonic part of
the internal potential function. The eigenelements of this Hamiltonian $H_{\text{vib}}$
have previously been obtained from the \textit{ab initio} calculations of
Martin \textit{et al.} \cite{jmlmtjlprt92}.

Let $E_{v_\nu}$ and $|...v_\nu...\rangle$ be the eigenvalues and eigenvectors of
the \textit{v}th level associated with the $\nu$th vibrational mode.

Thus, a first order perturbation treatment allows us to determine the frequency
shift $\Delta \omega _\nu$ associated with the vibrational fundamental transition
$\mid...0_\nu ...\rangle \rightarrow \mid ...1_\nu...\rangle $ for
the trapped molecule from the equation

\begin{equation}
\Delta \omega _\nu = \hslash^{-1}\left[ \langle ...1_\nu...\mid W_{\text{M}}(\{Q\})
\mid ...1_\nu...\rangle -
\langle ...0_\nu ...\mid W_{\text{M}}(\{Q\})\mid ...0_\nu...\rangle \right]; \label{eq14}
\end{equation}

\noindent in which all other modes remain in their fundamental states
and $\hslash = h/2\pi$. The computed values are given in Table 3. We note that each
of the vibrational modes are blue shifted with splitting of the degenerate
levels of $\nu_3$ mode (strong) and $\nu_4$ mode (weak) modelled from static effect.
Whereas for degenerate $\nu_2$ mode the blue shift of 47.0 cm$^{-1}$ from static
effect, is modified to red shift when dynamic effect on 1D tunneling motion
of the N atom is considered as discussed above.

\begin{table}[h!]
\caption{Vibrational frequency shifts (cm$^{-1}$) of NH$_3$ molecule
trapped in the C$_{60}$ cage, together with the gas-phase frequencies. }
\begin{tabular}{ccc}
\hline\hline
\begin{tabular}{c} Vibrational modes \end{tabular} & \begin{tabular}{c} Frequency
shifts (cm$^{-1}$) \end{tabular} & \begin{tabular}{c} Gas-phase  frequencies (cm$^{-1}$)
\end{tabular} \\
\hline
$\nu_1$ & \multicolumn{1}{c}{133.0} & \multicolumn{1}{c}{3337} \\
$\nu_2$ & \multicolumn{1}{c}{47.0} & \multicolumn{1}{c}{950} \\
$\nu_{3a}$ & \multicolumn{1}{c}{171.0} & \multicolumn{1}{c}{3448} \\
$\nu_{3b}$ & \multicolumn{1}{c}{69.6} & \multicolumn{1}{c}{} \\
$\nu_{4a}$ & \multicolumn{1}{c}{47.1} & \multicolumn{1}{c}{1627} \\
$\nu_{4b}$ & \multicolumn{1}{c}{34.0} & \multicolumn{1}{c}{} \\
\hline
\end{tabular}
\label{Table 2}
\end{table}

\section{Infrared absorption spectra}

\subsection{General}

The infrared absorption coefficient for the optically active ammonia molecule
trapped in a C$_{60}$ fullerene molecule at temperature \textit{T}, is defined
as the real part of the spectral density, \textit{i.e.}, the Fourier transform of the
time-dependent autocorrelation function $\phi (t)$

\begin{eqnarray}
I\left( \omega \right) &=&\frac{4\pi \mathcal{N}\omega }{3hc}\mathrm{Re
\int\limits_0^\infty dte^{i\omega t}\phi (t), \nonumber \\
\phi (t) &=&\text{Tr}\left[ \rho (0)\mathbf{\mu }_A(0)\mathbf{\mu
_A(t)\right], \label{eq15}
\end{eqnarray}

\noindent where $\omega $ is a running frequency variable
and $c$ the speed of light, $\rho (0)$
characterizes the initial canonical density matrix operator of the system, and
$\mu _A$ is the molecular dipole moment operator defined in the absolute frame
(see Appendix B). The trace (Tr) operation means an average over the initial
conditions at time $t=0$ and over all the possible evolutions of the system
between times $0$ and $t$.

However, since this work is devoted to constructing the infrared bar-spectra
of the trapped ammonia molecule, we assume the \textit{initial chaos hypothesis} for the
density matrix to be valid. This allows us to write the optical wave functions
$\mid v_\nu JMK\rangle$ as products of the renormalized vibrational and
orientational ones $\mid v_\nu \rangle \otimes |JMK\rangle_{v_{\nu}}$.
It must be noticed that the orientational states parametrically depend on the
vibrational states through the moments of inertia. This dependence will be
ignored below.

\subsection{Near infrared bar-spectrum}

Within this framework, the infrared bar-spectrum connected to the $\nu$th
fundamental vibrational mode transition
$\mid 0_\nu\rangle \longrightarrow \mid 1_\nu\rangle $
(all other modes being in their fundamental states) is written in the Lorentzian form as

\begin{eqnarray}
I_\nu (\omega ) = \frac{8\pi ^2\mathcal{N}}{3hc}\omega
\left| \left\langle 0_\nu \left| Q_\nu \right|
1_\nu\right\rangle \right| ^2
\sum\limits_{ {JMK}{J^{\prime }M^{\prime }K^{\prime }}}\frac
e^{-\beta E_{0JMK}}-e^{-\beta E_{1J^{\prime }M^{\prime }K^{\prime}}}}Z
\left| \left\langle JMK\left| \frac{\partial \mathbf{\mu }_
}{\partial Q_\nu }\right| J^{\prime }M^{\prime }K^{\prime
\right\rangle \right| ^2 \nonumber\\
\times \delta \left(\omega -\omega _{0JMK\rightarrow
1J^{\prime }M^{\prime }K^{\prime }}-\Delta \omega _{\nu}\right).
\label{eq16}
\end{eqnarray}

\noindent In this equation the $\left\langle ...\right\rangle $ brackets
refer to vibrational transition elements of the normal coordinate $Q_\nu $
and to the orientational transition elements of the first derivative of the
molecular dipole moment with respect to this coordinate. The numerical values
connected to all the NH$_3$ modes are given in Appendix B. Moreover, $Z$ defines the
orientation canonical partition function associated with the fundamental
vibrational level at temperature $T$, $\beta=(k_{B}T)^{-1}$ where $k_B$ is the
Boltzmann's constant, $\delta$ is the Dirac's function, and $\omega
_{0JMK\rightarrow 1J^{\prime }M^{\prime }K^{\prime }}=\hbar
^{-1}\left( E_{1J^{\prime }M^{\prime }K^{\prime }}-E_{0JMK}\right)$
is the transition frequency. Note that we use below the usual spectroscopic
nomenclature to designate the line transitions connected to the rotational
selection rules, that is, Q($J_K$) for $\Delta J = 0$, R($J_K$) for $\Delta J = +1$,
and P($J_K$) for $\Delta J = -1$.

\subsection{$\nu_2$ mode results and discussion}

In Figure 5, are given the near infrared bar-spectra of NH$_3$ trapped in fullerene
in the $\nu_2$ mode absorption region for three different temperatures, 10 K, 30 K
and 100 K. The spectra are calculated by taking into account the frequency shift and
allowed rovibrational transition elements following selection rules.

For the latter, as the first derivative of the molecular dipole moment with respect
to the $Q_2$ coordinate does not depend on the rotational spinning angle $\chi$,
$\Delta K= 0$ for the $K$ quantum number.

The infrared bar-spectrum consists of a Q line located at the pure
vibration-inversion frequency 751.6 cm$^{-1}$, and a set of
rotation-vibration-inversion lines P and R 14.8 cm$^{-1}$ apart, corresponding
to 2$B^{\text{eff}}$, on both sides of the Q branch.

\begin{figure}[h!]
\begin{center}
\includegraphics[width=15cm]{Figure5a.eps}
\end{center}
\begin{center}
\includegraphics[width=15cm]{Figure5b.eps}
\end{center}
\begin{center}
\includegraphics[width=15cm]{Figure5c.eps}
\end{center}
\caption{Near infrared bar-spectrum in the $\nu_2$ vibration-inversion frequency region for
NH$_3$ trapped in the fullerene nano-cage at temperatures 10, 30 and 100 K.}
\label{FIG. 5}
\end{figure}

At $T$ = 10 K, only the lowest few rotational states are populated. As a
result, the bar-spectrum shown in Figure 5a exhibits only three lines:
\textit{i}) a Q line at 751.6 cm$^{-1}$ which is in fact a superposition of
the three transitions Q($1_1$), Q($2_1$) and Q($2_2$),
\textit{ii}) a R($0_0$) line at 766.4 cm$^{-1}$
and \textit{iii}) a R line, superposition of R($1_0$) and R($1_1$) lines
at 781.3 cm$^{-1}$. There is only one P($1_0$) line
located at 736.8 cm$^{-1}$.

As the temperature increases, more and more rotational states
are populated and as result, more lines appear in the spectra with a different
distribution in terms of intensities. At 30 K, the Q line is the most intense
consisting of a superposition of Q($J_K$) lines for $J$ = 1 to 4
and $1\leq K \leq J$. Moreover, the intensity of lines R($0_0$)
and R($1_K$) decreases leading to the increase of the R($2_K$) line which is located
at 796.1 cm$^{-1}$. More lines appear in the R and P branches
(R($3_K$) and P($2_K$)).

At 100 K, the Q branch consists of a superposition of 6 lines for $J$ varying
from 1 to 6. The R and P branches are broader with lines of weaker intensities
than those at 10 and 30 K in particular for the maxima.

The modification of the spectra near $\nu_2$ with temperature shows that the spectral
signature of the trapped molecule can be used to sense the environmental
temperature of fullerene C$_{60}$. As the coupling with the vibrational modes
of C$_{60}$ is weak, the signal to noise ratio may be high enough for all the
lines to be observed. In our opinion this result can be extrapolated to other
types of simple molecules trapped in C$_{60}$.

Finally, the dynamical coupling of the vibration-rotation states of the NH$_3$
molecule with the vibrational modes of fullerene C$_{60}$ involves well separated
states for the latter. As a result, this coupling is weak and the energy from
local translational modes cannot redistribute itself in the thermal bath of C$_{60}$.
The linewidths are then negligible.

In conclusion, this work shows that the simulation of the spectra of NH$_3$ can
be used to probe the temperature of the media in which fullerene is observed.
This result can probably be extended to other molecules trapped in fullerene, in
particular when the coupling with the vibration modes of fullerene is weak, because
signal to noise ratio should then be high enough for the rotation lines to be observed.

\newpage

\textbf{ACKNOWLEDGEMENTS}

The authors thank the Regional Councels of Franche-Comt\'{e} and Ile de France,
the General Councel of the 78$^{\text{th}}$ district of Yvelines, the DGCIS
from Ministry of Industry in France, the MOVEO cluster for financing the FUI
project MEMOIRE under the label of cluster MOVEO, CNRS for financing support under
the Interdisciplinary Research Program of the INSU \textquotedblleft Planetary
Environments and Life Origins (EPOV)\textquotedblright, and also Camille Lakhlifi
for figures help.


text:
\section{INTRODUCTION}

In anomaly detection, we seek to differentiate normal behavior from anomalous behavior - essentially anything that deviates away from the nominal model. This task is easy if our confidence is high in the nominal model but as uncertainty rises, distinguishing normal from anomalous becomes more challenging. From a control theory perspective, one of the most fundamental sources of uncertainty is measurement noise and in this work we present an intuitive, but to-date-unexplored, idea of combining traditional model-based detection schemes with low-pass filtering to reduce the impact of noise on our ability to detect anomalies in control systems. Beyond demonstrating that this approach enables better detection, the novelty lies in characterizing this performance boost analytically, in characterizing how this effects the impact an attacker can have on the system, and in understanding what new kind of attacks are possible when filtering is involved.

There is a deep literature from Fault Detection that leverages model-based detectors to identify the occurrence of faults \cite{Patton_1}. More recently these tools have been reused in identifying the presence of attacks in control systems \cite{Cardenas}\nocite{mo2014detecting,Mo_3,Guo2016,milovsevic2018quantifying,umsonst2018anomaly,RuthsACC2018_Hidden,RuthsACC2018_Windowed,RuthsACC2018_Comparison}-\cite{RuthsACC2019_nongaussian}. Attacks, in particular their strategic and exploitative nature, offer new challenges to both control theory and anomaly detection. While many research groups have rallied behind the banner of analyzing attacks in control systems and the challenges they raise, in what follows we describe the main representative contributions of these groups that relate to characterizing the tuning and performance of model-based detectors.

While all papers \cite{Cardenas}-\cite{RuthsACC2019_nongaussian} investigate the use of model-based detectors, and \cite{Cardenas}-\cite{RuthsACC2018_Comparison} consider various types of attacks, it has only been recently that the impact of these attacks has been evaluated \cite{milovsevic2018quantifying}-\cite{RuthsACC2018_Comparison}. Of these \cite{Mo_3,RuthsACC2018_Hidden,RuthsACC2018_Comparison} use the reachable set to characterize attack impact and \cite{milovsevic2018quantifying,umsonst2018anomaly,RuthsACC2018_Windowed} use various norms of system state or the estimation error covariance \cite{Guo2016}. To a large extent this work on attacker impact has been facilitated by analytic methods to tune model-based detectors, as provided in \cite{Mo_3,umsonst2018anomaly,RuthsACC2018_Hidden,RuthsACC2018_Windowed,RuthsACC2018_Comparison,RuthsACC2019_nongaussian} and for work to define various types of worst-case attacks \cite{mo2014detecting}, \cite{Guo2016}-\cite{RuthsACC2018_Comparison}. In all of these cases Gaussian noise(s) are assumed on the measurements of the system (and possibly the system state), except in \cite{RuthsACC2019_nongaussian} which tunes a detector for generalized noise distributions.

While much has been done with model-based detectors, no study has integrated a low-pass filter to attenuate the effect of the noise on the attacker's capabilities to change the system state. The recent work to to tune classical detectors to desired levels of performance (desired false alarm rate) and to characterize the attacker's impact on the system state (given the attacker's desire to remain stealthy to classical detectors) positions us well to now add another layer, i.e., filtering, to the detection scheme. In Section \ref{sec:residual_detection}, we characterize the improvement gained by adding the filter (in terms of sensitivity to attacks and reduction in attacker capabilities) traded off with the new attack definitions that are stealthy to the introduced \textit{low-pass chi-squared detector}.

This notion of filtering is not exclusive to retrofit conventional residual-based detectors, but generalizes in compelling ways to other types of model-based observers. In particular, in Section \ref{sec:discontinuous_observer}, we review how the discontinuous observer, an estimator that uses a sliding mode to yield finite-time convergence, produces a discontinuous term that can be filtered to produce an estimate of the disturbance that enters a system \cite{robustobs}. In this case, we reframe the disturbance estimation problem as a method to approximate the anomaly (e.g., failure and/or attack), which provides an attractive alternative to residual-based detectors. In this case the disturbance estimate not only provides a way to detect the presence of an anomaly, but also find the form of the anomaly so that advanced mitigation strategies can be employed - such as using the estimated anomaly as feedback to avoid the estimate to deviate away from the true state value.

In Section \ref{sec:system} we introduce the system and our notation as well as the attack context; in Section \ref{sec:example} we provide an example that combines and compares both of these tools.

\section{SYSTEM \& CONTEXT} \label{sec:system}
We consider a continuous, linear, time-invariant control system
 \begin{equation} \label{eq:dLTIsystem}
 \left\{\begin{aligned}
 	\dot{x}(t) &= Ax(t) + Bu(t),\\
 	y(t) &= Cx(t) + \eta(t), \\
 \end{aligned}\right.
 \end{equation}
in which $x(t)\in\mathbb{R}^n$ is the system state, $u(t)\in\mathbb{R}^m$ is the input, A$\in\mathbb{R}^{n\times n}$ is the system matrix, $B\in\mathbb{R}^{n\times m}$ is the input matrix, $y\in\mathbb{R}^p$ is the system output, and $\eta(t) \in \mathbb{R}^p$ is zero mean Gaussian additive measurement noise with covariance $R$, i.e., $\eta(t) \sim \mathcal{N}(0,R)$. The pair $(A,C)$ is detectable and the pair $(A,B)$ is stablizable.

We use the output $y$ and an estimator to produce an estimate of the state, $\hat{x}(t)$. In this work we consider two different choices of estimators. We use this estimate for to construct the feedback law,
\begin{equation}
    u(t)=K\hat{x}(t).
\end{equation}
The process measurements are potentially corrupted by an attacker, who is able to read and arbitrarily change the sensor measurements. We model this as an additive attack $\delta(t) \in \mathbb{R}^p$ such that $y(t)$ is changed to $\bar{y}(t)$,
 \begin{equation}\label{outputb}
 	\bar{y}(t) = y(t) + \delta(t) =  Cx(t) + \eta(t) + \delta(t).
 \end{equation}

\subsection{Attacker Capabilities}
In this work we assume that the attacker has access to all system information including, e.g., dynamics, states, estimator states, and detector parameters. In particular, the attacker can view (disclosure) and edit (disruption) sensor measurements.

\section{RESIDUAL-BASED ANOMALY DETECTION} \label{sec:residual_detection}

The residual quantifies the difference between what we receive from the sensor, $\bar{y}(t)$ and what we expect based on an estimator, $\hat{y}(t)=C\hat{x}(t)$,
\begin{equation}
	r(t) = \bar{y}(t) - C\hat{x}(t),
\end{equation}
Anomaly detectors use this quantity to make real-time choices in the state of the control system they monitor.

To build an estimate of the system state we use a Luenberger observer with observer gain $L \in \mathbb{R}^{n \times p}$,
\begin{equation}
\begin{aligned}
	\dot{\hat{x}}(t) =& A\hat{x}(t) + Bu(t) + L(\bar{y}(t) - C\hat{x}(t)).
\end{aligned}
\end{equation}
Armed with this estimator, we can express the difference between the state and expected state as the estimation error $e(t)=x(t)-\hat{x}(t)$, which leads to the following coupled closed-loop equations,
\begin{equation} \label{eq:esterror}
\left\{ \begin{aligned}
    \dot{x}(t) &= \big( A + BK \big) x(t) - BKe(t), \\
	\dot{e}(t) &= \big( A - LC \big) e(t) - L\eta(t) - L \delta(t), \\
	r(t) &= Ce(t) + \eta(t) + \delta(t).
\end{aligned}\right.
\end{equation}
The steady state Kalman filter is the standard choice of estimator used when employing residual based detectors. This allows us to calculate the steady state covariance of the estimation error from the stochastic dynamics of $e(t)$ using the following Riccati equation,
\begin{equation}
\dot{P}=0=(A-LC) P + P(A-LC)'+LR L'.
\end{equation}

Because the detection process is inherently a sampled approach we discretize \eqref{eq:esterror} with uniform step size $\tau$ such that $\xi(k\tau)=\xi_k$,
\begin{equation} \label{eq:discrete}
\left\{ \begin{aligned}
    x_{k+1} &= F x_k + G e_k, \\
	e_{k+1} &= H e_k - L_d(\eta_k + \delta_k), \\
	r_k &= Ce_k + \eta_k + \delta_k.
\end{aligned}\right.
\end{equation}
Residual-based detection relies on quantifying the distribution of the residual under nominal operation, i.e., without attacks/anomalies. These statistics form a one-sided hypothesis test that is either accepted (no anomalies) or rejected (alarm raised). The discrete covariance of $\eta_k$ is $\Sigma_\eta = \tau R$ and the discrete covariance of the estimation error is $\Sigma_e = \tau P$.
In the absence of attack the distribution of the residual follows a zero-mean Gaussian distribution with the covariance \cite{Carlos_Justin2},
\begin{equation} \label{eq:residual covariance}
	\Sigma_r = E[r_kr'_k] = C\Sigma_e C'+\Sigma_\eta.
\end{equation}
In the following subsections we will review the conventional chi-squared detector and subsequently extend this detector to enhance its performance.

\subsection{Chi-Squared Detector}
The role of the detector is to create a non-negative scalar-valued random variable from the residual that can be easily compared with a threshold. One of the most widely used approaches to do this is given by the chi-squared detector which introduces a distance measure
\begin{equation} \label{eq:distancemeasure}
	z(t)= r'(t)\Sigma_r^{-1}r(t).
\end{equation}
Since the residual is normally distributed, the distance measure, as a sum of squared Gaussian random variables, is chi-squared distributed (with $p$ degrees of freedom, where $p$ is the number of the sensors). By normalizing the residual covariance, we eliminate system dependence. The detection procedure for the chi-squared detector is summarized as follows. For given a threshold $\alpha \in \mathbb{R}_{>0}$ and the distance measure in \eqref{eq:distancemeasure}
\begin{equation} \label{detector}
\left\{\begin{aligned}
	z(t) \leq \alpha &\quad\longrightarrow\quad \text{no alarm}, \\
	z(t) > \alpha &\quad\longrightarrow\quad \text{alarm}.
\end{aligned}\right.
\end{equation}
Since the chi-squared distribution has $p$ degrees of freedom the distance measure has mean $E[z(t)]=p$ and covariance $E[z^2]=2p$ during nominal (no anomaly) operation. In the presence of attacks/anomalies the distance measure will, in general, not fall according to a chi-squared distribution with mean value $p$ and covariance $2p$. As indicated in \eqref{detector}, the detector is required to make a decision to raise or not raise an alarm at every time instant. Therefore, it is impractical to build a sample distribution of ``current behavior'' to decide if the distance measure was still distributed in chi-squared fashion - this would take too long. Instead, the alarms, and more specifically, the rate of alarms can be used as a metric for how deviated the system is from nominal behavior.

Here \textit{alarm rate} is defined as the rate of generation of alarms by the detector, which empirically is the fraction of time instants in which an alarm was raised. We expect alarms to be raised regularly, even under normal operation, because of the infinite support of the measurement noise. We can predict the false alarm rate $\mathcal{A}$ - the alarm rate under normal operation - from the distance measure distribution, $\textbf{Pr}(z > \alpha)= \mathcal{A}$. More importantly, we can set the false alarm rate to a desired value $\mathcal{A}^*$ by selecting the threshold $\alpha$ appropriately.
\begin{lemma} \cite{Carlos_Justin2}. \label{lem:chisquared_tuning}
Consider a stochastic LTI system under normal operation and a chi-squared detector, with threshold $\alpha \in \mathbb{R}_{>0}$, $r_k \sim \mathcal{N}(0,\Sigma)$. Let $\alpha = \alpha^{*}:=2P^{-1}(1-\mathcal A^{*}, \frac{p}{2})$, where $P^{-1}(\cdot,\cdot)$ denotes the inverse regularized lower incomplete gamma function, then ${\mathcal A}= {\mathcal A}^{*}$.
\end{lemma}
If the attacker wants to remain stealthy (undetected) to the detector, the threshold of the chi-squared detector puts a limit on the distance measure. It, therefore, limits the capabilities of stealthy attacks to change the system behavior. We derive worst-case stealthy attacks by assuming a powerful, knowledgeable attacker that knows the noisy measurement $Cx+\eta$ and also knows the state estimate $\hat{x}$, the attack can compensate for the terms in the residual
\begin{equation}\label{eq:deltabar}
	\delta_k = - \bar{y}_k+C\hat{x}_k+\Sigma_r^{\frac{1}{2}}\bar{\delta}_k,
\end{equation}
such that
\begin{equation}
    z_k=\bar{\delta}'_k\bar{\delta}_k.
\end{equation}
By increasing the norm of $\bar{\delta}_k$ the attacker can increase the impact of the attack up to the point where the attack is identified by the detector; the ``direction'' of the vector $\bar{\delta}_k$ indicates the extent to which the attack impacts each sensor.

In \textit{\underline{zero alarm attacks}} the attacker takes the perspective that raising no alarms is the way to avoid detection, essentially keeping the distance measure below the threshold,
\begin{equation}
    z_k=\bar{\delta}'_k\bar{\delta}_k \leq \alpha,
\end{equation}
and hence $\| \bar{\delta}_k \| \leq \sqrt{\alpha}$.

In \textit{\underline{hidden attacks}}, the attacker realizes that under normal operation alarms are generated at the rate of $\mathcal{A}^*$, hence, it is reasonable to generate an attack sequence that mimics the false alarm rate,
\begin{equation} \label{eq:hiddendef}
	\textbf{Pr}(z_k > \alpha) = \textbf{Pr}\big(\|\bar\delta_k \| > \sqrt{\alpha}\big) = \mathcal{A}^*.
\end{equation}
Hidden attacks are inherently more potent than zero-alarm attacks, especially when keeping in mind that the norm of the attack can be arbitrarily large at the time instants in which alarms are raised.

For more detail on these attacks and ways in which to evaluate the impact of these attacks through state deviation or the induced reachable set see \cite{RuthsACC2018_Windowed,RuthsACC2018_Comparison}.

\subsection{Filtered Chi-Squared Detector}

In this section we introduce a modification to the chi-squared detector to leverage apriori knowledge that the noise is typically composed of substantially higher frequencies than the disturbances or many types of anomalies. By filtering the residual signal, we aim to reduce the covariance of the statistic that forms the distance measure so that attacks are more readily apparent. Exploiting frequency domain information has not been combined with attack detection and here we introduce the approach and characterize its performance.

We filter the residual element-wise using a bank of $p$ identical Butterworth filters (of second order) such that $\xi_i(t)=[\xi_{1i}(t),\xi_{2i}(t)]'$ and
\begin{equation} \label{eq:filter_dynamics}
    \left\{\begin{aligned}
    \dot\xi_i(t) &= \Phi \xi_i(t)+ \Psi r_i(t),\\
    \rho_i(t) &= [1\ 0]\xi_i(t)
    \end{aligned}\right.
\end{equation}
for $i=1,\dots,p$ and where
\begin{equation} \label{eq:filter_matrices}
        \Phi=\begin{bmatrix}0&1\\-\omega_c^2&-\sqrt{2}\omega_c\end{bmatrix},\ \  \Psi=\begin{bmatrix}0\\\omega_c^2\end{bmatrix}.
\end{equation}
and we use the filtered output $\rho(t)$ to construct a new distance measure. To characterize the distribution of this distance measure, first characterize the  covariance of $\rho_i(t)$.
\begin{theorem}
Given a residual signal $r_k\sim \mathcal{N}(0,\Sigma_r)$ filtered through the Butterworth filter in \eqref{eq:filter_dynamics}-\eqref{eq:filter_matrices} with bandwidth $\omega_c$, the discrete-time filter with sampling rate $\tau$ output $\rho_k$ is zero mean Gaussian with covariance,
\begin{equation}\label{eq:filtered_covariance}
    \Sigma_\rho = \frac{\tau \omega_c}{2\sqrt{2}}\Sigma_r.
\end{equation}
\end{theorem}
\vspace{2mm}

\begin{proof}
We first discretize the filter with time step $\tau$
\begin{equation} \label{eq:filter_dynamics_discrete}
    \left\{\begin{aligned}
    \xi_k^{(i)} &= \Phi_d \xi^{(i)}_k+ \Psi_d r^{(i)}_k,\\
    \rho_k^{(i)} &= [1\ 0]\xi^{(i)}_k
    \end{aligned}\right.
\end{equation}
where $\zeta^{(i)}_k = \zeta_i(k\tau)$ and for small enough $\tau$, $\Phi_d=I+\Phi\tau$ and $\Psi_d=\Psi\tau$. To find the covariance of the filtered residual vector $\rho_k$ we collect all filters together,
\begin{equation}\label{eq:combinedfilter}
\dot{\xi}_k = \widetilde{\Phi} \xi_k + \widetilde{\Psi} r_k,
\end{equation}
with $\widetilde{\Phi}=\Phi_d\otimes I_p$ and $\widetilde{\Psi}=\Psi_d\otimes I_p$ and where $\xi_k \in \mathbb{R}^{2p}$, $\widetilde{\Phi} \in \mathbb{R}^{2p\times2p}$ and $\widetilde{\Psi} \in \mathbb{R}^{2p\times p}$. We can define the covariance of $\xi_k$ as $\mathcal{P}=[\mathcal{P}_{i,j}]$, $i=1,\dots,p$ and $j=1,\dots,p$, where $\mathcal{P}_{i,j} \in \mathbb{R}^{2 \times 2}$. This total covariance is supplied by the discrete-time Riccati equation (in steady state),
\begin{equation}\label{eq:combinedfilter}
    0 = \widetilde{\Phi}\mathcal{P}\widetilde{\Phi}'-\mathcal{P} + \widetilde{\Psi}\Sigma_r\widetilde{\Psi}',
\end{equation}
and following some algebra manipulation reveals that we can split this into Riccati equations on each subblock
\begin{equation}\label{eq:splittedfilter}
\Phi_d \mathcal{P}_{i,j} \Phi_d' - \mathcal{P}_{i,j} + \Psi_d[\Sigma_r]_{i,j} \Psi_d'=0.
\end{equation}
From \eqref{eq:splittedfilter} we see that $\mathcal{P}_{i,j}=\mathcal{P}_{j,i}$ because $[\Sigma_r]_{i,j}=[\Sigma_r]_{j,i}$ and in addition based on \eqref{eq:combinedfilter}, we know that $\mathcal{P}_{i,j}=\mathcal{P}_{j,i}'$, therefore matrix $\mathcal{P}_{i,j}$ is a symmetric square matrix. The covariance of $\rho(t)$ can be computed element-wise using the structure of $\Phi_d$ and $\Psi_d$ leaving,
\begin{align}
    [\Sigma_\rho]&_{i,j} = [\mathcal{P}_{i,j}]_{11}\\ &= -\tau\omega_c {\tiny \frac{(\tau\omega_c)^2-\sqrt{2}\tau\omega_c+2}{(\tau\omega_c)^3-3\sqrt{2}(\tau\omega_c)^2+8\tau\omega_c-4\sqrt{2}}} \Sigma_r \nonumber
\end{align}
If the sampling time is taken small enough the higher order terms vanish leaving
\begin{equation}
    [\Sigma_{\rho}]_{i,j}=\frac{\tau\omega_c}{2\sqrt{2}} [\Sigma_r]_{i,j},
\end{equation}
which leads directly to \eqref{eq:filtered_covariance}.
\end{proof}

\begin{remark}
Note that $\tau$ is the sampling time of the system, which highlights that the filtered covariance is expected to be significantly smaller than the covariance of the original residual. By reducing the covariance of the nominal behavior using this filtering technique the aim is to make it easier to distinguish attacks and anomalies (especially those with low frequency components).
\end{remark}

Using this covariance, we can create a new normalized distance measure
\begin{equation}
    z_k=\rho'_k \Sigma_\rho^{-1} \rho_k = \frac{2\sqrt{2}}{\tau\omega_c} \rho'_k \Sigma_r^{-1} \rho_k.
\end{equation}
Because $\rho_k$ is a Gaussian random variable, $z(t)$ is also chi-squared. In fact, because of the normalization by the covariance, this functions exactly like the conventional chi-squared detector (including using Lemma \ref{lem:chisquared_tuning} to select the threshold). The major difference lies in the scale of the normalizing covariance; here being smaller provides an advantage to identify smaller attacks. To execute these attacks, attackers will also need to know the cut-off frequency of the filter. Analogous to \eqref{eq:deltabar}, we define
\begin{equation}
    \delta_k=-y_k+C\hat{x}_k+\left(\frac{\tau\omega_c}{2\sqrt{2}}\Sigma_r\right)^{\frac{1}{2}}\ \bar{\delta}_k,
\end{equation}
where $\bar{\delta}_k$ is an independent random variable that attackers may select to define different classes of stealthy attacks.

As before the \textit{\underline{zero alarm attack}} is generated in a way so that no alarm is generated by the detector. The precise characterization of the sequence $\bar{\delta}_k$ is a topic for future work, however, intuitively the fact that the attack signal is also filtered by the low-pass filter, it follows that the frequency content (as quantified by a Fourier Transform) and not just the amplitude (norm) of the attack vector plays a role in quantifying zero alarm attacks. To start, note that zero alarm attacks of the conventional chi-squared detector will also be zero alarm attacks of the filtered chi-squared detector. On top of this, there is an opportunity to inject high frequency content into the attack, which would be attenuated by the action of filter. If the filtered chi-squared used an \textit{ideal} low pass filter, it would be possible to inject any signal that had frequency content higher than the cut-off frequency of the filter. The second order filter requires some adjustment to balance the small contributions from the stopband frequencies with that of the passband frequencies.

Similar to the unfiltered case, we employ the zero alarm attack constraint in a probabilistic fashion for the corresponding \textit{hidden attack}. In this case $\bar{\delta}_k$ is selected such that the Fourier frequency content is now bounded in probability. Again the full characterization of these attacks are for future work.

\section{DISCONTINUOUS OBSERVER} \label{sec:discontinuous_observer}
In this section, we present an anomaly detector for \textit{smooth} attacks $\delta(t)$, based on a robust observer and a low pass filter, which is applicable  for a particular class of second order systems, i.e., systems of the form (\ref{eq:dLTIsystem}) with $x=[x_{1}\ x_{2}]^{T}$, $A\in\mathbb{R}^{2\times2}$, $B\in\mathbb{R}^{2}$, and $C\in\mathbb{R}^{1\times2}$. Note, however, that any physical plant that can be partitioned into a collection of second order systems can be captured by an extended version of the presented methods.

The observer is given by  \cite{robustobs}
\begin{align}
\dot{\hat{x}}&=A\hat{x}+Bu+\Gamma e+Bc_{3}\hbox{sign}(\bar{y}-\hat{y}),\label{obsi1}\\
\hat{y}&=C\hat{x},\label{obsi2}
\end{align}
where  $\hat{x}=[\hat{x}_{1}\ \hat{x}_{2}]^{T}$, $e=[(\bar{y}-\hat{y})\ (x_{2}-\hat{x}_{2})]^{T}$,  $\hat{x}_{i}\in\mathbb{R}$, and the diagonal matrix $\Gamma=\mbox{diag}(c_{1},c_{2})$, with $c_{1},c_{2}>0$. Finally, the gain $c_{3}$ of the discontinuous term can be chosen as described in \cite{robustobs}.

An interesting feature of this observer is that, by filtering the discontinuous term, it is possible to identify disturbances/anomalies. This is shown as follows. First, assume that the measurement is free of noise, i.e., $\eta(t)=0$ in (\ref{outputb}), and define the observation errors
\begin{equation}
e_{1}=\bar{y}-\hat{y},\  \  \  \mbox{and}\ \ \ \ e_{2}=x_{2}-\hat{x}_{2}.
\end{equation}
Then, by considering the output vector $C=[1\ 0]$, the corresponding observation error dynamics is described by
\begin{align}
\dot{e}_{1}&=e_{2}-c_{1}e_{1}+\dot{\delta},\label{ob1}\\
\dot{e}_{2}&=-(a+c_{2})e_{1}-be_{2}-c_{3}\mbox{sign}(e_{1})+a\delta.\label{ob2}
\end{align}
When the system reaches the discontinuous surface, it follows that $e_{1}=\dot{e}_{1}=0$. Substituting this into (\ref{ob1}) it follows that on the discontinuous surface, the observation error $e_{2}$ satisfies
\begin{equation}\label{e21}
    e_{2}=-\dot{\delta},
\end{equation}
and consequently, it also holds that
\begin{equation}\label{e22}
    \dot{e}_{2}=-\ddot{\delta}.
\end{equation}

Next, by using the equivalent control method \cite{utkin}, it follows from (\ref{ob2})-(\ref{e22}) that, on the discontinuous surface, the filtered version of the discontinuous term $c_{3}\mbox{sign}(e_{1})$ satisfies
\begin{equation}
    \overbar{c_{3}\mbox{sign}(e_{1})}=\ddot{\delta}+a\delta+b\dot{\delta},
\end{equation}
where the over bar denotes filtering.

In order to obtain $\overbar{c_{3}\mbox{sign}(e_{1})}$, we use the low-pass Butterworth filter, written in state-space form
\begin{align}
    \dot{x}_{f}&=A_{f}x_{f}+B_{f}c_{3}\mbox{sign}(e_{1}),\label{fil1}\\
    y_{f}&=C_{f}x_{f},\label{fil2}
\end{align}
where $x_{f}=[x_{f1}\ x_{f2}]^{T}$ and
\begin{align}
     A_{f}&=\left[\begin{array}{cc}0&1\\-\omega_{c}^{2}&-1.4142\omega_{c}\end{array}\right],\ \ B_{f}=\left[\begin{array}{c}0\\\omega_{c}^{2}\end{array}\right], \label{fil3} \\ C_{f}&=\left[\begin{array}{cc}1&0\end{array}\right], \label{fil4}
\end{align}
where $\omega_{c}$ is the cut-off frequency of the filter. Then, by choosing a $\omega_{c}$ that minimizes the phase delay, it is possible to assume that
\begin{equation}\label{filto}
    y_{f}\approx\ddot{\delta}+a\delta+b\dot{\delta}.
\end{equation}

\subsection{Anomaly output detector with robust observer and filter}\label{filin}
In order to design the anomaly output  detector, it is worth nothing that the output $y_{f}$ of the filter, see (\ref{filto}), can indeed be seen as a \textit{residual} because it contains information about the  attack/anomaly in the output of the plant.

Then, similar to (\ref{detector}), we define the following detector
\begin{equation} \label{detector2}
\left\{\begin{aligned}
	\left|y_{f}(t)\right| \leq \alpha_{f} &\longrightarrow \text{no alarm} \quad\longrightarrow\quad alarm=0, \\
	\left|y_{f}(t)\right| > \alpha_{f} &\quad\longrightarrow \text{alarm}\quad\longrightarrow\quad alarm=1.
\end{aligned}\right.
\end{equation}
For the time being, we do not have a formal procedure for determining the upper threshold $\alpha_{f}$. Instead, the value of $\alpha_{f}$ is numerically determined  by directly measuring the output of the filter, for sufficiently long time, and then the value of $\alpha_{f}$ will correspond to the $L_{\infty}$-norm computed from the obtained measurement.
\begin{remark}
Note that for the free-noise case, i.e., $\eta=0$, and considering that the attack is constant, i.e., $\delta:=\delta_{0}$, with $\delta_{0}\in\mathbb{R}$, then the output of the filter, see (\ref{filto}), satisfies
\begin{equation}\label{fre}
   y_{f}\approx a\delta_0.
\end{equation}
Hence, in this case, it is possible  to reconstruct the attack. Note that this fact opens the possibility of using the estimated attack in the control in order to make the system/plant immune to constant attacks by preventing the drift of the estimation error. However, it should be noted that in order to estimate the attack $\delta_{0}$, it is necessary to have knowledge of parameter $a$, which is a parameter of the plant.
\end{remark}

\section{EXAMPLE} \label{sec:example}
This section presents a numerical example, which illustrates the performance of the detectors presented in Sections \ref{sec:residual_detection} and \ref{sec:discontinuous_observer}. In particular, we consider a second order system of the form (\ref{eq:dLTIsystem}) with
\begin{equation}\label{22}
A=\left[\begin{array}{cc}0&1\\-4&-20\end{array}\right], \ B=\left[\begin{array}{c}0\\1\end{array}\right],\ C=\left[\begin{array}{cc}1&0\end{array}\right].
\end{equation}
Furthermore it is assumed that the output is influenced by zero-mean Gaussian noise $\eta(t) \sim \mathcal{N}(0,R)$.

\subsection{Residual-Based Methods}

We now demonstrate these tools and provide a comparison between the filtered and unfiltered version of the chi-squared detector. The system matrices are provided above; here we use a control gain $K=[1,1]$ and an observer gain $L=[0,2]'$. The sampling time is $T=0.001$ and the corresponding measurement noise covariance is $\Sigma_\eta=2$.

The idea behind this comparison is that filtration makes the covariance of the filtered residual smaller and removes the high frequency content of the measurement noise, therefore, the attack will be more identifiable when filter is applied. This ease of detecting the attack is apparent when the alarm rate corresponding to the filtered detector deviates more significantly from the desired/tuned false alarm rate, indicating the detector's sensitivity to the attack is higher. We consider a simple attack to demonstrate this point,
\begin{equation}
    \delta(t)=1.
\end{equation}

\begin{figure}[t]
\centering
\includegraphics[width=0.79\linewidth]{z_CDF_combined.eps}
\caption{This figure shows the performance of the proposed filter-based detector. Using a threshold of $\alpha=3.84$ both detectors have been tuned to provide a $\mathcal{A}^*=0.05$ (5\%) false alarm rate. The green curve plots the cumulative distribution function of the distance measure under attack when computed from the filtered residual $\rho_k$ and exhibits a dramatic difference from what the distribution looks like under normal operation (red). This is quantified by a higher than expected false alarm rate $\mathcal{A}=0.55$ (55\%). In contrast, the blue curve plots the cumulative distribution function of the distance measure under attack when computed from the residual without filtering. Not only does the distribution stay quite similar under attack, but the alarm rate is $7\%$ which is only $2\%$ bigger than the false alarm rate.}\label{fig:cdf}
\end{figure}

\begin{figure}[hb!]
\centering
\includegraphics[width=0.73\linewidth]{r_pdf_combined.eps}
\caption{The probability density functions, since they are not normalized, aptly demonstrate the role of shrinking covariance in the detection problem. The difference between the attacked and nominal behavior is more dramatic in the filtered residual $\rho_k$ (lower panel) than in the residual $r_k$ (upper panel).}\label{fig:pdf}
\end{figure}

Figure \ref{fig:cdf} shows the cumulative distribution of distance measure calculated from $r_k$ under attack (blue) compared with the cumulative distribution of the distance measure under no attacks (red). Because the distributions are closely similar, the conventional chi-squared detector is not able to easily distinguish the attack. This figure also compares the cumulative distribution of the distance measure calculated from the filtered residual $\rho_k$ (green). In this case the small attack significantly changes the distribution of the distance measure and detection is easy. This is captured quantitatively by comparing the false alarm rates. Here the detectors were tuned to provide $\mathcal{A}^*=0.05$ (5\%) false alarms (i.e., under no attacks) by selecting the threshold as $\alpha=3.84$.

Under attack the conventional detector's alarm rate was 7\% whereas the filtered detector's alarm rate was 55\%, demonstrating the dramatic difference in sensitivity provided by the additional filtration.

Because the residual and filtered residual are not themselves normalized (as opposed to the distance measure which is normalized), the role that covariance plays is quite clear in Fig. \ref{fig:pdf}. Here the reduced covariance of the filtered residual clearly distinguishes the attack scenario from the normal operation. This distinction is less obvious in the unfiltered residual.

\subsection{Discontinuous Observer}
For the case of the observer-based detector discussed in Section \ref{sec:discontinuous_observer}, we consider observer (\ref{obsi1})-(\ref{obsi2}) with $A$, $B$, and $C$ as given in (\ref{22}) and the gains $c_{1}=c_{2}=5$, and $c_{3}=12$. For this choice of $c_{i}$, $i=1,\ldots,3$, the observation error is globally asymptotically stable for the noise free ($\eta=0$) and attack free ($\delta=0$) case, see \cite{robustobs}.  Furthermore, the chosen cut-off frequency $\omega_{c}$ for the low-pass filter (\ref{fil1})-(\ref{fil4}) is $\omega_{c}=12$ [rad/s].

\begin{figure}[t]
	\centering
		\includegraphics[scale=0.52]{attacks.eps}
	\caption{Residual signal for the noise free case. a) Residual for the time varying attack (\ref{varying}). Black line: measured residual. Gray line: predicted residual from (\ref{filto}). b) Residual for the constant attack (\ref{cons}). Black line: measured residual. Gray line: predicted residual from Eq. (\ref{fre}).}
	\label{fig:attacks}
\end{figure}
Next, the performance of the observed-based detector is evaluated for two different attacks, namely a time varying attack and a constant attack. The time varying attack is described by
\begin{equation}\label{varying}
\delta(t)=\left\{\begin{array}{ccc}0&\mbox{if}&t\leq15,\\
                                    0.1\sin{t}&\mbox{otherwise},&\end{array}\right.
\end{equation}
whereas the constant attack is given by
\begin{equation}\label{cons}
\delta(t)=\left\{\begin{array}{ccc}0&\mbox{if}&t\leq15,\\
                                   0.1&\mbox{otherwise}.&\end{array}\right.
\end{equation}
Figure \ref{fig:attacks} shows the numerical and theoretical residuals $y_{f}$. Clearly, in both cases, the measured residual corresponds to the predicted residuals (\ref{filto}) and (\ref{fre}). The numerical results have been obtained  by using Simulink (MATLAB) and the solver Runge-Kutta with fixed step size of 0.001.

On the other hand, for the case of noisy output, it is assumed that the output contains zero-mean Gaussian noise $\eta$ with covariance $\Sigma_\eta=0.001$. In this case, by using the empirical procedure mentioned in Subsection \ref{filin}, it has been obtained that the threshold value of the detector (\ref{detector2}) is $\alpha_{f}=1.55$. It can be shown that if the attacks (\ref{varying}) and (\ref{cons}) are added to the noisy output, then the detector raises an alarm, i.e., those attacks can be identified.

\section{CONCLUSIONS}

In this paper we've introduced two techniques to exploit the a priori knowledge that the noise is a high frequency disturbance to make the attack/anomaly detection task easier. In the first case, we retrofit the conventional chi-squared detector with a filter to reduce the covariance of the detector's distance measure statistic. In the second case, we adopted a robust observer with a discontinuous term that reveals the presence of an attack or anomaly in the system. As fault and anomaly detection tools are poised to bridge to industry, such practical techniques are compelling ways to boost performance and constrain potential attackers even further.

\addtolength{\textheight}{-12cm}

\bibliographystyle{IEEEtran}

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction}
Recall that a {\it M\"obius} manifold is an oriented $(G,X)$-manifold (see e.g., \cite[3.1.2]{thurston}) of dimension $n$ with $X=S^n$ and $G$ the group $\mbox{M{\"o}b}_n$ of orientation-preserving M\"obius transformations (see e.g., \cite{hert} for details). We suppose that the charts participating in the definition of the M\"obius structure are orientation preserving. Two M\"obius structures on the same manifold are said to be {\it equivalent}
if one of them is induced from the other under some diffeomorphism homotopic to the identity.

For a M\"obius manifold $M$, we set $-M$ to be the M\"obius manifold ``inverse to $M$'', i.e., as a smooth oriented manifold, $-M$ is equal to $M$ with the orientation reversed, and the charts that define the M\"obius structure on $-M$ are the compositions
of the charts for $M$ and an orientation-reversing M\"obius transformation. In the sequel, all M\"obius diffeomorphisms are assumed to be orientation preserving. A diffeomorphism $M_1\to M_2$ between two M\"obius manifolds is said to be {\it anti-M\"obius}, if the composition
$-M_1\stackrel{\mathrm{Id}}\longrightarrow M_1\to M_2$ is M\"obius.

The definition of M\"obius structures can be naturally extended to manifolds with boundary.
Let $M$ is a M{\"o}bius manifold, and let $p:N_1\to N_2$ is an anti-M{\"o}bius diffeomorphism between two disjoint components of $\partial M$. Then the
manifold $M'$ obtained from $M$ by identifying $x$ and $p(x)$ for all
$x\in N_1$ can be equipped with a canonical M{\"o}bius structure.

Hence, one can define the {\it cobordism group $\Omega_n^{\mbox{\scriptsize \rm M{\"o}b}}$ of $n$-dimensional M\"obius manifolds} as follows.
We set the elements of $\Omega_n^{\mbox{\scriptsize \rm M{\"o}b}}$ to be $n$-dimensional closed M\"obius manifolds considered modulo
the equivalence relation $\sim$ with $M_1\sim M_2$, iff there exists an $n+1$-dimensional compact M\"obius manifold $N$ such that
$\partial N$ is M\"obius diffeomorphic to $M_1\sqcup -M_2$. The group operation
is induced by taking the disjoint union. The main result of this note is the following theorem.

\begin{theorem}
\label{nqs}$\Omega_1^{\mbox{\scriptsize \rm M{\"o}b}}=0$.
\end{theorem}

The proof is straightforward: we describe in section \ref{sec2} the M\"obius structures on the circle, and for each of those we construct in section \ref{sec3} a bounding surface. In section \ref{sec2}, we also show (lemma \ref{projective} and remark 2) that some cases are missing from N. Kuiper's classification \cite{kuiper} of projective structures
on the circle, and we complete that classification.

{\bf Remark 1.} One can define for any $n\geq 1$ the cobordism category $\EuScript{MC}_n$ of $n$-dimensional M\"obius manifolds as follows. We set the objects of $\EuScript{MC}_n$ to be compact M\"obius manifolds of dimension $n$,
and for any two objects $N_1$ and $N_2$ of $\EuScript{MC}_n$, morphisms from $N_1$ to $N_2$
are defined to be $n+1$-dimensional compact M\"obius manifolds $N$ such that $\partial N=N_1\sqcup -N_2$ (modulo M\"obius diffeomorphisms that are identical on $\partial N$). The categories $\EuScript{MC}_n$ are strict symmetric monoidal, and
the above theorem \ref{nqs} can be reformulated as $K_0(\EuScript{MC}_1)=0$. It would be interesting to learn more about $\EuScript{MC}_1$, e.g., a presentation by generators and relations etc.

The author is grateful to P. Vogel for proposing the problem and useful discussions.
\section{M{\"o}bius structures on the oriented circle}
\label{sec2}
In this section we recall the classification of M\"obius structures on the circle $S^1$.
A M{\"o}bius structure on a 1-dimensional manifold can
be viewed as a $(G,X)$-structure, where
$X$ is any circle $C$ in ${\mathbb C} P^1$,
and $G$ is the group of automorphisms of ${\mathbb C} P^1$ that
preserve both disks bounded by $C$. We shall often use the projective model, i.e., we shall take the circle ${\mathbb R} P^1\subset{\mathbb C} P^1$, so that $G$ becomes the group of orientation-preserving projective automorphisms of ${\mathbb R} P^1$.
Given a M\"obius structure $S$ on $S^1$, we denote by $[S]$ the equivalence class of $S$, and we define $\EuScript{M}(S^1)$ to be the set of all equivalence classes of M\"obius structures on $S^1$.

Let us identify the universal cover of $S^1$ with the (canonically oriented) ${\mathbb R}$.
Any M\"obius structure $S$ on $S^1$ constructed below will be defined by
an orientation-preserving developing map
$ F_S$ from ${\mathbb R}$ to a model circle $C\subset{\mathbb C} P^1$ and a holonomy homomorphism $f_S$ from ${\mathbb Z}$ to the subgroup of $\mathrm{PSL}_2({\mathbb C})$ corresponding to $C$;
the maps $F_S$ and $f_S$ should satisfy the condition
\begin{equation}
\label{eqqq}
 F_S(x+1)= f_S(1)\cdot F_S(x).
\end{equation}

Let $A_S$ be a matrix
that represents $ f_S(1)$ (this matrix is defined up to a
sign). We
shall distinguish the following cases:

\begin{enumerate}
\item $|tr(A_S)|\leq 2$, and $A_S$ is diagonalisable.
\item $|tr(A_S)|=2$, $A_S$ is nondiagonalisable.
\item $|tr(A_S)|>2$.
\end{enumerate}

The structures that verify (1), (2), (3), will be called respectively
{\it elliptic, parabolic} and {\it hyperbolic}.

For any $\alpha>0$ and $n>0$ let $\psi_{n,\alpha}$ be a continuous map $[0,1]\to{\mathbb R} P^1$ such that 1. $\psi_{n,\alpha}=\mathrm{e}^{\alpha x}-1,$ when $x$ is close to 0,
$\psi_{n,\alpha}=\mathrm{e}^{\alpha}(\mathrm{e}^{\alpha x}-1),$ when $x$ is close to 1, 2. $\psi^{-1}_{n,\alpha}(0)$ consists of $n+1$ elements, and 3.
$\psi_{n,\alpha}|(0,1)$ is an orientation-preserving local diffeomorphism.
Analogously, for any couple $(n,\varepsilon)$ (where $n>0$ is an integer, and $\varepsilon=\pm 1$)
let $\xi_{n,\varepsilon}$ be a continuous map $[0,1]\to{\mathbb R} P^1$ such that 1. $\xi_{n,\varepsilon}=-\frac1x,$ when $x$ is close to 0,
$\xi_{n,\varepsilon}=\varepsilon+\frac{1}{1-x},$ when $x$ is close to 1, 2. $\xi^{-1}_{n,\varepsilon}(\infty)$ consists of $n+1$ elements, and 3.
$\xi_{n,\varepsilon}|(0,1)$ is an orientation-preserving local diffeomorphism.
The following theorem gives an explicit representative for each element of $\EuScript M (S^1)$ (cf. \cite{kuiper}).

\begin{theorem}
\label{circle}
\begin{enumerate}
\item
If two M{\"o}bius structures on $S^1$ belong to different types, they are
not equivalent.
\item Equivalence classes of hyperbolic M{\"o}bius structures are
parametrised by couples $(n,\alpha)$, where $\alpha$ is a positive real number, and $n$ is a nonnegative integer.

The corresponding maps $ F_S:{\mathbb R}\to {\mathbb R} P^1,  f_S:{\mathbb Z}\to \mathrm{PSL}_2({\mathbb R})$ can be chosen as follows. Let $f_S(1)$ be the map $z\mapsto \mathrm{e}^{\alpha}z$.
If $n=0$, set $ F_S(x)=\mathrm{e}^{\alpha x}$, otherwise set $ F_S(x)=\psi_{n,\alpha} (x)$ for $x\in [0,1]$ and extend $ F_S$ to ${\mathbb R}$
using (\ref{eqqq}). Denote the equivalence class of hyperbolic M\"obius structures
corresponding to the couple $(n,\alpha)$ by
$H_{n,\alpha}$.
\item
Equivalence classes of parabolic M{\"o}bius
structures are parametrised by elements of the set
$(\{\mbox{nonnegative integers}\}\times\{\pm 1\})\setminus\{(0,-1)\}$.

The equivalence class that corresponds to the couple $(0,1)$ is represented by the M\"obius structure defined by the following maps
$ F_S:{\mathbb R}\to {\mathbb R} P^1,  f_S:{\mathbb Z}\to \mathrm{PSL}_2({\mathbb R})$: set $ F_S(x)=x,  f_S(1)(z)=z+1$.

The equivalence class corresponding to the couple $(n,\varepsilon)$ (where $n$ is a positive integer, and $\varepsilon$ is 1 or $-1$)
is represented by the M\"obius structure defined by the maps $ F_S:{\mathbb R}\to {\mathbb R} P^1,  F_S:{\mathbb Z}\to \mathrm{PSL}_2({\mathbb R})$ that
can be constructed as follows: set $ f_S(1)(z)=z+\varepsilon$, set $ F_S(x)=\xi_{n,\varepsilon} (x)$ for $x\in [0,1]$ and extend $ F_S$ to ${\mathbb R}$
using (\ref{eqqq}). Denote the equivalence class of parabolic M\"obius structures
corresponding to the couple $(n,\varepsilon)$ by
$P_{n,\varepsilon}$.
\item Equivalence classes of elliptic M{\"o}bius structures are
parametrised by positive real numbers. The equivalence class corresponding to $\alpha>0$ can be represented by the M\"obius structure defined by the maps $F_S:{\mathbb R}\to\{\mbox{complex numbers of modulus 1}\}$ and $f_S(1):{\mathbb Z}\to\mathrm{PU}_{1,1}$ given by $F_S(x)=\mathrm{e}^{i\alpha x},f_S(1)(z)=\mathrm{e}^{i\alpha}z$.
Denote this equivalence class by $E_\alpha$.
\end{enumerate}
\end{theorem}

{\bf Proof of theorem \ref{circle}.}
The argument is analogous to the one used in \cite{kuiper} to classify projective circles modulo
projective diffeomorphisms (but see remark 2 below).

First, let us show that any M\"obius structure $S$ on $S^1$ is equivalent to some structure described in theorem \ref{circle}.
Suppose, e.g., that $S$ is hyperbolic. Choose a base point in $S^1$, and denote by
$ F_S$ the corresponding developing map ${\mathbb R}\to{\mathbb R} P^1$; let $ f_S:{\mathbb Z}\to \mathrm{PSL}_2({\mathbb R})$ be the homomorphism
such that $ F_S$ and $ f_S$ satisfy (\ref{eqqq}). Then, replacing $S$ by an equivalent structure, we can assume that
$ f_S$ is given by the formula
$ f_S(n)(z)=\mathrm{e}^{\alpha n}z$ for some $\alpha>0$. The image of $ F_S$ is either the whole ${\mathbb R} P^1$ or not; we can assume (changing the base
point, if necessary) that $ F_S(0)$ is 1 in the first case and 0 in the second.
Now it is easy to see that $S$ is equivalent to some structure from the second
assertion of theorem \ref{circle}. The parabolic and elliptic cases are considered in an analogous way.

Now, let us prove that the structures introduced in theorem \ref{circle} are pairwise nonequivalent. Let $S$ be a M\"obius structure on $S^1$.
Suppose that $S$ is
defined by an orientation-preserving local diffeomorphism
$ F_S:{\mathbb R}\to {\mathbb R} P^1$ and a homomorphism $ f_S:{\mathbb Z}\to \mathrm{PSL}_2({\mathbb R})$ that satisfy (\ref{eqqq}).
We can associate the following invariants to $S$: the conjugacy class of
$f_S(1)$ in $\mathrm{PSL}_2({\mathbb R})$ and the number
$$\max_{x\in{\mathbb R},y\in{\mathbb R} P^1}\#( F_S^{-1}(y)\cap [x,x+1]).$$
It can be easily checked that these invariants are sufficient to distinguish any two different structures defined in theorem \ref{circle}.$\diamondsuit$

A {\it projective} structure on the circle is an $(\mathrm{PGL}_2({\mathbb R}),{\mathbb R} P^1)$-structure. Every M\"obius circle is a projective circle; conversely, one can orient a projective circle so that it will become a M\"obius circle.
\begin{lemma}
\label{projective}
Any two different M\"obius structures on $S^1$ defined in theorem \ref{circle} are not equivalent as projective structures.
\end{lemma}
{\bf Proof of Lemma \ref{projective}.}
We proceed as in the proof of theorem \ref{circle}. Suppose that the projective structure $S$ on $S^1$ is defined by
a couple $( F_S, f_S)$, where $ F_S$ is a local diffeomorphism ${\mathbb R}\to{\mathbb R} P^1$ and $f_S:{\mathbb Z}\mapsto\mathrm{PSL}_2({\mathbb R})$ is a homomorphism that satisfy (\ref{eqqq}).
It can be easily checked that the conjugacy class of $ F_S(1)$ in $\mathrm{PGL}_2({\mathbb R})$ and the number
$$\max_{x\in{\mathbb R},y\in{\mathbb R} P^1}\#( F_S^{-1}(y)\cap [x,x+2])$$
depend only on the equivalence class of $S$.
These invariants of projective structures distinguish any two nonequivalent M\"obius structures on $S^1$.$\diamondsuit$

{\bf Remark 2.} Due to Lemma \ref{projective} there exists a bijection $\EuScript M(S^1)
\leftrightarrow
($projective circles modulo projective diffeomorphisms). Hence, a classification of projective circles can be obtained from theorem \ref{circle}; the classification
given in \cite{kuiper} in not quite correct: if $n>0$, then
$S^1$ provided
with a M\"obius structure of class $P_{n,1}$ and $S^1$ provided
with a M\"obius structure of class $P_{n,-1}$ are not isomorphic as projective circles.

\label{anti}

{\bf Remark 3.} The interpretation in terms of developing maps allows us to introduce a topology on $\EuScript M(S^1)$. However, the resulting topological space
is nasty. The space $\EuScript M'(S^1)=\EuScript M(S^1)\setminus\{E_{2\pi k}|\mbox{$k$ is an integer $>0$}\}$
is the non-Hausdorff topological 1-manifold obtained by identifying the thin lines on Figure \ref{struct}; $\EuScript M'(S^1)$ is open in $\EuScript M(S^1)$,
and for any $k>0$ we can take the system of sets of the form
$\{E_{2\pi k}\}\cup (\mbox{a neighbourhood of $P_{k,1}$ in $\EuScript M'(S^1)$})\cup(\mbox{a neighbourhood of $P_{k,-1}$ in $\EuScript M'(S^1)$})$
as a local neighbourhood basis at $E_{2\pi k}$.
\begin{figure}\centering
\epsfbox{struct2.eps}
\caption{$\EuScript M'(S^1)
$}\label{struct}
\end{figure}

\section{Proof of theorem \ref{nqs}.}\label{sec3}

In order to prove theorem \ref{nqs}, it is enough
to show that any element of $\EuScript M(S^1)$ is equal to $[S]$, where $S$ is
the M\"obius structure on the boundary of some 2-dimensional M{\"o}bius surface.
Notice that a M\"obius structure on a 2-dimensional manifold
is the same as a $(\mathrm{PSL}_2({\mathbb C}),{\mathbb C} P^1)$-structure.

Let $N$ be the sphere minus three disjoint open disks (``the pants'').
The manifold $\tilde N$ is represented on Figure \ref{uncov1}.

\begin{figure}\centering
\epsfbox{uncov1.eps}
\caption{}\label{uncov1}
\end{figure}

Denote by
$\mathbf{F}(A,B)$ the free group on the generators $A,B$.
Suppose that $\mathbf{F}(A,B)$ acts on $\tilde N$ as shown on
Figure \ref{uncov} (the octagon in the middle
is a fundamental domain; the right and the left octagons are the
images of the fundamental domain under $A$ and $B$ respectively).

\begin{figure}\centering
\epsfbox{uncov.eps}
\caption{}\label{uncov}
\end{figure}

We shall equip $N$ (and hence, $\partial N$) with a M\"obius structure by constructing maps
$F:\tilde N\to {\mathbb C} P^1, f:\mathbf{F}(A,B)\to \mathrm{PSL}_2({\mathbb C})$ that
satisfy $$F(g\cdot x)=f(g)\cdot x,g\in\mathbf{F}(A,B),x\in\tilde N.$$ This can be done by choosing two
elements $a,b\in \mathrm{PSL}_2({\mathbb C})$ and an embedded rectangular octagon
in ${\mathbb C}$ with vertices $X_1,X_2,\ldots, X_8$ (see Figure \ref{fund};
arrows indicate the
orientations of the boundary)
such that the following conditions are satisfied:
\begin{description}\label{cond1}
\item[\rm (C1)] All segments $X_1X_2,\ldots,X_8X_1$ are arcs of circles;
denote the corresponding circles by $C_{X_1X_2},\ldots,C_{X_8X_1}$.
\item[\rm (C2)] $a(X_1)=X_4, a(X_2)=X_3, b(X_5)=X_8, b(X_6)=X_7$.
\item[\rm (C3)] $a(C_{X_2X_3})=C_{X_2X_3}, a(C_{X_8X_1})=C_{X_4X_5},
b(C_{X_6X_7})=C_{X_6X_7}, b(C_{X_4X_5})=C_{X_8X_1}$.
\end{description}
\begin{figure}\centering
\epsfbox{fund1.eps}
\caption{}\label{fund}
\end{figure}

\begin{figure}\centering
\epsfbox{oct1.eps}
\caption{}\label{oct}
\end{figure}

Let $x$ be a positive real number. Define the elements
$a_x, b\in \mathrm{PSL}_2({\mathbb C})$ respectively by $z\mapsto z+x, z\mapsto\frac{z}{1-z}$. The transformation $b\circ a_x$ is
represented by a matrix with trace $2-x$. Choose
the octagon as shown in Figure \ref{oct}. The horizontal line
at the bottom is the real axis, the circles $C_1$ and $C_2$ on
Figure \ref{oct} are defined respectively
by the equations $(\mathop{\rm Re}\nolimits z)^2-2\mathop{\rm Re}\nolimits z+(\mathop{\rm Im}\nolimits z)^2=0$
and $(\mathop{\rm Re}\nolimits z)^2+2\mathop{\rm Re}\nolimits z+(\mathop{\rm Im}\nolimits z)^2=0$.
The circles $C_1'$ and $C_2'$ are chosen in such a way that 1. $b(C_1')=C_2',a_x(C_2')=C_1'$,
2. $C_1'$ is orthogonal to $C_1$ and depends continuously on $x$, 3.
$C_1'\subset\{z|\mathop{\rm Re}\nolimits z\geq 0\}$, and the imaginary part of
the center of $C'_1$ is $>0$. Note that there are many ways to choose the circles $C_1'$ and $C_2''$ that satisfy these conditions.
The arc $X_6X_7$ is an arc of a sufficiently small circle in the upper half-plane that is tangent to the real axis at $0$.

Note that we have $b(z)=-\bar z$ for $z=1+\mathrm{e}^{it}$, i.e., the
restriction of $b$ to $C_1$ is the symmetry
with respect to the imaginary axis. This implies easily that
for any $x>0$, the octagon on Figure \ref{oct} and the
maps $a=a_x$ and $b$ satisfy the conditions (C1)-(C3).

For any $x>0$ we obtain M{\"o}bius structures on $N$ and all
components of $\partial N$. The components of $\partial N$ that correspond to $X_2X_3$ and $X_6X_7$
will always have parabolic structures that belong to $P_{0,1}$,
and the equivalence class of the structure on the third component changes as we change $x$.
Denote this structure by $S_x$; $S_x$ is elliptic for $x<4$,
parabolic for $x=4$ and hyperbolic for $x>4$.
Denote by $N_1(x)$ the M\"obius surface obtained by gluing together two parabolic components of
$\partial N$ (note that due to  Lemma \ref{projective}, for any M\"obius circle $C$ there exists an anti-M\"obius diffeomorphism $C\to C$).

Let us determine the equivalence class of $S_x$. Let $ F_{S_x}$ be a developing map ${\mathbb R}\to{\mathbb R} P^1$ that corresponds to $S_x$.

Note that $[S_4]=P_{0,1}$, which implies that the image of $ F_{S_4}$ is not the whole ${\mathbb R} P^1$.
Hence, if $D>0$, then we have $ F_{S_x}([0,D])\neq{\mathbb R} P^1$ for any $x$ sufficiently close to 4. This
implies the following lemma:

\begin{lemma}
\label{xx}
For any sufficiently
small $\alpha>0$ there exist $x_1, x_2$ such that the M\"obius structure on $\partial N_1(x_1)$ (respectively, $\partial N_1(x_2)$) belongs to
$H_{0,\alpha}$ (respectively, to $E_{\alpha}$).
\end{lemma}

Let us define the action of ${\mathbb N}$ on $\EuScript M(S^1)$ as follows: for any $k\in{\mathbb N}, \alpha>0,n\geq 0$ and
$\varepsilon\in\{\pm 1\}$ set
$kH_{n,\alpha}=H_{kn,k\alpha}, kP_{n,\varepsilon}=P_{kn,\varepsilon}, kE_{\alpha}=E_{k\alpha}$.

\begin{lemma}
\label{x}
Let $C'$ and $C''$ be oriented circles. Suppose
that $C''$ is equipped with a M\"obius structure, and denote this structure by $S''$. Let $p:C'\to C''$ be an orientation-preserving $k$-sheeted
covering. Denote by $S'$ the inverse image of $S''$ with respect to $p$. We have $[S']=k[S'']$.
\end{lemma}
$\diamondsuit$
\begin{lemma}
\label{xxxx}
Let $N$ be a compact M\"obius surface with one boundary component and nontrivial fundamental group, and let $k>0$ be an integer.
Denote by $S$ the M\"obius structure on $\partial N$.
There exists
a compact M\"obius surface $N'$ such that $\partial N'$ is a circle, and the M\"obius structure on $\partial N'$ belongs to $k[S]$.
\end{lemma}
{\bf Proof of Lemma \ref{xxxx}.} Let $\mathfrak S_k$ be the symmetric group on $k$ elements. If $k$ is odd, then there exists a homomorphism
$\pi_1(N)\to \mathfrak S_k$ that takes a generator of $\pi_1(\partial N)$ to some cycle of length $k$. Indeed, we can choose
a system $v_1,w_1,\ldots, v_g,w_g$ of free generators of $\pi_1 (N)$ so that $\pi_1 (\partial N)$ will be spanned by
$[v_1,w_1]\cdots [v_g,w_g]$. A required homomorphism can be constructed using the fact that in $\mathfrak{S}_k$ a cycle of length $k$ is a commutator
(e.g., we have
$$(1\ldots k)=\left(1\frac{k+3}{2}\ldots k\right)\left(12\ldots\frac{k+1}{2}\right)=\sigma\left(12\ldots\frac{k+1}{2}\right)^{-1}\sigma^{-1}\left(12\ldots\frac{k+1}{2}\right)$$ for some $\sigma\in\mathfrak{S}_k$).

Hence, for $k$ odd, there exists a surface
$N'$ with one boundary component and a $k$-sheeted covering $N'\to N$.
This, together with Lemmas  \ref{xx} and \ref{x}, allows us to construct for any $\alpha >0$ a M\"obius surface, whose boundary carries
a M\"obius structure that belongs to $E_{\alpha}$.

Now cut out a small M\"obius disk from $N$. The resulting M\"obius surface has two boundary components; one of these
components has the structure $S$, and the structure on the other one belongs to $E_{2\pi}$.
Note that if $M$ is a compact surface with two boundary components, then a generator
of one of these components can be taken as a free generator of $\pi_1(M)$; hence, for any $k>0$ there exists a $k$-sheeted covering $M'\to M$
that is cyclic over one of the components of $\partial M$. We have already constructed for any integer $l>0$ a M\"obius surface, whose boundary
carries a structure from $E_{2\pi l}$, and the lemma follows.
$\diamondsuit$

Lemmas \ref{xx}, \ref{x}, \ref{xxxx} imply that any of the classes
$P_{0,1}, H_{0,\alpha}$ or $E_{\alpha}, \alpha>0$ is the equivalence class of the M\"obius structure of the boundary of
some M\"obius surface. In order to complete the proof of theorem \ref{nqs}, we can proceed as follows. Suppose that $x>2$ is a real number and consider
the rectangular octagon represented on Figure \ref{oct2}.

\begin{figure}\centering
\epsfbox{oct2.eps}
\caption{}\label{oct2}
\end{figure}

Here the circles $C_1$ and $C_2$ are the same as in the proof of Lemma \ref{xx}. Suppose
that the centers of $C'_1$ and $C'_2$ are the points $\pm \frac x2 +i$, and the radii of these circles are $\frac x2-1$. As above, the octagon
and the maps $a_x$ and $b$ give us for any $x>2$ a M\"obius structure on $N$ (the pants). The induced M\"obius structure
on the component of $\partial N$ that corresponds to
the arc $X_2 X_3$ (respectively, the arc $X_6 X_7$) belongs to the class $P_{0,1}$ (respectively, $P_{1,-1}$). Denote by $S'_x$ the M\"obius structure on the
third component of $\partial N$. The structure $S'_x$ is elliptic for $x<4$; by attaching M\"obius surfaces to the elliptic component of $\partial N$ and to
the $P_{0,1}$-component, we obtain a M\"obius surface, whose boundary has a M\"obius structure of class $P_{1,-1}$, which means that
we can eliminate the $P_{1,-1}$-component as well, i.e., for any $x>2$ there exists a M\"obius surface, whose boundary carries the M\"obius structure
$S'_x$.

Let us determine the equivalence class of $S'_x,x\geq 4$. Note that we can choose a developing map $ F_{S'_x}:{\mathbb R}\to C'_1$ for the structure $S'_x$
so that $ F_{S'_x}$ sends the segment $[0,1/2]$ to the arc $X_5X_4$ and the segment $[1/2, 1]$ to the arc $a_x(X_1 X_8)$. It is easy to check that for $x\geq 4$,
the fixed points of $a_x\circ b$ are exactly the intersection points of $C'_1$ and the real axis.

Hence, for any $x>4$ the set $ F_{S'_x}^{-1}(\mbox{the fixed points of $a_x\circ b$})\cap [0,1)$ consists of two elements, which implies that for any $x>4$ we have
$S'_x\in H_{1,\alpha}$ with $\alpha=2\mathop{\rm arccosh}\nolimits (x/2-1)$.

An analogous argument shows that $S'_4\in P_{1,1}$ (note that $a_4\circ b$ acts on $C'_1$ ``clockwise''). The proof of theorem \ref{nqs} is easily completed
using Lemmas \ref{x} and \ref{xxxx}. $\diamondsuit$