text:
\begin{document}

\maketitle

\begin{abstract}
\noindent A translation surface $\S$ is said to have the finite
blocking property if for every pair $(O,A)$ of points in $\S$
there exists a finite number of ``blocking'' points $B_1, \dots ,
B_n$ such that every geodesic from $O$ to $A$ meets one of the
$B_i$'s. $\S$ is said to be purely periodic if the directional
flow is periodic in each direction whose directional flow contains
a periodic trajectory (this implies that $\S$ admits a cylinder
decomposition in such directions). We will prove that finite
blocking property implies pure periodicity. We will also classify
the surfaces that have the finite blocking property in genus $2$:
such surfaces are exactly the torus branched coverings.
Moreover, we prove that in every stratum, such surfaces form a
set of null measure. In the Appendix, we prove that completely periodic translation surfaces form a set of null measure in every stratum.\\

{\em \noindent Keywords: blocking property, polygonal billiards,
translation surfaces, Veech surfaces, torus branched covering,
complete periodicity, illumination, quadratic differentials.}\\

{\em \noindent AMS classification:
37D50,
37C27,
57M12,
51E21,
37E35,
32G15.
}

\end{abstract}

\clearpage

\section*{Introduction}

A {\em translation surface} is a triple $(\S, \Sigma, \omega)$
such that $\S$ is a topological compact connected surface,
$\Sigma$ is a finite subset of $\S$ (whose elements are called
{\em singularities}) and $\omega = (U_i,\phi_i)_{i\in I}$ is an
atlas of $\S \setminus \Sigma$ (consistent with the topological
structure on $\S$) such that the transition maps (i.e. the $\phi_j
\circ \phi_i^{-1} : \phi_i(U_i\cap U_j) \rightarrow \phi_j(U_i\cap
U_j)$ for $(i,j)\in I^2$) are translations. This atlas gives to
$\S \setminus \Sigma$ a Riemannian structure; therefore, we have
notions of length, angle, measure, geodesic... Moreover, we assume
that $\S$ is the completion of $\S \setminus \Sigma$ for this
metric.

A translation surface can also be seen as a holomorphic
differential $h$ on a Riemann surface (the singularities
correspond to the zeroes of the differential, and in an admissible
atlas $\omega$, $h$ is of the form $h=dz$). This point of view is
useful to give coordinates (and therefore a topology and a
measure) to the moduli space of all translation surfaces,
while the first one allows to make pictures in the plane.

Translation surfaces provide one of the main tool for the study of rational polygonal billiards.\\ \\

Since the unit tangent bundle of $\S$ enjoys a canonical global
decomposition $U\S=\S \times \mathbb{S}^1$
, the study of the geodesic flow on $\S$ can be done
through two points of view depending on whether the variable is
the first or the second projection:

\begin{enumerate}

\item We can fix one particular direction $\theta\in\mathbb{S}^1$:
this corresponds to the study of the {\em directional flow}
$\phi_{\theta} : \S\times \R \rightarrow \S $. We are usually
interested in the ergodic properties of this flow depending on
whether $\theta$ is a saddle connection direction or not, or on
the way $\theta$ can be approximated by such saddle connections
(see \cite{KMS}, \cite{Ve}, \cite{MT}, \cite{Vo}, \cite{Ma}, \cite{MS}, \cite{Ch}).
This class of problems is called {\em dynamical problems}.
This point of view is the most studied by mathematicians.

\item We can also fix one point in $x\in\S$: this corresponds to
the study of the {\em exponential flow} $exp_x : \mathbb{S}^1
\times \R \rightarrow \S $. We are usually interested in looking
at which points we can reach when we let the second variable vary,
or at the statistics of the return to (or near) $x$ (see \cite{Kl},

\cite{To}, \cite{Mo-HS}, \cite{Mo-V}, \cite{BGZ}, \cite{CG}).
This class of problems is called {\em illumination problems}.
This point of view is connected with computer science
(computational geometry, ray tracing, art gallery).

\end{enumerate}

Note that those two flows are defined on a dense $G_\delta$-set of
full measure but not everywhere because the singularities don't
allow some (few) geodesics to be defined on $\R$.\\

We are interested here by the connection between an illumination
problem called the finite blocking property, and a dynamical one
named the pure periodicity.\\

A translation surface $\S$ is said to have the {\em finite
blocking property} if for every pair  $(O,A)$ of points in $\S$,
there exists a finite number of points $B_1, \dots , B_n$
(different from $O$ and $A$) such that every geodesic from $O$ to
$A$ meets one of the $B_i$'s.\\

A translation surface $\S$ is said to be {\em purely periodic} if
for any $\theta\in\mathbb{S}^1$, the existence of a (non-singular) periodic orbit
in the direction $\theta$ implies that the directional flow
$\phi_{\theta}$ is periodic (i.e. there exists $T>0$ such that
$\phi_{\theta}^T=Id_{\S}$ a.e.).
In particular, the widths of the cylinders in such a direction $\theta$ are commensurable.
This notion is a stronger version of the complete periodicity defined in \cite{Cal}.\\

The paper is organized as follows:

In the first section, we recall some background and known results
about the finite blocking property. The second one is devoted to
the proof of the following result:\\

{\bf Theorem \ref{th-periodicity}.} {\em Let $\S$ be a translation
surface with the finite blocking property. Then $\S$ is purely periodic.}\\

The third section is devoted to some applications of this theorem.
First, we give an improvement of \cite{Mo-V}, Theorem 3:\\

{\bf Theorem \ref{th-measure}.} {\em In genus $g\geq 2$, the set
of translation surfaces that fail the finite blocking property is
of full measure in every stratum.}\\

Then, we give the complete classification in genus $2$:\\

{\bf Theorem \ref{th-genus2}.} {\em Let $\S$ be a translation
surface of genus $2$. Then $\S$ has the finite blocking property
if and only if $\S$ is a torus branched covering.}\\

In the Appendix, we prove the following known result which is needed in the proof of Theorem \ref{th-measure}:\\

{\bf Theorem \ref{th-mesure-zero-completement-periodiques}.}
{\em In genus $g\geq 2$, the set of completely periodic translation surfaces has measure zero in every stratum.}\\

\clearpage

\section{Background}\label{section-background}
In this section, we present some results about the finite blocking property
that were proved in \cite{Mo-V}.

\subsection{Stability}
The finite blocking property is stable under some classical
operations on translation surfaces:

\begin{description}

\item[Singularities]

