2D ArcPIC Code Description: Description of Methods and User ...

4 downloads 9309 Views 2MB Size Report
Sep 1, 2014 - 2D ArcPIC Code Description. Description of Methods and. User / Developer Manual. (second edition). 2D ArcPIC v 3.00. September 2014.
CERN – EUROPEAN ORGANIZATION FOR NUCLEAR RESEARCH

CLIC – Note – 1032

2D ARCPIC CODE DESCRIPTION: DESCRIPTION OF METHODS AND USER / DEVELOPER MANUAL (SECOND EDITION) Kyrre Sjobak, CERN, Geneva, Switzerland and University of Oslo, Norway Helga Timkó, CERN, Geneva, Switzerland and Helsinki Institute of Physics, Finland

01/09/2014

CERN-OPEN-2014-048

Abstract Vacuum discharges are one of the main limiting factors for future linear collider designs such as that of the Compact LInear Collider (CLIC). To optimize machine efficiency, maintaining the highest feasible accelerating gradient below a certain breakdown rate is desirable; understanding breakdowns can therefore help us to achieve this goal. As a part of ongoing theoretical research on vacuum discharges at the Helsinki Institute of Physics, the build-up of plasma can be investigated through the particle-in-cell method. For this purpose, we have developed the 2D ArcPIC code introduced here. We present an exhaustive description of the 2D ArcPIC code in several parts. In the first chapter, we introduce the particle-in-cell method in general and detail the techniques used in the code. In the second chapter, we describe the code and provide a documentation and derivation of the key equations occurring in it. In the third chapter, we describe utilities for running the code and analyzing the results. The last chapter contains suggestions for viable and useful avenues for further work on the code. The code and documentation are original work of the authors, written in 2010–2014, and is therefore under the copyright of the authors.

Geneva, Switzerland September 2014

2D ArcPIC Code Description Description of Methods and User / Developer Manual (second edition)

2D ArcPIC v 3.00 September 2014

Kyrre Sjobak1,3 , MSc and Helga Timkó1,2 , PhD

1 CERN, 2 Helsinki

Institute of Physics and 3 University of Oslo

Release Notes ArcPIC v2.07 (2011-03) Tested, debugged version. This version is described in CLIC Note 872 [1]. ArcPIC v2.08 (2012-07) Magnetic field changed from (Bz , Br , 0) to (Bz , 0, Bt ). Ion masses slightly corrected. ArcPIC v3.00 (2014-01) Added a modular system for particle boundary conditions, circuit model, and initial particle distributions. Added new particle boundary conditions and circuit models. Re-wrote the particle sorting algorithm, greatly reducing the memory footprint. Improved the Fowler-Nordheim model for electron emission from flat cathode surface, and also the field solver boundary conditions. Updated the documentation (this document).

iii

iv

Abstract Vacuum discharges are one of the main limiting factors for future linear collider designs such as that of the Compact LInear Collider (CLIC). To optimize machine efficiency, maintaining the highest feasible accelerating gradient below a certain breakdown rate is desirable; understanding breakdowns can therefore help us to achieve this goal. As a part of ongoing theoretical research on vacuum discharges at the Helsinki Institute of Physics, the build-up of plasma can be investigated through the particle-in-cell method. For this purpose, we have developed the 2D ArcPIC code introduced here. We present an exhaustive description of the 2D ArcPIC code in several parts. In the first chapter, we introduce the particle-in-cell method in general and detail the techniques used in the code. In the second chapter, we describe the code and provide a documentation and derivation of the key equations occurring in it. In the third chapter, we describe utilities for running the code and analyzing the results. The last chapter contains suggestions for viable and useful avenues for further work on the code. The code and documentation are original work of the authors, written in 2010–2014, and is therefore under the copyright of the authors.

Acknowledgments The development of the code has been carried out in frame of a collaboration between the Helsinki Institute of Physics and CERN. The professional guidance of Prof. R. Schneider and PhD K. Matyash from the Max-Planck-Institut für Plasmaphysik is gratefully acknowledged. We thank Lotta Mether (Helsinki Institute of Physics and CERN) for many insights into the properties of the models and the simulations, for playing an important role in debugging the 2D ArcPIC code, and for helping to write this manual. The research has received funding from the European Commission under the FP7 Research Infrastructures project EuCARD, grant agreement no. 227579. For further details, see http://cern.ch/eucard. We acknowledge CSC – IT Center for Science Ltd. in Espoo, Finland for the allocation of computational resources

Contents Release Notes

iii

Abstract

iv

Acknowledgments

iv

Contents

v

1

2

Methods 1.1 The particle-in-cell method . . . . . . . . 1.1.1 Superparticles . . . . . . . . . . 1.1.2 Field solver . . . . . . . . . . . . 1.1.3 Particle mover . . . . . . . . . . 1.1.4 Particle and field weighting . . . 1.1.5 Collisions . . . . . . . . . . . . . 1.1.6 Accuracy and stability conditions 1.2 The 2D ArcPIC model . . . . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

1 1 1 2 3 3 4 4 5

Code documentation 2.1 Structure of the code . . . . . . . . . . . . . . . . . . . . . 2.2 Scaling in 2D PIC . . . . . . . . . . . . . . . . . . . . . . . 2.2.1 The first equation of motion – particle displacement . 2.2.2 The second equation of motion – acceleration . . . . 2.2.3 The electric field . . . . . . . . . . . . . . . . . . . 2.2.4 Poisson’s equation . . . . . . . . . . . . . . . . . . 2.3 Electric field calculation . . . . . . . . . . . . . . . . . . . 2.4 Potential and electric field boundary conditions . . . . . . . 2.5 Interpolation scheme . . . . . . . . . . . . . . . . . . . . . 2.6 Particle pusher equations . . . . . . . . . . . . . . . . . . . 2.6.1 Electrostatic case . . . . . . . . . . . . . . . . . . . 2.6.2 Magnetic case . . . . . . . . . . . . . . . . . . . . . 2.7 Collision routines . . . . . . . . . . . . . . . . . . . . . . . 2.7.1 Coulomb collisions for e- –e- and Cu+ –Cu+ pairs . . 2.7.2 Other collisions . . . . . . . . . . . . . . . . . . . . 2.7.3 Parallelization using OpenMP . . . . . . . . . . . . 2.8 Circuit models . . . . . . . . . . . . . . . . . . . . . . . . . 2.8.1 Shockley-Ramo gap current calculation . . . . . . . 2.8.2 Rescaling to dimensionless units . . . . . . . . . . . 2.8.3 The output file circuit.dat . . . . . . . . . . . . 2.8.4 Models . . . . . . . . . . . . . . . . . . . . . . . . FixedVoltage_resistorCapacitor . . . . . . . . . . . FixedVoltage_resistor . . . . . . . . . . . . . . . . . TwoCapacitors . . . . . . . . . . . . . . . . . . . . 2.9 Surface models for particle boundary conditions . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . .

7 7 8 8 9 10 11 14 14 18 19 19 19 22 22 23 23 24 24 25 26 26 27 27 28 28

v

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

CONTENTS

vi

2.9.1 2.9.2 2.9.3 2.9.4

2.10

2.11 2.12 2.13

2.14 3

4

Output files . . . . . . . . . . . . . . . . . . . . . . . . Injection of particles . . . . . . . . . . . . . . . . . . . Fowler-Nordheim field emission . . . . . . . . . . . . . Models . . . . . . . . . . . . . . . . . . . . . . . . . . ArcSimple . . . . . . . . . . . . . . . . . . . . . . . . ArcBoundsOriginal . . . . . . . . . . . . . . . . . . . . ArcBoundsOriginalNewHS . . . . . . . . . . . . . . . FlexFN models . . . . . . . . . . . . . . . . . . . . . . Initial particle distribution . . . . . . . . . . . . . . . . . . . . 2.10.1 Generating a uniform density in 2D PIC . . . . . . . . . 2.10.2 The UniformRestricted initial particle distribution . . Density calculation . . . . . . . . . . . . . . . . . . . . . . . . Monitoring energy conservation . . . . . . . . . . . . . . . . . Output files . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.13.1 Output from main() . . . . . . . . . . . . . . . . . . . 2.13.2 Outputting using outputz.cpp, moms.cpp and vdf.cpp Control routines, stability of PIC . . . . . . . . . . . . . . . . .

Running 2D ArcPIC, utilities and analysis system 3.1 Supporting libraries . . . . . . . . . . . . . . 3.2 makeRundir.py . . . . . . . . . . . . . . . . 3.3 calcScaling.py . . . . . . . . . . . . . . . . . 3.4 Makefile . . . . . . . . . . . . . . . . . . . . 3.5 Analyses . . . . . . . . . . . . . . . . . . . . Recommended directions for future development

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . . . . . . . . . . . . . .

. . . . .

. . . . . . . . . . . . . . . . .

. . . . .

. . . . . . . . . . . . . . . . .

. . . . .

. . . . . . . . . . . . . . . . .

. . . . .

. . . . . . . . . . . . . . . . .

. . . . .

. . . . . . . . . . . . . . . . .

. . . . .

. . . . . . . . . . . . . . . . .

. . . . .

. . . . . . . . . . . . . . . . .

29 29 30 30 31 31 33 33 33 34 34 35 35 37 37 38 41

. . . . .

43 44 44 44 44 45 47

A List of files in the source code

49

B Example input.txt

53

C Miscellaneous useful formulas

57

Bibliography

59

1 Methods 1.1

The particle-in-cell method

As an introduction to the code, we shall first briefly review the particle-in-cell (PIC) method, a method that has been developed in the 1950’s by Buneman [2], Dawson [3], and others. A more thorough description of it can be found in the classic textbooks of Hockney & Eastwood [4] or of Birdsall & Langdon [5]. The PIC method is a technique widely applied to all sorts of plasma simulations, since it combines the description of individual “components” with the solution of partial differential equations (PDE). In plasma physics, these individual components translate to particles or fluid elements, while the set of PDEs translates to Maxwell’s equations; and thus the PIC method allows for a kinetic description of a plasma system. Why the method is called particle-in-cell will become clear somewhat later in conjunction with the weighting scheme. In PIC, particles and plasma macro-quantities, such as density, current density, or potential are treated in a different way. Particles are followed in a continuous phase space, in a Lagrangian (moving) frame, while macro-quantities are discretized onto mesh points and described in a Eulerian (stationary) frame of reference. The temporal description is the same for both particles and macro-quantities, discretised into time steps; however time steps for different particle species may vary. In that case, the fastest species and the macro-quantities should have the same temporal resolution. The whole simulation procedure is summarised schematically in Fig. 1.1. First, the simulation is set up by choosing simulation parameters (grid spacing, time step, superparticle weighting, etc.) adequate for the physical problem. By “adequate” we mean that parameters should be “well-guessed” in order to ensure the stability and accuracy of the code; we will come back to the criteria at the end of this section. Once the simulation is set up, the problem solving consists of a loop over following steps: 1. Solving Maxwell’s equations. 2. Obtaining fields and forces acting on the particles. 3. Moving particles by solving the equations of motion (EOM). 4. Applying collisions, boundary conditions for particles (injection, absorption, etc.). 5. Updating macro-quantities that are sources to Maxwell’s equations. We shall not describe all of these points from a general point of view. Instead, below we will focus on some important details of a PIC simulation. Methods specific to the 2D ArcPIC code will be detailed in Sec. 1.2.

1.1.1

Superparticles

Since usually the number of particles that would have to be simulated to model a real system exceeds our computational capacities, the PIC method makes use of so-called superparticles. These superparticles represent many real particles at the same time, that is, the real number of particles in the system is “rescaled” to a smaller number of computational particles in the

1

CHAPTER 1. METHODS

2

Figure 1.1: Schematic of the particle-in-cell simulation loop.

model. We may do so, because the particles move according to the Lorentz force, which depends only on the charge-to-mass ratio, and the charge-to-mass ratio is the same for both real and superparticles.

1.1.2

Field solver

Maxwell’s equations are solved in the field solver part of PIC. The most common approaches for solving PDEs use one of the following three methods: Finite difference method (FDM) replaces the continuous domain with a discrete grid of points; electromagnetic fields are calculated then on this grid. Derivatives are approximated with differences between neighbouring grid-point values and thus the FDM turns PDEs into algebraic equations. Finite element method (FEM) divides the continuous domain into a discrete mesh of elements. Each element contains a set of basis functions, which together form a basis for the whole function space in which the numerical solution to the PDEs exist. The Galerkin method is then applied to the PDEs, leading to a finite set of algebraic equations which can be solved numerically, typically using a library for linear algebra. Spectral methods such as the fast Fourier transform. Also in this case, the solution of the PDEs are written as a linear combination of basis functions, however, this time the basis functions are high order and defined globally over the whole domain. The domain itself is therefore not discretised, but remains continuous.

1.1. THE PARTICLE-IN-CELL METHOD

1.1.3

3

Particle mover

The equations of motion (EOM) for the particles are solved in the so-called particle mover or pusher, which is often the most time consuming part of PIC. It is therefore required to be of high accuracy and speed. The schemes used for the particle mover can be divided into implicit and explicit solvers. While implicit solvers calculate the particle velocity from the already updated fields, explicit solvers use only the old force from the previous time step, and are therefore simpler and faster, but require a smaller time step. Two frequently used schemes are the leapfrog method [5], which is an explicit solver, and the Boris scheme [6], which is implicit. For plasma applications, the leapfrog method takes the following form: xk+1 − xk = vk+1/2 , ∆t   vk+1/2 − vk−1/2 vk+1/2 + vk−1/2 q = × Bk , Ek + ∆t m 2

(1.1) (1.2)

where the subscript k refers to “old” quantities from the previous time step, k + 1 to updated quantities from the next time step (i.e. tk+1 = tk +∆t), and velocities are calculated in-between the usual time steps tk . The bolding of quantities in the equation indicate that these are 3vectors. This notation will be used throughout the document. In comparison, the equations of the Boris scheme are: xk+1 = xk + ∆tvk+1/2 , 0

0

vk+1/2 = u + q Ek ,

(1.3) (1.4)

where u0 = u + (u + (u × h)) × s , 0

(1.5)

u = vk−1/2 + q Ek ,

(1.6)

h = q 0 Bk ,

(1.7) 2

s = 2h/(1 + h ) , ∆t q and q 0 = · . 2 m

(1.8) (1.9)

We see that if B = 0, the leapfrog- and Boris scheme reduce to the same set of equations.

1.1.4

Particle and field weighting

The origin of the name “particle-in-cell” comes from the way plasma macro-quantities are assigned to simulation particles, that is, the particle weighting. Particles can be situated anywhere in the continuous domain, but macro-quantities (including the fields) are calculated only on the mesh points. To obtain the macro-quantities, one assumes that the particles have a given “shape” determined by the shape function S(x − X) , (1.10) where x is the coordinate of the particle and X the observation point (usually the grid point xi ). Perhaps the easiest and most used choice for the shape function is the so-called cloud-incell (CIC) scheme, which is a first order (linear) weighting scheme. In any case, whatever the

CHAPTER 1. METHODS

4

scheme is, the shape function has to satisfy the following conditions [7]: Space isotropy, charge conservation, and increasing accuracy (convergence) for higher-order terms. The fields obtained from the field solver are determined only on the grid points and can’t be used directly in the particle mover to calculate the force acting on particles. They thus have to be interpolated via the field weighting X E(x) = Ei S(xi − x) , (1.11) i

where the subscript i labels the grid point. To ensure that the forces acting on particles are selfconsistently obtained, the way of calculating macro-quantities from particle positions on the grid points and interpolating fields from grid points to particle positions has to be consistent, since they both appear in Maxwell’s equations. Above all, the field interpolation scheme should conserve momentum. This can be achieved by choosing the same weighting scheme for particles and fields, and by ensuring the appropriate space symmetry (i.e. no self-force and fulfilling the action-reaction law) of the field solver at the same time [7].

1.1.5

Collisions

As the field solver is required to be free of self-forces, inside a cell the field generated by a particle must decrease with decreasing distance from the particle, and hence inter-particle forces inside the cells are underestimated. This can be balanced with the aid of Coulomb collisions between charged particles. Simulating the interaction for every pair in a large system would be too computationally expensive, so several Monte Carlo methods have been developed instead. A widely used method is the binary collision model [8], in which particles are grouped according to their cell, then these particles are paired randomly, and finally the pairs are collided. In a real plasma, many other reactions may play a role, ranging from elastic collisions, such as collisions between charged and neutral particles, over inelastic collisions, such as electronneutral ionization collisions, to chemical reactions; each of them requiring separate treatment. Most of the collision models handling charged-neutral collisions use either the direct MonteCarlo scheme, in which all particles carry information about their collision probability, or the null-collision scheme [9, 10], which does not analyse all particles but uses the maximum collision probability for each charged species instead.

1.1.6

Accuracy and stability conditions

As in every simulation method, also in PIC, the time step and the grid spacing must be well chosen, so that the shortest time- and length scale phenomena are properly resolved in the problem. In addition, time step and grid spacing have an impact on the speed and accuracy of the code. As a rule of thumb, the following important conditions regarding the grid spacing ∆x and the time step ∆t should be fulfilled in order to ensure the stability of the solution: ∆x < 3.4λDb

(1.12)

−1 0.2ωpe

(1.13)

∆t ≤

These can be derived considering the harmonic oscillations of a one-dimensional unmagnetised plasma [11]. Perhaps not surprisingly, the natural length scale in the plasma is given by the −1 (Eq. C.2). Debye length λDb (Eq. C.1), and time scale by the inverse plasma frequency ωpe

1.2. THE 2D ARCPIC MODEL

5

Figure 1.2: Illustration of the (r, z)-geometry of the code, which can be interpreted as an infinitesimal slice of a three-dimensional cylindrically symmetric system.

