path: root/Docs/source/latex_theory/Boosted_frame/Boosted_frame.tex

% Copyright 2017-2019 Burlen Loring, Remi Lehe
%
% This file is part of WarpX.
%
% License: BSD-3-Clause-LBNL

\input{newcommands}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Moving window and optimal Lorentz boosted frame}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The simulations of plasma accelerators from first principles are extremely computationally intensive, due to the need to resolve the evolution of a driver (laser or particle beam) and an accelerated particle beam into a plasma structure that is orders of magnitude longer and wider than the accelerated beam. As is customary in the modeling of particle beam dynamics in standard particle accelerators, a moving window is commonly used to follow the driver, the wake and the accelerated beam. This results in huge savings, by avoiding the meshing of the entire plasma that is orders of magnitude longer than the other length scales of interest.

\begin{figure}
%\begin{centering}
\includegraphics[scale=0.6]{Boosted_frame.png}
%\par\end{centering}
\caption{\label{fig:PIC} A first principle simulation of a short driver beam (laser or charged particles) propagating through a plasma that is orders of magnitude longer necessitates a very large number of time steps. Recasting the simulation in a frame of reference that is moving close to the speed of light in the direction of the driver beam leads to simulating a driver beam that appears longer propagating through a plasma that appears shorter than in the laboratory. Thus, this relativistic transformation of space and time reduces the disparity of scales, and thereby the number of time steps to complete the simulation, by orders of magnitude.}
\end{figure}

Even using a moving window, however, a full PIC simulation of a plasma accelerator can be extraordinarily demanding computationally, as many time steps are needed to resolve the crossing of the short driver beam with the plasma column. As it turns out, choosing an optimal frame of reference that travels close to the speed of light in the direction of the laser or particle beam (as opposed to the usual choice of the laboratory frame) enables speedups by orders of magnitude \cite{Vayprl07,Vaypop2011}. This is a result of the properties of Lorentz contraction and dilation of space and time. In the frame of the laboratory, a very short driver (laser or particle) beam propagates through a much longer plasma column, necessitating millions to tens of millions of time steps for parameters in the range of the BELLA or FACET-II experiments. As sketched in Fig. \ref{fig:PIC}, in a frame moving with the driver beam in the plasma at velocity $v=\beta c$ (where $c$ is the speed of light in vacuum), the beam length is now elongated by $\approx(1+\beta)\gamma$ while the plasma contracts by $\gamma$ (where $\gamma=1/\sqrt{1-\beta^2}$ is the relativistic factor associated with the frame velocity). The number of time steps that is needed to simulate a ``longer'' beam through a ``shorter'' plasma is now reduced by up to $\approx(1+\beta) \gamma^2$ (a detailed derivation of the speedup is given below).

The modeling of a plasma acceleration stage in a boosted frame
involves the fully electromagnetic modeling of a plasma propagating at near the speed of light, for which Numerical Cerenkov
\cite{Borisjcp73,Habericnsp73} is a potential issue, as explained in more details below.
In addition, for a frame of reference moving in the direction of the accelerated beam (or equivalently the wake of the laser),
waves emitted by the plasma in the forward direction expand
while the ones emitted in the backward direction contract, following the properties of the Lorentz transformation.
If one had to resolve both forward and backward propagating
waves emitted from the plasma, there would be no gain in selecting a frame different from the laboratory frame. However,
the physics of interest for a laser wakefield is the laser driving the wake, the wake, and the accelerated beam.
Backscatter is weak in the short-pulse regime, and does not
interact as strongly with the beam as do the forward propagating waves
which stay in phase for a long period. It is thus often assumed that the backward propagating waves
can be neglected in the modeling of plasma accelerator stages. The accuracy  of this assumption has been demonstrated by
comparison between explicit codes which include both forward and backward waves and envelope or quasistatic codes which neglect backward waves
\cite{Geddesjp08,Geddespac09,Cowanaac08}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Theoretical speedup dependency with the frame boost}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The derivation that is given here reproduces the one given in \cite{Vaypop2011}, where the obtainable speedup is derived as an extension of the formula that was derived earlier\cite{Vayprl07}, taking in addition into account the group velocity of the laser as it traverses the plasma.

