Interactive Introduction to Plasma Physics Part II: Hands-on exercises Adam Obrusník, Lenka Zajíčková Copyright © 2013 Masaryk University Licensed under the Creative Commons Attribution-NonCommercial 3.0 Unported License (the "License"). You may not use this file except in compliance with the License. You may obtain a copy of the License at http: //creativecommons .org/licenses/by-nc/3.0. Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License. First printing, September 2014 1 Introduction................................................. 5 1.1 What this text is 5 1.2 What this text is not 5 1.3 Connection to the lecture 6 1.3.1 Weeks 2 & 3: Getting ready....................................... 6 1.3.2 Weeks 4 to 6: Particles in fields..................................... 6 1.3.3 Weeks 7 & 8: Data processing..................................... 6 1.3.4 Weeks 9 & 10: Distribution functions ................................ 6 1.3.5 Weeks 11 to 13: Final assignment .................................. 6 2 Matlab fundamentals........................................ 7 2.1 Installing Matlab at the Masaryk University 7 2.1.1 MUNI License.................................................. 7 2.1.2 Getting the installation files....................................... 7 2.1.3 Installation .................................................... 8 2.1.4 Running Matlab................................................ 9 2.2 Matlab fundamentals 9 2.2.1 Variables and indexing ......................................... 10 2.2.2 M-files....................................................... 11 2.2.3 Where to go from here ......................................... 13 3 Motion of particles in E and B fields .......................... 15 3.1 Particles in constant E and B fields 15 3.1.1 Underlying equations........................................... 15 3.1.2 Solving the system of ODEs in Matlab.............................. 17 3.1.3 The Implementation............................................ 17 3.1.4 Exercises..................................................... 19 3.2 Van Allen radiation belt 21 3.2.1 Underlying equations........................................... 21 3.2.2 Implementation............................................... 22 4 Data processing............................................ 29 4.1 Introduction 29 4.2 Loops and conditions in Matlab 29 4.2.1 For loop ..................................................... 30 4.2.2 If... elif... else conditioning ..................................... 32 5 Particle interactions in plasma............................... 35 5.1 Maxwell-Boltzmann distribution 35 5.1.1 Mean speeds................................................. 35 5.1.2 Shape of the Maxwell-Boltzmann distribution ........................ 38 5.2 Time evolution of a distribution function 38 5.3 Rate coefficients 42 6 Particle balance in plasma.................................. 45 6.1 Formulating the problem 45 6.1.1 The source term............................................... 45 6.1.2 The reaction scheme........................................... 47 6.2 Implementation 48 6.3 Parametric study 52 6.4 Conclusion 53 Bibliography ............................................... 55 Books 55 Articles 55 Online resources 55 1.1 What this text is This is a supplementary text to the exercises to the lecture F5170: Introduction to Plasma Physics. The lecture itself contains many theoretical derivations and considerations and the aim of this text is to mediate hands-on experience to students who come in contact with plasma physics for the first time. Throughout the text, you, as a student, will be presented to a number of problems, most of which have to be solved with the help of computers. You will, among other things, learn to solve systems of ordinary differential equations, to batch-process data or to numerically integrate data, which can not be described by an analytical function. You will find most of the skills and numerical techniques that you will learn throughout the semester useful also in other disciplines of physics. The hands-on examples are designed to be as straightforward and illustrative as possible. It has been taken into account that most of the students may initially lack the programming skills, required for completing this course. For this reason, all the code is richly commented and documented. The idea is that you will learn the basics of Matlab programming on these real-life commented examples. 1.2 What this text is not First and foremost, this text does not aim to be a comprehensive guide to Matlab programming and scripting. There are countless books and manuals on Matlab and its frameworks. The most relevant and accessible ones are listed in section 2.2.3. Be prepared that the text may occasionaly be too concise, the language may be sloppy and the logic not obvious. That is because at this point, the text is not designed to work on its own but rather in combination with the exercises to F5170: Introduction to Plasma Physics. This text is also not a reference that the user would be advised to cite in scientific publications or theses of any kind. It is always best to refer directly to the books or articles that this text is built on. Although obtaining the original references just to make sure a particular formula or figure is correct may seem tedious and pointless, it is a skill that you will certainly find useful in your future scientific work. Finally, this text is not a text on numerical methods. Various differential equation solvers 6 Chapter 1. Introduction are considered black boxes and their inner workings are not analysed in detail. The authors avoid non-dimensionalization of all the equations that are solved because using meaningful physical units has the advantage that you will not only understand the phenomena qualitatively but you should also get an idea about the magnitudes of various quantities in plasma physics. 1.3 Connection to the lecture The lecture spans across 13 weeks of a Fall semester at the Masaryk University. This text reflects the contents of the lecture and the idea is that you will progress by one chapter every two weeks. There is usually no exercise in the first week of the semester. 1.3.1 Weeks 2 & 3: Getting ready Since the lecture has to have a bit of a head start, you will have whole two weeks to get the required software to work on your computer. Both Windows and Linux users will have to install the university-licensed Matlab software (see chapter 2). Students should also study the second part of this chapter to learn the absolute basics of Matlab code and logic. 1.3.2 Weeks 4 to 6: Particles in fields In weeks 4 to 6, you should understand all the Matlab programs presented in chapter 3 and by the end of week 5, you should complete all the exercises in the respective chapter. The exercises focus on movement of particles in electric and magnetic fields and should help you to understand the nature of various particle drifts in plasmas. 1.3.3 Weeks 7 & 8: Data processing You will learn basics of batch data processing using examples from plasma physics. This includes interpolation and integration of discrete data in Matlab as well as data conversion. Although you will work with data relevant for plasma physics, you will find the skills obtained in this section useful in any other discipline of physics. 1.3.4 Weeks 9 & 10: Distribution functions This short chapter should help you to understand distribution functions. You will be provided with programs that can visualize the time development of a distribution function. You will also analyze the shape of the Maxwellian distribution for various conditions. 1.3.5 Weeks 11 to 13: Final assignment In the final assignment, you will utilize the knowledge obtained in all the previous sections. You will have to write or complete a program that solves particle balance in a laboratory argon plasma. You will also learn how the composition of the plasma changes with the pressure or the temperature of electrons and heavy species. Installing Matlab at the Masaryk University MUNI License Getting the installation files Installation Running Matlab Matlab fundamentals Variables and indexing M-files Where to go from here 2. Matlab fundamentals 2.1 2.1.1 Matlab (originally Matrix Laboratory) is commercial software and a scripting language being developed by the company MathWorks. Unlike many traditional programming languages, Matlab allows the user to manipulate matrices and vectors in an intuitive way. It allows the users to define their own functions which can have virtually any input or output. Installing Matlab at the Masaryk University MUNI License As of 2014, Masaryk University has the newest license of Matlab with a license pool of at least 250 licenses. This means, that 250 students or employees of MUNI can run Matlab simultaneously. In order to work properly, the MUNI license requires the user to be connected to the Internet and to MUNI's virtual private network (VPN). The license is limited to research and non-commercial usage. Getting the installation files The installation files for Matlab 2014b are available from the intranet of the Masaryk University. Please keep in mind that the Matlab installation file is currently approximately 7.5 GB large but inside Masaryk University's network, you will be able to download the installation files extremely fast. Starting from version 2013, Matlab is available only for 64bit computers and operating systems. If you want to use Matlab on a 32bit architecture, please ask the advisor. Howto 2.1.1 — Downloading Matlab files. Follow these steps to get all the necessary files 1. Open inet.muni.cz in your web browser 2. Log in using your university number (UCO) and your primary password 3. Go to Software —> Nabídka softwaru 4. Locate Matlab in the software list and click Získat 5. After agreeing with the License agreement, download the following three to your computer • The license .dat file • The installation image, R201xx_Windows. iso for Windows machines or R201xx_UNIX. iso for Linux-based operating systems 8 Chapter 2. Matlab fundamentals • The Autorizační kód - store it in a text file since the installer will ask you for it. 2.1.3 Installation Having downloaded the ISO installation images, you have to mount them and run the setup. Mounting an ISO image means creating a virtual DVD drive, containing the files from the ISO. Alternatively, you could burn the ISO image onto a DVD and install it from there. Howto 2.1.2 — Mounting the ISO on Windows. The most widely used software for mounting ISO files (i.e. creating a virtual DVD drive with the ISO files as a content) is probably Daemon Tools Lite. Despite being the best software for mounting virtual drives, Daemon tools Lite forces a lot of unwanted software on the user. If you want to avoid an unwanted toolbar in your browser, go through the installation carefully, unchecking the corresponding option. After the installation of Daemon Tools, mount the Matlab DVD: 1. Right-click the Daemon tools icon in your system tray 2. Go to Virtual DVD —> Device 0: —> Mount Image 3. Locate the ISO image and confirm 4. You will now see the Matlab DVD in My Computer Howto 2.1.3 — Mounting the ISO on Linux. Many modern Linux distributions have built-in functionality for mounting ISO images. You can usually access this functionality by right-clicking the ISO file and choosing the Mount option. If you cannot find a built-in mounting function in your distro, you can use AcetonelSO software, which is available in repositories of all major Linux distributions. If you prefer a command line solution, you can type the following to the terminal, without having to install anything mount -o exec R201xx_UNIX.iso /mnt/disk with /mnt/disk being the directory where you want to mount the ISO On Windows, the installer is started by the setup.exe file, while on Linux, you have to execute the ./install binary. In both cases, the files are located in the top-lever directory of the installation DVD. When the installer asks you whether you want to Activate Matlab, click Yes. On Linux, you may want to install Matlab as a super-user. Otherwise, a link in /usr/local/bin will not be created and you will not be able to run it simply by typing matlab to the command line. The Matlab installation itself is straightforward and user-friendly so follow the installation guide. You will be asked to provide first the Autorizační kód, then choose the installation location and then provide the license file. At one point, you will be asked about installing toolboxes. The toolboxes are additional function libraries, extending the functionality of Matlab. If you want to save some hard-drive space, you do not have to install toolboxes from these categories • Test and measurement • Computational finance • Computational Biology 2.2 Matlab fundamentals 9 • Code generation 2.1.4 Running Matlab Windows users will find a link to Matlab on their desktop and in their Start menu. On Linux, you can run Matlab by the matlab command (if you installed it under superuser privileges) or from its installation directory, /path/to/Matlab/R20xx/bin/matlab l After starting Matlab, the main window will open. It should be identical on all operating systems. MATLAB R2014a - + X HOME PLOTS APRS J L£i _ .©©J P m £ □ New New Open Script ^| Find Files Ä ^Compare Import Data BjJ, New Variable ^ Analyze Code |EH| ppj ©Preferences tr> Open Variable - g> Run and Time _ @ Set Path PESoupces Save Simulmk Layout Workspace ^ Clear Workspace ~* Clear Commands "» Library » |||| Parallel » file VARIABLE code simulink environment