1.2

The 2D ArcPIC model

The 2D ArcPIC code can be classified as a 2D3V PIC-MCC code. The “MCC” stands for the code being equipped with Monte-Carlo collision routines. The 2D3V means that particles are described in a five-dimensional phase space with two coordinate and three velocity components. The geometry of the coordinate space is (r, z), assuming cylindrical symmetry of the simulated system. This two-dimensional (r, z)-space can be pictured as an “infinitesimal” slice of a full three-dimensional cylindrically symmetric system (see Fig. 1.2). Indeed, in the pusher of the code (Sec. 2.6), particles are moved according to their three velocity components in three dimensions, and are then rotated back into the plane of reference (i.e. the simulation space). A uniform grid is used throughout all the simulation domain, with equal grid spacing in both coordinates ∆z = ∆r. Phenomena are described electrostatically only, that is, instead of solving the full set of Maxwell’s equations, only Poisson’s equation is taken into account (Sec. 2.2). The field solver transforms the Poisson equation into a matrix equation using the FDM, and uses the external package SuperLU [12] for inverting the matrix with the so-called “lower-upper factorization method”. The particle pusher (Sec. 2.6) uses the Boris method. Although the field solver is purely electrostatic, the code supports an optional external magnetic field. Also in this case, the EOMs are solved with the Boris method. For inter- or extrapolation between particle positions and grid points, the CIC weighting scheme is used (Secs. 2.5 & 2.11). Some special features of the code related to vacuum discharges are certain boundary conditions. This includes circuit models which matches the DC spark experiment [13] (Sec. 2.8), several surface models (Sec. 2.9) including sputtering and Fowler-Nordheim field emission of electrons (Sec. 2.9.3), and collisions (Sec. 2.7) which play a key role in the process. For optimization reasons, all the physical variables appearing in the code are rescaled to dimensionless units, so that unnecessary multiplications with constants are avoided. In Chapter 2 we will introduce the method of rescaling and derive all key equations in the form they appear in the code. However, may it be pointed out here that since only dimensionless quantities are used in the code, results could in fact apply to “any” time- and length scale. Only in two places in the code, namely in the collision (Sec. 2.7) and particle boundary condition (Sec. 2.9) routines, are we bound to certain dimensional values, which determines what is the “physical” scale of

6

CHAPTER 1. METHODS

the end result. The code is written in C++ language, and most of the supporting tools in Python. The development and running environment is provided by the Linux operating system and standard tools. The use of C++ allows for modularization using object-oriented techniques, creating a clear interface between the general parts of the code dealing with particle pushing, field calculation and collisions, and the parts implementing circuit models (Sec. 2.8), surface models (Sec. 2.9), and initial particle distributions (Sec. 2.10). In the future (Chapter 4), this modularization may be extended to other parts of the program.

2 Code documentation This chapter documents the code compiled into the 2D ArcPIC executable. It first describes the general structure of the code and the rescaling to dimensionless units. After this comes a description of the general parts of the code, dealing with particle pushing, field calculation, and collisions. Then comes a description of the modular parts of the code, which are the surface models (particle injection / removal), circuit models, and initial particle distributions. Finally, the output and monitoring subroutines are described.

2.1

Structure of the code

The 2D ArcPIC program is built for simulation of copper vacuum arcs in a 2D3V phase space, and uses three types of particles: Electrons (e- ), singly charged copper ions (Cu+ ), and neutral copper atoms (Cu). The geometry of the simulated domain and the coordinate system is illustrated in Figure 1.2. Here a flat cathode is assumed to be present at z = 0, and a flat anode at z = Zmax . A voltage is applied to these electrodes by a circuit model, and particles are injected/removed from the volume at its boundaries by a surface model. As described in Chapter 1, the PIC method is structured around a main time-stepping loop, illustrated in Figure 1.1. In 2D ArcPIC, each of the steps shown in Figure 1.1 are broken down into several sub-steps, resulting in the layout illustrated in Figure 2.1. In this loop, the position and velocity of the electrons are updated every time step of the simulation, as is the electric potential and field. The particle pusher is called at longer time intervals (ion time steps) for the slower ions and neutrals. The particle pusher is described in Sec. 2.2.1, 2.2.2 and 2.6, while

not output step

Timestep

Output and check of stability criterion Move ions & neutrals yes

Move e− Remove e−

Inject ions Ion step?

no

Timestep circuit & particle bound. cond.

Remove ions & neutrals

Inject e− not neutral inj. step

not e− inj. step

Update potential

Inject neutrals

Other collisions

not other collision step

finished

Input & initialize

not ion inj. step not e− collision step

start

e− +e− elastic collisions

Order particle arrays

Figure 2.1: The main sequence of the 2D ArcPIC code. The green rectangles indicates modular pieces, i.e. the currently loaded implementation of Circuit and ArcBounds.

7

CHAPTER 2. CODE DOCUMENTATION

8

the field solver is described in Sec. 2.2.3, 2.2.4, 2.3, and 2.4. After the particle pusher has moved the particles, their positions are checked to ensure that they are inside of the simulation domain. The particles which are not inside the simulation domain are then removed by the surface models. These models may use the position and energy of the removed particles as input to the algorithms handling the injection of new particles (sputtering etc.), as described in Sec. 2.9. When the particle positions have been resolved by the pusher and surface models, the particles are collided using Monte-Carlo collision techniques. This updates the velocity vectors of the particles. The collision algorithms are applied separately to the particles inside of each simulation grid cell, requiring the particle arrays to be sorted according to which grid cell the particle is in. Finally, after the field solver has ran, output of quantities for later analysis and check of stability criteria is performed, before the loop repeats.

2.2

Scaling in 2D PIC

To build-up and derive a fully dimensionless set of variables, we only need to rescale hierarchically the main equations that determine the evolution of the system: The first and second equation of motion together with the electric field calculation and the Poisson equation. Here is how we shall proceed. In each case we start from the dimensional equation (in SI-units) and multiply then by a rescaling factor. This rescaling factor allows us to define dimensionless quantities (denoted by a tilde) which are then used in the code and enable us to establish the connection between dimensional and dimensionless variables.

2.2.1

The first equation of motion – particle displacement

The first equation of motion determines the displacement of electrons (or ions, i.e. non-electron species) during one electron (ion) time step. This may simply be written as dz = v dt, where z denotes the particle’s coordinate, v its velocity and t time. In the Boris scheme, this equation takes the form1   1 (ion) (ion) dz = zi+1 − zi−1 = vz, i+ 1 + vz, i− 1 ∆t(ion) . × (2.1) ∆z 2 2 Defining dimensionless coordinates and velocities using the dimensional grid spacing and time step as rescaling factors yields z ∆z ∆t v˜z ≡vz ∆z ∆tion v˜z, ion ≡vz, ion . ∆z z˜ ≡

(2.2) (2.3) (2.4)

Further, the dimensionless time step and grid spacing2 can be related to the Debye length λDb and plasma frequency ωpe at reference density nref and temperature Tref as  f ≡ωpe ∆t = 0.2 = ω ∆t ˜ pe (2.5) 1 As a shorthand notation, the operations and factors that rescale dimensional equations to dimensionless ones are written down on the right hand side of the equations, marked by a slash. 2 This is based on the stability requirements in Sec. 1.1.6, but is stricter for ∆x.

2.2. SCALING IN 2D PIC

9

f ≡ ∆z = 0.5 ∆z λDb

  1 = . ˜ Db λ

(2.6)

Using this, Eq. (2.1) can be rescaled as: d˜ z = z˜i+1 − z˜i−1 =

∆t (v ˜z, i+ 1 + v˜z, i− 1 . 1 + v z, i− 12 ) = v 2 2 ∆z z, i+ 2

(2.7)

The above equations apply also for r and vr in a similar fashion; the code uses a homogeneous grid with ∆z = ∆r. Dimensional variables can be converted back from dimensionless variables as follows: zREAL f z˜ = ∆z λDb f vREAL ∆z vREAL = = ⇒ v˜ f λDb ωpe vt, e ∆t f v ion v˜ ∆z ⇒ REAL = f dt_ion vt, e ∆t ⇒

zREAL = ∆z z˜ vREAL =

∆z v˜ ∆t

ion vREAL =

∆z v˜ ∆tion

(2.8) (2.9) (2.10)

where dt_ion = ∆tion /∆t is an integer3 . In the special case of the thermal velocity (Eq. C.5) we obtain f f 1 ∆t ∆t ∆t . (2.11) = vt, e = v˜te = vt, e f ωpe λDb f ∆z ∆z ∆z | {z } vt, e

2.2.2

The second equation of motion – acceleration

The second equation of motion describes particle acceleration. Firstly, for electrons the dimensional equation (z-component) is: ∆t e vz, i+ 1 − vz, i− 1 = − Ez, i ∆t , × (2.12) 2 2 ∆z me with e being the elementary charge, me the electron mass, and E the electric field. Define the dimensionless electric field as ˜z or r ≡ E

e ∆t2 Ez or r 2me ∆z

(2.13)

to obtain the dimensionless equation  v˜z, i+ 1 − v˜z, i− 1 = −2 2

2

e ∆t2 2me ∆z



˜z, i . Ez, i = −2E

(2.14)

Note the convention of dividing the electric field contribution into two parts, hence the factor of two, which is convenient in the Boris scheme. So, in the pusher we update velocities as (when there is no magnetic field): ˜. ˜=v ˜ old − 2E v (2.15) 3 Throughout the whole document, a “variable” typeset as variable always denotes a certain variable appearing in the code.

CHAPTER 2. CODE DOCUMENTATION

10

For this to remain true also for ions, the electric field is rescaled for ions with the ion time step dt_ion = ∆tion /∆t and the ion mass Mi in e_ion.cpp to an ‘ion electric field’   dt_ion P ˜ Ei  dt_ion  me X ˜  2 me  i=1 ˜ Eion ≡ − dt_ion Ei = −dt_ion   Mi Mi  dt_ion  (2.16) i=1 = − dt_ion2

me ˜ hEi idt_ion . Mi

The summation is due to the fact that for ions the pusher is called only every ion step ∆tion instead of every time step ∆t. The ions velocity has to be rescaled with ∆tion as well, ∆t e ion ion ion vz, i+ 1 − vz, i− 1 = Ez, i ∆tion (2.17) × 2 2 Mi ∆z to obtain the same form of dimensionless equation as for electrons    e ∆t2 ion ion 2 me v˜z,i+ 1 − v˜z, i− 1 = −2 −dt_ion Ez, i . 2 2 Mi 2me ∆z

(2.18)

˜ion is first summed over dt_ion steps, To average out numerical fluctuations, E ˜ ion = E

dt_ion X

˜i E electron

(2.19)

i=1 e and then rescaled with −dt_ion m Mi to obtain

˜ ion = −dt_ion me E Mi

2.2.3

dt_ion X i=1

2 me ˜i hEion idt_ion . E electron = −dt_ion Mi

(2.20)

The electric field

~ takes the dimensional form Deriving the electric field from the potential ϕ: E = −∇ϕ ϕi−1 − ϕi+1 e ∆t2 Ez, i = . × (2.21) 2me ∆z 2∆z Defining then the dimensionless potential as ϕ˜ ≡

e ∆t2 ϕ me ∆z 2

results in the dimensionless equation   1 e ∆t2 1 ˜ Ez, i = (ϕi−1 − ϕi+1 ) = (ϕ˜i−1 − ϕ˜i+1 ) . 2 4 me ∆z 4

(2.22)

(2.23)

Dimensional variables are obtained by following rescaling factors: EREAL =

f 2me ∆z ˜ 2me ∆z 2 ˜ E = λDb ωpe E e ∆t2 e ∆t f2

(2.24)

2.2. SCALING IN 2D PIC

11

ϕREAL =

f2 me ∆z me ∆z 2 2 2 ϕ ˜ = ˜ 2 λDb ωpe ϕ e ∆t2 e ∆t f

(2.25)

The latter equation can be rewritten in terms of Te , since 2 λ2Db ωpe =

and therefore

2.2.4

Tref ε0 Tref e2 nref = = vt,2 e 2 e nref ε0 me me

(2.26)

f2 f 2  Tref  me ∆z ∆z 2 2 λDb ωpe = e ∆t e f2 f2 ∆t

(2.27)

f2 ϕREAL [V ] ∆z ⇒ = ϕ˜ . Tref [eV ] f2 ∆t

(2.28)

Poisson’s equation

Poisson’s equation is in its original form is ∇2 ϕ = − ερ0 = − εe0 (ne + ni ), with ε0 being the vacuum permittivity and using the convention ne 0

(2.30)

for the electron and ion number densities ne and ni , respectively. First, we shall derive the scaling that applies to ne . Written in cylindrical coordinates,  2  ∂ 1 1 ∂2 e + + ϕ = − ne . (2.31) ∂r2 r ∂r ∂z 2 ε0 With the finite differences scheme, for inner mesh points we get the so-called “five-point formula”: X ci ϕi (j, k) ≡aϕ(j − 1, k) + bϕ(j + 1, k) + cϕ(j, k − 1) + dϕ(j, k + 1) + eϕ(j, k) 5−p.

=−

e ne (j, k) ε0 (2.32)

where j ∈ [0, nr] and k ∈ [0, nz] enumerate mesh points, and where the factors ci are rescaled as follows:   1 1 2 a= − + rj+1 − rj−1 rj rj − rj−1   (2.33) 1 1 1 2 a ˜ = 2 − + ≡ ∆r r˜j+1 − r˜j−1 r˜j r˜j − r˜j−1 ∆z 2 So for a uniform grid r˜j = j and for inner mesh points we have:  f2  ∆z 1 a ˜=a ˜(j) = 1− 2j f2 ∆r

(2.34)

CHAPTER 2. CODE DOCUMENTATION

12

 f2  ˜b = ˜b(j) = ∆z 1 + 1 2j f2 ∆r c˜ = d˜ = 1

(2.35) (2.36)

2

e˜ = −2

f ∆z −2 f2 ∆r

(2.37)

f 2 = ∆z f 2 . With the aid of the above formulae, we can then rescale where now of course ∆r Poisson’s equation X e ∆t2 e (2.38) ci ϕi (j, k) = − ne . · ∆z 2 ε0 me ∆z 2 5−p.

Here we choose a scheme in which the volume of the cell is explicitly taken into account when density is calculated. That is, a particle in an inner cell will correspond to a higher density than a particle in an outer cell, since the outer cell has bigger volume. In the code, two routines are coupled for the calculation of potential and density: density_2D(...) and potential_backsolve_2D(...) or potential_backsolve_BC23(...). We are searching for a dimensionless charge q˜ such that the dimensionless density dens(j, k) at a grid point (j, k) would fulfil X

c˜i ϕ˜i (j, k) = −

5−p.

˜e, cell (j, k) N e2 ! ∆t2 ne (j, k) = − q˜ ≡ −g n_e(j, k) ≡ −dens(j, k) , ε0 m e V˜cell (j)

(2.39) ˜e, cell (j, k) is the number of particles “projected” onto the grid point (j, k) (see Eq. 2.41) where N and V˜cell (j) is the dimensionless volume of the cell (Eq. 2.46). Note the convention of the minus sign on the right-hand side. We can simplify this equation to ˜e, cell (j, k) e2 ne (j, k) ! N 2 ∆t2 ne (j, k) = ωpe ∆t2 = q˜ ε0 me | {z } nref V˜cell (j, k)

(2.40)

g2 ∆t 2

2 = e nref , with n since ωpe ref being a reference density chosen initially. In all these equations, ε0 me ˜ Ne, cell (j, k) is obtained as follows: It is the number of particles calculated on the grid point (j, k) by summing over the four neighbouring cell’s particle distribution f (rp ) and weighting it with the particle weighting function W (rp − rcell ), that is

X 5−p.

g2 ne (j, k) =! −dens(j, k) = −˜ q c˜i ϕ˜i (j, k) = −∆t nref

X  W (rp − rcell ) 4−c.

V˜cell (j)

 f (rp ) . (2.41)

In general, the volume of the cell calculated with the Verboncoeur method [14] on the grid point is independent of k and given as:  π   3 (r1 − r0 )(2r0 + r1 )∆z, (real) Vcell (j) = π3 [rj+1 (rj+1 + rj ) − rj−1 (rj + rj−1 )] ∆z,  π 3 (rnr − rnr−1 )(rnr−1 + 2rnr )∆z,

j=0 j ∈ (0, nr) j = nr

(2.42)

2.2. SCALING IN 2D PIC

13

For a homogeneous grid with rj = j∆r and ∆r = ∆z, this reduces to  π 2  j=0  3 1 · 1∆r ∆z, (real) π 2 Vcell (j) = 3 [(j + 1)(2j + 1) − (j − 1)(2j − 1))] ∆r ∆z, j ∈ (0, nr)  π 2 j = nr 3 1 · (3nr − 1)∆r ∆z,  π 2  j=0  3 ∆r ∆z, 2 = 2πj∆r ∆z, j ∈ (0, nr) .   1 2 π(nr − 3 )∆r ∆z, j = nr More generally, for rmin 6= 0 we would have   3rmin π 2  j=0  3 ∆r ∆z ∆r + 1 , (real) rmin 2 Vcell (j) = 2π∆r ∆z ∆r + j , j ∈ (0, nr) .    rmin 1 2 π∆r ∆z ∆r + (nr − 3 ) , j = nr

(2.43)

(2.44)

(2.45)

In its dimensionless form Vcell (j) = ∆z 3 V˜cell (j)  π  j=0 3, ˜ Vcell (j) = 2πj, j ∈ (0, nr) .   1 π(nr − 3 ), j = nr

(2.46)

