SciELO - Scientific Electronic Library Online

 
vol.63 issue1Relatividad general; su presencia temprana en MéxicoInvestigación en enseñanza de la física experimental en el siglo XXI author indexsubject indexsearch form
Home Pagealphabetic serial listing  

Services on Demand

Journal

Article

Indicators

Related links

  • Have no similar articlesSimilars in SciELO

Share


Revista mexicana de física E

Print version ISSN 1870-3542

Rev. mex. fís. E vol.63 n.1 México Jan./Jun. 2017

 

Education

On the accuracy of the Debye shielding model

M.A. Martínez-Fuentesa  * 

J.J.E. Herrera-Velázqueza 

aInstituto de Ciencias Nucleares, Universidad Nacional Autónoma de México, Apartado Postal 70-543, Ciudad Universitaria, 04511 México, D.F., México.


ABSTRACT

The expression for the Debye shielding in plasma physics is usually derived under the assumptions that the plasma particles are weakly coupled, so that their total kinetic energy is much greater than their electrostatic interaction energy, and that the velocity distributions of the plasma species are Maxwellian. The first assumption also establishes that the number of particles within a sphere with a Debye radius, known as the plasma parameter ND, should be significantly greater than 1, and determines the difference between weakly and strongly coupled plasmas. Under such assumptions, Poisson’s equation can be linearized, and a simple analytic expression is obtained for the electrostatic potential. However, textbooks rarely discuss the accuracy of this approximation. In this work we compare the linearized solution with a more precise numerical (or “exact”) solution, and show that the linearization, which underestimates the “exact” solution, is reasonably good even for ND ∼ 40. We give quantitative criteria to set the limit of the approximation when the number of particles is very small, or the distance to the test charge too short.

Keywords: Plasma physics; debye shielding; electrostatics

PACS: 52.27.Gr; 41.20.Cv; 02.30.Hq

1. Introduction

One of the most elementary concepts in classical plasma physics is that of Debye shielding, which establishes the characteristic distance λD in which the electrostatic field from a charged test particle can be shielded by particles of the opposite sign. If a test charge is placed in a homogenous plasma, in which the unperturbed electron and ion densities are the same (i.e.: satisfy the quasi-neutrality condition), it repels the charges of its own sign, while it “dresses” itself with charges of the opposite sign, resulting in an electrostatic shielding. The derivation of λD has been well presented with similar arguments in many textbooks, such as Refs. 1 to 7. In doing so, two assumptions are generally made: (i) The total kinetic energy of the particles is much greater than the electrostatic interaction energy between them, which means that the plasma is weakly coupled, and allows the linearization of the problem, (ii) Each particle species is in thermodynamic equilibrium, so their velocity distributions are Maxwellian. The latter assumption requires the density to be large enough to allow collisions to an equilibrium state, and does not apply for collision-less, non-equilibrium plasmas, such as those often found in space plasmas.

Modern textbooks discuss the case of strongly coupled plasmas, for which the first constraint does not apply 3,5,8, although they usually ignore the fact that many collision-less plasmas, particularly in space physics, even in steady state, are non-equilibrium systems, and therefore not Maxwellian. The case of shielding in quantum plasmas, in which Fermi statistics should be considered, has been presented in Ref. 8, but it will not concern us in this paper.

The question on the relaxation of assumption (ii) has been addressed by Bryant 9, who studied the shielding produced when the species are described by the so called kappadistributions for non-equilibrium plasmas. This is a more general family of distributions that include the MaxwellBoltzmann one as a special case, and is motivated by observations in space plasmas. While Bryant’s conclusion is that the shielding λD is shorter than for the Maxwellian case, given the same number density and mean energy, the weakly coupled assumption is kept. More recently, Livadiotis and Comas have reviewed the subject, providing a better theoretical ground for the kappa-distributions, under the point of view of non-extensive entropy 10. At a didactical level MeyerVernet 11 studied the case in which a Boltzmann equilibrium is not assumed, from a different point of view.

On the other hand, the weakly coupled plasma assumption, also poses the requirement that the number of particles ND within a sphere of radius λD be much greater than 1, so the case of classical strongly coupled plasmas, in which ND ∼ 1 (in the opposite limit of that of quantum plasmas, for which ND ≫ 1) is left open. It is also important to note that the linearized solution for the electrostatic potential, which results from the weakly coupled plasma assumption, is also limited in space, and breaks down for r < λD. The relaxation of the weakly coupled plasma assumption has been addressed by Mak 12, who compared the usual linearized solution for the electrostatic potential with the exact solution of Poisson’s equation for several cases in which ND À 1. While the conclusion was that the approximate solution is good beyond expectation for the cases studied, a few questions went unanswered. Particularly, when does the approximate solution breaks down as ND and the distance r are reduced?

