Full text of DOI:10.1002/1521-3986(200101)41:1<15::AID-CTPP15>3.0.CO;2-1

Contrib. Plasma Phys. 41 (2001) 1, 15−25
Effective Potentials, Energies, and Pair-distribution
Functions of Plasmas by Monte-Carlo Simulations
H. Wagenknecht, W. Ebeling, A. F o¨ rster
Institut f u¨ r Physik, Humboldt-Universit a¨ t zu Berlin,
Invalidenstraße 110, D-10115Berlin, Germany
e-mail: ebeling@physik.hu-berlin.de
Received 17 April 2000, in final form 2 October 2000
Abstract
The one-component plasma (OCP), an electron gas on a positive background, and the subsys-
tem of free charges in a two-component plasma (TCP) is investigated using a quasi-classical
approach. Quantum effects such as the Heisenberg uncertainty principle and the Pauli ex-
clusion principle are incorporated into our model by effective potentials. In case of the OCP
the symmetry effects are of special interest whereas in case of TCP the Heisenberg effects
play the significant role. The method of Slater sums is used to obtain effective potentials.
Monte Carlo simulations for both OCP and TCP were carried out to calculate the Coulombic
energy, the effective energy, and the pair-distribution functions.
1
Introduction:systems of charged fermions
The description of many-particle systems of charged fermions belongs to the most
interesting problems in modern physics. Such a description can be made either in
the physical or in the chemical picture. In the physical picture the system consists of
electrons and nuclei only and it is completely determined by the Coulomb interaction
of these point-like particles and their statistics.
In the chemical picture the description of the system is more complex. It is composed
of free electrons, free nuclei, and different types of bound states such as ions, atoms,
and molecules. All these particles are introduced as individual ”chemical species”
which interact by effective pair potentials. Although the treatment of a plasma within
the chemical picture yields a good description of the system it is essential to carry
out a careful quantum-mechanical calculation of all different types of effective pair
potentials. This way a double counting of effects can be avoided.
We investigate a plasma in the region of partial degeneration and strong coupling.
The next section gives an introduction of effective pair potentials. In section 3 we
present results for the simulations of an electron gas which is of interest in solid state
physics and which has been investigated by many researchers, e.g. [1]-[4]. In section
4
results for the mass-symmetrical plasma are given. Such systems may be observed
in semiconductors [5] or in astrophysics [6]. At the temperatures of the simulations
the plasma is almost fully ionised and hence interactions between bound states and
free charges are thermodynamically insignificant. For our investigation we restrict
ourselves to a system of free electrons (index e) and free positrons (index p).

16
Contrib. Plasma Phys. 41 (2001) 1
2
2
Effective potentials
.1 Coordinate-dependent potentials
Effective potentials can be obtained by the method of Slater sums which was already
used by Morita [7] and was further developed by Kelbg [8], Ebeling et al. [9] and
Rohde et al. [10]. All quantum effects are included into an effective pair potential by
Coul
modification of the original Coulomb energy U
to an effective energy:
eiej
eiej
→
F(r , m , T) ,
(1)
ij
ij
4
πε0rij
4πε r
0 ij
m m
i
j
where
r = | rꢀ − rꢀ |
and
m =
.
ij
i
j
ij
m + m
i
j
The effective potential of the particle i with respect to the position of the particle j
is obtained by dividing the effective energy by the charge e . The effective potential
j
depends not only on the distance between particles like the Coulomb potential but it
depends also on the temperature T and on the masses of particles m , m . At high
i
j
temperatures perturbation theory yields in first order [8]
ꢀ
ꢁ
ꢂ
ꢃ ꢄ
ꢂ
ꢃꢅ
2
r
√
rij
rij
ij
F(r , m , T) = 1 − exp −
+ π
1 − erf
,
(2)
ij
ij
2
λ
λij
λij
ij
ꢆ
where λ denotes the thermal de Broglie wavelength given by λ = ꢀ/ 2m k T, ꢀ
ij
ij
ij B
is the Planck constant and k is the Boltzmann constant. This way we find the Kelbg
B
K
interaction energy U
ij
eie
j
K
U =
F(r , m , T) .
(3)
ij
ij
ij
4
πε0rij
We assume that for the two-component plasma (TCP) of free charges the quasi-
classical Hamiltonian has the form
N
ꢇ
2
ꢇ
p
i
K
H =
+
U .
(4)
ij
2
mi
i=1
i<j
Here i, j = 1...N are the indices of the particles and N is the total number of particles.
4
We restrict our investigation to temperatures T ≥ 10 K, consider only the free charges,
and we assume that the real interactions have the same functional form as the Kelbg
potential with an additional free parameter γ . The function F(r , m , T) in Eqs.(3)
ij
ij
ij
is replaced by F(r , m , T, γ ) which is defined as
ij
ij
ij
ꢀ
ꢁ
ꢂ
ꢃ ꢄ
ꢂ
ꢃꢅ
2
r
√
rij
rij
ij
F(r , m , T, γ ) = 1 − exp −
+ π
1 − erf γ
.
(5)
ij
ij
ij
ij
2
λ
λ γ
ij ij
λij
ij

