J. Plasma Fusion Res. SERIES, Vol.7 (2006) 110-113

Monte Carlo Simulation of Cylindrical Langmuir Probe Sheaths
GOTO Makoto, KONDOH Yoshiomi1 , SHENG Min1 , TAKAHASHI Toshiki1 ,IWASAWA Naotaka1 and OKADA Tomio2
Kanto Polytechnic College, Oyama 323-0813, Japan
1 Gunma University, Kiryu 376-8515, Japan
2 Maebashi Institute of Technology, Maebashi 371-0816, Japan（Received: 5 October 2004 / Accepted: 21 december 2005）

Abstract
A self-consistent Monte Carlo modeling technique has been developed to study the sheath structure andproperty of cylindrical Langmuir probes. The introduced technique takes into account low-density and highenergy particles as well as enables the reduction of the calculation time of excitation collision cross sections.
The electron and ion current-voltage characteristics show good agreement with the probe theory and experimental results. The sheath structure and wide-range electron energy distributions are shown. The sweep timedependence on current-voltage characteristics is explained by the structure.
Keywords:
Monte Carlo simulation, cylindrical Langmuir probe, current-voltage characteristic, sheath structure,electron energy distribution
1．Introduction
An electric probe referred to as the Langmuir probe
is commonly used to diagnose plasma parameters (e.g.,electron temperature) for steady state, isotropic, and unmagnetized plasmas. Besides the measurement of quiet
plasmas, the requirement for observing the plasma parameters of RF glow discharges for etching, sputtering,and CVD applications have increased recently. Thisis due to the increasing role of plasma processing inmodern technologies. The results of simulation studies on RF glow discharges satisfactorily predict the experimental characteristics and provide useful information such as the spatial variation of plasma parametersin sheaths and plasma and their temporal evolutions.
This study is a benchmark to confirm that the simulation results show good agreement with the experimental current-voltage characteristics of the Langmuirprobe because it is established that the requirement for
further applications of simulation codes is inevitable.
In the present study, the particle simulations are carried out for the cylindrical Langmuir probe in argon
plasma. A fine subslab technique and the weight probability method are used [1, 2] and modified particularlyin the plasma source region (object of measurements)
so that high-energy particles are generated and simulated over a simulation space and time. Furthermore,some new techniques are introduced in order to reducethe calculation time.

The obtained electron current-voltage characteristics
are in good agreement with the Langmuir probe theory
and experiments. The electron current shows an exponential curve over three orders in the simulation resultsof a Maxwellian plasma source using 2000 particles in aring slab. The electron energy distributions are obtained
within a range of over four orders of magnitude. The
spatial variations of the electron and ion densities andpotential are shown, and the diﬀerence in the ion sheathand electron sheath is clarified. The sweep time dependence on current-voltage characteristics is explained by
the structure.
2．Simulation model and method
An infinitely long cylindrical Langmuir probe withradius ri (0.5 mm) is considered to be in the center of
an infinitely long cylindrical argon plasma with radiusro (5.5 mm), as shown in Fig. 1. The probe potential V p
with respect to the plasma end is swept from −5 V to 2 V
and from 2 V to −5 V with speeds of 0.025 V/µ s and0.2 V/µ s. The simulation particles of the electrons andions (Ar+ ) are traced under the influence of a radial selfconsistent electric field recalculated at the end of every
time step. The time interval dT s for the calculation ofthe electric field is 1.0 ns.
The simulation space is divided into Mm (= 50) ringslabs. A ring slab is divided into 10 subslabs. TheC 2006 by The Japan Society of Plasma

Science and Nuclear Fusion ResearchCorresponding author’s e-mail: goto@oyama-pc.ac.jp

110Monte Carlo Simulation of Cylindrical Langmuir Probe Sheaths

The radial electric field E r is calculated selfconsistently from the charge density distribution ρ byusing the following equation:  r  
 r0r01 11dξ