The purpose of this paper is to provide a better insight into the subject, using the same normalised approach as in Ref. 12, for comparison. In Sec. 2 we go through the usual discussion of the problem, reviewing what has been presented in textbooks, for the sake of completeness and in order to establish the notation. In Sec. 3 we compare the approximate analytical solution to the linearized problem with the exact numerical one, and in Sec. 4 we summarise the main conclusions. The paper is written with sufficient detail, for easy reading by an advanced undergraduate or graduate physics student.

2. Linearized Debye shielding

2.1. The Poisson equation

Let us consider a homogeneous plasma, in which the unperturbed large-scale average electron and ion densities neo, and nio, are equal: neo = nio = n. If a point test charge is introduced at the origin, it will create an electrostatic potential Φ(r), which only depends on the radial coordinate r, and which must satisfy Poisson’s equation

1r2ddrr2dΦ(r)dr=-ρc(r)ε0, (1)

where ρc(r) is the charge density distribution. We make the usual assumptions (i) and (ii) mentioned in the introduction, and following Ref. 12, we make an additional assumption: (iii) The positive ions are protons with infinite inertia, so they form a uniform background of density n. Since the electron density would then be given by ne(r) = nexp(eΦ(r)/kT), where −e is the charge of the electron, k is Boltzmann’s constant, and T is the electron temperature, then ρc (r) = ne[1−exp(eΦ(r)/kT)], and Eq (1) can be rewritten as

d2dr2[rΦ(r)]=-nerε0[1-expeΦrkT]. (2)

The numerical solution to this equation is what we shall call the “exact solution”.

2.2. The limit between weakly and strongly coupled plasmas

In this subsection we follow the discussion of Refs. 6 and 7, with some adaptations. The boundary between weakly and strongly coupled plasmas can be established when the electrostatic energy between the plasma particles is equal to the kinetic energy between them, so weakly coupled plasmas are those for which the former is smaller than the latter, and the strongly coupled plasmas are represented by the opposite case. If we consider an electron with velocity v and charge −e in the presence of an ion with charge Ze at rest, the distance of maximum approach rc will be defined by

12mv2-Ze24πε0rc=0 (3)

Taking v as the thermal velocity vth, in the case of three dimensions for an isotropic plasma, we can write (1/2)mvth2=(3/2)kT. Since Φ(r) = Ze/4πε0rc,

eΦ(r)kT=32, (4)

which can be taken as the limiting case in Poisson’s equation, such that weakly coupled plasmas happen for the right hand side is smaller than 3/2, and the opposite case for strongly coupled plasmas. For a weakly coupled plasma,

eΦ(r)/kT1, (5)

in which we may take 1 instead of 3/2. Let us now define the Debye length λDε0kT/ne2 1-7, whose significance will become clear in the linearized case. By taking the average distance between particles rd = n−1/3, (3) translates into

rdrc=12πZλD2rd2 (6)

from which it becomes clear that weakly coupled plasmas are those for which both rd and λD are greater than the maximum approach distance rc. Using this result, it is easy to find that, if we define Λ = λD/rc, as the ratio between the Debye length (the distance in linearized theory at which the electrostatic potential decays to e−1 in the linear theory studied in the next subsection,) and the distance of maximum approach, we find that in the limit between weakly and strongly coupled plasmas, Λ = (3/Z)ND, where 3ND4πnλD3 is the number of particles within a sphere defined by the Debye length (This definition for ND is for the sake of comparison with Ref. 12.) The familiar statement that in weakly coupled plasmasND > 1 , while the opposite is true for strongly coupled plasmas, is recovered, as well as (Z/3)Λ > 1, for which the limit increases linearly with Z. All follow from assumption (5).

2.3. The weakly coupled limit

The linearization for the weakly coupled plasmas arises after considering Φ(r)/kT ≪ 1, which states the potential energy between the particles is much smaller than their kinetic energy, and is equivalent to NDnλD31. This allows an expansion of the exponential in (2). Keeping the first two terms this yields for the approximate potential Φa(r)

