# «Association EURATOM / IPP.CR I N S T I T U T E O F P L A S M A P H Y S I C S, v.v. i. ACADEMY OF SCIENCES OF THE CZECH REPUBLIC ANNUAL REPORT ...»

We thus address here a persistent problem in the interpretation of probe data, which is the inconsistency between theoretical calculations of fast electron distributions that should give extremely large sheath potentials (~kV) and measurements, which only give 0-100 V. The reason might be due to secondary emission - the s.e.e. coefficient is of the order of 1 for energies in the 400 eV and higher range. That means that fast electrons do not contribute to the net electrical current, while only the thermal ones do. So instead of a very large potential drop with a very energetic electron flux reflected from the sheath back towards the SOL (what we get now), we would expect to have a reasonable potential drop with cold electrons sent back into the SOL.

As a basis for our kinetic treatment of secondary electron emission we adopt the semi-empirical form for the secondary electron emission – the s.e.e. - coefficient δ [1,2]

where δmax and Emax are respectively the maximum value of δ and the energy value where this maximum occurs. Further, Ekin is the electron kinetic energy, Vw0 is its potential energy, and E ⊥ its energy of gyromotion. The s.e.e coefficient thus normalized is shown in Fig.1. and selected values of δmax and Emax are shown in Table 1, both taken from Ref. [1].

ANNUAL REPORT 2007 ASSOCIATION EURATOM/IPP.CR Part II - PHYSICS The δ-curve decreases for EEmax owing to deeper penetration of the incident electron into the target. It can be seen that for most tokamak edge electron energies (i.e. thermals and LHgenerated fast electrons) and edge component materials of interest, the s.e.e coefficient δ lies between 0.5 and 2. In the context of a particle-in-cell (PIC or QPIC) code this is to be interpreted as follows. Suppose that during a particular time step ∆t, which in a QPIC code is limited in magnitude by the Courant condition ∆t∆z/vethermal, there are N electrons incident upon a material target. Clearly, we would like the number N to be very large with a well-defined energy distribution so that we could intepret δ ∆N as the number of secondary electrons generated per time step by the ∆N electrons falling into the energy range (E,E+∆E). However, ∆t is typically small and so is the corresponding number N of incident electrons per time step. Thus we cannot obtain significant ensemble statistics in a time step. On the other hand, in a typical QPIC simulation proceeding on the ion thermal transit time scale, the number of time steps is very large, typically around 105.

The probabilistic meaning of the δ-coefficient can be therefore retained by working with time – rather than ensemble – averages, adopting an electron-by-electron Monte-Carlo technique. This is done as follows. Take a particular electron which has reached the material target. (i.e. we assume that its energy exceeds in magnitude the self-consistent wall potential). This electron is characterized by its energy E and the corresponding δ in Eq. (1) or Fig. 1. We furthermore assign the electron a uniformly distributed random number RAND from the interval (0,1). For the sake of illustration, let us consider the following simple emission algorithm: if the s.e.e probability is high enough, i.e. if δ1 or RANDδ1, then a secondary electron is emitted and the wall global charge is unchanged. If, however, δRAND, then secondary emission does not occur and the wall charge state changes as the incident electron neutralizes one of the incident ions. In the long run, for N primary electrons in the energy range (E,E+∆E), this algorithm clearly generates Nδ(E) secondary electrons.

We limit ourselves here to intermediate values of δ such that no more than one secondary electron can be emitted. This has far-reaching consequences for the QPIC code and simulation, which thereby remains relatively simple to modify and execute. The essential simplification, as was already mentioned, is that a single electron emission event does not change the wall charge state. Hence the ion wall charge throughout the simulation only depends on ions reaching the walls. The only complication which can arise (and ocassionally it does) is that with many secondary events we might in a particular time step run out of electrons capable of maintaining global wall charge balance since the corresponding primary electrons do not contribute to neutralizing the wall ion charge. Now imagine the much more severe consequences of a multiple secondary event. First, with multiple secondary events the wall positive charge increases and many more electrons than in single events are required for wall neutralization.In addition, the positive electron „holes“ arising in multiple seconary events must be included in global wall charge balance.

In the present work we therefore only deal with single secondary events and discontinue secondary electron production when in a time step the number of secondary events plus wall ion charge exceeds the number of electrons incoming to the target..

