- Email: [email protected]
Multidimensional fluid–particle modelling technique in low-temperature argon plasma at low pressure
ARTICLE IN PRESS
Vacuum 82 (2008) 220–223 www.elsevier.com/locate/vacuum
Multidimensional fluid–particle modelling technique in low-temperature argon plasma at low pressure Petr Bartos˘ , Rudolf Hrach, Petr Jelı´ nek Department of Surface and Plasma Science, Faculty of Mathematics and Physics, Charles University, V Holesˇovicˇka´ch 2, 180 00, Prague 2, Czech Republic
Abstract In the contribution there is presented a new computer modelling technique, which can be used in multidimensional computer simulations of plasma at low pressure. The technique is based on the fluid–particle approach and can be used in the study of the range of physical phenomena, e.g. plasma behaviour in the vicinity of substrates, probe diagnostics, plasma technology, etc. The computer model is useable for the solution of three-dimensional problems and it minimizes the restrictions arising in complicated model geometry. The algorithm has been tested during the study of plasma sheath formation in the vicinity of probes with various geometries immersed into low-temperature argon plasma. Some results of these simulations are presented in the last section. r 2007 Elsevier Ltd. All rights reserved. PACS: 52.65.y; 52.40.Kh Keywords: Fluid–particle simulations; Plasma sheath; Argon plasma
1. Introduction
2. Model description
In recent years plasma is used widely in science, industry and technology. Scientific teams endeavour to optimize the plasma parameters of devices and so achieve the best performance of production lines. Unfortunately, this is often a very expensive and time-consuming activity. In the last decade the performance of computer systems has increased rapidly and therefore a new and less timeconsuming approach, computer modelling, can be used. The plasma parameters in technological devices can be predicted by this way with high precision, faster and more cheaply. On our department an effective computer model of low-temperature plasma has been developed. This model combines the advantages of both the fluid modelling technique (fast computation and therefore the possibility of the computer modelling of multidimensional problems) and particle modelling approach (precision in the description of collision processes in the plasma).
2.1. Model philosophy
Corresponding author.
E-mail address: [email protected] (P. Bartos˘ ). 0042-207X/$ - see front matter r 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.vacuum.2007.07.006
Models based on the fluid modelling technique are often used in recent studies. Unfortunately, the precision of these models depends on many parameters (diffusion coefficients and coefficients of mobility, collision frequencies, coefficient rates, etc.), which are very often complicated functions of energy. These dependences can be determined by experimental measurements and expressed as semiempirical formulas, e.g. Ref. [1], but this approach is too complicated. Alternatively, fluid–particle modelling techniques can be used, e.g. Refs. [2,3]. We developed a model of such type, the model arrangement is presented in Fig. 1. The algorithm consists of four separate parts. In the first step, there is solved a simplified fluid model that provides a first approximation of the electric potential distribution. Steps 2–4 are used in the for-loop and they are described in detail in next sections. The computer model has been written in the MATLAB programming language and the solvers of differential equations of the COMSOL Multiphysics package were used to find the solution.
ARTICLE IN PRESS P. Bartos˘ et al. / Vacuum 82 (2008) 220–223
221
the electron energy density we are the dependent variables the system of Eqs. (1–4) is solved for. The intensity of ~ ¼ rj: Vectors J~e electric field is given by the equation E ~ and J i denote the particle flux that can be expressed (in drift-diffusion approximation) as a superposition of a drift under the influence of electric field and in the concentration gradient
Fig. 1. Model scheme.
The assumptions of our model can be summarized into several points: (1) A low-temperature argon plasma is considered, which consists of argon neutrals, positively charged ions and electrons. (2) Magnetic field effects are not considered. (3) The ion and neutral temperature is assumed to be constant. This kind of particles is much heavier than electrons and the energy loss in the collision with electron is negligible. (4) Elastic scattering, ionization of neutrals by fast electrons and excitation of neutrals are considered. Other collisions processes were found to be negligible, although they are playing an important role at higher pressure and higher energy of particles. 2.2. Fluid model The system of differential equations consists (see [4,5]) of continuity equations for both electrons and positively charged argon ions, Poisson’s equation and the energy balance for electrons qni þ r J~i ¼ kion ni nn , qt
(1)
qne þ r J~e ¼ kion ne nn , qt
(2)
eðni ne Þ , 0 qwe 5 ~ 5 De rwe þ r me we E 3 3 qt ~ ~ þ ej e E þ C ela þ H ion Rion þ H exc Rexc ¼ 0.
Dj ¼
~ De rne , J~e ¼ ne me E
(5)
~ Di rni . J~i ¼ ni mi E
(6)
The diffusion coefficients and coefficients of mobility are given in the form Dj ¼ ðkB T j Þ=ðmj vj;n Þ; mj ¼ e=ðmj vj;n Þ; j ¼ e; i. e denotes the elementary charge, e0 the permitivity in vacuum, kB is the Boltzmann constant, mj denotes the mass of particle, Tj its temperature and vj;n ¼ nn sj;n vj the mean collision frequency of particles with neutrals. The analytic formula of total cross-section sj,n can be found in [6]. Hexc ¼ 11.5 eV and Hion ¼ 15.8 eV are the energy losses in considered non-elastic collisions. Reaction rates on the left-hand side of the Eq. (4) are defined by the relations Rion ¼ njnnkion and Rexc ¼ njnnkexc, where Z 1 km ¼ f ðe Þsm ðe Þvðe Þde m ¼ ion; exc. (7) 0
The integration in the Eq. (7) have to be done over all electron energies ee with the energy distribution function f(ee) (EEDF). The energy loss in elastic electron-neutral collisions can be expressed by the Krook’s approximation [7] 2me 3 ve;n ne e kB T n . C ela ¼ (8) 2 mn Initial and boundary conditions are an essential part of the problem formulation. In our case, initial conditions were equal to the solution in the previous step. The boundary conditions are:
(3)
ð4Þ
In the above equations, subscripts e, i and n denote electrons, positive ions and neutrals, respectively. The concentrations of particles n, the electric potential j and
Probe surface: electric potential j ¼ 10 V, concentration of positive ions ni ¼ 0 m3, flux of electrons ~ Ge ~ n¼ ~~ me ne E n þ vth ne =4 (see [8]), electron energy density we ¼ 3ne kB T e =2 þ ej, where ~ n denotes an unit outward vector. Undisturbed plasma: electric potential j ¼ 0V, concentration of charged particles ni ¼ ne ¼ 1 1015 m3, electron energy density we ¼ 3nekBTe/2.
2.3. Particle model In the second step of the algorithm, there is solved a particle model in order to obtain the EEDF as a function of space coordinates. The distribution of electric field near the substrate is given externally as an output of the first part of the computation and, therefore, the particle model can be of the non self-consistent type. This approach allows us to obtain the EEDF quickly (the computation takes
ARTICLE IN PRESS 222
P. Bartos˘ et al. / Vacuum 82 (2008) 220–223
approximately one tenth of the time that is necessary in the self-consistent approach [9]) and all important collision processes can be described with the highest precision. The source of particles is imaginarily placed in the region, where the intensity of electric field is supposed to be zero. The motion of particles has been realized by solving of Newton’s equations by the Euler’s algorithm. The collision processes are calculated by the Monte Carlo (MC) method. The random free path is generated according to the equation [10] lrand ¼ ltot ln g at the beginning of the computation, ltot is the mean free path, gA(0, 1) is a uniformly distributed random number. Because the collision cross-sections are functions of particle energy, the null collision method (NCL) [11] has to be used. The collision arises, when the particle travels the distance equal to the random free path lrand. When the steady state is reached, the working area is divided into subdomains and the magnitudes of particle velocities are saved in regular time-steps. Finally, the EEDF can be determined as a function of space coordinates.
Table 1 Time expenditure of single algorithm steps (in s) Step
1D
2D
3D
Simplified fluid model Non self-consistent model Computation of coefficients Improved fluid model
6 E1500 51 8
56 E1800 53 67
540 E2000 62 660
10 8 6 4 2 0 2
2.4. Determination of coefficients In the third step of the algorithm there are determined the diffusion coefficients De and coefficients of mobility me, reaction rates kion and kexc and energy gain Ce,n in elastic collisions of electrons with neutrals. These values are functions of the mean collision frequency vj;n ¼ nn sj;n v: Therefore, the mean value sj;n v that is equal to the integral (Eq. (7)), has to be calculated. The EEDF is given as a look-up table on each subdomain and, therefore, the integral (Eq. 7) can be computed neither analytically nor numerically by any classical differential scheme. But the MC method can be used effectively in the following way: N random numbers v(r) with the EEDF are generated on (r) each subdomain and corresponding cross-sections sðrÞ j;n (v ) are determined. the value sj;n v can be computed: PFinally, ðrÞ sj;n v ¼ kj ¼ N1 r sðrÞ j;n v . In our case, N ¼ 500,000. 3. Results The presented computer model has been tested during computer simulations of argon plasma behaviour in the vicinity of planar (1D), cylindrical (2D) and spherical (3D) probes. There were used plasma parameters corresponding to the pressure 133 Pa, electron temperature Te ¼ 23, 200 K, neutral and positive ion temperature Tg ¼ Ti ¼ 300 K. The concentration of charged particles in the area of undisturbed plasma has been set to 1 1015 m3. Time consumption of these calculations is presented in Table 1. Some results are presented in Figs. 2 and 3 for illustration. Although the results compare well with those obtained by the computation of particle models, the presented model provides the solution faster, especially in 2D and 3D problems.
0
× 103
–2
–2
–1
0
1
2 × 103
Fig. 2. Distribution of electric potential in the vicinity of spherical probe (3D geometry, cross-section via the probe centre).
× 10–19 14
1D 2D 3D
12
10
8
6
4 0
0.2
0.4
0.6
0.8
1
Fig. 3. Distribution of mean electron energy—dependence on the distance from the planar probe surface (1D), from the axis of cylindrical probe (2D) and from the centre of spherical probe (3D).
4. Discussion and conclusion The presented computer model can be used in various branches of science. Although plasma is here presented as a continuum of three components (electrons, positive ions,
ARTICLE IN PRESS P. Bartos˘ et al. / Vacuum 82 (2008) 220–223
neutrals), other components (negatively charged ions, excited atoms, etc.) can be included easily, if necessary. Although only elastic scattering, excitation and ionization of neutral argon atoms are taken into account in the presented model, other collision processes do not occur so often and they do not affect the electron energy to a large extent (see [12]). The precision of calculations might be relatively high, because all important parameters are recalculated iteratively during the computation. In contrast, lots of authors determine the coefficients in the pre-computation step and they do not change them during the computation. In comparison with particle models, the presented algorithm provides the results much faster. Moreover, it allows us to study multidimensional problems with a complicated non-symmetric model geometry, where particle models cannot be used. Acknowledgement The work is a part of the research plan MSM 0021620834 financed by the Ministry of Education of the
223
Czech Republic. Authors P. Bartosˇ and P. Jelı´ nek thank to the financial support of the grant KAN 101 120 701.
References [1] Ellis HW, Pai RY, McDaniel EW, Mason EA, Viehland LA. At Data Nucl Data Tables 1976;17:177–210. [2] Bogaerts A, Gijbels R. J Anal At Spectrom 2000;15:1191–201. [3] Bogaerts A, Gijbels R, Goedheer WJ. J Appl Phys 1995;78: 2233–41. [4] Chen G, Raja L L. J Appl Phys 2004;96:6073–81. [5] Passchier JDP, Goedheer WJ. J Appl Phys 1993;74:3744–51. [6] Bogaerts A, Gijbels R, Goedheer W. J Appl Phys 1999;38: 4404–15. [7] Mitchner M, Kruger CH. Partially ionized gases. New York: Wiley; 1973. [8] Hagelaar GJM, Hoog FJ, Kroesen GMW. Phys Rew E 2000;62: 1452–4. [9] Hrach R. Czech J Phys 2001;51:557–66. [10] Sˇimek J, Hrach R. Czech J Phys 2006;56:B1086–90. [11] Skulerud HR. Null collision method. Technical report EIP 72-1, Trondheim, 1972. [12] Bogaerts A, Gijbels R. J Appl Phys 2002;92:6408–21.
Vacuum 82 (2008) 220–223 www.elsevier.com/locate/vacuum
Multidimensional fluid–particle modelling technique in low-temperature argon plasma at low pressure Petr Bartos˘ , Rudolf Hrach, Petr Jelı´ nek Department of Surface and Plasma Science, Faculty of Mathematics and Physics, Charles University, V Holesˇovicˇka´ch 2, 180 00, Prague 2, Czech Republic
Abstract In the contribution there is presented a new computer modelling technique, which can be used in multidimensional computer simulations of plasma at low pressure. The technique is based on the fluid–particle approach and can be used in the study of the range of physical phenomena, e.g. plasma behaviour in the vicinity of substrates, probe diagnostics, plasma technology, etc. The computer model is useable for the solution of three-dimensional problems and it minimizes the restrictions arising in complicated model geometry. The algorithm has been tested during the study of plasma sheath formation in the vicinity of probes with various geometries immersed into low-temperature argon plasma. Some results of these simulations are presented in the last section. r 2007 Elsevier Ltd. All rights reserved. PACS: 52.65.y; 52.40.Kh Keywords: Fluid–particle simulations; Plasma sheath; Argon plasma
1. Introduction
2. Model description
In recent years plasma is used widely in science, industry and technology. Scientific teams endeavour to optimize the plasma parameters of devices and so achieve the best performance of production lines. Unfortunately, this is often a very expensive and time-consuming activity. In the last decade the performance of computer systems has increased rapidly and therefore a new and less timeconsuming approach, computer modelling, can be used. The plasma parameters in technological devices can be predicted by this way with high precision, faster and more cheaply. On our department an effective computer model of low-temperature plasma has been developed. This model combines the advantages of both the fluid modelling technique (fast computation and therefore the possibility of the computer modelling of multidimensional problems) and particle modelling approach (precision in the description of collision processes in the plasma).
2.1. Model philosophy
Corresponding author.
E-mail address: [email protected] (P. Bartos˘ ). 0042-207X/$ - see front matter r 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.vacuum.2007.07.006
Models based on the fluid modelling technique are often used in recent studies. Unfortunately, the precision of these models depends on many parameters (diffusion coefficients and coefficients of mobility, collision frequencies, coefficient rates, etc.), which are very often complicated functions of energy. These dependences can be determined by experimental measurements and expressed as semiempirical formulas, e.g. Ref. [1], but this approach is too complicated. Alternatively, fluid–particle modelling techniques can be used, e.g. Refs. [2,3]. We developed a model of such type, the model arrangement is presented in Fig. 1. The algorithm consists of four separate parts. In the first step, there is solved a simplified fluid model that provides a first approximation of the electric potential distribution. Steps 2–4 are used in the for-loop and they are described in detail in next sections. The computer model has been written in the MATLAB programming language and the solvers of differential equations of the COMSOL Multiphysics package were used to find the solution.
ARTICLE IN PRESS P. Bartos˘ et al. / Vacuum 82 (2008) 220–223
221
the electron energy density we are the dependent variables the system of Eqs. (1–4) is solved for. The intensity of ~ ¼ rj: Vectors J~e electric field is given by the equation E ~ and J i denote the particle flux that can be expressed (in drift-diffusion approximation) as a superposition of a drift under the influence of electric field and in the concentration gradient
Fig. 1. Model scheme.
The assumptions of our model can be summarized into several points: (1) A low-temperature argon plasma is considered, which consists of argon neutrals, positively charged ions and electrons. (2) Magnetic field effects are not considered. (3) The ion and neutral temperature is assumed to be constant. This kind of particles is much heavier than electrons and the energy loss in the collision with electron is negligible. (4) Elastic scattering, ionization of neutrals by fast electrons and excitation of neutrals are considered. Other collisions processes were found to be negligible, although they are playing an important role at higher pressure and higher energy of particles. 2.2. Fluid model The system of differential equations consists (see [4,5]) of continuity equations for both electrons and positively charged argon ions, Poisson’s equation and the energy balance for electrons qni þ r J~i ¼ kion ni nn , qt
(1)
qne þ r J~e ¼ kion ne nn , qt
(2)
eðni ne Þ , 0 qwe 5 ~ 5 De rwe þ r me we E 3 3 qt ~ ~ þ ej e E þ C ela þ H ion Rion þ H exc Rexc ¼ 0.
Dj ¼
~ De rne , J~e ¼ ne me E
(5)
~ Di rni . J~i ¼ ni mi E
(6)
The diffusion coefficients and coefficients of mobility are given in the form Dj ¼ ðkB T j Þ=ðmj vj;n Þ; mj ¼ e=ðmj vj;n Þ; j ¼ e; i. e denotes the elementary charge, e0 the permitivity in vacuum, kB is the Boltzmann constant, mj denotes the mass of particle, Tj its temperature and vj;n ¼ nn sj;n vj the mean collision frequency of particles with neutrals. The analytic formula of total cross-section sj,n can be found in [6]. Hexc ¼ 11.5 eV and Hion ¼ 15.8 eV are the energy losses in considered non-elastic collisions. Reaction rates on the left-hand side of the Eq. (4) are defined by the relations Rion ¼ njnnkion and Rexc ¼ njnnkexc, where Z 1 km ¼ f ðe Þsm ðe Þvðe Þde m ¼ ion; exc. (7) 0
The integration in the Eq. (7) have to be done over all electron energies ee with the energy distribution function f(ee) (EEDF). The energy loss in elastic electron-neutral collisions can be expressed by the Krook’s approximation [7] 2me 3 ve;n ne e kB T n . C ela ¼ (8) 2 mn Initial and boundary conditions are an essential part of the problem formulation. In our case, initial conditions were equal to the solution in the previous step. The boundary conditions are:
(3)
ð4Þ
In the above equations, subscripts e, i and n denote electrons, positive ions and neutrals, respectively. The concentrations of particles n, the electric potential j and
Probe surface: electric potential j ¼ 10 V, concentration of positive ions ni ¼ 0 m3, flux of electrons ~ Ge ~ n¼ ~~ me ne E n þ vth ne =4 (see [8]), electron energy density we ¼ 3ne kB T e =2 þ ej, where ~ n denotes an unit outward vector. Undisturbed plasma: electric potential j ¼ 0V, concentration of charged particles ni ¼ ne ¼ 1 1015 m3, electron energy density we ¼ 3nekBTe/2.
2.3. Particle model In the second step of the algorithm, there is solved a particle model in order to obtain the EEDF as a function of space coordinates. The distribution of electric field near the substrate is given externally as an output of the first part of the computation and, therefore, the particle model can be of the non self-consistent type. This approach allows us to obtain the EEDF quickly (the computation takes
ARTICLE IN PRESS 222
P. Bartos˘ et al. / Vacuum 82 (2008) 220–223
approximately one tenth of the time that is necessary in the self-consistent approach [9]) and all important collision processes can be described with the highest precision. The source of particles is imaginarily placed in the region, where the intensity of electric field is supposed to be zero. The motion of particles has been realized by solving of Newton’s equations by the Euler’s algorithm. The collision processes are calculated by the Monte Carlo (MC) method. The random free path is generated according to the equation [10] lrand ¼ ltot ln g at the beginning of the computation, ltot is the mean free path, gA(0, 1) is a uniformly distributed random number. Because the collision cross-sections are functions of particle energy, the null collision method (NCL) [11] has to be used. The collision arises, when the particle travels the distance equal to the random free path lrand. When the steady state is reached, the working area is divided into subdomains and the magnitudes of particle velocities are saved in regular time-steps. Finally, the EEDF can be determined as a function of space coordinates.
Table 1 Time expenditure of single algorithm steps (in s) Step
1D
2D
3D
Simplified fluid model Non self-consistent model Computation of coefficients Improved fluid model
6 E1500 51 8
56 E1800 53 67
540 E2000 62 660
10 8 6 4 2 0 2
2.4. Determination of coefficients In the third step of the algorithm there are determined the diffusion coefficients De and coefficients of mobility me, reaction rates kion and kexc and energy gain Ce,n in elastic collisions of electrons with neutrals. These values are functions of the mean collision frequency vj;n ¼ nn sj;n v: Therefore, the mean value sj;n v that is equal to the integral (Eq. (7)), has to be calculated. The EEDF is given as a look-up table on each subdomain and, therefore, the integral (Eq. 7) can be computed neither analytically nor numerically by any classical differential scheme. But the MC method can be used effectively in the following way: N random numbers v(r) with the EEDF are generated on (r) each subdomain and corresponding cross-sections sðrÞ j;n (v ) are determined. the value sj;n v can be computed: PFinally, ðrÞ sj;n v ¼ kj ¼ N1 r sðrÞ j;n v . In our case, N ¼ 500,000. 3. Results The presented computer model has been tested during computer simulations of argon plasma behaviour in the vicinity of planar (1D), cylindrical (2D) and spherical (3D) probes. There were used plasma parameters corresponding to the pressure 133 Pa, electron temperature Te ¼ 23, 200 K, neutral and positive ion temperature Tg ¼ Ti ¼ 300 K. The concentration of charged particles in the area of undisturbed plasma has been set to 1 1015 m3. Time consumption of these calculations is presented in Table 1. Some results are presented in Figs. 2 and 3 for illustration. Although the results compare well with those obtained by the computation of particle models, the presented model provides the solution faster, especially in 2D and 3D problems.
0
× 103
–2
–2
–1
0
1
2 × 103
Fig. 2. Distribution of electric potential in the vicinity of spherical probe (3D geometry, cross-section via the probe centre).
× 10–19 14
1D 2D 3D
12
10
8
6
4 0
0.2
0.4
0.6
0.8
1
Fig. 3. Distribution of mean electron energy—dependence on the distance from the planar probe surface (1D), from the axis of cylindrical probe (2D) and from the centre of spherical probe (3D).
4. Discussion and conclusion The presented computer model can be used in various branches of science. Although plasma is here presented as a continuum of three components (electrons, positive ions,
ARTICLE IN PRESS P. Bartos˘ et al. / Vacuum 82 (2008) 220–223
neutrals), other components (negatively charged ions, excited atoms, etc.) can be included easily, if necessary. Although only elastic scattering, excitation and ionization of neutral argon atoms are taken into account in the presented model, other collision processes do not occur so often and they do not affect the electron energy to a large extent (see [12]). The precision of calculations might be relatively high, because all important parameters are recalculated iteratively during the computation. In contrast, lots of authors determine the coefficients in the pre-computation step and they do not change them during the computation. In comparison with particle models, the presented algorithm provides the results much faster. Moreover, it allows us to study multidimensional problems with a complicated non-symmetric model geometry, where particle models cannot be used. Acknowledgement The work is a part of the research plan MSM 0021620834 financed by the Ministry of Education of the
223
Czech Republic. Authors P. Bartosˇ and P. Jelı´ nek thank to the financial support of the grant KAN 101 120 701.
References [1] Ellis HW, Pai RY, McDaniel EW, Mason EA, Viehland LA. At Data Nucl Data Tables 1976;17:177–210. [2] Bogaerts A, Gijbels R. J Anal At Spectrom 2000;15:1191–201. [3] Bogaerts A, Gijbels R, Goedheer WJ. J Appl Phys 1995;78: 2233–41. [4] Chen G, Raja L L. J Appl Phys 2004;96:6073–81. [5] Passchier JDP, Goedheer WJ. J Appl Phys 1993;74:3744–51. [6] Bogaerts A, Gijbels R, Goedheer W. J Appl Phys 1999;38: 4404–15. [7] Mitchner M, Kruger CH. Partially ionized gases. New York: Wiley; 1973. [8] Hagelaar GJM, Hoog FJ, Kroesen GMW. Phys Rew E 2000;62: 1452–4. [9] Hrach R. Czech J Phys 2001;51:557–66. [10] Sˇimek J, Hrach R. Czech J Phys 2006;56:B1086–90. [11] Skulerud HR. Null collision method. Technical report EIP 72-1, Trondheim, 1972. [12] Bogaerts A, Gijbels R. J Appl Phys 2002;92:6408–21.