\section{Introduction}

If an ink drop falls into a pond, the stain diffuses away until no trace of
it remains whatsoever. One may naturally say that such a process is in fact
irreversible. In the microscopic world, similar phenomena are also
ubiquitous. For instance, let us consider a many-spin quantum system in
thermal equilibrium where a local polarization excess is injected. Then,
this excitation would spread all over as consequence of spin-spin
interactions. Such an apparently irreversible process is known as \textit{%
spin diffusion} \cite{BlumeHubbard1970,forsterbook} and it can lead the
system back to equilibrium. However, this naive picture has its limitations.
On the one hand, spreading is not always the rule as there are physical
situations where the initial excitation does not vanish. This is the case of
Anderson localization \cite{altshuler2006,Zangara2015PIP} or when the
excitation remains topologically protected \cite{Wen_book}. On the other
hand, even in the cases where the system seems to have reached an
equilibrium state, the unitarity of quantum dynamics ensures a precise
memory of the non-equilibrium initial condition. Then, if some experimental
protocol could reverse the many-body dynamics, it would drive the system
back to the initial non-equilibrium state \cite{Hahn_atomicMemory}. Such a
general idea defines the Loschmidt echo (LE), which embodies the various
time-reversal procedures implemented in nuclear magnetic resonance (NMR)
\cite{prosen,Jacquod, scholarpedia}.

The first NMR time-reversal experiment was performed by Hahn in the 1950's
\cite{Hahn1950}. The procedure, known as spin echo, reverses the precession
dynamics of each independent spin around its local magnetic field by
inverting the sign of the Zeeman energy. However, the\ sign of the energy
associated to the spin-spin interactions is not inverted and, accordingly,
the echo signal is degraded. Such a decay occurs within the time scale $%
T_{2} $ that characterizes the spin-spin interactions. Indeed, these
interactions determine the survival of a spin excitation at short times as $%
\sim 1-\left( t/T_{2}\right) ^{2}$ and its later complex dynamics generating
a diffusive spreading. By the early 1970's, Kessemeier, Rhim, Pines and
Waugh implemented the reversal of the dynamics induced by the spin-spin
dipolar interaction \cite{Kessemeier1971,Rhim1971}. This results in the
\textquotedblleft magic echo\textquotedblright\, which indicates the
recovery of a global polarization state. Two decades later, Ernst and
collaborators introduced the \textquotedblleft polarization
echo\textquotedblright\ \cite{Ernst1992}. There, a local excitation injected
in a many-spin system is let to evolve, then time-reversed and finally
detected locally at the initial spot. While the success of these time
reversal echoes unambiguously evidenced the deterministic nature of
spin-dynamics in NMR, it is clear that the reversal is unavoidably degraded
by uncontrolled internal or environmental degrees of freedom or by
imperfections in the pulse sequences. Furthermore, the degradation seems to
occur in a time scale, say $T_{3}$, much shorter than a naive estimation of
the characteristic scale of these perturbations, say $\tau _{\Sigma }$.
Then, the question that arises is whether the complexity inherent to a large
number of correlated spins would enhance the fragility of the procedure
under perturbations.

A next generation of experiments in organic crystals \cite%
{patricia98,usaj-physicaA,MolPhys} seemed to confirm that the experimental $%
T_{3}$ never exceeds a few times $T_{2}$. In other words, $T_{3}$ keeps tied
to the time scale that characterizes the reversible many-body interaction.
This led to postulate that in an infinite many-spin system the complex
dynamics could favor the action of any small non-inverted interaction that
perturbs the reversal procedure. Thus, reversible interactions become
determinant for the irreversibility rate. This constitutes our \textit{%
Central Hypothesis of Irreversibility}. Such a wisdom is further reinforced
by the natural association of many-body complexity with a form of chaos \cite%
{Flambaum2000,Flambaum2001d} and\ the confirmation that quantum dynamics of
classically chaotic systems should manifest a dynamical instability \cite%
{peres1984} which leads to an \textit{environment-independent decoherence
rate}\cite{jalpa,Cookprb2004}.

During the last decade, solid-state NMR kept providing a versatile testing
bench to study time reversal in large spin arrays \cite%
{Suter2004,Boutis2012,Cappellaro2013,Claudia2014,Alvarez2014experimental}.
As a matter of fact, a standard experiment involves a crystalline sample
with an infinitely large number of spins. In contrast, the numerical test of
many spin dynamics has to be restricted to strictly finite systems \cite%
{Zangara2012,Fine2014}. While this appears to be a major limitation, it
allows the analysis of a situation that the experiments cannot achieve:
moving progressively from small systems to larger ones with a fully
controlled perturbation. The expectancy is that a sort of finite size
scaling may allow to identify an emergent mechanism that rules reversibility
in the thermodynamic limit \cite{Zangara2015PRA}. As in the experiments, the
witness for such a transition should be the LE as measured by a single spin
autocorrelation function $M_{1,1}$, i.e. a local polarization. For short we
call $M_{1,1}$\ the \textit{local LE}. It is not difficult to probe that $%
\Pi _{1,1}\equiv (M_{1,1}+1)/2$ is the probability that a given spin, say
the $1^{\mathrm{st}}$, remains up after the whole procedure. Besides, in a
case of a spin excitation in a 1D chain with $XY$ interactions \cite%
{mesoECO-PRL1995,mesoECO-exp,ERNST-CPL-1997} $M_{1,1}$ precisely coincides
\cite{Zangara2012} with the global overlap of two one-body wave functions as
defined for semiclassical systems \cite%
{jalpa,Jacquod-Beenakker-2001,Jacquod2002}. The square of the overlap
between the initial and final many-body wave functions, $M_{MB}$, defines a
\textit{global} or \textit{many-body LE}. It is important to notice that $%
M_{MB}$ has not been addressed experimentally, but nevertheless it is a
natural magnitude in theory \cite{prosen2007_rev,Santos2014a,Santos2014b}.
Thus, we are left without a precise relation between the object of
theoretical studies and experimental ones, i.e. $M_{MB}$\ and $M_{1,1}$\
respectively. This missing link is the central question we address in this
paper.

Here, we consider a system of $N$ spins whose initial state is given by a
local excitation injected in an high temperature state. Firstly, we discuss
the formal relation between $M_{MB}$\ and $M_{1,1}$, which is derived
exactly at least for very short times. In particular, we assess how the $N$%
-dependence or extensivity of $M_{MB}$ is evidenced in the time scales
involved. This leads us to hint that the revival of a many-spin state
results from the recovery of each individual spin configuration, much as if
they were statistically independent events. Since in the initial high
temperature state there are $N/2$ spins up, their rough statistical
independence would lead to a behavior of the sort of $M_{MB}\sim \left( \Pi
_{1,1}\right) ^{N/2}\sim \left( M_{1,1}\right) ^{N/4}$. This is confirmed by
the numerical evaluation of the LE in a specific spin model.

Furthermore, we discuss a general picture beyond the short-time regime,
where the decay of $M_{1,1}$ results from the interplay between the time
scale that characterizes the reversible interactions ($T_{2}$) and that of
the perturbation ($\tau _{\Sigma }$). This would ultimately lead to the
experimentally observed $T_{3}$. In such a sense, our analysis provides a
conceptual hinge between the theoretical and the experimental realms.

The paper is organized as follows. In Sec. \ref{Sec_formulation} we
introduce the LE framework:\ the initial state and the time-reversal
procedure. Here, we define both the local and global LE. In Sec. \ref%
{Sec_locglob} we compute the short time expansions for the local LE and its
non-local contributions (in particular, the many-body LE). This allows us to
discuss a general picture of the LE decay in terms of the times scales that
characterize $M_{MB}$\ and $M_{1,1}$. The dependence with $N$ is discussed
in terms of the extensivity of $M_{MB} $ and statistical independence of the
local autocorrelations. In Sec. \ref{Sec_1Dmodel} we assess our expectancies
by a numerical evaluation of the LE in a spin system. Section \ref%
{Sec_conclu} summarizes our main conclusions and some of the important open
questions in the field.

\section{The Loschmidt echo in spin systems.\label{Sec_formulation}}

Let us first specify the initial condition of a \textquotedblleft local
excitation in many-spin system\textquotedblright . We consider $N$ spins $%
1/2 $ in an infinite temperature state, i.e. completely depolarized mixture,
plus a locally injected polarization,

\begin{equation}
\hat{\rho}_{0}=\frac{1}{2^{N}}(\mathbf{\hat{I}}+2\hat{S}_{1}^{z}).
\label{rho_inicial}
\end{equation}%
Here, the spin $1$ is polarized while the others are not, i.e. $tr[\hat{S}%
_{i}^{z}\hat{\rho}_{0}]=\frac{1}{2}\delta _{i,1}$. Such an initial state can
be experimentally implemented not only in NMR \cite{Cappellaro2014_review}
but also in cold atoms \cite{Bloch2013}.

As in the early LE experiments \cite{patricia98,usaj-physicaA,MolPhys}, our
numerical evaluation focuses on an imperfect time reversed evolution of the
excitation, followed by a local measurement. The procedure is depicted in
Fig. \ref{Fig_LE_standard}. A many-spin Hamiltonian $\hat{H}_{0}$ rules the
\textit{forward} evolution of the system up to a certain time $t_{R}$. At
that moment, an inversion of the sign of $\hat{H}_{0}$ is performed, leading
to a symmetric \textit{backward }evolution. Nevertheless, there are
unavoidable perturbations, denoted by $\hat{\Sigma}$, that could arise from
the incomplete control of the Hamiltonian, acting on both periods. Thus,
evolution operators for these $t_{R}$-periods are $\hat{U}%
_{+}^{{}}(t_{R})=\exp [-\frac{\mathrm{i}}{\hbar }(\hat{H}_{0}+\hat{\Sigma}%
)t_{R}]$ and $\hat{U}_{-}^{{}}(t_{R})=\exp [-\frac{\mathrm{i}}{\hbar }(-\hat{%
H}_{0}+\hat{\Sigma})t_{R}]$ respectively. It is quite practical to define
the LE\ operator as:

\begin{equation}
\hat{U}_{LE}^{{}}(2t_{R})=\hat{U}_{-}^{{}}(t_{R})\hat{U}_{+}^{{}}(t_{R}),
\label{Ule}
\end{equation}%
which produces an imperfect refocusing at time $2t_{R}$. A local measurement
of the polarization, performed at site $1$, defines the local LE:

\begin{equation}
M_{1,1}(t)=2tr[\hat{S}_{1}^{z}\hat{U}_{LE}^{{}}(t)\hat{\rho}_{0}\hat{U}%
_{LE}^{\dag }(t)]=2tr[\hat{S}_{1}^{z}\hat{\rho}_{t}].
\label{autocorrelacion1}
\end{equation}%
Here, we choose as free variable $t=2t_{R}$, the total elapsed time in the
presence of the perturbation. The time dependence of $\hat{\rho}_{t}$ in the
Schr\"{o}dinger picture is,

\begin{equation}
\hat{\rho}_{t}=\hat{U}_{LE}^{{}}(t)\hat{\rho}_{0}\hat{U}_{LE}^{\dag }(t).
\label{autocorrelacion11}
\end{equation}

\begin{figure}[!h]
\centering\includegraphics[width=5.0in]{fig1.eps}
\caption{The scheme of the local LE, $M_{1,1}(t)$. (a) The initial state is
given by a local excitation in a high temperature spin system, as stated in
equation (\protect\ref{rho_inicial}). The system is let to evolve under the
Hamiltonian $\hat{H}_{0}+\hat{\Sigma}$ and the excitation diffuses all over
until a time $t=t_{R}$ (b). At that time, a time reversal is performed,
leading to a backward evolution ruled by $-\hat{H}_{0}+\hat{\Sigma}$. At
time $t=2t_{R}$ (c), a local measurement is performed at the initial spot.}
\label{Fig_LE_standard}
\end{figure}

Using equation (\ref{rho_inicial}), and after some algebraic manipulation, the LE
can be explicitly written as a correlation function:

\begin{equation}
M_{1,1}(t)=\frac{1}{2^{N-2}}tr[\hat{U}_{LE}^{\dag }(t)\hat{S}_{1}^{z}(0)\hat{%
U}_{LE}^{{}}(t)\hat{S}_{1}^{z}(0)]=\frac{tr[\hat{S}_{1}^{z}(t)\hat{S}%
_{1}^{z}(0)]}{tr[\hat{S}_{1}^{z}(0)\hat{S}_{1}^{z}(0)]}.
\label{autocorrelacion2}
\end{equation}%
Here, the time dependence is written according to the Heisenberg picture,

\begin{equation}
\hat{S}_{1}^{z}(t)=\hat{U}_{LE}^{\dag }(t)\hat{S}_{1}^{z}(0)\hat{U}%
_{LE}^{{}}(t).  \label{autocorrelacion21}
\end{equation}%
Notice that equation (\ref{autocorrelacion2}) is an explicit correlation function
at the same site but different times, i.e. an \textit{auto}correlation. This
kind of a correlation has been recently employed to address localization
phenomena in spin systems \cite{Zangara2013PRB,knap2015,Zangara2015PIP} and
it generalizes the standard one employed to assess spin diffusion \cite%
{LevPasCalvo1989}. In terms of the Hilbert--Schmidt inner product between
the initial and the time evolved density matrices, i.e. equations (\ref%
{rho_inicial}) and (\ref{autocorrelacion11}) respectively, the LE can be
written as \cite{usaj-physicaA,PazZurek2003,Bendersky2013}:

\begin{equation}
M_{1,1}(t)=2^{N}tr[\hat{\rho}_{0}\hat{\rho}_{t}]-1=2\frac{tr[\hat{\rho}_{0}%
\hat{\rho}_{t}]}{tr[\hat{\rho}_{0}\hat{\rho}_{0}]}-1,
\label{autocorrelacion3}
\end{equation}%
which, in the present case, progressively decays from $1$ to $0$ as it
occurs with the statistical overlap between two wave packets in the standard
LE definition \cite{jalpa}.

Equivalent expressions for the LE can be derived decomposing the statistical
state in a simpler basis. In order to proceed with the pure state
decomposition of $\hat{\rho}_{0}$, we consider the computational Ising basis
$\left\{ \left\vert \beta _{i}\right\rangle \right\} $, also known as $S^{z}$%
-decoupled basis. Additionally, we define the set $\mathcal{A}$ of indexes $j
$ that label basis states which have the $1^{\mathrm{st}}$ spin pointing up,
i.e. $j\in \mathcal{A}$ $\Leftrightarrow \hat{S}_{1}^{z}\left\vert \beta
_{j}\right\rangle =+\frac{1}{2}\left\vert \beta _{j}\right\rangle $. It is
straightforward to verify that $\hat{\rho}_{0}=\sum_{j\in \mathcal{A}%
}2^{-(N-1)}\left\vert \beta _{j}\right\rangle \left\langle \beta
_{j}\right\vert $. Then, as introduced in Ref. \cite{mesoECO-PRL1995},

\begin{equation}
M_{1,1}(t)=2\left[ \sum_{i\in \mathcal{A}}\sum_{j\in \mathcal{A}}\frac{1}{%
2^{N-1}}\left\vert \left\langle \beta _{j}\right\vert \hat{U}%
_{LE}^{{}}(t)\left\vert \beta _{i}\right\rangle \right\vert ^{2}-\frac{1}{2}%
\right] =2\left[ \Pi _{1,1}(t)-\frac{1}{2}\right] .  \label{autoco2}
\end{equation}%
Here, $\Pi _{1,1}(t)$ denotes the \textit{probability} that the $1^{st}$
spin keeps pointing up after a time $t$. After some manipulation,

\begin{equation}
\begin{split}
M_{1,1}(t)& =2\left[ \sum_{i\in \mathcal{A}}\sum_{j\in \mathcal{A}}\frac{1}{%
2^{N-1}}\left\vert \left\langle \beta _{j}\right\vert \hat{U}%
_{LE}^{{}}(t)\left\vert \beta _{i}\right\rangle \right\vert ^{2}-\frac{1}{2}%
\right] \\
& =\left[ \sum_{i\in \mathcal{A}}\frac{1}{2^{N-1}}\left( \vphantom{\sum_{i%
\in \mathcal{A}}}\left\vert \left\langle \beta _{i}\right\vert \hat{U}%
_{LE}^{{}}(t)\left\vert \beta _{i}\right\rangle \right\vert ^{2}+\right.
\right. \\
& \ \left. \left. \qquad \qquad \qquad \quad +\sum_{j\in \mathcal{A}\text{ (}%
j\neq i)}\left\vert \left\langle \beta _{j}\right\vert \hat{U}%
_{LE}^{{}}(t)\left\vert \beta _{i}\right\rangle \right\vert ^{2}-\sum_{j\in
\mathcal{B}}\left\vert \left\langle \beta _{j}\right\vert \hat{U}%
_{LE}^{{}}(t)\left\vert \beta _{i}\right\rangle \right\vert ^{2}\right) %
\right] .
\end{split}
\label{contribuciones}
\end{equation}%
Here, $\mathcal{B}$ stands for the complement of $\mathcal{A}$, i.e. $j\in
\mathcal{B}$ $\Leftrightarrow \hat{S}_{1}^{z}\left\vert \beta
_{j}\right\rangle =-\frac{1}{2}\left\vert \beta _{j}\right\rangle $. One can
naturally identify and define the two terms that contribute to the local
polarization $M_{1,1}(t)$. The first sum in equation (\ref{contribuciones})
stands for the average probability of revival of the many-body states,
denoted by $M_{MB}(t)$,

\begin{equation}
M_{MB}(t)=\sum_{i\in \mathcal{A}}\frac{1}{2^{N-1}}\left\vert \left\langle
\beta _{i}\right\vert \hat{U}_{LE}^{{}}(t)\left\vert \beta _{i}\right\rangle
\right\vert ^{2}.  \label{autocorrelacion51}
\end{equation}%
The second sum in equation (\ref{contribuciones}) represents the average
probability of changing the configuration of any spin except the $1^{\mathrm{%
st}}.$ The third sum stands for the average probability that the $1^{\mathrm{%
st}}$ spin has actually flipped, i.e. of all those processes that do not
contribute to $M_{1,1}(t)$. Then, the processes that contribute to $%
M_{1,1}(t)$ but not to $M_{MB}(t)$ are denoted as:

\begin{equation}
M_{X}(t)=\sum_{i\in \mathcal{A}}\frac{1}{2^{N-1}}\left( \sum_{j\in \mathcal{A%
}\text{ (}j\neq i)}\left\vert \left\langle \beta _{j}\right\vert \hat{U}%
_{LE}^{{}}(t)\left\vert \beta _{i}\right\rangle \right\vert ^{2}-\sum_{j\in
\mathcal{B}}\left\vert \left\langle \beta _{j}\right\vert \hat{U}%
_{LE}^{{}}(t)\left\vert \beta _{i}\right\rangle \right\vert ^{2}\right).
\label{autocorrelacion52}
\end{equation}%
This balance of probabilities leads to the appropriate asymptotic behavior
of $M_{1,1}(t)$ according to the symmetries that constrain the evolution.
The identification

\begin{equation}
M_{1,1}(t)=M_{MB}(t)+M_{X}(t)  \label{LE_contributions}
\end{equation}%
is a crucial step for the following discussions.

If we use the identity $\hat{S}_{1}^{z}=\hat{S}_{1}^{+}\hat{S}_{1}^{-}-\frac{%
1}{2}\mathbf{\hat{I}}$ in equation (\ref{autocorrelacion2}), the invariance of
the trace under cyclic permutations ensures that $tr[\hat{S}_{1}^{z}(t)\hat{S%
}_{1}^{z}(0)]=tr[\hat{S}_{1}^{-}(0)\hat{S}_{1}^{z}(t)\hat{S}_{1}^{+}(0)]-%
\frac{1}{2}tr[\hat{S}_{1}^{z}(t)]$. Since $tr[\hat{S}_{1}^{z}(t)]=tr[\hat{S}%
_{1}^{z}(0)]=0$, then:

\begin{eqnarray}
M_{1,1}(t) &=&2\sum_{i}\frac{1}{2^{N-1}}\left\langle \beta _{i}\right\vert
\hat{S}_{1}^{-}(0)\hat{U}_{LE}^{\dag }(t)\hat{S}_{1}^{z}(0)\hat{U}%
_{LE}^{{}}(t)\hat{S}_{1}^{+}(0)\left\vert \beta _{i}\right\rangle  \notag \\
&=&2\sum_{i\in \mathcal{A}}\frac{1}{2^{N-1}}\left\langle \beta
_{i}\right\vert \hat{U}_{LE}^{\dag }(t)\hat{S}_{1}^{z}\hat{U}%
_{LE}^{{}}(t)\left\vert \beta _{i}\right\rangle ,
\label{M11_promedioEnsamble}
\end{eqnarray}%
which is indeed an explicit way to rewrite equation (\ref{autocorrelacion1}) in
the form of an ensemble average. Remarkably, since $\hat{S}_{1}^{z}$ is a
local (\textquotedblleft one-body\textquotedblright ) operator, its
evaluation in equation (\ref{M11_promedioEnsamble}) can be replaced by the
expectation value in a single superposition state \cite{Alv-parallelism},

\begin{equation}
M_{1,1}(t)=2\left\langle \Psi _{neq}\right\vert \hat{U}_{LE}^{\dag }(t)\hat{S%
}_{1}^{z}\hat{U}_{LE}^{{}}(t)\left\vert \Psi _{neq}\right\rangle ,
\label{autocorrelacion7}
\end{equation}%
where:
\begin{equation}
\left\vert \Psi _{neq}\right\rangle =\sum\limits_{i\in \mathcal{A}}\frac{1}{%
\sqrt{2^{N-1}}}e^{\mathrm{i}\varphi _{i}}\text{\ }\left\vert \beta
_{i}\right\rangle .  \label{neqsup}
\end{equation}%
Here, $\varphi _{i}^{{}}$ is a random phase uniformly distributed in $%
[0,2\pi )$. As a matter of fact, the state defined in equation (\ref{neqsup}) is
a random superposition that can successfully mimic the dynamics of ensemble
calculations and provides a quadratic speedup of computational efforts \cite%
{Alv-parallelism,Fine2013,Pineda2014}.

\section{The Loschmidt echo dynamics\label{Sec_locglob}}

\subsection{Short time expansions and beyond.}

In order to analyze the $N$-dependence of the LE and its time scales, we
compute here the short time expansion of the magnitudes $M_{1,1}(t)$, $%
M_{MB}(t)$ and $M_{X}(t)$. Up to $2^{nd}$ order in time,

\begin{eqnarray}
M_{1,1}(t=2t_{R}) &=&2\sum_{i\in \mathcal{A}}\frac{1}{2^{N-1}}\left\langle
\beta _{i}\right\vert \hat{U}_{LE}^{\dag }(t)\hat{S}_{1}^{z}\hat{U}%
_{LE}^{{}}(t)\left\vert \beta _{i}\right\rangle  \notag \\
&=&1-\left( t/\hbar \right) ^{2}\sum_{i\in \mathcal{A}}\frac{1}{2^{N-1}}%
\left( \left\langle \beta _{i}\right\vert \hat{\Sigma}^{2}\left\vert \beta
_{i}\right\rangle -2\left\langle \beta _{i}\right\vert \hat{\Sigma}\hat{S}%
_{1}^{z}\hat{\Sigma}\left\vert \beta _{i}\right\rangle \right) +\mathcal{O}%
\left( \left( t/\hbar \right) ^{3}\right) .  \label{chplg_st1}
\end{eqnarray}%
Similarly, the leading contributions to $M_{MB}(t)$ and $M_{X}(t)$ are:

\begin{eqnarray}
M_{MB}(t) &=&\sum_{i\in \mathcal{A}}\frac{1}{2^{N-1}}\left\vert \left\langle
\beta _{i}\right\vert \hat{U}_{LE}^{{}}(t)\left\vert \beta _{i}\right\rangle
\right\vert ^{2}  \notag \\
&=&1-\left( t/\hbar \right) ^{2}\sum_{i\in \mathcal{A}}\frac{1}{2^{N-1}}%
\left( \left\langle \beta _{i}\right\vert \hat{\Sigma}^{2}\left\vert \beta
_{i}\right\rangle -\left\langle \beta _{i}\right\vert \hat{\Sigma}\left\vert
\beta _{i}\right\rangle ^{2}\right) +\mathcal{O}\left( \left( t/\hbar
\right) ^{3}\right) ,  \label{chplg_st2}
\end{eqnarray}%
and%
\begin{eqnarray}
M_{X}(t) &=&\sum_{i\in \mathcal{A}}\frac{1}{2^{N-1}}\left( \sum_{j\in
\mathcal{A}\text{ (}j\neq i)}\left\vert \left\langle \beta _{j}\right\vert
\hat{U}_{LE}^{{}}(t)\left\vert \beta _{i}\right\rangle \right\vert
^{2}-\sum_{j\in \mathcal{B}}\left\vert \left\langle \beta _{j}\right\vert
\hat{U}_{LE}^{{}}(t)\left\vert \beta _{i}\right\rangle \right\vert
^{2}\right)  \notag \\
&=&\left( t/\hbar \right) ^{2}\sum_{i\in \mathcal{A}}\frac{1}{2^{N-1}}\left(
2\left\langle \beta _{i}\right\vert \hat{\Sigma}\hat{S}_{1}^{z}\hat{\Sigma}%
\left\vert \beta _{i}\right\rangle -\left\langle \beta _{i}\right\vert \hat{%
\Sigma}\left\vert \beta _{i}\right\rangle ^{2}\right) +\mathcal{O}\left(
\left( t/\hbar \right) ^{3}\right) .  \label{chplg_st3}
\end{eqnarray}

