\section{Introduction}
Intrinsically localized modes, or discrete solitons (breathers), appear as
generic excitations in a large variety of physical systems \cite{rep2}, where
spatial periodicity (or discreteness) provides gaps in the linear dispersion
relation and nonlinearity allows for detuning the oscillation frequencies
into
these gaps, resulting in spatial localization. Particularly important
applications
of discrete nonlinear systems  are nonlinear optical waveguide arrays
\cite{rep1},  with possibilities to change and control all essential
parameters, such as geometry, dimensionality, nonlinearity, beam angle, etc.
For weakly coupled waveguides, the relevant mathematical description, derived
via coupled mode theory \cite{cj88,rep1}, is a Hamiltonian lattice equation
of the discrete nonlinear
Schr\"odinger (DNLS) type \cite{rep2,rep1,cj88,ej03}, where the particular
type of nonlinearity
depends on the nonlinear response of the medium.  The most studied case is a
DNLS equation with cubic on-site nonlinearity\cite{ej03}, corresponding to a
Kerr medium \cite{rep1}, which also appears generically as a modulational
equation for the small-amplitude dynamics in chains of coupled anharmonic
oscillators \cite{rep2}.   However, photorefractive media also enables the
generation of discrete spatial solitons \cite{rep1}, and the corresponding
lattice model is a DNLS
equation with saturable nonlinearity (s-DNLS) which was studied in a number
of theoretical works
\cite{prlkip,jpakhare,plaeil,prlmel,ob07,pre73,rnkf09}.

A particularly interesting property of the one-dimensional (1D) s-DNLS model
\cite{prlkip,jpakhare,plaeil,prlmel,ob07}
is the existence of certain ``sliding velocities'', where  localized  discrete
solitons may travel in the lattice without radiation. This behavior was
connected to the existence of ``transparent points'' associated with the vanishing
of a so-called Peierls-Nabarro (PN) potential barrier \cite{pre48}, usually
defined as the difference in energy (Hamiltonian) at constant power (norm)
between the two fundamental localized stationary solutions centered at one
site (odd mode) and symmetrically in between two sites (even mode). Close to
the points of vanishing of this energy difference are regions of stability
exchange between the even and odd solutions, associated in general with the
appearance of a family
of intermediate, asymmetric stationary solutions, connecting both types of
symmetric solutions at the bifurcation points \cite{pdaub}. The existence of
regimes of enhanced mobility close to such bifurcation points is also
well-known from studies of other one-dimensional lattice models
\cite{pdaub,pdcret,oje03}.

On the other hand, it is still an open issue whether moving discrete solitons
may exist as localized, radiationless modes also in two-dimensional (2D)
lattices. As was shown numerically in \cite{pre73}, a scenario with exchange
of stability through bifurcations with asymmetric stationary solutions
appears
also for the 2D saturable model in a square (isotropic) lattice, involving in
this case three different types of fundamental solutions \cite{krb00}:
one-site (odd-odd, OO), two-site (odd-even, OE), and four-site
(even-even, EE) modes. It was also shown numerically
in \cite{pre73}, that solutions with good (but generally not radiationless)
mobility in the axial directions may exist in these regimes, and that the
necessary energy needed for rendering a mobile stable stationary solution
agreed well with the concept of PN-barrier, if its definition was
extended to take into account also the energy for the relevant intermediate
stationary solution
(an analogous situation is well-known for the PN potential of
kinks \cite{pr82}). Very similar observations were also made later for a
2D DNLS model with cubic-quintic nonlinearity  \cite{ccmk09}. Moreover, in
the
low-power regime of the 2D s-DNLS equations, good mobility was also observed
in diagonal directions \cite{pre73}.

However, the relation between existence of regions of stability exchange,
small PN-barrier and mobile localized solutions in 2D is not that trivial, as
shown, e.g., for a model with cubic inter-site nonlinearities in \cite{oj09}:
even in regimes with small PN barrier and existence of intermediate
solutions,
the mobility may be very poor, if there is no continuous path in phase space
passing close to the relevant stationary modes. On the other hand, as was
observed for the low-power
(i.e., close to continuum limit) regime of a 2D lattice with quadratic
nonlinearity in \cite{susanto07}, the effective Peierls-Nabarro potential may
in some situations be weak enough to allow mobility in {\em arbitrary}
directions, without any direct connection to bifurcations and symmetry-broken
stationary solutions.

So there is clearly need for a better understanding of the conditions for
mobility in 2D lattices. It is the purpose of the present paper to generalize
the concept of PN-barrier as discussed above, and introduce a full 2D PN
potential surface describing the  pseudo-potential landscape in-between all
stationary modes. We will use a numerical constrained Newton-Raphson (NR)
method,
previously applied to 1D lattices in \cite{rodrigo1d}, to explicitly
construct
these surfaces for the
2D saturable model from \cite{pre73}, and show how parameter regimes and
directions of good
mobility may be immediately identified from smooth, flat parts of these
surfaces. We will also illustrate how the interplay between
translational motion in one lattice direction and oscillatory motion in the
orthogonal direction can be intuitively understood from the topology of the
corresponding PN surfaces.

The structure of the paper is as follows. In Sec.\ \ref{model}, we describe
the
2D s-DNLS model, its stationary localized solutions (Sec.\ \ref{stat}),
and the explicit implementation of the numerical constraint method used to
calculate the PN potential surfaces (Sec.\ \ref{constraint}).
We point out some essential differences between the 2D implementation and
earlier implementations for the 1D case. In Sec.\ \ref{surfaces} we identify
the regimes of parameter space where smooth surfaces can be found for the
isotropic lattice, and make a classification of appearing surfaces with
different topologies. In Sec.\ \ref{dynamics} we illustrate with direct
dynamical simulations various types of behavior for discrete solitons moving
slowly in different directions, corresponding to the surfaces discussed in
Sec.\ \ref{surfaces}.
We also here show examples on how simultaneous excitation of transverse
oscillation modes may affect
discrete solitons moving slowly in axial directions.
In Sec.\ \ref{anisotropy}, we analyze effects of adding a weak anisotropy to
the lattices in the parameter regimes used in previous sections, in order
mainly to investigate whether this may promote mobility in directions
different from the axial or diagonal ones (which presumably would be
associated with the appearance of additional "valleys" in the PN potential
surface). Finally, in Sec.\ \ref{conclusions} we make some concluding
remarks,
summarizing our results and pointing out some relevant directions for future
research.

\section{Model}\label{model}

We consider the following general form of the 2D s-DNLS equation in a
rectangular geometry~\cite{pre73}:
\begin{equation}
i\frac{\partial U_{n,m}}{\partial z}+\Delta_{n,m}^\alpha U_{n,m}
-\frac{\gamma U_{n,m}}{1+|U_{n,m}|^2} = 0\ ,
\label{sdnls}
\end{equation}
where $z$ corresponds to the normalized propagation coordinate,
$\gamma$ to an effective nonlinear parameter, and $U_{n,m}$ represents the
(complex) electric field amplitude at site $\{n,m\}$ in an array of
$N\times M$ sites.
$\Delta_{n,m}^\alpha U_{n,m}\equiv (U_{n+1,m}+U_{n-1,m})+\alpha(U_{n,m+1}+U_{n,m-1})$
defines the linear interaction between nearest-neighbor waveguides,
where $\alpha$
is an anisotropy parameter allowing us to study different 2D scenarios
from a completely isotropic lattice ($\alpha=1$) to a completely decoupled
set
of 1D arrays ($\alpha=0$). Model (\ref{sdnls}) possesses two conserved
quantities, the power (norm) defined as
\begin{equation}
P=\sum_{nm} |U_{n,m}|^2 \ ,
\label{pow}
\end{equation}
and the Hamiltonian (energy) defined as
\begin{eqnarray}
H=-\sum_{nm} \left[ \left(U_{n+1,m}+\alpha U_{n,m+1}\right)U_{n,m}^*
\right.\nonumber
\\ \left. -\frac{\gamma}{2}\ln (1+|U_{n,m}|^2)+\text{c.c.}\right] .
\label{ham}
\end{eqnarray}
\subsection{Stationary solutions}
\label{stat}

Stationary solutions to Eq.\ (\ref{sdnls}) are of the form
$U_{n,m}(z)=u_{n,m}\exp{(i \lambda z)}$, where $u_{n,m}$ is a $z$-independent
(generally complex) amplitude and $\lambda$ corresponds to the propagation
constant or
frequency \cite{pre73}.
Extended stationary solutions with constant amplitude (plane waves)
exist in frequency bands, whose edges depend on the amplitude $|u_{n,m}|$.
The linear band region is obtained for low-amplitude plane waves
($|u_{n,m}| \rightarrow 0$), and
it is easy to show that they will exist in the region
$\lambda \in [-\gamma-2(1+\alpha),-\gamma+2(1+\alpha)]$.
Moreover, for high-amplitude plane waves ($|u_{n,m}| \rightarrow \infty$),
the nonlinear term saturates and the corresponding interval is
$[-2(1+\alpha),2(1+\alpha)]$.

\begin{figure}[htbp]
\centering
\includegraphics[width=0.48\textwidth]{perfilesc.pdf}
\caption{(Color online) Examples of spatial profiles for (a) one-site (OO), (b) two-site horizontal (EO), (c) two-site vertical (OE), (d) four-site solutions (EE), (e) IS1, (f) IS2, (g) IS3-vertical, and (h) diagonal solutions, respectively, illustrated for an isotropic lattice ($\alpha = 1$).}
\label{profs}
\end{figure}

The fundamental localized stationary solutions to Eq.\ (\ref{sdnls}),
defined as solutions with real amplitudes $u_{n,m}$ having a single
maximum distributed on the sites in one unit-cell of the lattice, were
described for the isotropic case in \cite{pre73}. In general, we may
identify one-site (OO), two-site horizontal (EO), two-site vertical (OE),
and four-site (EE) solutions.
Figs.\ \ref{profs}(a)-(d) show typical profiles for these types of excitations.
As discussed in \cite{pre73},
these in-phase localized solutions bifurcate from the fundamental linear mode
(upper band edge) in the small-amplitude limit, and merge into the
corresponding high-amplitude plane-wave mode in the infinite-power limit.
As a consequence, they exist in the region
$\lambda\in [-\gamma+2(1+\alpha),2(1+\alpha)]$, so that the size of the
existence region will be just $\gamma$.
[A two-site diagonal solution \cite{Kalosakas} [see Fig.\ \ref{profs}(h)] could also be excited but only
in frequency regimes where solutions are very localized with
$\lambda \lesssim 0$, requiring large values of $\gamma $
(for $\alpha=1$, $\gamma \gtrsim 8.2$)].
In addition, in certain parameter regimes there are also asymmetric,
intermediate stationary solutions, associated with exchange of stability
between the symmetric ones \cite{pre73}.
We may identify three types of such solutions, termed henceforth
intermediate 1 (IS1), intermediate 2 (IS2),
and intermediate 3 (IS3), respectively. IS1 and IS3 both connect
a stable two-site mode, OE/EO, with another simultaneously stable solution:
the one-site mode OO (IS1), or the four-site mode EE (IS3) [see Figs.\ \ref{profs}(e) and (g)].
Thus, in these cases the unstable intermediate solutions act as carriers of instability between the
corresponding fundamental modes. The unstable IS2 solution exists when the two-site solutions are stable. It connects the unstable one-site solution OO with the likewise unstable four-site solution EE, and stabilizes the latter mode when reaching it [see Fig.\ \ref{profs}(f)].