H. Wagenknecht et al., Monte Carlo Simulations
17
For γ = 1 Eq.(5) coincides with Eq. (2). Now we use the method of Slater sums to
ij
correct the height of the Kelbg potential at zero-point distance [9]
e e
i j
F(r = 0, m , T, γ ) = −k T ln [S (r = 0)]
(6)
(7)
ij
ij
ij
B
ij ij
4
πε0rij
and simultaneously to conserve the correct first derivative for r → 0
ꢄ
ꢅ
eiej
πε0
1
2
γij
e e 1
i j
−
= −
4πε λ2
0
ij
.
2
4
λ
λ γ λ
ij ij ij
ij
In Eq. (6) S denotes the binary Slater sum of the particles i and j at zero distance
ij
including also symmetry effects coming from the different spin directions. It is given
by
∞
ꢊ
√
ꢈ
ꢉ
√
2
x dx
3
ep
2
ep
−x
e
S (r = 0) =
ep ep
πξ ζ(3) + ξ ζ(5) + 4 πξ
ꢋ
ꢌ , (8)
ep
πξep
1
− exp − x
0
∞
ꢊ
√
x2
x dx
−
S (r = 0) = S (r = 0) = 2 πξ
ee
e
ꢍ
ꢎ .
(9)
ee ee
pp pp
πξee
1
− exp − x
0
Values for the logarithm of the binary Slater sums for different temperatures are given
in Tab. 1. The interaction parameter ξ is defined by ξ = −e e /(4πε λ k T). In
ij
ij
i
j
0 ij B
Eq.(8) for the bound states the Brillouin-Planck-Larkin convention is used. ζ(k) de-
notes Riemann’s Zeta-function. In Eq.(9) the symmetry effects coming from different
spin directions of electrons and positrons are already included.
From Eq. (6) the parameter γ can be calculated for r = 0. It is given by
ij
ij
√π
e e
i j
γij = −
.
(10)
λ 4πε k T ln [S (r = 0)]
ij
0
B
ij ij
Eq. (7) is identically fulfilled which means that the first derivative of the effective
potential energy Eqs. (4,5) is correct for r = 0. This identity is fulfilled only for the
ij
given argument in the error-function in Eq.(5).
The advantage of effective potentials derived from Slater sums is that the pair-
correlations and the thermodynamic functions can be described correctly. On the
other hand the momentum distributions in computer simulations with these effective
potentials are always Maxwell distributions. This may be corrected by another type
of potential - the momentum-dependent potential which will be introduced in the next
subsection.