Assuming that the simulation box is a fixed number of plasma periods long, which implies the use (which is standard) of a moving window following
the wake and accelerated beam, the speedup is given by the ratio of the time taken by the laser pulse and the plasma to cross each other, divided by the shortest time scale of interest, that is the laser period. To first order, the wake velocity $v_w$ is set by the 1D group velocity of the laser driver, which in the linear (low intensity) limit, is given by \cite{Esareyrmp09}:

%
\begin{equation}
v_w/c=\beta_w=\left(1-\frac{\omega_p^2}{\omega^2}\right)^{1/2}
\end{equation}
%
where $\omega_p=\sqrt{(n_e e^2)/(\epsilon_0 m_e)}$ is the plasma frequency, $\omega=2\pi c/\lambda$ is the laser frequency, $n_e$ is the plasma density, $\lambda$ is the laser wavelength in vacuum, $\epsilon_0$ is the permittivity of vacuum, $c$ is the speed of light in vacuum, and $e$ and $m_e$ are respectively the charge and mass of the electron.

In practice, the runs are typically stopped when the last electron beam macro-particle exits the plasma, and a measure of the total time of the simulation is then given by
%
\begin{equation}
T=\frac{L+\eta \lambda_p}{v_w-v_p}
\end{equation}
%
where $\lambda_p\approx 2\pi c/\omega_p$ is the wake wavelength, $L$ is the plasma length, $v_w$ and $v_p=\beta_p c$ are respectively the velocity of the wake and of the plasma relative to the frame of reference, and $\eta$ is an adjustable parameter for taking into account the fraction of the wake which exited the plasma at the end of the simulation.
For a beam injected into the $n^{th}$ bucket, $\eta$ would be set to $n-1/2$. If positrons were considered, they would be injected half a wake period ahead of the location of the electrons injection position for a given period, and one would have $\eta=n-1$. The numerical cost $R_t$ scales as the ratio of the total time to the shortest timescale of interest, which is the inverse of the laser frequency, and is thus given by
%
\begin{equation}
R_t=\frac{T c}{\lambda}=\frac{\left(L+\eta \lambda_p\right)}{\left(\beta_w-\beta_p\right) \lambda}
\end{equation}
%
In the laboratory, $v_p=0$ and the expression simplifies to
%
\begin{equation}
R_{lab}=\frac{T c}{\lambda}=\frac{\left(L+\eta \lambda_p\right)}{\beta_w \lambda}
\end{equation}
%
In a frame moving at $\beta c$, the quantities become
\begin{eqnarray}
\lambda_p^*&=&\lambda_p/\left[\gamma \left(1-\beta_w \beta\right)\right] \\
L^*&=&L/\gamma \\
\lambda^*&=& \gamma\left(1+\beta\right) \lambda\\
\beta_w^*&=&\left(\beta_w-\beta\right)/\left(1-\beta_w\beta\right) \\
v_p^*&=&-\beta c \\
T^*&=&\frac{L^*+\eta \lambda_p^*}{v_w^*-v_p^*} \\
R_t^*&=&\frac{T^* c}{\lambda^*} = \frac{\left(L^*+\eta \lambda_p^*\right)}{\left(\beta_w^*+\beta\right) \lambda^*}
\end{eqnarray}
where $\gamma=1/\sqrt{1-\beta^2}$.

The expected speedup from performing the simulation in a boosted frame is given by the ratio of $R_{lab}$ and $R_t^*$

\begin{equation}
S=\frac{R_{lab}}{R_t^*}=\frac{\left(1+\beta\right)\left(L+\eta \lambda_p\right)}{\left(1-\beta\beta_w\right)L+\eta \lambda_p}
\label{Eq_scaling1d0}
\end{equation}

We note that assuming that $\beta_w\approx1$ (which is a valid approximation for most practical cases of interest) and that $\gamma<<\gamma_w$, this expression is consistent with the expression derived earlier \cite{Vayprl07} for the laser-plasma acceleration case, which states that $R_t^*=\alpha R_t/\left(1+\beta\right)$ with $\alpha=\left(1-\beta+l/L\right)/\left(1+l/L\right)$, where $l$ is the laser length which is generally proportional to $\eta \lambda_p$, and $S=R_t/R_T^*$. However, higher values of $\gamma$ are of interest for maximum speedup, as shown below.