\begin{figure}[htbp]
\centering
\includegraphics[width=0.46\textwidth]{PLc.pdf}\vspace{0.1cm}
\includegraphics[width=0.46\textwidth]{phb.pdf}\vspace{0.05cm}
\includegraphics[width=0.46\textwidth]{pgc.pdf}
\caption{(Color online) Properties of fundamental stationary solutions for
$\gamma=4$ and $\alpha=1$:
(a) Power versus frequency; (a)-inset: $\Delta P/P$ versus $\lambda$; (b) $\Delta H$ versus power; (c) Stability versus power.
Different solutions are plotted with different line-types:
one-site (thick solid blue), two-site (thick dotted blue),
four-site (thick dashed blue), IS1 (thin solid black),
IS2 (thin dotted black), and IS3 (thin dashed black), respectively.
}\label{PL}
\end{figure}

Figs.~\ref{PL} shows fundamental properties of different stationary solutions
for an isotropic lattice with $\gamma = 4$ (similar curves were plotted and
discussed in
\cite{pre73}, but we here include some additional examples to facilitate
the comparison with the full PN surfaces in Sec.\ \ref{surfaces}). Power versus frequency curve is shown in Fig.~\ref{PL} (a), where different solutions repeatedly cross each other in
the whole range of parameters. In Fig.~\ref{PL} (a)-inset we plot $\Delta P/P\equiv (P_i-P_{OO})/P_{OO}$ as a function of $\lambda$, for $i=OO, EO/OE, EE, IS1, IS2, IS3$. This figure shows clearer how solutions cross each other including many intermediate solutions appearing in the region $\lambda\sim 1.3$.
Moreover, we also observe these crossings when plotting $\Delta H\equiv H_i-H_{OO}$ versus power [Fig.~\ref{PL} (b)]. (We plot $\Delta H$ instead of $H$ because Hamiltonian differences between stationary solutions of the same power are generally quite small, which is indeed a
favorable scenario for moving solutions).
Inset in Fig.~\ref{PL} (b) shows a zoom in the region $P\sim 9.8$ where the
Hamiltonian of different fundamental solutions matches at some points.
This was initially interpreted as a vanishing PN barrier~\cite{prlkip},
which is actually not the case due to the intermediate solutions appearing
in-between, raising the effective energy barrier~\cite{pre73} [see Fig.~\ref{PL} (b)-inset when thick blue lines coincide]. However, good
mobility would be expected close to these regions but, still, it will be
strongly determined by the specific kick or perturbation given to the
solution
in order to put it in movement.

In order to study the stability of stationary localized solutions we
implement
a standard method \cite{jpakhare,pre73}, obtaining all the stable/unstable
linear
perturbation-modes and their corresponding instability measure
$g=\max_j[\Im \{\omega_j\}]$,
where $\omega_j$ are the oscillation frequencies
of the linear eigenmodes.  Thus, if $g=0$
the solution is (at least marginally) stable and if $g>0$ it is unstable.
Stability versus power for the same solutions as in Figs.\ \ref{PL} (a), (b)
is
shown in Fig.\ \ref{PL} (c). (Note that in the corresponding figures 2-3 of
\cite{pre73}, the quantity $G=-g^2$ was used as instability measure.)
From this figure, we can see how unstable
intermediate solutions IS1, IS2, IS3
appear when two or three solutions
(regarding OE and EO as different solutions)
are simultaneously stable.
Note that we never observed regions for simultaneously
stable one and four-site solutions for  isotropic coupling. However,
as will be discussed further in Sec.\ \ref{anisotropy}, in the anisotropic
case
such regions exist.

\subsection{Constraint method}
\label{constraint}

The constraint method allows us to construct energy surfaces connecting
stationary solutions for a given value of power. In that sense, it helps us
to effectively predict and interpret
the dynamics across the lattice. Critical points will
represent stationary solutions, and a coherent
movement across the lattice should transform one solution to the other by
keeping the power constant~\cite{pre48}.
The method was originally introduced by Aubry and
Cretegny~\cite{pdaub} and then implemented by Savin et al.~\cite{savin} to study the mobility of kinks in nonlinear Klein-Gordon lattices (see also Ref.\ \cite{Sepulchre}
for a related approach to analyze travelling
breathers in 1D oscillator chains in terms of an effective Hamiltonian.)
Lately, it was numerically
implemented to analyze surface states
in one-dimensional semi-infinite systems~\cite{rodrigo1d}.
By adiabatically changing the amplitude in
one
particular site, specifically chosen as the one after the main peak
(the peak is at $n_c$ and the constrained amplitude at $n_c+1$),
the one- and the two-site solutions could be connected. Ending the sweep
when $u_{n_c}=u_{n_c+1}$ and the center of mass of the constrained solution,
$X$, has varied from $n_c$ to $n_{c}+0.5$,
a one-dimensional energy surface, $H$ vs $X$, can be sketched.
Technically speaking, the method used in Ref.\ \cite{rodrigo1d} consists on
eliminating one equation from the Newton-Raphson problem, the one of the
constrained amplitude which is not anymore an unknown variable. However, as
the power is kept constant, an equation for $P$ is added and the
frequency $\lambda$ becomes a variable completing the variable-equations set.

In the present work, for the 2D lattice we implement
a more sophisticated method
where we explicitly vary the center of mass instead of the amplitude.
Since we are interested in the energy landscape around the fundamental
stationary solutions, we will assume also the amplitude  $u_{n,m}$ of the
constrained solutions to be real
and positive on the constraint sites.
Then, from the definition of the center of mass coordinates
in the horizontal and vertical directions,
\begin{equation}
X\equiv\frac{\sum_{nm} n |u_{n,m}|^2}{P}\ \ \ \ \ \
\text{and}
\ \ \ \ \ \ Y\equiv\frac{\sum_{nm} m |u_{n,m}|^2}{P}\ ,
\label{xy}
\end{equation}
we can solve for any pair of amplitudes $A \equiv  u_{n_A,m_A} $ and
$B\equiv  u_{n_B,m_B } $ at places $\{n_A,m_A\}$ and $\{n_B,m_B\}$:
\begin{eqnarray}
A&=&\sqrt{\frac{m_B S_n-n_B S_m +P(n_B Y-m_BX)}{n_Bm_A-n_Am_B}}\nonumber\\
&&\label{cmab}\\
B&=&\sqrt{\frac{n_A S_m -m_A S_n +P(m_A X-n_A Y)}{n_Bm_A-n_Am_B}}\nonumber
\label{cma}
\end{eqnarray}
where $S_j\equiv\sum_{nm}(1-\delta_{n,n_A}\delta_{m,m_A}-
\delta_{n,n_B}\delta_{m,m_B}) j  |u_{n,m}|^2$ with $j=n,m$.
Thus, as before the actual
constraints will be in the amplitudes $A$ and $B$, but now we can tune the
center of mass
as wished, from a given
stationary solution towards any other.
In general, the constrained solutions obtained in this way will of course
not be stationary solutions of the full system. To
identify any stationary solution - including the intermediate ones -  we
check whether the value of $\lambda$ obtained from the constrained
NR scheme coincides with the frequency of a hypothetical
true stationary solution to
the full  Eq.\ (\ref{sdnls}) with the computed amplitude profile.
Furthermore, constrained solutions allow us
to calculate their Hamiltonian, and therefore to construct an effective
energy landscape.

\begin{figure}[htbp]
\centering
\includegraphics[width=0.45\textwidth]{scheme.pdf}
\caption{(Color online)
Array scheme showing the constraint locations used to obtain
the energy surfaces below. }
\label{scheme}
\end{figure}
The choice for positions $\{n_A,m_A\}$ and $\{n_B,m_B\}$ is evidently not
unique, but in order to most efficiently trace out a smooth energy landscape
connecting the fundamental solutions (if it exists), the constraint
sites should preferably be chosen within the unit cell where the
amplitudes are large.
It turned out that, starting from a stationary
one-site solution centered at
$\{n_c,m_c\}$, in most cases
the best option is to choose the first constraint at
site $\{n_c+1,m_c\}$ or $\{n_{c},m_c+1\}$, and the other one at site
$\{n_c+1,m_c+1\}$, as sketched in Fig.\ \ref{scheme}. The results shown
in the following sections are obtained using these constraint sites. We
also tried using constraint sites at $\{n_c+1,m_c\}$  and $\{n_{c},m_c+1\}$
(i.e., along a diagonal); however, with this choice
we typically were not able to find the four-site EE solution when starting
from the one-site OO. Instead, the NR method with constraints on these sites
generally converges,
at $X=n_c+\frac{1}{2}, Y=m_c+\frac{1}{2}$, to a two-site diagonal solution
(which as remarked in Sec.\ \ref{stat} is an exact solution to the
unconstrained equation only for large $\gamma$), which is not expected to
be relevant for the mobility process.

\section{Isotropic lattice}

\subsection{Description of PN surfaces}\label{surfaces}

\begin{figure*}[h]
\includegraphics[width=0.45\textwidth]{g4a1p5.pdf}
\includegraphics[width=0.45\textwidth]{g4a1p945.pdf}
\includegraphics[width=0.45\textwidth]{g4a1p10.pdf}
\includegraphics[width=0.45\textwidth]{g4a1p12.pdf}
\includegraphics[width=0.45\textwidth]{g4a1p355.pdf}
\caption{(Color online)
Energy surfaces for $\gamma=4$, $\alpha=1$, in the five different
power regimes discussed in the text : (a) $P=5$; (b) $P=9.45$; (c) $P=10$;
(d) $P=12$; (e) $P=35.5$. The center of mass $\{X,Y\}$ for the four stationary
solutions are:
$\{8,8\}$ (OO), $\{8.5,8\}$ (EO), $\{8,8.5\}$ (OE), and $\{8.5,8.5\}$ (EE).
White dots denote local extrema corresponding to intermediate solutions.
(System size $N=M=15$ and fixed boundary conditions were used.)
}
\label{pot1}
\end{figure*}

By exploring the parameter-space $(\gamma, P)$, we have identified regions
where we were able to
compute complete two-dimensional energy surfaces (using constraint as shown in
Fig.\ \ref{scheme}) and regions where we cannot. In general,
complete surfaces could not be obtained for large values of
the nonlinearity constant, $\gamma \gtrsim 6$, when solutions get strongly
localized. In some cases we could trace some specific mobility directions
(e.g., along lattice axes or  diagonal) connecting stationary solutions
with almost equal Hamiltonian, but not
the full two-dimensional scenario involving all fundamental solutions.
Note also that, as remarked in \cite{pre73}, for large $\gamma$ and small
power the s-DNLS model essentially behaves as a cubic DNLS model with
effective nonlinearity $\gamma P$.
Consequently, we could not find complete
energy surfaces with all fundamental solutions for the cubic DNLS either
(although surfaces involving the two-site diagonal solution in place of the
EE solution could be obtained as mentioned in Sec.\ \ref{constraint};
for the cubic DNLS this two-site
solution does exist as a true stationary solution above some threshold
level of power \cite{Kalosakas}).