Let us consider a generic \textit{secular }(i.e. polarization conserving)
perturbation $\hat{\Sigma}$ given by a Hamiltonian with an arbitrary
anisotropy $\alpha $,%
\begin{equation}
\hat{\Sigma}=\sum_{i,j}^{N}(J_{\Sigma})_{ij}\left[ 2\alpha \hat{S}_{i}^{z}%
\hat{S}_{j}^{z}-\left( \hat{S}_{i}^{x}\hat{S}_{j}^{x}+\hat{S}_{i}^{y}\hat{S}%
_{j}^{y}\right) \right] .  \label{chplg_pert}
\end{equation}%
This is still quite general since even a double quantum perturbation ($\hat{S%
}_{i}^{+}\hat{S}_{j}^{+}+\hat{S}_{i}^{-}\hat{S}_{j}^{-}$, which does not
conserve polarization) can be reduced to a secular one by the truncating
effects of radiofrequency fields \cite{Zangara2015PRA}. In addition, we do
not consider here the case $\left[ \hat{\Sigma},\hat{S}_{1}^{z}\right] =0$
(e.g. pure Ising or on-site diagonal disorder), since in such a condition
the first non-trivial order in time is the $4^{th}$ (see Appendix). Then,
the following identities hold:%
\begin{eqnarray}
\sum_{i\in \mathcal{A}}\frac{1}{2^{N-1}}\left\langle \beta _{i}\right\vert
\hat{\Sigma}^{2}\left\vert \beta _{i}\right\rangle &=&2N\sigma ^{2}\left(
\frac{\alpha ^{2}}{4}+\frac{1}{8}\right) ,  \label{chplg_st4} \\
\sum_{i\in \mathcal{A}}\frac{1}{2^{N-1}}\left\langle \beta _{i}\right\vert
\hat{\Sigma}\hat{S}_{1}^{z}\hat{\Sigma}\left\vert \beta _{i}\right\rangle
&=&2N\sigma ^{2}\left( \frac{\alpha ^{2}}{8}+\frac{1}{16}\right) -\frac{1}{2}%
\sigma ^{2},  \label{chplg_st5} \\
\sum_{i\in \mathcal{A}}\frac{1}{2^{N-1}}\left\langle \beta _{i}\right\vert
\hat{\Sigma}\left\vert \beta _{i}\right\rangle ^{2} &=&2N\sigma ^{2}\frac{%
\alpha ^{2}}{4}.  \label{chplg_st6}
\end{eqnarray}%
Here, $\sigma ^{2}$ stands for the average local second moment of $\hat{%
\Sigma}$,%
\begin{equation}
\sigma _{{}}^{2}=\frac{1}{N}\sum_{i=1}^{N}\sigma _{i}^{2}=\frac{1}{N}%
\sum_{i=1}^{N}\left[ \sum_{j(\neq i)}^{N}\left( \frac{(J_{\Sigma})_{ij}}{2}%
\right) ^{2}\right] .  \label{chplg_sigmalocal}
\end{equation}%
Being a \textit{perturbation}, we notice that $\sigma ^{2}$ is much smaller
than the average local second moment $\sigma _{0}^{2}$ of the unperturbed
Hamiltonian $\hat{H}_{0}$. In terms of time scales,
\begin{equation}
T_{2}=\hbar /\sqrt{\sigma _{0}^{2}}\ll \hbar /\sqrt{\sigma _{{}}^{2}}=\tau
_{\Sigma }.  \label{time-comparison}
\end{equation}%
The identities in equations (\ref{chplg_st4}), (\ref{chplg_st5}) and (\ref%
{chplg_st6}) lead to

\begin{equation}
M_{1,1}(t)=1-\left( t/\tau _{\Sigma }\right) ^{2}+\mathcal{O}\left( \left(
t/\hbar \right) ^{3}\right) ,  \label{chplg_st7}
\end{equation}%
and,

\begin{eqnarray}
M_{MB}(t) &=&1-\frac{1}{4}N\left( t/\tau _{\Sigma }\right) ^{2}+\mathcal{O}%
\left( \left( t/\hbar \right) ^{3}\right) ,  \label{chplg_st8} \\
M_{X}(t) &=&\left( \frac{N-4}{4}\right) \left( t/\tau _{\Sigma }\right) ^{2}+%
\mathcal{O}\left( \left( t/\hbar \right) ^{3}\right) .  \label{chplg_st9}
\end{eqnarray}%
These expansions hold for $t<\left( \tau _{\Sigma }/N\right) $. Beyond such
a very short time regime, a general term in the expansion of $M_{1,1}(t)$
will be of the form

\begin{equation}
c_{(N,n)}t^{n}/\left( \tau _{\Sigma }^{k}T_{2}^{n-k}\right)
\label{termino_gral}
\end{equation}%
with $k\geq 2$ and the coefficient $c_{(N,n)}$ described by combinatorial
numbers of increasing size that depend on the topology of the interactions
(e.g. see \cite{Claudia2014,Zamar2015}). Since the experimental set up
corresponds to the limit described by equation (\ref{time-comparison}), this
expansion will be dominated by terms with the lowest possible order in the
weak interaction, i.e. $k=2$:%
\begin{equation}
(t/\tau _{\Sigma })^{2}\left[ 1+\sum_{n}c_{(N,n)}(t/T_{2}^{{}})^{n-2}\right]
.  \label{orden_dominante}
\end{equation}%
Equation (\ref{orden_dominante}) indicates that beyond the very short time
expansion, i.e. $(\tau _{\Sigma }/N)<t<\tau _{\Sigma }$, the dependence on $%
\tau _{\Sigma }^{{}}$ becomes superseded by the diverging terms in the scale
$T_{2}^{{}}$. This could lead to the new time scale $T_{3}$ which was seen
experimentally to be tied to $T_{2}$ as

\begin{equation}
T_{2}\lesssim T_{3}\ll \tau _{\Sigma }.  \label{comparison2}
\end{equation}%
In that sense, $T_{3}$ becomes characteristic of the complexity or
\textquotedblleft chaos\textquotedblright\ of the many spin system that
amplifies the small effect of the perturbation. In addition, it is important
to stress that, being an experimental fact, equation (\ref{comparison2})
corresponds to a system composed by infinitely many interacting spins. In
other words, equation (\ref{comparison2}) stands for the relations of time scales
in the thermodynamic limit. Quite on the contrary, any numerical simulation
involves a finite, very small indeed, number of spins where the
irreversibility rate $T_{3}$ would be essentially given by $\tau _{\Sigma}$%
. Then, the LE decay rate evaluated in a finite system would ultimately be
perturbation-dependent \cite{Fine2014}. Thus, our \textit{Central Hypothesis
of Irreversibility} would mean that equation (\ref{comparison2}) is an \textit{%
emergent} property. It should rely on the thermodynamic limit, which\
implies taking the limit $N\rightarrow \infty $ first, and then $\tau_{\Sigma}\rightarrow \infty $. The non-uniformity of these limits plays a crucial role to yield quantum phase transitions, as discussed in the context
of Anderson localization \cite{AndersonRMP-1978,Rodrigues1986,hmp-physB}.

The physical picture described above is schematically represented in Fig. %
\ref{dynamicalScheme}. There, we show the expected interplay between $%
M_{MB}(t)$ and $M_{X}(t)$ leading to $M_{1,1}(t)$. Indeed, as stated in
equations (\ref{chplg_st8}) and (\ref{chplg_st9}), the very short time
dependence of both contributions is extensive in $N$: $M_{MB}(t)$ decreases
as $1-N\sigma _{{}}^{2}t^{2}/4$ and $M_{X}(t)$ increases as $(N-4)\sigma
_{{}}^{2}t^{2}/4$. Such a precise balance provides for the short time decay
of $M_{1,1}(t)$ given by equation (\ref{chplg_st7}), i.e. $1-\sigma
_{{}}^{2}t^{2} $. Notice that there is no reason to assume that the decay of
$M_{MB}(t)$ would remain ruled by $\tau _{\Sigma }$. Beyond the very short
times, we expect that the time scale $T_{3}$ should also show up as $%
T_{3}/N^{\nu }$ with $\nu \sim 1$ in the decay of $M_{MB}(t)$(see below).
Furthermore, while $M_{MB}(t)$ goes monotonically to zero, $M_{X}(t)$
displays a highly non-trivial behavior. Indeed, $M_{X}(t)$ first increases
by feeding from the decay of $M_{MB}(t)$ until it reaches a maximum. This
growth indicates a progressive divergence of long-range correlations.
Thereafter, $M_{X}(t)$ should decay accounting for the fact that the state
remains properly normalized. This is precisely what $M_{1,1}(t)$ measures: a
conserved polarization that ultimately distributes uniformly within the spin
system. In an isolated finite system this implies the asymptotic plateau $%
M_{\infty }\sim 1/N$. As pointed above, the decay of both $M_{1,1}(t)$ and $%
M_{X}(t)$ occurs in a time scale $T_{3}$, which according to equation (\ref%
{orden_dominante}), is somewhat longer but close to the \textquotedblleft
diffusion\textquotedblright\ time $T_{2}$. This is the regime captured
experimentally.

\begin{figure}[h]
\centering\includegraphics[width=4.5in]{fig2.eps}
\caption{A pictorial scheme of the time dependence of $M_{1,1}(t)$ (dotted line) and its contributions, i.e. $M_{MB}(t)$ and $M_{X}(t)$ (solid lines). Their short-time expansions, as stated in equations  (\protect\ref{chplg_st7}), (\protect\ref{chplg_st8}) and (\protect\ref{chplg_st9}), are indicated with arrows. In particular, the expansion corresponding to the short-time behavior of $M_{1,1}(t)$ is plotted with a dashed line. (Online version in colour.)}
\label{dynamicalScheme}
\end{figure}

\subsection{The extensive decay hypothesis\label{Sec_scaling}}

The previous short-time expansions provide a hint on the scaling relation
between the local LE, $M_{1,1}(t),$ and the global one as embodied by $%
M_{MB}(t)$. In particular, let us first compare the probability of
refocusing the configuration (up or down) of a single spin, i.e. $\Pi
_{1,1}(t)$, and the probability of refocusing a complete many-spin state $%
M_{MB}(t)$. If the refocusing of each individual spin could be treated as an
independent event, then the scaling between $\Pi _{1,1}$ and $M_{MB}$ would
be extensive in $N$,

\begin{equation}
\left( \Pi _{1,1}(t)\right) ^{N/2}\simeq M_{MB}(t).
\label{scaling_locglobal}
\end{equation}%
Here, the factor $1/2$ in the exponent comes from equation (\ref{rho_inicial}),
i.e. the initial high temperature state, where basically half of the spins
point up, and half of them point down.\textbf{\ }Then, one can resort to the
picture of a lattice gas where\textbf{\ }$N/2$ particles jump among $N$
lattice sites. As in the well known Jordan-Wigner transformation \cite%
{Lieb-Mattis1961}, a fermion is associated to a spin pointing up and a
vacancy corresponds to a spin pointing down. Thus, the microstate of the gas
is completely described by the position of $N/2$ particles.

Strictly speaking, the notion of \textit{extensiveness} corresponds to standard thermodynamic quantities such as the entropy of the system. In addition, as discussed in \cite{MolPhys}, $S=-ln(M_{1,1}(t))$ is precisely a measure of the entropy. Then, the validity of equation (\ref{scaling_locglobal}) implies an extensivity relation between the entropy per spin and the total entropy of the system.

According to equations (\ref{autoco2}) and (\ref{chplg_st7}),

\begin{equation}
\Pi _{1,1}(t)=1-\frac{1}{2}\left( t/\tau _{\Sigma }\right) ^{2}+\mathcal{O}%
\left( \left( t/\hbar \right) ^{3}\right) ,  \label{scaling_aux}
\end{equation}%
which in turn, up to $2^{\mathtt{nd}}$ order in time, implies $\Pi _{1,1}(t)\simeq
\left( M_{1,1}(t)\right) ^{1/2}$. Thus, equation (\ref{scaling_locglobal}) yields

\begin{equation}
\left( M_{1,1}(t)\right) ^{N/4}\simeq M_{MB}(t).  \label{scaling_locglobal2}
\end{equation}%
This is precisely the relation verified between equations (\ref{chplg_st7}) and (%
\ref{chplg_st8}).

One might expect that beyond the very short-time decay, individual spin
autocorrelations deviate from the statistical independence. However, this
deviation will still have a local nature and therefore the $N$-extensivity
would remain valid. Indeed, we propose%
\begin{equation}
\left( M_{1,1}(t)\right) ^{\eta }\simeq M_{MB}(t),
\label{scaling_locglobal3}
\end{equation}%
where the exponent $\eta $ would be some appropriate function $\eta =\eta
(N,t)$. Our \textquotedblleft extensive decay hypothesis\textquotedblright\
implies that $\eta $ factorizes:%
\begin{equation}
\eta (N,t)=N\times f(t),  \label{factoriza_exponente}
\end{equation}%
where $f(t)$ stands for a function that encloses information of the
correlations originated by the system dynamics. Additionally,

\begin{equation}
\lim_{t\rightarrow 0^{+}}f(t)=\frac{1}{4}.  \label{limite_ft}
\end{equation}%
is required in order to recover equation (\ref{scaling_locglobal2}), i.e. the
statistical independence.

\section{A 1D model.\label{Sec_1Dmodel}}

The physical picture described above is discussed here under the light of a
specific model. In particular, we assess the validity of equations (\ref%
{scaling_locglobal2}) and (\ref{factoriza_exponente}). We consider a 1-D
spin chain with an anisotropic interaction described by:

\begin{equation}
\hat{H}_{0}=\sum_{i=1}^{N-1}J_{0}\left( \frac{1}{2}\hat{S}_{i}^{z}\hat{S}%
_{i+1}^{z}+\hat{S}_{i}^{x}\hat{S}_{i+1}^{x}+\hat{S}_{i}^{y}\hat{S}%
_{i+1}^{y}\right) ,  \label{chplg_H1vecinos}
\end{equation}%
with periodic boundary conditions, i.e. a ring configuration. Here, $J_{0}$
stands for the natural units of the spin-spin interaction energy. As a
perturbation $\hat{\Sigma}$ we choose a next nearest neighbors interaction
described by:

\begin{equation}
\hat{\Sigma}=\sum_{i=1}^{N-2}J_{\Sigma}\left( \frac{1}{2}\hat{S}_{i}^{z}\hat{%
S}_{i+2}^{z}+\hat{S}_{i}^{x}\hat{S}_{i+2}^{x}+\hat{S}_{i}^{y}\hat{S}%
_{i+2}^{y}\right) .  \label{chplg_H2vecinos}
\end{equation}
Such a perturbation appears naturally when one attempts to build an
effective one-body dynamics from linear crystals with dipolar interactions
\cite{elenaMQC}. This is also the case in a regular crystal, when the
natural non-secular dipole-dipole terms are truncated by the Zeeman energy
of the radiofrequency irradiation, which ultimately leads to effective
secular two-body next nearest neighbors interactions \cite{Zangara2015PRA}.

