Two-dimensional particle-in-cell simulations of transport in a ...
JOURNAL OF APPLIED PHYSICS 108, 103305 共2010兲
Two-dimensional particle-in-cell simulations of transport in a magnetized electronegative plasma E. Kawamura,a兲 A. J. Lichtenberg, and M. A. Lieberman Department of Electrical Engineering, University of California, Berkeley, California 94720-1770, USA
共Received 27 July 2010; accepted 27 September 2010; published online 19 November 2010兲 Particle transport in a uniformly magnetized electronegative plasma is studied in two-dimensional 共2D兲 geometry with insulating 共dielectric兲 boundaries. A 2D particle-in-cell 共PIC兲 code is employed, with the results compared to analytic one-dimensional models that approximate the end losses as volume losses. A modified oxygen reaction set is used to scale to the low densities used in PIC codes and also to approximately model other gases. The principal study is the limiting of the transverse electron flow due to strong electron magnetization. The plasma in the PIC calculation is maintained by axial currents that vary across the transverse dimension. For a cosine current profile nearly uniform electron temperature is obtained, which at the B-fields studied 共600–1200 G兲 give a small but significant fraction 共0.25 or less兲 of electron to negative ion transverse loss. For a more transverse-confined current, and approximating the higher mass and attachment reaction rate of iodine, the fraction of electron to negative ion transverse loss can be made very small. The models which have been constructed reasonably approximate the PIC results and indicate that the cross-field transport is nearly classical. © 2010 American Institute of Physics. 关doi:10.1063/1.3506518兴 I. INTRODUCTION
Ion-ion plasmas composed of positive and negative ions, with a negligible fraction of electrons, have potential applications for microelectronics etching1 and plasma thrusters for space missions.2 A method of producing such plasmas, in which the ion-ion plasma is formed at the periphery of a magnetized electronegative plasma core, has been explored in several experiments.3–7 Recently, a one-dimensional 共1D兲 fluid model, together with an analytic theory, was developed8 to approximate the configuration of one of these experiments.6,7 In that model, the two-dimensional 共2D兲 nature of the device was approximated, by including a volume axial loss frequency L that corresponds to the end losses. Rectangular, rather than cylindrical, coordinates were used, for simplicity, and the value of the magnetic field B0 was chosen such that the electrons were strongly magnetized, while the ions were only weakly magnetized. To make the numerical integration tractable, the electron temperature Te was assumed uniform. One then looked for special solutions, at a transverse boundary y = R, for which ne共R兲 = 0,
冏 冏 dne dy
= 0,
共1兲
y=R
as a function of Te, the ratio of central electron density to neutral density ne0 / ng, and the ratio of central negative ion to electron density ␣0 = nn0 / ne0. The boundary conditions imply that a pure ion-ion plasma exists at y = R. Solutions were found in the ne0 / ng − ␣0 plane for various values of R and length L, holding RL = constant. The latter constraint was useful in the numerics because the approximation used for L was a function of that constraint. The various assumptions a兲
Electronic mail: [email protected].
0021-8979/2010/108共10兲/103305/11/$30.00
that were made in that work remain open to investigation. The solutions do not imply the existence of a physical transverse wall at the position for which the specified electron boundary conditions are satisfied. Furthermore, at a physical wall where the particular boundary conditions are satisfied, it does not ensure that the special solution, rather than another solution at a different Te would be favored. It is also necessary to determine the validity of the expression used for L, and the conditions under which Te is relatively uniform over the cross-section. To examine the approximations of the numerical calculations and accompanying model, and to investigate other configurations of interest, a 2D particle-in-cell 共PIC兲 code is employed using an oxygenlike feedstock gas. PIC parameters are scaled to higher density, using a rescaled oxygen cross-section set, and the dimensions and pressure are adjusted to keep pd = const. Since the pressure is adjusted, the magnetic field is also adjusted to preserve the relative strength of the magnetization. This will become clearer from the analytic models which we develop to compare to the PIC simulations. In Sec. II, we describe the 2D PIC code and present some typical PIC simulation results for a uniform structure with nearly uniform Te, related to our previous work. In Sec. III, a model is developed to approximate the electron and ion transverse profiles observed in the PIC simulations, and the model results are compared with PIC simulations over a range of parameters. Some implications of various configurations for producing ion-ion plasmas are discussed. In Sec. IV, we apply these ideas to iodinelike parameters with a modified heating configuration that produces nearly electronfree ion-ion edge plasma. A two-region model is developed and compared to the PIC simulation. In Sec. V, a more complicated, but more practical, geometry is introduced, with PIC simulations showing that an edge ion-ion plasma is ob-
108, 103305-1
© 2010 American Institute of Physics
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jap.aip.org/jap/copyright.jsp
103305-2
J. Appl. Phys. 108, 103305 共2010兲
Kawamura, Lichtenberg, and Lieberman
tained downstream from a narrower source region. These results are compared to a coupled two-region global model. The results are summarized in Sec VI. II. 2D PIC SIMULATIONS
We use a 2D3v rectangular PIC code PDP2 共Ref. 9兲 to simulate an oxygen discharge in a rectangular metal chamber, with insulating walls. In this section, as well as Secs. III and IV, the chamber is 20 cm long with a parallel magnetic field B0 in the x-direction and 3 cm wide in the y-direction. A shorter device of 10 cm length was also explored. Because of the insulating walls of thickness 1/2 cm in the x direction and 1/4 cm in the y direction, the actual discharge dimensions were length L = 19 cm and width 2R = 2.5 cm, for the longer device. For the shorter system, the insulating walls were 1/4 cm in all directions, resulting in L = 9.5 cm and 2R = 2.5 cm. By 2D3v, we mean 2D displacement with three velocity components 共3v兲. In our simulation, each computer particle represents a cluster of about 106 real particles 共electrons or ions兲. The simulations are run with a sufficient number of particles to minimize the discrete particle noise. For each particle, PDP2 tracks the two displacement, and the three velocity components. Collisions between the particles and the background gas are included by using a Monte Carlo collision 共MCC兲 algorithm. A detailed description of the PIC-MCC model used is given in reference.10 In a typical electrostatic PIC simulation, for each timestep: 共i兲 the charge density 共x , y兲 is obtained by a linear weighting of the particles to the spatial grid; 共ii兲 共x , y兲 is used in Poisson’s equation to solve for the electric field E共x , y兲; 共iii兲 the force F = qE is used to advance the particles to new x and v; 共iv兲 the boundaries are checked, and out-ofbounds particles are removed; and 共v兲 a MCC handler checks for collisions and adjusts the particle velocities and numbers accordingly. PIC simulations obtain the fields, as well as the particle fluxes and densities, self-consistently, without making any assumptions about the particle temperatures or velocity distributions.
from 300 to 600 G. Following the experiments which used oxygen as the feedstock gas,4,6 the parameters for oxygen were employed.10,11 However, to scale the low PIC density to a higher physical density, the positive-negative ion recombination rate coefficient Krec was scaled by a factor Krecfac, since the Krecne product occurs in the equations 共see Ref. 11兲. For example, a scaling factor Krecfac = 20 is used to scale to a factor of 20 higher density than that used in the PIC simulation. This still results in a relatively low recombination rate as was investigated in most experiments and previous theory. Higher recombination rates are also explored as described below. Because the mass of oxygen atoms and molecules is not very large, the oxygen ions were significantly magnetized. To avoid this effect in understanding the magnetization of only the electrons, the ions were not magnetized in the fluid calculation. It is straightforward to investigate PIC results, and also the accompanying models, with both magnetized and unmagnetized ions, and this has been done in the results presented below. In the following sections, we also investigate the effect of higher masses, of heating that does not produce uniform transverse temperatures, and configurations that are not axially uniform. The plasma production and heating in an experiment is with a complex inductive field that probably involves the excitation of waves.4,6 In the PIC calculations a simpler form of plasma heating and maintenance is used: an rf current density directed parallel to the B-field. The local plasma conductivity acting on the current produces an electric field that accelerates and thus heats the plasma electrons. In the first set of PIC simulations, at 20 mTorr, Te is equilibrated along the field lines, but not necessarily across the field lines. Since our previous equilibrium and model calculations assumed uniform Te, for the PIC simulations in this section we designed a current profile to approximately achieve this result. Since, from the form of the plasma conductivity, the rf electric field was expected to vary as E ⬀ J / n, where J is the local current density, to achieve E ⬇ constant, transversely, we explored current profiles that might approximate J ⬀ n. The two transverse profiles used are
A. The PIC requirements
For stability and accuracy, 2D PIC simulations require a large number of particles per Debye length squared. Thus PIC codes run faster when simulating small plasma reactors with lower densities 共larger Debye lengths兲. Consequently, the PIC reactor was chosen to have a width of 3 cm and a length of 10–20 cm, about a factor of two to three smaller than typical experimental devices to which it can be compared. The average electron density in the device was below 1014 m−3, well below experimental parameters. These modifications of dimensions and density to accommodate PIC requirements necessitated rescaling of various parameters, as described below. The PIC dimensions, which are approximately half of the dimensions previously used in fluid equations8 to study the experiment,4,6,7 require increasing the pressure in the base case from 10 mTorr in typical experiments to 20 mTorr, to keep pd = const. To keep the degree of electron magnetization constant 共see Ref. 5兲 the magnetic field was also increased by a factor of 2, in the base case
and
冉 冊 冉 冊 4x y sin 2R L
共2兲
冉 冊 冉 冊
共3兲
Jx = J0 cos2
Jx = J0 cos
4x y sin . 2R L
The former approximates the profiles of the special solutions, satisfying Eq. 共1兲, while the latter is a profile found for equilibrium electron densities not dominated by either magnetization or negative ion effects. Since electrons move relatively freely along the B-field, the axial form of J was not considered important and a double sine with zero at the ends was used, to avoid excessive heating at the axial edges. Values of J0 in the range of 1 – 10 A / m2 were tried in the calculations. The higher values were initially considered to produce electron densities closer to the central range of the densities explored in experiments3–7 and the previously reported fluid calculations.8 However, an unexpected result
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jap.aip.org/jap/copyright.jsp
103305-3
J. Appl. Phys. 108, 103305 共2010兲
Kawamura, Lichtenberg, and Lieberman
(b)
2 1.5 1 0.5
1×10
0.005
0.01
1.5 1
0.015
y (m)
0.02
0
0.025
10×Ο −
O−
1×10
13
e− 0.01
0.015
y (m)
0.01
0.005
0.015
y (m)
0.02
0.025
0.02
0.025
O +2 O−
14
5×10
0.005
10×Ο +2
(d)
O2+
14
0
2
0.5
(c)
5×10
e−
2.5
Transverse T (V)
2.5
0
Transverse Density (m-3)
3
(a)
Axial Density (m-3)
Transverse Potential (V)
3
13
e− 0
0.025
0.05
0.075
0.1
x (m)
0.125
0.15
0.175
FIG. 1. PIC results for oxygen reaction rates and unmagnetized ions with Krecfac = 20, B0 = 600 G, L = 19 cm, R = 1.25 cm, p = 20 mTorr, and J0 = 1 A / m2.
was that the combination of the applied frequency of 13.56 MHz and the higher electron densities that were obtained at the higher current densities induced edge electric fields. This in turn produced higher electron temperatures near the transverse edges of the device. The lower current densities, resulting in lower electron densities, minimized this effect, and consequently, were used in most simulations. Because electrons move relatively freely along the field lines, large potentials build up at the device ends to keep the electron and ion flow ambipolar. Conducting transverse walls would result in potentials that inhibit negative ions from reaching them. Consequently, the transverse walls are insulated to obtain the desired type of transverse ion flows. The end walls could be either insulating or conducting. Both were tried in the simulations with some differences found due to the nonuniformity of the transverse electron density. We use insulated ends, which prevented nonzero net current densities flowing to the end walls, for the results of this section, but also used grounded metal end walls for other configurations. We found no significant difference in endloss flux between insulating and conducting end walls. B. PIC simulation results
Figure 1 shows PIC results for a 2D oxygen reactor with all insulating boundaries and unmagnetized ions. The discharge length L = 19 cm 共along B0兲 and transverse width 2R = 2.5 cm. Both the 1/4 cm dielectric in the transverse y direction and the 1/2 cm dielectric in the axial x direction have relative permittivity of four to separate the plasma from conducting boundaries. The pressure is p = 20 mTorr 共ng = 6.4⫻ 1020 m−3兲, the magnetic field is B0 = 600 G, and Krecfac = 20. The current profile has the form 共3兲 with J0 = 1 A / m2. The transverse potential is shown in Fig. 1共a兲,
transverse Te in Fig. 1共b兲, transverse densities in Fig. 1共c兲, and axial densities in Fig. 1共d兲. As expected, the potential profile, governed essentially by the ion energy, is nearly flat. The temperature Te is fairly constant in y, with slight increases near the transverse edges from the induced fields, as explained earlier. Both the electron density profile and the negative ion profile are nearly cosine functions. A selfconsistent determination of these profiles would be quite difficult, so we alternately accept them in our modeling. We note in this case that ␣0 ⬅ n−0 / ne0 ⬇ 4, and that ␣共y兲 ⬅ n−共y兲 / ne共y兲 ⬇ ␣0. Increasing the recombination to Krecfac = 100, with roughly the same electron density in the PIC simulation, the recombination has become significant. However, the only significant difference from Fig. 1 is a reduction in ␣0 to ␣0 ⬇ 3.0, due to the volume recombination. There are some other smaller differences which we will examine in the next section, after developing a model to compare to the simulations. If the ions are magnetized, this results in only a minor increase in O− due to slightly better confinement because the magnetizing term reduces the diffusion coefficient by only about a factor of two 关see 共7兲兴. Leaving the ions unmagnetized, but increasing the magnetic field by a factor of two to B0 = 1200 G reduces ␣0 somewhat due to some suppression of electron losses. The main effect is a very significant reduction in the radial electron loss, as will be seen in Sec. III. III. ANALYTIC MODEL AND COMPARISON WITH SIMULATIONS A. The model
We construct a 1D analytic model to compare to the PIC results. The key result that we take from the simulations is the observation that the transverse profiles of electrons and
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jap.aip.org/jap/copyright.jsp
J. Appl. Phys. 108, 103305 共2010兲
Kawamura, Lichtenberg, and Lieberman
with
L =
共4兲
冉 冊
d 2n e = 共iz − att兲ne − L共1 + ␣0兲ne . dy 2
共6兲
ⴱ is the transverse electron diffusion coefficient obHere, D⬜e tained from the complete ambipolar relations, as derived in ⴱ ⬇ D⬜e, the Appendix. Note that within a small factor D⬜e where for all species the transverse diffusion coefficients without ambipolar fields are11
De,i 2 2 , 1 + ce,i /e,i
De,i =
eTe,i . me,ie,i
共7兲
The solution to Eq. 共6兲 is obtained, using the observation from the PIC results that ne ⬇ 0 at y = ⫾ R 共note a shift in the y-coordinate from Fig. 1兲 ne共y兲 = ne0 cos
y , 2R
共8兲
such that Eq. 共6兲 yields ⴱ D⬜e
冉 冊 2R
2
ne0 = 共iz − att兲ne0 − L共1 + ␣0兲ne0 .
共9兲
Factoring out ne0, Eq. 共9兲 determines Te for given values of ␣0 and the other parameters. The negative ion balance equation is ⴱ − D⬜n
100Krecfac
0.2
d 2n − = attne − Krecn−共ne + n−兲, dy 2
共10兲
ⴱ with D⬜n calculated in the Appendix and ne as in Eq. 共8兲, the integration over the transverse dimension gives
4R ⴱ 2 ␣0ne0 = attne0 − Krecne0 ␣0共1 + ␣0兲R. D⬜n R
共11兲
Together with Eq. 共9兲, this gives Te and ␣0, for a given feedstock gas, in terms of device dimensions L and R, pressure p, the magnetic field B0, and the product Krecne0. For the model computations, we use a base case of L = 19 cm, R = 1.25 cm, and p = 20 mTorr 共ng = 6.4 ⫻ 1020 m−3兲, as in the PIC simulations. We magnetize the electrons with a 600 G axial magnetic field, and compare results with magnetized and unmagnetized ions, using oxygen parameters. We use the same cross-section set in the PIC
0.6 0.4
Irec /I⊥n Iz /I⊥p
0.2
20Krecfac 1
2
α0
3
4
5
0 0
1
2
α0
3
4
5
(c)
共5兲
In Eq. 共4兲, the flux and ion densities are taken to be surface- and volume-averages. Using Eq. 共5兲, approximating ni as ne共1 + ␣0兲, and recognizing that the electron and ion end losses must be equal if the sidewalls are insulated, the electron diffusion equation is
D⬜e,i =
0.6 0.4
(b)
I⊥e /I⊥p
0.8
1/2
.
0.8
(a)
0 0
2Di Te 1+ RL Ti
ⴱ − D⬜e
x 10
Current Ratios
2R ⌫z , Ln i ⬇ RL
−5
1
K recfac ne0 /n g
ions are quite similar, such that ␣0 ⬇ constant. An approximation taken from our previous work,8 is that the ion end loss can be related to the volume ion loss through an ion end loss frequency L
0.25
D (m 2/s)
103305-4
D⊥n ∗ D⊥n
0.2 0.15 0.1
D⊥p ∗ D⊥e
0.05 0 0
1
2
α0
3
4
5
FIG. 2. Theoretical results for oxygen reaction rates with Krecfac = 20 and 100, B0 = 600 G, L = 19 cm, R = 1.25 cm, p = 20 mTorr, and J0 = 1 A / m2.
simulations and the model. We plot results against the more physically intuitive quantity ␣0, rather than Krecfacne0 / ng. The two are uniquely related, as we see in Fig. 2共a兲. Note that the PIC simulation corresponding to this case has a density of 2.5⫻ 1013 m−3, such that the horizontal dashed and dotdashed lines correspond to PIC simulations with Krecfac = 20 and 100. They also correspond in the model to Krecfac = 1 and higher densities. Figure 2共b兲 gives important current ratios. Figure 2共c兲 gives the cross-field diffusion coefficients, showing that the ambipolar corrections to the diffusion coefficients are small 共see Appendix兲 and that the negative ion and electron diffusion coefficients are roughly equal. Because the recombination is relatively small compared to the negative ion flux to the wall, ␣0 does not change significantly for the two values of Krecfacne0 / ng, which differ by a factor of 5 关see Fig. 2共a兲兴. From the model equations, and from the results for the specific examples, we can draw some scaling relations which may be of use to designers. Using the approximation that the fluxes are uniform over their respective axial dimensions, the transverse electron and negative ion wall currents are obtained from Eqs. 共9兲 and 共11兲 as ⴱ L ne0 , I⬜e = 2D⬜e R
共12兲
ⴱ L ␣0ne0 . I⬜n = 2D⬜n R
共13兲
The axial loss current is obtained from Eq. 共4兲 as Iz = LRL共1 + ␣0兲ne0 .
共14兲
Substituting L from Eq. 共5兲 into Eq. 共14兲, we find
冉 冊
Iz = 2Di 1 +
Te Ti
1/2
共1 + ␣0兲ne0 .
共15兲
From the negative ion balance Eq. 共11兲, ne0 and ␣0 are related by
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jap.aip.org/jap/copyright.jsp
103305-5
J. Appl. Phys. 108, 103305 共2010兲
Kawamura, Lichtenberg, and Lieberman
ne0 =
ⴱ 4Ratt/ − D⬜n ␣0/R . Krec␣0共1 + ␣0兲R
共16兲
Equation 共16兲 determines ␣0 for any choice of ne0 关see Fig. 2共a兲兴. As can be seen from Eq. 共16兲, ␣0 must lie in the range ⴱ . The designer will generally want to 0 ⬍ ␣0 ⬍ 4attR2 / 2D⬜n achieve: 共a兲 high I⬜n 共high thrust兲, 共b兲 ne0 / ng not too small 共high propellant utilization兲, 共c兲 low Irec and low Iz 共high power efficiency兲, and 共d兲 low I⬜e 共small transverse electron effects兲, by proper choice of the system parameters R, L, mi, ng, B0, and ne0. From the above relations we identify the following relevant ratios in terms of physical parameters. From Eqs. 共12兲 and 共13兲, assuming electrons fully magnetized and ions essentially unmagnetized, we obtain ⴱ n2m1/2 1 I⬜e D⬜e ⬇ ⴱ ⬀ g2 i , I⬜n D⬜n ␣0 B 0␣ 0
共17兲
which is relevant to condition 共d兲. From Eq. 共11兲 or 共16兲, assuming reasonably high ␣0, we have Irec Krec␣0ne0R ⬀ Krec␣0R2nengm1/2 ⬇ i , ⴱ I⬜n D⬜n /R
共18兲
which is relevant to 共c兲. From Eqs. 共15兲 and 共13兲, assuming ⴱ ⬇ Di, we obtain Te / Ti Ⰷ 1 and D⬜n Iz I⬜n
⬇
冉冊
1 Te Ti
1/2
R , L
共19兲
which is also relevant to 共c兲. Finally, ␣0 can be solved from Eq. 共16兲 to give the limiting cases
␣0 ⬇
4R2 Kattng ⬀ R2Kattn2gm1/2 i , ⴱ 2 D⬜n
共20兲
for Irec / I⬜n small, and
␣0 ⬇
冉
4Kattng Krecne0
冊
1/2
,
共21兲
for Irec / I⬜n large. We note that the development in this section is effectively a global model once the approximation of ␣0 = const has been used. We will return to a global formalism in the more complicated geometries of subsequent sections.
TABLE I. Comparison of PIC with theory for oxygen reaction rates and magnetized ions with Krecfac = 20 and 100, B0 = 600 G, L = 19 cm, R = 1.25 cm, p = 20 mTorr, and J0 = 1 A / m2. 共a兲 Krecfac = 20
ne0 共m−3兲 Iatt I⬜e / I⬜n Irec / I⬜n I⬜n / Iatt ␣0 Te 共V兲 Iiz 共s−1兲 Iz / I⬜p ⴱ D⬜e 共m2 / s兲 ⴱ D⬜n 共m2 / s兲 L 共s−1兲
PIC
Theory
PIC
Theory
2.4⫻ 1013 1.1⫻ 1013 0.24 0.11 0.90 4.8 2.3 1.7⫻ 1013 0.38 0.25 0.29 1.5⫻ 103
2.4⫻ 1013 0.93⫻ 1013 0.21 0.11 0.90 3.9 2.1 1.5⫻ 1013 0.36 0.15 0.18 1.0⫻ 103
2.4⫻ 1013 1.0⫻ 1013 0.29 0.41 0.72 3.3 2.3 1.8⫻ 1013 0.43 0.19 0.20 2.0⫻ 103
2.4⫻ 1013 0.93⫻ 1013 0.28 0.46 0.68 3.0 2.1 1.4⫻ 1013 0.36 0.20 0.18 0.99⫻ 103
ues of L from the model to the simulations 共with both insulating and conducting end walls兲 indicates that the formula given in Eq. 共5兲 is reasonable. We also find that the ratio I⬜e / I⬜n can be reasonably well determined from the average diffusion coefficients and the assumption of a constant ␣0 I⬜e/I⬜n =
ⴱ D⬜e
ⴱ ␣0D⬜n
共22兲
.
Using the PIC results for the table in Eq. 共22兲 for case 1a and 1b, we obtain I⬜e / I⬜n = 0.18 and 0.29 respectively, in good agreement with the ratios in the table. Comparing case 1b to 1a, ␣0 has been significantly decreased both in the model and the simulation, as expected from increasing the recombination factor Krecfac. Irec / I⬜n has increased in the PIC results by a factor of 3.7, in good agreement with the factor of 3.4 increase from the scaling 共⬀Krec␣0兲 in Eq. 共18兲. In Table II, we present PIC results from varying important quantities which can be compared with the scalings from the model. In case 2a for unmagnetized ions, the theory preTABLE II. PIC results for oxygen reaction rates and both magnetized and unmagnetized ions with Krecfac = 20, R = 1.25 cm, p = 20 mTorr, and J0 = 1 A / m2. Unless otherwise noted, B0 = 600 G, L = 19 cm, and Kattfac = 1.
B. Comparison of model with PIC simulations
In Table I, we compare the results of two cases we have examined in Fig. 2 with the values obtained from PIC simulations. In the theory, we use the same density ne0 for the comparisons, but have multiplied the ion diffusion coefficients by a factor of 0.5 to make the ␣0’s similar in all comparisons. Such an adjustment could be due to profile differences between the model and simulation, or to the need for more accurate kinetic calculations of the diffusion. The results mainly speak for themselves, but we point out a few salient features. We see good agreement between PIC and theory for the important ratio of electron to ion radial flux I⬜e / I⬜n. One of the motivations of the simulations is to verify that the formula for L, which involves a number of approximations, is reasonably accurate. Comparing the val-
共b兲 Krecfac = 100
Unmag. ions
ne0 共m−3兲 Iatt I⬜e / I⬜n Irec / I⬜n I⬜n / Iatt ␣0 Te 共V兲 Iiz 共s−1兲 Iz / I⬜p ⴱ D⬜e 共m2 / s兲 ⴱ D⬜n 共m2 / s兲 L 共s−1兲
Mag. ions
共a兲 B0 = 600 G
共b兲 B0 = 1200 G
共c兲 L = 9.5 cm
共d兲 Kattfac = 3
2.5⫻ 1013 1.1⫻ 1013 0.23 0.077 0.92 3.9 2.3 1.7⫻ 1013 0.37 0.31 0.36 1.7⫻ 103
2.7⫻ 1013 1.1⫻ 1013 0.052 0.063 0.95 3.2 2.0 1.6⫻ 1013 0.40 0.05 0.28 1.8⫻ 103
2.1⫻ 1013 4.9⫻ 1012 0.27 0.14 0.89 7.5 2.8 1.2⫻ 1013 1.0 0.30 0.21 4.1⫻ 103
2.1⫻ 1013 3.2⫻ 1013 0.10 0.23 0.79 13 2.8 4.1⫻ 1013 0.26 0.24 0.29 1.0⫻ 103
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jap.aip.org/jap/copyright.jsp
103305-6
J. Appl. Phys. 108, 103305 共2010兲
Kawamura, Lichtenberg, and Lieberman 4
Transverse T (V)
1×10
e−
(a)
15
(b)
3
8×10
2
n (m -3)
6×10
1
0 0
0.005
0.01
10×I −
4×10
10×I +
2×10
0.015
y (m)
α
14
1000
I+ I−
14
α
14
100
14
10×e − 0.02
0.025
0.03
0 0
0.005
0.01
0.015
y (m)
0.02
0.025
10 0.03
FIG. 3. PIC results of particle 共a兲 temperature 共eV兲 and 共b兲 density profiles for iodinelike discharge and axially confined heating currents with Krecfac = 20, Kattfac = 5 共no threshold兲, B0 = 600 G, L = 19 cm, R = 1.25 cm, p = 20 mTorr, and J0 = 1 A / m2. The ␣ values are shown with a logarithmic scale.
dicts an increase in D⬜n by a factor of two from case 1a, the magnetized-ion case in Table I. The increase in the PIC diffusion coefficient is somewhat less, indicating other factors also enter, but rather weakly. In case 2b we increase the magnetic field by a factor of two, from the base case, 共but keeping the ions unmagnetized to correspond to a heavier mass in which 2ci / 2ci Ⰶ 1兲. Comparing case 2b with case 2a, we see that the ratio I⬜e / I⬜n is decreased approximately by the factor of 4 predicted from the scaling in Eq. 共22兲. The electron profile in this case lies between the model cos共y / 2R兲 profile and the cos2共y / 2R兲 profile which produced the special solution. In case 2c, with magnetized ions, we reduce the length by a factor of two to L = 10 cm 共9.5 cm when the dielectrics are subtracted兲. From 共4兲, the theory predicts that L is increased by approximately a factor of 2 from the value in case 1a. In the PIC simulation, the increase in L is by a larger factor of 2.7, indicating the importance of the end effects for short systems. Note that for case 2c, the axial positive ion loss is approximately equal to the radial positive ion loss 共Iz / I⬜p ⬇ 1兲, such that power would be used rather inefficiently. In case 2d, we investigate the effect of a larger attachment coefficient by increasing Katt by a factor of three from its value in case 1a. Since the recombination is small compared to the attachment 共less than 20% in this case兲 we expect, from 共11兲, that ␣0 would increase essentially by the same factor of three, which we confirm from the simulation by comparing the value in case 2d with that in case 1a. C. Discussion of uniform Te results
From the previous PIC simulations and accompanying analysis, we have found that it is quite difficult to achieve a result approximating the special solution 共1兲, for which the electron flux at the radial walls approaches zero, leaving equal numbers of positive and negative ions to be extracted. This solution can be approached at high magnetic field strengths, provided a heavy ion mass is employed such that the ions are not strongly magnetized. However, high magnetic fields 共B0 = 1200 G in the small device we considered兲 may not be practical for a particular application. In the PIC calculations and accompanying model, we have considered relatively uniform transverse electron temperatures 共exactly uniform in the model兲. Such a distribution is not necessary if the electron transport is significantly inhibited across the field. Achieving cooler edge electrons requires both a struc-
ture of the heating fields that heat primarily away from transverse walls, and also an appropriate geometry, magnetic field strength, and other parameters favoring axial electron loss over radial electron loss. Another factor that enters into the effectiveness of reducing the electron flux to the wall is the energy dependence of the attachment cross-section. If attachment is efficient at low energies, then a reduction in temperature near the transverse walls will locally favor attachment over ionization, and, therefore, deplete the electron population. Oxygen has an energy threshold for attachment of 4.2 eV, and thus at the generally fairly low electron energies existing in magnetically confined devices does not efficiently accomplish attachment electron absorption. IV. NONUNIFORM Te AND INCREASED Katt
The above discussion suggests studying a feedstock gas resembling iodine rather than oxygen, and heating the discharge nonuniformly, such that the edge electrons are cooler than those near the axis. Iodine has not been well characterized, but some information on attachment is available. We use mostly the well-known oxygen cross-sections except that we increase the attachment cross-section by a factor of 5 and remove the attachment threshold, as suggested by the literature.12 As in Sec. III, we account for the lower PIC density by using an enhanced recombination factor Krecfac = 20. The heavier mass of iodine is used, such that at the base case of B0 = 600 G, the ions are essentially unmagnetized but also diffuse more slowly. For the previous geometry to primarily heat the electrons away from the transverse walls, we use for 兩y兩 ⬍ R / 2 a transverse current profile J共y兲 = J0 cos2共y / R兲, with J0 = 1 A / m2 and R = 1.25 cm, with J共y兲 = 0 for 兩y兩 ⬎ R / 2. In Fig. 3共a兲, we give Te共y兲, and in Fig. 3共b兲 the densities for nI+ 共solid兲, nI− 共dashed兲 and 10ne 共dotted兲, to graph it on the same scale as the ion densities. The ␣ = nI− / ne values 共dot-dashed兲 are also shown in Fig. 3共b兲 but with a logarithmic scale. The central ␣0 ⬇ 27 is approximately a factor of 5 higher than for the oxygen case, reflecting the increase in Katt by 5, and the absence of a cutoff. The lowering of Te, mainly outside of the central half of the device, where the current is applied, has greatly decreased the ionization in these outer regions, so decreasing the electron density. We find the edge ␣共R兲 ⬇ 2000, which is quite high, resulting in I⬜e / I⬜n ⬇ 4 ⫻ 10−4, i.e., an ion-ion plasma. As one also notes, the electron density profile now approximates the electron profile
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jap.aip.org/jap/copyright.jsp
103305-7
J. Appl. Phys. 108, 103305 共2010兲
Kawamura, Lichtenberg, and Lieberman
with negligible wall electrons, studied previously,8 but is due to decreasing Te near the walls, rather than to a special uniform Te solution. Introducing a threshold energy of 4.2 eV for the attachment leads to a modest decrease in the central electronegativity to ␣0 ⬇ 15. The major change is that the edge electron density does not fall to such a low value. The resulting edge ␣共R兲 ⬇ 200 and ratio I⬜e / I⬜n ⬇ 1.5⫻ 10−3 are still in a useful range. The PIC iodinelike results show a high ionization frequency iz inside an active core of thickness 2R1 ⬇ 1.25 cm. In the edge region iz Ⰶ att, such that a tworegion diffusion model with no ionization in the edge region can be used to determine the equilibrium parameters. In the cases considered, the recombination is small compared to the attachment, so we also neglect recombination to make the analytic model more tractable. The electron and negative ion balance relations are then described by −
2 ⴱ d ne D⬜e 2
dy
ⴱ − D⬜n
= 共iz − att − L兲ne − Lnn ,
d 2n n = attne , dy 2
共23兲
共24兲
where we approximate iz to be constant in the core and zero in the edge region. The diffusion coefficients are given in Appendix. Differentiating Eq. 共23兲 twice and substituting Eq. 共24兲 to eliminate nn, we obtain 4 d2ne Latt ⴱ d ne − D⬜e = 共 − − 兲 + ⴱ ne . iz att L dy 4 dy 2 D⬜n
共25兲
In the core we can substitute a symmetric eigenfunction cos ky into Eq. 共25兲 to obtain the characteristic equation ⴱ 4 k − 共iz − att − L兲k2 + D⬜e
Latt = 0, ⴱ D⬜n
共26兲
which has two real positive roots k1 and k2 for 共iz − att ⴱ ⴱ − L兲2 ⬎ 4D⬜e Latt / D⬜n . The core solution can therefore be written as ne = A1 cos k1y + A2 cos k2y.
共27兲
In the edge, with iz = 0, we substitute the eigenfunctions either sinh 共R0 − y兲 or cosh 共R0 − y兲 into Eq. 共25兲 to obtain ⴱ 4 − 共att + L兲2 + D⬜e
Latt = 0, ⴱ D⬜n
共28兲
which yields two positive real roots 1 and 2 for 共att ⴱ ⴱ + L兲2 ⬎ 4D⬜e Latt / D⬜n . Guided by the PIC results, we use an edge solution for ne as in Eq. 共1兲 having both ne = 0 and dne / dy = 0 at the outer edge to obtain ne = A3关sinh 1共R0 − y兲 − 共1/2兲sinh 2共R0 − y兲兴 + A4关cosh 1共R0 − y兲 − cosh 2共R0 − y兲兴.
共29兲
Matching the core and edge ne’s and their first three derivatives at y = R1, we obtain four linear equations M · A = 0, where A is a column vector of the four unknown coefficients A1, A2, A3, and A4, and M is the corresponding 4 ⫻ 4 coeffi-
1×10 8×10 6×10
15
I+ I−
14
α 1000
14
n (m -3)
α
14
100
4×10 2×10
14
0 0
10×e − 0.005
0.01
0.015
y (m)
0.02
0.025
10 0.03
FIG. 4. Profiles of the ion densities nI+ 共solid兲, nI− 共dashed兲, and 10⫻ ne 共dotted兲 for the two-region model. The ␣ values 共dot-dashed兲 are shown with a logarithmic scale.
cient matrix. Setting det M = 0 yields the eigenvalue Te, the core electron temperature 共as well as all temperatureⴱ , etc.兲. The ratios of the A’s dependent quantities iz , att , D⬜e are then determined by the ratios of the cofactors of any row of M. Hence, ne is determined up to a constant factor. With ne共y兲 known, ni can be determined using the positive ion balance relation ⴱ − D⬜i
d 2n i = izne − Lni , dy 2
共30兲
again with iz a 共now-known兲 constant in the core and zero in the edge region. The solution to 共30兲 in the core is ni = B1 cos k1y + B2 cos k2y + B3 cosh Ly, ⴱ 2 k1,2 + L兲 B1,2 = A1,2iz / 共D⬜i
共31兲
ⴱ 1/2 L = 共L / D⬜i 兲 .
and The where first two terms are a particular solution of Eq. 共30兲 in the core and the third term is the symmetric homogenous solution. In the edge region, again guided by the PIC solution, we choose ⴱ dni / dy = niuBi at the outer the boundary condition that −D⬜i edge y = R0, where uBi = 共eTi / mi兲1/2 is approximately the Bohm velocity at the 共highly electronegative兲 edge. Then from Eq. 共30兲 with iz = 0, we obtain ⴱ ni = B4关sinh L共R0 − y兲 + 共D⬜i L/uBi兲cosh L共R0 − y兲兴.
共32兲 Matching the core and edge densities Eqs. 共31兲 and 共32兲, and their first derivatives, at y = R1, we obtain two linear equations to determine the two unknown coefficients B3 and B4. Hence, ne and ni 共as well as nn by quasineutrality兲 are determined up to a constant factor. Figure 4 shows the calculated density profiles versus transverse coordinate y for 10⫻ ne 共dotted兲, nI+ 共solid兲, and nI− 共dashed兲 from this two-region model for the iodinelike discharge shown in Fig. 3. The simulation values of R1 ⬇ 0.5 cm and the central electron density ne0 are used in the model. The results are plotted with the same scales as the PIC results in Fig. 3共b兲 to facilitate a comparison. There is reasonable agreement for both the shapes of the transverse density and ␣ profiles, and the central values of ␣. The shapes of the ion and electron radial profiles, both for the PIC simulation and for the two-region theoretical
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jap.aip.org/jap/copyright.jsp
103305-8
J. Appl. Phys. 108, 103305 共2010兲
Kawamura, Lichtenberg, and Lieberman
y (cm)
8 1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
5 mm
1
1
6 4
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
6 cm
1
1 1
1
1
1
1
1
1
1
1
1
5
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
4 cm
2 0
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1 1
1 1
1
15
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
6 cm 1
1
1
1
1
1
1 1
10
1
1
1
1
1
1
1
20
1
1
1
1
25
1
x (cm)
30
FIG. 5. Configuration with source region in center connected to two wider downstream chambers.
model, indicate that a global model can be constructed which will give reasonably accurately the important parameters, and furthermore, can be used for scaling. As will be done for a more complicated geometry in the following section, we use a cosine profile for the electrons over the central region, and a cosine profile for the negative ions over the entire diameter. The basic equations are the electron balance relation
iz = att + L共1 + ␣0兲
R0 R1
共33兲
and the negative ion balance
att =
2 R0 ⴱ Krecne0␣0共1 + ␣0兲 + D⬜n ␣0 , 4 R1 4R0R1
共34兲
where ␣0 is the on-axis electronegativity. Comparisons with the two-region model give good agreement, indicating that global models can be constructed in geometries with a hot electron core and a cold electron periphery.
V. MORE COMPLICATED GEOMETRIES
In experiments,6,7 an ion-ion plasma was produced near a boundary by having a step in the diameter with the main ionization source for electrons 共the heating region兲 confined to the narrower region. The cooler electron temperature near the outer part of the wider downstream region removed electrons by attachment, as described in Sec. IV. In the 2D rectangular approximation the PIC boundary conditions can readily be modified to a step geometry with the source heating confined to the narrower region. The configuration is shown in Fig. 5. For the cases investigated, the dielectric transverse 共radial兲 walls in the central source region were separated by 4 cm, and the length of the region was 10 cm. Symmetric downstream regions connected to the source region were each also 10 cm long, with grounded conducting ends. The transverse 共radial兲 walls were dielectric, backed by grounded conductors. Downstream conductor separations of 7 and 10 cm were investigated in simulations, with 1/2 cm dielectrics inside the conductors. The heating by an axial current in the central source region was a cosine in both directions, having the form Jx = J0 cos共y / 2R0兲cos共x / L0兲, where R0 is the half-width of the central region and L0 the length of the heating zone. As previously, J0 = 1 A / m2, x is the axial direction of the applied magnetic field B0 = 600 G, and y is the cross-field direction. As in Sec. IV, an iodinelike
discharge is simulated by using the base oxygen crosssection set, but increasing the attachment cross-section by a factor of five and removing the threshold. In Fig. 6, the results are given for what we shall call a base case, with Krecfac = 20 and p = 20 mTorr. The metal walls in the downstream regions are separated by 7 cm. Figures 6共a兲 and 6共b兲 give the radial densities at the center of the source region and the center of one downstream end region, respectively. We see that the electrons are well confined in the end region, while the positive and negative ion densities fall-off to the outer walls. In Figs. 6共c兲 and 6共d兲, the axial potential and axial densities, respectively, are shown. The potential fall in the end regions of approximately ⌬V ⬇ 1.5 V leads to rather large ion density decreases due to their low temperatures, but very little change in the on-axis density of the considerably higher temperature electrons. In Fig. 6共e兲, the axial Te is given, and in Fig. 6共f兲 the corresponding axial reaction rates Kiz and Katt. We note that the axial potential fall-off 关in Fig. 6共c兲兴 corresponds mainly to a decreasing electron temperature away from the heating zone 关in Fig. 6共e兲兴. The strong drop in Kiz in the end regions corresponds to the drop in Te, while the axial Katt is relatively unchanged. For the development of a theoretical treatment in Sec. VI, we note that the potential variation, which is a significant fraction of Te, leads to electric field driven positive ions flowing downstream and negative ions flowing upstream. The simulation results indicate that the particle balances in the central and end sections are coupled by mobilitydriven axial 共x兲 fluxes of positive and negative ions. At the junction between the central and right-hand end sections, ⌫i = i0niE and ⌫n = −n0nnE, where the electric field is estimated as E = Te0 / lE, with lE the scale length of a potential with characteristic value Te0. By flux neutrality, the corresponding electron flux is ⌫i − ⌫n. Accounting for these fluxes and with reasonable choices of the transverse density profiles, a simple global 共volume-averaged兲 model of the particle balances in the central and end sections can be used to describe the equilibrium of the coupled system. We let R0 and R1, and L0 and L1 be the transverse half-lengths and axial lengths of the central and end sections, and assume that all central cell transverse density profiles vary as cos共y / 2R0兲. Then for a unit length in the ignorable 共z兲 direction, the electron balance in the central section can be expressed as
iz0 = att0 + 2关i0共1 + ␣0兲 + n0␣0兴
2 Te0 ⴱ + D⬜e0 , l EL 0 4R20 共35兲
where ␣0 = nn / ne is the 共constant兲 electronegativity and the subscripts “0” denote quantities evaluated at the central cell temperatures. The last term gives the transverse electron loss to the insulating walls and the factor of 2 multiplying the bracket indicates a symmetric system with positive ions leaving and negative ions entering from both ends. The central negative ion balance relation can be similarly written as
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jap.aip.org/jap/copyright.jsp
J. Appl. Phys. 108, 103305 共2010兲
Kawamura, Lichtenberg, and Lieberman
1×10
15
I+ I−
(a)
End Cell Transverse n (m-3)
Center Cell Transverse n (m-3)
103305-9
6×10
4×10
5×10
14
2×10
20×e − 0 0
0.02
0.04
I+ I−
(b)
14
20×e −
14
0 0
0.06
y (m)
14
0.01
0.02
0.03
15
Axial Density n (m -3)
Axial Potential (V)
(c) 10
5
0 0
0.05
0.1
0.15
x (m)
0.2
1×10
5×10
15
5
8×10
Axial K (m 3 /s)
Axial T (V)
6×10
e−
4×10
2 1 0 0
2×10
5×I − 5×I + 0.05
0.1
0.15
x (m)
0.2
0.05
0.1
0.15
x (m)
0.2
0.25
-16
iz
-16
-16
att -16
0 0
0.25
0.07
20×e −
(f)
3
0.06
14
(e) 4
0.05
I+ I−
(d)
0 0
0.25
0.04
y (m)
0.05
0.1
0.15
x (m)
0.2
0.25
FIG. 6. PIC results for iodinelike discharge in configuration of Fig. 5 with Krecfac = 20, B0 = 600 G, p = 20 mTorr, and J0 = 1 A / m2.
att0 + 2n0␣0
Te0 2 ⴱ = Krec0ne0␣0共1 + ␣0兲 + D⬜n0 ␣0 2 , l EL 0 4 4R0 共36兲
where the last two terms give the volume recombination loss and the transverse negative ion loss. In a downstream end section, guided by the PIC simulation results, we choose the same on-axis density ne0 and transverse electron profile cos共y / 2R0兲 as in the central cell. The PIC results also indicate a negative ion profile that extends over the entire transverse dimension of the end cell. We model this using a simple cos共y / 2R1兲 transverse distribution and discuss this further in the conclusions. Since the electrons are confined to the smaller transverse dimension, there is a negligible transverse electron loss. However, there is a considerable axial loss, equal to the positive ion axial loss, which we model using the axial loss frequency L given from Eq. 共5兲
冉
L = Di1 1 +
Te1 Ti1
冊
1/2
2 , R 1L 1
共37兲
where we take the losses to be over the full radius on both sides of each end of the downstream sections, and the sub-
script 1 refers to downstream quantities. Equation 共37兲 somewhat overestimates L for this geometry. The electron balance relation in one end cell is
iz1 + 关i0共1 + ␣0兲 + n0␣0兴
Te0 R1 = att1 + L共1 + ␣1兲 , l EL 1 R0 共38兲
where ␣1 = nn共y = 0兲 / ne共y = 0兲 is the on-axis electronegativity. The negative ion balance in an end cell can be written as
att1 = n0␣0
Te0 R1 + Krec1ne0␣1共1 + ␣1兲 l EL 1 4 R0
ⴱ + D⬜n1 ␣1
2 . 4R0R1
共39兲
For a given electron density ne0 and specifying the heavy particle temperature as a fixed ratio of the electron temperature in each of the sections, then Eqs. 共35兲–共39兲 can be solved numerically to determine Te0, Te1, ␣0, and ␣1. Once these quantities are known, then, all of the particle generation and loss terms in the central and end sections can be evaluated.
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jap.aip.org/jap/copyright.jsp
103305-10
J. Appl. Phys. 108, 103305 共2010兲
Kawamura, Lichtenberg, and Lieberman
TABLE III. Comparison of PIC with theory for iodinelike discharge in configuration of Fig. 5 with p = 20 or 10 mTorr, Krecfac = 20, B0 = 600 G, and J 0 = 1 A / m 2. 20 mTorr Central cell ne0 共m−3兲 Iatt I⬜n / Iatt Irec / I⬜n Inin / I⬜n ␣0 Te 共V兲
PIC
Theory
10 mTorr PIC
Theory
End cells
1.0⫻ 1013 1.0⫻ 1013 5.5⫻ 1014 5.0⫻ 1014 0.68 0.62 1.4 1.3 0.73 0.43 105 98 2.8 3.6 20 mTorr PIC Theory
1.0⫻ 1013 1.0⫻ 1013 3.1⫻ 1014 2.5⫻ 1014 1.1 1.2 0.44 0.34 0.62 0.34 54 45 2.9 3.6 10 mTorr PIC Theory
Iatt1 I⬜n1 / Iatt1 Irec1 / I⬜n1 ␣01 Te1 共V兲 Iz / I⬜p1 L 共s−1兲
9.8⫻ 1014 0.13 3.5 58 1.8 1.4 270
6.7⫻ 1014 0.32 1.2 37 2.1 1.4 660
12⫻ 1014 0.27 1.9 62 3.1 0.79 750
5.0⫻ 1014 0.52 0.52 31 3.3 0.38 1540
The negative ion balance equations are insensitive to the electron temperature so that approximate negative ion results can be obtained, assuming nominal temperatures, without solving for Te0 and Te1 from the electron balance equations. However, the scale length lE is sensitive to the source and downstream temperatures, so approximate average values, taken from a few PIC cases, are used. In Table III, we compare PIC with theory for negative ion parameters, in two typical cases. The axial source electron density ne0, taken from the PIC results, is specified in the theory, for the comparison. A nominal value of lE = 0.18 m consistent with simulations at p = 20 mTorr is used in the theory at this pressure. For p = 10 mTorr, the drop in downstream Te is smaller, so that the effective lE is longer. For this pressure, we use lE = 0.36 m. The results in Table III give generally reasonable agreement between PIC simulations and theory. As with the simpler geometries the theoretical L factors work quite well. The approximate reduction in the ratio Irec / I⬜n from p = 20 mTorr to p = 10 mTorr, is consistent with the general scaling relations from Sec. III. The weaker expected drops in ␣0 are also found. In downstream regions the radial electron current is negligible, and therefore has not been listed. We have also investigated 10 cm width end regions, but the main effect was an increase in the overall recombination flux, which is not desirable, so we omit details of the results. VI. CONCLUSIONS AND FURTHER DISCUSSION
We have used a 2D PIC code to investigate the equilibrium of a device with magnetized electrons but only weakly magnetized ions, in an electronegative plasma. The purpose of the simulations is to study the possibility of extracting positive and negative ions in nearly equal numbers. In Secs. II and III, to conform to experiments and previous analysis, an oxygen feedstock gas was used, together with a heating
current that produced a relatively constant Te across the magnetic field lines. A device of 19 cm length and 2.5 cm transverse dimension made the ratio of end loss to transverse loss of the positive ions relatively small, allowing a 1D theoretical approximation to be made. The walls were insulating as is required to extract negative ions. A pressure p = 20 mTorr and an axial magnetic field B0 = 600 G were used in most of the simulations. At these parameters, the ions were significantly magnetized; therefore, to investigate the effect of ion magnetization, simulations were also performed with the ions unmagnetized. Because the PIC simulations require low electron densities, and the relative recombination rate is determined by the product Krecne, an enhanced Krec was used to rescale the recombination to normal values, designated as a Krecfac. For oxygen, at p = 20 mTorr, setting this factor to 20 gave relatively low recombination, while a factor of 100 gave significant recombination. Both were used in the simulations. A current directed along the magnetic field with a transverse cosine profile resulted in a relatively constant transverse Te. The results of the simulations indicated that the electronegativity ␣ = n− / ne was reasonably constant transversely. For fairly small ␣ this resulted in a significant ratio of electron to ion current I⬜e / I⬜n at the transverse wall, which we desired to be small. Higher magnetic fields were found to proportionally reduce this ratio. Because ␣ ⬇ ␣0 a simple model having this property could be developed and results of the model were in good agreement with the simulations. The model also allowed various scaling with parameters to be explored analytically. This, in turn, suggested the PIC parameters to be explored. The results from the simulations and the accompanying model indicated that a significant reduction in the ratio of electron to ion fluxes would be obtained with higher mass ions that suppressed ion magnetization, with higher ␣’s that would result from higher Katt, and with lower transverse edge electron temperature. The first two improvements would be achieved with a feedstock gas such as iodine. The third and probably most important effect can be obtained by limiting the plasma excitation and heating to the central portion of the transverse dimension. In Sec. IV, based on these insights, PIC simulations were performed with more iodinelike parameters of background neutrals, and with the excitation concentrated in the central half of the transverse dimension. The results of two representative simulations indicated excellent suppression of transverse electron flux at a pressure of p = 20 mTorr and a modest magnetic field strength of B0 = 600 G. A two-region model was developed that well- approximated the PIC results. Although parameter scaling can easily be explored, numerically, with the model, there are more parameters to study, such that more extensive PIC simulations would be required to validate the model over a wide range of parameters. Consequently, a global model was developed from which scaling could be obtained with results that were consistent with the two-region model. In Sec. V, a final set of simulations were performed to conform more closely to an experimental geometry which had a source region and wider downstream regions. If the electrons are significantly better confined, in the transverse direction, than the ions, this leads to cooler electrons at trans-
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jap.aip.org/jap/copyright.jsp
103305-11
J. Appl. Phys. 108, 103305 共2010兲
Kawamura, Lichtenberg, and Lieberman
verse positions beyond the transverse boundaries of the source region, and therefore naturally favors the production of ion-ion edge plasmas. This was, in fact, achieved, and coupled source-downstream global models were constructed, which gave good agreement with the PIC simulations. One of the results from all of the modeling was comparisons of the classical diffusion, both radially and axially, to the PIC results. We found that the PIC radial diffusion was approximately classical for the ion species. The axial diffusion coefficients controlled the L’s, which were found to agree reasonably well between theoretical calculations and PIC simulations.
dnn dne + 共D⬜i − D⬜e兲 dy dy . 共⬜i + ⬜n兲nn + 共⬜i + ⬜e兲ne
共D⬜i − D⬜n兲 E=
For a constant transverse electronegativity, nn共y兲 = ␣0ne共y兲,
共A6兲
we substitute Eq. 共A6兲 into Eq. 共A5兲 to eliminate the dependence of E on dne / dy to obtain n nE =
␣0共D⬜i − D⬜n兲 + D⬜i − D⬜e dnn . ␣0共⬜i + ⬜n兲 + ⬜i + ⬜e dy
ⴱ ⌫⬜n = − D⬜n
dnn , dy
ⴱ = D⬜n + ⬜n D⬜n
共A1兲
where the ⌫’s are given by the usual drift-diffusion relations ⌫⬜i = − D⬜i
dni + ⬜iniE, dy
共A2兲
⌫⬜n = − D⬜n
dnn − ⬜nnnE, dy
共A3兲
⌫⬜e = − D⬜e
dne − ⬜eneE. dy
共A4兲
Substituting Eqs. 共A2兲–共A4兲 into Eq. 共A1兲, with ni = nn + ne from quasineutrality, and solving for the electric field, we obtain
␣0共D⬜i − D⬜n兲 + D⬜i − D⬜e . ␣0共⬜i + ⬜n兲 + ⬜i + ⬜e
共A9兲
A similar substitution to eliminate the dependence of E on dnn / dy yields
␣0共D⬜i − D⬜n兲 + D⬜i − D⬜e . ␣0共⬜i + ⬜n兲 + ⬜i + ⬜e
共A10兲
A similar derivation for the positive ions yields
The effective cross-field diffusion coefficients are found from the ambipolar flux balance condition ⌫⬜i = ⌫⬜n + ⌫⬜e ,
共A8兲
then comparing Eq. 共A8兲 to Eq. 共A3兲, we obtain
ⴱ = D⬜e + ⬜e D⬜e
APPENDIX: EFFECTIVE CROSS-FIELD DIFFUSION COEFFICIENTS
共A7兲
Substituting Eq. 共A7兲 into Eq. 共A3兲 and introducing the efⴱ by fective diffusion coefficient D⬜n
ACKNOWLEDGMENTS
E.K., A.J.L., and M.A.L. acknowledge the support of California Industries and University of California Discovery Grant No. ele07-10283 under the IMPACT program and the Department of Energy Office of Fusion Energy Science Contract No. DE-SC0001939. The work is closely related to our previous collaboration with Gary Leray and Pascal Chabert at Ecole Polytechnique, France.
共A5兲
ⴱ = D⬜i − ⬜i D⬜i
␣0共D⬜i − D⬜n兲 + D⬜i − D⬜e . ␣0共⬜i + ⬜n兲 + ⬜i + ⬜e
1
共A11兲
S. K. Kanakasabapathy, L. J. Overzet, V. Midha, and D. J. Economou, Appl. Phys. Lett. 78, 22 共2001兲. A. Meige, G. Leray, J.-L. Raimbault, and P. Chabert, Appl. Phys. Lett. 92, 061501 共2008兲. 3 R. Kawai and T. Mieno, Jpn. J. Appl. Phys., Part 2 36, L1123 共1997兲. 4 P. Chabert, T. E. Sheridan, R. W. Boswell, and J. Perrin, Plasma Sources Sci. Technol. 8, 561 共1999兲. 5 S. G. Walton, D. Leonhardt, R. F. Fernsler, and R. A. Meger, Appl. Phys. Lett. 81, 987 共2002兲. 6 C. S. Corr, N. Plihon, and P. Chabert, J. Appl. Phys. 99, 103302 共2006兲. 7 A. Aanesland, A. Meige, and P. Chabert, J. Phys.: Conf. Ser. 162, 012009 共2009兲. 8 G. Leray, P. Chabert, A. J. Lichtenberg, and M. A. Lieberman, J. Phys. D: Appl. Phys. 42, 194020 共2009兲. 9 V. Vahedi and G. DiPeso, J. Comput. Phys. 131, 149 共1997兲. 10 V. Vahedi and M. Surendra, Comput. Phys. Commun. 87, 179 共1995兲. 11 M. A. Lieberman and A. J. Lichtenberg, Principles of Plasma Discharges and Materials Processing, 2nd ed. 共Wiley-Interscience, New York, 2005兲. 12 G. E. Caledonia, Chem. Rev. 75, 333 共1975兲. 2
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jap.aip.org/jap/copyright.jsp
Two-dimensional particle-in-cell simulations of transport in a magnetized electronegative plasma E. Kawamura,a兲 A. J. Lichtenberg, and M. A. Lieberman Department of Electrical Engineering, University of California, Berkeley, California 94720-1770, USA
共Received 27 July 2010; accepted 27 September 2010; published online 19 November 2010兲 Particle transport in a uniformly magnetized electronegative plasma is studied in two-dimensional 共2D兲 geometry with insulating 共dielectric兲 boundaries. A 2D particle-in-cell 共PIC兲 code is employed, with the results compared to analytic one-dimensional models that approximate the end losses as volume losses. A modified oxygen reaction set is used to scale to the low densities used in PIC codes and also to approximately model other gases. The principal study is the limiting of the transverse electron flow due to strong electron magnetization. The plasma in the PIC calculation is maintained by axial currents that vary across the transverse dimension. For a cosine current profile nearly uniform electron temperature is obtained, which at the B-fields studied 共600–1200 G兲 give a small but significant fraction 共0.25 or less兲 of electron to negative ion transverse loss. For a more transverse-confined current, and approximating the higher mass and attachment reaction rate of iodine, the fraction of electron to negative ion transverse loss can be made very small. The models which have been constructed reasonably approximate the PIC results and indicate that the cross-field transport is nearly classical. © 2010 American Institute of Physics. 关doi:10.1063/1.3506518兴 I. INTRODUCTION
Ion-ion plasmas composed of positive and negative ions, with a negligible fraction of electrons, have potential applications for microelectronics etching1 and plasma thrusters for space missions.2 A method of producing such plasmas, in which the ion-ion plasma is formed at the periphery of a magnetized electronegative plasma core, has been explored in several experiments.3–7 Recently, a one-dimensional 共1D兲 fluid model, together with an analytic theory, was developed8 to approximate the configuration of one of these experiments.6,7 In that model, the two-dimensional 共2D兲 nature of the device was approximated, by including a volume axial loss frequency L that corresponds to the end losses. Rectangular, rather than cylindrical, coordinates were used, for simplicity, and the value of the magnetic field B0 was chosen such that the electrons were strongly magnetized, while the ions were only weakly magnetized. To make the numerical integration tractable, the electron temperature Te was assumed uniform. One then looked for special solutions, at a transverse boundary y = R, for which ne共R兲 = 0,
冏 冏 dne dy
= 0,
共1兲
y=R
as a function of Te, the ratio of central electron density to neutral density ne0 / ng, and the ratio of central negative ion to electron density ␣0 = nn0 / ne0. The boundary conditions imply that a pure ion-ion plasma exists at y = R. Solutions were found in the ne0 / ng − ␣0 plane for various values of R and length L, holding RL = constant. The latter constraint was useful in the numerics because the approximation used for L was a function of that constraint. The various assumptions a兲
Electronic mail: [email protected].
0021-8979/2010/108共10兲/103305/11/$30.00
that were made in that work remain open to investigation. The solutions do not imply the existence of a physical transverse wall at the position for which the specified electron boundary conditions are satisfied. Furthermore, at a physical wall where the particular boundary conditions are satisfied, it does not ensure that the special solution, rather than another solution at a different Te would be favored. It is also necessary to determine the validity of the expression used for L, and the conditions under which Te is relatively uniform over the cross-section. To examine the approximations of the numerical calculations and accompanying model, and to investigate other configurations of interest, a 2D particle-in-cell 共PIC兲 code is employed using an oxygenlike feedstock gas. PIC parameters are scaled to higher density, using a rescaled oxygen cross-section set, and the dimensions and pressure are adjusted to keep pd = const. Since the pressure is adjusted, the magnetic field is also adjusted to preserve the relative strength of the magnetization. This will become clearer from the analytic models which we develop to compare to the PIC simulations. In Sec. II, we describe the 2D PIC code and present some typical PIC simulation results for a uniform structure with nearly uniform Te, related to our previous work. In Sec. III, a model is developed to approximate the electron and ion transverse profiles observed in the PIC simulations, and the model results are compared with PIC simulations over a range of parameters. Some implications of various configurations for producing ion-ion plasmas are discussed. In Sec. IV, we apply these ideas to iodinelike parameters with a modified heating configuration that produces nearly electronfree ion-ion edge plasma. A two-region model is developed and compared to the PIC simulation. In Sec. V, a more complicated, but more practical, geometry is introduced, with PIC simulations showing that an edge ion-ion plasma is ob-
108, 103305-1
© 2010 American Institute of Physics
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jap.aip.org/jap/copyright.jsp
103305-2
J. Appl. Phys. 108, 103305 共2010兲
Kawamura, Lichtenberg, and Lieberman
tained downstream from a narrower source region. These results are compared to a coupled two-region global model. The results are summarized in Sec VI. II. 2D PIC SIMULATIONS
We use a 2D3v rectangular PIC code PDP2 共Ref. 9兲 to simulate an oxygen discharge in a rectangular metal chamber, with insulating walls. In this section, as well as Secs. III and IV, the chamber is 20 cm long with a parallel magnetic field B0 in the x-direction and 3 cm wide in the y-direction. A shorter device of 10 cm length was also explored. Because of the insulating walls of thickness 1/2 cm in the x direction and 1/4 cm in the y direction, the actual discharge dimensions were length L = 19 cm and width 2R = 2.5 cm, for the longer device. For the shorter system, the insulating walls were 1/4 cm in all directions, resulting in L = 9.5 cm and 2R = 2.5 cm. By 2D3v, we mean 2D displacement with three velocity components 共3v兲. In our simulation, each computer particle represents a cluster of about 106 real particles 共electrons or ions兲. The simulations are run with a sufficient number of particles to minimize the discrete particle noise. For each particle, PDP2 tracks the two displacement, and the three velocity components. Collisions between the particles and the background gas are included by using a Monte Carlo collision 共MCC兲 algorithm. A detailed description of the PIC-MCC model used is given in reference.10 In a typical electrostatic PIC simulation, for each timestep: 共i兲 the charge density 共x , y兲 is obtained by a linear weighting of the particles to the spatial grid; 共ii兲 共x , y兲 is used in Poisson’s equation to solve for the electric field E共x , y兲; 共iii兲 the force F = qE is used to advance the particles to new x and v; 共iv兲 the boundaries are checked, and out-ofbounds particles are removed; and 共v兲 a MCC handler checks for collisions and adjusts the particle velocities and numbers accordingly. PIC simulations obtain the fields, as well as the particle fluxes and densities, self-consistently, without making any assumptions about the particle temperatures or velocity distributions.
from 300 to 600 G. Following the experiments which used oxygen as the feedstock gas,4,6 the parameters for oxygen were employed.10,11 However, to scale the low PIC density to a higher physical density, the positive-negative ion recombination rate coefficient Krec was scaled by a factor Krecfac, since the Krecne product occurs in the equations 共see Ref. 11兲. For example, a scaling factor Krecfac = 20 is used to scale to a factor of 20 higher density than that used in the PIC simulation. This still results in a relatively low recombination rate as was investigated in most experiments and previous theory. Higher recombination rates are also explored as described below. Because the mass of oxygen atoms and molecules is not very large, the oxygen ions were significantly magnetized. To avoid this effect in understanding the magnetization of only the electrons, the ions were not magnetized in the fluid calculation. It is straightforward to investigate PIC results, and also the accompanying models, with both magnetized and unmagnetized ions, and this has been done in the results presented below. In the following sections, we also investigate the effect of higher masses, of heating that does not produce uniform transverse temperatures, and configurations that are not axially uniform. The plasma production and heating in an experiment is with a complex inductive field that probably involves the excitation of waves.4,6 In the PIC calculations a simpler form of plasma heating and maintenance is used: an rf current density directed parallel to the B-field. The local plasma conductivity acting on the current produces an electric field that accelerates and thus heats the plasma electrons. In the first set of PIC simulations, at 20 mTorr, Te is equilibrated along the field lines, but not necessarily across the field lines. Since our previous equilibrium and model calculations assumed uniform Te, for the PIC simulations in this section we designed a current profile to approximately achieve this result. Since, from the form of the plasma conductivity, the rf electric field was expected to vary as E ⬀ J / n, where J is the local current density, to achieve E ⬇ constant, transversely, we explored current profiles that might approximate J ⬀ n. The two transverse profiles used are
A. The PIC requirements
For stability and accuracy, 2D PIC simulations require a large number of particles per Debye length squared. Thus PIC codes run faster when simulating small plasma reactors with lower densities 共larger Debye lengths兲. Consequently, the PIC reactor was chosen to have a width of 3 cm and a length of 10–20 cm, about a factor of two to three smaller than typical experimental devices to which it can be compared. The average electron density in the device was below 1014 m−3, well below experimental parameters. These modifications of dimensions and density to accommodate PIC requirements necessitated rescaling of various parameters, as described below. The PIC dimensions, which are approximately half of the dimensions previously used in fluid equations8 to study the experiment,4,6,7 require increasing the pressure in the base case from 10 mTorr in typical experiments to 20 mTorr, to keep pd = const. To keep the degree of electron magnetization constant 共see Ref. 5兲 the magnetic field was also increased by a factor of 2, in the base case
and
冉 冊 冉 冊 4x y sin 2R L
共2兲
冉 冊 冉 冊
共3兲
Jx = J0 cos2
Jx = J0 cos
4x y sin . 2R L
The former approximates the profiles of the special solutions, satisfying Eq. 共1兲, while the latter is a profile found for equilibrium electron densities not dominated by either magnetization or negative ion effects. Since electrons move relatively freely along the B-field, the axial form of J was not considered important and a double sine with zero at the ends was used, to avoid excessive heating at the axial edges. Values of J0 in the range of 1 – 10 A / m2 were tried in the calculations. The higher values were initially considered to produce electron densities closer to the central range of the densities explored in experiments3–7 and the previously reported fluid calculations.8 However, an unexpected result
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jap.aip.org/jap/copyright.jsp
103305-3
J. Appl. Phys. 108, 103305 共2010兲
Kawamura, Lichtenberg, and Lieberman
(b)
2 1.5 1 0.5
1×10
0.005
0.01
1.5 1
0.015
y (m)
0.02
0
0.025
10×Ο −
O−
1×10
13
e− 0.01
0.015
y (m)
0.01
0.005
0.015
y (m)
0.02
0.025
0.02
0.025
O +2 O−
14
5×10
0.005
10×Ο +2
(d)
O2+
14
0
2
0.5
(c)
5×10
e−
2.5
Transverse T (V)
2.5
0
Transverse Density (m-3)
3
(a)
Axial Density (m-3)
Transverse Potential (V)
3
13
e− 0
0.025
0.05
0.075
0.1
x (m)
0.125
0.15
0.175
FIG. 1. PIC results for oxygen reaction rates and unmagnetized ions with Krecfac = 20, B0 = 600 G, L = 19 cm, R = 1.25 cm, p = 20 mTorr, and J0 = 1 A / m2.
was that the combination of the applied frequency of 13.56 MHz and the higher electron densities that were obtained at the higher current densities induced edge electric fields. This in turn produced higher electron temperatures near the transverse edges of the device. The lower current densities, resulting in lower electron densities, minimized this effect, and consequently, were used in most simulations. Because electrons move relatively freely along the field lines, large potentials build up at the device ends to keep the electron and ion flow ambipolar. Conducting transverse walls would result in potentials that inhibit negative ions from reaching them. Consequently, the transverse walls are insulated to obtain the desired type of transverse ion flows. The end walls could be either insulating or conducting. Both were tried in the simulations with some differences found due to the nonuniformity of the transverse electron density. We use insulated ends, which prevented nonzero net current densities flowing to the end walls, for the results of this section, but also used grounded metal end walls for other configurations. We found no significant difference in endloss flux between insulating and conducting end walls. B. PIC simulation results
Figure 1 shows PIC results for a 2D oxygen reactor with all insulating boundaries and unmagnetized ions. The discharge length L = 19 cm 共along B0兲 and transverse width 2R = 2.5 cm. Both the 1/4 cm dielectric in the transverse y direction and the 1/2 cm dielectric in the axial x direction have relative permittivity of four to separate the plasma from conducting boundaries. The pressure is p = 20 mTorr 共ng = 6.4⫻ 1020 m−3兲, the magnetic field is B0 = 600 G, and Krecfac = 20. The current profile has the form 共3兲 with J0 = 1 A / m2. The transverse potential is shown in Fig. 1共a兲,
transverse Te in Fig. 1共b兲, transverse densities in Fig. 1共c兲, and axial densities in Fig. 1共d兲. As expected, the potential profile, governed essentially by the ion energy, is nearly flat. The temperature Te is fairly constant in y, with slight increases near the transverse edges from the induced fields, as explained earlier. Both the electron density profile and the negative ion profile are nearly cosine functions. A selfconsistent determination of these profiles would be quite difficult, so we alternately accept them in our modeling. We note in this case that ␣0 ⬅ n−0 / ne0 ⬇ 4, and that ␣共y兲 ⬅ n−共y兲 / ne共y兲 ⬇ ␣0. Increasing the recombination to Krecfac = 100, with roughly the same electron density in the PIC simulation, the recombination has become significant. However, the only significant difference from Fig. 1 is a reduction in ␣0 to ␣0 ⬇ 3.0, due to the volume recombination. There are some other smaller differences which we will examine in the next section, after developing a model to compare to the simulations. If the ions are magnetized, this results in only a minor increase in O− due to slightly better confinement because the magnetizing term reduces the diffusion coefficient by only about a factor of two 关see 共7兲兴. Leaving the ions unmagnetized, but increasing the magnetic field by a factor of two to B0 = 1200 G reduces ␣0 somewhat due to some suppression of electron losses. The main effect is a very significant reduction in the radial electron loss, as will be seen in Sec. III. III. ANALYTIC MODEL AND COMPARISON WITH SIMULATIONS A. The model
We construct a 1D analytic model to compare to the PIC results. The key result that we take from the simulations is the observation that the transverse profiles of electrons and
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jap.aip.org/jap/copyright.jsp
J. Appl. Phys. 108, 103305 共2010兲
Kawamura, Lichtenberg, and Lieberman
with
L =
共4兲
冉 冊
d 2n e = 共iz − att兲ne − L共1 + ␣0兲ne . dy 2
共6兲
ⴱ is the transverse electron diffusion coefficient obHere, D⬜e tained from the complete ambipolar relations, as derived in ⴱ ⬇ D⬜e, the Appendix. Note that within a small factor D⬜e where for all species the transverse diffusion coefficients without ambipolar fields are11
De,i 2 2 , 1 + ce,i /e,i
De,i =
eTe,i . me,ie,i
共7兲
The solution to Eq. 共6兲 is obtained, using the observation from the PIC results that ne ⬇ 0 at y = ⫾ R 共note a shift in the y-coordinate from Fig. 1兲 ne共y兲 = ne0 cos
y , 2R
共8兲
such that Eq. 共6兲 yields ⴱ D⬜e
冉 冊 2R
2
ne0 = 共iz − att兲ne0 − L共1 + ␣0兲ne0 .
共9兲
Factoring out ne0, Eq. 共9兲 determines Te for given values of ␣0 and the other parameters. The negative ion balance equation is ⴱ − D⬜n
100Krecfac
0.2
d 2n − = attne − Krecn−共ne + n−兲, dy 2
共10兲
ⴱ with D⬜n calculated in the Appendix and ne as in Eq. 共8兲, the integration over the transverse dimension gives
4R ⴱ 2 ␣0ne0 = attne0 − Krecne0 ␣0共1 + ␣0兲R. D⬜n R
共11兲
Together with Eq. 共9兲, this gives Te and ␣0, for a given feedstock gas, in terms of device dimensions L and R, pressure p, the magnetic field B0, and the product Krecne0. For the model computations, we use a base case of L = 19 cm, R = 1.25 cm, and p = 20 mTorr 共ng = 6.4 ⫻ 1020 m−3兲, as in the PIC simulations. We magnetize the electrons with a 600 G axial magnetic field, and compare results with magnetized and unmagnetized ions, using oxygen parameters. We use the same cross-section set in the PIC
0.6 0.4
Irec /I⊥n Iz /I⊥p
0.2
20Krecfac 1
2
α0
3
4
5
0 0
1
2
α0
3
4
5
(c)
共5兲
In Eq. 共4兲, the flux and ion densities are taken to be surface- and volume-averages. Using Eq. 共5兲, approximating ni as ne共1 + ␣0兲, and recognizing that the electron and ion end losses must be equal if the sidewalls are insulated, the electron diffusion equation is
D⬜e,i =
0.6 0.4
(b)
I⊥e /I⊥p
0.8
1/2
.
0.8
(a)
0 0
2Di Te 1+ RL Ti
ⴱ − D⬜e
x 10
Current Ratios
2R ⌫z , Ln i ⬇ RL
−5
1
K recfac ne0 /n g
ions are quite similar, such that ␣0 ⬇ constant. An approximation taken from our previous work,8 is that the ion end loss can be related to the volume ion loss through an ion end loss frequency L
0.25
D (m 2/s)
103305-4
D⊥n ∗ D⊥n
0.2 0.15 0.1
D⊥p ∗ D⊥e
0.05 0 0
1
2
α0
3
4
5
FIG. 2. Theoretical results for oxygen reaction rates with Krecfac = 20 and 100, B0 = 600 G, L = 19 cm, R = 1.25 cm, p = 20 mTorr, and J0 = 1 A / m2.
simulations and the model. We plot results against the more physically intuitive quantity ␣0, rather than Krecfacne0 / ng. The two are uniquely related, as we see in Fig. 2共a兲. Note that the PIC simulation corresponding to this case has a density of 2.5⫻ 1013 m−3, such that the horizontal dashed and dotdashed lines correspond to PIC simulations with Krecfac = 20 and 100. They also correspond in the model to Krecfac = 1 and higher densities. Figure 2共b兲 gives important current ratios. Figure 2共c兲 gives the cross-field diffusion coefficients, showing that the ambipolar corrections to the diffusion coefficients are small 共see Appendix兲 and that the negative ion and electron diffusion coefficients are roughly equal. Because the recombination is relatively small compared to the negative ion flux to the wall, ␣0 does not change significantly for the two values of Krecfacne0 / ng, which differ by a factor of 5 关see Fig. 2共a兲兴. From the model equations, and from the results for the specific examples, we can draw some scaling relations which may be of use to designers. Using the approximation that the fluxes are uniform over their respective axial dimensions, the transverse electron and negative ion wall currents are obtained from Eqs. 共9兲 and 共11兲 as ⴱ L ne0 , I⬜e = 2D⬜e R
共12兲
ⴱ L ␣0ne0 . I⬜n = 2D⬜n R
共13兲
The axial loss current is obtained from Eq. 共4兲 as Iz = LRL共1 + ␣0兲ne0 .
共14兲
Substituting L from Eq. 共5兲 into Eq. 共14兲, we find
冉 冊
Iz = 2Di 1 +
Te Ti
1/2
共1 + ␣0兲ne0 .
共15兲
From the negative ion balance Eq. 共11兲, ne0 and ␣0 are related by
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jap.aip.org/jap/copyright.jsp
103305-5
J. Appl. Phys. 108, 103305 共2010兲
Kawamura, Lichtenberg, and Lieberman
ne0 =
ⴱ 4Ratt/ − D⬜n ␣0/R . Krec␣0共1 + ␣0兲R
共16兲
Equation 共16兲 determines ␣0 for any choice of ne0 关see Fig. 2共a兲兴. As can be seen from Eq. 共16兲, ␣0 must lie in the range ⴱ . The designer will generally want to 0 ⬍ ␣0 ⬍ 4attR2 / 2D⬜n achieve: 共a兲 high I⬜n 共high thrust兲, 共b兲 ne0 / ng not too small 共high propellant utilization兲, 共c兲 low Irec and low Iz 共high power efficiency兲, and 共d兲 low I⬜e 共small transverse electron effects兲, by proper choice of the system parameters R, L, mi, ng, B0, and ne0. From the above relations we identify the following relevant ratios in terms of physical parameters. From Eqs. 共12兲 and 共13兲, assuming electrons fully magnetized and ions essentially unmagnetized, we obtain ⴱ n2m1/2 1 I⬜e D⬜e ⬇ ⴱ ⬀ g2 i , I⬜n D⬜n ␣0 B 0␣ 0
共17兲
which is relevant to condition 共d兲. From Eq. 共11兲 or 共16兲, assuming reasonably high ␣0, we have Irec Krec␣0ne0R ⬀ Krec␣0R2nengm1/2 ⬇ i , ⴱ I⬜n D⬜n /R
共18兲
which is relevant to 共c兲. From Eqs. 共15兲 and 共13兲, assuming ⴱ ⬇ Di, we obtain Te / Ti Ⰷ 1 and D⬜n Iz I⬜n
⬇
冉冊
1 Te Ti
1/2
R , L
共19兲
which is also relevant to 共c兲. Finally, ␣0 can be solved from Eq. 共16兲 to give the limiting cases
␣0 ⬇
4R2 Kattng ⬀ R2Kattn2gm1/2 i , ⴱ 2 D⬜n
共20兲
for Irec / I⬜n small, and
␣0 ⬇
冉
4Kattng Krecne0
冊
1/2
,
共21兲
for Irec / I⬜n large. We note that the development in this section is effectively a global model once the approximation of ␣0 = const has been used. We will return to a global formalism in the more complicated geometries of subsequent sections.
TABLE I. Comparison of PIC with theory for oxygen reaction rates and magnetized ions with Krecfac = 20 and 100, B0 = 600 G, L = 19 cm, R = 1.25 cm, p = 20 mTorr, and J0 = 1 A / m2. 共a兲 Krecfac = 20
ne0 共m−3兲 Iatt I⬜e / I⬜n Irec / I⬜n I⬜n / Iatt ␣0 Te 共V兲 Iiz 共s−1兲 Iz / I⬜p ⴱ D⬜e 共m2 / s兲 ⴱ D⬜n 共m2 / s兲 L 共s−1兲
PIC
Theory
PIC
Theory
2.4⫻ 1013 1.1⫻ 1013 0.24 0.11 0.90 4.8 2.3 1.7⫻ 1013 0.38 0.25 0.29 1.5⫻ 103
2.4⫻ 1013 0.93⫻ 1013 0.21 0.11 0.90 3.9 2.1 1.5⫻ 1013 0.36 0.15 0.18 1.0⫻ 103
2.4⫻ 1013 1.0⫻ 1013 0.29 0.41 0.72 3.3 2.3 1.8⫻ 1013 0.43 0.19 0.20 2.0⫻ 103
2.4⫻ 1013 0.93⫻ 1013 0.28 0.46 0.68 3.0 2.1 1.4⫻ 1013 0.36 0.20 0.18 0.99⫻ 103
ues of L from the model to the simulations 共with both insulating and conducting end walls兲 indicates that the formula given in Eq. 共5兲 is reasonable. We also find that the ratio I⬜e / I⬜n can be reasonably well determined from the average diffusion coefficients and the assumption of a constant ␣0 I⬜e/I⬜n =
ⴱ D⬜e
ⴱ ␣0D⬜n
共22兲
.
Using the PIC results for the table in Eq. 共22兲 for case 1a and 1b, we obtain I⬜e / I⬜n = 0.18 and 0.29 respectively, in good agreement with the ratios in the table. Comparing case 1b to 1a, ␣0 has been significantly decreased both in the model and the simulation, as expected from increasing the recombination factor Krecfac. Irec / I⬜n has increased in the PIC results by a factor of 3.7, in good agreement with the factor of 3.4 increase from the scaling 共⬀Krec␣0兲 in Eq. 共18兲. In Table II, we present PIC results from varying important quantities which can be compared with the scalings from the model. In case 2a for unmagnetized ions, the theory preTABLE II. PIC results for oxygen reaction rates and both magnetized and unmagnetized ions with Krecfac = 20, R = 1.25 cm, p = 20 mTorr, and J0 = 1 A / m2. Unless otherwise noted, B0 = 600 G, L = 19 cm, and Kattfac = 1.
B. Comparison of model with PIC simulations
In Table I, we compare the results of two cases we have examined in Fig. 2 with the values obtained from PIC simulations. In the theory, we use the same density ne0 for the comparisons, but have multiplied the ion diffusion coefficients by a factor of 0.5 to make the ␣0’s similar in all comparisons. Such an adjustment could be due to profile differences between the model and simulation, or to the need for more accurate kinetic calculations of the diffusion. The results mainly speak for themselves, but we point out a few salient features. We see good agreement between PIC and theory for the important ratio of electron to ion radial flux I⬜e / I⬜n. One of the motivations of the simulations is to verify that the formula for L, which involves a number of approximations, is reasonably accurate. Comparing the val-
共b兲 Krecfac = 100
Unmag. ions
ne0 共m−3兲 Iatt I⬜e / I⬜n Irec / I⬜n I⬜n / Iatt ␣0 Te 共V兲 Iiz 共s−1兲 Iz / I⬜p ⴱ D⬜e 共m2 / s兲 ⴱ D⬜n 共m2 / s兲 L 共s−1兲
Mag. ions
共a兲 B0 = 600 G
共b兲 B0 = 1200 G
共c兲 L = 9.5 cm
共d兲 Kattfac = 3
2.5⫻ 1013 1.1⫻ 1013 0.23 0.077 0.92 3.9 2.3 1.7⫻ 1013 0.37 0.31 0.36 1.7⫻ 103
2.7⫻ 1013 1.1⫻ 1013 0.052 0.063 0.95 3.2 2.0 1.6⫻ 1013 0.40 0.05 0.28 1.8⫻ 103
2.1⫻ 1013 4.9⫻ 1012 0.27 0.14 0.89 7.5 2.8 1.2⫻ 1013 1.0 0.30 0.21 4.1⫻ 103
2.1⫻ 1013 3.2⫻ 1013 0.10 0.23 0.79 13 2.8 4.1⫻ 1013 0.26 0.24 0.29 1.0⫻ 103
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jap.aip.org/jap/copyright.jsp
103305-6
J. Appl. Phys. 108, 103305 共2010兲
Kawamura, Lichtenberg, and Lieberman 4
Transverse T (V)
1×10
e−
(a)
15
(b)
3
8×10
2
n (m -3)
6×10
1
0 0
0.005
0.01
10×I −
4×10
10×I +
2×10
0.015
y (m)
α
14
1000
I+ I−
14
α
14
100
14
10×e − 0.02
0.025
0.03
0 0
0.005
0.01
0.015
y (m)
0.02
0.025
10 0.03
FIG. 3. PIC results of particle 共a兲 temperature 共eV兲 and 共b兲 density profiles for iodinelike discharge and axially confined heating currents with Krecfac = 20, Kattfac = 5 共no threshold兲, B0 = 600 G, L = 19 cm, R = 1.25 cm, p = 20 mTorr, and J0 = 1 A / m2. The ␣ values are shown with a logarithmic scale.
dicts an increase in D⬜n by a factor of two from case 1a, the magnetized-ion case in Table I. The increase in the PIC diffusion coefficient is somewhat less, indicating other factors also enter, but rather weakly. In case 2b we increase the magnetic field by a factor of two, from the base case, 共but keeping the ions unmagnetized to correspond to a heavier mass in which 2ci / 2ci Ⰶ 1兲. Comparing case 2b with case 2a, we see that the ratio I⬜e / I⬜n is decreased approximately by the factor of 4 predicted from the scaling in Eq. 共22兲. The electron profile in this case lies between the model cos共y / 2R兲 profile and the cos2共y / 2R兲 profile which produced the special solution. In case 2c, with magnetized ions, we reduce the length by a factor of two to L = 10 cm 共9.5 cm when the dielectrics are subtracted兲. From 共4兲, the theory predicts that L is increased by approximately a factor of 2 from the value in case 1a. In the PIC simulation, the increase in L is by a larger factor of 2.7, indicating the importance of the end effects for short systems. Note that for case 2c, the axial positive ion loss is approximately equal to the radial positive ion loss 共Iz / I⬜p ⬇ 1兲, such that power would be used rather inefficiently. In case 2d, we investigate the effect of a larger attachment coefficient by increasing Katt by a factor of three from its value in case 1a. Since the recombination is small compared to the attachment 共less than 20% in this case兲 we expect, from 共11兲, that ␣0 would increase essentially by the same factor of three, which we confirm from the simulation by comparing the value in case 2d with that in case 1a. C. Discussion of uniform Te results
From the previous PIC simulations and accompanying analysis, we have found that it is quite difficult to achieve a result approximating the special solution 共1兲, for which the electron flux at the radial walls approaches zero, leaving equal numbers of positive and negative ions to be extracted. This solution can be approached at high magnetic field strengths, provided a heavy ion mass is employed such that the ions are not strongly magnetized. However, high magnetic fields 共B0 = 1200 G in the small device we considered兲 may not be practical for a particular application. In the PIC calculations and accompanying model, we have considered relatively uniform transverse electron temperatures 共exactly uniform in the model兲. Such a distribution is not necessary if the electron transport is significantly inhibited across the field. Achieving cooler edge electrons requires both a struc-
ture of the heating fields that heat primarily away from transverse walls, and also an appropriate geometry, magnetic field strength, and other parameters favoring axial electron loss over radial electron loss. Another factor that enters into the effectiveness of reducing the electron flux to the wall is the energy dependence of the attachment cross-section. If attachment is efficient at low energies, then a reduction in temperature near the transverse walls will locally favor attachment over ionization, and, therefore, deplete the electron population. Oxygen has an energy threshold for attachment of 4.2 eV, and thus at the generally fairly low electron energies existing in magnetically confined devices does not efficiently accomplish attachment electron absorption. IV. NONUNIFORM Te AND INCREASED Katt
The above discussion suggests studying a feedstock gas resembling iodine rather than oxygen, and heating the discharge nonuniformly, such that the edge electrons are cooler than those near the axis. Iodine has not been well characterized, but some information on attachment is available. We use mostly the well-known oxygen cross-sections except that we increase the attachment cross-section by a factor of 5 and remove the attachment threshold, as suggested by the literature.12 As in Sec. III, we account for the lower PIC density by using an enhanced recombination factor Krecfac = 20. The heavier mass of iodine is used, such that at the base case of B0 = 600 G, the ions are essentially unmagnetized but also diffuse more slowly. For the previous geometry to primarily heat the electrons away from the transverse walls, we use for 兩y兩 ⬍ R / 2 a transverse current profile J共y兲 = J0 cos2共y / R兲, with J0 = 1 A / m2 and R = 1.25 cm, with J共y兲 = 0 for 兩y兩 ⬎ R / 2. In Fig. 3共a兲, we give Te共y兲, and in Fig. 3共b兲 the densities for nI+ 共solid兲, nI− 共dashed兲 and 10ne 共dotted兲, to graph it on the same scale as the ion densities. The ␣ = nI− / ne values 共dot-dashed兲 are also shown in Fig. 3共b兲 but with a logarithmic scale. The central ␣0 ⬇ 27 is approximately a factor of 5 higher than for the oxygen case, reflecting the increase in Katt by 5, and the absence of a cutoff. The lowering of Te, mainly outside of the central half of the device, where the current is applied, has greatly decreased the ionization in these outer regions, so decreasing the electron density. We find the edge ␣共R兲 ⬇ 2000, which is quite high, resulting in I⬜e / I⬜n ⬇ 4 ⫻ 10−4, i.e., an ion-ion plasma. As one also notes, the electron density profile now approximates the electron profile
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jap.aip.org/jap/copyright.jsp
103305-7
J. Appl. Phys. 108, 103305 共2010兲
Kawamura, Lichtenberg, and Lieberman
with negligible wall electrons, studied previously,8 but is due to decreasing Te near the walls, rather than to a special uniform Te solution. Introducing a threshold energy of 4.2 eV for the attachment leads to a modest decrease in the central electronegativity to ␣0 ⬇ 15. The major change is that the edge electron density does not fall to such a low value. The resulting edge ␣共R兲 ⬇ 200 and ratio I⬜e / I⬜n ⬇ 1.5⫻ 10−3 are still in a useful range. The PIC iodinelike results show a high ionization frequency iz inside an active core of thickness 2R1 ⬇ 1.25 cm. In the edge region iz Ⰶ att, such that a tworegion diffusion model with no ionization in the edge region can be used to determine the equilibrium parameters. In the cases considered, the recombination is small compared to the attachment, so we also neglect recombination to make the analytic model more tractable. The electron and negative ion balance relations are then described by −
2 ⴱ d ne D⬜e 2
dy
ⴱ − D⬜n
= 共iz − att − L兲ne − Lnn ,
d 2n n = attne , dy 2
共23兲
共24兲
where we approximate iz to be constant in the core and zero in the edge region. The diffusion coefficients are given in Appendix. Differentiating Eq. 共23兲 twice and substituting Eq. 共24兲 to eliminate nn, we obtain 4 d2ne Latt ⴱ d ne − D⬜e = 共 − − 兲 + ⴱ ne . iz att L dy 4 dy 2 D⬜n
共25兲
In the core we can substitute a symmetric eigenfunction cos ky into Eq. 共25兲 to obtain the characteristic equation ⴱ 4 k − 共iz − att − L兲k2 + D⬜e
Latt = 0, ⴱ D⬜n
共26兲
which has two real positive roots k1 and k2 for 共iz − att ⴱ ⴱ − L兲2 ⬎ 4D⬜e Latt / D⬜n . The core solution can therefore be written as ne = A1 cos k1y + A2 cos k2y.
共27兲
In the edge, with iz = 0, we substitute the eigenfunctions either sinh 共R0 − y兲 or cosh 共R0 − y兲 into Eq. 共25兲 to obtain ⴱ 4 − 共att + L兲2 + D⬜e
Latt = 0, ⴱ D⬜n
共28兲
which yields two positive real roots 1 and 2 for 共att ⴱ ⴱ + L兲2 ⬎ 4D⬜e Latt / D⬜n . Guided by the PIC results, we use an edge solution for ne as in Eq. 共1兲 having both ne = 0 and dne / dy = 0 at the outer edge to obtain ne = A3关sinh 1共R0 − y兲 − 共1/2兲sinh 2共R0 − y兲兴 + A4关cosh 1共R0 − y兲 − cosh 2共R0 − y兲兴.
共29兲
Matching the core and edge ne’s and their first three derivatives at y = R1, we obtain four linear equations M · A = 0, where A is a column vector of the four unknown coefficients A1, A2, A3, and A4, and M is the corresponding 4 ⫻ 4 coeffi-
1×10 8×10 6×10
15
I+ I−
14
α 1000
14
n (m -3)
α
14
100
4×10 2×10
14
0 0
10×e − 0.005
0.01
0.015
y (m)
0.02
0.025
10 0.03
FIG. 4. Profiles of the ion densities nI+ 共solid兲, nI− 共dashed兲, and 10⫻ ne 共dotted兲 for the two-region model. The ␣ values 共dot-dashed兲 are shown with a logarithmic scale.
cient matrix. Setting det M = 0 yields the eigenvalue Te, the core electron temperature 共as well as all temperatureⴱ , etc.兲. The ratios of the A’s dependent quantities iz , att , D⬜e are then determined by the ratios of the cofactors of any row of M. Hence, ne is determined up to a constant factor. With ne共y兲 known, ni can be determined using the positive ion balance relation ⴱ − D⬜i
d 2n i = izne − Lni , dy 2
共30兲
again with iz a 共now-known兲 constant in the core and zero in the edge region. The solution to 共30兲 in the core is ni = B1 cos k1y + B2 cos k2y + B3 cosh Ly, ⴱ 2 k1,2 + L兲 B1,2 = A1,2iz / 共D⬜i
共31兲
ⴱ 1/2 L = 共L / D⬜i 兲 .
and The where first two terms are a particular solution of Eq. 共30兲 in the core and the third term is the symmetric homogenous solution. In the edge region, again guided by the PIC solution, we choose ⴱ dni / dy = niuBi at the outer the boundary condition that −D⬜i edge y = R0, where uBi = 共eTi / mi兲1/2 is approximately the Bohm velocity at the 共highly electronegative兲 edge. Then from Eq. 共30兲 with iz = 0, we obtain ⴱ ni = B4关sinh L共R0 − y兲 + 共D⬜i L/uBi兲cosh L共R0 − y兲兴.
共32兲 Matching the core and edge densities Eqs. 共31兲 and 共32兲, and their first derivatives, at y = R1, we obtain two linear equations to determine the two unknown coefficients B3 and B4. Hence, ne and ni 共as well as nn by quasineutrality兲 are determined up to a constant factor. Figure 4 shows the calculated density profiles versus transverse coordinate y for 10⫻ ne 共dotted兲, nI+ 共solid兲, and nI− 共dashed兲 from this two-region model for the iodinelike discharge shown in Fig. 3. The simulation values of R1 ⬇ 0.5 cm and the central electron density ne0 are used in the model. The results are plotted with the same scales as the PIC results in Fig. 3共b兲 to facilitate a comparison. There is reasonable agreement for both the shapes of the transverse density and ␣ profiles, and the central values of ␣. The shapes of the ion and electron radial profiles, both for the PIC simulation and for the two-region theoretical
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jap.aip.org/jap/copyright.jsp
103305-8
J. Appl. Phys. 108, 103305 共2010兲
Kawamura, Lichtenberg, and Lieberman
y (cm)
8 1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
5 mm
1
1
6 4
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
6 cm
1
1 1
1
1
1
1
1
1
1
1
1
5
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
4 cm
2 0
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1 1
1 1
1
15
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
6 cm 1
1
1
1
1
1
1 1
10
1
1
1
1
1
1
1
20
1
1
1
1
25
1
x (cm)
30
FIG. 5. Configuration with source region in center connected to two wider downstream chambers.
model, indicate that a global model can be constructed which will give reasonably accurately the important parameters, and furthermore, can be used for scaling. As will be done for a more complicated geometry in the following section, we use a cosine profile for the electrons over the central region, and a cosine profile for the negative ions over the entire diameter. The basic equations are the electron balance relation
iz = att + L共1 + ␣0兲
R0 R1
共33兲
and the negative ion balance
att =
2 R0 ⴱ Krecne0␣0共1 + ␣0兲 + D⬜n ␣0 , 4 R1 4R0R1
共34兲
where ␣0 is the on-axis electronegativity. Comparisons with the two-region model give good agreement, indicating that global models can be constructed in geometries with a hot electron core and a cold electron periphery.
V. MORE COMPLICATED GEOMETRIES
In experiments,6,7 an ion-ion plasma was produced near a boundary by having a step in the diameter with the main ionization source for electrons 共the heating region兲 confined to the narrower region. The cooler electron temperature near the outer part of the wider downstream region removed electrons by attachment, as described in Sec. IV. In the 2D rectangular approximation the PIC boundary conditions can readily be modified to a step geometry with the source heating confined to the narrower region. The configuration is shown in Fig. 5. For the cases investigated, the dielectric transverse 共radial兲 walls in the central source region were separated by 4 cm, and the length of the region was 10 cm. Symmetric downstream regions connected to the source region were each also 10 cm long, with grounded conducting ends. The transverse 共radial兲 walls were dielectric, backed by grounded conductors. Downstream conductor separations of 7 and 10 cm were investigated in simulations, with 1/2 cm dielectrics inside the conductors. The heating by an axial current in the central source region was a cosine in both directions, having the form Jx = J0 cos共y / 2R0兲cos共x / L0兲, where R0 is the half-width of the central region and L0 the length of the heating zone. As previously, J0 = 1 A / m2, x is the axial direction of the applied magnetic field B0 = 600 G, and y is the cross-field direction. As in Sec. IV, an iodinelike
discharge is simulated by using the base oxygen crosssection set, but increasing the attachment cross-section by a factor of five and removing the threshold. In Fig. 6, the results are given for what we shall call a base case, with Krecfac = 20 and p = 20 mTorr. The metal walls in the downstream regions are separated by 7 cm. Figures 6共a兲 and 6共b兲 give the radial densities at the center of the source region and the center of one downstream end region, respectively. We see that the electrons are well confined in the end region, while the positive and negative ion densities fall-off to the outer walls. In Figs. 6共c兲 and 6共d兲, the axial potential and axial densities, respectively, are shown. The potential fall in the end regions of approximately ⌬V ⬇ 1.5 V leads to rather large ion density decreases due to their low temperatures, but very little change in the on-axis density of the considerably higher temperature electrons. In Fig. 6共e兲, the axial Te is given, and in Fig. 6共f兲 the corresponding axial reaction rates Kiz and Katt. We note that the axial potential fall-off 关in Fig. 6共c兲兴 corresponds mainly to a decreasing electron temperature away from the heating zone 关in Fig. 6共e兲兴. The strong drop in Kiz in the end regions corresponds to the drop in Te, while the axial Katt is relatively unchanged. For the development of a theoretical treatment in Sec. VI, we note that the potential variation, which is a significant fraction of Te, leads to electric field driven positive ions flowing downstream and negative ions flowing upstream. The simulation results indicate that the particle balances in the central and end sections are coupled by mobilitydriven axial 共x兲 fluxes of positive and negative ions. At the junction between the central and right-hand end sections, ⌫i = i0niE and ⌫n = −n0nnE, where the electric field is estimated as E = Te0 / lE, with lE the scale length of a potential with characteristic value Te0. By flux neutrality, the corresponding electron flux is ⌫i − ⌫n. Accounting for these fluxes and with reasonable choices of the transverse density profiles, a simple global 共volume-averaged兲 model of the particle balances in the central and end sections can be used to describe the equilibrium of the coupled system. We let R0 and R1, and L0 and L1 be the transverse half-lengths and axial lengths of the central and end sections, and assume that all central cell transverse density profiles vary as cos共y / 2R0兲. Then for a unit length in the ignorable 共z兲 direction, the electron balance in the central section can be expressed as
iz0 = att0 + 2关i0共1 + ␣0兲 + n0␣0兴
2 Te0 ⴱ + D⬜e0 , l EL 0 4R20 共35兲
where ␣0 = nn / ne is the 共constant兲 electronegativity and the subscripts “0” denote quantities evaluated at the central cell temperatures. The last term gives the transverse electron loss to the insulating walls and the factor of 2 multiplying the bracket indicates a symmetric system with positive ions leaving and negative ions entering from both ends. The central negative ion balance relation can be similarly written as
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jap.aip.org/jap/copyright.jsp
J. Appl. Phys. 108, 103305 共2010兲
Kawamura, Lichtenberg, and Lieberman
1×10
15
I+ I−
(a)
End Cell Transverse n (m-3)
Center Cell Transverse n (m-3)
103305-9
6×10
4×10
5×10
14
2×10
20×e − 0 0
0.02
0.04
I+ I−
(b)
14
20×e −
14
0 0
0.06
y (m)
14
0.01
0.02
0.03
15
Axial Density n (m -3)
Axial Potential (V)
(c) 10
5
0 0
0.05
0.1
0.15
x (m)
0.2
1×10
5×10
15
5
8×10
Axial K (m 3 /s)
Axial T (V)
6×10
e−
4×10
2 1 0 0
2×10
5×I − 5×I + 0.05
0.1
0.15
x (m)
0.2
0.05
0.1
0.15
x (m)
0.2
0.25
-16
iz
-16
-16
att -16
0 0
0.25
0.07
20×e −
(f)
3
0.06
14
(e) 4
0.05
I+ I−
(d)
0 0
0.25
0.04
y (m)
0.05
0.1
0.15
x (m)
0.2
0.25
FIG. 6. PIC results for iodinelike discharge in configuration of Fig. 5 with Krecfac = 20, B0 = 600 G, p = 20 mTorr, and J0 = 1 A / m2.
att0 + 2n0␣0
Te0 2 ⴱ = Krec0ne0␣0共1 + ␣0兲 + D⬜n0 ␣0 2 , l EL 0 4 4R0 共36兲
where the last two terms give the volume recombination loss and the transverse negative ion loss. In a downstream end section, guided by the PIC simulation results, we choose the same on-axis density ne0 and transverse electron profile cos共y / 2R0兲 as in the central cell. The PIC results also indicate a negative ion profile that extends over the entire transverse dimension of the end cell. We model this using a simple cos共y / 2R1兲 transverse distribution and discuss this further in the conclusions. Since the electrons are confined to the smaller transverse dimension, there is a negligible transverse electron loss. However, there is a considerable axial loss, equal to the positive ion axial loss, which we model using the axial loss frequency L given from Eq. 共5兲
冉
L = Di1 1 +
Te1 Ti1
冊
1/2
2 , R 1L 1
共37兲
where we take the losses to be over the full radius on both sides of each end of the downstream sections, and the sub-
script 1 refers to downstream quantities. Equation 共37兲 somewhat overestimates L for this geometry. The electron balance relation in one end cell is
iz1 + 关i0共1 + ␣0兲 + n0␣0兴
Te0 R1 = att1 + L共1 + ␣1兲 , l EL 1 R0 共38兲
where ␣1 = nn共y = 0兲 / ne共y = 0兲 is the on-axis electronegativity. The negative ion balance in an end cell can be written as
att1 = n0␣0
Te0 R1 + Krec1ne0␣1共1 + ␣1兲 l EL 1 4 R0
ⴱ + D⬜n1 ␣1
2 . 4R0R1
共39兲
For a given electron density ne0 and specifying the heavy particle temperature as a fixed ratio of the electron temperature in each of the sections, then Eqs. 共35兲–共39兲 can be solved numerically to determine Te0, Te1, ␣0, and ␣1. Once these quantities are known, then, all of the particle generation and loss terms in the central and end sections can be evaluated.
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jap.aip.org/jap/copyright.jsp
103305-10
J. Appl. Phys. 108, 103305 共2010兲
Kawamura, Lichtenberg, and Lieberman
TABLE III. Comparison of PIC with theory for iodinelike discharge in configuration of Fig. 5 with p = 20 or 10 mTorr, Krecfac = 20, B0 = 600 G, and J 0 = 1 A / m 2. 20 mTorr Central cell ne0 共m−3兲 Iatt I⬜n / Iatt Irec / I⬜n Inin / I⬜n ␣0 Te 共V兲
PIC
Theory
10 mTorr PIC
Theory
End cells
1.0⫻ 1013 1.0⫻ 1013 5.5⫻ 1014 5.0⫻ 1014 0.68 0.62 1.4 1.3 0.73 0.43 105 98 2.8 3.6 20 mTorr PIC Theory
1.0⫻ 1013 1.0⫻ 1013 3.1⫻ 1014 2.5⫻ 1014 1.1 1.2 0.44 0.34 0.62 0.34 54 45 2.9 3.6 10 mTorr PIC Theory
Iatt1 I⬜n1 / Iatt1 Irec1 / I⬜n1 ␣01 Te1 共V兲 Iz / I⬜p1 L 共s−1兲
9.8⫻ 1014 0.13 3.5 58 1.8 1.4 270
6.7⫻ 1014 0.32 1.2 37 2.1 1.4 660
12⫻ 1014 0.27 1.9 62 3.1 0.79 750
5.0⫻ 1014 0.52 0.52 31 3.3 0.38 1540
The negative ion balance equations are insensitive to the electron temperature so that approximate negative ion results can be obtained, assuming nominal temperatures, without solving for Te0 and Te1 from the electron balance equations. However, the scale length lE is sensitive to the source and downstream temperatures, so approximate average values, taken from a few PIC cases, are used. In Table III, we compare PIC with theory for negative ion parameters, in two typical cases. The axial source electron density ne0, taken from the PIC results, is specified in the theory, for the comparison. A nominal value of lE = 0.18 m consistent with simulations at p = 20 mTorr is used in the theory at this pressure. For p = 10 mTorr, the drop in downstream Te is smaller, so that the effective lE is longer. For this pressure, we use lE = 0.36 m. The results in Table III give generally reasonable agreement between PIC simulations and theory. As with the simpler geometries the theoretical L factors work quite well. The approximate reduction in the ratio Irec / I⬜n from p = 20 mTorr to p = 10 mTorr, is consistent with the general scaling relations from Sec. III. The weaker expected drops in ␣0 are also found. In downstream regions the radial electron current is negligible, and therefore has not been listed. We have also investigated 10 cm width end regions, but the main effect was an increase in the overall recombination flux, which is not desirable, so we omit details of the results. VI. CONCLUSIONS AND FURTHER DISCUSSION
We have used a 2D PIC code to investigate the equilibrium of a device with magnetized electrons but only weakly magnetized ions, in an electronegative plasma. The purpose of the simulations is to study the possibility of extracting positive and negative ions in nearly equal numbers. In Secs. II and III, to conform to experiments and previous analysis, an oxygen feedstock gas was used, together with a heating
current that produced a relatively constant Te across the magnetic field lines. A device of 19 cm length and 2.5 cm transverse dimension made the ratio of end loss to transverse loss of the positive ions relatively small, allowing a 1D theoretical approximation to be made. The walls were insulating as is required to extract negative ions. A pressure p = 20 mTorr and an axial magnetic field B0 = 600 G were used in most of the simulations. At these parameters, the ions were significantly magnetized; therefore, to investigate the effect of ion magnetization, simulations were also performed with the ions unmagnetized. Because the PIC simulations require low electron densities, and the relative recombination rate is determined by the product Krecne, an enhanced Krec was used to rescale the recombination to normal values, designated as a Krecfac. For oxygen, at p = 20 mTorr, setting this factor to 20 gave relatively low recombination, while a factor of 100 gave significant recombination. Both were used in the simulations. A current directed along the magnetic field with a transverse cosine profile resulted in a relatively constant transverse Te. The results of the simulations indicated that the electronegativity ␣ = n− / ne was reasonably constant transversely. For fairly small ␣ this resulted in a significant ratio of electron to ion current I⬜e / I⬜n at the transverse wall, which we desired to be small. Higher magnetic fields were found to proportionally reduce this ratio. Because ␣ ⬇ ␣0 a simple model having this property could be developed and results of the model were in good agreement with the simulations. The model also allowed various scaling with parameters to be explored analytically. This, in turn, suggested the PIC parameters to be explored. The results from the simulations and the accompanying model indicated that a significant reduction in the ratio of electron to ion fluxes would be obtained with higher mass ions that suppressed ion magnetization, with higher ␣’s that would result from higher Katt, and with lower transverse edge electron temperature. The first two improvements would be achieved with a feedstock gas such as iodine. The third and probably most important effect can be obtained by limiting the plasma excitation and heating to the central portion of the transverse dimension. In Sec. IV, based on these insights, PIC simulations were performed with more iodinelike parameters of background neutrals, and with the excitation concentrated in the central half of the transverse dimension. The results of two representative simulations indicated excellent suppression of transverse electron flux at a pressure of p = 20 mTorr and a modest magnetic field strength of B0 = 600 G. A two-region model was developed that well- approximated the PIC results. Although parameter scaling can easily be explored, numerically, with the model, there are more parameters to study, such that more extensive PIC simulations would be required to validate the model over a wide range of parameters. Consequently, a global model was developed from which scaling could be obtained with results that were consistent with the two-region model. In Sec. V, a final set of simulations were performed to conform more closely to an experimental geometry which had a source region and wider downstream regions. If the electrons are significantly better confined, in the transverse direction, than the ions, this leads to cooler electrons at trans-
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jap.aip.org/jap/copyright.jsp
103305-11
J. Appl. Phys. 108, 103305 共2010兲
Kawamura, Lichtenberg, and Lieberman
verse positions beyond the transverse boundaries of the source region, and therefore naturally favors the production of ion-ion edge plasmas. This was, in fact, achieved, and coupled source-downstream global models were constructed, which gave good agreement with the PIC simulations. One of the results from all of the modeling was comparisons of the classical diffusion, both radially and axially, to the PIC results. We found that the PIC radial diffusion was approximately classical for the ion species. The axial diffusion coefficients controlled the L’s, which were found to agree reasonably well between theoretical calculations and PIC simulations.
dnn dne + 共D⬜i − D⬜e兲 dy dy . 共⬜i + ⬜n兲nn + 共⬜i + ⬜e兲ne
共D⬜i − D⬜n兲 E=
For a constant transverse electronegativity, nn共y兲 = ␣0ne共y兲,
共A6兲
we substitute Eq. 共A6兲 into Eq. 共A5兲 to eliminate the dependence of E on dne / dy to obtain n nE =
␣0共D⬜i − D⬜n兲 + D⬜i − D⬜e dnn . ␣0共⬜i + ⬜n兲 + ⬜i + ⬜e dy
ⴱ ⌫⬜n = − D⬜n
dnn , dy
ⴱ = D⬜n + ⬜n D⬜n
共A1兲
where the ⌫’s are given by the usual drift-diffusion relations ⌫⬜i = − D⬜i
dni + ⬜iniE, dy
共A2兲
⌫⬜n = − D⬜n
dnn − ⬜nnnE, dy
共A3兲
⌫⬜e = − D⬜e
dne − ⬜eneE. dy
共A4兲
Substituting Eqs. 共A2兲–共A4兲 into Eq. 共A1兲, with ni = nn + ne from quasineutrality, and solving for the electric field, we obtain
␣0共D⬜i − D⬜n兲 + D⬜i − D⬜e . ␣0共⬜i + ⬜n兲 + ⬜i + ⬜e
共A9兲
A similar substitution to eliminate the dependence of E on dnn / dy yields
␣0共D⬜i − D⬜n兲 + D⬜i − D⬜e . ␣0共⬜i + ⬜n兲 + ⬜i + ⬜e
共A10兲
A similar derivation for the positive ions yields
The effective cross-field diffusion coefficients are found from the ambipolar flux balance condition ⌫⬜i = ⌫⬜n + ⌫⬜e ,
共A8兲
then comparing Eq. 共A8兲 to Eq. 共A3兲, we obtain
ⴱ = D⬜e + ⬜e D⬜e
APPENDIX: EFFECTIVE CROSS-FIELD DIFFUSION COEFFICIENTS
共A7兲
Substituting Eq. 共A7兲 into Eq. 共A3兲 and introducing the efⴱ by fective diffusion coefficient D⬜n
ACKNOWLEDGMENTS
E.K., A.J.L., and M.A.L. acknowledge the support of California Industries and University of California Discovery Grant No. ele07-10283 under the IMPACT program and the Department of Energy Office of Fusion Energy Science Contract No. DE-SC0001939. The work is closely related to our previous collaboration with Gary Leray and Pascal Chabert at Ecole Polytechnique, France.
共A5兲
ⴱ = D⬜i − ⬜i D⬜i
␣0共D⬜i − D⬜n兲 + D⬜i − D⬜e . ␣0共⬜i + ⬜n兲 + ⬜i + ⬜e
1
共A11兲
S. K. Kanakasabapathy, L. J. Overzet, V. Midha, and D. J. Economou, Appl. Phys. Lett. 78, 22 共2001兲. A. Meige, G. Leray, J.-L. Raimbault, and P. Chabert, Appl. Phys. Lett. 92, 061501 共2008兲. 3 R. Kawai and T. Mieno, Jpn. J. Appl. Phys., Part 2 36, L1123 共1997兲. 4 P. Chabert, T. E. Sheridan, R. W. Boswell, and J. Perrin, Plasma Sources Sci. Technol. 8, 561 共1999兲. 5 S. G. Walton, D. Leonhardt, R. F. Fernsler, and R. A. Meger, Appl. Phys. Lett. 81, 987 共2002兲. 6 C. S. Corr, N. Plihon, and P. Chabert, J. Appl. Phys. 99, 103302 共2006兲. 7 A. Aanesland, A. Meige, and P. Chabert, J. Phys.: Conf. Ser. 162, 012009 共2009兲. 8 G. Leray, P. Chabert, A. J. Lichtenberg, and M. A. Lieberman, J. Phys. D: Appl. Phys. 42, 194020 共2009兲. 9 V. Vahedi and G. DiPeso, J. Comput. Phys. 131, 149 共1997兲. 10 V. Vahedi and M. Surendra, Comput. Phys. Commun. 87, 179 共1995兲. 11 M. A. Lieberman and A. J. Lichtenberg, Principles of Plasma Discharges and Materials Processing, 2nd ed. 共Wiley-Interscience, New York, 2005兲. 12 G. E. Caledonia, Chem. Rev. 75, 333 共1975兲. 2
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jap.aip.org/jap/copyright.jsp