On the other hand, for very small values of $\gamma$, the computation of
energy surfaces becomes difficult for technical reasons: due to the widening
of the solutions, considerably larger lattices are needed to remove the
influence from boundary effects. Therefore, in
the following, we will present the main phenomenology for the isotropic
case ($\alpha=1$) found for intermediate values of $\gamma$,
$3 \lesssim \gamma \lesssim 5$, where complete
surfaces involving all four fundamental solutions
are obtained for all values of the power. We will show
results for the particular value $\gamma=4$, but
the scenario is found to be qualitatively the same for all $\gamma$ in this
interval.

From Fig.\ \ref{PL} (c), we may identify essentially
five different regimes where energy surfaces of qualitatively different
nature
should be expected, depending on the level of power.
The first
one is for low power, where (similarly to the cubic DNLS)
the one-site solution is always stable
and the other fundamental stationary solutions are all
unstable. The
corresponding energy surface is illustrated in  Fig.\ \ref{pot1} (a), where
the one-site solution yields the energy minimum, the two-site solutions
saddle points and the four-site solution the maximum. Note that, in this
low-power regime, the surfaces for this value of $\gamma$  are
still rather flat, and therefore some mobility may result if the one-site
solution is kicked to overcome the barriers, in the axial as well as in
the diagonal directions, as illustrated by Figs.\ 4 (a) and (c)
in Ref.\ \cite{pre73}. (Another example of mobility in this regime is discussed
in Sec.\ \ref{dynamics} below.)

The saturable nature of the system becomes evident at higher powers. In the
second regime, appearing for the first time when
$9.27 \lesssim P \lesssim 9.55$,  the one- and the two-site
solutions are stable simultaneously
and, as shown in
Fig.\ \ref{pot1} (b), these three points all correspond to local
minima of the surface. Intermediate solutions IS1 connecting the one- with the
two-site solutions in the
horizontal and vertical directions appear
as saddle points [white dots in
Fig.\ \ref{pot1} (b)].
Note that the energy landscape for the parameter values
in Fig.\ \ref{pot1} (b) is almost flat between the one-site and
two-site solutions in the axial directions,
leading to the very good axial mobility shown  in Fig.\ 4 (d)
of Ref.\ \cite{pre73}, while the maximum corresponding to the
unstable four-site solution creates a
too large effective barrier to overcome in the diagonal direction.

The third  power region is the one in which only the two-site solutions
are stable.
It  appears for the first time when $9.55 \lesssim P \lesssim 10.04$,
and is illustrated
in  Fig.\ \ref{pot1} (c). Here, the stable two-site solutions
correspond to two local minima of the surface, and the unstable one- and
four-site solutions both to local maxima. The two unstable solutions are
connected by the intermediate solution IS2 (white dot at the surface),
corresponding to a saddle point which for symmetry reasons
(for the isotropic lattice) will lie
along the diagonal connecting the unstable solutions.
As will be illustrated in
Sec.\ \ref{dynamics}, the easiest mobility  in this case is expected to occur
in a diagonal direction, connecting the two stable stationary solutions.

The fourth regime corresponds to simultaneously stable two- and four-site
solutions.  Such a regime appears for the first
time when
$10.04 \lesssim P \lesssim 10.32$.
If we further
increase the power, we will find a similar region with these three
simultaneously stable solutions for
$34.5\lesssim P \lesssim 36.25$. In this regime, as illustrated in
Fig.\ \ref{pot1} (e),
the two- and four-site solutions all correspond to minima,
the unstable
intermediate solutions IS3 to saddle points,  and
the one-site solution to a maximum.
Thus, comparing Figs.\ \ref{pot1} (b) and (e),
we see that the structure of the potential has been completely inverted.
Now, a big ``hill'' is located at the one-site
solution, and as a consequence no
mobility is expected involving this solution. As will
be shown below in  Sec.\ \ref{dynamics}, the simplest mobility
scenario will now be the one in which the two-site solutions travel through
the
four-site one across the lattice.

The fifth regime, with the four-site solution being the only stable one,
occurs for the first time when
$10.32 \lesssim P \lesssim 34.5$.
As illustrated in Fig.\ \ref{pot1} (d),  the four-site solution
is now a  minimum of the PN potential,
while the other three unstable
solutions correspond to saddle points (EO and OE) and maximum (OO),
respectively. (There are no intermediate solutions in this regime.)
Thus, by
increasing the power
we have now
reverted the
surface compared to the low-power regime in Fig.\ \ref{pot1} (a).

Further increasing of power shows for $36.25 \lesssim P \lesssim 43.6$ a
new region of stable two-site solutions, corresponding to the third regime.
For $43.6 \lesssim P \lesssim44.2$  the one- and two-site solutions are
simultaneously stable, so the scenario is equivalent to the second regime,
followed again for  $P\gtrsim 44.2$ by the first regime,
respectively. This second
complete inversion thus corresponds to a regain of the low-power
characteristics  of the surface.
Furthermore, repetition of these scenarios can be found for $P\gtrsim 118$,
and we confirmed, e.g., the existence of a complete, smooth surface, analogous
to the one shown in Fig.\ \ref{pot1} (d) for the fifth regime, for $P=120$.

\subsection{Mobility dynamics for isotropic lattice}\label{dynamics}

To explicitly show the connection between the different types of
energy surfaces described in Sec.\ \ref{surfaces} and the mobility of
localized solutions,
we numerically integrate model\ (\ref{sdnls}) by taking as initial condition a
stationary solution perturbed with a small
kick: $U_{n,m}(0)=u_{n,m} \exp[i k_x(n-n_c)+i k_y(m-m_c)]$,
where $k_x$ and $k_y$ correspond to the kick strength
in the horizontal and vertical
directions, respectively.
If the surfaces would be completely flat,
solutions should move even with an infinitesimally small kick.
However, as we have seen, in general surfaces are not flat due to
the discreteness and the self-induced PN potential, and although the PN
barrier can be very small in certain directions, it is generally non-zero.
Therefore, in order to put a localized solution in movement,
some amount of kinetic energy (represented by this kick) should be given to
effectively overcome the energy barriers.

\begin{figure*}[htbp]
\includegraphics[width=0.45\textwidth]{mobp5sc.pdf}
\includegraphics[width=0.45\textwidth]{mobp10.pdf}
\includegraphics[width=0.45\textwidth]{mobp12.pdf}
\includegraphics[width=0.45\textwidth]{mobp355.pdf}
\caption{(Color online)
Examples of mobility dynamics in the propagation direction $z$,
corresponding to surfaces shown in Fig.\ \ref{pot1} (a), (c), (d), and (e),
respectively. In each subfigure, top figures show profiles movement [(a) shows a colormap where colors were normalized to the maximum amplitude of each plot; this colormap also applies to (b)-(d)],
and bottom figures center of mass evolution for $X$ (full line) and $Y$
(dashed line). (a) $P=5$, $k_x=k_y=0.038$; (b) $P=10$, $-k_x=k_y=0.018$;
(c) $P=12$, $k_x=k_y=0.015$; (d) $P=35.5$, $k_x=0.015$ and $k_y=0$. (Other
parameter values are the same as in Fig.\ \ref{pot1}.)}
\label{dina1}
\end{figure*}
We first discuss an example from the first, low-power, regime
corresponding to the energy surface in Fig.~\ref{pot1} (a)
(as mentioned above, different mobility examples from this regime were
shown in Figs.\ 4 (a) and (c) in Ref.\ \cite{pre73}).
Here, we took the (unstable) EO solution and kicked it with a very small value
in
both directions, yielding the dynamics
numerically observed in Fig.~\ref{dina1} (a).
In the  $y$-direction, the kinetic energy is not sufficient to
overcome the barrier created by
the four-site solution, while in the $x$-direction
it can move towards the minimum corresponding to the one-site solution.
The resulting dynamics for the center-of-mass positions shows
how the movement in the horizontal direction gets combined with
oscillations in the vertical one.
The dependence $X$ vs $z$ also clearly shows how the
solution feels the potential in terms of its velocity: maximum velocities
occur around integer lattices sites (where the potential is a minimum) while
the
minimum velocities occur close to middle points
[see full line in Fig.~\ref{dina1} (a); compare also with the analogous
scenarios for the kicked OO mode in  Figs.\ 4 (a) and (c) of
Ref.\ \cite{pre73}].

In the second power regime, with energy surfaces as in Fig.~\ref{pot1} (b),
good mobility can be expected only in axial directions, with effective energy
barrier determined by the intermediate solution between the one- and the
two-site
solutions as discussed in Ref.\ \cite{pre73}. The dynamics
in this power regime, with a very small kick applied only in the axial
direction to one of the stable stationary solutions,
will be as already illustrated in Fig.\ 4 (d) of Ref.\ \cite{pre73}:
the solution moves very slowly and adiabatically traces the shape of the
potential with a minimal velocity at the places of the intermediate solutions.

In the third power regime, with surfaces as in Fig.~\ref{pot1} (c),
a very interesting kind of mobility is observed: a
diagonal
mobility between the (stable) horizontal and  vertical two-site solutions
as illustrated in  Fig.~\ref{dina1} (b). The initial
EO solution, kicked equally
in the $(-x)$- and $y$-directions,
gets sufficient kinetic energy to pass over the small barrier created
by the intermediate IS2 solution. It then continues through the OE solution,
passes another IS2 barrier, and then to the other EO
solution shifted by one site in both directions.
Although
the potential connecting these two solutions is not completely flat,
there is a very good
transport of energy in this direction, allowing mobility for more than one
lattice diagonal in the considered $15 \times 15$ lattice.

If we take a look at the surface in the fifth power regime
[Fig.~\ref{pot1} (d)], we
realize that an initial (unstable)
one-site solution may move in any direction by
slightly kicking it since it corresponds to a maximum.
Fig.~\ref{dina1} (c) shows an example for kicking the
one-site solution symmetrically in both directions in order to make it move
passing through the four-site, i.e., a diagonal movement.
The velocity has maximum in the minima
of the potential (corresponding to the EE solution) and minimum in the
potential maxima (OO solution) (no intermediate solutions appear in this
regime).

Finally, in the fourth power regime,
we see from the surface [Fig.~\ref{pot1} (e)] that the
two- and four-site solutions are stable simultaneously, presenting an
intermediate solution in-between that will define the effective energy barrier.
This barrier is very small and, therefore, a very small kick is also needed.
Figure~\ref{dina1} (d) shows the evolution starting from a two-site vertical
solution, passing through the intermediate one and arriving to the four-site
solution.

As illustrated by the example in Fig.\ \ref{dina1}(a), the
energy surfaces also provide nice intuitive interpretations to the dynamics
of moving discrete solitons with additional perturbations
{\em transverse} to the direction of motion.
In this
example, the dynamics in the two orthogonal directions
appear essentially independent of each other
(propagation in $X$ while oscillating in $Y$).
However, as we will illustrate
with another example below, there are situations where the particular
surface topologies may lead to a more intricate interplay between the
oscillatory and translational dynamics.