To determine now the dimensionless charge q˜, we shall consider the simple case of a uniform density in the whole system (see Sec. 2.10.1). Choosing nref as the uniform density, f 3 nr2 nz NDb particles have to be distributed like ∝ r over the whole system, that is nr×nz π ∆z f 3 nr2 NDb cells. In z-direction, the particle distribution is homogeneous, so in r-direction π ∆z R nr 2 particles have to be distributed like ∝ r, over an integrated volume of 0 r dr = nr2 . Let’s assume we place all particles directly onto grid points 1, 2, ..., (nr − 1). In this case ˜e, cell (j, k) = 2π ∆z f 3 j NDb , with4 j = we do not need to care about the weighting and N R ˜e, cell (1, k) = 2π ∆z f 3 NDb , 1, 2, ..., (nr−1) and j dj = nr2 /2. For the first cell, for instance, N ∀k, while V˜cell (1) = 2π. Substituting then back to Eq. 2.41, we obtain: ˜ g2 =! q˜Ne, cell (1, k) = q˜∆z f 3 NDb , ∆t ˜ Vcell (1)

(2.47)

f2 ∆t q˜ = . f 3 NDb ∆z

(2.48)

from which we can extract

The dimensionless charges defined in init.cpp and used in density.cpp are therefore: q˜e = −

g2 ∆t , g3 NDb ∆z

q˜i = −˜ qe . 4

⇒ ne < 0

(2.49)

⇒ ni > 0

(2.50)

Why we calculate here the integral and not the sum becomes clear in Sec. 2.10.1.

CHAPTER 2. CODE DOCUMENTATION

14

2.3

Electric field calculation

According to Eq. (2.23), for all inner cells we calculate the electric field from the potential as follows, ˜r (j, k) = 1 (ϕ(j E ˜ − 1, k) − ϕ(j ˜ + 1, k)) 4 ˜z (j, k) = 1 (ϕ(j, E ˜ k − 1) − ϕ(j, ˜ k + 1)) . 4

(2.51) (2.52)

At the boundaries, these equations have to be modified such that they are in harmony with the boundary conditions for the potential. This is discussed in Sec. 2.4.

2.4

Potential and electric field boundary conditions

Five different potential and electric field boundary conditions are available in 2D ArcPIC. The user can switch between the different boundary conditions offered by changing the variable BC, which is set in input.txt. The meaning of the different options are listed below, and the boundary conditions are illustrated in Fig. 2.2. In this section, tildes are consequently omitted in order to simplify the notation, but we are still referring to the dimensionless quantities. For breakdown simulations, one would normally use BC = 4. BC 0 for ϕ|z=0 = ϕC , ϕ|z=nz = ϕA ,



∂ϕ ∂r r=0

= 0, ϕ|r=nr = 0:

• Cathode: k = 0, j ∈ [0, nr] ϕ = ϕC

(2.53)

Er (j, 0) = 0 3 1 Ez (j, 0) = ϕ(j, 0) − ϕ(j, 1) + ϕ(j, 2) 4 4

(2.54) (2.55)

• Anode: k = nz, j ∈ [0, nr] ϕ = ϕA

(2.56)

Er (j, nz) = 0

(2.57)

1 3 Ez (j, nz) = − ϕ(j, nz) + ϕ(j, nz − 1) − ϕ(j, nz − 2) 4 4

(2.58)

• Symmetry axis: j = 0, k ∈ (0, nz) ϕ(0, k − 1) − 6ϕ(0, k) + ϕ(0, k + 1) + 4ϕ(1, k) = −ne (0, k) − ni (0, k) (2.59) ∂ϕ = Er (0, k) = 0 (symmetry) ∂r r=0 1 Ez (0, k) = (ϕ(0, k − 1) − ϕ(0, k + 1)) 4

(2.60) (2.61)

2.4. POTENTIAL AND ELECTRIC FIELD BOUNDARY CONDITIONS

15

(a) BC = 0, 1 and 4: Boundary conditions with fixed potential on Cathode and Anode.

(b) BC = 2 and 3: Boundary conditions with symmetries on “Cathode” and “Anode”. Red markers corresponds to BC = 3.

Figure 2.2: Boundary conditions for field solver, and how they are applied to different grid points (Sec. 2.4).

CHAPTER 2. CODE DOCUMENTATION

16

• ‘Infinity’: j = nr, k ∈ (0, nz) Note that in this special case Ez (nr, k) = 0 ∀k, so it should be made sure the boundary conditions altogether make sense (for instance, this BC can not be applied for a vacuum boundary in j = nr, since in case of the vacuum solution for ϕC 6= ϕA we would get Ez (nr, k) = const). This is just a special case of a fixed potential at the wall j = nr, and therefore we can not require either that far away from the center ∂ϕ = 0. ∂r r=∞

ϕ(nr, k) = 0

(2.62)

3 1 Er (nr, k) = − ϕ(nr, k) +ϕ(nr − 1, k) − ϕ(nr − 2, k) 4 | {z } 4

(2.63)

=0

BC 1 for ϕ|z=0

1 Ez (nr, k) = (ϕ(nr, k − 1) − ϕ(nr, k + 1)) = 0 4 ∂ϕ = ϕC , ϕ|z=nz = ϕA , ∂ϕ = 0, =0: ∂r ∂r r=0

(2.64)

r=nr

This model is suitable for breakdown simulation, but it is an old implementation which is included in order to be able to reproduce old results. New simulations should rather use BC=4, which is more accurate. Here we apply exactly the same BCs as for BC = 0, save for j = nr: • ‘Infinity’: j = nr, k ∈ (0, nz) ∂ϕ 3 1 = − ϕ(nr, k) + ϕ(nr − 1, k) − ϕ(nr − 2, k) = 0 ∂r j=nr 4 4 Er (nr, k) = 0 1 Ez (nr, k) = (ϕ(nr, k − 1) − ϕ(nr, k + 1)) 4

(2.65) (2.66) (2.67)

BC 2 for ϕ = const at r-boundaries, ∂ϕ/∂z = 0 at z-boundaries: • Cathode: k = 0, j ∈ (0, nr) 3 ∂ϕ 1 = ϕ(j, 0) − ϕ(j, 1) + ϕ(j, 2) = 0 ∂z k=0 4 4 1 Er (j, 0) = (ϕ(j − 1, 0) + ϕ(j + 1, 0)) 4 Ez (j, 0) = 0 • Anode: k = nz, j ∈ (0, nr) ∂ϕ 3 1 = − ϕ(j, nz) + ϕ(j, nz − 1) − ϕ(j, nz − 2) = 0 ∂z k=nz 4 4 1 Er (j, nz) = (ϕ(j − 1, nz) + ϕ(j + 1, nz)) 4 Ez (j, nz) = 0

(2.68) (2.69) (2.70)

(2.71) (2.72) (2.73)

2.4. POTENTIAL AND ELECTRIC FIELD BOUNDARY CONDITIONS

17

• Rmin : j = 0, k ∈ [0, nz] ϕ(0, k) = ϕmin 3 1 Er (0, k) = ϕ(0, k) − ϕ(1, k) + ϕ(2, k) 4 4 Ez (0, k) = 0 (no gradient along z)

(2.74) (2.75) (2.76)

• Rmax : j = nr, k ∈ [0, nz] ϕ(nr, k) = ϕmax 1 3 Er (nr, k) = − ϕ(nr, k) + ϕ(nr − 1, k) − ϕ(nr − 2, k) 4 4 Ez (nr, k) = 0

(2.77) (2.78) (2.79)

BC 3 for a periodic ϕ along z: Here boundary grid points are treated as inner grid points, and implicit periodicity is used. • Cathode: k = 0, j ∈ (0, nr)   1 1− ϕ(j − 1, 0) + ϕ(j, nz − 1) − 4ϕ(j, 0) + ϕ(j, 1)+ (2.80) 2j   1 ϕ(j + 1, 0) = −{ne (j, 0) + ne (j, nz)} − {ni (j, 0) + ni (j, nz)} + 1+ 2j 1 Er (j, 0) = (ϕ(j − 1, 0) − ϕ(j + 1, 0)) (2.81) 4 1 Ez (j, 0) = (ϕ(j, nz − 1) − ϕ(j, 1)) (2.82) 4 • Anode: k = nz, j ∈ (0, nr)   1 1− ϕ(j − 1, nz) + ϕ(j, nz − 1) − 4ϕ(j, nz) + ϕ(j, 1)+ (2.83) 2j   1 + 1+ ϕ(j + 1, nz) = −{ne (j, 0) + ne (j, nz)} − {ni (j, 0) + ni (j, nz)} 2j 1 (2.84) Er (j, nz) = (ϕ(j − 1, nz) − ϕ(j + 1, nz)) 4 1 Ez (j, nz) = (ϕ(j, nz − 1) − ϕ(j, 1)) (2.85) 4 • Symmetry axis: j = 0, k ∈ [0, nz] ϕ(0, k − 1) − 6ϕ(0, k) + ϕ(0, k + 1) + 4ϕ(1, k) = −ne (0, k) − ni (0, k) (2.86) Er (0, k) = 0 (symmetry) (2.87) 1 Ez (0, k) = (ϕ(0, k − 1) − ϕ(0, k + 1)) (2.88) 4 In the special case of k = 0 and k = nz we replace k − 1 with nz − 1 and k + 1 with 1 whenever k − 1 and k + 1 is not defined and modify the right-hand side of Eq. 2.86 as follows: −{ne (j, 0) + ne (j, nz)} − {ni (j, 0) + ni (j, nz)}.

CHAPTER 2. CODE DOCUMENTATION

18

• ‘Infinity’: j = nr, k ∈ [0, nz] 3 ∂ϕ 1 = − ϕ(nr, k) + ϕ(nr − 1, k) − ϕ(nr − 2, k) = 0 ∂r j=nr 4 4 Er (nr, k) = 0 1 Ez (nr, k) = (ϕ(nr, k − 1) − ϕ(nr, k + 1)) 4

(2.89) (2.90) (2.91)

Also here, in the special case of k = 0 and k = nz we replace k − 1 with nz − 1 and k + 1 with 1 whenever k − 1 and k + 1 is not defined. ∂ϕ BC 4 for ϕ|z=0 = ϕC , ϕ|z=nz = ϕA , ∂ϕ = 0, = 0: ∂r ∂r r=0

r=nr

This represents the same boundary conditions as BC = 1, however the implementation of the boundary conditions to Poisson’s equation along r = 0 and r = nr is changed. This leads to a smaller field error, especially near r = 0. • Symmetry axis: j = 0, k ∈ (0, nz) 2ϕ(1, k) + ϕ(0, k − 1) + ϕ(0, k + 1) − 4ϕ(0, k) = −ne (0, k) − ni (0, k) (2.92) • ‘Infinity’: j = nr, k ∈ (0, nz) 2ϕ(nr − 1, k) + ϕ(nr, k − 1) + ϕ(nr, k + 1) − 4ϕ(nr, k) = −ne (0, k) − ni (0, k) (2.93)

2.5

Interpolation scheme

Since in PIC macro-quantities are calculated on the grid but particles move inside the cells, it is often necessary to interpolate a macro-quantity from the grid points to a particle position. As an example, we here consider the interpolation of the electric field. However, the same scheme is used everywhere in the code when interpolation is required for a physical quantity (e.g. for the potential). Let a particle number n have the position (pa[n].p.r, pa[n].p.z). Define the integer and fractional parts of this position as (see Fig. 2.3): j = int (pa[n].p.r)

(2.94)

k = int (pa[n].p.z)

(2.95)

hr = frac (pa[n].p.r)

(2.96)

hz = frac (pa[n].p.z)

(2.97)

The i component of the electric field at the particle position is then given by the linear interpolation of the four surrounding grid point contributions: ˜i (r, z) = (1 − hr) · (1 − hz) · E ˜i (j, k) + (1 − hr) · hz · E ˜i (j, k + 1) E ˜i (j + 1, k) + hr · hz · E ˜i (j + 1, k + 1). + hr · (1 − hz) · E

(2.98) (2.99)

2.6. PARTICLE PUSHER EQUATIONS

19

Figure 2.3: Schematic representation of the interpolation from grid points to particle position.

2.6 2.6.1

Particle pusher equations Electrostatic case

To push a particle, the electric field is first interpolated to the particle position, as discussed in Sec. 2.5. The particle is then accelerated due to the electric part of the Lorentz force, ˜r , v˜r0 = v˜r − 2E ˜z , v˜z = v˜z − 2E

(2.100) (2.101)