ρξdξ −
ρξ ln
Er =r ε0ln(r/r)ξ0iririVp+,(1)ln(r0 /ri )where ε0 is the permittivity in vacuum.
The electron collisions (elastic, excitation, and ionization) and ion collisions (elastic and charge exchange)are considered only in neutral argon gas. The electron
and ion collision cross sections that are used are provided by the JILA of Colorado University [3].
In order to reduce the calculation time of the excitation collision procedure, we introduced a combined
cross-section system. We combined 25 excitation collision cross sections Qi (U) with 25 diﬀerent thresholdenergies U i of argon. The combined excitation crosssection Qex (U) and energy loss function E loss (U) for anelectron with kinetic energy U are calculated as follows:

Fig. 1 Schematic of structure and circuit of simulation
model.

space resolution is restricted due to the separation ofsubslabs. The weight probability method and fine subslab system are used [1, 2]. Thus, all the simulation
particles have not only three-dimensional positions andvelocities but also individual weights. The standard particle number N s of any electron and ion in a ring slab
is 2000. The plasma source consists of two outer ringslabs and three neighboring ring slabs. The test particles in the two outer ring slabs are supplied (renewed)
during each dT s . The plasma source density N p is set
at 1.0 × 1016 m−3 , and the electron energy kT e and ion
energy kT i of the two outer ring slabs are 2/3 eV and1/37 eV, respectively, where k, T e , and T i are the Boltzmann constant, electron temperature, and ion temperature, respectively.
At the end of every time step, the new electron velocities of the two outer ring slabs are determined by random values between 0 and 4(kT e /m)1/2 , where m is themass of an electron. The weights of the simulation particles are calculated under the condition that the velocity distribution is Maxwellian and density is N p . Thus, asuﬃcient number of high-energy electrons can be simulated. In this procedure, the highest energy of an electron is below the ionization energy. Thus, the ionization
collisions are negligibly small. On the other hand, theion velocities of the two outer ring slabs are randomly
selected from a group of 10000 velocities that comprise
the Maxwell distribution. Both the outer ring slabs haveequal ion weights. The electron and ion weights in thethree neighboring ring slabs of the two outer ring slabs
are regulated at the end of every time step in order tomaintain a constant plasma density N p .

Qex (U) =25

Qi (U),(2)

i=1

Eloss (U) =25

i=1

Qi (U) · Ui

25
Qi (U).

(3)i=1

By using the combined cross-section system and energy
loss function, we can reduce the time required to accurately calculate the electron energy after many excitation collisions in a statistical manner.
The ion and electron simulation times are individually set, as shown in Fig. 2. The acceleration factor
Ac is introduced to reduce the calculation time. In every time step dT s , ions are simulated throughout; however, electrons are simulated only during dT s /Ac . The
dT s is determined such that the ion plasma oscillation
time is longer than dT s and the change in ion density
in dT s is small. The electron simulation time dT s /Ac
is determined such that the electron plasma oscillation

Fig. 2 Time ﬂow chart of electron and ion simulation.
111

Monte Carlo Simulation of Cylindrical Langmuir Probe Sheathstime exceeds dT s /Ac , and the motion of electrons during dT s /Ac is faster than that of the ions during dT s .
Our simulations are performed under the condition thatno ionization collision occurs. We obtain the electron
current by multiplying the electron charge flow in theprobe during dT s /Ac by Ac . When Ac is below 2000, it
is confirmed that the current-voltage characteristics arein good agreement with the Langmuir probe theory andexperiments. In this study, Ac is equal to 500.
3．Results and discussion
We first tested the simulation code for a quasi-idealcondition. Figure 3 shows the voltage characteristicsof the electron and ion current densities under low argon pressure (reduced pressure p0 = 13.3 Pa). Figure 3shows the two types of sweep directions. The increasing and decreasing probe potential are swept at the same
speed 0.025 V/µs. These two curves are almost in goodagreement. For this condition, the mean free path ofthe electrons and ions are 19 mm and 0.38 mm, respectively, and the Debye length is 0.06 mm. In this case,
it is considered that the sheath is thin and does not undergo collision. The electron energy kT e obtained fromthe electron current curve shown in Fig. 3 is 0.67 eV.
The plasma density obtained from the electron saturation current is 1.0 × 1016 m−3 . These are in good agreement with a parameter of the plasma source. The obtained floating potential is −3.4 V, and the ratio of thesaturation currents between electrons and ions is around200. These are almost in agreement with the value calculated from the Bohm criterion. These results may
provide some amount of evidence for validating thissimulation.
The ion current appears to be somewhat small, especially in the low negative- and zero-voltage regions. It