For intense lasers ($a\sim 1$) typically used for acceleration, the energy gain is limited by dephasing \cite{Schroederprl2011}, which occurs over a scale length $L_d \sim \lambda_p^3/2\lambda^2$.
Acceleration is compromised beyond $L_d$ and in practice, the plasma length is proportional to the dephasing length, i.e. $L= \xi L_d$. In most cases, $\gamma_w^2>>1$, which allows the approximations $\beta_w\approx1-\lambda^2/2\lambda_p^2$, and $L=\xi \lambda_p^3/2\lambda^2\approx \xi \gamma_w^2 \lambda_p/2>>\eta \lambda_p$, so that Eq.(\ref{Eq_scaling1d0}) becomes
%
\begin{equation}
S=\left(1+\beta\right)^2\gamma^2\frac{\xi\gamma_w^2}{\xi\gamma_w^2+\left(1+\beta\right)\gamma^2\left(\xi\beta/2+2\eta\right)}
\label{Eq_scaling1d}
\end{equation}
%
For low values of $\gamma$, i.e. when $\gamma<<\gamma_w$, Eq.(\ref{Eq_scaling1d}) reduces to
%
\begin{equation}
S_{\gamma<<\gamma_w}=\left(1+\beta\right)^2\gamma^2
\label{Eq_scaling1d_simpl2}
\end{equation}
%
Conversely, if $\gamma\rightarrow\infty$, Eq.(\ref{Eq_scaling1d}) becomes
%
\begin{equation}
S_{\gamma\rightarrow\infty}=\frac{4}{1+4\eta/\xi}\gamma_w^2
\label{Eq_scaling_gamma_inf}
\end{equation}
%
Finally, in the frame of the wake, i.e. when $\gamma=\gamma_w$, assuming that $\beta_w\approx1$, Eq.(\ref{Eq_scaling1d}) gives
%
\begin{equation}
S_{\gamma=\gamma_w}\approx\frac{2}{1+2\eta/\xi}\gamma_w^2
\label{Eq_scaling_gamma_wake}
\end{equation}
Since $\eta$ and $\xi$ are of order unity, and the practical regimes of most interest satisfy $\gamma_w^2>>1$, the speedup that is obtained by using the frame of the wake will be near the maximum obtainable value given by Eq.(\ref{Eq_scaling_gamma_inf}).

Note that without the use of a moving window, the relativistic effects that are at play in the time domain would also be at play in the spatial domain \cite{Vayprl07}, and the $\gamma^2$ scaling would transform to $\gamma^4$. Hence, it is important to use a moving window even in simulations in a Lorentz boosted frame. For very high values of the boosted frame, the optimal velocity of the moving window may vanish (i.e. no moving window) or even reverse.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Numerical Stability and alternate formulation in a Galilean frame}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The numerical Cherenkov instability (NCI) \cite{Godfreyjcp74}
is the most serious numerical instability affecting multidimensional
PIC simulations of relativistic particle beams and streaming plasmas
\cite{Martinscpc10,VayAAC2010,Vayjcp2011,Spitkovsky:Icnsp2011,GodfreyJCP2013,XuJCP2013}.
It arises from coupling between possibly numerically distorted electromagnetic modes and spurious
beam modes, the latter due to the mismatch between the Lagrangian
treatment of particles and the Eulerian treatment of fields \cite{Godfreyjcp75}.

In recent papers the electromagnetic dispersion
relations for the numerical Cherenkov instability were derived and solved for both FDTD \cite{GodfreyJCP2013,GodfreyJCP2014_FDTD}
and PSATD \cite{GodfreyJCP2014_PSATD,GodfreyIEEE2014} algorithms.

