\subsection{Notations}

For any $n \in \mathbb{N}$ let $\mathbb{Y}_n$ denote the \emph{set of Young diagrams} with $n$ boxes.
By $\square_N\in\mathbb{Y}_{N^2}$ we~denote \emph{the square Young diagram} of side $N$.
For a Young diagram $\lambda \in \mathbb{Y}_n$ let $|\lambda| := n$ be the \emph{size} of diagram $\lambda$
and by $\mathcal{T}_\lambda$ denote the \emph{set of standard
	Young tableaux of~shape~$\lambda$}.
On~the~set~$\mathcal{T}_{\square_N}$ of square
tableaux with $N^2$ boxes we consider \emph{the uniform probability
	measure}~$\mathbb{P}_N$.

We draw Young tableaux in the \emph{French notation}.
For a tableau $T$, by $\operatorname{pos}_k(T) = (x_k, y_k)$ we denote the
\emph{Cartesian coordinates of entry~$k$}, i.e.~$x_k$ is its column and $y_k$
is its row. The~difference $u_k^T:=x_k- y_k$ will be called \emph{the
    $u$-coordinate} of the box $k$ (in the literature the name \emph{content} is
also used).

The irreducible representation of the symmetric group $\mathfrak{S}_n$ which corresponds
to~\mbox{$\lambda\in\mathbb{Y}_n$} will be denoted by $V^\lambda$
and its character by $\chi^\lambda$.

\subsection{Jeu de taquin, jeu de taquin path, evacuation path}

\bigskip

\textbf{\emph{Jeu de taquin}} is a transformation of tableaux which was
introduced by Sch\"utzenberger~\cite{Schuetzenberger1977}. It is a fundamental
tool for studying the combinatorics of Young tableaux,
Robinson--Schensted--Knuth (RSK) algorithm, and Littlewood--Richardson
coefficients
\cite{Fulton1997}.

Jeu de taquin acts as follows (see~\cref{fig:jdt}\subref{fig:jdtA} and
\subref{fig:jdtB}): we start with the bottom-left corner box of the given
tableau $T$, we erase it and obtain a~\emph{hole} in its place. Then we look
at the neigboring two boxes: the one to the right and the one  above the hole,
and choose the smaller one. We slide the selected box into the hole; as a
result the hole moves one box to the right or up. We continue sliding  as long as there
is some box to right or up to the hole. The resulting tableau will be denoted
by $j(T)$. The zig-zag path traversed by~the~hole will be called
\textbf{\emph{jeu de taquin path}}.

For a given tableau $T\in \mathcal{T}_\lambda$ with $n=|\lambda|$ boxes we may
iterate jeu de taquin $n$ times until we end with the empty tableau. During
each iteration the box with the biggest number $n$ either moves one node left
or down, or stays put. Its trajectory
\begin{equation}
\label{eq:evacuation}
 \on{pos}_n(T),\
\on{pos}_n\big( j(T) \big),\
\on{pos}_n\big( j^2(T) \big),\ \dots,\
\on{pos}_n\big( j^{n-1}(T) \big)
\end{equation}
will be called \textbf{\emph{evacuation path}}.

\begin{figure}[b]

    \centering
    \subfloat[]{
			\begin{tikzpicture}[scale=0.65]
					\fill[blue!20] (0,0) rectangle +(1,1);
					\fill[blue!20] (0,1) rectangle +(1,1);
					\fill[blue!20] (1,1) rectangle +(1,1);
					\fill[blue!20] (2,1) rectangle +(1,1);
					\fill[blue!20] (2,2) rectangle +(1,1);

					\begin{scope}
					\clip[](5,0) -- (5,1) -- (4,1) -- (4,2) -- (3,2) -- (3,3) -- (1,3) -- (1,4) -- (0,4) -- (0,0);
					\draw (0,0) grid (10,10);
					\end{scope}
					\draw[thick] (5,0) -- (5,1) -- (4,1) -- (4,2) -- (3,2) -- (3,3) -- (1,3) -- (1,4) -- (0,4) -- (0,0) -- cycle;
					\draw (0.5,0.5) node {$1$};
					\draw (0.5,1.5) node {$2$};
					\draw (1.5,0.5) node {$3$};
					\draw (2.5,0.5) node {$7$};
					\draw (3.5,0.5) node {$10$};
					\draw (4.5,0.5) node {$13$};
					\draw (1.5,1.5) node {$4$};
					\draw (0.5,2.5) node {$5$};
					\draw (2.5,1.5) node {$6$};
					\draw (1.5,2.5) node {$8$};
					\draw (0.5,3.5) node {$11$};
					\draw (2.5,2.5) node {$9$};
					\draw (3.5, 1.5) node {$12$};

					\end{tikzpicture}
\label{fig:jdtA}
}
\hfill
\subfloat[]{
			\begin{tikzpicture}[scale=0.65]
					\fill[blue!20] (0,0) rectangle +(1,1);
					\fill[blue!20] (0,1) rectangle +(1,1);
					\fill[blue!20] (1,1) rectangle +(1,1);
					\fill[blue!20] (2,1) rectangle +(1,1);

					\begin{scope}
					\clip[](5,0) -- (5,1) -- (4,1) -- (4,2) -- (2,2) -- (2,3) -- (1,3) -- (1,4) -- (0,4) -- (0,0);
					\draw (0,0) grid (10,10);
					\end{scope}
					\draw[thick] (5,0) -- (5,1) -- (4,1) -- (4,2) -- (2,2) -- (2,3) -- (1,3) -- (1,4) -- (0,4) -- (0,0) -- cycle;
					\draw (0.5,0.5) node {$2$};
					\draw (0.5,1.5) node {$4$};
					\draw (1.5,0.5) node {$3$};
					\draw (2.5,0.5) node {$7$};
					\draw (3.5,0.5) node {$10$};
					\draw (4.5,0.5) node {$13$};
					\draw (1.5,1.5) node {$6$};
					\draw (0.5,2.5) node {$5$};
					\draw (2.5,1.5) node {$9$};
					\draw (1.5,2.5) node {$8$};
					\draw (0.5,3.5) node {$11$};
					\draw (3.5, 1.5) node {$12$};
			\end{tikzpicture}
			\label{fig:jdtB}
		}
\hfill
\subfloat[]{
	\begin{tikzpicture}[scale=0.65]
	\fill[blue!20] (0,0) rectangle +(1,1);
	\fill[blue!20] (0,1) rectangle +(1,1);
	\fill[blue!20] (1,1) rectangle +(1,1);
	\fill[blue!20] (2,1) rectangle +(1,1);
	\fill[blue!40, green!60] (2,2) rectangle +(1,1);

	\begin{scope}
	\clip[](5,0) -- (5,1) -- (4,1) -- (4,2) -- (3,2) -- (3,3) -- (1,3) -- (1,4) -- (0,4) -- (0,0);
	\draw (0,0) grid (10,10);
	\end{scope}
	\draw[thick] (5,0) -- (5,1) -- (4,1) -- (4,2) -- (3,2) -- (3,3) -- (1,3) -- (1,4) -- (0,4) -- (0,0) -- cycle;
	\draw (0.5,0.5) node {$1$};
	\draw (0.5,1.5) node {$3$};
	\draw (1.5,0.5) node {$2$};
	\draw (2.5,0.5) node {$6$};
	\draw (3.5,0.5) node {$9$};
	\draw (4.5,0.5) node {$12$};
	\draw (1.5,1.5) node {$5$};
	\draw (0.5,2.5) node {$4$};
	\draw (2.5,1.5) node {$8$};
	\draw (1.5,2.5) node {$7$};
	\draw (0.5,3.5) node {$10$};
	\draw (3.5, 1.5) node {$11$};
	\draw (2.5, 2.5) node {$13$};
	\end{tikzpicture}
	\label{fig:jdtC} } \caption{\protect\subref{fig:jdtA} Standard Young tableau
$T$ of shape $\lambda$. The highlighted boxes form the~\emph{jeu de taquin
    path}. \protect\subref{fig:jdtB} The outcome $j(T)$ of jeu de taquin
transformation. \protect\subref{fig:jdtC} A bijection on $\mathcal{T}_\lambda$,
\emph{modified jeu de taquin $J(T)$}, is defined by creating an
additional box with~entry~\mbox{$|\lambda|+1$} at~the~final position of the
hole and then decreasing all entries by 1.} \label{fig:jdt}
\end{figure}

\subsection{The problem}
\label{sec:problems}

If we draw the boxes of a given square tableau $T_N\in \mathcal{T}_{\square_N}$ as little squares of
size $\frac{1}{N}$ then both the corresponding jeu de taquin path as well as the
evacuation path is a zig-zag line connecting the opposite corners of the unit square $[0,1]^2$.
We consider the case when $T_N\in \mathcal{T}_{\square_N}$ is a \emph{random} square tableau
sampled with the uniform distribution~$\mathbb{P}_N$.
\textbf{The~goal of the current
	paper is to find asymptotics of the corresponding scaled random zig-zag
	lines in the limit as $N\to\infty$, see \cref{fig:simulation-jdt}.}