We have used the convention that a geodesic stops when it reaches
a singularity. This convention is not restrictive for the study of
the finite blocking property since we can add those singularities
in any blocking configuration. Moreover the finite blocking
property is not perturbed when we add or remove removable
singularities (a singularity $\sigma \in \Sigma$ is said to be
{\em removable} if there exists an atlas $\omega' \supset \omega$
such that $(\S, \Sigma \setminus \{ \sigma \}, \omega')$ is a
translation surface).

\item[Branched covering]

A {\em branched translation covering} between two translation surfaces is a
map $\pi : (\S,\Sigma) \rightarrow (\S',\Sigma')$ that is
a topological branched covering that locally preserves the
translation structure. A {\em torus branched covering} is a
translation surface $\S$ such that there exists a branched translation covering
from $\S$ to a flat torus.

\begin{prop}\label{prop-covering}
Let $\pi : (\S,\Sigma) \rightarrow (\S',\Sigma')$ be a translation covering branched on a finite set
$\mathcal{R}' \subset \S'$. Then $\S$ has the finite blocking property if and
only if $\S'$ has.
\end{prop}

\item[Rational billiards]

Let $\P$ denote a polygon in $\mathbb{R}^2$. Let $\Gamma \subset
O(2,\R)$ be the group generated by the linear parts of the
reflections in the sides of $\P$. When $\Gamma$ is finite, we say
that $\P$ is a {\em rational polygonal billiard table}. When $\P$ is
simply connected, $\P$ is rational if and only if all the angles
between edges are rational multiples of $\pi$. A classical
construction due to Zemljakov and Katok (see \cite{ZK}, \cite{MT})
allows us to associate to each rational billiard table $\P$ a
translation surface $ZK(\P)$
(as an example, the translation surface associated to the square
billiard table is the torus $\mathbb{T}^2 = \R^2 / \Z^2 $).

\begin{prop}\label{prop-ZK}
Let $\P$ be a rational polygonal billiard table. Then $ZK(\P)$ has the
finite blocking property if and only if $\P$ has.
\end{prop}

\item[Action of \mathversion{bold}$GL(2,\R)$ \mathversion{normal}]

If $A\in GL(2,\R)$, we can define the translation surface $$A .
(\S, \Sigma,(U_i,\phi_i)_{i\in I}) \defi (\S, \Sigma, (U_i,A \circ
\phi_i)_{i\in I});$$ hence we have an action of $GL(2,\R)$ on the
class of translation surfaces. We classically consider only
elements of $SL(2,\R)$ (see \cite{MT}), but area is not assumed to be preserved here.

\begin{prop}\label{prop-GL2R}
Let $\S$ be a translation surface and $A$ be in $GL(2,\R)$. Then
$\S$ has the finite blocking property if and only if $A.\S$ has.
\end{prop}

\end{description}

\subsection{Sufficient conditions}

\begin{prop}[Fomin \cite{Fomin}]\label{prop-Fomin}
The square billiard table has the finite blocking property.
\end{prop}

Hence, combining propositions \ref{prop-covering}, \ref{prop-ZK}, \ref{prop-GL2R}
and \ref{prop-Fomin}, we have:

\begin{prop} \label{prop-torus-branched-covering}
Any torus branched covering has the finite blocking property.
\end{prop}

\subsection{Necessary conditions}
We first give some definitions and conventions concerning isometries and cylinders.
In this paper, ``isometry'' means ``orientation-preserving local isometry''.

Since an isometry $i$ from an open subset $U$ of $\R^2$ to a translation surface $\S$ can be uniquely extended by continuity to
$\overline{U}$ (it is uniformly continuous), we will also denote by $i$ this extension (that is also an isometry).

In this process, we may lose the injectivity of $i$.\\

A {\em subcylinder} $\C$ is an isometric copy of $\R / w \Z \times
(0,h)$ in a translation surface $\S$ ($w>0$, $h>0$). $w$ and $h$
are unique and called the {\em width} and the {\em height} of $\C$.
The {\em direction} of $\C$ is the direction of the image of $\R / w \Z \times \{ h /2 \}$
(which is considered here as an oriented closed geodesic); it is defined modulo $2\pi$.
The images in $\S$ of $\R / w \Z \times \{0\}$ and  $\R / w \Z \times \{h\}$

are called the {\em sides} of $\C$.
The side $\R / w \Z \times \{h\}$ is called the {\em upper side}.

A {\em cylinder} is a maximal subcylinder (for the inclusion).
Any (regular) periodic trajectory for $\phi_\theta$ can be
thickened to obtain a cylinder in the direction $\theta$.
Whenever $\S$ is not a torus, each side of a cylinder contains at least one
singularity and is a finite union of saddle connections
(a {\em saddle connection} is a geodesic joining two singularties).
We say that $\S$ admits a {\em cylinder decomposition} in the direction $\theta$ if
the (necessarily finite) union of cylinders whose direction is $\theta$ is dense in $\S$
(then, the complement of this union is the union of the sides of those cylinders).\\

Let $\C_1$ and $\C_2$ be two cylinders with width $w_1$ (\resp
$w_2$) and height $h_1$ (resp. $h_2$) in a translation surface
$\S$. We suppose that $\C_1$ and $\C_2$ are two different
cylinders in the same direction whose closures have a nontrivial
(i.e. not reduced to a finite set) intersection. The situation is described in Figure \ref{2cylinder.eps}.

\begin{figure}[!h]
\begin{center}
\psfrag{w1}{$w_1$} \psfrag{h1}{$h_1$} \psfrag{w2}{$w_2$}
\psfrag{h2}{$h_2$} \psfrag{C1}{$\C_1$} \psfrag{C2}{$\C_2$}
\psfrag{O}{$(0,0)$} \psfrag{l}{$l$}
\includegraphics[scale = 0.56]{2cylinder.eps}
\caption{\label{2cylinder.eps} The opposite vertical sides are
identified by translation while the dotted horizontal ones are
glued with the rest of the surface; $l>0$.}
\end{center}
\end{figure}

\begin{lem}\ \label{lem2cyln}
Let $\S$ be a translation surface that contains two different
cylinders $\C_1$ and $\C_2$ in the same direction whose closures
have a nontrivial intersection. If their widths are
uncommensurable, then $\S$ fails the finite blocking property.
\end{lem}

This lemma implies all the remaining results of section
\ref{section-background} and will be useful in the proof of
Theorem \ref{th-periodicity}.

\begin{prop}
The only regular polygons (considered as billiard tables) with the
finite blocking property are the equilateral triangle, the square
and the regular hexagon.
\end{prop}

\begin{prop}\label{prop-Veech}
A Veech surface has the finite blocking property if and only if it
is a torus covering, branched over only one point.
\end{prop}

In other terms, arithmetic surfaces are exactly the Veech surfaces
with the finite blocking property.
Note that proposition \ref{prop-Veech} was proven
independently
by Eugene Gutkin, and is the main result of
\cite{Gutkin}.

\subsection{Global point of view: moduli spaces of holomorphic forms}\label{subsection-moduli-space}

A singularity $\sigma \in \Sigma$ has a cone angle of the form
$2(k+1)\pi$, with $k\geq 0$; we say that $\sigma$ is of {\em
multiplicity} $k$. Removable singularities are the singularities of multiplicity $0$.
In terms of holomorphic differentials, it is equivalent to say that
there is a chart around $\sigma$ such that $h=z^{k}dz$ (i.e. $\sigma$ is a zero of order $k$ of $h$).

If $1\leq k_1\leq k_2 \leq \dots \leq k_n$ is
a sequence of integers whose sum is even, we denote by $\H(k_1, k_2,
\dots, k_n)$ the {\em stratum} of translation surfaces with
exactly $n$ singularities whose multiplicities are
$k_1, k_2, \dots, k_n$ (we consider only surfaces without removable singularities). A translation surface in $\H(k_1, k_2,
\dots, k_n)$ has genus $g=1+(k_1+k_2+ \dots +k_n)/2$.

Each stratum carries a natural topology and a $SL(2,\R)$-invariant measure that are
defined in \cite{Masur-ei} and \cite{Veech-ei}. Let $\S$ be an element of some $\H(k_1, k_2, \dots, k_n)$ and
let $\B=\{\gamma_1, \dots , \gamma_{2g+n-1}\}$ be a basis of the {\em relative} homology of $\S$:
it is just the concatenation of a basis of the first topological homology group $H_1(\S,\Z)$
with a set of $n-1$ curves from a singularity to the other ones.
If $\S'$ is another translation surface (built on the same topological surface),
let us denote by $hol_{\S'}$ the associated (translational) holonomy.
The map
$$\Phi \defi \ba \H(k_1, k_2,\dots, k_n)  \dans \mathbb{C}^{2g+n-1} \simeq (\R^2)^{2g+n-1} \\
\S'  \donne (hol_{\S'}(\gamma_1) , \dots , hol_{\S'}(\gamma_{2g+n-1}))  \ea$$
is named the {\em period map} and is a local homeomorphism in a neighbourhood of $\S$
and in this system of coordinates, the former measure is absolutely continuous relatively to Lebesgue.\\

\begin{prop}\label{prop-dense}
In genus $g\geq 2$, the set of translation surfaces that fail the
finite blocking property is dense in every stratum.
\end{prop}

In fact, we proved in \cite{Mo-V} that the hypotheses of Lemma \ref{lem2cyln} are satisfied on a dense subset of each stratum.
Unfortunately, those hypotheses are valid only on a set of zero measure.
Proposition \ref{prop-dense} will be considerably improved by Theorem \ref{th-measure}.

\clearpage

\section{Main result}

A translation surface is said to be {\em completely periodic} if the
following property holds (see \cite{Cal}): if there exists a periodic
orbit in the direction $\theta$, then the surface admits a cylinder
decomposition in the direction $\theta$.

\begin{prop}\label{prop-complete-periodicity}
Let $\S$ be a translation surface with the finite blocking
property. Then $\S$ is completely periodic.
\end{prop}

\begin{demo} Let $\theta$ be a direction that contains a periodic
orbit. Up to a rotation, we assume
that this direction is the horizontal one (defined by $\theta=0 \mod 2\pi$). We can thicken all of
such horizontal periodic orbits to obtain a finite number of cylinders. Assume by
contradiction, that $\S$ is not fully decomposed into cylinders in
the horizontal direction: there exists a cylinder $\C$ and singularities $\sigma$ and $\sigma'$ (possibly equal) that lie on the upper side of $\C$
such that the horizontal geodesic $\gamma$ (whose length is denoted by $l(\gamma)$ and image by $|\gamma|$)
from $\sigma$ to $\sigma'$ exists (the singularities are consecutive for the cyclic order on the upper side of $\C$)
and is not a part of a side of any other cylinder.
To be more precise, there exists $s>0$ (by compactness of $|\gamma|$ and by finiteness of the number of horizontal cylinders)
such that there is an injective isometry $i: (0,l(\gamma)) \times  (0,s) \rightarrow \S$ such that $i([0,l(\gamma)]\times\{0\})=|\gamma|$

and $i( (0,l(\gamma))\times (0,s))$ doesn't meet any horizontal cylinder.
Up to a dilatation (that belongs to $GL^{+}(2,\R)$), we assume that the width of $\C$ is $1$ and that its height is $2$.
The situation is described in Figure \ref{cyln-minim.eps}.

\begin{figure}[!h]
\begin{center}
\psfrag{g}{$|\gamma|$}
\psfrag{1}{$1$}
\psfrag{2}{$2$}
\psfrag{lg}{$l(\gamma)$}
\psfrag{C}{$\C$}
\psfrag{s}{$s$}
\psfrag{S}{$\sigma$}
\psfrag{S'}{$\sigma'$}
\includegraphics[scale = 0.58]{cyln-minim.eps}
\caption{\label{cyln-minim.eps} The opposite vertical black sides are
identified by translation; the dotted ones are glued with the rest of the surface;
the dark zone doesn't meet any horizontal cylinder.}
\end{center}
\end{figure}

The heart of the proof is organised into three steps. A summary
concludes each of them.

\begin{description}

\item[Step 1] {\em Construction of certain thin rectangles in $\S$ that contain a common point $A$.}

Let $\varepsilon > 0$ be smaller than $s$.

Let $T_{\varepsilon}$ be the supremum of all $t>0$ such that there exists an injective isometry
$j_t : (0,t) \times (0,\varepsilon) \rightarrow \S$ that coincides with $\gamma$ on $(0,\min(l(\gamma),t)) \times\{0\}$.
$T_{\varepsilon}$ is greater than or equal to $l(\gamma)$ and therefore strictly positive.
$\S$ has finite area hence $T_{\varepsilon}$ is finite.
Moreover $T_{\varepsilon}$ is a maximum (indeed, if $0 < t \leq t' < T_{\varepsilon}$, then $j_{t'}$ extends $j_t$):
there exists an injective isometry $j_{T_\varepsilon} : (0,T_{\varepsilon}) \times (0,\varepsilon) \rightarrow \S$
that coincide with $\gamma$ on $(0,l(\gamma)) \times\{0\}$.

There are only three possibilities depending on what happens on $j_{T_\varepsilon}(\{T_{\varepsilon}\}\times (0,\varepsilon))$
(if $j_{T_\varepsilon}(\{T_{\varepsilon}\}\times(0,\varepsilon))$ does not meet a singularity,
then $j_{T_\varepsilon}$ can be extended to an isometry that fails to be injective,
hence $I_{\varepsilon}=j_{T_\varepsilon}(\{T_{\varepsilon}\}\times (0,\varepsilon))$
and $J_{\varepsilon}=j_{T_\varepsilon}(\{0\}\times (0,\varepsilon))$ intersect each other):

\begin{enumerate}

\item $I_{\varepsilon}$ meets a singularity denoted by $\sigma_{\varepsilon}$.

\item $I_{\varepsilon} = J_{\varepsilon} $:
we just have constructed a cylinder, which is forbidden by the assumption, so this case does not occur.

\item The vertical intervals $I_{\varepsilon}$ and $J_{\varepsilon}$
intersect each other but we are neither in the first nor in the second case. Since $\sigma=j_{T_\varepsilon}(0,0)$ is a singularity,
$I_\varepsilon \cap J_{\varepsilon}$ is a subinterval of $J_{\varepsilon}$ that hits the extremity of $J_{\varepsilon}$ that is not $\sigma$
(i.e. $j_{T_\varepsilon}(0,\varepsilon)$) (see Figure \ref{enroulement.eps}).

\begin{figure}[!h]
\begin{center}
\psfrag{g}{$|\gamma|$}
\psfrag{1}{$1$}
\psfrag{2}{$2$}
\psfrag{e}{$\varepsilon$}
\psfrag{C}{$\C$}
\psfrag{s}{$s$}
\psfrag{S}{$\sigma$}
\psfrag{S'}{$\sigma'$}
\psfrag{O}{$O$}
\psfrag{A}{$A$}
\psfrag{h}{$h$}
\psfrag{l}{$l$}
\psfrag{e'}{$\varepsilon'$}
\psfrag{Ie}{$I_{\varepsilon}$}
\psfrag{Je}{$J_{\varepsilon}$}
\includegraphics[scale = 0.58]{enroulement.eps}
\caption{\label{enroulement.eps} The dark zone is the image of $j_{T_\varepsilon}$.}
\end{center}
\end{figure}

So, there exists $0 < \varepsilon' \leq \varepsilon$ such that $J_{\varepsilon} \setminus I_{\varepsilon} =  j_{T_\varepsilon}(\{0\}\times (0,\varepsilon'])$.

