# Efficient start-to-end 3D envelope modeling for two-stage laser wakefield
acceleration experiments

F. Massimo1 A. Beck1 J. Derouillat2 M. Grech3 M. Lobet2 F. Pérez3 I. Zemzemi1
and A. Specka1

(1 Laboratoire Leprince-Ringuet – École polytechnique, CNRS-IN2P3, Palaiseau
91128, France
2 Maison de la Simulation, CEA, CNRS, Université Paris-Sud, UVSQ, Université
Paris-Saclay, F-91191 Gif-sur-Yvette, France
3 Laboratoire d’Utilisation des Lasers Intenses, CNRS, École polytechnique,
CEA, Université Paris-Saclay, UPMC Université Paris 06: Sorbonne Universités,
F-91128 Palaiseau Cedex, France )

###### Abstract

Three dimensional Particle in Cell simulations of Laser Wakefield Acceleration
require a considerable amount of resources but are necessary to have realistic
predictions and to design future experiments. The planned experiments for the
Apollon laser also include two stages of plasma acceleration, for a total
plasma length of the order of tens of millimeters or centimeters. In this
context, where traditional 3D numerical simulations would be unfeasible, we
present the results of the application of a recently proposed envelope method,
to describe the laser pulse ant its interaction with the plasma without the
need to resolve its high frequency oscillations. The implementation of this
model in the code Smilei is described, as well as the results of benchmark
simulations against standard laser simulations and applications for the design
of two stage Apollon experiments.

††*Corresponding author. E-mail address<EMAIL_ADDRESS>
##  1 Introduction

The maximum electric field sustainable by an accelerating cavity determines
the minimum size of a particle accelerator. The breakdown limits of the
metallic accelerating cavities in conventional accelerators motivated the
accelerator community to find alternative technologies to achieve higher
accelerating gradients and thus smaller particle accelerators. One of the most
promising of these technologies to accelerate electrons is the Laser Wakefield
Acceleration (LWFA), i.e. their acceleration by plasma waves generated in the
wake of an intense laser pulse propagating in an under dense plasma [1, 2, 3,
4]. The future realization of the Apollon laser [5], will pave the way to
novel LWFA experiments with high laser power.

The importance of modeling in the LWFA field can hardly be overestimated,
since it represents a tool of experimental design, physical prediction,
analysis and understanding of the involved phenomena. Nonetheless, LWFA
modeling implies sensitive numerical challenges, which require the use of High
Performance Computing (HPC) techniques. The state-of-the-art type of codes
used to model LWFA is the Particle in Cell (PIC) [6] code, which samples the
plasma distribution function through macro-particles (MP), pushed by and
generating electromagnetic (EM) fields defined on a mesh. Vlasov’s equation is
thus solved following its characteristics, i.e. the MP equations of motion.
The PIC method self-consistently updates at each iteration the MP positions
and momenta and the EM fields they generate, as well as external EM fields
such as a laser pulse. In LWFA simulations, the largest scales to consider are
defined by the size of the plasma accelerator. For centimeters long
accelerating stages, the total propagation time is of the order of hundreds of
picoseconds. On the other hand, the smallest scales are usually defined by the
laser. Indeed, the standard PIC method imposes to resolve the laser wavelength
$\lambda_{0}$ and period $2\pi/\omega_{0}$ which are respectively of the order
of the micron and the femtoseconds. The discrepancy in scales quickly
demonstrates that realistic 3D or quasi-3D [7] simulations of this nature are
very costly. They even are beyond reach of standard laser techniques for
accelerating stages reaching tens of centimeters as is planned for multi-stage
LWFA. Henceforth, we will call standard laser techniques, standard laser
simulations those not using physical approximations to speed-up the
calculations. Examples of these techniques include but are not limited to
quasi-static approximation [8], Azimuthal Fourier decomposition [9], boosted
frame techniques [10] or hybrid models [11, 12, 13].

A possible technique to reduce the computing cost of LWFA simulations consists
in using a description of the laser pulse that takes into account only its
complex envelope, with length and transverse size of the order of the plasma
wavelength $\lambda_{p}\gg\lambda_{0}$, without the need to resolve each
optical cycle of the laser in time and space [8, 14, 15, 16, 12, 17]. The code
WAKE was one of the first codes with a laser envelope model to simulate laser
wakefield scenarios [8]. However, the quasi-static approximation (QSA) of the
code prevented the simulation of electrons injection. Time explicit (i.e.
without QSA) envelope codes were developed, with implicit schemes to solve the
evolution equation of the laser envelope [14, 16, 12, 17]. The use of implicit
solvers to solve an envelope equation requires non-trivial parallelization
techniques to be used in 3D, as they often require the inversion of a matrix
representing a transverse differential operator acting on the envelope [16,
12, 17, 18]. Recently, a time explicit 3D envelope code with an easily
parallelizable explicit solver for the envelope equation has been developed in
the PIC code ALaDyn [19].

In this paper, we present applications of this technique now also implemented
in the PIC code Smilei [20, 21], to demonstrate its suitability for studies
oriented to the LWFA experiments planned for Apollon.

In the second section, we briefly review the envelope model’s approximations
and equations. In the third section a numerical simulation of a second stage
experiment is discussed, comparing its results with a standard laser
simulation. In the fourth section, five simulations of possible working points
for the first injecting stage for Apollon are presented. In the Appendix, the
envelope model equations and their solution in Smilei are summarized.

##  2 Review of the envelope model

In many typical situations for LWFA, the spatial and temporal scales of
interest (e.g. the wavelength of the accelerating plasma wave, scaling as
$\lambda_{p}$) are significantly larger than the scales related to the laser
central wavelength $\lambda_{0}$. In these cases, if the laser pulse is much
longer than $\lambda_{0}$, the computation time can be substantially reduced
if one may sample only the laser envelope characteristic length (which
normally scales as $\lambda_{p}$ as well in typical LWFA setups) instead of
$\lambda_{0}$, as depicted in Fig. 1. With this lower resolution, the region
of interest for the electron acceleration, i.e. the end of the plasma “bubble”
in the first period of the plasma wave behind the laser, would still be
correctly described in most of the situations of interest for LWFA. Indeed,
the bubble formation is triggered by an averaged interaction of the particles
with the laser, i.e. the ponderomotive force [22], which derives from the
laser envelope.

Figure 1: Blue line: laser vector potential component $\hat{A}$ along the
polarization direction. Red line: the module of its complex envelope
$|\tilde{A}|$. The envelope is sampled by a number of points smaller by a
factor ten compared to $\hat{A}$.

The envelope model implemented in Smilei is similar to the one first
demonstrated in the PIC code ALaDyn, including the same solver for the
envelope equation in laboratory frame coordinates and the ponderomotive solver
for the particles’ equations of motion presented in [16, 19]. In the
following, the equations of this model are presented. Various numerical
schemes for their solution are detailed in [19] and those implemented in
Smilei are reviewed in A. Henceforth, normalized units will be used for all
quantities (choosing $k_{0}^{-1}$ as the normalized unit length, $c$ as the
normalized unit velocity, etc.).

The fundamental assumption of the model is the description of the laser pulse
vector potential in the complex polarization direction $\hat{A}(\mathbf{x},t)$
as a slowly varying envelope $\tilde{A}(\mathbf{x},t)$ modulated by fast
oscillations at wavelength $\lambda_{0}$:

$\hat{A}(\mathbf{x},t)=\textrm{Re}\left[\tilde{A}(\mathbf{x},t)e^{i(x-t)}\right],$
(1)

Thus, in general any physical quantity $A$ will be therefore given by the
summation of a slowly varying part $\bar{A}$ and a fast oscillating part
$\hat{A}$, i.e. $A=\bar{A}+\hat{A}$, where $\hat{A}$ has the same structure as
in Eq. 1 (we are using the same notation used in [16]). The laser vector
potential in vacuum would have only the fast oscillating part ($\bar{A}=0$).
In this context, “slowly varying” means that the space-temporal variations of
$\bar{A}$ and of the envelope $\tilde{A}$ of the fast oscillating part are
small enough to be treated perturbatively with respect to the ratio
$\epsilon=\lambda_{0}/\lambda_{p}$, as described in detail in [8, 22, 16]. The
laser envelope transverse size $R$ and longitudinal size $L$ are thus assumed
to scale as $R\approx L\approx\lambda_{0}/\epsilon$ [8, 22]. As described
thoroughly in the same references, the action of the laser envelope on the
plasma particles is modeled through the addition of a ponderomotive force term
in the particles equations of motion. This term, not representing a real
force, is a term rising from an averaging process in the perturbative
treatment of the particles motion over the laser optical cycles. The particles
equations of motions (in this case for electrons) thus read:

$\displaystyle\frac{d\mathbf{\bar{x}}}{dt}$ $\displaystyle=$
$\displaystyle\frac{\bar{\mathbf{p}}}{\bar{\gamma}}$ (2)
$\displaystyle\frac{d\mathbf{\bar{p}}}{dt}$ $\displaystyle=$
$\displaystyle-\left(\mathbf{\bar{E}}+\frac{\bar{\mathbf{p}}}{\bar{\gamma}}\times\mathbf{\bar{B}}\right)-\frac{1}{\bar{\gamma}}\nabla\Phi$
(3) $\displaystyle\bar{\gamma}$ $\displaystyle=$
$\displaystyle\sqrt{1+|\mathbf{\bar{p}}|^{2}+\Phi}$ (4)

where $\mathbf{\bar{x}}$, $\mathbf{\bar{p}}$, $\mathbf{\bar{E}}$,
$\mathbf{\bar{B}}$ are the slow-varying parts of the particles positions and
momenta and of the electric and magnetic field. The ponderomotive potential is
defined as $\Phi=|A|^{2}/2$. The Lorentz factor in the usual motion equations
is replaced by the ponderomotive Lorentz factor $\bar{\gamma}$ [8, 22]. The
quantities stored in an envelope simulation in Smilei are the slowly varying
MP positions and momenta $\mathbf{\bar{x}}$, $\mathbf{\bar{p}}$, the slowly
varying electromagnetic fields $\mathbf{\bar{E}}$, $\mathbf{\bar{B}}$, the
envelope $\tilde{A}$, the ponderomotive potential $\Phi$ and its gradient
$\nabla\Phi$. In a standard laser simulation, the complete (high frequency
part and low frequency part) MP positions and momenta $\mathbf{x}$,
$\mathbf{p}$ and electromagnetic fields $\mathbf{E}$, $\mathbf{B}$ are stored.

The evolution of the laser pulse envelope is derived combining d’Alembert’s
inhomogeneous equation and Eq. 1:

$\nabla^{2}\tilde{A}+2i\left(\partial_{x}\tilde{A}+\partial_{t}\tilde{A}\right)-\partial^{2}_{t}\tilde{A}=\chi\tilde{A}.$
(5)

The function $\chi$ represents the plasma susceptibility, which in a PIC code
can be computed similarly to the charge density [12, 19]. This term takes into
account the effect of the plasma on the laser propagation, and is necessary to
model phenomena of self-focusing [23]. The derivation of Eq. 5 assumes that
the high frequency contribution of the scalar potential can be ignored [3,
16].

In Smilei, no further assumption on the envelope is made and the full form of
Eq. 5 is solved. Thus, Eq. 5, solved in ALaDyn and Smilei, is physically
equivalent to the envelope equation in the code INF&RNO (Eq. 1 in [12, 18]),
where it is written in comoving coordinates and solved with an implicit
numerical scheme detailed in [18]. As explained in [12], we remark that
retaining all the terms in the derivation of the envelope equation allows in
principle to use this model in conjunction with the Lorentz boosted frame
method [10], since the form of Eq. 5 is invariant under Lorentz
transformations.