We present in
Fig.\ \ref{uns} a case with $P=44$.
With the notation from Sec.\ \ref{surfaces}, this value
of the power belongs to the second power regime, and
the structure of the surface is phenomenologically
identical to the one
sketched in Fig.\ \ref{pot1}(b), but for a higher level of power.
A picture of this energy landscape, periodically extended along
the $X$-direction,
is shown in Fig.\ \ref{uns} (a).
In this figure, the darkest regions
correspond to the positions for the minima (the OO and EO/OE modes),
while the brightest regions correspond to
maxima (EE modes). Note also the positions of the saddle points
(IS1, marked with arrows) in-between the OO and EO/OE modes.

\begin{figure}[htbp]
\centering
\hspace{0.58cm}
\includegraphics[width=0.47\textwidth]{unstaa2.pdf}
\includegraphics[width=0.47\textwidth]{unstab2.pdf}
\includegraphics[width=0.47\textwidth]{unstac2.pdf}
\includegraphics[width=0.47\textwidth]{unstad2.pdf}
\caption{(Color online)
Dynamics of an OO mode with $P=44$ ($\gamma=4, \alpha=1$),
kicked with $k_y=0.001$ and varying $k_x$.
(a) Extended energy surface. Darker (brighter) regions
correspond to lower (higher) values of the energy $H$.
(b) $Y_9,\ Y_{10},\ Y_{11},\ Y_{12}$ versus $k_x$ represented by (black)
filled circles,
(grey) squares, (blue) diamonds and (red) triangles, respectively.
(c) $Y(z)$ versus $X(z)$ for
$k_x=0.0111,\ 0.0141,\ 0.021,\ 0.049$ and $0.083$ represented by
thick (blue) solid,
thin (black) dashed, (blue) dotted, thin (black) solid and
thick (blue) dashed lines, respectively.
(d) Dashed (blue) and solid (black) lines correspond to $k_x=0.018$ and $k_x=0.035$,
respectively, for $N=M=41$.}
\label{uns}
\end{figure}

The kick we choose for the vertical direction is very small,
$k_y=0.001$, in order to disturb the horizontal movement only weakly.
Since the potential shape in the $Y$-direction - in a first approximation - will
correspond to a harmonic potential, the profile will essentially
perform harmonic oscillations in the $Y$-direction, with a frequency (determined by the curvature of the surface)
that will vary only slightly as the mode
translates in the $X$-direction.
To look for possible resonances between the oscillatory and translational
dynamics,
we implement the following scheme:
For fixed $k_x$ and
$k_y$ we numerically integrate model (\ref{sdnls}), starting with the OO
stationary solution centered at position $\{8,8\}$, and measure the z-values for
which $X$ is equal to $9, 10, 11$ and $12$, respectively.
For these z-values we compute
$Y_m\equiv Y(X=m)$, and plot in Fig.\ \ref{uns}(b) as a function of the
horizontal kick strength $k_x$.

The particular values of $k_x$ where symbols in Fig.\ \ref{uns}(b)
coincide for $Y=8$ correspond
to kick strengths where the natural frequency of the surface potential
is commensurate with the translational velocity.
In Fig.\ \ref{uns}(c) we show
some examples for these particular points.
For $k_x=0.0111$ the solution makes two cycles
before arriving to the next integer position in $X$. For $k_x=0.0141$ the
solution makes one and a half cycles, for $k_x=0.021$ one cycle, for
$k_x=0.049$ one half of cycle and for $k_x=0.083$ one quarter of cycle, before
arriving to the next integer position in $X$. As can be seen, around these
points the oscillations apparently remain bounded and the dynamics stable,
at least for long times. (Note that we do not see any visible phase-locking
effects between the translational and oscillatory motion here: the symbols in
Fig.\ \ref{uns}(b) seem to coincide only at isolated points and not in
intervals. This may be an indication that the dissipative effects on the
effective dynamics are very weak.)

However, in other regimes we observe that the amplitude of the oscillations
in the $Y$-direction rapidly grows. Two examples are illustrated in
Fig.\ \ref{uns}(d), corresponding to values of  $k_x$
where the $Y$-displacement in Fig.\ \ref{uns}(b) has maxima.
(Here, we have used a larger system  of $41\times 41$ sites,
in order to clearly identify the nature of this particular dynamics.)
The initial increase of the $Y$-amplitude suggests a kind
of unstable dynamics. However, as is clear from the longer simulation
in  Fig.\ \ref{uns}(d), the oscillation amplitude in the Y-direction remains
bounded and thus the trajectory does not escape in the vertical direction but
remains trapped by the barrier.
The typical dynamics, thus, exhibits amplitude-modulated oscillations
(beatings), with the largest amplitudes corresponding to maxima close to
integer $X$ where the surfaces are flatter.
Thus, this kind of mobility is a signature of the 2D topology of the
potential and it also validates the surfaces computed with our method.
(We performed several computations in different $\{P,k_x\}$-regimes finding
similar results as the ones presented in Fig.\ \ref{uns}).

\section{Effects of weak anisotropy}\label{anisotropy}

\begin{figure}[htbp]
\centering
\includegraphics[width=0.48\textwidth]{sta.pdf}
\caption{(Color online)
Density plots of $g$ versus $P$ and $\alpha$ for $\gamma =4$, for
the (a) OO; (b) EO; (c) OE; and (d) EE solutions, respectively.
Black means $g=0$ and lighter colors imply an increasing instability.}
\label{stabi}
\end{figure}
Considering weak anisotropy by putting
$\alpha < 1$ in Eq.\ (\ref{sdnls})
(implying a stronger coupling in the horizontal than in the vertical
direction),
the main qualitative modifications as concerns the
stationary solutions are shifts of the
stability regions for the different solution types (evidently, the
two-site horizontal, EO, and vertical, OE, solutions are now non-equivalent).
Focusing our discussion on the parameter value  $\gamma =4$ analyzed
in detail for the isotropic case in the previous sections,
the results for the stability analysis are presented in Fig.\ \ref{stabi}.
We find that, for fixed $\alpha < 1$ and increasing power,
the first stability inversion appears between the one-site OO
and the horizontal two-site EO solutions, and it occurs for lower values
of the power than in the isotropic case.
The first interval of stability of the four-site EE
solution gets narrower when $\alpha$ decreases, mainly because its lower limit
increases (with
a corresponding increase of the upper stability limit for the EO solution).
The first
stability region of the vertical two-site OE solution disappears for
$\alpha \lesssim 0.96$. On the other hand, its second stability regime for
higher powers
is enhanced, whereas the corresponding stability regime for its horizontal
counterpart gets narrower and disappears for $\alpha \lesssim 0.89$.

\begin{figure}
\includegraphics[width=0.45\textwidth]{g4a72p195.pdf}
\includegraphics[width=0.45\textwidth]{g4a07p20.pdf}
\caption{(Color online)
Energy surfaces for anisotropic lattices with $\gamma=4$:
(a) $\alpha=0.72$, $P=19.5$; (b) $\alpha=0.7$, $P=20$.}
\label{anpot1}
\end{figure}
With stronger anisotropy, some qualitatively new regimes not present in the
isotropic case will appear. For $\alpha \lesssim 0.76$, a new region of
stability of the OO mode appears around  $P\approx 20$,
leading to a small window of multistability where the
OO, EO and EE modes are stable simultaneously. An example of an energy
surface from this regime
is shown in Fig.~\ref{anpot1} (a),
where also two intermediate solutions [white dots] can be
found on the
edges. There is also a regime with  bistability of
the OO and EE modes only, as illustrated by
the pseudo-potential surface in Fig. \ref{anpot1} (b).
Note that, in this case, the intermediate
solution does not leave the edge, in contrast to
the IS2 mode existing on the diagonal between
the one- and four- site solutions for the isotropic case of simultaneously
stable EO/OE solution [Fig. \ref{pot1} (c)].
Since anisotropy destroys the symmetry between the OE and EO solutions,
bistability of both 2-site solutions
disappears already for a small amount of anisotropy, and for smaller
$\alpha$ the IS position can only
be found on an edge of the potential surface.

Thus, since no additional ``valleys'' were found to develop in the PN surfaces
(all minima appear on the edges), the best mobility for anisotropic lattices
should generally be expected
to occur along lattice directions (at least in cases where diagonal
coupling terms can be neglected, as assumed throughout this work).

\section{Conclusions}\label{conclusions}

In conclusion, we have studied in a deeper way the problem of
mobility of localized modes
in two-dimensional saturable discrete systems.
We numerically implemented
a constrained Newton-Raphson method
to construct full Peierls-Nabarro energy
surfaces, which appeared as very useful tools for predicting the dynamical
properties of localized excitations. Although these surfaces were never found
to be completely flat
(and therefore the corresponding Peierls-Nabarro barriers is strictly never zero),
parameter regimes and directions of good mobility were seen to correspond
to smooth, flat parts of the surfaces.

For the isotropic saturable model, five different surface topologies
could be identified in different power regimes,
depending on the
stability properties of the different fundamental stationary solitons.
By numerically studying the dynamics of perturbed stationary
solutions, we showed how these different topologies yielded qualitatively
different kinds of optimal mobility,
generally in axial or diagonal directions and with velocities varying
according to the shape of the potential
while the profile is propagating across the lattice.
The energy surfaces were also found to be very useful for interpreting the
dynamics resulting from the interplay between translational and oscillatory
motion in orthogonal directions.

We also studied the effect
of weak anisotropy, where the breaking of the symmetry of the system
also yields energy surfaces with lower symmetry. In addition, we
showed how this symmetry breaking could lead to a
change of stability properties for the fundamental
stationary solutions, and as a result two new types of surface topologies
were identified. From the shape of these surfaces, we predicted that the
best mobility in the anisotropic case should generally occur along
lattice directions.

Thus, the method developed here appears to be a very good tool for the
understanding of the mobility dynamics of localized
excitations in higher-dimensional lattices,
clearly showing the ``real'' energy barriers that localized solutions
experience.
We hope that this approach
will be useful also for other types of nonlinar lattice models.

Finally, from an application point of view, we have given further examples
on how the saturable nature of nonlinearity promotes the
generation of many different stability and mobility scenarios. The example
presented in section \ref{dynamics} shows the possibility to steer the
dynamics of the solution by using
the natural frequencies of its own self-induced potential. As an application
of this, we propose a discriminative optical switch which, for
example, will only be activated when $X$ and $Y$ positions are -
simultaneously -
an integer number. This is one more example of nonlinear lattices being the
key to success for controlling the propagation of light in future all-optical
technologies.

\begin{acknowledgments}

We are grateful to Serge Aubry and Sergej Flach who both, independently, suggested this approach to analyze 2D mobility, and also for discussions and comments on the manuscript. We thank Mario I. Molina for useful discussions. The authors acknowledge financial support from different sources: FONDECYT 1070897 and 7080001, FB0824/2008, and CONICYT fellowship. Authors acknowledge the hospitality of the Max Planck Institute for the Physics of Complex Systems (MPIPKS), Dresden, where this work was in part developed. M.J. also acknowledges support from the Swedish Research Council.
\end{acknowledgments}


\section{Introduction}