18
Contrib. Plasma Phys. 41 (2001) 1
T/K
ln(S (r = 0))
ee
ln(S (r = 0))
ep
1
0000
0000
0000
0000
0000
0000
0000
0000
0000
-9.04071
-6.84912
-5.82477
-5.19576
-4.75786
-4.42984
-4.17197
-3.96216
-3.78702
-3.63788
7.99727
6.4018
2
3
4
5
6
7
8
9
5.51923
4.92503
4.48641
4.14441
3.86776
3.63795
3.44311
3.27526
Tab. 1: Logarithm of the binary
Slater sums for different tempera-
tures
100000
2.2 Momentum-dependent potentials
Momentum-dependent effective potentials may be obtained from the Gaussian wave-
packets (GWP) approach. In this case the quasi-classical Hamiltonian has the follow-
ing structure
ꢂ
ꢃ
ꢂ
ꢃ
N
ꢇ
2
ꢇ
ꢇ
p
rij pij
,
r0 p0
eie
j
rij
i
P
H =
+
U
+
F
,
(11)
2
mi
4πε0
r0
i=1
i<j
i<j
with p = | pꢀ − pꢀ |. The interaction between two particles in a surrounding plasma is
ij
i
j
modelled in two ways, first by the Pauli energy U and second by an effective Coulomb
P
energy that is determined by a certain function F(r /r ). Starting with a GWP
ij
0
ꢏ
ꢐ
ꢂ
ꢃ
3
/4
2
1
πσ2
( ꢀx − ꢀr )
i ꢀp . ( ꢀx − ꢀr )
+
Ψ ( ꢀx ) =
0
exp −
(12)
2
2σ
ꢀ
effective expressions of type Eq.(11) are obtained by averaging the Hamilton operator
with respect to Ψ ( ꢀx ),
0
ꢊ
∗
ꢑ
H = d ꢀx Ψ ( ꢀx )HΨ ( ꢀx ) .
(13)
0
0
An effective momentum-dependent Pauli energy can be obtained when averaging the
operator of the kinetic energy with respect to an antisymmetrised two-particle wave
function of identical particles. The Pauli energy has the form
ꢂ
ꢃ
ꢂ
ꢃ
ꢂ
ꢃ
2
2
2
∆
rij
pij
P
ij
2
ij
U = U exp −
with ∆ =
+
(14)
0
2
r0
p0
and was already suggested by Dorso et al. [11] and was used by Klakow et al. [12] and
Ebeling et al. [13]. Here p , r , and U are parameters which satisfy the conditions
0
0
0
2
2
0
p r = ꢀ and U = ꢀ /(4m r ) = c k T. m is the electron mass. The constant c can
0
0
0
e
B
e
be obtained from the calculation of the pair-correlation function of the ideal gas
ꢊ
ꢀ
ꢁ
1
H
g(r) =
d pꢀ d pꢀ exp −
(15)
1
2
(
2πmk T)3
k T
B
B

H. Wagenknecht et al., Monte Carlo Simulations
19
and fitting g(r = 0) at zero distance to the exact value g(r = 0) = 1/2 for the ideal
Fermi gas
∞
ꢂ
ꢃ
ꢊ
ꢀ
ꢂ
ꢃ ꢂ
ꢃꢁ
3
/2
2
1
2
1
8U0
2U0
exp(p ) − 0.5
2
2
gr=0
=
= √
dp p exp −p
.
(16)
2
2 π k T
k T
B
exp(p ) − 1
B
0
This condition yields the value c = 1.015819. Furthermore, by the averaging procedure
Eq.(13) an effective Coulomb energy of the form
ꢂ
ꢃ
ꢂ
ꢃ
eie
j
rij
e e
i j
rij
KTR
Uij
=
F
=
erf √
(17)
4πε0
r0
4πε r
2πλij
0
ij
is obtained which was already used by Klakow, Toepffer, and Reinhard [12] and Ebeling
et al. [13]. The KTR potential does not provide a correct description of the spatial
correlations at small distances. This deficit may be corrected by a combination of
the KTR potential with the method of Slater sums. To receive the correct potential
energy at zero-point distance we propose a modification of the argument of the error
function in Eq.(17). The modified effective Coulomb energy reads
ꢂ
ꢃ
ꢂ
ꢃ
eie
j
rij
eie
j
rij
KTR,corr
Uij
=
F
=
erf
.
(18)
(19)
∗
4πε0
r0
4πε r
R
0
ij
ij
∗
In this expression R is given by
ij
2
eiej
1
∗
R = −√
ꢈ
ꢉ
ij
P
0
ij
π 4πε k T ln S (r = 0) + U (r = 0, p˜ /p )
0
B
ij
ij
ij
0
√
P
where for the argument in U we took the most probable momentum p˜ = m k T of
the pair of particles i, j. In the binary Slater sum for zero-particle distance S (r = 0)
ij
e B
0
ij ij
we only took the Coulombic effects into account and left out the symmetry effects
for identical particles since the symmetry effects are already included in the term
P
U (r = 0, p˜ /p ). Therefore we get
ij
ij
0
∞
ꢊ
√
2
x dx
0
0
−x
e
S (r = 0) = S (r = 0) = 4 πξ
ꢍ
ꢎ .
(20)
ee ee
pp pp
ee
πξee
1
− exp − x
0
0
The expression for S (r = 0) is identical to S (r = 0) which is given in Eq.(8) since
ep ep
ep ep
there are no exchange effects between the two kinds of particles involved. In Fig. 1 we
compare the shape of different potentials for the electron-electron interaction and for
the electron-positron interaction. Here, all effective energies and the Coulomb energy
are divided by the charge e in order to obtain the potential of the i-th particle with
i
respect to the particle j. For simplicity, in the following the index i is skipped for all
K
potentials. The original Kelbg potential V shows the right first derivative but does
KTR,corr
not have the correct height at r = 0. The corrected KTR potential V
has the
ij