The slowly varying electromagnetic fields $\bar{\mathbf{E}}$,
$\bar{\mathbf{B}}$ are updated solving Maxwell’s equations, with the slowly
varying current densities $\bar{\mathbf{J}}$ as source terms. Since the form
of Maxwell’s equations is unaltered, a Finite Difference Time Domain (FDTD)
scheme [24] is used to evolve these fields in the simulations described
hereafter. The deposition of $\bar{\mathbf{J}}$ is performed with the charge
conserving scheme by Esirkepov [25].

##  3 Benchmark case study: second stage simulations

An envisioned experiment for Apollon is the injection of electrons from a
first nonlinear plasma stage and their transport line towards a second plasma
stage, where they will be accelerated by plasma waves in the weakly nonlinear
regime. The lower accelerating gradients of weakly nonlinear regimes imply the
requirement of a large accelerating length for the second plasma stage to have
a sensitive energy gain. This represents a considerable challenge from the
experimental point of view, but also for the numerical simulation, which would
be unfeasible with a standard laser simulation techniques. In order to
demonstrate the suitability of the envelope model for second stage
simulations, we report a comparison between a standard laser and an envelope
simulation. The physical setup is given by a laser pulse injected in a plasma
and an idealized Gaussian electron beam injected from outside the plasma, as
for an external injection experiment into a second stage. The ideal laser,
plasma and beam parameters are the same as in [26] and are briefly recalled in
the following.

The Gaussian laser pulse, linearly polarized in the $y$ direction, has a waist
size $w_{0}=45$ $\mu$m. Its Gaussian temporal profile has an initial FWHM
duration $\tau_{0}=108$ fs in intensity and peak normalized field amplitude
$a_{0}=\sqrt{2}$. The idealized plasma profile is given by a transversely
parabolic profile $n_{e}(r)=n_{0}\left(1+\frac{\Delta
n}{n_{0}}\frac{r^{2}}{r^{2}_{0}}\right)$ (with $r$ the distance to axis,
$n_{0}=1.5\cdot 10^{17}$ cm-3, $\frac{\Delta n}{n_{0}}=0.25$, $r_{0}=45$
$\mu$m), starting from the initial right boundary of the moving window, which
has longitudinal and transverse dimensions $L_{x}=1400$ $k_{0}^{-1}$ and
$L_{y}=L_{z}=1600$ $k_{0}^{-1}$. The laser is focused at the beginning of the
plasma profile, and its initial center is chosen at a distance
$2\sqrt{2}\thinspace\tau_{0}$ from the plasma. The longitudinal cell size for
the standard laser simulation is $\Delta x_{\rm laser}=0.196$ $k_{0}^{-1}$ and
the transverse cell size is $\Delta y=\Delta z=3$ $k_{0}^{-1}$. The time step
has been chosen as $\Delta t_{\rm laser}=0.95\thinspace\Delta x_{\rm laser}$.
For the envelope simulation, the longitudinal grid cell size and the
integration time step have been set to $\Delta x_{\rm
envelope}=16\thinspace\Delta x_{\rm laser}$, $\Delta t_{\rm
envelope}=0.8\Delta x_{\rm envelope}$ respectively. In both the simulations, 8
particles per cell have been used.

According to our tests, the larger longitudinal cell size and integration
timestep in the envelope simulation seed the numerical Cherenkov radiation
(NCR) more quickly than in the standard laser simulations, as expected from
the underlying theory [27]. This numerical artifact has detrimental effects on
the beam emittance after long propagation distances. To partially cope with
the NCR, we use a binomial filter (2 passes) on the current densities at each
iteration [28]. An advantage of the envelope model, as well as of boosted
frame simulations as explained in [28], is that the laser is not modeled with
high frequency oscillations, so low-pass filtering does not risk to damp the
high frequency phenomena close to the laser as it would do in a standard laser
simulation [28]. Therefore, in the standard laser simulation no low-pass
filter is used.

The electron beam injected after the laser has an ideal Gaussian density
profile in all the phase space sub-planes and a total charge of $Q=30$ pC. Its
mean energy is $150$ MeV, with a $0.5\%$ rms energy spread. The beam
transverse rms size is $1.3$ $\mu$m, its longitudinal rms size $2$ $\mu$m,
with a transverse normalized emittance of $1$ mm-mrad. The beam initial
position is at a distance $0.75\lambda_{p}$ after the laser pulse, in a phase
with both a focusing and an accelerating field in the wake of the laser in the
plasma. All the electron beam MP carry the same electric charge, given by
$Q/N_{p}$, where $N_{p}=10^{6}$ is the number of MP which sample the beam. To
self-consistently initialize the EM fields of an electron beam that is already
relativistic at the beginning of the simulation, the procedure described in
[29, 30] was implemented in Smilei. To summarize, at the beginning of the
simulation the “relativistic Poisson’s equation” is solved, once the initial
beam charge density $\bar{\rho}$ is computed:

$\left(\frac{1}{\gamma^{2}_{0}}\partial^{2}_{x}+\nabla_{\perp}^{2}\right)\bar{\Phi}=-\bar{\rho}.$
(6)

Then, the low frequency EM fields $\mathbf{\bar{E}}$, $\mathbf{\bar{B}}$ at
the same time step $t=0$ are computed:

$\displaystyle\mathbf{\bar{E}}$ $\displaystyle=$
$\displaystyle\left(-\frac{1}{\gamma_{0}^{2}}\partial_{x},-\partial_{y},-\partial_{z}\right)\bar{\Phi},$
(7) $\displaystyle\mathbf{\bar{B}}$ $\displaystyle=$
$\displaystyle\beta_{0}\mathbf{\hat{x}}\times\mathbf{\bar{E}}.$ (8)

The quantities $\gamma_{0}=1/\sqrt{1-\beta^{2}_{0}}$ and $\beta_{0}$ represent
the initial mean beam Lorentz factor and initial normalized speed
respectively. The quantities $\bar{\rho}$ and $\bar{\Phi}$ are centered on the
primal grid in all the $x$, $y$, $z$ directions. The derivatives in the
previous equations are computed through finite differences. In order to
provide the properly space-centered fields for the FDTD scheme, the magnetic
field is then spatially interpolated in the Yee cell. Besides, in order to
provide the properly time-centered initial conditions for the FDTD scheme, the
magnetic field at time $-\Delta t/2$ is found using a backward FDTD “advance”
by $-\Delta t/2$.

Table 1 reports the electron beam parameters obtained by the two simulations
after a propagation distance of $15$ mm. We note a very good agreement for
certain parameters as the mean energy, the transported charge and the beam
duration after this distance. The most striking differences are found in the
energy spread and in the transverse plane parameters as the rms sizes
$\sigma_{i}$ ($i=y,\thinspace z$) and normalized emittances
$\varepsilon_{n,i}$. The normalized emittance is defined as
$\varepsilon_{n,i}=\sqrt{\sigma_{i}^{2}\sigma_{p_{i}}^{2}-\sigma_{ip_{i}}^{2}}$
($i=y,\thinspace z$), where $\sigma_{p_{i}}$ and $\sigma_{ip_{i}}$ are the
beam rms momentum spread and covariance between the coordinate $i$ and the
momentum $p_{i}$ in the transverse planes $y-p_{y}$ and $z-p_{z}$.

The differences in the beam emittance, which is conserved along the
propagation in the envelope simulation, is mainly due to the mitigation of NCR
thanks to the smoothing operated on the current density. The possibility of
low-pass filtering, and thus reduction of the effects of this numerical
artifact plaguing standard LWFA PIC simulations, without compromising the
laser propagation represents an important advantage of envelope simulations.
PIC simulations of LWFA not using envelope need more advanced techniques to
cope with NCR, like ad-hoc solvers for Maxwell’s Equations [31, 27, 32]
instead of the FDTD solver, or pseudo-spectral electromagnetic solvers [33]
coupled to the Galilean frame method [34].

The differences in energy can be explained noting that the envelope solver has
a minor numerical dispersion compared to a standard FDTD scheme [19]. It was
already shown in [32] that differences in the dispersion relation in the
laser/envelope propagation scheme can lead to significant differences in the
longitudinal phase space evolution of laser wakefield-accelerated beams.
Indeed, as explained in [32], less accurate description of the laser
propagation speed because of numerical dispersion can lead to an
overestimation of the dephasing effect hence an underestimation of the
accelerating field. From Fig. 2, where the longitudinal electric field on axis
is compared at $15$ mm, it is evident that the laser pulse is moving faster in
the envelope simulation, and after a long simulation time the difference in
the traveled distance becomes sensitive. Instead, the electron beam moves at
the same velocity (close to $c$) in both simulations and is found at the same
position (see Fig. 3) at the same simulation time. Thus, the electron beam is
subject to a numerically exaggerated dephasing in the standard laser
simulation. As can be seen from Fig. 2, the beam loading is sensitive to the
different wakefield phases where the electron beam is found in the two
simulations, where the laser propagates ad different speed. These differences
in beam loading and in dephasing yield a lower energy and a lower energy
spread in the standard laser simulation, clearly visible in the longitudinal
phase space distribution (Fig. 4).

To summarize, in these particular conditions the simulation of the second
stage can greatly benefit from a modeling based on the envelope approximation.
Compared to standard laser simulations, calculations using laser envelopes can
use digital filters to cope with NCR without damping the relevant physics near
the laser and the reduced numerical dispersion can yield more accurate
results. Besides, the envelope simulations only need a fraction ( 5% in this
case) of the computing resources required by a standard laser simulation. The
necessary computing resources for this simulation up to 15 mm with a standard
laser are not easily available normally and were provided by the Grand
Challenge ”Irene” 2018 project of GENCI.

| Initial values | Standard Laser (15 mm) | Envelope (15 mm)
---|---|---|---
$Q$($pC$) | 30 | 29.98 | 29.94
$\sigma_{x}$ ($\mu$m) | 2.0 | 2.0 | 2.0
$\sigma_{y}$ ($\mu$m) | 1.3 | 1.5 | 1.0
$\sigma_{z}$ ($\mu$m) | 1.3 | 1.4 | 1.0
$\varepsilon_{n,y}$ (mm-mrad) | 1.0 | 2.0 | 1.0
$\varepsilon_{n,z}$ (mm-mrad) | 1.0 | 2.1 | 1.0
$E$ (MeV) | 150 | 427 | 438
$\sigma_{E}/E$ (%) | 0.5 | 4.7 | 6.4
Table 1: Beam parameters: charge $Q$, rms sizes $\sigma_{i}$ ($i=x,y,z$),
normalized emittances $\varepsilon_{n,i}$ ($i=y,z$), mean energy $E$, rms
energy spread $\sigma_{E}/E$. First column: beam parameters at the beginning
of the simulation. Second and third columns: beam parameters after 15 mm of
propagation in the standard laser simulation and envelope simulation. Figure
2: Longitudinal electric field $E_{x}$ on axis after $15$ mm of propagation,
comparing the results obtained with the envelope model and the results
obtained through a standard laser simulation.  Figure 3: Plasma density after
$15$ mm of propagation. Top panel: results obtained with the envelope model.
Bottom panel: results obtained through a standard laser simulation. Figure 4:
Normalized longitudinal phase space distribution of the externally injected
and LWFA-accelerated electron beam after a propagation distance of $15$ mm.

##  4 Single stage simulations for Apollon

This Section illustrates the interest of the envelope model to run simulation
campaigns applied to real experimental setups.

In Section 3, we have compared the results of a standard laser simulation and
an envelope simulation for a long distance LWFA benchmark in a weakly
nonlinear regime. Examples of the use of time explicit envelope models for
highly nonlinear regimes with self-injection, like those discussed in this
Section, can be found for example in [11, 16].

