Comparison of Analytical and Numerical Solutions and PIC

Nuclear Energy for New Europe 2006

Portorož, Slovenia, September 18-21, 2006 http://www.djs.si/port2006

Potential Formation in a Plasma Diode Containing Two-electron Temperature Plasma – Comparison of Analytical and Numerical

Solutions and PIC Simulations

Tomaž Gyergyek^1,2, Milan Čerček^2,3, Borut Jurčič-Zlobec¹

1) University of Ljubljana, Faculty of electrical engineering, Tržaška 25, SI-1000 Ljubljana, Slovenia

2) “Jožef Stefan” Institute Jamova 39, SI-1000 Ljubljana, Slovenia

3) University of Maribor, Faculty of civil engineering, Smetanova 17, SI-2000 Maribor, Slovenia

Tomaz.gyergyek@fe.uni-lj.si, Milan.Cercek@ijs.si ABSTRACT

Potential formation in front of a negatively biased electron emitting electrode (collector) is studied by a fully kinetic one dimensional model of a plasma diode that contains electrons with a two-temperature maxwellian velocity distribution function. Analytical and numerical results of the model are compared to the results of the PIC simulation.

1 INTRODUCTION

In the edge plasmas of fusion devices energetic electron populations are readily produced. The presence of energetic electrons has a remarkable effect on the potential formation in the plasma and consequently on particle losses to the wall. Also electron emission from the walls that limit the plasma is a common occurrence in fusion devices. In this work we present a fully kinetic model of plasma potential formation in a bounded plasma system containing energetic and emitted electrons and compare the results of the model with numerical solutions of the Poisson equation and with the computer simulation.

2 MODEL

We consider a large planar electrode (collector) with its surface perpendicular to the x axis of the coordinate system. The collector is located at x = 0. This electrode absorbs all the particles that hit it. On the other hand it may also emit electrons. This electron emission may be thermal or secondary, triggered by the impact of incoming electrons and/or ions. The details of the emission mechanism are not essential for the model.

An infinitely large planar plasma source also has its surface perpendicular to the x axis.

The source is located at a certain distance x = L from the collector. The distance L is not essential for the model and may be taken arbitrarily. The source injects 3 groups of charged particles into the system: singly charged positive ions (index i), the cool electrons (index 1), the hot electrons (index 2). The electrons that are emitted from the collector have index 3. The particles i, 1 and 2 are injected from the source with half-maxwellian velocity distribution function with respective temperatures Ti, T1 and T2. The emitted electrons also have a half- maxwellian velocity distribution function at the collector with the temperature T3. We assume

that T2 > T1 and T1 >> T3. The ion temperature Ti can in principle be arbitrary, but is usually taken smaller than or equal to T1.

The collector is biased to a certain (negative!) potential ΦC. The potential an the electric field at the source are set to zero. This imposes the following boundary conditions:

( ) 0, d ( ) 0,

x L x L

dx

Φ = = Φ = = (1)

for the Poisson equation:

2

(

0

1 2 3

2

0

( ) ( ) ( ) ( ) .

i

e

d n x n x n x n x

dx ε

Φ= − − − −

)

⁽²⁾

where e0 is the elementary charge and nj(x) are particle densities. We are looking only for such solutions of the Poisson equation (2), where Φ(x) is a monotonically decreasing function in the direction towards the collector.

An ion that has left the source with a negligibly small initial velocity has at the distance x < L from the collector the velocity:

2 0 ( )

mi .

i

e x

v m

= − Φ (3)

The distribution function for the ions can therefore be written in the following way:

2

0 ( )

exp exp ( ).

2 2

i i

i i mi

i i i

m e x m v

f n

kT kT kT

π

⎛ Φ ⎞ ⎛ ⎞

= ⎜− ⎟ ⎜− ⎟

⎝ ⎠ ⎝ ⎠H v v− (4)

Here ni is the ion density at the source. Velocity direction towards the collector is positive and away from the collector is negative.

An electron that has left the collector source with a negligibly small initial velocity has at the distance x from the collector the velocity:

( )

2 0 ( )

C .

me

e

e x

v m

Φ − Φ

= − (5)

So the distribution functions for the electrons are written in the following way:

2 0

1 1

1 1 1

exp ( ) exp ( ),

2 2

e e

me