20
Contrib. Plasma Phys. 41 (2001) 1
12.00
6
.00
.00
0
K
VKTR,corr
VK,corr
VH
VCoul
Fig. 1: Effective Potentials at
−
6.00
T = 10000 K; Kelbg potential
K
V , corrected KTR potential
KT R, corr
V
corrected Kelbg
K, corr
potential V
, Hansen po-
V
H
tential V , and Coulomb po-
−
12.00
Coul
tential V
0
.00
1.00
2.00
r/λij
correct height but not the right shape as r → 0. Also shown is the corrected Kelbg
ij
K,corr
potential V
Eqs.(3,5,10) which we use in our TCP simulations, the potential of
H
Hansen et al. V ([14], see also work of C. Deutsch et al. [15]) and the Coulomb
Coul
potential V = e/4πε r.
0
Note that all effective potentials show the same long-range behaviour but differ for
interesting short-range distances where correlations appear.
3
Monte Carlo simulations for the electron gas
xc
Monte Carlo simulations were used to calculate the excess internal energy U of a
system consisting of N electrons as a function of Γ in the quasi-classical region (Fig. 2).
2
1/3
Γ denotes the coupling parameter and is given by Γ = e (4πn /3) /(4πε k T), n
e
0
B
e
is the electron density. The selected range of the coupling parameter Γ covers the
−2
region of weak coupling (Γ ≤ 10 ) up to the strongly coupled regime (Γ = 30).
For the classical OCP there exist analytical results for U (Γ) for Γ ≤ 0.4 [16] and
xc
very accurate MC calculations in the region of strong coupling 1 ≤ Γ ≤ 200 [4]. We
compare our results to a Pad e´ formula given by Kahlbaum [17] which covers the
whole OCP fluid range 0 ≤ Γ ≤ 200. Note that this Pad ´e formula is limited to the
pure classical plasma whereas our calculations include quantum effects through the
effective potentials.
In Figs. 3 and 4 the pair-distribution function for a temperature of 20000 K and
for different values of Γ is shown. While for Γ = 0.5 the distribution function increases
monotonously to 1, for the strong coupled region (Γ = 10) it displays already a next-
neighbour peak. The distances in the figures are given relatively the Bohr radius a0.