\begin{figure}
\centering
\hfill
\subfloat[]{
\subfile{evacuation-figures/kwadrat-jdt-duzo-100.tex}

\label{fig:simulation-jdt}
}
\hfill
\subfloat[]{
	\begin{tikzpicture}[scale=0.45]
	\tiny
	\clip[](10,0) -- (10,10) -- (0,10) -- (0,0);
	\draw[black!20] (0,0) grid (10,10);
	\draw[Set1-3-3,ultra thick] plot[smooth,scale=10] file {circle025.txt};
	\draw[Set1-3-2,ultra thick] (0,10) -- (10,0);
	\draw[Set1-3-1,ultra thick] plot[smooth,scale=10] file {circle075.txt};

	\draw[line width=1mm,Set1-3-1,dashed,rounded corners](10,4) -- (8,4) -- (8,5) -- (7,5) -- (7,7) -- (6,7) -- (6,8) -- (5,8) -- (5,9) -- (2,9) -- (2,10);
	\draw[line width=1mm,Set1-3-2,dashed,rounded corners](10,2) -- (8,2) -- (8,3) -- (7,3) -- (7,4) -- (5,4) -- (5,5) -- (4,5) -- (4,6) -- (2,6) -- (2,8) -- (1,8) -- (1,10);
	\draw[line width=1mm,Set1-3-3,dashed,rounded corners](9,0) -- (9,1) -- (6,1) -- (6,2) -- (3,2) -- (3,4) -- (2,4) -- (2,5) -- (1,5) -- (1,7) -- (0,7);
	\begin{scope}
	\end{scope}
	\draw (0.5,0.5) node {$1$};
	\draw (1.5,0.5) node {$2$};
	\draw (2.5,0.5) node {$3$};
	\draw (0.5,1.5) node {$4$};
	\draw (1.5,1.5) node {$5$};
	\draw (3.5,0.5) node {$6$};
	\draw (0.5,2.5) node {$7$};
	\draw (4.5,0.5) node {$8$};
	\draw (5.5,0.5) node {$9$};
	\draw (0.5,3.5) node {$10$};
	\draw (1.5,2.5) node {$11$};
	\draw (2.5,1.5) node {$12$};
	\draw (2.5,2.5) node {$13$};
	\draw (1.5,3.5) node {$14$};
	\draw (2.5,3.5) node {$15$};
	\draw (3.5,1.5) node {$16$};
	\draw (0.5,4.5) node {$17$};
	\draw (4.5,1.5) node {$18$};
	\draw (5.5,1.5) node {$19$};
	\draw (1.5,4.5) node {$20$};
	\draw (6.5,0.5) node {$21$};
	\draw (7.5,0.5) node {$22$};
	\draw (0.5,5.5) node {$23$};
	\draw (8.5,0.5) node {$24$};
	\draw (0.5,6.5) node {$25$};
	\draw (0.5,7.5) node {$26$};
	\draw (3.5,2.5) node {$27$};
	\draw (4.5,2.5) node {$28$};
	\draw (0.5,8.5) node {$29$};
	\draw (3.5,3.5) node {$30$};
	\draw (6.5,1.5) node {$31$};
	\draw (2.5,4.5) node {$32$};
	\draw (4.5,3.5) node {$33$};
	\draw (5.5,2.5) node {$34$};
	\draw (6.5,2.5) node {$35$};
	\draw (1.5,5.5) node {$36$};
	\draw (1.5,6.5) node {$37$};
	\draw (3.5,4.5) node {$38$};
	\draw (2.5,5.5) node {$39$};
	\draw (3.5,5.5) node {$40$};
	\draw (5.5,3.5) node {$41$};
	\draw (1.5,7.5) node {$42$};
	\draw (7.5,1.5) node {$43$};
	\draw (0.5,9.5) node {$44$};
	\draw (8.5,1.5) node {$45$};
	\draw (9.5,0.5) node {$46$};
	\draw (6.5,3.5) node {$47$};
	\draw (4.5,4.5) node {$48$};
	\draw (7.5,2.5) node {$49$};
	\draw (9.5,1.5) node {$50$};
	\draw (4.5,5.5) node {$51$};
	\draw (8.5,2.5) node {$52$};
	\draw (2.5,6.5) node {$53$};
	\draw (9.5,2.5) node {$54$};
	\draw (5.5,4.5) node {$55$};
	\draw (2.5,7.5) node {$56$};
	\draw (3.5,6.5) node {$57$};
	\draw (7.5,3.5) node {$58$};
	\draw (5.5,5.5) node {$59$};
	\draw (4.5,6.5) node {$60$};
	\draw (6.5,4.5) node {$61$};
	\draw (8.5,3.5) node {$62$};
	\draw (1.5,8.5) node {$63$};
	\draw (5.5,6.5) node {$64$};
	\draw (3.5,7.5) node {$65$};
	\draw (4.5,7.5) node {$66$};
	\draw (2.5,8.5) node {$67$};
	\draw (3.5,8.5) node {$68$};
	\draw (1.5,9.5) node {$69$};
	\draw (4.5,8.5) node {$70$};
	\draw (5.5,7.5) node {$71$};
	\draw (6.5,5.5) node {$72$};
	\draw (9.5,3.5) node {$73$};
	\draw (6.5,6.5) node {$74$};
	\draw (7.5,4.5) node {$75$};
	\draw (7.5,5.5) node {$76$};
	\draw (8.5,4.5) node {$77$};
	\draw (7.5,6.5) node {$78$};
	\draw (8.5,5.5) node {$79$};
	\draw (5.5,8.5) node {$80$};
	\draw (6.5,7.5) node {$81$};
	\draw (2.5,9.5) node {$82$};
	\draw (9.5,4.5) node {$83$};
	\draw (8.5,6.5) node {$84$};
	\draw (7.5,7.5) node {$85$};
	\draw (6.5,8.5) node {$86$};
	\draw (3.5,9.5) node {$87$};
	\draw (9.5,5.5) node {$88$};
	\draw (9.5,6.5) node {$89$};
	\draw (4.5,9.5) node {$90$};
	\draw (7.5,8.5) node {$91$};
	\draw (8.5,7.5) node {$92$};
	\draw (5.5,9.5) node {$93$};
	\draw (9.5,7.5) node {$94$};
	\draw (8.5,8.5) node {$95$};
	\draw (6.5,9.5) node {$96$};
	\draw (7.5,9.5) node {$97$};
	\draw (8.5,9.5) node {$98$};
	\draw (9.5,8.5) node {$99$};
	\draw (9.5,9.5) node {$100$};

	\draw[ultra thick](0,0) -- (0,10) -- (10,10) -- (10,0) -- cycle ;
	\end{tikzpicture}
\label{fig:level-curves}
}
\hfill \

	\caption{ \protect\subref{fig:simulation-jdt} Zigzag lines are sample jeu de
    taquin paths for random square tableaux of size $N=100$,
    selected so they cross the anti-diagonal near the corresponding
    meridians (smooth thick curves) with the longitudes
    $\psi\in\{\nicefrac{1}{10},\nicefrac{2}{10},\dots,\nicefrac{9}{10}\}$.
    The grey lines are the meridians with the longitudes
    $\psi\in\{\nicefrac{2}{100},\nicefrac{4}{100},\dots,\nicefrac{98}{100}\}$. \protect\subref{fig:level-curves} Sample random square tableau of size $N=10$. Zigzag lines are the level curves for $\alpha\in\{ \nicefrac{1}{4} ,\nicefrac{1}{2}, \nicefrac{3}{4} \}$. The smooth lines are the corresponding circles of~latitude~$g_\alpha$.}
\end{figure}

\section{Asymptotics of evacuation paths}\label{sec:circles}

It is easy to change the definition of jeu de taquin so that
it becomes a permutation $J\colon \mathcal{T}_\lambda\to \mathcal{T}_\lambda$ on the set
of standard tableaux of prescribed shape, see \cref{fig:jdtC}.

We keep notations from \cref{sec:problems} and consider the evacuation path
\eqref{eq:evacuation} for $T=T_N$ and $n=N^2$. The position of each of
the boxes in the evacuation path coincides with the~position of a specific box
in the standard Young tableau obtained by iterating $J$:
\begin{equation}
\label{eq:happy-box}
\on{pos}_{N^2}\big( j^{i}(T_N)\big)  = \on{pos}_{N^2-i} \big( J^{i}(T_N)\big).
\end{equation}
The latter standard tableaux $J^i(T_N)$ is a
uniformly random square Young tableau.

Some light on the problem of finding the right-hand side of
\eqref{eq:happy-box} is provided by the~work of Pittel and~Romik; one of their results
\cite[Theorem~2]{Pittel2007} gives an explicit family of~curves~$(g_\alpha)$
indexed by $\alpha\in[0,1]$ which fit inside the square $[0,1]^2$, see
\cref{fig:level-curves}.
Each~curve~$g_\alpha$ turns out to be  one of the level
curves of (scaled down) random standard Young tableau of shape $\square_N$ which
separates the boxes with entries $\leq \alpha N^2$ from the~boxes with entries
bigger than this threshold
(asymptotically, as $N\to\infty$, except for an~event of negligible probability).
We will refer to the curves $(g_\alpha)$ as \textbf{\emph{circles of latitude}}.

It follows that \textbf{asymptotically the (scaled down) position of the box
    \eqref{eq:happy-box} in the evacuation path is very close to some point on the
    circle of latitude $g_{1-t}$ where $t = \nicefrac{i}{N^2}$. \emph{The
        remaining difficulty is to pinpoint a specific location of this point on the
        curve}. For~this purpose we need a convenient parametrization on each circle
    of~latitude~$g_\alpha$.}

\subsection{Longitude. Geographic coordinates on the square}

The work of Pittel and Romik \cite[Theorem~2]{Pittel2007} additionally gives the
limit distribution of~the~(scaled down) position of the~box~$\lfloor \alpha
N^2\rfloor$ in a random square tableau with $N^2$ entries.
This probability distribution is supported on the curve $g_\alpha$ and will be denoted by $\mu_\alpha$.
We will use $\mu_\alpha$ to construct \emph{geographical coordinate system} on the unit square $[0,1]^2$.

For any point $p \in [0,1]^2$ there is exactly one curve $g_\alpha$ on which it lies. We say that
\emph{the~latitude of $p$ is equal to $\alpha$}. The \textbf{\emph{longitude} of $p$}, denoted by $\psi(p)\in
[0,1]$, is defined as the~measure~$\mu_\alpha$ of the set of points on the curve
$g_\alpha$ which have their $u$-coordinate smaller than the~\mbox{$u$-coordinate} of~$p$.
For given $\alpha,\psi\in [0,1]$ we denote by $P_{\alpha, \psi}$
the unique point of~the~square~$[0,1]^2$ with the appropriate latitude and
longitude.

\subsection{The first main result: asymptotics of evacuation paths}

For a given tableau $T_N\in \mathcal{T}_{\square_N}$ and $t\in[0,1]$ we denote
by
\[X_t= X_t(T_N) = \frac{1}{N} \on{pos}_{N^2}\left( j^{\lfloor t N^2 \rfloor}(T_N) \right) \in [0,1]^2 \]
the scaled down point from the evacuation path, cf.~\eqref{eq:evacuation}.

Our first main result states that, asymptotically, \emph{the scaled evacuation
	path in a random square tableau is a random meridian}, i.e.~a curve which
consists of points in the square $[0,1]^2$ with equal longitude $\psi$.
The probability distribution of this longitude $\psi$ is the uniform
distribution on the interval $[0,1]$.

\begin{thm}
\label{mainthm}

For each $N\in \mathbb{N}$ there exists a map~$\Psi_N\colon \mathcal{T}_{\square_N} \to [0,1]$ such
that for each $\varepsilon > 0$
\begin{equation}
\label{eq:main-theorem}
\lim_{N \to \infty} \mathbb{P}_N \Big\{ T_N \in \mathcal{T}_{\square_N}: \sup_{t \in [0,1]}
			\left| X_t(T_N) -
			P_{1-t, \Psi_N(T_N)}  \right| > \varepsilon \Big\} = 0.
\end{equation}
The probability distribution of the random variable $\Psi_N$ converges,
as $N\to\infty$, to the uniform distribution on the unit interval $[0,1]$.
\end{thm}

\section{Asymptotics of jeu de taquin paths}

Let $T$ be a standard Young tableau with $n$ boxes. We define
$\mathbf{q}(T)=(\mathbf{q}_1,\dots,\mathbf{q}_n) \subset \mathbb{N}^2$ to be the corresponding \emph{jeu de taquin
	path in the lazy parametrization}. More specifically, $\mathbf{q}_i$~is~defined as~the~last box
	along the jeu de taquin path corresponding to $T$
(cf.~\cref{fig:jdtA}) which contains a number $\leq i$.

Our second main result states that, asymptotically,
\emph{the scaled jeu de taquin path in~a~random square tableau is a random meridian}.

\begin{thm}\label{jdtthm}
For each $N\in \mathbb{N}$ there exists a map~$\widetilde{\Psi}_N\colon
\mathcal{T}_{\square_N} \to [0,1]$ such that for each $\varepsilon > 0$
\[\lim_{N \to \infty} \mathbb{P}_N \Big\{ T_N \in \mathcal{T}_{\square_N}: \sup_{t \in [0,1]}
\left|
 \frac{1}{N} \mathbf{q}_{\lceil t N^2\rceil}(T_N)  -
P_{t,\widetilde{\Psi}(T_N)} \right| > \varepsilon \Big\} = 0. \]
The probability distribution of the random variable $\widetilde{\Psi}_N$
converges, as $N\to\infty$, to the uniform distribution on the unit interval
$[0,1]$.
\end{thm}

An analogous problem was studied by Romik and \'Sniady \cite{Romik2015a} for
random tableaux obtained by applying Robinson--Schensted correspondence to a
random permutation (\emph{``Plancherel measure''}). The proof
presented there is not applicable in our context.

\begin{proof}
We show equivalence of Theorem~\ref{mainthm} and Theorem~\ref{jdtthm} by
comparing positions of~boxes in lazy jeu de taquin path and the evacuation path
in pairs of corresponding tableaux.