The local second moments $\sigma ^{2}$ and $\sigma _{0}^{2}$ of $\hat{\Sigma}
$ and $\hat{H}_{0}$ respectively can be evaluated as in equation (\ref%
{chplg_sigmalocal}):

\begin{eqnarray}
\sigma ^{2} &=&\frac{1}{2}\left( J_{\Sigma} \right) ^{2},
\label{chplg_sigmaModelo} \\
\sigma _{0}^{2} &=&\frac{1}{2}\left( J_{0}^{{}}\right) ^{2},
\label{chplg_sigmaH0}
\end{eqnarray}%
and constitute the main energy scales of our problem.

In Fig. \ref{Ecos_dinamica} we plot $M_{1,1}(t)$, $M_{MB}(t)$ and $M_{X}(t)$
for the particular choice $J_{\Sigma }=0.1J_{0}^{{}}$. Short-time expansions
given in equations (\ref{chplg_st7}), (\ref{chplg_st8}) and (\ref{chplg_st9})
are evaluated according to equation (\ref{chplg_sigmaModelo}). It is observed
that $M_{MB}(t)$ vanishes for long times.\ Actually, a close observation
shows that $M_{MB}(t\rightarrow \infty )\sim \mathcal{O}(2^{-N})$ (data not
shown). In addition, notice that $M_{X}(t\rightarrow \infty )\sim 1/N$. Such
an asymptotic contribution provides for the equidistribution of the spin
polarization $M_{1,1}(t\rightarrow \infty )\sim 1/N$. This long-time
saturation corresponds to the equilibration of a finite system.

In contrast to our schematic plot in Fig. \ref{dynamicalScheme}, here $%
M_{X}(t)$ does not get too close to $1$ and $M_{MB}(t)$ does not decay much
faster than $M_{1,1}(t)$. Since $M_{X}(t)$ provides for the whole $%
M_{1,1}(t) $ once that $M_{MB}(t)$ has fully decayed, the contribution of $%
M_{X}(t)$ is considerable only at long times. These effects are a
consequence of the relatively small size of the system considered. Indeed,
the case in Fig. \ref{Ecos_dinamica} corresponds to $N=14$ spins, and thus
the exponent that relates $M_{1,1}(t)$ and $M_{MB}(t)$ is quite small $%
(N/4)=3.5$. The need for larger systems indicates that revealing the
dominant orders in equation (\ref{orden_dominante}) is a major numerical
challenge that may go beyond the state-of-the-art techniques \cite%
{Dente2013CPC}.

\begin{figure}[h]
\centering\includegraphics[width=3.8in]{fig3.eps}
\caption{The local LE and its non local contributions. $M_{1,1}(t)$, $M_{MB}(t)$ and $M_{X}(t)$ correspond to the solid lines as indicated by the labels in the figure. $N=14$, $J_{\Sigma }=0.1J_{0}^{}$. The short-time expansions given in equations (\protect\ref{chplg_st7}), (\protect\ref{chplg_st8}) and (\protect\ref{chplg_st9}) are shown in dashed lines. (Online version in colour.)}
\label{Ecos_dinamica}
\end{figure}

In order to assess the accuracy of the \textquotedblleft extensive decay
hypothesis\textquotedblright , in Fig. \ref{Fig_factorizaexp} we address the
scaling relation between $M_{1,1}(t)$ and $M_{MB}(t)$ discussed in Sec. \ref%
{Sec_scaling}. In particular, we try out the factorization stated in equation (%
\ref{factoriza_exponente}). By plotting $\log (M_{MB}(t))/\left( \log
(M_{1,1}(t))N\right) $ as a function of time, we observe a unique function
which does not depend on $N$ or $J_{\Sigma }$, but it has a weak dependence
on time. Such a unique curve is indeed $f(t)$ as defined in equation (\ref%
{factoriza_exponente}). This means that the extensivity relation between $%
M_{1,1}(t)$ and $M_{MB}(t)$ is confirmed. The statistical independence, in
turn, fails progressively once $f(t)$ departs from the $1/4$ factor of the
ideal relation in equations (\ref{scaling_locglobal2}) and (\ref{limite_ft}).
Since beyond the short-time regime $f(t)$ decreases with time, we conclude
that the recovery of a single spin is tied to the recovery of its neighbors.
Thus, the spins are positively correlated and the revival probability of the
complete $N$-spin state is enhanced. This argument is particularly relevant
in 1D systems.

After the onset of the saturation regime, where $M_{1,1}\sim 1/N$ and $%
M_{MB}\sim \mathcal{O}(2^{-N})$, the universal scaling naturally becomes
noisy and the curves for different $N$ and $J_{\Sigma }$ separate each
other. Since the decay is faster for larger perturbations, the appearance of
such a spurious behavior is observed to occur first for the largest value of
$J_{\Sigma }$ considered ($J_{\Sigma }=0.3J_{0}^{{}}$, plus signs and triangles).

\begin{figure}[h]
\centering\includegraphics[trim=0 0 0 0 ,clip,width=4.5in]{fig4.eps}
\caption{The relation $\log (M_{MB}(t))/\left( N\log (M_{1,1}(t))\right)$ as a function of time. For $J_{\Sigma }=0.1J_{0}^{}$ the sizes plotted are: $N=10$ (circles), $N=12$ (squares), $N=14$ (diamonds), and $N=16$ (stars). For $J_{\Sigma }=0.2J_{0}^{}$ the sizes plotted are: $N=10$ (solid line), $N=12$ (dashed line), $N=14$ (dash-dot line), and $N=16$ (dotted line). For $J_{\Sigma }=0.3J_{0}^{}$ the sizes plotted are: $N=10$ (up triangles), $N=12$ (down triangles), $N=14$ (plus signs), and $N=16$ (right triangles). (Online version in colour.) }
\label{Fig_factorizaexp}
\end{figure}

\section{Conclusion \label{Sec_conclu}}

We presented a detailed analysis of the LE in interacting spin systems. As
in the NMR experiments, a local version of the LE, $M_{1,1}$, is defined as
a single spin autocorrelation function. Simultaneously, we define a global
LE, $M_{MB}$, as the average of the square overlap between many-body wave
functions that evolved under perturbed Hamiltonians. While the former
constitutes a specific experimental observable, the latter has only been
assessed theoretically. Here, we showed the formal relation between both
magnitudes, as far as their characteristic time scales and $N$-dependence
are concerned.

By analyzing a short-time expansion of $M_{1,1}$and $M_{MB}$ we derived a
precise the relation between their time scales. In this regime, the decay of
$M_{1,1}$ is given by the average local second moment of the perturbation ($%
\hbar /\tau _{\Sigma }=\sqrt{\sigma _{{}}^{2}}$), and the decay of $M_{MB}$\
by $N$ times the local scale ($N\hbar /\tau _{\Sigma }.$). This relation
hints a scaling law $M_{MB}\sim \left( M_{1,1}\right) ^{N/4}$ that accounts
for the extensivity of $M_{MB}$. In such a case, the recovery of a many-spin
state results from the recovery of each individual spin, much as if they
were independent events. The numerical evaluation in a specific spin model
shows that the exponent slightly diminishes with time, starting from the
initial $N/4$. This means that the recovery of a single spin is positively
correlated with the probability of recovery of its neighbors, and thus it
improves the probability of the revival of the complete $N$-spin state. A
precise control of these correlations may hint an experimental access to the
global autocorrelation, i.e. $M_{MB}$, just by measuring a single spin
(local) autocorrelation $M_{1,1}$. This would require an experimental
protocol capable to encode a local excitation into a correlated many-spin
state.

In addition, we discussed a general dynamical picture beyond the very
short-time regime. There, the decay of $M_{1,1}$ results from the interplay
between the time scale that characterizes the reversible interactions ($%
T_{2} $) and that of the perturbation ($\tau _{\Sigma }$). This would
ultimately lead to the experimentally observed $T_{3}$, which was found to
be roughly independent of $\tau _{\Sigma }$ but closely related to $T_{2}$.
The theoretical quest for the emergent $T_{3}$ time scale remains open and
it may be out of the reach of current numerical approaches. Assessing a fair
estimate analytically would require a detailed account of the higher order
processes that dress the quadratic term in the perturbative expansion.

Notice that our discussion lead us to identify $T_{3}$, and hence the
spin-spin interaction time $T_{2}$, as the time scales characterizing the
complexity or many-spin chaos. As such, they show up not only in the decay
of $M_{1,1}$ and $M_{MB}$, but also in the growth of $M_{X}=M_{1,1}-M_{MB}.$
Indeed, in the field of AdS/CFT there is an increasing interest in
characterizing the role of chaos in quantum dynamics \cite%
{Stanford_butterfly,Susskind2008,Susskind2015,Kitaev_charla}. There, chaos
manifests in the growth of four-body correlation functions, following an
early suggestion by Larkin and Ovchinnikov \cite{larkin1969}. They employed
semiclassical arguments to address disordered superconductors and probed
that the square dispersion of momentum should grow exponentially in a time
scale determined by the collisions with impurities, i.e. with the
unperturbed Hamiltonian (in our physical picture, $T_{2}$). Similarly, our
average multispin correlation $M_{X}$ would ultimately diverge within a time
scale $T_{3}/N$, i.e. independent of the perturbation. This is indeed a
measure of the decoherence, and hence of irreversibility, induced by
many-spin chaos. Of course, we do not have a precise characterization of
this time scale or the specific mathematical dependence on time. Thus, this
is a puzzling issue to explore in the field of many-body chaos. Besides the
obvious relevance for statistical mechanics and experimental physics, this
might also contribute to a possible pathway between quantum mechanics and
gravity.

\section*{Acknowledgments}
This work benefited from discussions with A.D. Dente and F. Pastawski. HMP
greatly acknowledges hospitality of A. Kitaev at Caltech, P.A. Lee at MIT
and V. Oganesyan at CUNY, where the issues discussed in this paper acquired
certain maturity. PRZ acknowledges M.C. Ba\~{n}uls and J.I. Cirac for their
kind hospitality at MPQ in Garching. We acknowledge financial support from CONICET, ANPCyT, SeCyT-UNC and
MinCyT-Cor. This work used computational resources from CCAD -- Universidad
Nacional de C\'{o}rdoba (http://ccad.unc.edu.ar/), in particular the
Mendieta Cluster, which is also part of SNCAD -- MinCyT, Rep\'{u}blica
Argentina.

\section{Appendix}

It is worthy to mention that very short-time expansions in equations (\ref%
{chplg_st7}), (\ref{chplg_st8}) and (\ref{chplg_st9}) do not depend on the
anisotropy $\alpha $ of the perturbation. In general, it can be proved that
if $\left[ \hat{\Sigma},\hat{S}_{1}^{z}\right] =0$ then:

\begin{equation}
M_{1,1}(t)=1-\frac{\left( t/\hbar \right)^{4}}{2^{N+3}}\sum_{i\in \mathcal{A}%
}\left( 2\left\langle \beta _{i}\right\vert \left[ \hat{\Sigma},\hat{H}_{0}%
\right] \hat{S}_{1}^{z}\left[ \hat{\Sigma},\hat{H}_{0}\right] \left\vert
\beta _{i}\right\rangle -\left\langle \beta _{i}\right\vert \left[ \hat{%
\Sigma},\hat{H}_{0}\right] ^{2}\left\vert \beta _{i}\right\rangle \right)+%
\mathcal{O}\left( \left( t/\hbar \right) ^{5}\right).  \label{a15}
\end{equation}%
This is precisely the case of a perturbation $\hat{\Sigma}$ enclosing
Anderson disorder and Ising interactions \cite{Zangara2013PRB} or
interactions with a fluctuating field \cite{LucasF}.


\section{Linear And Nonlinear Adaptive Filter}
The term \emph{filter} usually refers to a system that is designed to extract information about a prescribed quantity of interest from noisy data. An \emph{adaptive filter} is a filter that self-adjusts its input output mapping according to an optimization algorithm like predictive coding~\cite{Huang11} mostly driven by an error signal. Because of the complexity of the optimization algorithms, most adaptive filters are digital filters. With the available processing capabilities of current digital signal processors, adaptive filters have become much more popular and are now widely used in various fields such as sound wave based communication devices~\cite{Huang15}, face extraction~\cite{wu2011learning,wu2014learning}, camcorders and digital cameras, information retrieval~\cite{Nan15} and medical monitoring equipments.

\subsection{Linear Adaptive Filter}
In most adaptive signal processing applications, system linearity is assumed and adaptive linear filters are thus used. The traditional class of supervised adaptive filters rely on error-correction learning for their adaptive capability. Therefore, the error is the necessary element of a cost function, which is a criterion for optimum performance of the filter.

The linear filter includes a set of adaptively adjustable parameters (also known as weights), which is marked as $\bm \omega(n-1)$, where $n$ denotes discrete time. The input signal $\bm u(n)$ applied to the filter at time $n$, produces the actual response $y(n)$ via
\begin{equation}
y(n) = \bm \omega(n-1)^T \bm u(n)  \nonumber
\end{equation}
Then this actual response is compared with the corresponding desired response $d(n)$ to produce the error signal $e(n)$. The error signal, in turn, acts as a guide to adjust the weights $\bm \omega(n-1)$ by an incremental value denoted by $\triangle \bm \omega(n)$. On the next iteration, $\bm \omega(n)$ becomes the latest value of the weights to be updated. The adaptive filtering process is repeated continuously until the filter reaches a stop condition, which normally is that the weights adjustment is small enough.