d2dr2[rΦa(r)]=rΦa(r)λD2, (7)

whose solution is the well known Yukawa potential:

Φa(r)=Aexp(-r/λD)r. (8)

By taking the integration constant A = Ze/(4πεo), this allows the usual interpretation of the potential produced by a point charge of an ion Ze, with atomic number Z, at the origin (ρc(r) = Zeδ(r)), where δ(r) is the Dirac delta function), shielded by the electrons with a characteristic decay length λD. Once the electrostatic potential is known, the charge distribution for the linearized solution can be obtained from Eq. (1)

ρc(r)=Ze4π4πδ(r)-exp(-r/λD)rλD2. (9)

The subtlety in obtaining (9) is discussed in Appendix. The effect of the shielding can be appreciated by computing the charge within a sphere of radius r:

Q(r)=Ze[1-0rr'2ρρc(r')dr']=Ze(1+r/λD)exp(-r/λD). (10)

Thus, the charge is Ze at r = 0, and decays to 0 for large r. At the Debye radius the charge is Q(λD ) = 2exp(−1)Ze ∼ 0.736Ze, and falls to half its value at the origin, 0.5Ze, at r/λD ∼ 1.68, well beyond the Debye radius. This result is interesting, when considering λD as the cut-off impact parameter, when Coulomb collisions are considered.

3. Comparison between an accurate (numerical) and the approximate (analytical) solution

For the sake of comparison with Ref. 12, let us normalise the distance r and the potential Φ(r) in terms of the Debye length, changing to the following dimensionless variables:

ρ= r/λD, (11)

Ψρ= 4πε0λDΦ(r)/Ze=(eΦ(r)/kT)ND, (12)

Thus, Poisson’s Eq. (2) for the potential can be rewritten as

d2dρ2[ρΨ(ρ)]=-NDρ[1-exp(Ψ(ρ)/ND)], (13)

whose numerical solution we shall call the “exact” solution, while the normalized approximate solution would then be

Ψa(ρ)=exp(-ρ)ρ, (14)

valid for

eΦ(r)kT=1NDρexp(ρ)1. (15)

The main purpose of this paper is to compare the numerical solution Ψ(ρ) to Eq. (13), with the analytic solution (14) to the linearized equation.

Indeed, the approximation should be good for large distances, but Eq. (13) also tells us that it will break down for short distances (exp(−ρ) ∼ 1) when NDρ ≤ 1. This will happen if ND< 1, in very diluted plasmas, or for short distances, when r < λD. Also, remembering Subsec. 2.1, note that NDρ = 4πnλ2Dr = (Z/3)(r/rc), which means that for Z = 1, NDρ > 1 for r > 3rc. Therefore, the weak plasma approximation is expected to break around

r3rcexp(-r/λD)1 (16)

for short distances.

However, in order to fully understand what all this means, it is important to quantify these statements. The main purpose of this paper is to compare the numerical solution Ψ(ρ) to Eq. (13), with the analytic solution Ψa (ρ), (14), to the linearized equation.

The numerical solutions are obtained by means of a fourth-order Runge-Kutta routine 13, which is started at the tail of the solution, and integrated backwards, from larger to smaller distances. The initial values ρ0 for the integration can be chosen to be smaller for larger values of ND, because the validity of the approximation breaks at shorter distances. The way in which they were chosen was such that in the first step of integration we satisfy |Ψa −Ψ| < 10−11. It was found that if this condition was not met, the integration would not be independent of the initial condition, which is something else we sought. This of course, also depends on the step of integration. Another criterion used in order to choose the step of integration was based on the tolerance at which the approximation would break. For this purpose, three different tolerances, τ, were tried: 1, 5, and 10%. By trial and error, the integration steps were chosen in such a way that the distance ρd at which (|Ψa − Ψ|/Ψ) × 100 = τ occurred, would be such that τ could be determined within a 0.1% uncertainty.

A first result of this calculation we found, is that the approximate solution underestimates the exact one. Table I shows, for several values of ND , the initial values ρ0 of the integration, and the values ρd for which tolerances of 1, 5 and 10% between the two solutions fail. The same results are plotted in logarithmic scale in Fig. 1. For classical, weakly coupled plasmas, relevant ND values run from 40, for the Solar atmosphere, some gas discharges, and laser produced plasmas, to 108 for the interstellar gas, while for the Solar corona and thermonuclear plasmas, they are around 10714. In addition, we have explored smaller values of ND, including some smaller than one, for which a strongly coupled plasma is expected.