In the proof we use:
1. a~\emph{Robinson-Schensted-Knuth algorithm} which is a bijection which applied to a permutation $\sigma = (\sigma_1, \ldots, \sigma_n) \in \mathfrak{S}_n$
returns a pair $(P(\sigma),Q(\sigma))$ of~tableaux of~the~same shape;
2. a~\emph{Sch\"utzenberger involution} which is a~map
\[\mathfrak{S}_n \ni \sigma \mapsto \varepsilon^*(\sigma):=(n+1-\sigma_n, \ldots, n+1-\sigma_1)\]
and 3. the following identity~\cite[Section~2]{Romik2015a} in which \mbox{$s(\sigma) = s(\sigma_1, \sigma_2, \ldots, \sigma_n) := (\sigma_2, \ldots, \sigma_n)$}
\begin{equation}\label{shift}
\left(j \circ Q\right)(\sigma) = \left(Q \circ \; s\right) (\sigma), \ \ \ \sigma = (\sigma_1, \ldots, \sigma_n) \in \mathfrak{S}_n
\end{equation}
(\eqref{shift} holds true up to renumbering of the boxes on the left-hand side.)

Let $\sigma \in \mathfrak{S}_{N^2}$. In the lazy jeu de taquin path $\mathbf{q}(T) := \mathbf{q}(Q(\sigma))$ we have by \eqref{shift}:
\[\mathbf{q}_n(T) = Q\left(\sigma_1, \sigma_2 \ldots, \sigma_{n} \right) \setminus Q\left(\sigma_2, \ldots, \sigma_{n} \right).\]
On the other hand, using \eqref{shift} many times, in the evacuation path
\eqref{eq:evacuation} in the tableau $T^* := Q(\varepsilon^*(\sigma))$,
the~position of~the~box with entry $N^2$ in $j^{N^2-p}(T^*)$ is
\[Q\left(N^2+1-\sigma_{p}, \ldots, N^2+1-\sigma_2, N^2+1-\sigma_{1}\right) \setminus Q\left(N^2+1-\sigma_{p}, \ldots, N^2+1-\sigma_{2}\right).\]
By Greene theorem, this position is equal to $Q\left(\sigma_1, \sigma_2 \ldots, \sigma_{p} \right) \setminus Q\left(\sigma_2, \ldots, \sigma_{p} \right) = \mathbf{q}_p(T)$.
 \end{proof}

\section{Sketch of proof of \cref{mainthm}}

\begin{proof}[Proof of \cref{mainthm}]
In order to investigate the point $X_t$
we will pass to its geographic coordinates $\alpha(X_t)$ and $\psi(X_t)$;
our goal is to show that for each $c>0$ and $\varepsilon > 0$
\begin{equation}
\label{eq:latitudelongitude}
\sup_{t \in [0,1]}  \big|
\alpha(X_t)- (1-t)
\big|  < \varepsilon
\qquad \text{and} \qquad
\sup_{t \in [c,1-c]}  \big|
\psi(X_t) - \Psi_N(T_N)
\big|  < \varepsilon
\end{equation}
hold true except for $T_N$ in a set which has asymptotically negligible
probability $o(1)$.

The discussion from the beginning of \cref{sec:circles} shows
that  the first statement in
\eqref{eq:latitudelongitude} indeed holds true
(except for an event of negligible
probability, for $N\to\infty$).

If the second statement in \eqref{eq:latitudelongitude} instead of a supremum
would involve only a fixed value of $t_0\in (0,1)$, we could simply define
$\Psi_N(T_N):=\psi(X_{t_0})$ to be longitude of $X_{t_0}$.
Unfortunately, we need to justify that this choice of $\Psi_N(T_N)$ is
also good for other values of $t\neq t_0$. This kind of result is
provided by \cref{prop:all-the-same}.
\end{proof}

\begin{proposition}
\label{prop:all-the-same}
Assume that $0<t_1 < t_2 <1$. Then  for each $\varepsilon>0$
\[ \lim_{N\to\infty}
\mathbb{P}_N \left\{ T_N \in\mathcal{T}_{\square_N} :
    \big| \psi(X_{t_1}) - \psi(X_{t_2}) \big| > \varepsilon \right\} = 0.\]

\end{proposition}

\section{Sketch of proof of \cref{prop:all-the-same}. Surfers on the sink}

\cref{prop:all-the-same} is equivalent to the conjuction of the following two statements for $\varepsilon>0$:
\begin{align}
\lim_{N\to\infty} \mathbb{P}_N \left\{ T_N \in\mathcal{T}_{\square_N} :
\psi(X_{t_1}) - \psi(X_{t_2})  > \varepsilon \right\} &= 0,
\label{eq:sink-A}
\\
\lim_{N\to\infty} \mathbb{P}_N \left\{ T_N \in\mathcal{T}_{\square_N}
:  \psi(X_{t_1}) - \psi(X_{t_2})  < -\varepsilon \right\} &=
0. \label{eq:sink-B}
\end{align}
We present the proof of \eqref{eq:sink-A}; the other claim
\eqref{eq:sink-B} can be proved in an analogous way.

\subsection{Single surfer scenario}
\label{sec:single}

We present the story of the evacuation path in a different light. We will speak
about a~\emph{square pool of side $N$} (=the square Young diagram $\square_N$)
filled with \emph{$N^2-1$ particles of~water} (=the Young tableau $T_N$ with
the largest entry removed), a passive \emph{surfer} (=the box with the biggest
entry) and its behaviour when the pool is drained (=jeu de taquin). Our~goal~in
\cref{mainthm} is to show that, when pool is big enough, the surfer has some
typical paths along which he/she moves.

\textbf{In our proof of \cref{prop:all-the-same} we start our analysis at time $t_1$
when jeu de taquin was already applied $m=\lfloor t_1 N^2 \rfloor$ times.}
Our starting point is therefore the tableau
\[ T'_N:=j^m(T_N)\]
with $N^2-m$ boxes. Let us renumber the entries of this tableau so it becomes
standard. We denote $w=N^2-m-1$. In this way the boxes with numbers $1,\dots,w$
correspond to~the~\emph{water} and the box with the number $w+1$ to the
\emph{surfer}. By removing the box with the surfer
\[ W_N':= T_N' \setminus \{w+1\}\]
we get a standard Young tableau which encodes \emph{the configuration of the water}.

\subsection{Pieri tableaux. Multisurfer scenario}
\label{sec:multi}

Let $k=\lfloor \sqrt[4]{N} \rfloor$.
Let $M$ be a tableau in which the $k$ largest entries are numbered
by~consecutive integers $l+1,\dots,l+k$. We say that $M$ is a \emph{Pieri
    tableau} if these $k$ largest boxes are in \emph{an increasing order} from
north-west to south-east:
\[ u^{M}_{l+1} < \cdots < u^{M}_{l+k}.\]
By $\widetilde{\mathcal{T}}_{\square_N}$ we denote the set of standard Young tableaux of
shape $\square_N$ which are Pieri, and~by~$\widetilde{\mathbb{P}}_N$~the~uniform distribution
on~the~set~$\widetilde{\mathcal{T}}_{\square_N}$.

It is easy to check that if $M$ has at least $k+1$ boxes then $M$ is~a~Pieri tableau
if and only if $j(M)$ is a Pieri tableau.

\medskip

To pinpoint the position of the surfer in random water, we will need a point of~reference.
We introduce \textbf{the multisurfer scenario} in which
we consider the square pool filled with $N^2-k$ \emph{particles of water} on
which are surfing $k$ \emph{surfers} (=$k$~boxes with the biggest entries). We
assume that the multisurfers are in an increasing order; it~follows that this
scenario corresponds to a Pieri tableau $M_N\in
\widetilde{\mathcal{T}}_{\square_N}$ We assume that $M_N$ is a random
tableau sampled with the probability distribution $\widetilde{\mathbb{P}}_N$.

\textbf{We start our analysis when jeu de taquin was already applied $m+1-k$ times.}
Our~starting point is therefore the tableau
\[ M'_N:=j^{m+1-k}(M_N)\]
with $N^2+k-m-1=w+k$ boxes. Let us renumber the entries of this tableau so
it~becomes  standard. In this way, just as in the single surfer scenario, the
boxes with numbers $1,\dots,w$ correspond to the \emph{water}. On the
other hand, the boxes with the numbers $w+1,\dots,w+k$ correspond to
the \emph{multisurfers}. By removing the multisurfers
\begin{equation}
\label{eq:water} \widetilde{W}_N':= M_N' \setminus
\{w+1,\dots,w+k\} = J^{m+1-k} ( M_n ) \setminus
\{w+1,\dots, N^2\}
\end{equation}
we get a standard Young tableau which encodes the configuration of the water.

\subsection{Counting multisurfers gives longitude}
\label{multiplesketch-0}

For $0\leq i \leq N^2-k$ we consider the situation after jeu de taquin was
applied $i$ times in~the~multisurfer scenario, i.e.~$M'_N =
j^i(M_N)$, and for $u\in [-1,1]$ we define the random variable
$\overline{\Psi}_i(u)$ to~be~the~fraction of the number of multisurfers
which have scaled \mbox{$u$-coordinates} smaller than~$u$,~i.e.:~
\begin{equation}
\label{eq:psi-via-multi}
	\overline{\Psi}_i(u) := \frac{1}{k} \max\left\{ p\in\{1,\dots,k\}  :
	\frac{1}{N} u^{M_N'}_{w+p} \leq u \right\}.  \end{equation}

The following result states that $\overline{\Psi}_i$ gives a good
approximation of the longitude.

\begin{proposition}
\label{prop:longitude-experimental}
Let $\alpha\in(0,1)$ be fixed.
For a given $\psi\in[0,1]$ we set  $(x,y)=P_{\alpha,\psi}$ to
be the~Cartesian coordinates of the point with geographic coordinates
$\alpha,\psi$ and define $U_{\alpha,\psi}:=x-y$.

We set $i:=\lfloor (1-\alpha) N^2\rfloor$. The random variable
$\overline{\Psi}_{i}\left( U_{\alpha,\psi} \right)$ converges in
probability to $\psi$ and~this convergence is uniform over $\psi$; in
other words for each $\varepsilon>0$
\[
\widetilde{\mathbb{P}}_N\Big\{
M_N\in\widetilde{\mathcal{T}}_{\square_N} :
\sup_{\psi\in [0,1]}
\left| \psi - \overline{\Psi}_{i} \big( U_{\alpha,\psi} \big) \right| > \varepsilon
	   \Big\}  = O\left( \frac{1}{\sqrt{N}} \right). \]
\end{proposition}
\begin{proof}
	The scaled $u$-coordinates of the multisurfers may be encoded by a random
probability measure on $\mathbb{R}$
	\[ \mu := \frac{1}{k} \sum_{1\leq p\leq k} \delta\left(\frac{1}{N}\ u_{w+p}^{M_N'}  \right),  \]
where $\delta(x)$ denotes the point measure concentrated in $x$. Our goal is to
prove that the cumulative probability function of $\mu$ converges uniformly to
some explicit limit; the~convergence should hold true in probability.

In order to achieve this goal we investigate the \emph{moments} of the measure $\mu$:
\[ \mathcal{M}_r = \int_{\mathbb{R}} z^r \dif \mu =
\frac{1}{k} \sum_{1\leq p\leq k} \left(\frac{1}{N}\ u_{w+p}^{M_N'}  \right)^r \qquad \text{for } r\in\{1,2,\dots\}.
\]
By applying \cref{expvalue} for some specific symmetric polynomials $W$ we are
able to find explicit values for the mean value and the variance of the random
variable $\mathcal{M}_r$. It turns out that, asymptotically for $N\to\infty$,
the mean value $\mathbb{E} \mathcal{M}_r$ converges to the `right' value and the
variance $\operatorname{Var} \mathcal{M}_r$ tends to zero which implies that
the convergence of measures \emph{in the sense of moments} indeed holds true in
probability.

Since the limit measure is compatly supported, the convergence
\emph{in the sense of moments} implies also convergence of measures \emph{in
	the weak topology of probability measures}, as~required.
\end{proof}

\textbf{\cref{expvalue} provides a link between the statistical properties of
	the multisurfers and~the~representation theory of the symmetric groups.} The
problem of understanding (the~$u$-coordinates of) the positions of the
multisurfers after $i$ applications of jeu de taquin is equivalent to finding
(the  $u$-coordinates~of) the boxes with numbers from the~set~$M:=\{a,\dots,b\}:=\{N^2+1-k-i,\dots,N^2-i\}$
in the random tableau $S:=J^i(M_N)$.

We will view the symmetric group $\mathfrak{S}_k\subset \mathfrak{S}_{N^2}$ as the group of
permutations of~the~set~$M$ and define an element of the symmetric group
algebra \mbox{$p_{\mathfrak{S}_k} = \frac{1}{k!} \sum_{\sigma \in \mathfrak{S}_k} \sigma \in
\mathbb{C} \mathfrak{S}_{N^2}$}. We~also consider Jucys--Murphy elements $Z_s :=
\sum_{i=1}^{s-1} (i, s)\in\mathbb{C} \mathfrak{S}_{N^2}$   (see~\cite{Jucys1974}).

\begin{lem}\label{expvalue}
	Let $W(x_1, \dots, x_k)$ be a symmetric polynomial in $k$ variables. Let $S$
be a random element (sampled with the uniform distribution) of the set of
tableaux $T\in\mathcal{T}_{\square_N}$ with the additional property that
\[ u_{a}^T < \cdots < u_{b}^T.\]

Then
	\begin{equation}\label{exp}
	\mathbb{E}\ W\left(u^S_{a},\ \dots,\ u^S_{b}\right) = \frac{
\chi^{\square_N}\!\big(p_{\mathfrak{S}_k} \cdot W \left( Z_{a},\ \dots,\ Z_{b}\right) \big)}{ \chi^{\square_N}\!\left(p_{\mathfrak{S}_k} \right) } .
	\end{equation}
\end{lem}
\begin{proof}
The proof has the following ingredients.

Firstly, the observation that the vector space $V^{\square_N}$
has a linear basis $(e_{T})$ indexed by~standard Young tableaux
$T\in\mathcal{T}_{\square_N}$ in which the action of Jucys--Murphy elements is~diagonal,
with the eigenvalue equal to the $u$-coordinate of the appropriate
box:
\[ Z_s\  e_T = u^{T}_{s}\ e_T.\]

Secondly, with the help of Littlewood--Richardson coefficients, we investigate
the~decomposition of~the~restricted representation $V^{\square_N}
\big\downarrow^{\mathfrak{S}_{N^2}}_{\mathfrak{S}_{N^2-i-k} \times \mathfrak{S}_{k} \times \mathfrak{S}_i}$
into irreducible components
\end{proof}

\subsection{Presence of multisurfers does not influence the water}
\label{podobienstwowody}

We will show that the probability distribution of water in the single surfer
scenario is~very close to its multisurfer  counterpart.

\begin{lem}
\label{probdiv}
For any standard Young tableau $S$ with shape $\lambda:=\on{sh} S$ and $w$ boxes
\begin{equation}
\label{Pp=Pt}
\frac{\mathbb{P}_N\Big\{ T_N\in\mathcal{T}_{\square_N} : W_N'= S \Big\} }%
{\widetilde{\mathbb{P}}_N\left\{ M_N\in\widetilde{\mathcal{T}}_{\square_N} : \widetilde{W}_N' = S \right\} }
=
\frac{\mathbb{P}_N\Big\{ T_N\in\mathcal{T}_{\square_N} : \on{sh} W_N'= \lambda \Big\} }%
{\widetilde{\mathbb{P}}_N\left\{ M_N\in\widetilde{\mathcal{T}}_{\square_N} : \on{sh} \widetilde{W}_N' = \lambda \right\} }
= 1 + O\left(\frac{1}{\sqrt{N}}\right).
\end{equation}
\end{lem}
\begin{proof}
The first equality follows from the observation that the conditional
probability distribution of $W_N'$ under the event that $\on{sh} W_N'=\lambda$ is
the uniform measure on $\mathcal{T}_\lambda$ and~the~analogous result for
$\widetilde{W}_N'$.

Since
$\widetilde{\mathbb{P}}_N(\cdot) = \mathbb{P}_N(\cdot | \widetilde{\mathcal{T}}_{\sq_N})$ is a conditional probability,
Bayes rule can be applied and~the~middle part of \eqref{Pp=Pt} is equal to
\begin{equation}
\label{eq:Bayes}
\frac{\mathbb{P}_N\left(\widetilde{\mathcal{T}}_{\sq_N}\right)}{\mathbb{P}_N\left(\widetilde{\mathcal{T}}_{\sq_N} \big| \on{sh} W_N' = \lambda \right) }.
\end{equation}

We will show the following estimate for the denominator:
\begin{equation}
\label{eq:cond}
\mathbb{P}_N\left(\widetilde{\mathcal{T}}_{\sq_N} \big| \on{sh} W_N' = \lambda \right) = \frac{1 + O \left( \frac{1}{\sqrt{N}} \right)}{k!}.
\end{equation}
Equation \eqref{eq:water} shows that under condition $\on{sh} W_N' = \lambda$ the
conditional probability distribution of $U:=J^{m+1-k}
(M_N ) \setminus \{1,\dots,w\}$ coincides with the uniformly
random tableau with the skew shape $\square_N \setminus \lambda$. The conditional
probability \eqref{eq:cond} is thus equal to the probability that the smallest
$k$ boxes of $U$ (or, equivalently, the smallest $k$ boxes of the rectified
tableau $V:=\operatorname{rect} U$) are ordered from north-west to south-east
\emph{\'a la Pieri}. We denote by~$\mu$ the~shape of $V$. It follows that
\begin{equation}
\label{eq:average}
 \mathbb{P}_N\left(\widetilde{\mathcal{T}}_{\sq_N} \big| \on{sh} W_N' = \lambda \right) =
\sum_{{\mu\in\mathbb{Y}_{N^2-w} \text{ such that } \mu \subset \square_N }}
\mathbb{P}_N\left\{ T_N \in \mathcal{T}_{\square_N} :V = \mu \right\} \cdot P_\mu,
\end{equation}
where $P_\mu$ denotes the probability that $V$ (viewed as a uniformly random
tableau sampled from $\mathcal{T}_\mu$) has entries $1,\dots,k$ ordered from
north-west to south-east or, equivalently,
the shape of the restricted tableau $\operatorname{sh} \left( V|_{1,\dots,k} \right) = (k)$
is the one-row Young diagram with~$k$~boxes.

The link between standard Young tableaux and combinatorics of  irreducible
representations of the symmetric groups implies that
\begin{multline}
\label{eq:ortho}
P_\mu = \frac{\text{multiplicity of $V^{(k)}$ in $\left( V^{\mu}\big\downarrow^{\mathfrak{S}_{N^2-w}}_{\mathfrak{S}_k}\right)$}}%
{\operatorname{dim} V^{\mu}} =
\frac{
	\left\langle \chi^{\mu} \big\downarrow^{\mathfrak{S}_{N^2-w}}_{\mathfrak{S}_k}, \chi^{(k)} \right\rangle
}{
 \operatorname{dim} V^{\mu}
}
	= \\  \frac{1}{k!} \left( 1 +
\sum_{\pi\in \mathfrak{S}_k \text{ and } \pi \neq \operatorname{id}}
\frac{ \chi^{\mu}(\pi) }{
			\operatorname{dim} V^{\mu}
		} \right),
\end{multline}
where the second equality follows from the orthogonality of irreducible characters.

Let $C>0$ be fixed. We say that a Young diagram $\mu$ is $C$-balanced
\label{balanced}
if both
its number of~rows and its number of columns are $\leq C\sqrt{|\lambda|}$.
Each tableau $\mu$ which contributes to~\eqref{eq:average} is $C$-balanced for
$C=\frac{N}{\sqrt{t_1 N^2}}=\frac{1}{\sqrt{t_1}}$.
The
result of Fer\'ay and \'Sniady \cite[Theorem~1]{Feray2011} gives a good upper
bound for the character ratios on the right-hand side of \eqref{eq:ortho} for
$\mu$ in~the~class of~$C$-balanced diagrams which implies that
\[ P_\mu = \frac{1 + O \left( \frac{1}{\sqrt{N}} \right)}{k!} .\]
The latter estimate combined with \eqref{eq:average} implies that \eqref{eq:cond} indeed holds true.

\medskip

An analogous (but simpler) reasoning shows that the numerator  from \eqref{eq:Bayes} fulfills
\begin{equation}
\label{eq:cond2}
\mathbb{P}_N\left(\widetilde{\mathcal{T}}_{\sq_N}\right) = \frac{1+O\left( \frac{1}{\sqrt{N}} \right)}{k!}.
\end{equation}

The estimates \eqref{eq:cond} and \eqref{eq:cond2} for the terms in \eqref{eq:Bayes} complete the proof.
\end{proof}

\subsection{Single- and multisurfer scenario on the same water}
\label{eq:single-multi}

\textbf{In \cref{sec:single,sec:multi} we considered two random tableaux: $T_N'$ and
$M_N'$ defined on~two~different probability spaces. Our goal is
to define variants of these two random tableaux $\mathbf{T}'_N$ and $\mathbf{M}'_N$ on
\emph{the same} probability space in such a way that the corresponding
configurations of water would coincide: $\mathbf{W}'_N:=W_N'= \widetilde{W}_N'$.}

In order to achieve this goal we start by randomly sampling the tableau $\mathbf{W}'_N$
(=\emph{water}) with the probability distribution of $W_N'$ from
\cref{sec:single}. %
Once $\mathbf{W}'_N$ is selected we randomly choose the tableaux $\mathbf{T}'_N$ (=\emph{single
	surfer}) and the Pieri tableau $\mathbf{M}'_N$ (=\emph{multisurfers}) according to
the conditional probability distributions (conditioned by choosing $\mathbf{W}'_N$)
\begin{align*}
\mathbb{P}(\mathbf{T}'_N=S)&:= \mathbb{P}_N\left( T_N'=S | W_N'=\mathbf{W}'_N \right), \\
\mathbb{P}(\mathbf{M}'_N=S )&:= \widetilde{\mathbb{P}}_N\left( M_N'=S | \widetilde{W}_N'=\mathbf{W}'_N  \right)
\end{align*}
for an arbitrary standard tableau $S$.

The probability distribution of $\mathbf{T}'_N$ coincides with the distribution of
$T'_N$ from \cref{sec:single}. On the other hand \cref{probdiv} shows that
the total variation distance between the distribution of $\mathbf{M}'_N$ and the
distribution of $M'_N$ from \cref{sec:multi} is small, of~order~$O\left( \frac{1}{\sqrt{N}} \right)$.

\subsection{Ghosts of multisurfers do not overtake the surfer from right to left}

Let us fix a pair of tableaux $\mathbf{T}'_N$ and $\mathbf{M}'_N$ given by the above
construction. For any $0\leq q \leq w$
we will iteratively apply jeu de taquin $q$ times to both tableaux and compare the outcomes.
(Informally speaking, with the notations of \cref{multiplesketch-0}, this means that jeu de taquin was applied $i=\lfloor t_1 N^2\rfloor+q$ times to the tableau $T_N$.)

A simple inductive argument shows that after $q$ steps
the configurations of water in~$j^q(\mathbf{T}'_N)$ and~$j^q(\mathbf{M}'_N)$ coincide:
\[ j^q(\mathbf{T}'_N) \setminus \{w+1\} = j^q(\mathbf{M}'_N) \setminus\{w+1,\dots,w+k\}.\]
\textbf{In other words, we are dealing with two alternative universes in which
	the history of~water is the same, and differ only by the number of surfers on
	top. Since the~dynamics of the surfer and the multisurfers depends only on the
	water, there is no interaction between the surfer and the
	multisurfers (they are like ghosts).}

\medskip

We define $\widetilde{\psi}_q$ as the fraction of the multisurfers which are to
the left of the surfer:
\[ \widetilde{\psi}_q := \frac{1}{k} \max\left\{ p\in\{1,\dots,k\}  :
 u^{j^q(\mathbf{M}'_N)}_{w+p} \leq u^{j^q(\mathbf{T}'_N)}_{w+1} \right\}.  \]
This quantity is a modification of \eqref{eq:psi-via-multi} in which the role
of~$j^i(M_N)$ is played by~$j^q(\mathbf{M}'_N)$ evaluated at $u=\frac{1}{N}
u^{j^q(\mathbf{T}'_N)}_{w+1}$. The discussion from \cref{eq:single-multi} shows
that for the conclusion of \cref{prop:longitude-experimental} remains valid
also for this modification. It follows that both at time $t_1$ as~well
as~at~time $t_2$ we have
\begin{align}
\label{eq:almostA}
\left| \psi(X_{t_1}) - \widetilde{\psi}_0 \right| &< \nicefrac{\varepsilon}{2}, \\
\label{eq:almostB}
\left| \psi(X_{t_2}) - \widetilde{\psi}_{\lfloor t_2 N^2\rfloor - \lfloor t_1 N^2\rfloor} \right| &< \nicefrac{\varepsilon}{2}.
\end{align}
hold true, except for an event of negligible probability.

By comparing the action of jeu de taquin on $j^q(\mathbf{M}'_N)$ with its action on $j^q(\mathbf{T}'_N)$ it~follows that
the sequence $(\widetilde{\psi}_q)$ is weakly decreasing:
\begin{equation}
\label{eq:decreasing}
\widetilde{\psi}_0 \geq \cdots \geq \widetilde{\psi}_{w};
\end{equation}
in other words, \textbf{if we compare the \emph{relative} positions of~the~surfer
and the ghosts of~the~multisurfers in the alternative universe, the
	multisurfers can only move from the left of the surfer to the right, but not in
	the opposite direction}.

By combining \eqref{eq:almostA}, \eqref{eq:almostB} and \eqref{eq:decreasing}
it follows that \eqref{eq:sink-A} indeed holds true, as we claimed. This
concludes the proof of \cref{prop:all-the-same}.

\printbibliography

\end{document}


\section*{Appendices}
In these appendices we provide supporting documents for some of the statements in the paper. Appendix \ref{refA} shows that the problem of $N$ multi-channel Kondo impurities immersed in 2D and 3D can be reduced to a set of one-dimensional problems, including the superconducting proximity pairing term. Appendix \ref{refB} summarizes our convention for the Klein factors that are used in bosonization.
\subsection{Dimensional reduction}\label{refA}
For completeness, here we show how a problem of $M$ impurities in two dimension can be reduced to a sum of 1D metallic systems terminated at the spin impurities, appropriate for field theory analysis. This is a generalization of $M=2$ case from Ref.\,\cite{Affleck95}. We start from a two-dimensional system described by $H=H_0+H_K+H_\Delta$ defined in Eqs.\,\ref{eq1} and \ref{eq2}, and generalize it to the case that there are $M$ Kondo impurities located at positions $\vec d_n$ coupled to the conduction band.

\subsubsection{Kondo coupling}
We consider $M$ impurities located on a 2D plane at positions $\vec d_n$ shown in Fig.\,3. We have
\begin{equation}
H_K=\int{\frac{d^2kd^2k'}{(2\pi)^4}}\psi^{\dagger}_{\vec k}\frac{\vec\sigma}{2}\psi^{\vphantom{\dagger}}_{k'}\cdot\sum_{n=1}^MJ_n\vec S_n e^{-i(\vec k-\vec k').\vec d_n}.
\end{equation}
We start by defining the fermions
\begin{equation}
\psi_{nE}=\int{\frac{d^2k}{(2\pi)^2}}\delta(\epsilon_k-E)e^{i\vec k\cdot\vec d_n}\psi^{\vphantom{\dagger}}_{\vec k},
\end{equation}
in terms of which the Kondo Hamiltonian is
\begin{equation}
H_K=\int{dEdE'}\sum_{n=1}^MJ_n\psi^{\dagger}_{nE}\frac{\vec\sigma}{2}\psi^{\vphantom{\dagger}}_{nE'}\cdot\vec S_n.
\end{equation}
The problem is that $\psi_{nE}$ do not obey standard anticommutation relations. Rather
\begin{equation}
\acom{\psi^{\vphantom{\dagger}}_{nE}}{\psi^{\dagger}_{mE'}}=g_{nm}(E)\delta(E-E'),
\end{equation}
where
\begin{equation}
g_{nm}(E)=\int{\frac{d^2k}{(2\pi)^2}}\delta(\epsilon_k-E)e^{i\vec k\cdot(\vec d_n-\vec d_m)}.
\end{equation}
Representing the vector $\vec d_{nm}\equiv \vec d_n-\vec d_m$ in polar coordinates by $\vec d_{nm}=d_{nm}(\cos\phi_{nm},\sin\phi_{nm})$ we can write
\begin{eqnarray}
g_{nm}(E)&=&\frac{k_E}{2\pi\partial_k\epsilon_E}\int_0^{2\pi}\frac{d\phi}{2\pi}e^{ik_Ed_{nm}\cos(\phi-\phi_{nm})}\nonumber\\
&=&\frac{k_E}{2\pi\partial_k\epsilon_E}\sum_pe^{-ip\phi_{nm}}J_p(k_Ed_{nm})i^p\int_0^{2\pi}\frac{d\phi}{2\pi}e^{ip\phi}\nonumber\\
&=&g_EJ_0(k_Ed_{nm}),
\end{eqnarray}
where
\begin{equation}
g_E=\frac{k_E}{2\pi\partial_k\epsilon_E}.
\end{equation}
Remarkably, $g_{nm}$ depends only on mutual distances of the impurities, measured by the corresponding wavelength $k_Ed_{nm}$. The matrix $g_{nm}$ is real and symmetric and it has
real eigenvalues $\lambda_n$ and can be diagonalized by orthogonal eigenvectors $\vec u_n$ where $g\vec u_n=\lambda_n\vec u_n$. So we can write $g_{nm}=\sum_p u_{np}\lambda_{p}u^*_{mp}$ and the orthogonality is $\sum_pu_{np}u^*_{mp}=\delta_{nm}$. Defining new operators with
\begin{equation}
\psi_{nE}=\sum_i u_{ni}(E)\sqrt{\lambda_i(E)}\tilde{\psi}_{iE},
\end{equation}
and
\begin{equation}
\tilde{\psi}_{iE}\equiv\frac{1}{\sqrt{\lambda_i(E)}}\sum_n{u^*_{ni}}(E)\psi^{\vphantom{\dagger}}_{nE}
\end{equation}
we find that they are orthonormal
\begin{equation} \acom{\tilde{\psi}_{iE}}{\tilde{\psi}_{jE'}}=\delta_{ij}\delta({E-E'}),
\end{equation}
and the Kondo Hamiltonian becomes
\begin{equation}
H_K=\int{dEdE'}\sum_{nij}J_{nij}(E,E')\vec S_n\cdot \tilde{\psi}^{\dagger}_{iE}\frac{\vec\sigma}{2}\tilde{\psi}^{\vphantom{\dagger}}_{jE'},
\end{equation}
where
\begin{equation}
J_{nij}(E,E')=J_n\sqrt{\lambda_i(E)\lambda_j(E')}u_{ni}(E)u^*_{nj}(E').
\end{equation}
Although the $J_{nij}$ couplings are complex in general, the relation $J^*_{nij}(E,E')=J_{nji}(E',E)$ ensures the hermiticity of the Hamiltonian. Whether or not a spin $n$ couples the channels $i$ and $j$ depend on the product of wavefunctions of $i$ and $j$ at the site $n$, which can be tuned by moving potential scatterings induced by scanning tips (Fig.\,3). Next, we do a Taylor expansion of $J_{nij}(E,E')$ function around Fermi energy and keep only the leading relevant term.
$H_K$ can be written in the matrix form $[\bb J_n]_{ij}=J_{nij}$,
\begin{equation}
H_K=\sum_n  \Psi^{\dagger}(0)[\bb J_n\vec S_n\cdot\frac{\vec \sigma}{2}]\Psi^{\vphantom{\dagger}}(0).
\end{equation}
where $\Psi$ is a vector in $m=1...M$ index and the spin is implicit. As we see the impurities talk to all the channels and scatter electron between all the channels. A result of the Taylor expansoin of $J_{nij}(E,E')$ is that RKKY interactions will be induced between the spins
\begin{equation}
H_K\to H_K+\sum_{ij}J_H^{ij}\vec S_i\cdot\vec S_j.
\end{equation}
\subsubsection{Mapping to 1D}
The mapping to (unfolded) left-movers follows Ref.\,\cite{Affleck95}
\begin{equation}
\psi_L(x,t)=\frac{1}{\sqrt v}\int_{-D}^D dE e^{-iE(t+x/v)}\psi_E,
\end{equation}
for $x\in(-\infty,+\infty)$
with commutation relations
\begin{equation}
\{\psi^{\dagger}_L(x),\psi^{\vphantom{\dagger}}_L(y)\}=2\pi\delta(x-y).
\end{equation}
Alternatively, we can work in the folded space of left and right-movers defined as
\begin{equation}
\psi_R(x,t)=\psi_L(-x,t), \qquad (x>0).
\end{equation}
\subsubsection{Kinetic term}
As we saw above, the Kondo interaction involves $M$ 1D conduction bands pulled out of the 2D conduction band. The natural guess for their kinetic energy  is
\begin{eqnarray}
H_0(\vec d)&\equiv&\int{dE E\psi^{\dagger}_{nE}\psi^{\vphantom{\dagger}}_{nE}}\\
&=&\int\frac{d^2kd^2k'}{(2\pi)^4}\epsilon_k\delta(\epsilon_k-\epsilon_k')e^{i(\vec k'-\vec k)\cdot\vec  d}\psi^{\dagger}_{\vec k}\psi^{\vphantom{\dagger}}_{\vec k'}.\nonumber
\end{eqnarray}
Summing over all the positions we have
\begin{eqnarray}
\int{d^2r}H_0(\vec r)
&=&\delta(0)\int\frac{d^2kd^2k'}{(2\pi)^2}\epsilon_k\psi^{\dagger}_{\vec k}\psi^{\vphantom{\dagger}}_{\vec k},
\end{eqnarray}
which is the total Hamiltonian up to a $\delta(0)$, i.e. the total volume. Transformation from $\psi_{nE}$ to $\tilde{\psi}_{iE}$ is a unitary transformation at each energy and it doesn't change the kinetic part. Therefore, we have
\begin{equation}
H_0\to\sum_{n=1}^M\int{dE E\tilde{\psi}^{\dagger}_{nE}\tilde{\psi}^{\vphantom{\dagger}}_{nE}}
\end{equation}
\subsubsection{Pairing term}
Without loss of generality, we consider a channel-diagonal s-wave singlet pairing
\begin{equation}
H_\Delta=\int{\frac{d^2k}{(2\pi)^2}}(\Delta\psi^{\dagger}_{k\uparrow}\psi^{\dagger}_{-k\downarrow}+h.c.).
\end{equation}
Let us look at
\begin{eqnarray}
H_\Delta^+(\vec d_n)&\equiv&\Delta\int{dE}\psi^{\dagger}_{nE\uparrow}\psi^{\dagger}_{n,E\downarrow}\nonumber\\
&&\hspace{-1.75cm}=\Delta\int{dE}\int{\frac{d^2kd^2k'}{(2\pi)^4}}\delta(\epsilon_k-E)\delta(\epsilon_{k'}-E)\psi^{\dagger}_{k\uparrow}\psi^{\dagger}_{q\downarrow}e^{-i(\vec k+\vec k')\cdot \vec d_n}\nonumber
\end{eqnarray}
Summing over all the points gives a $(2\pi)^2\delta(\vec k+\vec k')$ term that  gives
\begin{equation}
\int{d^2r}H^+_\Delta(\vec r)=\delta(0)\Delta\int{\frac{d^2k}{(2\pi)^2}}\psi^{\dagger}_{k\uparrow}\psi^{\dagger}_{-k\downarrow},
\end{equation}
which is our initial Hamiltonian up to the system volume. Therefore,
\begin{eqnarray}
H_{\Delta}&\to &\Delta\sum_{n=1}^M\int{dE}(\psi^{\dagger}_{nE\uparrow}\psi^{\dagger}_{nE\downarrow}+h.c.)\nonumber\\
&=&\Delta\sum_{ij}\int{dE}(\tilde{\psi}^{\dagger}_{iE\uparrow}\tilde{\psi}^{\dagger}_{jE\downarrow}+h.c.).
\end{eqnarray}
\subsection{Klein factors}\label{refB}
Here, we briefly review the treatment of Klein-factors following Ref.\,\cite{Zarand00}. Denoting the total number of fermions $\mu$ with
\begin{equation}
N_\mu=\intoinf{dx}c^{\dagger}_{\mu}(x)c^{\vphantom{\dagger}}_{\mu}(x), \quad \mu\in\{1\uparrow,1\downarrow,2\uparrow,2\downarrow\},
\end{equation}
the Klein factors have the property (no summation):
\begin{eqnarray}
&&F^{\vphantom{\dagger}}_\mu F^{\dagger}_\mu=F^{\dagger}_\mu F^{\vphantom{\dagger}}_\mu=1\\
&&\{F^{\vphantom{\dagger}}_\mu,F^{\dagger}_{\nu}\}=2\delta_{\mu\nu}\\
&&\{F^{\vphantom{\dagger}}_\mu,F^{\vphantom{\dagger}}_\nu\}=0, \quad \forall \mu\neq \nu\\
&&\com{F_\mu}{N_\nu}=\delta_{\mu\nu}F_\nu\\
&&\com{F_\mu}{\phi_\nu(x)}=\com{F_\mu}{\theta_\nu(x)}=0.
\end{eqnarray}
In the paper we transform into collective boson basis:
\begin{eqnarray}
\mat{{\cal N}_c \\ {\cal N}_s \\ {\cal N}_f \\ {\cal N}_{sf}}=\frac{1}{2}\matc{rrrr}
{1 & 1 & 1 & 1 \\
1 & -1 & 1 & -1 \\
1 & 1 & -1 & -1 \\
1 & -1 & -1 & 1}
\mat{N_{1\uparrow} \\ N_{1\downarrow} \\N_{2\uparrow} \\ N_{2\downarrow}}.
\end{eqnarray}
In accord with this transformation, we define new Klein factors ${\cal F}_\mu$ for $\mu=c,s,f,sf$ that satisfy
\begin{eqnarray}
&&{\cal F}^{\vphantom{\dagger}}_\mu {\cal F}^{\dagger}_\mu={\cal F}^{\dagger}_\mu {\cal F}^{\vphantom{\dagger}}_\mu=1\\
&&\{{\cal F}^{\vphantom{\dagger}}_\mu,{\cal F}^{\dagger}_{\nu}\}=2\delta_{\mu\nu}\\
&&\{{\cal F}^{\vphantom{\dagger}}_\mu,{\cal F}^{\vphantom{\dagger}}_\nu\}=0, \quad \forall \mu\neq \nu\\
&&\com{{\cal F}_\mu}{{\cal N}_\nu}=\delta_{\mu\nu}{\cal F}_\nu\\
&&\com{{\cal F}_\mu}{\phi_\nu(x)}=\com{{\cal F}_\mu}{\theta_\nu(x)}=0.
\end{eqnarray}
It can be seen that these any off-diagonal relation combination, like $F^{\dagger}_\mu F^{\vphantom{\dagger}}_\nu$ change the population of the collective modes ${\cal N}_\zeta$ by integers and thus can be represented by a combination of new Klein factors\,\cite{Zarand00}. For example,
\begin{equation}
F^{\dagger}_{1\uparrow}F^{\vphantom{\dagger}}_{1\downarrow}={\cal F}^{\dagger}_{sf}{\cal F}^{\dagger}_s, \qquad F^{\dagger}_{2\uparrow}F^{\vphantom{\dagger}}_{2\downarrow}={\cal F}_{sf}{\cal F}^{\dagger}_s.
\end{equation}
In the paper, we repeatedly use these and similar relations.


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

A large number of papers and textbooks are devoted to the nonlinear theory of
the Langmuir waves and Landau damping~\cite{landau46_vib_elec_plasma}.
In addition to the standard way of explanation of physics, a large number of diverse approaches and interpretations have been published~\cite{Weitzner_1963, Weiland_1981, Li_1991, Sagan_1994, Soshnikov_1996, Brodin_1997, Neil_1965, Shatashvili_1982}, which are more advanced than
the standard one.

In these papers, the solution of the linearized Vlasov equation is
used. After finding the zero-order solution of the main equation, the authors
construct approximations of any higher order~\cite{Tsintsadze_2009, Mazitov1965, bernstein57:_exact_nonlin_plasm_oscil, Drummond_2004, Drummond_2005}. In the present
paper, however a new approach is presented, allowing to analyze the problem
self-consistently in an arbitrary order of nonlinear approximation.

The exact solution for Vlasov equation is found for the first time. Using
this solution, it is shown that the waves with finite amplitude are not
exposed to damping. Only waves with small amplitudes, when the
oscillation frequency of captured (in the wave-well) electrons is smaller
than the damping rate, can damp~\cite{Kadomtzev_1978, Shafranov_1965}. It is found that on the
fulfillment of the definite resonance condition, the waves with finite
amplitudes are unstable. For this necessitates, the fulfillment of the
definite resonance conditions, which are similar to conditions with the
parametric resonance.

\section{\label{sec:2}Exact solution of the Vlasov equation}

We start as usual from the collisionless Vlasov equation and Poisson's
equation,

\begin{eqnarray}
\label{eq1}
\frac{df}{dt}&=\frac{\partial f}{\partial t}+v\frac{\partial f}{\partial
z}-E\left( {z,t} \right)\frac{\partial f}{\partial v}=0,\\
\label{eq2}
\frac{\partial E}{\partial z}&=1-n,
\quad
E=-\frac{\partial \phi }{\partial z},
\end{eqnarray}

where the dimensionless values for the time, the coordinate, the velocity
and the electric potential are used

\begin{equation}
\label{eq3}
\left( {\omega_{pe} t} \right)\to t, \left( {z/\lambda_{De} } \right)\to z,
\left( {v/v_{Te} } \right)\to v, \left( {e\Phi /T_{e} } \right)\to \phi,
\end{equation}

$\omega_{pe} $ is the electron plasma frequency, $\lambda_{De} $ is the
electron Debye length and $v_{Te} $ is the electron thermal velocity. The
dimensionless electric field $E$ and the electron number density $n$ are
defined as follows
\begin{equation}
\label{eq4}
\frac{E}{\sqrt {4\pi \cdot n_{0} T_{e} } }\to E,
\quad
\frac{n}{n_{0} }\to n=\int {dv\cdot f\left( {z,t,v} \right)} .
\end{equation}
Here $T_{e} $ is the electron temperature in the energetic units. It is
assumed that ions stay in the equilibrium with the density $n_{0} $, which
results in the first term ("1") on the right-hand side (rhs) of
Poisson's Eq.~\eqref{eq2}.

The characteristic equations for Eq.~\eqref{eq1} reads
\begin{equation}
\label{eq5}
\frac{dz}{dt}=v,
\quad
\frac{dv}{dt}=-E\left( {z,t} \right).
\end{equation}

For constants of integrals $R$ and $U$ ($dR/dt=0$, $dU/dt=0$ ) we find

\begin{eqnarray}
\label{eq6}
R&=z-\int\limits^{t} {dt'\cdot {H}' \left[ {v+\int\limits_{{t}'}^t
{d{t}''\cdot E\left\{ {z\left( {{t}''} \right),{t}''} \right\}} } \right]},
\\
\label{eq7}
U&=v+\int\limits^t {d{t}'\cdot E\left\{ {z\left( {{t}'} \right),{t}'}
\right\}},
\end{eqnarray}

where the functions $z\left( {{t}'} \right)$ are defined
with the expressions:

\begin{eqnarray}
\label{eq8}
z\left( {{t}'} \right)&=z-\int\limits_{{t}'}^t {d{t}''\cdot {H}'\left[
{v+\int\limits_{{t}''}^t {d{t}'''\cdot E\left\{ {z\left( {{t}'''}
\right),{t}'''} \right\}} } \right]},
\\
\label{eq9}
z\left( {{t}''} \right)&=z-\int\limits_{{t}''}^t {d{t}'''\cdot {H}'\left[
{v+\int\limits_{{t}'''}^t {dt^{IV}\cdot E\left\{ {z\left( {t^{IV}}
\right),t^{IV}} \right\}} } \right]} ,
\\
\label{eq10}
z\left( {{t}'''} \right)&=z-\int\limits_{{t}'''}^t {dt^{IV}\cdot {H}'\left[
{v+\int\limits_{t^{IV}}^t {dt^{V}\cdot E\left\{ {z\left( {t^{V}}
\right),t^{V}} \right\}} } \right]}, \\
&\mathellipsis . \nonumber
\end{eqnarray}

The chain~\eqref{eq8}, \eqref{eq9}, \eqref{eq10}, \textellipsis, can be continued. Here the
values $z\left( {{t}'} \right)$, $z\left( {{t}''} \right)$, $z\left( {{t}'''}
\right)$, \textellipsis, must be substituted into the arguments of the
electric fields' expressions $E\left\{ {z\left( {{t}'} \right),{t}'}
\right\}$, $E\left\{ {z\left( {{t}''} \right),{t}''} \right\}$, $E\left\{
{z\left( {{t}'''} \right),{t}'''} \right\}$, \textellipsis, etc.
In Eqs.~\eqref{eq6}-\eqref{eq8} the function $H\left( x \right)$ is defined as follows

\begin{equation}
\label{eq11}
H\left( x \right)=\frac{1}{2}x^{2},
\quad
{H}'\left( x \right)=x.
\end{equation}

In the following the upper dashes in ${H}'\left( x \right)$ will denote the
derivative with the whole argument of the function $H\left( x \right)$ and
in ${E}'\left\{ {z\left( t \right),t} \right\}$ will denote a derivative only
with respect to $z\left( t \right), \quad {E}'\left\{ {z\left( t \right),t}
\right\}=\partial E\left\{ {z\left( t \right),t} \right\}/\partial
z\left( t \right)$. The solution of Eq.~\eqref{eq1} can be represented in the form

\begin{equation}
\label{eq12}
f=f\left\{ {R\left( {v,z,t} \right),U\left( {v,z,t} \right)} \right\}.
\end{equation}
Substituting Eq.~\eqref{eq12} into Eq.~\eqref{eq1}, using the definitions~\eqref{eq6} and~\eqref{eq7}, and
successively carrying out the derivatives we obtain

\begin{equation}
\label{eq13}
\eqalign{
\frac{df}{dt}=\frac{\partial f}{\partial U}\cdot
\frac{dU}{dt}+\frac{\partial f}{\partial R}\cdot \frac{dR}{dt}= \cr
=\left\{ {\frac{\partial f}{\partial U}\cdot \int\limits^t {d{t}'\cdot
{E}'\left\{ {z\left( {{t}'} \right),{t}'} \right\}+\frac{\partial
f}{\partial R}} } \right\}\cdot \left( {-{H}'\left[ v \right]+v} \right)=0.}
\end{equation}

In Eq.~\eqref{eq13} all other terms cancel each other. According to the second
relation from~\eqref{eq11} ${H}'\left( v \right)=v$ and hence the solution ~\eqref{eq12} satisfies the kinetic equation~\eqref{eq1}.
Defining the initial distribution function we assume that at the initial
moment $t_{0} \to -\infty $ the electron distribution function depends only
on the velocity

\begin{equation}
\label{eq14}
f\left\{ {R\left( {v,z,t_{0} } \right),U\left( {v,z,t_{0} } \right)}
\right\}=\left. {f_{0} \left( v \right)} \right|_{t_{0\to -\infty } } .
\end{equation}

As $f_{0} \left( v \right)$ we can choose the Maxwell distribution function
with a normalizing coefficient, $\left( {1/\sqrt {2\pi } } \right)\cdot \exp
\left( {-v^{2}/2} \right)$. From the relations at the initial time $t_{0} $,

\begin{equation}
\label{eq15}
R\left( {v,z,t_{0} } \right)=R \quad \textrm{and} \quad U\left( {v,z,t_{0} } \right)=U
\end{equation}

and definitions given by Eqs.~\eqref{eq6} and \eqref{eq7} we can find expressions for the velocity $v$ and the coordinate $z$:

\begin{eqnarray}
\label{eq16}
v&=U-\int\limits^{t_{0} } {d{t}'\cdot E\left\{ {z\left( {{t}'} \right),{t}'}
\right\}},
\\
\label{eq17}
z&=R+\int\limits^{t_{0} } {d{t}'\cdot {H}'} \left[
{v+\int\limits_{{t}'}^{t_{0} } {d{t}''\cdot E\left\{ {z\left( {{t}''}
\right),{t}''} \right\}} } \right].
\end{eqnarray}

Into Eqs.~\eqref{eq16} and \eqref{eq17} the explicit expressions \eqref{eq6} and \eqref{eq7} can be again substituted. Hence using the initial condition \eqref{eq14} for the electron distribution function we find

\begin{equation}
\label{eq18}
f_{0} =f_{0} \left\{ {v+\int\limits_{t_{0} \to \infty }^t {d{t}'\cdot
E\left\{ {z\left( {{t}'} \right),{t}'} \right\}} } \right\},
\end{equation}

where values $z\left( {{t}'} \right)$, $z\left( {{t}'''} \right)$,
\textellipsis, are defined by Eqs.~\eqref{eq8},~\eqref{eq9}, \textellipsis. We can represent Poisson's Eq.~\eqref{eq2}  in the following form,

\begin{equation}
\label{eq19}
\frac{\partial E\left( {z,t} \right)}{\partial z}=1-\int\limits_{-\infty
}^\infty {dv\cdot f_{0} \left\{ {v+\int\limits_{t_{0} \to -\infty }^t
{d{t}'\cdot E\left\{ {z\left( {{t}'} \right),{t}'} \right\}} } \right\}} .
\end{equation}

It is convenient to introduce a new variable - s, defined by the relation:

\begin{equation}
\label{eq20}
v+\int\limits_{-\infty }^t {d{t}'\cdot E\left\{ {z\left( {{t}'}
\right),{t}'} \right\}} =s,
\end{equation}

which allows to simplify the argument of $f_{0} $. Then for the explicit
expression for the velocity $v$ we have

\begin{equation}
\label{eq21}
\eqalign {v=s-\int\limits_{t_{0} \to -\infty }^t {d{t}'\cdot E\left[ {z\left( {{t}',s}
\right) } \right] } = s-\int\limits_{t_{0\to .\infty } }^t {d{t}'\cdot E \left[ {z-s\left( {t-{t}'} \right)}+ \right.}\nonumber \cr
+\int\limits_{{t}'}^t {d{t}''\int\limits_{t_{0} }^{{t}''} {d{t}'''\cdot E\left\{ {z-s\left( {t-{t}''} \right)} \right.} } + \nonumber \cr
+\int\limits_{{t}'''}^t {dt^{IV}\int\limits_{t_{0} }^{t^{IV}} {dt^{V}\cdot
E\left\langle {z-s\left( {t-t^{V}} \right)+....} \right\rangle ,} } \left.
{{t}'''} \right\} \left. {{t}''} \right],}
\end{equation}

where by analogy to Eqs.~\eqref{eq8}, \eqref{eq9}, \textellipsis, we have introduced the relations:

\begin{equation}
\label{eq22}
z\left( {{t}',s} \right)=z-s\left\{ {t-{t}'} \right\}+\int\limits_{{t}'}^t
{d{t}''\int\limits_{t_{0} }^{{t}''} {d{t}'''\cdot E\left[ {z\left(
{{t}''',s} \right),{t}'''} \right]} },
\end{equation}

\begin{equation}
\label{eq23}
\eqalign{
z\left( {{t}''',s} \right)=z-s\left\{ {t-{t}'''}
\right\}+\int\limits_{{t}'''}^t {dt^{IV}\int\limits_{t_{0} }^{t^{IV}}
{dt^{V}\cdot E\left[ {z\left( {t^{V},s} \right),t^{V}} \right]} }\cr
\mathellipsis}
\end{equation}

Poisson's equation can be represented in the form

\begin{eqnarray}
\label{eq24}
\frac{\partial E\left( {z,t} \right)}{\partial z}=1-\int\limits_{-\infty
}^\infty {ds\cdot \frac{dv\left( s \right)}{ds}f_{0} \left( s \right)}
=1+\int\limits_{-\infty }^\infty {ds\cdot \frac{df_{0} \left( s
\right)}{ds}\cdot v\left( s \right)}.
\end{eqnarray}

Substituting Eq.~\eqref{eq21} into Eq.~\eqref{eq24} for the Maxwell
distribution function $f_{0} \left( s \right)$ the first term in~\eqref{eq24} is
compensated by the positive charge of ions:

\begin{equation}
\label{eq25}
\frac{\partial E}{\partial z}=-\int\limits_{-\infty }^{+\infty } {ds\cdot
\frac{\partial f_{0} \left( s \right)}{\partial s}\cdot \int\limits_{t_{0}
\to -\infty }^t {d{t}'\cdot E\left[ {z\left( {{t}'s} \right),{t}'} \right]}
} .
\end{equation}

Taking the first derivative with time and restricting ourselves to keeping
the terms up to second powers of the electric strength (such a restriction
will keep everywhere throughout the calculations) we find

\begin{equation}
\label{eq26}
\eqalign{
\frac{\partial }{\partial t}\frac{\partial E}{\partial
z}=-\int\limits_{-\infty }^{+\infty } {ds} \cdot \frac{\partial f_{0} \left(
s \right)}{\partial s}\left\{ {-s\int\limits_{t_{0} \to -\infty }^t
{d{t}'\cdot {E}'\left[ {z\left( {{t}',s} \right),{t}'} \right]} } \right. \cr
+ \int\limits_{t_{0} \to -\infty }^t {d{t}'\cdot {E}'\left[ {z\left( {{t}',s}
\right),{t}'} \right]}
\int\limits_{t_{0} \to -\infty }^t {d{t}'''\cdot E\left[ {z\left( {{t}''',s}
\right),{t}'''} \right]}
\cr
{-s\int\limits_{t_{0} \to -\infty }^t {d{t}'\cdot {E}'\left[ {z\left(
{{t}',s} \right),{t}'} \right]}} \cdot \cr
\left.\cdot\int\limits_{{t}'}^t
{d{t}''\int\limits_{t_{0} \to -\infty }^{{t}''} {d{t}'''\cdot {E}'\left[
{z\left( {{t}''',s} \right),{t}'''} \right]} }   \right\}.}
\end{equation}

In the last term in the rhs of this equation we can change the ordering of
integrals in the following manner

\begin{equation}
\label{eq27}
\int\limits_{t_{0} \to -\infty }^t {d{t}'\int\limits_{{t}'}^t {d{t}''\cdot
Q\left( {{t}',{t}''} \right)=\int\limits_{t_{0} \to -\infty }^t
{d{t}''\int\limits_{t_{0} \to -\infty }^{{t}''} {d{t}'} } } } \cdot Q\left(
{{t}',{t}''} \right).
\end{equation}

Then the second derivative with time of Eq.~\eqref{eq26} gives

\begin{equation}
\label{eq28}
\eqalign{
\frac{\partial^{2}}{\partial t^{2}}\frac{\partial E}{\partial
z}=-\frac{\partial E}{\partial z}-\int\limits_{-\infty }^{+\infty } {ds}
\cdot \frac{\partial f_{0} \left( s \right)}{\partial s}
\cdot \left\{ {s^{2}\cdot
\int\limits_{t_{0} \to -\infty }^t {d{t}'\cdot {E}''\left[ {z\left( {{t}',s}
\right),{t}'} \right]} } -\right.\cr
-s\cdot \frac{\partial }{\partial z}\cdot \int\limits_{t_{0} \to -\infty }^t
{d{t}'\cdot {E}'\left\{ {z-s\left( {t-{t}'} \right),{t}'} \right\}} \cr
\cdot \int\limits_{t_{0} \to -\infty }^t {d{t}'''\cdot E\left\{ {z-s\left(
{t-{t}'''} \right),{t}'''} \right\}}\cr
+\frac{\partial }{\partial t}\cdot \int\limits_{t_{0} \to -\infty }^t
{d{t}'\cdot {E}'\left\{ {z-s\left( {t-{t}'} \right),{t}'} \right\}}
\cdot \left.\int\limits_{t_{0} \to -\infty }^t {d{t}'''\cdot E\left\{ {z-s\left(
{t-{t}'''} \right),{t}'''} \right\}}\right\}.}
\end{equation}

In Eq.~\eqref{eq28} we can transform the last term (with the first derivative with
time) as follows

\begin{equation}
\label{eq29}
\eqalign{
\frac{\partial }{\partial t}\cdot \int\limits_{t_{0} \to -\infty }^t
{d{t}'\cdot {E}'\left[ {z\left( {{t}',s} \right),{t}'} \right]}
\int\limits_{t_{0} \to -\infty }^t {d{t}'''\cdot E\left[ {z\left( {{t}''',s}
\right),{t}'''} \right]} \cong
\cr
\cong \frac{\partial }{\partial z}\cdot E\left( {z,t}
\right)\int\limits_{t_{0} \to -\infty }^t {d{t}'\cdot E\left\{ {z-s\left(
{t-{t}'} \right),{t}'} \right\}}-s\cdot \frac{\partial }{\partial z}\cdot
\cr
\cdot\int\limits_{t_{0} \to -\infty }^t
{d{t}'\cdot E\left\{ {z-s\left( {t-{t}''} \right),{t}'} \right\}}\cdot \cr
\cdot\int\limits_{t_{0} \to -\infty }^t {d{t}'''\cdot E} \left\{ {z-s\left(
{t-{t}'''} \right),{t}'''} \right\} ,}
\end{equation}

In the first term on the right-hand side of Eq.~\eqref{eq29}, containing the squared
electric field, we can, following Eq.~\eqref{eq25}, use a simplified linearized
expression

\begin{equation}
\label{eq30}
\frac{\partial E}{\partial z}\cong -\int\limits_{-\infty }^{+\infty }
{ds\cdot \frac{\partial f_{0} \left( s \right)}{\partial
s}\int\limits_{t_{0} \to -\infty }^t {d{t}'\cdot E\left\{ {z-s\left(
{t-{t}'} \right),{t}'} \right\}} } .
\end{equation}

Using the relations \eqref{eq29} and \eqref{eq30} the fourth derivative with time of Eq.~\eqref{eq28} we can represent in the form

\begin{equation}
\label{eq31}
\eqalign{
\frac{\partial^{4}E}{\partial t^{4}}+\frac{\partial^{2}E}{\partial
t^{2}}+3\frac{\partial^{2}E}{\partial z^{2}}+\frac{\partial^{3}}{\partial
z^{3}}\cdot \cr
\cdot\int\limits_{-\infty }^{+\infty } {ds} \cdot s^{4}\cdot
\frac{\partial f_{0} \left( s \right)}{\partial s}\int\limits_{t_{0} \to
-\infty }^t {d{t}'\cdot E\left\{ {z-s\left( {t-{t}'} \right),{t}'} \right\}+}
\cr
+\frac{\partial }{\partial z}\frac{\partial^{2}}{\partial t^{2}}\left\{
{\int\limits_{t_{0} \to -\infty }^t {d{t}'\cdot E\left( {z,{t}'}
\right)}}\cdot \right. \cr
\left.\cdot \int\limits_{t_{0} \to -\infty }^t {d{t}'\cdot E\left( {z,{t}'}
\right)}  -\frac{1}{2}\cdot E\left( {z,t} \right)\cdot E\left( {z,t}
\right) \right\}=0.}
\end{equation}

On obtaining Eq.~\eqref{eq31}, the first derivative $\left( {\partial /\partial z}
\right)$ in front of every term has been canceled and the explicit
expressions for the integrals from the Maxwell distribution function
$-f_{0} \left( s \right)$ is used.

\section{\label{sec:level3}Waves in the weak nonlinear case}

The first three terms of Eq.~\eqref{eq31} should describe Landau damping
in the linear approximation. By means of the Fourier expansion of these three terms for the
electric amplitudes we find

\begin{equation}
\label{eq32}
E\left( {\omega ,k} \right)=\frac{1}{2\pi }\int\limits_{-\infty }^{+\infty }
{dt\int\limits_{-\infty }^{+\infty } {dz} \cdot E\left( {z,t} \right)\cdot
\exp \left\{ {i\left( {\omega t-kz} \right)} \right\}} ,
\end{equation}

which results in the following dispersion relation

\begin{equation}
\label{eq33}
\omega^{4}-\omega^{2}-3k^{2}+\int\limits_{-\infty }^{+\infty } {ds\cdot
\frac{k^{3}s^{3}}{\omega -ks}\cdot s\frac{\partial f_{0} \left( s
\right)}{\partial s}} =0.
\end{equation}

A, simple transformation of the last equation gives

\begin{equation}
\label{eq34}
\omega^{4}+\omega^{3}\int\limits_{-\infty }^{+\infty } {ds\frac{s}{\omega
-ks}\frac{\partial f_{0} \left( s \right)}{\partial s}=0} .
\end{equation}

After expanding the denominator of \eqref{eq34} in powers of $\left( {ks/\omega }
\right)$ for the real and the imaginary parts of the frequency, $\omega
=\omega_{0} +i\gamma $ , we obtain the equalities

\begin{equation}
\label{eq35}
\omega_{0}^{2} =1+3k^{2} \textrm{and} \gamma =-\sqrt {\frac{\pi }{8}}
\frac{1}{k^{3}}\exp \left\{ {-\frac{1}{2k^{2}}-\frac{3}{2}} \right\},
\end{equation}

which determine the frequency and the Landau damping of the high-frequency
Langmuir waves~\cite{alexandrov84:_princ_plasm_elect}.

In fact Landau damping describes the initial stage of the electrons
capturing by the wave cavity (the minimum region of the wave), herewith
the amplitude must be very small, smaller than the value proportional to
$\gamma $, namely$\sqrt {eE_{0} /T_{e} \cdot k} \ll\gamma /k$~\cite{Kadomtzev_1978} (here
$E_{0} $ is the amplitude of the wave). In other words, this inequality means
that the oscillation frequency of captured electrons in the wave-well, must
be much smaller than the damping rate of the wave $-\gamma $; the electrons
are pushed by the back-side wall of the wave-cavity and during the
time-interval of passing the cavity width, the wave should be
damped~\cite{Kadomtzev_1978, Shafranov_1965}.

Below we consider the case when the inverse inequality is fulfilled, that
means the predominance of the electrons' oscillation frequency in the well
over the damping rate $-\gamma$. Then the captured electrons are reflected
many times from the cavity walls (getting and losing the energy) and on
average the wave keeps its energy -- hence at $\sqrt{eE_{0} /T_{e} \cdot k} \cdot k\gg\gamma$ the damping of the wave does not take place~\cite{Kadomtzev_1978, Shafranov_1965}.

To simplify further calculations, it is convenient to introduce the
following auxiliary value,

\begin{equation}
\label{eq36}
I\left( {z,t} \right)=\int\limits_{t_{0} \to -\infty }^t {d{t}'\cdot E\left(
{z,{t}'} \right)},
\end{equation}

for which from Eq.~\eqref{eq31} we can obtain the following equation, neglecting the
term corresponding to the wave damping:

\begin{equation}
\label{eq37}
\frac{\partial^{4}I}{\partial t^{4}}+\frac{\partial^{2}I}{\partial
t^{2}}+3\frac{\partial^{2}I}{\partial z^{2}}+
\frac{\partial }{\partial z}\frac{\partial }{\partial t}\left\{ {I^{2}\left(
{z,t} \right)-\frac{1}{2}\left( {\frac{\partial I\left( {z,t}
\right)}{\partial t}} \right)^{2}} \right\}=0.
\end{equation}

By means of Eq.~\eqref{eq37} we
can construct the expression for the value $\frac{1}{2}\left(
{\frac{\partial I\left( {z,t} \right)}{\partial t}} \right)^{2}$. The
straightforward calculations give

\begin{equation}
\label{eq38}
\eqalign{
\frac{1}{2}\left( {\frac{\partial I\left( {z,t} \right)}{\partial t}}
\right)^{2}
\cong -\frac{1}{2}I^{2}\left( {z,t} \right)-3\int\limits_{-\infty
}^t {d{t}'\frac{\partial I\left( {z,{t}'} \right)}{\partial {t}'}}
\frac{\partial^{2}}{\partial z^{2}}\cdot\cr
\cdot \int\limits_{-\infty }^{{t}'}
{d{t}''\int\limits_{-\infty }^{{t}''} {d{t}'''\cdot I\left( {z,{t}'''}
\right)} } -\int\limits_{-\infty }^t {d{t}'\frac{\partial I\left( {z,{t}'}
\right)}{\partial {t}'}\frac{\partial }{\partial z}}
\cdot\int\limits_{-\infty
}^{{t}'} {d{t}''\cdot I^{2}\left( {z,{t}''} \right)+}  \mathellipsis .}
\end{equation}

Further we a) hold on to the approximation usually used in the theory of
Landau damping -- the assumption of the smooth dependence of the
electric field on the spatial coordinate. Therefore in Eq.~\eqref{eq37} the terms only up to the second derivative with the spatial $z-$coordinate are kept, and b) use
the relation $\left( {\partial I/\partial t} \right)=E\left( {z,t}
\right)=-\partial \phi /\partial z$, which follows from Eq.~\eqref{eq36} and
restrict with the quadratic term $I^{2}$ in the last term in the
rhs of Eq.~\eqref{eq38}. Assuming the dependence of all unknown values in the
argument $\xi =z-V\cdot t$, we find that $I\left( \xi \right)=\left( {1/V}
\right)\cdot \phi \left( \xi \right)$ , where $V$ is the dimensionless
velocity of the waves, which is assumed to be large, $V\gg 1$. Under
restrictions a) and b) the last two terms in the rhs of Eq.~\eqref{eq38} give
negligibly small contributions. Substituting the remaining first term in the
rhs into Eq.~\eqref{eq37} after two-times integration we obtain

\begin{equation}
\label{eq39}
V^{4}\frac{\partial^{2}\phi }{\partial \xi^{2}}+\left( {V^{2}+3}
\right)\cdot \phi -\frac{3}{2}\phi^{2}=\frac{1}{2}C_{1} .
\end{equation}

Here the constant $C_{1} $does not depend on time. In the rhs of Eq.~\eqref{eq39} we
have neglected the term proportional to $\xi $, leading to the nonphysical
result at $\xi \to \infty $. Multiplying Eq.~\eqref{eq39} by $\left( {\partial \phi
/\partial \xi } \right)$ after integration we find:

\begin{equation}
\label{eq40}
V^{4}\left( {\frac{\partial \phi }{\partial \xi }} \right)^{2}=\phi
^{3}-\left( {V^{2}+3} \right)\cdot \phi^{2}+C_{1} \cdot \phi +C_{2} ,
\end{equation}

where the constants of integrations $C_{1} $ and $C_{2} $ can be expressed
in terms of the minimum $- \quad \phi_{n} $ and the maximum $- \quad \phi_{m} $
values of the potential, $\left. {\left( {\partial \phi /\partial \xi }
\right)} \right|_{\phi_{m} ,\phi_{n} } =0$. As a result we obtain

\begin{equation}
\label{eq41}
V^{4}\left( {\frac{\partial \phi }{\partial \xi }} \right)^{2}=(\phi_{m}
-\phi )\cdot \left( {\phi_{s} -\phi } \right)\cdot \left( {\phi -\phi_{n}
} \right).
\end{equation}

The constant $\phi_{s} $ (together with $\phi_{m} $ and $\phi_{n} )$
defines the wave velocity $V$:

\begin{equation}
\label{eq42}
V^{2}=\phi_{m} +\phi_{s} +\phi_{n} -3\quad .
\end{equation}

At $\phi_{m} >\phi_{s} >\phi_{n} $ the solution of Eq.~\eqref{eq41} represents a
periodic wave, described by the expression (see the Fig.~\ref{fig:1})

\begin{equation}
\label{eq43}
\phi =\phi_{m} -\left( {\phi_{m} -\phi_{n} } \right)\cdot dn^{2}\left\{
{x,s} \right\},
\end{equation}

where $dn\left\{ {x,s} \right\}$ is a Jacobi elliptic function with the
module $-s$~\cite{Jahnke_1945},

\begin{equation}
\label{eq44}
dn\left\{ {x,s} \right\}=\frac{\pi }{2\cdot K\left( s \right)}+\frac{2\pi
}{K\left( s \right)}\sum\nolimits_{n=1}^\infty {\frac{q^{n}}{1+q^{2n}}\cos
\left\{ {\frac{\pi \cdot n\cdot x}{K\left( s \right)}} \right\}} ,
\end{equation}

\begin{equation}
\label{eq45}
\eqalign{
q=\exp \left[ {-\pi \frac{K\left( {{s}'} \right)}{K\left( s \right)}}
\right],\quad {s}'=\sqrt {1-s^{2}}, \cr
x=\frac{\sqrt {\phi_{m} -\phi_{n} } }{2\cdot V^{2}}\cdot \xi,
\quad s^{2}=\frac{\phi_{s} -\phi_{n} }{\phi_{m} -\phi_{n}}.}
\end{equation}

The function $dn^{2}\left\{ {x,s} \right\}$ is a periodic function with the
period $2K\left( s \right)$, where $K\left( s \right)$ is the complete
elliptic integral of the first kind,

\begin{equation}
\label{eq46}
K\left( s \right)=\int\limits_0^{\pi /2} {\frac{d\alpha }{\sqrt
{1-s^{2}\cdot \sin \alpha } }} ,
\end{equation}

therefore the wave length of the periodic solution \eqref{eq43} can be defined
according to the relation

\begin{equation}
\label{eq47}
\lambda =\frac{4}{\sqrt {\phi_{m} -\phi_{n} } }\cdot V^{2}K\left( s
\right)
\end{equation}

\begin{figure}[htpb]
\centering
\includegraphics{NonlinearLandau_fig1}
\caption{Sketch of the dependence of $dn^{2}\left\{ {x,s} \right\}$ on $x$ at $s=const; \Delta \phi =\left( {1-a^{2}} \right)\cdot \left( {\phi_{m} -\phi_{n} } \right)$, $a$ is the minimum value of the curve $dn\left\{ {x,s} \right\}$}
\label{fig:1}
\end{figure}

\section{Instability of waves}
We can now analyze the stability of the waves \eqref{eq43}. Starting from Eq.~\eqref{eq39}
we introduce the potential perturbation $\delta \phi $ according to the
relation $\phi \to \phi +\delta \phi $ . For the unperturbed part of the
potential$-\phi $ remains the expression \eqref{eq43} and for $\delta \phi $ we
obtain

\begin{equation}
\label{eq48}
\frac{\partial^{2}\delta \phi }{\partial \xi^{2}}+\frac{1}{V^{2}}\left[
{1+\frac{3-\phi }{V^{2}}+\frac{3}{V^{2}}\left\{ {1-dn^{2}\left( {z,s^{2}}
\right)} \right\}} \right]\cdot \delta \phi =0.
\end{equation}

An explicit solution of this equation with the periodic coefficient can be
found applying a Hill's method~\cite{Whittaker_1973}. This method is rather cumbersome with
the long and difficult calculations. Here we simplify the equation recalling
the condition used above, $V \gg 1$, and assuming that
the parameter $s$ is small, $s \ll 1$. Then Eq.~\eqref{eq48} can be reduced to
Mathieu's equation~\cite{Whittaker_1973, landau05:_mech}, which describes the phenomenon known as a parametric resonance,

\begin{equation}
\label{eq49}
\frac{\partial^{2}\delta \phi }{\partial \xi^{2}}+\frac{1}{V^{2}}\left\{
{1+h\cdot \cos k\xi } \right\}\cdot \delta \phi =0 ,
\end{equation}

where $h\cong \frac{12}{V^{2}}\frac{\phi_{m} -\phi_{n} }{Ch\pi }$ and
$k\cong \frac{\sqrt {\phi_{m} -\phi_{n} } }{V^{2}}$. Applying a standard
method at the fulfilment of the resonance condition $k=2/V$ Eq.~\eqref{eq49} gives
the following expression for the rate of the instability $\bar{{\gamma
}}=\frac{1}{4}\frac{h}{V}$~\cite{landau05:_mech}.

\section{Summary}

Finding the exact analytic solution for the collisionless Vlasov's equation,
the nonlinear stage of the Langmuir waves is analyzed in the first
non-vanishing (quadratic) nonlinear approximation. The Langmuir waves with
the finite amplitude, and with the oscillation frequency (of electrons in
the wave-well) larger than the damping rate (found in the linear
approximation), do not damp and tend to keep the periodic structure. On
the fulfilment of the definite resonance conditions, the waves are unstable.
The finding of the corresponding rate is quite similar to the procedure applied in
the investigation of the parametric instability.

\section{Acknowledgment}

This work has been partially supported by the grants P2-0073 and P2-0256 of the Slovenian Research Agency. IV would like to acknowledge the financial support from the EUROfusion Consortium under grant agreement No 633053  (WP-EDU). The views and opinions expressed herein do not necessarily reflect those of the European Commission.

\newcommand{}{}
\bibliographystyle{unsrtnat}
\providecommand{\noopsort}[1]{}\providecommand{\singleletter}[1]{#1}%