5

THERMODYNAMICS OF PLASMAS WITH DEVELOPED

IONIZATION

5.1One–component plasma with the neutralizing charge background

We shall begin the description of the thermodynamics of strongly ionized plasma by discussing the most popular and well–studied model of the one–component plasma. A one–component plasma (OCP) represents a system of point ions placed in a homogeneous medium of charges of opposite sign (Brush et al. 1966; Hansen 1973; Pollock and Hansen 1973; Ng 1974; Lieb and Narnhofer 1975; DeWitt and Hubbard 1976; Galam and Hansen 1976; Zamalin et al. 1977; DeWitt 1977; Gann et al. 1979; Fisher et al. 1979; Baus and Hansen 1980; Slattery et al. 1980; Ichimaru 1982, 1992; Dubin and O’Neil 1999). Such a model serves as a good approximation for the plasma at ultrahigh pressure realized in the center of white dwarfs and heavy planets of Jupiter type (see Table 5.1). In these cases, matter is ionized under the e ect of pressure and degenerate electrons have su cient

kinetic energy, εkin ≈ εF = (3π2)2/3 n2e/3/(2m), to produce an almost uniform background density distribution. Due to the small electron mass, the kinetic energy of electrons at high density (rs → 0) is εF kT , and its pressure is much higher than that of an ion subsystem. In fact, two systems are realized in this case: the Coulomb system of point nuclei described by the Boltzmann statistics, and the quantum electron fluid. The weak interaction between these components results in an insignificant increase of the electron density in the neighborhood of nuclei (polarization), and most attention is concentrated on the analysis of the Coulomb internuclear interaction. The OCP model is the simplest, nontrivial model of a plasma, since there is no doubt concerning the form of the interaction potential, and the absence of quantum e ects enables one to exclude from the treatment the formation of bound states (of molecules, atoms, and ions) (Rushbrooke 1968; Kovalenko and Fisher 1973) and the e ect of degeneracy and interference (Rosenfeld and Ashcroft 1979; Springer et al. 1973; DeWitt 1976; Gorobchenko and Maksimov 1980; Ebeling et al. 1991). Therefore, the OCP model has been comprehensively studied, both analytically and numerically, in a broad range of nonideality parameters.

5.1.1Monte Carlo method

A great deal of results have been obtained by Zamalin, Norman, and Filinov (1977) in the framework of the OCP model by using the Monte Carlo method. This enables one to perform “computer” experiments with a plasma of any den-

169

170 THERMODYNAMICS OF PLASMAS WITH DEVELOPED IONIZATION

Table 5.1 One–component plasma of astrophysical objects

Parameter		Jupiter		White dwarf		Neutron star
Z		1 (H)		6 (C)		26 (Fe)
n	, cm−3	6	1024	5	1030	1032
i			·4		·8	10	8
T , K		10		10
γ		50		10–200		870
rs, a0		0.65		0.4 · 10−2		0.8 · 10−3

sity. The method employs the first principles of statistical physics and is based on direct computer calculations of mean thermodynamic values,

F = Q−1(V, N, T )	. . . F (q) exp{−βUN (q)}dN q,	(5.1)

where Q(V, N, T ) = . . . exp{−βUN (q)}dN q is the configuration integral defining the equilibrium properties of the thermodynamic system and q ≡ q1, . . ., qN are the coordinates of the particles. The interparticle interaction potential is assumed to be pre–assigned and, in most calculations, binary,

UN (q) =	Φab(rij ), rij = qi − qj .	(5.2)

a,b,i<j

By using a pseudo–random number generator, the Monte Carlo method in this case is reduced to the generation of a Markov chain which represents a set of space configurations, A1, A2, . . ., AN , with probability ωij of the Ai → Aj transition,

ωij = exp{−β[U (Aj ) − U (Ai)]}.

(5.3)

In doing so, the thermodynamic functions are averaged along the obtained Markov chain, which is equivalent to averaging over the canonical ensemble.

Direct calculations are usually performed for a single cubic elementary cell containing N particles. In order to take into account the interaction of these particles with those of the neighboring cells, the Ewald procedure and periodic boundary conditions were employed.

In the pioneering study by Brush et al. (1966), which involved the use of the Monte Carlo method in the OCP, the values of energy, heat capacity cV , and binary correlation function g(r) were calculated as the mean over α 105 configurations, with the number of particles in a cell being N 108 (for tests, N varied from 32 to 500).

Having the dense plasma of Jupiter in mind, DeWitt and Hubbard (1976) performed numerical calculations of a hydrogen plasma as well as a plasma of a mixture of light elements, taking into account the polarization of the background (rs = 0) by the method of linear response theory.

The calculations of Brush et al. (1966) describe OCP in a liquid state in the range 0.05 γ 100 (from 32 to 500 particles per cell). The calculations of

ONE–COMPONENT PLASMA

171

.

Fig. 5.1. The binary correlation function of the OCP for di erent values of the nonideality

parameter γ (Ichimaru 1982).

Hansen (1973) cover the range 1 γ 160 (from 16 to 250 particles per cell). The results of these studies are in good agreement, except for the region of the highest values of γ which is described more accurately by Hansen (1973). OCP in a crystalline state was used by Pollock and Hansen (1973) for the range 140 γ 300 (128 and 250 particles per cell). Detailed calculations of OCP properties were performed also by Ng (1974), Slattery et al. (1980), and Stringfellow et al. (1990).

5.1.2Results of calculation

The binary correlation function g(r) is presented in Fig. 5.1. Figure 5.2 demonstrates the behavior of the static structure factor,

