The simulation model assumes that the beam electron dynamics is fully relativistic. The equations of motion of the particles are thus written as
where 
is the particle momentum,
is the relativistic factor, m and q are
respectively the mass and the charge of the simulation particle
and c is the velocity of light.
In the electrostatic approximation, the electric field is expressed
in terms of an electric potential 
and 
is
a static magnetic field imposed externally. 
can, for
example, be computed from a set of coils, using the Biot-Savart equations.
Each particle is pushed using Eq. (1) which is time-discretized
by the usual second order leapfrog scheme [5, 6].
Using the axisymmetry of the geometry and cylindrical coordinates in
which the azimuthal coordinate 
is an ignorable variable, the particle
phase space coordinates are 
and the charge density can be constructed from the ensemble of 
electrons as
The variational form of the Poisson equation, using the charge density as defined by Eq. (2), is then
where the integration is performed over the 2D computational domain
, and 
is an arbitrary
weight function from the same function
space as 
. The Finite Element discretization of
the Poisson equation is then completely defined by choosing
the bilinear basis functions
, 
on a (r,z) mesh of
N grid points [7]. The discretized potential can thus be
expressed as 
,
where the coefficients 
are obtained by
solving the following linear matrix equation:
Dirichlet boundary conditions, 
, are specified on
perfectly conducting walls by modifying the matrix 
and the charges 
in a standard manner.
Natural boundary conditions (normal derivative 
)
are automatically included in the matrix formulation above. The latter
conditions are utilized for the field at the symmetry axis r=0 and
to model the open sides of the boundaries.
The electric field at the particle position is finally calculated from
Noe that, by the choice of the bilinear basis functions ,
the charge
deposited on the grid , as defined in Eq. (4b), is
computed by an linear assignment, while the electric field interpolation
on the particle is stepwise constant (linear) in r and linear (constant)
in z for 
( 
). This scheme is similar to the
energy-conserving scheme described in [5].
At every simulation time step,
the beam is continuously injected at the left end of
the computational domain, z=0, with a prescribed constant current I,
a given number of injected particles 
and a distribution
function of guiding centers 
,
where 
are respectively the momenta perpendicular and
parallel to ,
and F is uniform in the pitch angles 
.
Finite spreads in the guiding center radius 
as well as in the
velocities can be prescribed. The charge of these
particles is given by 
while their mass m is
such that the charge to mass ratio is that of an electron
. At both ends of the domain,
free boundary conditions are applied
for particles with 
at z=0 and 
at 
.
The total number of particles is thus a time varying function with
. It increases linearly with time up to approximately
one transit of the beam electrons across the computational domain.
Electrons of a non-drifting plasma could also be
simulated together with the beam electrons. The
simulated charge of these 
plasma electrons is then
, where 
and 
is the constant plasma density. At both ends of the domain,
either periodic or reflecting boundary conditions are applied to these
plasma electrons so as to keep constant.
The equations (1) and (4) form a complete electrostatic model for the self-consistent evolution of the electron beam. It has been implemented on a massively parallel computer, using the domain decomposition technique to parallelize the particles as well as the electrostatic field. The parallel direct matrix solver described in [8] is utilized for solving the discretized Poisson equation (4).
The extension of this model to include electromagnetic effects consists
of solving the Maxwell equations for the fields 
and 
. In the axisymmetric case considered here, the
electromagnetic fields can be decomposed into the decoupled TM polarization
and the TE polarization 
.
The particle equations (1) remain unchanged but the
current density has to be calculated:
The detailed description of this electromagnetic model and its implementation is given in [9].