m e x m v

f n

kT kT kT

π

⎛ Φ ⎞ ⎛ ⎞

= ⎜ ⎟ ⎜− ⎟

⎝ ⎠ ⎝ ⎠H v v− (6)

2 0

2 2

2 2 2

exp ( ) exp ( ),

2 2

e e

me

m e x m v

f n

kT kT kT

π

⎛ Φ ⎞ ⎛ ⎞

= ⎜ ⎟ ⎜− ⎟

⎝ ⎠ ⎝ ⎠H v v− (7)

( )

²

0

3 3

3 3 3

exp ( ) exp ( ).

2 2

e C

me

e x

m m_ev

f n H

kT kT kT

π

Φ − Φ

⎛ ⎞ ⎛ ⎞

= ⎜ ⎟ ⎜ ⎟

⎝ ⎠

⎝ ⎠ − v −v (8)

Here n1, n2, and n3 are the respective electron densities at the source and at the collector.

We now introduce the following variables:

3 0 0 0

2

1 1 1 1 1

3 0

2 1

0 2

1 1 1 0 1 0

( ) ( 0)

, , , , , ,

, , , 2 , , , .

e i C

C i

i

D

e D

m T T T e x e x e

m T T T kT kT kT

n n n kT v x

v u z

n n n m v n e

μ τ σ

α β ε λ ε

λ

1

kT1

Φ Φ = Φ

= = Θ = = Ψ = Ψ = =

= = = = = = =

(9)

With these variables the distribution functions are written in the following way:

(

^,

)

^{exp exp 2}

(

i

F u α u H u μ

τ μτ

π τ μ

⎛ ⎞

⎛ Ψ⎞

)

^,

Ψ = ⎜− ⎟ ⎜− ⎟ −

⎝ ⎠ ⎝ ⎠ − Ψ (10)

( ) ( ) ( )

²

( )

1

, 1 exp exp _C

F u u H u

Ψ = π Ψ − + Ψ − Ψ , (11)

( )

²

(

2 , exp exp u _C

F u β H u

π

⎛ ⎞

⎛Ψ⎞

Ψ = Θ ⎜⎝Θ⎟⎠ ⎜⎝−Θ ⎟⎠ + Ψ − Ψ

)

^{, (12)}

( )

²

(

3 , exp ^C exp u _C

)

.

F u ε

σ σ

π σ

⎛ ⎞

Ψ − Ψ

⎛ ⎞

Ψ = ⎜⎝ ⎟⎠ ⎜⎝− ⎟⎠H − Ψ − Ψ −u (13)

The n-th moment of the velocity distribution is obtained in the following way. The distribution function is multiplied by the selected power of the velocity uⁿ and then integrated over the velocity space. In this way the moment of the distribution function is obtained as a function of the potential Ψ. The zero moments (n = 0) are particle densities and the first moments (n = 1) are fluxes. For the distribution functions (10) - (13) the first and zero moments are the following:

( ) ( )

( ) ( ) ( ) ^{( ) ( )}

( ) ( )

( )

1 1

2 2

3

exp Erfc , ,

2 2

1 1

exp 1 Erf , exp ,

2 2

exp 1 Erf , exp ,

2 2

exp Erfc

2

i i

C C

C C

C C

N J

N J

N J

N

α α

τ τ π

π

β β

π ε

σ σ

⎛ ⎞

Ψ Ψ

⎛ ⎞

Ψ = ⎜⎝− ⎟⎠ ⎜⎜⎝ − ⎟⎟⎠ Ψ =

Ψ = Ψ + Ψ − Ψ Ψ = Ψ

⎛ Ψ − Ψ ⎞ Ψ

Ψ Θ ⎛ ⎞

⎛ ⎞

Ψ = ⎜⎝Θ⎟⎠⎜⎜⎝ + Θ ⎟⎟⎠ Ψ = ⎜⎝ ⎟⎠

⎛ ⎞

Ψ − Ψ Ψ − Ψ

⎛ ⎞

Ψ = ⎜⎝ ⎟⎠ ⎜⎜⎝ ⎟⎟⎠

( )

μτ

Θ

, 3 ,

J ε σ2

Ψ = − π

(14)

where

( ) ( )

²

^{( )} ( )

0

2 2

Erf exp , Erfc exp .

x