S(q) = 1 + ni	dr[g(r) − 1]e−iqr.		(5.4)
The static dielectric permeability ε(q,		0) is given in Fig. 5.3,
ε(q, 0) = [1 − qD2 /q2S(q)]−1		, qD2 = 4πni(Ze)2/kT.	(5.5)

For γ 1, these functions vary monotonically and are described by the linearized Debye approximation,

SD(q) = q2(q2 + qD2 )−1, ε(q, 0) = (q2 + qD2 )/q2.

(5.6)

The pair correlation function is monotonic up to γ 2.5 and is in qualitative agreement with the Debye approximation,

g(r) = 1 − Γ	rD	exp(−r/rD),
	r		(5.7)

Γ = (eZ)2/(rDkT ), rD = (4πZ2e2ni/kT )−1/2.

172 THERMODYNAMICS OF PLASMAS WITH DEVELOPED IONIZATION

S

Fig. 5.2. The structure factor S(q) of the OCP for di erent values of the nonideality parameter γ (Ichimaru 1982).

.

Fig. 5.3. The static dielectric function of the OCP for di erent values of the nonideality parameter γ (Ichimaru 1982).

At γ 2.5, oscillations of g(r) appear, suggesting the emergence of short– range order: the system changes from the ideal gas to a liquid state. With an increase of γ, the oscillations increase and an e ective hard core is formed.

Given g(r), one can calculate the internal energy,

ONE–COMPONENT PLASMA

173

u = 3nikT /2 + uex,

uex/(nikT ) = (ni/2kT )

(5.8) dr(Z2e2/r)[g(r) − 1].

In the presence of the neutralizing background, the potential energy is not equal to zero solely because of the correlations in the charge positions. Therefore, the mean potential (excess) energy of interaction, uex, is often referred to as the “correlation energy”. For Γ 1, by substituting (5.7) in (5.8), one can easily obtain the classic a Debye–H¨uckel result,

√

uex/(nikT ) = −Γ2 = − 23 γ3/2.

Although the potential of the interparticle interaction is repulsive, the correlation energy is negative because of the neutralizing background. A “correlation hole” is formed around each charge, where the presence of another charge is improbable because of the repulsion. Then the major contribution to the potential energy is due to the interaction of the point–like charge with the background within the hole, which obviously corresponds to attraction.

The results of Monte Carlo calculations for 1 γ 160 were approximated by Slattery et al. (1980) with an accuracy of 3 ·10−5 by the following expression:

uex/(nikT ) = aγ + bγ1/4 + cγ−1/4 + d,

(5.9)

where a = −0.897 52, b = 0.945 44, c = 0.179 54, d = −0.800 49. In subsequent work by Stringfellow et al. (1990) and Dubin and O’Neil (1999) the accuracy increased substantially, which caused some changes in the approximations. By integrating Eq. (5.9) over γ, one can obtain the free energy of the OCP,

f = F/(nikT ) = aγ +4(bγ1/4 −cγ−1/4)+(d+3) ln γ −(a+4b−4c+1.135). (5.10)

The referencing of f in Eq. (5.10) to the point γ = 1 was done by Slattery et al. (1980) by using expressions for uex valid at γ 1.

Note that in accordance with the virial theorem, pV = u/3, a negative value of u at large γ leads to the negative pressure of the ion component. This, however, does not make the OCP unstable – the total pressure is, of course, positive because of the high value of the electron background pressure. The isothermal compressibility, χT , behaves similarly,

nikT χT	= 1 + 3	nikT		+ 9 dγ nikT .		(5.11)
1	1		uex(γ)		γ d uex(γ)

The results of the calculation for solidified OCP (Fig. 5.4) were given for the bcc lattice, which has the lowest energy. Slattery et al. (1980) approximated the excess energy at 160 γ 300 by the expression

	3
uex/nikT = abccγ +	2 bγ−2 = 0.895 929γ + 1.5 + 2980/γ2.	(5.12)

174 THERMODYNAMICS OF PLASMAS WITH DEVELOPED IONIZATION

Fig. 5.4. Binary correlation function for γ = 160. Solid line is for the solid phase (Pollock and Hansen 1973), broken line is for the liquid phase (Hansen 1973), dots are for the HMC approximation (Ng 1974).

For the physical interpretation of the individual terms in the approximate expressions, such as Eqs. (5.9) and (5.10), simplified models were used by Lieb and Narnhofer (1975), Ichimaru (1977), and Dubin and O’Neil (1999). We discuss them below.

5.1.3Ion–sphere and hard–sphere models for OCP

The ion–sphere model (valid for γ 1), or the cell model, breaks the system up into spheres with radius rs. Each sphere (called the Wigner–Seitz cell) accommodates an ion surrounded by a homogeneous cloud of neutralizing electrons the density 3Ze/4πrs3. Since the cells are not overlapped and the total charge is equal to zero, the interaction between them is absent and the potential energy is simply equal to the sum the energies of individual cells. Within this model, one can easily calculate the mean excess energy,

uex/(kniT ) = −	9	γ +	3	−	9	γ.	(5.13)
	10		2		10

Note that Eq. (5.13) satisfies one of the global inequalities imposed on OCP by the singularity of the Coulomb potential (Lieb and Narnhofer 1975). Equation (5.13) should be compared with Eq. (5.10), as well as with Eq. (5.12). In the latter, abcc is the Madelung constant for the bcc lattice. Therefore, the dominant term in Eq. (5.12) is static, and the other terms represent the “thermal” part associated with harmonic oscillations of ions. The contribution of the thermal part to uex, p, and χ proves to be less than 3% at γ 100.

As we can see, the leading term in expressions (5.9) and (5.10) is close to the Madelung term for the bcc lattice of Coulomb charges. OCP yields the bcc