v ... •••/ (3.15) The column vector q contains the values of our independent variable (time) while the columns of the matrix q contain the values of our dependent variables (positions and velocities) at given times. 3.1.3 The Implementation We are going to need two files to solve the motion equation. The first one, odef un .m will contain the function evaluating the r.h.s. of the system of differential equation while the file solve .m will contain the initial solution itself and plotting of the result. These two files have to be placed in the same directory and only the latter will be executed directly from Matlab. Let us first look at the file odef un. m. Program 1: odefun.m function dqdt = odefun(t,q), l °/0 This is the input function for the 2 °/0 ordinary differential equation solver ode45 . 3 °/0 Input: independent variable t(time), 6-dim vector q 4 °/o q = [x ,y , z , v_x , v_y , v_z] 5 °/0 Output: time derivative of q, dqdt = dq/dt 6 7 qe = -1.602e-19; % Elementary charge in 8 Coulomb 18 Chapter 3. Motion of particles in E and B fields m = 9.109e-31; '/, Particle mass in kg 9 10 °/0 Now, we define position and velocity variables n °/0 using the vector of dependent variables, q. 12 % This would not be necessary but using 13 °/0 "x" instead of q(l), etc.. makes the code 14 °/0 below more readable 15 x = q( 1) ; 16 y = q(2) ; 17 z = q(3) ; 18 vx = q(4) ; 19 vy = q(5) ; 20 vz = q (6) ; 21 22 Ex = 0 ; 23 Ey = 2e-6; 24 Ez = 0 ; 25 Bx = 0 ; 26 By = 0 ; 27 Bz = le-7; 28 29 °/0 Now, we calculate the components of the 6-dim 30 % vector dq/dt . 31 dxdt = vx; °/0 dx/dt = v_x 32 dydt = vy; °/0 dy/dt = v_y 33 dzdt = vz; °/0 dz/dt = v_z 34 dvxdt = qe/m*(Ex + vy*Bz - vz*By) ; °/0 dv_x/dt 35 dvydt = qe/m*(Ey + vz*Bx - vx*Bz); °/0 dv_y/dt 36 dvzdt = qe/m*(Ez + vx*By - vy*Bx) ; °/0 dv_z/dt 37 % Finally , we have to assign the calculated 38 component s % to the vector variable dqdt , which is 39 % the output of the function 40 dqdt = [dxdt; dydt; dzdt; dvxdt; dvydt; dvzdt]; 41 42 end % end of function 43 Although the source code is richly commented, let us now look closer at some of its important parts. Apparently, the file odef un .m is a Matlab function file. Line 1 tells the computer that this file contains a function called odefun which takes the independent variable t and the vector of dependent variables q as the input and it returns its derivative q (written as dqdt in ASCII). On lines 8 and 9, the elementary charge and particle mass are defined. The new set of variables which is defined on lines 16 to 22 is not necessary for the correct operation of the whole program but it is much more illustrative to use x, y, etc ... on lines 35 to 37, where the derivatives are expressed. On lines 23 to 28, the electric and magnetic fields are defined. Defining the constants qO and m inside the function is not a good practice but we do it here to make the source code as clear as possible. Since the function odef un() will be called many times by the ODE solver, it would be better to define these constants as global variables. 3.1 Particles in constant E and B fields 19 Having denned the function, describing the right hand side of the system of ODEs (3.13), we can finally solve it and plot the results. The relevant commands are included in the file solve .m. Program 1: solve.m % This matlab script solves the equation of motion l % for a charged particle in presence of electric 2 % and magnetic fields. 3 4 °/0 First , we need to define the initial conditions . 5 °/0 Remember, q = [x , y , z , v_x , v_y , v_z] 6 qO = [0;0 ; 0;le2 ; 0;0]; i 8 °/0 Now we define the time interval , on which 9 °/t we want to solve the ODE . 10 ti = 0; n tf = 2.5e-3; 12 N = 1000; 13 timespan = linspace (ti , tf , N) ; 14 % The function linspace (ti, tf , N) creates 15 % a linearly spaced vector beginning at "ti" 16 % ending at "tf" with "N" steps/components. 17 18 °/t And now for the solving itself 19 °/t The notation @odefun tells the program 20 °/0 that the first argument is a function. 21 [t, q] = ode45(@odefun , timespan, qO) ; 22 °/0 The result will be a vector t and a matrix q with N rows. 23 24 % The following commands display the calculated trajectory 25 mesh([q(: ,1) q(:,l)], [q(:,2) q(:,2)], [q(:,3) q(:,3)], [t ( :) 26 t (:)],' EdgeColor ' , 'interp', 'FaceColor', 'none'); °/. plot view(2) °/0 sets the plot to top view 27 set(gca,'FontSize',16,'fontWeight','bold ') % sets font size 28 xlabel('x [m] ') °/0 sets title of x-axis 29 ylabel('y [m] ') °/0 sets title of y-axis 30 zlabel('z [m] ') °/0 sets title of z-axis 31 title(['Particle trajectory, ti=', num2str(ti), ' s, tf=', 32 num2str(tf), ' s' ]) % sets figure title 33 print -depsc '0utput.eps' % prints the figure to file 34 3.1.4 Exercises Exercise 3.1 Copy the source code presented in the previous section to your computer or use m-files provided in the Progl_Particle_motion directory. • Run the program with pre-defined values you should get a curve similar to figure 3.2. • What kind of drift is observed? • What is the direction of the drift velocity for an electron and a positron? 20 Chapter 3. Motion of particles in E and B fields Particle trajectory, ti=0 s, tf=0.0025 s 0.01 r „0.006- 0.002 - \ /\ y y~y \ \ \ \ ! \ 1 \l i I j i \ \ \ 1 \ \ \ / / / / i \ j \ \ / \ / v y v \ \ \ 0.01 0.02 0.03 0.04 0.05 Figure 3.2: A trajectory of a particle in E = (0,1(T6,0) V/m, B = (0,0,1(T7)T andq(0) = (0,0,0, lOOm/s, 0,0,0) Exercise 3.2 Change the charge qe and the mass m of the particle to those of a proton. • How many times do you have to increase the time scale in order to see the characteristic trajectory of the drift? • Compare amplitude of the oscillations of electrons and protons and the magnitude of their drift velocities (read these data from the plots). Exercise 3.3 Modify the Matlab program from exercise 3.1 in such a way that the motion equation is solved without electric field but with magnetic field quadratically growing in the direction of y E = (0,0,0) (3.16) B = (0,0,fiz) (3.17) \ 2 y Bz=B0^j-J (3.18) qo = (0,0,0,V2v0,V2v0,0) (3.19) In the equation above, yo is the characteristic distance that your particle travels. You can either find this empirically or you can estimate it from the theoretical velocity magnitude of your particle and the time scale tf of your simulation. Do you observe any drift? For what values of Bq, yo and vo did you observe it? What is the direction of the drift velocity for an electron and for a positron? ■ Advanced exercise 3.1 Modify the program from this section so that the electric field changes harmonically with time. We recommend using Ex = E0-cos(oit), (3.20) Ey = 0, (3.21) Ez = 0, (3.22) 3.2 Van Allen radiation belt 21 and setting the initial velocity to vx = 0, (3.23) vy = v0, (3.24) vz = 0. (3.25) Now, examine for several frequencies ft) (for example 106, 107, 108 and 109 Hz) the motion of a proton and an electron. • How do they react to the field? • Compare how easily the electron and proton are influenced by the external field for various frequencies. 3.2 Van Allen radiation belt The discovery of van Allen radiation belts is often seen as the first big discovery of the space age. Their existence was prediced by James van Allen and confirmed by the Explorer I satellite in 1958. The belts are a direct consequence of the Earth's magnetic field, which traps high-energy particles and prevents them from reaching the surface of the Earth. Since no shielding is ever perfect, there is always some flux of particles towards the surface of the Earth. These particles then ionize and excite the gas molecules in the lower layers of the atmosphere, which we observe as the aurora (pictured in the header of this chapter). There are two stable van Allen belts, the first is located at 1000 to 6000 km above the surface while the second is located at 13000 to 60000 km [Bakl2]. The third van Allen belt has been discovered very recently and it seems to appear and disappear over time [Scil3; Shp+13]. Since you should now have a very good understanding of charged particle motion in electric and magnetic fields, saying that the energetic particles are simply trapped in the magnetic field is perhaps too vague. In this section, we will modify the Program 1 we used before and use to calculate a trajectory of a high-energy proton in the Earth's magnetic field. In the process, we will learn what global variables are and how to use them in Matlab. You will also become familiar with several advanced plotting commands. This section is inspired by Lecture notes of S. Markidis [Marl3]. 3.2.1 Underlying equations The equations that will have to be solved and the method of their solution are identical to the previous section 3.1. The only principal difference is that the electric field will be zero and the magnetic field will be a strong function of position and will be approximated by the field of a magnetic dipole with constant magnetic dipole moment m. As mentioned above, we approximate the geomagnetic field by a field of a dipole. This is a very rough approximation which is valid only up to the distance of several Re (Earth radii). Above that the magnetic field is strongly non-symmetrical due to its interactions with the solar wind. Most of the scientific activity in this field, both experimental and theoretical, is endorsed by NASA under the Living with a star program (lws.gsfc.nasa.gov). Some very advanced theoretical models of the geomagnetic field and space weather are available at ccmc.gsfc.nasa.gov. In the vector form, the magnetic field of a dipole is Ho f 3r(m • r) m 4n \ r5 r3 ,(r) = $C^£l_-*| (3.26) 22 Chapter 3. Motion of particles in E and B fields where m is the magnetic moment, [xq the vacuum permeability and r = {x,y,z), (3.27 r=|r|. (3.28 Exercise 3.4 We will work in a cartesian coordinate system with the origin at the center of the Earth and the z-axis identical with the Earth's axis. In this system, the magnetic moment of the planet is m = (0,0,Af). (3.29) • Express components of the magnetic field B = (Bx,By,Bz) in cartesian coordinates. • What is the value of M if the geomagnetic field at the equator is 3.12 • 10~5 T? 3.2.2 Implementation Just like in the previous case, we have to define the function which will return the right hand sides of our system of ODEs. The corresponding function m-file, odef un .m looks very similar to the previous cases. Program 2: odefun.m function dqdt = odefun(t,q), l °/0 This is the input function for the 2 °/0 ordinary differential equation solver ode45 . 3 °/0 Input: independent variable t(time), 6-dim vector q 4 °/o q = [x ,y , z , v_x , v_y , v_z] 5 °/0 Output: time derivative of q, dqdt = dq/dt 6 7 °/0 In solve.m, we define some global variables. Now, 8 we need to tell °/0 Matlab , that we will use "m" and "qe" 9 global m; global qe ; 10 n % Re-naming position and velocity variables 12 x = q( 1) ; 13 y = q(2) ; 14 z = q(3) ; 15 vx = q(4) ; 16 vy = q(5) ; 17 vz = q(6) ; 18 19 °/0 The electric field is zero here 20 Ex = 0 ; 21 Ey = 0 ; 22 Ez = 0 ; 23 24 °/t The magnetic field is now calculated using another 25 function, bfield(r) which is a function of position. The function is stored in bfield.m 3.2 Van Allen radiation belt 23 rO = [q(l), q(2) , q(3)]; '/, The position vector 26 B = bfield(rO); % the function returns a vector B 27 % ... and we use this vector to define our component- 28 wise B-field variables Bx = B (1) ; 29 By = B (2) ; 30 Bz = B (3) ; 31 32 °/t From here on, it is~similar to Program 1 33 dxdt = vx; °/0 dx/dt = v_x 34 dydt = vy; °/0 dy/dt = v_y 35 dzdt = vz; °/0 dz/dt = v_z 36 dvxdt = qe/m*(Ex + vy*Bz - vz*By) ; °/0 dv_x/dt 37 dvydt = qe/m*(Ey + vz*Bx - vx*Bz); °/0 dv_y/dt 38 dvzdt = qe/m*(Ez + vx*By - vy*Bx) ; °/0 dv_z/dt 39 dqdt = [dxdt; dydt; dzdt; dvxdt; dvydt; dvzdt]; 40 41 end °/0 end of function 42 The first difference is found on line 10. Instead of setting the variables m and qe to some particular value, we defined them as global variables. If you are not familiar with the concept of global and local variables, see the remark below. A global variable is a variable that is available in all functions and programs where it is declared global. Normal variables (local) are not shared (inherited) between functions and programs. Imagine that you want to use the variable m both in solve .m and odefun.m • If you use local variables, you have to assign the value of the variable, i.e. m = 1.627e-27;, both in odefun.m and in solve.m. Therefore, you have to write the numerical value into every file. • If you use global variables, you can set the value m = 1.627e-27; only in solve, m and you declare the variable as global by writing global m; both to odefun.m and in solve .m. Therefore, the numerical value is written in one file only. The advantage of the second approach is apparent. If you want to change the mass m, you only have to re-write it once. The second difference in odef un. m is around line 26. You can see that the magnetic field is not defined explicitly but it is obtained from an external function bf ield(). The function is stored in the file bf ield. m, it also uses some global variables (also defined in solve. m) and you have to complete it by writing correct expressions for Bx, By, Bz. Program 2: bfield.m function B = bfield(rO), l % The function calculates the geomagnetic field B at 2 the position rO. It works with cartesian coordinates. % The function uses the following global variables: 3 global Re; global q; global m; global muO; global M; 4 5 % If you did the exercises , the lines below should be 6 clear 24 Chapter 3. Motion of particles in E and B fields x = rO(l) ; 7 y = rO(2) ; 8 z = rO(3) ; 9 r = sqrt (x~2+y~2+z~2) ; 10 n % Complete the following based on the exercise above 12 Bx = XXX; 13 By = YYY; 14 Bz = ZZZ; 15 16 B = [Bx; By; Bz] ; 17 end 18 Finally, we need a script file, which includes all the necessary variables and calls the ode45 solver. Program 2: solve.m % This matlab script solves the equation of motion 1 % for a charged particle in Earth's magnetic field. 2 3 % To make the program run faster , we will define all our 4 % constants as global variables. 5 global Re; global m; global qe; global muO; global M; 6 Re = 6378137; % Earth radius in meters 7 m = 1.627e-27; % Particle mass in kg 8 qe = 1.602e-19; % particle charge in Coulomb 9 esc = 1; X Earth scale - for plotting 10 c = 299792458; °/. Speed of light in m/s 11 mu0 = 4*pi*le-7; % Vacuum permeability in V*s/(A*m) 12 M = 7.94e22; °/0 Earth's magnetic moment in H/m 13 14 °/0 Initial conditions for position 15 °/0 Defined in spherical coordinates for convenience ... 16 rO = 3*Re; 17 phiO = 0; is thetaO = pi/2; 19 °/0 ... and transformed to cartesian. 20 xO = rO*sin(thetaO)* cos(phi 0) 21 yO = r0*sin(theta0)*sin(phi0) 22 zO = r0*cos (thetaO) 23 24 °/0 Initial conditions for velocity 25 °/0 defined using particle energy and two angles 26 Ek_eV = 5e7; % energy in eV 27 Ek = Ek_eV*l.602e-19; '/, energy in Joule 28 v_r0 = c*(l+m*c~2/Ek)~-0.5 % relativistic velocity magnitude 29 v_phi0 = 0; °/0 azimuthal angle in rad 30 v_theta0 = pi/4; °/0 polar angle in rad 31 32 °/0 Coordinate transformation for velocity. 33 3.2 Van Allen radiation belt 25 vxO = v_r 0* sin (v_thetaO) * cos (v_phi 0) 34 vyO = v_r 0* sin (v_thetaO) * sin (v_phi 0) 35 vzO = v_rO*cos (v_thetaO) 36 37 °/0 Defining the vector of initial conditions 38 qO = [xO ; yO ; zO ; vxO ; vyO ; vzO] ; 39 40 % Defining the time interval 41 ti = 0; % initial time in seconds 42 tf = 20; % final time in seconds 43 N = 10000; °/0 number of steps 44 timespan = linspace (ti , tf , N) ; 45 46 % Solving the ODE. 47 % Please note that odefun(t,q) is different from Program 1 48 [t, q] = ode45(@odefun , timespan, qO) ; 49 50 51 °/0 Having solved the equation, the data has to be plotted. 52 % It is not necessary to understand what each of the commands 53 % below does . 54 55 h = figure; % Creates a new figure 56 57 °/0 The following set of commands plots the magnetic lines of 58 force. [Adapted from MagLForce script by A. Abokhodair] n = 4; 59 d2r = pi/180; r2d=l/d2r; 60 tht=d2r*(0:5:360) ' ; 61 phi=d2r*(0:ceil(180/n):180); 62 hh=phi*r2d; 63 A = r0; 64 r=A*sin(tht) .~2; 65 rho=r . * sin (tht) ; 66 x=rho*cos(phi); y = rho * sin (phi ) ; 67 [nR , nC] =size (x) ; 68 u=ones (1 ,nC) ; z=r . * cos (tht) *u ; 69 plot3(x/Re,y/Re,z/Re,'r','LineWidth ' ,1 .0) ; 70 71 hold on; °/0 This tells Matlab that whatever you are going to 72 plot next will be included in the same figure ! ! ! 73 % Plotting an ellipsoid, representing the Earth. 74 [u , v , w] = sphere (30) ; 75 surf(esc*u, esc*v, esc*0.9*w); 76 colormap ('default') ; 77 camlight right ; 78 lighting phong ; 79 axis equal % This forces the axes to have the same scale! so 26 Chapter 3. Motion of particles in E and B fields % Plotting the particle trajectory and add labels to axes. plot3(q(:,1)/Re, q(:,2)/Re, q(:,3)/Re,'LineWidth',1.0) set(gca,'FontSize ' , 18) % sets font size xlabel('x [R_e]') °/» sets title of x-axis ylabel('y [R_e]') °/» sets title of y-axis zlabel('z [R_e]') °/» sets title of z-axis title([' Proton trajectory, E = 50 MeV, t_i=', num2str(ti) , ' s, t_f=', num2str(tf), ' s']) % sets figure title print -dpng -r200 ' ProtonMagField . png' °/0 prints the figure to file using the png format with 200 dpi resolution. The solve .m file is quite similar to the one we used in the previous one. On lines 6 to 13, global variables are defined. It is best practice that global variables declared global and assigned their values at the beginning of the first file that is executed. On lines 15 to 39, the initial conditions are defined. To make the them more intuitive, the initial conditions are defined in spherical coordinates and then transformed to Cartesian ones. The differential equation is solved on line 49 and from line 57 onwards, the plotting commands are listed. Please note the usage of the hold on command, which tells Matlab that all the next plots will be added to the same figure. Also note that all data that is plotted is scaled by the Earth's radius Re. If you completed the expression for the magnetic field correctly, running the program with the pre-defined values should give a plot similar to figure 3.3. Proton trajectory, E = 50 MeV, ti=0 s, tf=20 s Exercise 3.6 Even though our model is quite simple, it can be helpful in understanding the basic structures of the radiation belts. Try changing the initial position of the proton. • What is the maximum initial distance rO, for which the proton still has a stable 3.2 Van Allen radiation belt 27 trajectory? What is the minimum initial distance, for which the proton does not hit the surface of the Earth? Advanced exercise 3.2 Replace the proton in the previous model with an electron and try to find the initial conditions (especially the distance from the Earth), for which the electron will have stable trajectory. Think before you code, will electrons require higher or lower magnetic field to be confined? How does the drift of an electron differ from that of a proton? 4.1 Introduction The Matlab skills that you will learn in this chapter will be useful also outside plasma physics. In particular, the chapter focuses on data processing, including interpolation and integration of discrete data. Furthermore, you will learn to use so-called loops and conditions in order to batch-process large sets of data. 4.2 Loops and conditions in Matlab Let us first look at two very important elements of Matlab programming and programming in general, loops and conditions. If you have previous experience with programming, this section will probably seem quite trivial but if you do not, you will learn skills that you will almost certainly find useful in your future career of a scientist. Loops are stuctures that allow you to repeat a certain sequence of instruction multiple times. Such repetitive tasks in physics may include • plotting multiple data files in the same format, • conveting large amounts of data, • solving a differential equation for multiple sets of initial/boundary conditions. • ... There are two types of loops in Matlab, the for loop and the while loop. The former executes the specified sequence of commands N times (e.g. "integrate the data in first 50 files") while the latter keeps executing the sequence as long as some expression is true (e.g. "keep integrating the data until you reach an integral larger than 10"). However, the range of applications of the while loop is very narrow and it is quite easy to make it endless. For this reason, we will only focus on the for loop, which is completely sufficient for purposes of most physicists. Conditions are used when the sequence of commands that you want to execute depends on whether some statement is true or false. They are specially useful in combination with loops. An example you can think about is automated plotting; you can create a program that will plot the data in logartihmic scale if the difference between the lowest and the highest values is larger than two orders of magnitude. 30 Chapter 4. Data processing 4.2.1 For loop As already mentioned, the for loops makes it possible to perform a certain operation N times. The loop is used as follows. for j=l:10, l °/0 your commands 2 end 3 The loop above will run ten times. The variable j is called the loop variable and you can use it inside the loop. Let us now write a program that actually does something at least a little useful. In particular, we will create plots ofy = xn with n £ (1,6) and x £ (0,1). figure; % creates figure l colors = ['r', 'g', 'b', 'm', 'k', 'y', 'c', 'r']; % this is 2 a vector of strings (texts) for j=l:8, X will run for 8 values of j 3 x = linspace(0 ,1,50); % we saw this before, creates a 4 linearly spaced vector from 0 to 1 with 50 values, color = colors(j); % selects j-th element from colors 5 plot(x, x.~j, color); 6 hold on; °/0 we want everything in the same figure i end 8 Exercise 4.1 Modify the program above so that it plots siní nx) y= Kx ' forxG (0,27T)and«G (1,5) (4.1) Do not forget that x is a vector and you have to use element-wise operations. ■ In plasma physics as well as in other disciplines of physics, one often has to work with data obtained from literature and often it happens that the data have a different unit than you need. The next exercise program serves as a batch-converter of collisional cross-sections which have been specified in the unit of cm2 to the unit of m2. In order to write the program, we are going to need two more useful functions, loadO and dlmwrite (). The load function is very intuitive, it loads data into a matlab variable from a file. The input data file can have several formats (tab-separated data, comma-separated data, etc..) and it is loaded to a variable by the command matrix = load('datafile.dat '); i On the other hand, the function for saving data to an ASCII file is equally straightforward, dlmwrite('new_datafile.dat', matrix, '\t'); i As you can see, it takes three parameters, namely the output file, the variable that is saved to the file and the delimiter delimiting columns in the output file. Typically used delimiters are ',' (comma), ';' (semicolon) and ' \t' (tab). In the next exercise, you also need to know how to work with strings in Matlab. String is merely a programming term for a text. Since string in Matlab is basically a vector of letters, you can use the same syntax with the exception that strings have to be enclosed in quotation marks ('). 4.2 Loops and conditions in Matlab 31 Let us demonstrate how to join strings on a simple example. The following script will open and plot the content of files cskl.dat to csk3.dat (these files can be found in the Prog3_cross_section_convert directory). colors = ['r', 'g', 'b', 'm', 'k', 'y', 'c', 'r']; figure; for j=l:3, color = colors(j); matrix = load(['csk', num2str(j), '.dat']); xdata = matrix(:,l); °/0 first column ydata = matrix(:,2); °/0 second column plot(xdata, ydata, color); hold on; end xlim ( [0 , le6]) ; % set x-axis limits The data you just plotted are not random curves, they are cross sections of collisions of electrons with argon atoms. The x-axis is the electron temperature in Kelvin while the v-axis is the col-lisional cross-section in m2. In particular, the cskl .dat file contains the data for an elastic collision e + Ar-^e + Ar, (4.2) csk2 .dat the data for argon excitation e + Ar^e + Ar*, (4.3) and csk3 .dat the data for ionization e + Ar->2e + Ar+. (4.4) Exercise 4.2 The program above produces a plot, which is not very intuitive mainly because the magnitudes of the plots are very different. Modify the program so that • the electron energy is expressed in the unit of eV, • the v-axis scale is logarithmic. The resulting plot should look comparable to this one -18.5 "0 10 20 30 40 50 Electron energy [eV] with the red line corresponding to the elastic collision, the green line to excitation and the blue line to ionization. • Look at the picture, what are the excitation and ionization thresholds? 32 Chapter 4. Data processing Exercise 4.3 Modify the plotting program so that it converts the x-data from Kelvin to elec-tronvolt and saves each cross-section to tab-delimited csevN. dat with N being the corresponding file number. ■ 4.2.2 If... elif... else conditioning In the previous section, you learned about for loops which allow you to perform the same task repetitively. You can further step up your programs by using conditions. By including conditions in your programs, you can a certain set of commands if a statement is true and an another set of commands if it is false. Let's assume that you are developing a very complicated program, in which you often need to find inverse matrices. As you know, the inverse matrix A-1 to matrix A is defined as A-1 A I (4.5) where / is the identity matrix. Exercise 4.4 The Matlab function which allows you to calculate the inverse matrix is called inv() and its only argument is the matrix that you want to invert. • Define a 3x3 singular matrix A (determinant is zero) and try to invert it using inv (). What does the resulting matrix look like? As you saw in the exercise above, trying to calculate inverse matrix to a singular matrix does not lead to any useful result. It would be much more practical to define a custom function which will perform the inversion only if the input matrix is not singular. If the matrix is singular, your function will display an error and stop your program. The function we just described could look something like this function inverted = myinv(matrix), rows = length(matrix(:,1)); °/0 calculates the length of the first column, i.e. the number of rows cols = length (matrix (1 ,:)) ; '/, calculates the length of the first row, i.e. the number of colums if (rows ~= cols), °/. if rows DOES NOT EQUAL cols disp('The supplied matrix is not a square matrix') return °/0 terminates the script elseif(det(matrix) == 0); % if determinant is zero disp('The matrix is singular.') return °/0 terminates the script else, % if both the conditions above are false inverted = inv(matrix); °/0 performs inversion and returns the inverted matrix end end The syntax of the if . . . elseif. . . else conditioning should now be evident from the code above. You can see that the statement that needs to be evaluated is written in round brackets after if or elseif. Matlab tests the conditions from up to bottom, i.e. 1. The shape of the matrix is tested first. If not square, the error is displayed. 10 11 12 13 14 4.2 Loops and conditions in Matlab 33 2. If the matrix is square, the determinant is tested. If zero, the error is displayed 3. Only if the matrix is square and non singular, the inversion is performed. Exercise 4.5 Verify that the function myintO behaves correctly, try entering several non-square matrices and several singular matrices. Add another else if statement which will display a warning if your matrix rank is larger than 10. Test if it works. ■ Advanced exercise 4.1 Look into the directory Prog4_cross_section_convert_if which contains three files with collisional cross-sections. Two of these files include electron temperature in eV while one of them includes electron energy in Kelvin. • Modify the program that we used in exercise 4.2 using a reasonable if. . .else condition so that it decides which files should be converted and which not. • A comparable simple piece of code could be very useful when compiling a collisional-radiative plasma model with thousands of similar files that need to be preprocessed. However, what are the limitations of your program? Maxwell-Boltzmann distribution Mean speeds Shape of the Maxwell-Boltzmann distribution Time evolution of a distribution function Rate coefficients 5. Particle interactions in plasma This chapter illustrates the differences between equilibrium and non-equlibrium distribution functions and shows the impact of the shape of the distribution function on macroscopic variables. In this chapter, you will learn basics of numerical integration and interpolation in Matlab. 5.1 Maxwell-Boltzmann distribution 5.1.1 Mean speeds In the lecture, you learned that under equilibrium conditions, the velocity of particles is governed by the Maxwell-Boltzmann distribution. In particular, the distribution function for the velocity magnitude F(v) is 9 / m \3/2 /— mv2\ Distribution functions, especially those of electrons, can tell you a lot about your plasma, if you know how to read and interpret them. In particular, it is very important to understand what influence do distribution functions of different shapes have on macroscopic variables that everyone can intuitively understand. The macroscopic variables are very often given by quite simple integrals of distribution functions. Therefore, this section will teach you how to numerically integrate functions in Matlab. Some of the integrals that we will evaluate numerically in this section also have an ana-— lytical solution. However, in most laboratory plasmas, the distribution functions of some particles (especially electrons) can not be described by an analytical expression. Hence we will prefer numerical integration over analytical. Again, knowing how to integrate data numerically will be useful for you even in other disciplines of physics. There are several numerical techniques for numerical integration, the simplest one is the trapezoidal rule that you may remember from mathematical analysis lectures. In the trapezoidal rule, the function is approximated by trapezoids and the total surface area of the trapezoids gives you the estimate of the integral of the function between points a and b. The trapezoidal rule is schematically illustrated in figure 5.1(a). A more advanced method is the Simpson's rule which approximates the function f(x) to be integrated by several consecutive polynomials Pt(x) 36 Chapter 5. Particle interactions in plasma (see figure 5.1(b)). Therefore, the Simpson's rule can achieve better accuracy, especially with fast-varying functions. 9 9-l ^2 (a) Trapezoidal rule P2(x) 9 9-l 9 2 (b) Simpson's rule Figure 5.1: An illustriation of the trapezoidal rule for numerical integration (left) and the Simpson rule (right). In the Simpson rule, the upper side of the "trapezoids" is not linear but rather given by a polynomial, which can be integrated analytically. Figure 5.1(a) also suggests, that you can make the trapezoidal rule arbitrarily accurate by decreasing the size of the trapezoids. Therefore, the trapezoidal rule will be sufficient for our pourposes. The trapezoidal rule integration algorithm is implemented in Matlab under the name trapz(). The function takes two vectors as arguments. The first contains the x-data while the second contains f(x). x = linspace(0, pi, 100); fx = sin(x . ~2) ./log(x) ; I = trapz (x , fx) ; You do not have to specify the lower and upper bounds of the definite integral because these are given by the minimum and maximum values of the vector x. In other words, the two lines above calculate the follwing definite integral ™x sinx ■áx - 71 sin(x2) ■ dx (5.2) In the following exercise, you will learn to use the trapz () function and evaluate its accuracy on something that can be verified analytically. Exercise 5.1 Program 5 below is a template for plotting the Maxwell-Boltzmann distribution and for calculating and plotting the mean speeds of the distribution function. Complete the program by • filling in the expression for F(v) on line 7 (speed distribution function in 3D), • calculating the mean speed on line 12 and mean quadratic speed on line 13 (both using the trapezoidal rule), • providing analytical expressions for mean speed, mean quadratic speed and most probable speed on lines 16 to 18. Once you complete the program, answer the following questions. • Do the numerical results agree with the analytical ones? • Look at the plot. Which of the three speeds is the lowest, which is the highest? Does 5.1 Maxwell-Boltzmann distribution 37 their order change when you increase/decrease the temperature? • Try changing setting the number of steps in linspace () to 10, 50, 100, 500. What can you say about the agreement of numerical and analytical results? Program 5: Maxwell-Boltzmann distribution and speed mO = 1.6605e-27; % a.m.u. in kilograms l m = 40*m0; °/» particle mass in kilograms 2 kB = 1.380e-23; '/, Boltzmann constant in m~2 kg s~-2 K~-l 3 T = 1000; % temperature in Kelvin 4 5 v = linspace (0 ,4000 ,500) ; °/. x data 6 Fv = . . .; X distribution function i plot(v , Fv) ; 8 hold on; 9 10 % numerical expressions n v_mean = . . . ; '/, mean speed 12 v_sq_mean = . . . ; '/, mean quadratic speed 13 14 % analytical expressions 15 v_mean_an = . . . ; '/, mean speed 16 v_sq_mean_an = . . . ; '/, mean quadratic speed 17 v_mp_an = . . . ; '/, most probable speed 18 19 lineht = max(Fv)* 1.1; % line height, only for plotting 20 plot ([v_mean , v.mean] , [0, lineht], 'r'); °/» mean speed 21 (red) plot ([ v_sq_mean , v_sq_mean] , [0, lineht], 'g'); °/» mean 22 quadratic speed (green) plot ( [ v_mp_an, v_mp_an] , [0, lineht], 'm'); % most 23 probable speed (pink) hold off 24 Advanced exercise 5.1 In the program above, we calculated the mean velocity and the mean quadratic velocity numerically while the most probable velocity was only calculated analytically. To calculate it numerically, you would have to find a maximum of our distribution function (which is given by vectors v and Fv in our program). Use the internet to figure out a way to do it. ■ As you know from the lecture, the distribution function is normalized so that F(v)dv = n. (5.3) 0 Therefore, when you integrate the distribution function from some velocity vo to 00, you obtain the number of particles which have velocity greater or equal to vo- If you integrate F(v) from vo to v\, you obtain the number of particles with velocities in the interval (vo, v\). With this in mind, try to complete the exercise 5.2. 38 Chapter 5. Particle interactions in plasma Exercise 5.2 The mass of nitrogen molecules is 28 a.m.u. and their number density at ambient conditions is approximately 1.7 • 1025 m~3. How many nitrogen molecules in your room are faster than 50, 500, 1000, 2500, 5000 and 10000 m/s? ■ 5.1.2 Shape of the Maxwell-Boltzmann distribution You have probably noticed that the equilibrium Maxwell-Boltzmann distribution function is given only by two macroscopic parameters, in particular the mass of a particle m and the temperature T. The exercise below does not require any programming but if you do not know the answer right away, you will find the program in the previous section very helpful. Exercise 5.3 The following figure shows three equilibrium distributions of speed with unknown parameters m and T. • If the particle mass m is constant, which of the distribution functions will corresponds to lowest temperature and which to highest? • If, on the other hand, the temperature T is constant, which of the distribution functions corresponds to lowest and highest particle mass? 0 500 1000 1500 2000 2500 3000 Velocity magnitude, v [m/s] 5.2 Time evolution of a distribution function In the previous section, we thoroughly analyzed the key properties of the Maxwell-Boltzmann distribution. However, the Maxwell-Boltzmann distribution is only achieved in equilibrium and prior to achieving equilibrium, the distribution function can be almost arbitrary. In addition, we only took into account homogeneous and isotropic distribution functions in the previous section (i.e. depends only on velocity magnitude). In this section, we will still assume that the space is homogeneous (i.e. the distribution function does not depend on the coordinate r) but now, the distribution function can also be anisotropic (i.e. can depend on the velocity vector v, not only on its magnitude). The program that you will use for analyzing the time evolution of the distribution function is not very complicated. It is based on the simplest non-equilibrium formulation of the Boltz-mann kinetic equation. In this approximation, we assume that the external force acting on the particles is zero, F = 0, (5.4) 5.2 Time evolution of a distribution function 39 and the collision term is expressed in the Krook form df\ f-fo dt J -vm •(/-/<)) (5.5) coll where /o is the equilibrium distribution function, t is the relaxation time and vm is the collision frequency. The whole Boltzmann equation, therefore, reduces to ^|^ = -vm-(/(v,0-/o(v)) (5.6) which has an analytical solution of f(y,t) =/0(v) + [/(v,0)-/0(v)]e-v™f. (5.7) Therefore, the program actually does not solve the differential equation but only plots the solution in an intuitive way. The source code with extensive comments is listed below. The following lines are crucial for you when using the program • Line 5 sets the time interval, on which the solution is plotted • Line 6 sets the velocity interval, similar to the previous program • Lines 11 to 13 set the initial distribution function. • Line 32 changes how the plots scale. If it is commented out, the plots will have variable scale, if you uncomment it, the scale will be constant. Program 6: Plotting time evolution of a distribution function mO = 1.6605e-27; % a.m.u. in kilograms m = 40*m0; °/» particle mass in kilograms kB = 1.380e-23; '/, Boltzmann constant in m~2 kg s~-2 K~-l Vm = 5e8; % collision frequency, Hz time = linspace(0, le-8, 100); % time interval v = linspace (-2000 ,2000 ,200) ; °/. velocity interval % initial distribution function, f_0 sigma = 10; % Initial width vO = 500; % initial velocity fx_init = l/(sqrt(2*pi)*sigma)*exp(-(v-vO)."2/(2*sigma~2)); fy_init = v*0; fz_init = v*0; °/0 Calculating the equilibrium temperature , the following follows from energy conservation. v_sq = sqrt(trapz(v, (fx_init + fy_init + fz_init) . *v . ~2)) ; Teq = l/3*m*v_sq~2/kB; % equilibrium distribution function fx_eq = sqrt(m/(2*pi*kB*Teq)) * exp(-m*v.~2/(2*kB*Teq)); fy_eq = sqrt(m/(2*pi*kB*Teq)) * exp(-m*v.~2/(2*kB*Teq)); fz_eq = sqrt(m/(2*pi*kB*Teq)) * exp(-m*v.~2/(2*kB*Teq)); fx = fy = fz = fx_init; fy_init; fz_init ; 40 Chapter 5. Particle interactions in plasma hFig = figure; % creating the figure, setting dimensions set(hFig, 'Position', [100 300 1000 300]); for j=1:length(time), % for loop making the time steps absmax = max([fx, fy, fz])*l.l; % variable y-axis scale % absmax = max([max(fx_init),max(fy_init),max(fz_init)]); % uncomment this line for fixed y-axis scale subplot(1 , 3,1); % first subplot - x plot(v, fx, 'b', 'LineWidth', 2); ylabel('f(v_x) ') ; xlabel( ' v_x'); ylim([0, absmax]); whitebg('white') subplot (1 , 3 , 2) ; °/0 second subplot - y plot(v, fy, 'r', 'LineWidth', 2); ylabel('f(v_y) ') ; xlabel( ' v_y'); ylim( [0 , absmax]) subplot(1 , 3 , 3); % third subplot - z plot(v, fz, 'm', 'LineWidth', 2); ylabel('f(v_z) ') ; xlabel( ' v_z'); ylim([0, absmax]); °/0 setting figure title - time ax = axes('Units','Normal','Position' , [.075 .075 .85 .85],' Visible','off ') ; set(get(ax,'Title '),'Visible' , 'on') title(['t = ', num2str(time(j)) , ' s']); °/0 calculating new values of fx, f y , fz fx = fx_eq + (fx_init - fx_eq)*exp(-time(j)*Vm); fy = fy_eq + (fy_init - fy_eq)*exp(-time(j)*Vm); fz = fz_eq + (fy_init - fz_eq)*exp(-time(j)*Vm); M(j) = getframe(gcf); % appends a frame to variable M end movie2avi(M, 'out.avi', 'fps', 7); % saves the movie Exercise 5.4 Run the program above, watch the output animation and answer the following questions • What is the physical meaning of the initial condition for the distribution function? • What kind of particles could the distribution functions describe? • The collision frequency is 5 • 108 Hz, which is a resonable value. What is the time necessary for reaching the equilibrium? • Try increasing and decreasing the collision frequency in Program 6. What happens with the time necessary for reaching equilibrium and why? 5.2 Time evolution of a distribution function 41 Here is an example of what the output should look like at various times. Please note that the y-axis scale is changing. 0.03 0.025 „ 0.02 e- 0.015 0.01 0.005 01— -2000 x 10" 6 t = 0 ns 0.03 0.025 0.02 0.015 0.01 0.005 2000 -2000 x 10" 0 v y 3 ns 6 >-, £ 4 0.03 0.025 0.02 0.015 0.01 0.005 2000 01— -2000 x 10" 6 £ 4 2000 -2000 2000 -2000 2000 -2000 2000 x 10" 1.5 —x 1 0.5 -2000 1.5 x 10" 0.5 t = 6.5 ns x10~3 1.5 —>. 1 0.5 2000 -2000 1.5 t = 10 ns x10~3 0.5 x 10" 1.5 —M 1 > 0.5 2000 -2000 1.5 x 10" 0.5 2000 -2000 2000 -2000 2000 -2000 2000 42 Chapter 5. Particle interactions in plasma Rate coefficients Finally, let us demonstrate how different distribution functions can influence macroscopic properties, such as the rate coefficients of plasmachemical reactions. The rate coefficient of a reaction is calculated using the collisional cross section ar as p co kr = (ffr(v)v) = / ffr(v)vF(v)dv. (5.8) Jo As you can see, we will again be working only with the distribution function for speed because it is reasonable to assume that the cross-section of a particular reaction does not depend on the whole velocity vector but only on its magnitude. When calculating the integral (5.8), there is one technical issue that has to be overcome, the interpolation of the collisional cross-section <7r(v). Cross-sections of most reactions can not be expressed analytically and are available only as tabulated data. The interpolation is necessary beacause you need to compute the element-wise product ar(v)vF(v) and <7r(v) and F(v) do, generally, have different data sampling. In Matlab, you can interpolate ID data using the interpl () function. Let's say that you have a function given by vectors xdata (sample points) and ydata (values) and you want to obtain function values with much finer sampling. The syntax would be as follows xdata = [1 , 5 , 10 , 15, 20]; ydata = [1, 35, 110, 235, 410]; xinterp = linspace(8 ,17 , 0.5) ; yinterp = interp1(xdata, ydata, xinterp, 'pchip', 0); You can see that the interpl () function has five parameters, the vector of sample points (xdata), the vector of function values (ydata), the vector of new sample points (xinterp) and the interpolation method. The fifth parameter (zero) means that no extrapolation is performed and for all points outside the interval given by xdata, the function value is zero. With interpl(), you can choose from several interpolation methods. The most popular choices are 'nearest' which interpolates the data based on the nearest neighbour, 'linear' which performs linear interpolation, 'pchip' which performs cubic interpolation and ' spline' which interpolates the data using piecewise-continuous polynomials. Let us now examine, what will happen if we try to calculate the rate coefficient for electron impact ionization of argon e + Ar^2e + Ar+. (5.9) The calculation will be performed for two distribution functions, the Maxwell-Stefan distribution function and for a nearly monoenergetic beam with the same mean speed. The program below should not contain anything we have not encountered before. Its output are two numbers, the rate coefficient for Maxwell-Bolrzmann distribution and the rate coefficient for the monoenergetic beam. "Program 7: Electron impact ionization" m = 9.109e-31; % electron mass kB = 1.380e-23; '/, Boltzmann constant in m~2 kg s~-2 K~-l TeV = 15; % electron temperature in eV T = TeV*11604; °/0 electron temperature in Kelvin 5.3 Rate coefficients 43 5 % loading the cross - section 6 sigdata = load('sigmaion.dat'); 7 xsig = sigdata(:,1) ; 8 ysig = sigdata(:,2); 9 10 v = linspace(0,2e7,le6); % sampling for speed n sigma = interpl(xsig, ysig, v, 'pchip', 0); % interpolating 12 the cross - section 13 °/t Maxwell -Boltzmann distribution 14 Fv = 4*pi*v•"2 . *(m/(2*pi*kB*T))~1•5.*exp(-m*v."2/(2*kB*T)); % 15 distribution function kr_MB = trapz(v, v . *Fv .* sigma) 16 vmean = trapz(v, Fv.*v); 17 18 % Distribution function of a nearly monoenergetic beam 19 vO = vmean ; 20 width = vmean/100; 21 Fv2 = l/(sqrt(2*pi)* width) * exp ( - (v-vO) . ~ 2/(2* width "2)) ; °/0 22 distribution function kr_beam = trapz(v, v . * Fv2 . * sigma) 23 Exercise 5.5 Run the program and compare the two rate coefficients. You will notice that they are not the same, although the mean velocity is the same in both cases. Provide an explanation. ■ Exercise 5.6 Modify the program so that it plots a>(v) and answer the following questions. • What is the ionization threshold in electronvolts? • Where does the cross-section reach the maximum value. Exercise 5.7 Run the program for several values of electron temperature (defined in elec-tronvolt on line 3). • What happens with the rate coefficients with increasing electron temperature? • Is there electron temperature for which the rate coefficient for the nearly monoenergetic beam exceeds the Maxwell-Boltzmann coefficient? Provide an explanation why this is/is not possible. Advanced exercise 5.2 Rewrite Program 7 so that it uses electron energy rather than electron velocity for the calculations. Using electron energy is much more common in plasma physics. ■ There are several publicly available databases, from which you can obtain collisional cross-sections. The most comprehensive one is probably LxCat, available at lxcat.net. Another quite extensive database is the ALADDIN database, available at www-amdis.iaea.org/ALADDIN/. As you already know, plasmas consist of many types of particles (species). Apart from neutral ground-state atoms and molecules, there can be various excited species (as discussed in Interactive Introduction to Plasma Physics: part I), positive or negative ionic species and, of course, the electrons by which the plasma is sustained. Furthermore, each species can generally have different temperature. Knowing the temperature of individual species and their concentration is very important for the applications of plasmas. For example, in biological applications of laboratory plasmas, it is crucial to mainain high electron temperature while keeping the temperature of heavy particles low. 6.1 Formulating the problem In this last chapter of this text, you will implement a program that calculates the concentrations of individual particle species in atmospheric-pressure argon plasma for a given electron temperature Te and gas temperature Tg. The model will be again implemented in Matlab and it will be zero-dimensional. Practcally, this means that the number densities of each particle species are only a function of time, not position, na=na{t). (6.1) 6.1.1 The source term Under the zero-dimensional approximation that we made, the continuity equation for species a, simplifies to an ordinary differential equation dna V (6.2) dt This equation looks relatively simple but it has to be pointed out that the source term Sa is a complicated function of temperatures (which are constant in our model) and number denisties of other particle species. Generally speaking, the source term for particle species a is a sum of contributions from all reactions that include the species a, Sa=^Sa,r- (6-3) 46 Chapter 6. Particle balance in plasma The contribution Sa,r represents how many particles of species a are produced or consumed in reaction r per unit time. The source term depends on the order of the reaction. In our model, we will consider two orders of reactions, two-body reactions and three-body reactions. As the names suggest, two-body reactions are reactions between two particles and three body reactions are reactions between three particles at the same time. However, keep in mind that the order of the reaction tells you only how many reactants enter the reaction, it does not say anything about the number of products. In reality, no actual three-body collisions take place but plasma physicists often use this approximation at high pressures (including atmospheric pressure). A three-body reaction is actually given by two consecutive two-body reactions A + B—>C + D, (6.4) C + E—>F. (6.5) If the second reaction is faster and more probable than the first reaction, we can write the process as a three-body reaction A + B + E—>D + F. (6.6) An example of a two body reaction is the following electron-impact ionization of argon (later, we will designate this reaction R4). e + Ar* —> 2e + Ar+ (R4) (6.7) In this reaction, an electron collides with an excited argon atom and provides the energy necessary for ionization. The reaction will, contribute to source terms of all the reactants and products, in particular Se,R4 = +h(Te,Tg) -ne-nAr*, (6.8) >$Ar+,R4 = +h {Te, Tg) ■ ne ■ «Ar*, (6.9) >W,R4 = -h(Te,Tg)-ne-nAr*. (6.10) You can see, that the source term contribution contains the rate coefficient, which can depend on the temperatures of electrons and heavy particles. Our model is a two-temperature model, in which the electron temperature is different but all the other types of particles (ground-state, ions, excited species) have the same temperature Tg. You can also see that the contributions differ by the ± sign, which is also quite straightforward. Since reaction (R4) produces one new electron and one new ion, the sign of the corresponding source terms is positive. On the other hand, the reaction consumes one excited argon atom, therefore the corresponding source term contribution has a negative sign. An example of a three-body reaction is for example the formation of argon molecular ion. You may find the existence of this particle quite surprising, given that argon is an inert gas, but in high-pressure plasmas, they play quite an important role. The reaction leading to the formation of this ion is Ar+ + 2Ar —> Ar+ + Ar. (R8) (6.11) Analogical to the two-body reactions, the three-body reaction contributes to the source terms of individual species the following way ^Ar+.RS = +ks(Tg) ■ "Ar • «Ar • «Ar+, (6-12) ^Ar+,R8 = -h(Tg) ■ «Ar • «Ar • «Ar+> (6.13) ^Ar,R8 = -h(Tg) ■ nAr ■ nAr ■ nAr+. (6.14) 6.1 Formulating the problem 47 This time, the rate coefficient depends only on the temperature of heavy species Tg because no electrons enter the reaction. In addition, the rate coefficient will have a different unit (see exercise below). Exercise 6.1 In our formulation, the source term always has the unit of l/(m3 -s), what is the physical interpretation? What is the unit of the rate coefficient kr for two-body and three-body reactions? ■ 6.1.2 The reaction scheme The reaction scheme implemented in this work is adapted from an article by Baeva et. al. which focuses on numerical simulations of an atmospheric-pressure microwave plasma jet operating in argon [Bae+12]. The model includes the particle species listed in table 6.1. It should Table 6.1: The list of all species included in the argon plasma model Ar ground-state argon atom (concentration obtained from the state equation) Ar* excited argon atom (groups resonant and metastable 4s states of argon) Ar+ argon atomic ion ArJ argon molecular ion e electron be pointed out that the continuity equation (6.2) is solved only for four species, Ar*, Ar+, ArJ and e. The number density of ground-state argon can be determined from the state equation for the ideal gas, which works very good for monoatomic gases at high pressures. If the argon gas was not ionized, the number density of argon would simply be P nKx = —— (6.15) kBTg where p is the pressure, k# the Boltzmann constant and Tg the temperature of heavy particles. However, the equation above is not valid if a part of the argon gas is ionized or excited. Exercise 6.2 What is the actual number density of ground-state argon in plasma if the number density of excited atoms is «Ar*. the number density of atomic ions is «Ar+ and the number density of molecular ions is «Ar+ ? ■ The particle species listed in table 6.1 above can undergo 11 reactions in total. The complete set of these reactions is listed in table 6.2. In this reaction scheme, the authors assumed that the energy distribution functions of electrons and heavy particles are both Maxwellian. This, combined with advanced fitting, allowed them to express the rate coefficients analytically. Exercise 6.3 Reactions (Rl), (R3) and (R4) in table 6.2 all include an exponential part exp(—C/Te) which is zero at low electron temperatures and approaches one at higher electron temperatures. • What do you think is the unit and the physical meaning of the constants 11.65, 15.76 and 4.11? • Looking at the rate coefficients of these three reactions, which one will probably be the most important ionization channel? 48 Chapter 6. Particle balance in plasma Table 6.2: The reaction scheme for high-pressure argon plasma (adapted from [Bae+12]). The electron temperature is always in eV while the gas temperature is in K. reaction rate coefficient unit (Rl) e + Ar—>e + Ar* 4.9- 1(T15 • v%-exp(-11.65/re) m3/s (R2) e + Ar*—>e + Ar 4.8 • 10~16 • v% m3/s (R3) e + Ar—>2e + Ar+ 1.27 • 1(T14 • v%-exp (-15.76/7;) m3/s (R4) e + Ar*—>2e + Ar+ 1.37 • 10~13 • v%-exp (-4.11/Te) m3/s (R5) e + e + Ar+—>e + Ar 8.75 • 10~39 • T~4-5 m6/s (R6) e + Ar+^Ar + Ar* 1.04 • 10'12 • (re/0.026)-°-67 • t^TexpMWr.) mVs (R7) e + Ar+—>e + Ar + Ar+ 1.11 • 10 12 • exp (-2.94 - 3 m3/s (R8) Ar++2Ar—>Ar + Ar+ 2.25 • 10-43(7g/300)-°-4 m6/s (R9) Ar + Ar+—>Ar++2Ar 0.522 • lO-15?^1 exp(-15131/7g) m6/s (RIO) Ar*+Ar*—>e + Ar+ + Ar 6.2-10~16 m3/s (Rll) Ar*+Ar—>Ar + Ar 3.0-10~21 m3/s IExercise 6.4 Express the source terms of the following four species: e, Ar*, Ar+, ArJ using rate coefficients k\ tok\\ and the species' denisties «Ar,«e,"Ar*,«Ar+ >nArJ■ ■ As you can see, even for the relatively simple argon atom, the number of possible reactions is quite high. The number of reactions taking place in plasmas rises dramatically with the complexity of the gas in which plasma is ignited. With molecular gases, various rotational and vibrational excitations have to be taken into account, as well as dissociation and re-association of the molecule. Just to get the idea of the complexity, you would need several dozen reactions to provide a minimum description of a plasma in an O2/N2 mixture. In order to describe plasmas in argon with ambient air (including humidity and other admixtures), the number of reactions rises to more than 4000 [GB13]. Implementation Now we know everything in order to solve the system of ordinary differential equations in the following form d_ di ( "e \ «Ar* «Ar+ V "AT+ / Sat* \ sa4 J (6.16) We have solved a system of mathematically similar ordinary differential equations (ODEs) in chapter 3. Therefore, if you are unsure about anything regarding the code, try looking back at the exercises where we dealt with motion of particles in E and B fields. Since this is the final chapter of this text, you will only be provided with a template of the whole program and you have to fill in the missing pieces. The main file containing the solution looks like this: "Program 8: solve.m" % setting global variables global p; 6.2 Implementation 49 global kB; global Tg; global Te ; p = le5; % pressure in Pa kB = 1.38e-23; % Boltzmann constant in m~2 kg s~-2 K~-l Tg = 400; % Gas temperature in Kelvin Te = 2; °/0 Electron temperature in eV tsteps = 1000; % number of time steps tspan = logspace(-11, -6, tsteps); % this time, we do not use linear spacing for the time but rather logalithmic % Initial number densities in m~-3 nArs_init = lel2; °/, Ar* nArp_init = lel2; °/, Ar + nAr2p_init = lel2; °/, Ar_2 + ne_init = nArp_init + nAr2p_init ; °/0 electrons , follows from global neutrality °/t The vector of initial densities initial = [nArs_init , nArp_init , nAr2p_init, ne_init]; % Solving the system of (DDEs [t, sol] = ode45( ' odefun ' , tspan, initial); % Plotting the data close all figure ; hold on; nAr = p./(kB*Tg)-sol( ,2) -0 plot(loglO (t) plot (loglO (t) plot(loglO (t) plot (loglO (t) plot(loglO (t) ylim([12, 26]) legend('Ar', loglO(nAr) loglO (sol ( loglO(sol ( loglO (sol ( loglO(sol ( 5*sol( c>) ; ,2)) ,3)) ,4)) ,3)-sol(:, ) ) ; ); ); ); Ar ~* ' 'northwest') 1 Ar~ + Ar 2~ + electrons 'Location xlabel('Time, log_{10} [s]'); ylabel('Number density, log_{10} [m~{-3}]'); title(['p=', num2str(p), ' Pa, T_g=' , num2str(Tg), ' K, T_e=' , num2str(Te), ' eV ']) ; set(gca,'fontsize ' , 16); The program above should not surprise you now. The biggest difference with our previous programs solving ordinary differential equations lies on line 12. In particular, we do not use a linearly-spaced vector for time and we use a logarithmically spaced vector instead. This is because the processess in plasma are often logarithmic with respect to time. Similar to chapter 3, the solve. m program uses the ode45 () solved to calculate the time evolution of our system. The function evaluating the right-hand side of our system of equations (6.16) is called odef un() and is stored in a Matlab function file with a similar name. 50 Chapter 6. Particle balance in plasma "Program 8: odefun.m" function dqdt = odefun(t, q), % We will use global variables defined in solve.m global p; global kB; global Tg; global Te; % the vector q contains the densities of Ar*, Ar+, Ar_2+ and electrons nArs = q(1); nArp = q(2); nAr2p = q(3); ne = q(4) ; % the density of argon is calculated using the state equation and the densities of other heavy species nAr = p ./(kB*Tg)-nArs-2*nAr2p-nArp; % The rate coefficients are all stored in separate function files kl = f_kl(Te , Tg) k2 = f_k2(Te, Tg) k3 = f_k3(Te , Tg) k4 = f_k4(Te , Tg) k5 = f_k5(Te , Tg) k6 = f_k6(Te , Tg) k7 = f_k7(Te , Tg) k8 = f_k8(Te , Tg) k9 = f_k9(Te , Tg) klO = f_klO(Te, Tg); kll = f_kll(Te , Tg); % Now, it is necessary to define the source terms SArs = +kl*ne*nAr -k2*ne*nArs . -k4*ne*nArs . +k6*ne*nAr2p -2*klO*nArs*nArs . . . +kll*nArs*nAr; SArp = . . . ; SAr2p = . . . ; Se = ... ; % and finally , the derivative of the input vector is returned dqdt = [SArs; SArp; SAr2p; Se] ; end In the odef un() function, the source terms for individual particle species have to be calculated. The source term for Ar* (variable SArs) has already been pre filled but the source terms for Ar+ (variable SArp), ArJ (variable SAr2p) and electrons (variable Se) have to be completed based on exercise 6.4. Furthermore, the rate coefficients kl to kll are defined using eleven external functions f_kl() to f_kll() which take the gas temperature Tg in Kelvin and electron temperature Te in eV as input. Of course, not all the rate coefficients will depend on both Tg 6.2 Implementation 51 and re but for consistency, both these temperatures are passed to the functions. You will have to create these 11 simple functions based on the reactions and rate coefficients in table 6.2. Exercise 6.5 Complete the program, the template for which was provided above. In particular, you will have to 1. Writing the expressions for particle source terms in file odef un. m. 2. Creating 11 function filenení vyloučeno, že autoři neudělali chybus, f _kl .mtof _kll .m. Each of these files will calculate the rate coefficient of the corresponding reaction. Please note that the functions must be able to work with vector arguments, i.e. you have to use element-wise operations with Tg and Te. 3. Run the code with pre-defined values of pressure p, Tg and Te. Check if the output looks similar to the figuře below. p=100000 Pa, T =400 K, T =2eV g e 26 r ? 24-E Time, log1Q[s] Please take into consideration that this is the first time this text has been used in class and the authors may have also made a mistake in the source terms. If your output looks a little different and you are sure that your source terms are correct, do not worry. Now, when your program works, you can try changing various parameters in order to get a qualitative and quantitative idea how the argon plasma works. Exercise 6.6 Try increasing and decreasing the electron temperature and answer the following questions. • How does the equilibrium electron density change and why? • How does the ignition time change? • What is the dominant ion in the ignition phase and what is the dominant ion when the plasma stabilizes? Does the dominance of the two ions change with electron temperature? Exercise 6.7 Run the program for your chosen value of electron temperature and for pressures of 103, 104, 105 and 106 Pa and answer the following questions. • How does the equilibrium electron density change and why? • How does the ignition time change? 52 Chapter 6. Particle balance in plasma Parametric study Programs which are in principle similar to the program that we developed in this section are often used in plasma physics for examining plasma properties at various conditions, although the number of reactions taken into account is typically much higher (several thousand). It is often desirable to know how the plasma composition changes at various electron temperatures, neutral gas temperatures or pressures. Advanced exercise 6.1 Modify the program in the previous section so that it solves the system of equations for several values of TJp and plots the steady-state number densities as functions of TJp. You can do this by adding a for loop. What makes this exercise difficult is the fact that the ignition time changes quite quickly with electron temperature and pressure. Therefore, the time interval has to be updated in each iteration according to the current value of Te and p, otherwise the solution will take very long. ■ When solving similar sets of equations for a large number of conditions, advanced numerical techniques have to be employed to achieve reasonable computation times. Sometimes, the equations are solved in a logarithmic form, i.e. not for the number densities of particles but for their logarithms, lj = \nnj. (6.17) This trick can make the convergence much faster because the values of rij vary by many orders of magnitude during the solution while values of lj vary only within an order of magnitude. The exercise above may be quite difficult and unnecessarily time-consuming for students who do not intend to use Matlab on regular basis. However, even if you didn't implement the program in the advanced exercise above, try to complete the following one. 6.4 Conclusion 53 Exercise 6.8 As already mentioned, the composition of plasma often changes dramatically with electron temperature. The figure below shows the number densities of various particle species as a function of electron temperature at the constant pressure of p = 105 Pa and gas temperature Tg = 400 K. Think about the following questions and answer them. • If you wanted to use the plasma as a light source, what temperature would you use and why? • What is a simple explanation for the decrease in ArJ with increasing Te? p=100000 Pa, T =400 K 26r_, 60 0.5 1 1.5 2 2.5 3 Electron temperature, [eV] Conclusion This brings us to the end of the study material. If you have come thus far and completed all the exercises successfully, accept our congratulations. The authors sincerely hope that the exercises and demonstrations in this study text helped you to deepen the understanding of the complex processes taking part in plasmas. We also hope that the Matlab skills that you learned here will be useful to you, no matter what discipline of physics you end up in. Fields gauged by Poritv more: pore * Sociology is Pí^kůwgv is biduwy is wch is nust oh. hd; i didn't TuST^PPLJED Jtór APPLIED JUST APPLIED APFUED PHY9CS. 5ee YOO WSPiU- PSVtHOLOGy BIOLOGY. CHEM1STW its MICE TO the W\V OVER THERE. SXIOOjlSnS PSYCHOLOGISTS GlW-QGlSTS CHEMISTS PHYSICISTS MATHErWiCIArJS Figure 6.1: Downloaded from xkcd.org Books [Bit04] [GK08] J. A. Bittencourt. Fundamentals of Plasma Physics. 3rd edition. New York: Springer, 2004 (cited on page 12). Dan M. Goebel and Ira Katz. Fundamentals of Electric Propulsion: Ion and Hall Thrusters. 1st edition. Pasadena: California Institute of Technology, 2008 (cited on page 16). Articles [Bae+12] [Bakl2] [GBl 3] [Shp+13] M. Baeva et al. "Modeling of microwave-induced plasma in argon at atmospheric pressure". In: Physical Review E 85.5 (May 2012), page 056404. ISSN: 1539-3755. DOl:10.1103/PhysRevE.85.056404. URL: http://link.aps.org/doi/10. 1103/PhysRevE. 85.056404 (cited on pages 47, 48). Daniel N Baker. "New Twists in Earth's radiation belts". In: American Scientist 102 (2012) (cited on page 21). W Van Gaens and A Bogaerts. "Kinetic modelling for an atmospheric pressure argon plasma jet in humid air". In: Journal of Physics D: Applied Physics 46.27 (2013), page 275201. URL: http: //stacks. iop. org/0022- 3727/46/i=27/a=275201 (cited on page 48). Yuri Y Shprits et al. "Unusual stable trapping of the ultrarelativistic electrons in the Van Allen radiation belts". In: Nature Physics 9.11 (2013), pages 699-703. ISSN: 1745-2473. DOl: 10 .1038/nphys2760. URL: http://www.nature.com/ doif inder/10.1038/nphys2760 (cited on page 21). Online resources [FEI10] FEI. Introduction to Electron Microscopy. July 2010. URL: http: //www. f ei . com/ documents/introduction-to-microscopy-document/ (cited on page 15). 56 Chapter 6. Particle balance in plasma [Marl3] Stefano Markidis. Lectures on Computational Plasma Physics. Sept. 2013. URL: https : / /www. pdc . kth . se/education/computational-plasma-physics/ computational-plasma-physics (cited on page 21). [Scil3] Scitechdaily. Scientists Explain the Formation of the Unusual Third Van Allen Radiation Ring. Sept. 2013. URL: http : / /scitechdaily . com/scientists-explain - format ion - unusual - third - van - alien - radiat ion - ring/ (cited on page 21).