x

x t dt x t dt2

π π

∞

=

∫

− =

∫

− ⁽¹⁵⁾

The total current density Jt to the collector is then given by:

( )

3 1 2

| | 1 exp exp .

2 2 2 2

C

t i C

J J J J J α μτ ε σ β

π π π π

Ψ

Θ ⎛

= + − − = + − Ψ − ⎜ Θ⎝ ⎠

⎞⎟ (16) With the variables (9) and particle densities given by (14) the Poisson equation (2) is written in the following form:

( ) ( )

( ( ) ) ( ( ) ) ^{( ) ( )}

( ) ( )

2

2 exp Erfc

+exp 1 Erf exp 1 Erf

exp Erfc .

C C

C C

z z

d dz

z z

z z

z z

α τ τ

β

ε σ σ

⎛ ⎞

Ψ Ψ

⎛ ⎞

Ψ = − ⎜⎝− ⎟⎠ ⎜⎜⎝ − ⎟⎟⎠+

⎡ ⎛ ⎞⎤

Ψ Ψ −

⎛ ⎞

⎡ ⎤ Ψ

⎢ ⎜ ⎟⎥

Ψ ⎢⎣ + Ψ − Ψ ⎥⎦+ ⎜ Θ ⎟ + ⎜ Θ ⎟ +

⎢ ⎥

⎝ ⎠⎣ ⎝ ⎠⎦

⎛ ⎞

Ψ − Ψ Ψ − Ψ

⎛ ⎞

⎜ ⎟

+ ⎜⎝ ⎟⎠ ⎜⎝ ⎟⎠

(17)

At a certain distance z = z0 from the collector the plasma is neutral and the potential at that point has the value Ψ(z = z0) = ΨP. So the neutrality condition can be obtained directly from (17):

( ) ( )

0

2

2 exp Erfc

+exp 1 Erf exp 1 Erf

exp Erfc 0.

P P

z z

P C

P

P P C

P C P C

d

dz α

τ τ

β

ε σ σ

=

⎛ ⎞

⎛ Ψ⎞ = − ⎛−Ψ ⎞ ⎜ −Ψ ⎟+

⎜ ⎟ ⎜⎝ ⎟⎠ ⎜ ⎟

⎝ ⎠ ⎝ ⎠

⎡ ⎛ Ψ − Ψ ⎞⎤

⎛Ψ ⎞

⎡ ⎤

Ψ ⎣ + Ψ − Ψ ⎦+ ⎜⎝ Θ ⎟⎠⎢⎢⎣ + ⎜⎜⎝ Θ ⎟⎟⎠⎥⎥⎦+

⎛ ⎞

Ψ − Ψ Ψ − Ψ

⎛ ⎞

+ ⎜⎝ ⎟⎠ ⎜⎜⎝ ⎟⎟⎠=

(18)

At z = z0 the potential Ψ(z) obviously has an inflection point. Since:

2 2

2

1 ,

2

d d d d

dz dz dz dz

Ψ Ψ Ψ

⎛ ⎞ =

⎜ ⎟

⎝ ⎠

the Poisson equation (17) is multiplied by dΨ/dz and integrated once over Ψ from Ψ = 0 (at the source) to Ψ = ΨP (at the inflection point). The following equation is obtained:

( ) ( ) ( )

1 2 exp Erfc

+exp 1 Erf 2 exp 1 Erf

exp 1 Erf 2 exp 1 Erf

2

P P P

P P C C P C C C

P C C P C C C

P

P C

ατ π τ τ τ

π

β π

εσ π σ

⎡ Ψ ⎛ Ψ ⎞ Ψ ⎤

− ⎢⎣ − − − ⎜⎝− ⎟⎠ − ⎥⎦+

⎡ ⎤ ⎡ ⎤

Ψ ⎣ + Ψ − Ψ ⎦− Ψ Ψ − Ψ − −Ψ − +⎣ −Ψ ⎦+

⎡ ⎛Ψ ⎞⎛ Ψ − Ψ ⎞ ⎛Ψ ⎞⎛ Ψ − Ψ −Ψ ⎞ ⎛ −Ψ ⎤ + Θ⎢⎢⎣ ⎜⎝ Θ ⎟⎠⎜⎜⎝ + Θ ⎟⎟⎠− ⎜⎝ Θ ⎟⎠⎜⎜⎝ Θ − Θ ⎟ ⎜⎟ ⎜⎠ ⎝− + Θ ⎥⎥⎦+