We now turn to the question of global charge conservation and the determination of selfconsistent wall potentials in the presence of single secondary electrons ( i.e. the global positive wall charge per time step q comes only from incident ions). Even without secondaries this is a complicated non-linear problem involving global charge conservation together with Kirchoff's second law for the tokamak circuit with the discharge. We thus need do solve the system of two equations

where the elements of the potential arrays VL,R(1:q) are the kinetic energies of the highest-energy incoming electrons. A solution of these two equations exists if the maximum V on one side is not smaller than the minimum V on the other, and at best it is only approximate, since the charges q are integers and the potentials Vl,R take on discrete values. The solutions are the wall potentials VwL,R, functions of the appropriate wall electron charges qeL,R. In practice we arrange the first q most energetic electrons on each side in descending order of kinetic energy, so that the first elements correspond to the electrons least likely to be reflected from, i.e. most likely to be "passing", the potential barrier. On each side we need q such electrons since we do not know beforehand which will be the combination of wall electron charges that will satisfy (2) and (3).

Let us now consider what happens when a secondary electron is emitted. As already mentioned, an electron which gives rise to a single secondary does not contribute to neutralizing the wall ion charge. Therefore another electron, less energetic than the ones already in the list, must be admitted. Thus for every secondary event, less energetic electrons take the place of the originally most energetic ones travelling to the target. It follows that the solutions of Eq. (3), the wall potentials VwL,R, decrease in magnitude in the presence of secondary electron emission. This tendency is naturally stronger with multiple secondary events.

The solution itself of Eqs (2,3) is quite complicated. Without secondaries, the algorithm is fairly straightforward, but with secondaries an additional major complication arises owing to the dependence of energy E in Eq. (1) on the unknown value Vw of the wall potential. With secondaries, the resulting non-linear trancendental equations are resolved by a bisection method, which will be presented elsewhere. Convergence of the method is fast, typically in less than 100 iterations.

Results for thermal electrons are obtained mostly without problems, but we have been experiencing difficulties resulting from too many secondaries associated with fast primaries generated when lower-hybrid power is switched on. As mentioned above, we dealt with this problem by disabling secondary electron production when a certain critical number of secondaries is exceeded. Selected results are presented for illustration below.

In conclusion, we introduced secondary electron production in the QPIC code, carried out necessary code modifications and debugging, and obtained some first exploratory results. The main difficulties with introducing secondary electrons have been identified, some of which (e.g.

Part II - PHYSICS multiple secondary events) will be dealt with in the continuation of this project. We consider the milestones set out for this part (B) of the mission to have been met.

Selected results Table II Some results of selected 6 QPIC simulations at the following conditions: "natural" Tore Supra SOL, Te0=Ti0=50 eV, n=5x1017 1/m3 ; initial number of particles per cell = 5000 ions, 5000 electrons; zero background drift velocity; Emax=300 eV. Further, E0 is the LH grill electric field strength, and we imposed a symmetric LH grill spectrum in order to test the expected symmetry of the simulation. The subscript "L" refers to the left wall, "R" to the right wall.

The wall potentials VwL,R given in Table II are time-averages taken over the period when the

**simulation has reached equilibrium, as indicated below in Fig. 2 :**

Fig 2. a) The number of ions and electrons in the simulation region as a function of time. The simulation is seen to reach equilibrium in about 3 ion transit times. b) The evolution of wall potentials VwL,R at the left and right targets, respectively. Results are obtained at simulation conditions of Table II.

ANNUAL REPORT 2007 ASSOCIATION EURATOM/IPP.CR We note a small effect of secondary electrons on Vwall as dmax increases in the three simulations at thermal conditions of Table II, but a relatively much larger effect at LH conditions. This indicates that the effect of secondary electrons should be taken into account in the intrepretation of probe diagnostics.

References [1] E. W. Thomas, Supp. Nucl. Fusion 1 (1991) 79.

[2] R. Kollath, Handbuch der Physik, Vol. 21, Springer Verlag, Berlin (1956) 232.

[3] S. Thomas and E. B. Pattison, J. Phys. D3 (1971) 349.

**In collaboration with:**