Several solutions have been proposed to mitigate the NCI \cite{GodfreyJCP2014,GodfreyIEEE2014,GodfreyJCP2014_PSATD,GodfreyCPC2015,YuCPC2015,YuCPC2015-Circ}. Although
these solutions efficiently reduce the numerical instability,
they typically introduce either strong smoothing of the currents and
fields, or arbitrary numerical corrections, which are
tuned specifically against the NCI and go beyond the
natural discretization of the underlying physical equation. Therefore,
it is sometimes unclear to what extent these added corrections could impact the
physics at stake for a given resolution.

For instance, NCI-specific corrections include periodically smoothing
the electromagnetic field components \cite{Martinscpc10},
using a special time step \cite{VayAAC2010,Vayjcp2011} or
applying a wide-band smoothing of the current components \cite{VayAAC2010,Vayjcp2011,VayPOPL2011}. Another set of mitigation methods
involve scaling the deposited
currents by a carefully-designed wavenumber-dependent factor
\cite{GodfreyJCP2014_FDTD,GodfreyIEEE2014} or slightly modifying the
ratio of electric and magnetic fields ($E/B$) before gathering their
value onto the macroparticles
\cite{GodfreyJCP2014_PSATD,GodfreyCPC2015}.
Yet another set of NCI-specific corrections
\cite{YuCPC2015,YuCPC2015-Circ} consists
in combining a small timestep $\Delta t$, a sharp low-pass spatial filter,
and a spectral or high-order scheme that is tuned so as to
create a small, artificial ``bump'' in the dispersion relation
\cite{YuCPC2015}. While most mitigation methods have only been applied
to Cartesian geometry, this last
set of methods (\cite{YuCPC2015,YuCPC2015-Circ})
has the remarkable property that it can be applied
\cite{YuCPC2015-Circ} to both Cartesian geometry and
quasi-cylindrical geometry (i.e. cylindrical geometry with
azimuthal Fourier decomposition \cite{LifschitzJCP2009,DavidsonJCP2015,Lehe2016}). However,
the use of a small timestep proportionally slows down the progress of
the simulation, and the artificial ``bump'' is again an arbitrary correction
that departs from the underlying physics.

A new scheme was recently proposed, in \cite{KirchenARXIV2016,LeheARXIV2016}, which
completely eliminates the NCI for a plasma drifting at a uniform relativistic velocity
-- with no arbitrary correction -- by simply integrating
the PIC equations in \emph{Galilean coordinates} (also known as
\emph{comoving coordinates}). More precisely, in the new
method, the Maxwell equations \emph{in Galilean coordinates} are integrated
analytically, using only natural hypotheses, within the PSATD
framework (Pseudo-Spectral-Analytical-Time-Domain \cite{Habericnsp73,VayJCP2013}).

The idea of the proposed scheme is to perform a Galilean change of
coordinates, and to carry out the simulation in the new coordinates:
\begin{equation}
\label{eq:change-var}
\vec{x}' = \vec{x} - \vgal t
\end{equation}
where $\vec{x} = x\,\vec{u}_x + y\,\vec{u}_y + z\,\vec{u}_z$ and
$\vec{x}' = x'\,\vec{u}_x + y'\,\vec{u}_y + z'\,\vec{u}_z$ are the
position vectors in the standard and Galilean coordinates
respectively.

When choosing $\vgal= \vec{v}_0$, where
$\vec{v}_0$ is the speed of the bulk of the relativistic
plasma, the plasma does not move with respect to the grid in the Galilean
coordinates $\vec{x}'$ -- or, equivalently, in the standard
coordinates $\vec{x}$, the grid moves along with the plasma. The heuristic intuition behind this scheme
is that these coordinates should prevent the discrepancy between the Lagrangian and
Eulerian point of view, which gives rise to the NCI \cite{Godfreyjcp75}.

An important remark is that the Galilean change of
coordinates (\ref{eq:change-var}) is a simple translation. Thus, when used in
the context of Lorentz-boosted simulations, it does
of course preserve the relativistic dilatation of space and time which gives rise to the
characteristic computational speedup of the boosted-frame technique.