Ψ − Ψ −Ψ

+ −

⎞⎟⎟

⎠

exp Erfc exp Erfc 0.

C P C P C C C

σ σ σ σ σ

⎡ ⎛ ⎞+ ⎛Ψ − Ψ ⎞ Ψ − Ψ − ⎛−Ψ ⎞ Ψ ⎤

⎢ ⎜⎜ ⎟⎟ ⎜⎝ ⎟⎠ ⎜⎝ ⎟⎠ ⎥

⎢ ⎝ ⎠ ⎥

⎣ − ⎦=

(19) This equation is called the zero electric field condition at the inflection point.

If the electron emission from the collector increases, at a certain level of emission a negative space charge may start to accumulate at the collector surface. At certain critical electron emission the electric field that accelerates the negative electrons away from the collector becomes equal to zero. The zero electric field condition at the collector can be derived in a very similar way as the zero electric field condition at the inflection point (19).

Again the Poisson equation (17) is multiplied by dΨ/dz and integrated once over Ψ, only now the integration limits go from Ψ = ΨP (at the inflection point) to Ψ = ΨC (at the collector).

The following equation is obtained:

( ) ( )

exp Erfc exp Erfc 2

+exp 1 2 exp 1 Erf

exp 1 2 exp 1 Erf

C C

P P P

C P C P P C

C P C P P C

ατ τ τ τ τ π τ τ

π

β π

⎡ ⎛ Ψ ⎞ ⎛ Ψ ⎞ ⎛ Ψ ⎞ ⎛ Ψ ⎞ ⎛ Ψ Ψ ⎞⎤

− ⎢⎢⎣ ⎜⎝− ⎟⎠ ⎜⎜⎝ − ⎟⎟⎠− ⎜⎝− ⎟⎠ ⎜⎜⎝ − ⎟⎟⎠+ ⎜⎜⎝ − − − ⎟⎟⎠⎥⎥⎦

⎛ ⎞ ⎡ ⎤

Ψ ⎜⎝ + Ψ − Ψ ⎟⎠− Ψ ⎣ + Ψ − Ψ ⎦+

⎡ ⎛Ψ ⎞⎛ Ψ − Ψ ⎞ ⎛Ψ ⎞⎡ Ψ − Ψ ⎤ + Θ⎢⎢⎣ ⎜⎝ Θ ⎟⎜⎠⎜⎝ + Θ ⎟⎟⎠− ⎜⎝ Θ ⎟⎠⎢⎣ + Θ ⎥⎦

C +

1 2 ^{P C} exp ^{P C} Erfc ^{P C} 0.

εσ π σ σ σ

⎤+

⎥⎥⎦

⎡ Ψ − Ψ ⎛Ψ − Ψ ⎞ Ψ − Ψ ⎤

+ ⎢ − − ⎜ ⎟ ⎥=

⎝ ⎠

⎣ ⎦

(20)

In an experiment the parameters like μ, Θ, τ, σ and β are determined by the method of the plasma production. The potential of the collector ΨC can be determined by the external power supply connected to it. If the collector is floating, then the potential ΨC is determined by the floating condition of zero total current to the collector. This means that Jt given by (16) is set to zero. The electron emission from the collector (parameter ε) can be considered as given by the experimental set-up, if the emission is thermal and temperature limited. If it is space charge limited, then the zero electric field condition at the collector (20) is fulfilled. To summarize: if the collector is floating and the emission is critical, equations (16) (with Jt = 0), (18), (19) and (20) form a system of 4 equations for 4 unknown quantities: α, ΨP, ΨC and the critical value of ε for a given set of the parameters μ, Θ, τ, σ and β. If the collector potential ΨC is also given, we have to solve the system of 3 equations (18), (19) and (20) for 3 unknown quantities: α, ΨP and the critical value of ε. If also the emission ε is determined (e.

g. by the collector temperature) and it is below the space charge limit then we only have to solve 2 equations (18) and (19) for α and ΨP. The potential profile Ψ(z) can be found by solving the Poisson equation (17) for a given set of the parameters α, Θ, τ, ε, σ and β. This can only be done numerically. The collector potential ΨC and the derivative dΨ/dz (electric field at the collector) must in this case be given as the boundary conditions for the Poisson equation (17).

