praticals_CR.tex

\documentclass[11pt,a4paper]{article}

% ADDED PACKAGE
\usepackage{latexsym}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{graphicx}   %  include figure files
\usepackage{epsfig}       %  files .eps
\usepackage{dcolumn}   %  align tble columns on decimal points
\usepackage{bm}           %  bold math
\usepackage{color}        %  textcolor
\usepackage{ulem}        %  sout (strike out)

\usepackage{tabularx}
\usepackage[T1]{fontenc}
\usepackage{hyperref}
\PassOptionsToPackage{unicode}{hyperref}
\PassOptionsToPackage{naturalnames}{hyperref}

\usepackage{times}

\usepackage{wrapfig}

\usepackage[noanswer]{exercise}

\newcommand{\rom}[1]{\uppercase\expandafter{\romannumeral #1\relax}}

\renewcommand{\QuestionNB}{\arabic{Exercise}.\arabic{Question} }

\renewcommand{\abstractname}{Summary}

\renewcommand{\ExerciseListName}{Exercise~}

\renewcommand\ExerciseListHeader{%
    \bfseries\large\ExerciseListName\ExerciseHeaderNB : \ExerciseHeaderTitle\ExerciseHeaderOrigin%
\addcontentsline{toc}{subsubsection}{\ExerciseListName\ExerciseHeaderNB : \ExerciseHeaderTitle\ExerciseHeaderOrigin}
}


\usepackage[top=40pt,bottom=50pt,left=40pt,right=40pt]{geometry}

\usepackage{filecontents,listings}

 
\makeatletter
\def\lst@MSkipToFirst{%
    \global\advance\lst@lineno\@ne
    \ifnum \lst@lineno=\lst@firstline
        \def\lst@next{\lst@LeaveMode \global\lst@newlines\z@
        \lst@OnceAtEOL \global\let\lst@OnceAtEOL\@empty
        \lst@InitLstNumber % Added to work with modified \lsthk@PreInit.
        \lsthk@InitVarsBOL
        \c@lstnumber=\numexpr-1+\lst@lineno % this enforces the displayed line numbers to always be the input line numbers
        \lst@BOLGobble}%
        \expandafter\lst@next
    \fi}
\makeatother
%\let\oldlstinputlisting\lstinputlisting
%\renewcommand{\lstinputlisting}[2][]{%
%\oldlstinputlisting[caption={file \texttt{\detokenize{#2}}},#1]{#2}%
%}


\usepackage{color}

\definecolor{mygreen}{rgb}{0,0.6,0}
\definecolor{mygray}{rgb}{0.5,0.5,0.5}
\definecolor{mymauve}{rgb}{0.58,0,0.82}
\definecolor{backcolour}{rgb}{0.95,0.95,0.92}
\definecolor{classcolour}{rgb}{0.2,0.2,0.9}
\definecolor{paramcolour}{rgb}{0.9,0.2,0.2}


\lstset{ %
  backgroundcolor=\color{backcolour},   
  basicstyle=\linespread{0.85}\ttfamily\footnotesize,        % the size of the fonts that are used for the code
  breakatwhitespace=false,         % sets if automatic breaks should only happen at whitespace
  breaklines=true,                 % sets automatic line breaking
  captionpos=b,                    % sets the caption-position to bottom
  commentstyle=\color{mygreen},    % comment style
  escapeinside={\%*}{*)},          % if you want to add LaTeX within your code
  extendedchars=true,              % lets you use non-ASCII characters; for 8-bits encodings only, does not work with UTF-8
  keepspaces=true,                 % keeps spaces in text, useful for keeping indentation of code (possibly needs columns=flexible)
  keywordstyle=\color{blue},       % keyword style
  language=Python,                 % the language of the code
  otherkeywords={*,...},            % if you want to add more keywords to the set
  numbers=left,                    % where to put the line-numbers; possible values are (none, left, right)
  numbersep=5pt,                   % how far the line-numbers are from the code
  numberstyle=\tiny\color{black}, % the style that is used for the line-numbers
  rulecolor=\color{black},         % if not set, the frame-color may be changed on line-breaks within not-black text (e.g. comments (green here))
  showspaces=false,                % show spaces everywhere adding particular underscores; it overrides 'showstringspaces'
  showstringspaces=false,          % underline spaces within strings only
  showtabs=false,                  % show tabs within strings adding particular underscores
  stepnumber=1,                    % the step between two line-numbers. If it's 1, each line will be numbered
  stringstyle=\color{mymauve},     % string literal style
  tabsize=2,	                   % sets default tabsize to 2 spaces
  morekeywords={Main,Species,DiagScalar,DiagParticleBinning,DiagFields},
  keywordstyle=\color{magenta},
  keywords=[2]{time_frozen, grid_length, simulation_time, number_density,cell_length,axes, particles_per_cell, timestep},
  keywordstyle=[2]\color{classcolour},
  keywords=[3]{xnumber,xamplitude},
  keywordstyle=[3]\color{paramcolour}
}

\newcommand{\class}[1] {{\color{magenta}\texttt{#1}}}
\newcommand{\code}[1] {{\color{classcolour}\texttt{#1}}}
\newcommand{\param}[1] {{\color{paramcolour}\texttt{#1}}}


% NEW COMMAND
\newcommand{\vE}{\mathbf{E}}
\newcommand{\vB}{\mathbf{B}}
\newcommand{\vJ}{\mathbf{J}}
\newcommand{\vx}{\mathbf{x}}
\newcommand{\vp}{\mathbf{p}}
\newcommand{\vv}{\mathbf{v}}
\newcommand{\vu}{\mathbf{u}}
\newcommand{\vF}{\mathbf{F}}

\usepackage[multiple]{footmisc}


\makeatletter
\let\@fnsymbol\@arabic
\makeatother

\usepackage{etoolbox}
\makeatletter
\patchcmd{\@maketitle}{\LARGE \@title}{\fontsize{24}{19.2}\selectfont\textbf{\@title}}{}{}
\makeatother

\renewcommand{\contentsname}{Summary}

\newcommand{\addsection}[1] {\newpage
\section*{#1}
\addcontentsline{toc}{section}{#1}}

\newcommand{\addsubsection}[1] {
\subsection*{#1}
\addcontentsline{toc}{subsection}{#1}}



\begin{document}
%-------
% Title
%-------
\title{UPMC - TP plasma: \\ \vspace{10pt} Introduction to plasma physics \\ \vspace{10pt} with Particle-In-Cell simulations}


\author{Caterina~Riconda \\ {\tiny\url{caterina.riconda@upmc.fr}} \and Tommaso Vinci \\ {\tiny\url{tommaso.vinci@polytechnique.edu}} \and Mickael~Grech \\ {\tiny\url{mickael.grech@polytechnique.edu}}}

\date{}

\maketitle              

\hrule
\vspace*{20pt}

These exercises are  intended as an introduction to plasma physics through Particle-in-cell (PIC) simulations.

The code used is Smilei and is currently developed at LULI in collaboration with various laboratories 


%\begin{center}
%\includegraphics[width=0.5\textwidth]{SmileiLogo.pdf}
%\end{center}

\tableofcontents
\subsection*{Links:}
\begin{tabularx}{0.95\textwidth}{rl}
These notes: & \url{https://github.com/SmileiPIC/TPUPMC} \\
Smilei homepage: & \url{https://Smileipic.github.io/Smilei}\\
Smilei download: &\url{https://github.com/SmileiPIC/Smilei}
\end{tabularx}


%%%%%%%%%%%%%%%%%%
%  INTRODUCTION  %
%%%%%%%%%%%%%%%%%%
\newpage
\addsection{Introduction}\label{intro}

\addsubsection{The Particle-In-Cell code Smilei}

Smilei is a Particle-In-Cell (PIC) code initially developed to study laser-plasma interaction at extreme intensities. 
The PIC method aims at solving the Maxwell-Vlasov system of  equations~\cite{birdsall_langdon}.

The electromagnetic fields ($\vE,\vB$) are described on a (here cartesian) grid, and Maxwell's equations:
\begin{eqnarray}\label{eq_Maxwell}
\nabla \cdot \vE &=& \frac{\rho}{\epsilon_0} \,, \nonumber \\
\nabla \cdot \vB &=& 0 \,,\\
\nabla \times \vE &=& -\partial_t \vB \,, \nonumber\\
\nabla \times \vB &=& \mu_0\, \vJ + \mu_0 \epsilon_0\,\partial_t \vE \,, \nonumber
\end{eqnarray}
are solved using the finite-domain time-difference (FDTD) method~\cite{taflove_2005}.
In Eqs.~\eqref{eq_Maxwell}, $\epsilon_0$ and $\mu_0$ are the vacuum permittivity and permeability, respectively. 
As for $\rho$ and $\vJ$, they denote the total charge density and total current, respectively.

Both quantities $\rho$ and $\vJ$ are due to the presence of a plasma.
The kinetic description of a collisionless (fully ionised) plasma is encompassed by the (here relativistic) Vlasov equation, which for each species $s$ of the plasma, reads:
\begin{eqnarray}\label{eq_Vlasov}
\left(\partial_t  + \frac{\vp}{m_s \gamma} \cdot \nabla + \vF_L \cdot \nabla_{\vp} \right) f_s = 0\,,
\end{eqnarray}
where $f_s(t,\vx,\vv)$ is the species distribution function, $\vF_L$ is the external force acting on the species (in our case the Lorentz force), $m_s$ is the mass of one (physical) element of species $s$ (i.e. the mass of one ion or one electron of the plasma) and $\gamma = \sqrt{1+\vp^2/(m_s^2\,c^2)}$ is the so-called Lorentz factor of an element of species $s$ with momentum~$\vp$.

The PIC method consists in solving the Vlasov Eq.~\eqref{eq_Vlasov} by discretizing the distribution function as a sum of $N_s$ so-called macro-particles:
\begin{eqnarray}\label{eq_macrop}
f_\alpha (t,\vx,\vv) = \sum_{j=1}^{N_s}\,\alpha_s^j\,\,S(\vx-\vx^j)\,\delta(\vv-\vv^j)\,,
\end{eqnarray}
where $\vx^j$ and $\vp^j$ are the position and momentum of the $j^{th}$ macro-particle, respectively, $\alpha^j$ is the so-called particle-weight, $S(\vx)$ is the so-called shape of the macro-particle [note that $S(\vx)$ is defined such that $S(-\vx)=S(\vx)$ and $\int d\vx\,S(\vx)=1$], and $\delta(\vp)$ is the Dirac-distribution. Introducing Eq.~\eqref{eq_macrop} in Vlasov Eq.~\eqref{eq_Vlasov} simplifies the problem to solving, for each macro-particle of each species, its relativistic equations of motion:
\begin{eqnarray*}
\frac{d\vx}{dt}  &=& \frac{\vp}{m_s \gamma}\,,\\
\frac{d\vp}{dt}  &=&  q\,\left( \vE^j + \frac{\vp}{m_s \gamma} \times \vB^j \right)\,,
\end{eqnarray*}
where we have explicitly written the Lorentz force acting on the $j^{th}$ particle with charge $q$ and mass $m_s$, and $\vE^j$ and $\vB^j$ are the electric and magnetic fields seen by this particle, respectively.

\addsubsection{Fluid equations}

The fluid equations studied in the course can be deduced by the kinetic description, and are thus contained in the kinetic description 
in the limit where kinetic effects are not dominant, and the particle distribution function can be described as a local Maxwellian. 

The connection between the kinetic description and the fluid description is given by 
\begin{eqnarray*}
n_s &=& \int f_s d^3p \\
vu_s &=&\int \vv f_s d^3p \\
k_B T_s &=&\int \frac{p^2}{2m_s}  f_s d^3p
\end{eqnarray*}

We remind their form here, for  $s=i,e$ :

\begin{eqnarray*}
\frac{\partial  n_s}{\partial t}+
 {\bf{\nabla}} \cdot (n_s \vu_s)  &=&  0\, \\
m_s n_s \left[\frac{\partial \vu_s}{\partial t}+\vu_s \cdot {\bf{\nabla}} \vu_s
\right] &=& q_s n_s \left( \vE + \vu_s \times \vB \right) - {\bf{\nabla}} p_s \,  \\
p_s &=& n_s T_s \, \\  
T_s &=& C_s n_s^{\gamma_s-1}
\end{eqnarray*}

\addsubsection{Normalizations}

The PIC code Smilei is an electromagnetic PIC code. As such, it solves Maxwell's equations in which a characteristic velocity, that of light $c = (\epsilon_0 \mu_0)^{-1/2}$ appears. It is therefore natural to normalise velocities in the code to $c$. Similarly, it is interesting to normalise the charge of all particles in the plasma to $e$ (with $-e$ the charge of the electron) and all masses to $m_e$ the electron mass.

To proceed with our normalization, we now need to define a characteristic time. When considering a plasma density with constant density $n_{e0} = Z\,n_{i0}$ ($Z$ being the charge state of one ion in the plasma), a characteristic frequency arises that is the so-called electron-plasma frequency $\omega_{pe} = \sqrt{e^2 n_{e0}/(\epsilon_0 m_e)}$. Hence, we will now normalise times to $\omega_{pe}^{-1}$.

\addsubsection{Technical information: running the code}

In this practicals, you will need to run some computer simulations. The code Smilei is located in the home ($\sim$\texttt{/Smilei}). It has already been compiled and installed on the computer you're using. The documentation has been built. This means you don't need to concentrate on the code itself but on the physics.

The code input lists for these 2 exercises and this text is located in \texttt{$\sim$/TPUPMC}. Symbolic link to these files are located on the desktop as well as the code documentation (file \texttt{index.html})

To analyse the data and plot the results you will get from the simulations we provide a processing interface written in python (\texttt{$\sim$/TPUPMC/smileiQt.py}). Again, concentrate on the physics!

To run the simulation and automatically launch the post processing interface, we provide a script \texttt{tpupmc.sh}. To do so, you just need to type:

\texttt{tpupmc.sh $\sim$/TPUPMC/TP\_proj1.py}
or 
\texttt{tpupmc.sh $\sim$/TPUPMC/TP\_proj2.py}

The code will run and you will see the text output and then the post processing window will appear ad you'll be able to analyse your data. 

\begin{ExerciseList}
\vspace{0.5cm}
  \Exercise[title={Normalizations},label={ex1}] \leavevmode

    \ExeText{Considering the normalizations discussed above,}
    \Question{show that distances are now expressed in units of $c/\omega_{pe}$}
    \Question{show that momentum are in units of $m_e c$ and energies in units of $m_e c^2$}
    \Question{show that electric and magnetic fields are in units of $m_e\,c\,\omega_{pe}/e$ and $m_e\,\omega_{pe}/e$, respectively}
    \Question{show that all densities (total density, electron and ion densities) are in units of $n_e$}
\end{ExerciseList}


%%%%%%%%%%%%%%%%%%
%        Project 1            %
%%%%%%%%%%%%%%%%%%
\newpage
\section{Electron plasma oscillations}\label{proj1}

This first project aims at illustrating a central mechanism of plasma physics, the so-called electron plasma oscillations.

The code input file for this simulation are available in the file \texttt{$\sim$/TPUPMC/TP\_proj1.py} and reproduced in annex of this document.
Here, we just reproduce the sections relating to the definition of the main simulation parameters the species an their density profiles.

\addsubsection{Initial set-up}
To do so, we consider a plasma with an initial  constant density $n_{e0} = Z\,n_{i0}$ (equilibrium solution).

We also fix $n_{e0}=1$ such that times are now normalised to $\omega_{pe}^{-1}$, distances to $c/\omega_{pe}$ and so forth.

In addition, we will consider the simplest case of a one-dimensional in space (1D) infinite -~and initially cold~- plasma.
In Smilei, periodic boundary conditions are used to simulate the infinite plasma.

Those values are defined in the \class{Main} part of the namelist:
\lstinputlisting[linerange=4-14]{TP_proj1.py}

The plasma length is \code{grid\_length}  (in units of $c/\omega_{pe}$) and the time over which the simulation is run (typically a few tens of $\omega_{pe}^{-1}$) is given by \code{simulation\_time}. 

First, the ions are defined as a neutralising back ground and will not move during the simulation (i.e. their mass can be considered as infinite: \code{time\_frozen} > \code{simulation\_time}) and is defined in the \class{Species} (here for the ions):
\lstinputlisting[linerange=16-27]{TP_proj1.py}

Now, for electrons, we apply a small perturbation on their density (\code{number\_density}), and will investigate what is the plasma response to such a perturbation.
Here, we consider a periodic (cosine) perturbation of the (normalised) electron density:
\begin{eqnarray}\label{eq_pert}
n_e(t=0,x) = n_{e0}+n_{e1}(x) = 1 + \delta n\,\cos(k\,x)\,,
\end{eqnarray}
where $k$ is the mode of the parturbation and $\delta n$ its amplitude (\param{xamplitude}). 
In the simulation, what we actually fix are the simulation length $L_{\rm sim}$ (\code{grid\_length}) and the number of modes (see below) $N$ so that $k=2\pi\,N/L_{\rm sim}$.

Which translates in the code below ($N$=\param{xnumber}):
\lstinputlisting[linerange={28-38}]{TP_proj1.py}

\newpage 

\begin{ExerciseList}
\vspace{0.5cm}
  \Exercise[title={Simulations of electron-plasma oscillations},label={ex2}] \leavevmode
    \Question{\label{q2_1}In the limit of small perturbations (i.e.\ $\delta n \ll 1$ in Eq.~\ref{eq_pert}), derive analytically the equation for $n_e(t,x)$ for the case of a cold plasma with immobile ions.   What is the theoretical frequency of plasma oscillations in normalised units? Write the solution for the density perturbation $n_{e1}(x,t)$ in the form of a stationary wave with an arbitrary phase}
    \Question{\label{q2_2}Derive the relationship between $n_{e1}(t,x)$, $u_{e1,x}(t,x)$ and $E_{1,x}(t,x)$ from the linearised  fluid equation system, for a cold plasma, including Maxwell's equations.}
    \Question{In the initial set up the plasma macro-particles  can be uniformly distributed over the whole plasma length. Let us say that the interval between two near neighbours is given by  $x_{0a}-x_{0b}$, then   the equilibrium density is given by  $n_{e0} = \frac{1}{(x_{0a}-x_{0b})}$. A small perturbation $\delta x$  corresponds to a particle displacement and to a new position  $x = x_0 +  \delta x$. At time $t=0$ we initialize the system with a small perturbation and we can deduce the corresponding initial density  $n_e(x,t=0)= \frac{1}{(x_{a}-x_{b})}$. For a small perturbation we can write  $n_e(x,t=0)=n_{e0}+n_{e1}(x)$. Derive  $n_{e1}(x)$ in the limit  $\delta x << x_{0a}-x_{0b}$. If we take  $\delta x = x_1 \sin(kx +\varphi)$, what is the form of $n_{e1}(x)$?}

    \Question{\label{q2_5}Edit the code namelist  \texttt{xed $\sim$/TPUPMC/TP\_proj1.py}  and add the missing diagnostic electric field (\texttt{'Ex'}) in the \class{DiagFields}}

    \Question{To run the PIC simulation \texttt{TP\_proj1.py}. Then, in the post processing window that appears select diagnostics and choose \texttt{Uelm}, \texttt{Ex}, \texttt{Rho\_eon}, \texttt{jx}, and the phase space diagnostic (in the code this can be achieved with a \class{DiagParticleBinning}):

\lstinputlisting[linerange=45-53]{TP_proj1.py}

Where \code{axes} is in this case a two component list with 4 values:  \texttt{"name"}, minimum, maximum and number of bins. In the above case we're binning the x-px, thus the first component refers to space coordinate, and second to momentum coordinate. Typically the normalised momentum will be same order of magnitude of the normalised density fluctuations (the quantity \param{xamplitude} in \code{number\_density} of \class{Species}).

This will allow you to visualise the electrostatic field energy averaged over space as a function of time,  the electric field, the electron density and the electron current as a function of space and the momentum and position of the particles in the code (phase space) by clicking on the \texttt{plot} tab}.

    \Question{How is perturbed electron current $j_{1,x}$ related to the variable $u_{e1,x}(t,x)$?}

Notice that we are exciting a stationary wave. 

    \Question{From the simulations, deduce the explicit form of $n_{e1}(x,t)$ (i.e. by considering  $\delta n = \delta n_0 \sin (kx+\phi) \cos (\omega t)$) equation and the simulation results deduce the actual value of the arbitrary phase introduced in question (\ref{ex2}.\ref{q2_1})) and discuss the wave evolution in time. What is the characteristic time for the electric field oscillations (in normalised as well as physical units)? What about density oscillations? What is the relative phase between the electric field and density oscillations?  Are results consistent with theoretical predictions from exercise (\ref{ex2}.\ref{q2_1}) and (\ref{ex2}.\ref{q2_2})?}
    \Question{Compute analytically the time evolution of the total electrostatic field energy (variable  \texttt{Uelm} in the code):
$$U_{\rm es}(t) = \int_{0}^{L_{\rm sim}} dx E_x^2(t,x)/2$$
Plot the corresponding energy from the PIC simulation. Discuss the results: is the oscillation frequency of the energy consistent with the theoretical expectation?}
    \Question{Increase the number of modes to 2 and 3 in the input file (\param{xnumber} in \code{number\_density} of \class{Species}). Discuss the resulting changes in terms of the characteristic oscillation time.}
    \Question{Go back to mode 1. Increase the number of macro-particles per-cell used for the electron species \\(\param{particles\_per\_cell}). What does it change in your simulation? Why?}

\end{ExerciseList}

\newpage

\addsubsection{Effect of the temporal and spatial resolution}

As we are dealing with numerical simulation, the algorithms used require to be stable to fulfil some well-defined conditions:\\

\begin{itemize}
\item The FDTD method used to solve Maxwell's equations in Smilei requires the Courant-Friedrichs-Lewy (CFL) condition to be fulfilled. Basically, this condition states that the time-step $\Delta t$ should be short enough for the information on one cell not to be transported (at the light velocity) further away than one spatial-step $\Delta x$. In 1D, this condition simply requires that $c\,\Delta t < \Delta x$ (i.e. one has to take a value of \code{timestep} in our simulation (normalised) input file always smaller than \code{cell\_length}). Keep this in mind when changing the time-step in the simulations.
\item  The leap-frog scheme that is used to solve the particle equations of motion also requires the time-step to fulfill the condition: $\Delta t < 2$. If this conditions is not fulfilled, the leap-frog scheme becomes unstable.
Moreover even in the limit of small (non normalised) $\Delta t^\prime$ we introduce an error in the frequency calculation, given by 
$$\omega = \omega_0  [1+ (\omega_0 \Delta t^\prime/2 )^2/6+... ]$$
where we have indicated the theoretical value of the frequency  by $\omega_0$
\end{itemize}

In addition to these stability conditions, the time-step and cell-length have to be small enough to describe the physics of interest. For instance, to a given time-step $\Delta t$ correspond a maximum frequency, the so-called Nyquist (critical) frequency $\nu_c =(2\,\Delta t)^{-1}$, beyond which oscillations can not be resolved. Similarly, for a given cell-length $\Delta x$, the modes $k>(\pi/\Delta x)$ can not be resolved by the code.  The phenomenon of aliasing appears in this case, that is the code couples high frequencies(wave vectors)  with low frequencies(wave vectors), i.e. interprets high frequencies (wave vectors) as low frequencies (wave vectors), and fills the low frequency(wave vectors) space. 

The next two exercises (\ref{ex3} and \ref{ex4}) aim at investigating the problem of time (then spatial) resolution. 

As discussed above, the choice of the time interval can affect the precision and the performance of the code. We will thus take different time intervals. 


\begin{ExerciseList}
\vspace{0.5cm}
  \Exercise[title={Effect of the temporal resolution},label={ex3}] \leavevmode
   \Question{What are the units of $\Delta t$?}

   \Question{How many points describe a period of oscillation for a temporal resolution of $\Delta t = 0.2$}
    \Question{Verify that you are back to 10 particles per cell. Keeping $\Delta t = 0.2$, calculate the correction to the frequency from the formula above.  
    
    Take in the code $\code{timestep} = 0.2$, $\code{cell\_length} = 0.24$, $\code{grid\_length}= 4.8$, and in the \\\class{DiagParticleBinning} for the second component of \code{axes} \texttt{["px", -0.12, 0.12,100]} (extremes values for the momentum). For the duration of the simulation, it is useful not to have more than 20 oscillations. Deduce the oscillation frequency form the time evolution of the electrostatic energy in the code and compare with the theoretical result. Comment the result.}
    \Question{Calculate the value of $\Delta t $ such that the Nyquist frequency is equal to the frequency of oscillation of the electrostatic energy. N.B. $\omega_c = 2 \pi \nu_c$. Why do we calculate this value?}
    \Question{Measure in the code  the time of oscillation of the electrostatic energy for different value of $\Delta t$ (\code{timestep}) and modify \code{cell\_length} accordingly: you need to take care to fulfil the CFL condition in your choice of parameters, e.g. by using $\Delta x = 1.2\,\Delta t$ (note that you also need to modify $\code{grid\_length}=20~\Delta x$). Deduce the oscillation time of the plasma wave, the oscillation frequency of the electrostatic energy, and the oscillation frequency of the of the plasma wave, you can write your results in the table. The maximum value of $\Delta t $ that we will take is $\Delta t = 2.1$ and $\code{simulation\_time}=100*2*\texttt{math.pi}$. What happens for this value? Take a longer duration of the simulation if needed.}
    \Question{Make a plot of the oscillation frequency of the plasma wave as a function of $\Delta t$, and discuss the result. }
\end{ExerciseList}


\vspace{10pt}
\begin{tabularx}{0.95\textwidth}{X|X|X|X|X}
\centering\textbf{$\Delta t$} & \centering\textbf{$\tau_{E,n}$} & \centering\textbf{$\tau_{P,n}$} &  \centering\textbf{$\omega_{E,n}$} & \centering\arraybackslash\textbf{$\omega_{P,n}$} \\
\hline
 & & & & \\  & & & & \\ \hline  & & & & \\  & & & & \\ \hline  & & & & \\  & & & & \\ \hline  & & & & \\  & & & & \\ \hline  & & & & \\  & & & & \\ \hline  & & & & \\  & & & & \\ \hline  & & & & \\  & & & & \\ 
\end{tabularx}

\begin{ExerciseList}
\vspace{0.5cm}
  \Exercise[title={Effect of the spatial resolution and aliasing effect},label={ex4}] \leavevmode

    \ExeText{In this exercise we will fix $\Delta t = 0.01$, the plasma size to $1.$, the cell size to $0.02$, 
 in the \class{DiagParticleBinning} take the second value of  the first component of \code{axes} = 1. (maximum value for the position/plasma length), and for the the second component of \code{axes} : \texttt{["px", -0.02, 0.02, 100]} (extremes values for the momentum). 
Vary the wave vector $k$ by varying the quantity \param{xnumber} in \code{number\_density}, since the  wave-number $k$ is defined as :
\begin{equation*}
k = {\rm \param{xnumber}} \frac{2\pi}{\rm \code{grid\_length}}\,.
\end{equation*}
}

    \Question{Theoretically does the cold plasma oscillations frequency depend on the wave vector?}
    \Question{If we modify the wave vector without modifying the plasma size and the cell size (\code{cell\_length}), are we still resolving correctly one wavelength? (You can look at the code results to help you answer this question).}
    \Question{What is the limiting value of the wave vector for the given cell size? This is the equivalent of the Nyquist frequency, the theoretical limit of sampling.}
    \Question{Vary the wave vector $k$ up to the theoretical limiting value  and further. For each $k$ value considered measure the period of the electrostatic fields energy oscillation, deduce the period of the plasma oscillations and the frequency for the energy and the fields. 
You can use the table to note your results. Plot the oscillation frequency as a function of the wave vector. Comment the plot.}
\end{ExerciseList}

\vspace{10pt}
\begin{tabularx}{0.95\textwidth}{X|X|X|X|X}
\centering\large\textbf{$k$} & \centering\large\textbf{$\tau_{E,n}$} & \centering\large\textbf{$\tau_{P,n}$} &  \centering\large\textbf{$\omega_{E,n}$} & \centering\arraybackslash\large\textbf{$\omega_{P,n}$} \\
\hline
 & & & & \\  & & & & \\ \hline  & & & & \\  & & & & \\ \hline  & & & & \\  & & & & \\ \hline  & & & & \\  & & & & \\ \hline  & & & & \\  & & & & \\ \hline  & & & & \\  & & & & \\ \hline  & & & & \\  & & & & \\  
\end{tabularx}

\newpage 
\addsubsection{Kinetic (wave-breaking) effects}

We will now investigate the changes brought by taking a large amplitude for the initial density perturbation ($\delta~n~\lesssim~1$). To do so, it is necessary to take sufficiently small time-step and cell-length (take care to fulfil the CFL condition in your choice of parameters, e.g. by using $\Delta x = 2\,\Delta t $).

\begin{ExerciseList}
\vspace{0.5cm}
  \Exercise[title={Kinetic effects for large amplitude plasma oscillations},label={ex5}] \leavevmode

    \Question{In the conditions of exercise \ref{ex2}, (e.g. \code{grid\_length}=1., \code{cell\_length}=0.02, \code{timestep}=0.01 and $\delta n = 0.01$) plot the electron phase-space ($x,p_x$). Does this bring any additional information as compared to the fluid description? Justify your answer.}
    \Question{Now, increase the amplitude of the density perturbation to $\delta n = 0.5$ and reduce the time step to $\Delta t =  10^{-3}$ and the cell size accordingly. Check  the electron phase-space ($x,p_x$). (Note that you may have to increase the minimum and maximum value of the momentum when defining the phase-space diagnostics in the input file, parameters of the maximum and minimum momentum of the binning, see Exercise \ref{ex2}.\ref{q2_5}). Discuss what you observed: 
    \begin{itemize}
    \item does the kinetic description brings more information that the linearised fluid approach? 
    \item What about the temporal evolution of the electrostatic field energy?
    \end{itemize}}
    
    \Question{Now increase $\delta n$ to 0.9. Check the  electron phase-space ($x,p_x$) and the temporal evolution of the electrostatic field energy and of the particles kinetic energy. 
    \begin{itemize}
    \item How does it differ from what was observed in the limit $\delta n \ll 1$ (exercise \ref{ex2})? 
    \item Can you explain such differences? 
    \item What about the phase-space? 
    \item Does this bring any additional information as compared to the fluid description? Justify your answer.
    \end{itemize}}
\end{ExerciseList}

%%%%%%%%%%%%%%%%%%
%        Project 2            %
%%%%%%%%%%%%%%%%%%
\newpage
\section{Two-stream instability}\label{proj2}

\addsubsection{Initial set-up and theory}

This second project aims at illustrating the so-called two-stream instability.
To do so, we consider a 1D infinite plasma with averaged constant density $n_{e0} = 1$. 
In Smilei, this is simulated by considering periodic boundary conditions in the x-direction. 
To create the two streams, we consider two electron species with initially opposite mean-velocities $v_{\pm} = \pm v_0$ (with $v_0 = 0.1$ in units of the light velocity) in the x-direction.
Each species (labelled $+$ and $-$ depending on the direction of their mean-velocity) is then initialized such that:
\begin{eqnarray*}
n_e(t=0,x) = n_e^+(x) + n_e^-(x)\,,
\end{eqnarray*}
where:
\begin{eqnarray*}
n_e^+(t=0,x) = n_e^-(t=0,x) = \frac{n_{e0}}{2}\,\left[1 + \delta n\,\cos(k\,x)\right]\,.
\end{eqnarray*}

\begin{wrapfigure}[15]{r}{0.4\textwidth}
\def\svgwidth{\linewidth}
\input{Fig_2stream.pdf_tex}
\caption{Two stream instability}
\label{Fig_2stream}
\end{wrapfigure}

Here, $\delta n = 10^{-2}$ is the amplitude of the perturbation and $k$ its mode.
As periodic conditions allow us for simulating an infinite plasma, and taking $n=1$ mode  (in the input file, see Proj.~\ref{proj1}), the wave-number $k$ is directly linked to the length of the simulation box by:
\begin{eqnarray}\label{eq_k}
k = \frac{2\pi\,n}{\code{grid\_length}}\,.
\end{eqnarray}


Now, we recall the dispersion relation for such a system. In normalised units, it reads:
\begin{eqnarray}\label{eq_disprel_2stream}
\frac{1}{(\omega - k\,v_0)^2} + \frac{1}{(\omega + k\,v_0)^2} = 2\,,
\end{eqnarray}
where $k$ is in unit of $\omega_{p0}/c$ and $\omega$ is in units of $\omega_{p0}$.


Note that here, $\omega_{p0}$ denotes the plasma frequency at the total (background) electron density $n_0$, while the beam (half-density) $n_0/2 = \langle n_e^{\pm}(t=0) \rangle$ was used in the lecture.


\begin{ExerciseList} 
\vspace{0.5cm}
  \Exercise[title={Dispersion relation, characteristic frequencies and growth rate},label={ex6}] \leavevmode

    \Question{For a given value of $v_0 k$, show that Eq.~\eqref{eq_disprel_2stream} is equivalent to the biquadratic equation on $\omega$:
\begin{eqnarray}\label{eq_biquad}
\omega^4 - \left[1+2 (v_0 k)^2\right]\,\omega^2 + (v_0 k)^2\,\left[(v_0 k)^2-1\right] = 0\,\,.
\end{eqnarray}}
    \Question{Solve Eq.~\eqref{eq_biquad} to obtain:
\begin{eqnarray}\label{eqomega}
\omega^2 = \frac{1}{2}\,\left[ \left(1+2 (v_0 k)^2\right) \pm \sqrt{1+8 (v_0 k)^2} \right]
\end{eqnarray}
Show that it admits imaginary solution for $v_0 k<1$. To what does this imaginary solution correspond to?}
    \Question{Plot, as a function of $v_0 k$ (e.g. in the range $0 - 5$),  the solution for $\omega$ (use dashed lines for the imaginary solution). [For the current exercise, you can directly use the plot given in Figure~\ref{Fig_2stream}]}
    \Question{Considering that the dispersion relation is obtained assuming all perturbed quantities vary as $\mathrm{e}^{i(\omega t - k\,x)}$, show that, for $v_0 k>1$, the energy in the electrostatic field varies as 
\begin{equation*}
A + B\,\cos(2\,\omega_a\,t) + C\,\cos(2\,\omega_b\,t) + D\,\cos\left[(\omega_a-\omega_b)\,t\right] + E\,\cos\left[(\omega_a+\omega_b)\,t\right]
\end{equation*}
 with $A$, $B$, $C$, $D$ and $E$ some numerical constants. What do $\omega_a$ and $\omega_b$ stand for?}
    \Question{Similarly, show that, for $v_0 k<1$, the energy in the electrostatic field goes as  
\begin{equation*}
A + B\,\cos(2\omega t) + C\,\cos(\omega t)\,\mathrm{e}^{\gamma t}+D\,\mathrm{e}^{2\gamma t}
\end{equation*}

What do $\omega$ and $\gamma$ stand for? }
\end{ExerciseList}



\addsubsection{Numerical simulation of the 2-stream instability}

We will now perform numerical simulation using the PIC code Smilei for different values of $v_0 k$.
As discussed above, we take $v_0=10^{-1}$ which ensures that relativistic effects can be neglected.
In addition, we use periodic boundary conditions, so that $k$ is uniquely defined by the simulation length using Eq.~\eqref{eq_k}.
The input file to start these simulation is given in the file \texttt{TPUPMC/TP\_proj2.py}.

If not, set the density perturbation $\delta n = 0.001$.

\begin{ExerciseList}
\vspace{0.5cm}
  \Exercise[title={Simulation for $v_0 k > 1$},label={ex7}] \leavevmode

    \Question{Compute the value of \code{grid\_length} corresponding to the limit between the stable and unstable solutions.}
     
    \Question{Compute the values of $v_0 k$ corresponding to different values of $k$ using  \code{grid\_length} $= 0.62, 0.60, 0.31, 0.16$. Calculate the corresponding values of $\omega$ from Eq. \ref{eqomega}.}
    \Question{Run the corresponding simulations.  For the phase-space diagnostics, in the input file \texttt{TP\_proj2.py}. For simulations with \code{grid\_length} $= 0.62$ and $0.60$, extract from the field energy temporal evolution the characteristic frequencies (try to estimate error bars for your measurements). Report them on the analytical predictions (from Exercise \ref{ex6}). Why does it become more difficult to extract the characteristic frequencies for \code{grid\_length}~$= 0.31 - 0.16$ ?}
    \Question{ Discuss the general behaviour of 
        \begin{enumerate}
    \item the energy in the electrostatic field,
    \item the behaviour in phase-space of the electron species
    \end{enumerate} }
\end{ExerciseList}


\begin{ExerciseList}
\vspace{0.5cm}
  \Exercise[title={Simulation for {$v_0 k < 1$}},label={ex8}] \leavevmode

    \Question{Compute the values of $v_0 k$ corresponding to \code{grid\_length} $= 0.69 - 1.03 - 1.68$. Calculate the corresponding values of $\omega$ and $\gamma$ from Eq. \ref{eqomega}.}
    \Question{Run the corresponding simulations. For the phase-space diagnostics, in the input file \texttt{TP\_proj2.py},  take for the extremes values of the momentum in \class{DiagParticleBinning} section (i.e.\ the  \texttt{"px"}  part of the \code{axes})   to $\pm 0.4$, respectively. For each simulation (i.e. for different value of the parameter $v_0 k$), extract from the field energy temporal evolution the characteristic frequency and growth rate of the instability (it is convenient to consider the energy plot in log-scale. To do so, just type on the corresponding panel the letter L). Report them on the analytical predictions (from Exercise \ref{ex6}).}
    \Question{For each simulation, discuss the general behaviour of 
    \begin{enumerate}
    \item the energy in the electrostatic field,
    \item the behaviour in phase-space of the electron species
    \end{enumerate} }
    \Question{For these simulations, using the temporal evolution of the electrostatic energy, discuss over what characteristic time the instability is in its linear phase.}
    \Question{Now, considering the phase-space distribution of the electron species, discuss what happens in the nonlinear phase of the instability.}
\end{ExerciseList}

\vspace{10pt}
\hrule

\addcontentsline{toc}{section}{References}
\begin{thebibliography}{2}
\bibitem{birdsall_langdon} C.~K.~Birdsall and A.~B.~Langdon, {\it Plasma physics via computer simulation}, McGraw-Hill, New York (1985).
\bibitem{taflove_2005} A.~Taflove and S.~C.~Hagness, {\it Computational Electrodynamics: The Finite-Difference Time-Domain Method}, 3rd~ed. Norwood, MA: Artech House (2005).
\bibitem{boris_1970} J.~P.~Boris, {\it Relativistic plasma simulations - Optimization of a hybrid code}, Proc. 4th Conf. Num. Sim. of Plasmas {\bf 3}, 1970.
\end{thebibliography}

%%%%%%%%%%%%%%%%
%  REFERENCES  %
%%%%%%%%%%%%%%%%
\newpage
\addsection{Smilei namelists}

\addsubsection{\small\texttt{TP\_proj1.py}}
\lstinputlisting{TP_proj1.py}

\newpage
\addsubsection{\small\texttt{TP\_proj2.py}}
\lstinputlisting{TP_proj2.py}






\end{document}