Time series are sequences of values measured over time. Examples include
financial data, such as stock prices,  or medical data, such as blood sugar
levels. Classifying time series is an important class of problems
which is applicable to music classification~\citep{springerlink:10.1007/s11634-007-0016-x},
medical diagnostic~\citep{Sternickel2002109} or bioinformatics~\citep{Legrand2008215}.

Nearest Neighbor (NN) methods classify time series efficiently and accurately~\citep{1454226}.
The 1-NN method is especially simple: we merely have to find the nearest labeled instance.

We need to specify a distance: the Euclidean and Dynamic Time Warping~\citep{1163055} distances are popular choices.
However, we can also learn a distance based on some training data~\citep{yang2006distance,weinberger2009distance}.
Given the training data set made of classes of time series instances,
we can either learn a single (global) distance function, or learn one distance
function per class~\citep{csatari2010,paredes2000class,1221358}. That is, to compute the distance between a test element and an instance of class~$j$, we use
a distance function specific to class~$j$.

Because the Euclidean distance is popular for NN classification, it is tempting to consider
generalized ellipsoid distances~\citep{ishikawa1998mindreader}, that is, distances of the form
\begin{eqnarray*}
D(x,y) = (x-y)^{\top} M (x-y).
\end{eqnarray*}
When $M$ is a positive semi-definite matrix, the square root of this distance is a pseudometric:
it is symmetric, non-negative and it satisfies the triangle inequality ($\sqrt{D(x,y)}+\sqrt{D(y,z)}\leq \sqrt{D(x,z)}$).
Further, when $M$ is a positive definite matrix, then the square root becomes a
metric because $D(x,x)=0 \Rightarrow x=0$.
When the matrix $M$ is the identity matrix, we recover the (squared) Euclidean distance.
We get the  Mahalanobis distance when solving for the matrix $M$ minimizing the sum of distances
$\sum_{x,y}D(x,y)$ (see~\S~\ref{sec:mathstuff}).
We can require $M$ to be diagonal, thus defining the diagonal Mahalanobis distance~\citep{1221358,ishikawa1998mindreader}.
In the more general case, solving
for $M$ can be difficult, as it often involves inverting a low-rank matrix.
Perhaps partly due to these mathematical difficulties, there has been no attempt to use the general
Mahalanobis distance to classify time series (to our knowledge). Thus, for the first time,
we apply the full-matrix Mahalanobis distance for time series classification. To solve
the mathematical difficulties, we use both an approach based on a pseudoinverse
and on covariance shrinkage. With both the diagonal
and the covariance shrinkage approaches, the
square root of the distance $D$ is a metric (under mild assumptions) whereas our approach based on a  pseudoinverse merely
generates a pseudo-metric.

While we find that  the full-matrix Mahalanobis distance is not competitive when
relying on a pseudoinverse, we get good results with covariance shrinkage
or using the diagonal Mahalanobis distance.
Moreover, we find that the class-based Mahalanobis distance
is preferable to the global Mahalanobis distance.

\section{Related Works}

Several distance functions are used for time series classification, such as
\begin{itemize}
\item Dynamic Time Warping (DTW)~\citep{ratanamahatana2004making},
\item DISSIM~\citep{frentzos2007index},
\item  Threshold
Queries~\citep{abfalg2006similarity},
\item Edit distances~\citep{chen:mln,chen2005robust},
\item Longest Common Subsequences (LCSS)~\citep{878994},
\item Swale~\citep{1247544},
\item  SpADe~\citep{chen2007spade},
\item and Cluster, Then Classify (CTC)~\citep{keogh1998enhanced}.
\end{itemize}
\cite{1454226} presented an extensive comparison of these distance functions and
concluded that DTW is among the best measures and that the accuracy of the Euclidean distance
converges to DTW as the size of the training set increases.

In a general Machine Learning setting, \cite{paredes2000class,1221358} compared Euclidean distance with the conventional and class-based Mahalanobis distances.
 One of our contribution is to validate these generic results on time series: instead of tens of features, we have hundreds or even thousands of values which makes the problem mathematically more challenging: the rank of our covariance matrices are
 often tiny compared to their sizes.

More generally, distance metric learning has an extensive literature~\citep{255644,hastie1996discriminant,1598381,short1980new}.
We refer the reader to \cite{weinberger2009distance} for a review.

A conventional distance-learning approach is to find an optimal generalized ellipsoid
distance with respect to a specific loss function.
The LMNN algorithm proposed by \cite{weinberger2009distance} takes a different approach. It seeks to force nearest neighbors to
belong to the same class and it separates instances from different classes by a large margin.
LMNN can be
formulated as a semi-definite programming problem.
They also propose a modification
which they call multiple metrics LMNN as it learns different distances for each class.

There are many extensions and alternatives to NN classification. For example,
\cite{1558510} use instance weights to improve classification.
Meanwhile, \cite{1553530} learn a distance per instance.

\subsection{Dynamic Time Warping (DTW)}

As already stated, The DTW is one of the most accurate distance function for time-series classification.
The DTW was invented
to recognize spoken words~\citep{sakoe1978dpa}, but it is also used
for various problems such as
handwriting recognition~\citep{bahlmann2004wio,niels2005udt},
 chromosome classification~\citep{legrand2007ccu},
 networking~\citep{springerlink:10.1007/978-3-642-15880-3_8}
  or
  shape retrieval~\citep{TPAMI.2005.21,marzal2006cbs}.

For simplicity, consider two time series $x$ and $y$ of equal length ($n$).
  Recall that the $l_p$ distance between two
  time series is \begin{eqnarray*}\sqrt[p]{\sum_{i=1}^n |x_i - y_i|^p}\end{eqnarray*} if $p$ is a positive
  integer or $\max_{i=1,\ldots,n} |x_i - y_i|$ if $p=\infty$.
 The $l_1$ distance is also called the Manhattan distance whereas
 the $l_2$ distance is the Euclidean distance.

The intuition behind the DTW is that one could speak the same sentence by
speeding up and slowing down without changing the meaning of
the sounds.
 The DTW is a generalization of the $l_p$ distance which allows the
 data to be realigned.