A set of four 3D simulations is performed in order to probe the parameter
space accessible to the first Apollon shots. The very first experiments will
be done with increasing laser energy starting from 5 J and up to 15 J. We
expect the laser pulse duration to be between 20 and 30 fs. In order to fix
parameters for a simple single stage LWFA experiment within this window of
laser pulse energy and length, one can use the criteria defined in [35]. The
plasma density window is limited by the diffraction-limited propagation on one
side $\omega_{p}\tau_{L}\geq 1$, and the depletion-limited propagation on the
other $\omega_{p}\tau_{L}\leq\pi/2$. And once a density is chosen, the minimum
power is $P_{-}=10P_{cr}$ in order for the relativistic self-focusing to
counter balance diffraction. We remind that $P_{cr}=16.2(n_{c}/n_{0})$ GW,
where $n_{0}$ and $n_{c}$ are respectively the electron and critical density
respectively [23]. The laser wavelength is $\rm\lambda_{0}=0.8\ \mu m$. The
maximum power is $P_{+}=25P_{cr}$ in order to stay in the self-guided
propagation without critical oscillations of the laser spot size. The laser
pulse length being already fixed, the criteria on $P$ translates directly into
a condition on the energy. This is illustrated in Fig. 5 which shows, for
different pulse lengths, the minimum and maximum energy as a function of
density within the prescribed density window. It is worthwhile to note that,
as the pulse becomes shorter in time, the energy window narrows and the
optimal energy decreases. Future high energy LWFA experiments based on a
similar self-guided laser pulse willing to reach high charges or energy should
therefore make use of longer laser pulses.

Figure 5: Density window prescribed by pulse length. Minimum energy $E_{-}$
and maximum energy $E_{+}$ in order to guarantee proper self-focusing and
propagation are plotted as functions of the plasma density. The plasma density
ranges between $n_{-}$, in order to have self-focusing effectively counter-
balancing diffraction, and $n_{+}$ above which laser depletion makes the
acceleration inefficient and unstable. Optimum parameters are located at the
center of the delimited zone. Accessible parameters with energies of 5, 10 and
15 J are plotted as gray dashed lines. Colored dots mark the parameters chosen
for the scan simulations.

The waist size has been chosen as $40$ $\mu$m and focused at the start of the
plasma in all the simulations in order to better compare the various working
points. The plasma density of all the simulations was defined with a $300$
$\mu$m upramp from $0$ to the plateau density value $n_{0}$, reported in Table
2. After the upramp, the plasma density remains constant. Besides, the plasma
density is chosen uniform in the transverse direction, to rely only on
relativistic self-focusing for the laser guiding [23]. In all the scan
simulations, 8 particles per cell have been used. The integration timestep was
always set to $\Delta t=0.9\thinspace\Delta x$, where $\Delta x$ is the
longitudinal mesh cell size. The transverse mesh cell size has been set to
$\Delta y=\Delta z=3.5$ $k_{0}^{-1}$ for all the simulations. As in the
simulations reported in Section 3, we used a binomial filter on the current
densities at each iteration (2 passes). Again, this smoothing aims at limiting
the NCR impact on the simulation.

We remark that a similar set of equivalent 3D simulations would have been
unfeasible with standard laser simulations. For reader’s convenience, the
parameters of the four simulations have been summarized in table 2, including
the longitudinal mesh cell size.

Simulation | Pulse duration [fs] | $a_{0}$ |  $n_{0}$ $[10^{18}$ cm${}^{-3}]$ | $P/P_{cr}$ | Laser Energy [J] |  $\Delta x$ [$k_{0}^{-1}$]
---|---|---|---|---|---|---
1 | 20 | 2.96 | 1.1 | 18.3 | 10 | 0.5
2 | 25 | 2.64 | 0.9 | 12.0 | 10 | 0.667
3 | 25 | 3.24 | 1.1 | 22.0 | 15 | 0.667
4 | 30 | 2.96 | 0.8 | 13.3 | 15 | 0.8
Table 2: Laser and Plasma Parameters for the scan simulations. The waist size
has been chosen as $40$ $\mu$m for all the four simulations. The laser
strength parameter[3] is denoted by $a_{0}$. The longitudinal mesh cell size
$\Delta x$ is also reported. The transverse mesh cell size has been set to
$\Delta y=\Delta z=3.5$ $k_{0}^{-1}$ and the integration timestep to $\Delta
t=0.9\thinspace\Delta x$.

Figure 6 (top left panel) depicts the evolution of the peak absolute value of
the electric field envelope $|\tilde{E}|$ for the four simulations. From the
definition of the envelope $\tilde{A}$ of the vector potential $\hat{A}$ along
the polarization direction (Eq. 1), this quantity is defined as the amplitude
of the envelope of $\hat{E}=-\partial_{t}\hat{A}$, i.e.
$\tilde{E}=-(\partial_{t}-i)\tilde{A}$. As desired, the ratios $P/P_{cr}\gg 1$
lead to relativistic self-focusing [23], which suddenly enlarges the bubble
behind the laser, triggering electron injection. The importance of dark
current (low energy tail in electrons energy distribution) depends strongly on
the set of parameters used. An ad hoc definition of the beam is used here in
order to separate it from the dark current when necessary.

After injection, the energy distribution of the electrons shows a distinct
peak at $E_{spectrum\thinspace peak}$ and a full width half maximum $\Delta
E$. We define the injected beam, at a given time, as the electrons having an
energy within an interval of width $2\Delta E$ around $E_{spectrum\thinspace
peak}$.

Figure 6 also shows the evolution of the total charge $Q_{tot}$ above $300$
MeV (bottom left panel) and the energy of the spectrum peak
$E_{spectrum\thinspace peak}$ (top right panel). As expected, injection occurs
when the laser self-focuses [35]. Since secondary injections can create
additional peaks in the electron spectrum, the energy of the peak
$E_{spectrum\thinspace peak}$ can vary suddenly when a peak overtakes another
one and becomes the new maximum of the distribution (see red line or green
line of the top right panel in Figure 6).

Figure 6:  Top left panel: Evolution of the absolute value of the laser peak
electric field envelope $|\tilde{E}|$. Top right panel: evolution of the
energy of the electron spectrum peak $E_{spectrum\thinspace peak}$ for
particles with energy above $300$ MeV. Bottom left panel: Evolution of the the
total charge $Q_{tot}$ above $300$ MeV. In all the panels, the colored markers
highlight the “optimum” laser propagation distance, before the electron beam
starts its rotation in the longitudinal phase space

To provide interesting working points for the Apollon experiments, in Table 3
we report the injected beam parameters after a certain laser propagation
distance, different for each simulation. This optimal distance has been chosen
by carefully following the evolution of the injected beam in the longitudinal
phase space $x-p_{x}$. At that point, the electron beam reaches the maximum
spectral density. It corresponds to the moment when phase space rotation
minimizes the beam energy spread [36]. As seen on Fig. 6, the acceleration
process goes on past this point but at the cost of a reduced beam quality and
possible additional dark current. For all the simulations, this optimal
propagation distance is found around 6-7 mm. In all these working points, the
spectrum peak energy is around 1 GeV, with an energy spread lower than $10\%$.
Fig. 7 reports the electron energy spectrum for all the simulations at this
“optimum” laser propagation distance.

Figure 7: Energy spectrum of the single stage simulations at the “optimum” laser propagation distance. The cut-off energy is 300 MeV.  Simulation | 1 | 2 | 3 | 4
---|---|---|---|---
$n_{0}$ [$10^{18}$ cm-3] | 1.1 | 0.9 | 1.1 | 0.8
$a_{0}$ | 2.96 | 2.64 | 3.24 | 2.96
$P/P_{cr}$ | 18.3 | 12.0 | 22.0 | 13.3
Laser Energy [J] | 10 | 10 | 15 | 15
Pulse duration [fs] | 20 | 25 | 25 | 30
Laser propagation distance [mm] | 6.3 | 7.2 | 6.5 | 7.6
$Q_{beam}$ [pC] | 263 | 48 | 543 | 24
$Q_{total,\thinspace>300\thinspace MeV}$ [pC] | 426 | 72 | 729 | 24
$E_{spectrum\thinspace peak}$ [MeV] | 870 | 740 | 930 | 1130
$\Delta E/E_{spectrum\thinspace peak}$ [%] | 8.3 | 3.2 | 6.4 | 2.0
$2\sigma_{y}$ [$\mu$m] | 3.0 | 2.3 | 2.9 | 0.5
$2\sigma_{z}$ [$\mu$m] | 2.8 | 2.3 | 3.1 | 0.5
$\varepsilon_{n,y}$ [mm-mrad] | 14.5 | 1.3 | 12.3 | 0.4
$\varepsilon_{n,z}$ [mm-mrad] | 11.8 | 1.3 | 13.1 | 0.4
Table 3: Resume of the plasma plateau density $n_{0}$, initial laser pulse
parameters used in the set of first stage simulations. For the electron beam,
the peak in the energy spectrum $E_{spectrum\thinspace peak}$, the relative
energy spread $\Delta E/E_{spectrum\thinspace peak}$ , rms size (2 standard
deviations $\sigma$) and the normalized emittances $\varepsilon_{n,i}$
($i=y,\thinspace z$) in the transverse planes are reported. These parameters
are measured after a certain laser propagation distance, before the phase
space rotation, reported in the table.

Even in the relatively well constrained range of parameters, significant
differences are observed. First, there is a factor 20 between the minimum and
maximum beam charges. This parameter appears to be controlled primarily by the
$P/P_{\rm cr}$ ratio. Lasers close to the $P/P_{\rm cr}=10$ are subject to
milder self-focusing and trigger less self-injection (simulations 2 and 4).
Conversely, for higher $P/P_{\rm cr}$, self-injection is dramatically enhanced
(simulations 1 and 3). This dependency is expected [35] and provides a
practical way to optically control the injected charge even though the
sensitivity is such that it might be difficult to achieve a very good
accuracy.

The energy evolution is roughly similar in all cases with the exception of
simulation 2. This evolution is governed by the accelerating field $E_{\rm
acc}$. In the bubble regime, $E_{\rm acc}$ is proportional to
$\sqrt{a_{0}n_{0}}$[3]. In the simulations, $n_{0}$ is constant and therefore,
$a_{0}$ evolution dictates the evolution of the accelerating field. This
evolution is shown as the normalized peak $|\tilde{E}|$ in figure 6.
Simulation 2 has a small density and the weakest $a_{0}$ throughout the run
because of a low inital $a_{0}$ but also the weakest self-focusing. This
explains the small acceleration observed in that case with respect to the
other simulations. Simulation 4, performs quite well in terms of acceleration
in spite of having the smallest density and only a medium initial $a_{0}$
thanks to a superior laser guiding. Simulation 3, in spite of having the
highest $a_{0}$ and density do not outperform other simulations in terms of
energy because in that case the beam loading effect is not negligible any
more[37].

The computing time needed by the 3D simulations for 1 mm is of the order of
30-60 kh. We remark that a parameter scan of four 12 mm long 3D simulations as
the one presented in this work would have been extremely costly with standard
laser simulations, which are typically slower by at least a factor 20. This
justifies the interest in time-explicit envelope models [12, 16], especially
for preliminary studies for experiments like those planned for Apollon.

##  5 Conclusions

A recently published, time explicit, easily parallelizable envelope method for
PIC codes has the potential to model laser-plasma interaction involved in many
LWFA setups with considerable speedups compared to standard laser simulations.
The implementation of this time-explicit model in Smilei has been described.
The suitability of this envelope model for single-stage with self-injection
and second-stage LWFA simulations has been explored. The laser and plasma
parameters chosen for these studies lie within the intervals of interest for
the upcoming Apollon LWFA experiments.

The envelope model used in this work has been benchmarked against a standard
laser simulation in a scenario with external injection of a witness electron
beam, yielding a speedup of 20 compared to a standard laser simulation and a
good agreement even at a distance of 15 mm. Besides, the envelope model allows
to model more accurately the longitudinal and transverse phase space evolution
of the injected beam. Indeed, in the longitudinal direction a more accurate
prediction of the laser propagation speed by the envelope equation reduces the
numerical dephasing of the electron beam caused by FDTD solvers of Maxwell’s
equations. In the transverse direction, the envelope simulations allow for
low-pass filtering which in this case considerably reduced the growth of NCR,
conserving the beam emittance.