seems that ion current decay with a low negative probe
potential is primarily due to the collision of ions andalso the growth of a thin and weak reverse field surrounding the probe surface to a certain extent. It is confirmed that the ion current decreases significantly as thegas pressure increases. For the ions, there exists a potential wall with a value of around −kT i when V p is 0 V.
Although this may require more precise consideration,
it beyond the scope of this study.
Figure 4 shows the current density characteristics at a
sweep speed of 0.2 V/µs. These curves exhibit hysteresis characteristics. The electron current for a decreasing
sweep voltage from 0 V to −1 V tends to remain at thesaturation value. The two ion current curves are different over a wide range. The ion current for the decreasing voltage is less than that for the increasing voltage near a space potential. The relation varies inversely
over a deep negative probe potential. These may be explained qualitatively by considering the density distributions and mass of charged particles.
The radial distributions of the plasma parameters areobtained in the simulation shown in Fig. 3. We present
some typical results of the simulation. Figure 5 shows
the potential distributions of four probe potentials (−5,−2, 0, and 2 V). Figures 6 and 7 show the distributionsof electron density Ne and ion density Ni , respectively.
All the results for this situation contain spatiotemporal fluctuations. For example, the space potential in theplasma region varies with a magnitude of around kT e .
The time averaged value is considered.
The negative probe potential have a long-range influence up to the plasma source. When a negative potential
is applied, there is a growth in a wide sheath and presheath. In the long-range pre-sheath, this weak field
accelerates the ions. In order to sustain an eﬀectively

Fig. 3 Electron and ion current density characteristics at
a sweep speed of 0.025 V/µs.

Fig. 4 Electron and ion current density characteristics at
a sweep speed of 0.2 V/µs.

112

Monte Carlo Simulation of Cylindrical Langmuir Probe Sheaths

Fig. 5 Potential distributions.

Fig. 8 Electric energy distributions at Vp = –5 V.

sheath obtained at a negative probe voltage, while apositive voltage leads to the shallow sheath and prevents plasma density decay. A high plasma density near
the probe causes high electron and ion currents beginning from the positive probe potential; however, a large
mass of ions leads to a delay in the increase of ion current, thereby resulting in the hysteresis curve shown in
Fig. 4.
The electron energy distributions of four positions ata negative potential (V p = −5 V) are shown in Fig. 8.
The wide-energy distributions (over four orders of magnitude) are obtained. By employing the weight probability method, the supplied energy distribution in theplasma source has a wide range. Furthermore, the energy loss due to a collision is negligibly small. Thus,
the Boltzmann relation is confirmed, as shown in Fig.
8. On the other hand, the ion distributions in the sheath
are not Maxwellian. The accelerated ions have wider
band distributions.

Fig. 6 Electron density distributions.

4．Concluding remarks
A cylindrical Langmuir probe simulation has beendeveloped. The fine subslab technique and weight probability method are used and modified particularly in theplasma source region. A combined cross section system
is introduced to reduce the simulation time. The obtained electron and ion current-voltage characteristics
show good agreement with the probe theory and experimental results. The profiles of the potential and electronand ion densities are clarified, and wide-range electronenergy distributions are shown.

Fig. 7 Ion density distributions.

large ion flow, not only diﬀusion flow but also drift flowis required in the plasma. The electric field for gathering ions is governed by the electrons. The field is
controlled within a value where electrons play the roleof charge neutrality. On the other hand, the positive potential is perfectly shielded within a short region. Whenthe probe potential is positive or zero, the sheath is narrow and prevents plasma density decay near the sheath.
When the ion current is low or zero, a very weak field
can compensate the ion diﬀusion of the density gradient. The electrons are accelerated by a positive potential through the short region. These are verified fromthe mean energy distributions of the electrons and ions.
The features comprise the wide sheath and pre-

[1]
[2]
[3]

113

References
M. Goto et al., Jpn. J. Appl. Phys. 36, 4815
(1997).
M. Goto and Y. Kondoh, Jpn. J. Appl. Phys. 37,308 (1998).
JILA of Colorado Univ. home page:jilawww.colorado.edu/research/colldata.html.