To compute the DTW between $x$ and $y$, you must find a
many-to-many matching between the data points in  $x$ and the data points
in $y$. That is each data point from one series must be matched with at least
one data point with the other series. One such matching is the trivial one, which maps the
first data point from $x$ to the first data point in $y$,
the second data point in $x$ to the second data point in $y$, and so on.
Write the
set matches $(i,j)$ as $\Gamma$ so that the trivial matching
is just $\Gamma=\{(1,1), (2,2), \ldots, (n,n) \}$.
The $l_p$ cost corresponding to a matching is defined as
\begin{eqnarray*}\sqrt[p]{\sum_{(i,j)\in \Gamma} |x_i - y_j|^p}\end{eqnarray*} if $p$ is
in an integer or
$\max_{(i,j)\in \Gamma} |x_i - y_j|$ if $p=\infty$.
Typically, $p$ is either 1 or 2: for our purposes we
assume $p=2$.
For a given $p$, the DTW is defined as the minimal cost
over all allowed matchings $\Gamma$.
Typically, we require matchings to be monotonic:
if both $(i,j)$ and $(i+1,j')$ are in $\Gamma$ then $j'\geq j$,
that is, we cannot warp back in time.
Moreover, some matches might be forbidden, maybe because
the data points are too far apart~\citep{itakura1975mpr,sakoe1978dpa}.
\cite{Yu20112787} has proposed learning this warping constraint
from the data.
Except when $p= \infty$, the DTW fails to satisfy the
triangle inequality:
the DTW is not a metric distance.

The computational cost of the DTW is sometimes a challenge~\citep{Salvador2007561}.
To alleviate this problem, several strategies have been proposed
including lower bounds and R*-tree indexes~\citep{ratanamahatana2005tmd,Lemire:2009:FRT:1542560.1542887,Ouyang:2010:HDS:1884499.1884520}.

\cite{Gaudin:2006:ATW:1193211.1193725} proposed a weighted
version of the DTW called Adaptable Time Warping. Instead of
computing $\sum_{(i,j)\in \Gamma} |x_i - y_j|^p$, it computes
$\sum_{(i,j)\in \Gamma} M_{i,j}|x_i - y_j|^p$ where $M$ is some matrix.
Unfortunately, finding the optimal matrix $M$ can be a challenge.
\cite{Jeong2010} investigated another form of weighted DTW
where you seek the minimize
\begin{eqnarray*}\sqrt[p]{\sum_{(i,j)\in \Gamma} w_{|i-j|}|x_i - y_j|^p}\end{eqnarray*}
where $w$ is some weight vector.
Many other variations on the DTW  distance have been
proposed, e.g.,~\cite{springerlink:10.1007/s11634-006-0004-6}.

\section{Mahalanobis distance}
\label{sec:mathstuff}

 For completeness, we derive the Mahalanobis distance~\citep{mahalanobis1936generalized} as an optimal form of generalized ellipsoid distance.
We seek $M$ minimizing
\begin{eqnarray*}
\sum_{x,y \in S}{(x-y)^{\top} M (x-y)} = \sum_{x,y \in S}{\left(\sum_{k=1}^n\sum_{l=1}^n{(x_k-y_k) m_{kl}
(x_l-y_l)}\right)}
\end{eqnarray*}
where $S$ is some class of time series. We add a regularization constraint on the determinant ($\det(M) =
1$).

We solve the minimization
problem by the Lagrange's multiplier method with the Lagrangian
\begin{eqnarray*}
L(M,\lambda) &=& \sum_{x,y \in S}\left(\sum_{k=1}^n\sum_{l=1}^n{(x_k-y_k) m_{k,l} (x_l-y_l)}\right) - \lambda (\det M -1).
\end{eqnarray*}

We want to compute the derivative of the
Lagrangian, and of $\det(M)$, with respect
to $m_{k,l}$. By Laplace expansion, we have for all $l$ that
\begin{eqnarray*}\det(M) = \sum_{k=1}^n{(-1)^{k+l} m_{k,l} \det(M_{k,l})} = 1\end{eqnarray*} where $M_{k,l}$ is the $(k,l)$ minor of $M$: an $(n-1)\times(n-1)$ matrix
obtained by deleting $k$-th row and $l$-th column of $M$.
Thus, we have ${\partial \det(M)}/{\partial m_{k,l}} = (-1)^{k+l} \det(M_{k,l}) $.
Setting the derivatives to zero, we get
\begin{eqnarray*}
\frac{\partial L(M,\lambda)}{\partial m_{k,l}} = \sum_{x,y \in S}{(x_k-y_k)(x_l-y_l)} -
\lambda (-1)^{k+l} \det(M_{k,l})  = 0
\end{eqnarray*}
and therefore
\begin{eqnarray*}
\det(M_{k,l}) = \frac{\sum_{x,y \in S}{(x_k-y_k)(x_l-y_l)}}{\lambda (-1)^{k+l}}.
\end{eqnarray*}
Because $\det(M)=1$ and using Cramer's rule, the inverse matrix $M^{-1}$ can be represented as
\begin{eqnarray*}
m_{k,l}^{-1} = \frac{(-1)^{k+l} \det(M_{k,l})}{\det(M)} = (-1)^{k+l} \det(M_{k,l}).
\end{eqnarray*}
Hence, we have
\begin{eqnarray*}
m_{k,l}^{-1} = \frac{\sum_{x,y \in S}{(x_k-y_k)(x_l-y_l)}}{\lambda }.
\end{eqnarray*}
Thus, we have that
$M^{-1} \propto C$ where $C$ is the covariance matrix. Because we require $M$ to be positive definite and to satisfy $\det(M)=1$, we set $M = (\det(C))^{\frac{1}{n}} C^{-1}$ which produces the Mahalanobis distance\footnote{The original Mahalanobis distance is defined with $M=C^{-1}$}.

Our derivation assumes that the covariance matrix is (numerically) invertible.
This fails in practice.
In \S~\ref{sec:computing}, we review some solutions.

\section{Computing Mahalanobis distances for time series}
\label{sec:computing}

 As a rule of thumb, the covariance matrix becomes singular when the number of instances is smaller or about the same as the number of attributes.
 This is a common problem with time series: whereas individual time series might have thousands of samples, there may only be a few labeled time series in each class.

The most straight-forward solution is to limit the matrix $M$ to its diagonal--thus producing a weighted Euclidean distance.
Revisiting the derivation of \S~\ref{sec:mathstuff} where we require $m_{k,l}=0$ for $k\neq l$, we get that the inverse of the Mahalanobis matrix $M$ must be equal to the inverse
of the diagonal of the covariance matrix:
$M^{-1} \propto \mathrm{diag} ( C )$. As long as the variance of
each attribute in our training sets is different from zero -- a condition
satisfied in practice in our experiments, the problem
is well posed and the result is a positive-definite matrix.
In such a diagonal case, the number of parameters
to learn grows only linearly with the number of attributes in the time series.
In contrast, the number of elements in the full covariance matrix grows
quadratically. The speed of the computation of the distances also depends
on the number of non-zero elements in the Mahalanobis matrix $M$.

 Alas, the diagonal Mahalanobis  distance fails to take into account the information
 off the diagonal in the covariance matrix.
See Figure \ref{fig:signcov} for the covariance matrix of a class of time series. It is clear from the figure that the covariance matrix has significant values off the diagonal.
There are even block-like patterns in the matrix corresponding to specific time intervals.

\begin{figure}
\centering

\subfloat[Sample of time series]{

	\includegraphics[width=0.69\columnwidth]{cylindersignal.pdf}

	}
\subfloat[Sample covariance]{

	\includegraphics[width=0.29\columnwidth]{cylindercovariance.pdf}

	}
\caption{\label{fig:signcov} The Cylinder class from the CBF data set and its sample covariance.
Higher absolute values in the matrix are presented using darker colors.
}
\end{figure}

 Could it be that non-diagonal Mahalanobis  distance could be superior or at least
 competitive with the diagonal Mahalanobis  distance?
 It is tempting to use banded matrices, but the restriction of a positive definite
matrix to a band may fail to be positive definite. Block-diagonal matrices~\citep{Matton20101303} can preserve positive definiteness, but learning which blocks
to use in the context of time series might be difficult.
Instead, we propose two approaches: one is based on the widely used Moore-Penrose pseudoinverse, and the other is covariance shrinkage.
See Figure \ref{fig:signmatrices} for the three different covariance estimates of the same class.

\begin{figure}
\centering

\subfloat[Sample covariance]{

	\includegraphics[width=0.3\columnwidth]{funnel-samplecov-crop.pdf}

	}
\subfloat[Shrinked covariance]{

	\includegraphics[width=0.3\columnwidth]{funnel-shrinkage-crop.pdf}

	}
\subfloat[Diagonal covariance]{

	\includegraphics[width=0.3\columnwidth]{funnel-diagonal-crop.pdf}

	}
\caption{\label{fig:signmatrices}
The covariance estimates of the Funnel class in the CBF data set.
Large absolute values are in darker colors.
Both the shrinked
and diagonal covariances are positive definite whereas
the sample covariance matrix is singular. }
\end{figure}

The approach based on the pseudoinverse is based on the singular value decomposition (SVD). We write the SVD as $C= U\Sigma V^*$ where $\Sigma$ is a diagonal matrix with eigenvalues $\lambda_1, \lambda_2,\ldots$ and $U$ and $V$ are unitary matrices.
The Moore-Penrose pseudoinverse is given by $V\Sigma^+ U^*$ where $\Sigma^+$
is the diagonal matrix made of the eigenvalues $1/\lambda_1, 1/\lambda_2,\ldots$ with
the convention that $1/0=0$. The pseudo-determinant is the product of the non-zero eigenvalues of $\Sigma$.
We set $M$ equal to the
pseudoinverse of the covariance matrix---normalized so that it has a pseudo-determinant
of one. This solution is equivalent to projecting the time series data on the
subspace corresponding to the non-zero eigenvalues of $\Sigma$. That is, the matrix $M$ will be singular.

Covariance shrinkage is an estimation method for problems with small number of instances and large number of attributes~\citep{stein}. It has better theoretical and practical properties for such data sets as the estimated covariance matrix is guaranteed to be non-singular.
Let $x_1, x_2, \ldots, x_n$ be a set of time series. We write $x_{i}= (x_{i1},x_{i2},\ldots,x_{ip}) $ where $x_{ik}$ is the $k^{\mathrm{th}}$ value (attribute) of the time series $x_i$.
One element of the (sample) covariance matrix $S$ is\begin{eqnarray*}
s_{ij} = \frac{1}{n-1} \sum_{k=1}^n (x_{ki} - \bar x_i )(x_{kj} - \bar x_j )
 \end{eqnarray*}
 where $\bar x_i = \frac{1}{n} \sum_{k=1}^n x_{ki}$.
It is positive semi-definite but can be singular.
To fix the problem that $S$ is singular, we can replace it
with  an estimation of the form \begin{eqnarray*}C^*=\lambda T + (1- \lambda) S\end{eqnarray*} for some suitably chosen
target $T$. If $T$ is a  positive definite matrix and $\lambda\in (0,1]$, we have that  $\lambda T + (1- \lambda) S$
must be positive definite. Moreover,
the smallest eigenvalue of $\lambda T + (1- \lambda) S$ must be at least as large as $\lambda$ times the smallest eigenvalue of $T$.
We have used the target recommended by \cite{schafer2005shrinkage} which is the diagonal of the unrestricted covariance estimate, $T = \mathrm{diag} ( C )$.
It is positive definite in our examples. For $\lambda$, we use the estimation proposed by \cite{schafer2005shrinkage}: it is computationally inexpensive.

Thus, finally, we consider six types of Mahalanobis distances: two localities (global or class-based) and three estimators (pseudoinverse, shrinkage, or diagonal).

\section{Experiments}

The main goal of our experiments is to evaluate Mahalanobis distances and the class-based approach on time series.
A secondary goal is to evaluate the LMNN method.

We begin all tests with training data set made of classes of instances. When applicable,
distances are learned from this data set. We then attempt to classify some test
data using 1-NN and we report the classification error.

The code for the experiments is available online~\citep{ourmatlabcode} with instructions on how the results can
be reproduced. For LMNN, we use the source code provided by \cite{lmnncode} for the experiments
with default parameters.

\subsection{Data sets}
We use the UCR time series classification benchmark~\citep{keoghbenchmark} for our experiments as it includes
diverse time series data sets from many domains.
It has predefined training-test splits for the experiments,
so the results can be compared across different papers.
We removed the two data that are not z-normalized by default (Beef and Coffee).
Indeed, z-normalization improves substantially the classification accuracy---irrespective of the chosen distance.
Thus, for fair results, we should z-normalize them, but this may create confusion with previously
reported numbers.
We also removed the Wafer data set as all distances classify it nearly perfectly.
The remaining 17~data sets were used for the comparison of different methods.

\subsection{Best Mahalanobis distance for 1-NN accuracy}

We compare the various Mahalanobis distances in Table~\ref{tab:ourresultsforMahalanobis}. We have left out the Moore-Penrose pseudoinverse, because its error rates were twice as high on average compared to the other variants.
What is immediately apparent is that the class-based metrics give better classification results.

The diagonal Mahalanobis is somewhat better and they are also considerably faster computationally, but the shrinkage estimate yields significantly better results for several data sets (e.g.~Adiac and Face~(four)).
Thus, out of the six variations, we recommend the class-based shrinkage estimate and the class-based diagonal Mahalanobis distance.

\begin{table}
\caption{\label{tab:ourresultsforMahalanobis}Classification error for the various Mahalanobis distances.}
\centering
\begin{tabular}{l|cccc}\hline
\multirow{2}{*}{Data set}           & \multicolumn{2}{c}{Shrinkage}  & \multicolumn{2}{c}{Diagonal}              \\
             & global &  class-based         & global & class-based \\\hline
50 words			& 0.36		& 0.71	&  0.34			&  \textbf{0.32}\\
Adiac               & 0.33     & \textbf{0.28}          &  0.37         &  0.36\\
CBF                 &  0.52     & \textbf{0.04} &  0.16         &   0.05 \\
ECG                 &  0.13     & 0.09          &0.10           & \textbf{0.08}     \\
Fish                &  0.33     & \textbf{0.15}          & 0.19 & 0.18  \\
Face (all)          &  0.31     & 0.27	&  0.32         & \textbf{0.25}\\
Face (four)         &  0.45     & \textbf{0.10}          & 0.16         &0.17  \\
Gun-Point           &  \textbf{0.06}     & 0.10          & 0.10          & 0.11 \\
Lighting-2          &  0.49     & 0.31 &  \textbf{0.25}          &  \textbf{0.25}  \\
Lighting-7          &  0.59     & 0.32          & 0.36          &  \textbf{0.23} \\
OSU Leaf            &  0.69     & 0.69 & \textbf{0.46}          & \textbf{0.46}    \\
OliveOil          &0.17			& 0.17  & 0.17          &\textbf{0.13}  \\
Swedish Leaf        & 0.25     & \textbf{0.14}          & 0.21          & 0.18 \\
Trace               &  0.27     & 0.09 & 0.21          & \textbf{0.07}   \\
Two Patterns        &  \textbf{0.10}     & \textbf{0.10} & 0.12          &  0.12\\
Synthetic Control   &  0.23     & 0.10 & 0.13          &  \textbf{0.09}\\
Yoga              & 0.24 		& 0.21	& \textbf{0.17}         & \textbf{0.17}	 \\\hline
\# of best errors   & 2 & 6 & 3 & 10 \\ \hline
\end{tabular}
\end{table}

\subsection{Comparing competitive distances}

How do the class-based Mahalanobis distances fare compared to competitive distances? Computationally,
the diagonal Mahalanobis is inexpensive compared to schemes such as the DTW or LMNN\@.
Regarding the 1-NN classification error rate, we give the results in Table~\ref{tab:ourresults}.
As expected~\citep{1454226}, no distance is better on all data sets. However, because the diagonal Mahalanobis distance is closely related
to the Euclidean distance, we compare their classification accuracy.
In two data sets, the Euclidean distance outperformed the class-based Mahalanobis distance and only by small differences (0.09 versus 0.10-0.12). Meanwhile, the class-based diagonal Mahalanobis outperformed the Euclidean distance 12 times, and sometimes by large margins (0.07 versus 0.24 and 0.05 versus 0.15).

The LMNN is also competitive: its classification error is sometimes half that of the Euclidean distance.
However, the class-based LMNN gets the best result among all methods only twice as opposed to
five~times for the global LMNN. Moreover, the global LMNN significantly outperforms the class-based
LMNN on the Two Patterns data set (0.05 versus 0.24). For time series data sets, the class-based LMNN is not an improvement over the global LMNN.

We have to note that DTW has the lowest error rates and provides best results for half of the data sets, but it is much slower than Mahalanobis distances.

\begin{table}
\caption{\label{tab:ourresults}Classification errors for some competitive schemes.  We use class-based Mahalanobis distances. For the 50~words data set, the LMNN computation fails because it has a class with only one instance.}
\centering
\begin{tabular}{l|ccccccc}\hline
\multirow{1}{*}{Data set}           & \multirow{1}{*}{Euclidean} & \multirow{1}{*}{DTW}           & \multicolumn{2}{c}{C.-b. Mahalanobis} & \multicolumn{2}{c}{LMNN}  &  \\
             &  &     &   shrink.   & diag.  &    & c.-b. \\\hline
50 words			&  0.37		& \textbf{0.31}	& 0.71	&   0.32				& ---			 & ---\\
Adiac               &  0.39     & 0.40          & 0.28         &   0.36               & \textbf{0.23}      & 0.32\\
CBF                 &  0.15     & \textbf{0.00} & 0.04         &   0.05               &  0.15           &  0.15 \\
ECG                 &  0.12     & 0.23          & 0.09           & 0.08        	&  0.10          & \textbf{0.07}\\
Fish                &  0.22     & 0.17      	& 0.15	     & 0.18                 &  \textbf{0.13}  & 0.14\\
Face (all)          &  0.29     & 0.19			& 0.27        & 0.25         &  \textbf{0.16}           & 0.20\\
Face (four)         &  0.22     & 0.17          & \textbf{0.10}     & 0.17                  & 0.16   &  0.16\\
Gun-Point           &  0.09     & 0.09          & 0.10        & 0.11                 &\textbf{0.05}      & 0.09\\
Lighting-2          &  0.25     & \textbf{0.13} & 0.31          & 0.25                 &  0.41            & 0.34\\
Lighting-7          &  0.42     & 0.27          & 0.32         &  \textbf{0.23}       &  0.51            &  0.48\\
OSU Leaf            &  0.48     & \textbf{0.41} & 0.69         & 0.46                 & 0.57             & 0.54\\
OliveOil          &\textbf{0.13}&\textbf{0.13}  & 0.17         &\textbf{0.13}         & \textbf{0.13}    & \textbf{0.13}\\
Swedish Leaf        &  0.21     & 0.21          & \textbf{0.14}   & 0.18        & 0.21            & 0.19\\
Trace               &  0.24     & \textbf{0.00} & 0.09          & 0.07                 & 0.20            & 0.20\\
Two Patterns        &  0.09     & \textbf{0.00} & 0.10          &  0.12                & 0.05		  	& 0.24\\
Synthetic Control   &  0.12     & \textbf{0.01} & 0.10          & 0.09                 & 0.03          & 0.09\\
Yoga              & 0.17 		& \textbf{0.16}	& 0.21          & 0.17					& 0.18		 & 0.18 \\
\hline
\# of best errors   & 1 & 9 & 2 & 2 & 5 & 2\\ \hline
\end{tabular}
\end{table}

\subsection{Effect of the number of instances per  class}

Whereas Table~\ref{tab:ourresults} shows that the Mahalanobis distances are far superior to the Euclidean distance on some data sets, this result
 is linked to the number of instances
per class. For example, on the Wafer data set (which we removed), there are many instances per class (500), and correspondingly, all distances give a negligible classification error.

Thus, we considered three different synthetic time-series data sets with varying numbers of instances per class:
Cylinder-Bell-Funnel (CBF)~\citep{921732},
  Control Charts (CC)~\citep{pham1998ccp} and
 Waveform~\citep{breiman1998car}.
Test sets have 1\,000~instances per class whereas training sets have between 10 to 1\,000 instances.
We repeated each test ten times, with different training sets.
Fig.~\ref{fig:synthdatasets} shows that whereas the class-based diagonal Mahalanobis
is superior to the Euclidean distance when there are few instances, this benefit
is less significant as the number of instances increases. Indeed, the classification
accuracy of the Euclidean distance grows closer to perfection and it becomes more difficult
for alternatives to be far superior.

\begin{figure}
\centering\includegraphics[width=0.8\columnwidth]{synthdatasets}
\caption{\label{fig:synthdatasets}Ratios of the 1-NN classification accuracies using the class-based diagonal Mahalanobis and Euclidean distances}
\end{figure}

\section{Conclusion}
The Mahalanobis distances have received little attention for time series classification and we are not
surprised given their poor performance
as a 1-NN classifier when used in a straight-forward manner. However, by learning one Mahalanobis distance per class we get a competitive classifier when using either covariance shrinkage or a diagonal approach.
Moreover, the diagonal Mahalanobis distance is particularly appealing computationally: we only
need to compute the variances of the components.
Meanwhile, we get good results with the LMNN on time series data, though it is more expensive.
The DTW is superior, but computationally expensive.

\begin{acknowledgements}
This work is supported by NSERC grant 261437.
\end{acknowledgements}

\bibliographystyle{spbasic}


\section{Preamble: Direct detection of binarity with high spatial resolution imaging}
A major breakthrough on the study on asymmetrical Planetary Nebulae (PNs) has been the recognition by a large part of our community of the growing importance of binary systems as main shaping agent of the bipolar and asymmetrical PNs. Companions encompassing a large range of mass, from the stellar objects to jovian-mass planets are suspected to deeply influence the ejecta when the star reaches the AGBs or even as early as the RGB \citep{2009PASP..121..316D}. This interaction can potentially influence dramatically the fate of the star, leading to poorly known evolutionary paths, influencing deeply the time scales of the different evolutionary stages, eventually bearing only little similarities compared to the time scales involved for the evolution of a single, naked star \citep{2007BaltA..16..104F,2007BaltA..16..112V}. Nevertheless, it is still far from excluded that a single star may also provide the conditions for the shaping of an asymmetrical nebula, if it expels at one time or another some mass with a significant pole-to-equator density gradient \citep{2002ARA&A..40..439B}.

Many detections of companion of PNs central stars were reported during this conference \citep[De Marco, Hrivnak, Hajduk, Miszalski...these proceedings]{2010NewA...15..483B, 2009A&A...496..813M} some based on difficult long-term spectral monitoring, and many on the recent extensive photometric campaigns from automated telescopes. The high spatial resolution techniques may at first sight represent a large potential for detecting some well-separated  (2-100\,milli-arcsecond, hereafter mas) companions, but they suffer from intrinsic constraints that do not allow them to currently play a significant role in that domain. Adaptive optics on 8m-class telescopes are mostly limited in terms of spatial resolution, whilst optical interferometry has strong limitations in terms of contrast and imaging capabilities. The brighter the circumstellar environment and the central source, the harder will be the detectability of a small point-like structure in the close vicinity. The best case is the detection of low $Teff$ companions ($\Delta \sim$5mag, d$\sim$10-100mas) around hot, and preferentially naked sources for which it is easy to separate the SED of the stellar components. For a review of the methods used for detecting binarity, the reader is referred to \citet{2008AIPC.1057....1J}.

\section{Defining the different kind of equatorial over-densities}
Whatever the origin of the shaping process of an asymmetrical planetary nebula, an equatorial over-density of material is often involved in the models, and a growing number of observations unveil their presence at many stages of the star evolution. A wealth of new high spatial resolution techniques are now routinely available for the observer, namely the adaptive optics techniques, the speckle and lucky imaging techniques (often called the 'burst modes') in the optical and the infrared, and also the interferometry in the optical, infrared, millimetric and radio wavelengths. All these techniques have their own spatial resolution, sensitivity, contrast and astrometric constraints to the point that it is often quite difficult to compare and put into context the outcome of these observations.
Hence, the vocabulary can often be somehow misleading. The detection of disks around evolved sources, as in  post-AGB stars, may be closely associated with a high probability to see a binary system. However, this claim is highly dependent on the kind of structure encompassed in the 'disk' denotation, and one shall in the following sections precise further the equatorial overdensities, dividing them into the 'torus' (or 'outflowing wind'), and the 'stratified disks' families. The detection of compact, hot, dust-less accretion disks whose SEDs peak in the UV/B and that are hardly detectable by the above mentioned infrared techniques is out of the scope of this review.

\subsection{Dusty torii: wind related equatorial overdensities}
An equatorial overdensity is a region of higher density and lower expansion speed compared to the polar circumstellar regions. The kinematics of such a structure is dominantly radial, and the total angular momentum carried by the structure is limited. Equatorial overdensities are short-term structures (these are deflected winds), e.g. if the mass-loss ceases, then the fate of the material is to rapidly expand and vanish. The increase of density toward the equatorial plane can be ascribed easily by a dependence on the co-latitude of the star in spherical coordinates, and there is no mass nor energy storage in the structure.

To date, the best examples of such torii originate from millimetric interferometry that associates a good spatial resolution ($\sim$1") and an excellent spectral resolution (i.e. 1\,km.s$^{-1}$). One can cite as good example of expanding torii, \citet{2008ApJ...673..934D} or \citet{2007A&A...473..207P}. In the near- or mid-IR, it is more difficult to access the density distribution of the circumstellar material due to opacity effects. Nevertheless, if the scale height of the structure is comparable to its radial extension, without clear sign of marked density increase toward the equatorial plane, the torus hypothesis can be favored \citep{2009A&A...503..837V,2006A&A...448..203L}. Long slit spectroscopy can provide further evidence for this classification when a significant expansion velocity (i.e. $\sim$15-40 km.s$^{-1}$) of the structure can be measured \citep{2010MNRAS.401..405J,2010arXiv1006.5873J}.

\subsection{Stratified disks}
A disk exhibits a clear vertical stratification whose scale height is governed by the gas pressure only. Its aperture angle is very small (i.e. less than $\sim$10 degrees) and its kinematics is dominated by (quasi) Keplerian velocities, with a small expansion component ($\leq$10\,km.s$^{-1}$).  Thus its lifetime is much longer than structures blown within a wind, competing or even exceeding the reference time scale of a PN lifetime of a few tens of thousand years \citep{2003ARA&A..41..391V}.

The density structure of a Keplerian disk is best ascribed in cylindrical coordinates with a radial law in the equatorial plane, and a vertical stratification perpendicular to this plane, for the hydrostatic equilibrium.

Such a model is used in a wealth of astronomical contexts, and in particular for the formation of Young Stellar Objects (YSOs). Some examples of application of a similar model can be found in this list: \citet{2009A&A...500.1065D,2008A&A...491..809F,2003ApJ...588..373W,2003A&A...398..607D}.

The density law used in this model is
\begin{equation}
	\rho\left(r,z\right)=\rho_{0}\left(\frac{R_{*}}{r}\right)^\alpha \exp{\left[-\frac{1}{2}\left(\frac{z}{h\left(r\right)}\right)^2\right]}
\end{equation}
\begin{equation}
	h\left(r\right)=h_{0}\left(\frac{r}{R_{*}}\right)^\beta
\end{equation}
\noindent
where {\it r} is the radial distance in the disk's mid-plane, $R_{*}$ is the stellar radius, $\beta$ defines the flaring of the disk, $\alpha$ defines the density law in the mid-plane and $h_{0}$ is the scale height at a reference distance (often 100 AU)  from the star.
Until recently, such Keplerian disks could only be studied by millimetric interferometry at the highest possible spatial resolution and for the closest targets. One of the best example of such disks is the Red Rectangle \citep{2009ApJ...693.1946W,2005A&A...441.1031B}. Again, the great advantage of such a technique is to provide both the spatial distribution and the kinematics of the source \citep{2007A&A...468L..45B}.

In absence of any information on the kinematics in the infrared due to the limited spectral resolution of optical interferometric observations, one has to rely on some indirect evidence to claim for the discovery of a stratified dusty disk, namely one has to prove at least that the disk opening angle is small, or better, that the density law follows radially and vertically the model described above. This is a possible task when the disk is seen close to edge-on. In this case, the radial and vertical directions are well separated on the sky.
The best example is the discovery of the edge-on disk in the core of Menzel 3, very well fitted by a stratified disk model \citep{2007A&A...473L..29C}. However, our group has conducted some tests trying to invert the disks parameters from artificial interferometric datasets (Niccolini et al., A\&A, accepted). When the disk is seen at low inclination, many degeneracies appear between the parameters of the vertical density law (such as the flaring parameter $\beta$), and the parameters of the radial law. Moreover, the size of the dust grains and and their composition affect also critically the fits.

To improve the constraints on the disk temperature law, a good approach is to perform near-IR and mid-IR observations using the AMBER and MIDI instruments of the VLTI respectively. The post-AGB binary IRAS 08544-4431 was studied this way \citep{2007A&A...474L..45D} and the same strategy is used for YSOs \citep{2008ApJ...676..490K}. It is more difficult to model the near-IR interferometric data due to the potential spatial complexity of the dusty disk's inner rim, and an intense theoretical and observational effort is currently performed, mainly in the YSOs community to better understand the so-called 'puffed-up inner rims' \citep{2010arXiv1006.3485D,2009A&A...506.1199K,2007ApJ...661..374T}. A better determination of the spatial properties and fine chemical content of the dust forming region is a challenge for the future. Other interesting targets are the double-chemistry sources such as BM\,Gem, in which a companion has recently been discovered \citep{2008ApJ...682..499I, 2008A&A...490..173O, 2006A&A...445.1015O}. These sources seem to harbor systemically an equatorially enhanced circumstellar environment, but there is no firm confirmation yet that the structure is best described by a torus or a disk.

\section{Disk evolution}
\subsection{Stratified disk and binarity}

Stratified disks detected in evolved systems are potentially highly correlated with binaries as demonstrated by Van Winckel and collaborators on the environment of binary post-AGBs \citep[Gielen et al. these proceedings]{2008A&A...490..725G,2007BaltA..16..112V,2006MmSAI..77..943V,2006A&A...448..641D}. Grain growth, settling, radial mixing and crystallization are efficient in such an environment and the circumbinary disc of these sources seems to be governed by the same physical processes that govern the proto-planetary discs around young stellar objects. It seems that another distinctive character of these long-lived stratified disks as seen in the infrared would be their content in highly processed grains \citep{2010A&A...510A..30M,2008A&A...490..725G, 2006A&A...451..951I}.
The key point for the stabilization of the disk is to provide enough angular momentum \citep{2008NewA...13..157A, 1997ApJS..112..487S}. A star may (hardly) supply this angular momentum via magnetic fields \citep{2004ApJ...614..737F, 2002ApJ...576..413U}. But even in this case the formation and stabilization of a Keplerian disk remains a challenge. Of course, this argument does not apply to the accretion disks encountered around YSOs, for which the difficulty in the contrary is to understand how the excess of angular momentum is dissipated.
The angular momentum provided by a low-mass stellar or even sub-stellar companion may potentially have a dramatic influence on a growing RGB or AGB star \citep{2010arXiv1002.2216N, 2010NewA...15..483B, 2007MNRAS.376..599N, 2006MNRAS.370.2004N,2000ApJ...538..241S}. The presence of a stratified disk, and the associated Lindblad resonances, seem to be a key ingredient in the orbital evolution of the binary system. Without this ingredient, it is difficult to reproduce the observed morphology of the eccentricity - period diagram (Dermine T. et al., in press)

{\it My personal opinion is that the discovery of a stratified disk with proven Keplerian kinematics is directly connected to the influence of a companion, albeit the few exceptions presented above, namely the Young Stellar Objects or the critical velocity rotating massive sources. This hypothesis must be confirmed by further observations. }

\subsection{Time scale of formation}
When dealing with the theoretical building-up of a disk, time matters, and any constraint on the time-scale on which the observed structures were formed is of importance. An interesting study was presented by \citet{2007ApJ...663..342H} on the close time-scale connection between jets and torii, based on many kinematical studies. $\pi^1$ Gruis is a good example of a recently formed structure \citep{2008A&A...482..561S, 2006ApJ...645..605C}. OH231.8+4.2 is another example of a similar process caught in the act in IRAS16342-3814 \citep{2010arXiv1004.2659E, 2006ApJ...646L.123M,2007A&A...468..189D}.

Some recent examples show how fast may be the building up of an equatorial overdensity.
The born-again stars, V605\,Aql \citep{2008A&A...479..817H} and the Sakurai's Object are surrounded by an equatorial overdensity with a large scale height that may be described as a torus \citep{2009A&A...493L..17C}. In the case of the Sakurai's object, the torus was detected in 2007, about 10 years after the dust began to form. A dense equatorial structure was formed even faster, in less than 2yrs, around the slow dust-forming nova V1280 Sco \citep[Chesneau et al. in preparation]{2008A&A...487..223C}. There is no doubt in this case that the slow ($\sim$500km.s$^{-1}$) ejecta from the outbursting nova were deeply affected by the common envelope phase that lasted more than tens of orbital period of the companion, leading to the fast formation of a bipolar nebula (see also Evans, these proceedings).

These examples show that as soon as the primary in a binary system gets larger and increases its mass-loss rate when evolving, the influence of a companion can rapidly focus the ejected material onto the equatorial plane of the system, leading to an equatorial overdensity, as already investigated theoretically \citep{1998ApJ...497..303M,1993MNRAS.265..946T}. This physical process depends on many parameters (mass ratio, orbital parameters, mass-loss rate of the primary...), but the efficiency is such that a large number of targets might be affected at one stage or another \citep{2010NewA...15..483B}.

\subsection{Time scale of dissipation and fate}
The observations of the inner circumstellar structures around evolved sources have been to date too scarce to put them into an evolutionary sequence. As written above, the torii are supposed to expand and dissipate much faster than the Keplerian disks. One may even consider the extreme case in which a binary system surrounded by a stable circumbinary disk continuously replenished by the interaction of the stars may remain virtually unchanged for time scale as long as many 10$^5$-10$^6$years, as proposed for the Red Rectangle. A proposed unified picture of the such an interacting binary evolution involving the presence of a torus is presented in \citet{2007BaltA..16..104F}.

What is the fate of a stratified disk? An 'old' dissipating disk can see its density and kinematical structure deeply affected by the fast evolving wind of the central star. \citet{2010A&A...514A..54G} performed an in-depth kinematical study of the dissipating disk in the core of the M2-29 nebula \citep{2008A&A...490L...7H}. The torus found around Hen2-113 might also be dissipating under the influence of the fast radiative wind emitted by the Wolf-Rayet star in the core of the nebula \citep{2006A&A...448..203L}.

\section{Stratified disks around massive stars}
There is now a large bunch of evidence that the disks encountered around Be stars and at least some B[e] supergiants are rotating close to Keplerian velocities \citep{2007A&A...464...59M}. The Be stars are proven statistically to represent the tail of the fastest rotators among B stars, and their disks of plasma may be explained without invoking the influence of a companion, by a combination of an extreme centrifugal force at the equator and some pulsation properties of the star which leads to erratic ejection of material with high angular momentum \citep{2009ApJ...701..396C, 2009A&A...506...95H}. By contrast, there is still no consensus for the more massive supergiant counterparts, the B[e] stars that are surrounded by dense disks of plasma {\it and} dust. The dust survives much closer to this hot star than expected so far  \citep{2010A&A...512A..73M, 2009A&A...507..317M, 2007A&A...464...81D}.The B[e] supergiant's rotation rate is strongly decreased by the increase of their radius while leaving the main sequence, and the rotation alone is probably far from being sufficient to explain the formation of a disk without invoking the influence of a close companion \citep{2007ApJ...667..497M}. Note also, that dusty circumbinary disk are commonly encountered around interacting binary systems \citep{2009A&A...499..827N}.

The example of the A[e] supergiant HD\,62623 is informative in this context \citep[Millour, Meilland et al., submitted]{2010A&A...512A..73M}. HD\,62623 is an A supergiant showing the characteristics of the 'B[e]' spectral type, namely a spectrum dominated by strong emission lines and a large infrared excess. Spectrally and spatially resolved observations of AMBER/VLTI in the Br$\gamma$ line have shown that the supergiant lies in a cavity, and is surrounded by a dense disk of plasma. The Br$\gamma$ line in the location of the central star is {\it in absorption} showing that the star is not different from a normal member of its class such as Deneb (A3Ia), albeit with a significantly large v$sini$ of about 50km.s$^{-1}$. By contrast, the Balmer and Bracket lines are wider ($vsini$ of about 120km.s$^{-1}$), and the AMBER observations demonstrated that they originate from a disk of plasma, most probably in Keplerian rotation (Millour, Meilland et al. submitted). In absence of any proof of binarity, it is often difficult to understand how such a dense equatorial disk could have been generated. However, HD62623 is a known binary with a stellar companion that orbits close to the supergiant with a period of about 136 days \citep{1995A&A...293..363P}. The mass ratio inferred is very large, and the companion is probably a solar mass star. \citet{1995A&A...293..363P} proposed that an efficient angular momentum transfer occurs near the L2 Lagrangian point of the system, propelling the mass lost from the supergiant by its radiative wind and probably also by strong tides into a stable dense circumbinary disk \citep[in the context of an AGB star]{2009A&A...507..891D}. A similar idea was proposed by \citet{2009A&A...507..891D} in relation with radiation pressure acting on the wind of AGB stars and modifying the Roche lobe geometry, therefore probably easing the formation of a circumbinary disk.

The comparison between the disks encountered around low and intermediate mass stars and those observed around the B[e] supergiants, a rare spectral type among the zoo of massive stars might not appear relevant at first sight. Yet, recent Spitzer observations of 9 LMC B[e] stars showed an interesting homogeneity of their spectra, and a great similarity with the spectra of post-AGBs harboring dense dusty disks \citep{2010AJ....139.1993K}. This tightens further the connection between B[e] stars and binarity.

\acknowledgements This article benefited from fruitful discussions with many persons, and in particular Alain Jorissen.