Table I Initial value ρ0 at which the integration of (13) was started, and distance ρd at which tolerances of 1, 5, and 10% between the approximate and exact solutions are reached, for different values of ND

ρ0 ρd
ND 1% 5% 10%
0.1 13 3.569 2.314 1.828
0.3 11 2.686 1.510 1.077
1 11 1.775 0.734 0.401
2 11 1.288 0.368 0.137
5 7 0.708 5.11×10−2 1.94×10−2
10 5 0.340 8.43×10−3 7.22×10−3
40 5 1.83×10−3 1.40×10−3 1.34×10−3
1×102 5 5.51×10−4 4.86×10−4 4.71×10−4
5×102 4 8.70×10−5 8.05×10−5 7.87×10−5
1×103 2 4.02×10−5 3.75×10−5 3.68×10−5
1×104 2 3.27×10−6 3.10 ×10−6 3.06×10−6
1×105 0.1 2.77×10−7 2.66×10−7 2.63×10−7
1×106 0.01 2.41×10−8 2.33 ×10−8 2.31×10−8
1×107 0.001 2.14×10−9 2.08×10−9 2.07×10−9

Figure 1 (In colour in the on-line version) Distances ρd at which the approximate solution breaks, for different values of N D. The blue dots are for the case in which a 1% tolerance is taken, the green squares for 5%, and the red diamonds for 10%. The blue line (below) is the fit for results that go from ND = 40 to 1×107. The red line (above) is the curve log10ρ = −log10ND

For large values of ND, the distance values at which the approximation breaks down are very similar regardless of the tolerance. In the plot in Fig. 1, they overlap down to ND = 40. As ND is further reduced, the distance at which the approximation breaks down increases, and is larger for smaller tolerances. However, it is interesting to observe that the differences obtained for different tolerances tend to diminish for fractional values of ND. While we have followed the calculations into very small values of ND, it should be noted that the description of the charge density, used in the Poisson’s equation, in terms of a Boltzmann distribution for such a case may lose its meaning, due to the increasing importance of the temporal fluctuations about the average. Indeed, the use of the Maxwell-Boltzmann distribution assumes the particle species are in equilibrium, which requires the density to be large enough for collisions to allow relaxation to such a state.

Thus, the linear approximation is very good when ND ranges between 40 and 107, and it is possible to fit the results to the curve log10ρd = 1.18828−1.07493ND (lower line in Fig. 1), which yields the empirical law

ND1.07493ρd=0.06482. (17)

This should be compared to the curve log10ρ=−log10 ND (upper line in Fig. 1), which stands for the limit NDρ = 1. Therefore, our result gives a quantitative meaning to the statement NDρ ≪ 1. We should note that this result differs with that found in Ref. 12:

NDρ=0.01 (18)

Our work goes further than that of Ref. 12, in that we study the cases of smaller values of ND. Besides, only the case in which the tolerance is 10% was reported. We did not explore ND as large as in that previous work (1015), because it is to be expected that the approximation will improve as ND increases, while we are more interested in exploring when the approximation breaks down.

4. Conclusions

In order to study the Debye shielding, the numerical (exact) solution for the electrostatic potential Ψ(ρ) was obtained, from Eq. (13), and it was compared to the approximate solution, given by the Yukawa potential (14), found for the linearized approximation. The latter underestimates the exact solution, but is a good approximation for values of ND as small as 40. Even for ND = 5, the approximation is still good, down to r = 0.05λD if one requires only a 5% tolerance. The distances ρd, at which the approximation breaks, given fixed values of ND, for tolerances of 1, 5 and 10% were obtained, and it was observed that they are practically the same, regardless of the tolerance, for ND > 40. Using these data, the empirical law (17) was found, which gives a clear quantitative meaning to the statement NDρ ≪ 1, for the validity of the linear approximation.

From the didactical point of view, this result is illustrative of the validity of the usual approximation, which is never quantitatively explained in textbooks, even when a Maxwell Boltzmann equilibrium distribution is assumed, which may not necessarily be right in many cases.