An important issue in the adaptive design, no matter linear or nonlinear adaptive filter, is to ensure the learning curve is convergent with an increasing number of iterations. Under this condition, we define the system is in a steady-state.

\subsection{Nonlinear Adaptive Filter}
Even though linear adaptive filtering can approximate non-linearity, the performance of adaptive linear filters is not satisfactory in applications where nonlinearities are significant. Hence more advanced nonlinear models are required. At present, Neural Networks and \emph{Kernel Adaptive Filters} are popular nonlinear models.  By providing linearity in a high dimension feature space, Reproduce Kernel Hilbert Space (RKHS), and universal approximation in Euclidean space with universal kernels\cite{kernelbook}, the kernel adaptive filters are attracting more attention. Through a reproducing kernel, kernel adaptive filters map data from an input space to RKHS, where appropriate linear methods are applied to the transformed data. This procedure implements a nonlinear treatment for the data in the input space. Comparing with other nonlinear techniques, kernel adaptive filters have the following features:
\begin{itemize}
  \item They can be universal approximations whenever the kernel is universal.
  \item They have no local minima with the squared error cost function.
  \item They have moderate complexity in terms of computation and memory.
  \item They belong to online learning method~\cite{Huang14} and have good tracking ability to handle nonstationary conditions.
\end{itemize}
The details about kernel adaptive filters are introduced next.

\section{Reproducing Kernel Hilbert Space }
Reproducing Kernel Hilbert Space, RKHS for short, is a complete inner product space associate with a Mercer kernel.  A Mercer kernel is a continuous, symmetric and positive definite function $\kappa$ : $\mathbb U \times \mathbb U \rightarrow \mathbb R $, where $\mathbb U$ is the input domain in Euclidean space $\mathbb R^L$ ($L$ is the input order). The commonly used kernels includes the Gaussian kernel (\ref{equ:gaussianKernel}) and the polynomial kernel (\ref{equ:gpolyKernel})  \cite{kernelbook}.
\begin{equation}
\kappa(\bm u, \bm u^{\prime}) = \text{exp}(-\frac{\| \bm u - \bm u^{\prime} \|^2}{\sigma^2})
\label{equ:gaussianKernel}
\end{equation}
\begin{equation}
\label{equ:gpolyKernel}
\kappa(\bm u, \bm u^{\prime})= (\bm u^T \bm u^{\prime}+1)^p
\end{equation}
If only one free parameter is inserted into the kernel function, $\kappa(\bm u, \cdot)$ is expressed as a transformed feature vector  $\varphi(\bm u)$ through a mapping $\bm\varphi: \mathbb U \rightarrow \mathbb H $, where $\mathbb H$ is the RKHS. Therefore $\bm\varphi(\bm u) = \kappa(\bm u, \cdot) $.  One of the most important properties of RKHS for practical applications is the one, called the ``kernel trick''
\begin{equation}
\bm \varphi(\bm u)^T \bm \varphi(\bm u^{\prime})= \kappa(\bm u, \bm u^{\prime})
\end{equation}
that allows computing the inner products between two RKHS functions as a scalar evaluation in the input space by the kernel.

Besides ``kernel trick'', some other properties of RKHS related to this work are as follows. Assume $\mathbb H$ be any RKHS of all real-valued functions of $\bm u$ that are generated by the $\kappa(\bm u, \cdot)$. Suppose now two functions $h(\cdot)$ and $g(\cdot)$ are picked from the space $\mathbb H$ that are respectively represented by
\begin{equation}
h(\cdot) = \sum_{i=1}^l a_i\kappa(\bm c_i,\cdot)= \sum_{i=1}^l a_i \bm \varphi (\bm c_i) \nonumber
\end{equation}
\begin{equation}
g(\cdot) = \sum_{j=1}^m b_j\kappa(\tilde{\bm c}_j,\cdot)= \sum_{j=1}^m b_j \bm \varphi (\tilde{\bm c}_j)\nonumber
\end{equation}
where $a_i$ and $b_j$ are the expansion coefficients and both $\bm c_i$ and $\tilde{\bm c}_j \in \mathbb U$ for all $i$ and $j$.
\begin{enumerate}
  \item Symmetry
   \begin{equation}
   <h,g> = <g,h> \nonumber
   \end{equation}
  \item Scaling and distributive property
    \begin{equation}
    <(cf+dg),h> = c<f,h>+d<g,h> \nonumber
    \end{equation}
  \item Squared norm
    \begin{equation}
    \|h\|^2 = <h,h> \geq 0 \nonumber
    \end{equation}
\end{enumerate}
\section{Kernel Adaptive Filtering}
The kernel method is a powerful nonparametric modeling tool for pattern analysis and statistical signal processing. Through a nonlinear mapping, kernel methods transform the data into a set of points in a RKHS. Then various methods are utilized to find relationship between the data. There are many successful examples of this methodology including support vector machines (SVM) \cite{svm}, kernel regularization network \cite{krn}, kernel principal component analysis (kernel PCA) \cite{kcomponent} and kernel fisher discriminant analysis \cite{kernelfisher}. Kernel adaptive filters is a class of kernel methods. As a member of kernel adaptive filter, kernel affine projection algorithms (KAPA) include the kernel least mean square (KLMS) \cite{klms} as the simplest element and kernel recursive least squares (KRLS) \cite{KRLS} as the most computationally demanding.

The main idea of kernel adaptive filtering can be summarized as follows: Transform the input data into a high-dimension feature space $\mathbb F$, via a Mercer kernel. Then appropriate linear methods are subsequently applied on the transformed data. As long as a linear method in the feature space can be formulated in terms of inner products, there is no need to do computation in high-dimension space basing on ``kernel trick''. It has been proved that the kernel adaptive filters with universal kernel has universal approximation property. i.e. for any continuous input-output mapping $f: \mathbb{U}\rightarrow \mathbb{R}$, $\forall \varsigma >0, \exists \{ \bm u(i)\}_{i\in N} \in \mathbb{U}$ and real number ${c_i}_{i\in N}$, such that $\| f-\sum_i c_i\kappa(\bm u(i), \cdot)\|_2 < \varsigma$. This universal approximation property guarantees that the kernel method is capable of superior performance in nonlinear tasks. If we express a vector in $\mathbb F$ as
\begin{equation}
\bm \Omega=\sum_i c_i \bm \varphi(\bm u(i))
\end{equation}
we obtain
\begin{equation}
\|f-\Omega^T \bm \varphi\|_2 < \varsigma
\end{equation}

Furthermore, in the view of supervised learning, which requires the availability of a collection of desired responses, error cost functions (or error criteria) play significant role. Kernel adaptive filters provide a generalization of linear adaptive filters because these are a special case of the former when expressed in the dual space. Kernel adaptive filters exhibit a growing radial basis function network, learning the network topology and adaptive the free parameters directly from the data at the same time \cite{kernelbook}. In the following of this section, KRLS and KLMS are introduced respectively.

Consider the learning of a nonlinear function $f: \mathbb{U}\rightarrow \mathbb{R}$ based on a known sequence $(\bm{u}(1), d(1)), (\bm{u}(2), d(2)),\ldots, (\bm{u}(n), d(n))$, where $\mathbb{U}\in \mathbb{R}^L $ is the input space, $\bm{u}(i)$, $i= 1,...,n$ is the system input at sample time $i$, and $d(i)$ is the corresponding desired response.
\subsection{Kernel Recursive Least Square Algorithm}
The KRLS is actually the recursive least squares algorithm (RLS) algorithm in RKHS. At each iteration, one needs to solve the regularized least squares regression to obtain $f$:
\begin{equation}
\min_{f \in \mathbb{H}} \sum_{i=1}^n \|f(\bm u(i))-d(i)\|^2 + \lambda \|f\|_\mathbb{H}^2
\end{equation}
where $\lambda$ is the regularization term and $\|\cdot\|^2_\mathbb{H}$ denotes the norm in $\mathbb{H}$.

This problem can alternatively be solved in a feature space and results in KRLS. The learning problem of KRLS with regularization in feature space $\mathbb{F}$ is to find a high-dimensional weight vector $\bm \Omega \in \mathbb{F}$ that minimizes
\begin{equation}
\label{Equ:KRLSRegu}
\min_{\bm \Omega \in \mathbb{F}} \sum_{i=1}^n \|\bm \Omega^T \bm \varphi(\bm u(i))-d(i)\|^2 + \lambda \| \bm \Omega\|^2_\mathbb{F}
\end{equation}
where $\|\cdot\|^2_\mathbb{F}$ denotes the norm in $\mathbb{F}$.
\subsubsection{Approximate linear dependency for sparsification}
Similar to RLS, KRLS achieves high accuracy and has fast convergence rate in stationary scenario. However the good performance is obtained at a  cost of high computation complexity, $O(n^2)$, $n$ is the number of processed sample.

In order to decrease the computational complexity of KRLS, sparisification techniques are adopted. Specifically, sparsification in kernel methods is related to less computational complexity and memory consumption. Furthermore, the system's generalization ability is also influenced by sparsification in machine-learning algorithms and signal processing \cite{kernelbook}. Novelty Criterion, Approximate Linear Dependency (ALD) \cite{KRLS} , Prediction Various, Surprise and Quantized techniques are common strategies for sparisification. Among them, ALD is an effective sparisification technique for KRLS because it is solved in the feature space, unlike most of the other techniques.  Before introducing ALD, several concepts should be interpreted:

Network size: the number of data utilized to describe  the system model.

Center $\bm c$:  the data utilized to build the system. Therefore, the number of centers equals to network size.

Center dictionary $\mathcal{C}$: the set of all centers $\bm c$.

Suppose after having observed $n-1$ training samples, we have established a center dictionary $\mathcal{C}(n-1) = \{\bm{ c}_i\}_{i=1}^{K(n-1)}$, where $K(n-1)$ is the cardinality of $\mathcal{C}(n-1)$. When a new sample $\{\bm u(n), d(n)\}$ is presented, ALD tests whether there exists a coefficient $\bm a(n) =(a_1, \ldots, a_{K(n-1)})$ satisfying
\begin{equation}
d_2 = \min_{\bm a}\left\Vert \sum_{i=1}^{K(n-1)}a_{ni}\bm\varphi(\bm c_i) - \bm\varphi(\bm u(n))\right \Vert^2 \leq \delta
\end{equation}
where $\delta $, called the approximation level, is the threshold determining the sparisification level as well as system accuracy. If this condition is satisfied, the new sample in feature space could be linearly approximated by the centers in $\mathcal{C}(n-1)$. Therefore the effect of this sample to the mapping can be expressed through existing centers in center dictionary and there is no need to augment the center dictionary. Otherwise, the new sample whose feature vector is not approximately dependent on the samples should be added into dictionary and, consequently, a new coefficient corresponding to this center will be included. Through straight calculation it is easy to obtain,
\begin{equation}
\bm a(n) = \bm{\tilde K}(n-1)^{-1} \bm h(n)
\end{equation}
\begin{equation}
d_2 = \kappa(\bm u(n), \bm u(n)) - \bm h(n)^T \bm a(n)
\end{equation}
where $\bm h(n) = [\kappa(\bm c_1, \bm u(n)), \ldots,\kappa(\bm c_{K(n-1)}, \bm u(n))]^T$, the matrix $[\bm{\tilde K}(n-1)]_{ij} = \kappa(\bm c_i,\bm c_j) $. ALD not only is an effective approach to sparsification but also improves the overall stability of the algorithm because of its relation with the eigenvalue of $\bm{\tilde K}$ \cite{kernelbook}.
\subsubsection{Regularized KRLS-ALD}
Now they are two simply way to deal with KRLS-ALD: 1) Setting regularization term equal to 0 \cite{KRLS}. 2) Discarding the samples outside of the center dictionary and ignoring the influence of these samples \cite{kernelbook}. However, both of these two strategies have disadvantages. With the sparsification, the probability of overfitting decreases, even not happens when data number is small enough. This is the inspiration that method 1 used to set the regularization term equal to 0. Unfortunately, the success of this method depends on  the extent of sparsification. If the approximation level is not large enough, it doesn't mitigate overfitting. Owning to discarding useful information, the convergency speed and accuracy performance of the second method may be not satisfactory. In order to overcome these drawbacks, I proposed a general structure of regularized KRLS-ALD.