In the single-stage simulations, four potential working points in the regions
of parameters of interest have been tested, finding the injection of $\approx
1$ GeV electron beams with energy spread lower than $10\%$. These working
points have been found examining the evolution of the injection and
acceleration process for 12 mm in four 3D simulation, a propagation distance
which would need considerable computing resources to be simulated with a
standard 3D PIC code.

The use of this envelope model may thus represent an important investigation
tool for the study of LWFA and parameter explorations with reduced resource
requirements. Further improvements could include spectral solvers for the
envelope equation and Maxwell’s equations to cope with numerical Cherenkov
radiation and the use of cylindrical symmetry to speed up even further the
simulations (albeit losing the full 3D characteristics of the phenomena).

##  Appendix A Envelope PIC loop

We report a brief summary of the numerical solution of the equations involved
in the envelope model used for the simulations of this work. More details on
the derivation of the equations and the numerical schemes used to solve them
can be found in [19]. The envelope equation that Smilei solves is

$\nabla^{2}\tilde{A}+2i\left(\partial_{x}\tilde{A}+\partial_{t}\tilde{A}\right)-\partial^{2}_{t}\tilde{A}=\chi\tilde{A},\quad$
(9)

in which the the susceptibility $\chi$ is defined as

$\chi(\mathbf{x})=\sum_{s}\,r_{s}^{2}\,\sum_{p}\,\frac{w_{p}}{\bar{\gamma}_{p}}\,S\big{(}\mathbf{x}-\mathbf{\bar{x}}_{p}\big{)}\
,$ (10)

where $\bar{\gamma}_{p}$ is the averaged or ponderomotive Lorentz factor of
the macroparticle $p$ and $r_{s}=q_{s}/m_{s}$ is the charge to mass ratio
($q_{s}$ and $m_{s}$ are respectively the particle species $s$ charge and mass
normalized by the elementary charge $e$ and electron mass $m_{e}$). The MP
weight and shape factor are denoted by $w_{p}$ and $S(\mathbf{x})$
respectively (see [6] for the definition of shape factor). The averaged
Lorentz factor is defined from the averaged particle momentum
$\mathbf{\bar{u}}_{p}=\mathbf{\bar{p}}_{p}/m_{s}$ and the ponderomotive
potential $\Phi=|\tilde{A}|^{2}/2$ [22, 8, 16]:

$\bar{\gamma}_{p}=\sqrt{1+\mathbf{\bar{u}}_{p}\cdot\mathbf{\bar{u}}_{p}+r_{s}^{2}\Phi(\mathbf{\bar{x}}_{p})}.$
(11)

Maxwell’s equations retain their form, except for substituting the
electromagnetic fields and the source terms with their respective low
frequency components (denoted with a bar):

$\displaystyle\partial_{t}\mathbf{\bar{E}}$ $\displaystyle=$
$\displaystyle\nabla\times\mathbf{\bar{B}}-\mathbf{\bar{J}}$ (12)
$\displaystyle\partial_{t}\mathbf{\bar{B}}$ $\displaystyle=$
$\displaystyle-\nabla\times\mathbf{\bar{E}}.$ (13)

The low frequency current density $\mathbf{\bar{J}}$ is projected on the grid
through Esirkepov’s method [25].

Also the MP equations of motion remain similar to the usual ones, but with a
crucial difference: the low frequency components of the position and momentum
are pushed by both the low frequency electromagnetic fields (the Lorentz force
term) and the ponderomotive force
$\mathbf{F}_{pond}=-r^{2}_{s}\thinspace\frac{1}{4\bar{\gamma}_{p}}\nabla\Phi$.
This term takes into account the averaged effect of the laser on the particles
and can be computed from the envelope $\tilde{A}$. The resulting equations of
motions for the particle of index $p$ are:

$\displaystyle\frac{d\mathbf{\bar{x}}_{p}}{dt}$ $\displaystyle=$
$\displaystyle\frac{\mathbf{\bar{u}}_{p}}{\bar{\gamma}_{p}}\,$ (14)
$\displaystyle\frac{d\mathbf{\bar{u}}_{p}}{dt}$ $\displaystyle=$
$\displaystyle
r_{s}\,\left(\mathbf{\bar{E}}_{p}+\frac{\mathbf{\bar{u}}_{p}}{\bar{\gamma}_{p}}\times\mathbf{\bar{B}}_{p}\right)-r^{2}_{s}\thinspace\frac{1}{4\bar{\gamma}_{p}}\nabla\Phi_{p}.$
(15)

The nonlinearities introduced in the envelope source term (Eqs. 10,11 ) and
the ponderomotive force by the terms $\bar{\gamma}$ (Eq. 11) and $\nabla\Phi$
require a modification of the standard PIC algorithm. The PIC loop is changed
to implement the solution of the ponderomotive equations (Eqs. 9,14), adding
the envelope equation solver, the susceptibility deposition and a
ponderomotive particle push, as depicted in Fig.8.

Figure 8: Ponderomotive Particle in Cell loop. The definition of ponderomotive
potential $\Phi=\frac{|\tilde{A}|^{2}}{2}$ is used.

The interpolation of the fields to the MP positions and the deposition of the
current density and susceptibility on the grid are implemented as in standard
PIC codes [6]. In the simulations presented in this work, we used an order 2
shape function, without loss of generality.

As explained in detail in [19], Eq. 5 can be discretized using finite
differences both in time and space. The second order centered finite
differences yield an explicit scheme:

$\displaystyle\tilde{A}_{ijk}^{n+1}=\frac{1+i\Delta t}{1+\Delta
t^{2}}\Bigg{[}2\tilde{A}_{ijk}^{n}-(1+i\Delta
t)\tilde{A}_{ijk}^{n-1}+\Bigg{(}\nabla^{2}\tilde{A}_{ijk}\bigg{|}_{ijk}^{n}-\chi^{n}_{ijk}\tilde{A}_{ijk}^{n}+2i\Delta
t^{2}\frac{\tilde{A}_{i+1\thinspace jk}^{n}-\tilde{A}_{i-1\thinspace
jk}^{n}}{2\Delta z}\Bigg{)}\Bigg{]}\quad$ (16)

where

$\displaystyle\nabla^{2}\tilde{A}\bigg{|}_{ijk}^{n}=$
$\displaystyle\frac{\tilde{A}_{i+1\thinspace
jk}^{n}-2\tilde{A}_{ijk}^{n}+\tilde{A}_{i-1\thinspace jk}^{n}}{\Delta
x^{2}}+\frac{\tilde{A}_{ij-1\thinspace
k}^{n}-2\tilde{A}_{ijk}^{n}+\tilde{A}_{ij+1\thinspace k}^{n}}{\Delta
y^{2}}+\frac{\tilde{A}_{ijk-1}^{n}-2\tilde{A}_{ijk}^{n}+\tilde{A}_{ijk+1}^{n}}{\Delta
z^{2}}.$ (17)

The indices $i$,$j$,$k$ refer to the mesh cell indices along the $x$, $y$, $z$
directions. The envelope, the ponderomotive potential and the susceptibility
are centered in space as the charge density in all directions and time-
centered as the electric field in the Yee scheme [24] .

Once the susceptibility term at time step $n$ is known, the past and present
envelope fields ($\tilde{A}^{n-1}$, $\tilde{A}^{n}$ respectively) can be used
to compute $\tilde{A}^{n+1}$ using Eq. 16. The stencil of this explicit
envelope solver contains three points in space for each direction.

To project the susceptibility on the grid and to update the MP momenta , as
explained in [19], the ponderomotive Lorentz factor at the timestep $n$ is
necessary. This quantity cannot be directly computed from Eq. 11, since the MP
momentum is known only at the timestep $n-1/2$. Thus, a two step computation
(derived in [19]) is used to obtain an approximation of the desired Lorentz
factor $\bar{\gamma}_{p}^{n}$:

$\displaystyle\gamma_{0p}^{2}$ $\displaystyle=$ $\displaystyle
1+r_{s}^{2}\Phi_{p}^{n}+\mathbf{\bar{u}}_{p}^{n-1/2}\cdot\mathbf{\bar{u}}_{p}^{n-1/2},$
(18) $\displaystyle\bar{\gamma}_{p}^{n}$ $\displaystyle=$
$\displaystyle\gamma_{0p}+\frac{1}{2\gamma_{0p}^{2}}\Bigg{(}\gamma_{0p}\frac{\Delta
t}{2}r_{s}\mathbf{\bar{E}}^{n}-\frac{\Delta
t}{4}r_{s}^{2}\nabla\Phi_{p}^{n}\Bigg{)}\cdot\mathbf{\bar{u}}_{p}^{n-1/2}.$
(19)

This ponderomotive Lorentz factor $\bar{\gamma}_{p}^{n}$ is used in the
deposition of susceptibility, following Eq. 10.

A modified Boris pusher [6] can be used to update the MP momenta. The only
modifications to the standard scheme are the use of
$\left(r_{s}\mathbf{\bar{E}}_{p}^{n}-r_{s}^{2}\frac{\nabla\Phi^{n}_{p}}{2\bar{\gamma}_{p}^{n}}\right)$
instead of $r_{s}\mathbf{\bar{E}}_{p}^{n}$ as the source term changing the
particle $p$ energy and the term
$r_{s}\frac{\mathbf{\bar{B}}^{n}_{p}}{\bar{\gamma}_{p}^{n}}$ for its momentum
rotation under the effect of the magnetic field. The ponderomotive Lorentz
factor $\bar{\gamma}_{p}^{n}$ used by the pusher is the one from Eq. 19.

Again, to update the particle $p$ position through a leapfrog scheme, the
necessary ponderomotive Lorentz factor at timestep $n+1/2$ cannot be computed
directly from Eq. 11, since the ponderomotive potential $\Phi$ and its
gradient can be interpolated only at the known MP position at timestep $n$.
Hence, a two-step approximation of the ponderomotive Lorentz factor
$\bar{\gamma}_{p}^{n+1/2}$ (derived in [19]) is used:

$\displaystyle\gamma_{0p}^{2}$ $\displaystyle=$ $\displaystyle
1+\mathbf{\bar{u}}_{p}^{n+1/2}\cdot\mathbf{\bar{u}}_{p}^{n+1/2}+r_{s}^{2}\Phi_{p}^{n+1/2}(\mathbf{\bar{x}}_{p}^{n}),$
(20) $\displaystyle\bar{\gamma}^{n+1/2}_{p}$ $\displaystyle=$
$\displaystyle\gamma_{0p}+\frac{\Delta
t}{4\gamma^{2}_{0p}}\mathbf{\bar{u}}_{p}^{n+1/2}\cdot
r_{s}^{2}\nabla\Phi_{p}^{n+1/2}(\mathbf{\bar{x}}_{p}^{n}).$ (21)

Note that the ponderomotive potential and its gradient in Eqs. 20, 21 are
defined at the timestep $n+1/2$. Since these quantities known at the timesteps
$n$ and $n+1$, their value at $n+1/2$ can be obtained through linear
interpolation. After this interpolation in time, they are interpolated at the
known particle position at the timestep $n$ and then used to compute
$\bar{\gamma}_{p}^{n+1/2}$ through Eqs. 20, 21.

Finally, using $\bar{\gamma}_{p}^{n+1/2}$ of Eq. 21, the updated position can
be computed:

$\mathbf{\bar{x}}_{p}^{n+1}=\mathbf{\bar{x}}_{p}^{n}+\frac{\mathbf{\bar{u}}_{p}^{n+1/2}}{\bar{\gamma}^{n+1/2}_{p}}\Delta
t.$ (22)

## Acknowledgments

