E L E C T R O N A N D ION R E L A X A T I O N IN A T O M I C A N D M O L E C U L A R MODERATORS By K i Yee Leung B. Sc. (Chemistry), Queen's University, 1987 M . Sc. (Chemistry), University of British Columbia, 1990 A T H E S I S S U B M I T T E D T H E IN P A R T I A L R E Q U I R E M E N T S D O C T O R F O R O F T H E F U L F I L L M E N T D E G R E E O F P H I L O S O P H Y in T H E F A C U L T Y O F G R A D U A T E S T U D I E S (Department of Chemistry) We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y O F BRITISH C O L U M B I A July 1997 Â© K i Yee Leung, 1997 O F In presenting degree freely at the available copying of department publication this of in partial fulfilment University of British Columbia, for this or thesis reference thesis by this for his thesis and scholarly or for her The University of British Columbia Vancouver, Canada I I further purposes gain the shall requirements agree that agree may representatives. financial permission. DE-6 (2/88) study. of be It not that the Library by understood be an advanced shall permission for granted is for allowed the that without make it extensive head of copying my my or written Abstract The standard technique for the solution of kinetic theory and transport problems usually involves the expansion of the velocity distribution function ( V D F ) in some suitable set of basis functions. For steady problems, the Boltzmann equation (BE) for the distribution function is then reduced to algebraic equations for the expansion coefficients. For time dependent problems, the B E is reduced to ordinary differential equations. The main objective with this solution technique is to calculate the transport coefficients which can be expressed in terms of the lower order expansion coefficients. Recent experiments on different physical systems are capable of probing the details of the velocity distribution functions. The main objective of this thesis is to develop accurate and efficient numerical methods for the B E with the aim to calculate the distribution function. For an ensemble of ions dilutely dispersed in a background of neutrals under the influence of an applied electric field, the exact solution of the Bhatnagar, Gross, and Krook ( B G K ) collision model is used to study different numerical solutions of the B E . In this study, the main theme of the thesis is presented. In place of the traditional expansion techniques, a discretization of the V D F is proposed. The discretization is based upon a set of polynomial basis functions that yield the optimum convergence of the solution. For electrons dilutely dispersed in a background of neutrals, the B E is approximated by a Fokker-Planck equation (FPE). The one-dimensional F P E has been used to describe many different phenomena in chemistry, physics and engineering. A detailed comparison of the eigenvalues and eigenfunctions of the F P operator is carried out with different numerical methods including a finite difference (FD) scheme, Lagrange interpolation n (LI) method, the Quadrature Discretization Method (QDM) with speed points and the Q D M with Davydov points. The Q D M involves the discretization of the distribution at a set of points that coincide with the points in a quadrature that is based upon some weight function. The Davydov distribution is the steady distribution for electrons in a gas moderator at some given electric field strength. The time dependent solutions of the FP are obtained with the different energy discretizations applied in conjunction with an appropriate time stepping procedure. The Q D M with Davydov points provides a more rapid convergence of the eigenvalues and eigenfunctions relative to the other methods used. The relaxation of a nonequilibrium distribution of electrons in a mixture of CCI4 with either A r or Ne is studied. The electron collisions in this analysis include elastic collisions between electrons and CCI4 and the inert gas, vibrationally inelastic collisions between electrons and C C I 4 , as well as the electron attachment reaction with C C I 4 . The time dependent electron energy distribution function is determined from the B E and the energy relaxation times are determined. The coupling of the thermalization process and the attachment process is discussed in detail. The results from the calculations are compared with experimental studies, and the methodology of the experimental reduction of the data is examined in detail. In addition, the nature of the spectrum of the Boltzmann collision operator for realistic ion-atom interactions is studied with a discretization of the integral collision operator. The singular nature of the collision operator is considered in detail with a novel interpolation technique. The eigenvalues and eigenfunctions for the Henyey-Greenstein (HG) model differential cross section are calculated. The relaxation of anisotropic ion distributions with and without an applied electric field is studied. 111 Table of Contents Abstract ii Table of Contents vi List of Tables vii List of Figures ix Acknowledgement xiii 1 Introduction 1 1.1 General methods for ion transport system 2 1.2 General methods for electron transport system 6 1.3 Theme and outline of the thesis 7 2 The Generation of Basis Functions for Kinetic Theory Problems 10 2.1 Introduction 10 2.2 Ion velocity distribution function for the B G K model . 11 2.3 The expansion of B G K solution in Hermite polynomials 14 2.4 Gautschi Stieltjes procedure for the construction of non-classical basis functions 19 2.5 Quadrature Discretization Method (QDM) 29 2.6 Ion velocity distribution function for the time dependent B G K model . . 31 iv 3 Numerical Methods of Solution of the F P E for Electron Thermalization 35 3.1 Introduction 35 3.2 The Lorentz-Fokker-Planck equation; electron thermalization 39 3.3 Numerical Solution of the Fokker-Planck Eigenvalue Problem 42 3.3.1 Quadrature Discretization Method (QDM) 42 3.3.2 Lagrange Interpolation 44 3.3.3 Finite difference scheme 45 3.4 3.5 4 5 Time Dependent Solutions of the Fokker-Planck Equation 46 3.4.1 Finite Difference Time Discretization 47 3.4.2 Q D M and Lagrange Interpolation Results and discussion â€¢ â€¢ â€¢ 48 . . . 48 3.5.1 Calculation of eigenvalues and eigenfunctions 48 3.5.2 Calculation of electron velocity distribution function 62 Electron Thermalization and Attachment in C C L j / A r and C C l / N e 78 4.1 Introduction 78 4.2 The Boltzmann Equation 81 4.2.1 Inelastic collision operator 81 4.2.2 Cross Section Data Set 84 4 4.3 Calculations and Results 88 4.4 The steady electron distribution at long times 105 Relaxation of Ion Distribution Functions by Ion-Neutral Collisions 116 5.1 Introduction 116 5.2 Boltzmann collision operator 121 5.3 Relaxation of anisotropic ion distribution function 124 5.3.1 124 Quadrature Discretization Method (QDM) v 5.4 6 5.3.2 The Multi-Domain (MD) approach 127 5.3.3 Results and discussion 129 Numerical solution of the time dependent B E with applied electric field . 5.4.1 Numerical algorithm 5.4.2 Results and discussion 136 â€¢ â€¢ â€¢ 136 140 S u m m a r y and O u t l o o k 144 Bibliography 148 vi List of Tables 2.1 Comparison of the convergence of with Gautschi Procedure with drift term, D' = 1 2.2 23 Comparison of the convergence of ct^ with Gautschi Procedure with drift term D' = 3 2.3 24 Comparison of the convergence of with Danloy's Procedure with drift term D' = 1 2.4 25 Comparison of the convergence of with Danloy's Procedure with drift term D' = 3 26 2.5 Convergence of the expansion of the B G K solution with different D' values 28 3.1 Comparison of the convergence of eigenvalues with Q D M , LI, and F D approaches for electrons in helium at E/n = 0.0 T d (10~ Vcm ) 17 3.2 2 Comparison of the convergence of eigenvalues with Q D M , LI, and F D approaches for electrons in argon at E/n = 0.0 T d 3.3 55 Comparison of the convergence of eigenvalues with Q D M , L I , and F D approaches for electrons in argon at E/n = 0.25 T d 3.6 54 Comparison of the convergence of eigenvalues with Q D M , LI, and F D approaches for electrons in helium at E/n = 0.50 T d 3.5 51 Comparison of the convergence of eigenvalues with Q D M , L I , and F D approaches for electrons in helium at E/n = 0.25 T d 3.4 50 56 Comparison of the convergence of eigenvalues with Q D M , LI, and F D approaches for electrons in argon at E/n = 0.50 T d Vll 57 3.7 Comparison of the convergence of eigenvalues with Q D M , LI, and F D approaches for electrons in xenon at E/n = 0.25 T d 3.8 59 Comparison of the convergence of eigenvalues with Q D M , LI, and F D approaches for electrons in xenon at E/n = 0.50 T d 4.1 60 Electron energy relaxation times pr^! and ratio between E^/Eth for CCI4 in A r at different concentration 4.2 96 Electron energy relaxation times pri.i and ratio between E^/E h t for CC1 4 in Ne at different concentration 5.1 97 Convergence properties of eigenvalues of the collision operator (only with H = 0) with Q D M for H G fit cross sections 5.2 Convergence properties of eigenvalues of the collision operator (only with Â£ = 0) with M D approach for H G fit cross sections 5.3 132 133 Reduced drift velocity for hard sphere cross section with unit mass ratio at different E/n 142 vm List of Figures 2.1 Ratio of the expansion of the B G K solution in Hermite Polynomials to the actual solution with different D' values 17 2.2 Convergence of the expansion of the B G K solution in Hermite polynomials 18 2.3 Time evolution of the B G K ion V D F with different D' at various initial conditions 3.1 34 Interpolation error for the LI scheme versus N with equi-spaced grid and Legendre quadrature points 3.2 53 Davydov distributions for electrons in helium and argon at different E/n values 3.3 61 Eigenfunction ^ 3 for F P operator with helium cross section at E/n = 0.25 T d and N = 5, 8, 10 and 15 3.4 63 Eigenfunction ip3 for F P operator with helium cross section at E/n = 0.50 T d and N = 5, 8, 10 and 15 3.5 64 Eigenfunction ipr for F P operator with helium cross section at E/N â€” 0.25 T d and N = 8, 10, 15 and 30 3.6 65 Eigenfunction ^7 for F P operator with helium cross section at E/n = 0.50 T d and N = 8, 10, 15 and 30 3.7 66 Eigenfunction ip for F P operator with argon cross section at E/n = 0.25 3 T d and N = 8, 10, 15 and 30 3.8 67 Eigenfunction V> for F P operator with argon cross section at E/n = 0.50 3 T d and N = 10, 15, 20 and 30 . ix 68 3.9 Eigenfunction tp for F P operator with argon cross section at E/n = 0.25 7 Td and N = 10, 15, 20 and 30 69 3.10 Eigenfunction ^ 3 for F P operator with argon cross section at E/n = 0.50 Td and N = 10, 15, 20 and 30 70 3.11 Exact solution of the time evolution of the V D F for electrons in xenon at E/n = 0.25 and 0.5 T d calculated with F D scheme 3.12 Time evolution of the V D F for electrons in xenon at E/n 72 = 0.25 T d calculated with various schemes (QDM, LI and FD) 3.13 Time evolution of the V D F for electrons in xenon at E/n 73 = 0.50 T d calculated with various schemes (QDM, LI and FD) 75 3.14 Variation of the maximum eigenvalue AN for the F P operator versus no. of grids, N 4.1 77 Energy dependence of momentum transfer cross sections for CCI4, Argon and Neon 4.2 85 Energy dependence of the four vibrationally inelastic cross sections for electron-CCU collisions 4.3 87 Energy dependence of the attachment cross sections for electron-CCLi collisions 89 4.4 Electron relaxation in C C L i / A r and CC-14/Ne mixtures without attachment 90 4.5 Electron relaxation in C C L i / A r with attachment 93 4.6 Electron relaxation in C C U / N e with attachment 94 4.7 The attachment heating effect for electrons in C C L i / A r and C C l / N e 4.8 The variation of 4 1/pTi.i . . 99 versus mole fraction CCI4 in C C ^ / N e mixtures at extremely small concentration 100 x 4.9 The variation of 1/pri.! versus mole fraction CCI4 for C C l / A r mixtures 4 and CCLj/Ne mixtures 101 4.10 The variation of the attachment rate coefficient versus time and the extraction of energy relaxation rate coefficients 103 4.11 Time dependence of the electron distribution function for electrons in a C C U / A r mixture 106 4.12 Electron attachment cross section, e-CCl Cross section used by Kowari 4 and Shizgal and Cross section used in the present work 107 4.13 Time variation of the average energy for CCI4 in Argon with all collision processes 108 4.14 Time variation of the average energy for CCI4 in Neon without vibrationally inelastic collisions Ill 4.15 The attachment heating and cooling effects for electrons in CCLj/Ne mixtures without vibrationally inelastic collisions 113 4.16 Steady electron distribution functions as the lowest eigenfunction of the Fokker-Planck operator; CCI4 in Neon without vibrationally inelastic collisions . . . 114 5.1 Differential cross section for (LI ,He) collisions with energy (eV) . . . . . 5.2 The isotropic kernel K~a(x, y) for two Henyey-Greenstein (HG) model cross + sections 5.3 130 Relaxation of the anisotropic ion distribution function for H G model cross section with e = 0.9 5.4 125 134 Relaxation of the anisotropic ion distribution function for H G model cross section with e = 0.1 135 xi 5.5 Time variation of the average ion energy and directed ion velocity with different H G model cross sections 5.6 137 Relaxation of the ion distribution for a hard sphere cross section with different E/n values 5.7 141 Time variation of the average reduced drift velocity with different E/n at various initial conditions 143 xii Acknowledgement I wish to express my sincere thanks to my research supervisor, Dr. B . Shizgal. It has been a pleasure to have worked with him, and his support, direction and assistance will always be greatly appreciated. Thanks are also due to the various members of Dr. Shizgal's research group for contributing to an enjoyable working environment. Special thanks go to Dr. Ken-ichi Kowari for the permission to use his program for the electron/CCU attachment problem. Finally, I wish to express my great appreciation to my parents and my wife Christina for all their patience and encouragement. This thesis is dedicated to them. xm Chapter 1 Introduction There are numerous physical situations in gaseous systems for which a theoretical development is based upon kinetic equations such as the Boltzmann equation (BE) and the Fokker-Planck equation (FPE) neutral and ionized systems [1] â€” [4]. These include transport properties of [5] â€” [10], reactive systems [11], plasmas [12], neutron thermalization [13,14], hot atom chemistry, aerosol dynamics [15,16] and space physics [17] â€” [20]. Accurate and efficient methods of solution for these kinetic equations are important in order to compare theory with experiment in a variety of different physical situations. The objective of this work is to develop new numerical methods for the solution of the B E and F P E for the velocity distribution functions ( V D F ) of different neutral or charged species. The macroscopic properties such as the drift velocity, diffusion coefficient and other transport properties are the appropriate averages of the velocity distribution. There are experimental techniques primarily Doppler spectroscopy, which have been used to probe the details of velocity distribution functions [21] â€” [25]. The distribution functions are probed by using single frequency laser-induced fluorescence to scan the Doppler profile at specific detection times relative to some initial pulse. The Doppler spectroscopy has been used for measuring the V D F for systems which are of practical importance for the study of the transport processes both in the laboratory [5,6], [26][28] and in the upper atmosphere [29]- [32]. For example, Bastian et al. measured the velocity distributions for B a ions drifting in A r under the influence of an external + 1 Chapter 1. Introduction 2 electric field [21]; Nan and Houston studied the relaxation of energetic sulphur atoms in He and A r gases [24]; and Park et al. studied the time evolution of the V D F of energetic hydrogen atoms in 0 , N and rare gases [25]. Doppler spectroscopic techniques are 2 2 useful in comparison with theoretical calculations [33, 34]. However both theoretical and experimental knowledge of the ion or neutral velocity distributions is still quite limited. The standard method for the solution of the kinetic equations is to expand the V D F in a suitable set of basis functions [1, 5, 6], [35] â€” [38]. The B E and/or F P E are reduced to ordinary differential equation and/or matrix equations for the expansion coefficients. The transport quantities such as, the drift velocity and diffusion coefficient, can be expressed in terms of the lower order coefficients. Consequently, traditional transport theory [1,5,6], [35] is not concerned with the calculation of the V D F . The V D F is required for any comparison with the above mentioned Doppler type measurements. The present thesis is directed towards the development of numerical methods that would provide accurate VDFs. The present work considers two major applications. The first considers the transport of electrons in atomic and molecular moderators. The second part considers the transport of ions in similar atomic moderators. Since the electron/moderator and ion/moderator mass ratios are very different, the elastic e-neutral and ion-neutral collision integrals in the B E are treated differently. This thesis provides a simple theoretical methodology for both systems. The systems studied in this thesis are spatially homogeneous and the V D F depends only upon the velocity v and the time t. 1.1 General methods for ion transport system To motivate the theme of this thesis, a classic kinetic theory problem is considered. The calculation of the mobility of ions drifting through a background neutral gas species Chapter 1. Introduction 3 under the influence of an external uniform electric field has been studied for over 50 years [1,5,6], [35] â€” [38]. There is a continuing interest in the calculation of the variation of the mobility with the density reduced electric field E/n [26,33,34]. The main objective of these workers was directed primarily towards the calculation of the ion mobility and not the ion V D F . However, during the past few years, the development of Doppler spectroscopy [21] â€” [23], which can probe the details of the V D F has stimulated the need for an efficient methodology for the calculation of ion distribution functions. In general, the ion V D F under the influence of an external electric field is very anisotropic and non-Maxwellian. The difficulty in the calculation of the anisotropic distribution functions from the B E is due to the large departures from Maxwellian. The kinetic theory that has been developed to date is based upon a moment method in which the ion distribution function is expanded about some zero order estimate which serves as the weight function for a basis set of orthogonal polynomials. The weight function used in the earliest studies of ion mobilities was a Gaussian, w(x) = e~ . x2 It is hoped that the weight function is close to the actual solution so that only a small number of additional terms have to be calculated from the B E . A thorough review of the details of this approach has been given by Mason and McDaniel [6]. The V D F of an ensemble of ions, dilutely dispersed in a background of neutral species, f(v,t) is given by the Boltzmann equation, where e is the ion charge, E is the electric field and m is the mass of the ion species. J is the collision operator which describes the effects of ion-neutral collisions on the velocity distribution function and is given by J[f] = 2 v r | J[ff[ M - ff Hg,x)gsm d dy , M 1 X X 1 (1.1.2) Chapter 1. Introduction 4 where / / ^ ( v i ) is the distribution function of the neutrals assumed to be Maxwellian, Â°"(<7) x) l s the differential elastic cross section for ion -neutral collisions, g = v â€” v i is the relative ion - neutral velocity, and x is the scattering angle in the center of mass frame. The quantities v and v i are the velocities before a collision and v ' and v i ' are the velocities after a collision. The system is assumed to be spatially uniform. Since we consider systems subjected only to an external uniform electric field E , the velocity v is specified by its magnitude v and orientation relative to the E field direction given by the angle 0. The traditional approach, known as the one temperature theory [5,9], is the expansion f{v,e,t)^Y Y, n{t)L ^{x )p (p), a l (1.1.3) 2 t l n Â£ where x = rnv /2kTb is the reduced energy, k is the Boltzmann constant, T& is the neutral 2 2 gas temperature and p = cos 9, L (% ) are associated Laguerre polynomials and Pt(p) 2 n 2 are Legendre polynomials. Laguerre polynomials satisfy an orthogonality relation with a Gaussian weight function fÂ°Â° e-^L^(x )L ^(x )x dx Jo 2 f + n 2 2(+1 2 =6 . nnl (1.1.4) The product of Legendre and Laguerre polynomials are referred to as Burnett functions. The choice of Legendre polynomials in the angle variable is made because the collision operator in the B E is diagonal with this basis set [39]. The choice of Laguerre polynomials is made for two reasons [6]. First, they are the basis functions orthogonal with the equilibrium Maxwellian distribution as the weight function. Second, they are the eigenfunctions of the collision operator in the B E for the pair potential, V(r) = Vo/r , and 4 so are thought to be useful whenever the collision frequency is slowly varying. For other potentials, there is no compelling reason why this traditional choice of basis functions is the best. The criteria for the choice of basis functions is the ease of their use in a specific Chapter 1. Introduction 5 application and the rate of convergence of the solution with the size of the basis set. However, it has been shown that for large E fields, the one temperature theory does not converge rapidly. For large fields, the anisotropy of the V D F is large and many terms in the expansion in Legendre polynomials are required. The choice of the Burnett functions as basis functions leads to results for the mobility which are a power series expansion in [5]. E/n The calculation of ion mobilities as a function of E/n improves by replacing Th with an ion basis temperature T in the definition of reduced energy, that is, x = 2 rnv /2k,T. 2 This approach is referred to as the two temperature theory [6]. With the proper choice of T , the two temperature theory gives better results for the calculation of ion mobilities than the one temperature theory, but still converges poorly at high E/n values. To permit the accurate calculation of ion mobilities at high E fields and to account for the large, anisotropy, Mason and Viehland [6] chose a drifting bi-Maxwellian weight function characterized by the temperature parameters Tji and Tj_, given by m â€¢WBM(vx,v ,v ) y z = exp ( l v + l) v rn(v - w) z 2kTÂ±_ 2kT\\ (1.1.5) where w, the unknown drift velocity and the temperature parameters Tj|, Tj. are to be determined. The basis functions for this weight function are Hermite polynomials in v ,v , x y and v . The orthogonality relation in each velocity direction is z / ' e-^H (i )H ,{( )di n where ^ = rnv /2kTÂ±. 2 x n x x = 6 ,, nn (1.1.6) This choice of basis set obscures the useful symmetry property that occurs with one temperature basis functions, where Legendre polynomials are used. For the bi-Maxwellian weight function, the calculation of the matrix elements of the collision operator becomes complicated . The use of a bi-Maxwellian weight function permits an accurate calculation of ion mobilities at higher E/n values where the two temperature Chapter 1. Introduction 6 theory fails. Since the ion speed is scaled with a third basis temperature parameter, the above method is referred as the three temperature theory [6]. A disadvantage of this approach is that the representation of the collision operator is cumbersome and the use of large basis sets becomes difficult. In addition, the three temperature method was designed to provide accurate values of the mobility but not the ion V D F . 1.2 General methods for electron transport system In this case, we are interested in the V D F of an ensemble of electrons dilutely dispersed in a background of atoms. This is very similar to the ion system of section 1.1. The important difference between the two systems is that for electrons the small electron/moderator mass ratio leads to some simplification. The theory has been considered by several researchers and developed via different methods. These include (i) the displaced pseudo Maxwellian approximation introduced by Mozumder [40]; (ii) the Monte Carlo simulations used by Koura [41]; (iii) a variety of techniques, such as moment methods, finite difference and Monte Carlo simulations used by Braglia [42] in order to solve the F P E ; (iv) the approach based on the Q D M for solving the F P E , developed by Shizgal and coworkers [43]- [45]. For a spatially uniform electron transport system, the V D F is given by a linear FokkerPlanck equation derived from the spatially independent B E [46] â€” [48]. Similar to the ion case, the orientation of the velocity vector v relative to the direction of the E field is given by angle 0. The electron V D F , / ( v , i ) , is expanded in Legendre polynomials P\(p)- With the substitution of this expansion into the F P E , an infinite set of coupled differential equations is obtained for the expansion coefficients f)(v, t) [39], [43] â€” [45], [48] as discussed in Chapter 3. In order to solve the infinite set of equations, the expansion is truncated at a finite Chapter 1. Introduction 7 number of terms. Owing to the small electron and rare gas moderator mass ratio, the randomization of the electron velocity direction will occur much faster than the rate of energy exchange. The electron V D F will therefore become isotropic in a very short time compared with the energy exchange. Consequently, only the two terms with / = 0 and 1 are retained, and the set of coupled equations reduces to the two equations for fo and f\. This procedure is referred to as a two-term approximation [39,43,48]. If inelastic collisions are important, then the anisotropy of the electron V D F is larger and more terms in the expansion of the distribution are required. The use of the two-term approximation has been shown by Lin et al. [49] to be excellent approximation for elastic electron-moderator collisions. Furthermore for a steady field, the relaxation of f\ to its steady value occurs on a time scale much shorter than the time scale for f to approach equilibrium. Consequently, we 0 set dfi/dt = 0 and express f\ in terms of fo. With the replacement of f\ for / , a single 0 partial differential equation for the isotropic component of the V D F is obtained. The steady solutions at infinite time f (v,oo) 0 can be calculated from this F P E by setting dfo/dt = 0 and are referred to as the steady electron velocity distribution function D(x). In the absence of electric field, it is a Maxwellian distribution, whereas for a steady d.c. electric field, it is a Davydov distribution [45]. This Davydov distribution can be easily calculated and involves simple integrals over the electron-moderator collision cross-section. It is useful to note that for the "ion-problem" discussed in Section 1.2, the two term approximation is useless since the anisotropics are large. 1.3 T h e m e and outline of the thesis The main objective of this thesis is the development of discretization procedures for the solution of kinetic theory problems. We consider the representation of the V D F by Chapter 1. Introduction 8 a set of grid points rather than by an expansion in a basis set. A numerical method known as the Quadrature Discretization Method (QDM) [50], which is based upon the representation of the V D F by its values at a set of quadrature points, is employed for solving the B E and F P E . The solutions are determined at a set of points which coincide with the points of a quadrature. In most cases, the quadrature is based upon a set of orthogonal polynomials with an arbitrary weight function. In all cases where the generalized orthogonal polynomials expansion has been applied, the success of a given expansion solution in evaluating the actual distribution function depends upon a good choice of the weight function. The optimum choice of the weight function is one which gives a good approximation to the actual distribution with a small number of polynomials. In general, the weight function that provides rapid convergence of the solution to the B E and F P E is the steady distribution at infinite time. The orthogonal polynomials for such nonclassical weight functions are constructed with the Gautschi discretized Stieljes procedure [51,52]. Chapter 2 considers the ion problem of Section 1.1 with a simple model for the collision operator developed by Bhatnagar, Gross, and Krook ( B G K ) [53]. This model provides an analytic solution for the ion V D F that allows us to study the convergence properties of the solution of the B E with different weight functions and different polynomial basis sets. The method introduced by Gautschi for the generation of a basis set, and the Q D M are discussed in this chapter as well. Chapter 3 outlines the Fokker-Planck theory and the two-term approximation for electron neutral systems. The eigenvalue spectrum of the F P operator and the time evolution of the electron V D F are studied by different numerical methods, including finite-difference (FD), Lagrange interpolation (LI) and the Q D M . Again, this chapter demonstrates the significant of the choice of weight function in obtaining faster convergence of the electron V D F . Chapter 4 considers the coupling between the electron thermalization and electron attachment in CCLj/Ar and C C U / N e Chapter 1. Introduction 9 mixtures. In these systems, the collisional processes that are included are, elastic e~~inert gas, and e - C C L ; collisions, e~-CCl vibrationally inelastic and e ~ - C C l attachment - 4 4 collisions. The e~-CCl4 inelastic collisions are described by nonlocal difference operators and a finite difference technique is used to solve the B E . A finite difference approach is also used for the time derivative and the resulting set of linear equations is solve with the successive over-relaxation (SOR) technique. The theoretical time dependence of the electron density and temperature, as well as electron relaxation times are compared with experimental data of Warman and Sauer [54] and Shimamori and Sunagawa [55]. Chapter 5 describes the relaxation of the ion V D F for ion neutral systems. A study of the relaxation of some initial anisotropic non-Maxwellian ion distribution for a model based on a realistic differential cross section is considered. The ion V D F is determined by the Q D M , which is essentially the one discussed in Chapters 2 and 3, employed there for differential operators and here for integral operators. In order to simplify the calculation, the differential cross section is fitted with a Henyey-Greenstein function [56] which duplicates the forward scattering very well. A study of the ion V D F under the influence of an applied electric field is also carried out. The ion V D F is written in terms of the longitudinal velocity component along the E direction and the velocity component perpendicular to the E direction. The collision operator of the B E is written as an integral operator with a well defined kernel, and the integrals involved are evaluated with a Gaussian quadrature rule. This is combined with a F D scheme for both time and speed derivatives so as to determine the time evolution of the ion V D F . Chapter 2 The Generation of Basis Functions for Kinetic Theory Problems 2.1 Introduction The V D F of an ion moving through a background of neutrals in the presence of an external uniform electric field is given by the solution of the B E . As mentioned in Section 1.1, the standard technique for the solution of this kinetic equation involves the expansion of the ion V D F in some suitable set of basis functions, which is constructed from a weight function, w(x) [5,9,35,36]. A Gaussian weight function is generally chosen and the polynomials orthogonal with this weight function are products of Laguerre and Legendre polynomials. Under the influence of an applied E field, the ion V D F , in addition to being highly anisotropic, decays very slowly as the velocity component in the direction of the applied electric field increases. The solution of the B E decays very slowly in the large speed region. Hence, the ion V D F defined with a Gaussian weight can diverge [57,58]. This has been discussed by Skullerud for the one temperature approach [57] and it is possible that this is also true in both the two and three temperature theories. In this chapter, we demonstrate the importance of the choice of basis set by considering a very simple model. We approximate the collision operator with the model developed by Bhatnagar, Gross, and Krook ( B G K ) [53]. The B G K model provides a useful analytic solution for the ion V D F in order to study the convergence properties of numerical solutions of the B E . We first study the convergence properties of the ion V D F with Hermite polynomials which are basis functions of a Gaussian weight. Then we study 10 Chapter 2. The Generation of Basis Functions for Kinetic Theory Problems 11 the convergence properties of the ion V D F with an alternate polynomial set defined by a weight function which is based upon the B G K solution. The choice of B G K solution, which includes the drift velocity, as the weight function is anticipated to give a better approximation to the ion velocity distribution function. The polynomials orthogonal with respect to B G K weight function are generated with a procedure developed by Gautschi [51,52]. The use of polynomials generated from various weight functions enhances the accuracy and flexibility of the expansion technique for kinetic theory. The organization of the chapter is as follows. Section 2.2 summarizes the derivation of the analytic solution of the B G K model first introduced by Whealton and Woo [59]. The B G K model, although not a realistic model, provides some insights into the nature of the ion V D F . Section 2.3 discusses the expansion of ion V D F in Hermite polynomials, that is, with a Gaussian weight function [60, 61]. A detailed discussion of the Gautschi procedure and the convergence properties of the solution of the B G K basis set are presented in Section 2.4. Section 2.5 outlines the details of Q D M numerical method. Lastly, Section 2.6 studies the time evolution of an ion V D F for the B G K model by using both finite difference (FD) and pseudospectral (PS) techniques. The difference between the two methods is that they employ different techniques in the treatment of the speed derivative. In addition, the steady ion V D F s obtained from these calculations are also compared to the analytic B G K solution. 2.2 Ion velocity distribution function for the B G K model The B G K model involves the approximation of the collision operator J with the simple form given by Chapter 2. The Generation of Basis Functions for Kinetic Theory Problems 12 The velocity independent collision frequency v represents the rate of collisions in some approximate way and f is the bi-Maxwellian distribution defined as M where vÂ± is the speed of the ion in a plane perpendicular to the applied uniform electric field E and v\\ is the speed of the ion along the E field direction. Eq. (2.2.1) was first introduced by Whealton and Woo [59] for equal ion and neutral masses. This was later modified by Makabe et al. [62] for unequal ion and neutral masses. The time and spatially independent B G K Boltzmann equation for equal ion and neutral masses is \v mv J ov\, 0 where v'^ = i>||/i>o a n d 'Â± = vÂ±/v 11 are the dimensionless speed, and v = \j2kTk/m is v 0 0 the thermal speed. Since there is no correlation occurring between the two velocity components t>|| and vÂ±, the ion V D F in a plane perpendicular to the electric field remains Maxwellian and unaffected by the drift. The effects of the uniform E field is only represented by the change of the ion V D F in the ||-direction. When we multiply Eq. (2.2.3) by v'j_ and integrate over v'^, the resulting equation is fb\df(v') U - s t with D = â€” , N = ' 1 ' W + = * ' " " â€¢ â€¢ < 2 - 2 - 4 > r- and we have used poo 2TT / where f(v' ,v'^) = f(v'^)f (v' ), M Â± L / M Â« K c f o l = 1, and f (v' ) M L = 'fcr ~"^â€¢ ^ e 27r fc (2.2.5) w e define a reduced drift velocity D' = D/v , and a dimensionless speed in the moving reference frame given by 0 c[| = â€” D', x = c||/D', we get ^ l ox + f( ) x = Ne- ' ^ . D 2 2 (2.2.6) Chapter 2. The Generation of Basis Functions for Kinetic Theory Problems 13 In order to solve the above differential equation, we multiply it with a factor e and x integrate over x. With the condition of f(x) = 0 at x = â€” oo, the general solution of Eq. (2.2.6) can be written as f(x) = Ne- [ x e - ' ^ dx. X x D 2 (2.2.7) 2 By completing square in the term (x â€” D' (l + x) ) in Eq. (2.2.7) and with the definitions 2 B = JJ577-I and A = ^ 7 â€” 2 D', we now have f(x) = Ne e~ B x f e- '^ -^dx. J â€” OO D x (2.2.1 With the transformation y = D'x â€” A , Eq. (2.2.8) can be reduced to B -x Ne /(*) e f D l x D' _ A J- e y dy. (2.2.9) In order to do the integral in Eq. (2.2.9), we first separate the integral into two different integrals (2.2.10) e dy + e dy D U-oo Jo The first integral in Eq. (2.2.10) is elementary and the second can be expressed in terms f( ) x = F v â€” / y y of the error function so that Q f(x) = N^ -^â€” 2D e [l - erf (A - D'x)], (2.2.11) where (2.2.12) erf(x) = â€”-â€” / e ' dt. 7T Jo With the definitions of the complementary error function erfc(s) = 1 â€” erf(x) and c[| = D'x, the final result for the normalized ion V D F in the ||-direction is / ( c ii ) = 2Vr C I | / Â°' e r f c ( A " ii c : (2.2.13) Chapter 2. The Generation of Basis Functions for Kinetic Theory Problems 14 Eq. (2.2.13) is the B G K analytic solution for the ion V D F along the applied field direction and is controlled by the reduced drift velocity D'. The B G K model provides a known exact solution which enables us to study the convergence and approximation of the orthogonal polynomial solution to the B E . In the following section, we will use the B G K solution to examine the convergence properties of Hermite polynomials at different D' values. As mentioned earlier, they are basis polynomials of a Gaussian weight function which generally serves as the zero order approximation for the Maxwellian distribution. 2.3 The expansion of B G K solution in Hermite polynomials A standard technique in kinetic theory is to represent the distribution function as a series expansion in suitable basis functions [5, 9,10] in the velocity space of the form N n where f\ is the N th order approximation to the velocity distribution function / , a n are expansion coefficients, and w is a weight function which can be considered as a zero" order approximation to the distribution function / . The zero t/l 1 order approximation can be obtained by solving some simple related problems, by making an educated guess based upon physical arguments, or even just an inspired guess. Once the weight function is chosen, a complete set of orthogonal polynomials (f> (v) are defined and satisfy an n orthogonality relation given by J w(v)4> (v)4> (v)dv = 6 , m where 6 mn n mn (2.3.2) is the Kronecker delta. Putting Eqs. (2.3.1) and (2.3.2) together, we obtain directly the expansion coefficients a in the form, n a = J f(v)(j) (v)dv. n n (2.3.3) Chapter 2. The Generation of Basis Functions for Kinetic Theory Problems Making use of Eqs. (2.3.1), (2.3.2) and (2.3.3), we realize that f 15 will converge towards A the actual distribution function / only if the integral J w f dv l 2 exists (2.3.4) 7^ oo [57,60]. Eq. (2.3.4) reduces the choice of the zero th order approximation and serves as a qualitative criterion for selecting it. Therefore, the choice of the weight function determines whether the series expansion will be convergent, and if convergent whether a good approximation can be obtained with a limited number of polynomials. To illustrate the convergence difficulties that can occur, we study the expansion of the B G K solution, Eq. (2.2.13) in a set of normalized Hermite polynomials H n [60], that is f ( 4 ) ~ e- il X> r7 ( f|), (2-3.5) s2c nJ n 5C where w = e~ U and s is a scaling factor for the reduced velocity variable analogous to s C the basis temperature parameter in the two and three temperature theories. Using Eq. (2.3.3), the expansion coefficients are given as -c'/D' evfc(A - C||)-# (sc||)d(.sc||), n (2.3.6) where Eq. (1.1.6) has been used. With an integration by parts and the use of the relation dH (x)ldx n = y/2nH -i(x)i we find n n = (D's) 1 7y^Tj7 / #n(-scj|)e ( c i l + Â£ ) ' dc'^ + \/2naâ€ž-i^ (2.3.7) If y = cj| + D', Eq. (2.3.7) can be written as / \ 1 1 fÂ°Â° i iâ€” The integral in Eq. (2.3.8) can be determined with Gauss-Hermite quadratures. {2.3.1 Chapter 2. The Generation of Basis Functions for Kinetic Theory Problems 16 Figs. 2.1 A and 2.IB show the variation of the partial sums of the expansion series of Eq. (2.3.5), for selected values of cj| versus the number of terms N in the expansion with two different values of D'. Shown in the figure is the ratio of the approximation in the RHS of Eq. (2.3.5) and the B G K analytic solution given in Eq. (2.2.13). It is clear from the results of Figs. 2.1 A and 2.IB that the series diverges. Increasing the reduced drift velocity D' tends to drive the B G K distribution towards the large c|| region and the distribution decays much slower than the Gaussian weight function. As a result, the integral in Eq. (2.3.4) fails to exist and the expansion series of Eq. (2.3.5) diverges. The present results agree with those of Skullerud [57] and Ferrari [63] who suggested that this should be a general behavior for the solution of the Boltzmann equation based upon a Gaussian weight function. Also shown in Figs. 2.1C and 2.ID is the behavior of the partial sums versus N for a weight function made wider by choosing s < 1. The series expansion might behave better if the weight functions were made wider since a larger portion of the actual distribution will be encompassed. However, this does not happen and is clearly seen from the figures. In Figs. 2.1C and 2.ID, the convergence can be improved but remains nonuniform in speed and diverges if a sufficient number of terms is retained. In Fig. 2.2, the exact B G K distribution, /(cj|) is shown (dashed curve) and compared with the result from Eq. (2.3.5) (solid curve) for fixed and different s. Provided that N is not too big, an optimum value of s can be chosen to give a reasonable agreement with the exact B G K solution. Widening the Hermite weight function with smaller s increases the best value of A^. This eventually leads to a somewhat better approximation, but does not solve the convergence problem. Thus, it is important to find another weight function which is close to the actual distribution that satisfies the convergence condition as well. In the next section, we will study the convergence properties of a set of basis functions by choosing the B G K analytic solution as the weight function. The major problem Chapter 2. The Generation of Basis Functions for Kinetic Theory Problems 17 0 20 10 N 30 12 16 N Figure 2.1: Convergence of the expansion of the B G K solution in Hermite Polynomials (A) D' = 0.25, cj, = 0.8, s = 0.943; (B) D' = 0.40, cj, = 0.5, s = 0.943; (C) D' = 0.60, cf = 0.8, s = 0.315; (D) D' = 0.60, cf, = 2.0, s = 0.315. Chapter 2. The Generation of Basis Functions for Kinetic Theory Problems 18 Figure 2.2: Convergence of the expansion of the B G K solution in Hermite polynomials; variation with the scale factor s: D' = 0.6, N = 50, s is equal to (A) 0.1, (B) 0.2, (C) 0.3 and (D) 0.315 ( B G K solution and Expansion solution). Chapter 2. The Generation of Basis Functions for Kinetic Theory Problems 19 associated with the selected B G K weight function is to generate a set of orthogonal polynomials. Fortunately, this is overcome with the use of a numerical procedure developed by Gautschi. 2.4 Gautschi Stieltjes procedure for the construction of non-classical basis functions The generation of a set of polynomials orthonormal with respect to some weight function w(x) has been discussed in many texts [64] and papers [50,65,66]. For a non-classical weight function such as the B G K distribution, the general approach of constructing the associated orthogonal polynomials basis set is based upon a three term recurrence relation [65,67]. For a set of polynomials {Qk(x)} with Qq = 1, and Q - i = 0 the three term recurrence relation is given by Qk+i{x) = (x - a )Q (x) k - /3 Qk-i{x). k k (2.4.1) The polynomials generated by Eq. (2.4.1) are not normalized. The corresponding normalized polynomials Pk(x) = Qk{ )/^hk x xP (x) = Jp\^P (x) k k+1 satisfy the recurrence relation + a P {x) + yJ%P ^(x), k k k (2.4.2) where the normalization factors are given by, h = J w(x)Q (x)dx. 2 k k For the weight function w(x) = x e~ , 2 relations for the coefficients a k and x2 (2.4.3) Shizgal [50,65,66] derived a set of recurrence These recurrence relations were numerically unstable and extended precision arithmetic was used. Clarke and Shizgal [68] studied analytic techniques to overcome the numerical instabilities. Alternatively, Gautschi's Chapter 2. The Generation of Basis Functions for Kinetic Theory Problems 20 Stieltjes procedure can be used to generate the cxk and fa coefficients. The recurrence coefficients in Eq. (2.4.1) are given by, and fa = With Q Q = 1, we can use Eqs. (2.4.4) and (2.4.5) to obtain a? and we set /3 = 0. 0 0 Applying Eq. (2.4.1), with k = 0 gives us Q i , and together with Qo, allows us to calculate ct\ and fa from Eqs. (2.4.4) and (2.4.5). Knowing ct\ and fa, we can use again Eq.(2.4.1), with k = 1 to determine Q , which together with Q\ gives a and fa from Eqs. (2.4.4) and 2 2 (2.4.5) etc. However, with this alternating Stieltjes procedure, the normalization factors hk are found to become increasingly small as k increases. Consequently, the calculation of the integrals in Eqs. (2.4.4) and (2.4.5) can fail due to the accumulation of round off errors [51,65,68,69]. Therefore, the associated orthogonal polynomials are difficult to construct because of the loss of accuracy. Gautschi's Stieltjes procedure involves the accurate calculation of Eqs. (2.4.4) and (2.4.5) by subdividing the domain of interest into many subdomains and evaluating the contribution from each subdomain with a high order quadrature. Gautschi discusses the use of Gauss-Legendre and Feijer quadrature rules, and in this way, all the recurrence coefficients o^, fa and orthogonal polynomials Pk{x) can be calculated. Another numerical method to construct orthogonal polynomials for a certain class of weight functions was provided by Danloy [70]. Danloy's method requires the reduction of integrals in Eqs. (2.4.4) and (2.4.5) to forms that permit their exact evaluation with known quadratures. This is only possible with certain weight functions. We use Danloy's Chapter 2. The Generation of Basis Functions for Kinetic Theory Problems 21 procedure, which is exact, to test the accuracy of Gautschi's method which is numerical. As a result, the coefficients a , fa and Qk(x) are determined through the Stieltjes k procedure by a three term recurrence relation. Danloy's method exactly computes the normalizations on the interval [a,6], h = / w(x)Ql(x)dx. Ja (2.4.6) k We are now interested in constructing the polynomials orthogonal with respect to the weight function equal to the B G K solution of the ion problem, that is given by Eq. (2.2.13). Thus if x = cf^/D', we introduce the B G K analytic solution as w(x) â€” ^-[e~ eric(A x - D'x)] into Eq. (2.4.6) and let x = cj|//_)', we then have pB roo h =â€” e- edc{A - D'x)Q (x)dx, 2 Jâ€”oo x (2.4.7) 2 k and with the definition of the error function erf(x), Eq. (2.4.7) can be rewritten as h =(-=) e~ Vy7T/ J-oo J-oo e~ Q (x)dydx. x y k (2.4.8) 2 k With an elementary transformation, Eq. (2.4.8) yields ( ^ ) /-~ = h ~ e X Q l { x ) d X -' (2 4 9) - If we let x = z + ^nr, we obtai am D' II By completing square in e~^ {^y-^-w) ^ y r )d e ~ Z Q l { z ^ - + v )dz ( 2 - 4 J 0 ) ~D ' and with the definition of t = (y + jjjr) , we get y +M hk= ' e {y2+ r ) 2 f ^ d t J Â°Â°e-Ql(z+ ~^ t o + i )dz. A (2.4.11) Eq. (2.4.11) involves a t-integration with weight function e ' , and a ^-integration with - 2 weight function e~ . Hence hk can be determined exactly with Gauss-Hermite and Gaussz Laguerre quadratures, since the integrand is a polynomial of finite degree. If we apply M Chapter 2. The Generation of Basis Functions for Kinetic Theory Problems 22 Gauss-Hermite points tj and weights vfp, and M Gauss-Laguerre points Zi and weights (*) to these integrals, we can rewrite Eq. (2.4.11) in the form of discrete sums of positive quantities, (2.4.12) As both Gauss-Hermite and Gauss-Laguerre quadratures are known, we can use Eq. (2.4.12) to obtain accurate numerical results for hk with the application of these Gaussian quadrature rules. Similarly, we can also get accurate a k and (3 values with the above k procedure. The recurrence coefficient a of the B G K weight function from Eq. (2.2.13) k at different D' values are determined by both Gautschi and Danloy methods. Results obtained with the two methods agree very well. Tables 2.1 and 2.2 show the convergence of a k calculated from the Stieltjes proce- dure developed by Gautschi at various set of intervals (NINT) with different number of quadratures (NP) at each interval for two different D' values. The results in Tables 2.1 (D' = 1) and 2.2 (D' = 3) show that the Gautschi's method provides convergence for a^o up to 12 significant figures with about 3000 and 4000 quadrature points, respectively. In Tables 2.3 and 2.4, different orders of a are calculated from Danloy's method at different k M (or NP) values. The results in Table 2.3 and 2.4 show that Danloy's method only requires about 2500 quadrature points to obtain the same kind of accuracy provided by Gautschi method for a 50 at two different D' values. Apparently, Danloy's method is a better approach because it calculates the integrals in Eqs. (2.4.4) and (2.4.5) exactly. However, unlike the Gautschi procedure, Danloy's method only applies to certain weight functions. The Stieltjes procedure developed by Gautschi works very well as long as a suitable number of intervals and the number of quadratures in each interval are chosen. Once a , f3 and orthogonal polynomials Q {x) are determined, we can study the k k k convergence properties of different basis functions. Chapter 2. The Generation of Basis Functions for Kinetic Theory Problems NP "20 "10 ^30 a 40 23 Â«50 NINT = 10 40 80 19.0617299873 17.7217300809 36 .2554698794 36 .8920237438 160 17.7215980596 17.7215980604 17.7215980604 240 320 84 .2633453871 36 .9032746136 62.5735518539 56.1400578847 56.3308057551 75. .8762724767 96.0311801782 94.1488732114 95.4932161492 36 .9032745655 36 .9032745655 56.3308052313 56.3308052313 75 .8762714465 75, .8762714465 95.4932154027 95.4932154027 78, .4902833734 75, .8769499168 75. .8762714447 75. .8762714465 75, .8762714465 100.7775740118 95.4929467558 95.4932154376 95.4932154027 95.4932154027 75. .4132972179 NINT = 20 20 40 17.6396442331 17.7216040721 37. .5882236120 36. .9029808777 54.2811644182 56.3310510424 80 120 160 17.7215980604 17.7215980604 36. .9032745656 36, .9032745655 36. .9032745655 56.3308052303 56.3308052313 56.3308052313 NINT = 40 10 20 40 60 80 16.9686353663 17.7216916603 17.7215980692 17.7215980604 17.7215980604 39, .2284549699 36, .9054123821 55.7642613262 56.3192737267 74, .7763171614 75, .8572879458 91.9129964908 95.5527404516 36, .9032742423 36, ,9032745655 36, .9032745655 56.3308073579 75, .8762786630 56.3308052313 â€¢ 75,,8762714465 56.3308052313 75, .8762714465 95.4932167379 95.4932154026 95.4932154027 NINT = 80 5 18.8206313603 36. ,0578063491 52.8041521410 76. ,7706550337 92.8757274616 10 17.7182624160 17.7215993441 36. ,9212299396 36. ,9032557564 56.3750994916 75, ,8039235730 95.2602788422 36, ,9032745662 56.3307690123 56.3308052313 75. ,8763159143 75. ,8762714060 95.4933597269 95.4932152013 40 17.7215980603 17.7215980604 36, ,9032745655 56.3308052313 75, .8762714465 50 17.7215980604 36. .9032745655 56.3308052313 75. .8762714465 95.4932154027 95.4932154027 20 30 Table 2.1: Comparison of the convergence of ak with Gautschi Procedure with weight function, w(x) = ^-[e~ er{c(A â€” D'x)] where x = cUD' and D' = 1. x Chapter 2. The Generation of Basis Functions for Kinetic Theory Problems NP Â«10 Â«20 Â«30 24 Â«40 NINT = 10 40 80 160 240 320 400 480 55 9673124102 58 3185681997 58 4380134395 58 4380203142 58 4380203283 58 4380203283 126 7567945475 118 9929772156 117 8449552491 117 8449851025 117 8449849259 117 8449849258 117 8449849259 167.2887790643 176.0167149040 177.4342304191 177.4342336438 177.4342335767 177.4342335780 177.4342335780 257 2182913803 240 4006993513 237 1104153198 237 1096887356 237 1096898804 237 1096898771 237 1096898771 323 7652176000 292 9028115171 296 8341061426 296 8374700349 296 8374679621 296 8374679686 296 8374679686 253 3578647033 239 4648832496 237 0996426620 237 1096798288 237 1096902996 237 1096898771 237 1096898771 293 3879345320 300 6626159328 296 8682519833 296 8374746497 296 8374667967 296 8374679686 296 8374679686 242 7247111181 235 3324636915 237 1044085963 237 1096432274 237 1096902522 237 1096898767 237 1096898771 237 1096898771 300 299 296 296 296 296 296 296 2977224689 1049686090 8515435791 8375657911 8374666273 8374679704 8374679685 8374679685 245 2960170641 237 9596957329 237 0541548134 237 1096998527 237 1096971286 237 1096897950 237 1096898772 237 1096898771 237 1096898771 311 300 296 296 296 296 296 296 296 8839480405 5275197324 9944856309 8374465837 8374494427 8374682542 8374679677 8374679685 8374679685 NINT = 20 20 40 80 120 160 240 280 64 2302754349 58 5543029931 58 4379941037 58 4380199748 58 4380203296 58 4380203283 58 4380203283 126 2738244097 116 4023290944 117 8421476041 117 8449842918 117 8449849386 117 8449849258 117 8449849258 190.7797216213 180.7031762677 177.4427783951 177.4342389863 177.4342334598 177.4342335780 177.4342335780 NINT = 40 10 20 40 60 80 100 120 140 59 8626678367 58 6868582785 58 4379414495 58 4380205121 58 4380203284 58 4380203283 58 4380203283 121 3430434513 116 9143294247 117 8449069994 117 8449803284 117 8449849336 117 8449849258 117 8449849258 177.5965335054 179.4528793494 177.4354014843 177.4342509407 177.4342335019 177.4342335780 177.4342335780 NINT = 80 5 10 20 30 40 50 60 70 80 61 0394442631 58 9387884686 58 4379805739 58 4380285834 58 4380203304 58 4380203281 58 4380203283 58 4380203283 124 0784118829 118 8582090897 117 8407231169 117 8449622327 117 8449855760 117 8449849246 117 8449849258 117 8449849258 174.0288987309 178.9208522358 177.4527737307 177.4342369026 177.4342309038 177.4342335952 177.4342335779 177.4342335780 177.4342335780 Table 2.2: Comparison of the convergence of a with Gautschi Procedure with weight function, w(x) = ^-[e 'erfc(/4 â€” D'x)] where x = cf^/D' and D' = 3. k -x Chapter 2. The Generation of Basis Functions for Kinetic Theory Problems NP aio 8 9 10 11 12 15 18 19 20 21 22 25 28 29 30 31 32 35 38 39 40 41 42 45 48 49 50 51 52 12.4694098442 8.8160419155 10.7075494377 17.7215980604 17.7215980604 Â«20 15.6149462720 20.6081722870 34.9302544554 17.3300426017 27.4683433515 36.9032745655 36.9032745655 Â«40 0:30 36.6200998624 39.9522627074 55.8167913273 26.4196629359 45.7867337115 56.3308052313 56.3308052313 25 Â«50 50.6957080108 59.3122582486 73.5084103965 36.9383337780 64.8437346262 75.8762714465 75.8762714465 77.9419102966 87.81935045186 48.82540989136 84.3111664044 95.4932154026 95.4932154026 Table 2.3: Comparison of the convergence of cx with Danloy's Procedure with weight function, w(x) = ^-[e~ evfc(A â€” D'x)] where x = cj|/D' and D' = 1. k x Chapter 2. The Generation of Basis Functions for Kinetic Theory Problems NP Â«10 8 9 10 11 12 15 18 19 20 21 22 25 28 29 30 31 32 35 38 39 40 41 42 45 48 49 50 51 52 33.9476503884 45.1698092913 9.7021016138 "20 Â«30 a 40 26 "50 58.4380203283 58.4380203283 60.3023397311 90.2530155902 99.5023907415 27.96034764289 117.8449849258 117.8449849258 99.5638142705 120.4371532797 145.2981528741 150.1401712948 49.9556444731 177.4342335780 177.4342335780 144.5418063253 145.6396387102 181.4659877401 199.9120480977 196.7300661542 74.5589346091 237.1096898771 237.1096898771 193.7698832080 203.1484304226 242.7475295431 255.0301105417 239.1545229850 101.1964767242 296.8374679686 296.8374679686 Table 2.4: Comparison of the convergence of a with Danloy's Procedure with weight function, w(x) = ^-[e~ evic(A â€” D'x)] where x â€” cUD' and D' = 3. k x Chapter 2. The Generation of Basis Functions for Kinetic Theory Problems 27 To study the convergence properties of basis functions with respect to the B G K weight function, we expand the distribution from Eq. (2.2.13) in this set of basis functions as follows, /(*, D') Â» w (x, BGK D'J Â£ a Pi \x, D'J, BGK n (2.4.13) n=0 D'J â€” ^Â§r[e~ ' ' erfc(A where wbgk{ , x P( \x, BGK x D w â€” x)] is the B G K weight function, x = c\\ and D'J are orthogonal polynomials with respect to weight function wbgk B G K weight function, wbgk â€¢ The has the same form as in the B G K distribution / in Eq. (2.2.13), with the exception that D' was replaced by a new parameter D' . Further, the w expansion coefficients in Eq. (2.4.13) are given by oo f(x,D')P^ \x,D'Jdx â€¢ (2.4.14) BGK / -oo Making use of Eqs. (2.4.13) and (2.4.14), we can determine the N th order approximation to / . The integral in Eq. (2.4.14) can be evaluated by quadratures based upon the wbgk [50, 64]. Table 2.5 shows the convergence of the expansion f(x, D') with D' in the range from 0.4 to 1.2 in B G K polynomials p( ) BGK with D' = 0.6. The solutions of the w B G K series expansion at both x values converge remarkably well except with D' = 1.2. For D' < D' , the distribution vanishes in the large speed region decays faster than w the weight function. The approximate solution to the distribution therefore remains finite and converges reasonably rapidly. However, when D' > D' , the approximate w solution converges only when the condition in Eq. (2.3.4) is satisfied. For D' = 1.2, the distribution vanishes in the large speed region decays much slower than the weight function and therefore the approximate solution to the distribution diverges. In Sections 2.3 and 2.4, we have demonstrated the significance of the choice of the weight function in a given application. It appears that the capability of generating nonclassical orthogonal polynomials with arbitrary weight functions increases the accuracy and flexibility of the expansion technique in kinetic theory. Chapter 2. The Generation of Basis Functions for Kinetic Theory Problems N D' = 0.40 0.50 0.70 0.80 1.0 1.2 0..3437 0 .3440 0 .3434 0..3433 0..3434 0.2801 0.2968 0.2997 0.2991 0.2991 0..2023 0,.2390 0,.2933 0,.2975 0..2914 0..2875 0..2411 0,.2863 0,.0281 0..0278 0..0281 0..0282 0..0282 0.0416 0.0356 0.0370 0.0370 0.0370 0..0581 o.,0376 0.,0420 0.,0510 0.,0480 0.,0394 0.,0415 28 x = 0.33 2 5 10 15 20 30 36 40 0,.4286 0..4403 0,.4442 0..4449 0..4453 0..4451 0..4450 0 .4459 0..4675 0..4763 0 .4799 0 .4804 0 .4801 0 .4800 0.3703 0.3684 0.3677 0.3677 x = 2.55 2 5 10 15 20 24 30 o..0101 0..0081 0..0073 0..0074 0..0074 0..0135 0..0126 0,.0122 0..0122 0.0225 0.0228 0.0230 0.0230 â€¢ Table 2.5: Convergence of the expansion of the B G K solution with different D' values in N B G K polynomials with weight function, w(x) = fjy-[e~ / â„¢erfc(A â€” x)] where x = cj| and D' = 0.6. x D w Chapter 2. The Generation of Basis Functions for Kinetic Theory Problems 29 In Section 2.6, we describe a numerical scheme for the solution of the time dependent B E equation with a B G K collision operator in Eq. (2.2.1). The numerical method will include the Q D M and the non-classical B G K weight function wbgk- The Q D M has been described by Shizgal and Blackmore [50,65,66] and used with great success for electron transport [43] â€” [45], optically bistability [71], and quantum mechanics [72] . A n outline of the Q D M is provided in the following section. 2.5 Quadrature Discretization Method (QDM) The Q D M is based upon the discretization of the distribution function on a grid of points rather than the expansion in a set of basis functions [50]. The main theme of the Q D M is that it employs non-classical weight functions (such as the B G K weight function ) and the associated polynomials in the discretization procedure. We consider a set of orthogonal polynomials, {R (x)}, defined by a weight function n w(x) on the interval [a, b], such that, / w(x)R (x)R (x)dx Ja n Associated with the roots, m of the N ih =S . (2.5.1) nm order polynomial, R^(x), / w(x)f(x)dx is a quadrature rule, Â« ^tu,/(i ), *=i a (2.5.2) t where W{ are the weights. The calculation of the points, and weights, W{, is described in previous papers by Shizgal et al. [50,65,66], and in standard references discrete representation of a function f(x) points, that is, f(xi). [64]. The is the values of the function at the discrete This representation is entirely equivalent to the polynomial basis representation which is given in terms of the coefficients a , in the polynomial expansion n N-l f{x) = J2 "nRn(x)n=0 (2.5.3) Chapter 2. The Generation of Basis Functions for Kinetic Theory Problems 30 The expansion coefficients a in Eq. (2.5.3) are given by n a = / w(x)R (x)f(x)dx, Ja n (2.5.4) n and Eq. (2.5.4) can be evaluated numerically with the quadrature rule, Eq. (2.5.2), that is, N a = ^2w R (x )f(x,). n i n If we substitute Eq. (2.5.5) into Eq. (2.5.3), fix) the N values f(xi), (2.5.5) i can be calculated at any point from that is, N f(x) = Y;ii(x)f(x ). (2.5.6) i The interpolating polynomial, T(x) is given by N-l Ii(x) = Wi Rn(Xi)Rn(x), (2.5.7) n=0 and satisfies the completeness relation (2.5.1 Ii{ j) â€” x The derivative of the interpolation operator provides the algorithm for numerical differentiation [73], defined by, df dx N X=X t The matrix D,- w. /wjl'j(xi) (2.5.9) Wi. J=l is the derivative matrix, and N-l Dij = y/WiWj (2.5.10) Rn{ i) n{x ). X R 3 n=0 where R' is the first order derivative to the n th n order polynomial R n with respect to x. The derivative operator is represented as a finite matrix of dimension equal to the number of grid points. The algorithms, Eqs. (2.5.6) and (2.5.7), can be used to reduce a general Chapter 2. The Generation of Basis Functions for Kinetic Theory Problems 31 differential operator equation to a set of coupled algebraic equations in the function evaluated at the quadrature points. This approach is referred to as a pseudospectral method (PS) whereas the expansion in a basis set is termed a spectral method. In the next section, we will study the time dependent B G K Boltzmann equation with the application of F D and PS method. The time evolution profiles of the ion velocity distribution will be presented and the steady solutions at infinite time will also be compared to the B G K solution in Eq. (2.2.13). 2.6 Ion velocity distribution function for the time dependent B G K model The objective of this section is to develop a numerical method for the solution of the transient V D F using a nonclassical polynomial set. The thermalization of ions is studied because there has been growing interest in the determination of the V D F both experimentally [21] â€” [25] and theoretically [33,34]. We employ here the B G K model, with and x = cf, = v' â€” D'. The B G K time dependent B E for the definition of D' = u VQvm II ' II equal ion and neutral masses is then ^ = /"(,)-/(,,,), + (2.6.1) where f(x,t') is the time dependent ion V D F along the applied electric field, f (x) M -(x+D') ^ 2 e _ ^ - jv g n e â€” dimensionless time and v is the time and speed independent collision frequency. One of the simplest methods for solving Eq. (2.6.1) is the Euler explicit method. The scheme employs a forward difference for the time derivative and backward difference for the speed derivative. Applying this method to Eq. (2.6.1) gives Ii Ji i r y J i Ji â€” 1 7 ^ + D' ' â€” = f M" At' A J 1 A where t' = n'At', n X i = iAx and f? = f(x ,t' ). t n t The values of ft, (2.6.2) of the ion V D F on the speed grid at the time t' + A t ' are calculated with an iterative procedure. Chapter 2. The Generation of Basis Functions for Kinetic Theory Problems A n alternative approach is a pseudospectral method 32 [74] which applies a forward difference for the time derivative and a pseudospectral method for the speed derivative. The pseudospectral method we employed here is the Q D M discussed in Section 2.5. For the Q D M , the algorithm for differentiation of the time dependent distribution in the discrete velocity space is written as ^ f e Q = ^;b^(x i > 0, . (2-6.3) 3 where f(x,t') and g{x,f) = w (x)g(x,t') BGK operator DÂ»j for g(x,t') The derivative = J2n=o a (t)P^ \x). BGK n is derived from Eq. (2.5.10), that is D ^ ^ ^ P i ^ ^ x . O P ^ ^ x O , (2.6.4) n=0 where Wj and Xj are quadrature weights and points associated with B G K weight function We now first substitute WBGR'(x)g(x,t') wbgk{x). for f(x,t') into the space derivative of Eq. (2.6.1), and employ a forward difference for the time derivative. We have f n + (x)-r( ) i x + D' WBGi<(x)g n At' where w' (x), (x) + g (x)w' (x) = f n BGK 0 ) - P(x), (2.6.5) the first derivative of wbgk with respect to x is given by BGK BGK\ ) W 2 -(A-x)* /5F e = X 2D' erfc(A-x) D' (2.6.6) w and g , the first derivative of g with respect to x is determined with D ; J in Eq. (2.6.3). n n The final equation of the B G K time-dependent B E with the pseudospectral scheme in discrete space is f n + \ x ) - r ( t At' X i ) W GI<(xi) B E ^i39 ( j) n X + g (xi)w'g (xi n = f M n (x,)-P( X i ). GK (2.6.7) W i t h a given initial condition, we can evaluate the RHS and the second term of the LHS in Eq. (2.6.7). Then we multiply these two terms by A t ' and add to the corresponding Chapter 2. The Generation of Basis Functions for Kinetic Theory Problems 33 distribution array element. The result is a new distribution function at time t' + A t ' . Repeating the above procedure with the newly acquired distribution, we can obtain the distribution function at any time t'. f (xi) n+1 However, for the PS method, the calculation of is accurate only if a small number of terms is needed in the expansion. In order to satisfy the above condition, we select the wbgk as the initial condition and choose the applied drift D' close to D' . We now compare the results of the ion V D F obtained with w the F D algorithm, Eq. (2.6.2) and the PS algorithm, Eq. (2.6.7). The time evolution of f(x,t') at D' = 1.0 and 2.0 calculated with the F D and PS methods are shown in Figs. 2.3A and 2.3B. In Fig. 2.3A, we choose the B G K weight function wbgk with D' â€” 0.6 as the initial condition. As time proceeds, the ion V D F w decays much slower in the cj| direction than the initial condition. For longer time, it reaches a steady state which matches with the analytic B G K solution in Eq. (2.2.13). Results from both the F D and PS methods agree very well. In Fig 2.3B, we choose another B G K weight function wbgk with D' w = 1.5 as the initial condition. Similar relaxation behavior of the ion V D F is observed. Again, results from the F D and PS methods agree very well and the steady solutions at long times agree with the B G K analytic solution. Although there are limitations in the solution of ion V D F for the PS method, the studies of the time-dependent B G K model act as a guideline to the methodology of solving the B E which has a more complicated collision operator. By using this method, we can calculate the relaxation of the ion V D F under the influence of an applied electric field as in Ch. 5. Chapter 2. The Generation of Basis Functions for Kinetic Theory Problems -4 -2 0 2 4 34 6 Figure 2.3: Time evolution of the B G K ion V D F with (A) D' = 1.0 and D' = 0.6 and with (B) D' = 2.0 and D' = 1.5. For both (A) and (B), t' is equal to (a) 0.0, (b) 1.5 and (c) 6.6. (The solid lines are the distribution function calculated with the F D method with 100 points where A t ' and Ac[| = 0.03. The triangles and the solid dots are the distribution function calculated with the pseudospectral method with 50 Gaussian-BGK points where A t ' = 0.03). The weight function wbgk(c\\) â€” [ e V ' e r f c ( A â€” cj|)] is chosen as the initial condition. w w - Â£> 0 Chapter 3 Numerical Methods of Solution of the F P E for Electron Thermalization 3.1 Introduction The study of electron thermalization and transport in both atomic [55], [75]- [82] and molecular moderators [83]- [92] continues to be an active field of study. The subject has important applications to radiation chemistry and physics, discharge devices, ionospheric applications, plasma processing of materials and plasma chemistry. The thermalization of energetic electrons in atomic and/or molecular moderators proceeds owing to collisions between the electrons and the constituents of the moderator assumed to be present in large excess. Electron-electron collisions are assumed not to occur and are not included in the theoretical analysis. A complete treatment would require the inclusion of elastic and inelastic collisions, and chemically reactive processes such as ionization and attachment. The main interest is to determine the time dependence of the electron energy distribution function starting from some initial arbitrary distribution function. The time scale for the relaxation process and the nature of the time dependent nonequilibrium distribution function are of importance. The theoretical approach is based upon the B E which for electrons, with small electron to moderator mass ratio, can be approximated by a F P E . Also, the anisotropy of the electron distribution function is expected to be small and two terms in the Legendre expansion of the distribution function is sufficient. A standard procedure for the definition of the electron relaxation time is to determine the time required to reach 10% or 1% of the equilibrium energy. Detailed reviews of the current 35 Chapter 3. Numerical Methods of Solution of the FPE for Electron Thermalization 36 status of the field were presented by Braglia [42] and by Shizgal et al [93]. This chapter considers a detailed comparison of several different techniques in the study of electron thermalization in inert gases for low initial energies on the order of several electron volts. Inelastic collisions are not included and only electron-atom elastic collisions are required. The initial energy loss by inelastic collisions is extremely fast and it is the much slower energy exchange by elastic collisions that determines the measured relaxation time. Although these atomic systems are not complicated, there is still considerable improvement to be made in the agreement between calculations and measurements [55]. Also, several interesting phenomena have been predicted theoretically. Two unexpected effects have been reported in the course of these studies. One is the "negative mobility effect" for which the transient electron mobility can be negative for some initial transient, predicted first by Braglia and Ferrari [94] and later by Shizgal and McMahon [44], and also confirmed with a Monte Carlo simulation by Koura [41]. Experimental observations of the effect were reported by Warman et al. [95]. Another effect is the "negative differential conductivity effect" which refers to the decrease in the electron drift velocity with increasing electric field strength over a particular field range. It was long recognized as arising from inelastic collisions and therefore should only occur in molecular systems. Calculations of the drift velocity of electrons in He-Xe and Ffe-Kr mixtures by Shizgal [96] demonstrated that this negative differential conductivity effect can occur for particular concentrations of these mixtures. The effect was later confirmed by Nagpal and Garscadden [97] unaware of the earlier work. There appears to be no experimental confirmation of this effect to date. There have been numerous theoretical approaches to the problem of electron relaxation in the inert gases. The basis for the theoretical approach is the B E with realistic data for the electron-atom momentum transfer cross section. The Boltzmann equation is Chapter 3. Numerical Methods of Solution of the FPE for Electron Thermalization 37 approximated by a Fokker-Planck equation and only two terms in the Legendre expansion of the angular part of the distribution function are included. The earliest treatment was carried out by Mozumder [40], and Tembe and Mozumder [98] who assumed that the electron distribution function remains a local Maxwellian with a time dependent temperature. They derived expressions for the rate of change of the electron temperature from the B E but the approach did not involve a solution of the B E . Knierem et al. [99] adapted a moment method of solution of the B E used previously in the study of ion transport. Because of their choice of Burnett functions as basis functions for the expansion of the distribution function, their solutions did not converge well. Koura [100] developed a Monte-Carlo simulation to determine the electron velocity distribution function during the relaxation. Because of the small mass ratio, the computations were very time consuming. Braglia and co-workers [42] employed a variety of techniques including moment methods, finite difference schemes and Monte-Carlo simulations to solve the F P E . Shizgal and co-workers [44,45], [101]- [103] introduced the Q D M for the solution of the F P E . The Q D M was developed by Shizgal [65] for the solution of the integral B E and Shizgal and Blackmore [104] for the solution of differential equations and in the study of nonequilibrium systems with bistable states [66]. It has been used extensively in the study of electron relaxation problems [101]- [103], in the solution of the Schroedinger equation [72] and the Navier-Stokes equation [105]. The Q D M is a discretization method based on the evaluation of the solution at a set of quadrature points defined by the roots of the N th order polynomial of a suitable set orthogonal with respect to some weight function in a specified interval. The main feature of the method is the use of nonclassical polynomials based on a weight function which closely approximates the solution and hence provides for rapid convergence. In this sense it differs from the usual spectral methods [74,106] which generally use Fourier series or Chebyshev polynomials. Chapter 3. Numerical Methods of Solution of the FPE for Electron Thermalization 38 There have been several numerical methods proposed for the solution of such FokkerPlanck equations. The most widely used being a finite difference method developed some time ago by Chang and Cooper [107] and recently refined by Larsen et al. [108] and Epperlein [109]. Robson et al. [110] and Robson and Pritz [111] proposed a derivative matrix technique, based upon a Lagrange interpolation, for the solution of a differential equation. A very useful discussion of the numerical solution of Fokker-Planck equations was recently published by Park and Petrosian [112]. The chapter is thus directed towards a detailed comparison of several different numerical methods in studying electron thermalization in Ar. In view of recent experiments, there is a need to reexamine the relaxation of electrons in argon and study the dependence on the form and energy of the initial distribution function. A comparison of several different methods of solution of the relevant F P E is carried out. In particular, the Q D M is employed with a quadrature based on the steady Davydov distribution function as weight function. The Davydov distribution is characterized by the electron-atom momentum transfer cross section and the electric field strength. The Gautschi Stieltjes procedure [51], described in Section 2.4, is used to generate the orthogonal polynomials and the associated quadrature weights and points. We compare this new Q D M quadrature with the speed polynomials, defined on [0, oo] with weight function w(x) = x exp(â€”x ) used 2 2 in most of the previous works. We also consider solutions obtained with a finite difference technique and the differential quadrature based on Lagrange interpolation, discussed by Robson et al. [110]. These authors suggested that their technique based on a Lagrange interpolation is superior to the Q D M . A comparison of the different discretization schemes is considered in the first instance with regard to the eigenvalue problem, (3.1.1) Chapter 3. Numerical Methods of Solution of the FPE for Electron Thermalization 39 where L is the Fokker-Planck (FP) operator defined explicitly later, and A and <f> are n n the eigenvalues and eigenfunctions, respectively. We then consider the time evolution of f(x,t) with some initial condition, / ( x , 0 ) , with the discretization of the time derivative and advancing the solution in time with either an implicit or explicit scheme. The previous work by Shizgal and co-workers [45,101-103] considered the time dependent solution expressed in the eigenfunctions of the F P operator. A study is carried out of the stability of the discretized spatial operator in the different schemes and related to the eigenvalue spectrum. Section 3.2 defines the different physical or model systems considered. The different discretization schemes of the F P operator are outlined in Section 3.3. The time evolution of the solution obtained with several different time discretizations is presented in Section 3.4 for the spatial discretizations in Section 3.3. A stability analysis is discussed in Section 3.5 based upon the eigenvalue spectrum obtained in Section 3.3. 3.2 The Lorentz-Fokker-Planck equation; electron thermalization The motivation for the present study arises from the work of Shizgal and co-workers on the transient behavior of the distribution function for electrons dilutely dispersed in a large excess of an inert gas [42,44,45]. As mentioned in the Section 3.1, there is a continuing interest in the calculation and measurement of electron relaxation times. Also, these systems provide useful benchmarks for different numerical methods. We introduce a new discretization which provides for rapid convergence for these systems. The time dependent, anisotropic, spatially uniform electron V D F is expanded in Legendre polynomials, that is, oo (3.2.1) where 6 is the angle between v and the polar axis chosen in the direction of the electric Chapter 3. Numerical Methods of Solution of the FPE for Electron Thermalization 40 field. For inert gases at low electron energies, for which only elastic collisions need to be included, only the terms in I = 0 and Â£ =1 are retained. With the use of Eq. (3.2.1) we get the usual two-term approximation [113], df eE d 2 ^at + T (-K+ )^ 6m ov v 0 L 1 m d = -m\v ^ i -Ov L kT d +â€” rnvi~) ov Jfo, b i (3-2.2) df\ eE df ^ 7 + â€” â€” = -vh. Ot m ov , 0 (3.2.3) In E q . (3.2.2), m\ is the moderator mass, E is the electric field strength, k is the Boltzmann constant, Tj, is the temperature of the moderator, and the collision frequency v(v) = nva(v), where cr(y) is the momentum transfer cross section and n is the number density of the moderator. As discussed in a previous paper [45], we set dj\/dt = 0 and substitute fx from Eq. (3.2.3) into Eq. (3.2.2). We define a dimensionless time t' = t / r where 1 nm 2kTh - = 7â€”<r \ r lm\ V "2 with CTQ some convenient hard sphere cross section, and define a = ctJctq. 3.2.4 0 introduce the definition of the dimensionless speed variable x = yrnv /2kT, 2 temperature T differs from the bath temperature Tf,. The quantity s 2 We also where the = T/Tb is a parameter used to scale the points in the Q D M . In this way, we get the equation for fo given by, fo s d dF d r 2 = T^[ 2x x dfo] 2 4 a f o + S ^ B { x ) ^ zoom ( 3 - 2 5 ) where B(x) = xa(x) + ^T %. xa(x) 7 (3.2.6) In Eq. (3.2.6), the quantity a is a field-strength parameter given by 2 M f eE \ 6m \nkTi,aoJ 2 , â€ž â€ž â€žx Chapter 3. Numerical Methods of Solution of the FPE for Electron Thermalization In Eq. (3.2.7), the electric field occurs as the density reduced field, E/n, 10~ 17 41 in units of V c m which is a Townsend (Td). From Eq. (3.2.5), the steady solution attained 2 at f = oo is given by ,A2; fo(x,oo) = D(x) = Cexp -2s L 2 f Jo (3.2.8) ^f-^dx' rf(x') where C is a normalization constant. The function D(x) is the steady electron distribution at equilibrium referred to as the Davydov distribution. In the absence of an electric field, B(x) is equal to xcr and D{x) is a Maxwellian. If we write f (x,t') 0 = D(x)g(x,t'), then the F P E is of the form [44], dg(x,t') (3.2.9) -Lg(x,t') dt' w here If d d 2 (3.2.10) dx' with A(x) = 2B(x) 2s x a( 2 2 B'(x). (3.2.1i; The F P E for this problem can also be rewritten in the usual form (3.2.12) or where J(x,t') x D(x). 2 = -[A(x)P(x,t') dP(x,t') dJ_ dt' dx + dB(x)P(x,t')]/sdx) If we set P(x, t ) = x D(x)g(x, 1 2 dg(x,t') dt' 1 sx D(x) 2 (3.2.13) is flux and vanishes for P(x,t') = t'), then the F P E is of the form, d dx x D(x)B(x 2 dg(x,t') dx (3.2.14) Chapter 3. Numerical Methods of Solution of the FPE for Electron Thermalization 42 The method of solution employed in the previous work [42,44,45,42] involved the formal solution of Eq. (3.2.9) for g(x,t') as given by, g(x,t') = e- 'g(x,0). Lt (3.2.15) If the initial distribution g(x,0) is expanded in the eigenfunctions of L, that is, oo d( ,0) = zZ n<t>n{x), x a (3.2.16) n=0 then the time evolution is given by, oo g(x,t') = Za <f> (x)e- ', Xnt z n n (3.2.17) 71 = 0 where a are determined from the initial condition and n L(f> (x) = \ <l> (x). n n n (3.2.18) Average transport properties can be evaluated from these expansions. For example, the electron energy relative to the thermal energy at equilibrium is given by, oo E(f) Eh = E^" ' t 3.3 3.3.1 U Tl=0 (3.2.19) Numerical Solution of the Fokker-Planck Eigenvalue Problem Quadrature Discretization Method (QDM) The solution of differential equations by differential quadrature methods is based upon the representation of the derivative operator, d/dx, in a discrete basis. The method has been discussed in Section 2.5 [44,45,66,72]. One important aspect of this discretization procedure is to preserve the self-adjoint property of the F P operator, that is, the discretized representation of the F P operator should be symmetric. Another important consideration is the choice of mesh points (or equivalently basis functions) determined by the quadrature weight function, w(x). Chapter 3. Numerical Methods of Solution of the FPE for Electron Thermalization 43 In this section, we present an outline of the derivation of the Q D M representation of the Fokker-Planck operator. The F P operator, L defined in Eq. (3.2.9) is self-adjoint with respect to the equilibrium distribution, Pq(x) = x D(x) 2 polynomial set {R (x)} as weight function. If the is orthogonal with respect to weight function w(x), we define the n set of polynomials {Q {x)} n orthogonal with respect to Po(x) as weight function, that is, w(x) Qn{x) = Po(x) 1 R (x) (3.3.1) n which satisfies the orthogonality relation / 1/2 roo Po(x)Q (x)Q {x)dx n =6 . m (3.3.2). nm JO The matrix elements of the operator L in this basis set are given by, roo Lnm = / Po(x)Q (x)LQ (x)dx, Jo n (3.3.3) m which with an integration by parts becomes, L = - nm Jo P {x)B(x)Q {x)Q (x)dx. 0 n (3.3.4) m The,boundary conditions which correspond to zero flux at both boundaries, so that, x D(x)B(x)â€”â€”â€”\ = x D(x)B(x)â€”â€”â€”\ . x=0 (3.3.5) x=00 However, as mentioned earlier, L is self adjoint with respect to the equilibrium function Po(x), and not the weight function w(x). [P (x)/w(x)] L[w(x)/P (x)] 2 can be constructed by evaluating the matrix representa- 2 0 A symmetric representation of L = 0 tion in the basis set of functions with respect to w(x), and transforming to the discrete representation with the appropriate unitary transformation. In terms of the polynomials {R (x)}, Eq. (3.3.1), and weight function w(x) we have n with Eq. (3.3.4), Lnm = ~ j roo ^ ( ^dx~ B W X d d + Kx)}Rn{x)[â€” + h(x)]R (x) m (3.3.6) Chapter 3. Numerical Methods of Solution of the FPE for Electron Thermalization 44 where the function h(x) is , . : w'(x) \x D(x)Y 2 which arises owing to the choice of a weight function w(x) differing from the equilibrium distribution Pq(x). If we now use the quadrature rule for the weight w(x), we have that, N = - E B(x )w [R (x ) + h(x )R (x )][R! (x ) + h(x )R (x )]. (3.3.8) k=l The derivatives are evaluated with the Q D M algorithm Eq. (2.5.10) which gives L nm L k k m k k N N = ~lZ\/^Rrn,{Xj)lZ^JyTiRn{Xi)YZ j=l i=l m k n k k n k N (3.3.9) B mn fc=l We then transform the F P operator in Eq. (3.3.9) from the basis set representation to the discrete representation of the function by using the unitary matrix T, with elements Tin = ^/w~iR {xi). n Finally, we have the Q D M representative in the discrete basis given by, Lij = - N B{x ) [D,- + hiS ] [Dj + hjSjk]. fc=i k fc ik (3.3.10) k Since Lij given by Eq. (3.3.10) is clearly symmetric, the eigenvalues A are real. The n details of the derivations in this section have appeared in previous papers [45,72,104,66]. 3.3.2 Lagrange Interpolation In a recent paper, Robson et al. [110] suggested that the use of a Lagrange interpolation algorithm with a uniform grid in place of an interpolation based on a set of orthogonal polynomials would be more efficient. Their result, previously obtained by other workers, [73] is based upon the Lagrange interpolation polynomial given by, * w = ffi(Â»~ 11j^i=\\ i x * r ^j) ( 3 - 3 1 1 ) Chapter 3. Numerical Methods of Solution of the FPE for Electron Thermalization 45 which also satisfies Eq. (2.5.8). The points are chosen anywhere on the interval of interest, and define the roots of the polynomial N-l P {x) N where is the coefficient of x . J J (x - xj), = k N (3.3.12) It therefore follows that N *(*)=, P n [x [ X } (3.3.13) Xijr^yxi) p ( V and the derivative matrix is Dij =Â£-(x ). (3.3.14) J We use the form given by Eqs. (3.3.13) and (3.3.14), and write the representation of the F P operator in this form, + sffisQjfoQDffl [x B(x )D( )]% 2 r t L l3 Xi - 3 -jjj-- â€”, . (3.3.15) which is not symmetric as is the Q D M representation, Eq. (3.3.10). Robson et al. [110] suggested that this is preferred since the evaluation of the derivative matrix is faster with Eqs. (3.3.13) and (3.3.14) rather than with Eq. (2.5.10), and the free choice of the grid points is advantageous. They evaluated the accuracy of the Q D M and Lagrange differentiation algorithms by considering the calculation of the second derivative of some arbitrary function. In a subsequent paper, Robson and Pritz [111] applied their approach to the solution of several applied problems. 3.3.3 Finite difference scheme The self-adjoint form of the F P operator was shown to be consistent with zero flux at the boundaries, Eq. (3.3.5). This boundary condition is also related to particle conservation dt - D(x)g(x, t')x dx = x D(x)B(x) 2 2 ^ } |~ 0 = 0. (3.3.16) Chapter 3. Numerical Methods of Solution of the FPE for Electron Thermalization 46 A n y useful discretization would have to ensure particle conservation which yields, Ao = 0. W e discretize the speed variable according to 0 = x\ < x < x 2 Â£,â€¢+1 = X{ + Ax, and A x = x /(N â€” 1) where x max 3 < x^ â€” x with max is the speed point chosen large max enough so that the flux boundary condition is satisfied. We also introduce a shifted grid defined b y Xi -[/ = x; + Ax/2. W i t h a centered difference for the derivative, + 2 dg(x,t') _ g(x i,t') l dx -g(xj,t') i+ k + 1 ' ~ Ax 2 (ooi \ 7 ' [ 6 - 6 A i ) the finite difference representation of the eigenvalue problem is, N 2)2Lij(f>n(xj) - \ <l> (x ) 3=1 n n (3.3.18) t where = [ A 1 , (Ax) xfDiBi + \ y xj,,Di B , +1 1 2 x 2 2 + 1 / + A+i/ 2 1 2 2 x | _l + 2 + 1 1 A =- , , * T ^ (Ax) 2 i+1 Â» = 1, A, 3.3.19 2 x DB (Ax) x H^ + x 2 + 1 / 2 / 2 + i 2 A+i/ i+ ^ \ A+i/2 2 i = l, N - l , (3.3.21) w i t h the understanding that the first term i n the fraction on the R H S of E q . (3.3.19) vanishes for i = 1 and the second term vanishes for i = N i n order to enforce the boundary conditions. 3.4 Time Dependent Solutions of the Fokker-Planck Equation T h e previous work on the thermalization of electrons employed the eigenfunction expansion solution of the F P E . One of the difficulties w i t h this approach is that the calculation converges poorly for very small times and large initial energies for an initial delta function distribution. T h e reason for this is that the expansion coefficients i n E q . (3.2.16) required the calculation of the eigenfunctions at large argument and these are not accurately calculated. T h i s aspect of the eigenfunction expansion was discussed at length Chapter 3. Numerical Methods of Solution of the FPE for Electron Thermalization 47 in previous papers. In this section we compare the eigenfunction approach with a time discretization using the F D method was well as the Q D M . 3.4.1 Finite Difference Time Discretization Following Larsen et al [108], we discretize the velocity variable as described earlier in Section 3.3.3. With t' = nAt' for the discrete time variable using a backward Euler difference for the time derivative, that is, dg(x,t>) _g? -g? +l dt' where gâ„¢ = g(xi i/ ,nAt'), + dg(x,t') _ g f 8V 2 + 1 At' ' { the F D version for the time dependent F P E becomes - g? _ xl,D B (g^ z+1 - g f ) - x D B^ + 1 l+1 2 +1 t At' (Ax)2x 2 + 1 / 2 A+i/2 - g^ ' 1 ' ' j which after rearranging terms can be written as E i=i = Hi (3-4-3) where V is a tridiagonal matrix whose elements are given by V = l + AtL u u = AtLi^i Vi, i+1 = AtL i = l,.....N, (3.4.4) i = 2, N, (3.4.5) N - 1. (3.4.6) i = 1, hl+1 The new values of the distribution, c7â„¢ , i = 1, ...N are then obtained from Eq. (3.4.3) +1 by inverting V . Notice that the matrix V is diagonally dominant i.e. \Vu\ > IV^-I] + lVi_i4.il and that the off-diagonal elements are negative; therefore, if gf > 0 for all i , then gâ„¢ > 0 for all i and the positivity of the distribution is preserved in time as mentioned +1 by Larsen et al. [108]. Chapter 3. Numerical Methods of Solution of the FPE for Electron Thermalization 48 The choice of a backward Euler difference for the time derivative leading to an implicit scheme is dictated by stability arguments. For an explicit scheme, it can be easily seen that the matrix elements are not uniformly bounded functions of x near x = 0 when an electric field is present (a ^ 0). We also observed that for an explicit scheme, the zero eigenvalue is inaccurate and hence an explicit scheme would lead to unphysical behavior and the relaxation to equilibrium is not guaranteed. 3.4.2 Q D M and Lagrange Interpolation As with the F D scheme, the solution of Eqs. (3.2.9) and (3.2.10) was a first order Euler method in the time is used. The discretized version of Eq. (3.2.9) becomes, 9 t A / ' =Uigt\ (3-4-7) where Lij is the matrix representation of the Fokker-Planck operator. This results in the following implicit scheme for the time evolution of the distribution function, given by, ^^D^-At'^-]- ^. 1 â€¢ (3-4.8) i=l The discretized form of Lij for the two cases are as given previously. 3.5 3.5.1 Results and discussion Calculation of eigenvalues and eigenfunctions In this section, we have carried out a comparison of the accuracy of several different discretization procedures in the evaluation of the eigenvalues and eigenfunctions of the F P operator defined by Eqs. (3.2.10) and (3.2.18). These include the Q D M with speed and Davydov quadratures, F D scheme, and LI procedure with uniform grids and grids based on Legendre polynomials. The eigenvalues and eigenfunctions are determined by Chapter 3. Numerical Methods of Solution of the FPE for Electron Thermalization 49 the numerical diagonalization of the F P operator. In the latter case, the representation of the F P operator is not symmetric and the eigenvalues can be complex. Both F D and LI schemes require specifications of some maximum speed x max so as to define the right hand side boundary. Tables 3.1 to 3.8 show the convergence of the lower order eigenvalues with different discretization methods for the F P operator defined for collisions of electrons with either helium, argon or xenon. The momentum transfer cross sections are those reported by Nesbet [114] for helium and by Mozumder [40] for argon and xenon. The cross section for helium is approximately independent of the electron energy, while argon and xenon have cross sections that vary rapidly with electron energy and have a distinctive minimum at low energy. The minimum in the collision cross section for the heavier gases (Ar and Xe) is known as the Ramsauer Townsend minimum. As mentioned earlier, in the absence of electric field as shown in Tables 3.1 and 3.2, the Davydov weight function is the speed weight function and there is no difference between the Davydov and speed quadrature points. Table 3.1 shows the convergence of the eigenvalues for electrons in helium with the F D scheme, Q D M with speed points, and LI procedure with uniform and Legendre grids. It is clear from results in the table that the Q D M with speed points gives the most .rapid convergence. It is marginally more rapid than the LI procedure with Legendre quadrature points. The other methods give good results with larger number of points. Table 3.2 shows similar results for electrons in argon. The rate of convergence is much slower than in Table 3.1 because of the strong energy variation of the argon cross section and the Ramsauer Townsend minimum that occurs. These results show, as in Table 3.1, that the Q D M with speed points gives the most rapid convergence, the F D approach gives relatively poor results. Results with LI uniform grid procedure could not be obtained for argon due to spurious oscillations of the Chapter 3. Numerical Methods of Solution of the FPE for Electron Thermalization N 10 20 30 50 80 100 5 10 15 20 25 30 35 5 10 15 20 25 30 5 8 10 12 15 20 25 Ai A 2 A 3 FD 78.130 82.682 89.155 90.497 90.569 90.550 A 5 24.643 25.645 25.639 25.613 25.599 25.596 43.770 54.600 55.388 55.444 55.398 55.380 21.213 25.562 25.586 25.590 25.592 25.607 25.672 LI (uniform g rids) 60.437 170.05 54.400 86.648 210.00 55.332 90.309 174.18 55.341 90.473 175.32 55.322 90.541 175.54 55.226 90.726 175.01 54.881 92.436 162.37 23.252 25.593 25.589 25.589 25.595 25.589 25.589 249.26 151.24 155.41 171.82 175.10 175.43 LI (Legendre Â£frids) 62.513 171.63 55.598 87.369 189.93 55.341 90.464 176.67 55.342 90.466 175.49 55.342 90.465 175.49 Q D M (speed) 55.978 99.925 55.345 90.668 193.80 55.342 90.473 178.67 55.342 90.466 175.82 90.465 175.50 175.49 A 50 7 684.93 262.12 246.95 254.74 273.40 276.14 468.43 303.16 278.22 277.74 277.18 251.72 468.36 258.06 278.22 277.90 277.90 439.57 329.35 293.14 279.12 277.90 277.90 Table 3.1: Comparison of the convergence of eigenvalues (in units of r ) with Q D M , LI, and F D approaches for electrons in helium (E/n = 0.0 Td). a Chapter 3. Numerical Methods of Solution of the FPE for Electron Thermalization 51 N. A 2 A A 3 FD 8.4724 10.142 10.333 10.413 10.438 10.444 A 5 7 ( a ) 10 20 30 50 80 100 4.4437 3.9691 3.8970 3.8622 3.8507 3.8481 7.4067 6.9381 6.9367 6.9347 6.9339 6.9337 LL 10.247 10.104 10.459 10.456 10.453 10.453 22 438 17 898 19 547 20 140 20 325 20 366 35 166 27 338 30 559 32 540 33 130 33 262 25 582 17 336 19 810 20 459 20 441 20 436 20 436 53 700 38 754 28 511 32 931 33 536 33 473 33 469 33 470 {b) 20 30 40 50 60 70 80 100 3.6771 3.8543 3.8442 3.8435 3.8434 7.5207 6.9791 6.9287 6.9334 6.9334 Q D M (speed) (=) 10 20 30 40 5.0 7.1889 3.8576 3.8438 3.8437 14.929 6.9760 6.9336 6.9335 29.098 10.550 10.454 10.454 81.894 20.326 20.437 20.436 184.54 40.282 33.614 33.473 33.470 (a) F D with x = 5.0. (6) LI with Legendre grids and x = 8.0. (c) Speed with scale factor, s = 0.336. max max 2 Table 3.2: Comparison of the convergence of eigenvalues (in units of r ) with Q D M , LI, and F D approaches for electrons in argon (E/n â€” 0.0 Td). _ 1 Chapter 3. Numerical Methods of Solution of the FPE for Electron Thermalization 52 interpolating polynomials near the edges, as discussed by Fornberg in another application [73]. The interpolation error Ek \f {x ) exact k - f interpolation^^ i .calculated for the s LI scheme with both an equi-spaced grid and Legendre quadrature points. The test function / used here is the second eigenfunction of the F P operator in Eq. (3.3.10) for electrons in argon with no applied electric field. For equi-spaced grid, the interpolation error, which is shown in in Fig. 3.1A, increases with increasing N. The approximate function does not converge. By contrast, the interpolation error for Legendre quadrature points is shown in Fig. 3.IB and provides a converged accurate interpolation versus N. This agrees with the work by Fornberg that an equi-spaced grid yields oscillations in the interpolating polynomials near the edges, and considerable roundoff error as N increases. On the contrary, Legendre quadrature points seem to suppress these oscillations and ensures converged value with larger N as clearly shown in Fig. 3.1. We have also considered a similar study for the F P operator with two electric field strengths 0.25 T d and 0.5 T d . The density reduced electric field values, E/n, Townsends (Td) where T d = 1 0 - 1 7 Vcm . 2 are in Again, it was not possible to obtain useful results from the LI uniform grid procedure with an applied electric field due to roundoff error. Tables 3.3 and 3.4 illustrate convergence of the lower order eigenvalues for electrons in helium with an applied electric field. The results in Table 3.3 show that the Q D M with Davydov points provides the convergence of the first nine eigenvalues up to 5 significant figures within 25 points. It is slightly faster than the Q D M with speed and the Lf scheme with Legendre quadrature points, whereas the F D method gives rather poor convergence. W i t h a higher electric field strength E/n = 0.5 T d as shown in Table 3.4, similar results to Table 3.3 are observed. The convergence of the eigenvalues for electrons in argon with E/n = 0.25 T d and 0.5 T d is shown in Tables 3.5 and 3.6 respectively. With the strong energy variation argon cross section, we notice from the results in both Tables 3.5 and 3.6 that the Q D M with Davydov points gives more rapid convergence than the Q D M with Chapter 3. Numerical Methods of Solution of the FPE for Electron Thermalization 53 i- o fcl so o L_J â€¢16 0 10 20 30 L_J I . I 40 50 60 70 40 50 60 70 N o e o 0 10 20 30 N Figure 3.1: Interpolation error Â£ \ f { x ) - f interpolated ^ ^ us N with (A) equi-spaced grid (B) Legendre quadrature points. f (x) is the second eigenfunction of the Fokker-Planck operator using the argon cross section at E/n â€” 0. e x a c t f c k x for t exact h e L I s c n e m e v e r S Chapter 3. Numerical Methods of Solution of the FPE for Electron N Ai A 2 A A 3 10 30 50 70 100 133.77 136.35 136.45 136.48 136.49 284.62 308.66 309.40 309.58 309.67 FD 414.73 502.19 504.58 505.13 505.40 5 10 15 20 30 131.77 136.72 136.50 347.03 310.64 309.76 309.75 LI<Â» 776.29 480.85 505.32 505.64 A 5 Thermalization 54 9 ( a ) 1238.8 927.72 938.69 941.05 942.19 22983 1849.9 1918.4 1932.2 1938.8 1628.2 949.17 943.18 943.21 2036.6 1900.6 1943.3 Q D M (speed) 332.10 981.58 309.81 508.44 1043.3 309.77 505.68 944.19 505.67 943.26 943.25 (c) 5 10 15 20 30 5 10 15 20 25 138.15 136.51 Q D M (Davyd ov) 311.98 535.42 309.77 505.68 944.38 505.67 943.25 136.55 136.51 89621 2332.8 1961.2 1943.3 52260 1966.8 1943.4 1943.3 (a) F D with x = 8.1. (6) LI with Legendre grids and maximum x (c) Speed with scale factor, s = 2.0. max max = 9.0. 2 Table 3.3: Comparison of the convergence of eigenvalues (in units of r and F D approaches for electrons in helium (E/n = 0.25 Td). _ 1 ) with Q D M , LI, Chapter 3. Numerical Methods N of Solution Ai A 2 of the FPE for Electron A A 3 10 30 50 70 100 200.02 202.95 203.06 203.08 203.09 427.30 458.13 458.85 459.01 459.11 FD 659.90 744.38 746.78 747.28 747.53 10 15 20 30 40 203.18 203.11 459.85 459.18 459.16 750.01 747.73 747.75 A 5 Thermalization 55 9 ( a ) 1452.0 1374.3 1386.2 1389.4 1389.7 1244.7 1376.0 1390.1 1390.3 Q D M (speed) 459.60 774.13 2248.0 459.18 747.78 1395.1 747.77 1390.4 44911 2734.6 2815.8 2836.2 2838.1 2537.0 2513.9 2837.0 2837.5 (c) 10 15 20 30 5 10 15 20 203.13 203.11 Q D M (Davyd ov) 461.29 822.97 459.18 747.78 1394.5 747.77 1390.3 1390.4 203.14 203.11 5056.0 2885.0 2837.6 93053 2847.6 2837.6 (a) F D with x = 10.1. (b) LI with Legendre grids and x = 11.0. (c) Speed with scale factor, s â€” 2.0. max max 2 Table 3.4: Comparison of the convergence of eigenvalues (in units of r ) with Q D M , Lf, and F D approaches for electrons in helium (E/n = 0.50 Td). 1 Chapter 3. Numerical Methods of Solution of the FPE for Electron Thermalization 56 N Ai A 2 10 30 50 70 100 290.74 295.24 295.35 295.37 295.38 684.07 761.75 763.02 763.26 763.37 15 20 25 30 40 295.53 295.40 295.40 763.20 763.45 763.45 A A 3 FD>) 1306.2 1331.8 1337.9 1339.0 1339.4 1328.6 1339.1 1339.8 1339.8 A 5 4237.0 2658.7 2707.1 2714.4 2717.4 9 5424.8 5978.4 6055.6 6091.9 3035.5 2748.9 5887.3 2719.6 6476.3 2719.6 '6098.5 6094.2 Q D M (speed) 766.30 1419.9 4286.9 763.56 1341.5 2809.5 763.51 1340.0 2730.1 763.50 1339.9 2719.8 2719.8 (c) 20 25 30 40 50 60 5 10 15 20 30 295.60 295.43 295.41 Q D M (Davyd ov) 767.20 6392.7 763.50 1340.4 2863.3 1339.9 2720.6 2719.8 295.84 295.41 7119.7 6149.9 6095.6 6094.5 59671 10633 6123.6 6094.5 (a) F D with x = 13.5. (b) LI with Legendre grids and x = 14.0. (c) Speed with scale factor, s = 2.65. max max 2 Table 3.5: Comparison of the convergence of eigenvalues (in units of r and F D approaches for electrons in argon (E/n = 0.25 Td). _ 1 ) with Q D M , LI, Chapter 3. Numerical Methods of Solution of the FPE for Electron Thermalization 57 N Ai A 2 A A 3 10 30 50 70 100 495.99 504.82 505.01 505.06 505.08 1184.8 1308.0 1309.9 1310.3 1310.5 FD 2296.5 2291.3 2300.0 2301.5 2302.1 15 20 30 40 50 508.12 505.10 1402.3 1310.5 1310.6 2186.4 2302.4 2302.5 A 5 9 ( a ) 20 30 35 40 50 60 5 10 20 30 6393.1 4583.3 4652.2 4662.0 4665.4 4676.7 4667.7 4667.5 Q D M (speed) W 2268.8 4962.9 1311.4 2304.0 5021.4 1310.7 2302.8 4686.5 2320.7 4668.0 4667.9 735.22 505.21 505.13 Q D M (Davydov) 1317.4 1310.7 2303.2 4867.0 2302.7 4667.8 4667.9 505.44 505.13 9404.9 10230 10341 10388 10812 10283 10380 10379 .21776 12702 10749 10383 10378 10432 10379 (a) F D with x = 15.8. (b) LI with Legendre grids and x = 17.0. (c) Speed with scale factor, s = 2.65. max max 2 Table 3.6: Comparison of the convergence of eigenvalues (in units of r ) with Q D M , LI, and F D approaches for electrons in argon (E/n = 0.50 Td). 1 Chapter 3. Numerical Methods of Solution of the FPE for Electron Thermalization 58 speed and the LI with Legendre quadrature points. This can also be seen for the case of electrons in xenon shown in Tables 3.7 and 3.8. Again, the F D method gives rather poor convergence in all cases. The results in these tables illustrate that the Q D M quadratures based upon the Davydov weight function converge with the smallest number of points among the four discretization methods, although with a sufficient number of quadrature points, all approaches give converged results. The Davydov points are comparable and superior to other methods, especially at higher electric fields and for cross sections that vary strongly with energy. This is reasonable since the quadrature points determined with the Davydov distribution as the weight function depend upon the cross section and the electric field strength. The Davydov distribution is the steady state distribution when an electric field is imposed. The electric field heats up the electrons and drives them to higher energy regions. The effect of the applied electric field strength on the Davydov distributions for two different collision cross sections helium and argon is shown in Fig. 3.2. We also show the eigenfunctions i/> and ipr for electrons in helium and electrons 3 in argon with an applied electric field. The peculiar shape of the eigenfunctions can be understood by considering the potential that occurs in the Schroedinger equation corresponding to the F P E . With the definition of y = JQ r~T a n < ^ ^{x) = e ^^ (a;), c n we reduce the F P operator into a Schroedinger equation [44,45,93,102]of the form (3.5.1) where (3.5.2) and (3.5.3) Chapter 3. Numerical Methods of Solution of the FPE for Electron Thermalization 59 iV Ai A 2 A A 3 A 5 9 FJj(a) 10 30 50 70 100 747.21 737.91 737.23 737.03 736.93 1036.2 1034.8 1035.2 1035.3 1035.4 1814.6 2213.0 2221.3 2223.2 2224.2 3810.1 4226.9 4285.3 4298.8 4305.6 664550 9923.8 10559 10704 10782 15 20 30 40 50 774.07 746.93 736.82 736.85 1052.6 1035.9 1035.4 2197.6 2224.3 2225.1 2225.0 4148.1 4308.7 4311.6 4311.6 11436 10490 10823 10809 10810 10 20 30 35 40 50 60 70 Q D M (speed; 2345.1 3668.6 9103.1 1040.4 2263.1 4426.1 1035.5 2225.7 4324.7 1035.5 2225.3 4314.9 2225.2 4312.2 4311.9 4311.8 891.68 740.13 736.90 736.89 QDM 1035.3 1035.5 1035.5 10 738.22 â€¢ 15 736,96 20 736.89 25 (Davydov) 2247.8 4266.7 2225.1 4315.6 2225.2 4312.0 4311.8 79266 14297 11027 10902 10854 10815 10811 10810 20562 11429 10800 10810 (a) F D with x â€” 10.7. (b) LI with Legendre grids and x â€” 12.1. (c) Speed with scale factor, s = 2.65. max max 2 Table 3.7: Comparison of the convergence of eigenvalues (in units of r ) with Q D M , LI, and F D approaches for electrons in xenon (E/n = 0.25 Td). 1 Chapter 3. Numerical Methods of Solution of the FPE for Electron Thermalization 60 N Ai A 2 A A 3 A 5 9 F_p_(Â«) 10 30 50 70 100 1374.6 1403.5 1404.4 1404.6 1404.7 3592.2 3619.4 3615.5 3614.1 3613.4 4400.4 4619.4 4629.6 4631.2 4631.9 8208.1 10675 10726 .10717 10707 20 30 40 50 60 1402.8 1404.9 1404.8 3615.9 3613.2 3612.7 LTW 4625.6 4632.6 4632.5 10197 10699 10695 20 30 35 40 50 60 70 10 20 30 35 Q D M (speed) 1422.5. 3627.1 4693.4 1405.1 3614.2 4637.8 1404.9 3613.2 4633.8 3613.0 4632.9 4632.8 1406.5 1404.9 11137 10722 10705 10698 10696 Q D M (Davydov) 3620.9 4674.3 10505 3613.0 4632.9 10696 4632.8 10696 20081 20881 20990 21052 17867 20360 20668 21016 21018 36966 21641 21170 20930 21024 21019 21018 58687 21098 21019 21018 (a) F D with x = 12.8. (b) LI with Legendre grids and x = 13.6. (c) Speed with scale factor, s = 3.0. max max 2 Table 3.8: Comparison of the convergence of eigenvalues (in units of r *) with Q D M , L I , and F D approaches for electrons in xenon (E/n = 0.50 Td). Chapter 3. Numerical Methods of Solution of the FPE for Electron Thermalization 61 Figure 3.2: Davydov distributions for electrons in helium and argon at E/n equal to (a) 0.0, (b) 0.5 and (c) 1.0 T d ( 1 0 - V c m ) . 17 2 Chapter 3. Numerical Methods of Solution of the FPE for Electron Thermalization 62 Fig. 3.3 shows the convergence of i/> for electrons in helium versus the number of 3 points N used for each method with E/nâ€”0.25. The solid curve calculated with 70 Davydov points is considered to be the exact solution. It can be seen that the Q D M with speed or Davydov points and LI scheme give converged results with just N=10, whereas F D scheme does not. The results for a higher electric field E/nâ€”0.5 T d are shown in Fig. 3.4. In this case, the Q D M with speed points and LI results are not completely converged at 10 points. The increase in electric field strength increases the electron average energy and shifts the eigenfunctions to higher energy as well. Consequently, it requires more points to give converged results. Similar results are obtained for ip shown in Figs. 3.5 7 and 3.6. A l l methods have converged results by N=30. The eigenfunctions for electrons in argon with E/n=0.25 and 0.5 T d are also determined. In each plot the solid curves are the exact results obtained with 70 Davydov points to ensure convergence. Figs. 3.7 and 3.8 show the convergence behavior with F D scheme, Q D M with speed and Davydov points and LI procedure. At N=10, both F D , Q D M with speed points and LI have not converged, but all four methods have converged by N=30. The convergence behavior is similar for if>7 as shown in Figs. 3.9 and 3.10. Note that F D scheme, Q D M with speed points and LI procedure give inaccurate results for N=15, but have not completely converged by N=30. Summing up results from the tables and figures, we can conclude that the Q D M with Davydov points is more efficient than either the F D scheme, Q D M with speed points or LI scheme. 3.5.2 Calculation of electron velocity distribution function In addition to the determination of the eigenvalues and eigenfunctions, we are also interested in the study of the time relaxation of the electron V D F . In previous papers by Shizgal et al. [43]- [45], the time dependent electron distribution is calculated by the expansion in the eigenfunctions of the F P operator defined in Eqs. (3.2.15) to (3.2.17). Chapter 3. Numerical Methods of Solution of the FPE for Electron Thermalization 63 Figure 3.3: Eigenfunction tps for Fokker-Planck operator with helium cross section at E/n = 0.25 T d and N = 5, 8, 10 and 15, where (â€¢ Q D M with Davydov points, + Q D M with speed points, o F D scheme, * LI points based on Legendre polynomials). Chapter 3. Numerical Methods of Solution of the FPE for Electron Thermalization 64 Figure 3.4: Eigenfunction xps for Fokker-Planck operator with helium cross section at E/n = 0.50 T d and N = 5, 8, 10 and 15, where (â€¢ Q D M with Davydov pointts, + Q D M with speed pts, o F D scheme, * LI points based on Legendre polynomials). Chapter 3. Numerical Methods of Solution of the FPE for Electron Thermalization 65 Figure 3.5: Eigenfunction %pj for Fokker-Planck operator with helium cross section at E/N = 0.25 T d and N = 8, 10, 15 and 30 , where (â€¢ Q D M with Davydov points, + Q D M with speed points, o F D scheme, * LI points based on Legendre polynomials).. Chapter 3. Numerical Methods of Solution of the FPE for Electron Thermalization 66 Figure 3.6: Eigenfunction fa for Fokker-Planck operator with helium cross section at Ein = 0.50 T d and N = 8, 10, 15 and 30 , where (â€¢ Q D M with Davydov points, + Q D M with speed points, o F D scheme, * LI points based on Legendre polynomials). Chapter 3. Numerical Methods of Solution of the FPE for Electron Thermalization 67 Figure 3.7: Eigenfunction ^ 3 f Â° Fokker-Planck operator with argon cross section at E/n = 0.25 T d and N = 8, 10, 15 and 30 , where (â€¢ Q D M with Davydov points, + Q D M with speed points, o F D scheme, * LI points based on Legendre polynomials). r Chapter 3. Numerical Methods of Solution of the FPE for Electron Thermalization 68 Figure 3.8: Eigenfunction fa for Fokker-Planck operator with argon cross section at E/n = 0.50 T d and N = 10, 15, 20 and 30 , where (â€¢ Q D M with Davydov points, + Q D M with speed points, o F D scheme, * LI points based on Legendre polynomials). Chapter 3. Numerical Methods of Solution of the FPE for Electron Thermalization 69 Figure 3.9: Eigenfunction ipj for Fokker-Planck operator with argon cross section at E/n = 0.25 T d and N = 10, 15, 20 and 30 , where (â€¢ Q D M with Davydov points, + Q D M with speed points, o F D scheme, * LI points based on Legendre polynomials). Chapter 3. Numerical Methods of Solution of the FPE for Electron Thermalization 70 Figure 3.10: Eigenfunction fa for Fokker-Planck operator with argon cross section at E/n = 0.50 T d and N = 10, 15, 20 and 30 , where (â€¢ Q D M with Davydov points, + Q D M with speed points, o F D scheme, * LI points based on Legendre polynomials). Chapter 3. Numerical Methods of Solution of the FPE for Electron Thermalization 71 The expansion coefficients a are determined from the initial g(x, 0) and are given by n roo Â°n(0) = / w(x)g(x,0)R (x)dx. (3.5.4) n J0 The Q D M works well for long time but can fail for very short times. In this section, we employ an alternative approach which uses the direct time integration of the discretized equations both for F D , Q D M and LI schemes. A n implicit backward Euler method which provides stable solutions is used to advance the distribution in time. This analysis is useful with regards to the application of pseudospectral methods to diffusion equations and fluid flow problems. In Fig.. 3.11, we show the time evolution of the V D F for electrons in xenon with E/n = 0.25 and 0.5 T d . The results are obtained with the F D scheme with 100 points for an initial Gaussian distribution, that is ' x D(x)g(x,0) 2 = x exp(-4(x 2 - x ) ), 2 0 (3.5.5) where XQ is the initial peak of the distribution. We have used the F D scheme with a large number of points as the "exact" distribution owing to the positivity constraint as discussed by Larsen et al [10.8]. The distribution function remains positive for all times due to the nature of the V matrix with elements defined in Eqs. (3.4.4), (3.4.5) and (3.4.6). When the electric field is applied, the steady state distribution peaks at a reduced speed x that is well past the Ramsauer Townsend minimum. Consequently, only the momentum transfer cross section at the higher energy is sampled and a similar relaxation behavior towards the steady state is seen in both Figs. 3.11A and 3.11B. Results from both Q D M with Davydov points, F D scheme and LI method based on Legendre polynomials with E/n = 0.25 and 0.5 T d are also presented in Figs. 3.12 and 3.13 respectively. A comparison of the accuracy between these schemes is performed. The distribution for t' = 0.0005 is shown in Figs 3.12A and 3.12B and for t' = 0.002 is shown in Figs 3.12C and 3.12D for N = 15 and N = 30 respectively. The initial condition Chapter 3. Numerical Methods of Solution of the FPE for Electron Thermalization 72 x Figure 3.11: Time evolution of the V D F , f(x, t') for electrons in xenon at (A) E/n = 0.25 Td with t' equal to (a) 0.0, (b) 0.0005, (c) 0.001, (d) 0.002 and (e) co; (B) E/n = 0.50 T d with t' equal to (a) 0.0, (b) 0.0002, (c) 0.0005, (d) 0.001 and (e) oo. The initial condition is a Gaussian distribution peaked at x = 9.5 and 11.2 respectively and T d in units of ( 1 0 V c m ) . The distribution functions are obtained with the F D method with 100 points. 0 _17 2 Chapter 3. Numerical Methods of Solution of the FPE for Electron Thermalization 73 0.8 r - ^ - r Figure 3.12: Electron velocity distribution function, f(x,t') for electrons in xenon at E/n = 0.25 T d ( 1 0 - V c m ) and At' = 0.000001. t' = 0.0005 and N is equal to (A) 15, (B) 30; t' = 0.002 and N is equal to (C) 15, (D) 30, where (â€¢ Q D M with Davydov points, o F D scheme, * LI points based on Legendre polynomials). 17 2 Chapter 3. Numerical Methods of Solution of the FPE for Electron Thermalization 74 is a Gaussian distribution with a initial peak at x = 9.5. The solid line is the "exact 0 solution" calculated with the F D method with 100 points. The time step for all methods in this case is At' = 0.000001. In Fig. 3.12A with N =15, only the Q D M solutions gives good results. The F D scheme has not completely converged and the LI procedure exhibits oscillations and departs significantly from the exact distribution. With N = 30, all the methods are essentially exact as can be seen in Fig. 3.12B. A t t' â€” 0.002, all the transients have decayed and both Q D M and F D methods yield accurate results with as few as 15 points. The LI solutions are now positive but totally inaccurate as shown in Fig. 3.12C. As we go from N = 15 to N = 30, all methods give accurate results. In Fig. 3.13, the results of a similar study is shown for E/n = 0.50 T d . As shown in Fig. 3.11, the steady state distribution is shifted outwards with respect to the case in Fig. 3.12; therefore, we have chosen the initial Gaussian with a initial peak at larger x â€” 11.2. In this case, the distribution for t' = 0.0002 is shown in Figs. 3.13A and 0 3.13B and for t' = 0.001 is shown in Figs. 3.13C and 3.13D for N = 20 and N = 50 respectively. In order to acquire accurate results for such time scales, we pick a smaller time step A t ' = 0.0000005. At t' = 0.0002, the Q D M gives good results for N = 20, while the solutions of the F D and LI schemes have not converged. Even at the longer time, t' = 0.001, only the Q D M provides converged solutions for N = 20. With N = 50, all methods give fully converged results for both times. A n explicit forward Euler time integrating scheme has also been studied and compared to the implicit one. Both explicit and implicit schemes agree very well with the above given time steps. However, if we choose a bigger time step than the given one, only the implicit Euler method will give stable solutions. With an explicit scheme, the initial errors would not be damped, but rather amplified as time proceeds and the time integration becomes unstable. The choice of the time step is dictated by the accuracy with which one wants to represent the initial transients of the distribution. For a given N , each method of discretizing Chapter 3. Numerical Methods of Solution of the FPE for Electron Thermalization 75 Figure 3.13: Electron velocity distribution function, f(x,t') for electrons in xenon at E/n = 0.5 T d ( 1 0 - V c m ) and At' = 0.0000005. t' = 0.0002 and N is equal to (A) 20, (B) 50; t' = 0.001 and N is equal to (C) 20, (D) 50, where (â€¢ Q D M with Davydov points, o F D scheme, * LI points based on Legendre polynomials). 17 2 Chapter 3. Numerical Methods of Solution of the FPE for Electron Thermalization 76 the F P operator gives a set of N eigenvalues, X , n = 1 n N of which AN is the largest. The spectrum of eigenvalues determines both the relaxation behavior as given by Eq. (3.2.17), and the stability and time step restrictions for the time integration of the discretized equations. Clearly, each term in the sum in Eq. (3.2.f 7) decays on a time scale of the order A " and the fast transients will be well represented by the numerical solution 1 if AN A t ~ 1. For a first order explicit method used here, the smallest time step size is 0 ( 1 / A N ) . In Fig. 3.14A, we show AN as a function of N for the hard sphere cross section at E/n = 0 for both F D and Q D M with Davydov points. From the results, we see that the discretization with Davydov points is less restrictive with regards to time step in the explicit Euler scheme. However, this condition is not true in the case of E/n ^ 0. Fig. 3.14B shows that the largest eigenvalues AN are very sensitive to the applied electric field at large N . As a result, for large N , the Q D M with Davydov points requires a smaller time step to achieve stability than the F D scheme. Fortunately, for the time integration of the Q D M discretized equation, only a small number of Q D M Davydov points is needed in order to obtain converged results. With results from tables and figures, we conclude that the Q D M with Davydov points provides faster convergence and more accurate solutions than other discretization methods in all aspects that we have studied. Chapter 3. Numerical Methods of Solution of the FPE for Electron Thermalization 77 4000 0 10 20 30 40 50 60 N Figure 3.14: Variation of the maximum eigenvalue AN for the Fokker-Planck operator with hard sphere cross section versus N ( A ) AN verse N at zero field for both curves (a) F D scheme and (b) Q D M with Davydov points. (B) The variation of AN with N for Q D M with Davydov points at various field strengths, E/n (a) 0, (b) 0.00025, (c) 0.0005, (d) 0.001 T d 2 Chapter 4 Electron Thermalization and Attachment in C C l / A r and C C l / N e 4 4.1 4 Introduction The sttidy of the thermalization of electrons in atomic and molecular moderators is an important field of study and has a long history. Detailed reviews of the present status of the field that provide a bibliography of previous works were presented by Braglia [42] and by Shizgal et al [93]. The nature of the time dependent electron distribution function and average properties are important aspects of gas discharges [115], radiation chemistry and physics [116], the development of gaseous insulators and switching devices [117], the interpretation of electron swarm experiments [118], environmental concerns and kinetic processes [119,120]. The main objective is to consider an initial nonequilibrium electron distribution function and to determine the time dependent distribution function and the nature of the steady distribution at long times. The thermalization of energetic electrons in atomic and/or molecular moderators proceeds owing to collisions between the electrons and the constituents of the moderator assumed to be present in large excess. Electron-electron collisions are assumed not to occur and are not included in the theoretical analysis. A complete treatment would require the inclusion of elastic and inelastic collisions, and chemically reactive processes such as ionization and attachment. Recent papers reported continuing studies of electron thermalization in atomic moderators [75]- [77] in comparison with the earlier work by Shizgal and co-workers [43]- [45], [96], [101]- [103], [121], Mozumder [40,98] and Koura [41,100]. 78 Chapter 4. Electron Thermalization and Attachment in CCU/Ar and CCU/Ne 79 There have also been several papers on the thermalization of electrons in molecular moderators such as C H [77,86], S i H [122], H [85] and C H / A r mixtures [77]. In addition 4 4 2 4 to these calculations, it is of considerable interest and importance to consider electron thermalization in electronegative gases for which it is important to include the attachment of electrons. The attachment of electrons to electronegative gases in swarm type experiments is an important process in connection with gas discharges, plasma etching, radiation chemistry and physics. The main objective of this chapter is the study of the time evolution of the electron nonequilibrium distribution function from some initial distribution under the influence of elastic, and inelastic collisions of electrons with the ambient gases, and electron attachment to some electronegative gas. We consider a mixture of either A r or Ne and CC1 so that the electron collisions to be considered are elastic electron-inert gas and 4 electron-CC1 , and vibrationally inelastic electron-CCl collisions and electron attach4 4 ment to CC1 . The rate constant for the attachment of electrons to CC1 is very large. 4 4 The inert gas moderators were chosen so as to have an example of a gas (Ar) with a momentum transfer cross section characterized by a Ramsauer-Townsend minimum, and a second moderator (Ne) for which the momentum transfer cross section is nearly constant with energy. It is of considerable interest to contrast the relaxation behavior for these different moderators. We consider a numerical solution of the time dependent B E for the electron energy distribution function, and calculate the time variation of the electron number density and electron average energy. These quantities are characterized by time dependent energy relaxation rate coefficients and the attachment rate coefficient. Relaxation times are calculated in terms of the time required for the electron energy to reach within 10% of the steady average energy. The steady average energy of the electrons can be above the ambient temperature of the moderators because of the effect of attachment heating. This Chapter 4. Electron Thermalization and Attachment in CCl /Ar 4 and CCl /Ne 4 80 effect was discussed in the previous papers by Shizgal [123]- [125] except that in these works inelastic collisions were neglected. The previous works by Shizgal showed that there can be a steady state distribution even though electron attachment involves a loss of electrons and a steady state would not be expected to occur. The basis for this theoretical work was the F P E for elastic electron inert gas collisions and a term due to electron attachment. However, if there is a separation of relaxation times such that the energy thermalization occurs on a much faster time scale than attachment, then there can be a pseudo-steady distribution. The energy per particle attains a limiting value with the electron energy and the electron number density varying slowing with time on the time scale of the attachment process. In a recent paper, Kowari and Shizgal [126] extended this previous work to include vibrationally inelastic collisions as well. The main thrust of this previous paper was with regard to the existence of a steady distribution in the presence of the attachment process and a loss of electrons. Unfortunately, the low energy portion of the attachment cross section used in this previous paper, the portion below O.lOeV, was not the correct cross section for the system studied. The calculations are repeated here for the correct cross section and the results are reinterpreted. Warman and Sauer [54] reported an experimental technique to determine the energy relaxation rate coefficients and the relaxation times in pure inert gases and other moderators by following the reaction kinetics of electron attachment to C C I 4 . In this chapter, we consider the experimental method used by Warman and Sauer to extract relaxation times from the electron kinetics associated with attachment. We compare the relaxation times extracted from the electron kinetics with those determined directly from the relaxation of the average energy. We also compare with the very recent measurements by Shimamori and Sunagawa [55]. The formalism of this chapter follows the earlier papers on electron relaxation with the Chapter 4. Electron Thermalization and Attachment in CCl /Ar 4 and CCl /Ne 4 81 inclusion of vibrationally inelastic collisions but without electron attachment [85,86,122], and the recent paper [126] that also includes electron attachment. In Section 4.2, the B E for the electron energy distribution function is considered, and the numerical solution is discussed. The results and their discussion is presented in Section 4.3. A reinterpretation of the work of Kowari and Shizgal [126] with regard to the existence of a steady electron distribution is presented in Section 4.4. 4.2 The Boltzmann Equation The theoretical method employed here follows closely that used previously [85,86] which is based upon the B E with collision terms for elastic, vibrationally inelastic collisions. The changes from this previous work due to the inclusion of the attachment to the molecular gas was discussed by Kowari and Shizgal [126]. Because there is no electric field applied to the system, it is reasonable to consider only the isotropic portion of the electron distribution function. We also assume that the electron distribution function is spatially homogeneous. 4.2.1 Inelastic collision operator We write the B E with the collision term which explicitly includes electron attachment as well as elastic and inelastic collisions, as given by a/ J(f) dt Mf) where / = f(v,t). In Eq. (4.2.1), J i(f) e sion operators, respectively, and J {f) a + Mf) + and Ji (f) n Uf), (4.2.1) are the elastic and inelastic colli- is the electron-attachment collision operator. Eq. (4.2.1) no longer conserves the number of electrons because of electron attachment. We Chapter 4. Electron Thermalization and Attachment in CCl /Ar 4 and CCU/Ne 82 denote the mass of an electron by m and that of a molecule by M. For the small mass ratio m/M the elastic collision operator is very well approximated by the differential F P operator [43]- [45], [96], [101]- [103], [121] m AL d ~^ MyV z r , ... ov kTh d ,, â€ž. , mv ov where AL, is the number density of moderator 7 , k is the Boltzmann constant, T is the temperature of the moderator, v is the speed of an electron, and cr^(v) is the momentum transfer cross section for collisions of electrons with moderator 7 . We take 7 = 1 for CCI4 and 7 = 2 for either A r or Ne. The inelastic collision operator can be written as the sum over all the inelastic processes involved and can be cast in the form of a difference operator. The derivation of the inelastic collision operator is given in the Appendix of Ref. [86]. If the total inelastic cross section between molecular states i and j is c;j(t>), then the inelastic collision operator is given by Jin^^v^nâ€”aijivy^-niaijiv)/^ hj where v' = \Jv â€” (4.2.3) and e;j = tj â€” t{ is the excitation energy transfer for internal 2 molecular states i and j . The collision is inelastic if > 0 and j > i, and it is super elastic if e j < 0 and j < i: The population density of the molecular state is denoted by 2 n i â€” J2.e~ i/ b e kT ' a ^ i U i ~ ' The electron attachment collision operator is J = â€”Niva (v)f(v), a section for electron attachment to CCI4 is a (v). a a The B E for the study of the relaxation of electrons is then, df ^ + miV-y d 4 { -.^ , kT d b NivY,[niâ€”o-ij(v')f(v')-niaij(v)f(v)] v 12 where the cross Chapter 4. Electron Thermalization and Attachment in CCU/Ar - and CCU/Ne Niv<r {v)f(v). 83 (4.2.4) a In radiation physics and chemistry, it is customary to use an energy-dependent density p(c,t) = 27r(^)2 y/ef(v, t), where e = and the track-length distribution (or incremen- tal degradation spectrum) z(e,t) = vp(e,t). With the definitions of p(e,t) and z(e,t) in Eq. (4.2.4), we can obtain the expression for the time evolution of z(e,t), that is, ldz(e,t) f 2 V 3 / d 2 + A i ^ [ n i e r ^ e ' ^ e ' , * ) - ni<Tij(e)z(e, t)] - N va (e)z(e,t), 1 (4.2.5) a where e' = |mt>' = e + e,;j. 2 We consider the relaxation of electrons in a gas due to elastic and inelastic collisions. The inelastic collisions that may be important in the low energy regime are rotational and vibrational collisions; however, the present calculations do not involve rotationally inelastic collision processes explicitly because we use the momentum transfer cross section that is vibrationally elastic. The cross section for vibrational excitation from a vibrational state V to V in the vibrational mode v is denoted by a^ ' v E^ v and the energy transfer by . The population density of the vibrational state V in the vibrational mode v is denoted by . It is necessary to replace the summation J2ij = J2i 12 j with an explicit reference to the vibrational states, so that 12i 12 j â€” 12v 12v 12vThe B E with the differential operator for elastic collisions and the difference operator Chapter 4. Electron Thermalization and Attachment in CCU/Ar and CCl /Ne 4 84 for vibrationally inelastic collisions is expressed by 1 dz(e,t) dt = J( ) + S(e~e )S{t) z 0 d =E . d ( \ 5 -tcr^(e)z + â€” (e o- (e)J fcT â€” 2 7 + e ay(e)kT}) (z 6 2 + de 2 ve, Ai v V V - EE-rE-r(^,o . v + V N T (e)z(e,t) l( a V (4.2.6) S(e-e )6(t). 0 In Eq. (4.2.6), the delta functions on the right-hand side show explicitly the initial electron distribution as a delta function at energy erj. The electron number density and energy are given by n(t) = j f(v,t)dv (4.2.7) O TKIV / f(v,t)â€”dv, respectively. The average electron energy per electron is then E(t) = 4.2.2 Etotal (4.2.8) (4.2.9) n (*) Cross Section Data Set The momentum transfer cross sections for elastic collisions with either A r , Ne or C C I 4 are shown in Fig. 4.1. The momentum transfer cross section for A r used in the present calculations is constructed as follows. We obtain the data below 1 eV by using the modified effective range theory (MERT) procedure given by Haddad and O'Malley [127]. We calculate the M E R T procedure for the energy range of 0 eV and 0.32 eV with the M E R T parameters [127] for the upper limit of 0.5 eV and for the energy range of 0.32 eV Figure 4.1: Energy dependence of momentum transfer cross sections. cr (E) is in units of 1 0 c m for (a) CC1 , (b) Argon and (c) Neon. m - 1 6 2 4 Chapter 4. Electron Thermalization and Attachment in CCl /Ar and CCl /Ne 4 4 86 and 1 eV with that for the upper limit of 1 eV. We use Table 1 of Ref. [128] between 1 eV and 4 eV and Table 2 of Ref. [129] above 4 eV. The momentum transfer cross section of Ne is taken from O'Malley and Crompton [130], and the one for CC1 is from Hayashi [131]. 4 The cross section for A r and C C I 4 are characterized by Ramsauer-Townsend minima whereas the one for Ne is not. CCI4 belongs to T d symmetry and has nine normal modes of which only four, denoted by vl â€” z/4, are distinguishable because of degeneracy. The excitation energies from the vibrational ground state to the first excited vibrational state for each mode are 0.0569 eV, 0.0269 eV, 0.0962 eV, and 0.0389 eV, respectively [130]. The four vibrationally inelastic cross sections characterized by these threshold energies are shown in Fig. 4.2. Because the population densities of some of the highly excited vibrational states in each normal mode are not negligible (for example, nj, is 9.8 percent, n* is 22.8 percent, n* a 2 3 is 2.4 percent and n* is 17.3 percent), we take account of the transitions up to V = 3 , 4 8, 2, and 5 for vl â€” v4, respectively; all with A V = 1 . In order to take account of these transitions, we need cross sections from vibrationally excited states to the upper states of those states in each vibrational mode. However, available cross sections on vibrationally inelastic collisions for both modes are those for V=0â€”>1, and the scaling laws Ref. [132] to construct vibrationally inelastic cross sections for vibrationally excited states are not applicable to the present case. Therefore, we assume the cross sections with the transition from V to V + l where V>1 for the four normal modes are the same as those from V=0 to 1. The cross sections with the transitions from V + l to V for the four normal modes are obtained using the relation for microscopic reversibility, e<7;j(e) = (e â€” eij)cfji(e â€” e j). S We use Hayashi's electron attachment cross-section [131] for electron attachment to CCI4. Hayashi's electron attachment cross section is consistent with that of Chutjian and Alajajian [133] at low energies, and we use the analytical expression of Chutjian and Chapter 4. Electron Thermalization and Attachment in CCU/Ar 0.01 0.1 1 and CCU/Ne 87 10 E (eV) Figure 4.2: Energy dependence of the four vibrationally inelastic cross sections for electron-CCU collisions. The threshold energies are, 0.0269, 0.0389, 0.0569 and 0.0962 eV, respectively. Chapter 4. Electron Thermalization and Attachment in CCl /Ar and-CCl /Ne 4 88 4 Alajajian cr(e) = 5 4 5 { ^ ^ e x p [ - ( â€” ^ - ) 1 + e x p ( â€”)} . 0.005 0.056 (4.2.10) 2 ; V ; J F V J; v ; below 0.1 eV. The thermal attachment rate coefficient at 300K is given by K h = t 2.94xl0 cm s _7 3 _1 in excellent accord with the value of 2.96xl0 cm s~ _7 3 1 at 296K re- ported by Blaustien and Christophorou [135], and in good agreement with the value of 3.5xl0 cm s _7 3 _1 at (298Â± 3K). The attachment cross section is shown in Fig. 4.3. The dashed curve in the figure is the power law fit to the high energy portion of the attachment cross section. In the previous papers [124]- [126], this analytic expression was incorrectly reported with power of e as | rather than as In the first two papers cited, this was just a typographical error and the calculations were reported for the correct cross section. However, in the most recent work by Kowari and Shizgal [126], an error in a computer code resulted in the use of an incorrect cross section. This is discussed in detail in Section 4.4. 4.3 Calculations and Results The B E for the specific cross sections is solved with a finite difference technique for both energy and time as discussed in previous papers [85,122]. The initial distribution function is a delta function at 1.0 eV. We use a temperature of 300 K for the background gases CCI4 and either A r or Ne with the total number density of 2 . 6 9 x l 0 c m . At each time 19 -3 step we solve the finite difference equation with the successive over-relaxation method [136]. In theoretical calculations, unlike in experiments, it is possible to study the effect of the omission or inclusion of different cross sections to understand the role of different processes. In Fig. 4.4, we show the relaxation of the electron energy in C C L t / A r (Fig. 4.4A) and C C L i / N e (Fig. 4.4C) in the absence of the electron attachment process. The Chapter 4. Electron Thermalization and Attachment in CCl /Ar 4 and CCl /Ne 89 4 Figure 4.3: Energy dependence of the attachment cross sections for electron-CCL, collisions. The solid line is the energy dependence of the electron-CCL, attachment cross section. The dashed line is the approximate power law cross section given by a (E) = 0.5178.E- - , with E in eV and a in 1(T cm . a 1 88 16 a 2 Chapter 4. Electron Thermalization and Attachment in CCU/Ar and CCU/Ne 90 Figure 4.4: Electron relaxation in C C U / A r and C C U / N e mixtures without attachment. (A) Time variation of the average energy in C C U / A r for different C C U mole fractions equal to (a) 3xl0~ (b) 10~ , (c) 5xl0~ , (d) 3 x l 0 - and (e) IO" . (B) Time variation of the energy relaxation coefficient in C C U / A r for different C C U mole fractions equal to (a) 1 0 , (b) 5 x l 0 and (c) 3 x l 0 . (C) Time variation of the average energy in C C U / N e for different CC1 mole fractions equal to (a) 10" , (b) 6 x l 0 - , (c) 3 x l 0 " and (d) 10~ . ( D ) Time variation of the energy relaxation coefficient in C C U / N e for different C C U mole fractions equal to (a) IO" , (b) 6 x l Q - , (c) 3 x l Q - and (d) IO" . 4 -4 4 5 - 5 5 5 - 5 4 5 5 4 5 4 5 5 5 Chapter 4. Electron Thermalization and Attachment in CCU/Ar and CCU/Ne 91 different curves are for different CC1 mole fractions. The electron energy decreases with 4 time and attains the average energy E h = Â§ ArTJ,, which for a moderator temperature at t 300K is 0.0388 eV, so that E(0)/E = 25.79. The graphs on the right-hand-side in Fig. th 4.4 show the time variation of the energy relaxation rate coefficient defined by, K (t) E E(t) = - l d E(t) Eth th The behavior of the electron energy relaxation with A r or Ne as moderators is similar. The rapid increase in the rate of energy relaxation with the addition of C C I 4 is due to the inclusion of the vibrationally inelastic collisions. If Kg in Eq. (4.3.1) is assumed to be constant, the integration of Eq. (4.3.1) leads to an exponential decay of the average energy, that is, . â€¢ . (m_^) where E(oo)/E h t = (m_*MyÂ« ,. (4 2) = 1 if there is no attachment. The time variation of the energy rate coefficients is indicative of the fact that the energy relaxation is not a pure exponential decay. This is anticipated since for this linear problem, the time dependence can be expressed as a sum of exponential terms with each term characterized by an eigenvalue of the collision operator. This was the approach adopted with analysis that did not include the'vibrationally inelastic collisions [123]- [125]. With the inclusion of inelastic collisions, the eigenvalue spectrum of the operator should be modified by the occurrence of very large eigenvalues that characterize the time scale for vibrationally inelastic collisions. Semilogarithmic plots of E(t)/E h t versus time that are not shown here suggest that the energy relaxation could be fitted to a sum of three exponential terms. This observation is similar to that reported recently by Krajcar-Bronic and Kimura [77] who studied electron energy relaxation in pure C H , and for which they demonstrated a relaxation 4 occurring over two distinct time domains. It would be of considerable interest to extract and compare from both the theoretical and experimental data, the 2 or 3 dominant Chapter 4. Electron Thermalization and Attachment in CCl /Ar 4 and CCU/Ne 92 exponential terms [137]. If Eq. (4.3.2) were valid then there is a simple relationship between K# and pr^j given by [54], 2.303 r i , . E(0) . , where N = 3.22xl0 c u T is the density at 1 torr and 300K. 16 3 In Figs. 4.5 and 4.6, we show the time variation of the electron energy and density with the inclusion of the attachment process for C C l / A r and C C l / N e mixtures, respectively. 4 4 The main difference compared with the behavior without attachment is that the electron number density decays rapidly and the average electron energy attains a steady value above the thermal energy of the moderators. This is the attachment heating effect. The effect increases with increasing concentration of C C I 4 as shown in Fig. 4.7 and in Tables 4.1 and 4.2. In Figs. 4.5B and 4.6B, we also show the energy relaxation rate coefficients, and in Figs. 4.5D and 4.6D the time dependence of the attachment rate coefficient defined by, T r l K(i) . 1 dn = - - T t . . . (4.3.4) The energy relaxation rate coefficients shown in Figs. 4.5B and 4.6B vary with time indicative of a multi-exponential decay, as was the case without attachment. The inclusion of attachment appears to have reduced the K# values by roughly a factor of 2 in the case of A r as moderator. The reduction in K g with Ne as the moderator is somewhat less. The attachment rate coefficients increase rapidly and attain a steady value which is significantly less than the thermal values at 300K. The steady electron distribution function is presumably considerably perturbed from a Maxwellian at 300K. This nonequilibrium effect appears to be somewhat larger for A r than for Ne, perhaps due to the overall smaller momentum transfer cross section for A r . The main quantity of interest, in particular when comparing with experiment, is the relaxation time for the energy relaxation. Since the relaxation of the average energy Chapter 4. Electron Thermalization and Attachment in CCU/Ar and CCU/Ne 93 Figure 4.5: Electron relaxation in C C L i / A r with attachment. (A) Time variation of the average energy for different CCL4 mole fractions equal to (a) 1 . 5 x l 0 , (b) 1 0 , (c) 5 x l 0 ~ , (d) 3 x l 0 and (e) 1 0 . (B) Time variation of the energy relaxation coefficient for different C C I 4 mole fractions equal to (a) 1 0 , (b) 5xl0~ and (c) 3 x l 0 . (C) Time variation of the number density for different C C I 4 mole fractions equal to (a) 1 . 5 x l 0 , (b) 1 0 , (c) 5 x l 0 , (d) 3 x l 0 and (e) 1 0 . (D) Time variation of the attachment rate coefficient for different CC1 mole fractions equal to (a) 10~ , (b) 5xl0~ , (c) 3 x l 0 ~ , (d) lCT and (e) 5 x l Q - . -3 4 - 4 - 3 -4 -4 5 - 5 -3 -3 - 4 - 4 -4 3 4 4 5 4 4 Chapter 4. Electron Thermalization and Attachment in CCU/Ar and CCU/Ne 94 Figure 4.6: Electron relaxation in C C U / N e with attachment. (A) Time variation of the average energy for different C C I 4 mole fractions equal to (a) 3 x l 0 , (b) 1 0 , (c) 3 x l 0 and (d) 1 0 . (B) Time variation of the energy relaxation coefficient for different C C I 4 mole fractions equal to (a) 1 0 , (b) 3 x l 0 and (c) 1 0 . (C) Time variation of the number density for different C C I 4 mole fractions equal to (a) 3 x l 0 , (b) 1 0 , ( c ) 3 x l 0 and (d) 1 0 . (D) Time variation of the attachment rate coefficient for different C C I 4 mole fractions equal to (a) 5xlCT , (b) I O , (c) 5 x l ( T and (d) IO" . - 4 - 4 - 5 -5 -4 - 5 -5 - 4 -4 -5 4 - 4 5 5 -5 ChapterA. Electron Thermalization and Attachment in CCl /Ar 4 and CCl /Ne 95 4 (and the number density) is not a pure exponential, there is no unique definition of a relaxation time. It has become customary to define a relaxation time as the time required to approach within 10% of the thermal energy, denoted by T ^ . The variation of this thermalization time, with and without attachment, versus CC1 mole fraction is 4 shown in Tables 4.1 and 4.2. Also shown in the tables is the steady electron energy relative to the thermal energy. This illustrates the attachment heating effect with increasing CC1 concentration. The results in Table 4.1 show that the thermalization time 4 decreases rapidly with an increase in the CC1 concentration. This is due primarily to 4 the increased energy loss from vibrationally inelastic collisions. The electron attachment process increases the relaxation time at all concentrations shown in Table 4.1. This is consistent with the decrease in the KE values shown previously in Figs. 4.4B and 4.5B. With increasing C C I 4 concentration, the nonequilibrium attachment heating effect increases strongly as seen from the rapid increase in E(oo)/E ht The asterisks in the table indicate that a relaxation time could not be defined owing to a minimum in the time variation of E(t)/E h. t The results for Ne as moderator in Table 4.2 are similar in some respects and different in others. The relaxation times decrease with increasing C C I 4 concentration but the inclusion of the attachment process decreases the relaxation times for small concentrations and increases the relaxation time for the larger concentrations. It is important to note that the results for very small concentrations that are shown in Table 4.2 were not included in Table 4.1. The determination of the final steady electron energy for the smallest concentrations requires considerable computation and their determination is difficult. We did not obtain reliable results for A r at the smallest concentrations (Xcci ~ 1 0 ) , -8 4 and for Ne only the relaxation times without attachment could be obtained. In this case, the final energy is the thermal energy, E(oo)/E h t = 1, and the time integrations do not have to be carried out to a steady state. The asterisks in Table 4.2 also indicate that Chapter 4. Electron Thermalization and Attachment in CCU/Ar Xcci (a) 4 1.0( -7) 3.0( -7) 5.0( -7) 1.0( -6) 3.0( -6) 5.0( -6) 8.0( -6) 1.0( -5) 2.0( -5) 3.0( -5) 4.0( -5) 5.0( -5) 6.0( -5) 7.0( -5) 8.0( -5) 1.0( -4) 3.0( -4) 5.0( -4) 1.0( -3) (*>) PTl.l PT-l.l 5091 2942 2171 1404 683.5 481.4 343.6 291.1 169.1 120.6 94.37 77.48 66.01 57.44 50.84 41.30 14.52 8.90 4.59 5416 4219 3141 1769 812.7 572.6 416.7 354.5 220.5 164.6 133.3 112.2 96.29 84.03 73.49 59.33 * * * and CCl /Ne 4 96 Eoo/Eth 1.445 1.741 1.982 2.344 2.735 2.903 3.055 3.127 3.357 3.495 3.597 3.683 3.760 3.827 3.896 3.982 4.830 5.219 5.829 C C I 4 in A r without attachment. in A r with attachment. PT1.1 = 5570 ps torr for Ar; ref [43]. ^ CCI4 Table 4.1: Electron energy relaxation times pri.i (in units of //sec-torr) and ratio between Eoo/Eth for C G I 4 in A r at different concentration. Chapter 4. Electron Thermalization and Attachment in CCl /Ar 4 (a) XcCL, P^l.l 1.0( -8) 3.0( -8) 5.0( -8) 6.0( -8) 8.0( -8) 1.0( -7) 2.0( -7) 3.0( -7) 4.0( -7) 6.0( -7) 8.0( -7) 1.0( -6) 3.0( -6) 1.0( -5) 3.0( -5) 1.0( -4) 3.0( -4 2944 2842 2748 2704 2620 2543 2224 1988 1804 1537 1349 1209 638.1 271.0 112.8 39.67 14.30 0) PTl.l * * * * 2339 2083 1802 1583 1311 1035 873.6 769.2 448.2 248.2 126.9 52.44 13.12 and CCU/Ne 97 Eooj Eth * * * 1.087 1.104 1.134 1.200 1.270 1.323 1.411 1.476 1.526 1.812 2.187 2.762 3.475 4.502 in Ne without attachment. C C I 4 in Ne with attachment. pr .i = 2960 /is torr for Ne; ref [43]. ( ) CCI4 a ( i ) a Table 4.2: Electron energy relaxation times pri.i (in units of /isec-torr) and ratio between Eoo/Eth for C C I 4 in Ne at different concentration. Chapter 4. Electron Thermalization and Attachment in CCU/Ar and CCU/Ne 98 reliable results could not be obtained for the entries at the lowest concentrations. In Figure 4.7, the solid curves are for C C l / N e and the dashed curves are for C C l / A r . The 4 4 curves labelled (a) correspond to the heating in the absence of inelastic collisions whereas the curves labelled (b) are for the situation with inelastic collisions. The effect of the inelastic collisions to cool the electrons is clear. In Fig. 4.8, we show the results of calculations without the attachment process for C C l / N e mixtures that are also included in 4 Table 4.2. The dots show the results of the calculations at five CC1 concentrations. The 4 solid line through the points is a linear fit the slope of which gives TI.I for CC1 without 4 attachment. We find from the slope that pTi.i = 1.68xl0 -3 ps torr without attachment. Shimamori and Sunagawa reported the experimental value of 3xl0~ ps torr, which in3 cludes the effect of attachment. The lower solid curve is the result without inelastics and with attachment. The intercept in the figure provides the relaxation time for pure Ne and is 2994 ps.torr in agreement with the previous result of 2960 ^s.torr [43]- [45]. The concentration dependence of the relaxation times for both C C l / A r and C C l / N e 4 4 for much larger concentrations is shown in Fig. 4.9. The curve labelled (a) is with inelastics and without attachment. The curves labelled (b) include both inelastics and attachment. The lower solid curve labelled (c) does not include inelastic collisions but does include attachment. This illustrates the lengthening of the relaxation time due to the attachment process at relatively large CC1 concentrations. Notice that for A r as 4 moderator the curve labelled (b) lies above the curve labelled (a) at the lowest concentrations. This shows the decrease in the relaxation time due to attachment as previously discussed in connection with the results in Table 4.2. In Fig. 4.10, we carry out a calculation of the energy relaxation rate coefficient, assumed independent of time, from the time variation of the electron density as previously done by Warman and Sauer [54] in the interpretation of their experimental results. Their work is based on the power law approximation of the attachment cross section versus Chapter 4. Electron Thermalization and Attachment in CCU/Ar x and CCU/Ne 99 cci ( " ) 10 4 4 Figure 4 . 7 : The attachment heating effect for electrons in C C U / A r (dashed curves) and C C U / N e (solid curves) mixtures, (a) without inelastic collisions, (b) with inelastic collisions. Chapter 4. Electron Thermalization and Attachment in CCl /Ar 4 and CCl /Ne 4 100 0.40 Q 32 0.0 1 1 1 0.2 1 1 1 1 0.4 x 0.6 1 !â€” 1 1 0.8 1.0 cci ( " ) 10 7 4 Figure 4.8: The variation of'l/pri.i versus mole fraction C C I 4 in C C l / N e mixtures. The solid symbols are the results of the calculation and the solid line through the points is a linear fit. The line below without symbols is the result without inelastic collisions. Note that the mole fraction of C C I 4 is extremely small. 4 Chapter 4. Electron Thermalization and Attachment in CCU/Ar 0.0 0.2 0.4 0.6 0.8 1.0 0.8 1.0 and CCU/Ne 101 Xcci/ " ) 1 0 0.0 0.2 0.4 4 0.6 Xcci/ - ) 1 0 4 Figure 4.9: The variation of 1/pri.i versus mole fraction C C I 4 . (a) with inelastic collisions and without attachment, (b) with inelastic collisions and attachment and (c) without inelastic collisions and with attachment. (A) C C U / A r mixtures, (B) C C U / N e mixtures. Chapter 4. Electron Thermalization and Attachment in CCl /Ar 4 and CCl /Ne 4 102 energy, that is, Â°[ ) E = Y> -- (4 3 5) for which the equilibrium attachment rate coefficient varies with temperature as given by, K{T(t)) = \^a [kT(t)]-v ? r(2 +1 - s), 2 0 (4.3.6) where T(x) is the Gamma function. If one sets, E(t) = |fcT(t) and E h = |A;Ti(0), then t this gives the relation between the time dependent attachment rate coefficient and the time dependent average electron energy, that is, K(t) E{t) K YE,ih -p (4.3.7) J th where p = s â€” | . Eq. (4.3.5) assumes that the electron velocity distribution function remains a Maxwellian with a time dependent temperature. This concept was used previously by Mozumder [40] in his study of electron thermalization and a similar approach was also considered in the treatment of hot atom reactions [91], [138]- [140]. If Eq. (4.3.6) is combined with Eq. (4.3.2), we get that In ( K ( t ) y v K th J 1 / v â€¢ -K t. E (4.3. This is the analysis by Warman and Sauer by which they extracted values of Kg from the kinetics of the attachment process. We have used the calculated n(t) and K(t) data and extracted KE values by fitting the results to Eq. (4.3.8). This is shown in Figs. 4.10A and 4.10C for C C l / A r and C C l / N e mixtures, respectively. We have chosen to fit 4 4 the curves between (approximately) the times for which n(t)/n(0) = 0.8 and n(i)/n(0) = 0.05. The dashed lines in Figs. 4.10A and 4.10C are the linear fits to the calculated K(t) curves and the slopes are related to the KE coefficients as given by Eq. (4.3.8). It is important to emphasize that Eq. (4.3.7) is valid for a power law attachment cross section, Chapter 4. Electron Thermalization and Attachment in CCl /Ar and CCl /Ne 4 0 20 40 60 t(ns) 80 100 0.0 0.2 103 4 0.4 0.6 0.8 X (10- ) 1.0 4 ccl4 Figure 4.10: The variation of the attachment rate coefficient versus time and the extraction of energy relaxation rate coefficients. (A) Variation of \o[(K(t)lKth) ^ â€” 1] versus time for C C l / A r mixtures. The dashed lines are linear fits to the calculations. The slopes yield approximate energy relaxation rate coefficients analogous to those determined experimentally. The mole fraction of CC1 are (a) 1 0 , (b) 5 x l 0 and (c) 3 x l 0 . (B) The concentration dependence of the energy relaxation rate coefficients from (A). The triangles are for p = 0.726 and the solid dots are for p = 1.38. (C) Variation of \n[(K(t)/Kth) ^ â€” 1] versus time for C C l / N e mixtures. The dashed lines are linear fits to the calculations. The slopes yield approximate energy relaxation rate coefficients analogous to those determined experimentally. The mole fraction of CC1 are (a) 1 0 , (b) 4xl0~ and (c) 2 x l 0 . (D) The concentration dependence of the energy relaxation rate coefficients from (C). The triangles are for p = 0.726 and the solid dots are for p = 1.38. -1 4 -4 - 5 4 - 5 -1 4 4 - 4 5 - 5 Chapter 4. Electron Thermalization and Attachment in CCl /Ar and CCl /Ne 4 104 4 an electron distribution that remains Maxwellian and time independent K E coefficients. From the results presented here, these assumptions are not rigorously valid. The variation of the KE coefficients with CC1 concentration determined from the 4 slopes of the fits in Figs. 4.10A and 4.10C is shown in Figs. 4.10B and 4.10D for C C l / A r 4 and C C l / N e mixtures, respectively. The symbols are the actual results and the lines 4 are the linear fits to the data points. Our analysis is not rigorous and the "data" points in these figures show some scatter analogous to experimental data. The variation of the energy exchange rate coefficient with CC1 concentration is clear. The results indicated 4 by triangles are for the choice p = 0.726, which is the value used by Warman and Sauer [54]. The results for the circles are the present choice, p = 1.38, which corresponds to the power law fit of the attachment cross section used here. The intercepts of the straight lines in the Figs. 4.10B and 4.10D coincide reasonably well with the K values of pure E A r and pure Ne, respectively. The energy relaxation rate coefficient for electrons in a pure inert gas can be obtained from the F P operator and is given by, 167rAL?7l R E ' = 0M2 where x 0.80xl0 â€” -12 E/kTh cm s 3 _1 values of 0 . 1 3 x l 0 / roo j2KT /m / Jo b (4.3.9) exp{-x )x a {x)dx, 2 5 m is the reduced speed. We find that K is 0 . 8 7 x l 0 ~ c m s and 12 E 3 _1 for A r and Ne, respectively. Warman and Sauer [54] reported the _12 cm s give the value of ( 1 0 3 -13 _1 and 0 . 2 5 x l 0 ~ c m s 12 - 10~ ) c m s 14 3 3 -1 whereas Shirnamori and Sunagawa [55] for A r . The experimental results are somewhat _1 below the theoretical estimates. The intercepts in Fig. 4.10B are somewhat larger for A r whereas the intercepts in Fig. 4.10D for Ne are in good agreement with Eq. (4.3.9). The slopes of the straight lines give K E values for CC1 , and these depend upon the 4 assumed power law for the attachment cross section, Eq. (4.3.4). The slopes of the lines for both A r and Ne gives for pure CC1 the K 4 cm s 3 _ 1 E values of 5 x l 0 - 8 cm s 3 _1 and 3 . 8 x l 0 -8 for p = 0.726 and 1.38, respectively. Shirnamori and Sunagawa recently reported Chapter 4. Electron Thermalization the experimental estimate of 1 0 - 7 and Attachment cm s 3 - 1 in CCl /Ar 4 and CCU/Ne 105 for KE for pure C C I 4 , which is in reasonable agreement with the estimates from Fig. 4.10. The electron distribution functions for C C L i / A r and XQCI* = 1 0 - 4 are shown for different times in Figs. 4.11A - 4.11C. Fig. 4.11D is the steady distribution without attachment. The effect of the vibrationally inelastic collisions is clearly evident in these graphs especially at short times. The momentum transfer collisions serve at later times to broaden out the distribution function. This behavior is analogous to what was observed in previous applications. The effect of the attachment process is to remove low energy electrons and create a distribution function with a pronounced peak at high energies. This is also clearly seen in these graphs. 4.4 The steady electron distribution at long times The main interest in this chapter is the time scale for the approach of the electron distribution to a steady state as governed by elastic, vibrationally inelastic and electron attachment collisions. The results are summarized in Tables 4.1 and 4.2 in terms of the T u relaxation times. In this section, we reconsider the calculations in the recent paper by Kowari and Shizgal [126] regarding the establishment of a steady electron distribution. The conclusions stated in this previous paper were based on calculations with an incorrect cross section. In Fig. 4.12, the e-GCU attachment cross section in the low energy region is shown. The curve labelled (b) is the correct cross section as discussed in Section 4.3. The curve labelled (a) is the incorrect cross section which is a linear extrapolation of the actual cross section below 0.10 eV. Fig. 4.13 shows the time dependence of the average electron energy in an C C U / A r mixture with X c c i 4 = 5xl0 - 6 for both attachment cross sections shown in Fig. 4.12. The correct cross section, (b), gives a steady distribution whereas the one used previously, (a) appears, at first sight Chapter 4. Electron Thermalization and Attachment in CCU/Ar and CCl /Ne 4 106 Figure 4.11: Time dependence of the electron distribution function for electrons in a C C l / A r mixture. X c i = 1 0 . The time is from (A) 0.32 ns to 1.28 ns in intervals of 0.32 ns, (B) 1.6 ns to 6.4 ns in intervals of 1.6 ns, (C) 8 ns to 40 ns in intervals of 8 ns, (D) Same as (C) but without attachment. -4 4 C 4 Chapter 4. Electron Thermalization and Attachment in CCl /Ar and CCl /Ne 4 4 600 01 0.00 i I I i 0.05 0.10 i I 0.15 E(eV) Figure 4.12: Electron attachment cross section, e-CCl : (a) Cross section used by Kowari and Shizgal, (b) Cross section used in the present work. 4 107 Chapter 4. Electron Thermalization and Attachment in CCk/Ar and CCl /Ne 4 0 300 600 900 1200 1500 1800 t(ns) Figure 4.13: Time variation of the average energy, C C I 4 in Argon with all collision processes included, X c c i = 5 x l 0 . (a) Results for the attachment cross section used by Kowari and Shizgal, Fig. 4.12(a), (b) Results for the attachment cross section used in the present work, Fig. 4.f 2(b). - 6 4 108 Chapter 4. Electron Thermalization and Attachment in CCU/Ar and CCU/Ne not to yield a steady state. We here elaborate on the interpretation of the different behavior obtained with these two attachment cross sections. In fact, both cross sections give a steady electron distribution; one leads, (b), to a high steady state temperature with T t ady > Tb, that is, attachment heating occurs, whereas the other cross section, s e (a), leads to a low steady state temperature with T dy stea < Tb, that is electron cooling occurs. Since the electron density decays to zero, there is no rigorous steady state. However, previous studies [123]- [125], which neglected vibrationally inelastic collisions, demonstrated that a quasi-steady state can occur associated with the phenomenon of attachment heating. The phenomena of attachment heating as well as cooling have been discussed at length by Ness and Robson [91], Robson [138] and by Shizgal [123]- [125]. This is very similar to the discussions of high and low temperature steady states in hot atom chemistry [139], [140]. Although the attachment process removes electrons from the system, the electron distribution can attain a steady value with a steady average energy per electron. The analysis of these papers is based on the solution of the linear F P E for the problem. The electron V D F can be written as a sum of the eigenfunctions of the F P operator [123]- [125], L , and is given by, (4.4.1) n where L%l) = A ^ , x = \JEjkTb is the reduced speed, and the f n n n n coefficients are determined from the initial distribution function. The electron number density and the energy are written in the form, oo â€”X t n n=0 (4.4.3) 109 Chapter 4. Electron Thermalization and Attachment in CCU/Ar and CCU/Ne 110 The average energy per electron is then E{t) . _ E E (t) total (4.4.4) n(t)E th th a e~ + ... Xot 0 With the assumption that the electron density decays to zero faster than the energy attains some steady state, the long time dependence of the average energy per electron is, M Â« *1 + * L - < ^ > ' , Eh t (4.4.5) e Â«o Â«o provided that the higher order exponential terms (A â€” A ) are all much larger than the 0 n lowest one shown in Eq. (4.4.3). In Fig. 4.14, we show the time dependence of the average electron energy and without vibrationally inelastic collisions. Curves (a) and (b) are for a C C U / N e mixture for the attachment cross sections (a) and (b), respectively. Cross section (b) gives a steady state characterized by attachment heating whereas cross section (a) gives a low temperature steady state by attachment cooling effect. Figure 4.14A illustrates the situation with a relatively small C C U concentration. Curve (b) in the figure is for the cross section (b) and attachment heating occurs, with the asymptotic steady electron energy given by, E(oo)IE h t = 1-69. Curve (a) is for cross section (a) and there is a low temperature steady state (attachment cooling) such that E(oo)/'Eth = 0.34. In Fig. 4.14B curve (c) is for a C C U / A r mixture and yields a high temperature steady state for both attachment cross sections. At first sight, one would be inclined to interpret the results as there being a steady electron distribution for C C U / A r for both attachment cross sections and for C C U / N e only for cross section (b). However, this interpretation is not correct as there are steady distributions for all three cases in Fig. 4.14B. The steady values for the average energy, E(oo)/E h t are 5.184 and 0.00359 for curves (b) and (a), respectively. Figure 4.14: Time variation of the average energy, without vibrationally inelastic collisions. (A) CCI4 in Neon, X c c i = 4xl0~ ; (a) Results for the attachment cross section used by Kowari and Shizgal, Fig. 4.12(a), showing attachment cooling, (b) Results for the attachment cross section used in the present work, Fig. 4.12(b), showing attachment heating. (B) X c c i = 2xl0~ (a) CC1 in Neon; results for the attachment cross section used by Kowari and Shizgal, Fig. 4.12(a), with considerable attachment heating, (b) CCI4 in Neon; results for the attachment cross section used in the present work, Fig. 4.12(b). (c) CCI4 in Argon with both cross sections in Fig. 4.12. 7 4 5 4 4 Chapter 4. Electron Thermalization and Attachment in CCU/Ar and CCU/Ne 112 Fig. 4.15 shows the change in the average electron energy for C C l / N e mixtures versus 4 the mole fraction CC1 for both cross sections. For cross section (a) there is clearly a 4 cooling of the electrons whereas as for cross section (b) there is a heating. Cross section (b) leads to a hot steady state situation or attachment heating, whereas cross section (a) yields a cold steady state or electron cooling. The steady distribution is given by the eigenfunction of the Fokker-Planck operator with the lowest eigenvalue. These eigenfunctions can be interpreted as eigenfunctions of a Schroedinger equation with a potential function determined from the F P operator. This potential which determines the nature of the spectrum of the operator, depends upon the energy dependence of the momentum transfer and attachment cross sections as well as the mass of the inert gas moderator. For C C U / A r , there is only electron heating for both cross sections and we get the same results with either one. The low energy portion of the attachment cross section is not sampled. For C C U / N e , it appears that the cross section (a) gives a low temperature steady state. The time dependence of the average electron energy for cross section (a) shown in Fig. 4.14 depends upon the eigenvalues and the density and energy coefficients, a and n 6 , which are determined by the eigenfunctions. If there is a high temperature steady n state (attachment heating) the lowest eigenfunction, i/>o(x), the steady distribution is displaced to higher energies with increasing CC1 concentrations. This is shown in Fig. 4 4.16B for cross section (b). On the other hand, if there is a low temperature steady state (attachment cooling) the lowest eigenfunction, the steady distribution, is displaced towards lower energies with increasing C C U concentrations as shown in Fig. 4.16A. The accompanying density and energy coefficients can be small. Several of the higher eigenfunctions can also be concentrated at lower energies. The time variation of the average energy for Fig 4.14B, curve (a), can be interpreted as follows. With increasing C C U , the eigenvalue Ao approaches A i , so that after an initial transient the average Chapter 4. Electron Thermalization and Attachment in CCU/Ar and CCU/Ne 113 5 Figure 4.15: The attachment heating and cooling effects for electrons in C C U / N e mixtures without vibrationally inelastic collisions, (a) Results for the attachment cross section used by Kowari and Shizgal, Fig. 4.12(a). (b) Results for the attachment cross section used in the present work, Fig. 4.12(b). Chapter 4. Electron Thermalization and Attachment in CCl /Ar and CCU/Ne 114 4 x Figure 4.16: Steady electron distribution functions as the lowest eigenfunction of the Fokker-Planck operator; C C I 4 in Neon without vibrationally inelastic collisions (A) Results for the attachment cross section used by Kowari and Shizgal, Fig. 4.12(a). X c c i equal to (a) 0, (b) 10~ , (c) 2 x l 0 and (d) 4 x l 0 . (B) Results for the attachment cross section used in the present work, Fig. 4.12(b). X c c i equal to (a) 0, (b) 5 x l 0 , (c) 1 0 and (d) 5 x l 0 ~ . 4 7 - 7 - 7 - 7 4 6 - 6 Chapter 4. Electron Thermalization and Attachment in CCU/Ar and CCU/Ne 115 electron energy attains a nearly time independent value given by the ratio (&o + &i)/ao in Eq. (4.4.3). The exponential term over this time period from about 0.4 ps to almost 1.5 ps is essentially unity so that E(t)/E h t = 5.18 corresponding the the plateau value in Fig. 4.14. However, for times greater than about 1.5 ps the exponential decays rapidly and E(t)jE h t = b /a 0 0 = 0.00359 . In terms of the time dependent distribution function, one can say that it gets trapped briefly in an unstable high temperature state and then decays into the actual low temperature stable state at very long times. The preceding discussion is based on the F P E without vibrationally inelastic collision terms. The B E with these inelastic collisions is also linear and a discussion in terms of the eigenvalues and eigenfunction is possible analogous to the one just presented. The results shown in Fig. 4.13 for the system with inelastic collisions could presumably be interpreted similarly. The absence of a plateau for curve (a) is presumably due to the influence of the inelastic terms in shortening the relaxation time and hence the difference in the two lowest eigenvalues is much larger in this case than for the model problem considered. Chapter 5 Relaxation of Ion Distribution Functions by Ion-Neutral Collisions 5.1 Introduction The transport processes associated with an ensemble of ions or atoms dilutely dispersed in a background of neutral species, is an important problem in chemistry and physics [6,9] with applications to aeronomy, astrophysics, and other areas [141] â€” [149]. The approach to equilibrium of the V D F is a classical nonequilibrium kinetic theory problem which is governed by the linear B E . One method of solution of the B E , analogous to the solution of the F P E in Chapter 3, is the expansion of the V D F in the eigenfunctions of the Boltzmann collision operator. The nature of the spectrum of the collision operator in the Boltzmann equation, that is the eigenvalues and associated eigenfunctions, is of considerable interest in many aspects of theoretical gas dynamics. However, very little is known for realistic collision cross sections, and results are only available for the hard sphere (HS) cross section and the polarization interactions (where the potential varies as 1/r , and r is the ion-neutral atom separation). In this chapter, we consider a preliminary 4 study of the spectrum of the collision operator for realistic collision cross sections. We chose ion-neutral systems since these systems are of practical importance for a study of ion transport processes in the laboratory, as well as in the terrestrial ionosphere [6, 9,141]. Also, the interaction potentials for several ion-neutral systems have been published and the differential scattering cross sections can be determined. A major objective of this chapter is to consider the relaxation of an initial nonequilibrium ion V D F for a realistic 116 Chapter 5. Relaxation of Ion Distribution Functions by Ion-Neutral Collisions 117 system. Experimental techniques, primarily Doppler spectroscopy, have been used recently to probe the details of velocity distribution functions [21] â€” [24]. The distribution functions are generally probed using single frequency laser-induced fluorescence (LIF) to scan the Doppler profile at specific detection times relative to some initial pulse. Bastian et al. measured the velocity distributions for B a ions drifting in A r under the influence of an + external electric field [21]; Nan and Houston studied the relaxation of hot sulphur atoms in He and A r gases [24]; and Park, Shafer and Bersohn studied the time evolution of the V D F of hydrogen atoms in O2, N and rare gases [25]. These techniques of deducing the 2 projection of the three dimensional ion V D F in a given direction are useful in comparison with theoretical calculations [33,34]. Both theoretical and experimental knowledge of the ion velocity distributions is still quite limited. Therefore, we are interested in how the spectrum of the collision operator, and hence the nature of the ion relaxation to equilibrium, depends upon the differential cross section. Given a specific pair interaction potential, the cross section can be evaluated either with a quantum calculation or with a semi classical approach [146]. The eigenfunctions of the Boltzmann collision operator for such systems can be determined from their expansion in a suitable basis set and the numerical diagonalization of the matrix representative of the collision operator in this basis set. This requires the determination of the matrix elements of the collision operator in the chosen basis set. The evaluation of these matrix elements has been carried out for a limited number of basis functions, the most popular being the Burnett functions [6,9,150]. As men- tioned in Section 1.1, the Burnett functions are products of Legendre and associated Laguerre polynomials. The eigenvalues of the collision operator are known for only a few interatomic potentials, in particular for V(r) = Vo/r , which is a model for ion-atom 4 interactions, referred to as "Maxwell-molecules" [150]. The eigenvalues are also known Chapter 5. Relaxation of Ion Distribution Functions by Ion-Neutral Collisions 118 for hard spheres. The calculation of the matrix elements of Boltzmann collision operator with realistic cross sections for arbitrary basis functions is not straight-forward, ft is also known that these calculations are subject to considerable numerical round off errors. Alternatively, the collision operator can be written explicitly as an integral operator and reduced to matrix form with a suitable quadrature. This latter procedure is adopted in the present work. It is essentially the Quadrature Discretization Method (QDM) discussed in Chapters 2 and 3, employed there for differential operators and here for integral operators. This approach was used by Shizgal [65,151], Shizgal and Blackmore [152], Lindenfeld, Reeves and Shizgal [153] for a hard sphere differential cross section. Clarke and Shizgal [146] considered accurate differential cross sections for ( H , H ) collisions. They also studied the relaxation of an initial nonequilibrium H + + distribution function. In this work, the eigenvalues of the collision operator were calculated with the Q D M discussed in Section 2.5. However, owing to the occurrence of a sharp cusp in the kernel of the collision operator for highly anisotropic differential cross sections, an alternative numerical method referred to as the Multi-Domain (MD) method [154] was employed. The main feature of this approach is to divide the semi-infinite interval, in the first instance, into two domains where the position of the cusp defines the boundary of the two domains. To improve the accuracy of the calculation, additional domains are introduced as discussed in detail later. The convergence properties of low order eigenvalues with various differential cross sections are studied in detail. In addition, we also consider the relaxation to equilibrium of an initial non-Maxwellian anisotropic ion distribution. We are interested in the time scale for the establishment of an isotropic distribution and for energy transfer. A classical kinetic theory problem that has received considerable attention for over 50 years, is the determination of the steady ion V D F for ions dilutely dispersed in a neutral gas and under the influence of an external electric field. This problem was briefly Chapter 5. Relaxation of Ion Distribution Functions by Ion-Neutral Collisions 119 discussed in Chapters 1 and 2. It has been developed principally by Mason and Viehland [5,6], as well as by several other researchers [35,36,57,58,62]. The main difficulty for these systems is that, for large electric fields, the ion V D F is very anisotropic. The expansion in Burnett functions, referred to as the one temperature theory converges only for the smallest fields. This was later modified by choosing the temperature T that defines the reduced energy to be different from the actual bath temperature Tf, [5,9]. This approach introduces a parameter s = T/Tf, analogous to the scaling procedure in 2 Chapter 3 in the study of electron relaxation. Although this modification improves the one temperature results, the limitation of the Burnett functions in accurately representing the anisotropy of the distribution remains. The method that is presently routinely used is the three temperature approach which employs a drifting bi-Maxwellian, Eq. (1.1.5), as weight function. The basis functions for the velocity in cartesian coordinates are the Hermite polynomials. With this choice, we lose the spherical symmetry of the collision operator, and the calculation of the matrix elements associated with these basis functions becomes very complicated. Moreover, both the two and three temperature methods may not converge for large speeds because the ion V D F vanishes slower than the Gaussian weight function [57, 63]. Viehland has also introduced the Gram-Charlier weight function and associated basis functions to improve the convergence of both the two and three temperature methods [155]. The purpose of this approach is to overcome the deficiency of these previous expansion methods by building anisotropy into the weight function. However, a high price is paid both analytically and computationally in determining the matrix elements; see Ref. [155]-Appendix. The main objectives of the work described in this chapter are twofold. In the first instance, we are concerned with a study of the relaxation of some initial nonequilibrium anisotropic ion distribution. We employ a model based upon the quantum mechanical differential cross section for a specific ion-neutral potential. This.is analogous to the work Chapter 5. Relaxation of Ion Distribution Functions by Ion-Neutral Collisions by Clarke and Shizgal [146] on the relaxation of H + 120 in a background of H . The B E for this application is written explicitly in terms of the kernel of the integral operator. A quadrature procedure for integration, as discussed in Chapter 2, is used to reduce the B E to matrix form. In this way, the distribution function is determined at a discrete set of points rather than as an expansion in a basis set. However, this procedure gives poor results when the kernel is sharply peaked. In such a case, the alternative M D method is employed. The second application is the study of the relaxation of the ion V D F under the influence of an external uniform electric field. The method we employed here is similar to the one by Drallos and Wedehra [82] for the solution of the electron V D F . The ion distribution function is written in terms of vy, the longitudinal velocity component along the electric field direction, and vj_, the velocity component perpendicular to the electric field direction. Again, we use the kernel form of the collision operator and the integrals are evaluated with Gaussian quadratures. This is combined with a F D scheme for the drift term so as to determine the time evolution of the ion V D F : This algorithm is simple and straight-forward. A simple model consisting of a HS differential cross section, and with the ion and neutral mass ratio equal to unity is studied. We summarize this chapter as follows: Section 5.2 describes the relaxation dynamics in terms of the B E . The B E collision operator is given as an integral operator with a well defined kernel. The calculation of the differential cross section is also briefly discussed. Section 5.3 describes the numerical methods for the solution of the time dependent B E which includes both the Q D M and the M D approach. Section 5.4 introduces the numerical scheme for the study of ion relaxation with an applied electric field. Chapter 5. Relaxation of Ion Distribution Functions by Ion-Neutral Collisions 5.2 121 Boltzmann collision operator We consider a nonequilibrium system of ions with mass m, dilutely dispersed in a bath of atoms denoted by a subscript 1, mass rri\. We assume that the number density of atoms n is much larger than the number of density of ions, so that only ion-atom collisions need to be considered and the atom distribution is a Maxwellian at the temperature TJ,, that is /j (v) = n(m /2TrkT f e- ' M /2 l . The B E for the time evolution of the ion velocity 2kTb m v2 b 1 distribution function is given by / / [/(v')/i (vi') M where cr(<7,x) M 1 1 X (5.2.1) the differential elastic cross section, g = v â€” v i is the relative collision IS velocity, and x f^)f ^ )Hg,x)9smxd dvr IS the scattering angle in the center of mass frame. The quantities v and V i are the velocities before a collision and v ' and v^ are the velocities after a collision. The time dependent B E can be rewritten in the form (5.2.2) where the collision frequency Z(y) is given by (5.2.3) The kernel, K(y,v'), / dxsm csc (^)\v tf(v.v') x exp is explicitly given by [156] 4 X Jo Z -v' m (1 - 2 M ) 2kT M 2 2 [ 2 (5.2.4) where I is the modified Bessel function of order zero, Mi = m/(m + m-i) and M = Q 2 Chapter 5. Relaxation of Ion Distribution Functions by Ion-Neutral Collisions m i / ( m + T O i ) . In general, K(v,v') 122 depends only upon the magnitudes of the velocity v, v', and p = cos 9, where 9 is the angle between v and v'. We rewrite K(v, v') as K~(v, v', p). If we define dimensionless speeds x = V/VQ and y = v'/vo, where v = y(2kT/m), 0 can derive the dimensionless kernel K(x,y,p) r / \ / i+ r\ 3 2 from Eq. (5.2.4) [157]: 4 exp X V 7 o X s m X c s c 2r 4r 2TT / (77) e CTO 2 JO x (i + r) P 4r <7n ^ o t <!Â» x/o(( i y ^ ) o i / 7 c o t ( | ; (5.2.5) x where T = mi/m, we A â€” 4n<7oy27rA;T/m, cr = l A , 7 = x + y â€” 2xyp, a = 2xy(l â€” p.)2., 2 2 0 and E = (&T7csc (^))(l + r ) / ( 4 r ) . We can transform from x to u with the relationships 2 ,2 =_ 7cot (^) and -udu = â€” d sinxcsc (|). 2 Eq. (5.2.5) becomes 4 X K(x,y,p) = ( i - ry exp < â€” [x + 4r 8TTT / 5 2 I X / where s g(E,x) 4r r 2 - i 2r i + r / ((-^)au)du, 0 (5.2.6) 0o = 2 t a n - ( 7 / / u ) and E = kT(<y + u ) ( l + T)/(4r). If we choose T = 1, we can 1 X uexp ;i + r) + y) + 2 1 2 2 further reduce the kernel to a simpler expression, K(x,y,p) = A(â€”â€”] exp(â€”x ) / uexp(â€”u )â€”-â€”''-^-Io(au)du. 2 XTT^/j/ 2 JO (5.2.7) (To The ion V D F can be written as an expansion in Legendre polynomials [152] f(v,t) = ZMv,t)P (p), z 1=0 i (5.2.8) Chapter 5. Relaxation of Ion Distribution Functions by Ion-Neutral Collisions 123 where p = cos 6 and 8 is the polar angle between v and some fixed axis in velocity space. We can decompose the dimensionless kernel K(x,y,p) by the Â£ th component K~i(x,y), where K (x,y) = 2TT j\ K(x,y,p))P (p)dp. (5.2.9) 2 t e We can introduce the dimensionless time t' = tA and reduce the time dependent B E to a set of uncoupled differential equations, that is roo 1_ dft / A Jo dt' .K (x,y)f ( ,t')dy-Z(x)f (x,t>) Â£ i y (5.2.10) Â£ We consider (Li ,He) as a typical system. The interaction potential proposed by Kout+ selos, Manson, and Viehland ( K M V ) [158] is of the form V(r) = v A\ exp(- b\r, B\ exp( 0 P -h cu (ci .+cÂ£ ) nd {Cj + iap nd ct ) p ' h 1 = r > 1.28r [ exp[â€”(1.28r â€” l ) ] 2 m m (5.2.11) otherwise. In atomic units, the parameters in Eq. (5.2.11) are v = 0.1135, A\ = 146.98, a\ = 1.5024, 0 5a = 70.198, b = 1.4041, p = 0.7185, C = 0.6916, C l nd x Cl nd = 5.307, Ct 4 p = 2.004, and r m 6 = 1.2217, Cf sv = 0.2696, = 3.64. The differential cross sections, a(E,x), were determined from a quantum mechanical calculation of the phase shifts for this pair interaction. For the scattering angle Xi the scattering amplitude F(x) c a n expressed in terms of the phase shift Sg of the partial D e waves of the angular momentum quantum number t oo (X) F = \ E( ^ + !)^( 2 cos x)e iSt sin S lt (5.2.12) i=o where k is the wave number and Sg are calculated by the Wentzel, Kramers, and Brillouin fc ( W K B ) method [159]. The differential cross sections at energy E are given by <r(E,x) = \F(x)\\ (5.2.13) Chapter 5. Relaxation of Ion Distribution Functions by Ion-Neutral Collisions 124 where the scattering amplitude depends upon the interaction energies between collision pairs. The total cross section <r (E) is defined as Jot v (E) tot = 2TT Â£ <x(E, ) sin dxX (5.2.14) X We use the Ffenyey-Greenstein (HG) [56] function given by <7(E,x) = <T i o ( (E)(l-/)/(l+^-2^cosx) 3 / 2 (5.2.15) with g â€” 1 â€” e to fit the quantum differential cross sections. For a given energy, as e â€”> 1, the anisotropy of the differential cross section decreases and tends to an energy dependent "hard sphere" cross section. Unlike the classical case, the quantum differential cross sections are finite at zero scattering angle. Varying the e values enable us to study the effects of the differential cross section on the eigenvalue spectrum of the collision operator. The differential cross section for (Li ,He) collision with the K M V interaction poten+ tials in Eq. (5.2.14) for two energies are shown in Fig 5.1. The dashed curve is the H G fit given by Eq. (5.2.15). The parameters for the energies 0.01 and 0.1 eV of Figs. 5.1A and 5.IB are a (E) tot = 235A and 108A, and e = 0.0057 and 0.0039, respectively. In the preliminary study of the role of the differential cross sections on the spectrum of the collision operator and on the relaxation behavior, it was sufficient to use the simple H G fit to the cross sections which duplicates the small angle scattering very well. The detailed oscillations of the quantum scattering cross sections are not expected to be important. 5.3 5.3.1 Relaxation of anisotropic ion distribution function Quadrature Discretization Method (QDM) We solve the Boltzmann equation, Eq. (5.2.10), numerically using the Q D M . This method has been described in detailin Section 2.5 [65,146], so we only give an outline here. It is Chapter 5. Relaxation of Ion Distribution Functions by Ion-Neutral Collisions 125 LU LU Figure 5.1: Differential cross section for (LI ,He) collision with energy (eV) (A) 0.01 and (B) 0.1 ( quantum mechanical cross section and H G fit cross section). + Chapter 5. Relaxation of Ion Distribution Functions by Ion-Neutral Collisions 126 convenient to consider a reduced speed defined by x = The where v â€” (2kT/mi) ! . 1 V/VQ 0 2 basis of the QDM is to use a set of convenient orthogonal polynomials to approximate the distribution function. A l l quantities are then evaluated at discrete points (the zeros of the N th order polynomial), and integrals are evaluated using Gaussian quadratures. Here we use the "speed" polynomials, defined as the set of orthogonal polynomials R {x) n which satisfy f Â°Â° 0 w(x)R (x)R (x)dx n = 6 m where w(x) = x e~ . 2 nm x2 With the change of variable to reduced speed, the integration over y on the right hand side of Eq. (5.2.10) can be performed with a quadrature of the form / Â°Â° w(y)g(y)dy ~ 0 J2f=i j9( j) w x where Xj and WJ are the corresponding quadrature points and weights. Eq. (5.2.10), written as a quadrature sum is then of the form [65] QÂ£ (5.3.1) = 2^ ijM J> )> X M T where the subscript i denotes evaluation at the point X{ and I<e(xi, XJ) - Z(xi)S /w(Xj). i: â€¢ (5.3.2) In order to conserve particle number accurately, we also calculate the collision frequency Z(x{) with the same quadrature [65], that is, Z( ) Xi = iZwjK(x^Xj)f (x )/(w(Xj)f (x )). M M 3 t (5.3.3) 3=1 The u integral of the kernel in Eq. (5.2.6) is evaluated by using a standard Simpson's rule integration and the p integration in Eq. (5.2.9) is performed by using a Legendre quadrature. The solution of Eq. (5.3.1) is N /H*M0 = exp(-AfV) I > i f Â« ? (5.3.4) 3=0 with = f;bW)(- ))l / (x ,0), 1 i k=0 / jk j ,(5.3.5) Chapter 5. Relaxation of Ion Distribution Functions by Ion-Neutral where U is a matrix which diagonalizes M , that is, U _ 1 Collisions 127 M U = A (where A is the diagonal matrix of eigenvalues with A Â°' = 0 < AJÂ°' < A Â°' < . . . ) . The term uffl is the 0 i th 2 component of the eigenfunction with eigenvalue A Â°' = 0, i.e. a Maxwellian. Thus, 0 / decays to a Maxwellian as t' â€”> co with a relaxation time of the order of l / A ^ L The spectrum of the collision operator determines the detailed time evolution of the distribution function. We consider the relaxation of an arbitrary initial anisotropic, non-Maxwellian ion distribution function / ( j , ,o A c ) = (A + i)'Â«p[-a(g,-Â»')]_ ( 5 . 3 6 ) where p is a small integer that determines the extent of the initial anisotropy; EQ is the initial reduced energy; a is a factor that specifies the width of the distribution, C is a normalization factor such that 27r / Â°Â° j 0 \ x f(x, 2 p)dpdx = 1 and p = cos (9, where 8 is the angle between x and z direction. Since the distribution function is expanded in Legendre polynomials, we can describe the initial coefficients fg(xi,0) in Eq. (5.3.5) as N fi(xi,0) = Y^Wif(xi,pi,0)Pe(pi), (5.3.7) where W{ and pi are the weights and points of the Legendre integration. Eq. (5.3.4) together with Eq. (5.2.8) provides the solution for the time dependent Boltzmann equation. The coefficients a!p in Eq. (5.3.5) are evaluated from Eqs. (5.3.6) and (5.3.7). 5.3.2 The Multi-Domain (MD) approach ft is well known that the numerical solution of Eq. (5.2.10) is made difficult by the occurrence of a narrow cusp in the kernel at x = y. This was discussed at length by Lindenfeld and Shizgal [154] where a Multi-Domain (MD) approach was introduced. In Chapter 5. Relaxation of Ion Distribution Functions by Ion-Neutral Collisions 128 order to accurately integrate over the cusp, for fixed x the semi-infinite interval in y is divided into two domains [0,x] and [x,oo]. This ensures that one grid point coincides with the position of the cusp. Since the kernel is very strongly peaked, additional subdomains are introduced on either side of x. In each subdomain, a Legendre quadrature is used. The accuracy of this procedure is checked by calculating the collision frequency from Eq. (5.2.3) and also directly in terms of the total cross section. A fine grid is chosen for fixed x and an interpolation is then used to relate the fine grid to the speed points. The M D method divides the integration range into many separate intervals with one of the interval end points at x. With the choice of the subintervals within the integration range, a set of M integration points Zk(xi) and weights Wk(xi) are defined for each speed X{. At each of the Nxi points, the kernel integration over z leads to a quadrature approximation of the form M K (x, z)ft(z)dz = w (xi)Kt(xi, e z (xi))ft(z {xi)). k " k (5.3.8) k =i fc This leads to MN unknown quantities ft(z (xi)) k instead of just N. In order to reduce Eq. (5.3.8) to a set of N equations in N unknown quantities fe(xj), an n th Lagrange interpolation procedure is employed, that is, N )( )l n M*k(xi)) = E iM*i))M*i), 3 D (5-3-9) i=i where Df\z {x ))= k ft i=i & % * - X l Z X ' {Xl)) ' (5-3.10) J X With Eqs. (5.3.8), (5.3.9), and (5.3.10), we can replace the elements M~\f in Eq. (5.3.2) with = k (xi,xj)-Z(xi)8ij, 1 e where M Ke(xi,Xj) = Yj k{xi)Kt(xi,z (xi))Df\z (xi)), k=i w k k (5.3.11) (5.3.12) Chapter 5. Relaxation of Ion Distribution Functions by Ion-Neutral Collisions 129 and Z(xi) is the collision frequency determined with M D so as to maintain detailed balance. Although the M D method provides a better approximation to the kernel integration, it yields a nonsymmetric matrix M which gives complex eigenvalues and eigenfunctions. As a result, we can not apply the M D procedure for solving the time dependent B E as we did with the Q D M . However, when the kernel is severely peaked, the Q D M fails to approximate the kernel integration properly. 5.3.3 Results and discussion There appears to be only one previous study of the relaxation of ions in a bath of neutral atoms. This was the study of ( H , H ) relaxation by Clarke and Shizgal [146]. The system + chosen for study in the present work (Li ,He) is more classical than ( H , H ) in terms of + + the nature of the differential cross section. In particular, the scattering cross section for (Li ,He) is far more strongly peaked in the forward direction than ( H , H ) . The present + + work is directed towards a study of the effect of the nature of the differential cross section on the relaxation behavior. We therefore consider the H G fit to (Li ,He) cross section + chosen as a model system and vary the parameter e which controls the cross section at small angles. In order to simplify the problem, we choose the mass ratio as unity. We have used both speed quadrature (QDM) and the M D approach to calculate the eigenvalues of the collision operator (with only I = 0) for an H G model cross section with e in the range 0.01 to 0.9. The results are summarized in Tables 5.1 and 5.2. For e â€”> 0.01, the differential cross section becomes sharply peaked in the forward direction and the kernel is characterized by a narrow cusp. The isotropic kernels, KQ, with e = 0.9 and 0.1 are plotted in Fig. 5.2. Decreasing e increases the sharpness of peak at the cusp where x = y. The results in Table 5.1 show that the Q D M with speed quadrature points provides convergence of the first seven eigenvalues up to 3 significant figures, Figure 5.2: The isotropic kernel K (x,y) with e equal to (A) 0.9 and (B) 0.1, for the H G cross section. 0 Chapter 5. Relaxation of Ion Distribution Functions by Ion-Neutral Collisions for e = 0.9 and e = 0.5. 131 For these values of e, the differential cross section is not strongly peaked in the forward direction and the kernel is also not strongly singular. However, the convergence begins to deteriorate for e < 0.1 and is rather poor for Ay ' and 0 e = 0.01. B y contrast, the M D results in Table 2 provides remarkable convergence for the larger eigenvalues and larger e. For example, with e = 0.05, the first three eigenvalues are converged to 3 significant figures with the M D and 20 quadrature points, whereas the speed points fail to give convergence to 2 significant figures. Moreover, with e as specified in Table 2, most eigenvalues converge up to 5 significant figures with M D and 40 quadrature points. We have also considered the relaxation of an initially anisotropic non-Maxwellian ion distribution function given by Eq. (5.3.6) with p = 2, E = 700 and & = 0.00045. We 0 show the relaxation of the anisotropic distribution in Figs. 5.3 and 5.4 for the H G fit cross sections with e = 0.9 and 0.1, and <r (E) = 4 9 . 8 E ioi - 03 3 7 and a mass ratio of unity. The total cross section a , is determined by a power law fit to the total cross sections tot which is calculated from the realistic differential cross sections. In the figures, x is the reduced speed and \x = jx defined in Section 5.3.1. Since the kernel is not severely peaked for these values of e, we employ the speed quadrature for this calculation by using 80 speed and 80 Legendre quadrature points. For e = 0.9 in Fig. 5.3, the relaxation of the isotropic and anisotropic portions of the distribution appear to occur on the same time scale. Apparently, for this e, there is a large exchange of energy on collision which results in a depletion of the high energy particles and a growth of distribution at low energy. The growth of the distribution at low speed gives a bimodal form during relaxation. This is similar to the results by Blackmore and Shizgal [152] on the relaxation of V D F for HS differential cross section for unit mass ratios. For a smaller e = 0.1 in Fig. 5.4, bimodal distributions do not occur. There is a rapid degrading of speed followed by a slow removal of the anisotropy. The results show persistence of the anisotropy due to the Chapter 5. Relaxation of Ion Distribution Functions by Ion-Neutral Collisions N â€¢ Ar 4Â°Â» \(Â°) Ai 58.1287 58.2494 58.2649 58.2691 60.1844 60.3132 60.3286 60.3327 61.6099 61.7571 61.7724 61.7761 41.0835 41.7793 41.8906 41.9241 44.9949 45.9119 46.0567 46.1002 50.5318 51.8977 52.1031 52.1643 8.71527 9.71481 10.0636 10.2276 9.79740 11.1277 11.5921 11.8120 11.3940 13.4828 14.2016 14.5467 4.22622 4.77534 4.99171 5.10468 4.75274 5.47923 5.76478 5.91426 5.52435 6.65524 7.09148 7.32205 0.81706 0.92788 0.97379 0.99906 0.91894 1.06522 1.12558 1.15878 1.06787 1.29479 1.26215 1.43719 0) 132 A<Â°> e = 0.9 20 40 60 80 34.1333 34.1647 34.1688 34.1700 47.3967 47.4681 47.4777 47.4804 54.2967 54.3995 54.4132 54.4170 e = 0.5 20 40 60 80 18.9174 19.0311 19.0486 19.0537 29.1069 29.3911 29.4364 29.4499 36.0013 36.4836 36.5612 36.5844 e = 0.1 20 40 60 80 3.48729 3.67725 3.74150 3.77123 5.71013 6.13674 6.28342 6.35178 7.38403 8.08196 8.32430 8.43758 e = 0.05 20 40 60 80 1.68379 1.79152 1.83375 1.85616 2.76307 3.00147 3.09513 3.14441 3.57783 3.96409 4.11618 4.19569 e = 0.01 20 40 60 80 0.32511 0.34709 0.35623 0.36146 0.53385 0.58226 0.60237 0.61362 0.69155 0.76971 0.80215 0.82008 Table 5.1: Convergence properties of eigenvalues (in units 4ncr (27rA;T/ra) / ) of the collision operator (only with Â£ = 0) with Q D M for H G fit cross sections where e is range from 0.01 to 0.9. 1 0 2 Chapter 5. Relaxation of Ion Distribution Functions by Ion-Neutral Collisions N A<Â°> \(0) 2 \(o) 133 x(o) A e = 0.9 20 30 40 34.1710 34.1710 34.1710 47.4828 47.4828 47.4828 54.4202 54.4203 54.4203 58.2727 58.2727 58.2727 60.3369 60.3361 60.3361 61.7410 61.7744 6.1.7783 41.9553 41.9554 41.9554 46.1411 46.1406 46.1406 52.2159 52.2201 52.2203 10.5250 10.5288 10.5288 12.2128 12.2170 12.2169 15.1977 15.1803 15.1825 5.39926 5.41289 5.41300 6.30600 6.32672 6.32637 7.84520 7.95187 7.96425 1.10040 1.11286 1.11334 1.29370 1.31101 1.30963 1.54155 1.61815 1.66209 e = 0.5 20 30 40 19.0583 19.0583 19.0583 29.4625 29.4625 29.4625 36.6062 36.6062 36.6062 e = 0.1 20 30 40 3.81969 3.81971 3.81971 6.47212 6.47204 6.47204 8.64385 8.64331 8.64331 e = 0.05 20 30 40 1.90961 1.90966 1.90966 3.27136 3.27110 3.27109 4.41038 4.40834 4.40834 e = 0.01 20 30 40 0.38258 0.38263 0.38263 0.66259 0.66233 0.66232 0.90223 0.90019 0.90019 Table 5.2: Convergence properties of eigenvalues (in units 4n,0- (27rÂ£:T/m) ' ) of the collision operator( only with Â£ = 0) with M D approach for H G fit cross sections where e is range from 0.01 to 0.9. 1/ 2 o Figure 5.3: Relaxation of the anisotropic ion distribution function for Henyey-Greenstein model cross section (e = 0.9) with t' equal to (A) 0.0, (B) 0.02, (C) 0.03 and (D) 0.05. Figure 5.4: Relaxation of the anisotropic ion distribution function for Henyey-Greenstein model cross section (e = 0.1) with t' equal to (A) 0.0, (B) 0.1, (C) 0.2 and (D) 0.4. Chapter 5. Relaxation of Ion Distribution Functions by Ion-Neutral Collisions 136 strong forward scattering. The time evolution of the reduced energy is given by E(t') = / (5.3.13) f(v,t')(mv /2)dv 2 and since it is determined by the isotropic portion of the distribution, only \ \ contribute to the relaxation. B y contrast, the decay of the directed velocity given by, (5.3.14) depends upon the anisotropic portion of the distribution and A- (Â£ 7^ 0) eigenvalues t contribute to the relaxation. The decay of these average quantities in Fig. 5.5 with several e illustrate the effects of e on the relaxation process. It is clear that both the average energy and directed velocity relax much faster with increasing e, that is, with an decrease of the forward scattering. The dependence of the relaxation behavior with e is similar to the variation of mass ratio [152], for the HS cross section. 5.4 5.4.1 Numerical solution of the time dependent B E with applied electric field Numerical algorithm As discussed in Section 5.1 and briefly in Chapter 2, the determination of the ion distribution for ions in a heat bath of atoms under the influence of a steady electric field is a classic kinetic problem. For the most part the interest is towards the determination of the ion mobility and not the distribution function. Here we are interested in the development of new numerical methods towards the determination of the ion V D F which would also yield the mobility. The main problem is the need to determine distribution functions far from Maxwellian both in terms of the energy as well as the large anisotropics that can occur. As mentioned in Chapter 2, under the influence of an applied electric field, the tail Figure 5.5: (A) Time variation of the average ion energy. (B) Time variation of the directed ion velocity with e equal to (a) 0.1, (b) 0.5 and (c) 0.9. Chapter 5. Relaxation of Ion Distribution Functions by Ion-Neutral Collisions 138 of ion V D F tends to fall off much slower than any Gaussian [57]. Hence the distribution function defined with a Gaussian weight diverges. In addition, it is very difficult to find a weight function which is not only suitable for the anisotropic ion V D F , but also simple enough for determining the matrix elements of the collision operator. We now consider an applied electric field E along the z axis. The presence of the electric field accelerates the ions and creates an anisotropy in that direction. In order to deal with the highly anisotropic ion V D F , we rewrite the kernel in terms of V||, the velocity component || to the electric field, and vj_, the velocity component _L to the electric field [82,160]. W i t h the introduction of the reduced speed xy = v^/v and xÂ± = vÂ±/v 0 0 where VQ â€” (2kTlm\Y^ , we write the time dependent B E as 2 df(x ,x ,t') l} L + dt' mv Aj\n) dx~, 0 = (5.4.1) J where t' = tA is the dimensionless time, n is the density of the neutrals, and J' = J A is the dimensionless collision operator given by fOO OO / / -oo J0 K(x\\,x y\\, yÂ±)fi{y\\,yÂ±,t)y dy dy\\ u x - Z(x\\, x )h{x\\, x x x ,t). x (5.4.2) The kernel K(x\\,xj_,y\\,yÂ±) is derived from Eq. (5.2.6) and described as â€¢ ^_ A(i + ry ^ H ' ^ l b ^ 8*T^(xsine )(y in9 )Js 1 x exp 1 + rS 2 1 - 4T X 2 2T 2 foo S 2r l + rv l + r 2 ( i + \g(E,x) / uexp â€”â€”u Jo V 41in the/ form â€¢ crof and the collision frequency is written 0 oo / -oo roo j JO K(x\\,x ,y ,y ) Â± {] _ x + (l _ ( f l l Â± ^ 2 r i+r ^ + \j io â€” - â€” a u i d u , r T 2 ds ^ 0 + .2 ^ V 2 jM ly\ * Jl \X) (z ^ A (5.4.3) / y dyÂ±dy Â± h (5.4.4) Chapter 5. Relaxation of Ion Distribution Functions by Ion-Neutral Collisions where fN is a Maxwellian; a = 2xy(\- i )2/S; u = S cot ( /2); S = 2 2 t 2 2 139 (x +y -2xyp) 2 2 X 2 ] So and Si are the limits at \i = 0,1 respectively; where [i is the cosine of the angle between v and v'; x = x\ + xjj and 61 is the angle between x and xy; y = y\ + y and 9 is the 2 2 2 2 angle between y and vn; a = x Â° J J II' ^. \ 0 y ) cot#i cot6> ; and a = ~â€”â€¢ \ â€¢ a 2 2xy s i n Â» i sin#2 1 2 ' a â€¢ The S 2a;i/sinffisin#2 integral in Eq. (5.4.3) is performed by using a Chebyshev quadrature. We reduce the integro-differential equation in Eq. (5.4.1), by introducing two sets of grids for xy and x\_. We choose Hermite and speed quadrature weights and points for and respectively. The integrations over these variables in the collision operator in Eqs. (5.4.2) and (5.4.4) can be represented by quadrature sums. The derivative with respect to x\\ in the electric field term in Eq. (5.4.1) is evaluated as a finite difference with the Hermite quadrature points. The time derivative is approximated by a forward finite difference with a uniform grid for the time. This scheme yields the discretized B E given by rn+l _ fn / pF \ f n â€” f N n N (5.4.5) where w(x) = x e~ 2 x2 and w'(x) â€” e~ . x2 The subscripts i, k denote the || velocity compo- nent represented by Hermite quadrature whereas the subscripts j, I denote the J_ velocity component represented by speed points. The discretized representation of the B E given by Eq. (5.4.5) is clearly very straight-forward to obtain. It simply requires the evaluation of the kernel at the chosen set of grid points. The choice of grid points is important and different choices can be easily implemented. A n alternate choice to Hermite quadrature could for example be the quadrature based on the B G K solution discussed in Chapter 2. B y contrast, the moment methods by Viehland and Mason [9,6] and by Viehland [155] are restricted to particular basis functions and a change in the basis sets entails a Chapter 5. Relaxation of Ion Distribution Functions by Ion-Neutral Collisions 140 reevaluation of the matrix elements. However, the evaluation of the Kij,k,i i Eq. (5.4.5) n requires a lot of computational time and memory space. Fortunately, we only need to do this once. After the Ki,j k,i is calculated, we store it as a four dimensional array and ; reuse it for each time iterative step. In order to solve the kernel integration with the Gaussian quadrature, we are forced to use the non-uniform Hermite points as the grid for the speed derivative. This will cause some discrepancies in the calculation, in particular when the distribution is highly anisotropic. 5.4.2 Results and discussion The development of the Doppler spectroscopy has stimulated the need for an efficient methodology for the calculation of the ion distribution function. However, as mention earlier, most of the works are directed primarily towards the ion mobility but not the ion distribution. Here, we have developed a simple method in which ion V D F can be calculated. With a given initial distribution, we can evaluate the collision terms in Eq. (5.4.5) with Gaussian quadrature. We then multiply each of the collision terms by A t ' and add to the corresponding distribution function array element. The result is a new distribution function at time t' + A t ' . Repeating the above procedure with the newly calculated distribution, we can obtain the distribution function at any t'. A simple model with the HS differential cross section and unit mass ratio is used to study the ion relaxation with an external applied electric field, E/n. Although this is a simple model, it provides a preliminary study for the ion relaxation process with an applied electric field. Fig. 5.6 shows the time evolution of the ion velocity distribution with an initial bi-Maxwellian function e~ ^e~ ^ for E/n x x = 0.5 T d ( 1 0 -17 V c m ) . This 2 calculation is performed by using 30 Hermite, 30 speed and 60 Chebyshev points. The uniform time step A t ' is chosen as 0.03. At the beginning, the distribution function is Chapter 5. Relaxation of Ion Distribution Functions by Ion-Neutral Collisions 141 A outermost contour line is at 0.05. The second outermost contour line is at 0.1. After that each contour line is separated equally by 0.1). Chapter 5. Relaxation of Ion Distribution Functions by Ion-Neutral Collisions E/n (Td) 0.10 0.25 0.50 0.75 1.00 w(oo) 142 % Deviation Present results Viehland [9] 0.600 0.151 0.303 0.448 0.579 0.615 0.153 0.299 0.435 0.559 2.3 1.3 1.3 3.2 3.1 Table 5.3: Reduced thermal drift velocity w(oo) for hard sphere cross section with unit 1/2 mass ratio at different E/n, where w = w/vo and vo = (2kT/m) symmetric on both positive and negative values of #y. However, it becomes anisotropic and shifts gradually towards higher xy as time proceeds. For a longer time, the ion distribution function reaches a steady state. The relaxation of the average drift velocity w(t') = / fvdv is also plotted in Fig. 5.7A in reduced units, w = w/vo, for different initial energies. It is clear that all curves in Fig. 5.7A reach the same thermal drift velocity and independent of the initial energy. Moreover, the relaxation of the drift velocity for field strengths up to 1 T d are also shown in Fig. 5.7B. Increasing the field strength allows the ions to drift faster along that direction and to reach the thermal drift velocity in a shorter time. From this figure, all curves approach nicely towards the thermal drift velocity and agree very well with the results reported by Viehland as shown in Table 5.3 [9]. Calculation with this model gave bimodal distributions for E/n > 0.5 T d which appears to be nonphysical. At present, we do not have an explanation of this effect and suspect that it is due to some numerical aspects of the method. Chapter 5. Relaxation of Ion Distribution Functions by Ion-Neutral Collisions 143 Figure 5.7: (A) Time variation of the average reduced drift velocity with E/n = 0.5 T d (10~ Vcm ) for different initial velocity equal to (a) 2, (b) 1 and (c) 0. (B) Time variation of the average reduced drift velocity for several E/n in T d (a) 1, (b) 0.75 and (c) 0.50. 17 2 Chapter 6 Summary and Outlook The traditional method of solution of the B E of kinetic theory involves the expansion of the V D F in a basis set of orthogonal polynomials. In this thesis, we have introduced an alternate approach which involves the representation of the distribution function on a grid of quadrature points defined by a weight function. With the Stieltjes procedure developed by Gautschi, quadrature points based upon basis sets for arbitrary weight functions can be calculated. The criteria for the choice of the basis functions is the ease of their use in a particular application and the rate of convergence of the solution with the size of the basis set. In Chapter 2, we studied the convergence properties of the solution of the B E for a simple model of the collision operator developed by Bhatnagar, Gross, and Krook [53] and used Hermite polynomials as basis functions. This described the distribution function of ions dilutely dispersed in a background of neutrals under the influence of an applied electric field. In the presence of an applied electric field, the results obtained from the expansion of the ion V D F in a set of Hermite polynomials diverges. These agree with those by Skullerud [57] and Ferrari [63] who suggested that this should be a general behavior for the solution of the B E based upon a Gaussian weight function. The convergence properties of the ion V D F with an alternate set of polynomials based upon the B G K solution of the B E as weight function. The results obtained from the expansion of the ion V D F in these polynomials converge very well. This study serves to demonstrate the significance of the choice of the weight function in a given application. The polynomial set associated with the non-classical weight functions based upon the 144 Chapter 6. Summary and Outlook 145 B G K solution is determined with the Stieltjes procedure developed by Gautschi [51,52]. It has been shown that Gautschi's method works very well as long as a suitable number of intervals and quadrature points is chosen. Moreover, Gautschi's method was also used for the generation of polynomials with respect to some weight function in electron-neutral collisions. In Chapter 3, we reconsider the relaxation of electrons in rare gas moderators under the influence of an applied electric field. The B E is approximated by the F P E and only two terms in the Legendre expansion of the angular part of the distribution function are sufficient. A numerical method known as the Q D M is employed for solving the F P E eigenvalue problem. The previous work by Shizgal and coworkers [43]- [45] used the speed points and weights defined by the weight function w(x) â€” x e~ , whereas in the 2 x2 present work, the quadrature weights and points are defined in terms of the Davydov distribution as the weight function. The Davydov distribution is the steady distribution for electrons in a gas moderator at some given electric field strength. In this way, the discretization depends upon the electric field and the electron moderator momentum transfer cross section. The orthogonal polynomials and associated weights and points with respect to the Davydov distribution are generated with Gautschi's procedure. A detailed comparison is made with the Q D M with speed points, F D technique and LI approach with the Legendre polynomials to define the derivative operator in the discrete representation. The present results show that the Q D M with Davydov points provides the most rapid convergence of all the methods. In Chapter 4, we consider the study of electron-inert gas elastic collisions. We also consider the study of a mixture of either A r or Ne and C C U , including elastic electroninert gas and electron-CCI4, vibrationally inelastic electron-CCI4 collisions, and electron attachment to C C I 4 . A finite difference technique for both energy and time is employed Chapter 6. Summary and Outlook 146 for the solution of the B E . The time dependent electron number density, electron average energy, energy relaxation rate coefficients, attachment rate coefficients and energy relaxation times are calculated with a wide range of C C I 4 concentrations. Comparison of the theoretical results with experiment is made. The nature of the spectrum of the collision operator in the B E , that is the eigenvalues and associated eigenfunctions, is of considerable importance in many aspects of theoretical gas dynamics. One of our interests in this thesis is to study the eigenvalue spectrum of the collision operator for a realistic cross section of ion-neutral collisions. With the interaction potential for a given collision pair, the differential cross section can be evaluated from a quantum or from a semi-classical calculation. In this preliminary study on the relaxation behavior, we use a simple H G fit to the cross section which duplicates the small angle scattering very well. We have used both a Q D M and M D approach with speed points to calculate the eigenvalues of the collision operator. The present results indicate that the M D approach provides extremely rapid convergence of the eigenvalues, while the Q D M works well only when the model cross section is not strongly peaked. With the Q D M , we are able to study the time relaxation of an initially anisotropic, non-Maxwellian distribution for different H G model cross sections. For a differential cross section that is not strongly peaked in the forward direction, the relaxation of the isotropic and anisotropic portions of the V D F appears to occur on the same time scale which gradually gives a bimodal form during relaxation. By contrast, for strong forward scattering, the bimodal distribution does not occur and there is a persistence of the anisotropy. In addition, we also consider the relaxation of the ion V D F under the influence of an external electric field. The ion V D F is written in terms of the longitudinal velocity component along the E direction and the velocity component perpendicular to the E direction. The collision operator of the B E is described in kernel form and the integrals Chapter 6. Summary and Outlook 147 are evaluated by a Gaussian quadrature rule. This is combined with a F D scheme for both time and speed derivatives so as to determine the time evolution of the ion V D F . However, this model only works well for a small range of electric field. The present work has concentrated on the development of kinetic models for the solution of the V D F for both electron/neutral and ion/neutral systems. However, some of these models are very simple and thus this work is in some instances a preliminary study towards a treatment of real systems, in particular for ion-neutral collisions. For a complete understanding of the relaxation of the V D F , future work will be mainly concerned with: (i) a detailed analysis of the relaxation of the electron V D F for larger field strengths. Owing to the restriction of the momentum transfer cross section, the Q D M with Davydov only works well up to 0.5 Td; (ii) a detailed analysis of the relaxation of the electron V D F for small time. As mention earlier, the Q D M works well for the long times, but not for short times; (iii) a detailed analysis of the steady state at long times for the electron-CCU attachment process; (iv) a detailed analysis of the relaxation of the ion V D F with more realistic differential cross section models as well as the exact quantum cross section; (v) a detailed analysis of the relaxation of the ion V D F under the influence of the applied electric field, for large electric field strengths. Bibliography S. Chapman and T . G . Cowling, The Mathematical (Cambridge Univ. Press, New York, 1964). Ferziger J.H. and Kaper H.G., Mathematical (North Holland, Amsterdam, 1964). Theory of Non-uniform Theory of Transport Gases Processes in Gases of Ions in Gases in Gases (John I. Kuscer and N . Corngold, Phys. Rev. A 139, 981 (1965). I. Kuscer and N . Corngold, Phys. Rev. A 140, 5 (1965). E. W . McDaniel and E. A . Mason, The Mobility (John Wiley, New York, 1978). E. A . Mason and E. W . McDaniel, Transport Wiley, New York, 1988). and Diffusion Properties of Ions G. H . Wannier, Bell System Tech. J. 32, 170-254 (1953). T. Kihara, Rev. Mod. Phys. 25, 844-852 (1953). L. A . Viehland and E. A . Mason, Ann. Phys. 91, 499-533 (1975). L. A . Viehland and E. A . Mason, Ann. Phys. 110, 287-328 (1978). B. D. Shizgal and D. Napier, Physica A 223, 50-86 (1996). J. C. Light and J . Ross and K . E. Shuler, Kinetic Processes ed. A. Hochstim (Academic Press, New York, 1969). M . William, Mathematical don, 1971). Methods in Particle K . M . Case and P.F. Zweifel, Linear ma, 1967). Transport Transport Theory in Gases Theory and Plasmas (Butternorth, Lon- (Addison - Wesley, Reading S. Loyalka, Phys. Fluid A 1, 612 (1989). J. Barrett and B . Shizgal, J. Chem. Phys. 91, 6505 (1990). R. W . Schunk, and J . C. G. Walker, Planet. Space Sci. 20, 2175-2191 (1972). 148 Bibliography 149 [18] J. St.-Maurice and R. W. Schunk, Planet. Space Sci. 21, 1115-1130 (1973). [19] J. St.-Maurice and R. W . Schunk, Planet. Space Sci. 22, 1-18 (1974). [20] J. St.-Maurice and R. W . Schunk, Planet. Space Sci. 25, 243-260 (1977). [21] M . J . Bastian, C. P. Lauenstein, V . M . Bierbaum, and S. R. Leone, J . Chem. Phys. 98, 9496-9512 (1994). [22] Y . Matsumi, M . Shamsuddin, Y . Sato, and M . Kawasali, J. Chem. Phys. 101, 96109618 (1994). [23] Y . Matsumi, A . M . Sarwaruddin Chowdhury, J. Chem. Phys. 104, 7036-7044 (1996). [24] G. Nan and P. L. Houston, J. Chem. Phys. 97, 7865-7872, (1992). [25] J . Park, N . Shafer, and R. Bersohn, J. Chem. Phys. 91, 7861-7871 (1989). [26] M . Hogan and P. Ong, J. Phys. B: At. Mol. Opt. Phys 21, 1417-1428 (1988). [27] L. A . Viehland and A . S. Dickinson, Chem. Phys. 193, 255-286 (1995). [28] L. A . Viehland, A . S. Dickinson and R. Maclagan, Chem. Phys. 211, 1-15 (1996). [29] Lie-Svendsen, M . H . Rees and K . Stamnes, Planet. Space Sci., 39, 929-943 (1991). [30]' T. L. Killen and P. B . Hays, J. Geophys. Res. 88, 10163-10169 (1983). [31] R. D. Sharma, V . Kharchenko, Y . Sun, and A . Dalgarno, J . Geophys. Res. 101, 275-281 (1996). [32] V . Kharchenko, J. Tharamel and A . Dalgarno, J . atmos. terr. Phys. 59, 107-115 (1997). [33] P. Ong, M . Hogan, K . Lam and L. A . Viehland, Phys. Rev. A 45, 3997 (1992). [34] L. A . Viehland, and D. Hampt, J. Chem. Phys. 97, 4964 (1993). [35] S. L . L i n , L. A . Viehland and E. A . Mason, Chem. Phys. 37, 411-424 (1979). [36] T. Y . Wu, Kinetic 1966). Equations of Gases and Plasma (Addison-Wesley, Reading ma, [37] Grad H . , Commun. pure appl. Math 2, 331-407 (1949). [38] Grad. H . , Rarefied 1960). Gas Dynamics ed F. Devienne 100-138 (London, Pengamon, Bibliography 150 [39] K . Kumar and R. E. Robson, Aust. J . Phys. 26 157 (1973). [40] A . Mozumder, J. Chem. Phys. 72, 1657, 6289 (1980); ibid 74, 6911 (1981); ibid 76, 3277 (1982). [41] K . Koura, J . Chem. Phys. 87, 6481 (1987). [42] G. L. Braglia, Riv. Nuovo Cimento 18, 1 (1995). [43] B . Shizgal and D. R. A . McMahon, J . Chem. Phys. 84, 4854 (1984). [44] D. R. A . McMahon and B. Shizgal, Phys. Rev. A 31, 1894 (1985). [45] B . Shizgal and D. R. A . McMahon, Phys. Rev. A 32, 3669 (1985). [46] J. Duderstadt and W.R. Martin, Transport [47] H . Risken, The Fokker - Planck Equation, [48] M . Mitchner and C. H . Kruger, Partially Theory (Wiley, New York, 1979). 2nd ed (Springer, Berlin, 1979). Ionized Gases (Wiley, New York, 1973). [49] S. L. L i n , R. E. Robson and E. A . Mason, J. Chem. Phys. 71, 3483 (1979). [50] B . Shizgal and R. Blackmore, J. Computat. Phys. 55, 313-327 (1984). [51] W . Gautschi, J . Comput. Appl. Math. 12, 61-76 (1985). [52] W . Gautschi, Trans. Math. Software 20, 21-62 (1994). [53] P. L. Bahatnager, E. P. Gross, and M . Krook, Phys. Rev. 94, 511-525 (1954). [54] J . M . Warman and M . C. Sauer Jr., J . Chem. Phys. 62, 1971 (1975). [55] H . Shimamori and T. Sunagawa, Chem. Phys. Letters 267, 334 (1997). [56] K . Stamnes, Planet. Space Sci., 28, 427 (1980). [57] H . R. Skullerud, J Phys B : Atom Mol. Phys. 17, 913-929 (1984). [58] K . Kumar, H . R. Skullerud and R. E. Robson, Aust. J. Phys 33, 343 (1980). [59] J . H . Whealton and B . S. Woo, Phys. Rev. A 6 , 2319-2325 (1971). [60] D. Hubert, J. Geophys. Res. 87, 8255-8262 (1982). [61] D. Hubert, J . Phys D: Appl. Phys. 18, 1521-1532 (1985). [62] T. Makabe, K . Misawa, and T. Mori, J. Phys. D : Appl. Phys. 14, 199-206 (1981). Bibliography 151 L. Ferrari, Physica 93A 531-552 (1978). P. Nevai, in Orthogonal Polynomials: Theory and Practice (Kluwer Academic, Norwell, M A , 1990); M . J . Davis and P. Rabinowitz, Methods of Numerical Integration (Academic Press, New York, 1975); G. Freud, Orthogonal Polynomials (Pergamon Press, New York, 1971); T.S. Chihara, A n Introduction to Orthogonal Polynomials (Gordon and Breach, New York, 1978). B. Shizgal, J. Computat. Phys. 41, 309-329 (1981). R. Blackmore and B . Shizgal, Phys. Rev. A 31, 1855 (1985). H. S. Wilf, Mathematical Methods of Digital Computation (Wiley, New York, 1963). A. S. Clarke and B . Shizgal, J . Comput. Phys. 104, 140 (1993). J. D. Jackson, Mathematics York, 1975). for Quantum Mechanic, 93-94 (Academic Press, New B. Danloy, Maths. Comput. 27, 861-879 (1973). R. Blackmore, U . Weinent and B . Shizgal, Transport 15, 181 (vol.1), 120 (vol.2) (1986). Theory and Statistical Physics B. Shizgal and H . Chen, J. Chem. Phys. 104, 4137 (1996). B. Fornberg, A Practical Guide to Pseudospectral Press, Cambridge, 1996). Methods (Cambridge University C. Canuto, M . Y . Hussaini, A . Quarteroni and T . A . Zang, Spectral Method in Dynamics (Springer, 1987). Fluid I. Krajcar-Bronic, M . Kimura and M . Inokuti, J. Chem. Phys. 102, 6552 (1996). M . Dillon and M . Kimura, Phys. Rev. A 52, 1178 (1995). I. Krajcar Bronic and M . Kimura, J. Chem. Phys. 103, 7104 (1995). G. L. Braglia and P. Minari, Contr. Plasma Phys. 34, 711 (1994). K . Kondo and H . Tagashira, J. Phys. D: Appl. Phys. 26, 1948 (1993). R. Winkler, J. Wilhelm, G. L. Braglia and M . Diligenti, Nuovo Cimento, D, 12, 975 (1990). [81] G. L. Braglia, M . Diligenti, J. Wilhelm and R. Winkler, Nuovo Cimento D, 12, 257 (1990). Bibliography 152 [82] P. J . Drallos and J. M . Wadehra, Phys. Rev. A 40, 1967 (1989). [83] M . Yumoto, Y . Kuroda and T. Sakai, J. Phys. D: Appl. Phys. 24, 1594 (1991). [84] M . Yousfi, A . Hennad and A . Alkan, Phys. Rev. E 49, 3264 (1994). [85] K . Kowari, Phys. Rev. A 53, 853 (1996). [86] K . Kowari, L. Demeio and B . D. Shizgal, J. Chem. Phys. 97, 2061 (1992). [87] S. Yachi, H . Date. K . Kitamori and H . Tagashira, J . Phys. D: Appl. Phys. 24, 573 (1991) . [88] H . Date, S. Yachi, K . Kondo and H . Tagashira, J. Phys. D: Appl. Phys. 25, 42 (1992) . [89] S. Ushiroda, S. Kajita and Y . Kondo, J . Phys. D: Appl. Phys. 23, 47 (1990). [90] N . Goto and T. Makabe, J . Phys. D: Appl. Phys. 23, 686 (1990). [91] K . F . Ness and R. Robson, Phys. Rev. A 34, 2185 (1986). [92] R. Robson and K . F. Ness, Phys. Rev. A 33, 2068 (1986). [93] B . Shizgal, D. A . R. McMahon and L. A . Viehland, Radiat. Phys. and Chem. 34, 35 (1989). [94] G L . Braglia and L. Ferrari, Nuovo Cimento B 4, 245 (1971). [95] J . M . Warman, U . Sowada and M . P. deHaas, Phys. Rev. A 31, 1974 (1985). [96] B . Shizgal, Chem. Phys. 147, 271 (1990). [97] R. Nagpal and A . Carscadden, Phys. Rev. Letters 73, 1598 (1994). [98] B . L . Tembe and A . Mozumder, J . Chem. Phys. 78, 2030 (1983); ibid 81, 2492 (1984); Phys Rev. 27, 3274 (1983). [99] K . D. Knierim, M . Waldman and E. A . Mason, J. Chem. Phys. 77, 943 (1982). [100] K . Koura, J . Chem. Phys. 84, 6227 (1986); ibid 82, 2566 (1985); ibid 81, 4180 (1984). [101] K . Leung, J . C. Barrett and B . Shizgal, Chem. Phys. 147, 257 (1990). [102] L. A . Viehland, S. Ranganathan and B. Shizgal, J . Chem. Phys. 88, 362 (1988). Bibliography 153 [103] D. R. A . McMahon, K . Ness and B. Shizgal, J. Phys. B: At. Mol. Phys. .19, 2759 (1986). [104] R. Blackmore and B . Shizgal, J. Computat. Phys. 55, 313 (1984). [105] G. Mansell, W . Merryfield, B . Shizgal and U . Weinert, Comput. Meth. Appl. Mech. Engrg. 104 295 (1993); H . H . Yang, B . R. Seymour and B . D. Shizgal, Comp. Fluids 23, 829 (1994); H . H . Yang and B. Shizgal, Comput. Meth. Appl. Mech. Engrg. 118 47 (1994). [106] D. Gottlieb and S. Orzag, Numerical Analysis of Spectral Methods ( S I A M - C B M S , Philadelphia, 1977). [107] J . S. Chang and G. Cooper, J. Computat. Phys. 6, 1 (1970). [108] E . W . Larsen, C. D. Livermore, G. C. Pomraning and J . G. Sanderson, J. Computat. Phys. 61, 359 (1985). [109] E. M . Epperlein, Laser and Particle Beams 12, 257 (1994); J. Computat. Phys. 112, 291 (1994). [110] R. E. Robson, K . F. Ness, G. E. Sneddon and L. A . Viehland, J. Computat. Phys. 92, 213 (1991). [Ill] R. E. Robson and A . Pritz, Aust. J. Phys. 46, 465 (1993). [112] B . T. Park and V . Petrosian, Astrophys. J. Suppl. Series 103, 255 (1996); Astrophys. J . 446, 699 (1995). [113] A . Gilardini, Low Energy Electron Collisions in Gases (Wiley, New York, 1974); L. G. H . Huxley and R. W . Crompton The Diffusion and Drift of Electrons in Gases (Wiley, New York, 1974). [114] R. K . Nesbet, Phys. Rev. A 20, 58 (1979). [115] S. R. Hunter, J.Carter and L. G. Christophorou, J. Chem. Phys. 80, 6150 (1984). [116] A . Pagamenta, M . Kimura, M . Inokuti and K . Kowari, J. Chem. Phys. 89, 6220 (1988). [117] L . G. Christophorou, Nuc. Instrum. Methods Phys. Res. A268, 424 (1988); L. G. Christophorou and R. J . Van Brunt, I E E E Trans. Dielectrics Elect. Insul. 2, 952 (1995).. Bibliography 154 [118] H . Itoh Y . Miura, N . Ikuta, Y . Nakao and H . Tagashira. J. Phys. D: Appl. Phys. 21, 922 (1988); K . Satoh, H . Itoh, Y . Nakao and H . Tagashira. J. Phys. D: Appl. Phys. 21, 931 (1988); [119] L. G. Christophorou and S. R. Hunter, in Electron-Molecule Interactions and their Applications, vol 2, edited, by L. G. Christophorou (Acedemic Press, New York, 1984). - . ' . [120] L. G. Christophorou and M . 0 . Pace, Gaseous Dielectrics IV, (Pergamon Press, New York, 1984). [121] T. Nishigori and B . Shizgal, J . Chem. Phys. 89, 3275 (1988); B . Shizgal and T. Nishigori, Chem. Phys. Letters 171, 493 (1990). [122] K . Kowari and B . Shizgal, Chem. Phys. 185 1, (1984). [123] B . Shizgal, Chem. Phys. Letters 138, 65 (1987). [124] B . Shizgal, J . Phys. B: At. Mol. Opt. Phys. 21, 1699 (1988). [125] B . Shizgal, J . Phys. B: At. Mol. Opt. Phys. 24, 2909 (1991). [126] K . Kowari and B . Shizgal, Chem. Phys. Letters 260, 365 (1996). [127] G. N . Haddad and T. F. O'Malley, Aust. J. Phys., 35, 35 (1982). [128] H . B . Milloy, R. W . Crompton, J. A . Rees, and A . G. Robertson, Aust. J . Phys., 30, 61 (1977). [129] Isao Shimamura, Scientific Papers of the Institute of Physical and Chemical Research, 82, 1 (1989). [130] T. F. O'Malley and R. W . Crompton, J. Phys. B: Atom. Molec. P h y s B l 3 , 3451 (1980). [131] M . Hayashi, Swarm Studies and Inelastic Electron Molecule collisions (Springer, New York, 1987), 167. private communication [132] I. Shimamura, Phys. Rev. A 46, 1394 (1992). [133] A . Chutjian and S. H . Alajajian, Phys. Rev., A 3 1 , 2885 (1985). [134] T. Shimanouchi, Tables of Molecular Vibrational Frequencies, Consolidated Volume I, Nat. Stand. Ref. Data. Ser., Nat. Bur. Stand (US), 39, (1972). [135] R. P. Blaustein and L. G. Christophorou, J. Chem. Phys. 49, 1526 (1968). Bibliography 155 [136 W. H . Press, B . P. Flannery, S. A . Teukolsky, and W . T. Vetterling, Recipes 2nd edition (Cambridge, 1989). Numerical [137; M . Garcia-Gracia, J . Letosa, M . Artacho and J. M . Fornies-Marquina, I E E E Trans, on Magnetics 33, 1464 (1997); M . R. Smith and S. T. Nichols, Nucl. Instrum. Method 205, 479 (1983); M . R. Smith, S. Cohn-Sfetcu and H . A . Buckmaster, Tecnometrics 18, 467 (1976). [138 R. E. Robson, J. Chem. Phys. 85, 4486 (1986). [139 J. Keizer, J. Chem. Phys. 58, 4524 (1973). [140 B. Shizgal, J. Chem. Phys. 69, 5355 (1978). [141 M . H . Rees, Physics and Chemistry University Press, Cambridge, 1989). of the Upper Atmosphere, ch. 3, (Cambridge [142 Z. Alfassi and S^Amiel, J. Chem. Phys. 57, 5085-5087 (1972). [143 K . Nanbu and Y . Watanabe, Rep. Inst. High Speed Mech. 34, 21-54 (1981). [144; J. R. Stallcop and H . Partridge, Phys. Rev. 32, 639 (1985). [145 R. G. Burnside, C. A . Tepley and V . B. Wickwar, Annal. Geophys. 5A, 343 (1987). [146 A . S. Clarke and B . Shizgal, Phys. Rev. E49, 347 (1994). [147 B . Shizgal and D. Hubert, Rarefied Gas Dynamics 16, 3 (1989). [148 K . Waldeer and H . Urbassek, Physica A 176, 325-344 (1991). [149 J. E. Salah, Geophys. Res. Lett. 20, 1543 ; J. R. Stallcop et. al., J.Chem. Phys. 95, 6429 (1991). [150 Wang Chang C. S. and G. E. Uhlenbeck "The Kinetic Theory of Gases", Studies in Statistical Mechanics V.5 G. E. Uhlenbeck and De Boer Eds Elsiner, New York, 155-190 (1970). [151 B. Shizgal, Can. J. Phys. 42, 97 (1984). [152 B. Shizgal and R. Blackmore, Chem. Phys. 77, 417 (1983). [153 B . Shizgal, M . J. Lindenfeld and R. Reeves, Chem. Phys. 56, 249 (1981). [154 B . Shizgal and M . J. Lindenfeld, Planet. Space Sci. 27, 1321 (1979). [155 L. A . Viehland, Chem. Phys. 179, 71-92 (1994). Bibliography 156 [156] R. Kapral arid J . Ross, J. Chem. Phys. 52, 3, 1238 (1970). [157] L. Monchick and E. A . Mason, Phys. Fluid 10, 1377 (1967). [158] A . D. Koutselos, E. A . Mason and L. A . Viehland, J. Chem. Phys. 93, 7125 (1990). [159] R. B . Bernstein, J . Chem. Phys. 33, 795-805 (1960). [160] P. Drallos and M . E. Riley, Phys. Rev. A 44, 1405-1408 (1991).
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Electron and ion relaxation in atomic and molecular...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Electron and ion relaxation in atomic and molecular moderators Leung, Ki Yee 1997
pdf
Page Metadata
Item Metadata
Title | Electron and ion relaxation in atomic and molecular moderators |
Creator |
Leung, Ki Yee |
Date Issued | 1997 |
Description | The standard technique for the solution of kinetic theory and transport problems usually involves the expansion of the velocity distribution function (VDF) in some suitable set of basis functions. For steady problems, the Boltzmann equation (BE) for the distri-bution function is then reduced to algebraic equations for the expansion coefficients. For time dependent problems, the B E is reduced to ordinary differential equations. The main objective with this solution technique is to calculate the transport coefficients which can be expressed in terms of the lower order expansion coefficients. Recent experiments on different physical systems are capable of probing the details of the velocity distribution functions. The main objective of this thesis is to develop accurate and efficient numerical methods for the BE with the aim to calculate the distribution function. For an ensemble of ions dilutely dispersed in a background of neutrals under the influence of an applied electric field, the exact solution of the Bhatnagar, Gross, and Krook (BGK) collision model is used to study different numerical solutions of the BE. In this study, the main theme of the thesis is presented. In place of the traditional expansion techniques, a discretization of the VDF is proposed. The discretization is based upon a set of polynomial basis functions that yield the optimum convergence of the solution. For electrons dilutely dispersed in a background of neutrals, the BE is approximated by a Fokker-Planck equation (FPE). The one-dimensional F P E has been used to describe many different phenomena in chemistry, physics and engineering. A detailed comparison of the eigenvalues and eigenfunctions of the FP operator is carried out with different numerical methods including a finite difference (FD) scheme, Lagrange interpolation (LI) method, the Quadrature Discretization Method (QDM) with speed points and the QDM with Davydov points. The QDM involves the discretization of the distribution at a set of points that coincide with the points in a quadrature that is based upon some weight function. The Davydov distribution is the steady distribution for electrons in a gas moderator at some given electric field strength. The time dependent solutions of the FP are obtained with the different energy discretizations applied in conjunction with an appropriate time stepping procedure. The QDM with Davydov points provides a more rapid convergence of the eigenvalues and eigenfunctions relative to the other methods used. The relaxation of a nonequilibrium distribution of electrons in a mixture of CCI4 with either Ar or Ne is studied. The electron collisions in this analysis include elastic collisions between electrons and CCI4 and the inert gas, vibrationally inelastic collisions between electrons and CCI4, as well as the electron attachment reaction with CCI4. The time dependent electron energy distribution function is determined from the BE and the energy relaxation times are determined. The coupling of the thermalization process and the attachment process is discussed in detail. The results from the calculations are compared with experimental studies, and the methodology of the experimental reduction of the data is examined in detail. In addition, the nature of the spectrum of the Boltzmann collision operator for realistic ion-atom interactions is studied with a discretization of the integral collision operator. The singular nature of the collision operator is considered in detail with a novel interpolation technique. The eigenvalues and eigenfunctions for the Henyey-Greenstein (HG) model differential cross section are calculated. The relaxation of anisotropic ion distri-butions with and without an applied electric field is studied. |
Extent | 6856126 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-04-01 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0060250 |
URI | http://hdl.handle.net/2429/6710 |
Degree |
Doctor of Philosophy - PhD |
Program |
Chemistry |
Affiliation |
Science, Faculty of Chemistry, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 1997-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_1997-25092X.pdf [ 6.54MB ]
- Metadata
- JSON: 831-1.0060250.json
- JSON-LD: 831-1.0060250-ld.json
- RDF/XML (Pretty): 831-1.0060250-rdf.xml
- RDF/JSON: 831-1.0060250-rdf.json
- Turtle: 831-1.0060250-turtle.txt
- N-Triples: 831-1.0060250-rdf-ntriples.txt
- Original Record: 831-1.0060250-source.json
- Full Text
- 831-1.0060250-fulltext.txt
- Citation
- 831-1.0060250.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0060250/manifest