\section{#1}\setcounter{equation}{0}}
\renewcommand{\theequation}{\thesection.\arabic{equation}}
\newcommand{\begin{equation}}{\begin{equation}}
\newcommand{\end{equation}}{\end{equation}}
\newcommand{\begin{eqnarray}}{\begin{eqnarray}}
\newcommand{\end{eqnarray}}{\end{eqnarray}}
\newcommand{\nonumber}{\nonumber}
\newcommand{{\underline{1}}}{{\underline{1}}}
\newcommand{{\underline{2}}}{{\underline{2}}}
\newcommand{{\underline{3}}}{{\underline{3}}}
\newcommand{{\underline{4}}}{{\underline{4}}}
\newcommand{{\underline{a}}}{{\underline{a}}}
\newcommand{{\underline{b}}}{{\underline{b}}}
\newcommand{{\underline{c}}}{{\underline{c}}}
\newcommand{{\underline{d}}}{{\underline{d}}}
\newcommand{{\underline{i}}}{{\underline{i}}}
\newcommand{{\underline{j}}}{{\underline{j}}}
\newcommand{\ku}{{\underline{k}}}
\newcommand{{\underline{l}}}{{\underline{l}}}
\newcommand{{\underline{I}}}{{\underline{I}}}
\newcommand{{\underline{J}}}{{\underline{J}}}

\newcommand{{\mathbb R}}{{\mathbb R}}
\newcommand{{\mathbb C}}{{\mathbb C}}
\newcommand{{\mathbb Q}}{{\mathbb Q}}
\newcommand{{\mathbb Z}}{{\mathbb Z}}
\newcommand{{\mathbb N}}{{\mathbb N}}

\def\dt#1{{\buildrel {\hbox{\LARGE .}} \over {#1}}}

\newcommand{\bm}[1]{\mbox{\boldmath$#1$}}

\def\double #1{#1{\hbox{\kern-2pt $#1$}}}

\newcommand{{\hat{m}}}{{\hat{m}}}
\newcommand{{\hat{n}}}{{\hat{n}}}
\newcommand{{\hat{p}}}{{\hat{p}}}
\newcommand{{\hat{q}}}{{\hat{q}}}
\newcommand{{\hat{r}}}{{\hat{r}}}
\newcommand{{\hat{a}}}{{\hat{a}}}
\newcommand{{\hat{b}}}{{\hat{b}}}
\newcommand{{\hat{c}}}{{\hat{c}}}
\newcommand{{\hat{d}}}{{\hat{d}}}
\newcommand{{\hat{e}}}{{\hat{e}}}
\newcommand{{\hat{M}}}{{\hat{M}}}
\newcommand{{\hat{N}}}{{\hat{N}}}
\newcommand{{\hat{A}}}{{\hat{A}}}
\newcommand{{\hat{B}}}{{\hat{B}}}
\newcommand{{\hat{C}}}{{\hat{C}}}
\newcommand{{\hat{i}}}{{\hat{i}}}
\newcommand{{\hat{j}}}{{\hat{j}}}
\newcommand{{\hat{k}}}{{\hat{k}}}
\newcommand{{\hat{l}}}{{\hat{l}}}

\newcommand{{\hat{\alpha}}}{{\hat{\alpha}}}
\newcommand{{\hat{\beta}}}{{\hat{\beta}}}
\newcommand{{\hat{\gamma}}}{{\hat{\gamma}}}
\newcommand{{\hat{\delta}}}{{\hat{\delta}}}
\newcommand{{\hat{\rho}}}{{\hat{\rho}}}
\newcommand{{\hat{\tau}}}{{\hat{\tau}}}

\newcommand{{\dot\gamma}}{{\dot\gamma}}
\newcommand{{\dot\delta}}{{\dot\delta}}

\newcommand{{\tilde{\sigma}}}{{\tilde{\sigma}}}
\newcommand{{\tilde{\omega}}}{{\tilde{\omega}}}

\renewcommand{\Bar}{\overline}

\newcommand{{\underline{\alpha}}}{{\underline{\alpha}}}
\newcommand{{\underline{\beta}}}{{\underline{\beta}}}
\newcommand{{\underline{\gamma}}}{{\underline{\gamma}}}
\newcommand{{\underline{\delta}}}{{\underline{\delta}}}
\newcommand{{\underline{\rho}}}{{\underline{\rho}}}
\newcommand{{\underline{\tau}}}{{\underline{\tau}}}

\newcommand{{\underline{\ad}}}{{\underline{\ad}}}
\newcommand{{\underline{\bd}}}{{\underline{\bd}}}
\newcommand{{\underline{\gd}}}{{\underline{{\dot\gamma}}}}
\newcommand{{\underline{\dd}}}{{\underline{{\dot\delta}}}}
\newcommand{{\underline{\dot{\rho}}}}{{\underline{\dot{\rho}}}}
\newcommand{{\underline{\dot{\tau}}}}{{\underline{\dot{\tau}}}}

\newcommand{{\underline{\hal}}}{{\underline{{\hat{\alpha}}}}}
\newcommand{{\underline{\hbe}}}{{\underline{{\hat{\beta}}}}}
\newcommand{{\underline{\hga}}}{{\underline{{\hat{\gamma}}}}}
\newcommand{{\underline{\hde}}}{{\underline{{\hat{\delta}}}}}
\newcommand{{\underline{\hrh}}}{{\underline{{\hat{\rho}}}}}

\newcommand{{\nabla}}{{\nabla}}
\newcommand{{\bar{\nabla}}}{{\bar{\nabla}}}

\newcommand{{\bar{\sigma}}}{{\bar{\sigma}}}

\newcommand{{\theta}}{{\theta}}
\newcommand{{\bar{\theta}}}{{\bar{\theta}}}
\newcommand{{\bar{\theta}}}{{\bar{\theta}}}

\newcommand{{\dot{\rho}}}{{\dot{\rho}}}
\newcommand{{{\tau}}}{{{\tau}}}
\newcommand{{\dot{\ta}}}{{\dot{{{\tau}}}}}

\newcommand{{(u^+u^-)}}{{(u^+u^-)}}

\newcommand{{{\bar{\zeta}}}}{{{\bar{\zeta}}}}

\newcommand{{\bm L}}{{\bm L}}
\newcommand{{\bm R}}{{\bm R}}
\newcommand{{\boxplus}}{{\boxplus}}
\newcommand{{\boxminus}}{{\boxminus}}
\newcommand{{\oplus}}{{\oplus}}
\newcommand{{\ominus}}{{\ominus}}

\newcommand{{\overline{I}}}{{\overline{I}}}
\newcommand{{\overline{J}}}{{\overline{J}}}
\newcommand{{\overline{K}}}{{\overline{K}}}
\newcommand{{\overline{L}}}{{\overline{L}}}
\newcommand{{\overline{M}}}{{\overline{M}}}
\newcommand{{\overline{N}}}{{\overline{N}}}
\newcommand{{\overline{P}}}{{\overline{P}}}
\newcommand{{\overline{Q}}}{{\overline{Q}}}

\newcommand{{\underline{I}}}{{\underline{I}}}
\newcommand{{\underline{J}}}{{\underline{J}}}
\newcommand{{\underline{K}}}{{\underline{K}}}
\newcommand{{\underline{L}}}{{\underline{L}}}
\newcommand{{\underline{M}}}{{\underline{M}}}
\newcommand{{\underline{N}}}{{\underline{N}}}
\newcommand{{\underline{Q}}}{{\underline{Q}}}
\newcommand{{\underline{P}}}{{\underline{P}}}

\newcommand{{\mathbb D}}{{\mathbb D}}
\newcommand{{\mathbb \DB}}{{\mathbb \bar{D}}}
\newcommand{{\mathbb S}}{{\mathbb S}}

\newcommand{{\bf D}}{{\bf D}}
\newcommand{{\bar{\bfD}}}{{\bar{{\bf D}}}}

\newcommand{{\bm D}}{{\bm D}}
\newcommand{{\bar{\bmD}}}{{\bar{{\bm D}}}}

\newcommand{\bm{\nabla}}{\bm{\nabla}}
\newcommand{\bar{\bm{\nabla}}}{\bar{\bm{\nabla}}}

\newcommand{\underline}{\underline}
\newcommand{{\epsilon}}{{\epsilon}}
\newcommand{\bar{\sigma}}{\bar{\sigma}}
\newcommand{{\dot{\alpha}}}{{\dot{\alpha}}}
\newcommand{{\dot{\beta}}}{{\dot{\beta}}}
\newcommand{{\dot{\gamma}}}{{\dot{\gamma}}}
\newcommand{{\dot{\delta}}}{{\dot{\delta}}}
\newcommand{{\dot{\mu}}}{{\dot{\mu}}}
\newcommand{{\dot{\nu}}}{{\dot{\nu}}}
\newcommand{{\dot{\phi}}}{{\dot{\phi}}}
\newcommand{{\bar\theta}}{{\bar\theta}}
\newcommand{{\mathcal N}}{{\mathcal N}}

\def\mathsurround=0pt{\mathsurround=0pt}
\def\fracm#1#2{\hbox{\large{${\frac{{#1}}{{#2}}}$}}}

\def\eqalign#1{\,\vcenter{\openup2\jot \mathsurround=0pt
        \ialign{\strut \hfil$\displaystyle{##}$&$
        \displaystyle{{}##}$\hfil\crcr#1\crcr}}\,}
\newif\ifdtup
\def\panorama{\global\dtuptrue \openup2\jot \mathsurround=0pt
        \everycr{\noalign{\ifdtup \global\dtupfalse
        \vskip-\lineskiplimit \vskip\normallineskiplimit
        \else \penalty\interdisplaylinepenalty \fi}}}
\def\eqalignno#1{\panorama \tabskip=\humongous                  %
eqalignno
        \halign to\displaywidth{\hfil$\displaystyle{##}$
        \tabskip=0pt&$\displaystyle{{}##}$\hfil
        \tabskip=\humongous&\llap{$##$}\tabskip=0pt
        \crcr#1\crcr}}
\def\eqalignnotwo#1{\panorama \tabskip=\humongous
        \halign to\displaywidth{\hfil$\displaystyle{##}$
        \tabskip=0pt&$\displaystyle{{}##}$
        \tabskip=0pt&$\displaystyle{{}##}$\hfil
        \tabskip=\humongous&\llap{$##$}\tabskip=0pt
        \crcr#1\crcr}}

\def\de{{\nabla}}
\def{\bar{\nabla}}}                                         % \bar{del{{\bar{\nabla}}}
\def\sp#1{{}^{#1}}
\def\sb#1{{}_{#1}}
\def{\cal M}{{\cal M}}
\def{\cal R}{{\cal R}}
\def{\cal Y}{{\cal Y}}
\def{\cal F}{{\cal F}}

\def{\bar{\de}}{{\bar{\de}}}

\newcommand{\begin{subequations}}{\begin{subequations}}
\newcommand{\end{subequations}}{\end{subequations}}
\newcommand{\boxedalign}[1]{%
  \[\fbox{%
      \addtolength{\linewidth}{-2\fboxsep}%
      \addtolength{\linewidth}{-2\fboxrule}%
      \begin{minipage}{\linewidth}\vspace{-8pt}%
      \begin{align}#1\end{align}%
      \end{minipage}\nonumber%
    }\]
    }

\newcommand{{\bar i}}{{\bar i}}
\newcommand{{\bar j}}{{\bar j}}
\newcommand{{\bar k}}{{\bar k}}
\newcommand{{\bar l}}{{\bar l}}
\newcommand{{\bar p}}{{\bar p}}
\newcommand{{\bar q}}{{\bar q}}
\newcommand{{\bar 1}}{{\bar 1}}
\newcommand{{\bar 2}}{{\bar 2}}
\newcommand{{\bar 0}}{{\bar 0}}
\newcommand{{\bar n}}{{\bar n}}
\newcommand{{\bar m}}{{\bar m}}
\newcommand{{\bar 4}}{{\bar 4}}

\newcommand{{\rm L}}{{\rm L}}
\newcommand{{\rm R}}{{\rm R}}

\newcommand{{{\bm s}}}{{{\bm s}}}
\newcommand{{{\bar{\mu}}}}{{{\bar{\mu}}}}
\newcommand{{{\bm S}}}{{{\bm S}}}

\newcommand{\mathfrak F}{\mathfrak F}

\newcommand{{\bar{\varphi}}}{{\bar{\varphi}}}

\newcommand{{\bar{\xi}}}{{\bar{\xi}}}
\newcommand{{\bar{\lambda}}}{{\bar{\lambda}}}

\newcommand{{\bm{\xi}}}{{\bm{\xi}}}
\newcommand{{\bm{\xb}}}{{\bm{{\bar{\xi}}}}}

\newcommand{{\bm{\omega}}}{{\bm{\omega}}}

\newcommand{\eqinbox}[1]{%
  \[\fbox{%
      \addtolength{\linewidth}{-2\fboxsep}%
      \addtolength{\linewidth}{-2\fboxrule}%
      \begin{minipage}{\linewidth}\vspace{-8pt}%
      \begin{align}#1\end{align}%
      \end{minipage}\nonumber%
    }\]
    }

\newcommand{\notag \\}{\notag \\}
\newcommand{\varepsilon}{\varepsilon}

\numberwithin{equation}{section}
\newcommand{{\nabla}}{{\nabla}}
\newcommand{{\mathrm{c.c.}}}{{\mathrm{c.c.}}}

\newcommand{{\veps}}{{\varepsilon}}

\newcommand{\mathsf{Sp}}{\mathsf{Sp}}
\newcommand{\mathsf{SU}}{\mathsf{SU}}
\newcommand{\mathsf{SL}}{\mathsf{SL}}
\newcommand{\mathsf{GL}}{\mathsf{GL}}
\newcommand{\mathsf{SO}}{\mathsf{SO}}
\newcommand{\mathsf{O}}{\mathsf{O}}
\newcommand{\mathsf{U}}{\mathsf{U}}
\newcommand{\mathsf{S}}{\mathsf{S}}
\newcommand{\mathsf{PSU}}{\mathsf{PSU}}
\newcommand{\mathsf{PSL}}{\mathsf{PSL}}
\newcommand{\mathsf{OSp}}{\mathsf{OSp}}
\newcommand{\mathsf{Spin}}{\mathsf{Spin}}
\newcommand{\mathsf{Mat}}{\mathsf{Mat}}
\newcommand{\mathsf{ISO}}{\mathsf{ISO}}

\begin{document}

\begin{titlepage}
\begin{flushright}
October, 2021 \\
\end{flushright}
\vspace{5mm}

\begin{center}
{\Large \bf
Local supersymmetry: Variations on a theme by Volkov and Soroka}
\end{center}

\begin{center}

{\bf Sergei M. Kuzenko} \\
\vspace{5mm}

\footnotesize{
{\it Department of Physics M013, The University of Western Australia\\
35 Stirling Highway, Perth W.A. 6009, Australia}}
~\\
\vspace{2mm}
~\\
Email: \texttt{
sergei.kuzenko@uwa.edu.au}\\
\vspace{2mm}

\end{center}

\begin{abstract}
\baselineskip=14pt
As is well known, supergravity was discovered in 1976 by Ferrara, Freedman and van Nieuwenhuizen, and independently by Deser and Zumino. In 1994, Volkov, one of the co-creators of supersymmetry, posted two preprints to the hep-th archive containing ``supergravity before 1976'' in the title. He argued that the key ingredients of $\mathcal{N} = 1$ supergravity in four dimensions had appeared in his earlier work with Soroka on the super-Higgs effect. In order to better appreciate his point of view, in this paper we first revisit the Volkov-Soroka construction and review it. Then we demonstrate that, for specially chosen parameters of the theory, the Volkov-Soroka action is invariant under two different local supersymmetries. One of them is present for arbitrary values of the parameters and acts on the Goldstino, while the other supersymmetry emerges only in a special case and leaves the Goldstino invariant. The former can be used to gauge away the Goldstino, and then the resulting action coincides with that proposed by  Deser and Zumino  for consistent supergravity. In this  sense, pure $\mathcal{N} = 1$ supergravity is a special case of the Volkov-Soroka theory, and it could have been discovered before 1976.
\end{abstract}
\vspace{5mm}

\vfill

\vfill
\end{titlepage}

\newpage
\renewcommand{\thefootnote}{\arabic{footnote}}
\setcounter{footnote}{0}

\allowdisplaybreaks

\section{Introduction}

Dmitry V. Volkov was one of the co-discoverers of supersymmetry. His most prominent results on rigid and local supersymmetry are the Goldstino model (jointly with Akulov) \cite{VA,AV} and the super-Higgs mechanism (jointly with Soroka) \cite{VS,VS2}. His approach to nonlinear realisations of internal and space-time symmetries \cite{Volkov73}, which paved the way to \cite{VA,AV,VS,VS2},
has also been highly influential. However, it is less known that Ref. \cite{AV} also pioneered the following fundamental concepts of modern theoretical physics: (i)  the $\cN$-extended super-Poincar\'e group and, hence,  the $\cN$-extended super-Poincar\'e algebra\footnote{The $\cN=1$ super-Poincar\'e algebra
was discovered in 1971 by Golfand and Likhtman
\cite{GL}.} for $\cN>1$; and (ii) $\cN$-extended Minkowski superspace.\footnote{The Akulov-Volkov paper \cite{AV} was submitted to the journal {\it Theoretical and Mathematical Physics} on 8 January 1973 and published in January 1974, before the publication of the first paper on supersymmetry by Wess and Zumino \cite{WZ} and long before the work by Salam and Strathdee \cite{SS} devoted to the $\cN=1$ superspace approach. It remains largely unknown, perhaps because
it was published in a Russian journal. }
It is quite remarkable that all these results had appeared before
the first paper by Wess and Zumino on supersymmetry \cite{WZ} was published on 18 February, 1974.

As is well known, supergravity was discovered in 1976 \cite{FvNF,DZ}. In 1994, less than two years before his death, Volkov posted two preprints to the hep-th archive \cite{Volkov1,Volkov2}, both of which contained ``supergravity before 1976'' in the title.
He argued that the crucial ingredients of $\cN=1$ supergravity in four dimensions had appeared in \cite{VS,VS2}.\footnote{Similar thoughts were also expressed by Soroka a  few years later \cite{Soroka1,Soroka2}, see also \cite{Nurmagambetov, Dup}.} In order to better appreciate his point of view, in this paper we first provide a detailed review of the Volkov-Soroka construction. Then we demonstrate that for specially chosen parameters the theory possesses two local supersymmetries. One of them can be used to gauge away the Goldstino, and then the resulting action coincides with that proposed by  Deser and Zumino  to describe  consistent supergravity. This shows that  pure $\mathcal{N} = 1$ supergravity is a special case of the Volkov-Soroka theory, and it could have been discovered before 1976.

\section{A review of the Volkov-Soroka construction}

This section is devoted to a pedagogical review of the Volkov-Soroka construction
 \cite{VS,VS2}.

Let  ${\mathfrak P}{(4|\cN )}$ denote the four-dimensional
$\cN$-extended super-Poincar\'e group introduced in \cite{AV}.
Any element $g \in {\mathfrak  P}{(4 | \cN )}$ is a $(4+\cN)\times (4+\cN)$
{supermatrix} of the form\footnote{Our parametrisation of the elements of
${\mathfrak P}{(4|\cN )}$ follows \cite{Kuzenko2010}.}
\begin{subequations}\label{SP}
\begin{eqnarray}
g &=& S(b,  \ve) \,h (M,U)  \equiv S h ~,
\label{SP1} \\
S(b, \ve ) &:= &
\left(
\begin{array}{c | c ||c}
  \mathbbm{1}_2  ~& ~ 0&~{ 0}   \\
\hline
-{\rm i}\,\tilde{ b}_{(+)}  \phantom{\Big|}  ~& ~\mathbbm{1}_2 ~&~ {  2{ \epsilon}^\dagger} \\
\hline
\hline
{ 2\epsilon} ~& ~{ 0}&~\mathbbm{1}_\cN
\end{array}
\right)
=
\left(
\begin{array}{r | c ||c}
\delta_\alpha{}^\beta  ~& ~ 0&~{ 0}   \\
\hline
-{\rm i}\,b^{{\dot \alpha} \beta}_{(+)}  ~& ~\delta^{\dot \alpha}{}_{\dot \beta} &~ { 2{\bar \epsilon}^{{\dot \alpha} j}}
\label{SP-g} \\
\hline
\hline
{2\epsilon_i{}^\beta} ~& ~{ 0}&~\delta_i{}^j
\end{array}
\right) ~,
\label{SP2}  \\
{ h(M,U)} &:=&
\left(
\begin{array}{c | c ||c}
 M  ~& ~ 0&~{ 0}   \\
\hline
0  ~& ~(M^{-1})^\dagger&~ { 0} \\
\hline
\hline
{ 0} ~& ~{ 0}&~U
\end{array}
\right)
= \left(
\begin{array}{c | c| |c}
M_\alpha{}^\beta  ~& ~ 0&~{ 0}   \\
\hline
0  ~& ({\bar M}^{-1})_{\dot \beta}{}^{\dot \alpha}  &~ {  0} \\
\hline
\hline
{  0} ~& ~{  0}&~U_i{}^j
\end{array}
\right) ~,
\label{SP3}
\end{eqnarray}
\end{subequations}
where $M =(M_\alpha{}^\beta)  \in \mathsf{SL}(2,{\mathbb C})$,
$U =(U_i{}^j) \in \mathsf{U}(\cN)$, and
\begin{eqnarray}
\tilde{ b}_{(\pm)} = \tilde{ b} \pm 2{\rm i} \epsilon^\dagger \epsilon~.
\label{2.2}
\end{eqnarray}
The tilde notation in \eqref{SP2} and \eqref{2.2} reflects the fact that there are two types of relativistic Pauli matrices, $\sigma_a$ and $\tilde \sigma_a$, see the appendix.
The group element $S(b, \ve)$
is labelled  by
4 commuting real  parameters $b^a$ and
$4\cN $ anti-commuting
complex parameters $\ve = (\epsilon , \epsilon^\dagger)$, where
$\epsilon=( \epsilon_i{}^\alpha)$  and
$\epsilon^\dagger =  ({\bar \epsilon}^{{\dot \alpha} i} )$,
$  {\bar \epsilon}^{{\dot \alpha} i} :=
\overline{\, \epsilon_i{}^\alpha \,}$. In the vector notation, eq. \eqref{2.2} reads
\begin{eqnarray}
 b^a_{(\pm )}:= b^a \pm {\rm i} \,\epsilon_i \sigma^a {\bar \epsilon}^i
=b^a \pm {\rm i}\, \epsilon_i{}^\alpha (\sigma^a)_{\alpha \ad} {\bar \epsilon}^{{\dot \alpha} i}  ~.
\label{+-}
\end{eqnarray}

Introduce Goldstone fields
$Z^A (x) = \big(X^a (x) , \Theta_i{}^\alpha (x) , \bar \Theta^{\ad i} (x) \big)$
for spacetime translations $\big(X^a \big)$ and supersymmetry transformations
$\big( \Theta_i{}^\alpha  , \bar \Theta^{\ad i}  \big)$.
They parametrise the homogeneous space
($\cN$-extended Minkowski superspace)
\begin{eqnarray}
{\mathbb M}^{4|4\cN}
=\frac{ {\mathfrak P}(4|\cN) }{ \mathsf{SL}(2,{\mathbb C}) \times \mathsf{U}(\cN) }
\label{2.4}
\end{eqnarray}
according to the rule:
\begin{eqnarray}
{\mathfrak S}(Z ) &= &
\left(
\begin{array}{c | c ||c}
  \mathbbm{1}_2  ~& ~ 0&~{ 0}   \\
\hline
-{\rm i}\,\tilde{ X}_{(+)}  \phantom{\Big|}  ~& ~\mathbbm{1}_2 ~&~ {  2{ \Theta}^\dagger} \\
\hline
\hline
{ 2\Theta} ~& ~{ 0}&~\mathbbm{1}_\cN
\end{array}
\right)
~ \implies ~
{\mathfrak S}^{-1} (Z ) =
\left(
\begin{array}{c | c ||c}
  \mathbbm{1}_2  ~& ~ 0&~{ 0}   \\
\hline
{\rm i}\,\tilde{ X}_{(-)}  \phantom{\Big|}  ~& ~\mathbbm{1}_2 ~&~  - 2{ \Theta}^\dagger \\
\hline
\hline
- 2\Theta ~& ~{ 0}&~\mathbbm{1}_\cN
\end{array}
\right)~,~~~
\end{eqnarray}
where
\begin{eqnarray}
\tilde{ X}_{(\pm)} = \tilde{ X} \pm 2{\rm i} \Theta^\dagger \Theta~.
\end{eqnarray}

We define gauge super-Poincar\'e transformations by
\begin{eqnarray}
g(x): Z(x) \to Z'(x) ~, \qquad g {\mathfrak S}(Z) = {\mathfrak S}(Z') h~,
\end{eqnarray}
with $g=Sh$.
This is equivalent to the following transformations of the Goldstone fields:
\begin{subequations}\label{2.8}
\begin{eqnarray}
 S(b,  \ve ) :\qquad
 \tilde X' &=& \tilde X + \tilde b +2{\rm i} (\epsilon^\dagger \Theta- \Theta^\dagger \epsilon) ~, \label{2.8a}\\
 \Theta' &=& \Theta +\epsilon \label{2.8b}
 \end{eqnarray}
 \end{subequations}
and
\begin{subequations} \label{2.9}
\begin{eqnarray}
h(M, U): \qquad \tilde{X}' &=& = (M^\dagger)^{-1} \tilde{X} M^{-1} ~,\\
\Theta'&=& U \Theta M^{-1} ~.
\end{eqnarray}
\end{subequations}

Introduce a connection $ {\mathfrak A}  = {\rm d} x^m {\mathfrak A}_m $ taking its values in the super-Poincar\'e algebra,
\begin{eqnarray}
{\mathfrak A}  &:= &
\left(
\begin{array}{c | c ||c}
  \Omega  ~& ~ 0&~{ 0}   \\
\hline
-{\rm i}\,\tilde{ e}  \phantom{\Big|}  ~& ~-\Omega^\dagger ~&~ {  2{ \psi}^\dagger} \\
\hline
\hline
{ 2\psi} ~& ~{ 0}&~{\rm i} V
\end{array}
\right)
=
\left(
\begin{array}{r | c ||c}
\Omega_\alpha{}^\beta  ~& ~ 0&~{ 0}   \\
\hline
-{\rm i}\,e^{{\dot \alpha} \beta}  ~& ~ -\bar \Omega^\ad{}_\bd &~ { 2{\bar \psi}^{{\dot \alpha} j}} \\
\hline
\hline
{2\psi_i{}^\beta} ~& ~{ 0}&~ {\rm i} V_i{}^j
\end{array}
\right) ~,
\end{eqnarray}
and possessing the gauge transformation law
\begin{eqnarray}
 {\mathfrak A}' = g  {\mathfrak A} g^{-1} + g {\rm d} g^{-1}~.
 \end{eqnarray}
 Here the one-forms $\Omega_\alpha{}^\beta $ and $\bar \Omega^\ad{}_\bd $ are the spinor counterparts of the Lorentz connection $\Omega^{ab} = {\rm d} x^m \Omega_m{}^{ab} = - \Omega^{ba} $ such that
 \begin{eqnarray}
 \Omega_\alpha{}^\beta = \frac12 (\sigma^{ab})_\alpha{}^\beta\, \Omega_{ab} ~, \qquad
\bar  \Omega^\ad{}_\bd =- \frac12 (\tilde \sigma^{ab})^\ad{}_\bd \,\Omega_{ab} ~.
 \end{eqnarray}
 The Lorentz connection is an independent field. It is expressed in terms of the other fields on the mass shell.
 The one-form $e^{{\dot \alpha} \beta}$ is the spinor counterpart of the vierbein
 $e^a = {\rm d} x^m e_m{}^a$.  The fermionic one-forms $\psi_i{}^\beta $
 and ${\bar \psi}^{ \ad j} $ describe $\cN$ gravitini. Finally, the anti-Hermitian one-form
 ${\rm i} V$ is the $\mathsf{U}(\cN)$ gauge field.

Associated with ${\mathfrak S}$ and ${\mathfrak A}$ is another connection
\begin{eqnarray}
{\mathbb A} := {\mathfrak S}^{-1}{\mathfrak A} {\mathfrak S}
+ {\mathfrak S}^{-1} {\rm d} {\mathfrak S}
~,
\end{eqnarray}
which is characterised by the gauge transformation law
\begin{eqnarray}
{\mathbb A}' = h {\mathbb A} h^{-1} + h \,{\rm d} h^{-1}~,
\label{2.14}
\end{eqnarray}
for an arbitrary gauge parameter $g = Sh$.
This transformation law tells us that $\mathbb A$ is invariant under all gauge
transformations of the form $g = S(b, \ve)$ which describe local spacetime translations and supersymmetry transformations.
The connection ${\mathbb A} $ is the main object of the Volkov-Soroka construction.
It has the form
\begin{eqnarray}
{\mathbb A}  &:= &
\left(
\begin{array}{c | c ||c}
  \Omega  ~& ~ 0&~{ 0}   \\
\hline
-{\rm i}\,\tilde{ E}  \phantom{\Big|}  ~& ~-\Omega^\dagger ~&~ {  2{ \Psi}^\dagger} \\
\hline
\hline
{ 2\Psi} ~& ~{ 0}&~ {\rm i} V
\end{array}
\right)~,
\end{eqnarray}
where we have defined
\begin{subequations}\label{2.15}
\begin{eqnarray}
\Psi &:=& \psi +\cD \Theta~,  \label{2.15a} \\
\Psi^\dagger &:=& \psi^\dagger +\cD \Theta^\dagger~,
\label{2.15.b} \\
\tilde{ E} &:=& \tilde{ e} + \cD \tilde{ X} + 4{\rm i} (\Psi^\dagger - \frac12 \cD \Theta^\dagger) \Theta
- 4{\rm i} \Theta^\dagger ( \Psi -\frac12 \cD \Theta)~,
\label{2.15c}
\end{eqnarray}
\end{subequations}
and $\cD$ denotes the covariant derivative,
\begin{subequations}
\begin{eqnarray}
\cD \Theta &=& {\rm d} \Theta - \Theta \Omega + {\rm i} V \Theta~, \\
\cD \Theta^\dagger &=& {\rm d} \Theta^\dagger - \Omega^\dagger \Theta^\dagger - {\rm i}  \Theta^\dagger V ~,\\
\cD \tilde{ X} &=& {\rm d}  \tilde{ X} -\Omega^\dagger  \tilde{ X} -  \tilde{ X} \Omega~.
 \end{eqnarray}
 \end{subequations}

Equation \eqref{2.14} is equivalent to the following gauge transformation laws:
\begin{subequations}
\begin{eqnarray}
\Omega' &=&  M \Omega M^{-1}  +M {\rm d} M^{-1} ~,\\
{\rm i} V'&=& U ({\rm i} V )U^{-1}  +  U {\rm d} U^{-1}
\end{eqnarray}
\end{subequations}
and
\begin{subequations}
\begin{eqnarray}
\tilde{E}' &=&  (M^\dagger)^{-1} \tilde{E} M^{-1} ~,\\
\Psi'&=& U \Psi M^{-1} ~,
\end{eqnarray}
\end{subequations}
We see that the supersymmetric one-forms $E^a$ and $\Psi_i{}^\beta$ transform as tensors with respect to the Lorentz and $\mathsf{U}(\cN)$ gauge groups.

Making use of \eqref{2.8}, we deduce the local supersymmetry transformation
of the gravitini and the vielbein
\begin{subequations}
\begin{eqnarray}
\psi' &=& \psi - \cD \epsilon ~,\\
\tilde{e}' &=& = \tilde e + 4 {\rm i} (\epsilon^\dagger \psi - \psi^\dagger \epsilon)
+{\rm i} (\cD \epsilon^\dagger \epsilon - \epsilon^\dagger \cD \epsilon) ~.
\end{eqnarray}
\end{subequations}
In the infinitesimal case, this transformation can be rewritten in the form
\begin{eqnarray}
\delta_\ve \psi = - \cD \epsilon~, \qquad
\delta_\ve e^a =  2{\rm i} \, {\rm tr} \big[ \sigma^a (\psi^\dagger \epsilon  - \epsilon^\dagger  \psi ) \big] ~.
\label{2.21}
\end{eqnarray}
As pointed out by Volkov \cite{Volkov1}, these transformation laws  coincide with those used by Deser and Zumino in their construction of $\cN=1$ supergravity \cite{DZ}.
We should remark that the supersymmetry transformations of the Goldstone fields $X^a$ and $\Theta_i{}^\beta$ are given by the relations \eqref{2.8}.

Let us consider a local Poincar\'e translation, $S(b,0)$. It only acts on the Goldstone vector field $X^a$ and the vierbein $e^a$,
\begin{eqnarray}
X'{}^a = X^a + b^a~, \qquad e'{}^a = e^a - \cD b^a~.
\label{2.22}
\end{eqnarray}
We have two types of gauge transformations with vector-like parameters, the general coordinates transformations and the local Poincar\'e translations. The latter gauge freedom can be fixed by imposing the condition $X^a =0$, and then we stay only with the general coordinate invariance. However, in what follows we will keep $X^a$ intact.

The curvature tensor is given by
\begin{eqnarray}
{\mathbb R}={\rm d} {\mathbb A} - {\mathbb A} \wedge {\mathbb A}~, \qquad
{\mathbb R}' = h {\mathbb R} h^{-1}~.
\end{eqnarray}
Its explicit form is
\begin{eqnarray}
{\mathbb R}  &:= &
\left(
\begin{array}{c | c ||c}
  R  ~& ~ 0&~{ 0}   \\
\hline
-{\rm i}\,\tilde{ \mathbb T}  \phantom{\Big|}  ~& ~-R^\dagger ~&~   2 \cD \Psi^\dagger \\
\hline
\hline
 2\cD \Psi ~& ~{ 0}&~ {\rm i} F
\end{array}
\right)~,
\end{eqnarray}
where $R = (R_\alpha{}^\beta) $ and $R^\dagger = (\bar R^\ad{}_\bd )$
form the Lorentz curvature, $F=(F_i{}^j)$ is the
Yang-Mills  field strength,
\begin{subequations}
\begin{eqnarray}
\cD \Psi &=& {\rm d} \Psi - \Psi \wedge \Omega - {\rm i} V \wedge \Psi~,
\\
\cD \Psi^\dagger &=& {\rm d} \Psi^\dagger +\Omega^\dagger \wedge \Psi^\dagger
- {\rm i} \Psi^\dagger \wedge V
\end{eqnarray}
\end{subequations}
are the gravitino field strengths, and
\begin{eqnarray}
\tilde {\mathbb T} = {\rm d} \tilde E - \tilde E \wedge \Omega + \Omega^\dagger \wedge \tilde E
- 4{\rm i} \Psi^\dagger \wedge \Psi = \cD \tilde E - 4{\rm i} \Psi^\dagger \wedge \Psi
\end{eqnarray}
is the supersymmetric torsion tensor.  In the vector notation, the torsion tensor reads
\begin{eqnarray}
{\mathbb T}^a = \cD E^a + 2{\rm i} \Psi \wedge \sigma^a \bar \Psi~.
\label{torsion}
\end{eqnarray}

It should be pointed out that the exterior derivative is defined to obey the property
\begin{eqnarray}
{\rm d} \Big( \Sigma_p \wedge \Sigma_q \Big) = \Sigma_p \wedge {\rm d} \Sigma_q
+(-1)^q {\rm d} \Sigma_p \wedge \Sigma_q~,
\end{eqnarray}
which is used for superforms \cite{WB}.

The above results allow one to construct gauge-invariant functionals that can be used to construct a locally supersymmetric action. The invariants proposed in \cite{VS,VS2} are the following:
\begin{itemize}
\item The Einstein-Hilbert action
\begin{eqnarray}
S_{\rm EH} = \frac 14 \int \ve_{abcd} E^a \wedge E^b \wedge R^{cd} ~;
\label{EH}
\end{eqnarray}
\item The Rarita-Schwinger action
\begin{eqnarray}
S_{\rm RS} = \frac12 \int \Big( \Psi_i \wedge E^a \wedge \sigma_a \cD \bar \Psi^i - \cD \Psi_i \wedge e^a \wedge \sigma_a \bar \Psi^i \Big) ~;
\label{RS}
\end{eqnarray}
\item The cosmological term
\begin{eqnarray}
S_{\rm cosmological } = \frac{1}{24} \int \ve_{abcd} E^a \wedge E^b \wedge E^{c} \wedge E^d ~.
\label{2.31}
\end{eqnarray}
\item $\mathsf{O}(\cN)$-invariant mass term
\begin{eqnarray}
S_{\rm mass} = \frac{{\rm i}}{4}  \int  E^a \wedge E^b \wedge
\Big(\delta^{ij}  \Psi_i \wedge \sigma_{ab} \Psi_j  - \delta_{ij} \bar \Psi ^i\wedge \tilde{\sigma}_{ab} \bar \Psi^j \Big) ~.
\label{2.32}
\end{eqnarray}
\end{itemize}
The functionals $S_{\rm EH} $, $S_{\rm RS} $ and $S_{\rm cosmological}$ are $\mathsf{U}(\cN)$ invariant.  The cosmological term, eq. \eqref{2.31},  also  contains the kinetic term for the Goldstini \cite{VA,AV}.
The mass term \eqref{2.32} is invariant under local internal transformations only if the group $\mathsf{U}(\cN)$ is replaced with $\mathsf{O}(\cN)$, and the gauge connection ${\rm i} V$ takes its values in the Lie algebra $\mathfrak{so}(\cN)$. The Yang-Mills action associated with $V$ is obviously supersymmetric, but it will not be used  in what follows.

In the  $\cN=1$ case, the Volkov-Soroka theory is described by the general action
\begin{eqnarray}
S= S_{\rm EH} +4c \,S_{\rm RS} + 4m\, S_{\rm mass} + \lambda S_{\rm cosmological} ~,
\label{2.33}
\end{eqnarray}
with $c$, $m$ and $\lambda$ coupling constants.

\section{Second local supersymmetry}

In this section our consideration is restricted to the $\cN=1$ case, and  the $\mathsf{U}(1)$
gauge field is switched off, $V=0$. We are going to  demonstrate that the action
\begin{eqnarray}
S_{\rm SUGRA} = S_{\rm EH} + 4 S_{\rm RS}
\label{3.1}
\end{eqnarray}
is invariant under
a new local supersymmetry transformation described by the parameter  ${\bm \ve} = ({\bm \epsilon}^\alpha , \bar{\bm \epsilon}^\ad)$.
It acts on the composite fields
$\Psi^\alpha$ and $E^a$, defined by eq. \eqref{2.15}, and the Lorentz connection as follows
\begin{subequations} \label{3.2}
\begin{eqnarray}
\delta_{\bm \ve} \Psi^\alpha &=& - \cD {\bm \epsilon}^\alpha~, \qquad
\delta_{\bm \ve} E^a =  2{\rm i} \, (\Psi \sigma^a \bar {\bm \epsilon}   - {\bm \epsilon} \sigma^a  \Psi )  ~, \label{3.2a} \\
  \frac 14  \ve_{abcd} \delta_{\bm \ve} \Omega^{bc} \wedge E^d &=&
  -{\bm \epsilon} \sigma_a \cD \bar \Psi
  -  \cD \Psi \sigma_a \bar{\bm \epsilon} ~. \label{3.2b}
\end{eqnarray}
The Goldstone fields are inert under this transformation,
\begin{eqnarray}
\delta_{{\bm \ve}}  X^a &=& 0~, \qquad \qquad \delta_{ {\bm \ve}} \Theta^\alpha =0~. \label{3.2c}
\end{eqnarray}
The elementary field $\psi^\alpha $ and $e^a$ transform as follows:
\begin{eqnarray}
\delta_{\bm \ve} \psi^\alpha &=& - \cD {\bm \epsilon}^\alpha +\Theta^\beta \delta_{\bm \ve} \Omega_\beta{}^\alpha~, \label{3.2d}\\
\delta_{\bm \ve} e^a &=&  2{\rm i} \big( \Psi \sigma^a \bar {\bm \epsilon} - {\bm \epsilon} \sigma^a \bar \Psi \big)
 +2{\rm i} \big( \Theta \sigma^a \cD \bar {\bm \epsilon} - \cD {\bm \epsilon} \sigma^a \bar \Theta \big)
 -\delta_{\bm \ve} \Omega^a{}_b X^b\nonumber \\
 && + \frac12 \ve^{abcd} \delta_{\bm \ve} \Omega_{bc} \Theta\sigma_d \bar \Theta~.
 \label{3.2e}
\end{eqnarray}
\end{subequations}
The dependence on $\delta_{\bm \ve} \Omega$ in \eqref{3.2d} and \eqref{3.2e} is such that
the composite fields $\Psi^\alpha$ and $E^a$ remain unchanged when the connection gets the displacement $\Omega \to \Omega + \delta_{\bm \ve} \Omega$. We should point out  that the action \eqref{3.1}
involves only  the one-forms  $\Psi^\alpha$, $\bar \Psi^\ad$, $E^a$ and $\Omega^{ab}$ and their descendants. We also remark that the transformations \eqref{3.2a} and \eqref{3.2b} reduce to those given by Deser and Zumino \cite{DZ} if the Goldstone fields
$X^a$ and $\Theta^\alpha$ are switched off.

It should be pointed out that eq. \eqref{3.2b} uniquely determines
$\delta_{\bm \ve} \Omega^{bc}$. Indeed, given a vector-valued two-form
\begin{eqnarray}
\Sigma_a = \frac12 E^c \wedge E^b \Sigma_{a,bc} ~,
\end{eqnarray}
the following equation
\begin{eqnarray}
\frac12 \ve_{abcd}  \omega^{bc} \wedge E^d \equiv \widetilde{\omega}_{ab} \wedge E^b = \Sigma_a~,
\qquad  \widetilde{\omega}_{ab} = E^c \widetilde{\omega}_{c,ab}
\end{eqnarray}
on the one-form ${\omega}^{ab} = E^c {\omega}_c{}^{ab} $, which takes its values in the Lorentz algebra,
has the unique solution
\begin{eqnarray}
\widetilde{\omega}_{c,ab} = \frac12 \Big( \Sigma_{a,bc} - \Sigma_{b,ac} - \Sigma_{c,ab} \Big)~.
\end{eqnarray}

Now we turn to demonstrating that the action \eqref{3.1} is invariant under the supersymmetry transformation \eqref{3.2}. Let us first take into account the variations
in \eqref{3.2a}.
Varying the Einstein-Hilbert action \eqref{EH} for $\bm \epsilon \neq 0$ and $\bar {\bm \epsilon} =0$
gives
\begin{eqnarray}
\delta^{(1)}_{\bm \epsilon} S_{\rm EH} = -{\rm i}  \int \ve_{abcd}
R^{ab} \wedge E^c \wedge {\bm \epsilon} \sigma^d \bar \Psi   ~.
\label{3.3}
\end{eqnarray}
Varying the Rarita-Schwinger action \eqref{RS} gives
\begin{eqnarray}
\delta^{(1)}_{\bm \epsilon} S_{\rm RS} &=&  \int \Big\{
\frac12  \cD  E^a \wedge {\bm \epsilon} \sigma_a \cD \bar \Psi
+{\rm i} {\bm \epsilon}\sigma^a \bar \Psi \,\Psi \sigma_a \cD \bar \Psi
\Big\}
 \nonumber \\
&&
+\frac12 \int \Big\{ E^a \wedge {\bm \epsilon} \sigma_a R^\dagger \bar \Psi
-E^a \wedge {\bm \epsilon} R \wedge \sigma^a \bar \Psi
\Big\} ~,
\label{3.4}
\end{eqnarray}
where we have used the relations
\begin{eqnarray}
\cD \cD {\bm \epsilon} = - {\bm \epsilon} R~, \qquad
\cD \cD \Psi^\dagger = - R^\dagger \wedge \Psi^\dagger~.
\end{eqnarray}
Here $\delta^{(1)}$ means that we have taken into account only the variations \eqref{3.2a}.

Making use of the identities \eqref{A.14}, the curvature-dependent contributions in \eqref{3.4} can be rearranged as
\begin{eqnarray}
\frac12 \int \Big\{ E^a \wedge {\bm \epsilon} \sigma_a R^\dagger \bar \Psi
-E^a \wedge {\bm \epsilon} R \wedge \sigma^a \bar \Psi
\Big\} =
\frac{{\rm i}}{4}  \int \ve_{abcd}
R^{ab} \wedge E^c \wedge {\bm \epsilon} \sigma^d \bar \Psi  ~.
\label{3.6}
\end{eqnarray}
In order for the curvature contributions in \eqref{3.3} and \eqref{3.4}  to cancel each other, we must consider the action \eqref{3.1}, for which we obtain
\begin{eqnarray}
\delta^{(1)}_{\bm \epsilon} (S_{\rm EH} + 4 S_{\rm RS} )
=  2 \int \Big(
  \cD  E^a + 2{\rm i} \Psi \sigma^a \bar \Psi \Big)
  \wedge {\bm \epsilon} \sigma_a \cD \bar \Psi
  =2 \int  {\bm \epsilon} \sigma_a \cD \bar \Psi \wedge {\mathbb T}^a
  ~,
  \end{eqnarray}
with ${\mathbb T}^a$ being the supersymmetric torsion tensor
\eqref{torsion}. Adding the complex conjugate part gives
\begin{eqnarray}
\delta^{(1)}_{\bm \ve} (S_{\rm EH} + 4 S_{\rm RS} )
  = 2\int \Big(  {\bm \epsilon} \sigma_a \cD \bar \Psi
  +  \cD \Psi \sigma_a \bar{\bm \epsilon} \Big)\wedge {\mathbb T}^a
  ~.
  \label{3.8}
  \end{eqnarray}

Next, let us vary the action \eqref{3.1} with respect to the Lorentz connection $\Omega^{ab}$.
We give the Lorentz connection a small disturbance, $ \Omega \to \Omega + \delta_{\bm \ve}  \Omega $,
with $  \delta_{\bm \ve} \Omega $ to be determined below, and assume that the fields  $\psi^\alpha$ and $e^a$ also acquire $  \delta_{\bm \ve} \Omega $-dependent variations given in \eqref{3.2d} and \eqref{3.2e}.
We denote $\delta^{(2)}$ the corresponding variation.
Direct calculations give
\begin{eqnarray}
\delta^{(2)}_{\bm \ve} (S_{\rm EH} + 4 S_{\rm RS} )
=  \frac12 \int \ve_{abcd} \delta_{\bm \ve} \Omega^{bc} \wedge E^d \wedge {\mathbb T}^a~.
\label{3.9}
\end{eqnarray}

Combining the results \eqref{3.8} and \eqref{3.9}, we end up with
\begin{eqnarray}
\delta_{\bm \ve}  (S_{\rm EH} + 4 S_{\rm RS} )
=2\int \Big(  {\bm \epsilon} \sigma_a \cD \bar \Psi
  +  \cD \Psi \sigma_a \bar{\bm \epsilon}
  +\frac 14  \ve_{abcd} \delta_{\bm \ve} \Omega^{bc} \wedge E^d
  \Big)\wedge {\mathbb T}^a ~.
\end{eqnarray}
This variation vanishes if $\delta_{\bm \ve} \Omega$ is given by eq. \eqref{3.2b}.

We have demonstrated that the theory \eqref{3.1} has two types of local supersymmetry.
The original Volkov-Soroka supersymmetry is described by the relations
\begin{subequations} \label{3.11}
\begin{eqnarray}
\delta_{ \ve} \psi^\alpha &=& - \cD { \epsilon}^\alpha~, \qquad
\delta_{\ve} e^a =  2{\rm i} \, (\psi \sigma^a \bar { \epsilon}   - { \epsilon} \sigma^a \bar  \psi )  ~, \\
\delta_{ { \ve}} \Theta^\alpha &=&\epsilon^\alpha ~,\qquad
\delta_{{ \ve}}  X^a = {\rm i} (\Theta \sigma^a \bar \epsilon-  \epsilon \sigma^a \bar \Theta)~.
\end{eqnarray}
\end{subequations}
The second supersymmetry, introduced in this work,
 is given by eq. \eqref{3.2}. The gauge transformations
\eqref{2.22} and  \eqref{3.11} allow us to gauge away the Goldstone fields
$Z^A (x) = \big(X^a (x) , \Theta^\alpha (x) , \bar \Theta^{\ad } (x) \big)$
by
imposing the conditions
\begin{eqnarray}
X^a=0~, \qquad
\Theta^\alpha=0~.
\end{eqnarray}
As a result, the action \eqref{3.1} turns into the supergravity action proposed by Deser and Zumino \cite{DZ}, and the local supersymmetry transformation laws \eqref{3.2a} turn into those given in \cite{DZ}.

\section{Conclusion}

 A few years ago, it was shown \cite{BMST} that the Volkov-Soroka theory
 \eqref{2.33},
with non-vanishing parameters $c$, $m$ and $\lambda$,
is equivalent to spontaneously broken $\cN=1$ supergravity \cite{DZ2,LR}
which was
called de Sitter supergravity in \cite{BFKVP}.

In this paper we have demonstrated that the action \eqref{3.1} is a St\"uckelberg-type extension of the $\cN=1$ supergravity action proposed by Deser and Zumino.
Equivalently, the latter theory is a gauged-fixed version of \eqref{3.1}.
Therefore,  pure $\mathcal{N} = 1$ supergravity is a special case of the Volkov-Soroka theory, and it could have been discovered before 1976.

The new local supersymmetry \eqref{3.2} of the action \eqref{3.1} is the main original result of this paper.
\\

\noindent
{\bf Acknowledgements:}\\
I thank Evgeny Buchbinder, Ian McArthur and Dmitri Sorokin for discussions, and Darren Grasso and Michael Ponds for comments on the manuscript.
This work was supported in part by the Australian Government through
the Australian Research Council, project No. DP200101944.


\section{Introduction}\label{sec:Intro}

With constant increase of computational power data-driven methods have gained a significant popularity in the past few decades with applications spanning all fields of science \cite{Montans19}. The study of {\it intrinsic localized modes} (ILM) in crystal lattice models, also known as {\it discrete breathers} (DB), is an active field of research in mathematics and physics \cite{MacKay94,Flach08,ProcMica14,ram20} which relies heavily on obtained data from numerical simulations. Thus, emerging data-driven methods provide great opportunity to explore new techniques for understanding nonlinear wave phenomena in crystal lattice models.

Such localized wave solutions are of particular interest in material science where transport of charge in silicates by moving nonlinear localized modes is experimentally confirmed, the phenomenon called {\it hyperconductivity} \cite{ram20,Russell19}. In addition, localized modes provide mechanisms for the formation of long decorated dark lines in muscovite mica \cite{ProcMica14} and can have a functional role in proteins \cite{Senet15,Piazza08}. Strong evidence of localized phenomena in crystal lattice models is provided by numerical experiments \cite{Chris11,Marin01,Bajars15,Shepelev20,Chetverikov16,Bajars21}, to mention but a few.

There are already well established theoretical and numerical understandings of the existence of one-dimensional stationary DB \cite{MacKay94,Flach08,Aubry06}, i.e., spatially localized time-periodic excitations. The same may not be said about propagating DB \cite{Flach99} and, in particular, in higher dimensions. There are still open theoretical questions regarding the existence and properties of propagating DB in general nonintegrable lattices. Mobile discrete breathers are not expected to be long-lived, since the lattice models considered are likely to be nonintegrable and their lifespan is subject to interactions with lattice defects and the phonon background, which may be viewed as thermal noise. Thus, an obvious challenge is to study the existence and interactions of long-lived propagating discrete breather solutions \cite{Bajars15,Bajars21} to understand their role in the transport of energy in sputtering experiments \cite{Russell07}. Such demanding computations require good data-based numerical tools to detect and track breather solutions in space and time which serves as a motivation for our work.

In this work we develop classification algorithms to differentiate between nonlinear localized waves and nonlocalized linear waves in numerical simulations of one-dimensional crystal lattice model. Obtained classifiers rely on locally sampled data in opposed to nonlocal time series analysis methods, such as discrete Fourier and wavelet transforms. Such classifiers can be efficiently trained with different given labeled datasets which are not only limited to numerical simulation data. Trained classifiers then can be used to detect localization regions, e.g., in numerical simulations, which is the first step towards fully automated tool for quantitative data-based analysis of complex numerical experiments. Importantly, our analysis and methodology, which follows, extends, in general, to any one-dimensional crystal lattice models which support intrinsic localized mode solutions.

The paper is organized in the following way. In Section \ref{sec:MathModel} mathematical model of one-dimensional crystal lattice is formulated. Spectral properties of discrete lattice waves are discussed in Section \ref{sec:Spectrum}. Dimensionality reduction methods and classification algorithms are presented in Section \ref{sec:Classification}. In Section \ref{sec:Applications} application of the classification algorithms is demonstrated for detecting localized modes in numerical simulations. Discussion and conclusions are provided in Section \ref{sec:Conclusions}.

\section{Mathematical model}\label{sec:MathModel}

In this work we are concerned with one-dimensional crystal lattice models of classical molecular dynamics exhibiting ILM solutions. In particular, we consider Hamiltonian dynamics of $N$ particles arising from the dimensionless Hamiltonian
\begin{equation}\label{eq:Hamilt}
H = K + U + V
  = \sum_{n=1}^{N} \left( \frac{1}{2} \dot{q}_{n}^2 +
  U(q_{n}) + V(|q_{n+1}-q_{n}|) \right),
\end{equation}
where $K$ is the kinetic energy, $U$ is the on-site potential energy, $V$ is the interaction potential of particles, $q_n\in\mathbb{R}$ is the position of the $n^{th}$ particle, $\dot{q}_{n}$ is its time derivative and $|\cdot|$ is the Euclidean distance. In the dimensionless form \eqref{eq:Hamilt} particle equilibrium distance $\Delta=q_{n+1}-q_n$ as well as all masses are equal to one.

For example, such model \eqref{eq:Hamilt} arises as one-dimensional model of the K-K layer of layered silicate muscovite mica \cite{Chris11} or simplified model of two-dimensional hexagonal cristal lattice \cite{Bajars15,Bajars21}. In such case, the on-site potential $U$ models forces from atoms above and below the K-K sheet. Thus, the on-site potential is modeled as smoothed periodic function \cite{Chris11}:
\begin{equation}\label{eq:OnSite}
U(q_n) = 1-\cos(2 \pi q_n).
\end{equation}
Note that more realistic potentials for muscovite mica can also be considered \cite{Juan19}.

There are multiple options for modeling the interaction potential $V$ \cite{Allen89}. In this work we consider well known scaled Lennard-Jones potential in the following form:
\begin{equation}\label{eq:LJ}
V(r_n) = \epsilon \left( \left( \frac{1}{r_n} \right)^{12} - 2 \left(
    \frac{1}{r_n} \right)^{6} \right), \quad r_n = |q_{n+1}-q_n|,
\end{equation}
where the parameter $\epsilon>0$ describes the ratio of the atomic well depths associated to the interaction and the on-site potential. That is, large values of $\epsilon$ lead to a model describing one-dimensional Lennard-Jones fluid, while small values of $\epsilon$ lead to system of decoupled nonlinear oscillators. Further simplifications in \eqref{eq:Hamilt} have been made to only consider close neighbor interactions with periodic boundary conditions, i.e., ${q}_{N+1} = {q}_{1}$.

From the Hamiltonian \eqref{eq:Hamilt} obtained Hamiltonian equations:
\begin{align}
\dot{q}_n & = p_n, \label{eq:q}\\
\dot{p}_n & = -U'(q_n) + V'(|q_{n+1}-q_{n}|) - V'(|q_{n}-q_{n-1}|), \label{eq:p}
\end{align}
where $p_n$ is momentum, are solved numerically with the second order time reversible symplectic Verlet method \cite{Hairer}. In the following, without loss of generality, all analysis and numerical simulations are performed with time step $\tau=0.01$ and $\epsilon=0.05$.

To excite DB solutions we consider particles in their dynamical equilibrium states, i.e., $q^0_n=n-1$ and $p^0_n=0$, while exciting four neighboring particle momenta with the pattern
\begin{equation}\label{eq:patern}
\B{p}_0 = \gamma \, (-1, 2, -2, 1)^T, \quad \gamma > 0.
\end{equation}
Such excitation, due to nonlinearity, generates in our model stationary discrete breather. In the current study we only consider numerical data from stationary DB solutions while the extension to differentiate between stationary and moving breather solutions is left for future work.

In general, larger values of $\gamma$ generate localized waves with larger particle displacements. While such excitation \eqref{eq:patern} does not generate exact wave solutions the amount of generated linear phonon waves has been observed to be minimal. Not all values of $\gamma$ can be considered. Too small values will generate phonons, while too large values will force particles to move out from their potential wells or in addition to the stationary breather will also generate moving breathers. In addition, appropriate values of $\gamma$ depend on $\epsilon$. With spectral analysis of Section \ref{sec:Spectrum} and $\epsilon=0.05$ we find $\gamma$ lower bound $0.15$ for the initial momenta pattern \eqref{eq:patern} to generate stationary DB. On the other hand, to generate arbitrary phonon waves we can choose small random momenta values $p_n$ as initial conditions, as confirmed by the spectral analysis in Section \ref{sec:Spectrum}.

Both described initial conditions will be used to generate datasets in Section \ref{sec:Classification} to obtain classification and detection algorithms which are capable of distinguishing between localized and phonon waves. To see discrepancies between both types of waves in Figure \ref{fig:sol} we illustrate particle displacement function $u_n = q_n - q_n^0$ from two numerical simulations with $N=24$ at final computational time $T_{end}=5$. In Figure \ref{fig:Bsol} we show displacement function of a localized discrete breather solution computed with $\gamma=0.3$, while in Figure \ref{fig:Psol} we demonstrate particle displacements for phonon waves computed with random initial momentum values. It is easy to see larger displacement values for the breather solution and that it is spatially localized compared to the approximately linear waves. In addition, in Figure \ref{fig:Bsol} we have added markers to eight particle displacement values on which the breather is localized. Exactly eight data point values will be used in the classification algorithms of Section \ref{sec:Classification}.

\begin{figure}[t]
\centering
\subfigure[]{\label{fig:Bsol}
\includegraphics[trim=1cm 0.1cm 1cm 1cm,clip=true,width=0.48\textwidth]{Fig1a.jpg}}
\subfigure[]{\label{fig:Psol}
\includegraphics[trim=1cm 0.1cm 1cm 1cm,clip=true,width=0.48\textwidth]{Fig1b.jpg}}
\caption{Particle displacements from the numerical simulations at $T_{end}=5$ with $\epsilon=0.05$ and $N=24$. (a) stationary discrete breather solution initiated with $\gamma=0.3$ and markers indicating eight particles on which the breather is localized. (b) phonon wave solution which does not exhibit localization.}\label{fig:sol}
\end{figure}

\section{Spectral properties of lattice waves} \label{sec:Spectrum}

It is well known that linear phonon waves in crystal lattice models can be characterized by their dispersion relation \cite{Flach08}. For the model \eqref{eq:q}--\eqref{eq:p} the dispersion relation is
\begin{equation}\label{eq:dispersion}
\omega^2 = \omega_0^2 + 4 c^2 \sin^2 \left( \frac{k}{2} \right),
\end{equation}
where $\omega$ and $k$ are the phonon frequency and wave number, respectively, $\omega_0^2=\frac{\partial^2 U}{\partial q_n^2}(q_n^0)=4\pi^2$ is the oscillation frequency of the isolated oscillators when $\epsilon=0$ and $ c^2 = \frac{\partial^2 V}{\partial q_n^2}(q_n^0)= 72 \epsilon $ is the speed of sound.

In Figure \ref{fig:SP} we plot normalized two-dimensional spectrum of numerical simulation data in space and time, i.e., contours of 2D discrete Fourier transform. We consider a lattice of $N=128$ particles and integrate in time until $T_{end}=400$. The solid line indicates the dispersion relation \eqref{eq:dispersion}.

Excited stationary breather solution by \eqref{eq:patern} are optical breathers with particle frequencies above the phonon linear spectrum as can be seen in Figure \ref{fig:SPb}. Only small amount of energy is localized on the dispersion curve indicating minimal presence of phonons in numerical simulation.

In Figure \ref{fig:SPp} we considered random initial conditions for the momentum variable $p_n$ drawn from the normal distribution with zero mean and variance $\sigma^2=0.001$. Such initial condition leads to full spectrum of phonons which coincides with the dispersion curve \eqref{eq:dispersion}. Thus, this provides representative data of phonons for data analysis.

\begin{figure}[t]
\centering
\subfigure[]{\label{fig:SPb}
\includegraphics[trim=0.5cm 0.5cm 0.8cm 1cm,clip=true,width=0.48\textwidth]{Fig2a.jpg}}
\subfigure[]{\label{fig:SPp}
\includegraphics[trim=0.5cm 0.5cm 0.8cm 1cm,clip=true,width=0.48\textwidth]{Fig2b.jpg}}
\caption{Contours of normalized two-dimensional discrete Fourier transform in space and time of numerical simulation with $N=128$ and $T_{end}=400$. Dispersion relation \eqref{eq:dispersion} is indicated by the solid line. (a) stationary discrete breather spectrum computed with the initial pattern \eqref{eq:patern} and $\gamma=0.3$. (b) phonon wave spectrum computed with random initial condition and coinciding with the dispersion curve \eqref{eq:dispersion}.}\label{fig:SP}
\end{figure}

Spectral analysis of this section not only provides insights into the lattice wave properties but also allows to justify the methods for obtaining a so called labeled training and testing data for classification problem presented in the following section.

\section{Classification of lattice waves} \label{sec:Classification}

In this section we describe classification algorithms of crystal lattice waves based on locally sampled data. Machine learning techniques \cite{ML} are considered for data dimensionality reduction and building an efficient classifier. All calculations are performed in {\it Python} with open source {\it scikit-learn}\footnote{\url{https://scikit-learn.org/stable/}} library.

We differentiate localized discrete breather solutions, see Figure \ref{fig:Bsol}, from nonlocalized linear phonons, see Figure \ref{fig:Psol}. Training (including validation) and testing data is obtained from total of $N_{sim}=1000$ numerical simulations with $0.75$ localized and nonlocalized wave data proportion. Local data of particle displacements, momenta and energies are sampled from eight neighboring particles indicated by red dots in Figure \ref{fig:Bsol}. We consider three different datasets for classification: $X_{u}$, $X_{u,p}$ and $X_{u,p,E}$, where the dataset $X_{u}$ only contains particle displacement $u_n$ values, the dataset $X_{u,p}$ contains particle displacement $u_n$ and momenta $p_n$ values, while the dataset $X_{u,p,E}$ in addition to displacement and momenta values contains also particle energy density data defined by
\begin{equation}\label{eq:En}
E_n := \frac{1}{2} {p}_{n}^2 + U(q_{n}) + \frac{1}{2}\left( V(|q_{n}-q_{n-1}|) + V(|q_{n+1}-q_{n}|) \right).
\end{equation}

We excite discrete breather solutions with momenta pattern \eqref{eq:patern} where parameter $\gamma$ values are randomly chosen from the uniform distribution $\mathcal{U}(0.15,1.2)$. Linear phonon waves are generated with random initial particle displacement and momenta values from the uniform distribution $\mathcal{U}(-0.01,0.01)$.

\subsection{Dimensionality reduction}

To increase the efficiency of the classification algorithms we consider two dimensionality reduction techniques, i.e., {\it Principal Component Analysis} (PCA) and {\it Locally Linear Embedding} (LLE) \cite{ML}. PCA is a data projection algorithm onto so called principal components. In our experiments we find that with projection down to only two principal components we are able to preserve more than $95\%$ of datasets's variance. On the other hand, LLE is a dimensionality reduction technique that computes low dimensional, neighborhood preserving embeddings of high dimensional data. Performing LLE we considered twenty neighbors for each data point.

In Figure \ref{fig:dimred} we illustrate dimension reduction of all three datasets $X_{u}$, $X_{u,p}$ and $X_{u,p,E}$ applying both dimensionality reduction techniques PCA and LLE, respectively. Colorbars in Figure \ref{fig:dimred} indicate $\gamma$ values of the pattern \eqref{eq:patern} used to obtain particular DB solution data. Apart from the PCA projection of the dataset $X_{u,p}$, see Figure \ref{fig:PCAWP}, the $\gamma$ values can be traced along the continuous curve in two dimensions. Interestingly, in Figure \ref{fig:LLEX} the curve is a straight line parallel to $y$-axis.

In Figures \ref{fig:PCAW}, \ref{fig:PCAWP} and \ref{fig:PCAX} we projected datasets down onto two principal components $(PC1,PC2)$ while in Figures \ref{fig:LLEW}, \ref{fig:LLEWP} and \ref{fig:LLEX} we show local data linear embedding in two-dimensional $(LLE1,LLE2)$ coordinates. Notice that linear wave data is hardly separated from nonlinear localized wave data in Figures \ref{fig:PCAW}, \ref{fig:LLEW} and \ref{fig:PCAWP}. For the dataset $X_{u,p}$ Figure \ref{fig:LLEWP} suggests that LLE is superior over PCA for the classification problem. Locally linear embedding of the dataset $X_{u,p,E}$ is able to produce the largest separation between two types of wave data. Thus, information of particle displacements and momenta alone may not be sufficient to build an optimal classifier for detecting ILM, e.g., in numerical simulations, see Section \ref{sec:Applications}.

\begin{figure}[!ht]
\centering
\subfigure[]{\label{fig:PCAW}
\includegraphics[width=0.47\textwidth]{Fig3a.jpg}}
\subfigure[]{\label{fig:LLEW}
\includegraphics[width=0.47\textwidth]{Fig3b.jpg}}
\subfigure[]{\label{fig:PCAWP}
\includegraphics[width=0.47\textwidth]{Fig3c.jpg}}
\subfigure[]{\label{fig:LLEWP}
\includegraphics[width=0.47\textwidth]{Fig3d.jpg}}
\subfigure[]{\label{fig:PCAX}
\includegraphics[width=0.47\textwidth]{Fig3e.jpg}}
\subfigure[]{\label{fig:LLEX}
\includegraphics[width=0.47\textwidth]{Fig3f.jpg}}
\caption{Dimensionality reduction with PCA and LLE of three datasets $X_{u}$, $X_{u,p}$ and $X_{u,p,E}$ and decision boundaries of SVC. (a) PCA of $X_{u}$ with nonlinear SVC. (b) LLE of $X_{u}$ with nonlinear SVC. (c) PCA of $X_{u,p}$ with nonlinear SVC. (d) LLE of $X_{u,p}$ with nonlinear SVC. (e) PCA of $X_{u,p,E}$ with linear SVC. (f) LLE of $X_{u,p,E}$ with linear SVC.}\label{fig:dimred}
\end{figure}

\subsection{Classification algorithms}

Linear and nonlinear kernel {\it Support Vector Machine} \cite{ML} classifiers (SVC) are trained for lattice wave data classification. SVC are supervised learning algorithms that construct a decision boundary with the largest distance to the nearest training-data point of each data class. The Gaussian {\it Radial Basis Function} (RBF) kernel is used in nonlinear SVC which showed improved performance over other nonlinear kernels, such as polynomial and hyperbolic tangent. Soft margin SVC are considered where optimal hyperparameter values are found with {\it grid search} technique.

To evaluate the performance of linear and nonlinear classifiers for each dataset $\mathrm{K}$-fold cross-validation was considered with $0.7$ training and testing data proportion. In our experiments we considered $\mathrm{K}=5$. The {\it precision} and {\it recall} scores were recorded on cross-validation and testing datasets. The precision is a metric of accuracy of the positive (localized wave) predictions while the recall is a metric of positive instances that are correctly detected by the classifier.

Averaged precision and recall values over $100$ random {\it train-test-splits} of datasets $X_{u}$, $X_{u,p}$ and $X_{u,p,E}$ are listed in Table \ref{tab:Precision_recall}. For each train-test-split we compute precision and recall scores for the testing data and validation dataset obtained from the $\mathrm{K}$-fold of training data. Notice that in most cases recall values are smaller than precision values which is attributed to the fact that more localized wave data is classified as linear opposed to linear waves being classified as localized. In addition to mean values we also computed (not shown) variance values which in all cases, if not exactly zero, did not exceed $10^{-2}$ value.

Examining Table \ref{tab:Precision_recall} we observe exceptionally good results except for the PCA of the dataset $X_{u,p}$ when linear SVC is used. That is also evident in Figure \ref{fig:PCAWP}, where it is easy to see that nonlinear waves can not be separated from linear waves by a straight line. In that case the use of nonlinear SVC is required. Very good results are demonstrated by PCA and LLE dimensionality reduction of the dataset $X_{u,p,E}$, either with linear or nonlinear SVC. In this case for simplicity we advocate the use of linear SVC. Precision and recall values equal to one state that none of the data samples were incorrectly classified by SVC. As follows we consider linear SVC for the dataset $X_{u,p,E}$ and nonlinear classifier for the datasets $X_{u}$ and $X_{u,p}$, which is illustrated in Figure \ref{fig:dimred} where the prediction boundaries and contours of the decision function are drawn.

In summary, Table \ref{tab:Precision_recall} shows that the classifier is not sensitive to random train-test-split of the datasets and that for the dataset $X_{u,p,E}$ SVC is essentially a perfect classifier for both dimensionality reduction techniques. To further differentiate the classifiers we consider their classification performance on imperfect data, see the following section.

\begin{table}[t]
	\begin{center}
		\begin{tabular}{|c||c|c||c|c||c|}
			\hline
			\multirow{2}{*}{} & \multicolumn{2}{c||}{{Linear kernel}} &
			\multicolumn{2}{c||}{{RBF kernel}} & \\
			\cline{2-5}
			& {PCA} & {LLE} & {PCA} & {LLE} & \\
			\specialrule{.1em}{.05em}{.05em}
			{Precision} (validation set) & 0.9985 & 0.9623 & 1 & 0.9953
			& \multirow{4}{*}{$X_u$} \\
			{Recall} (validation set) & 0.9533 & 0.9478 & 0.9998 & 0.9354 & \\
			\cline{1-5}
			{Precision} (testing set) & 0.9984 & 0.9658 & 1 & 0.9965 & \\
			{Recall} (testing set) & 0.9516 & 0.9487 & 0.9992 & 0.9370 & \\
			\specialrule{.1em}{.05em}{.05em}
			{Precision} (validation set) & 0.7947 & 0.9927 & 1 & 1
			& \multirow{4}{*}{$X_{u,p}$} \\
			{Recall} (validation set) &
0.9544 & 0.9993 & 0.9754 & 0.9990 & \\
			\cline{1-5}
			{Precision} (testing set) & 0.7978 & 0.9928 & 1 & 0.9999 & \\
			{Recall} (testing set) &
0.9519 & 0.9999 & 0.9758 & 0.9991 & \\

			\specialrule{.1em}{.05em}{.05em}
			{Precision} (validation set) & 1 & 1 & 1 & 1 & \multirow{4}{*}{$X_{u,p,E}$} \\
			{Recall} (validation set) & 1 & 1 & 1 & 1 & \\
			\cline{1-5}
			{Precision} (testing set) & 1 & 1 & 1 & 1 & \\
			{Recall} (testing set) & 1 & 1 & 1 & 1 & \\
			\hline
		\end{tabular}
	\end{center}
	\caption{Precision and recall metrics averaged over $100$ random train-test-splits of different classification algorithms, dimensionality reduction techniques and datasets.}
	\label{tab:Precision_recall}
\end{table}

\subsection{Classification of imperfect data}

So far we have considered datasets of perfectly sampled data of linear and localized nonlinear waves. With application in mind, i.e., detecting ILM in numerical simulations, we consider additional datasets $X_{u,p,E}^{n_{off}}$ of localized discrete breather data from $N_{sim}=1000$ simulations with randomly chosen $\gamma$ values from the uniform distribution $\mathcal{U}(0.15,1.2)$. We define index $n_{off}$, where $n_{off}=0,1,\dots,8$, which indicates the shift in particles from which the data is sampled, i.e., the data is sampled from particles $q_{m+n_{off}}$, where particle $q_m$ is the first particle in Figure \ref{fig:Bsol} identified with red dot. Thus, the sampled datasets $X_{u,p,E}^{n_{off}}$ consist of partially linear and nonlinear wave data.

\begin{figure}[t]
\centering
\includegraphics[width=0.6\textwidth]{Fig4.jpg}
\caption{Percentages of classified imperfect localized wave data of the datasets $X_{u,p,E}^{n_{off}}$ by linear SVC with dimensionality reduction methods PCA and LLE.}
\label{fig:dataoff}
\end{figure}

We classify samples from the datasets $X_{u,p,E}^{n_{off}}$ using linear SVC and both dimensionality reduction methods PCA and LLE. In Figure \ref{fig:dataoff} we plot the percentage of wave data being classified as localized wave or linear wave depending on the shift index $n_{off}$. When $n_{off}=0$, the data consists of purely nonlinear wave data and we can see $100\%$ identification of it. On the other hand, when $n_{off}=8$ the dataset does not contain information of nonlinear waves and we see $100\%$ identification of linear wave data. Both classifiers demonstrate gradual decrease in localized wave identification as $n_{off}$ value increases. Notice that the curve of LLE drops faster to zero compared to PCA. This has an impact on detection region size, i.e., SVC with PCA produces wider detection regions of localization as indicated by numerical simulations of Section \ref{sec:Applications}.

\section{Applications}\label{sec:Applications}

In this section we apply classification algorithms developed above to detect localized DB of the one-dimensional crystal lattice model \eqref{eq:Hamilt}. We consider sliding window object detection method by sliding a window of eight particle size at each time step over the computational data. To detect region of a localized wave we define a discrete density function $\rho_n$, $n=1,\dots,N$, which takes a positive value if data at particular eight particles is classified as being of localized wave, otherwise the density value is equal to zero. In the following Figures \ref{fig:ApplA} and \ref{fig:ApplBC} we identify all particles by red dots where the localization is detected.

For the first example we consider a lattice of $N=64$ particles and excite two stationary breathers using momenta pattern \eqref{eq:patern} with $\gamma=0.55$ and $\gamma=0.35$ at particles with indices $(16,17,18,19)$ and $(44,45,46,47)$, respectively. We solve the Hamiltonian system \eqref{eq:q}--\eqref{eq:p} until time $T_{end}=1$ and detect the regions of localization, see Figure \ref{fig:ApplA}.

\begin{figure}[t]
\centering
\subfigure[]{\label{fig:ApplA_1}
\includegraphics[trim=1.1cm 0.5cm 2cm 3.5cm,clip=true,width=0.47\textwidth]{Fig5a.jpg}}
\subfigure[]{\label{fig:ApplA_2}
\includegraphics[trim=1.1cm 0.5cm 2cm 3.5cm,clip=true,width=0.47\textwidth]{Fig5b.jpg}}
\subfigure[]{\label{fig:ApplA_3}
\includegraphics[trim=1.1cm 0.5cm 2cm 3.5cm,clip=true,width=0.47\textwidth]{Fig5c.jpg}}
\subfigure[]{\label{fig:ApplA_4}
\includegraphics[trim=1.1cm 0.5cm 2cm 3.5cm,clip=true,width=0.47\textwidth]{Fig5d.jpg}}
\caption{Detection of numerically simulated two stationary breathers excited by momenta pattern \eqref{eq:patern} with $\gamma=0.55$ and $\gamma=0.35$ on the left and right, respectively. Detected localization region particles are identified by red dots. (a) linear SVC with PCA trained on the dataset $X_{u,p,E}$. (b) linear SVC with LLE trained on the dataset $X_{u,p,E}$. (c) nonlinear SVC with PCA trained on the dataset $X_{u}$. (d) nonlinear SVC with LLE trained on the dataset $X_{u}$.}\label{fig:ApplA}
\end{figure}

For this experiment we used four classifiers. In Figure \ref{fig:ApplA_1} we demonstrate results of linear SVC classifier with PCA, while in Figure \ref{fig:ApplA_2} we consider linear SVC classifier with LLE. Both classifiers are trained on the dataset $X_{u,p,E}$ and are able to capture regions of localization very well. As previously suggested, see Figure \ref{fig:dataoff}, on average SVC classifier with PCA produces wider detection regions, as can be seen from number of colored particles.

In comparison, nonlinear classifiers trained on the dataset $X_{u}$ containing only the information of particle displacements, see Figures \ref{fig:ApplA_3} and \ref{fig:ApplA_4}, are not always able to capture small displacement ILM. Similar issues (not shown) were observed applying nonlinear classifiers trained on the dataset $X_{u,p}$. Despite these shortcomings the classifiers can still be applied if additional information of the energy \eqref{eq:En} is unknown.

To further demonstrate methods' robustness and future research directions we consider two additional examples. In the first example we consider a stationary discrete breather in noisy background by initially randomly exciting all particles and their momenta with addition to the pattern \eqref{eq:patern} with $\gamma=0.5$. In this example we performed calculations until time $T_{end}=10$ and used linear SVC classifier with LLE trained on the dataset $X_{u,p,E}$, see Figure \ref{fig:ApplC}. Figure \ref{fig:ApplC} illustrates very well that the method was able to detect the regions of localization despite the presence of large amount of phonon waves. These results suggest that the classification algorithms can be applied to thermostatic dynamics, i.e., during experiments at fixed temperature.

\begin{figure}[t]
\centering
\subfigure[]{\label{fig:ApplC}
\includegraphics[trim=1.1cm 0.5cm 2cm 3.5cm,clip=true,width=0.47\textwidth]{Fig6a.jpg}}
\subfigure[]{\label{fig:ApplB}
\includegraphics[trim=1.1cm 0.5cm 2cm 3.5cm,clip=true,width=0.47\textwidth]{Fig6b.jpg}}
\caption{Detection of localization regions in numerical simulations with linear SVC and LLE trained on the dataset $X_{u,p,E}$. Particles of detected regions are identified by red dots. (a) simulation of stationary discrete breather ($\gamma=0.5$) in noisy background. (b) simulation of two mobile breather collision initiated by momenta pattern \eqref{eq:pattern2} with $\gamma=0.45$ (left) and $\gamma=-0.65$ (right).}\label{fig:ApplBC}
\end{figure}

In the second example we consider a simulation of two mobile breather collision initiated by three neighboring particle momenta with pattern
\begin{equation}\label{eq:pattern2}
\B{p}_0 = \gamma \, (-1, 2, -1)^T.
\end{equation}
The mobile breather propagating to right was excited with $\gamma=0.45$ while the mobile breather propagating to left was excited with $\gamma=-0.65$. Integration in time was performed until $T_{end}=100$. Around time $t=50$ both breathers collide and pass through or reflect from each other and continue to propagate, see Figure \ref{fig:ApplB}. Although the classifier was trained on data obtained from stationary discrete breathers, we are still able to detect the localization regions before, during and after the collision. As in Figure \ref{fig:ApplC} we used linear SVC with LLE trained on the dataset $X_{u,p,E}$. These satisfactory results show that the methodology can be used beyond just detection of stationary DB.

\section{Discussion and conclusions}\label{sec:Conclusions}

The results of Section \ref{sec:Applications} demonstrate that classification algorithms, such as SVC, either linear or nonlinear, can be used to detect localization regions in numerical simulations of ILM. The strength of the developed approach lies in the fact that localization regions can be detected from locally sampled data. In future research we will look for the method which makes predictions with the most optimal width of the localization region. For example, in this work we only considered data sampled from eight neighboring particles, other options could also be investigated.

Numerical examples demonstrated that our classifiers are also able to detect mobile DB but lack the ability to differentiate between these two types of localized solutions. To address this we will explore multilabel and multioutput classification algorithms in future work. In addition, we plan to extend our work to higher dimensional crystal lattice models and include kink solutions.

For efficiency and visualization purposes with high success we applied two dimensionality reduction methods before training SVC. The best results were obtained with LLE on the dataset $X_{u,p,E}$, such that we could use linear SVC and find more optimal localization regions in numerical experiments. In future work we plan to explore, in particular, more other manifold learning techniques.

To demonstrate versatility of our approach we considered three different datasets for training. We found that the best results are obtained if in addition to particle displacement values particle momenta and energy density values are also added to the dataset. Nevertheless, classifiers can be built on all three datasets and overall good performance can be achieved.

\section*{Acknowledgements}
J.~Baj\={a}rs acknowledges support from PostDocLatvia grant No.1.1.1.2/VIAA/4/20/617.

\bibliographystyle{unsrt}