F. Massimo was supported by P2IO LabEx (ANR-10-LABX-0038) in the framework
“Investissements d’Avenir” (ANR-11-IDEX-0003-01) managed by the Agence
Nationale de la Recherche (ANR, France). This work was granted access to the
HPC resources of TGCC/CINES under the allocation 2018-A0050510062 and Grand
Challenge “Irene” 2018 project gch0313 made by GENCI. The authors are grateful
to the TGCC and CINES engineers for their support. The authors thank the
engineers of the LLR HPC clusters for resources and help. The authors are
grateful to the ALaDyn development team for the help and discussions during
the development of the envelope model, in particular D. Terzani, A. Marocchino
and S. Sinigardi and to Gilles Maynard for fruitful discussions.

## Bibliography

## References

  * [1] T. Tajima and J. M. Dawson.  Laser electron accelerator.  Phys. Rev. Lett., 43:267–270, Jul 1979.
  * [2] V. Malka, S. Fritzler, E. Lefebvre, M.-M. Aleonard, F. Burgy, J.-P. Chambaret, J.-F. Chemin, K. Krushelnick, G. Malka, S. P. D. Mangles, Z. Najmudin, M. Pittman, J.-P. Rousseau, J.-N. Scheurer, B. Walton, and A. E. Dangor.  Electron acceleration by a wake field forced by an intense ultrashort laser pulse.  Science, 298(5598):1596–1600, 2002.
  * [3] E. Esarey, C. B. Schroeder, and W. P. Leemans.  Physics of laser-driven plasma-based electron accelerators.  Rev. Mod. Phys., 81:1229–1285, Aug 2009.
  * [4] V. Malka.  Laser plasma accelerators.  Physics of Plasmas, 19(5):055501, 2012.
  * [5] B. Cros, B.S. Paradkar, X. Davoine, A. Chancé, F.G. Desforges, S. Dobosz-Dufrénoy, N. Delerue, J. Ju, T.L. Audet, G. Maynard, M. Lobet, L. Gremillet, P. Mora, J. Schwindling, O. DelferriÃšre, C. Bruni, C. Rimbault, T. Vinatier, A. Di Piazza, M. Grech, C. Riconda, J.R. MarquÃšs, A. Beck, A. Specka, Ph. Martin, P. Monot, D. Normand, F. Mathieu, P. Audebert, and F. Amiranoff.  Laser plasma acceleration of electrons with multi-pw laser beams in the frame of cilex.  Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 740:27 – 33, 2014.  Proceedings of the first European Advanced Accelerator Concepts Workshop 2013.
  * [6] C. K. Birdsall and A. B. Langdon.  Plasma Physics via Computer Simulation.  Taylor and Francis Group, 2004.
  * [7] X. Davoine, E. Lefebvre, J. Faure, C. Rechatin, A. Lifschitz, and V. Malka.  Simulation of quasimonoenergetic electron beams produced by colliding pulse wakefield acceleration.  Physics of Plasmas, 15(11):113102, 2008.
  * [8] Patrick Mora and Jr. Thomas M. Antonsen.  Kinetic modeling of intense, short laser pulses propagating in tenuous plasmas.  Physics of Plasmas, 4(1):217–229, 1997.
  * [9] A. Lifschitz, X. Davoine, E. Lefebvre, J. Faure, C. Rechatin, and V. Malka.  Particle-in-Cell modelling of laser–plasma interaction using Fourier decomposition.  Journal of Computational Physics, 228(5):1803–1814, November 2008\.
  * [10] J.-L. Vay.  Noninvariance of space- and time-scale ranges under a lorentz transformation and the implications for the study of relativistic interactions.  Phys. Rev. Lett., 98:130405, Mar 2007.
  * [11] C. Benedetti, C. B. Schroeder, E. Esarey, C. G. R. Geddes, and W. P. Leemans.  Efficient modeling of laser-plasma accelerators with inf&rno.  AIP Conference Proceedings, 1299(1):250–255, 2010.
  * [12] Carlo Benedetti, Eric Esarey, Wim P. Leemans, and Carl B. Schroeder.  Efficient Modeling of Laser-plasma Accelerators Using the Ponderomotive-based Code INF&RNO.  In Proceedings, 11th International Computational Accelerator Physics Conference (ICAP 2012): Rostock-Warnemünde, Germany, August 19-24, 2012, page THAAI2, 2012.
  * [13] F. Massimo, S. Atzeni, and A. Marocchino.  Comparisons of time explicit hybrid kinetic-fluid code architect for plasma wakefield acceleration with a full pic code.  Journal of Computational Physics, 327(Supplement C):841 – 850, 2016\.
  * [14] D. F. Gordon, W. B. Mori, and T. M. Antonsen.  A ponderomotive guiding center particle-in-cell code for efficient modeling of laser-plasma interactions.  IEEE Transactions on Plasma Science, 28(4):1135–1143, Aug 2000\.
  * [15] C. Huang, V.K. Decyk, C. Ren, M. Zhou, W. Lu, W.B. Mori, J.H. Cooley, T.M. Antonsen, and T. Katsouleas.  Quickpic: A highly efficient particle-in-cell code for modeling wakefield acceleration in plasmas.  Journal of Computational Physics, 217(2):658 – 679, 2006.
  * [16] Benjamin M. Cowan, David L. Bruhwiler, Estelle Cormier-Michel, Eric Esarey, Cameron G.R. Geddes, Peter Messmer, and Kevin M. Paul.  Characteristics of an envelope model for laser–plasma accelerator simulation.  Journal of Computational Physics, 230(1):61 – 86, 2011.
  * [17] A. Helm, J. Vieira, L. Silva, and R. Fonseca.  Implementation of a 3D version of ponderomotive guiding center solver in particle-in-cell code OSIRIS.  In APS Meeting Abstracts, page GP10.011, October 2016.
  * [18] C Benedetti, C B Schroeder, C G R Geddes, E Esarey, and W P Leemans.  An accurate and efficient laser-envelope solver for the modeling of laser-plasma accelerators.  Plasma Physics and Controlled Fusion, 60(1):014002, 2018.
  * [19] Davide Terzani and Pasquale Londrillo.  A fast and accurate numerical implementation of the envelope model for laser-plasma dynamics.  Computer Physics Communications, 2019.
  * [20] J. Derouillat, A. Beck, F. Pérez, T. Vinci, M. Chiaramello, A. Grassi, M. Flé, G. Bouchard, I. Plotnikov, N. Aunai, J. Dargent, C. Riconda, and M. Grech.  Smilei : A collaborative, open-source, multi-purpose particle-in-cell code for plasma simulation.  Computer Physics Communications, 222:351 – 373, 2018.
  * [21] Arnaud Beck, Julien Dérouillat, Mathieu Lobet, Asma Farjallah, Francesco Massimo, Imen Zemzemi, Frédéric Perez, Tommaso Vinci, and Mickael Grech.  Adaptive SIMD optimizations in particle-in-cell codes with fine-grain particle sorting.  Computer Physics Communication, (in press) 2019.
  * [22] Brice Quesnel and Patrick Mora.  Theory and simulation of the interaction of ultraintense laser pulses with electrons in vacuum.  Phys. Rev. E, 58:3719–3732, Sep 1998.
  * [23] G. Z. Sun, E. Ott, Y. C. Lee, and P. Guzdar.  Self focusing of short intense pulses in plasmas.  The Physics of Fluids, 30(2):526–532, 1987.
  * [24] K. Yee.  Numerical solution of inital boundary value problems involving maxwell’s equations in isotropic media.  IEEE Transactions on Antennas and Propagation, 14:302–307, May 1966\.
  * [25] T.Zh. Esirkepov.  Exact charge conservation scheme for particle-in-cell simulation with an arbitrary form-factor.  Computer Physics Communications, 135(2):144 – 153, 2001.
  * [26] Xiangkun Li, Alban Mosnier, and Phu Anh Phi Nghiem.  Design of a 5 gev laser-plasma accelerating module in the quasi-linear regime.  Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 2018.
  * [27] R. Lehe, A. Lifschitz, C. Thaury, V. Malka, and X. Davoine.  Numerical growth of emittance in simulations of laser-wakefield acceleration.  Phys. Rev. ST Accel. Beams, 16(2):021301, 2013.
  * [28] J.-L. Vay, C.G.R. Geddes, E. Cormier-Michel, and D.P. Grote.  Numerical methods for instability mitigation in the modeling of laser wakefield accelerators in a lorentz-boosted frame.  Journal of Computational Physics, 230(15):5908 – 5929, 2011.
  * [29] J.-L. Vay.  Simulation of beams or plasmas crossing at relativistic velocity.  Physics of Plasmas, 15(5):056701, 2008.
  * [30] F. Massimo, A. Marocchino, and A.R. Rossi.  Electromagnetic self-consistent field initialization and fluid advance techniques for hybrid-kinetic pwfa code architect.  Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 829(Supplement C):378 – 382, 2016.  2nd European Advanced Accelerator Concepts Workshop - EAAC 2015.
  * [31] M. Karkkainen, E. Gjonaj, T. Lau, and T. Weiland.  Low-dispersion wake field calculation tools.  Proceedings of ICAP, (1):35.
  * [32] Benjamin M. Cowan, David L. Bruhwiler, John R. Cary, Estelle Cormier-Michel, and Cameron G. R. Geddes.  Generalized algorithm for control of numerical dispersion in explicit time-domain electromagnetic simulations.  Phys. Rev. ST Accel. Beams, 16:041303, Apr 2013.
  * [33] S. Jalas, I. Dornmair, R. Lehe, H. Vincenti, J.-L. Vay, M. Kirchen, and A. R. Maier.  Accurate modeling of plasma acceleration with arbitrary order pseudo-spectral particle-in-cell methods.  Physics of Plasmas, 24(3):033115, 2017.
  * [34] Remi Lehe, Manuel Kirchen, Brendan B. Godfrey, Andreas R. Maier, and Jean-Luc Vay.  Elimination of numerical cherenkov instability in flowing-plasma particle-in-cell simulations by using galilean coordinates.  Phys. Rev. E, 94:053305, Nov 2016.
  * [35] A. Beck, S.Y. Kalmykov, X. Davoine, A. Lifschitz, B.A. Shadwick, V. Malka, and A. Specka.  Physical processes at work in sub-30fs, pw laser pulse-driven plasma accelerators: Towards gev electron acceleration experiments at cilex facility.  Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 740:67 – 73, 2014.  Proceedings of the first European Advanced Accelerator Concepts Workshop 2013.
  * [36] W. Lu, M. Tzoufras, C. Joshi, F. S. Tsung, W. B. Mori, J. Vieira, R. A. Fonseca, and L. O. Silva.  Generating multi-gev electron bunches using single stage laser wakefield acceleration in a 3d nonlinear regime.  Phys. Rev. ST Accel. Beams, 10:061301, Jun 2007.
  * [37] M. Tzoufras, W. Lu, F. S. Tsung, C. Huang, W. B. Mori, Thomas C. Katsouleas, J. Vieira, R. A. Fonseca, and L. O. Silva.  Beam loading in the nonlinear regime of plasma-based acceleration.  Phys. Rev. Lett., 101:145002, 2008.


Permutability of Backlund Transformation for $N=1$ Supersymmetric Sinh-Gordon
111 key words:Backlund Transformation, Supersymmetric sinh-Gordon, non linear
superposition formula
Pacs 02.30.Ik

J.F. Gomes222corresponding author<EMAIL_ADDRESS>L.H. Ymai and A.H.
Zimerman

Instituto de Física Teórica-UNESP
Rua Pamplona 145
fax (55)11 31779080
01405-900 São Paulo, Brazil

###### Abstract

The permutability of two Backlund transformations is employed to construct a
non linear superposition formula to generate a class of solutions for the
$N=1$ super sinh-Gordon model.

Backlund transformations (BT) relating two different soliton solutions are
known to be characteristic of certain class of nonlinear equations. A
remarkable consequence is that from a particular soliton solution, a second
solution can be generated by Backlund transformation. This second solution, in
turn, generates a third one and such structure allows to construct conditions
for the permutability of two sequences of BT.

