2394 IEEE TRANSACTIONS ON PLASMA SCIENCE, VOL. 38, NO. 9, SEPTEMBER 2010
Time-Step Considerations in Particle Simulation
Algorithms for Coulomb Collisions in Plasmas
Bruce I. Cohen, Andris M. Dimits, Alex Friedman, and Russel E. Caflisch
Abstract—The accuracy of first-order Euler and higher-order
time-integration algorithms for grid-based Langevin equations
collision models in a specific relaxation test problem is assessed.
We show that statistical noise errors can overshadow time-step
errors and argue that statistical noise errors can be conflated
with time-step effects. Using a higher-order integration scheme
may not achieve any benefit in accuracy for examples of prac-
tical interest. We also investigate the collisional relaxation of an
initial electron-ion relative drift and the collisional relaxation to a
resistive steady-state in which a quasi-steady current is driven by
a constant applied electric field, as functions of the time step used
to resolve the collision processes using binary and grid-based, test-
particle Langevin equations models. We compare results from two
grid-based Langevin equations collision algorithms to results from
a binary collision algorithm for modeling electron-ion collisions.
Some guidance is provided on how large a time step can be used
compared to the inverse of the characteristic collision frequency
for specific relaxation processes.
Index Terms—Algorithms, collision processes, computer appli-
cations, numerical analysis, particle collisions, plasmas.
I. INTRODUCTION
T
HERE are two popular types of algorithms for Coulomb
collisions in particle simulations of plasmas using finite-
sized particles and deposition of charge and current densities
onto a grid (particle-in-cell simulation, i.e., PIC simulation).
In the binary algorithm, particles in a subdomain, e.g., a cell,
are grouped into discrete pairs of interacting particles such
that the relative velocity is scattered through an angle whose
statistical variance is dictated by the theory of Coulomb colli-
sions in a plasma in the Fokker–Planck limit [1], [2]. The post-
collision velocities of the interacting pair conserve momentum
and energy relative to the pre-collision velocities. In the second
type of algorithm, the collisions are modeled by defining test
and field particles. The test-particle velocity is subject to drag
and diffusion in three velocity dimensions using Langevin
equations whose drag and diffusion coefficients depend jointly
Manuscript received November 13, 2009; revised March 29, 2010; accepted
April 14, 2010. Date of publication June 1, 2010; date of current version
September 10, 2010. This work was performed under the auspices of the U.S.
Department of Energy by the Lawrence Livermore National Laboratory under
Contract DE-AC52-07NA27344 and under Grant DE-FG02-05ER-25710 at
University of California at Los Angesles under the Multiscale Initiative pro-
gram supported by the DOE Office of Scientific Computing Research.
B. I. Cohen, A. M. Dimits, and A. Friedman are with the Lawrence
Livermore National Laboratory, Livermore, CA 94551 USA (e-mail: cohen1@
llnl.gov; dimits1@llnl.gov; afriedman@lbl.gov).
R. E. Caflisch is with the Mathematics Department, University of California
at Los Angeles, Los Angeles, CA 90024 USA (e-mail: caflisch@math.ucla.edu).
Color versions of one or more of the figures in this paper are available online
at http://ieeexplore.ieee.org.
Digital Object Identifier 10.1109/TPS.2010.2049589
on the velocity of the test particle and the moments of the field-
particle velocity distribution deposited on the configuration-
space mesh [3]–[7]. The grid-based Langevin equations model
conserves the particle number trivially and conserves energy
and momentum approximately in a statistical sense after aver-
aging over many collisions and over the velocity distribution
functions, although energy and momentum conservation can
be repaired by scaling and shifting velocities after the Monte
Carlo collisions occur on each time step. The drag and diffusion
coefficients are derived from the classical theory of screened
Coulomb collisions in the Fokker–Planck limit [4], [8]–[10].
It is of practical interest to assess the accuracy of the time
integration of the collisional evolution of the plasma velocity
distribution using these algorithms.
We investigate the accuracy issues for first-order Euler and
higher-order time-integration algorithms for two grid-based
Langevin equations collision models for a specific relaxation
test problem. In an example of practical interest using nu-
merical parameters that are typical for plasma simulations, we
find that statistical noise errors can dominate systematic time-
step errors. We then argue that statistical noise errors can be
conflated with time-step effects. Moreover, we find that using
a higher-order integration scheme may achieve no benefit in
accuracy. We also find that when a higher-order Milstein cor-
rection [11], [12] to the Langevin equations model is included,
there is no significant change in the results of the collisional
relaxation process for the specific example considered. Results
using the binary collision algorithm for the same collisional
relaxation test problem have been reported by Wang et al. [13].
In the Wang et al. study, it was found that the mixing of
statistical errors with time-step effects made it difficult to obtain
unambiguous and clear scalings of the errors with respect to the
time step used. Here, we do not undertake a detailed conver-
gence analysis of the grid-based Langevin-equation collision
algorithm with respect to the time step and the particle number.
Instead, we limit ourselves to the examination of the influence
of changing the time step on the results of simulations using
two variants of the grid-based Langevin equations collision
algorithm for numbers of particles that are typical of those used
in well-resolved plasma physics studies.
We also investigate the collisional relaxation of an initial
electron-ion relative drift and the relaxation to a resistive
steady-state in which a quasi-steady current is driven by a
constant applied electric field, as functions of the time step
used to resolve the collision processes. We show that one of
the two grid-based Langevin equations models investigated has
an unfavorable mass-ratio scaling such that modeling electron-
ion collisions can require a much smaller time step than that
0093-3813/$26.00 © 2010 IEEE
COHEN et al.: PARTICLE SIMULATION ALGORITHMS FOR COULOMB COLLISIONS 2395
required using either the binary collision algorithm or a second
Langevin equations algorithm that employs spherical polar
velocity coordinates [6] in the example studied.
This paper is organized as follows. In Section II we pro-
vide brief overviews of the binary and test-particle, Langevin-
equations Coulomb collision algorithms used in particle codes.
We also provide some pertinent discussions of their properties.
In Section III, we present results from simulations using the
test-particle Langevin equations and the grid-based Coulomb
collision algorithms in which we have employed either a first-
order Euler integration, a higher-order predictor-corrector time
integration, or a first-order Euler integration including the
Milstein correction. In Section IV, we show simulation results
using the binary collision algorithm to study the collisional re-
laxation of a relative drift between electron and ion species. We
also present results for the collisional relaxation to a resistive
steady state given a constant electric field using both binary and
test-particle Langevin-equations Coulomb collision algorithms.
A brief summary is presented in Section V. The findings here
are intended to be of practical value to computational plasma
physicists undertaking kinetic simulations in which Coulomb
collisions are included. The findings here also intend to give
insight into the behavior of the collision algorithms consi-
dered here.
II. C
OLLISION ALGORITHMS IN PARTICLE CODES
The use of grid-based Langevin equations to model Coulomb
collisions is well established [3]–[7], [10]. The algorithms
are based on the classic theory describing screened Coulomb
collisions in the Fokker–Planck limit [8], [9], [14], which
yields the ensemble-averaged drag and diffusion coefficients
Δv/Δt and ΔvΔv/Δt. The algorithm for a velocity-
dependent Langevin equations collision operator in an isotropic
Maxwellian background plasma can be represented in a con-
venient approximation [7] as follows for a test particle with
velocity v in the local mean drift frame of the background field
particles:
Δv
z
= F
d
Δt + g(v)Δt
1/2
N
1
g(v)=
A
D
G
v
√
2v
th,f
/v
1/2
F
d
= −(m
f
/2T
f
)A
D
(1 + m
t
/m
f
)G
A
D
=8πn
f
q
2
t
q
2
f
ln Λ/m
2
t
G(u)/u ≈
1
2
1
u
3
+
3
√
π
4
u = v
t
/
√
2v
th,f
,v
th,f
=(T
f
/m
f
)
1/2
Δv
⊥1,2
=[A
D
(Φ − G)/v]
1/2
Δt
1/2
N
2,3
, Φ=erf(u)
N
1,2,3
≡Gaussian random nos.
N
i
=0,
N
2
i
=1 (1)
where ln Λ is the Coulomb logarithm, erf is the error function,
and the subscripts t and f are the test and field particles,
respectively. The z-axis is aligned with the velocity of the
test particle before the collision. Equation (1) describes the
velocity increments Δv
z
and Δv
⊥1,2
(in the z-direction and the
two binormal directions) acquired by the test particles due to
collisions with the field particles in the drift frame of the field
particles. The z-axis in (1) coincides with the velocity vector
of the test particle in the local mean drift frame of the field
particles before the test particle is collisionally scattered. The
scattered velocity vector is transformed back to the laboratory
Cartesian frame with the rotation matrix given in [1] or [6], and
the local mean drift of the field particles is added. Equation (1)
is a discretized solution of the Fokker–Planck equation for the
probability density of velocities f(v)
∂
∂t
f(v)=−
∂
∂v
· [F
d
(v)f(v)] +
1
2
∂
2
∂v∂v
· [D(v)f(v)]
F
d
= Δv/Δt, D = Δv Δv/Δt. (2)
In (1), we have corrected a minor typographical error that
appeared in the approximation to G, which was introduced by
Sherlock in (16)–(18) in his publication. [7] In the general cir-
cumstance in which the background field particle’s velocity dis-
tribution function is not an isotropic Maxwellian, Rosenbluth
potentials must be constructed; consequently, the collision op-
erator acquires more structure. [4], [14] In the simulations, the
field particles are composed of all of the particles of a specific
species, with the density, mean drift, and temperature moments
of the particle velocity distribution deposited on the configura-
tion space grid using an interpolation scheme (linear interpo-
lation in our simulations). The collision operator (1) conserves
total momentum and energy approximately if we average over
all of the particles and over an ideal distribution of random
numbers for all species present. In (19) of Manheimer et al.,
[4] a finite Δt correction to the drag F
d
is introduced to improve
energy conservation. Energy and momentum conservation can
be repaired by scaling and shifting velocities after the Monte
Carlo collisions occur on each time step [4], [6].
Although no magnetic field effects are included in the for-
mulation of the Fokker–Planck collisions presented here, the
collision formulation is applicable to magnetized plasmas in the
following sense. In the Fokker–Planck limit [8], [9], [14], many
infinitesimal small-angle collisions are assumed to occur with
individual collision events whose duration τ is arbitrarily short
so that the product of the acceleration due to electromagnetic
Newton–Lorentz forces and the time duration τ is arbitrarily
small compared to particle velocity. In this limit, the collision is
unaffected by any electromagnetic fields present, although the
overall particle trajectory is influenced by the electromagnetic
fields. From a computational perspective, this invites the use
of operator splitting to accommodate both the collisions and
the Newton–Lorentz forces due to electromagnetic fields. The
result of the collision should be insensitive to when during the
time step the collision event occurs.
A variant of the grid-based Langevin-equations Coulomb
collision operator has been introduced by Lemons et al. [6].
The methodology in [6] scatters the velocity vector of the
test particle using spherical polar coordinates (v, θ,φ).The
magnitude of the velocity is subject to both drag and diffusion,
and the polar angle is subject to diffusion. There is an advantage
in using spherical coordinates for scattering processes that are
2396 IEEE TRANSACTIONS ON PLASMA SCIENCE, VOL. 38, NO. 9, SEPTEMBER 2010
dominantly pitch-angle scattering at nearly constant energy of
the test particle, e.g., scattering of a light particle by a very
heavy particle. We will illustrate this in the following discussion
and example calculations. The Langevin equations in [6] are
repeated here for convenience
Δv = F
v
Δt + g(v)Δt
1/2
N
1
,
g(v) =
A
D
G
v
√
2v
th,f
/v
1/2
F
v
= −(A
D
/2v
2
)
×
(1 + m
t
/m
f
)v
2
/v
2
th,f
+1
G − Φ
,
A
D
=8πn
f
q
2
t
q
2
f
ln Λ/m
2
t
G(u)/u ≈
1
2
1
u
3
+
3
√
π
4
, Φ=erf(u),
u =v
t
/
√
2v
th,f
, v
th,f
=(T
f
/m
f
)
1/2
Δθ =
A
D
(Φ − G)/v
3
1/2
Δt
1/2
N
2
Δφ =2πU[0, 1]
N
1,2
≡Gaussian random nos.,
<N
i
> =0,<N
2
i
> =1,
U ≡uniform random no., 0 ≤ U ≤ 1. (3)
At low test-particle energies, F
v
and Δθ diverge in (3). To
resolve the divergence, Lemons et al. [6] retain the dominant
terms in the expression for the change in the speed; thus,
they derive v(t +Δt)=
v(t)
2
+(4/3
√
2π)A
D
Δt/v
th,f
for
v
2
≤ (4/3
√
2π)A
D
Δt/v
th,f
. However, Δθ
2
∝ 1/v
2
, Δv
2
⊥
∝
v
2
Δθ
2
is not divergent.
An algorithm for binary Coulomb collisions has been in-
troduced in the classic work of Takiziuka and Abe [1]. In
the binary algorithm, equally weighted particles in a cell are
paired and then the relative velocity vector of the two particles
is scattered through a random scattering angle with variance
dictated by the theory of screened Coulomb collisions in a
plasma [8], [9]. After the relative velocity vector is scattered,
the two scattered particle velocities are reconstructed such that
particle momentum and energy are conserved algebraically.
Particle number is conserved identically. Nanbu [2] extended
the algorithm of Takizuka and Abe [1] to allow for a larger time
step by aggregating multiple collisions. In Takizuka and Abe
[1], the relative velocity of a pair is scattered through an angle
Θ with variance related to
δ
2
=2πq
2
1
q
2
2
n
(1,2)<
ln ΛΔt/m
2
12
u
3
u = |v
1
− v
2
|,m
12
= m
1
m
2
/(m
1
+ m
2
),δ≡ tan
Θ
2
(4)
and through a random angle φ about the axis of the relative
velocity before the scattering event. The post-collision velocity
vectors of the scattered pair are constructed from the scattered
relative velocity vector. There is no separation of test and field
particles in the binary scheme, and there is no assumption that
the velocity distribution is isotropic and Maxwellian. However,
there is an implicit assumption that the value of ln Λ derived for
an abitrary velocity distribution (computed locally in a spatial
cell) deviates insignificantly from the value of ln Λ derived for
a Maxwellian. This method conserves particle number, energy,
and momentum.
Both the grid-based Langevin equations and binary collision
algorithms generalize to relativistic collisions. Both algorithms
are formally accurate through O(Δt
1/2
) and both produce drag
Δv/Δt and diffusion coefficients ΔvΔv/Δt that agree
with the Spitzer-Trubnikov theory through O(Δt) assuming
perfect statistics.
The accuracy of both algorithms requires that the velocity
change and the angle scattered by a test particle in one time step
must be small. In the binary collision algorithm, we can easily
deduce the scaling of the variance with respect to the charge
state and mass for electon-electron, electron-ion, and ion-ion
scattering. We note the reduced mass in (4), m
ee
∼ m
ei
/2=
m
ie
/2 m
ii
, s o that for low charge-state ions and T
i
∼ T
e
,
the time step for accurately resolving the electron-electron
binary collisions is comparable to that required for resolving
electron-ion collisions. Further, ion-ion collisions can be re-
solved with a substantially larger time step. However, the grid-
based Langevin equations algorithms have different numerical
characteristics. From (1), one can show that the magnitudes of
the drag and diffusion coefficients are monotonically increasing
as the test-particle speed v goes to zero. In the low-velocity
limit for test partcles, one can show that the perpendicular
and parallel velocity diffusion coefficients for scattering a test
particle on a field particle in (1) have the following scalings:
D
⊥,,ei
/D
∗
⊥,,ee
∼ Z
i
(Z
i
n
i
/n
e
)
× (T
e
/T
i
)
1/2
(m
i
/m
e
)
1/2
∝ (m
i
/m
e
)
1/2
D
⊥,,ii
/D
∗
⊥,,ee
∼ Z
3
i
(Z
i
n
i
/n
e
)
× (T
e
/T
i
)
1/2
(m
e
/m
i
)
3/2
∝ (m
e
/m
i
)
3/2
D
⊥,,ie
/D
∗
⊥,,ee
∼ Z
i
(Z
i
n
i
/n
e
)
× (m
e
/m
i
)
2
∝ (m
e
/m
i
)
2
(5)
where the ion charge is Z
i
e, D
⊥kl
= Δv
⊥
Δv
⊥
/Δt and
D
kl
= Δv
Δv
/Δt for species k scattering on species l.For
T
i
∼ T
e
, however, a thermal electron has a much larger velocity
than a thermal i on; thus, the drag and diffusion coefficients for
electron-ion scattering should be evaluated for v
t
v
th,f
.In
the large argument limit of (1)
D
,tf
∝
n
f
q
2
t
q
2
f
/m
2
t
v
3
v
2
th,f
D
⊥,tf
∝
n
f
q
2
t
q
2
f
/m
2
t
vv
2
th,f
v
2
th,f
F
d
∝−
n
f
q
2
t
q
2
f
/m
2
t
v
th,f
v
2
v
th,f
. (6)
The drag and diffusion coefficients go to zero at different rates
for large test-particle velocities, and D
⊥,tf
and the drag F
d
have no dependence on the mass of the field particles. How-
ever , D
,tf
∝ T
f
/m
f
, and thus the parallel velocity diffusion
COHEN et al.: PARTICLE SIMULATION ALGORITHMS FOR COULOMB COLLISIONS 2397
coefficient for electron-ion collisions is smaller than that for
electron-electron collisions by Z
i
m
e
/m
i
for T
e
∼ T
i
and de-
creases as 1/v
3
. We believe that the unfavorable mass-ratio
scaling for the ratio of the electron-ion collisional diffusion to
the electron-electron collisional diffusion in (5), Z
i
(m
i
/m
e
)
1/2
at low velocities leads to the required use of a significantly
smaller time step to accurately resolve the electron-ion colli-
sions for the grid-based Langevin equations collision algorithm
in (1), which is borne out in our simulation experience (illus-
trated in Section IV).
In the Lemons et al. algorithm for large test-particle
velocities,
D
v,tf
∝
n
f
q
2
t
q
2
f
/m
2
t
v
3
v
2
th,f
D
⊥,tf
∝
n
f
q
2
t
q
2
f
/m
2
t
vv
2
th,f
v
2
th,f
F
v
∝−(m
t
/m
f
)
n
f
q
2
t
q
2
f
/m
2
t
v
th,f
v
2
v
th,f
. (7)
Compared to the algorithm based on (1) in the large velocity
limit, the major difference is that the drag coefficient in the
Lemons et al., algorithm has F
v
/F
d
∝ (m
t
/m
f
); thus, the
drag for electron-ion collisions is much smaller compared to
the electron-electron collisions in the Lemons et al. algorithm.
Furthermore, in the Lemons et al. algorithm, the scattering of
the low-velocity test particles, where the drag dominates the
stochastic terms, is performed to remove divergences. [8] In
the limit that m
e
/m
i
→ 0, the scattering of electrons on ions
should only produce angle scattering of the electron without
changing the electron test-particle energy. The Lemons et al. [6]
Langevin equations algorithm in (3) for electron collisions on
ions is dominated by angle scattering. Moreover, the diffusion
D
v,ei
in the magnitude of the electron test-particle velocity
and the drag F
v,ei
are relatively weak. The slowing of the
electrons on the ions comes from the angle scattering. Thus, we
expect that (3) should be better able to accommodate electron-
ion scattering; and it should allow larger time steps than the
algorithm based on (1). An example of simulating electron-ion
collisions is presented in Section IV.
When the collisional scattering of the velocity vector in one
time step is too large using the grid-based Langevin equations
algorithm based on (1), relaxation rates and energy and momen-
tum conservation become inaccurate. If only the scattering an-
gle Δθ becomes too large in (3) when the time step is too large,
the angular diffusion rate becomes inaccurate, but the change in
the test-particle energy may s till be small. In contrast, when the
time step becomes too large in the binary collision algorithm,
conservation of energy and momentum is still preserved, and
although relaxation rates may not be reproduced accurately, the
binary collision algorithm fails gracefully.
In the work of Wang et al., [13] the convergence properties of
the Takizuka and Abe and the Nanbu binary collision operators
with respect to particle number and time step were s tudied. It
was found that the Nanbu collision algorithm achieved a factor-
of-two improvement in relative accuracy over the Takizuka and
Abe basic algorithm for the same time step. The underlying
properties of the Nanbu algorithm were studied analytically in
the work of Dimits et al. [15].
III. C
ORRECTIONS TO FIRST-ORDER EULER INTEGRATION
OF THE
GRID-BASED LANGEVIN EQUATIONS
COLLISION OPERATOR
In this section, we report the results of the studies of the grid-
based Langevin equations collision operator (1), attempting to
extend it to higher-order accuracy than first-order Euler. We
explore several time-discretization schemes for the grid-based
Langevin equations collision operator, which are represented
schematically in the following finite-difference Langevin equa-
tions fashioned after (1). In the cases considered in this section,
only ion-ion collisions are considered.
First-order Euler
y
n+1
= y
n
+Δt
dy
dt
(y
n
)+N
1
Δt
1/2
D(y
n
). (8)
Predictor-corrector (modified Euler)
predictor : y
n+1
= y
n
+Δt
dy
dt
(y
n
)+N
1
Δt
1/2
D(y
n
)
corrector : y
n+1
= y
n
+Δt
dy
dt
y
n+1
+ y
n
/2
+ N
2
Δt
1/2
D
y
n+1
+ y
n
/2
.
(9)
Predictor-corrector (two-step scheme)
predictor : y
n+1/2
= y
n
+
1
2
Δt
dy
dt
(y
n
)
+ N
1
1
2
Δt
1/2
D(y
n
)
corrector : y
n+1
= y
n
+Δt
dy
dt
(y
n+1/2
)
+ N
2
Δt
1/2
D(y
n+1/2
). (10)
Partial corrector
predictor : y
n+1
=y
n
+Δt
dy
dt
y
n
+y
n−1
/2
+N
1
Δt
1/2
D
y
n
+y
n−1
/2
partial corrector : y
n+1
=y
n
+Δt
dy
dt
y
n+1
+y
n
/2
+N
2
Δt
1/2
D
y
n+1
+y
n
/2
(11)
subject to the constraint relations that the ensemble averages
(y
n+1
− y
n
)/Δt = dy/dt|
y
n
and (y
n+1
− y
n
)
2
/Δt =
D(y
n
) (to lowest order in powers of Δt) are the drag and
diffusion coefficients, respectively, for the vector components
of the test-particle velocities given in (1). The partially
corrected Euler algorithm [16] is a special case in the family
of multistep Runge–Kutta methods and is closely related to the
two-stage explicit Adams method. It has advantages in that it
is easy to initialize and it can be represented with a reduced
2398 IEEE TRANSACTIONS ON PLASMA SCIENCE, VOL. 38, NO. 9, SEPTEMBER 2010
Fig. 1. BZOHAR simulations of like-species Coulomb collisions using the
grid-based Langevin equations algorithm (1) with first-order Euler time inte-
gration to study the relaxation of a weak temperature anisotropy: (a) simulated
temperature anisotropy versus time and (b) the difference of the simulated
temperature anisotropy from theory versus time for various values of the time
step (color online).
Fig. 2. BZOHAR simulations of like-species Coulomb collisions using the
grid-based Langevin equations algorithm (1) with partial corrector time inte-
gration to study the relaxation of a weak temperature anisotropy: (a) simulated
temperature anisotropy versus time and (b) the difference of the simulated
temperature anisotropy from theory versus time for various values of the time
step (color online).
Fig. 3. BZOHAR simulations of like-species Coulomb collisions using the
grid-based Langevin equations algorithm (1) with modified Euler predictor-
corrector time integration to study the relaxation of a weak temperature
anisotropy: (a) simulated temperature anisotropy versus time and (b) the
difference of the simulated temperature anisotropy from theory versus time for
various values of the time step (color online).
number of function evaluations: the corrector-step function
evaluations can be saved for reuse in the function evaluations
in the next predictor step. The method may be a natural
“upgrade” for some codes using predictor-corrector methods.
It is an upgrade in that it may be possible to nearly double the
efficiency of the calculation even with little rearrangement of
the code.
Figs. 1–3 show the results of a series of BZOHAR sim-
ulations [5] using the grid-based Langevin equations ion-ion
collision operator to study the collisional relaxation of a weak
ion temperature anisotropy (T
y
=0.95T
x
,T
x
= T
z
) in which
we vary the product of the characteristic ion-ion collision
frequency ν
∗
and the time step Δt for 1 000 particles per cell
and N
p
=2×10
6
particles in one spatial dimension, where
ν
∗
=(8/3
√
π)ν
0
, ν
0
=
√
2πq
4
n ln Λ/
√
mT
3/2
and ν
0
is the
Braginskii [20] characteristic collision frequency for like-
species collisions. In all of the BZOHAR simulations reported
here, the electrons are modeled as a fluid with a Boltzmann
response: the ions are particles; there are self-consistent electric
fields at thermal levels and particle advection; and the grid
cell size is chosen equal to the electron Debye length. The
value of ln Λ was scaled so that the value of ν
∗
Δt can be
set artificially. The exponential relaxation rate for a weak
temperature anisotropy is given by Trubnikov [9], ν
relax
=
(8/5
√
2π)ν
0
. Fig. 1 shows results for the relaxation of a
temperature anisotropy using the first-order Euler integration
scheme, (4). Figs. 2 and 3 show the corresponding results for
the first predictor-corrector algorithm (modified Euler) and the
partial corrector and (9) and (11), respectively. In each of these
figures, we plot the relaxation of the temperature anisotropy
normalized to the initial anisotropy as a function of time
and, separately, the difference of the normalized temperature
anisotropy with respect to the asymptotic theory exp(−ν
relax
t)
for several different values of ν
∗
Δt. The number of particles
per cell and the total number of particles used here are typical
of those commonly used in well-resolved particle simulations
of many plasma phenomena. In these simulations, we observe
the weakly anisotropic Maxwellian ion velocity distributions
relax by transferring energy to the colder velocity dimension
from the hotter dimensions. The lower energy ions being more
collisional tend to relax their temperature anisotropy faster than
the more energetic ions.
The object of the scan with respect to ν
∗
Δt is to assess
any trend in the convergence of the s imulation results to the
asymptotic theory. For completeness, we list several of the
sources that might contribute to the deviations of the simulation
results from the theory: i) the theory is asymptotic in the
small parameter (T
y
− T
x
)/T
x
, but a finite value of the initial
anisotropy is used; ii) the drag and diffusion coefficients used
here are calculated from an assumed isotropic Maxwellian
velocity distribution, which is only an approximation to the
actual weakly anisotropic distribution of test-particle velocities;
iii) errors are associated with the finite value of ν
∗
Δt used
in the discrete time integration; iv) a finite number of test-
particle velocities is used to resolve the velocity distribution
(and the associated random numbers used in initializing the
velocity distribution of the particles); and v) there is a deviation
from an ideal distribution of the finite set of random numbers
associated with the finite number of collisions during any
time interval. There are systematic errors associated with i),
ii), and iii), and random errors associated with iv) and v).
Moreover, the random errors in iv) and v) are i ndependent
of one another. However, the error analysis has an additional
complication. We note that in a fixed physical time interval τ ,
the number of collision events is determined by τ/Δt. Hence,
as we increase ν
∗
Δt, the number of collision events and the
number of random numbers N
r
associated with the collisions