and particles are moved in the direction of the velocity they have (∆ri = vi ∆t), z˜0 = z˜ + ∆˜ z = z˜ + v˜z , q q 0 2 2 r + ∆˜ rx ) + ∆˜ ry = (˜ r + v˜r )2 + v˜t2 . r˜ = (˜

(2.102) (2.103)

where z˜0 and r˜0 stand for the updated quantities. The change in r-coordinate is illustrated in Fig. 2.4a. Finally, since the code is two-dimensional, we are effectively simulating an infinitely thin “slice” of a cylinder in ϑ-direction. Therefore, after pushing the particles we must rotate them back into this slice, i.e. to a given constant ϑ such as ϑ = 0. This is illustrated Fig. 2.4b. r+∆r t Defining now sin α = ∆r r0 and cos α = r0 , we may perform the rotation by transforming velocities as

2.6.2

v˜r0 = cos α · v˜r + sin α · v˜t

(2.104)

v˜t0

(2.105)

= − sin α · v˜r + cos α · v˜t .

Magnetic case

In the code, it is possible to apply an external magnetic field along the z- and/or ϑ-direction Bext = (Bz , Br = 0, Bt ) in the right-handed coordinate system (z, r, t). The non-zero components of v × B are:       vz Bz vr Bt (2.106) v × B = vr  ×  0  = vt Bz − vz Bt  vt Bt −vr Bz

CHAPTER 2. CODE DOCUMENTATION

20

(a) Updating velocities

(b) Rotating the reference frame

Figure 2.4: Illustration of the particle pusher, where particles are treated in a Lagrangian frame. The coordinates change in the direction of the velocities, that are calculated from the electric field. The new reference frame is given by the updated coordinates (z 0 , r0 ), into which velocities that are calculated in the old frame have to be transformed.

Denoting the “old” velocity vector from the previous time step as v0 = (vz0 , vr0 , vt0 ) and the updated velocity vector from the next time step as v0 = (vz0 , vr0 , vt0 ), we can update the dimensional velocity vector according to the Boris method [6] (also discussed in Sec. 1.1.3) in the following four steps (treating first electrons): 1. Electric push (I) e∆t Ez 2me e∆t Er vra = vr0 − 2me

vza = vz0 −

(2.107) (2.108)

2. Magnetic push (I) e∆t a v Bt 2me r e∆t 0 vrb = vra − {v Bz − vza Bt } 2me t e∆t a vtb = vt0 + v Bz 2me r

vzb = vza −

(2.109) (2.110) (2.111)

3. Magnetic push (II) vzc = vza −

e∆t me

vrb Bt 2 e∆t 1 + 2me {Bz2 + Bt2 }

(2.112)

vrc = vra −

e∆t me

vtb Bz − vzb Bt  2 e∆t 1 + 2m {Bz2 + Bt2 } e

(2.113)



2.6. PARTICLE PUSHER EQUATIONS

vt0 = vt0 +

21

e∆t t me

vrb Bz 2 e∆t {Bz2 + Bt2 } 1 + 2m e 

(2.114)

4. Electric push (II) e∆t Er 2me e∆t vz0 = vzc − Ez 2me vr0 = vrc −

(2.115) (2.116)

Once the velocities are updated, particle positions are updated and velocities rotated in the same way as in the electrostatic case, see Eqs. 2.102, 2.103, 2.104, and 2.105. ˜ = e ∆t2 E as before, and analoTo obtain the dimensionless equations, we rescale E ˜= gously v˜B as

e ∆t2 2me ∆z vB.

2me ∆z

Since v˜ =

∆t ∆z v,

the dimensionless magnetic field has to be defined

˜ ≡ e∆t B . B 2me The dimensionless velocity update equations are then:

(2.117)

1. Electric push (I) ˜z v˜za = v˜z0 − E ˜r v˜ra = v˜r0 − E

(2.118) (2.119)

2. Magnetic push (I) ˜t v˜zb = v˜za − v˜ra B ˜z − v˜za B ˜t } v˜rb = v˜ra − {˜ vt0 B

(2.120)

˜z v˜ra B

(2.122)

v˜tb

=

v˜t0

+

(2.121)

3. Magnetic push (II) v˜zc = v˜za −

˜t 2˜ vrb B 2 ˜ +B ˜ 2} 1 + {B z

(2.123)

t

˜z − 2˜ ˜t 2˜ vtb B vzb B ˜2 + B ˜ 2} 1 + {B z t b ˜ 2˜ vr Bz v˜t0 = v˜t0 + ˜z2 + B ˜ 2} 1 + {B

v˜rc = v˜ra −

(2.124) (2.125)

t

4. Electric push (II) ˜z v˜z0 = v˜zc − E ˜r v˜r0 = v˜rc − E

(2.126) (2.127)

˜ by E ˜ion and multiply B ˜ by the factor To apply the same equations also for ions, replace E me −dt_ion Mi (one factor of dt_ion is contained already in the ion velocity).

CHAPTER 2. CODE DOCUMENTATION

22

2.7 2.7.1

Collision routines Coulomb collisions for e- –e- and Cu+ –Cu+ pairs

Coulomb collisions between the same species are treated with the Takizuka-Abe binary collision model [8] in the function coll_el_knm_2D(), with the implementation described in the PhD thesis of K. Matyash [15]. In the model, the changes in velocity due to Coulomb collisions are characterised by the angle ϑ, and the variable δ ≡ tan ϑ/2 is chosen randomly from a Gaussian distribution such that the variance5 hδ 2 i is given by hδ 2 i =

qα2 qβ2 nL ln Λ 3 8πε20 m2∗ vrel

∆tcoll

(2.128)

where α and β denotes the species, nL is the lowest of different species’ densities, ln Λ is the Coulomb logarithm, m∗ is the reduced mass, and vrel the relative velocity. In our case α = β = e− or Cu+ , therefore nL = nα = nβ , and m∗ = Mp /2, where Mp is the electron or ion mass. When rescaling the above expression to dimensionless quantities, just as before, vrel = ∆tcoll ≡ ∆t · ncoll_el, masses are expressed in me , and the density in a cell can be calculated based on the number of simulation particles in the cylindrical geometry as ∆z ˜(ion) , ∆t(ion) v

n=

NSP Np NSP Np , = Vcell π(2j + 1)∆z 3

(2.129)

where Np is the number of simulation particles in a given cell determined by points {(j, k), (j + 1, k + 1)}, recorded in the variable ordcount[j,k], and the volume of that cell6 is Vcell = 2 R k+1 R j+1 rj+1 −rj2 dz = 2π∆z 2πr dr = π(2j + 1)∆r2 ∆z. 2 k j Substituting these into the expression for hδ 2 i we arrive at   qα2 qβ2 nL ln Λ e4 me 2 2 hδ i = × 3 ∆tcoll =vel_scale 2πε2 m2 Mp 8πε20 m2∗ vrel 0 e Np ∆t4 NSP ln Λ 3 ncoll_el 3 , 3 ∆z π(2j + 1)∆z v˜rel

(2.130)

3 . Note that where vel_scale is 1 for electrons and dt_ion3 for ions due to the scaling of vrel in order to keep the collision frequency for simulated particles the same as for real particles, we used qα = ±e and mα = me− (Cu+ ) and did not multiply charges and masses by NSP .

Since

e4 ε20 m2e

=

4 ωpe 2 nref

f we can simplify the above to and ωpe ∆t = ∆t,

ln Λ hδ i = vel_scale 2 2π | 2



me Mp

2 f 4 ∆t NSP 1 Np 6 n2 λ6 ncoll_el 2j + 1 v 3 ˜rel f {z ∆z ref Db }

(2.131)

≡Acoll

which are exactly the factors appearing in the code. Note that hδi = 0. The cell volumes should not to be confused with the Verboncoeur volumes, which are calculated on the grid points, not between them. 5

6

2.7. COLLISION ROUTINES

2.7.2

23

Other collisions

In this subsection, we treat the amplitude scaling in collision routines which apply a known (obtained by measurement or fit) cross-section σ. Linear interpolation to cross-section data is carried out in init.cpp. The following routines in collisions.cpp all belong to this category: coll_ion_neutral_noSP_2D(...) , treating ion-neutral elastic collisions with charge exchange, Cu+ + Cu −→ Cu + Cu+ coll_n_n_2D(...) , handling neutral-neutral elastic collisions, Cu + Cu −→ Cu + Cu coll_el_all_fake_2D(...) , covering electron-neutral elastic collisions, e− + Cu −→ e− + Cu coll_el_neutrals_2D(...) , treating the impact ionisation process, e− + Cu −→ 2 e− + Cu+ . In all these routines, the probability that a collision occurs is approximated as Pi = 1 − e−σnCu vrel ∆tcoll ≈ σ nH vrel ∆tcoll ,

(2.132)

where nH is the number density of the most populous species in that cell. The requirement that this probability should remain unaltered when rescaling to dimensionless quantities determines the scaling of σ. As before, vrel = ∆t∆z v˜(ion) , ∆tcoll ≡ ∆t · ncoll_ion, and the dimensional (ion) density in a given cell is n = through demanding σ ˜ Np v˜ · 1 = CScoeff · |

NSP Np Vcell

NSP Np . π(2j+1)∆z 3

=

We then can determine CScoeff ≡ σ ˜ /σ

π(2j + 1)∆z 3 ∆t 1 · · (σ nH vrel ∆tcoll ) NSP ∆z ∆t ncoll_ion {z }

(2.133)

!

=1

nref λ3Db ncoll_ion NSP ncoll_ion = 2 g2 λ2 NDb π(2j + 1)∆z π(2j + 1)∆z Db nref λDb ncoll_ion 1 = g2 NDb 2j + 1 π ∆z

⇒ CScoeff =

(2.134)

In the code, the correction with the factor 1/(2j + 1) happens separately when calculating the probability in Eq. 2.132.

2.7.3

Parallelization using OpenMP

In order to speed up processing, the neutral-neutral collisions were parallelized using OpenMP [16], which is a compiler-based parallelization technique for SMP computer systems with shared memory. This allows for easy parallelization of existing codes, requiring minimal rewrites. The Neutral-neutral collision routine was targeted, as profiling revealed that the largest amount of CPU time was spent in the function coll_n_n_2D(...). This function was therefore adapted such that when the number of neutral superparticles in the system exceeds the limit

CHAPTER 2. CODE DOCUMENTATION

24

ϕsr,cathode

ϕsr

ϕsr,anode

+ + + + Cathode

I

+ + +

+ +

+

+ + +

z

z Cathode

Anode

(a) Induced charge on electrodes (red) by a negative charge in the gap (blue).

Anode

(b) Induced current in circuit due to movement of charge in gap.

Figure 2.5: Illustration of the Shockley-Ramo theorem as applied to two parallel plates with free charges between them.

set in COLL_N_N_2D_OMP_MINPARTICLES (a #defined constant in colls.h with default value 1000), the work is split across numParaThreads threads. The variable numParaThreads is set in input.txt. This splitting takes advantage of the fact that the particles in each cell are handled independently by the collision algorithm, by tasking each CPU with a subset of the total number of cells. The load is balanced by having each CPU handling an approximately equal number of particles. The parallelization speed-up was tested by pre-filling the gap with a large number of neutrals, using purely absorbing wall boundary conditions, and running the code with 1, 2, 3, and 4 CPUs. For this scenario, an almost linear increase in number of time steps per second was observed.

2.8

Circuit models

The electric potentials ϕA and ϕC on the anode and cathode are set by the currently loaded circuit model at the end of each time step. 2D ArcPIC contains several circuit models, all implemented in circuit.cpp. The circuit models are expected to provide the updated potentials, based on the gap current Igap (t) and the time since the simulation was started.

2.8.1

Shockley-Ramo gap current calculation

In order to accurately estimate the gap current, Igap is calculated from the flow of the surface charges on the electrodes induced by the charges in the gap (Fig. 2.5a). This is done using the Shockley-Ramo theorem [17, 18]. For the geometry used in this simulation, the induced charge on the cathode from the free charges in the gap is given as Qind.cat. = −

X i

qi ϕsr (zi ) where ϕsr (zi ) = 1 −

zi Zmax

,

(2.135)

2.8. CIRCUIT MODELS

25

the sum runs over all charged particles in the gap, zi is the position of particle i and qi is the charge. The induced current is then given as dQind.cat. , (2.136) dt where the sign comes from that the direction of positive current is defined such that electrons flowing from the cathode to the anode creates a positive current. In the code, the charge transport is calculated as ∆Qind.cat. , (2.137) Igap ≈ − ∆t where ∆Qind.cat. is the charge transported in the time step between t = ti−1 and t = ti . For this to work correctly when injecting and removing particles from the domain, the induced charge described in Eq. (2.135) must take into account that the Shockley-Ramo potential inside the cathode is ϕsr = 1. This means that charges which at t = ti are in the gap and are inducing a charge −qϕsr (z) on the cathode, were at t = ti−1 already inducing a charge −1q. In practice, this is handled by adding a term −Qinj. cat. (ti ) for charge transport through the cathode plane. The charge transport is thus given as Igap = −

∆Q(ti ) = −Qind. cat. (ti ) + Qind. cat. (ti−1 ) − Qinj. cat. (ti ) ,

(2.138)

 where Qinj. cat. (ti ) = −e Ne− inj. − Ne− removed − NCu+ inj. + NCu+ removed is the total charge of the particles injected and removed by the surface model at the cathode (see Sec. 2.9), and e the elementary charge. The use of this method to calculate the current has the effect of smoothing the current signal, compared to using the charge crossing the surfaces per time step. This can be important if using circuit a model with little or no gap capacitance, such as FixedVoltage_resistor. A possible improvement of the method, which eliminates the need to sum over all particles as in Eq. 2.135 by directly using the output of the field solver, is outlined in Ch. 4.

2.8.2

Rescaling to dimensionless units

In order to implement a circuit model, we need to rescale it to dimensionless units. To do this, the voltages and potentials are rescaled using Eq. 2.28, and currents to number of superparticles per time step, such that the total current is rescaled by the factor r f · NDb f · NDb I˜ [1] ∆t 1 me ∆t ∆t = = = 2.693026 × 105 , (2.139) 3/2 I [A] eNSP ε0 e (Tref /e) (Tref /e)3/2 where the reference temperature is given in units of [Tref /e] = V . Charge, capacitance, and resistance have to scale with them. With the reference density given in units of [nref ] = 1/cm3 , they scale as follows: √ √ ˜ [1] Q I˜ 1 1 e NDb · nref = = = 3/2 3/2 Q [C] I ∆t eNSP ε0 (Tref /e) √ 10 NDb · nref = 1.519260 × 10 , (2.140) (Tref /e)3/2 ˜U C˜ [1] Q 1 = = ˜ C [F ] QU eNSP

f ∆z f ∆t

!2

√ Tref e = 3/2 e ε0

f ∆z f ∆t

!2

√ NDb · nref p Tref /e

CHAPTER 2. CODE DOCUMENTATION

26

10

= 1.519260 × 10

˜ [1] ˜I R U = = R [Ω] U I˜

f ∆z f ∆t

!2

√ NDb · nref p , Tref /e

!2

p r f Tref /e e ∆t 1 eNSP = ε0 2 Tref /e ∆t me ∆z NDb f p f Tref /e ∆t . = 3.713294 × 10−6 2 NDb f ∆z

2.8.3

f ∆t f ∆z

(2.141)

(2.142)

The output file circuit.dat

The Circuit class writes some output every7 file_timestep number of time-steps into the ASCII file circuit.dat. By default, the output line contains the current time step, ∆Q as ˜ = ϕ˜A − ϕ˜C . In a number of superparticles since last output, and the current gap voltage U addition to the per-time-step data, header lines starting with “#” are written at the top of the file which describe the file format and name the model. Some of the circuit models override the default output, adding more columns. It is expected by the analysis tools that they still keep the default format for the first columns.

2.8.4

Models

Several classes inherits from and implements the virtual methods from the abstract class Circuit. The most used models corresponding to these classes are listed below. Which model is loaded is determined by input.txt, where the circuit model section is parsed by readInputSection(...). Calling the right circuit constructor is handled by the static method Circuit::LoadCircuit(...), meaning that all circuit models must be listed in this method to be usable. The most important functions which must be implemented are: Constructor A constructor with one input argument std::vector options should be provided. This is called from the static function Circuit::LoadCircuit(...), which uses readInputSection(...) to parse input.txt and load the correct model. The options argument contains a list of parameters loaded from input.txt, containing data such as initial potentials, parameters of the circuit elements and similar things. init() Re-scales parameters read from input.txt to correct units. timestep(double deltaQ, unsigned int nstep) Called from main() at every time step. The arguments are the number of superparticles deltaQ transported across the gap during the previous time step, and the current time step nstep. This function is expected to update the member variables U0 and UNz, which stores the cathode and anode potential. The most important circuit models are: 7

The variable file_timestep, which is set by the constructor of the Circuit-implementation, is normally set to 1. This means the output file is written to every time step of the simulation. The constructor often sets this variable from input.txt.

2.8. CIRCUIT MODELS

27

Igap Vgap

Anode PIC Cathode

Icirc

R Cgap

U

(a) Implemented equivalent circuit, where the PIC model is acts as a “black box” which sees some electric potential Vgap and creates some current Igap .

dQgap dt

Icirc. PIC Cgap =

Qgap Vgap

R

Vgap

Icap = Igap − Icirc = −

Igap

Anode

+Qgap

U

−Qgap Cathode (b) Physical process – the Igap is the flux of charged particles between the cathode and anode just around the spark, ignoring the displacement current in this area.

Figure 2.6: Circuit model FixedVoltage_resistorCapacitor.

FixedVoltage_resistorCapacitor This model, which is illustrated in Fig. 2.6, is the one most closely matching the current setup of the DC spark experiment. It consists of the local discharge gap capacitance Cgap , which is connected to an external voltage supply through a resistor. This circuit can be described by the equations Z t 1 Vgap (t) = Vt=0 + Icirc (t) − Igap (t) dt , (2.143) Cgap 0 where Icirc =

U −Vgap R

ϕA (ti ) =

and Vt=0 = U . In the code, this is implemented as

∆t (Icirc (ti−1 ) − Igap (ti )) + ϕA (ti−1 ) and Cgap

ϕC = 0 ,

(2.144)

and the current is then updated as Icirc (ti ) =

U − ϕA (ti ) . R

(2.145)

This model adds one column to circuit.dat, containing the circuit current Icirc in the same units as Igap . FixedVoltage_resistor This model is a simpler version of FixedVoltage_resistorCapacitor which lacks the gap capacitance. It is implemented as ϕA (ti ) = U − Igap (ti ) R

and

ϕC = 0 .

(2.146)

CHAPTER 2. CODE DOCUMENTATION

28

Figure 2.7: The discharge gap as part of an RC-circuit. The energy available for the discharge is stored in the external capacitor.

TwoCapacitors This model is the same as found as in v2.07, and used when comparing results to the old version. It is illustrated in Fig. 2.7, and is similar to the old CERN DC spark experiment circuit. Initially, the capacitor Cext is assumed to be fully charged with the initial voltage ϕA − U0 given by the user in input.txt, Q0 = (UN z − U0 )Cext . The capacitor charge Q(t) and the voltage over the discharge gap U (t) are then updated every time step as follows: Q(ti+1 ) = Q(ti ) − Igap (ti ) ∆t , Q(ti+1 ) U (ti+1 ) = − Rext Igap (ti ) , Cext

(2.147) (2.148)

with the current I(ti ) calculated as described in Sec. 2.8.1. After a time tchange which is set in input.txt, the capacitance is changed to a larger value (set in input.txt), emulating the time delay before the charges on an external capacitor becomes available. ˜ ext This model adds two columns to circuit.dat: The current rescaled capacitor charge Q and a flag indicating whether t > tchange or not.

2.9

Surface models for particle boundary conditions

As for the circuit models (Sec. 2.8), the particle boundary conditions are implemented using a modular system. In this case, each model is contained in a class, and these classes all inherit from and implement the virtual functions in the abstract class ArcBounds. The most important virtual functions the models classes must implement are remove_e (...) / remove_i (...) / remove_n (...) and inject_e (...) / inject_i (...) / inject_n (...). The system is designed to easily facilitate addition of new models for the surface physics, without needing to change other parts of the code. The particle removal functions are called just after the particle pushers (as described in Fig. 2.1), and remove all particles which have been pushed outside the edges of the simulation domain. These functions may be implemented directly by the model (for example if the removed particles affect the injection of particles through sputtering), or by inheriting from the class ArcRemover, which provides a minimal implementation. The particle injection functions are called on the injection time steps, as shown in Fig. 2.1. Their role is to add more particles to the particle arrays, using algorithms modeling FowlerNordheim field emission, sputtering of neutrals etc.

2.9. SURFACE MODELS FOR PARTICLE BOUNDARY CONDITIONS

2.9.1

29

Output files

The surface models generate several output files. Most of this output comes from ArcBounds directly, such that implementing classes only need to fill arrays and counters which are then read by the file writing functions in ArcBounds. arcbounds.dat contains information about how many superparticles of each type are injected and removed at each simulation domain border. It is written every file_timestep, which is an object variable set by the constructor of the implementing class, often using input.txt. It is usually set to 1. removed_i.dat lists the position and velocity of all ion (Cu+ ) superparticles which are removed from the simulation domain, and when they were removed. removed_n.dat lists the position and velocity of all neutral (Cu) superparticles which are removed from the simulation domain, and when they were removed. removed_e.dat lists the position and velocity of all electron superparticles which are removed from the simulation domain, and when they were removed. This file can become extremely large, and it may thus in many cases be desirable to skip downloading it if the simulation is ran on a remote machine (See Ch. 3). jhist_cathode.dat contains one line written every global output step, which stores the sum of the charge of the superparticles removed/injected from every cell on the cathode in the time interval between the output steps. This is useful for calculating the current density as a function of time and position on the electrode. jhist_anode.dat contains one line written every global output step, which contains the sum of the charge of the superparticles removed/injected from every cell on the anode in the time interval between the output steps. This is useful for calculating the current density as a function of time and position on the electrode. In addition, the classes may generate extra output files by implementing the virtual function writeFile_extras(...), by default defined as a virtual do-nothing method in ArcBounds.

2.9.2

Injection of particles

Most of the models described in Sec. 2.9.4 use the same injection scheme, which is described in this subsection. The initial velocity of the electrons is for all three pcomponents sampled from a Gaussian random distribution with σ = 0.01 vth , where vth = Tref /me . This distribution is truncated below σ = 5, such that the velocity away from the initial point |v|2 > 0.05vth . For the neutrals, the initial velocity is determined the same way as for the electrons, using a temperature Tion = Tref · Ti_over_Te, where the variable Ti_over_Te is set in input.txt. The initial z-position of the particles is normally given as z = vz ∆t , while the way the radial position is chosen is different for different emission models.

(2.149)

CHAPTER 2. CODE DOCUMENTATION

30

For field emitted electrons, it is in some models possible to randomly sample the injection time across the injection time step, and use this to determine the initial position (and update the initial velocity). To do this, it is tracked for a fractional time step using an Euler method, taking into account that the Boris method treats the velocity and position at different points in time. Using this scheme, the position of the particle is updated to r = r0 + (1 − R)∆tinj v0

(2.150)

where r0 is the original position (which is on the surface of the emitter), ∆tinj the injection time step, and R a uniformly distributed random number. The velocity is then updated as   ∆tinj 1 v = v0 + (1 − R) ∆t a , (2.151) − ∆t 2 where v0 is the initial velocity sampled as described above, R is the same random number as in qe Eq. 2.150, and the acceleration vector is given as a = m Ez eˆz , where the field Ez eˆz is sampled e on the nearest mesh grid on the surface.

2.9.3

Fowler-Nordheim field emission

When simulating vacuum arc discharges, field emission is an important process. The 2D ArcPIC code models this by injecting electrons according to the Fowler-Nordheim formula [19], using the Wang and Loew approximation [20]. The current density of the emitted electrons are thus given as ! 2 1.5414 × 10−6 Eloc 6.8309 · 109 φ3/2 v(s) jFN = exp , (2.152) φ · t2 (s) Eloc where the field emission current density is measured in [jFN ] = A/m2 , the local surface field in [Eloc ] = V/m, and the work function in [φ] = eV; the elliptical functions are approximated with √ Eloc 2 −5 t(s) = 1 and v(s) = 0.956 − 1.062 · s , where s = 3.7947 · 10 φ . Inserting the value φ = 4.5 eV and simplifying, we arrive at the form used in the code   62.338 2 jFN = 4.7133 × 109 Eloc · exp − , (2.153) Eloc where [jFN ] = A/cm2 and [Eloc ] = GV/m.

2.9.4

Models

Several models are provided with the code, all inheriting from the ArcBounds class. The most commonly used models are presented below. Which model is loaded is determined by input.txt, where the particle boundary parameters section is parsed by readInputSection(...). There is significant overlap in both functionality and code between the different models, as they are frequently created by branching off- and borrowing pieces from each other. This development method is well supported by the ArcBounds framework. As mentioned above, the class ArcRemover provides a simple “default” implementation of remove_e(...), remove_i(...) and remove_n(...). In this implementation, these functions only do the minimum required, which is to delete the particles which have left the domain,

2.9. SURFACE MODELS FOR PARTICLE BOUNDARY CONDITIONS

Remit Rtip

r

Secondary electron yield SEY = 0.5

Anode

Cathode

je (r) = jFN (βflat · Ez (r)) jCu (r) = rCu/e · je (r)/e

Sputtering of neutrals (low intensity model) Anode

Cathode

r

31

Enhanced sputtering

2 Itip = πRtip · jFN (βtip · Ez (0)) ICu = rCu/e · Ie /e

z (a) Field emission and neutral evaporation.

z (b) Sputtering and secondary electron yield.

Figure 2.8: The emission processes modeled in ArcBoundsOriginal and ArcBoundsOriginalNewHS.

and register this in various arrays which are used for writing output files and for current calculation. They thus also provide a template for writing such functions from scratch. To use this implementation, the new ArcBounds class should inherit from ArcRemover and not directly from ArcBounds. ArcSimple This very simple model injects a fixed number of electron, ion and neutral superparticles with zero velocity into a random location which is uniformly sampled in the range z ∈ (0, ∆z), r ∈ (0, ∆z) (not taking into account that dV ∝ r dr) at every injection time step. How many particles to create of each type is set in input.txt. This model is mainly intended for testing of the ArcBounds framework. Removal of particles leaving the simulation domain is provided by ArcRemover. ArcBoundsOriginal This model implements field emission from the cathode, single electron yield and sputtering, as illustrated in Fig. 2.8. It corresponds to setting PBC=0 and PIS=0 in v2.07, and was used for comparing results with previous versions while updating other parts of the code. The model ArcBoundsOriginalNewHS is very similar, as it is based on ArcBoundsOriginal. A cylindrical field emitter of radius Rtip is assumed to be located at the cathode, around the symmetry axis. The field emitter is assumed only implicitly through particle inlet conditions; no assumptions on a particular shape or other properties of the field emitter have been made. This field emitter is assigned a field enhancement factor βtip , and serves as a source of electron field emission. Since in reality surface roughness and smaller emitters are practically unavoidable, we assign also an average field enhancement factor βflat to the “flat” cathode surface, i.e. outside the central “tip”. This area also serves as a source of electron field emission.

32

CHAPTER 2. CODE DOCUMENTATION

For the current originating from the “tip”, we used Rtip to convert from injection current 2 j density to total current as Itip = πRtip FN (βtip Ez (r = 0, z = 0)), and a possibly slightly larger radius of Remit for the actual injection of particles. The current density jtip (Elocal ) is calculated as described in Sec. 2.9.3. This is done in order to properly resolve the density, avoiding numerical stability problems in the small, discretized region around the field emitter.. Injection is carried out with a uniform number density of injected particles above the injection area, such that more superparticles are injected into the simulation cells which are further away from the symmetry axis, as the size of the volume elements increase with the distance from the axis. The generated particles are then injected as described in Sec. 2.9.2, with or without a fractional time step push (to be chosen by the user in input.txt). In addition to the injection of particles emitted from the “tip”, we also inject particles emitted from the rest of the cathode. As for the tip, the current density is calculated around each mesh point on the cathode surface as described in Sec. 2.9.3, and then multiplied with the area of the emitting area to find the current. In this model, flat surface injection is applied for r > Remit . The generated particles are then injected as described in Sec. 2.9.2, with or without a fractional time step push. Since the field emitters are exposed to a very high Eloc = βEz and a current density jFN that can heat up the emission region significantly, the field emitter also serves as a source of neutrals and ions [21] through a mixture of thermal and field evaporation [22, 23]. We here assume a simple model in which the field-assisted thermal evaporation of Cu neutrals follows the field emission current in a constant ratio rCu/e . The injection of evaporated neutrals is handled as described in Sec. 2.9.2, and they are all injected from r < Remit . Besides field emission and neutral evaporation, several sputtering phenomena are taken into account: • Cu and Cu+ can sputter Cu at both anode and cathode with an experimentally measured, energy-dependent physical sputtering yield [24]. • High-energy ions from the plasma impinging at the cathode with a high flux can lead to heat spike sputtering with an enhanced sputtering yield previously calculated with MD simulations [25]. Above an ion flux of 6 × 1025 1/s/cm2 incident on the cathode surface, we inject a fixed number of neutral superparticles Yenh. = 1000 per time-step, regardless of the number of incident particles. This overrides any “normal” sputtering yield from ions, but not from the neutrals. • High-energy ions may cause the sputtering of electrons at the cathode. Above an ion impact energy of 100 eV, we apply a constant secondary electron yield (SEY) of 0.5. The particles are injected with the velocity distributions as described in Sec. 2.9.2, but no fractional time step push is applied. For low intensity sputtering and SEY, the injection position is either the same as the incoming particle causing the sputtering. In the case of heatspike sputtering, the position is sampled from a radial Gaussian distribution with σ = r˜heatspike · ∆r, where r˜heatspike is the outermost edge of the outermost cell where the ion flux was above the threshold. The ArcBoundsOriginal model generates one extra output file arcbounds_original.dat, which contains additional information such as how many particles was injected from the tip, r˜, Eloc etc. It is written to every file_timestep.

2.10. INITIAL PARTICLE DISTRIBUTION

33

ArcBoundsOriginalNewHS This model is mostly identical to ArcBoundsOriginal, the differences being in the flat emission, the enhanced sputtering process, and that more data is written to arcbounds_original.dat. The model is illustrated in Figure 2.8. The electron emission works in the same way as for ArcBoundsOriginal, except that the border Rborder from which the flat emission is applied can be set by the user in input.txt. Rborder may be set equal to Remit to produce the same behavior as in ArcBoundsOriginal, or to any other value. The standard setting is Rborder = Rtip , which produces an overlap of the injection areas for the tip and flat injection while conserving the total current. The flat emission also uses the same model as for ArcBoundsOriginal, except that it decouples the evaporation from the “tip” and “flat” surface. Here, the “neutral current density” jCu is determined by the electron current density je at the same point. Thus neutrals evaporated due to electron tip current is injected at r < Remit , while those evaporated due to the flat current is injected at r > Rborder . The ArcBoundsOriginalNewHS class implements three different high intensity sputtering models, selectable in input.txt. Of these, the one named “c” is almost always used, and is the one described here. Here, the enhanced sputtering process is modeled in a similar way as in ArcBoundsOriginal, finding the outermost edge r˜heatspike of the outermost cell where the flux of ions and neutrals with kinetic energy8 > 23.383 eV was greater than the threshold jenh. sput. . These particles are counted, and a sputtering yield Yenh. is applied to decide how many particles to inject. However, here the yield is a ratio determining how many new particles (on average) to inject per incoming particle. This sputtering is also applied in adittion to the “normal” low-intensity sputtering, for both ions and neutrals. The sputtered neutrals are injected with a Gaussian distribution in r˜, such that the σ = r˜heatspike for this distribution. FlexFN models The FlexFN- and FlexFN2 classes are not in themselves instantiable, but provide a framework for easily implementing field emission in a new surface model. They contain methods for calculating and injecting the field-emitted electrons, given the α and β for each grid point (FlexFN) or for an arbitrary large set of “rings” defined between two radii (FlexFN2). The factor α is a scaling factor for the current density, such that the injected current density is jinj. = α jFN (β Ez (r)).

2.10

Initial particle distribution

As for the circuit- (Sec. 2.8) and surface- (Sec. 2.9) models, there is also a modular system for injecting an initial particle distribution. In this case, the models are again contained in separate classes, which inherit from the abstract class InitialParticles and implement it’s virtual functions. There is currently only one implementation, UniformRestricted, which is described below. Whether to load this model or not, and the parameters for the model, are controlled through input.txt. Here readInputSection(...) is used to parse the parameters from the initial particle distribution section. 8

This is also the threshold energy for the low-intensity sputtering.

CHAPTER 2. CODE DOCUMENTATION

34

2.10.1

Generating a uniform density in 2D PIC

Over the whole domain (r ∈ [0, Rmax ], z ∈ [0, Zmax ]), we would like to distribute Z 2 NREAL = nref dV = πRmax Zmax nref

(2.154)

real particles with density nref , corresponding to nref NREAL 2 = πRmax Zmax NSP NSP (nr)2 ∆r2 nz ∆z f 3 nr2 nzNDb =π NDb = π ∆z 3 λDb

NSIM =

(2.155) (2.156)

simulation particles distributed with a uniform density, where NDb is defined as the number of superparticles in a Debye cube corresponding to the density nref . To obtain a uniform density, the particle distribution has to have the same dependency on the coordinates as the volume element: f (r, z) = c r (2.157) where c is a constant that is determined through the normalization requirement: Z nz Z nr c ! f 3 nr2 nzNDb dz dr c · r = nr2 nz = π ∆z 2 0 0 f 3 NDb . ⇒ c = 2π ∆z The probability distribution function of one simulation particle is then    2 2 1 p(r, z) = 2 r = r ≡ pf (r)pg (z) 2 nr nz nr nz

(2.158) (2.159)

(2.160)

To generate particles in such a distribution, we calculate the cumulative distribution functions Z r √ 2 0 0 r2 [0, 1] 3 F (r) = r dr = −→ F −1 (RAND) = nr · RAND (2.161) 2 2 nr 0 nr Z z z 1 [0, 1] 3 G(z) = dz 0 = −→ G−1 (RAND) = nz · RAND (2.162) nz nz 0 f 3 nr2 nzNDb particles with following positions: In conclusion, we need to generate π ∆z z˜ = RAND · nz √ r˜ = RAND · nr

2.10.2

(2.163) (2.164)

The UniformRestricted initial particle distribution

This class fills a given sub-volume between a minimum and a maximum r and z with ions, electrons an/or neutrals with a given density and temperature. The positions are generated as discussed inpSec. 2.10.1, while the velocity components arep generated p as N (µ = 0, σ = σv˜ ). Here σv˜ = Tinj /Tref v˜te for electrons, and σv˜ = dt_ion me /Mi Tinj /Tref v˜te for ions or neutrals.

2.11. DENSITY CALCULATION

35

Figure 2.9: Density calculation on a given grid point (marked with a cross). Contributions come from four neighbouring cells, taken into account with the same V˜cell (j). Vice versa, one particle contributes to the density of four surrounding grid points, with different V˜cell (j) depending on the j value of the grid points.

2.11

Density calculation

As it was mentioned in the introduction (Sec. 1.1), in order to ensure self-consistency, extrapolation of quantities from particle position to grid points has to be carried out in analogy with the interpolation from grid points to the particle position. We therefore define hr = frac (pa[n].p.r) and hz = frac (pa[n].p.z) as in Sec. 2.5. The contribution from each particle to the density of the four closest grid points is then (as illustrated in Fig. 2.9) (1 − hr) (1 − hz) V˜cell (j) (1 − hr) hz dens(j, k + 1)+ = V˜cell (j) hr (1 − hz) dens(j + 1, k)+ = V˜cell (j + 1) hr hz dens(j + 1, k + 1)+ = , ˜ Vcell (j + 1) dens(j, k)+ =

(2.165) (2.166) (2.167) (2.168)

where V˜cell (j) is the Verboncoeur volume of the PIC cell surrounding the grid point (see Sec. 2.2.4). In the end, according to Eq. 2.41, this is to be multiplied by the species’ dimensionless charge q˜p

2.12

Monitoring energy conservation

At any instant, the total energy in the system can be followed via energy_total_2D() in engy.cpp. The total energy is the sum of kinetic and potential energy, X X W =T +V = mi vi,2 n + qi ϕi, n (ri, n ) , (2.169) i,n

i, n

where we sum over species (i) and number of particles within that species (n). In the following, we will split the total energy into particle (e, n, i) and electrostatic potential energy (p) contributions: Wtot = We + Wn + Wi + Wp . (2.170)

CHAPTER 2. CODE DOCUMENTATION

36

By following the values of Wtot , energy conservation may be monitored directly in noncollisional systems, if boundary conditions are appropriate (e.g. reflection is energy conserving; for other boundary conditions, the energy injected or lost at the walls has to be calculated). In systems with external particle and/or energy sources, Wtot can for instance serve as an indicator for reaching steady-state, while changes in the ratio of kinetic and potential energy can for example indicate the build-up of a sheath. Note that since every simulation particle represents NSP real particles, the kinetic energy contributions will scale as ∝ mi ∝ NSP , whereas the electrostatic potential energy scales as 2 , and the electric field energy density scales as ∝ 1. To compare these quantities, ∝ qi qj ∝ NSP we have to make sure that they all scale in the same way with respect to NSP . In the following, W∗ (∗ = e, n, i, p) stands for the physical energy, a sum of energies over real particles. 1. Electrons (all velocities dimensionless) We v2 pa[].p.vz2 + pa[].p.vr2 + pa[].p.vt2 = 2 = NSP Tref 2vt, e 2vt,2 e

(2.171)

2. Neutrals (all velocities dimensionless) v2 pa[].p.vz2 + pa[].p.vr2 + pa[].p.vt2 Wn = 2 = me 2 NSP Tref 2cs, Cu 2M v Cu t, e

(2.172)

3. Ions (all velocities dimensionless) Wi pa[].p.vz2 + pa[].p.vr2 + pa[].p.vt2 v2 = = 2 NSP Tref 2cs, Cu+ 2 Mme+ vt,2 e

(2.173)

Cu

4. The electrostatic potential energy of a super-point-charge is V (r) = qϕ(r): g2 Wp qϕ(r) (2.28) q ∆z = ϕ(r) ˜ , = g2 NSP Tref Tref e ∆t |{z}

(2.174)

±1

where the sign is + for ions and − for electrons, and the potential is interpolated from the grid to the position of the charge, since that will give us the maximum of the potential energy. 5. Alternatively, the energy contained in the electric field (Wf ) may be calculated, too. However, the electric field is a secondary quantity, since it is derived from the potential, therefore, the error will be bigger in this case. To follow the same scaling as above, this term needs to be rescaled also to NSP P 1 2 Wf j,k (E(j, k) ∆V (j)) 2 ε0 = , (2.175) NSP Tref NSP Tref where E(j, k) is the electric field on the grid and V (j) is the corresponding volume of the cell, calculated on the grid. In terms of dimensionless variables: P 1 3 2˜ Wf j,k (E(j, k) Vcell (j)) 2 ε0 ∆z = (2.176) NSP Tref NSP Tref

2.13. OUTPUT FILES

37

1 f3 3 P i,j,k 2 ε0 ∆z λDb

=

h



i2 f 2me ∆z 2 E 2 (j, k) V ˜ ˜cell (j) λ ω e ∆t f 2 Db pe i nref λ3Db NDb Tref

X f5 1 ε0 4m2e ∆z 2 2 ˜i2 (j, k)V˜cell (j)) = NDb λ2Db ωpe ωpe (E 4 2 2 nref Tref e ∆t f | {z } i,j,k

(2.177)

(2.178)

=Tref /me

f 5 e2 n ∆z

X ε0 4me 1 ref ˜i2 (j, k)V˜cell (j)) (E = NDb 2 nref e2 ∆t f 4 ε0 me i,j,k

(2.179)

f5 X 4NDb ∆z ˜i2 (j, k)V˜cell (j)) . = (E 4 f 2∆t i,j,k

(2.180)

f ∆z, f we finally obtain Given that the dimensionless thermal velocity v˜te = ∆t/ f3

4NDb ∆z X X Wf f2 ∆t ˜i2 (j, k)V˜cell (j)) = 4 ˜ 2 V˜ ) . = ( E (E NSP Tref 2vt,2 e 2˜ q vt,2 e i,j,k

(2.181)

i,j,k

Note once more that energy conservation can only be demanded, if there are no sinks, no sources, no collisions, no externally pumped-in energy and so forth.

2.13

Output files

In addition to the output files written by the circuit models (Sec. 2.8) and surface models (Sec. 2.9), there are common outputs describing the fields, particle positions and densities, and other data. These files are described in this section. Most of these outputs are written only every output time step, and how often this happens is −1 set by the variable dt_out from input.txt in units of ωref . Also, some of the outputs are averaged across several time steps, and how many is set by the variable av_time from input.txt −1 in units of ωref . These are converted to number of time-steps in calc_parameters_2D() in init.cpp as   av_time nav_time = int , (2.182) ωref + 0.1 and similarly for nav_dt from dt_out.

2.13.1

Output from main()

Two output files are produced from the main() function: mainStats.dat lists, for every time step of the simulation, the step number and the corresponding time in nanoseconds, the wall time needed to advance from the previous time step, and the number of superparticles of each type which is currently tracked by the program timeIndex.dat provides a list of the main output steps, which are the time-steps when the fields etc. are written to files. Every line in the file lists the step number and the corresponding time in nanoseconds. This is useful for the analyses, which can use this to find the correct file names to load.

CHAPTER 2. CODE DOCUMENTATION

38

LOCKFILE is created when ArcPIC is started. The goal is to prevent accidental restarts of the executable overwriting previously produced simulation data. If the file is already present, ArcPIC will refuse to launch. It only contains a text with a short description of what it is used for.

2.13.2

Outputting using outputz.cpp, moms.cpp and vdf.cpp

All quantities in the output files are given in units of reference quantities (Te , nref , etc.). Therefore prior to outputting, variables used in the code have to be properly rescaled, which we discuss below. Outputting density We sum up the density (calculated as described in Sec. 2.11) on each grid point over av_time, and store it in the variable dens_av[j*NR+k]. The rescaled density in units of the reference density is then: (r, z) nREAL dens_av[j*NR+k] e , = sign nref f2 n_aver∆t

(2.183)

where sign is +1 for ions and −1 for electrons. The output files are called ne_XXXXXXXX.dat (for e- ), Cup_XXXXXXXX.dat (for Cu+ ), and Cu_XXXXXXXX.dat (for Cu), where XXXXXXXX is the time step number in the same format as listed in timeIndex.dat. The output is handled by the function out_dens_2D(...), and the file-names are stored in the global variables fn_e, fn_i and fn_n. Outputting coordinates We here output the position and velocity of all the particles in the system. This output is only active if the variable OUT_COORD is set to 0 in input.txt. The position is rescaled to units of λDb , such that rout =

r fr . = ∆z˜ λDb

(2.184)

−1 , and must thus be rescaled as Similarly, the velocities are given in units of λdb /ωpe

vout =

f v ∆z ˜, = v f dt_ion λdb ωpe ∆t

(2.185)

where dt_ion = 1 for electrons. The output files are called re_XXXXXXXX.dat (for e- ), rCup_XXXXXXXX.dat (for Cu+ ), and rCu_XXXXXXXX.dat (for Cu), where XXXXXXXX is the time step number in the same format as listed in timeIndex.dat. The output is handled by the function out_coords_2D(...), and the file-names are stored in the global variables fr_e, fr_i and fr_n. Outputting potential We here use directly the potential calculated in phi.cpp, averaged over nav_time steps. Therefore Eq. 2.28 still applies: f2 ϕREAL [V ] ∆z = ϕ˜av Tref [eV ] f2 n_aver ∆t

(2.186)

2.13. OUTPUT FILES

39

The output file is called phiXXXXXXXX.dat, where XXXXXXXX is the time step number in the same format as listed in timeIndex.dat. The output is handled by the function out_phi_2D(...), and the file-name stored in the global variable fphi. Outputting electric fields Similarly to the potential, the r- and z-components of the electric field are averaged over nav_time steps. The output is Eout =

f 2 ∆z f 2 n_aver ∆t

˜, E

(2.187)

and as Eq. 2.25 applies, the real field is E=

me 2 λDb ωpe Eout . e

(2.188)

The output files are called EzXXXXXXXX.dat (for the z-component) and ErXXXXXXXX.dat (for the r-component), where XXXXXXXX is the time step number in the same format as listed in timeIndex.dat. The output is handled by the function out_efield_2D(...), and the file-names are stored in the global variables fEz and fEr. This output file is only written if the variable OUT_EFIELD is set to 0 in input.txt. Outputting velocity In moms.cpp in the function aver_moments_2D(), velocities are summed over (i) averaging time steps and (ii) all particles in the given cell: mom[].ui = P pa[].p.vi . This is then rescaled to units of cs, Cu+ in out_velz_2D() from outputz.cpp by calculating mom[].ui /(u0 mom[].n), where (a) for electrons

s u0 = cs

M_ions[0] = cs, Cu+ , M_ions[1]

(2.189)

and (b) for ions s u0 = cs dt_ion

M_ions[0] = M_ions[1]

s

s Tref M_ions[0] dt_ion = cs, Cu+ dt_ion . mH+ M_ions[1] (2.190)

and the variable mom[].n is the number of super-particles in that cell. Finally we get exactly the average velocity in units of cs, Cu+ . For electrons, for instance (and similarly for ions) P mom[].ui hui i, n vi P = = . (2.191) cs, Cu+ u0 mom[].n cs, Cu+ i, n 1 P P Here the sums go over averaging time steps ( i ) and particles in the cell ( n ). The output files are called uezXXXXXXXX.dat (z-component), uerXXXXXXXX.dat (rcomponent) and uetXXXXXXXX.dat (t-component) for e- , and uizXXXXXXXX.dat (zcomponent) uirXXXXXXXX.dat (r-component) and uitXXXXXXXX.dat (t-component) for Cu+ . Here XXXXXXXX is the time step number in the same format as listed in timeIndex.dat. The file-names are stored in the global variables fv_ez, fv_er, fv_et, fv_iz, fv_ir and fv_it.

CHAPTER 2. CODE DOCUMENTATION

40

Outputting temperature Correspondingly, in aver_moments_2D(), temperatures are summed over (i) averaging time steps and (ii) all particles in the given cell: mom[].ti + = pa[].p.v2i . This is then rescaled to units of Te in outputz.cpp, out_tempz_2D(). Let’s assume a Maxwell-Boltzmann distribution9 ∝ e−

m(v−v0 )2 2T

(2.192)

in each of the directions separately, where, in general, there might be some stream velocity v0 in that direction. Then the temperature is nothing else but the standard deviation squared of this distribution:  #2  PN 2 " PN v i=1 vi  T = σ 2 = m(hv 2 i − hvi2 ) = m  i=1 i − (2.193) N N To output temperature in units of Te , we have to calculate: (i) for electrons: v2 T = 2 = Tref vt, e

P v2 Pi, n i i, n 1



hP

i vi 2 i, n 1

Pi, n

vt,2 e

so, in principle, T = Tref

mom_el[].t mom_el[].n





mom_el[].u mom_el[].n

,

(2.194)

.

(2.195)

2

vt,2 e

However, if out_velz_2D() has been called already before, the second term is already given in units of cs, Cu+ and we get T = Tref



mom_el[].t − (mom_el[].u)2 2 u0 mom_el[].n

with u0 = cs, Cu+ and fnorm = (ii) for ions:

 · fnorm ,

me mCu+ .

2 2 vion vion m + T = Cu = , Tref me (dt_ion)2 vt,2 e (dt_ion)2 c2s, Cu+

which implies T = Tref

(2.196)

mom_ion[].t mom_ion[].n



mom_ion[].u mom_ion[].n 2 (dt_ion) c2s, Cu+



(2.197)

2 .

(2.198)

However, in combination with out_velz_2D() we have to use the following expression in out_tempz_2D():   T mom_ion[].t 2 = − (mom_ion[].u) · fnorm , (2.199) Tref u20 mom_ion[].n with u0 = cs, Cu+ · dt_ion and fnorm = 1. 9

Note that discharge plasmas are usually non-thermal.

2.14. CONTROL ROUTINES, STABILITY OF PIC

41

The output files are called TezXXXXXXXX.dat (z-component), TerXXXXXXXX.dat (rcomponent) and TetXXXXXXXX.dat (t-component) for e- , and TizXXXXXXXX.dat (zcomponent) TirXXXXXXXX.dat (r-component) and TitXXXXXXXX.dat (t-component) for Cu+ . Here XXXXXXXX is the time step number in the same format as listed in timeIndex.dat. The file-names are stored in the global variables fT_ez, fT_er, fT_et, fT_iz, fT_ir and fT_it. Outputting velocity distribution If the variable OUT_VDF is set to 0 in input.txt, velocity distribution is outputted with a spatial resolution of 2 × 2 cells and a resolution of Nvdst equal parts in velocity space, centred around zero velocity. A given range of velocity can be sampled like this: For electrons, |v| ≤ Nvdst const vt, e , where we chose const = 200/3 and typically Nvdst = 401, resulting in a resolution of |v| ≤ 3vt, e . For ions, vt, e is replaced by cs, i . The velocity distributions are collected by the function vel_dst_along_2D(...), and the output is handled by out_fv_along_2D(...). The output files are called vdfDPXXXXXXXX.dat, where P is the particle (e, Cup or Cu), D the direction (z, r or abs) and XXXXXXXX is the time step number in the same format as listed in timeIndex.dat. The file names are stored in the global variables named fvdf_PD where P and D is the particle (this time e, i or n) and direction (z, r or abs). The velocity distributions are collected in the arrays vdf_PD (same naming scheme as for file name variables), which are only created on the heap if OUT_VDF is 0, as these arrays can be quite large.

2.14

Control routines, stability of PIC

To ensure the stability of the solution calculated with PIC, the smallest time- and length scale phenomena have to be well resolved, and therefore the following conditions have to be fulfilled all the time (cf. Sec. 1.1.6 and Sec. 2.2): −1 ∆t . 0.2 ωpe

(2.200)

∆z . 0.5 λDb

(2.201)

Since the computation happens in dimensionless quantities, the electron temperature and the −1 and λ electron density corresponding to the ωpe Db have to be initially guessed, these guessed values are the reference values Tref and nref , and iteratively adjusted until the guess is sufficiently close to reality. To have a feedback on whether the chosen scaling is suitable or not, a diagnostic routine is built in, which investigates the fulfilment of −1 ∆t = 0.2 ωpe

guessed

⇒√

1 nref

−1 . 0.2 ωpe 1 .√ ne

real

ne . nref

(2.202) (2.203) (2.204)

and similarly the fulfillment of ∆z = 0.5 (λDb )guessed . 0.5 (λDb )real r r Tref Te ⇒ . nref ne

(2.205) (2.206)

42

CHAPTER 2. CODE DOCUMENTATION

ne nref . Te Tref

(2.207)

In the code, a warning is written to the file error.log if ne > 5 nref in any cell, or if on average ne > nref in any row of cells between the cathode and the anode. Especially checked is that if on average ne > 5nref in the row of cells closest to the cathode, and breaking this stability requirement leads to program termination. ref Similarly, warnings are issued to the standard output if nTee > 2 nTref on average along any row of cells, or to error.log if it happens in any one single cell.

as

Note that here Te is here not defined as standard deviation of a Maxwellian distribution, but P 2 i |v|i P . (2.208) i1

This means that for instance the drift velocity is also taken into account (contrary to the temperature outputted, cf. Sec. 2.13.2), which is important when an external electric field is applied,.

3 Installing and running the 2D ArcPIC code, utilities and analysis system This chapter describes how to run the 2D ArcPIC code, and what functionality is provided by the accompanying tooling. The 2D ArcPIC code and the surrounding tooling is built to run for long periods of time on remote Linux machines, such as Hippu at CSC – IT center for Science Ltd. in Espoo, Finland. The first step for running the code is checking it out from the SVN repository and into your home directory. The newest version of the code is then contained in the subfolder pic2d/trunk, and code modifications should be done here so that they can easily be checked back into the SVN repository. This folder also contains the script makeRundir.py (Sec. 3.2), which installs the simulation program and the analysis scripts into a new folder (which is created). The install script also sets up the Makefile (Sec. 3.4) correctly, and creates log files which documents which SVN revision was installed, the output of svn diff etc. The RUNNAME referenced below is the name of the newly created folder. The user should then enter the newly created directory, and make any required changes to input.txt and dim.h here. The script calcScaling.py (Sec. 3.3) is useful for verifying that the parameters are correct. If necessary, configurations for using the correct compiler and external libraries should also be loaded. In case of long runs, it is recommended to use the terminal multiplexer screen, which allows the user to disconnect the SSH connection to the server without interrupting the simulation. The next step is then to compile the code using make and the Makefile which is in the src sub-directory (Sec. 3.4). This produces an executable ArcPIC_RUNNAME , which is stored in the root of the folder tree created by makeRundir. The simulation is normally started as ./ArcPIC_RUNNAME | tee readme.txt, which produces log file readme.txt containing all which has been written to the program’s standard output, as well as writing it directly to the terminal. This executable produces output to the out subfolder, as well as to the files which are saved to the same folder as the executable, i.e. circuit.dat, arcbounds.dat, mainStats.dat, LOCKFILE, error.dat and sometimes extra files from the currently loaded ArcBounds implementation such as arcbounds_original.dat. To download the output from the server to a local directory, the rsync utility is very useful. It also allows partial downloads, i.e. if the code is running, only download what has happened since last download. An example of rsync usage is given: rsync -rvvpzhltu -P machine :remote_path local_path --exclude="removed_e.dat"

This synchronizes the content of remote_path on machine to local_path on the current machine, while skipping files which has been altered on the local copy (typically analysis scripts), and skipping files named removed_e.dat1 . The --exclude option can use wildcards, as can also the remote_path . If remote_path does not contain a wildcard, it should typically be a folder name, omitting the final “/” on the end of the path to also get the folder and not just its contents. The local_path is typically just the current working directory “.”. If 2D ArcPIC has been started but the user wishes to start the run again (for example after realizing a mistake in input.txt), the target co in the Makefile deletes all output files, 1

These files are huge and usually not very useful – see Sec. 2.9.1 for more information

43

44

CHAPTER 3. RUNNING 2D ARCPIC, UTILITIES AND ANALYSIS SYSTEM

including LOCKFILE.

3.1

Supporting libraries

The 2D ArcPIC code uses two external libraries: SuperLU [12] and Gnu Scientific Library (GSL) [26]. Additionally, SuperLU depends on BLAS, which is typically provided by ATLAS [27] or by system libraries such as Intel Math Kernel Libraries (MKL). Of these, SuperLU and ATLAS are distributed together with ArcPIC, and are found in the folder support/pic2d/trunk in the SVN repository. Note that if compiling ATLAS from scratch, the machine must be completely unloaded and CPU throttling should be disabled. It is also possible to use a version of BLAS distributed together with SuperLU – this is often easier for debugging, but slower.

3.2

makeRundir.py

This Python script enables easy installation of the program into a new folder. It is run as follows: ./makeRundir.py new_path This creates the folder new_path (or exits with an error if it already exists), sets up the Makefile (Sec. 3.4), copies the analysis scripts (Sec. 3.5) and supporting files from the paths analysis/pic2d/trunk and support/analysis/pic2d/gle_template_redblue/trunk in the SVN working copy of the SVN repository into the folder GLEanalysis, sets the correct permission bits (chmod) and creates a log file detailing the current status of the SVN working copy.

3.3

calcScaling.py

This script parses input.txt, and calcuates the system size in physical units, gridsizes ∆z, particle weighting Nsp , field strengths etc. This is useful for checking that the system is specified correctly, and also provides scaling factors necessary for interpreting the output data and generating plots with physical units. It is therefore provided with two interfaces – it can either be ran as a stand-alone program on the command line, or it can be imported as a Python module providing a class for parsing input.txt. The user of the code will typically use the first interface, while the analysis scripts will use the second.

3.4

Makefile

The Makefile, which is found in the src-folder, is used to build the program. It contains the following targets: d runs gccmakedep, which checks dependencies for compilation and sets up the Makefile. This results in a large number of lines being written below the line #DO NOT DELETE in the Makefile. e compiles and links the main executable, ArcPIC_RUNNAME .

3.5. ANALYSES

45

tests This target compiles and links the test executables. c cleans all object files (but not executables). co deletes all output files. The normal procedure after installing 2D ArcPIC to a folder and having set up input.txt and dim.h is to run the targets c, d and e in that order. Inside the Makefile, the variable RUNNAME contains the name of the install folder, and WORK it’s complete path. If the folder is renamed or moved, these variables must be updated. These are normally set by makeRundir.py. There are also several different version of the variables specifying the supporting libraries and compiler flags. If debugging or running on a different type of machine where the operating system provides high-performance math libraries (such as MKL which is provided on hippu.csc.fi), it is often useful to modify these before compiling.

3.5

Analyses

Many different analysis scripts are distributed with the 2D ArcPIC code. These are written in Python or GLE [28] + Unix shell. Describing each these analyses are outside the scope of this manual, as they change more often than the main code. They are therefore considered to be documented by their own source code (incl. comments). As described in Sec. 3.2, the analysis scripts are installed to the folder GLEanalysis by makeRundir.py. These scripts reads the output files from the 2D ArcPIC code, does varying amounts of processing on the output, and produces plots and derived values. In general, the GLE analyses are considered to be depreciated and are not maintained. The actively developed analysis are thus the Python scripts. As described in Sec. 3.3, the analyses load calcScaling.py as a module in order to calculate the conversion factors from the output units to physical units.

4 Recommended directions for future development of the 2D ArcPIC code The 2D ArcPIC code has recently undergone several structural changes, factoring the logic into separate silos for the different pieces of the model, thus making the program more modular. This has enabled the development of different models for the surface (Sec. 2.9) and external circuit (Sec. 2.8). There are, however, still many things which are worth improving in the code itself, on top of improvement and experimentation with the physics models. This chapter lists some of these issues, for the benefit of future maintainers. Improved particle manager One basic piece of infrastructure which could benefit from rewriting is the particle manager. This could greatly improve the performance and flexibility of the code, while making it easier to understand and creating a solid platform for future developments. A modular approach where each class of particles are handled by a “particle list” object could simplify the code by eliminating several global variables used for storing the properties of the particle types, and make it easier to add more particle types or states. There are also possibilities for performance improvements by reducing the amount of data copying now done in the particle ordering algorithm, making it simpler and faster. A cleanup / rewrite of this code would touch most parts of the code, so it would be useful to plan it together with other changes, like the ones listed below. Dynamical particle weighting Today, the particle/superparticle ratio Nsp is the same for all particles throughout the whole simulation, even if the densities changes over many orders of magnitude in the course of a typical simulation. This leads to poor resolution of the particle distributions in phase-space in the initial stages of the breakdown, and a very large number of particles in the final stages. A possible solution to this problem is the dynamic particle weighting technique [29, 30], which works by using smaller particle weights at low densities, and higher weights when the density increases, combining or splitting super-particles according to changes in density. keeping the number of particles in a given part of the phase-space roughly constant. One drawback of this method is that it does not conserve both energy and momentum exactly when merging particles. Modular collisions handling The current collision system implemented in collisions.cpp replicates the same or very similar code in multiple functions. Unifying these functions and making them parallel (using for example a technique similar to the one described in Sec. 2.7.3) could significantly speed up the code and make it more transparent. Another benefit would be to make it easier to setup and configure new types of collisions, especially if adding more particle types or states to the code. Binary file outputs Currently, all output files from the code are ASCII tables, which makes them very large (on the order of several 100 GB for a long simulation), and also taking more time than necessary to read for the analyses. Converting to a binary format such as

47

48

CHAPTER 4. RECOMMENDED DIRECTIONS FOR FUTURE DEVELOPMENT

NetCDF [31] or HDF5 [32] for at least the larger output files would help these issues, at the cost of changing some parts of the output and analysis routines. Checkpointing Another avenue for improvement would be to re-activate the checkpointing procedure. It is likely beneficial to do this in conjunction with implementing binary file outputs, as a binary output file would ensure that exactly the same values are restored in memory as was written out. Faster calculation of Igap The Shockley-Ramo theorem is provides an effective way to calculate the induced charge on the electrodes without needing to run a field solver. However, 2D ArcPIC does run a field solver that includes the charges in the gap. This field can be directly used to calculate the induced charge on the cathode Qind.cat. due to the charges in the gap as Z R Ez r dr − Q0ind.cat. , (4.1) Qind.cat. = 0 2π 0

Q0ind.cat.

is the charge induced by the external electric field, i.e. the capacitance of where the simulation domain. This is given as Q0ind.cat. = 0 2πR2

U . Z

(4.2)

A List of files in the source code input.txt This is the standard version of the input file, which sets the system parameters. Note that if the number of mesh cells (nr and nz) are changed, the user must also update h/dim.h. An example of this file is shown in Appendix B. It is recommended to not change the parameters in this copy of input.txt before every run, but rather first install to a new folder with makeRundir.py and then change input.txt and h/dim.h there. makeRundir.py This Python script installs the sources and analysis scripts to a given path, and then sets up the Makefile. This is done so that several versions of the code can be installed at different locations and ran in parallel with different input.txt, without touching the “development” source folder controlled by SVN. See Sec. 3.2 for more information. calcScaling.py This Python script parses input.txt and converts the values from dimensionless units to normal units. It is very useful for checking the input file. It is also used by many of the analyses for converting the values in the output files into normal units. See Sec. 3.3 for more information. h This folder contains all the .h-files, describing the functions and classes defined in the .cpp-files. Here, only the header files which are implementing functionality and are thus not only listing function signatures in the corresponding .cpp files, are listed. dim.h contains #defined constants setting the system size, maximum number of superparticles in the system etc. Note that if the number of superparticles are set too high, the size of the particle arrays can exceed the capacity of the program’s data area, leading to compiler errors unless the capacity is increased (via compiler options). arrays1.h defines all global arrays used by the program for such purposes as holding the particles. var.h defines a set of global variables, which hold data necessary for scaling to/from dimensionless units etc. It also defines the pointers used to reference the Circuit, ArcBounds and InitialParticles object currently in use, as well as configuration for the Random Number Generator (RNG), parallelization, and flags indication which output files to produce. pic.h defines multiple data types used throughout the program, such as Particle and PhasespaceVar. mydef.h #defines useful constants such as PI. Also defined are macros for finding the maximum of two numbers and squaring a number. outp.h contains global variables used to hold the filenames for the next set of output files, as generated by filenames.cpp. input.h defines the functions in input.cpp and #defines some configurations for readInputSection(...).

49

50

APPENDIX A. LIST OF FILES IN THE SOURCE CODE

random.h defines the functions in random.cpp, which is a wrapper for the GNU Scientific Library (GSL) [26] random number generator. It also #defines a constant which sets which of the RNGs distributed with GSL to use, as well as the RAND macro which expands to a function call to Random(), which returns a uniformly distributed random number between 0 and 1. arcbounds.h defines the abstract class ArcBounds. The classes inheriting from this implement surface models and particle injection algorithms, as described in Sec. 2.9. circuit.h defines the abstract class Circuit. The classes inheriting from this control the potential difference between anode and the cathode as a function of current and time, as described in Sec. 2.8. initialParticles.h defines the abstract class InitialParticles. The classes inheriting from this can be used to pre-seed the volume with particles sampled from some distribution, as described in Sec. 2.10. colls.h defines the functions in collisions.cpp. It also #defines the constant COLL_N_N_2D_OMP_MINPARTICLES, which configures the parallelization as described in Sec. 2.7.3. src This folder contains all the .cpp files which are compiled to create the 2D ArcPIC executable. Only the most important files are listed here. main.cpp contains the main() loop, as illustrated in Fig. 1.1. It is also responsible for initializing the code, as well as writing of the LOCKFILE, out/timeIndex.dat and mainStats.dat output files (Sec. 2.13.1). init.cpp initializes the main global variables and collision cross-section data. input.cpp parses the input file input.txt and creates the output folder if it doesn’t exist. Some of the input file is parsed simply by reading lines with fscanf(...), while the Arcbounds, Circuit and InitialParticles classes are partially responsible for parsing their own input data. These classes use the readInputSection(...) method defined in input.cpp to do this arsing. print_par.cpp contains the function print_parameters_2D(), which prints the parameters read in from the input file during initialization, similarly to calcScaling.py. arcbounds.cpp defines the common services and default function implementations for ArcBounds classes, as described in Sec. 2.9. These classes implement surface physics models, and which class to load can be selected by the input file. Some classes inheriting from and implementing all or part of the virtual methods in ArcBounds are also defined in this file, such as ArcRemover, ArcDummy and ArcSimple. arcbound_original.cpp contains the same surface physics model as used in 2D ArcPIC v2.07 [1] with PIS and PBC both set to 0, except fixing a bug which disabled flat surface field emission in the older version. arcbound_original_newheatspike.cpp contains an updated version of the original model described above, using a new model for enhanced sputtering at high ion bombardment fluxes.

51

circuit.cpp defines the common services and default function implementations for Circuit classes, as described in Sec. 2.8. The file also contains all the classes which implements the abstract class Circuit. initialParticles.cpp defines the common services and default function implementations for InitialParticles classes, as described in Sec. 2.10. It also contains the class UniformRestricted, which is currently the only implementation of abstract class InitialParticles. collisions.cpp implements the collision routines described in Sec. 2.7. order.cpp implements a particle sorter, which sorts the particle arrays according to cell, as required by the collision routines. random.cpp implements a wrapper around GSL’s random number generator. It is capable of generating separate random number streams for parallel threads, resulting in a thread-safe running and repeatable result. The RNGs are initialized using a single seed, which during program initialization is used to seed a secondary generator which generates the seeds for the primary generators used during program execution. Functions for generating uniform- and normal-distributed random numbers are provided. phi.cpp implements the solver for the electric potential ϕ, ˜ as described in Sec. 2.2.4 and 2.4, using SuperLU [12]. It is split in two parts – one which factorizes the fieldsolver matrix using a LU-decomposition method, and is called during initialization of the program. The other part is called at every time step in order to update the potential, using the current charge distribution and anode/cathode potentials. efield.cpp finds the electric field at every grid point by calculating the derivative of the electric potential. See Sec. 2.3 and 2.4 for more information. e_ion.cpp contains functions for re-scaling the electric field using Eq. (2.20) for use by the particle pusher when acting on ions. push.cpp implements the particle pusher, as described in Sec. 2.6 and 1.1.3. outputz.cpp contains several functions for writing output files, as described in Sec. 2.13.2. moms.cpp contains helper functions for the output file writing functions in outputz.cpp. density.cpp contains the function density_2D(...), which calculates the particle density as described in Sec. 2.11. This is used for output to files, for the field solver, and for checking stability. vdf.cpp defines the function vel_dst_along_2D(...), which is used for output of the average and standard deviation (temperature) of the velocity distributions as explained in Sec. 2.13.2. If the flag OUT_VDF is 0 (controlled by input.txt), create velocity distributions to be written to file by functions in outputz.cpp. This is normally disabled, as generating these distributions requires large amounts of memory. filenames.cpp creates the filenames for the output files, which includes a time code. The filenames are stored in the global variables defined in outp.h.

52

APPENDIX A. LIST OF FILES IN THE SOURCE CODE

engy.cpp contains functions for calculating the kinetic and potential energy of the particles in the plasma, as discussed in Sec. 2.12. colltest.cpp compiles to an executable colltest, which is a test routine for checking that the collisions work properly. rngtest.cpp compiles to an executable rngtest, which is a test routine for checking that the RNG code works correctly. rngparatest.cpp compiles to an executable rngparatest, which is a test routine for checking that the parallelization of the RNG works correctly. picFieldTest.cpp compiles to an executable picFieldTest, which is a test routine calculating the error from estimating the field at r = z = 0 (position of “tip” emitter in many emission models) using PIC and the field solver instead of analytically. Makefile enables compilation of the code. See Sec. 3.4 for more information. inputdata folder contains data files which are read by some models for initialization. tests This folder contains Python scripts for analyzing the data from the test executables colltest, rngtest and picFieldTest. These executables are also compiled into this folder, with the exception of picFieldTest, which is compiled into the root folder as it is using some of the standard input and output routines. colltestAna.py This Python script parses the output from the colltest executable and produces plots. plotPicFieldTest.py This Python script parses the output from the picFieldTest executable and produces plots. rngtestAna.py This Python script parses the output from the rngtest executable and produces plots.

B Example input.txt 1 ∗∗∗ INPUT FILE FOR 2D ARC−PIC CODE ∗∗∗ 2 ∗∗∗ ( C ) HELGA TIMKO, 2010 ∗∗∗ 3 ∗∗∗ CERN and H e l s i n k i I n s t i t u t e o f P h y s i c s ∗∗∗ 4 ∗∗∗ ( C ) KYRRE SJOBAK , 2012−2013 ∗∗∗ 5 ∗∗∗ CERN and U n i v e r s i t y o f O s l o ∗∗∗ 6 ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗ 7 8 ∗ SCALING PARAMETERS ∗ 9 10 R e f e r e n c e d e n s i t y i n 1 / cm3 : 1 6 . 0 e18 11 R e f e r e n c e t e m p e r a t u r e i n eV : 2900. 12 P a r t i c l e s i n a Debye c u b e Ndb : 1500. 13 14 Number o f c e l l s nr , nz : 480 120 15 G r i d s i z e dz i n Debyes : .5 16 T i m e s t e p d t i n Omega_pe−s : .2 17 18 19 ∗ TIMESTEPS ∗ 20 21 I o n t i m e s t e p d t _ i o n i n d t−s : 5 22 C o l l i s i o n t i m e s t e p n c o l l _ e l i n d t : 5 23 C o l l i s i o n t i m e s t e p n c o l l _ i o n i n d t : 5 24 S t a b i l i t y d i a g n o s t i c s i n d t : 1000 25 26 S i m u l a t i o n t i m e i n O_pe−s : 3000000 27 O u t p u t t i n g t i m e d t _ o u t i n O_pe−s : 1000. 28 S t a r t o f a v e r a g i n g i n O_pe−s : 1. 29 D u r a t i o n o f a v e r a g i n g i n O_pe−s : 1. 30 31 32 ∗ FIELDS , PARTICLES AND BOUNDARY CONDITIONS ∗ 33 34 E x t e r n a l d i m e n s i o n l e s s B_z f i e l d : 0. 35 E x t e r n a l d i m e n s i o n l e s s B_t f i e l d : 0. 36 37 − I n j e c t i o n t i m e s t e p s − 38 I n j e c t i o n t i m e s t e p i n d t , e l e c t r o n s : 1 39 I n j e c t i o n t i m e s t e p i n d t , n e u t r a l s : 5 40 I n j e c t i o n t i m e s t e p i n d t , i o n s : 10 41 42 ∗ P a r t i c l e b o u n d a r y p a r a m e t e r s ∗ 43 ArcBoundName : ArcOriginalNewHS 44 ∗∗∗ ArcDummy 45 / / / 46 ∗∗∗ A r c S i m p l e 47 n e _ i n j e c t : 10 Number o f e l e c t r o n s t o i n j e c t / t i m e s t e p 48 n i _ i n j e c t : 10 Number o f i o n s ( t y p e 1 , Cu + ) t o i n j e c t / t i m e s t e p 49 n n _ i n j e c t : 10 Number o f n e u t r a l s t o i n j e c t / t i m e s t e p 50 f i l e _ t i m e s t e p : 1 How o f t e n t o w r i t e t o ’ a r c b o u n d s . d a t ’ [ dt ] 51 f i l e _ t i m e s t e p _ h i s t : 0 How o f t e n t o w r i t e c u r r e n t h i s t o g r a m s ( 0 = n o r m a l o u t p u t s t e p s ) 52 / / / 53 ∗∗∗ A r c O r i g i n a l 54 b e t a : 35. I n i t i a l Fowler−Nordheim b e t a o f t i p 55 b e t a _ f : 2. Fowler−Nordheim b e t a o f f l a t s u r f a c e 56 j _ m e l t : 5 0 . e8 C u r r e n t d e n s i t y where b e t a i s s e t e q u a l t o b e t a _ f [A / cm ^ 2 ] 57 R e m i s s i o n : 4. R a d i u s o f i n j e c t i o n o f FN e m i t t e d e− from t i p [ dz ] 58 R e m i s s i o n _ t h e o r : 1 . Radius of t i p , used f o r c u r r e n t c a l c [ dz ] 59 SEY : 0.5 Secondary e l e c t r o n y i e l d 60 r_Cu_e : 0 . 0 1 5 R a t i o o f num . Cu e v a p o r a t e d t o e l e c t r o n s FN−e m i t t e d 61 t i p P i : n I n c l u d e t h e m i s s i n g PI i n t h e t i p a r e a ? [y/n] 62 f r a c I n j e c t S t e p : n Inject electrons at fractional timestep [y/n] 63 a l p h a _ f l a t : 1.0 F r a c t i o n a l area of f l a t s u r f a c e e m i t t i n g 64 h e a t s p i k e _ t h r e s : 1 e7 T h r e s h o l d f o r i o n c u r r e n t t h r o u g h a c e l l t o t r i g g e r h e a t s p i k e [A / cm ^ 2 ] 65 h e a t s p i k e _ y i e l d : 2 1 3 5 2 . 1 7 8 1 0 6 6 9 9 9 9 7 Number o f i n j e c t e d Cu n e u t r a l s / s t a n d a r d t i m e s t e p ( ~ 1 . 7 7 f s ) 66 67 68 69 70 71 72 73

i n an h e a t s p i k e . D e f a u l t v a l u e = 1000 s u p e r p a r t i c l e s file_timestep : 1 How o f t e n t o w r i t e t o ’ a r c b o u n d s . d a t ’ /// ∗∗∗ ArcOriginalNewHS beta : 35. I n i t i a l Fowler−Nordheim b e t a o f t i p beta_f : 2. Fowler−Nordheim b e t a o f f l a t s u r f a c e j_melt : 5 0 . e8 C u r r e n t d e n s i t y where b e t a i s s e t e q u a l t o b e t a _ f do_erosion : n Does n e u t r a l e v a p o r a t i o n e r o d e t h e t i p , l o w e r i n g Remission : 8. R a d i u s o f i n j e c t i o n o f FN e m i t t e d e− from t i p

53

[ dt ]

[A / cm ^ 2 ] beta ? [y / n] [ dz ]

54

APPENDIX B. EXAMPLE INPUT.TXT

74 R e m i s s i o n _ t h e o r : 1 . 1 2 8 3 7 9 1 6 7 0 9 5 5 1 2 6 R a d i u s o f t i p , u s e d f o r c u r r e n t c a l c [ dz ] 75 R b o r d e r : 1 . 1 2 8 3 7 9 1 6 7 0 9 5 5 1 2 6 Where t o s t a r t f l a t s u r f a c e e m i s s i o n [ dz ] 76 SEY : 0.5 Secondary e l e c t r o n y i e l d 77 r_Cu_e : 0 . 0 1 5 R a t i o o f num . Cu e v a p o r a t e d t o e l e c t r o n s FN−e m i t t e d 78 r _ C u _ e _ f l a t : 0 . 0 1 5 R a t i o o f num . Cu e v a p o r a t e d t o e l e c t r o n s FN−e m i t t e d from f l a t s u r f a c e 79 e v a p _ f l a t _ c e n t e r : n E v a p o r a t e c o p p e r from f l a t s u r f a c e i n s i d e R e m i s s i o n ? [ y / n ] 80 f r a c I n j e c t S t e p : y Inject electrons at fractional timestep [y/n] 81 a l p h a _ f l a t : 1.0 F r a c t i o n a l area of f l a t s u r f a c e e m i t t i n g 82 h e a t s p i k e _ t h r e s : 1 e25 T h r e s h o l d f o r i o n c u r r e n t t h r o u g h a c e l l t o t r i g g e r h e a t s p i k e [ p a r t i c l e s / 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151

cm ^ 2 / s ] heatspike_YR : 1 Y i e l d r a t i o p = ( Number o f i n j e c t e d Cu n e u t r a l s ) / ( number o f i m p a c t o r s ) in a heatspike heatspike_model : c How t o d e t e r m i n e number o f i m p a c t o r s and t h e r a d i u s [ a / b / c ] file_timestep : 1 How o f t e n t o w r i t e t o ’ a r c b o u n d s . d a t ’ [ dt ] /// ∗∗∗ F l e x F N _ r i n g alpha_ring : 1. F r a c t i o n a l area of e m i t t i n g ring e m i t t i n g beta_ring : 35 Fowler−Nordheim b e t a o f t h e s e e m i t t e r s idx_ring1 : 2 Lowest mesh node i d x ( i n r ) e m i t t i n g idx_ring2 : 2 H i g h e s t mesh node i d x ( i n r ) e m i t t i n g file_timestep : 1 How o f t e n t o w r i t e t o ’ a r c b o u n d s . d a t ’ [ dt ] /// ∗∗∗ FlexFN_twoComp alpha1 : 0.1 Alpha f o r component 1 alpha2 : 1.0 Alpha f o r component 2 beta1 : 3 5 . 0 B e t a f o r component 1 beta2 : 1.0 B e t a f o r component 2 idx1 : 0 H i g h e s t mesh node i n d e x e m i t t i n g [ dr ] file_timestep : 1 How o f t e n t o w r i t e t o ’ a r c b o u n d s . d a t ’ [ dt ] /// ∗∗∗ F l e x F N _ t w o C o m p _ o r i g N e u t r a l s alpha1 : 0.1 Alpha f o r component 1 alpha2 : 1.0 Alpha f o r component 2 beta1 : 3 5 . 0 B e t a f o r component 1 beta2 : 1.0 B e t a f o r component 2 idx1 : 0 H i g h e s t mesh node i n d e x e m i t t i n g [ dr ] SEY : 0.5 Secondary e l e c t r o n y i e l d r_Cu_e : 0 . 0 1 5 R a t i o o f num . Cu e v a p o r a t e d t o e l e c t r o n s FN−e m i t t e d doHeatspike : y Enable / d i s a b l e h e a t s p i k e s p u t t e r i n g [y/n] file_timestep : 1 How o f t e n t o w r i t e t o ’ a r c b o u n d s . d a t ’ [ dt ] /// ∗∗∗ F l e x F N 2 _ r i n g ring_alpha : 1. F r a c t i o n a l area of e m i t t i n g ring e m i t t i n g ring_beta : 35 Fowler−Nordheim b e t a o f t h e s e e m i t t e r s ring_r1 : 0.0 Lowest p o i n t ( i n r ) e m i t t i n g [ dz ] ring_r2 : 2.0 Highest point ( in r ) emitting [ dz ] file_timestep : 1 How o f t e n t o w r i t e t o ’ a r c b o u n d s . d a t ’ [ dt ] /// ∗∗∗ F l e x F N 2 _ t a b l e N e u t r a l s emitter_radius : 1. R a d i u s o f e m i t t e r [ dz ] emitter_alpha : 1.0 F r a c t i o n a l area of e m i t t i n g area e m i t t i n g emitter_beta : 35 Fowler−Nordheim b e t a o f e m i t t e r flat_alpha : 1.0 F r a c t i o n a l area of f l a t s u r f a c e e m i t t i n g flat_beta : 1.0 Fowler−Nordheim b e t a o f f l a t s u r f a c e file_timestep : 1 How o f t e n t o w r i t e t o ’ a r c b o u n d s . d a t ’ [ dt ] /// ∗∗∗ SumTipFNField beta_tip : 35. Fowler−Nordheim b e t a o f t i p Remission : 4. R a d i u s o f i n j e c t i o n o f FN e m i t t e d e− from t i p [ dz ] R e m i s s i o n _ t h e o r : 0 . 5 6 4 1 8 9 5 8 3 5 4 7 7 5 6 2 8 R a d i u s o f t i p , u s e d f o r c u r r e n t c a l c [ dz ] file_timestep : 1 How o f t e n t o w r i t e t o ’ a r c b o u n d s . d a t ’ [ dt ] /// ∗∗∗ END ∗ Circuit parameters ∗ CircuitName : F i x e d V o l t a g e _ r e s i s t o r C a p a c i t o r ∗∗∗ F i x e d V o l t a g e U0 : 0. Ground p o t e n t i a l Zmin [ Te ] UNz : 2. P o t e n t i a l a t Zmax [ Te ] file_timestep : 1 How o f t e n t o w r i t e t o ’ c i r c u i t . d a t ’ /// ∗∗∗ F i x e d V o l t a g e _ r e s i s t o r UNz : 0.6 S u p p l y v o l t a g e & i n i t i a l p o t e n t i a l a t Zmax [ Te ] Rseries : 100.0 External s e r i e s r e s i s t a n c e [Ohm] file_timestep : 1 How o f t e n t o w r i t e t o ’ c i r c u i t . d a t ’ [ dt ] /// ∗∗∗ F i x e d V o l t a g e _ r e s i s t o r C a p a c i t o r UNz : 0.6 S u p p l y v o l t a g e & i n i t i a l p o t e n t i a l a t Zmax [ Te ] Rseries : 1000.0 External s e r i e s r e s i s t a n c e [Ohm] Cgap : 1 e−12 Gap c a p a c i t a n c e [F]

55

152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221

file_timestep : 1 How o f t e n t o w r i t e t o ’ c i r c u i t . d a t ’ [ dt ] /// ∗∗∗ F i x e d V o l t a g e _ r a m p U0 : 0. Ground p o t e n t i a l Zmin [ Te ] UNz : 2. P o t e n t i a l a t Zmax [ Te ] rampTime : 0 . 5 e−3 Time t o r e a c h f u l l p o t . [ n s ] file_timestep : 1 How o f t e n t o w r i t e t o ’ c i r c u i t . d a t ’ /// ∗∗∗ T w o C a p a c i t o r s U0 : 0. Ground p o t e n t i a l Zmin [ Te ] UNz : 2. P o t e n t i a l a t Zmax [ Te ] C_ext : 0 . 5 e−12 Gap c a p a c i t a n c e [F] C_ext2 : 1 . e−9 External capacitance [F] R_ext : 0. External resistance [Ohm] t_change : 50.0 When t o c h a n g e from i n t e r n a l t o e x t e r n a l c a p a c i t a n c e [ n s ] avoid_negativ : n Lock UNz > U0? [y/n] file_timestep : 1 How o f t e n t o w r i t e t o ’ c i r c u i t . d a t ’ /// ∗∗∗ F r e q u e n c y E x c i t a t i o n FreqGHz : 12.0 E x c i t a t i o n o s c i l l a t i o n f r e q u e n c y [ GHz ] UNz_ampl : 2.0 Amplitude of e x c i t a t i o n file_timestep : 1 How o f t e n t o w r i t e t o ’ c i r c u i t . d a t ’ /// ∗∗∗ END o f c i r c u i t p a r a m e t e r s ∗ I n i t i a l p a r t i c l e d i s t r i b u t i o n ( ’ None ’ i s an I n i t i a l P a r t i c l e T y p e : None ∗∗∗ U n i f o r m R e s t r i c t e d density : 1 e18 Particle density maxR : 5.0 Area t o f i l l maxZ : 5.0 −−− " " −−− minR : 0.0 −−− " " −−− minZ : 0.0 −−− " " −−− doInject_e : n Inject electrons ? doInject_i : n Inject ions ? doInject_n : y Inject neutrals ? Tinj : 300 Injection temperature /// ∗∗∗ END o f i n i t i a l p a r t i c l e d i s t r i b u t i o n

accepted option ) ∗

[ p a r t i c l e s / cm ^ 3 ] [ dz ]

[y/n] [y/n] [y/n] [ eV ]

∗ OPTIONS ∗ ( 0= y e s 1= no ) Outputting p a r t i c l e coordinates : Outputting e l e c t r i c field : O u t p u t t i n g VDF : M a g n e t i c f i e l d on :

0 0 1 1

( 0 = new run , o t h e r w i s e u s e n r o f p r e v i o u s r u n + 1 ) C o n t i n u e from which r u n : 0 ∗ MISCELLANEOUS ∗ Ti_over_Te : mi_over_me : Seed f o r random number g e n e r a t o r : F i e l d b o u n d a r y c o n d i t i o n BC : Number o f CPUs t o u s e :

.005 1836.15 1234 4 2

∗∗∗ INFO : ∗∗∗ I n c a s e o f r e−s t a r t , o n l y Ndb , d t _ o u t , R e m i s s i o n , R e m i s s i o n _ t h e o r , and CONTINUATION a r e a c t i v e . Remember t o u p d a t e t h e no . o f r u n . ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗

C Miscellaneous useful formulas Debye length r λDb =

ε0 Te = e 2 ne

r

552635 × Te [cm] ne

(C.1)

Plasma frequency s ωpe =

  √ 1 e2 ne = 56414.6 ne ε0 m e s

(C.2)

Gyroradius (Larmor radius)

me v⊥ (C.3) eB where B is the magnitude of a constant external magnetic field and v⊥ is the speed of the electron perpendicular to the magnetic field. rg =

Gyrofrequency νg =

eB 2πme

(C.4)

Thermal velocity r vt, e = λDb ωpe =

Te me

(C.5)

Sound velocity r cs, i =

me vt, e Mi

(C.6)

Number of superparticles nref λ3Db , (C.7) NDb where NDb is defined as the number of superparticles in a Debye cube corresponding to the density nref . NSP =

Ionisation rate Using uniform and monoenergetic electron and neutral beams in a head-on collision provides an easy way to cross-check the ionisation routine. The theoretically predicted ionisation rate is given by R = ne · nCu · vrel · σ(E),

(C.8)

where ne and nCu are the electron and neutral number densities, respectively, vrel is the relative speed of the beams and the total cross-section σ is calculated at the energy E = 1 1 2 2 m∗ vrel , with m∗ ≈ me being the reduced mass. The unit of R is [R] = cm3 ·s (or 1 equally well [R] = m3 ·s ). This rate is then to be compared with the simulated rate NSP · Γsim Γreal = , (C.9) V · ∆t V · ∆t where Γsim is the number of simulated ions yielded in the collisions during a timestep of ∆t in a volume of V , and where electrons and neutrals were initially uniformly distributed, e.g. V = nz π (2j + 1)∆r2 ∆z for particles distributed in the inner j cells. R=

57

APPENDIX C. MISCELLANEOUS USEFUL FORMULAS

58

Ion species Sorts[ ]

0 1 2

H+ (reference species) Cu+ Cu

Masses M_ions[0] = MH+ = 1836.15 me

(C.10)

M_ions[1] = MCu+ = (1836.15 × 63.546/1.0073 − 1)me

(C.11)

M_ions[2] = MCu = (1836.15 × 63.546/1.0073)me

(C.12)

The numerical factor 1836.15 is stored by the variable mi_over_me, which is set in input.txt.

Bibliography [1]

H. Timko, 2D Arc-PIC code description: methods and documentation, CLIC Note 872, CERN, 2011.

[2]

O. Buneman, Dissipation of Currents in Ionized Media, Phys. Rev. 115 (1959) 503.

[3]

J. Dawson, One-dimensional Plasma Model, Phys. Fluids 5 (1962) 445.

[4]

R. Hockney and J.W.Eastwood, Computer Simulation Using Particles, IOP Publishing Ltd, Bristol and Philadelphia, fifth edition, 1999, Chaps. 1.1-1.5, 2.1-2.4, 9.1-9.3.

[5]

C. Birdsall and A. Langdon, Plasma Physics via Computer Simulation, McGraw-Hill, Inc., United States of America, first edition, 1985.

[6]

J. Boris, Relativistic plasma simulation-optimization of a hybrid code, in Proceedings of the 4th Conference on Numerical Simulation of Plasmas, pages 3–67, Naval Res. Lab., Washington, D.C., 1970.

[7]

D. Tskhakaya, Chapter 6: The Particle-in-Cell Method, in Computational Many-Particle Physics, edited by H. Fehske, R. Schneider, and A. Weiße, pages 161–189, SpringerVerlag, Berlin Heidelberg, first edition, 2008.

[8]

T. Takizuka and H. Abe, A binary collision model for plasma simulation with a particle code, Journal of Computational Physics 25 (1977) 205.

[9]

C. Birdsall, Particle-in-cell charged-particle simulations, plus Monte Carlo collisions with neutral atoms, PIC-MCC, IEEE Transactions on Plasma Science 19 (1991) "65.

[10] V. Vahedi and M. Surendra, A Monte Carlo collision model for the particle-in-cell method: applications to argon and oxygen discharges, Computer Physics Communications 87 (1995) 179, Particle Simulation Methods. [11] D. Tskhakaya, K. Matyash, R. Schneider, and F. Taccogna, The Particle-In-Cell Method, Contrib. Plasm. Phys. 47 (2007) 563. [12] SuperLU general purpose library, http://crd.lbl.gov/~xiaoye/SuperLU/, project funded by DOE, NSF and DARPA. [13] N.Shipman, S. Calatroni, and W. Wuensch, MEASUREMENT OF THE DYNAMIC RESPONSE OF THE CERN DC SPARK SYSTEM AND PRELIMINARY ESTIMATES OF THE BREAKDOWN TURN-ON TIME, in Proceedings of IPAC2012, New Orleans, Louisiana, USA, 2012. [14] J. P. Verboncoeur, Symmetric Spline Weighting for Charge and Current Density in Particle Simulation, Journal of Computational Physics 174 (2001) 421. [15] K. Matyash, Kinetic Modeling of Multi-Component Edge Plasmas, PhD thesis, ErnstMoritz-Arndt-Universität Greifswald, 2003.

59

60

BIBLIOGRAPHY

[16] OpenMP Architecture Review Board, OpenMP Application Program Interface Version 3.0, 2008. [17] W. Shockley, Currents to Conductors Induced by a Moving Point Charge, Journal of Applied Physics 9 (1938) 635. [18] S. Ramo, Currents Induced by Electron Motion, Proceedings of the IRE 27 (1939) 584. [19] R. H. Fowler and L. Nordheim, Electron Emission in Intense Electric Fields, Royal Society of London Proceedings Series A 119 (1928) 173. [20] J. W. Wang and G. A. Loew, Field emission and rf breakdown in copper linac structures, Part. Accel. 30 (1989) 225. [21] S. Parviainen, F. Djurabekova, H. Timko, and K. Nordlund, Electronic processes in molecular dynamics simulations of nanoscale metal tips under electric fields, Comp. Mater. Sci. 50 (2011) 2075 . [22] R. G. Forbes, Field evaporation theory: a review of basic ideas, Appl. Surf. Sci. 87-88 (1995) 1 , Proceedings of the 41st International Field Emission Symposium. [23] D. R. Kingham, Model calculations of tunnelling and thermal evaporation rate constants relating to field evaporation, J. Phys. D: Appl. Phys. 15 (1982) 2537. [24] Y. Yamamura and H. Tawara, Energy dependence of ion-induced sputtering yields from monatomic solids at normal incidence, At. Data Nucl. Data 62 (March 1996) 149. [25] H. Timko, F. Djurabekova, K. Nordlund, L. Costelle, K. Matyash, R. Schneider, A. Toerklep, G. Arnau-Izquierdo, A. Descoeudres, S. Calatroni, M. Taborelli, and W. Wuensch, Mechanism of surface modification in the plasma-surface interaction in electrical arcs, Phys. Rev. B 81 (2010) 184109. [26] M. Galassi, J. Davies, J. Theiler, B. Gough, G. Jungman, P. Alken, M. Booth, and F. Rossi, GNU Scientific Library Reference Manual, Network Theory Ltd., 3rd edition, 2009. [27] See homepage for details, Automatically Tuned Linear Algebra Software (ATLAS) homepage, http://math-atlas.sourceforge.net/. [28] Graphics Layout Engine (GLE) homepage, http://glx.sourceforge.net/. [29] M. M. Hopkins, P. S. Crozier, J. J. Boerner, L. C. Musson, T. P. Hughes, E. V. B. R. Hooper, and M. T. B. S. N. Laboratories), Ingredients for 3D Vacuum Arc Discharge Simulation, Poster at MevArc’12, https://indico.cern.ch/event/208932/contribution/8/ material/slides/1.pdf, 2012. [30] J. Boerner, M. Hopkins, E. Barnat, P. Crozier, M. Bettencourt, C. Moore, R. Hooper, and L. Musson, The Path to Full Geometry 3D Vacuum Arc Simulation, Presentation at MevArc’12, https://indico.cern.ch/event/208932/contribution/9/ material/slides/1.pdf, 2012. [31] Network Common Data Form (NetCDF), http://www.unidata.ucar.edu/software/ netcdf/. [32] Hierarchical Data Format 5, http://www.hdfgroup.org/HDF5/.