3 RESULTS

In this paper we are only interested in a floating collector with a given (temperature limited) electron emission. Our main objective is to compare numerical solutions of the Poisson equation (17) with the analytical predictions of the model and with the results of the PIC simulations. We first select the following parameters: μ = 1/1836 (proton mass), τ = 0.1, β = 0.31, Θ = 20 and σ = 0.01 and then we solve the system of equations (16) (with Jt = 0), (18), (19) and (20). In this case the system of equation (16), (18), (19) and (20) has 3 solutions which are the following: 1) ΨP = -1.75731, ΨC = -22.8964, α = 7.0144, ε = 3.89487; 2) ΨP = -7.62295, ΨC = -28.0594, α = 6.67414, ε = 2.91603; 3) ΨP = -13.2845, ΨC = -27.6732, α = 6.69849, ε = 2.98069. The solutions are most distinguished by the value of ΨP. We call the solution with the ΨP closest to zero the low solution. The solution with the most negative ΨP

is called the high solution and solution with the intermediate ΨP is called the middle solution.

The nature of these solutions is discussed elsewhere [1,2]. Let us just mention that the low and the high solution correspond to the cases, when the plasma potential is determined by the cool or by the hot electron population, while the middle solution is not a physical one. We now focus on the comparison with numerical solutions. We proceed in the following way: the same parameters as before are selected: μ = 1/1836, τ = 0.1, β = 0.31, Θ = 20 and σ = 0.01.

The emission coefficient ε is now also selected and gradually increased, but always kept below the critical value, found in the above solutions. The system of equations (16), (18) and (19) is then solved for each set of parameters. The solutions are given in Table 1. In this way α, ΨP and ΨC are found. Then α is inserted into (17) together with other parameters and ΨC is used as the boundary condition. In order to solve the equation (17) we need an additional boundary condition, which is the value of the derivative dΨ/dz ≡ ΨC|. This boundary condition is found by guessing. Initial value for ΨC| is usually taken close to 1. Then the equation (17) is solved numerically using the 4-th order Runge-Kutta method and the graph of the solution is plotted. From the graph it can be seen easily, whether the selected value for ΨC|

is to large or to small and it is corrected. When gradually more and more decimal places of ΨC| are found the region where Ψ(z) has (almost) a constant value becomes longer and longer.

Table 1: Solutions of the system of equations (16), (18) and (19) for μ = 1/1836, τ = 0.1, β = 0.31, Θ = 20 and σ = 0.01, while ε is gradually increased.

ε(t) ΨC ΨC| ΨP(t) ΨP(n) α(t) α(n) β(n) ε(n) 0 -65.8719 3.42852260437 -1.7252 -1.73 6.9731 3.33 0.31 0 -66.9065 3.2098160816 -8.6125 -2.13 6.6215

0 -66.901 3.16628121 -10.9045 -10.9 6.6233

0.5 -52.2937 2.8041197712 -1.7259 -1.73 6.9740 3.31 0.31 0.25 0.5 -52.81 2.49531662 -8.55904 -11.0 6.62359

0.5 -52.8069 2.4963949 -10.9846 -10.99 6.6257

1.0 -44.2789 2.3574303849 -1.7276 -1.73 6.9763 3.36 0.31 0.5 1.0 -44.6204 2.00887947 -8.43068 -11.26 6.62882

1.0 -44.6173 2.010586639 -11.1898 -11.19 6.63193

1.5 -38.5718 1.98704831267 -1.7302 -1.73 6.9798 2.93 0.30 0.75 1.5 -38.8247 1.5927055 -8.25575 -11.61 6.63674

1.5 -38.821 1.59568307 -11.5009 -11.50 6.64171

2.0 -34.1362 1.6550131368 -1.73367 -1.73 6.98449 3.48 0.30 1.0 2.0 -34.335 1.197958588 -8.05042 -12.08 6.64734

2.0 -34.3302 1.203525831 -11.9245 -11.93 6.6555