Finally, it must be noted that the results reviewed in this paper are valid for classical plasmas, such that the distance of maximum approach rc, is much greater than the electron’s thermal wavelength λe=2mkT. If ND ≫ 1the electrons become degenerate, and the interference between wave functions must be considered for high temperatures. In this case, the classical model is no longer valid, and the shielding is provided by the Thomas-Fermi length rTF=(π/3ne)2/4me2 (Ref. 8).

REFERENCES

1. D.R. Nicholson, Introduction to Plasma Theory (New York: John Wiley & Sons, 1983) p. 1-5. [ Links ]

2. R.J. Goldston and P.H. Rutherford, Introduction to Plasma Theory (Bristol: IOP Publishing Ltd., Bristol, 1995) p. 13-16. [ Links ]

3. P.H. Bellan, Fundamentals of Plasma Physics (Cambridge: Cambridge University Press, 2003) p. 7-11. [ Links ]

4. J.A. Bittencourt, Fundamentals of Plasma Physics Third Edition (New York: Springer Verlag, 2004) p. 273-279. [ Links ]

5. A. Piel, Plasma Physics: An Introduction to Laboratory, Space and Fusion Plasmas (Berlin: Springer Verlag, 2010) p. 35-39. [ Links ]

6. R. Hazeltine and F. L. Waelbroeck, The Framework of Plasma Physics (Reading, Mass: Perseus Books, 1998) p. 7. R. [ Links ]

7. R. Fitzpatrick, Plasma Physics: An Introduction (Boca Raton: CRC Press, 2015) p. 6 [ Links ]

8. V. Fortov, I. Iakubov and A. Khrapak, Physics of Strongly Coupled Plasma (Oxford: Clarendon Press, 2006). [ Links ]

9. D.A. Bryant, J. Plasma Phys. 56 (1996) 87. [ Links ]

10. G. Livadiotis and D.J. Comas, J. Plasma Phys. 80 (1993) 341. [ Links ]

11. N. Meyer-Vernet, Am. J. Phys. 61 (1993) 249. [ Links ]

12. P.H. Mak, Plasma Physics and Controlled Fusion 34 (1992) 453. [ Links ]

13. S.E. Koonin and D.C. Meredith, Computational Physics (Westview Press, 1990) p. 32. [ Links ]

14. J.D. Huba, NRL Plasma Formulary (Washington, D.C: Naval Research Laboratory, 2016) p. 40. [ Links ]

15. W. Greiner, Classical Electrodynamics (New York: Springer Verlag, 1996) p. 34-36. [ Links ]

Appendix

A. The calculation of the charge density

The charge density distribution (9) can be obtained from Poisson’s equation, from the potential (8), by appropriately dealing with the singularity. As explained in Classical Electrodynamics textbooks (take for instance Ref. 13), the potential produced by a point charge Ze at the origin

Φp(r)=Ze4πε0r (A.1)

satisfies the Poisson equation

2Φp(r)=-Zeδ(r)ε0, (A.2)

where δ(r) is the Dirac delta function.

In the Debye shielding case, we provide two ways to proceed. The first one, suggested by Greiner in Ref. 15, is to sum and subtract the point charge potential in (8). Therefore,

2Φa(r)=Ze4πε02[1r+exp(-r/λD)-1r]=Ze4πε0[-4πδ(r)+exp(-r/λD)rλD2] (A.3)

which yields (9), using (1). The second way may be by using the identity of the Laplacian for a product of two scalar functions ∇2(fg) = ∇2(f)g + ∇f · ∇g + f2(g), where f = Ze/(4πε 0 r), and g = exp(−r/λ D ). By these means, we get

2Φa(r)=Ze4πε0[-4πδ(r)exp(-r/λD)+exp(-r/λD)rλD2], (A.4)

which is equivalent to (A.3), since when integration is made to obtain the charge within the sphere of radius r, it yields the same result (10) for the total charge in a sphere of radius r, due to the property of the Dirac delta function

f(r)δ(r-r0)dV=f(r0), (A.5)

when integrating over a domain containing r0.

In plasma physics textbooks, this basic electrostatics is sue is seldom discussed, although Refs. 3 to 5, and 7 do mention the first term in the charge density, but not as a result of the Laplacian of (8). We believe it is a good example to stress the importance of properly dealing with singularities.

Received: October 14, 2016; Accepted: December 05, 2016

Present address: Laboratorio Nacional de Fusion, Ciemat, Madrid, Spain. e-mail: herrera@nucleares.unam.mx

Creative Commons License This is an open-access article distributed under the terms of the Creative Commons Attribution License