The study of superposition principle for soliton solutions of the sine-Gordon
equation was employed to show that the order of two BT is, in fact, irrelevant
[1]. Such condition became known as the permutability theorem as was applied
to the KdV and for the $N=1$ super KdV equations in [2] and in [3]
respectively.

Soliton solutions for the $N=1$ super sine-Gordon were obtained in [4] from
the super Hirota’s formalism and in [5] using dressing transformation and
vertex operators. Backlund solutions for super mKdV were considered in [6]. In
this paper we use the superfield approach for the BT as proposed in [7] to
derive a closed algebraic superposition formula for soliton solutions of the
$N=1$ super sinh-Gordon model assuming that two sucessive BT commute.

The model is described by the following equation of motion written within the
superfield formalism [7]

$\displaystyle D_{x}D_{t}\Phi=2i\sinh\Phi,$ (1)

where the bosonic superfield $\Phi$ is given in components by

$\displaystyle\Phi=\phi+\theta_{1}\bar{\psi}+i\theta_{2}\psi-\theta_{1}\theta_{2}2i\sinh\phi,$
(2)

where $\theta_{1}$ and $\theta_{2}$ are Grassmann variables (i.e.
$\theta_{1}^{2}=\theta_{2}^{2}=0$ and
$\theta_{1}\theta_{2}+\theta_{2}\theta_{1}=0$ ) The superderivatives

$\displaystyle D_{x}=\partial_{\theta_{1}}+\theta_{1}\partial_{x},\qquad
D_{t}=\partial_{\theta_{2}}+\theta_{2}\partial_{t},$ (3)

satisfy

$\displaystyle D_{x}^{2}=\partial_{x},\qquad D_{t}^{2}=\partial_{t},\qquad
D_{x}D_{t}=-D_{t}D_{x}.$ (4)

The Backlund transformation for eqn. (1) is given by [7]

$\displaystyle D_{x}(\Phi_{0}-\Phi_{1})$ $\displaystyle=$
$\displaystyle-\frac{4i}{\beta_{1}}f_{0,1}\cosh\left(\frac{\Phi_{0}+\Phi_{1}}{2}\right),$
(5) $\displaystyle D_{t}(\Phi_{0}+\Phi_{1})$ $\displaystyle=$ $\displaystyle
2\beta_{1}f_{0,1}\cosh\left(\frac{\Phi_{0}-\Phi_{1}}{2}\right),$ (6)

where $\beta_{1}$ is an arbitrary parameter (spectral parameter) and the
auxiliary fermionic superfield $f_{0,1}$ (i.e. $f_{0,1}^{2}=0$)

$\displaystyle
f_{0,1}=f_{1}^{(0,1)}+\theta_{1}b_{1}^{(0,1)}+\theta_{2}b_{2}^{(0,1)}+\theta_{1}\theta_{2}f_{2}^{(0,1)},$
(7)

satisfy

$\displaystyle
D_{x}f_{0,1}=\frac{2i}{\beta_{1}}\sinh\left(\frac{\Phi_{0}+\Phi_{1}}{2}\right),\qquad
D_{t}f_{0,1}=\beta_{1}\sinh\left(\frac{\Phi_{0}-\Phi_{1}}{2}\right).$ (8)

Consider now two successive Backlund transformations. The first one involving
superfields $\Phi_{0}$ and $\Phi_{1}$ and the parameter $\beta_{1}$ whilst the
second involves $\Phi_{1}$ and $\Phi_{3}$ with $\beta_{2}$. The Permutability
theorem states that the order in which such Backlund transformations are
employed is irrelevant, i.e. we might as well consider the first involving
$\Phi_{0}$ and $\Phi_{2}$ with $\beta_{1}$ followed by a second, involving
$\Phi_{2}$ and $\Phi_{3}$ with $\beta_{2}$. Similar to (5) we have,

$\displaystyle D_{x}(\Phi_{0}-\Phi_{1})$ $\displaystyle=$
$\displaystyle-\frac{4i}{\beta_{1}}f_{0,1}\cosh\left(\frac{\Phi_{0}+\Phi_{1}}{2}\right),$
(9) $\displaystyle D_{x}(\Phi_{1}-\Phi_{3})$ $\displaystyle=$
$\displaystyle-\frac{4i}{\beta_{2}}f_{1,3}\cosh\left(\frac{\Phi_{1}+\Phi_{3}}{2}\right),$
(10) $\displaystyle D_{x}(\Phi_{0}-\Phi_{2})$ $\displaystyle=$
$\displaystyle-\frac{4i}{\beta_{2}}f_{0,2}\cosh\left(\frac{\Phi_{0}+\Phi_{2}}{2}\right),$
(11) $\displaystyle D_{x}(\Phi_{2}-\Phi_{3})$ $\displaystyle=$
$\displaystyle-\frac{4i}{\beta_{1}}f_{2,3}\cosh\left(\frac{\Phi_{2}+\Phi_{3}}{2}\right).$
(12)

and from (8),

$\displaystyle D_{x}f_{0,1}$ $\displaystyle=$
$\displaystyle\frac{2i}{\beta_{1}}\sinh\left(\frac{\Phi_{0}+\Phi_{1}}{2}\right),$
(13) $\displaystyle D_{x}f_{1,3}$ $\displaystyle=$
$\displaystyle\frac{2i}{\beta_{2}}\sinh\left(\frac{\Phi_{1}+\Phi_{3}}{2}\right),$
(14) $\displaystyle D_{x}f_{0,2}$ $\displaystyle=$
$\displaystyle\frac{2i}{\beta_{2}}\sinh\left(\frac{\Phi_{0}+\Phi_{2}}{2}\right),$
(15) $\displaystyle D_{x}f_{2,3}$ $\displaystyle=$
$\displaystyle\frac{2i}{\beta_{1}}\sinh\left(\frac{\Phi_{2}+\Phi_{3}}{2}\right).$
(16)

The equality of the sum of equations (9) and (10) with the sum of (11) and
(12) yields the following relation

$\displaystyle\frac{1}{\beta_{1}}f_{0,1}\cosh\left(\frac{\Phi_{0}+\Phi_{1}}{2}\right)+\frac{1}{\beta_{2}}f_{1,3}\cosh\left(\frac{\Phi_{1}+\Phi_{3}}{2}\right)$
$\displaystyle=\frac{1}{\beta_{2}}f_{0,2}\cosh\left(\frac{\Phi_{0}+\Phi_{2}}{2}\right)+\frac{1}{\beta_{1}}f_{2,3}\cosh\left(\frac{\Phi_{2}+\Phi_{3}}{2}\right).$
(17)

Analogously from (6) we find

$\displaystyle D_{t}(\Phi_{0}+\Phi_{1})$ $\displaystyle=$ $\displaystyle
2\beta_{1}f_{0,1}\cosh\left(\frac{\Phi_{0}-\Phi_{1}}{2}\right),$ (18)
$\displaystyle D_{t}(\Phi_{1}+\Phi_{3})$ $\displaystyle=$ $\displaystyle
2\beta_{2}f_{1,3}\cosh\left(\frac{\Phi_{1}-\Phi_{3}}{2}\right),$ (19)
$\displaystyle D_{t}(\Phi_{0}+\Phi_{2})$ $\displaystyle=$ $\displaystyle
2\beta_{2}f_{0,2}\cosh\left(\frac{\Phi_{0}-\Phi_{2}}{2}\right),$ (20)
$\displaystyle D_{t}(\Phi_{2}+\Phi_{3})$ $\displaystyle=$ $\displaystyle
2\beta_{1}f_{2,3}\cosh\left(\frac{\Phi_{2}-\Phi_{3}}{2}\right).$ (21)

Equating now the difference of the first two, (18) and (19) and the last two
equations, (20) and (21) we get

$\displaystyle\beta_{1}f_{0,1}\cosh\left(\frac{\Phi_{0}-\Phi_{1}}{2}\right)-\beta_{2}f_{1,3}\cosh\left(\frac{\Phi_{1}-\Phi_{3}}{2}\right)$
$\displaystyle=\beta_{2}f_{0,2}\cosh\left(\frac{\Phi_{0}-\Phi_{2}}{2}\right)-\beta_{1}f_{2,3}\cosh\left(\frac{\Phi_{2}-\Phi_{3}}{2}\right).$
(22)

Solving (17) and (22), for $f_{1,3}$ and $f_{2,3}$, we get

$\displaystyle f_{1,3}$ $\displaystyle=$
$\displaystyle\Lambda_{1,3}^{(1)}f_{0,1}+\Lambda_{1,3}^{(2)}f_{0,2},$ (23)
$\displaystyle f_{2,3}$ $\displaystyle=$
$\displaystyle\Lambda_{2,3}^{(1)}f_{0,1}+\Lambda_{2,3}^{(2)}f_{0,2},$ (24)

where the coefficients $\Lambda$ are given as

$\displaystyle\Lambda_{1,3}^{(1)}=-\beta_{1}\beta_{2}\frac{\left[\cosh\left(\frac{\Phi_{0}+\Phi_{1}}{2}\right)\cosh\left(\frac{\Phi_{2}-\Phi_{3}}{2}\right)+\cosh\left(\frac{\Phi_{0}-\Phi_{1}}{2}\right)\cosh\left(\frac{\Phi_{2}+\Phi_{3}}{2}\right)\right]}{\left[\cosh\left(\frac{\Phi_{2}-\Phi_{3}}{2}\right)\cosh\left(\frac{\Phi_{1}+\Phi_{3}}{2}\right)\beta_{1}^{2}-\cosh\left(\frac{\Phi_{1}-\Phi_{3}}{2}\right)\cosh\left(\frac{\Phi_{2}+\Phi_{3}}{2}\right)\beta_{2}^{2}\right]},$
$\displaystyle\Lambda_{1,3}^{(2)}=\frac{\left[\cosh\left(\frac{\Phi_{0}+\Phi_{2}}{2}\right)\cosh\left(\frac{\Phi_{2}-\Phi_{3}}{2}\right)\beta_{1}^{2}+\cosh\left(\frac{\Phi_{0}-\Phi_{2}}{2}\right)\cosh\left(\frac{\Phi_{2}+\Phi_{3}}{2}\right)\beta_{2}^{2}\right]}{\left[\cosh\left(\frac{\Phi_{2}-\Phi_{3}}{2}\right)\cosh\left(\frac{\Phi_{1}+\Phi_{3}}{2}\right)\beta_{1}^{2}-\cosh\left(\frac{\Phi_{1}-\Phi_{3}}{2}\right)\cosh\left(\frac{\Phi_{2}+\Phi_{3}}{2}\right)\beta_{2}^{2}\right]},$
$\displaystyle\Lambda_{2,3}^{(1)}=-\frac{\left[\cosh\left(\frac{\Phi_{0}-\Phi_{1}}{2}\right)\cosh\left(\frac{\Phi_{1}+\Phi_{3}}{2}\right)\beta_{1}^{2}+\cosh\left(\frac{\Phi_{0}+\Phi_{1}}{2}\right)\cosh\left(\frac{\Phi_{1}-\Phi_{3}}{2}\right)\beta_{2}^{2}\right]}{\left[\cosh\left(\frac{\Phi_{2}-\Phi_{3}}{2}\right)\cosh\left(\frac{\Phi_{1}+\Phi_{3}}{2}\right)\beta_{1}^{2}-\cosh\left(\frac{\Phi_{1}-\Phi_{3}}{2}\right)\cosh\left(\frac{\Phi_{2}+\Phi_{3}}{2}\right)\beta_{2}^{2}\right]},$
$\displaystyle\Lambda_{2,3}^{(2)}=\beta_{1}\beta_{2}\frac{\left[\cosh\left(\frac{\Phi_{0}+\Phi_{2}}{2}\right)\cosh\left(\frac{\Phi_{1}-\Phi_{3}}{2}\right)+\cosh\left(\frac{\Phi_{0}-\Phi_{2}}{2}\right)\cosh\left(\frac{\Phi_{1}+\Phi_{3}}{2}\right)\right]}{\left[\cosh\left(\frac{\Phi_{2}-\Phi_{3}}{2}\right)\cosh\left(\frac{\Phi_{1}+\Phi_{3}}{2}\right)\beta_{1}^{2}-\cosh\left(\frac{\Phi_{1}-\Phi_{3}}{2}\right)\cosh\left(\frac{\Phi_{2}+\Phi_{3}}{2}\right)\beta_{2}^{2}\right]},$
(25)