Figure 1: Numerical solutions Ψ(z) of the Poisson equation (17), the negative derivative (the electric field) -dΨ/dz, the particle densities and ΔN for the following parameters: μ = 1/1836, τ = 0.1, β = 0.31, Θ = 20, σ = 0.01 and ε = 2 are plotted versus z. In the left column the low solution is shown, in the middle column the middle solution is shown and in the right

column the high solution is shown.

In Fig. 1 we show a plot of Ψ(z) versus z of the numerical solutions that correspond to the low, middle and high solution of the system of equations (16), (18) and (19) for = 1/1836, τ = 0.1, β = 0.31, Θ = 20, σ = 0.01 and ε = 2. The corresponding values of ΨC, ΨC|, ΨP and α are also given in the Table. In the Table 1 the values of ΨP and α that are obtained from the system of equations (16), (18) and (19) are labelled ΨP(t) and α(t). In separate columns we give the values of ΨP(n) and α(n). These are obtained from the numerical solutions of the equation (17). Note that for the low and the high solution the theoretical value ΨP(t) and the numerical value ΨP(n), which is simply read from the constant portion of the Ψ(z) versus z graph (see Fig. 1) are in very good agreement. For the middle solution the difference between ΨP(t) and ΨP(n) is considerable. In Fig. 1 we also show the electric field profile -dΨ/dz, the density profiles of all 4 particle species found from (14) and ΔN = Ni – Ne1 – Ne2 – Ne3. The potential profile Ψ(z) has a local maximum value for the low, the middle and the high solution. But only for the low solution the value of this maximum comes very close to zero and therefore (almost) fulfills the boundary condition (1). So for the low solution this point could be identified as the position of the plasma source. This position therefore depends on the precision of ΨC|.

Figure 2: Comparison of PIC simulation (left) and numerical solution of (17) (right).

For the low solution the particle densities at the source can be determined and so also the ratios α(n) and β(n) and ε(n) can be determined and compared to the corresponding theoretical values α(t) and β(t) and ε(t). For β the agreement is very good, for ε the difference is exactly a factor of 2 and for α approximately a factor of 2.

In Fig. 2 we show a comparison between a result of a computer PIC simulation obtained by the XPDP1 code [3] and a numerical solution of (17). In the simulation τ and Θ at the

source and σ at the collector can be selected independently. The other independently selected input parameters are influxes of the particle species from the source and from the collector.

These influxes must be selected so, that the electric field at the source is zero. Then α, β and ε are obtained as results of the simulation, together with particle densities, collector floating potential and other parameters. In the simulation shown in Fig. 2 we have selected kT1 = 1 eV, kT2 = 20 eV, kTi = 0.1 eV, so that τ = 0.1, Θ = 20, ε = 0 and μ = 1/1838.132. Then we

have found that the collector floating potential is -67.3916 V, which gives directly ΨC = -67.3916. From the simulation we also find β = 0.3562. With these parameters we solve

the system of equations (16) and (18) to get α = 7.43032 and ΨP = -1.95799. The later result is in good agreement with simulation. Then we solve the equation (17) numerically for the parameters that τ = 0.1, Θ = 20, ε = 0, μ = 1/1838.132, ΨC = -67.3916 and α = 7.43032.

Qualitative and also quantitative agreement of the simulation and the numerical solution is very good. The ratio between the ion and the cool electron density at the source in the computer simulation is αsim = 2.97, while for the numerical solution it is α(n) = 2.92.

4 CONCLUSIONS

We have compared an analytical fully kinetic model of a bounded plasma system with numerical solutions of the corresponding Poisson equation and the results of a PIC simulation.

We found very good agreement in the predicted plasma potential. There is also a very good agreement between the numerical solution and the computer simulation in the predicted ion density at the source.

ACKNOWLEDGMENTS

This work been carried out within the Association EURATOM-MHEST. The content of the publication is the sole responsibility of its authors and it does not necessarily represent the view of the Commission or its services.

REFERENCES

[1] T. Gyergyek, M. Čerček, B. Jurčič-Zlobec, Czechoslovak journal of physics, vol. 56, (2006), suppl. B, pp. B733-B739J.

[2] T. Gyergyek, M. Čerček, Czechoslovak Journal of Physics, vol. 54, (2004), pp. 431-460 [3] J. P. Verbocoeur, M. V. Alves, V. Vahedi and C. K. Birdsall, J. Comput. Phys., vol. 104,

(1993), pp. 321-328