Define the matrices $\bm \Phi(n) = [\bm\varphi(\bm u(1)),\ldots,\bm\varphi(\bm u(n))]$, $\bm{\tilde{\Phi}}(n) = [\bm\varphi(\bm c_n),\ldots, \bm\varphi(\bm c_{K(n)})]$, according to \cite{KRLS},
\begin{equation}
\bm\Phi(n) = \bm{\tilde{\Phi}}(n) \bm A(n)^T + \bm \Phi(n)^{res}
\end{equation}
\begin{equation}
\bm \Omega(n+1) = \bm{\tilde\Phi}(n)\bm{\tilde \alpha}(n+1)
\end{equation}
where $\bm A(n) =[\bm a(1),...,\bm a(n)]^T$. Then, the cost function becomes
\begin{equation}
\label{Equ:costfunction}
\begin{split}
L(\bm{\tilde \alpha}(n+1)) &=\sum_{i=1}^n \|\bm \Omega(n)^T \bm \varphi(\bm u(i))-d(i)\|^2 + \lambda \| \bm \Omega(n)\|^2_\mathbb{F} \\
&= \|\bm \Phi(n)^T \bm{\tilde\Phi}(n) \bm{\tilde \alpha}(n+1) - \bm d(n)\|^2 + \lambda \|\bm \Omega(n) \|^2_\mathbb{F}  \\
& \approx \|\bm A(n) \bm {\tilde K}(n) \bm{\tilde \alpha}(n+1) - \bm d(n) \|^2 + \lambda\| \bm{\tilde\Phi}(n)\bm{\tilde \alpha}(n+1) \|^2_\mathbb{F} \\
\end{split}
\end{equation}
where $\bm d(n) = [d(1),...d(n)]^T$. In order to minimize $L(\bm{\tilde \alpha}(n))$, we take the derivative with respect to $\bm{\tilde \alpha}(n)$ and obtain
\begin{equation}
\label{Equ:derivative}
\frac{\partial L}{\partial \bm{\tilde \alpha}}(n+1) = 2(\bm A(n) \bm {\tilde K}(n))^T (\bm A(n) \bm {\tilde K}(n) \bm{\tilde \alpha}(n+1) - \bm d(n)) + 2\lambda  \bm{\tilde\Phi}(n)^T \bm{\tilde\Phi}(n)\bm{\tilde \alpha}(n+1)
\end{equation}
At the extremum the system  solution is:
\begin{equation}
\label{Equ:solution}
\bm{\tilde \alpha}(n+1)= [\bm A(n)^T \bm A(n) \bm{\tilde K}(n) + \lambda \bm I]^{-1}\bm A(n)^T \bm d(n)
\end{equation}
In the online scenario, at each time step, we are faced with either one of the following two cases.
\begin{enumerate}
  \item $ \bm \varphi(\bm x(n))$ could be approximated by $\mathcal C(n-1)$ according to ALD, that is $d_2 \leq \delta$. In this case, $\mathcal{C}(n) = \mathcal{C}(n-1)$.
  \item $d_2 > \delta$ and the new data $\bm \varphi(\bm x_n)$ is not ALD on $\mathcal{C}(n-1)$. Therefore, $\mathcal{C}(n) = \mathcal C(n-1) \bigcup \{\bm x(n)\}$.
\end{enumerate}
The key issue in this problem is how to design a iterative solution to obtain $\bm{\tilde \alpha}_n$. In the following, we denote $\bm P(n) = [\bm A(n)^T \bm A(n) \bm {\tilde K}(n)+ \lambda \bm I]^{-1}$ and derive the KRLS with regularization for each of these two cases.