H. Wagenknecht et al., Monte Carlo Simulations
21
0.0
−
−
−
10.0
20.0
30.0
MC simulations
Pade
formula
’
Fig. 2: Excess Internal Energy
of the OCP
1e−02
1e−01
1e+00
1e+01
Γ
θ is the degeneracy and is given by θ = T/T where T is the Fermi temperature.
F
F
4
Monte Carlo simulations for mass-symmetrical plasmas
The equilibrium state of an electron-positron plasma is simulated using the corrected
K,corr
Kelbg potential V
by a Monte Carlo run with 1600 particles. As a thermodynamic
characteristic of TCP we present in Fig. 5 the effective energy
ꢇ
eie
j
eff
U =
F(r , m , T, γ )
(21)
(22)
ij
ij
ij
4πε0
i<j
and the Coulombic energy
ꢇ
eie
j
Coul
U
=
4
πε r
0 ij
i<j
at a temperature of 10000 K as a function of Γ. The result for the corrected Kelbg
potential of Eqs.(5,10) is compared with molecular dynamics simulations by Valuev
eff,V
Coul,V
[
19, 20] (U
and U
), the classical Debye limit
√
Debye
3/2
= − 6Γ Nk T
U
(23)
B
Pad´e
with N the number of electrons, and a Pad ´e formula by Beule etal. [21] (U
)
√
6
Γ3/2
1/2
Pad´e
U
= −Nk T
(24)
(25)
B
1
4
1
+ (6π) Γ3/2τ1/2
Here τ is given by
2
2
4
τ = (2k Th (4πꢁ ) /m e ).
B
0
ep

22
Contrib. Plasma Phys. 41 (2001) 1
1
1
0
.5
.0
.5
1
.0
.5
.0
0
0
0.0
0.0
0.0
50.0
100.0
2.0
4.0
6.0
r/a0
r/a0
Fig. 3: g (r) for T=20000 K, Γ = 0.5, Fig. 4: g (r) for T=20000 K, Γ = 10,
ee
ee
1
9
−3
23
−3
n = 5.1 × 10 cm , θ = 34.3
n = 4.1 × 10 cm , θ = 0.08
e
e
Pair-correlation functions for the electron-electron g and the electron-positron g cor-
ee
ep
relations are shown in Figs. 6, 7 and 8 for three given temperatures: T = 45500K,
T = 98000K, T = 455000K. For equal particles the electron-electron and positron-
positron correlations are averaged over different spin directions. Since the masses of
electrons and positrons are equal it yields one single curve. The results are compared
with path integral Monte Carlo (PIMC) results [18] (gee,PIMC , gep,PIMC) and with molec-
ular dynamics (MD) simulations [19] (gee,MD , gep,MD). For the given temperature region
a qualitatively agreement of all three methods is observed. Only for a limited range
small divergences between the different results can be found.
5
Conclusion
An OCP and a mass symmetrical TCP was investigated by quasi-classical Monte
Carlo simulations. The interactions between particles were simulated by effective pair
potentials calculated by the method of Slater sums. We found the thermodynamic
functions from simulations in good agreement with other simulation methods and
with theoretical results. We could show that the pair-correlation functions are in good
agreement with PIMC and MD methods. In a forthcoming paper we will skip the
condition of equal masses and investigate the electron-proton plasma for a regime of
strong coupling and partial degeneration.
Acknowledgements
We would like to thank Burkhard Militzer from the University of Illinois at Urbana-
Champaign and Ilya Valuev from Moscow Institute of Physics and Technology for their
kindly providing the simulation data.

H. Wagenknecht et al., Monte Carlo Simulations
23
0.0
eff
U
U
U
U
U
U
−
2.0
Coul
eff,V
Coul,V
Debye
Fig. 5: Energies of TCP for
various Γ at T = 10000 K;
eff
effective
energy
U
U
,
,
Coul
Coulomb energy
effective energy of simu-
eff, V
lations of Valuev U
Pad
’
e
Coulomb energy of simu-
Coul, V
lations of Valuev U
,
classical Debye results, and
Pade’ formula.
−
4.0
0.0
1.0
2.0
Γ
4
3
2
1
0
.0
6.0
gee
gee
gep
.0
.0
.0
.0
gep
gee,PIMC
gep,PIMC
gee,MD
gep,MD
4
.0
.0
gee,PIMC
gep,PIMC
gee,MD
gep,MD
2
0.0
0.0
0
.0
2.0
4.0
2.0
4.0
r/a0
r/a0
Fig. 6: Pair-distribution functions for the Fig. 7: Pair-distribution functions for the
2
3
−3
22
−3
TCP at T = 45500 K, ne = 10 cm
Γ = 2.7, θ = 0.5
,
TCP at T = 98000 K, ne = 10 cm ,
Γ = 0.59, θ = 5