R. Pánek, Association EURATOM-IPP.CR The Faculty of Mathematics and Physics continued the work on development of hybrid code. The program under development presents a self-consistent model of central and border regions of this tokamak for different experimental conditions. The goal is to determine the anticipated range of plasma parameters in this system and to forecast the conditions for superior plasma confinement. Confinement like that is characterized by a low-temperature and low power load of the divertor. In addition the doctoral student Josef Havlíček participated in improving the EFIT 2006 program. This program enables calculating the current density, plasma shape and pressure and further plasma parameters. Our diploma, now PhD, student Michael Komm improved the SPICE code in such a manner that the program runs approximately 20 times faster.

The work in 2007 covered two main areas: (i) Numerical investigations of plasma parameters in COMPASS tokamak and (ii) Preparation of 3D particle code of magnetized plasmas. In frame of the first part numerical investigation of plasma parameters in COMPASS tokamak was performed and presented in the 34th EPS Plasma Conference on Plasma Physics in Warsaw, Poland [1]. The plasma parameters in the device were analyzed in the frame of the self-consistent model of central plasma and edge region. The possibility of achieving high recycling and detached regimes in COMPASS divertor was discussed (see Fig. 1).

** Figure 1. Total energy losses Ploss, energy confinement time τE, average plasma temperature Tcore, plasma temperature at the LCMS Tes, plasma temperature at the plate Tep and normalized parameter βn as functions of normalized volume average density ne/nG for different configurations.**

Scaling (B).

In frame of the task preparation of 3D particle code of magnetized plasmas there was developed during 2007 the self-consistent fully 3D Particle-In Cell code for modeling of plasma-solid interaction with new Poisson solver based on direct LU decomposition combined with multigrid approach. The model was tested on a cylindrical cavity geometry, i.e. on a hollow cylindrical ANNUAL REPORT 2007 ASSOCIATION EURATOM/IPP.CR chamber opened to the plasma, in both collisionless and collisional plasmas (see Fig. 2). The basic studied feature was the influence of non-axial orientation of magnetic field – this parameter is important for the analysis of experimental data from the Katsumata probe. Earlier version of the model was published in [2], the improved model developed in 2007 in [3,4].

** Figure 2. The angular (left) and axial (right) distributions of ions on the inner cylindrical wall, with the magnetic field inclination α and ambient gas pressure as parameters.**

**References:**

[1] E. Havlíčková, R. Zagórski, R. Pánek, Numerical Investigations of Plasma Parameters in COMPASS Tokamak, 34 EPS Conference on Plasma Physics, 2-6 July 2007, Warsaw, Poland, proceedings reference ECA Vol. 31F, P4-117 (2007).

[2] Z. Pekárek, R. Hrach, Influence of Non-Axial Magnetic Field in a 3D3V Paraticle-In-Cell Plasma Model, Vacuum 82 (2008), 244-247.

[3] Z. Pekárek, Š. Roučka, R. Hrach, 3D particle simulations of plasma-solid interaction:

Magnetized plasma and a cylindrical cavity, 17th Interational Vacuum Congress IVC-17, 2-6 July 2007, Stockholm, Sweden.

[4] Z. Pekárek, R. Hrach: Multigrid methods as a Poisson equation solver for 3D PIC models of plasma-solid interaction, Conference on Computational Physics CCP2007, 5-8 September 2007, Brussels, Belgium.

This contribution presents the Particle-in-Cell (PIC) modeling of the Ball-pen probe using the 2D XOOPIC code developed in Berkeley University. The ball-pen probe has been designed for the direct measurements of the plasma potential. The probe head consists of a movable conducting collector, which is partially shielded by the isolator. This construction allows modification of the ratio of the electron and ion saturation currents with respect to the collector position.

The basic idea of the direct plasma potential measurement by the ball-pen probe [1,2] is to adjust − + the ratio R = I sat I sat in Eq. (1) to be equal to one. When this is achieved, the floating potential of the probe Vfl is equal to the plasma potential Φ as follows from V fl = Φ − T e ln( R ) (1) Figure 1 shows the results of the PIC simulations (Te=Ti=20 eV, Φ=0, Bt=1T) of I-V characteristics of the ball-pen probe using 2D XOOPIC (developed in Berkeley University) [3].