Acting with $D_{x}$ in eqn. (9)-(12) and using (13)- (16), we find

$\displaystyle\partial_{x}(\Phi_{0}-\Phi_{1})$ $\displaystyle=$
$\displaystyle\frac{4}{\beta_{1}^{2}}\sinh(\Phi_{0}+\Phi_{1})+\frac{4i}{\beta_{1}}f_{0,1}D_{x}\left[\cosh\left(\frac{\Phi_{0}+\Phi_{1}}{2}\right)\right].$
(26) $\displaystyle\partial_{x}(\Phi_{1}-\Phi_{3})$ $\displaystyle=$
$\displaystyle\frac{4}{\beta_{2}^{2}}\sinh(\Phi_{1}+\Phi_{3})+\frac{4i}{\beta_{2}}f_{1,3}D_{x}\left[\cosh\left(\frac{\Phi_{1}+\Phi_{3}}{2}\right)\right].$
(27) $\displaystyle\partial_{x}(\Phi_{0}-\Phi_{2})$ $\displaystyle=$
$\displaystyle\frac{4}{\beta_{2}^{2}}\sinh(\Phi_{0}+\Phi_{2})+\frac{4i}{\beta_{2}}f_{0,2}D_{x}\left[\cosh\left(\frac{\Phi_{0}+\Phi_{2}}{2}\right)\right].$
(28) $\displaystyle\partial_{x}(\Phi_{2}-\Phi_{3})$ $\displaystyle=$
$\displaystyle\frac{4}{\beta_{1}^{2}}\sinh(\Phi_{2}+\Phi_{3})+\frac{4i}{\beta_{1}}f_{2,3}D_{x}\left[\cosh\left(\frac{\Phi_{2}+\Phi_{3}}{2}\right)\right].$
(29)

Notice that equations (26)-(29) correspond formally to the pure bosonic case
when the terms proportional to the fermionic superfields $f_{i,j}$ are
neglected. Equating the R.H.S. of the sum of eqns. (26) with (27) and (28)
with (29) we find

$\displaystyle\frac{1}{\beta_{1}^{2}}\left(\sinh(\Phi_{0}+\Phi_{1})-\sinh(\Phi_{2}+\Phi_{3})\right)+\frac{1}{\beta_{2}^{2}}\left(\sinh(\Phi_{1}+\Phi_{3})-\sinh(\Phi_{0}+\Phi_{2})\right)$
(30) $\displaystyle=$
$\displaystyle-\frac{i}{\beta_{1}}f_{0,1}D_{x}\left[\cosh\left(\frac{\Phi_{0}+\Phi_{1}}{2}\right)\right]-\frac{i}{\beta_{2}}f_{1,3}D_{x}\left[\cosh\left(\frac{\Phi_{1}+\Phi_{3}}{2}\right)\right]$
$\displaystyle+$
$\displaystyle\frac{i}{\beta_{2}}f_{0,2}D_{x}\left[\cosh\left(\frac{\Phi_{0}+\Phi_{2}}{2}\right)\right]+\frac{i}{\beta_{1}}f_{2,3}D_{x}\left[\cosh\left(\frac{\Phi_{2}+\Phi_{3}}{2}\right)\right].$

Factorizing the L.H.S. of (30)

$\displaystyle\frac{1}{\beta_{1}^{2}}\left(\sinh(\Phi_{0}+\Phi_{1})-\sinh(\Phi_{2}+\Phi_{3})\right)+\frac{1}{\beta_{2}^{2}}\left(\sinh(\Phi_{1}+\Phi_{3})-\sinh(\Phi_{0}+\Phi_{2})\right)$
$\displaystyle=$ $\displaystyle
2\cosh({{\Phi_{0}+\Phi_{1}+\Phi_{2}+\Phi_{3}}\over{2}})\left({{1}\over{\beta_{1}^{2}}}\sinh({{\Phi_{0}+\Phi_{1}-\Phi_{2}-\Phi_{3}}\over{2}})+{{1}\over{\beta_{2}^{2}}}\sinh({{-\Phi_{0}+\Phi_{1}-\Phi_{2}+\Phi_{3}}\over{2}})\right)$
$\displaystyle=2\cosh({{\Phi_{0}+\Phi_{1}+\Phi_{2}+\Phi_{3}}\over{2}})\left(\sinh({{\Phi_{0}-\Phi_{3}}\over{2}})\cosh({{\Phi_{1}-\Phi_{2}}\over{2}}){{(\beta_{1}^{2}-\beta_{2}^{2})}\over{\beta_{1}^{2}\beta_{2}^{2}}}\right.$
$\displaystyle\left.+\sinh({{\Phi_{1}-\Phi_{2}}\over{2}})\cosh({{\Phi_{0}-\Phi_{3}}\over{2}})({{{\beta_{1}^{2}+\beta_{2}^{2}}\over{\beta_{1}^{2}\beta_{2}^{2}}}})\right).$
(31)

The vanishing of this expression leads to

$\displaystyle\Phi_{3}-\Phi_{0}=2arctanh\left(\delta\;\tanh({{\Phi_{1}-\Phi_{2}}\over{2}})\right)\equiv\Gamma(\Phi_{1}-\Phi_{2}),\qquad\delta=-\left({{{\beta_{1}^{2}+\beta_{2}^{2}}\over{\beta_{1}^{2}-\beta_{2}^{2}}}}\right).$
(32)

For the more general case taking into account the fermionic superfields
$f_{i,j}$ we propose the following ansatz,

$\displaystyle\Phi_{3}=\Phi_{0}+\Gamma(\Phi_{1}-\Phi_{2})+\Delta,$ (33)

where $\Delta$ is a bosonic superfield proportional to the product
$f_{0,1}f_{0,2}$, i.e.,

$\displaystyle\Delta=\lambda
f_{0,1}f_{0,2},\qquad\lambda=\lambda(\Phi_{1}-\Phi_{2}).$ (34)

Due to the fact that $\Delta^{2}=0$, eqns. (25) take the general form

$\displaystyle\Lambda_{1,3}^{(1)}$ $\displaystyle=$
$\displaystyle-a+c_{1}f_{0,1}f_{0,2},$ $\displaystyle\Lambda_{1,3}^{(2)}$
$\displaystyle=$ $\displaystyle-b+c_{2}f_{0,1}f_{0,2},$
$\displaystyle\Lambda_{2,3}^{(1)}$ $\displaystyle=$ $\displaystyle
b+c_{3}f_{0,1}f_{0,2},$ $\displaystyle\Lambda_{2,3}^{(2)}$ $\displaystyle=$
$\displaystyle a+c_{4}f_{0,1}f_{0,2},$ (35)

where $c_{i}=c_{i}(\Phi_{0},\Phi_{1},\Phi_{2}),i=1,\cdots 4$ do not contribute
to eqns (23) and (24). Substituting (33) and (34) in (25) we obtain

$\displaystyle
a=\frac{\delta_{1}}{\sqrt{1-\delta^{2}\tanh^{2}\left(\frac{\Phi_{1}-\Phi_{2}}{2}\right)}},\qquad
b=\frac{\delta\,\mathrm{sech}\left(\frac{\Phi_{1}-\Phi_{2}}{2}\right)}{\sqrt{1-\delta^{2}\tanh^{2}\left(\frac{\Phi_{1}-\Phi_{2}}{2}\right)}}$
(36)

where $\delta_{1}=\frac{2\beta_{1}\beta_{2}}{(\beta_{1}^{2}-\beta_{2}^{2})}$.
Inserting (35) in (23) and (24) we find, since $f_{01}^{2}=f_{02}^{2}=0$,

$\displaystyle f_{1,3}$ $\displaystyle=$ $\displaystyle-af_{0,1}-bf_{0,2},$
$\displaystyle f_{2,3}$ $\displaystyle=$ $\displaystyle bf_{0,1}+af_{0,2}.$
(37)

From the fact that $\delta^{2}-\delta_{1}^{2}=1$, it follows that

$\displaystyle f_{1,3}f_{2,3}=f_{0,1}f_{0,2}.$

Adding (9) and (10) we find

$\displaystyle
D_{x}(\Phi_{3}-\Phi_{0})=\frac{4i}{\beta_{1}}f_{0,1}\cosh\left(\frac{\Phi_{0}+\Phi_{1}}{2}\right)+\frac{4i}{\beta_{2}}f_{1,3}\cosh\left(\frac{\Phi_{1}+\Phi_{3}}{2}\right).$
(38)

Substituting (33) and (37) in (38) we find

$\displaystyle
f_{0,1}\Sigma_{1}+f_{0,2}\Sigma_{2}+(D_{x}\lambda)f_{0,1}f_{0,2}=0,$ (39)

where

$\displaystyle\Sigma_{1}$ $\displaystyle=$
$\displaystyle\partial_{\xi}\Gamma_{|_{\xi=(\Phi_{1}-\Phi_{2})}}\frac{4i}{\beta_{1}}\cosh\left(\frac{\Phi_{0}+\Phi_{1}}{2}\right)-\lambda\frac{2i}{\beta_{2}}\sinh\left(\frac{\Phi_{0}+\Phi_{2}}{2}\right)$
$\displaystyle-\frac{4i}{\beta_{1}}\cosh\left(\frac{\Phi_{0}+\Phi_{1}}{2}\right)-\Lambda_{1,3}^{(1)}\frac{4i}{\beta_{1}}\cosh\left(\frac{\Phi_{0}+\Phi_{1}+\Gamma}{2}\right),$
$\displaystyle\Sigma_{2}$ $\displaystyle=$
$\displaystyle-\partial_{\xi}\Gamma_{|_{\xi=(\Phi_{1}-\Phi_{2})}}\frac{4i}{\beta_{2}}\cosh\left(\frac{\Phi_{0}+\Phi_{2}}{2}\right)+\lambda\frac{2i}{\beta_{1}}\sinh\left(\frac{\Phi_{0}+\Phi_{1}}{2}\right)$
$\displaystyle-\Lambda_{1,3}^{(2)}\frac{4i}{\beta_{1}}\cosh\left(\frac{\Phi_{0}+\Phi_{1}+\Gamma}{2}\right).$

The last term in (39) vanishes since

$\displaystyle
D_{x}\lambda=\partial_{\xi}\lambda_{|_{\xi=(\Phi_{1}-\Phi_{2})}}D_{x}(\Phi_{1}-\Phi_{2}),$

and $D_{x}(\Phi_{1}-\Phi_{2})$, from (9) and (11), is proportional to
$f_{0,1}$ and $f_{0,2}$. Since $f_{0,1}$ and $f_{0,2}$ are independent, (39)
yields a pair of algebraic equations for $\lambda$, i.e.
$\Sigma_{1}=\Sigma_{2}=0$ which are satisfied by

$\displaystyle\lambda=-\frac{4\sinh\left(\frac{\Phi_{1}-\Phi_{2}}{2}\right)\beta_{1}\beta_{2}(\beta_{1}^{2}+\beta_{2}^{2})}{\beta_{1}^{4}+\beta_{2}^{4}-2\cosh(\Phi_{1}-\Phi_{2})\beta_{1}^{2}\beta_{2}^{2}}.$
(40)

and therefore