Another important remark is that the Galilean scheme is \emph{not}
equivalent to a moving window (and in fact the Galilean scheme can be
independently \emph{combined} with a moving window). Whereas in a
moving window, gridpoints are added and removed so as to effectively
translate the boundaries, in the Galilean scheme the gridpoints
\emph{themselves} are not only translated but in this case, the physical equations
are modified accordingly. Most importantly, the assumed time evolution of
the current $\vec{J}$ within one timestep is different in a standard PSATD scheme with moving
window and in a Galilean PSATD scheme \cite{LeheARXIV2016}.

In the Galilean coordinates $\vec{x}'$, the equations of particle
motion and the Maxwell equations take the form
\begin{subequations}
\begin{align}
\frac{d\vec{x}'}{dt} &= \frac{\vec{p}}{\gamma m} - \vgal \label{eq:motion1} \\
\frac{d\vec{p}}{dt} &= q \left( \vec{E} +
\frac{\vec{p}}{\gamma m} \times \vec{B} \right) \label{eq:motion2}\\
\left( \Dt{\;} - \vgal\cdot\nab\right)\vec{B} &= -\nab\times\vec{E} \label{eq:maxwell1}\\
\frac{1}{c^2}\left( \Dt{\;} - \vgal\cdot\nab\right)\vec{E} &= \nab\times\vec{B} - \mu_0\vec{J} \label{eq:maxwell2}
\end{align}
\end{subequations}
where $\nab$ denotes a spatial derivative with respect to the
Galilean coordinates $\vec{x}'$.

Integrating these equations from $t=n\Delta
t$ to $t=(n+1)\Delta t$ results in the following update equations (see
\cite{LeheARXIV2016} for the details of the derivation):
%
\begin{subequations}
\begin{align}
\fb^{n+1} &= \theta^2 C \fb^n
 -\frac{\theta^2 S}{ck}i\vec{k}\times \fe^n \nonumber \\
& + \;\frac{\theta \chi_1}{\epsilon_0c^2k^2}\;i\vec{k} \times
                     \fj^{n+1/2} \label{eq:disc-maxwell1}\\
\fe^{n+1} &=  \theta^2 C  \fe^n
 +\frac{\theta^2 S}{k} \,c i\vec{k}\times \fb^n \nonumber \\
& +\frac{i\nu \theta \chi_1 - \theta^2S}{\epsilon_0 ck} \; \fj^{n+1/2}\nonumber \\
& - \frac{1}{\epsilon_0k^2}\left(\; \chi_2\;\mc{\rho}^{n+1} -
  \theta^2\chi_3\;\mc{\rho}^{n} \;\right) i\vec{k} \label{eq:disc-maxwell2}
\end{align}
\end{subequations}
%
where we used the short-hand notations $\fe^n \equiv
%
\fe(\vec{k}, n\Delta t)$, $\fb^n \equiv
\fb(\vec{k}, n\Delta t)$ as well as:
\begin{subequations}
\begin{align}
&C = \cos(ck\Delta t) \quad S = \sin(ck\Delta t) \quad k
= |\vec{k}| \label{eq:def-C-S}\\&
\nu = \frac{\vec{k}\cdot\vgal}{ck} \quad \theta =
  e^{i\vec{k}\cdot\vgal\Delta t/2} \quad \theta^* =
  e^{-i\vec{k}\cdot\vgal\Delta t/2} \label{eq:def-nu-theta}\\&
\chi_1 =  \frac{1}{1 -\nu^2} \left( \theta^* -  C \theta + i
  \nu \theta S \right) \label{eq:def-chi1}\\&
\chi_2 = \frac{\chi_1 - \theta(1-C)}{\theta^*-\theta} \quad
\chi_3 = \frac{\chi_1-\theta^*(1-C)}{\theta^*-\theta} \label{eq:def-chi23}
\end{align}
\end{subequations}
Note that, in the limit $\vgal=\vec{0}$,
(\ref{eq:disc-maxwell1}) and (\ref{eq:disc-maxwell2}) reduce to the standard PSATD
equations \cite{Habericnsp73}, as expected.
As shown in \cite{KirchenARXIV2016,LeheARXIV2016},
the elimination of the NCI with the new Galilean integration is verified empirically via PIC simulations of uniform drifting plasmas and laser-driven plasma acceleration stages, and confirmed by a theoretical analysis of the instability.