24
Contrib. Plasma Phys. 41 (2001) 1
6.0
4.0
2.0
0.0
gee
gep
gee,PIMC
gep,PIMC
gee,MD
gep,MD
Fig. 8: Pair-distribution
functions
for
the
TCP at T = 455000K,
2
3
−3
n = 10 cm , Γ = 0.275,
e
θ = 5
0.0
2.0
4.0
6.0
8.0
r/a0
References
[
1] Ichimaru, S., Rev. Mod. Phys. 54 (1982) 1017; Ichimaru, S., Hiroshi, H., Tanaka,
S., Phys. Rep. 149 (1987) 91
[
[
[
2] Perrot, F., Dharma-wardana, M.W.C., Phys. Rev. A 30 (1984) 2619
3] Stolzmann, W., Bl o¨ cker, T., Astron. & Astrophys. 314 (1996) 1024
4] DeWitt, H.E., Slattery, W.L., Stringfellow, G.S., in: Strongly Coupled
Plasma Physics, ed. Ichimaru, S., North-Holland, Amsterdam, 1990, p. 635
5] Zimmermann, R., Many-Particle Theory of Highly Excited Semiconductors,
TEUBNER-TEXTE zur Physik, Bd. 18, Teubner Verlagsgesellschaft, Leipzig, 1988
6] Blinnikov, S., Dunina-Barkovskaya, N.V., Nadyozhin, D.K., Astrophys. J.,
Suppl. Ser. 106 (1996) 171
[
[
[
[
[
7] Morita, T., Prog. Theor. Phys. 20 (1958) 920; 23 (1960) 829
8] Kelbg, G., Ann. Phys. 12 (1963) 219
9] Ebeling, W., Kelbg, G., Rohde, K., Ann. Phys. 21 (1968) 235
[
[
[
[
[
10] Rohde, K., Kelbg, G., Ebeling, W., Ann. Phys. 22 (1968) 1; 25 (1970) 80
11] Dorso, C., Duarte, S., Randrup, J., Phys. Lett. B 188 (1987) 287
12] Klakow, D., Toepffer, C., Reinhard, P.-G., Phys. Lett. A 192 (1994) 55
13] Ebeling, W., Schautz, F., Phys. Rev. E 56 (1996) 3498
14] Hansen, J.P., McDonald, I.R., Phys. Rev. A 23 (1981) 2041, formula (16); Phys.
Lett. 97 A (1983) 42; Bernu, B., Hansen, J.P., Mazighi, R., Phys. Lett. 100 A
(
1983) 28
[
15] Deutsch, C., Phys. Lett. 60 A (1977) 317; Deutsch, C., Gombert, M.M., Minoo,
H., Phys. Lett. 66 A (1978) 381; Minoo, H., Gombert, M.M., Deutsch, C., Phys.
Rev. A 23 (1981) 924

H. Wagenknecht et al., Monte Carlo Simulations
25
[
[
16] Cohen, E.G.D., Murphy, T.J., Phys. Fluids 12 (1969) 1404
17] Kahlbaum, T., in: Physics of Strongly Coupled Plasmas, ed. Kraeft, W.D.,
Schlanges, M., World Scientific, Singapore, 1995, p. 99
[
[
[
18] Ebeling, W., Militzer, B., Schautz, F., Contrib. Plasma Phys. 30 (1997) 553
19] Valuev, I.A., private communication (1999)
20] Ebeling, W., Norman, G.E., Valuev, A.A., Valuev, I.A., Contrib. Plasma Phys.
39 (1999) 61
[
21] Beule, D., Ebeling, W., F o¨ rster, A., Physica A 241 (1997) 719