\emph{Dictionary doesn't change:} In this case, $\bm {\tilde \Phi} (n)= \bm {\tilde \Phi}(n-1)$ and hence $\bm {\tilde K}(n) = \bm {\tilde K}(n-1)$. Only $\bm A$ changes between time steps: $\bm A(n)= [\bm A(n-1)^T, \bm a(n)]^T$. Therefore,
\begin{equation}
\bm A(n)^T \bm A(n) = \bm A(n-1)^T \bm A(n-1) + \bm a(n) \bm a(n)^T
\end{equation}
Such  matrix $\bm P(n)$ can be expressed as $\bm P(n) = [\bm P(n-1)^{-1} + \bm a(n) \bm a(n)^T \bm {\tilde K}(n)]^{-1}$. According to the matrix inversion lemma, assume $P(n-1) = A$, $\bm a(n) = B$, $\bm a(n)^T \bm {\tilde K}(n) = C$, $I=D$, yields,
\begin{equation}\label{Equ:PCenterUnchange}
\bm P(n)= \bm P(n-1)- \frac{\bm P(n-1)\bm a_n\cdot \bm a(n)^T \bm {\tilde K}(n) \cdot \bm P(n-1)}{1 + \bm a(n)^T \bm {\tilde K}(n) \cdot \bm  P(n-1)\bm a(n)}
\end{equation}
Defining $\bm q(n)=\bm P(n-1)\bm a_n/(1 + \bm a(n)^T \bm {\tilde K}(n) \cdot \bm  P(n-1)\bm a(n))$ and $\bm s(n)=\bm a(n)^T \bm {\tilde K}(n)$ ,The coefficient vector $\bm{\tilde \alpha}(n+1)$ could then be expressed as
\begin{equation}
\label{Equ:CoeffCenterUnchange}\begin{split}
\bm{\tilde \alpha}(n+1) & = \bm P(n) \bm A(n)^T \bm d(n) \\
& = \!\left[\!\bm P(n-1)\!-\! \frac{\bm P(n-1)\bm a(n)\cdot \bm a(n)^T \bm {\tilde K}(n) \cdot \bm P(n-1)}{1 + \bm a(n)^T \bm {\tilde K}(n) \cdot \bm P(n-1)\bm a(n)}\!\right]\!\! [\bm A(n-1)^T\bm d(n-1) \!+\! \bm a(n) d(n)] \\
& = \bm{\tilde \alpha}(n) + \frac{\bm P(n-1)\bm a(n)[d(n) - \bm a(n)^T \bm {\tilde K}(n) \bm{\tilde \alpha}(n)]}{1 + \bm a(n)^T \bm {\tilde K}(n) \cdot \bm P(n-1)\bm a(n)}\\
& = \bm{\tilde \alpha}(n)- \bm q(n)(d(n)-\bm s(n))\bm{\tilde \alpha}(n)\\
\end{split}
\end{equation}

\emph{The size of center dictionary increases:} In this condition, $\bm {\tilde \Phi}(n) = [\bm {\tilde \Phi}(n) , \bm \varphi(\bm u(n)) ]$. The matrix $\bm A$ changes to
\begin{equation}
\label{Equ:ACenterChange}
\bm A(n)  = \left[
  \begin{array}{cc}
    \bm A(n-1)  & \bm 0 \\
    \bm 0 & 1\\
  \end{array}
\right]
\end{equation}
Therefore,
\begin{equation}
\bm A(n) ^T \bm A(n)  = \left[
  \begin{array}{cc}
    \bm A(n-1) ^T \bm A(n-1) & \bm 0 \\
    \bm 0 & 1\\
  \end{array}
\right]
\end{equation}
\begin{equation}
\bm P(n)  = \left[
  \begin{array}{cc}
    \bm P(n-1) ^{-1} & \bm A(n-1) ^T \bm A(n-1) \bm h(n)  \\
    \bm h(n) ^T & \lambda + \kappa_{nn}\\
  \end{array}
  \right]^{-1}
\end{equation}
where $\bm h(n) = [\kappa(\bm c_1, \bm u(n) ), \ldots,\kappa(\bm c_{k_{n-1}}, \bm u(n))]^T$ and $\kappa_{nn}$ is the simplification of $\kappa(\bm u(n) , \bm u(n))$. Utilize the block matrix inversion identity, we obtain
\begin{equation}
\label{Equ:PCenterChange}
\bm P(n)  = \gamma(n) ^{-1} \left[
  \begin{array}{cc}
    \bm P(n-1) {-1}\gamma(n)  + \bm z_A(n)  \bm z(n) ^T & -\bm z_{A}(n)   \\
    -\bm z(n) ^T  & 1\\
  \end{array}
  \right]
\end{equation}
where $\gamma(n)  = \lambda + \kappa_{nn} - \bm h(n) ^T \bm z_A(n)  $, $\bm z_A(n)  = \bm P(n-1) \bm A(n-1) ^T \bm A(n)  \bm h(n) $, and $\bm z(n)  =\bm  P(n-1) ^T \bm h(n) $. Such the coefficient vector is updated as
\begin{equation}
\label{Equ:CoeffCenterChange}\begin{split}
\bm{\tilde \alpha}(n+1)  & = \bm P(n+1)  \bm A(n) ^T \bm d(n)  \\
& = \gamma(n) ^{-1} \left[
  \begin{array}{cc}
    \bm P(n-1) \gamma(n)  + \bm z_{A}(n)  \bm z(n) ^T & -\bm z_{A}(n)   \\
    -\bm z(n) ^T  & 1\\
  \end{array}
  \right] \left[
  \begin{array}{c}
  \bm A(n-1) ^T \bm d(n-1)  \\
  d_n \\
  \end{array}
  \right]   \\
  & =\left[\begin{array}{c}
  \bm{\tilde \alpha}(n) - \bm z_{A}(n)  \gamma(n) ^{-1} e(n)  \\
  \gamma(n) ^{-1}e(n)
  \end{array} \right]
\end{split}
\end{equation}
where $e(n) = d(n)  - \bm h_n^T \bm{\tilde \alpha }(n) $.

We now have obtained a recursive algorithm to solve the KRLS with regularization, which is referred to as  is described in pseudocode in Algorithm.\ref{algo_KRLS}.
\begin{algorithm}
\caption{Kernel RLS with regularization algorithm} \label{algo_KRLS}
\begin{algorithmic}
\STATE  \emph{\textbf{Initialization}:} Select the threshold $\delta > 0$ and the regularization parameter $\lambda > 0 $,
\STATE $\mathcal{C}(1) = \{\bm u(1)\}$, $P(1) = [\kappa_{11}+\lambda]^{-1}$, $\bm {\tilde K}(1) = \kappa_{11}$,
\STATE $\bm {\tilde K}(1)^{-1} = \kappa_{11}^{-1}$, $\bm A(1) = 1$, $\bm{\tilde \alpha}(1) = \bm P(1) y(1)$
\FOR {$n = 2,3...$}
\STATE Get the new sample: $(\bm u(n), d(n))$;
\STATE Compute $\bm h(n)$
\STATE ALD test: $\bm a(n) = \bm {\tilde K}(n-1)^{-1} \bm h(n)$
\STATE $d_2(n) = \kappa_{nn}-\bm h(n)^T \bm a(n)$
\IF { $d_2(n) \leq \delta^2$ }
\STATE  $\mathcal{C}(n) = \mathcal{C}(n-1)$, $\bm {\tilde K}(n) = \bm {\tilde K}(n-1)$, $\bm {\tilde K}(n)^{-1} = \bm {\tilde K}(n-1)^{-1}$
\STATE Compute $\bm{\tilde \alpha}(n+1)$  (Equ.\ref{Equ:CoeffCenterUnchange})
\STATE Update $\bm P(n)$ (Equ.\ref{Equ:PCenterUnchange})
\STATE Update $\bm A(n)$, $\bm A(n) = [\bm A(n-1)^T, \bm a(n)]^T$
\ELSE
\STATE $\mathcal{C}(n) = \mathcal{C}(n-1) \bigcup \{\bm u(n)\}$
\STATE Compute $\bm{\tilde \alpha}(n+1)$ (Equ.\ref{Equ:CoeffCenterChange})
\STATE Update $\bm P(n)$ (Equ.\ref{Equ:PCenterChange})
\STATE Update $\bm {\tilde K}(n)^{-1}$ and $\bm {\tilde K}(n)$
\STATE Update $\bm A(n)$ (Equ.\ref{Equ:ACenterChange})
\ENDIF
\ENDFOR
\STATE return $\bm{\tilde \alpha}(n+1)$, $\mathcal{C}(n)$
\end{algorithmic}
\end{algorithm}
Table \ref{tab:KRLSComplexity} summarizes the computational costs per iteration for KRLS-ALD with and without regularization.
\begin{table}[!t]
    \caption{\label{tab:KRLSComplexity}Computational costs per iteration}
    \centering
    {\small\begin{tabular} {cc}
    \hline
    \hline
    ALD test &  $O(K(n))^2$\\
    \hline
    Update $\bm P(n)$ & $O(K(n))^2$ \\
     \hline
    Update $\bm{\tilde \alpha}(n)$  & $O(K(n))$ \\
     \hline
    \end{tabular}}
\end{table}

\subsection{Kernel Least Mean Square Algorithm}
The KLMS utilizes the gradient descent techniques to search for the optimization solution of KRLS and is the least mean square algorithm (LMS) in RKHS. KLMS is obtained by minimizing the instantaneous cost function:
\begin{equation}\begin{split}
J(n)  & =\frac{1}{2}e(n)^2 \\
& = \frac{1}{2}\|\bm \Omega^T \bm \varphi(\bm u(n))-d(n)\|^2 \\
\end{split}\end{equation}
 Assume the initial condition of weight is $\bm \Omega(0)=0$, then use the LMS algorithm on RKHS, yields
\begin{equation}
\bm \Omega(n+1) = \bm \Omega(n)+\eta e(n)\bm \varphi(\bm u(n))
\end{equation}
However, the dimensionality of $\bm\varphi(.)$ is very high, so an alternative way is needed.
\begin{equation}\begin{split}
\bm \Omega(n+1) & = \bm \Omega(n)+ \eta e(n)\bm \varphi(\bm u(n)) \\
= & \bm \Omega(n-1) + \eta e(n-1)\bm \varphi(\bm u(n-1))+ \eta e(n)\bm \varphi(\bm u(n))\\
& \ldots \\
& = \eta \sum^{n}_{i=1}e(i)\bm \varphi(\bm u(i))\\
\end{split}\end{equation}
According to the ``kernel trick'' the output of system of the new input $\bm u(n)$ can be expressed as
\begin{equation}\begin{split}
\bm \Omega(n)^T \bm \varphi(\bm u(n)) & = [\eta \sum^{n-1}_{i=1}e(i)\bm \varphi(\bm u(i))]\bm \varphi(\bm u(n))\\
& =\eta \sum^{n-1}_{i=1}e(i)\kappa(\bm u(i),\bm u(n))\\
\end{split}\end{equation}
In conclusion, the learning rule of KLMS in the original space is as show in Algorithm. \ref{Kernel least mean square algorithm}.
\begin{algorithm}
\caption{Kernel least mean square algorithm}
\label{Kernel least mean square algorithm}
\begin{algorithmic}
 \STATE  \emph{\textbf{Initialization}:} Select the stepsize $\eta$,
 \STATE  $ e(1)= d(1) $;
 \STATE  $ y(1) = \eta e(1)$;
 \STATE  \emph{Computation}
 \WHILE {$\{\bm{u}(n), d(n)\}$ \emph{available}}
 \STATE $ y(n) =\eta\sum_{i=1}^{n-1}[e(i)\kappa(\bm{u}(i),\bm{u}(n))]$;
 \STATE $e(n) = d(n) - y(n)$;
 \ENDWHILE
\end{algorithmic}
\end{algorithm}
Among all of the kernel adaptive filters, KLMS is unique. It provides well-posedness solution with finite data \cite{klms} and naturally creates a growing radial-basis function (RBF) network \cite{kernelbook}. Moreover, as an online learning algorithm, KLMS is much simpler to implement, with respect to computational complexity and memory storage, than other batch-model kernel methods.
\subsection{Kernel Affine Projection Algorithm} The KLMS algorithm simply uses the current values to search the optimal solution. While the affine projection algorithm (APA) adopts better approximation by using the $K$ most recent samples and provides a trade-off between computation complexity and system performance. More interestingly, the KAPA provides a general framework for several exiting techniques, including KLMS, sliding window KRLS and kernel regularization networks \cite{kernelbook}.

\section{Discussion}
Even though kernel adaptive filters have advantages mentioned above and have been proved to be useful in complicated nonlinear regression and classification problems, they still have some drawbacks. For example, the conventional cost function, mean-square error (MSE) criterion, can not obtain the best performance in non-Gaussian situations. Moreover, the linear growing structure with each new sample leads to high computation burden and memory requirement. Therefore, computational complexity and memory increase linearly with time as more samples are processed, particulary for continuous scenarios which hinder so far online application, such as in DSP and FPGA implementations. Therefore to apply these powerful kernel adaptive filters in practice, we need to address two main issues first:
\begin{itemize}
  \item \textbf{How to obtain a robust performance in non-Gaussian situations?}
  \item \textbf{How to decrease the growing computational burden and memory requirement?}
\end{itemize}
These two issues are important aspects to judge system performance and are highly related. Improving robust performance normally results a the computational complexity increase. Therefore, appropriate processes, such as reducing the network size, to decrease computational complexity are necessary. On the other hand, approximation techniques to decrease computation complexity may suffer from system performance getting worse and could be compensated by improving system robust performance. Throughout the literature, similar problems have been studied from different perspectives leading to a myriad of techniques. Yet few techniques have been adopted to kernel adaptive filters. Our goal is to provide appropriate methods which take into account the intrinsic properties of kernel adaptive filters to solve these two problems.

There are various factors influencing kernel adaptive filters performance, such as cost function and kernel function. First, let us consider the cost function. MSE criterion is equivalent to maximum likelihood technique in linear and Gaussian condition  and obtains good performance, while it is not enough in other scenarios. The least mean p - power error (MPE) criterion, builds  a linear weighted combination of various powers of the prediction error, is investigated to solve this problem. However, many free parameters and prior knowledge requirements limit its wide application in practice. Recently, Information theoretic learning (ITL) has been proved more efficient to train adaptive systems. Different from conventional error cost function, ITL optimizes the information content of the prediction error to achieve the best performance in terms of information filtering. By taking into account the whole signal distribution, adaptive systems training through information theoretic criteria have better performance in various applications especially for which Gaussianity assumption is not valid. Furthermore, the nonparametric property in cost function and clear physical meaning also motivate us to introduce the ITL into the kernel adaptive filters to achieve a robust system and propose the kernel maximum correntropy algorithm (KMC).

Even though the accuracy of a system is improved, the redundant computational and memory burden is what we have to face with. Therefore the other important problem for the practical application of kernel adaptive filters is compacting network structure, which helps us decrease the total computational complexity. By including the important information and discarding relatively less useful data, growing and pruning techniques maintain the system network size into an acceptable range. Different from previous sparsification techniques, the quantized kernel least mean square utilized ``redundant'' data to update the network not purely discarded them to obtain a compact system with higher accuracy. As a complement of
spasification, pruning strategies makes the system structure more compact. Sensitive analysis based pruning strategies, discarding the units which make insignificant contribution to the overall network, are robust and widely utilized in various areas. Generally, omitting any information brings the accuracy down. Hence, the growing and pruning techniques is a trade off between system performance and compact structure.
Therefore, a measure of \emph{significance} is adopted in kernel adaptive filter to estimate the influence of process data and decide what information will be discarded.

The significance measure guides the system to fix the network size on a predefined threshold. However, this is not a real nonstationary learning. What we expected is that the network size should optimally be dictated by the complexity of the true system and the signal, while the system accuracy is acceptable. If we solve this problem, we have a truly online algorithm to handle real world nonstationary problems. This problem is transferred as how to obtain a trade-off between system complexity and accuracy performance. A variety of information criteria have been proposed so far to deal with this compromise problem. Among them, Minimal Description Length (MDL) has  great advantages of small computational costs and robust against noise, which have been proven in many applications particular in information learning area. MDL criterion has two formulations: batch model and online model.  Taking the approximation level selection in KRLS-ALD as an example,  the batch model MDL in kernel adaptive filters is illustrated firstly. Then we proposed an  KLMS sparsification algorithm  to explain the online model MDL. Owning to this proposed algorithm  separating the input (feature) space  with quantization techniques, it is called QKLMS-MDL.

This article mainly focus on the improvement of KLMS, the simplest of the kernel adaptive filters, but we believe, further extentions to other kernel adaptive filters, such as kernel recursive least square algorithm (KRLS) will broaden the the scope of the applications currently addressed in KLMS.

\bibliographystyle{unsrt}


\section{\label{}}

\section{INTRODUCTION}
The rigid body is an important idealization of moving and interacting parts and beings. The key feature is the fixed distance between two points in the body, since forces and torques cannot cause deformation of its shape~\citep{Stewart:2000}. The subject is relevant for the development of video-game physics engines since the simulation of motions and interactions of rigid bodies approximates reality fairly well~\citep{Garstenauer:2006}. Prosthetic limbs and robotic devices are other important application areas of this theme.

Here we offer a tool to stimulate engineering students to continue their learning path in the subject. We show mathematical formulations and numerical simulations of the interdependent motions.

\section{MATHEMATICAL FORMULATIONS}

The analysis has its origins in~\cite{Shames:1996}, where there is a solved example for the instantaneous position as shown in~Fig.~\ref{fig:subfig1}, and general guidelines for a computer project are proposed. In what follows, equations for the kinematics are derived showing angle and length constraints, as well as velocity and acceleration relations while the disk undergoes either a constant $2~rad/s$, or a variable $\left( 0.10~+~2.3\cdot\,\alpha \right)\;rad/s$, with $\alpha$ being the turn angle of the disk~[Fig.~\ref{fig:subfig2}].

\begin{figure*}[t]
\centering
\subfigure[Initial position]{
    \includegraphics[width=0.4\textwidth]{mec1.png}
    \label{fig:subfig1}
}
\subfigure[Intermediate position]{
    \includegraphics[width=0.4\textwidth]{mec2.png}
    \label{fig:subfig2}
}
\caption[Optional caption for list of figures]{Articulated mechanism at its \subref{fig:subfig1} initial and \subref{fig:subfig2} intermediate positions. Reproduced from~\cite{Shames:1996}.}
\label{fig:both}
\end{figure*}

\subsection{Angles and length constraints}

The three angles describing rotation of the articulated mechanism are $\alpha$, $\beta$, and $\gamma$. A one letter notation is introduced to write trigonometric functions as:

\begin{align}
      l&=sin\, \alpha&\ \ \ \ m&=cos\, \alpha \\
      n&=sin\, \beta&\ \ \ \ o&=cos\, \beta \\
      p&=sin\, \gamma&\ \ \ \ q&=cos\, \gamma
\end{align}

Since the length of the linking bar AB cannot change, the horizontal displacement of pin A must be equal to the horizontal displacement of pin B, or:

\begin{equation}
\overline{OA} \cdot\; l\; =\; \overline{BC}\; -\; \overline{BC} \cdot\; o
\label{eq1}
\end{equation}

An overbar means the length of a given straight line segment. It can be inferred from Fig.~\ref{fig:subfig1} that $\overline{OA}\,=\,\overline{BC}$. Simplifying:

\begin{equation}
l\; =\; 1\; -\; o
\label{eq2}
\end{equation}

A right triangle with $\overline{AB}$ as the hypothenuse will have, at its intermediate position in Fig.~\ref{fig:subfig2}, the following opposite side to angle $\gamma$:

\begin{equation}
\overline{BC} \cdot\; n\; +\; \overline{OA} \cdot\; m
\label{eq3}
\end{equation}

And the adjacent side:

\begin{equation}
\overline{OC}\; -\; \overline{OA} \cdot\; l\; -\; \overline{BC} \cdot\; o
\label{eq4}
\end{equation}

Applying the Pitagorean theorem in Fig.~\ref{fig:subfig1}:

\begin{equation}
\overline{AB}\; =\; \sqrt{\left(\overline{OC}\; -\; \overline{BC}\right)^{2}\; +\; \left(\overline{OA}\right)^{2}}
\label{eq5}
\end{equation}

These relations are needed to obtain $\beta$, $\gamma$, and the equations for velocities and accelerations in the following sections.

\subsection{Velocities}

With pin A viewed as belonging to the disk, its velocity is:

\begin{equation}
\vec{V}_A =  \vec{\omega} \times \vec{\rho}_{OA}
\label{eq6}
\end{equation}

Where $\vec{\omega}\ =\ 2\; \hat{k}\ rad/s$ for the constant rotation case, or $\vec{\omega}\ =\ \left(~0.10~+~2.3\cdot\,\alpha \right) \hat{k}\ rad/s$ for the variable rotation one. For $\alpha$ going from $0^o$ to $90^o$ counterclockwise starting at the lowest position, $\omega$ varies from $0.10\ rad/s$ to $3.7\ rad/s$ when $\omega$ varies linearly with $\alpha$. This accounts for a $1.9\ rad/s$ average angular velocity. We can thus write:

\begin{equation}
  \alpha(t)\ =\ 1.9\cdot\,t
\label{eq6a}
\end{equation}

Taking the time derivative of $\omega$ we get:

\begin{equation}
\dot{\omega}\ =\ 2.3\cdot\,\dot{\alpha}
\label{eq6b}
\end{equation}

And

\begin{equation}
  \dot{\omega}\ =\ \left(2.3 \right)\,\left(1.9 \right)\ =\ 4.4\,rad/s^2
\label{eq6c}
\end{equation}

With pin B belonging to the output bar BC, its velocity is:

\begin{equation}
\vec{V}_B = \vec{\omega}_{BC} \times \vec{\rho}_{CB}
\label{eq7}
\end{equation}

Where $\vec{\omega}_{BC}$ is the still unknown angular velocity vector of bar BC. The motion of linking bar AB can be described using Chasles theorem~\citep{Shames:1996}:

\begin{equation}
\vec{V}_A = \vec{V}_B + \vec{\omega}_{AB} \times \vec{\rho}_{BA}
\label{eq8}
\end{equation}

Where $\vec{\omega}_{AB}$ is the still to be determined angular velocity vector. $\vec{V}_A$ and $\vec{V}_B$ are always tangential to the circular trajectories of pins A and B. For the mechanism's intermediate position shown in Fig.~\ref{fig:subfig2}, $x$ and $y$ components for both velocities can be written as:

\begin{align}
      V_{Ax}&=m \cdot V_A&\ \ \ \ V_{Ay}&=l \cdot V_A \label{eq8a} \\
      V_{Bx}&=n \cdot V_B&\ \ \ \ V_{By}&=o \cdot V_B \label{eq8b}
\end{align}

The fixed vectors in the rigid bodies going from O to A $\left(\vec{\rho}_{OA}\right)$, C to B $\left(\vec{\rho}_{CB}\right)$, and B to A $\left(\vec{\rho}_{BA}\right)$ can also be expressed using their $x$ and $y$ components. An inspection of Fig.~\ref{fig:subfig2} allows one to write:

\begin{align}
      \rho_{OAx}&= + l \cdot \overline{OA} &\ \ \ \ \rho_{OAy}&= - m \cdot \overline{OA} \label{eq8c} \\
      \rho_{CBx}&= - o \cdot \overline{BC} &\ \ \ \ \rho_{CBy}&= + n \cdot \overline{BC} \label{eq8d} \\
      \rho_{BAx}&= - q \cdot \overline{AB} &\ \ \ \ \rho_{BAy}&= - p \cdot \overline{AB} \label{eq8e}
\end{align}

Substituting (\ref{eq8a}) -- (\ref{eq8e}) back in Eqs.~(\ref{eq6}),~(\ref{eq7}), and~(\ref{eq8}):

\begin{equation}
V_B = V_A \cdot \frac{m \cdot q + l \cdot p}{n \cdot q + o \cdot p}
\label{eq9}
\end{equation}

and

\begin{equation}
\omega_{AB} = \frac{V_A}{\overline{AB}}\cdot \left(\frac{m \cdot o - l \cdot n}{o \cdot p + n \cdot q}\right)
\label{eq10}
\end{equation}

With~(\ref{eq9}) in Eqs.~(\ref{eq7}) and~(\ref{eq8b}), $\omega_{BC}$ can be calculated.

\subsection{Accelerations}

We now need to distinguish between the \emph{constant} and the \emph{variable} $\omega$ scenarios.
Pin A only has radial acceleration with constant $\omega$, but it has both radial and tangential components otherwise. Pin B, as before, will have both radial and tangential components.

\subsubsection{Constant $\omega$}
Taking the derivative of Eq.~(\ref{eq8}) with respect to time:

\begin{equation}
\vec{a}_A = \vec{a}_{BR} + \vec{a}_{BT} + \dot{\vec{\omega}}_{AB} \times \vec{\rho}_{BA} + \vec{\omega}_{AB} \times \left( \vec{\omega}_{AB} \times \vec{\rho}_{BA} \right)
\label{eq11}
\end{equation}

The linear accelerations $\vec{a}_A$, $\vec{a}_{BR}$ and $\vec{a}_{BT}$ can be expressed in terms of unit vectors $\hat{i}$ and $\hat{j}$:

\begin{equation}
\vec{a}_A = \left( \omega^{2}\, \overline{OA} \right) \cdot \left(-l\, \hat{i} + m\, \hat{j} \right)
\label{eq12}
\end{equation}

\begin{equation}
\vec{a}_{BR} = \left( {\omega_{BC}}^{2}\, \overline{BC} \right) \cdot \left(o\, \hat{i} - n\, \hat{j} \right)
\label{eq13}
\end{equation}

\begin{equation}
\vec{a}_{BT} = - \left( \dot{\omega}_{BC} \overline{BC} \right) \cdot \left(n\, \hat{i} + o\, \hat{j} \right)
\label{eq14}
\end{equation}

Substituting back into Eq.~(\ref{eq11}) and working the vector algebra, one is able to obtain expressions for angular accelerations of bars AB and BC. As these expressions are long, they are grouped into four different terms as $AB_1$ through $AB_4$ for $\dot{\omega}_{AB}$, and $BC_1$ through $BC_4$ for $\dot{\omega}_{BC}$. It should be noted that these quantities point along the z-direction. The final result is:

\begin{align}
      AB_{1}&=\omega^2\;\overline{OA}\cdot \left( l o + m n \right)&\ \ \ \ AB_{2}&={\omega_{BC}}^2\;\overline{BC}\cdot \left(o^2 + n^2 \right) \label{eq15a}\\
      AB_{3}&={\omega_{AB}}^2\;\overline{AB}\cdot \left( o q - n p \right)&\ \ \ \ AB_{4}&=\overline{AB}\;\cdot \left(o p + n q \right) \label{eq15b}
\end{align}

\begin{equation}
\dot{\omega}_{AB} = - \frac{AB_1 + AB_2 + AB_3}{AB_4}
\label{eq15}
\end{equation}

\begin{align}
      BC_{1}&=\omega^2\;\overline{OA}\cdot \left( m p - l q \right)&\ \ \ \ BC_{2}&={\omega_{AB}}^2\ \overline{AB} \cdot \left(p^2 + q^2 \right) \\
      BC_{3}&={\omega_{BC}}^2\;\overline{BC}\cdot \left( n p - o q  \right)&\ \ \ \ BC_{4}&=\overline{BC}\;\cdot \left(n q + o p \right)
\end{align}

\begin{equation}
\dot{\omega}_{BC} = - \frac{BC_1 - BC_2 + BC_3}{BC_4}
\label{eq16}
\end{equation}

\subsubsection{Variable $\omega$}

Since the disk now has a variable rate of rotation $\dot{\omega}$, pin A's acceleration is:

\begin{equation}
\vec{a}_A = \vec{a}_{AR}~+~\vec{a}_{AT}
\label{eq17}
\end{equation}

Referring back to Fig.~\ref{fig:subfig2}, the tangential component can be written as:

\begin{equation}
  \vec{a}_{AT}~=~\dot{\omega}\,\overline{OA}\,\left( m \hat{i}\,+\,l \hat{j} \right)
\end{equation}

The radial component is still given by Eq.~(\ref{eq12}). As for the constant case, angular accelerations for bars AB and BC are written as:

\begin{align}
      AB_{1}&=\omega^2\;\overline{OA}\cdot \left( l o + m n \right)&\ \ \ \ AB_{2}&={\dot{\omega}}\;\overline{OA}\cdot \left(m o - l n \right) \label{eq18a}\\
      AB_{3}&={\omega_{BC}}^2\;\overline{BC}\cdot \left(o^2 + n^2 \right)&\ \ \ \ AB_{4}&={\omega_{AB}}^2\;\overline{AB}\cdot \left( o q - n p \right) \label{eq18b}\\
      AB_{5}&=\overline{AB}\;\cdot \left(o p + n q \right)&\ \ \ \ \label{eq18c}\\
\end{align}

\begin{equation}
\dot{\omega}_{AB} = \frac{-AB_1 + AB_2 - AB_3 - AB_4}{AB_5}
\label{eq18}
\end{equation}

And,

\begin{align}
      BC_{1}&=\omega^2\;\overline{OA}\cdot \left( m p - l q \right)&\ \ \ \ BC_{2}&={\dot{\omega}}^2\ \overline{OA} \cdot \left(m q + l p \right) \\
      BC_{3}&={\omega_{BC}}^2 \;\overline{BC}\cdot \left( n p - o q  \right)&\ \ \ \ BC_{4}&={\omega_{AB}}^2\ \overline{AB} \cdot \left(p^2 + q^2 \right) \\
      BC_{5}&=\overline{BC}\cdot \left( n q + o p  \right)&\ \ \ \
\end{align}

\begin{equation}
\dot{\omega}_{BC} = \frac{-BC_1 - BC_2 - BC_3 + BC_4}{BC_5}
\label{eq19}
\end{equation}

This concludes the kinematical formulation of the articulated mechanism that has been implemented numerically and placed under \emph{git} version control~\citep{Fleischfresser:2014}.

\section{NUMERICAL EXPERIMENTS}

\subsection{Constant $\omega$}
Geometrical and kinematics parameters for the $90^{o}$ turn of the disk are shown in Fig.~\ref{fig4}. Of note is the change of rotation sense for the linking bar AB from counterclockwise (positive $\omega_{AB}$) to clockwise (negative $\omega_{AB}$) when $\alpha = 30^{o}$ (top right panel). This feature is also evident when $\gamma$ reaches a maximum (top left panel), and $V_B$ and $\omega_{BC}$ remain nearly constant around $\alpha = 30^{o}$ (bottom left and top right panels). Angular accelerations, despite their intricate mathematical expressions [Eqs.~(\ref{eq15a}) through~(\ref{eq16})], are slowly varying functions after the start-up period.

\begin{figure*}[t]
\centering
\includegraphics[angle=0, scale=0.45]{todos.png}
\caption{\textit{Top left}: $\beta$ and $\gamma$ vs. $\alpha$. \textit{Bottom left}: linear velocities of pins A and B vs. $\alpha$. \textit{Top right}: angular velocities of bars AB and BC vs. $\alpha$. \textit{Bottom right}: angular accelerations of bars AB and BC vs. $\alpha$. Angular step size $\Delta \alpha\,=\,0.5^{o}$.}
\label{fig4}
\end{figure*}

\begin{table}[!ht]
\begin{center}
\caption{Head and tail of linear and angular velocities. See text for comments on cells with asterisks.}
\begin{tabular}{|c|c|c|c|c|c|c|c|}
\hline
$\alpha$ ($^{o}$) & $\beta$ ($^{o}$) & $\gamma$ ($^{o}$) & $V_A$ ($m/s$) & $V_B$ ($m/s$) & $\omega$ ($rad/s$) & $\omega_{AB}$ ($rad/s$) & $\omega_{BC}$ ($rad/s$) \\
\hline
 0 & 0 & 16.7 & 0.6 &  2* &  2* &  2* & -6.7\\
 \hline
 \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\
\hline
 90 & 90 & 16.7 &  0.6* & 0.18 & 2 & -0.6* & -0.6* \\
\hline
\end{tabular}
\label{tab1}
\end{center}
\end{table}

\pagebreak

Table~\ref{tab2} has the computed angular accelerations' initial and final values. The results for $\alpha\,=\,0^{o}$ serve to validate the numerical implementation since the computations match the textbook calculations in~\cite{Shames:1996}. The same is true for the angular velocities of Tab.~\ref{tab1}.

\begin{table}[!ht]
\begin{center}
\caption{Head and tail of angular accelerations.}
\begin{tabular}{|c|c|c|c|c|}
\hline
$\alpha$ ($^{o}$) & $\beta$ ($^{o}$) & $\gamma$ ($^{o}$) & $\dot{\omega}_{AB}$ ($rad/s^2$) & $\dot{\omega}_{BC}$ ($rad/s^2$)\\
\hline
 0 & 0 & 16.7 & -57.8 & 192.6\\
 \hline
 \vdots & \vdots & \vdots & \vdots & \vdots \\
\hline
 90 & 90 & 16.7 & 0.0 & 5.2 \\
\hline
\end{tabular}
\label{tab2}
\end{center}
\end{table}

\subsection{Variable $\omega$}

Results are now displayed in Fig.~\ref{fig5}. Note the linear increase of $V_A$ (bottom left panel), and how angular speeds and accelerations evolve without resemblance to the constant $\omega$ case. In particular, it is now easy to see when both bars start decelerating during the $90^{o}$ counterclockwise turn.

\begin{figure*}[h]
\centering
\includegraphics[angle=0, scale=0.45]{variavel.png}
\caption{\textit{Top left}: $\beta$ and $\gamma$ vs. $\alpha$. \textit{Bottom left}: linear velocities of pins A and B vs. $\alpha$. \textit{Top right}: angular velocities of bars AB and BC vs. $\alpha$. \textit{Bottom right}: angular accelerations of bars AB and BC vs. $\alpha$. Angular step size $\Delta \alpha\,=\,0.5^{o}$.}
\label{fig5}
\end{figure*}

\begin{table}[!ht]
\begin{center}
\caption{Head and tail of linear and angular velocities. See text for comments on cells with asterisks.}
\begin{tabular}{|c|c|c|c|c|c|c|c|}
\hline
$\alpha$ ($^{o}$) & $\beta$ ($^{o}$) & $\gamma$ ($^{o}$) & $V_A$ ($m/s$) & $V_B$ ($m/s$) & $\omega$ ($rad/s$) & $\omega_{AB}$ ($rad/s$) & $\omega_{BC}$ ($rad/s$) \\
\hline
 0 & 0 & 16.7 & 0.03 &  0.1* &  0.1* &  0.1* & -0.33\\
 \hline
 \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\
\hline
 90 & 90 & 16.7 &  2.1* & 0.63 & 7.0 & -2.1* & -2.1* \\
\hline
\end{tabular}
\label{tab3}
\end{center}
\end{table}

\begin{table}[!ht]
\begin{center}
\caption{Head and tail of angular accelerations.}
\begin{tabular}{|c|c|c|c|c|}
\hline
$\alpha$ ($^{o}$) & $\beta$ ($^{o}$) & $\gamma$ ($^{o}$) & $\dot{\omega}_{AB}$ ($rad/s^2$) & $\dot{\omega}_{BC}$ ($rad/s^2$)\\
\hline
 0 & 0 & 16.7 & 4.2 & -14.1\\
 \hline
 \vdots & \vdots & \vdots & \vdots & \vdots \\
\hline
 90 & 90 & 16.7 & -1.3 & 61.7 \\
\hline
\end{tabular}
\label{tab4}
\end{center}
\end{table}

\section{DISCUSSION AND CONCLUSIONS}
The relations needed to simulate the motion of an articulated mechanism were presented and the kinematics was analyzed numerically. Results for initial and final positions ($\alpha = 0^{o}$ and $90^{o}$ respectively) are shown in Tab.~\ref{tab1}, Tab.~\ref{tab2}, Tab.~\ref{tab3} and Tab.~\ref{tab4}. At the initial position ($\alpha = 0^{o}$) both $V_B$ and $\omega_{AB}$ have the same numerical values of the disk's rotation rate (cells with asterisks). Paper and pencil solutions give instantaneous values while numerical experiments capture the motion's evolution. Hence, if only a paper and pencil solution is attempted for $\alpha = 0^{o}$, one may be tempted to believe these parameters always have the same numerical values. The same would be true if one attempts a paper and pencil solution for the final position ($\alpha = 90^{o}$), since $\omega_{AB}$ and $\omega_{BC}$ have the same absolute values of the constant linear speed of a point at the edge of the disk (cells with asterisks).

The experiments also demonstrated the inversion of rotation for the linking bar under plane motion for two different scenarios: one where the rate of rotation of the disk is constant, and one where it is made to vary linearly with the angle of rotation. Analysis such as the one presented here complement routine paper and pencil solutions which alone may lead to erroneous conclusions about the motions. A design modification of the original mechanism is shown in Fig.~\ref{fig6}. The distance between the fixed centers of rotation is shorter than the original mechanism to allow full disk revolution while keeping the linking bar size unchanged, and the disk rotation is now clockwise. The full revolution is not possible in the original configuration (Fig.~\ref{fig:both}).

Despite the complexities of the rigid bodies' formulations described in this article, one should keep in mind that this is still a relatively simple mechanics problem. All principles applied are exact and the numerical implementations only required algebraic expressions. Complications arise if one attempts to formulate and model collisions between rigid bodies. Impulse-momentum principles are needed compared to the limited kinematics toolbox. Moreover, collision laws are not well understood and numerical implementations need to handle differential equations with discontinuities that may lead to energy conservation inconsistencies and other difficulties~(\citep{Stewart:2000},~\citep{Garstenauer:2006},~\citep{Chatterjee:1997}).

\begin{figure*}[h]
\centering
\includegraphics[angle=0, scale=0.7]{mecanismoArticulado.png}
\caption{Articulated mechanism design modification that allows for full disk revolution.}
\label{fig6}
\end{figure*}

\newpage