$\displaystyle\Phi_{3}=\Phi_{0}+2\,\mathrm{Arctanh}\left[\left(\frac{\beta_{2}^{2}+\beta_{1}^{2}}{\beta_{2}^{2}-\beta_{1}^{2}}\right)\tanh\left(\frac{\Phi_{1}-\Phi_{2}}{2}\right)e^{\Omega
f_{0,1}f_{0,2}}\right],$ (41)

where

$\displaystyle\Omega=\delta_{1}\mathrm{sech}\left(\frac{\Phi_{1}-\Phi_{2}}{2}\right).$

In order to write eqn. (41) in components, we need to specify the superfields
$f_{0,1},f_{0,2}$ in terms of the components of $\Phi_{0},\Phi_{1},\Phi_{2}$.
These are given by eqns. (5.91), (5.94), (5.96) and (5.98) in the the appendix
of ref. [8]. Introducing $\sigma_{k}=-\frac{2}{\beta_{k}^{2}}$ $(k=1,2)$, the
solution (41) in components according to (2) becomes

$\displaystyle\phi_{3}$ $\displaystyle=$
$\displaystyle\phi_{0}+2\,\mathrm{Arctanh}\left[\delta\tanh\left(\frac{\phi_{1}-\phi_{2}}{2}\right)\right]$
(42)
$\displaystyle-\frac{\Delta_{2}}{8\sqrt{\sigma_{1}\sigma_{2}}}\left[\frac{\bar{\psi}_{0}(\bar{\psi}_{1}-\bar{\psi}_{2})+\bar{\psi}_{1}\bar{\psi}_{2}}{\cosh\left(\frac{\phi_{0}+\phi_{1}}{2}\right)\cosh\left(\frac{\phi_{0}+\phi_{2}}{2}\right)}\right],$
$\displaystyle\bar{\psi}_{3}$ $\displaystyle=$
$\displaystyle\bar{\psi}_{0}+\Delta_{1}(\bar{\psi}_{1}-\bar{\psi}_{2})$ (43)
$\displaystyle-\frac{\Delta_{2}}{2}\left[\sqrt{\frac{\sigma_{2}}{\sigma_{1}}}\frac{\sinh\left(\frac{\phi_{0}+\phi_{2}}{2}\right)}{\cosh\left(\frac{\phi_{0}+\phi_{1}}{2}\right)}(\bar{\psi}_{0}-\bar{\psi}_{1})-\sqrt{\frac{\sigma_{1}}{\sigma_{2}}}\frac{\sinh\left(\frac{\phi_{0}+\phi_{1}}{2}\right)}{\cosh\left(\frac{\phi_{0}+\phi_{2}}{2}\right)}(\bar{\psi}_{0}-\bar{\psi}_{2})\right],$
$\displaystyle\psi_{3}$ $\displaystyle=$
$\displaystyle\psi_{0}+\Delta_{1}(\psi_{1}-\psi_{2})$ (44)
$\displaystyle-\frac{\Delta_{2}}{2}\left[\sqrt{\frac{\sigma_{2}}{\sigma_{1}}}\frac{\sinh\left(\frac{\phi_{0}-\phi_{1}}{2}\right)}{\cosh\left(\frac{\phi_{0}-\phi_{2}}{2}\right)}(\psi_{0}+\psi_{2})-\sqrt{\frac{\sigma_{1}}{\sigma_{2}}}\frac{\sinh\left(\frac{\phi_{0}-\phi_{2}}{2}\right)}{\cosh\left(\frac{\phi_{0}-\phi_{1}}{2}\right)}(\psi_{0}+\psi_{1})\right],$

where

$\displaystyle\Delta_{1}$ $\displaystyle=$
$\displaystyle\frac{2}{\sinh\left(\phi_{1}-\phi_{2}\right)}\left[\frac{\delta\tanh\left(\frac{\phi_{1}-\phi_{2}}{2}\right)}{1-\delta^{2}\tanh^{2}\left(\frac{\phi_{1}-\phi_{2}}{2}\right)}\right],$
$\displaystyle\Delta_{2}$ $\displaystyle=$
$\displaystyle\frac{A\sinh\left(\frac{\phi_{1}-\phi_{2}}{2}\right)}{B-\sinh^{2}\left(\frac{\phi_{1}-\phi_{2}}{2}\right)},$
$\displaystyle\delta$ $\displaystyle=$
$\displaystyle\frac{\sigma_{1}+\sigma_{2}}{\sigma_{1}-\sigma_{2}},\qquad
A=\frac{\sigma_{1}+\sigma_{2}}{\sqrt{\sigma_{1}\sigma_{2}}},\qquad
B=\frac{(\sigma_{1}-\sigma_{2})^{2}}{4\,\sigma_{1}\sigma_{2}}.$ (45)

One soliton Solution

The Backlund equations in components for $\Phi_{0}=0$ take the form )from (7):

$\displaystyle\partial_{x}\phi_{1}$ $\displaystyle=$ $\displaystyle
2\,\sigma_{1}\sinh\phi_{1},\qquad\partial_{t}\phi_{1}=\frac{2}{\sigma_{1}}\sinh\phi_{1},$
(46) $\displaystyle\bar{\psi}_{1}$ $\displaystyle=$ $\displaystyle
2\sqrt{2\sigma_{1}}\cosh\left(\phi_{1}/2\right)f_{1}^{(0,1)},\qquad\psi_{1}=2\sqrt{{{2}\over{\sigma_{1}}}}\cosh\left(\phi_{1}/2\right)f_{1}^{(0,1)}$
(47) $\displaystyle\partial_{x}f_{1}^{(0,1)}$ $\displaystyle=$
$\displaystyle\sqrt{{{\sigma_{1}}\over{2}}}\cosh\left(\phi_{1}/2\right)\bar{\psi}_{1},\qquad\partial_{t}f_{1}^{(0,1)}={{1}\over{\sqrt{2\sigma_{1}}}}\cosh\left(\phi_{1}/2\right)\psi_{1}$
(48)

By direct integration we obtain

$\displaystyle\phi_{1}$ $\displaystyle=$
$\displaystyle\ln\left(\frac{1+E_{1}}{1-E_{1}}\right),\qquad
E_{1}=b_{1}\exp\left(2\sigma_{1}x+2\sigma_{1}^{-1}t\right),$ (49)
$\displaystyle\bar{\psi}_{1}$ $\displaystyle=$
$\displaystyle\epsilon_{1}\frac{a_{1}}{b_{1}}E_{1}\left(\frac{1}{1+E_{1}}+\frac{1}{1-E_{1}}\right),\qquad\psi_{1}=\frac{\bar{\psi}_{1}}{\sigma_{1}},$
(50)

where $a_{1}$ and $b_{1}$ are arbitrary constants and $\epsilon_{1}$ is a
fermionic parameter.

Two soliton Solution

Choosing $\Phi_{0}=0$ and $\Phi_{1},\Phi_{2}$ as one soliton solutions with
components

$\displaystyle\phi_{k}$ $\displaystyle=$
$\displaystyle\ln\left(\frac{1+E_{k}}{1-E_{k}}\right),\qquad
E_{k}=b_{k}\exp\left(2\sigma_{k}x+2\sigma_{k}^{-1}t\right),$
$\displaystyle\bar{\psi}_{k}$ $\displaystyle=$
$\displaystyle\epsilon_{k}\frac{a_{k}}{b_{k}}E_{k}\left(\frac{1}{1+E_{k}}+\frac{1}{1-E_{k}}\right),$
$\displaystyle\psi_{k}$ $\displaystyle=$
$\displaystyle\frac{\bar{\psi}_{k}}{\sigma_{k}},\qquad k=1,2$

where $a_{k}$, $b_{k}$ are arbitrary constants and $\epsilon_{k}$ fermiônic
parameters we find from (41),

$\displaystyle\phi_{3}$ $\displaystyle=$ $\displaystyle
2\,\mathrm{Arctanh}\left[\delta\tanh\left(\frac{\phi_{1}-\phi_{2}}{2}\right)\right]-\frac{\Delta_{2}}{8\sqrt{\sigma_{1}\sigma_{2}}}\left[\frac{\bar{\psi}_{1}\bar{\psi}_{2}}{\cosh\left(\frac{\phi_{1}}{2}\right)\cosh\left(\frac{\phi_{2}}{2}\right)}\right],$
(51) $\displaystyle\bar{\psi}_{3}$ $\displaystyle=$
$\displaystyle\left[\Delta_{1}+\frac{\Delta_{2}}{2}\sqrt{\frac{\sigma_{2}}{\sigma_{1}}}\frac{\sinh\left(\frac{\phi_{2}}{2}\right)}{\cosh\left(\frac{\phi_{1}}{2}\right)}\right]\bar{\psi}_{1}-\left[\Delta_{1}+\frac{\Delta_{2}}{2}\sqrt{\frac{\sigma_{1}}{\sigma_{2}}}\frac{\sinh\left(\frac{\phi_{1}}{2}\right)}{\cosh\left(\frac{\phi_{2}}{2}\right)}\right]\bar{\psi}_{2},$
(52) $\displaystyle\psi_{3}$ $\displaystyle=$
$\displaystyle\left[\Delta_{1}-\frac{\Delta_{2}}{2}\sqrt{\frac{\sigma_{1}}{\sigma_{2}}}\frac{\sinh\left(\frac{\phi_{2}}{2}\right)}{\cosh\left(\frac{\phi_{1}}{2}\right)}\right]\psi_{1}-\left[\Delta_{1}-\frac{\Delta_{2}}{2}\sqrt{\frac{\sigma_{2}}{\sigma_{1}}}\frac{\sinh\left(\frac{\phi_{1}}{2}\right)}{\cosh\left(\frac{\phi_{2}}{2}\right)}\right]\psi_{2},$
(53)

By rescaling of parameters

$\displaystyle\sigma_{k}\to\gamma_{k},\qquad\epsilon_{k}\to c_{k}\qquad k=1,2$
$\displaystyle
b_{1}\to\frac{b_{1}}{2}\left(\frac{\gamma_{1}-\gamma_{2}}{\gamma_{1}+\gamma_{2}}\right),\qquad
b_{2}\to-\frac{b_{2}}{2}\left(\frac{\gamma_{1}-\gamma_{2}}{\gamma_{1}+\gamma_{2}}\right),$
$\displaystyle
a_{1}\to-\gamma_{1}\left(\frac{\gamma_{1}-\gamma_{2}}{\gamma_{1}+\gamma_{2}}\right),\qquad
a_{2}\to\gamma_{2}\left(\frac{\gamma_{1}-\gamma_{2}}{\gamma_{1}+\gamma_{2}}\right).$

we verify that solution (51) and (52) coincide precisely with solution
(3.28)-(3.29) of ref. [5] (after choosing $\gamma_{3}=-\gamma_{1}$ and
$\gamma_{4}=-\gamma_{2}$).

Acknowledgements
LHY acknowledges support from Fapesp, JFG and AHZ thank CNPq for partial
support.

## References

  * [1] C. Rogers, in “Soliton Theory:a survey of results”, Ed. A. Fordy, Manchester Univ. Press (1990)
  * [2] H. Wahlquist and F. Estabrook, Phys. Rev. Lett. 31 (1973) 1386
  * [3] Q.P. Liu and Y.F. Xie Phys. Lett. 325A (2004) 139; Q.P. Liu and Xing-Biao Hu J. Physics A38 (2005) 6371
  * [4] B. Grammaticos, A. Ramani and A.S. Carstea, J. Physics A34 (2001) 4881
  * [5] J.F. Gomes, L. H. Ymai and A.H. Zimerman, Phys. Lett. 359A (2006) 630, hep-th/0607107
  * [6] Q.P. Liu and Meng-Xia Zhang, Nonlinearity  18 (2005) 1597
  * [7] M. Chaichian and P. Kulish, Phys. Lett. 78B (1978) 413
  * [8] J.F. Gomes, L. H. Ymai and A.H. Zimerman, J. Physics A39 (2006) 7471, hep-th/0601014