Now, the same three possibilities are offered to $j_{T_{\varepsilon'}}$, but the third case can't happen. Indeed, otherwise there is $0<\varepsilon'' < \varepsilon'$ such that $j_{T_{\varepsilon'}}(0,\varepsilon')=j_{T_{\varepsilon'}}(T_{\varepsilon},0)=j_{T_{\varepsilon'}}(T_{\varepsilon'},\varepsilon'-\varepsilon'')$ is not a singularity but
coming out of it are at least 2 horizontal lines
(equivalently, this point has a conical angle greater than $2\pi$).

\end{enumerate}

So, up to reducing the size of $\varepsilon$ (by replacing $\varepsilon$ by $\varepsilon'$ whenever we are in the third case),
we can consider that we are in the first case, and,
since $\Sigma$ is finite, there exists a sequence $(\varepsilon_n)_{n\in\N}$ of strictly positive numbers that converges to $0$
and a singularity $\sigma_{\infty} \in \Sigma$ such that for all $n$ in $\N$, $\sigma_{\varepsilon_n}=\sigma_{\infty}$.
If $n\in \N$, there exists $0<y_n<\varepsilon_{n}$ such that $\sigma_{\infty}=j_{T_{\varepsilon_{n}}}(T_{\varepsilon_n}, y_n)$.
Let $A=j_{T_{\varepsilon_{n}}}(T_{\varepsilon_n} - l(\gamma)/2 , y_n)$
($A$ does not depend on $n\in\N$ and is in the image of every $j_{T_{\varepsilon_n}}$).

{\em
To sum up this step and fix notations, we proved that there exists
$A\in \S$ such that for all $\varepsilon > 0$ there exists
$l_{\varepsilon} \geq l(\gamma)$ and
$0 < h_{\varepsilon} \leq \varepsilon$
such that there exists an injective isometry $i_{\varepsilon}$ from
$(0,l_{\varepsilon}) \times (0,h_{\varepsilon})$ to $\S$
that coincides with $\gamma$ on $(0,l(\gamma)) \times\{0\}$
and whose image contains $A$.}

The situation is summarized in Figure \ref{situation-step-1.eps}. Let $O$ be the point drawn on this Figure.
We denote $A=i_{\varepsilon}(x_{\varepsilon}, y_{\varepsilon})$.

\begin{figure}[!h]
\begin{center}
\psfrag{g}{$|\gamma|$}
\psfrag{1}{$1$}

\psfrag{e}{$\varepsilon$}
\psfrag{C}{$\C$}
\psfrag{s}{$s$}
\psfrag{S}{$\sigma$}
\psfrag{S'}{$\sigma'$}
\psfrag{O}{$O$}
\psfrag{A}{$A$}
\psfrag{h}{$h_{\varepsilon}$}
\psfrag{l}{$l_{\varepsilon}$}
\psfrag{y}{$y_{\varepsilon}$}
\includegraphics[scale = 0.6]{situation-step-1.eps}
\caption{\label{situation-step-1.eps}

Summary of Step 1. }
\end{center}
\end{figure}

\item[Step 2] {\em Construction of ``teepees'' of more and more geodesics between $O$ and $A$.}

Let $\varepsilon>0$.
The unfolding of $\C$ allows us to define an extension of $i_{\varepsilon}$ as in Figure \ref{thales1.eps}:
we get an isometry $i_{\varepsilon} : (0,l_{\varepsilon}) \times [0,h_{\varepsilon}) \bigcup \R\times (-2,0) \rightarrow \S$.
In $\R^2$, we denote without ambiguity $\sigma=(0,0)$, $\sigma'=(l(\gamma),0)$ and $A=(x_{\varepsilon}, y_{\varepsilon})$.

\begin{figure}[!h]
\begin{center}
\psfrag{g}{$|\gamma|$}
\psfrag{1}{$1$}

\psfrag{e}{$\varepsilon$}
\psfrag{C}{$\C$}
\psfrag{s}{$s$}
\psfrag{S}{$\sigma$}
\psfrag{S'}{$\sigma'$}
\psfrag{s}{$S_{\varepsilon}$}
\psfrag{s'}{$S'_{\varepsilon}$}
\psfrag{O}{$O$}
\psfrag{A}{$A$}
\psfrag{h}{$h_{\varepsilon}$}
\psfrag{lg}{$l(\gamma)$}
\psfrag{d}{$d(S_{\varepsilon},S_{\varepsilon}')$}
\psfrag{y}{$y_{\varepsilon}$}
\includegraphics[scale = 0.44]{thales1.eps}
\caption{\label{thales1.eps} First (topological) application of the theorem of Thales.}
\end{center}
\end{figure}

Let $S_{\varepsilon}$ (resp. $S_{\varepsilon}'$) be the point of the line that passes through $A$
and $\sigma$ (resp. $\sigma'$) whose y-coordinate is $-1$.
Since $(\sigma, \sigma')$ and $(S_{\varepsilon}, S_{\varepsilon}')$ are parallel,
the theorem of Thales (about similar triangles) asserts that $$\frac{d(S_{\varepsilon},S_{\varepsilon}')}{ (1+y_{\varepsilon})} = \frac{l(\gamma)}{ y_{\varepsilon}} , $$
where $d$ stands for the euclidean distance of $\R^2$.
Let $\mathcal{O}_{\varepsilon} = i_{\varepsilon}^{-1}(\{O\})\cap (S_{\varepsilon}, S_{\varepsilon}')$.
We have $$\card \, \mathcal{O}_{\varepsilon} \geq d(S_{\varepsilon},S_{\varepsilon}') -1 = l(\gamma) (1+y_{\varepsilon}) / y_{\varepsilon} -1  \xrightarrow[ \varepsilon \rightarrow 0]{} \infty .$$

Since the open triangle $(A,S_{\varepsilon},S_{\varepsilon}')$ is included in the domain of $i_{\varepsilon}$, we define
a set $\T_{\varepsilon}$ (called a {\em teepee}) of geodesics in $\S$ from $O$ to $A$ by taking the images by $i_{\varepsilon}$
of the segments joining an element of $\mathcal{O}_{\varepsilon}$ to $A$.
Of course, $$\card \, \T_{\varepsilon} = \card \, \mathcal{O}_{\varepsilon}  \xrightarrow[ \varepsilon \rightarrow 0]{} \infty  . $$

{\em
To sum up this step, we constructed, for all $\varepsilon > 0$, a set $\T_{\varepsilon}$
of more and more geodesics from $O$ to $A$ whose images are included in
$\C \cup | \gamma | \cup i_{\varepsilon} ((0,l_{\varepsilon}) \times (0,h_{\varepsilon}))$.
}

\item[Step 3] {\em The density of the geodesics of a teepee that are blocked by a fixed point vanishes.}

Let $B$ be a fixed point in $\S$ that is different from $O$ and $A$.
For $\varepsilon >0$, let $\T_{\varepsilon} (B)$ be the set of elements of $\T_{\varepsilon}$ that are blocked by $B$
(i.e. that pass through $B$).

\begin{itemize}
\item Suppose that $B$ is not in $\C$. Let $\varepsilon >0$.
If $B$ is not in the image of $i_{\varepsilon}$, then $\T_{\varepsilon} (B)$ is empty.
Otherwise, $B$ is in $i_{\varepsilon}((0,l_{\varepsilon}) \times [0,h_{\varepsilon}))$,
then the cardinal of $\T_{\varepsilon} (B)$ is at most $1$
(indeed, the elements of $\T_{\varepsilon}$ are disjoint on $i_{\varepsilon}((0,l_{\varepsilon}) \times [0,h_{\varepsilon})) \setminus \{A\}$).

Hence, $\card \, \T_{\varepsilon}(B) / \card \, \T_{\varepsilon}  \xrightarrow[ \varepsilon \rightarrow 0]{} 0$.

\item Suppose that $B$ is in $\C$.
Let $-h$ be the y-coordinate of any preimage of $B$ by any $i_{\varepsilon}$. We can assume that $h\in (0,1)$,
otherwise $\T_{\varepsilon} (B)$ is empty for any $\varepsilon >0$.
Let $\varepsilon > 0$.

Let $\beta_1$ and $\beta_2$ be two distinct elements of $\T_{\varepsilon}(B)$:
for $i\in\{1,2\}$, there exists $O_i$ in $\mathcal{O}_{\varepsilon}$ and $B_i$ in $i_{\varepsilon}^{-1}(B)$ such that
$B_i$ lies in the segment $[O_i,A]$ whose image by $i_{\varepsilon}$ is $\beta_{i}$ (see Figure \ref{thales2.eps}).

\begin{figure}[!h]
\begin{center}
\psfrag{g}{}

\psfrag{1}{$1$}

\psfrag{e}{$\varepsilon$}
\psfrag{C}{$\C$}
\psfrag{s}{$s$}
\psfrag{S}{$\sigma$}
\psfrag{S'}{$\sigma'$}
\psfrag{s}{$S_{\varepsilon}$}
\psfrag{s'}{$S'_{\varepsilon}$}
\psfrag{O1}{$O_1$}
\psfrag{O2}{$O_2$}
\psfrag{B1}{$B_1$}
\psfrag{B2}{$B_2$}
\psfrag{A}{$A$}
\psfrag{h}{$h$}
\psfrag{pe}{$p_{\varepsilon}$}
\psfrag{qe}{$q_{\varepsilon}$}
\psfrag{y}{$y_{\varepsilon}$}
\includegraphics[scale = 0.44]{thales2.eps}
\caption{\label{thales2.eps} Second (arithmetic) application of the theorem of Thales.}
\end{center}
\end{figure}

We choose $\beta_1$ and $\beta_2$ such that $d(O_1,O_2)$ is minimal, and set $p_{\varepsilon}=d(B_1,B_2)$ and $q_{\varepsilon}=d(O_1,O_2)$.
By minimality, $$ \card \, \T_{\varepsilon}(B) \leq \frac{\card \, \T_{\varepsilon}}{ q_{\varepsilon}} +1 . $$

Since $i_{\varepsilon}(B_1)=i_{\varepsilon}(B_2)=B$ and $i_{\varepsilon}(O_1)=i_{\varepsilon}(O_2)=O$,
$p_{\varepsilon}$ and $q_{\varepsilon}$ are
integers.

Since $(B_1, B_2)$ and $(O_1,O_2)$ are parallel,
the theorem of Thales asserts that $$\frac{p_{\varepsilon}}{q_{\varepsilon}} = \frac{(h+y_{\varepsilon})}{ (1+y_{\varepsilon})} . $$

This quotient
converges to $h$ when $\varepsilon$ converges to $0$ but never takes the value $h$,
hence the denominator $ q_{\varepsilon}  \xrightarrow[ \varepsilon \rightarrow 0]{} \infty$

Hence, $$\frac{\card \, \T_{\varepsilon}(B) }{ \card \, \T_{\varepsilon} } \leq \frac{1 }{ q_{\varepsilon} } + \frac{1 }{  \card \, \T_{\varepsilon}}  \xrightarrow[ \varepsilon \rightarrow 0]{} 0  . $$

\end{itemize}

{\em
To sum up this step, we proved that for any point $B$ in $\S$ that is different from $O$ and $A$,
the density of the geodesics of $\T_{\varepsilon}$ that are blocked by $B$ tends to $0$ when $\varepsilon$ tends to $0$.
}

\end{description}

Now, we assumed that $\S$ has the finite blocking property,
so there exists a finite set $\{B_1,\dots,B_n\}$ in $\S \setminus\{O,A\}$ that blocks any geodesic from $O$ to $A$.

Let $\varepsilon > 0$ be small enough such that $\Sigma_{i=1}^n \card \, \T_{\varepsilon}(B_i) / \card \, \T_{\varepsilon} < 1 $:
there exists a geodesic in $\T_{\varepsilon} \setminus \bigcup_{i=1}^{n} \T_{\varepsilon}(B_i)$ from $O$ to $A$ that doesn't meet any $B_i$,
leading to a contradiction.
So, $\S$ is completely periodic.

\end{demo}

\begin{thm}\label{th-periodicity}
Let $\S$ be a translation surface with the finite blocking
property. Then $\S$ is purely periodic.
\end{thm}

\begin{demo} Let $\theta$ be a direction that contains a periodic trajectory. By virtue of Proposition \ref{prop-complete-periodicity},
$\S$ admits a cylinder decomposition $(\C_i)_{i=1}^n$ in the direction $\theta$.
Starting from one of those cylinders and applying Lemma \ref{lem2cyln} step by step,
we show that all the widths of those cylinders are commensurable (remember that $\S$ is connected).
Hence there exists $w>0$ such that the width of $\C_i$ is equal to $p_i w$ with $p_i\in \N^*$.
Setting $T= w p_1 \dots p_n > 0$, we have $\phi_{\theta}^T = Id_{\S}$ a.e.

\end{demo}

\clearpage

\section{Some applications}

\begin{thm}\label{th-measure}
In genus $g\geq 2$, the set of translation surfaces that fail the
finite blocking property is of full measure in every stratum.
\end{thm}

\begin{demo}
Thanks to Proposition~\ref{prop-complete-periodicity}, it suffices to prove that completely periodic translation surfaces form a set of null measure in every stratum (in genus $g\geq 2$). This fact is well known in the community and should be attributed to the folklore.
Since it is written nowhere, we will give a proof in the Appendix.
The author would like to thank Erwan Lanneau, Samuel Leli\`evre, Anton Zorich and the referee for helpful discussions or comments about this result.

Note that Theorem \ref{th-measure} will be considerably improved in \cite{Mo-homology}.

\end{demo}

\begin{thm}\label{th-genus2}
Let $\S$ be a translation surface of genus $2$. Then $\S$ has the
finite blocking property if and only if $\S$ is a torus branched
covering.
\end{thm}

\begin{demo}

By Proposition \ref{prop-torus-branched-covering}, if $\S$ is a torus
branched covering, then $\S$ has the finite blocking property.

Conversely, suppose that $\S$ has the finite blocking property. By
Proposition \ref{prop-complete-periodicity}, $\S$ is completely
periodic. Such surfaces of genus $2$ were classified in \cite{Cal}. If
$\S$ is in $\H(2)$, then $\S$ is a Veech surface and by Proposition
\ref{prop-Veech}, $\S$ is a torus branched covering. If $\S$ is in
$\H(1,1)$, we assume by contradiction that $\S$ is not a torus branched
covering. A direct consequence of Theorem 1.2 of \cite{Cal} is that after rescaling $\S$, there
exists a square-free integer $d>0$, four positive numbers
$w_1$, $w_2$, $s_1$, $s_2$ in
$\Q(\sqrt{d})$ and a direction $\theta \in \mathbb{S}^1$ such that:

\begin{itemize}
   \item $\S$ is decomposed into cylinders $\C_1$, $\C_2$ (and sometimes $\C_3$) in the direction $\theta$
   \item $w_1$ (\resp $w_2$) is the width of $\C_1$ (\resp $\C_2$)
   \item $$\frac{w_1}{w_2}  \overline{s}_1 + \overline{s}_2 = 0$$ \\
\end{itemize}

Since $w_1$, $w_2$, $s_1$, $s_2$ are in $\Q(\sqrt{d})$, there exists two rational
numbers $r_1$ and $r_2$ such that $$\frac{w_1}{w_2} s_1 + s_2 = r_1 + r_2 \sqrt{d}.$$

The addition of those two equalities gives $\frac{w_1}{w_2} (s_1 +
\overline{s}_1) + (s_2 + \overline{s}_2) = r_1 + r_2 \sqrt{d}$.
Theorem~\ref{th-periodicity} asserts that $\S$ is purely periodic, so
$\frac{w_1}{w_2}$ is rational, so $\frac{w_1}{w_2} (s_1 +
\overline{s}_1) + (s_2 + \overline{s}_2)$ is rational, so $r_2=0$.
The subtraction of the two equalities leads to $r_1=0$.

Hence, $\frac{w_1}{w_2} s_1 + s_2$ is both null and positive,
which is impossible: $\S$ is a torus branched covering.

\end{demo}

Of course,
the notion of finite blocking property and the results presented in this paper can be translated in the vocabulary of quadratic differentials.

\clearpage

\section*{Conclusion}
We conclude this paper by the following open question:
\begin{center} Does pure periodicity imply being a torus branched covering? \end{center}

It was positively answered in \cite{Mo-V} for Veech surfaces and in the present paper for surfaces of genus two.
Proving that in full generality would provide an equivalence between the following assertions:

\begin{enumerate}

              \item $\S$ is a torus branched covering
              \item there exists $A\in GL(2,\R)$ such that $hol(H_1(A.\S,\Z))=\Z + i \Z$
              \item $Vect_{\Q}(hol(H_1(\S,\Z)))$ has dimension $2$

              \item $\S$ has the finite blocking property
              \item $\S$ has the bounded blocking property: there exists $n\in\N$ such that for every pair $(O,A)$ of points in $\S$, there exists $n$ points $B_1, \dots , B_n$ (different from $O$ and $A$) such that every geodesic from $O$ to $A$ meets one of the $B_i$'s
              \item $\S$ has the middle blocking property: For every pair $(O,A)$ of points in $\S$, the set of the midpoints of the geodesics from $O$ to $A$ is finite
              \item $\S$ has the bounded middle blocking property: there exists $n\in\N$ such that the set of the midpoints of the geodesics from $O$ to $A$ has cardinal less than $n$

              \item $\S$ is purely periodic: in every direction $\theta$ that contains a periodic orbit, the directional flow $\phi_{\theta}$ is periodic (i.e. there exists $T>0$ such that $\phi_{\theta}^T=Id_{\S}$ a.e.)
              \item in every direction that contains a periodic orbit, $\S$ is decomposable into cylinders whose widths are commensurable

\end{enumerate}

Properties $1$, $2$ and $3$ are geometric, whereas properties $4$, $5$, $6$, $7$ are exponential and properties $8$ and $9$ are dynamic.\\

If $\C$ denotes the subgroup of $H_1(\S,\Z)$ generated by the periodic orbits,
it is already possible to prove that if $\S$ is purely periodic,
then  $Vect_{\Q}(hol(\C))$ has dimension $2$.
Consequently, when $\C=H_1(\S,\Z)$,
then the nine previous assertions are equivalent \cite{Mo-homology}.

\clearpage

config: arxiv


text:
\begin{document}

\begin{frontmatter}

\selectlanguage{english}
\title{WENO schemes applied to the quasi-relativistic Vlasov--Maxwell model for laser-plasma interaction}

\selectlanguage{english}
\author[authorlabel1]{Francesco Vecil},
\ead{francesco.vecil@gmail.com}
\author[authorlabel2]{Pep Mulet Mestre},
\ead{jose.mulet@uv.es}
\author[authorlabel3]{Simon Labrunie}
\ead{simon.labrunie@univ-lorraine.fr}

\address[authorlabel1]{Laboratoire de Math\'ematiques, Universit\'e Blaise Pascal (Clermont-Ferrand 2), UMR 6620, CNRS, Campus des C\'ezeaux B.P. 80026, 63171 Aubi\`ere (France)}
\address[authorlabel2]{Universitat de Val\`encia, Departament de Matem\`atica Aplicada, calle del Doctor Moliner 50, Burjassot 46100 (Spain)}
\address[authorlabel3]{  Universit\'e de Lorraine, Institut \'Elie Cartan de Lorraine, UMR 7502, 54506 Vand{\oe}uvre-l\`es-Nancy (France)\\
  CNRS, Institut \'Elie Cartan de Lorraine, UMR 7502, 54506 Vand{\oe}uvre-l\`es-Nancy (France)}

\date{February 7, 2014}

\medskip
\begin{center}
{\small Received *****; accepted after revision +++++\\
Presented by £££££}
\end{center}

\begin{abstract}
In this paper we focus on WENO-based methods for the simulation of the
1D Quasi-Relativistic Vlasov--Maxwell (QRVM) model used to describe how a laser wave
interacts with and heats a plasma by penetrating into it.
We propose several non-oscillatory methods based on either Runge--Kutta
(explicit) or Time-Splitting (implicit) time discretizations.
We then show preliminary numerical experiments.

\vskip 0.5\baselineskip

\selectlanguage{francais}

\noindent{\bf R\'esum\'e}
\vskip 0.5\baselineskip
\noindent
{\bf Sch\'emas WENO appliqu\'es au mod\`ele Vlasov--Maxwell quasi-relativiste pour l'interaction laser-plasma. }
Dans cet article, nous nous int\'eressons aux m\'ethodes de type WENO pour la simulation du mod\`ele
Vlasov--Maxwell quasi-relativiste (QRVM) 1D, utilis\'e pour d\'ecrire la fa\c{c}on dont
une onde laser interagit avec un plasma et le r\'echauffe en le p\'en\'etrant. Nous proposons
plusieurs m\'ethodes non oscillatoires fond\'ees sur
des discr\'etisations en temps soit Runge--Kutta (explicites) soit Time-Splitting (implicites).
Ensuite, nous pr\'esentons des exp\'eriences num\'eriques pr\'eliminaires.

\keyword{Vlasov--Maxwell; WENO; laser-plasma interaction; Runge--Kutta schemes; Strang splitting}
\vskip 0.5\baselineskip
\noindent{\small{\it Mots-cl\'es~:} Vlasov--Maxwell~; WENO~;
interaction laser-plasma~; sch\'emas de Runge--Kutta~; splitting de Strang}}

\end{abstract}
\end{frontmatter}

\selectlanguage{english}

\section{Introduction}
\noindent The object of our study is the dimensionless 1D quasi-relativistic Vlasov--Maxwell (QRVM) system:
\begin{equation}
\label{QRVM}
\left\{ \begin{array}{l}
\displaystyle \frac{\partial f}{\partial t} + v(p) \frac{\partial f}{\partial x}
+ F{(t,x)} \frac{\partial f}{\partial p} = 0 \quad \mbox{(collisionless Vlasov)},
\qquad F = - \left( E + A B \right)\\[4mm]
\displaystyle  \frac{\partial E}{\partial x} = \eta^{-2} \left( \varrho_{\mathrm{ext}}- \varrho \right) \qquad \mbox{(Poisson equation)}\\[4mm]
\displaystyle \frac{\partial A}{\partial t} = -\mathcal{E}, \qquad
\frac{\partial \mathcal{E}}{\partial t} = \eta^{-2} A\varrho -\frac{\partial B}{\partial x},\qquad
\frac{\partial B}{\partial t} = -\frac{\partial \mathcal{E}}{\partial x},\qquad
\qquad \mbox{(Maxwell equations)}\\[4mm]
\displaystyle v(p) = \frac{p}{\sqrt{1+p^2}}, \qquad \varrho = \int f \, \mathrm{d}p \qquad \mbox{(relativistic character)}
\end{array} \right.
\end{equation}
solved for $(t,x,p) \in [0,+\infty[ \times [0,1] \times \mathbb{R}$,
endowed with periodic boundary conditions in $x$.
Problem \eqref{QRVM} needs several initial conditions: one for the distribution function $f$; three for the magnetic potential $A$
and its derivatives, the magnetic field $B$ and the transverse electric field $\mathcal{E}$, which
are related by
\begin{equation}
\mathcal{E} = -\frac{\partial A}{\partial t}
\qquad \mbox{and} \qquad
B = \frac{\partial A}{\partial x}.
\label{BE}
\end{equation}
The quantity $\varrho_{\text{ext}}$ represents the immobile ion background which keeps the plasma neutral.
The interest in this Vlasov--Maxwell system is motivated by its importance in plasma physics:
\begingroup
it describes laser-plasma interaction, i.e. the action of a laser wave,
called \emph{pump}, penetrating into a plasma
and heating it, while interacting with electrostatic waves and accelerating the electrons.
This model, and its variants, have been long known in the plasma physics community~\cite[and references therein]{GBS+90,Huot2003512}.
Its derivation and a discussion about the global existence and uniqueness of classical solutions can be found in \cite{CarLab06}.

\smallbreak

In order to solve~\eqref{QRVM} numerically, one has to choose a time discretization method, a Vlasov solver and a Maxwell solver.
So far, characteristic solvers have been generally used for the Maxwell part, combined with various semi-Lagrangian methods~\cite{GBS+90,Huot2003512,BoCr09} for Vlasov, as well as wavelets~\cite{BLG+08}. Time-splitting methods were often used for the quasi-relativistic model, though they are unstable with a fully relativistic model~\cite{Huot2003512}.

\smallbreak

The goal of this article is to introduce several Weighted Essentially Non-Oscillatory (WENO) schemes for the QRVM model, and to
perform preliminary tests and comparisons, in order to decide which schemes are more suitable. In Table \ref{overall} we summarize all the combinations we have considered and tested.
\endgroup

\begin{table}[ht]
\centering
\begin{tabular}{|c||l||l|}
\hline
\textbf{time integration}              & RK                            & TS\\\hline
\textbf{Vlasov equation}               & FDWENO                        & \emph{DSLWENO}\\
                                       &                               & \emph{CSLWENO}\\\hline
\textbf{Maxwell equations}             & RK                            & RK\\
                                       & LF                            & LF\\\hline
\end{tabular}
\caption{\textbf{The overall integration strategy.} The schemes in italic are implicit.}
\label{overall}
\end{table}

\noindent \textbf{RK} refers to the Total-Variation-Diminishing \textbf{R}unge--\textbf{K}utta scheme \cite{CCGMS}.\\
\textbf{TS} refers to the \textbf{T}ime-\textbf{S}plitting (Strang) scheme \cite{strang,CarVec05}.\\
\textbf{FDWENO} stands for the \textbf{F}inite-\textbf{D}ifference \textbf{W}eighted \textbf{E}ssentially \textbf{N}on \textbf{O}scillatory
interpolator for the approximation of partial derivatives \cite{CCGMS}.\\
\textbf{DSLWENO} stands for the non-conservative \textbf{D}irect \textbf{S}emi-\textbf{L}agrangian scheme \cite{CarVec05},
coupled to the Point-Value \textbf{WENO} interpolator \cite{AraBaeBelMul2011,MR2607148,CarVec05}.\\
\textbf{CSLWENO} stands for the \textbf{C}onservative \textbf{S}emi-\textbf{L}agrangian scheme,
based on the Flux-Balance-Method (FBM) \cite{FilSonBer01}, coupled to the FBM\textbf{WENO} described later on.\\
\textbf{LF} stands for \textbf{L}eap-\textbf{F}rog scheme (aka \emph{Yee} scheme).

For the sake of clarity,
we shall make use of a three-word notation to describe the coupling:
\{\emph{time discretization}\}-\{\emph{Vlasov solver}\}-\{\emph{Maxwell solver}\},
e.g., TS-DSLWENO-LF.

\smallbreak

The outline of this paper is the following: in Section \ref{numericalschemes}
we describe the initial and boundary conditions and the discretization of the system;
in Section \ref{timeintegration} we describe the time-integration strategy;
in Section \ref{numericalresults} we show numerical experiments;
and we conclude in Section~\ref{weno-is-good-for-you}.

\section{Initial and boundary conditions, and discretization}
\label{numericalschemes}

\subsection{Initialization}
\label{initialize}
\noindent Problem \eqref{QRVM} needs two initializations:
one for the distribution function $f(t,x,p)$, and one for the
electro-magnetic variables $A(t,x)$, $B(t,x)$ and $\mathcal{E}(t,x)$.

\subsubsection{Initialization of the distribution function}
\noindent We suppose that a proportion $1-\alpha$ of the
electrons are thermalized at a (dimensionless) cold velocity $v_{\mathrm{cold}}$, while
the remaining proportion $\alpha$ are hot with (dimensionless) velocity $v_{\mathrm{hot}}$
\begin{equation*}
G(p) = (1-\alpha) \, G_{\mathrm{cold}}(p) + \alpha \, G_{\mathrm{hot}}(p),
\end{equation*}
where we have split the Maxwellian $G(p)$ into
a cold part $G_{\mathrm{cold}}(p)$, described by a classical Gaussian,
and a hot part $G_{\mathrm{hot}}(p)$, described by a J\"uttner distribution:
\begin{equation*}
\underbrace{G_{\mathrm{cold}}(p)
= \frac{\exp\left(-\frac{p^2}{2 \, v_{\mathrm{cold}}^2}\right)}{\sqrt{2 \pi} \, v_{\mathrm{cold}}}}_{\mbox{normalized classical Gaussian}}
\qquad \mbox{and} \qquad
\underbrace{G_{\mathrm{hot}}(p)
= \frac{\exp\left( -\frac {\sqrt{1+p^2}-1}{v_{\mathrm{hot}}^2} \right)}
{\mathop{\mathlarger{\int}}_{\!\!\!\mathbb{R}} \exp\left( -\frac {\sqrt{1+p^2}-1}{v_{\mathrm{hot}}^2} \right) \, \mathrm{d}p}}_{\mbox{normalized J\"uttner distribution}}.
\end{equation*}

\noindent We shall introduce a fluctuation
for the initial density
\begin{equation*}
\varrho(0,x) = 1 + \frac{\varepsilon \, v_{\mathrm{cold}} \left( 0.6 \, k_{\mathrm{pla}} \right)}
{ \sqrt{1 + 3 \, v_{\mathrm{cold}}^{2}  \left( 0.6 \, k_{\mathrm{pla}} \right)^{2}} }
\cos\left( 2 \pi \, k_{\mathrm{pla}} \, x \right),
\end{equation*}
for some spatial frequency $k_{\text{pla}}$.
Consequently, a fluctuation is also introduced for the Maxwellian, hence,
all in all, the initial distribution function reads
\begin{equation*}
f(0,x,p) = \varrho(0,x)
\cdot
G(p - \varepsilon \, v_{\mathrm{cold}} \cos\left( 2 \pi \, k_{\mathrm{pla}} \, x \right)).
\end{equation*}

\subsubsection{Initialization of the electro-magnetic field}
\noindent
The initial conditions for $A$, $B$ and $\mathcal{E}$ describe the pump wave
which is going to interact with the plasma wave due to the density fluctuations.

Depending on the coupling we choose between the Vlasov and the Maxwell solvers,
we shall need to set $A$, $B$ and $\mathcal{E}$ at different
initial times and positions, which is why we keep the maximum generality
by writing them as $(t,x)$-dependent:
\begin{align*}
A\left(t,x\right)           &= A_0 \sin \left(2\pi \, k_{\mathrm{pump}} \, x - \omega_{0} \, t\right), \qquad
B\left(t,x\right)            = 2\pi \, A_0 \, k_{\mathrm{pump}} \cos \left(2\pi \, k_{\mathrm{pump}} \, x - \omega_{0} \, t \right),\\[2mm]
\mathcal{E}\left(t,x\right) &= A_0 \, \omega_{0} \cos \left(2\pi \, k_{\mathrm{pump}} \, x - \omega_{0} \, t\right).
\end{align*}
One checks that the relation between $A$, $B$ and $\mathcal{E}$ at any chosen initial time is given by \eqref{BE}.

\subsubsection{Boundary conditions}
\label{sec_boundcond}
\noindent
Problem~\eqref{QRVM} is endowed with periodic boundary conditions in the $x$-dimension.
To keep the computational domain bounded and enforce mass conservation, we use
Neumann boundary conditions in the $p$-dimension.
Actually, if the size of the domain is properly chosen, no
electrons should reach the $p$-border. The boundary conditions are implemented as:
\begin{align*}
f_{-i,j}    &= f_{i+N_{x},j}      && \mbox{for } i=1,\dots,N_{\mathrm{ghp}}, \quad j=1,\dots,N_{p}-1\\
f_{N_{x}+i,j} &= f_{i           ,j} && \mbox{for } i=1,\dots,N_{\mathrm{ghp}},  \quad j=1,\dots,N_{p}-1\\
f_{i,N_{p}+j} &= f_{i,N_{p}-1}     && \mbox{for } i=0,\dots,N_{x}, \quad j=0,\dots,N_{\mathrm{ghp}}\\
f_{i,-j}    &= f_{i,1}          && \mbox{for } i=0,\dots,N_{x}, \quad j=0,\dots,N_{\mathrm{ghp}}.
\end{align*}

\subsection{Discretization}
\label{amr}
\noindent We mesh the computational domain $\displaystyle \Omega = [0,1] \times \left[-p_{\max}, p_{\max}\right]$
by uniform grids:
\begin{equation*}
\left( x_{i}, p_{j} \right) = \left( i \, \Delta x, j \, \Delta p \right),
\qquad \!\!\!\! \left( \Delta x, \Delta p \right) = \left( \dfrac{1}{N_{x}}, \dfrac{2 \, p_{\max}}{N_{p}} \right). \nonumber
\end{equation*}
In order to take into account the boundary conditions,
ghost points outside the physical domain are used.

\section{Time integration}
\label{timeintegration}
\noindent In this section, we take care of the time integration
for the Vlasov equation
\begin{equation}
\frac{\partial f}{\partial t} + v(p) \frac{\partial f}{\partial x} + F(t,x) \frac{\partial f}{\partial p} = 0,
\qquad F = -\left( E + A \,  B \right)
\label{vlasov}
\end{equation}
and for the set of Maxwell equations
\begin{equation}
\frac{\partial A}{\partial t} = -\mathcal{E}, \qquad
\frac{\partial \mathcal{E}}{\partial t} = \eta^{-2} A\varrho -\frac{\partial B}{\partial x},\qquad
\frac{\partial B}{\partial t} = -\frac{\partial \mathcal{E}}{\partial x}.
\label{Maxwell}
\end{equation}

\noindent As for the Poisson equation
\begin{equation*}
\frac{\partial E}{\partial x} = \eta^{-2} \left( \varrho_{\mathrm{ext}}- \varrho \right),
\end{equation*}
we use the fast, spectrally-accurate solver,
whose details can be found in \cite{MulVec2012}.

\noindent We wish to test two different integration strategies, which are summarized in Table \ref{overall}.

\noindent TS is \emph{implicit} in the sense that it generally uses
implicit schemes for advection, thus weakening the constraints on the time step;
on the other hand, RK is \emph{explicit}, thus it requires a CFL condition.

\noindent This section is organized as follows:
in Section \ref{RK-schemes} we introduce the Runge--Kutta based schemes;
in Section \ref{TS-schemes} we introduce the Strang-splitting based schemes;
in Section \ref{maxwelleqs} we introduce leap-frog and multi-stage schemes to integrate \eqref{Maxwell};
in Section \ref{checasino!} we summarize all the resulting schemes.

\subsection{RK-FDWENO scheme}
\label{RK-schemes}
\noindent The explicit third-order TVD \emph{Runge--Kutta} strategy consists in
integrating, from time $t^{n}$ to $t^{n+1}$, the Vlasov equation
\begin{equation*}
\frac{\partial f}{\partial t} = - v(p) \frac{\partial f}{\partial x} - F(t,x) \frac{\partial f}{\partial p} =: \mathcal{H}\left[ t, f \right]
\end{equation*}
as
\begin{align}
f^{n,1} &= f^{n} + \Delta t \, \mathcal{H}\left[ t^{n}, f^{n} \right],
\quad
f^{n,2} = \frac{3}{4} f^{n} + \frac{1}{4} f^{n,1} + \frac{1}{4} \Delta t \, \mathcal{H}\left[ t^{n}+\Delta t, f^{n,1} \right], \nonumber\\
f^{n+1} &= \frac{1}{3} f^{n} + \frac{2}{3} f^{n,2} + \frac{2}{3} \Delta t \, \mathcal{H}\left[ t^{n}+\frac{\Delta t}{2}, f^{n,2} \right].
\label{TVD-RK3}
\end{align}

\noindent The partial derivatives are approximated through the fifth-order FDWENO routine for
finite differences, whose details can be found, for instance,
in \cite{CGMS-JCP,CGMS-JCP2} and references therein.
As this scheme is quite classical, we believe it does not deserve further details here.
The scheme is subject to a CFL constraint for stability:
\begin{equation*}
\Delta t < \dfrac{1}{ \dfrac{\left\| v(p) \right\|_{\infty}}{\Delta x} + \dfrac{\left\| F \right\|_{\infty}}{\Delta p} }.
\end{equation*}
Remark that we have to use the correct upwinding and that,
with proper boundary conditions (see Section~\ref{sec_boundcond}), the scheme enforces mass conservation.

\noindent RK requires the calculation of the Lorentz force
at three different times
\begin{equation*}
F\left( t^{n} \right) =: F^{n}, \quad
F\left( t^{n}+\Delta t \right) =: F^{n+1}, \quad
F\left( t^{n}+\frac{1}{2} \Delta t \right) =: F^{n+1/2}.
\end{equation*}
Computing the electrostatic field $E(t)$ at the desired times is easy,
because it is consistent with the distribution function $f(t)$; conversely, obtaining
the magnetic variables $A(t)$ and  $B(t)$ is slightly more complicated, because they follow
their own evolution equations. In case the time integrator for the Maxwell equations
does not provide us with $A$ and $B$ at the desired times,
we can estimate them by interpolations.

\subsection{TS-DSLWENO and TS-CSLWENO schemes}
\label{TS-schemes}
\noindent The (Strang) \emph{Time-Splitting} strategy \cite{strang,CheKno76} approximates the integration of the Vlasov equation
\begin{equation*}
\frac{\partial f}{\partial t} + v(p) \frac{\partial f}{\partial x} + F{(t,x)} \frac{\partial f}{\partial p} = 0
\end{equation*}
as a combination of partial solutions along the $x$-dimension and the $p$-dimension:
\begin{equation}
\frac{\partial f}{\partial t} + v(p) \frac{\partial f}{\partial x} = 0
\qquad \mbox{and} \qquad
\frac{\partial f}{\partial t} + F{(t,x)} \frac{\partial f}{\partial p} = 0.
\label{both}
\end{equation}
We advect $f^{n} \longmapsto f^{n+1}$ by means of the advection
field evaluated at time $t^{n+1/2}$, a strategy called \emph{prediction/correction} \cite{JakMarOwr2000,MarOwr2001},
summarized on Figure \ref{predcorr},
which gives a scheme of order~2 in time as soon as $F(t^{n+1/2})$ is approximated at order~1.
\begin{figure}[ht]
\centering
\subfigure[Prediction.]{\includegraphics[height=4.5cm,keepaspectratio]{griglia2.pdf}\label{here}}
\hspace{1cm}
\subfigure[Correction.]{\includegraphics[height=4.5cm,keepaspectratio]{griglia1.pdf}\label{nomodifydensity}}
\caption{\textbf{Prediction/correction strategy.}}
\label{predcorr}
\end{figure}

\noindent In principle, the one-dimensional PDEs \eqref{both} can be solved by means of any time integrator;
here we propose a direct semi-Lagrangian (DSL) strategy (non-conservative),
fully described in \cite{MulVec2012},
and a conservative semi-Lagrangian (CSL) strategy,
described in Section \ref{cslweno};
\emph{semi-Lagrangian} means that the method is characteristics-based.

\subsubsection{CSL integration for 1D advection problems}
\label{cslweno}
\noindent The model equation which we solve is
\begin{equation*}
\frac{\partial u}{\partial t} + \frac{\partial}{\partial x} \left[ a(t,x) \, u \right] = 0,
\qquad u(t^{\star},x) = u^{\star}(x), \qquad (t,x) \in [0,+\infty[ \times I
\end{equation*}
(being $a : [0,+\infty[ \times I \to \mathbb{R}$ and $I \subseteq \mathbb{R}$ an interval)
by means of a semi-Lagrangian conservative method; this strategy
is taken from \cite{FilSonBer01}.
To this end, we evolve approximated cell averages
\begin{equation*}
u^{n+1}_{i} \approx {\frac1{\Delta x}\int_{x_{i-1/2}}^{x_{i+1/2}} } u\left( t^{n+1}, \xi \right) \, \mathrm{d}\xi
\end{equation*}
and use a semi-Lagrangian strategy by following the characteristics backward,
along which $J \, u$ is conserved,
\begin{equation}
\int_{x_{i-1/2}}^{x_{i+1/2}} u\left( t^{n+1}, \xi \right) \, \mathrm{d}\xi
= \int_{x_{i-1/2}}^{x_{i+1/2}} u\left( t^{n}, \mathcal{X}\left( t^{n}; t^{n+1}, \xi \right) \right) \,
J\left( t^{n}; t^{n+1}, \xi \right) \, \mathrm{d}\xi,
\label{fbm}
\end{equation}
with $\mathcal{X}(s;t,x)$ the characteristic and $J(s; t, x)$ its Jacobian:
\begin{equation*}
\frac{\mathrm{d} \mathcal{X}(s; t, x)}{\mathrm{d}s} = a\left(s, \mathcal{X}(s; t, x)\right), \qquad \mathcal{X}(t;t,x) = x,
\qquad J(s; t, x) := \det \dfrac{\partial \mathcal{X}(s; t, x)}{\partial x}.
\end{equation*}
If we change variables $\displaystyle \eta = \mathcal{X}\left( t^{n}; t^{n+1}, \xi \right)$
into \eqref{fbm}, we get:
\begin{equation}
{\frac1{\Delta x}\int_{x_{i-1/2}}^{x_{i+1/2}} } u\left( t^{n+1}, \xi \right) \, \mathrm{d}\xi
= \frac{1}{\Delta x} \int_{x^{\mathrm{back}}_{i-1/2}}^{x^{\mathrm{back}}_{i+1/2}}
u\left( t^{n}, \eta \right) \, \mathrm{d}\eta
= \frac{U^{n}\left( x^{\mathrm{back}}_{i+1/2} \right) - U^{n}\left( x^{\mathrm{back}}_{i-1/2} \right)}{\Delta x},
\label{fbmscheme0}
\end{equation}
where we have set $\displaystyle x^{\mathrm{back}} := \mathcal{X}\left( t^{n}; t^{n+1}, x \right)$
and  $U^{n}$ is a primitive of $u(t^{n}, \cdot)$. This gives the following scheme:
\begin{equation}
u_{i}^{n+1}
= \frac{\widetilde U^{n}\left( x^{\mathrm{back}}_{i+1/2} \right) -
  \widetilde U^{n}\left( x^{\mathrm{back}}_{i-1/2} \right)}{\Delta x},
\label{fbmscheme}
\end{equation}
where $\widetilde U^{n}$ is an approximation of $U^{n}$ based on values of $\left( u_{j}^{n} \right)_j$.
The scheme is conservative if $u$ is compactly supported or under periodic boundary conditions.
In our application, the computations are simplified by $a$ being a real constant.
\footnote{Recall that the advection field in the $x$-dimension is independent of~$x$, and similarly in the $p$-dimension; furthermore $F(t,x)$ is approximated by $F(t^{n+1/2},x)$ on the time interval $[t^n,t^{n+1}]$.}
Therefore, we have explicit characteristics $\mathcal{X}(s; t, x) = x + a(s-t)$, so
\begin{equation*}
u^{n+1}_{i}
= \frac{\widetilde U^{n}\left( x_{i+1/2} - a \, \Delta t \right) -
  \widetilde U^{n}\left( x_{i-1/2} - a \, \Delta t \right)}{\Delta x}.
\end{equation*}

\subsubsection{The WENO reconstruction for CSL (called FBMWENO)}
\label{cslweno2}
\noindent In order to set up the scheme \eqref{fbmscheme} we need an interpolator for the primitive $U$
(dropping the time-dependency notation from now on).
In the WENO fashion, we shall perform a convex combination of several Lagrange polynomials
interpolating $U$ at different substencils.
We can adjust two parameters in order to obtain all the possible combinations:
the degree $r_{\mathrm{tot}}$ of the Lagrange polynomial interpolating $U(x)$ in the whole stencil $\mathcal{S}$
(which thus contains $r_{\mathrm{tot}}+1$ points),
and the degree $r_{\mathrm{sub}}$ of the Lagrange polynomials
in the substencils (each substencil contains $r_{\mathrm{sub}}+1$ points).
Let us also introduce the number of substencils $N_{\mathrm{sub}} := r_{\mathrm{tot}}-r_{\mathrm{sub}}+1$.

\noindent
Let us denote $P_{\nu}^{r}(x)$ the Lagrange polynomial interpolating the point values of the primitive $U$
at points $\left\{ x_{\nu-r}, \dots, x_{\nu} \right\}$.
If $\displaystyle \mathcal{S} = \left\{ x_{\mathrm{left}} , \dots, x_{\mathrm{left}+r_{\mathrm{tot}}} =: x_{\mathrm{right}} \right\}$ is
the big stencil used to approximate $U(x)$, then
\begin{equation*}
U(x) \approx \widetilde U(x):=\sum_{\ell = 0}^{N_{\mathrm{sub}}-1} \omega_{\ell}(x) \, P^{r_{\mathrm{sub}}}_{\mathrm{right}-\ell}(x).
\end{equation*}
In order to define the weights
\begin{equation*}
\omega_{\ell}(x) := \frac{\tilde\omega_{\ell}(x)}{\sum_{\ell'=0}^{N_{\mathrm{sub}}-1} \tilde\omega_{\ell'}(x)}, \qquad
\tilde\omega_{\ell}(x) := \frac{C_{\ell}(x)}{\left( 10^{-6} + \sigma_{\ell} \right)^{2}}, \qquad \ell=0, \dots, N_{\mathrm{sub}}-1
\end{equation*}
we need two ingredients: the polynomials $\left\{ C_{\ell}(x) \right\}_{\ell=0}^{N_{\mathrm{sub}}-1}$
defined by the relation
\begin{equation*}
P_{\mathrm{right}}^{r_{\mathrm{tot}}}(x) = \sum_{\ell=0}^{N_{\mathrm{sub}}-1} C_{\ell}(x) \, P_{\mathrm{right} - \ell}^{r_{\mathrm{sub}}}(x),
\end{equation*}
and the \emph{smoothness indicators} $\left\{ \sigma_{\ell} \right\}_{\ell=0}^{N_{\mathrm{sub}}-1}$,
which we wish to define in such a way that the $u^{n+1}_{i}$ given by~\eqref{fbmscheme}
is not polluted by spurious oscillations.
To this end, we are not interested in the smoothness of $U$, rather in the smoothness of $u$.

\noindent Now, the derivative of $P^r_{\nu}(x)$ is a lower-order approximation to $u(x)$:
\begin{equation*}
\mathcal{P}^r_{\nu}(x) := \frac{\mathrm{d} P^r_{\nu}}{\mathrm{d} x} (x) \approx u(x),
\end{equation*}
in the sense that if $P^r_{\nu}(x)$ approximates $U(x)$ at order $r+1$,
$\mathcal{P}^r_{\nu}(x)$ approximates $u(x)$ at order $r$.
We now fix the interval $I := \left[ x_{i-1/2}, x_{i+1/2} \right]$
that contains the evaluation point and define the smoothness measurement as
in the Jiang--Shu fashion \cite{JS}: for $\ell = 0, \dots, N_{\mathrm{sub}}-1$
\begin{equation*}
\sigma_{\ell} :=
\sum_{k=1}^{r_{\mathrm{sub}}-1} \Delta x^{2k-1} \int_{I} \left[ \left( \mathcal{P}_{\mathrm{right} - \ell}^{r_{\mathrm{sub}}} \right)^{(k)} (\xi) \right]^{2} \, \mathrm{d}\xi.
\end{equation*}
\noindent The polynomials $C_\ell(x)$ and constants $\sigma_\ell$ for $\left( r_{\mathrm{tot}}, r_{\mathrm{sub}}
\right) = (5,3)$ are given in Appendix \ref{constants_fbmweno}.

\subsection{Integration of the Maxwell equations}
\label{maxwelleqs}
\noindent We test two strategies: a \emph{leap-frog}-type Yee scheme
and a Runge--Kutta scheme. The Yee scheme will be coupled to both
schemes for the Vlasov equation and the Runge--Kutta scheme will only
be coupled to the Runge--Kutta scheme for the Vlasov equation.

In any case, once we have updated the ponderomotive force
$\displaystyle \mathcal{F} = A \, B$
up to time $t^{n}$, we impose it has numerically zero average:
\begin{equation*}
\bar{\mathcal{F}}^{n} := \frac{1}{N_{x}} \sum_{i=0}^{N_{x}-1} \mathcal{F}^{n}_{i},
\qquad \mbox{then} \qquad \mathcal{F}^{n}_{i} \longmapsto \mathcal{F}^{n}_{i} - \bar{\mathcal{F}}^{n}.
\end{equation*}

The LF scheme that we use for the Maxwell equations is second-order accurate in both space and time,
and is known as the \emph{Yee} scheme. It is of the \emph{leap-frog} type
with half-shifted variables: see Figure \ref{leapfrog} for a sketch.
\begin{figure}
\centering
\subfigure[leap $A^{n}_{i} \mapsto A^{n+1}_{i}$]{\includegraphics[scale=.4]{YEE_2.pdf}\label{yee_1}}
\subfigure[leap $B^{n}_{i+1/2} \mapsto B^{n+1}_{i+1/2}$]{\includegraphics[scale=.4]{YEE_3.pdf}\label{yee_2}}
\subfigure[leap $\mathcal{E}^{n+1/2}_{i} \mapsto \mathcal{E}^{n+3/2}_{i}$]{\includegraphics[scale=.4]{YEE_4.pdf}\label{yee_3}}
\subfigure[at time $t^{n+1}$]{\includegraphics[scale=.4]{YEE_5.pdf}\label{yee_4}}
\caption{\textbf{\emph{Leap-frog} strategy.} The scheme is second-order in both time and space,
  because all the $t$- and $x$-derivatives are approximated by centered differences.
  Inside the figure $\displaystyle \mathcal{E}_{x} := \frac{\partial \mathcal{E}}{\partial x}$
  and $\displaystyle  B_{x} := \frac{\partial B}{\partial x}$.}
\label{leapfrog}
\end{figure}
Knowing $\varrho^{n+1}$, we advance in time
$\left( A^{n}, B^{n}, \mathcal{E}^{n+1/2} \right) \longmapsto \left( A^{n+1}, B^{n+1}, \mathcal{E}^{n+3/2} \right)$
by centered finite differences:

\begin{itemize}
\item The evolution of the vector potential $A$ (Figure \ref{yee_1})
\begin{equation*}
\frac{\partial A}{\partial t} = -\mathcal{E}
\qquad \mbox{gives} \qquad
A^{n+1}_{i}=A^{n}_{i} -\mathcal{E}^{n+1/2}_{i} \Delta t.
\end{equation*}

\item The evolution of the magnetic field $B$ (Figure \ref{yee_2}),
\begin{equation*}
\frac{\partial B}{\partial t} = -\frac{\partial \mathcal{E}}{\partial x}
\qquad
\mbox{gives}
\qquad
B^{n+1}_{i+1/2}=B^{n}_{i+1/2}-\frac{\Delta
t}{\Delta x}\left( \mathcal{E}^{n+1/2}_{i+1}-\mathcal{E}^{n+1/2}_{i} \right).
\end{equation*}

\item The evolution of the transverse electric field $\mathcal{E}$ (Figure \ref{yee_3}),
\begin{equation*}
\frac{\partial \mathcal{E}}{\partial t}
= \eta^{-2} A \, \varrho -\frac{\partial B}{\partial x}
\quad
\mbox{gives}
\quad
\mathcal{E}^{n+3/2}_{i} =
\mathcal{E}^{n+1/2}_{i}+\eta^{-2} A^{n+1}_{i} \varrho^{n+1}_{i} \Delta t -
\frac{\Delta t}{\Delta x} \left( B^{n+1}_{i+1/2}-B^{n+1}_{i-1/2} \right).
\end{equation*}

\end{itemize}

\subsection{Summary of the schemes}
\label{checasino!}
\noindent In order to construct the schemes resulting from the different choices for the time integrators of
the Vlasov and the Maxwell equations (see Table \ref{overall}), we have to be particularly careful
in order to fit each block properly within the coupling.

\subsubsection{TS-DSLWENO-LF and TS-CSLWENO-LF schemes}
\noindent The scheme  to advance
\begin{equation*}
\left( f^{n}, A^{n-1/2}, B^{n-1/2}, \mathcal{E}^{n} \right) \longmapsto \left( f^{n+1}, A^{n+1/2}, B^{n+1/2}, \mathcal{E}^{n+1} \right)
\end{equation*}
is sketched on Figure \ref{TS-DSLWENO-LF}.
Notice that the time indices of $A,\ B,\ \mathcal{E}$ have been shifted by one half w.r.t. Section~\ref{maxwelleqs}, so as to have the force at hand at time~$t^{n+1/2}$, as explained in Section~\ref{TS-schemes}. Thus, $\varrho$ and~$E$ must be available at~$t^{n+1/2}$. This is done by computing them after the first half-advection in~$x$~\cite{CheKno76}, see Figure~\ref{nomodifydensity}.

The difference between the two schemes is
how the steps in Figure \ref{init_TS_LF_1.pdf} and Figure \ref{init_TS_LF_3.pdf} are performed,
with a non-conservative method for DSLWENO and with a conservative one for CSLWENO.

\begin{figure}[ht]
\centering
\subfigure[Strang 1/3]{\includegraphics[scale=.34]{init_TS_LF_1.pdf}\label{init_TS_LF_1.pdf}}
\subfigure[leap-frog]{\includegraphics[scale=.34]{init_TS_LF_2.pdf}\label{init_TS_LF_2.pdf}}
\subfigure[Strang 2/3, 3/3]{\includegraphics[scale=.34]{init_TS_LF_3.pdf}\label{init_TS_LF_3.pdf}}
\subfigure[at time $t^{n+1}$]{\includegraphics[scale=.34]{init_TS_LF_4.pdf}\label{init_TS_LF_4.pdf}}
\caption{\textbf{TS-DSLWENO-LF and TS-CSLWENO-LF schemes.} The schemes differ in how the Strang stages are performed.}
\label{TS-DSLWENO-LF}
\end{figure}

\subsubsection{RK-FDWENO-RK scheme}
This scheme is obtained by applying the third-order TVD
Runge--Kutta ODE solver \eqref{TVD-RK3} to a  discretization in $x$ and $p$ of the
Vlasov--Maxwell equations
\begin{equation*}
\frac{\partial}{\partial t} \left( \begin{array}{c} f \\\hdashline A \\ B \\ \mathcal{E} \end{array} \right)
=
\left( \begin{array}{c} - v(p) \frac{\partial f}{\partial x} - (E+A \, B)(t, x) \frac{\partial f}{\partial p} \\\hdashline -\mathcal{E} \\ -\frac{\partial \mathcal{E}}{\partial x} \\ \eta^{-2} A \, \varrho -\frac{\partial B}{\partial x} \end{array} \right)
=: \mathcal{H} \left[ t, \left( \begin{array}{c} f \\\hdashline A \\ B \\ \mathcal{E} \end{array} \right) \right],
\end{equation*}
where, as mentioned in Section \ref{RK-schemes}, the $x$ and $p$ derivatives in the Vlasov
equation are discretized by WENO finite differences, the $x$ derivatives in
the Maxwell equations are discretized by linear finite differences;
$\varrho$ is discretized by the midpoint quadrature rule,
and $E$ is computed by the Poisson solver.

\subsubsection{RK-FDWENO-LF scheme}
\noindent The resulting scheme is depicted in Figure
\ref{RK-FDWENO-LF}.
Remark that the Yee scheme forces the time step $\Delta t$ to be kept fixed,
despite the adaptive character of the Runge--Kutta scheme.

\begin{figure}
\centering
\subfigure[RK 1/3]{\includegraphics[scale=.33]{init_RK_LF_1.pdf}\label{init_RK_LF_1.pdf}}
\subfigure[leap-frog]{\includegraphics[scale=.33]{init_RK_LF_2.pdf}\label{init_RK_LF_2.pdf}}
\subfigure[RK 2/3]{\includegraphics[scale=.33]{init_RK_LF_3.pdf}\label{init_RK_LF_3.pdf}}
\subfigure[RK 3/3]{\includegraphics[scale=.33]{init_RK_LF_4.pdf}\label{init_RK_LF_4.pdf}}
\caption{\textbf{RK-FDWENO-LF scheme.} This scheme is second-order in time (because of
the interpolation and the first-order approximation of $\varrho$ used to evolve $\mathcal{E}$) and second-order in space.}
\label{RK-FDWENO-LF}
\end{figure}

\section{Results for the quasi-relativistic Vlasov--Maxwell system}
\label{numericalresults}
\noindent No WENO-based scheme has yet been extensively tested on the QRVM problem. Therefore,
our first task is to decide which among the overall integration strategies introduced in Table~\ref{overall}
are suitable.

\subsection{Empirical stability results}
\noindent All the schemes proposed in this article seem stable
from empirical observation, but RK-FDWENO-LF requires
extremely small time steps in order not to blow up.
A summary is given in Table \ref{estabilidad}.

\begin{table}[ht]
\centering
\begin{small}
\begin{tabular}{|l||c|c|}
\hline
Vlasov $\downarrow$ \hspace{3cm} Maxwell $\rightarrow$  & \textbf{LF}                                   & \textbf{RK}\\\hline\hline
\textbf{RK}-\textbf{FDWENO}-                            & \scalebox{3}{\frownie}                        & \scalebox{3}{\smiley}\\\hline
\textbf{TS}-\textbf{DSLWENO}-                           & \scalebox{3}{\frownie}                        & not couplable\\\hline
\textbf{TS}-\textbf{CSLWENO}-                           & \scalebox{3}{\smiley}                         & not couplable\\\hline
\end{tabular}
\end{small}
\caption{\textbf{Quality of the results.}}
\label{estabilidad}
\end{table}

\noindent
The evolution equations for $B$ and $\mathcal{E}$ can be rewritten as
\begin{equation*}
\frac{\partial \left( B \pm \mathcal{E} \right)}{\partial t}
\pm \frac{\partial \left( B \pm \mathcal{E} \right)}{\partial x}
= \pm \eta^{-2} \, A \, \varrho,
\end{equation*}
therefore the condition
\begin{equation}
\Delta t < \Delta x
\label{perlastabilita}
\end{equation}
seems reasonable as constraint for stability of an explicit scheme.

\noindent If we take as reference a $400 \times 400$ mesh, $\Delta x$ would be equal to $0.0025$.
Notwithstanding, experiments suggest the threshold $\Delta t$ should be of order $10^{-5}$ for RK-FDWENO-LF.
In the other cases, the RK-FDWENO-RK scheme, the TS-DSLWENO-LF scheme and the TS-CSLWENO-LF scheme,
if the CFL parameter or the $\Delta t$
are adapted so as to fulfill \eqref{perlastabilita}, the simulations appear stable.

\subsection{Quality of the results}
\noindent
On Figure \ref{evol_compared} we compare at similar stages the evolution
computed by the three most stable schemes. The dynamic of laser-plasma interaction~\cite{GBS+90,Huot2003512,BLG+08} is precisely captured.
The plasma wave, initiated by the initial fluctuations of the electron density,
exchanges energy with the electrons and with the transverse electromagnetic wave.
Vortices appear in phase space, due to the particles
getting trapped by the plasma wave's potential well and bouncing on its separatrices. The vortices
show an oscillating behavior: they periodically inflate and deflate.
One observes the well-known ``filamentation'' phenomenon: thin
structures appear, then they are stretched thinner and folded, again
and again.

We see that in the short term both the RK-based and the TS-based schemes behave well,
but TS-CSLWENO-LF diffuses the microscopic details more than RK-FDWENO-RK, as the long-time behavior ($t = 300$) shows.

On Figure \ref{toutensemble} we plot the conservation properties:
the relative variation (w.r.t. time $t=0$) of the mass, of the $L^2$-norm and of the total energy
\begin{equation}
W(t)=\overbrace{\underbrace{\frac{1}{2} \int_0^1 A^2 \varrho \, \mathrm{d}x}_{\mathrm{WTK}(t) := \mathrm{kinetic}}
+ \underbrace{\frac{\eta^{2}}{2} \int_0^1 \left[ \mathcal{E}^2 + B^2 \right] \mathrm{d}x}_{\mathrm{WTP}(t) := \mathrm{potential}}}^{\mathrm{WT}(t) := \mathrm{transversal}}
+
\overbrace{\underbrace{\int_0^1 \int_{\mathbb{R}} \left( \sqrt{1+p^2} -1 \right) f \,\mathrm{d}p \, \mathrm{d}x}_{\mathrm{WLK}(t) := \mathrm{kinetic}} + \underbrace{\frac{\eta^{2}}{2} \int_0^1 E^2 \,\mathrm{d}x}_{\mathrm{WLP}(t) := \mathrm{potential}}}^{\mathrm{WL}(t) := \mathrm{longitudinal}},
\label{totalenergy}
\end{equation}
which is shown in~\cite{CarLab06} to be conserved by the system.
We observe that around time 300 TS-DSLWENO-LF has gained about 13 \% w.r.t. the normalized mass,
which means that the plasma is strongly non-neutral, hence even the integration of the
Poisson equation becomes meaningless because the periodicity is lost.
RK-FDWENO-RK conserves better the $L^{2}$-norm, i.e. the microscopic details inside the computational domain, and the total energy.

\begin{figure}[ht]
\centering{\hspace*{-4em}  \includegraphics[scale = 1]{evol_compared.pdf} }
  \caption{
    \textbf{Comparison.} Evolution of the system up to time $\approx 300$, for a $400 \times 400$ mesh, using three different schemes.
    Left column: the explicit conservative RK-FDWENO-RK.
    Central column: the implicit non-conservative TS-DSLWENO-LF.
    Right column: the implicit conservative TS-CSLWENO-LF.
  }
  \label{evol_compared}
\end{figure}

\begin{figure}[ht]
  \centering
  \includegraphics[scale = 1]{toutensemble.pdf}
  \caption{
    \textbf{Conservation properties.}
    Top: the relative variation of total mass w.r.t. the initial condition.
    Center: the relative variation of the $L^2$-norm w.r.t. the initial condition.
    Bottom: the relative variation of the total energy \eqref{totalenergy} w.r.t. the initial condition.
  }
  \label{toutensemble}
\end{figure}

\section{Conclusion}
\label{weno-is-good-for-you}
We have performed some preliminary tests of several WENO-based schemes to simulate the 1D quasi-relativistic Vlasov--Maxwell system, which models laser-plasma interaction.
WENO schemes, with their high accuracy and robustness to the steep gradients created by filamentation, are ideally suited to capture the dynamic of this interaction.
Indeed, our test cases have reproduced the qualitative behavior known from the literature since~\cite{GBS+90}.

To decide which schemes are more suitable for the simulation of the QRVM problem,
we tested the various combinations of Table~\ref{overall}.
Some of them immediately appear unsatisfactory, either because they require ridiculously small
time steps, or because they are strongly non-conservative.
The two strategies which show the best behavior are RK-FDWENO-RK and TS-CSLWENO-LF, which
are both conservative; the advantage of TS-CSLWENO-LF is its implicit character and weaker
constraints on the time step, while its drawback is that
in the long time it shows a more diffusive behavior.

{F}rom the computational point of view, WENO-based schemes have several other advantages.
They are easily parallelizable: see for instance \cite{ManCac09} for a parallel version of RK-FDWENO.
They can be made adaptive relatively easily:
see \cite{MulVec2012} for an AMR version of TS-DSLWENO, or
\cite[and references therein]{BarMarMul2011,BaeMul2006} for an AMR version of RK-FDWENO;
the built-in computation of smoothness indicators points to the regions which have to be refined (or de-refined).
This will be presented in a future publication.

\clearpage

\appendix

\section{Constants}
\label{Constants}
\noindent The constants involved in the dimensionless system are:
\begin{eqnarray*}
&& \eta = \frac{3}{10 \pi}, \qquad
k_{\mathrm{pump}} = 4, \qquad
k_{\mathrm{pla}} = 2, \qquad
\omega_{0} = \sqrt{\eta^{-2} +k_{\mathrm{pump}}^2 }, \qquad
A_0 = \frac{2.5}{\omega_0},\\
&& \alpha = 0.05, \qquad
v_{\mathrm{cold}} = \sqrt{\frac{15}{511}}, \qquad
v_{\mathrm{hot}} = \sqrt{\frac{100}{511}}, \qquad
\varepsilon = \frac{\sqrt{2}}{10}, \qquad
p_{\max} = 8.
\end{eqnarray*}

\section{Constants for FBMWENO}
\label{constants_fbmweno}
\noindent If we let $x \in \left] x_{i-1/2}, x_{i+1/2} \right[$ and the interpolant is centered in the stencil,
\begin{equation*}
\mathrm{left} = i-5/2, \qquad \mathrm{right} = i+5/2, \qquad N_{\mathrm{sub}} = 3.
\end{equation*}
The polynomials $\left\{ C_{\ell}(x) \right\}_{\ell=0}^{N_{\mathrm{sub}}-1}$ are
\begin{footnotesize}
\begin{equation*}
C_{0} =  \frac{\left( x - x_{i-5/2} \right)\left( x - x_{i-3/2} \right)}{20 \, \Delta x^{2}}, \quad
C_{1} = -\frac{\left( x - x_{i-5/2} \right)\left( x - x_{i+5/2} \right)}{10 \, \Delta x^{2}}, \quad
C_{2} =  \frac{\left( x - x_{i+5/2} \right)\left( x - x_{i+3/2} \right)}{20 \, \Delta x^{2}}.
\end{equation*}
\end{footnotesize}
The smoothness indicators are
\begin{footnotesize}
\begin{align*}
\sigma_{0} &=
\frac{10}{3} \left( U_{i-1/2} \right)^{2} - 17 \, U_{i-1/2} \, U_{i+1/2} + 14 \, U_{i-1/2} \, U_{i+3/2} -\frac{11}{3} U_{i-1/2} \, U_{i+5/2}
+ 22 \left( U_{i+1/2} \right)^{2} \\& - \frac{111}{3} U_{i+1/2} \, U_{i+3/2} + 10 \, U_{i+1/2} \, U_{i+5/2}
+ 16 \left( U_{i+3/2} \right)^{2} - 9 \, U_{i+3/2} \, U_{i+5/2}
 +  \frac{4}{3} \left( U_{i+5/2} \right)^{2},\\[3mm]
\sigma_{1} &=
  \frac{4}{3} \left( U_{i-3/2} \right)^{2} -7 \, U_{i-3/2} \, U_{i-1/2} + 6 \, U_{i-3/2} \, U_{i+1/2} - \frac{5}{3} U_{i-3/2} \, U_{i+3/2}
  + 10 \left( U_{i-1/2} \right)^{2} \\&- 19 \, U_{i-1/2} \, U_{i+1/2} + 6 \, U_{i-1/2} \, U_{i+3/2}
  + 10 \left( U_{i+1/2} \right)^{2} - 7 \, U_{i+1/2} \, U_{i+3/2} + \frac{4}{3} \left( U_{i+3/2} \right)^{2},\\[3mm]
\sigma_{2} &=
\frac{10}{3} \left( U_{i+1/2} \right)^{2}
- 17 \, U_{i+1/2} \, U_{i-1/2}
+ 14 \, U_{i+1/2} \, U_{i-3/2}
-\frac{11}{3} U_{i+1/2} \, U_{i-5/2}
+ 22 \left( U_{i-1/2} \right)^{2}
\\&- \frac{111}{3} U_{i-1/2} \, U_{i-3/2}
 + 10 \, U_{i-1/2} U_{i-5/2}
+ 16 \left( U_{i-3/2} \right)^{2} -9 \, U_{i-5/2} \, U_{i-3/2} +
  \frac{4}{3} \left( U_{i-5/2} \right)^{2}.
\end{align*}
\end{footnotesize}

\section*{Acknowledgments}
\noindent Francesco Vecil and Pep Mulet acknowledge financial support
from MINECO project MTM2011-22741.


text:
\section{Introduction}
We consider in this paper  nonlinear regression model of the form
\begin{equation}
\label{model}Y_n = f(X_n,\theta) + \varepsilon_n, \ \ n\in\bkN,
\end{equation}
where the observations $(X_n,Y_n)_n$ are independent random vectors in $\mathbb{R}^p\times \mathbb{R}$, $(\varepsilon_n)_n$ are independent, identically distributed zero-mean non observable random variables. Moreover, $f: \mathbb{R}^p\times \mathbb{R}^q \longrightarrow \mathbb{R}$
and $\theta \in \mathbb{R}^q$ is the unknown
parameter to estimate.
\vspace{1ex}

Nonlinear regression models are a standard tool for modeling real phenomena
with a complex dynamic.
It has a wide range of applications including maching learning (\cite{bach2014adaptivity}), ecology (\cite{komori2016asymmetric}), econometry (\cite{varian2014big}), pharmacokinetics (\cite{Bates88}), epidemiology (\cite{suarez2017}) or biology (\cite{Giurcuaneanu05}).
Most of the time, the parameter  $\theta$ is estimated using
the least squares method, thus it is estimated by

\[ \widehat\theta_n = \arg\min_{h \in \bkR^q} \sum_{j=1}^n (Y_j - f(X_j, h))^2 \]

Many authors have studied the asymptotic behavior of the least squares estimator,
under various assumptions with different methods.
For example, \cite{Jennrich69} consider the case when $\theta$ belongs to a compact set,
$f$ is a continuous function with a unique extremum.
\cite{Wu81} consider similar local assumptions.
\cite{lai94, Skouras, Geer90, Geer96, Yao00} generalized the consistency results.
Under a stronger assumption about the errors $\varepsilon$,
that is, considering that the errors are uniformly subgaussian,
\cite{Geer90} obtained sharper stochastic bounds using empirical process methods. Under second moment assumptions on the errors and regularity assumptions (lipschitz conditions) on $f$, \cite{Pollard06} established weak consistency and a central
limit theorem.
\cite{Yang} consider a compact set for $\theta$,
under strong regularity assumptions and the algorithm allows to reach a stationary point without certainty that it is the one we are looking for.
They built hypothesis tests and confidence intervals.
Finally, let also mention \cite{Yao00}
which considers the case of stable nonlinear autoregressive models
and establishes the strong consistency of the least squares estimator.
\vspace{1ex}

However, in practice, the calculation of the least square estimator is not explicit in most cases and therefore requires the implementation of a deterministic approximation algorithm. The second order Gauss-Newton algorithm is generally used  (or sometimes Gauss-Marquard algorithm to avoid inversion problems
(see for example \cite{Bates88}).

 These algorithms are  therefore not adapted to the case where the data are acquired sequentially and at high frequencies.
 In such a situation, stochastic algorithms
offer an interesting alternative.
One example is the stochastic gradient algorithm,
defined by
\begin{equation}
\theta_n = \theta_{n-1} - \gamma_n \nabla_\theta f(X_n, \theta_{n-1}) \pa{Y_n - f(X_n, \theta_{n-1})},
\end{equation}
where  $(\gamma_n)$ is a sequence of  positive real numbers decreasing towards 0. Thanks to its recursive nature,
this algorithm does not require to store all the data and can be updated automatically when the data sequentially arrive.

It is thus adapted to very large datasets of data arriving at high speed.
We refer to \cite{robbins1951} and to its averaged version \cite{PolyakJud92}
and \cite{ruppert1988efficient}.
However, even if it is often used in practice as in the case of neural networks,
the use of stochastic gradient algorithms can lead to non satistying results. Indeed, for such models, it amounts to take the same step sequence for each coordinates. Nevertheless, as explained in \cite{BGBP2019}, if the Hessian of the function we would like to minimize has eigenvalues at different scales,
this kind of "uniform" step sequences is not adapted.
\vspace{1ex}

In this paper, we propose an alternative strategy to stochastic gradient algorithms,
in the spirit of the Gauss-Newton algorithm.
It is defined by :
\begin{eqnarray}
\phi_{n} & = & \nabla_\theta f(X_{n}, \theta_n)
\label{GNun} \\
S_{n}^{-1} &= & S_{n-1}^{-1} - (1 + \phi_{n}^T  S_n^{-1} \phi_{n})^{-1}
S_{n-1}^{-1} \phi_{n} \phi_{n}^T S_{n-1}^{-1}
\label{GNdeux} \\
\theta_{n} &= &\theta_{n-1} + S_{n}^{-1}\phi_{n}\pa{Y_{n} -f(X_{n}, \theta_{n-1})}
\label{GNtrois}
\end{eqnarray}
where the initial value $\theta_0$ can be arbitrarily chosen and $S_0$ is a positive definite deterministic matrix,
typically $S_0 = I_q$ where $I_q$ denotes the identity matrix of order $q$.
Remark that thanks to Riccati's formula (see \cite[p. 96]{Duf97}, also called Sherman Morrison's formula (\cite{Sherman50}), $S_n^{-1}$ is the inverse of the matrix $S_{n}$ defined by

\[S_n = S_0 + \sum_{j=1}^{n}\phi_j \phi_j^T.\]

When the function $f$ is linear of the form $f(x,\theta) = \theta^T x$,
 algorithm (\ref{GNun})-(\ref{GNtrois})  rewrites as the standard recursive least-squares estimators
(see \cite{Duf97}) defined by:

\begin{eqnarray*}
S_n^{-1} &= & S_{n-1}^{-1} - (1 + X_n^T  S_{n-1}^{-1} X_n)^{-1}
S_{n-1}^{-1} X_nX_n^T S_{n-1}^{-1} \\
\theta_n &= &\theta_{n-1} + S_{n}^{-1}X_n\pa{Y_n - \theta_{n-1}^T X_n}
\end{eqnarray*}

This algorithm can be considered as a Newton stochastic algorithm since the matrix $n^{-1} S_n$ is an estimate of the Hessian matrix of the least squares criterion.
\vspace{1ex}

To the best of our knowledge and apart from the least squares estimate mentioned above,
second order stochastic algorithms are hardly ever used and studied since they often require
the inversion of a matrix at each step, which can be very expensive in term of time calculation.
To overcome this problem some authors (see for instance \cite{mokhtari2014res,lucchi2015variance,byrd2016stochastic})
use the BFGS (for \emph{Broyden-Fletcher-Goldfarb-Shanno}) algorithm
which is based on the recursive estimation of a matrix
whose behavior is closed to the one of the inverse of the Hessian matrix.
Nevertheless, this last estimate need a regularization of  the objective function,
leading to unsatisfactory estimation of the unknown parameter.
In a recent paper  dedicated to estimation of parameters in logistic regression models  \citep{BGBP2019},
the authors propose a truncated Stochastic Newton algorithm. This truncation opens the way for online stochastic Newton algorithm without necessity to penalize the objective function.

In the same spirit of this work, and to relax assumption on the function $f$,
we  consider in fact a modified version  of the
 stochastic Gauss-Newton algorithm defined by \eqref{GNtrois},
that enables us to obtain the asymptotic efficiency of the estimates
in a larger  area of assumptions.

In addition, we introduce the following new \emph{Averaged Stochastic Gauss-Newton} algorithm (ASN for short)
defined by

\begin{eqnarray}
\phi_n & = & \nabla_h f(X_n, \overline{\theta}_{n-1})
\label{AGNun} \\
S_n^{-1} &= & S_{n-1}^{-1} - (1 + \phi_n^T  S_{n-1}^{-1} \phi_n)^{-1}
S_{n-1}^{-1} \phi_n \phi_n^T S_{n-1}^{-1}
\label{AGNdeux} \\
\theta_n &= &\theta_{n-1} + n^\beta S_n^{-1}\phi_n\pa{Y_n -f(X_n, \theta_{n-1})}
\label{AGNtrois} \\
\overline{\theta}_n &=& \overline{\theta}_{n-1}
\label{AGNquatre}  + \dfrac{1}{n}\,\pa{\theta_n - \overline{\theta}_{n-1}}
\end{eqnarray}

where $\beta\in (0, 1/2)$, $\overline{\theta}_0 = 0$.

The introduction of the term $n^\beta$ before the term $S_n^{-1}$  in (\ref{AGNtrois}) allows the algorithm to move quickly which enables to reduce the sensibility to a bad initialization.
The averaging step allows to maintain an optimal asymptotic behavior.
Indeed, under assumptions, we first give the rate of convergence of the estimates, before proving their asymptotic efficiency.
\vspace{1ex}

The paper is organized as follows.
Framework and algorithms are introduced in Section~\ref{sec:framework}.
In Section~\ref{sectiontheo}, we give the almost sure rates of convergence of the estimates and establish their asymptotic normality. A simulation study illustrating the interest of averaging is presented in Section~\ref{sectionsimu}. Proofs are postponed in Section~\ref{sectionproof} while some general results used in the proofs on almost sure rates of convergence for martingales  are given in Section~\ref{sectionmartingales}.


TITLE: Is every point in the Cantor set rational?
QUESTION [0 upvotes]: I saw a post somewhere that stated that "each point in the Cantor set is the endpoint of some closed interval of the form $\big[\frac{n}{3^k}, \frac{n+1}{3^k}\big]$, with $n, k \in \mathbb{N}$, so each point in the Cantor set is rational." Is that true?

REPLY [1 votes]: Note that $$\frac{1}{4}\in C$$ is a rational but not an endpoint.
(plus those uncountably many irrationals hiding in the corners...)