UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Collisional relaxation in plasmas Crossman Statter, Gregory Christopher 1988

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-UBC_1989_A6_7 S72.pdf [ 5.07MB ]
Metadata
JSON: 831-1.0059427.json
JSON-LD: 831-1.0059427-ld.json
RDF/XML (Pretty): 831-1.0059427-rdf.xml
RDF/JSON: 831-1.0059427-rdf.json
Turtle: 831-1.0059427-turtle.txt
N-Triples: 831-1.0059427-rdf-ntriples.txt
Original Record: 831-1.0059427-source.json
Full Text
831-1.0059427-fulltext.txt
Citation
831-1.0059427.ris

Full Text

COLLISIONAL RELAXATION IN PLASMAS By GREGORY CHRISTOPHER CROSSMAN STATTER B.Sc, Sydney University, 1979 M.Sc, Sydney University, 1982 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS OF MASTER OF SCIENCE in THE FACULTY OF GRADUATE STUDIES Department of Chemistry We accept this thesis as conforming to the required standard UNIVERSITY OF BRITISH COLUMBIA October, 1988 © Gregory Christopher Statter, 1988. In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of The University of British Columbia Vancouver, Canada DE-6 (2/88) ii Abstract The energy relaxation of two types of fully ionized plasma systems are determined with the solution of the Fokker-Planck equation. In both cases the plasma constituents are treated as being point-like and structureless and the plasma relaxes collisionally in the absence of spatial gradients and external electric and magnetic fields. The first plasma system consists of one plasma species dilutely dispersed in a second plasma which acts as a heat bath at equilibrium. The initial energy distribution of the dilute constituent is chosen to be a delta function and the approach to a Maxwellian distribution at the heat bath temperature is determined. The second plasma system consists of just one plasma species which initially possesses a bi-Maxwellian ion velocity distribution function (VDF). The average of the kinetic energy, m<v,f >/2 for particle motions in some arbitrary z direction serves to define a temperature parameter T„ =m<v,?>/k for this degree of freedom. Similarly, the temperature Tx=m<v^>/2k, where v ±is the velocity component in the plane perpendicular to z, parametrizes the average energy for the other two translational degrees of freedom The relaxation of T„ and T± to a common temperature, T, via self-collisions of the plasma is studied. In both cases the collisional relaxation can be described by a linear collision operator and the expansion of the distribution function in the eigenfunctions of the Fokker-Planck operator is considered. The reciprocals of the corresponding eigenvalues are the characteristic relaxation times for the system. For the first plasma system, the temperature relaxation time determined with the solution of the Fokker-Planck equation is compared with the relaxation time calculated with the assumption that the distribution function remains Maxwellian all the time. For the second plasma system the relaxation time is the characteristic time for the relaxation of one of the temperature components and is compared iii with the relaxation time calculated with the assumption that the distribution function remains bi-Maxwellian all the time. The results are compared with other theoretical and experimental results. A study of the eigenvalue spectrum of the Fokker-Planck operator and the temporal approach to equilibrium is also emphasized. iv Table of Contents Abstract ii Table of Contents iii List of Tables v Captions to Figures vi Acknowledgements ix I Chapter 1 1 1.1 Introduction 1 1.2 Formal Solution 3 1.3 Numerical Methods 5 1.4 The Equations to be Solved 8 II Relaxation of a Charged Species in an Equilibrium Fully Ionized Plasma. 11 II. 1 Introduction 11 11.2 Time Evolution of the Distribution Function 12 11.3 Numerical Method 16 n.4 The Equivalent Schroedinger Eigenvalue Problem 17 11.5 Perturbation Expansion of L about y=0 19 11.6 Results and Discussion 21 11.7 Conclusions 49 HI Anisotropic Temperature Relaxation due to Coulomb Collisions in a Spatially Homogeneous One Component Plasma .50 III. 1 Introduction 50 III.2 Theoretical Approach 53 (a) Formulation of the Problem and the Linearization of the Fokker-Planck equation 53 (b) Relaxation of T,,(t) 58 (c) The Time-Decay of the Expansion Coefficients 60 V (d) The Initial Expansion Coefficients, un/(0) 64 111.3 Results and Discussion 66 111.4 Conclusions .' 90 IV Summary 91 Appendix 1 93 Appendix 2 104 REFERENCES 106 vi List of Tables Chapter I Chapter n Table n. 1 Convergence of Eigenvalues 27 Table n.2 Variation of the Eigenvalues with the scale factor, s 31 Table n.3 Comparison of Relaxation Times 47 & 48 Chapter m Table UL 1 Convergence of Eigenvalues 69 Table RI.2 Convergence of Relaxation times, tj, as a function of cc(0) and basis set size, N 71 Table in.3 Comparison of Relaxation Times 89 Chapter IV Appendix 1 Table Al Convergence of M Matrix Elements using Speed, Associated Laguerre (/=0) and Laguerre Quadratures 98 Appendix 2 REFERENCES vii CAPTIONS TO FIGURES II Figure n. 1: Potential functions in the Schroedinger equation equivalent to the Fokker-Planck eigenvalue equation 24 & 25 Figure II.2: Eigerrfunctions of the Fokker-Planck operator, curve a, and the first order perturbation approximation to the eigenfunctions, curve b (see equation (II.5.7)). T b = 300K; y = 0.04 (A) X3 = 10.102, (B) X4 = 12.659, (C) Xs = 14.775 and (D) X6 = 16.394. These eigenvalues are the quasi-bound, "discrete" states shown in the potential of Figure 1 and converge rapidly as shown in Table 111 for some of the eigenvalues 28 & 29 Figure TL3: Temperature relaxation (A) and the effective relaxation rate (B): T b = 300K, T(0)/Tb = 24 and (a) y = 0.04, (b) y = 0.09, (c) y = 0.16 and (d) y = 0.25 32 Figure LI.4: Variation of eigenvalues of the Fokker-Planck operator with mass ratio. The points are eigenvalues obtained from the Fokker-Planck operator, whilst the lines correspond to the first order perturbation approximation to the eigenvalues; see equation (U.5.6) of the text 34 Figure JJ.5: Ratio of the absolute values of the sum of the contributions from the continuum spectrum, Sc, to the discrete spectrum, Sd. T b = 300K, T(0)/Tb = 24 and (a) y = 0.04, (b) y = 0.0625, (c) y = 0.09, (d) y = 0.1225 and (e) y = 0.16; see text for definitions and discussion 35 Figure II.6: Variation of the energy expansion coefficients, c(A), with eigenvalue, A,, for the continuum portion of the spectrum Tb= 300K, T(0)/Tb = 24, c(K) in units of 103. The triangles are the numerical results for four different values of the scaling parameter s 38 Figure n.7: Variation of the energy expansion coefficients, c(X), with eigenvalue, X, for the continuum portion of the spectrum. T b = 300K; the triangles are the numerical results for four different values of the scaling parameter, s. T(0)/Tb is equal to (A) 10.667 and (B) 16.667, c(K) in units of 102 39 Figure n.8: Turning points in the potential function of the Schroedinger equation equivalent to the Fokker-Planck eigenvalue equation 40 viii Figure II.9: Figure n. 10: Figure 11.11: Figure UJ.1: Figure HI.2: Figure LLI.3: Figure m.4: Figure HI.5: Comparison of the continuum eigenfunctions determined from the diagonalization of the matrix representative of the Fokker-Planck operator and the JWKB approximate functions. T b = 300K, s = 0.6042, y = 0.16; (A) X = 4.855 and (B) X = 7.603 43 Effective relaxation rate (A) and temperature relaxation (B) Tb= 300K, T(0)/Tb = 24 and (a) y = 0.25, (b) y = 1 and (c) y = 25. x is the time in units of to for the temperature difference T(t')-Tb to decrease by 1/e of the initial value 44 Logarithmic plot of x versus y for the results in Table II.3. Lines a-c are the Fokker-Planck relaxation times, x, for T(0)/Tb = 2.667, 10.667 and 24 respectively, whilst lines e-f are the Spitzer relaxation times, Xj for these T(0)/Tb values 49 CAPTIONS TO FIGURES HI Eigenfunction plots. Number of basis functions used in A is 40, B 80. Curve 'a' consists of the lowest 15 eigenfunctions for A and the lowest 35 eigenfunctions for B. Curves b and c, in both figures, consist of eigenfunctions with larger eigenvalues than those used in curve a 70 Variation of the relaxation time, tj, with basis set size, N. Note that both figures show that tj possess a region of convergence. Most relaxation times are satisfactorily estimated with N=10. Figure B leads to the estimate that t x =1.45±0.02 73 Variation of the initial expansion coefficients, ur2(0) with r. Curves a-f correspond to a(0) = 2,3,4,5,6 and 10 respectively. Note the growth in the coefficients for a ( 0 ) £ 5 74 Temperature coefficients, b(2)[X(2)] as a function of X.(2). A contains plots of b(2)[Xi2h for a(0) = 0, 0.2, 0.4....2.0. B is for a(0) = 6.0. Both plots used 80 basis functions .....76 Relaxation of the anisotropic temperature parameters T„(t')/T (A and B), cc(t') (C and E) and e(t') (D and F), for initial temperature anisotropies as labelled 77-79 ix Figure LTJ.6: Theoretical and experimental results for relaxation of Tu(t') and T±(t'). The experimental points are the HDM data from their figure 2. The number densities, n, for the theoretical curves are as labelled 80 Figure LTJ.7: Theoretical and experimental rates of relaxation for various electron plasma densities and temperatures. The experimental points are the HDM points from their figure 3. The triangles refer to a(0)<l in their experiment and the circles to cc(0)>l. The upper line is the rate determined from equation (TJJ.3.5) for an electron plasma, the lower line is the result from the present study for a (0) = 8 81 Figure LTJ..8: The function h(T±/T) = (T/T e f f) 3 / 2 given in equation (LTI.3.7) plotted versus a. A is for a < 1, while B is for a > 1, corresponding to the upper and lower equations in equation (IH.3.7) 86 Figure LTJ.9: Ratio of the expansion coefficients, wr 2(t')/ur 2(t') for the initial bi-Maxwellian ion VDF with cc(0) = 2. The former coefficient, wr 2(t'), is the coefficient for the ion VDF always being bi-Maxwellian, and ur 2(t') is calculated from the Fokker Planck equation, and is obtained from 'equations (III.2.24 and 32) 91 X Acknowledgement I would like to thank Dr. B. Shizgal for his suggestion of these topics and his assistance and advice with many problems. I would like to thank Dr. M. Comisarow for his interest in these topics in connection with Ion Cyclotron Resonance Broadening and Drs. Kevin Ness and Jonathan Barrett whom were most helpful and patient in explaining many aspects related to this work. I am grateful to the Chemistry Department and to the university for providing me with a University Graduate Fellowship (U.G.F). 1 Chapter I 1.1 Introduction This thesis deals with the study of the collisional relaxation to equilibrium of non-equilibrium distributions of charged particles, or plasmas. The forces of interaction between these charged particles are Coulombic and the rate of relaxation of the plasma species, determined from the Fokker-Planck collision operator, is of fundamental concern to plasma collisional processes in general with numerous important practical applications. There have been many theoretical treatments of this subject^1-12) as well as numerous experimental studies^3-18). The rate at which charged particles lose energy via Coulomb collisions with the electrons and ions of the plasma is an important consideration with regard to a very broad branch of physics and chemistry. For example in the development of controlled thermonuclear fusionC19-26) the heating of a fusion plasma can be acheived by the injection of energetic charged particles. Coulomb collisions are also an important consideration in the solar wind(27-34), whose constituents are typically protons, helium nuclei and electrons possessing different temperatures. Furthermore, the velocity distribution of the constituents in the solar wind are frequently modelled as being anisotropic(27-34), that is, the average kinetic energy for motion in some direction differs form the average energy for motion in the perpendicular direction. These considerations are also very important in the ionosphere, aeronomy and the polar wind(33-36), where, for example, oxygen ions (0+) can be up to 10 times hotter along one direction than in the perpendicular direction^6). Also the operating conditions of many instruments such as the inductively coupled plasma(37-39) and the (Fourier Transform) ion cyclotron resonance mass spectrometerC40-43) may depend on the relaxation behaviour of such charged particle systems. For example charge particle interactions can lead to Coulomb broadening in ion 2 cyclotron resonance mass spectrometry^0-43). A final example of the importance of Coulomb collisions is in the evolution of star clustersC44-46). Since systems of stars are governed by the 1/r gravitational potential, the problem of the relaxation of star clusters and galaxies is analogous to the problem of the relaxation of a plasma species which is governed by the 1/r Coulomb potential. The study of the relaxation (or thermalization) of charged particles in ionized gases is clearly an important endeavor. The two problems studied in this thesis are done so in the absence of spatial gradients and electric and magnetic fields. They are fundamental and important to the study of collisional relaxation operative in many plasmas. Specifically these two problems are concerned with: (i) The relaxation of an isotropic system where charged 'test' particles of mass m are introduced into an equilibrium fully ionized plasma (or moderator) of mass M. The test particles are chosen to possess an initial isotropic delta function distribution of energies and the rate of equilibration of test particle energy (or temperature) with the moderator species is studied both as a function of moderator to test particle mass ratio, y=M/m, and as a function of the initial test particle energy. Since the moderator is at equilibrium and present in large excess the problem is linear. The second problem studied in chapter UL is: (ii) The relaxation of an initial bi-Maxwellian (anisotropic) single component plasma. This problem, unlike the former, is non-linear because the plasma test particles undergo collisional relaxation with themselves only (there is no moderator species present). For this problem the average kinetic energy, m<v2>/2 for particle motions along some 3 arbitrary axis, z, serves to define a temperature parameter TB =m<v,f>/k for this degree of freedom. Similarly the temperature T±=m< v 2 > 12k where v ± is the velocity component in the plane perpendicular to z, parametrizes the average energy for the other two translational degrees of freedom. The relaxation of T„ and T x to a common temperature, T, via self-collisions of the plasma is studied. 1.2 Formal Solution Both physical situations i.e (i) and (ii) above can be described by a Fokker-Planck equation of the form ^ = Lf(M) , (1.2.1) where Lis a linear operator and f(x,t) is the time dependent distribution function, where x is a dimensionless or reduced velocity, defined by x = 7m/2kT' v, where T' is a temperature parameter that may differ from the equilibrium temperature, T. In problem (ii) the operator L results from the linearization of a non-linear Fokker-Planck operator, by assuming that the ion velocity distribution function (VDF) is only slightly perturbed from its equilibrium value The specific forms for these operators are dealt with in chapters II and HI respectively and the details of the linearization involved with problem (ii) are presented in chapter HT of this thesis. The formal solution of equation (L2.1) is f (x,t) = e-Ltf(x,0) (1.2.2) 4 where f(x,0) is the initial VDF. If f(x,t) is expanded in the eigenfunctions, {<|>k}JL0, of L then we may write f(x,t) = X c k e A t V x ) • (L2'3> k where 1^= W> ^ being an eigenvalue of L, and c± an expansion coefficient of f(x,0) in the eigenfunction basis. Both and ^ are calculated numerically by diagonalizing the matrix representation of L in a suitable basis as discussed in section 1.3 on numerical methods. The E implies that the spectrum of L is discrete although this is not always the k case(47-51), and in general f(x,t) is expressed as a sum over the discrete eigenvalues and an integral over the continuum portion of the spectrum The nature of the eigenvalue spectrum of the collision operator determines the details of the time evolution of the distribution function and there have been many general discussions of this in the literature, particularly with regard to neutron transport theory^52) hot atom chemistry^53) and other kinetic theory problems( 4 7« 4 8). Much of the interest is with respect to the presence or absence of a continuous portion of the eigenvalue spectrum. It has generally been thought that the eigenfunction method of solution will not be valid or will converge slowly if the continuous portion of the spectrum is not included. However, it has recently been demonstrated that expansions with finite basis sets which do not include the continuum eigenfunctions in a rigorous fashion can yield converged solutions(49)- These considerations are also of considerable interest in connection with the present problems where the continuum eigenfunctions and eigenvalues often dominate the approach to equilibrium. One of the objectives of the present work is the study of this eigenvalue spectrum and the relationship 5 with the temporal approach to equilibrium. 1.3 Numerical Methods Two methods for detennining the eigenvalues and eigenfunctions of L are implemented in this thesis. The first method of solution, used in chapter II, is based on the Discrete Ordinate (DO) method by Shizgal and Blackmore(54'55), where the eigenfunctions are determined at a set of points which coincide with the points of a quadrature procedure based on the weight function w(x) = x e x . This particular choice of weight function, as discussed by ShizgaK56), is particularly well suited to many kinetic theory problems since it is based on the particle speed rather than the particle energy. The method consists of representing functions, in particular the eigenfunctions and distribution functions, as column vectors, f, whose N elements are f. =ywif(x.) . The points {x;} are the quadrature points, and {wj} are the quadrature weights corresponding to the integration rule, ~ N J dx w(x) f(x) « X w i f (xi> ' ( L 3 ' X ) o i = i where f(x) is typically some function of the reduced speed, x, defined on the interval (0,oo). Formally, the function f(x) is represented in terms of its N expansion coefficients, P, in the speed polynomial basis, {Bn(x)}~=Q. For a more detailed discussion of these polynomials, which are orthonormal with respect to the weight function w(x), the reader is referred to ShizgaK56). The transformation to the DO representation from the polynomial representation is then achieved through the unitary transformation matrix T whose elements 6 arc given by T n j = BJx^Jwy The crux of the DO method as developed by Shizgal and Blackmore(54»55) is the representation of the derivative operator d/dx in the discrete space defined by the points in the quadrature rule, equation (1.3.1). This procedure enables the linear Fokker-Planck operator, L, to be discretized in this same discrete space. Since this discretized operator L can be written in an explicit form (see for example chapter n equation (II.3.1) and Blackmore and ShizgaK54)) it is not necessary to consider the explicit expansion of the Fokker-Planck operator in the speed polynomial basisC54-56). The second method of solution, used in chapter in, is based on an explicit polynomial expansion of the linearized Fokker-Planck operator for self-collisions. A transformation to a DO basis set is not employed, although there would be little difficulty accomplishing this. Instead the eigenfunctions are determined in a polynomial basis and not at quadrature points. The Burnett functions(57_6°) which consist of a spherical harmonic, Y, m(6,<|>), for the angular dependence of the VDF multiplied by a generalized Laguerre (or Sonine(57-60)) polynomial, L / n + 1 / 2(x 2), for the energy variable are used as basis functions. The Burnett functions are the eigenfunctions of the Boltzmann collision operator for Maxwell molecules^9-62), (characterized by a 1 /r 4 interaction potential) and are often used in kinetic theory problems(63). For the azhnuthally (or cylindrically) symmetric anisotropic VDF for the plasma considered in chapter HI the Burnett functions can be written as, Vf( X) = x'l4+ i yV)P|(«»e) , 0-3.2) where m = 0 in the spherical harmonics reduces them to the Legendre polynomials, P^ cosO), where 8 is the polar angle describing the angle of inclination of the velocity, v, to 7 the symmetry or z axis of the VDR In chapter LTJ the matrix elements of the Fokker-Planck operator in the Burnett function basis are reduced to integrals involving functions of the reduced speed, x, w(x) and combinations of L / n + 1 / 2(x 2) (the angular integrals being readily performed). The technique of generating matrix elements from the generating function for the polynomial basis is well known and widely used in kinetic theoryC5 9 - 6 1'6 3 - 6 7) and is the method used to obtain the matrix elements in chapter ILL It is preferable to evaluating the individual matrix elements since the reduced speed variable, x, appears in exponential form instead of in polynomial formC64'65). Additionally, the use of the generating function leads to recursion relations for the matrix elements which would otherwise not be as easily derived. For example, that an element of a set of matrix elements is given by, Kl " < C V ) IL I L' s + 1 / 2(x 2)> , (1.3.3) where the angular brackets denote / dx f(0)(x), where f*°) is the equilibrium Maxwellian o VDF and L operates on the basis function to its right The generating function for the matrix elements, M ( / \ is then obtained from the generating function for L?n + 1 / 2(x2), which is(59-£l,63-68) G f t . x ) = ( l - 5 ) - ' - 3 ' 2 H t p ( ^ ) = ±CCU2V) • (1.3.4) 1 ^ m=0 Hence, 8 M ( ° = Y, Mr.? $ V = < G($> X ) 1 £ 1 G ( t 1' X ) > (L 3-5) r,s = 0 is the generating function for the M matrix elements, and the matrix element, , is i*S obtained from the coefficient of C/TJs in the expression forM ( , ). The numerical codes for both problems (i) and (ii) are very efficient and the computation times are typically cbtermined by the generation of the matrix elements in their pertinent bases and their diagonalization to yield the eigenvalues and eigenfunctions. For example, on a u-Vax H machine the total computation time even for large basis set sizes (e.g 80 basis functions) never exceeds two minutes for either code. 1.4 The Equations to be Solved The notation used in this thesis closely follows that due to Hinton(5) and it is assumed that, for the problems to be solved, the effects of charged particle interactions are adequately described by Fokker-Planck equations of the form*5), IT = • ( I A L A ) where CJW = • [ A . f a " i £ • ( D a 0 ] ("•">) is the Fokker-Planck collision term, and fa and fb are the VDFs for species a and b respectively. In equation (1.4. lb) A a is the dynamical friction vector and D a is the velocity diffusion tensor. They are written; 9 A a = "2 % [l + f dV fb(V) JL , (i.4.2a) ma b j u and D a = -2 ^ J dv' fb(v') [U2LJUU] , CL4.2b) where u is the relative velocity defined by u = v - v' and c a b = 2jce2eJlnA . (1.4.3) In equation (1.4.3), ea=zae and e^ =z^ e. are the charges on species a and b respectively and the terms with index b represent the effect of scattering of particles of species a by particles of species b. For the problem considered in chapter n, species b is the bath (or moderator) species, whilst in chapter HI, where like species scattering is considered, species b is identical to a (i.e a=b). In equation (I.4.2b) dyadic notation has been used for the second rank tensors, with 1 being the unit dyadic, equivalent to the Kronecker delta in subscript notation. The Coulomb logarithm, InA, is defined by Spitzer^1'69-71) as lnA = ln(XD/b0) , (1.4.4) where Xn = A - ^ - and bn = (1.4.5) are, respectively, the Debye shielding distance^1'69-71) (for species b shielded by electrons) and the distance of closest approach for a typical thermal particle. In equation (1.4.5) nb is 10 the number density of species b. The dimensionless parameter A is frequently referred to as the plasma parameter^9-71) and is the average number of particles present in a Debye sphere (a sphere of radius XQ). For most plasmas^1'69-71) A»10 4 so that at any one time a particle of species a is interacting with a very large number of particles of species b. For the isotropic linear test particle problem considered in chapter U, fb is always Maxwellian enabling the integrals appearing in the expressions for the friction and diffusion tensors (equations (L4.2a and b)) to be evaluated analytically, so that the resulting Fokker-Planck equation is of differential form only. For problem (ii) these integrals cannot be evaluated analytically and the resulting Fokker-Planck equation is integro-differential in form. The methodology employed for solving this latter integro-differential Fokker-Planck equation starts with the Landau form of the Fokker-Planck operator^ 5*72) and employs the generating function technique based on associated Laguerre polynomials, previously discussed in section 3. It is important to point out that the use of the Fokker-Planck operator is only approximate in that it neglects collisions with small impact parameters and large scattering angles. This has been discussed recently by several workers^73-75) in connection with applications to fusion devices, by ShoubC76) for plasma systems and gravitational systems, and by Retterert45'46) with regard to the application of kinetic theory to globular clusters. For more detail on the specific equations solved in this thesis the reader is referred to chapters II and DX The theoretical approach and methodology described above are now implemented in chapters II and ILL In chapter II we solve the linear test particle problem using the DO method (i.e problem (i) defined on page 2), while in chapter HI we solve the linearized non-linear Fokker-Planck operator for a one component plasma using an expansion in Burnett functions (i.e problem (ii) defined on page 2). 11 Chapter IT Relaxation of a Charged Species in an Equilibrium Fully Ionized Plasma  II.1 Introduction This chapter is concerned with the relaxation of an isotropic distribution of a charged species, or charged "test particles" dilutely and uniformly dispersed in a second charged species which serves as a heat bath at equilibrium. The extension of the present work to the relaxation in a multicomponent equilibrium plasma is staightforward. We also do not include the effects of external electromagnetic fields in this study. The kinetic equation for this physical problem is linear and should be contrasted with the situation for which none of the system species can be considered to be at equilibrium and the resulting kinetic equations are non-linear^2-4-29*30) as discussed in chapter I. The approach adopted in the present chapter is based on the Fokker-Planck equation (FPE) for the test particle distribution function where an expansion of the distribution function in the eigenfunctions of the Fokker-Planck operator is considered (see chapter 1.2). The corresponding eigenvalues are the reciprocals of the characteristic relaxation times for the approach of the system to equilibrium. A detailed study of the eigenvalue spectrum of the Fokker-Planck operator is considered in the context of the transformation of the FPE to the equivalent Schroedinger equation. The methodology of the present work follows closely with the work on the transient behaviour of electrons in neutral moderators(77»78) and numerical methods that have been developed recently for the solution of a large class of Fokker-Planck equationst54-56). An outline of the methodology was given in chapter 1.3. The emphasis of the present work is with regard to the time dependence of the energy (or equivalently, the temperature) and estimates of the energy relaxation times. We 12 consider in detail the eigenvalue spectrum of the Fokker-Planck operator versus the mass ratio, y=M/m, where M is the moderator or heat bath species' mass and m is the mass of the "test particle". The approach to equilibrium depends on the nature of this eigenvalue spectrum and in particular on the existence of discrete or continuum eigenvalues. For non-zero Y, we show that the eigenvalue spectrum consists of several "quasi-discrete" eigenvalues plus a continuum. The approach to equilibrium is governed by both types of eigenvalues depending on Y, the initial distribution and the time. The present approach to such problems differs from the finite difference methods employed by other workers<2.6.9.11,79-82). The theoretical approach that is adopted for the time evolution of the test particle (ion or electron) distribution function is presented in the next section, section LT.2. The numerical methods and the relationship with a Schroedinger equation are briefly discussed in sections II.3 and LL4, respectively, whilst section LI.5 deals with a perturbation expansion of the Fokker-Planck operator in powers of Y- The results and their discussion together with a comparison of the analytical work by CorngokK83'84) and by Garrett(85) is presented in section LI.6. LT.2 Time Evolution of the Distribution Function The time evolution of the spatially homogeneous velocity distribution function for the test particle discussed in Chapter I is, (H.2.1) where 13 G(v) = erf ( — ) - ( — ) exp^— ) (II.2.2) and A = 4 7 C n b f / z b m A mM (II.2.3) In equation (TJ.2.2) and (II.2.3); nb and T b refer respectively to the moderator (or bath) number density and temperature; InA is the Coulomb logarithm (see equations (1.4.4 and 5)); e is the electronic charge, and z and are the charges on the test and bath particles respectively. The function G(v) results from the evaluation of the integrals in the expression for the dynamical friction vector and the diffusion tensor with the equilibrium Maxwellian distribution for the moderator (see chapter 1.4 for the specific quantities and integrals). The distribution function f(v,t) is the isotropic distribution function, that is, the /=0 component of the expansion of f(v, t) in Legendre polynomials. We employ reduced units, t'=t/tQ and x for the time and speed respectively, where 3 I ' 2 / 2kT h \ 3 / 2 / M \ 1 / 2 ^^iVhr) m d X = W v • ( I L 2 - 4 ) In equation (H2.4) the reduced speed, x, has been defined with respect to T', which can differ from T b. With these definitions, the equation for the distribution function is, where the quantity s 2=T'/T b is the parameter used to scale the points in the discrete 14 ordinate method(54'55), and Gj(sx) = 1^j2 G(v). This particular definition of G1 is discussed later. The equilibrium solution to equation (II.2.5) is Maxwellian, that is, _ 2 2 _ 2 2 f (x,oo) = C e s x , where C is a normalization constant. If we set f(x, t1) = e s x g(x.t'), the Fokker-Planck equation for g(x,t') becomes = -Lg(x,t') , (II.2.6) where L = F[ D ( X ) !T- B ( X ) !P ( I L 2-7 ) is the Fokker-Planck operator, with D(x)=2s 2xB(x)-2Bix)-M(x) x dx (II.2.8) and G.(sx) s V (II.2.9) The method of solution employed in the present chapter is based on the expansion of g(x,t') in the eigenfunctions, (J^ , of L (see section L2), that is g(x,f) = Xv|>k(x)e-Xkt' , (H.2.10) k=0 15 where the eigenfunctions, fy., are defined by L<|>k=-\<|>k , (H.2.11) and satisfy the orthogonality relation OO J x V ' ^ t o ^ W d x = 5 j k . (n.2.12) With the form of the distribution function given in equation (1T.2.10), the time dependence of the average energy (or equivalently the temperature) is given by OO 1 M = I$s J e _ s 2 x 2 g ( x ' t ' ) x 4 d x ( l m 3 a ) o OO = X c k e ~ X k t ' (LTI.2.13b) k=0 where OO m 2 2 ck = | s ^ j e~s x x4<|)k(x) dx (H.2.14) is the expression for the energy expansion coefficients. With t=0 in equation (LT.2.10) {ak} are recognized as the expansion coefficients of g(x,0). This eigenfunction approach is generally considered to be valid if the spectrum of the Fokker-Planck operator is discrete. 16 JJ.3 Numerical Method The eigenfunctions and eigenvalues of the FPE are determined with the discrete ordinate (DO) method introduced by Shizgal and Blackmore(54«55) and discussed previously in chapter 1.3. The explicit expansion of the eigenfunctions in a basis set is not employed, rather the eigenfunctions are ctetermined at a set of points which coincide with the points of a quadrature procedure based on speed polynomials, Bn(x), orthogonal with respect to the weight function w(x) = x e x (for a good discussion on this particular choice of weight function and a method of generating the speed polynomials and speed quadrature points and weights via recursion relations see ShizgaK56)). The operator, L, is self-adjoint with respect to the equilibrium distribution, x2f(x,oo), and not with respect to the weight function w(x). As discussed at length by Shizgal and Blackmore(54>55), the symmetric representation of L can be constructed by evaluating the matrix representation in the basis set orthonormal with respect to x2f(x,<»), and transforming to the DO representation with the appropriate unitary transformations. The result is that the eigenvalues and the corresponding eigenfunctions, yk(xj) =•/ x? f (x.,oo) /w(x.) <t>k(Xj) are deteimined by diagonalizing the symmetric matrix, Ly, in the (finite) DO basis set given by, N L i j = - r £ B ( x k > [ D k i + 8 j ( k h ( x k ) ] . di.3.1) s k = l where, h(x) = - ['•fO^tf = x ( s 2 _ i) (II.3.2) 2w(x) 2x2f(x,oo) and D^ is the discrete ordinate representation of the derivative operator^54'55). The 17 parameter s is chosen so as to be able to scale the quadrature points appropriately. An important advantage of the DO method of solution is that the eigenfunctions are determined at the quadrature points appropriate for the integral evaluation of the coefficients in the calculation of average values. For the average energy E(t')/E(oo) given by equation (n.2.13), we have oo 0 where wi are the quadrature weights in the quadrature rule (see equation (1.3.1)). n.4 The Equivalent Schroedinger Eigenvalue Problem The transient behaviour of the charged particle system is determined with the expansion of the distribution function in the eigenfunctions of the Fokker-Planck operator. Since the time dependence of the distribution function is governed by the eigenvalues, [\) (see equation (1.2.3)), it is useful to consider the form of the eigenvalue spectrum. Consequently, we consider the nature of the 'potential' in the Schroedinger equation corresponding to the Fokker-Planck eigenvalue problem. If in equations (LU..2.7) and (LTJ.2.11) we change variables to X y(x) = f - ^ g/BOO (II.4.1) 18 together with the definition, <|>k(x) = ec*x**Fk(x), then we can shown that T ~ [V(y) - * k]*F k = 0 , (II.4.2) dy z where and X C[y(x)] = i j l g l d x ' + |lnB(x) . (II.4.4) With equations (U.4.2-4) and the definition given by equations (TL2.8 and 9) and (TJ..4.1) we find that the equivalent Schroedinger potential may be written in the form x dx 4 dx2 16B(x)v dx / • ^ " 15?) " 3 [ l + T - i ? ] e X P ( - Y x 2 ) - ^ 7 e x p < - 2 7 x 2 ) (LT..4.5), where the definition in equations (II.2.8 and 9) (with s=l) has been used. From the equation above and the definition of G l t we have that V(0) =-3(1+ y/2) and V(x)=l/x as x-*». Since the potential is finite at infinity, the discrete spectrum of the Fokker-Planck operator is expected to be empty^3-85). For y = 0, i.e for heavy test particles on light moderator particles, referred to as the 19 Rayleigh limit in kinetic theory<47»48), Gj(x) = x3, B(x) = 1, y = x and V(x) = x2 - 3. For this case, the eigenvalues of the Fokker-Planck equation are given by X =^4n (the Fokker-Planck equation in this limit reducing to the associated Laguerre differential equation, for which the eigenfunctions are the associated Laguerre polynomials, L n ). This result is also familiar to workers in kinetic theory and corresponds to the situation for which the cross-section, o, varies inversely with x(47*48), which corresponds to the Maxwell model collision cross section(47>48). In the study of electron thermalizationC77'78), the function B(x) = xc(x), where o is the electron-atom momentum transfer cross section. For the present study, it is useful to consider some effective "cross section" oeff(x) = B(x)/x, with B(x) given in equation (LI.2.9). n.5 Perturbation Expansion of L about y=0 Much of the discrete portion of the eigenvalue spectrum can be explained in terms of a simple perturbation expansion of the Fokker-Planck operator about the Rayleigh (or y=0) limit. Proceeding in the customary manner (as is detailed in most quantum mechanics texts, see for example P.W.Atkins^86)) we expand the Fokker-Planck operator (given by equation (LI.2.7) with s = 1) in powers of y, L = L ( 0 ) + Y L ( 1 ) + Y 2 L ( 2 ) + ... . (H.5.1) where L ( 0 ) is the Rayleigh limit Fokker-Planck operator and L ( 1 ) and L ( 2 ) are the first and second order expansions of L in y and y 2 respectively. Expanding equations (IJJ.2.8 and 9) as power series in y yields 20 BOO = 1 - f Y x 2 + £ y V + ...+ 2%i± + _ . (II.5.2) so that L W = 3(-)n n r/3 _ W + _ d l l n!(2n + 3) y L\2 'dy ' ( I L 5 , 3 ) where y = x . Since the Rayleigh limit eigenfunctions, ^ \y)=L*/2 (x2), are known it is possible to calculate the mass ratio dependence for the eigenvalues near the Rayleigh limit. For example, the first order correction to the eigenvalue is^86) C - f ^ e - C ^ C , (.1.5.4) and the first order correction to the eigenfunction ist86) - •« + L .(0) >(6) ( I L 5- 5) where X' denotes summation over all k except k=m. To second order in y we find that X„ , n [4 - M | ± i l T - 3 ^ - ^ - 5 9 ) ^ ] _ ^ 21 where the second and third terms correspond to the first and second order contributions (as powers of y) to the eigenvalue. It is also useful to note that the first order perturbation approximation to the eigenfunctions (unnormalized) is given by = L n 2 " f Y [n(n + 1) ll*l - (n - 1) (n + l/2)l£t _ £ ( n - 1) (n + 1/2)1^*x ], (II.5.7) which provides a useful means of comparison with the true Fokker-Planck discrete eigenfunctions. II.6 Results and Discussion The calculation of the relaxation of the test particle energy given by equation (LI.2.13) requires the eigenvalues, X^, and energy coefficients, c^ The eigenvalues and eigenfunctions are (tetermined by diagonalizing the DO matrix representative of the Fokker-Planck operator, that is, the matrix Ly given by equation (LT.3.1). The coefficients ck are determined for an initial delta function distribution at x0. A major interest of the present work is the examination of the energy (or temperature) relaxation obtained from the solution of the FPE in comparison with the temperature relaxation detennined with the assumption that the distribution function remains Maxwellian, that is dT(f) _ T(f) - T b df ' U i ° l ; where x , often referred to as the equilibration time, is given by 22 t = i ( l t m f (H.6.2) «i 4 v mTb ' Equation (LT.6.1) is expressed in the same t' units as used for the Fokker-Planck equation, equation (LT..2.6). An exponential relaxation time is defined as the time x such that T(x) - T H = [T(0) - T(oo)] / e . (IL6.3) For initial temperatures near T b , equation (TJ.6.2) provides an estimate for the exponential relaxation time x. = , where the temperature is replaced with the b initial value and the variation of the temperature with time is ignored. Many workers^1'8'35'45'46) have employed this approach to estimate the thermalization or equilibration times of charged species in ionized gases. We introduce the definitions, 1/2 a 2 = M + i and U(f) = ( l + ^Ef) (0.6.4) m v m l 7 D and solve equation (TI.6.1) to yield, f = a 3 l n ( [ a + U ( t ' ) 1 [ a - U ( 0 ) ] | + ^[U(0)-U(t')l - irU 3 (f)-U 3 (0)l 1 4 I [a-U(t')][a + U(0)] i 2 K n 6 L K > W J (II.6.5) Equation (U.6.5) is a transcendental equation and the relaxation time, %2 (see equation LT.6.3 and put x=x2), determined from it includes the variation of t „ with time. We return 23 to a comparison of the relaxation times Xj and x2 with the relaxation time x determined from the Fokker-Planck equation later in this section. As discussed in Section II.4 and 5, in the limit y=0, the eigenfunctions of the Fokker-Planck operator are the Laguerre polynomials, L n . Since the energy, x2, is 1/2 orthogonal to L n for n > 2, equation (L1.2.13b) reduces to, I f l ^ + c ^ ' 0 (IL6.6) E(<») 1 The result for the relaxation time in equation (JJ..6.6) agrees exactly with equation (LT.6.2) since in the Rayleigh limit x„ is independent of T(t') and is equal to 1 /4 (as is 1A,). We consider the nonequilibrium effects associated with the thermalization process as the mass ratio departs from zero. The transformation of the Fokker-Planck equation to the Schroedinger equation provides a useful interpretation of the variation of the eigenvalue spectrum with mass ratio. Figures II.1A-D show the potential function for four different mass ratios, and the departure from the quadratic form for zero mass ratio. The horizontal lines in the potential functions in Figures LI.1A-D show the positions of these eigenvalues or "states". With the asymptotic expansion of G(x) in equation (IJ.2.2) for large x, the potential V(x) given by equation (LI.4.5) varies as 1/x as x-»oo, and is clearly not bounded for non-zero mass ratios. Hence, the spectrum of the Fokker-Planck operator will have a continuum portion. Diagonalization of the DO representative of the Fokker-Planck operator with N quadrature points (i.e N basis functions) yields a finite number of eigenvalues some of which converge very rapidly as N is increased. The convergence of these eigenvalues versus N is shown in Table LT.l for several mass ratios, and the rapidity of the convergence is clear. 24 Figure II. 1 Potential functions in the Schroedinger equation equivalent to the Fokker-Planck eigenvalue equation. Figure II. 1 (continued) Potential functions in the Schroedinger equation equivalent to the Fokker-Planck eigenvalue equation. 26 Table III Convergence of Eigenvalues^ ) y 0.04 0.09 0.16 0.04 0.09 0.04 0.04 N *2 *6 2 4.257 4.079 4 3.761 3.467 3.774 7.368 6.657 6 3.760 3.462 3.136 7.134 6.076 16.41 8 3.057 7.131 6.022 13.10 30.44 10 3.044 6.017 12.71 19.99 15 3.041 6.016 12.66 16/70 20 16.41 25 16.39 (b) 3.760 3.460 3.040 7.136 6.056 12.74 16.80 (c) 3.761 3.462 3.048 7.132 6.038 12.69 16.62 (a) in units of tfl . (b) First order perturbation result for Xn from equation (JJ..5.6). (c) Second order perturbation result for Xn from equation (LL5.6). 27 We refer to these eigenvalues as discrete although they are not sharp in the quantum mechanical sense. Figures II.2A-D show the Fokker-Planck and perturbation eigenfunctions associated with the four uppermost eigenvalues in Figure LT.1B with Y=0.04. Near x=0 the agreement between the perturbation and Fokker-Planck eigenfunctions is very good, but deteriorates for eigenfunctions corresponding to eigenvalues that are near the maximum in the Schroedinger potential. The somewhat more diffuse nature of the eigenfunction in figure II.2D corresponding to X 6= 16.39 is associated with the topmost 'discrete' eigenvalue and the slower convergence in Table LT.l are indications of the quantum mechanical tunneling of this "state" through the potential. Indeed evidence of this tunneling can be seen in figure LI.2D where following the decay on penetrating the potential barrier (at about x=6.3) the eigenfunction undergoes a low amplitude oscillation on emerging through the barrier at about x=9.5). It is useful to note that the maximum attained by the potential for this mass ratio is about 17.16, just above this (discrete) eigenvalue; see figure II. IB. The remaining eigenvalues that result from the numerical diagonalization vary considerably with increasing N and display no comparable convergence. We consider these eigenvalues as belonging to the continuum. For fixed N, the quantity s = -JT/Tb is varied with the result that the "discrete eigenvalues" remain unchanged whereas the eigenvalues in the continuum vary as shown in Table II.2. Variation of the variable s causes the eigenfunctions associated with the discrete eigenvalues shown in figure H..2 to be either compressed closer to the origin or spread out towards larger x. The value of the corresponding eigenvalue remains unchanged unless the scaling is made too large or too small. It is interesting to note that CorngoldC84) considers the discretum (i.e the discrete portion of the spectrum of L) for non-zero mass ratios to be empty, which is rigorously correct although we prefer to divide the spectrum into a discrete portion (the "long-lived 28 0.4 1 i • i i — r — • — i — • — i — 0.3 -1 °*2 -1 0.1 -? n n OX) S -0.1 --0.2 B -- ° - 3 c • < • • . • . I 1 2 4 6 8 10 12 14 Figure II.2 Eigenfunctions of the Fokker-Planck operator, curve a, and the first order perturbation approximation to the eigenfunctions, curve b (see equation (IL5.7)). T b = 300K; Y=0.04 (A) X 3 = 10.102 and (B) X,4 = 12.659. These eigenvalues are the quasi-bound, "discrete" states shown in the potential of Figure II. IB and converge rapidly as shown in Table II. 1. for some of the eigenvalues. 29 0.3 0.2 J 0.1 1 O.Or 3 1 -O. lH 5 -0.2r--0.3 - 0 . 4 J 12 14 10 12 14 Figure TJ.2 (continued) Eigenfunctions of the Fokker-Planck operator, curve a, and the first order perturbation approximation to the eigenfunctions, curve b (see equation (II.5.7)). T b = 300K;y = 0.04 (C)X5= 14.775 and (D) X6 = 16.394. 30 Table IJ.2 Variation of the Eigenvalues with the scale factor, s X <b> X8 Xg 0.75 14.78 16.39 17.38 18.12 19.20 0.80 17.28 17.75 18.59 0.85 17.07 17.45 18.10 0.90 16.62 17.28 17.67 0.95 16.03 16.98 17.38 1.00 15.40 16.43 17.15 1.05 14.76 15.81 16.66 1.10 14.14 15.16 16.04 (a) N=50; y= 0.04; Xn in units of t 0 \ (b) Eigenvalues X5 and X6 are the two highest states shown in Figure 1 for this mass ratio. The value of the potential at the maximum is 17.10. 31 states" bound in the potential functions as shown in figure LT.2) and a continuum portion that extends from zero to infinity. The time dependence of the energy, or equivalently the temperature, is then given as a sum (actually an integral) over the continuum. These features of the eigenvalue spectrum and the non-convergence of the eigenvalues in the continuum would appear to invalidate the present approach based on a discrete basis set. We find, however, that the time dependence of the energy, E(t') /E(°°) is well converged over a large range of mass ratios and times. For most of the results presented in this chapter, 50-80 quadrature points were employed. For situations close to equilibrium, as few as 30 quadrature points are required. Figure II.3A shows the temperature relaxation for an initial delta function distribution for several different mass ratios with T(0)/Tb=24. Also shown in figure LI.3B is the effective eigenvalue for the relaxation defined by If the spectrum were entirely discrete, the final approach to equilibrium would be a pure exponential and governed by the smallest nonzero eigenvalue. For sufficiently long times, X,eff would coincide with this smallest non-zero eigenvalue. For three of the mass ratios considered in Figure LI.3B (curves a-c), the Fokker-Planck operator has at least one discrete level and X e f f tends to the lowest (non-zero) eigenvalue in each case. These are 3.760, 3.462 and 3.041 for y=0.04, 0.09 and 0.016 respectively. A comparison of Figures LT.3A and B shows clearly that the pure exponential approach to equilibrium occurs when the system is indeed very close to thermal equilibrium. For curve d in figure LT.3A (7=0.25), there are no discrete eigenvalues. A discussion of the results for this mass ratio (II.6.7) 32 Figure H3 Temperature relaxation (A) and the effective relaxation rate (B): T b = 300K, T(0)/T b = 24 and y equal to (a) 0.04, (b) 0.09, (c) 0.16 and (d) 0.25. 33 is presented later in comparison with results for larger mass ratios. The mass dependence of the discrete eigenvalues, relative to the values for zero mass ratio is shown in figure II.4 and they exhibit a strong linear dependence with mass ratio up to the mass ratio for which that discrete level disappears. The first order perturbation approach discussed previously in section II.5 describes this linear dependence very well (see the lines through the points) and little is gained with respect to convergence by going to the second order perturbation approximation as can be seen in Table II. 1. A careful study of the Fokker-Planck discrete eigenvalues versus the perturbation approximation reveals that the most severe discrepancies arise only near the maximum in the equivalent Schroedinger potential, or the last discrete eigenvalue. This is a result of the the failure of the perturbation expansion near the continuum The ratio of the absolute values of the sum over the continuum terms, Sc, to the sum over the discrete eigenvalues, Sd, (not including the XQ term) is shown in figure JJ.5 for several different mass ratios. As might be expected, the continuum makes an important contribution for very small times and decays quickly to zero, especially for the smaller mass ratios. The situation for curve d in Figure n.5 arises owing to the contribution from the discretum being initially negative and passing through zero. The sign of the contribution of the absolute value of the sum over the discretum appears to depend on whether there are an even or odd number of discrete eigenvalues. For the cases shown in figure n.5, the contribution from the continuum decreases with time, at least for the time scale shown in the figures. We will return to this aspect later in the chapter and consider why the present calculations give converged results for the temperature relaxation even though there are a large number of unconverged (continuum) eigenvalues. Strictly speaking equation (11.2.13b) should be written, 34 Figure II.4 Variation of the eigenvalues of the Fokker-Planck operator with mass ratio. The points are eigenvalues obtained from the Fokker-Planck operator, whilst the lines correspond to the first order perturbation approximation to the eigenvalues; see equation (LI.5.6) of the text 35 0 1 2 3 4 5 Figure II.5 Ratio of the absolute values of the sum of the contributions from the continuum spectrum, S c, to the discrete spectrum, Sd. T b = 300K, T(0)/Tb = 24 and (a) y = 0.04, (b) y = 0.0625, (c) Y = 0.09, (d) y = 0.1225 and (e) Y = 0.16; see text for definitions and discussion. 36 (II.6.8) where D is the number of discrete eigenvalues, b the upper bound in the potential of the Schroedinger equation equivalent to the Fokker-Planck equation below which the eigenvalues are discrete. The sum over the remaining N-D terms is replaced by the integral over the continuum eigenvalues which lie at and above b. The reason for the convergence of E(t')/E(<») is that for each different choice of the number of quadrature points, N, and scale factor, s, the eigenvalues and corresponding coefficients, although unconverged, appear to fall on a common c(X) versus X curve. Figures LL6 and 7 show this behaviour. Figure n.6 shows the results for two different mass ratios whereas figures IL7A and B are for Y=0.25 and two different initial speeds x0=4 and 5 respectively. The triangles in these graphs are the numerically determined coefficients, c^ , for N=100 and four different values of s. The solid curve is the line connecting the symbols. It is clear from these results that although the individual eigenvalues Xk and coefficients c± vary with s and are unconverged, the integral in equation (LI.6.8) converges for a large range of times, from very short to very long. One of the more important features is that the first large maximum in the c(X) versus X curve in Figures n.6 and n.7 occurs at a value of X that decreases with increasing mass ratio. This first maximum occurs when X ~ 7.45,4.33 and 2.45 for y = 0.09, 0.16 and 0.25, respectively. This corresponds remarkably well with the maximum value in the Schroedinger potential function, V,,^, for which the above mass ratios yield V m a x = 7.48, 4.10 and 2.55. A useful approach towards an understanding of the behaviour of c(X) versus X and 37 0 5 10 ^ 15 20 Figure II.6 Variation of the energy expansion coefficients, c(X), versus eigenvalue, X, for the continuum portion of the spectrum. Tb= 300K, T(0)/Tb = 24, c(X) in units of 103. The triangles are the numerical results for four different values of the scaling parameter, s. Figure n j Show the variation of the energy expansion coefficients, c(X), versus eigenvalue, X, for the continuum portion of the spectrum. The triangles are the numerical results for four different values of the scaling parameter, s. Tb= 300K, T(0)/Tb = 10.667 and 16.667 in A and B respectively, c(K) is in units of 102. 39 in turn the time dependence of the energy, E(t')/E(«>), is to consider a JWKB approximation to the equivalent Schroedinger eigenfunctions. For eigenfunctions corresponding to "k < V m a x , there are three turning points (i.e where V(x) - ^=0), one at X j = 0, another at X j on the increasing portion of the potential and a third at x3 far from the origin on the decreasing asymptotic portion of the potential (see figure LI.8 for the pictorial representation of these turning points). Most of the eigenfunctions with X < V m a x are concentrated in the region x>x3. Only for several eigenfunctions satisfying the semi-classical quantization condition (II.6.9) V(x) Figure LT.8 Turning points in the potential function of the Schroedinger equation equivalent to the Fokker-Planck eigenvalue equation. 40 are the eigenfunctions concentrated in the region of the innermost part of the potential (shown in figure II.8). In equation (LL6.9) above, S(x) is the phase function given as X J V B(x') S(x) = I / ~ ,\ ' dx' , (H.6.10) and the corresponding JWKB eigenfunction is, -1/4 4>(X;x) = F(X) sin[S(x)] [B(x)] (H.6.11) [X-V(x)]1 /4 x e- 2/2 When V(x) - Xk< 0 S(x) is imaginary and the eigenfunctions decay (i.e for the region between x2 and x3 in figure II.8), whilst when V(x)-X,k>0 S(x) is real and the eigenfunctions oscillate like sinusoids. This accounts for the qualitative differences between the discrete and continuum eigenfunctions shown respectively in figures LT.2 and 9. The numerically determined eigenfunctions are the eigenvectors of a (finite dimensional) symmetric matrix (equation (LT.3.1) and are (by definition) L2-integrable 2 „2_2 functions, with weight function x e . The JWKB function given by equations (LI.6.9) and (LI.6.10) are non-normalizable functions over the semi-infinite interval [0,<»]. However, a point by point comparison of the eigenfunctions that result from the diagonalization of the DO representative of the Fokker-Planck operator and the corresponding JWKB functions shows excellent agreement to within a normalization constant which is denoted by F(X) in equation (LI.6.11). Figure II.9 shows two such 41 0.3 0.3 8 12 ' X Figure II.9 Comparison of continuum eigenfunctions detenrtined from the diagonalization of the matrix representative of the Fokker-Planck operator and the JWKB approximate functions.Tb=300K, s = 0.6042, X = 0.4; (A) X = 4.855 and (B) X = 7.603. 42 continuum eigenfunctions for y=0.09 with eigenvalues equal to 4.86 and 7.60, somewhat above the maximum of the potential which is 4.10. The triangles in figure LT.9 are the results obtained with the diagonalization of the matrix representative of the Fokker-Planck operator, and the solid line is the corresponding JWKB function. The normalization of the JWKB function was chosen to give the agreement that is shown. If the present numerical quadrature procedure is used to determine the normalization of the JWKB eigenfunction, the norm that results is somewhat different from that obtained by the point-by-point comparison, owing to discrepancies in the last few quadrature points. Nevertheless, the excellent agreement between the numerically determined eigenfunctions and the JWKB functions provides added support for the success of the present numerical work. The r' normalizable eigenvectors determined from the diagonalization of the DO representative of the Fokker-Planck operator are excellent approximations to the continuum eigenfunctions. The approximation of continuum states with L2-integrable functions, has also been employed by LanghoffC87) and Reinhardt^ 88) in the quantum mechanical study of photoionization processes. For the initial delta function distribution we have that a k = <t>k(x0) in equation (JJ.2.10), and the oscillatory behaviour of the curves in figures n.6 and LT.7 reflect the oscillatory variation of these continuum eigenfunctions with X. The amplitude of the c(X) oscillations decreases with increasing X owing to the decrease in the integrals of w(x)x2<|>k(x) in equation (LL2.14). Of particular interest is the way in which the approach to equilibrium might be different for mass ratios where there is one or no discrete eigenvalue. As discussed earlier in this chapter, if there is one discrete eigenvalue, the approach to equilibrium will eventually become exponential, governed by this eigenvalue. Figure JJ.10A shows the 43 -9 • 1 « 1 1 • 1 0 2 t y T 4 6 Figure LI. 10 Effective relaxation rate (A) and temperature relaxation (B); Tb= 300K, T(0)/Tb=24 and X equal to (a) 0.25, (b) 1.0 and (c) 25.0. x is the time in units of tQ for the temperature difference T(t') - T b to decrease by 1/e of its initial value. 44 effective eigenvalue defined with equations (II.6.1) and the time dependence of log([T(t') - T^ ] / T ^ for the larger mass ratios y = 0.25,1 and 25 for which there are no discrete eigenvalues (except XQ). Since the time dependence is considerably different for these mass ratios (much slower for y=25 than for y=0.25), the reduced time is scaled with the exponential relaxation time x defined in equation (LT.6.3) so as to show the results for all three mass ratios on the same graph. The time scale is very large in these cases since at t'/t = 5, the values of (T(t) - / T b are of the order of 10"9 to 10"7. It is clear from the results in Figure II.10B that the approach to equilibrium is qualitatively close to exponential, and almost independent of mass ratio for this choice of time scale. This independence on y arises from the fact that the potential function V(x) or equivalently B(x) or G(x) attain large mass ratio asymptotic forms for mass ratios not much greater than 1. For large mass ratios, we have that G l W - | j | and so V(x) = ^ [ 1 - , (n.6.12) where we can show that V(x) has a maximum at x = 1/45/2 = 1.295 for which V m a x = 0.8212/y3/2. It happens that the exponential relaxation time x varies inversely with y 3 / 2 for y somewhat greater than unity. The present results with regard to the long time dependence remain qualitative and a more definitive interpretation would require extraction of the analytic form of the c(X) versus X curve from the JWKB approximation. It is however important to mention that from a practical point of view that this is a somewhat academic question since this asymptotic time dependence occurs when the system is extremely close to thermal equilibrium. For sufficiently long times, the final approach to equilibrium will be 45 determined by the small X portion of the c(X) versus X curve. This will also be the case for mass ratios for which the discretum is not empty. Calculations of the temperature relaxation for such mass ratios at extremely long times show that the contribution from the continuum eventually increases and becomes larger than the contribution from the discretum. For such long times the eigenvalues associated with eigenfunctions concentrated to the right of the rightmost turning point dominate the relaxation. The curves shown in figure 1T.3B extrapolated to very long times eventually begin to decrease from the flat portions, and the ratios Sg/Sjj shown in figure LI.5 also increase. Table LT.3 compares the relaxation times, x, calculated with solutions of the FPE with those calculated from equations.(II.6.2 and 5) denoted by %v and x2 and defined earlier in this section. There is complete agreement among all three relaxation times for the Rayleigh limit, (y=0), as discussed earlier. The departure of the nonequilibrium relaxation time x from the other two estimates increases with increasing mass ratio. x2/x attains a maximum value for mass ratios near unity for the smaller initial temperature ratio and at smaller mass ratios with increasing initial temperature ratio. For the larger mass ratios, x varies inversely as ,y 3 / 2, as discussed earlier and since for large T(0)/Tb, xl (and presumably x2) vary as y3?2 (equation (LI.6.2)), the ratio x2/x tends to a constant value independent of Y. It is important to note that in many applications it is xl that is used as an estimate of the relaxation time rather than x2. The departure of x from x2 is considerably greater than from x2. The graph of x versus Y, shown in figure 11.11, for the results in Table LT.3 has two distinct linear portions (on a logarithmic scale), one for y< 1 and one for Y> 1. For the larger mass ratios, the slope is very close to 3/2, whereas for the smaller mass ratios the slope depends on the initial temperature ratio. 46 Table n,3 Comparison of Relaxation Times^ Y *2 t2/x T(0)/TM = 2,667 0.01 0.254 0.260 0.258 1.016 0.04 0.265 0.291 0.281 1.060 0.09 0.284 0.345 0.323 1.137 0.16 0.313 0.426 0.383 1.224 0.25 0.354 0.538 0.466 1.316 1.00 0.858 1.755 1.389 1.619 4.00 5.483 9.962 7.037 1.283 25.00 85.28 139.2 95.39 1.119 100.0 682.3 1095 747 1.095 T(0VT(~) = 10.667 0.01 0.266 0.291 0.277 1.041 0.04 0.319 0.426 0.364 1.141 0.09 0.425 0.686 0.527 1.240 0.16 0.616 1.113 0.786 1.276 0.25 0.935 1.755 1.168 1.249 1.00 6.006 9.962 5.875 0.978 4.00 47.72 72.14 40.78 0.855 25.00 745.5 1095 610.8 0.819 100.0 5964 8722 4857 0.814 47 Table LI.3 (continued)  Comparison of Relaxation Times Y X x2lx T(0)/T(oo) = 24.0 0.01 0.287 0.345 0.311 1.084 0.04 0.424 0.686 0.518 1.222 0.09 0.755 1.404 0.937 1.241 0.16 1.435 2.662 1.649 1.149 0.25 2.620 4.630 2.745 1.048 1.00 20.28 31.25 17.21 0.849 4.00 162.1 238.8 128.7 0.794 25.00 2532 3683 1973 0.779 100.0 20260 29410 15740 0.777 (a) Relaxation times are in units of 1 0 and represent the times required for T(t') - T(«>) to decay to 1 / e of the initial value; x is determined from the solution to the Fokker-Planck equation with an initial delta function distribution; xl is the relaxation time given by equation (LI.6.2) with T = T(0); x2 is the relaxation time as determined with equation (II.6.5). 48 Y Figure 11.11 Logarithmic plot of x versus yfor the results in Table U.3. Lines a-c are the Fokker-Planck relaxation times, x, for T(0)/T b = 2.667, 10.667 and 24 respectively, whilst lines e-f are the Spitzer relaxation times, Xj, for these T(0)/Tb values. 49 II.7 Conclusions The results of the present study indicate that expansions in a finite polynomial basis can yield converged solutions to the energy (or temperature) relaxation despite the presence of a continuum eigenvalue spectrum. In particular the convergence of the energy relaxation was shown to be a result of the validity of the approximation of equation (LTJ.2.13b) to equation (II.6.8). The presence of a c(X) versus X curve (figures II.6 and 7) was considered in the light of a JWKB approximation to the Fokker-Planck continuum eigenfunctions and the eigenvalues were classified according to the form of the potential in the Schrc«dinger equation equivalent of the Fokker-Planck equation. The relaxation of the energy was shown to give good agreement with the Spitzer relaxation times for a wide range of mass ratios and initial energies. The long time approach of the system to equilibrium for mass ratios possessing no discrete eigenvalue spectrum was shown to be qualitatively close to exponential and was attributable to the dominant part of the c(X) versus X curve at long times being peaked in the vicinity of the maximum in the equivalent Schroedinger potential. In the next chapter we solve the non-linear "test particle" problem for an initially anisotropic one component plasma relaxing via self collisions. The spectral properties of the linearized operator for this problem are briefly considered and the convergence of the anisotropic energy (or temperature) is explored in more detail in light of the present findings. 50 Chapter JJJ Anisotropic Temperature Relaxation due to Coulomb CoUisions in a Spatially  Homogeneous One Component Plasma This chapter is concerned with the collisional relaxation of the anisotropic temperature components, T, =m<v2>/k andTx =m<v 2>/2k in an initial bi-Maxwellian single component plasma (the quantities <v2> and <v2> denote averages of v 2 and v 2 over the ion velocity distribution function (VDF)). The plasma is assumed to be spatially homogeneous and collisional relaxation occurs in the absence of external electric and magnetic fields. The Fokker-Planck equation (FPE) describing this system is non-linear and is linearized by expanding the ion VDF, f, about its equilibrium Maxwellian. Results from this study are compared with several recent experimental results as well as recent theoretical solutions of the full non-linear Fokker-Planck equation. The method employed here for solving the FPE is based on the generating function approach for the associated Laguerre polynomials introduced in chapter I and developed further in the work by Wong(67) and in appendix 1. m i Introduction The collisional relaxation to equilibrium of a spatially homogeneous single component anisotropic plasma in the absence of electric and magnetic fields is an important problem in plasma physics. Such anisotropic distributions are common to many trapped ion cells( 1 6' 1 8> 9°), occur in ion cyclotron resonance heatingO7-91-93), the solar wind(27-34), ionosphere and many aeronomical phenoma(33-3$. These plasmas are characterized by two temperature components so that the plasma species is hotter along one direction than in a 51 perpendicular direction. This non-equilibrium anisotropic distribution in the plasma relaxes in the presence of collisions owing to a transfer of energy from the translational degree of freedom defining the hotter plasma temperature to that defining the cooler temperature until a state of equilibrium is reached. A simple anisotropic single component plasma of physical interest which possesses azimuthal (or cylindrical) symmetry is characterized by two temperature components, T ( | =m<v2>/k and T ± = m< v2>/2k. We choose to describe this plasma by the bi-Maxwellian ion VDF, whose form is written explicitly in the next section. This plasma relaxes collisionally via energy exchange between the degrees of freedom in the plasma. The equuibrium temperature of the plasma, T, is given by, T = T" "^ 2 T x (III.l.l) The degree of anisotropy for the bi-Maxwellian ion VDF is most conveniently described in terms of the temperature parameter, a(t), given by a(t) = i|L = l/e(t) , (in.1.2) where the rate of relaxation of the initial bi-Maxwellian ion VDF depends on the initial temperature parameter, a(0). The present research parallels the work by Livi and Marsch(3°) and the work by Jorna and Wood^11) who employed the non-linear Fokker-Planck operator. For the theoretical study in this chapter the linearized approximation to the non-linear equation is used. This approximation is based on an expansion of the ion velocity distribution 52 function, f, about its equilibrium Maxwellian, f*°\ and is discussed more fully in the next section. It is important to point out that the linearization facilitates an expansion of the operator in terms of its eigenfunctions and eigenvalues. The procedure for this was introduced in chapter I of this thesis and required the operator to be linear. The remainder of this chapter is organized as follows: Section ITJ.2 Deals with the theoretical aspects and is divided into four parts. Part (a) focuses on the linearization of the Fokker-Planck equation. Part (b) establishes a connection between T„(t) and the time dependent expansion coefficients. Part (c) connects the time dependent expansion coefficients with the eigenvalues and eigenfunctions of the linearized Fokker-Planck operator as well as their initial values. Also included is the dependence on the initial expansion coefficients for which an expression is given in part (d). Section LII.3 Deals with the results and discussion and compares the present with other theoretical and experimental work. The velocity distribution function for the single component plasma studied in this chapter is expanded in a Burnett function basis (introduced in chapter I) about Hence the eigenfunctions of the linearized operator are linear combinations of these basis functions and the eigenvalues are calculated by diagonalizing the matrix representation of L in the Burnett function basis. The means for obtaining the matrix elements is based on the generating function approach introduced in chapter I and outlined in appendix 1 where some low order matrix elements are given. 53 HI 2 Theoretical Approach (a) Formulation of the Problem and the Linearization of the Fokker-Planck Equation The Fokker-Planck equation for a spatially uniform plasma undergoing self collisions in the absence of external electromagnetic fields may be written in the form - | f = C[f,f] , (in.2.1) where C[f,f] =l2LeVA a f d v . u 2 i - u u a t r |L _ f ar \ ( m 2 2 ) m2 ov J u 3 V dv dv'/ is the Landau form for the collision term(5'67'72), u =v-v' is the relative velocity, and InA is the Coulomb logarithm, previously introduced in equation (1.4.4). Note that f = f(v',t) and f = f(v,t) in equation (HI.2.2). The Landau form for the collision term is readily derived from equation (1.4. lb) by making use of the substitution JL - I J . f u2l-uiA  u 3 ~ 2 dy'' V u 3 ) in the expression for the dynamical friction vector, A (equation (I.4.2a)), and then integrating by parts. In this chapter we wish to solve equation (HL2.1) subject to the initial conditions 1/2 2 2 (m.2.3) 54 where; n is the number density of the plasma species, T,(0) and T±(0) are the initial parallel and perpendicular temperature components, and v and v± are the components of the projection of the ion velocity vector, v, onto the axes defining T„and T r We choose to define the symmetry axis of the VDF as the z axis which is also the axis defining T r Equation (JJI.2.3) is a bi-Maxwellian ion velocity distribution function, that can be described in terms of the initial temperature anisotropy, a(0), defined in equation (in. 1.2). It is useful to note that f(v,0) is cylindrically symmetric about the z axis defining TH(0). As it stands equation (LTJ.2.1) is non-linear in f and hence cannot be solved by the techniques of linear algebra. We therefore linearize it by expanding f as f = f ( 0 ) + f ( 1 ) , (IH.2.4) where f<0)« • Kw) 3 ' V r a v V 2 " <ra-2-5> is the equilibrium Maxwellian ion velocity distribution function, T the equilibrium temperature previously defined in equation (III. 1.1) and f^sf ^(v,t) is a time decaying perturbation describing the anisotropy in the ion velocity distribution function and is required to vanish at equilibrium. Substitution of equation (IJI.2.4) into equation (HI.2.1) yields i t = c[f(1!f(0)] + c[f(0!f(1)] + c[f(0!f(0)i + ctf(1!f(1)] . (in.2.6) dt The term C[&°\ f*0*] is equal to the rate of change of the equilibrium ion velocity 55 distribution function and is trivially zero, whilst the term Ctf .^f*1*] is the non-linear collision term, which is small for small departures from equilibrium. The linearization of equation (LTI.2.1) with (IJJ.2.4) simply works by discarding this term. Consequently, the linearization reduces equation (LLI.2.1) to W. = C [ f d ) fOTj + c [ f <P> fd)j = ( L d + £ l ) f = L f (in.2.7) where L is a linear operator (for the formal solution of equation (IJJ.2.7) see the discussion in chapter 1.2). In equation (HI.2.7) f may be replaced by f ( 1 ) since d&0)/dt = Itf® = 8. As noted by Hinton(5), the first term on the right hand side of equation (HI.2.7), C[&1\ f*0)] =L^^\ is a differential operator acting on f because the integral over v1 in equation (LII.2.2) is easily performed when the argument of f*0) is v'. In plasma kinetic theory these integrals, which appear in the expressions for the dynamical friction vector and the velocity diffusion tensor in equations (1.4.2), are typically formulated in terms of the Rosenbluth potentialst2-5'7), so named because they are solutions of Poisson's equation, like the electrostatic potentials. The details of these important potentials are not given here, instead the reader is referred to the literature by Rosenbluth et aK2'3), Hinton(5) and Shkarofsky et aK7). We note in passing that Lp (for isotropic velocity distributions of test particles) was the form of the Fokker-Planck operator considered in the last chapter. The second term in equation (HI.2.7) is an integral operator, Qf^.f*1)] =Ljf acting on f T h i s is because only the operations involving 9/9v in equation (HI.2.2) can be performed, leaving integrals over v'. Consequently equation (LTI.2.7) is an integro-differential equation whose solution is approximated by expanding f^ in Burnett 56 functions. For the initial bi-Maxwellian ion velocity distribution function considered here the Burnett functions are rewritten from equation (1.3.2), \|/^ (x) = x 4+1/2(x2) E (cos6) . (LTI.2.8a) These functions form an orthogonal set as follows, n-WJdxyyxivPwe-*2 = h / i n8 / r8 n n , where h = 2 r(/ + n + 3/2)  />n 2/+1 / I n ! and ~ < f t f ' is the reduced velocity. The statistical mechanical averages of over the ion velocity distribution function, f, are referred to as moments, and of particular interest in this context are the time dependent moments associated with We expand f ^  in Burnett functions as oo f ( 1 )(x,t) = <D(x,t)f(0)(x) = X u n / (t)¥2)(x)f(°)(x) . (ffl.2.10) nj = 0 where u j(t) is a time dependent expansion coefficient of the function f ( 1 l In order to (ffl.2.8b) (LTI.2.8c) (IJI.2.9) 57 follow the details of the evolution of the ion velocity distribution function from its initial bi-Maxwellian form to its equilibrium Maxwellian form, the time dependence of u .^(t) is required. In particular the time evolution of the n = 0,1 = 2 moment (i.e the moment of 1/2 (3 v,2- v2)), which is associated with Tu(t), requires the computation of the time evolution of UQ2(t). This is shown explicitly in part(b) of this section. Substitution of equation (HI.2.10) into equation (HJ.2.7) yields for the moment, OO J d V V- I = j d V ^ C [ f ' f ] S X Un./(t) Kn + < » ) (m.2.11) n = 0 where the = symbol has been introduced because of the linearization. The following definitions also apply, < » - J * *J? C [<r<° f<0,» f <»>] - Jdv i " v"' ^ v"' (in 2.12) C - jdv V « C[f<0,» ¥ » f ( ° » ] - Jdvf<»V«> I, V « (in.2.13) These definitions are equivalent to the definitions used by Wong(67). The computation of the M and N matrix elements is outlined in appendix 1 of this thesis and in the paper by Wong(67). The matrix elements are calculated using the generating function approach based on the associated Laguerre polynomials. This approach was introduced in chapter I and is also expanded upon in appendix 1 of this thesis. The following parts to this section deal with the relaxation of the anisotropic 58 temperature, Tn (t), and feature expressions facilitating the computation of u0 2(t). (V) Relaxation of T/O The objective here is to derive an expression for the time dependence of Tu(t) in terms of the time dependent expansion coefficients, u^ t^). The relaxation of Tu(t) and u (^t) are expected to depend on the eigenvalues and eigenfunctions of L as well as the initial conditions. This subsection deals with the explicit relation of T|((t) with u^ t^), whilst the subsequent subsections deal with the explicit dependence of ii^ Ct) on the eigenvalues and eigenfunctions of L as well as the initial conditions. The procedure through which T„(t) is derived starts from the expression for the average parallel component of energy, E„(t), which can be written in terms of T„(t) as, E^t) = |kT„ (t ) = Jdv imv,f f (v,t) , (IH.2.14) where v„ = v cos8 = vfj. and vL = v (1 - \i2)m (IH.2.15) are, respectively, the projections of v onto the axes defining T„ and T x in terms of |i, the cosine of the angle 6 between v and the axis defining Tu. Substituting these expressions as well as equations (UL2.4), (JH.2.5), (LTL2.9) and (IJJ.2.10) into equation (JJJ..2.14) yields dx e"x2 [ l + X u r /<l> Yr(,)(*)Vx2 , (IH.2.16) T,l = o where n^ 2 is simply a linear combination of the three Burnett functions, 59 u 2x 2 = £v®<x) + ^V 0 0 )(x) - jvj0 )(x) (ffl.2.17) so that upon utilizing the orthogonality of these functions (see equation (I.3.8b and c)) it can be shown that, m = l + u0>2(t) + u 0 0 ( t ) +u1>0(t) . (in.2.18) Equation (m.2.18) can be further simplified by invoking the additional constraints of conservation of particle number and energy (the violation of the conservation of either of these quantities can lead to changes in the rates of decay via changes in the particle number density (i.e n—>n(t)) and changes in the collision frequency due to a change in T). These conservation properties require that the moments of y*® andx}^ vanish so that uo o )^ = u i o^ ) = ® f° r ^  times throughout the relaxation. Hence equation (Ln.2.18) reduces to M2 = l+u 0 2(t) , (ffl.2.19) which requires the determination of the expansion coefficient u02(t). It is useful to note that at equilibrium T,(<»)/T = 1 implying that U Q 2 ( ° ° ) = 0. Additionally at equilibrium we can expect all expansion coefficients u^/00) - 0, so that f ^ v . o © ) =0 and the system is Maxwellian as required. In the ensuing part (c) to this section it is shown how the time dependence of the expansion coefficients (and hence that of Tu(t) and T±(t)) depend on the initial conditions 60 and the linearized Fokker-Planck operator, L, via its eigenvalues and eigenfunctions. (c) The Time Decay of the Expansion Coefficients Although the previous part to this section established a relation between Tu(t) and u02(t), it did not establish a connection between either of these quantities and the initial conditions or the eigenvalues and eigenfunctions of L It is the purpose of this section to establish such a connection. The time dependence of the expansion coefficient, un / (t) is obtained first by recognizing that 3f/dt = df (1)/3t so that substituting equation (1TI.2.10) into equation (IJI.2.11) yields, equation (1.3.4). It is convenient to introduce the matrix elements A whose components A ^ are defined by (III.2.20) where h / j n is a normalization factor for the Burnett functions previously introduced in A(/) _ t K ^ O m ' k nJh. h / l r * /,m /,k (in.2.2i) where 3/m (kT) 3/2 T = (IH.2.22) 61 is the self-collision momentum deflection time. It is important to note that n/i is the basic time scale associated with the Fokker-Planck collision term and that it appears explicitly in the expressions for the M and N matrix elements (see appendix 1). The above choice for the A matrix elements, which are dimensionless, necessitates the transformation to dimensionless time, t' = t/x. The definition in equation (ILT.2.21), enables equation (HJ.2.20) to be written in the form (') oo #^ - I ° > A « , (m.2.23) k = 0 where cmV) = V % „ u m > / ( t ) . (IH.2.24) Equation (HI.2.23) may also be written in matrix form as ± c(/)(f) = A ( / ) .c ( / )(f) . (IH.2.25) dt' In order to solve equation (HL2.25) we must transform to the a new set of coefficients, c', of the time derivative operator. These coefficients necessarily satisfy the eigenvalue equation c1 (,)(f) = - A ( / > . c' ( / )(f) , (IH.2.26) where A is the eigenvalue matrix of A . We can express A in terms of A via the 62 similarity transformation [ U ' V ' . A ^ U ^ - A ^ , (IH.2.27) where U ( / ) is the matrix of eigenvectors of A ( / )and [U ( /*J _ 1 the inverse of U ( / ) . This enables the establishment of the relation c , ( / )(f) = [U ( /l~ 1.c ( / )(f) • (in.2.28) Solving (IH.2.26) for c' yields c , ( 0(f) = e- A ( , ) tlc , ( / )(0) , (IH.2.29) where c'^ (0) is the initial new coefficient vector. Since the matrix elements of are Ap ^  = Xp * 8p equation (LTI.2.29) can be written in component form as c;(/)(t') = e ^ ' c ^ O ) . (IH.2.30) We regain the original coefficients via the inverse of 001.2.28) to yield c ( Z ) (f) = U ( 0 . e _ A llcr<o(0) (m.2.31a) = V^.e-^WY^iO) , (JJ.L2.31b) 63 or, in component form, r s r where S Equations (HI.2.32) and (HI.2.33) have exploited the diagonal symmetry in the matrix and [U^ l^ _ 1 = [ U ^ T . This relation between the inverse eigenvector matrix and its transpose arises only when A ^ is a real symmetric matrix. As a consequence of this the computation of the inverse matrix is unnecessary. The importance of equations (HI.2.32) and (HI.2.33) is that the time dependence of the c^ (t') coefficients (and hence the u^Ct) coefficients and T„(t)) depend on the initial coefficients and the eigenvalues and eigenfunctions of the linearized Fokker-Planck operator, L. It was the objective of this subsection to establish these connections. Since the eigenvalues and eigenfunctions of L are obtained by diagonalization of the A matrix introduced in equation (1TJ.2.21) (see appendix 1 for the calculation of the M and N matrix elements) all that is required now is the computation of the initial expansion coefficients u n j(0). This is done in part(d) which now follows. It should be noted that equations (m.2.32), (m.2.33), (ffl.2.19) and (ILL2.24) enable T„(t) to be written explicitly in terms of the eigenvalues and eigenfunctions of L as 64 — ' = 1 + U 0 , 2 ( T , ) = 1+Xbr2)e"^r 1 ' (m-2-34) where a(2> oa<2> b j 2 ) = -^==r = —7=r • (m.2.35) 1 J \ o / 3 (2) / (2)\ The functional form of b (X ) is of interest in connection with the eigenvalue spectrum of L and the long time approach to equilibrium of the system (d) The Initial Expansion Coefficients u^ CO) The initial expansion coefficients u f^O) are calculated from the expression for the bi-Maxwellian ion velocity distribution function given in equation (LTI.2.3). Transforming this equation from a velocity distribution function into a reduced velocity distribution function by using equations (LTI.2.15), (HJ.1.1) and (HI. 1.2) for t = 0, yields the reduced bi-Maxwellian velocity distribution function, (UI.2.36) Substitution of the equation for the expansion of f about its equilibrium Maxwellian solution (equation (LTI.2.4)) and equation (LTI.2.10) for the expansion of f^ 1^ gives for the moment of \x) (as a function of time), Jdx f(x, t) Vr°(x) = n $ Q8 Q + h, rur,(t) (IH.2.37) 65 where ht r is defined in equation (LTI.2.8c). When equation (LLT.2.36) is substituted into equation (LTL2.37) for t=0, the following expression for ur/(0) is obtained: Since the integration over the polar angle, (|>, gives 2n and because only ur 2(t) is required for the time evolution of T„(t) (see equation (LTI.2.19)) we have, oo 1 U r.2 ( 0, = F 7 | = J Jdxdtl/Lf(x2)P>,exp{-x2[^][^-a(0)(l-^)]} 2,r o -1 (m.2.39) Consequently the time dependence of the c (^t') coefficients can be determined from equation (LTI.2.32), where equation (LTL2.24) relates ^ (^0) to cj '(0). Equation (III.2.39) can be evaluated using a quadrature rule. For example the integration with respect to the x variable may be performed by using a speed quadrature 1 - X 2 rule (based on the weight function w(x) = x e ), whilst the integration over \i is most conveniently performed using a Legendre quadrature rule. It is, however, possible to express the expansion coefficients in the form 66 where q - 1 5 V (M(-3)m(m +l)! 1 } %2 Vm; (2m+5)!! (m.z.4ij m = 0 is the a-independent part of u ,(0). The double factorial function appearing above is defined by (2m + 5)!! = (2m + 5)(2m + 3)...3.1, and = (T_^y m , is the combinatorial function. The details of this derivation are given in appendix 2 of this thesis. Equations (LU.2.40 and 41) constitute a much more rapid and accurate method than the numerical integration of equation (UI.2.39) using speed and Legendre quadratures, since the q r 2's need only be generated once. However, for the temperature anisotropics studied in this chapter, where convergence is achieved for small basis set sizes, the use of a quadrature rule, such as that suggested for equation (LTJ.2.39) is just as satisfactory. We are now in a position to analyze the various features of the relaxation of T„(t) to equilibrium This is the purpose of the next section. IIL3 Results and Discussion The emphasis of this section is the comparison of the theoretical results, obtained from the prescription given above, with recent experimental results^6-18) and other theoretical results obtained by employing the full non-linear Fokker-Planck operator^  ^ O ). Li particular, this section concentrates on the relaxation times from the anisotropic electron plasma considered by Hyatt, Driscoll and Malmberg (HDM) in their magnetic mirror apparatus^ 6) as well as relaxation times for the non-linear Fokker-Planck operator studies by Joma and Wood(n) (JW) and Livi and Marsch(3°) (LM). Prior to a discussion of these important comparisons some of the more immediate features of the theoretical results from 67 this work are given. A numerical code for the solution of T, (t) as a function of the initial conditions and for the linearized Fokker-Planck operator formulated in equation (LTL2.7) was written. It was found that the diagonalization of the matrix A gave rise to eigenvalues that change with the basis set size, N, as shown in Table III.l. These eigenvalues show no convergence with increases in the basis set size leading to the conclusion that the eigenvalue spectrum is continuous. Further support for this conclusion is gained by an inspection of the eigenfunctions shown in figures LTXIA and B. These eigenfunctions have been chosen to be orthonormal on the interval [0,oo] with respect to unit weight function. What is important here is that in both figures the curves labelled 'a' consists of many different eigenfunctions that fall on the same curve as explained in the caption. These eigenfunctions correspond to different eigenvalues and are typical of continuous eigenfunctions, which possess many nodes and do not decay. Discrete eigenfunctions on the other hand possess an integral number of nodes and decay due to the penetration of a potential barrier. Despite the nonconvergent behaviour of this eigenvalue spectrum the macroscopic observables, such as the temperature relaxation, converge for comparatively small basis set sizes, N, (usually fewer than 10 Burnett functions seems to be sufficient for a converged estimate of the relaxation time). Table HL2 shows the convergence of the relaxation time, ^  versus N, where T„(tt) - T = [T„(0) - TJ/e defines tj. The results lead to the conclusion that tj is rapidly convergent (to four significant figure accuracy) except for large oc(0)'s where it diverges (eg.for a(0)=10 tj diverges for N>25). This is illustrated for a(0) = 8 and 10 in figures LTJ..2A and B respectively. Figure LTJ.2A is seen to possess a 'flat' region suggesting a region of convergence, whilst there is no 'flat' region in figure HJ.2B and tT quickly diverges. As a consequence of this, for a(0) = 10 is estimated as 1.45± 0.02. A likely explanation for the divergent behaviour can be found in figure HI.3 where 68 Table mTl Variation of eigenvalue. V . vs r. as a function of basis set size. N N 5 10 20 30 50 r = 0 0.2690 0.1495 0.0913 0.0705 0.0520 1 0.4548 0.1968 0.1055 0.0781 0.0557 2 0.7637 0.2621 0.1209 0.0857 0.0591 4 3.3712 0.5098 0.1600 0.1026 0.0660 9 8.2108 0.3897 0.1682 0.0464 14 1.4888 0.3218 0.1158 19 19.978 0.7883 0.1633 # In units of 1/x. 69 0.3 x Figure ITJ.1 Eigenfunction plots. Number of basis functions used in A is 40, B 80. Curve 'a' consists of the lowest 15 eigenfunctions for A and the lowest 35 eigenfunctions for B. Curves b and c, in both figures, consist of eigenfunctions with larger eigenvalues than those used in curve a. 70 Table IJJ.2 Convergence of Relaxation times. tr. as a function of a.(0) versus basis set size. N* N 1 2 3 5 10 20 a(0)=0.00 1.179 1.150 1.145 1.153 1.152 0.10 1.167 1.166 1.170 0.25 1.191 1.193 0.50 1.226 1.229 1.228 2.00 1.348 1.344 1.342 4.00 1.419 1.410 1.400 1.399 6.00 1.456 1.446 1.431 1.426 8.00 1.478 1.468 1.450 1.442 10.00 1.493 1.485 1.464 1.452 # In units of x. 71 N 50 60 Figure ILT.2 Variation of the relaxation time, t r with basis set size, N. Note that both figures show that possess a region of convergence. Most relaxation times are satisfactorily estimated with N=10. Figure B leads to the estimate that tT=1.45±0.02. 72 Figure IIJ.3 Variation of the initial expansion coefficients, u r 2(0) with r.Curves a-f correspond to a(0) = 2, 3,4,5, 6 and 10 respectively. Note the growth in the coefficients for oc(0) > 5. 73 the initial expansion coefficient, ur2(0), relevant to the relaxation of TH(t) is plotted as a function of r for six different values of a(0). The logarithmic scale used for the y axis indicates that values of ct(0) £ 5 are associated with expansion coefficients that grow rapidly with r, whilst values of a(0)^4 are associated with expansion coefficients that diminish in size with increasing r. This growth in the expansion coefficients for ce(0) £ 5 results in the non-convergent behaviour of tt shown in figures HL2A and B and the (2) / (2)\ oscillatory behaviour of the temperature coefficients, b (X ), shown in figure LTJ.4B. No such divergence in tt is found for a(0) <4, where the expansion coefficients quickly become small. In connection with the continuous eigenvalues spectrum and the aforementioned (2) temperature coefficients, br , appearing in equation (UI.2.34), we find it convenient to consider equation (IJJ.2.34) as being a discretization of ^=i+J>(x(V*%x<2> . (m.3.i) 0 Figures HJ.4A and B show the functional form for this temperature coefficient (2) ( (2)\ (2) b [X ) as a function of X . For small values of a(0), the coefficients trace out a smooth curve (see figure HI.4A) just as in the previous chapter (see figures n.6 and 7), although the qualitative features of the curve here are very different. For large values of rx(0) the smooth hump distribution deteriorates into the highly oscillatory plots shown in figures HL4B. The is probably attributable to the large values of the expansion coefficients, 74 Figure IIJ.4 Temperature coefficients, b(2)[A.(2)] as a function of X ( 2 ). A contains plots of bv W '] for a(0) = 0, 0.2, 0.4....2.0. B is for a(0) = 6.0. Both plots used 80 basis functions. 75 ur 2(0), previously mentioned in connection with the convergence of tT as a function of N. It is somewhat surprising that these temperature coefficients can yield converged values of Figures IJJ..5A-F show the relaxation of Tn(t')/T, a(t') and e(t') (defined in equation (JJI.1.2)) as a function of t' for the initial conditions as described in the figures. Note that the relaxation rates for a> 1 are smaller than the relaxation rates for rx< 1. Similar differences in the fundamental relaxation rates for the a(t') curves are evident from figures LU.5C and E. The explanation for these differences can be easily appreciated by considering the two extreme beam and disc-like distributions corresponding respectively to a=0 and rx=oo. The plasma corresponding to a=0 (i.e T n » T j ) has only one degree of freedom for exchange energy, whereas the plasma corresponding to a - o o has two degrees of freedom available for energy exchange, and therefore relaxes faster. Connected with this work are the experimental results of Hyatt, Driscoll and MalmbergO6) (HDM) concerning the relaxation of anisotropic temperature distributions of electrons in a magnetized pure electron plasma. The data from two of their figures (figures 2 and 3) are shown as the circles and triangles in figures III.6 and 7 and compared with the theoretical results from this study. Figure HT.6 compares the HDM experimental data obtained from their figure 2 for the relaxation of T,(t) and Tx(t) with results from this study. Unfortunately HDM make no reference to the number density used for their data (except for an earlier mention of electron densities in the range 3xl0 6<n£3.5xl0 7cm" 3 ). Note that the result n=5xl06cm~3 lies very close to the HDM data for T„(t), but lies significandy above their data for T±(t). The relaxation rates given by HDM for their data in figure LTI.6 are v„ =79 sec-1 and vx= 121 sec-1 corresponding to the reciprocal of the relaxation time for T„ and T± respectively. Theoretically, the two temperature components should possess the same relaxation rates, 76 Figure m.5 Relaxation of the anisotropic temperature parameters T|((t')/T (A and B). 77 Figure IH.5 (continued) Relaxation of the anisotropic temperature parameter ct(t'), in C and e(t') in D for initial temperature anisotropics as labelled. 78 Figure LTI.5 (continued') Relaxation of the anisotropic temperature parameter a(t'), in E and e(t') in F for initial temperature anisotropics as labelled. 79 time (ms) Figure ITI.6 Theoretical and experimental results for relaxation of T | ((t') and T ± ( t ' ) . The experimental points are the H D M data from their figure 2. The number densities, n, for the theoretical curves are as labelled. 80 Figures III.7 Theoretical and experimental rates of relaxation for various electron plasma densities and temperatures. The experimantal points are the HDM points from their figure 3. The triangles refer to a(0)<l in their experiment and the circles to a(0)>l. The upper line is the rate determined from equation (III.3.5) for an electron plasma, the lower line is the result from the present study for a(0) = 8. 81 the discrepancy arising in the HDM experiment because v„ is "some measure of the energetic tail of the parallel distribution". Since the Coulomb collision cross section decreases as (energy)-2 the collisional relaxation rate for v„ is expected to be less than that for v ± as stated by HDM and as verified by their data. These author's consequently believe that their data for v x is a much better estimate than that for vM. It is interesting to note that the theoretical number density corresponding to v ± (for which a(0) = 0.23) is n = 4.8xl06cm-3. This value is smaller than that expected from the theoretical curves in figure LTJ.6. It is consistent with the data points shown in figure LLI.7 and confirms that if the results of the HDM data from figure ILT.6 are indicated in figure III.7, then n « 5xl06cm-3 is very close to their experimental value. We note in passing that the value for the relaxation rate used by HDM in figure III.6 would be represented by a triangle in figure LTI.7 at virtually the same value of n as predicted from our results. The HDM work also suggests that the relaxation rate for their anisotropic pure electron plasma is comparable to the Ichimaru and Rosenbluth (LR) rate(10). This can be deduced from a comparison of the HDM data points in figure IJJ.7 with the upper line, which corresponds to the LR rate. The lower line (rate) in the figure corresponds to the relaxation rate calculated from the present theoretical work for rx(0) = 8. Unfortunately once again, the anisotropy parameters used by HDM for their data in figure ITI.7 are not known (except for the data from figure HI.6). Since the theoretical results from this work and the work by Jorna and Wood(11) suggest fundamental differences in the rate of relaxation when cc(0) < 1 than when a(0) > 1 it is suggested that experimental results confirning this, or otherwise, would be extremely important in assessing the validity of the Fokker-Planck operator. This is also of importance in connection with recent theoretical work which suggests that the use of the Fokker-Planck collision operator for 1/r type interaction 82 potentials is unsound, and that the full Boltzmann collision operator should be employed^76). It is interesting to note the HDM statement that the least squares line of best fit for their data in figure ffl.7 yields a line 5% below the IR line. Since the IR line corresponds to cc(0)»0.2 from the results of this work, it is not surprising that the HDM line of best fit lies below the IR line. The appearance of the HDM data in figure LTI.7 suggests that their rates for compression, a(0)>l, lie below their rates for expansion, a(0) >1, in accordance with that expected from this work. There are several other important points to make in connection with the HDM and IR results. First, an excellent approximation for the IR rate can be derived by considering the expansion of equation (UI.2.11) to first order in the temperature anisotropy (i.e n = 0, / = 2), so that for this latter moment we have j dv X|/02) | = Jdv V<2)C[f, f] - u 0 2 ( t ) (M™ + N<2>) (III.3.2) The eigenvalue and temperature coefficient are easily determined to be: X « - A « - ^ and tf> = V 2 (0 ) = 2 * $ = £ . (in.3.3) yielding to first order, the approximation M2 o 1+ 2 r a ^ - 1 1 e-*?l' = 1 + 2 r a ^ - 1 1 (III 3 4) T La(0) + 2-l La(t') + 2J ' (2) Consequently the rate constant, v=XI is the reiprocal of the relaxation time (in units of 83 sec-1) can be written as v = l / Z n e 4 - J n A ( m 3 5 ) 5 V m (kT) 3 / 2 which differs from the LR expression only by the replacement of T with the effective temperature, T e f f , defined by LR as i = J i f %tSLz!Ll _ (m.3.6) The integral can be performed analytically yielding, on multiplication by T* 12, 5 ra(t)+21 3' 2f 2a(t) + l wf/pfrPi +fih)) - 1 La(t)-1J l2ra(t) - 11 l o g ^ a w + 2 V a(t)-1 J 3/2_ • 2 / 3 L a ( 0 - U »-2[a(t)-l] " V T W * w / 2Va(t) <TJ ' 1 ^ _ra(t)-f2-i 3/ 2f i-2a(t) sm-i/rr^ tj - 2 r ^ T l 2/T L1 -a(t)J I 2[1 -a(t)]  i K ) 2Vl-a(t)J (LU.3.7) where the uppermost expression in the braces applies for a(t)£ 1, and the lower for a(t)< 1 (see equation (LTJ.1.1) for the definition of a(t)). HDM refer to (T/T e f f ) 3 / 2 as h(Tj/T) and a plot of this function is provided in figures LTI.8A and B. It deviates marginally (=5%) from unity in the range of the HDM experiment, but significantly when T x /T-> 0 and for large values of T ± / T , as shown. The important point to be made here is that the use of equation (LTJ.3.5) in place of the LR expression is justified for the HDM experimental range. 84 1 .3 I — 1 — i — ' — i — ' — i — < — i — • — i — • — i — < — r n q I • i i i i i i i i i . i i i i i i _ J 1 2 3 4 5 6 7 8 9 10 a Figure III.8 The function h (T ± /T) = (T/T e f f ) given in equation (ILI.6.8) plotted versus a . A is for a < 1, while B is for a > l , corresponding to the upper and lower equations in equation (III.6.8). 85 Although we have specifically ignored the effects of collisional relaxation in the presence of electric and magnetic fields, it is a relatively straightforward matter to include the effects of a strong homogeneous magnetic fiektf89). Indeed the agreement between the HDM experiment, where the electron plasma is confined in a magnetic bottle, and the IR theory depends upon the Montgomery, Joyce and Turner (MJT) approximation^ 89) that the charged particle's gyro-radius can function as an effective cut-off for the long range Coulomb potential. The physical conditions validating this assumption require the gyro-radius (or Larmor radius) to be significantly less than the Debye length, but significantly greater than the distance of closest approach. MJT state that these criteria should not be too difficult to satisfy for electrons in many plasmas although they are more stringent for ions. Note that the Larmor radius is given by TL=V±/(Q0, where vx is the perpendicular velocity component (it depends on the ratio W i n for an average plasma particle) and co0=qB / m, is the cyclotron frequency. Also of interest is the experimental work by Church(18), who studied the collisional relaxation between the axial and radial degrees of freedom for protons in a radio-frequency quadrupole trap. His results yield agreement to within 25% of the Spitzer relaxation times. Strictly speaking, however, the IR relaxation time or the reciprocal of the rate appearing in equation (JU.3.5) above should be used in place of the Spitzer relaxation time (which, for self colliding particles, is lower by a factor of 1.2(16)). This is because the Spitzer relaxation time applies to isotropic ion velocity distribution functions, whereas equation (JJI.3.5) applies to anisotropic ion velocity distribution functions. The consequences of this modification introduce a greater disparity between Church's relaxation times and the one obtained from equation (LTL3.5), the disparity now being larger than Church's stated 25%. It is also important to note that the trap technique used in this experiment does not require the presence of a magnetic field. 86 Another experimental study dealing with the effect of anisotropic relaxation is the work of Amemiya, Oyama and Sakamoto(17)who measured the rate of temperature relaxation in an ECR (electron cyclotron resonance plasma). These authors have also used the rate equation for isotropic relaxation given by Spitzer, but find difficulty explaining the relaxation in terms of Coulomb collisions only. Of equal significance in their experiment are ion-neutral collisions, accounting for 16 to 66% of the relaxation. These authors consider the possibility of a very large contribution to the relaxation rate from drift instability, arising due to the non-linear nature of their ECR plasmaC94). These authors consider this relaxation mechanism to be more significant in their experiment than that of Coulomb collisions. They also postulate the existence of several other energy relaxation mechanisms such as particle reflection at the mirror points (where there are significant magnetic field inhomogeneities). No relaxation times are given in this study. Apart from this present theoretical study one should consider also the work of Jorna and Wood(11) (JW) and the work of Livi and Marsch(3°) (LM), who solved the full non-linear Fokker-Planck equation for an initial bi-Maxwellian ion velocity distribution function. A problem arises in relating relaxation times from the full non-linear operator to relaxation times from the work presented here in that the results from the non-linear work have, at best a qualitative interpretation only. For example, both JW and LM only display plots of the temperature relaxation. These plots presumably are based on some kind of least squares curve fitting between various points and LM actually quote their results with 10% errors. One must also be careful with the many different definitions of the collision time scale. This study employs the definition of the Coulomb self collision time in equation (JJI.2.22), whilst LM employ a definition differing from ours by a factor of 2/3/7. Figure 5 of their work suggests that the relaxation time for their double beam distribution with a(0)« 2.5 is TL M=2.7±0.3, which translates to t a » l in our dimensionless time units 87 (xa refers to the 1/e collisional relaxation time for cc(t)). This result is in excellent agreement with the result from the present theoretical work (1.018) shown in Table LTJ.3. Another result from the LM work is that from the bi-Maxwellian velocity distribution function, for which they chose e(0) = T1(0)/T||(0) = 4. The relaxation time from their figure 11 is estimated to be 1.510.15 which translates into our units as xe«0.56±0.05. The result from this study, xe=0.51 is just within range of this LM value. Additionally, LM quote an 'average e-folding time' for their bi-Maxwellian study. With the assumtion that this is the l/2e relaxation time their x' L M=3.4±0.3 corresponds to x' e«1.3±0.1 in our units, comparable to the value from our Fokker-Planck study of x'£=1.38. The Jorna and Wood (JW) studyO1) features two results that enable comparison with this study. These results correspond to curves c and d of their figure 1, for which e(0)=10 and ct(0)=10 respectively. With the assumption that the time scale used by these authors is based on electron collisions (they never state what plasma system they are considering!) their time scale differs from ours by 3%. The results from their figure 1 as well as the Livi and Marsch results are shown in the Table LTJ.3. This table also includes many other results from this study, enabling comparison with future experimental and theoretical work. In particular we note that the JW results are about 40% smaller than our values. It is difficult to assess where this discrepancy might arise since the JW initial anisotropies are very large. The plot of the relaxation time for T,|(t')/T versus the basis set size in figure LTL2 suggested that there might be an error associated with cc(0) =10, but not of the order of 40%. There exists the possibility that the linearized operator is a poor approximation to the full non-linear operator for these large anisotropy parameters, but detailed results from non-linear Fokker-Planck studies are missing. The final plot appearing in figure ILT.9 is to draw attention to the fact that although the initial ion velocity distribution function chosen for this problem is bi-Maxwellian, at no 88 Table mt3  Comparison of Relaxation times* ct(0) 0 1.152 1.516 0.0003 1/10 1.170 1.486 0.232 0.14 (JW) 1/8 1.174 1.480 0.284 1/6 1.180 1.468 0.365 1/4 1.193 1.447 0.511 0.5610.05 (LM) 1/2 1.228 1.387 0.849 2 1.342 1.096 1.747 2.5 1.361 1.018 1.0±0.1 (LM) 1.894 4 1.399 0.833 2.185 6 1.426 0.666 2.406 8 1.443 0.552 2.540 10 1.45 0.47 0.27 (JW) 2.63 In units of x, where xjt, x a and x£ denote the relaxation times for T||(t')/T, ct(t') and e(t') respectively. 89 t' Figure IJI.9 Ratio of the expansion coefficients, w r 2(t')/ur 2(t') for the initial bi-Maxwellian ion VDF with a(0) = 2. The former coefficient, w r 2(t'), is the coefficient for the ion VDF always being bi-Maxwellian, and u 2(t') is calculated from the Fokker Planck equation, and is obtained from equations (III.2.24 and 32). 90 time other than at t=0 does it remain bi-Maxwellian. This is because the energetic 'tail' of the distribution function relaxes slower than the more energetically moderate parts of the velocity distribution function. The expansion coefficients assuming that the temperature relaxation is always bi-Maxwellian are denoted wr2(t') and the figure shows wr 2(t')/ur 2(t') versus r, where ur^ (t') is obtained from equations (JJI.2.24 and 32). III.4 Conclusions The present study based on the linearized Fokker-Planck equation for a single component anisotropic plasma has yielded excellent agreement with experiment The results from this study suggest where attention should be focused in anisotropic charged particle relaxation studies. In particular experimental and theoretical studies with the non-linear Fokker-Planck operator should be conducted over a larger range of temperature anisotropics than conducted so far. This would enable a quantitative assessment concerning how good an approximation the linearized Fokker-Planck operator is to the non-linear Fokker-Planck operator and, as well, how well both operators compare with experiment. It would be interesting to see whether or not a more thorough experimental study could yield differences in the rates of relaxation for a> 1 and ct< 1. The numerical code developed for this problem has application to the study of many other plasma systems such as the two component plasma. This problem, for isotropic ion velocity distribution functions, would be of interest in connection with the work on the linear isotropic test particle problem, because in the limit of one species having a much larger number density than the other, the two component plasma relaxation problem is physically the same as the linear test particle problem 91 Chapter IV  Summary The aim of this thesis has been to account for the collisional relaxation in two kinds of plasma systems. Both plasma systems were modeled as being spatially homogeneous and devoid of external electric and magnetic fields. The collisional relaxation was described in terms of a linear Fokker-Planck equation and the distribution function for the plasma species was expanded in the eigenfunctions of the Fokker-Planck operator. These eigenfunctions were represented by two different methods in this thesis; the first involved their representation at the quadrature points of a quadrature rule based on speed polynomials, Bn(x), orthogonal with respect to the weight function, w(x)=x e x , whilst the second method involved their representation in terms of associated Laguerre polynomials. Both methods yielded converged solutions for the relaxation time for isotropic and anisotropic velocity distributions, despite the representation of continuum eigenfunctions in a finite polynomial basis. This suggests that a large number of kinetic theory problems may be represented in terms of finite polynomial basis functions, although the spectrum of the collision operator may have a continuous portion. The problems studied in this thesis were compared with results from other Fokker-Planck studies and experimental work. Excellent agreement was found suggesting that the use of the Fokker-Planck collision operator is well founded for the present application. The readers attention is drawn to recent literature suggesting the contraryC76) and advocating the use of the full Boltzmann collision operator for a correct description of plasma collisional relaxation. The Boltzmann collision operator would predict significantly larger relaxation rates for both plasma systems considered in this thesis based on Shoub's fmdingsC76). The present theoretical work parallels closely with numerous other theoretical and experimental studies, particularly in connection with chapter HI where the linearized 92 Fokker-Planck equation is seen to be in excellent agreement with the full non-linear Fokker-Planck equationC11 )^ and the experimental work of Hyatt, Driscoll and Malmberg (HDM)(16) provided the departures from equilibrium are not too large. This statement requires further qualification and quantification in the form of more extensive studies on the non-linear Fokker-Planck operator. Comparison with the HDM results could be.made more complete if experimentalists turn their attention to looking for fundamental differences in the rates of relaxation for plasmas with T„ > T and T„<T. In connection with the findings of HDM the work of Shoub using the Boltzmann collision operatorC76) predicts significantly larger rates of relaxation. Consequently, further studies involving the Boltzmann collision operator would be of interest for the specific range of physical parameters considered in the HDM experiment Finally there are many aspects of the present work which can be extended. For example, the numerical code for the generation of the M and N matrix elements of chapter in based on appendix 1, can easily be applied to the problem of the two component or multispecies plasma relaxation problem. This problem can be used to study the importance of plasma species number density ratios on the rate of relaxation. Although, like the problem of chapter LTI, this problem is non-linear but can be linearized when one plasma species number density is very much larger than the other's. The problem reduces to the one considered in chapter H where the differential form of the Fokker-Planck operator, L^ , acts on the velocity distribution function of the rninor constituent A comparison of the two numerical methods used in the present work might also be of interest with respect to convergence of the energy relaxation. 9 3 Appendix 1 We begin this appendix by discussing methods for the computation of the M matrix elements. These are matrix elements of the opertor L Q as explained in chapter L T I . 2 (see equation ( U I . 2 . 1 2 ) ) . A formula due to Wong that enables the matrix elements to be computed via a quadrature rule is introduced and the results for several different quadrature rules are tabulated. Next, recursion relations for the M matrix elements are introduced. A numerical method for their computation is outlined and several analytical expressions for these matrix elements are given. Finally, the N matrix elements are discussed. These are matrix elements of the operator Lj (see chapter H T . 2 equation ( H T . 2 . 1 3 ) ) is derived via a generating function approach extending the work of Wong(67). An explicit expression can be written for these matrix elements. This appendix is general in the sense that multiple species collisions can be considered (i.e one can consider sums over the species, b), which gives rise to a dependence of the matrix elements on the mass ratio, p = - / m t / m a - The collision term given by equation (1.4.lb) is therefore of relevance to the following development The M Matrix Elements These matrix elements are defined by equation ( L U . 3 . 1 2 ) and are given by Wong(67) in the form, oo = f dxe^x 2 7 - 1 IViSSllV; + 2 / x 3 d G l(L' ML N + L M L ; ) m-1 \ B 2 / + 1 J L dx m n Hx 7 o + ( / 2 x d ^ + / ( / + i ) G ' ) L M L N ] , (Al.l) 94 where the integration over the Legendre polynomials has been performed and, Jill 3/nT (kT) T , = \— a ' ' (A1.2) * 472TeaeJnblnA is the two particle momentum deflection time. All other quantities in equation (A 1.1) are defined in chapter HI with 'a' and 'b' referring to the two plasma species and G(xb) is defined as, _ 2 GC^) = ( \ + l/2xb) erf (x^ ) + ^  , (A1.3) where xb= px and G' denotes differentiation of G with respect to xb. Also in equation / + 1/2 o 0 (A 1.1), L n=L„ (x ) and L' denotes differentiation of L with respect to x . (Note that the x5 term appearing in equation (A 1.1) appears incorrectly as x3 in Wong's paper. Equation (A 1.1) also makes apparent the symmetry, M ^ n = M ^ m for all p. Equation (A 1.1) has a form that enables numerical evaluation by several different quadrature rules. Such quadrature rules were defined in equation (1.3.1) of this thesis. For 2 _x2 example we may choosing w(x) = x e * for a quadrature rule based on speed points and weights, or w(y) = J~y e"y for a quadrature rule based on associated Laguerre points and weights (/ = 0) or w(y) = e_> for a quadrature rule based on Laguerre points and weights (y = x2). The performance of all three of these quadrature rules are tabulated in table Al , which suggests that the speed and associated Laguerre quadrature rules are the best for evaluating integrals of the kind in equation (A 1.1). They enable a rapid computation of any 95 Table Al Convergence of M Matrix elements using speed, associated Laguerre (/ =0) and Laguerre Quadratures. N Speed Ass. Laguerre Laguerre. M<2) 2 0.241 -0.191 0.521 6 0.410 0.540 0.367 10 0.0615 0.0825 -0.00380 15 0.0626 0.0626 0.0519 20 0.0589 25 0.0607 30 0.0615 40 0.0621 ^9,9 2 2.132 2.531 2.988 6 2.171 2.279 2.273 10 2.729 3.462 3.289 15 2.681 2.720 2.565 20 2.685 2.685 2.659 25 2.674 96 particular matrix element without the need for recursive schemes such as required by Wong's generating function technique. The generating function for the M matrix elements can be derived from the recipe outlined in equations (1.3.3-5) of chapter I and applied to equation (Al.l). Wong defines^ 67) the generating function for the M matrix elements as, oo M(%n) = X^C^V ' (A1.4) m,n = 0 and employs equation (1.3.4) and its first derivative with respect to xz to derive his system of equations (26) through (31). Note that the matrix element M ^ n is obtainable from the generating function expansion in (A 1.4) by picking out the coefficient of £ m T | n . Numerical methods based on (A 1.4) yield these coefficients by using power series expansions as discussed below. The recursive scheme for the evaluation of Wong's matrix elements is given in his equations (26) through (31) which we rewrite below as M«> = - i Htg®. {a<0"B, • . « > B , + I + « « > B , + 2 • a«>B, + , } , (A..5) ab V 7 1 where J3 / + 1 = [l-^T| + p 2 ( l - 0(l-Tl ) ] C / + 1 - Ijlic, (A1.6) and 97 C - 2 1 C< i 1 / + 1 2/-1 i-$n (i-^)2(i-n)2(i-^T|) l p [ l - ^ + P 2 ( l - ^ ) ( l - T l ) ] / + 1/2 (A 1.7) The above equations require C 0 = 0, and for the 'a' coefficients appearing in equation (A 1.5) Wong writes: (0 _ _ 3/(/-l) f A 1 8 , a0 " 2(2/-1) ' ( A 1 ' * a ) a; = * l 1 J / + 2(/+l)(c>n)-(5/ + 4)cX] , (A1.8b) i ® i 2(2/-1) 4° = f[- / (^ + T l ) -4^T 1 + (3/ + 4 ) ^ T | ( ^ + T 1)-4(/+ 1)^2T,2] , (A1.8c) 4° = | (2 / + 3 ^ T i [ l - ( ^ + Ti) + ^ T i ( ^ + n)-^ 2n 2] . (A1.8d) These equations completely describe a recursive evaluation of the generating function for the M matrix elements. The crux of the computation being equation (A 1.7) where the denominator must be expanded in a (double) power series. It is the second term on the right hand side of this equation which poses the most problems, because the product (or convolution) of two double power series must be computed for each value of /. This can be accomplish by employing the power series expansions given below: (1-cx)' = £ ( £ ) ( - c x ) k k=0 (A1.9a) 98 where = — ^ ^  is the combinatorial function. Equation (A1.9a) is the binomial expansion and its reciprocal is the infinite power series, 1 7 = X ( / + £-1)<cx>k . CA1.%) (1-cx)* kT and hence i = y n/+k-i/2) ( .k A 9 . ( i - c x ) ' - 1 / 2 h r(/-V2)ki ( c x ) ( A L 9 c ) for / integral. These expansions constitute the means for obtaining the matrix elements from the generating function. From equations (A 1.7 and 9b) it then follows that, ^•(.-tfa-U-w = • ( A U O a ) where A (£, rj) = A (T|, %) requires aij = a .^ For i £ j we can show that . q + i)q+2K3H + 3) ( A U O b ) *0 O which is an integer. Similarly we can show that 99 D(/)<^> = ~r H T7TT72 = 2 1 / + 1/2 (AM la) where D ( />(£,TJ ) = D ( / ) (n,£) requires d?? = d??. For i£ j we can show that, d( / ) = i ( V - i V ( o1 V ^ rg + k + i-n/2) / p4 \ k i.j j ! r ( / + l / 2 ) V ~ V p 2 ^ / g (k + i-j)! Vp 4-!^ (Al.llb) Consequently, the computation of the second term in equation (A 1.7) can be performed by writing, oo oo P(P + 1) i,j=U r.s=0 oo 1 oo s n , n a * i+i/i Z X Z X a i , d j - i , s - r ^ V (A1.12) P(P + 1) j=0 i=0 s=0 r=0 These last equations are at the heart of the numerical method but their usefulness in computing analytical expressions by hand for the M matrix elements is dubious. Instead it seems easier to employ equation (Al.l) using the integrals below, oo a(k,p) = Jdxe" xV k" 3 erf(px) (A1.13a) o = P g rdc- j -^Hk- i ) ! 2 / 7 p0 (p 2+l) f c-J-3 / 2(k-j-2)! and 100 P ( k i p ) = 2$= [dx xa<k-i>c-*»+i>* = pr(k-l/2) (A1.14) / 7 J (p2 + i ) k - 1 / 2 /F Hence we can show that /=0 yields, M ? > - — a n d M<°> = ~ 9 P 3 , / 2 A (A1.15) 1 , 1 (p 2 + l ) 3 / 2 U * 2(p2+l)5/2 *ab when /= 1, = ~P ,„ A and M£> = ~ 3 p ' A (A1.16) V ° 2(p2 +1)1 / 2 t A °- 1 4(p2+l)3/2 t a b and when /=2, M<2> = ~ 3P ( ? 3P tA> — (A1.17) 0 > 0 10(p 2+l) 3 / 2 * a b where the /= 1 terms differ from those listed in Hinton's article^ 5) by a factor of 2 and in the absence of the collision frequency, njx^ (see equations (264-74) of Hinton(5)). Hinton(5) has also made use of the equality p =Vmb/ma = v a/v b in his equations (270-74). As a conclusion to this section we note thai when p = 1, M*2' = -^ L-^ s. which is used to derive the expression for the collisional relaxation rate for the single component plasma studied in chapter JH (see equation (7H.3.5)). The N Matrix Elements We now consider the generating function for the N matrix elements and show that the expression given by Wong can be reduced further to an explicit analytical form for any 101 matrix element These matrix elements are matrix elements of the operator Lj in the Burnett function basis (see chapter HT.2 equation (HI.2 .13)). Wong's equations (44) through (46)(67) describe the generating function for these matrix elements and are written, N ( / ) r f m - n« H/+1/2) 3 { I2 , /(2/-1) f t . , r\\ . 4/2-1 p_l where (A1.18) D = p(l-c,) + 1^1 (A1.19) By employing equations (A1.9a-c) of this appendix we can, for example, show that the coefficient of £ m r | n in D"n'+ 1 / 2 is, Coefficient = , p/ + T"T i / . ^nV^'f <A1'20> (p2 + l) 7 r ( / - 1/2) m! n! with simple alterations leading to the coefficients of £ m T | n in D ~ / - 1 / 2 and D - / - 3 / 2 . Collecting the various coefficients in (A1.18) enables the coefficient of £ m T i n in N ( / ) to be determined. Carrying through the necessary algebra then yields a simple and compact expression for the N matrix elements, given below as, N (0 = 3nap f + 2 m - 1 r ( / + m + n-l/2)(/ + 2m) (7 + 2n) ( A 1 2 1 ) m , n 2 (21 + l)xAhfit m! n! (p2 + 1 ) l + m + n " 1 / 2 which, unlike Wong's equations (44) through (46) requires no recursive scheme at all. For 102 numerical purposes we can write the above equation in the form = ± p 2 m Sm> = p 2 ( m - n ) , (A1.22) m,n x m.n r n.m ' ab which makes explicit the appearance of a symmetrical part, S^n = S®m in the N matrix elements. We also note that equation (A1.22) satisfies the symmetry relation NmlntCJ = , (A1.23) where the bracketed collision term, C^, denotes a on b collisions and vice-versa for C^. That is, with respect to a on b collisions the quantity appearing on the left hand side of equation (A 1.23) is unchanged on interchanging a with b and m with n. This is also equivalent to N m > ] = O"1] • (A 1.24) Conservation Properties of a One Component Plasma The case of a one component plasma introduces many simplifications into the generating equations (especially for the M matrix elements). Apart from the various mathematical simplications that arise, the physical constraints for a non-reactive, isolated plasma system in the absence of electric and magnetic fields, require that particle number density, momentum and energy be conserved. Consideration of the moments of these conserved quantities (i.e \\f®\ Y j 0 ) and ) with the ion VDF in the Burnett function 103 basis (see equations (ILT.2.9) and (LTJ.2.15)) yields the following additional relations for a one component plasma, < - - " . ? . - » . - < > - < : • (A1.25) An important connection with classical kinetic theory is gained here where linear combinations of these conserved quantities are referred to as 'summational invariants^59'60). The conservation properties of the plasma considered here thereby account for the presence of 5 zero eigenvalues (1 for particle conservation, 3 for momentum conservation and 1 for energy conservation). In a similar fashion the conservation properties of multicomponent plasmas can also be established. For example in a two component plasma undergoing inelastic collisional relaxation, particle number density is conserved with respect to each component (2 zero eigenvalues), system energy is conserved (1 zero eigenvalue) and system momentum is conserved (3 zero eigenvalues) to give a total of 6 zero eigenvalues for this two component plasma. For a non-reactive k component plasma relaxing via elastic Coulomb collisions we can therefore expect k+4 zero eigenvalues. 104 Appendix 2 An expression for the Initial Expansion Coefficients, 2(0) In order to solve the linearized Fokker-Planck equation for the time dependence of the anisotropic temperature, Tu(t), the initial expansion coefficients, ur^ (0), are required. Although a quadrature rule is quite satisfactory for the evaluation of the double integral appearing in equation (JJI.2.39), a little more insight may be gained by evaluating the integrals analytically. A method for acheiving this is outlined in this appendix. First the integral over \i can be performed by noting that, l LJx) = 1 Jdn(3n 2-l)e-^ 2 - l (A2.1) where we have written a=a(0) and, Y=y(x) = xyp^o7) , P=-£3~? (A2'2) which reduces ur 2(0) to, 0 105 where, - - S C t t * ! ) 3 * (A2.4) 3V3ap(l-a) The expression in equation (A2.3) may be evaluated by the generating function technique previously introduced in equations (1.3.3-5) of chapter 1. We introduce, V . i £ i m f ^ (A2.5) so that on substituting equations (A2.3) and (1.3.4) into equation (A2.5) we can write, V r9^?3* m CA2.6) 4(a + 2 H l + a $ H l -2a$) The coefficient of £ r in equation (A2.6) is then related to ur 2(0) via equation (A2.5). Expanding the above equation as a power series in % yields for u^O), which is equivalent to the expression in equation (in.2.40). 106 REFERENCES 1. L.Spitzer. Phvsics of Fully Ionized Gases (John Wiley, New York, 1962). 2. W.M.MacDonald, M.N.Rosenbluth and W.Chuck, Phys.Rev., 1QL 350 (1957). 3. M.N.Rosenbluth, W.M.MacDonald and DUudd, Phys.Rev., 107, 1 (1957). 4. D.GMontgomery and D.A.Tidman, Plasma Kinetic Theory (McGraw-Hill, New York, 1964). 5. F.L.Hinton in Collisional Transport in Plasmas in Basic Plasma Physics I, ed. A.A.Galeev and R.N Sudan (North-Holland, Amsterdam, 1983). 6. T.H.Kho, Phys.Rev, A22, 666 (1985). 7. LP.Shkarofsky, T.WJohnston and M.P.Bachinski, The Particle Kinetics of  Plasmas (Addison-Wesley, Reading, Mass., 1966). 8. M.Mitchener and CJi.Kruger. Partially Ionized Gases (John Wiley, New York, 1973). 9. H.C.Kranzer, Phys.Fluids, 3J), 1340 (1987). 10. S.Ichimaru and M.N.Rosenbluth, Phys.Fluids, 13, 2778 (1970). 11. SJorna and L.Wood, Phys.Rev.A., 26, 397 (1987). 12. M.Eyni and A.S.Kaufman, J.Phys..A., iQ, L83 (1977) 13. R.J.Burke and R.F.Post, PhysJluids, 17_, 1422 (1974). 14. D.M.Cox, H.H.Brown, LKlavan and B.Bederson, Phys.Rev.A, 10, 1409 (1974) 15. M.Katsch and K.Wiesemann, Plasma Phys, 22,627 (1980). 16. A-WJiyatt, CRDriscoll and J.H.Malmberg, Phys.Rev.Lett., 5_9_, 2975 (1987). 17. RAmemiya, ROyama, Y.Sakamoto, J.Phys.Soc.(Japan), 5J>, 2401 (1987). 18. D.A.Church, Phys.Rev.A., 3J, 277 (1988). 19. RJ3Jiazeltine and F.LJiinton, Phys Jluids, 1& 1883 (1973). 20. H.W.Drawin in Reactions Under Plasma Conditions, ed. by M.Venugopalan (John Wiley, New York, 1971), Ch.4. 21. D.CKhandekar and D.C.Sahni, Phys.Fluids, 21 1464 (1980). 22. HTsuji, M.Katsurai, T.Sekiguchi and N.Nakano, NucLFusion, 16, 287 (1976). 107 23. N.Mizuno, H.Irie and M.Sato, J.Phys.Soc.(Japan), 54, 1830 (1985). 24. Yu.N.Dnestrovsky and D.P.Kostomarov, NucUusion, i i , 141 (1971). 25. E.B.Corman, W.E.Loewe, G.E.Cooper and A.M.Winslow, NucLFusion, 15_, 373 (1975). 26. P.A.Haldy and J.Ligou, NucUusion, 11, 6 (1977). 27. L.W.Klein, K.W.Ogilvie and L.F.Burlaga, J.Geophys.Res.^ X 7389 (1985). 28. R.Hernandez and E.Marsch, J.Geophys.Res., 2Q, 11062 (1985). 29. E.Marsch and S.Livi, Ann.Geophysica, 3_, 545 (1985). 30. S.Livi and E.Marsch, Ann.Geophysica, 4, 333 (1986) 31. E.Marsch and S.Livi, Phys.Fluids, 28., 1379 (1985). 32. J.CBrandt, Introduction to the Solar Wind. (W.H.Freeman, 1970). 33. A.R.Barakat and R.W.Schunk, J.Geophys.Res., £9_, 9771 (1984). 34. A.R.Barakat and R.W.Schunk, Plasma Phys., 24, 389 (1982). 35. P.M.Banks and GJ.Lewak, Phys.Fluids, LL 804 (1968). 36. T.E.Holzer, J.A.Fedder & P.M.Banks, J.Geophys.Res., 26, 2453 (1971). 37. Although there seems to be no specific reference to the importance of Coulomb collisions in inductively coupled plasmas we may conclude that their effects might be important via Stark broadening of excited state ions. Also the non-equilibrium distributions of ions in these plasmas would cool and isotropize via Coulomb collisions. See for example: J.W.Mills and CMJIieftje, Spectochim. Acta, 29JB, 859 (1984) and, 38. I.J.M.M.Raaijmakers, P.W.J.M.Bowmans, B. van der Sijde and D.C.Schram, Spectochim. Acta, 2813,697 (1983) 39. M.W.Blades, B.L.Caughlin, ZiLWalker and L.L.Burton, Prog, analyt. Spectrosc, 10_, 55 (1987) 40. T-C Lin Wang and A.G.Marshall, IntJ.Mass.Spectrom.Ion Processes 6j£, 287 (1986). 41. J.B.Jeffries, SJi.Barlow and G.H.Dunn, IntJ.Mass.Sprectom.Ion Processes, 54,169 (1983). 42. T.J.Francl, M.G.Sherman, R.L.Hunter, MJ.Locke, W.D.Bowers and R.T.McIver, IntJ.Mass.Sprectrom.Ion Processes, 54, 189 (1983). 108 43. M.B.Comisarow, private communications. 44. L.Spitzer. Physical Processes of the Interstellar Medium (John Wiley, New York, 1978), Ch. 2. 45. J.M.Retterer, Astron.J., £5_, 249 (1980) 46. J.M.Retterer, Astron.J., M , 370 (1979). 47. M.R.Hoare, Adv.Chem.Phys., 2Q, 135 (1971). 48. M.R.Hoare and CJi.KapUnsky, J.ChemPhys., 52, 3336 (1970). 49. B.Shizgal and R.Blackmore, CanJ.Phys., 42, 97 (1984). 50. R.JJlner and I.Kuscer, J.StatPhys., 20_, 303 (1979). 51. N.G.van Kampen, Stochastic Processes in Physics and Chemistry. (North-Holland publishing, Third edn.) 125 (1984). 52. N. Corngold, P. Michael and W. Wollman, Nuc. Sci. Eng. 15, 13 (1963). 53. B. Shizgal, J. Chem. Phys. 24, 1401 (1981). 54. R.Blackmore and B.Shizgal, Phys.Rev.A, 3_L 1855 (1985). 55. B.Shizgal and R.Blackmore, LComputatPhys., 55,313 (1984). 56. B.Shizgal, J.Computat.Phys, 41, 309 (1981). 57. D.Burnett, Proc.Lond.Math.Soc., 22, 385 (1935). 58. D.Burnett, Proc.Lond.Math.Soc., 40_, 382 (1935). 59. S.Chapman and T.G.Cowling, The Mathematical Theory of Non-uniform Gases (Cambridge Univ. Press, Third ed., London 1970) 60. J.Ferziger and H.Kaper, Mathematical Theory of Transport Processes in Gases (North-Holland, Amsterdam, 1970). 61. C.S.Wang-Chang, G.E.Uhlenbeck and J.de Boer, in: Studies in Statistical  Mechanics, eds. J.de Boer and G J£.Uhlenbeck, North Holland, Amsterdam, 1970). 62. L.A.Viehland and E.A.Mason, Ann.Phys., 91., 499 (1975). 63. K.Kumar, AustJ.Phys., 2Q, 205 (1967). 109 64. M.J.Lindenfeld and B.Shizgal, CheraPhys., 41, 81 (1979). 65. B.Shizgal and J.M.Fitzpatrick, Chem.Phys., & 54 (1974). 66. G.W.Ford, Phys.Fluids, H , 515 (1968). 67. S.K.Wong, Phys.Fluids, 2&, 1695 (1985). 68. G.Arfken, Mathematical Methods for Physicists (Academic Press, New York, second ed.,1970) 69. F.F.Chen, Plasma Physics and Controlled Fusion. (Plenum Press, New York and London, second ed., 1985). 70. D.R.Nicholson, Introduction to Plasma Theory. (John Wiley & Sons, 1983). 71. S.Ichimaru, Basic Principles of Plasma Phvsics. (Benjamin Inc., 1973). 72. L.Landau, Phys.Z.Sowjetunion, 1Q_, 154.(1936). 73. G.Kamelander, Nuc.Sci.Eng.8Ji, 355 (1984). 74. G.CPomraning, Nuc.Sci. Eng.85, 116 (1983). 75. M.Caro and J.Ligou, Nuc.Sci.Eng, £3, 242 (1983). 76. E.C.Shoub, Phys.Fluids 2Q, 1430 (1987). 77. B.Shizgal and D.R.A.McMahon, Phys.Rev.A, 21,1894 (1985). 78. D.R.A.McMahon and B.Shizgal, Phys.Rev.A, 22, 3669 (1985). 79. J.S.Chang and G.Cooper, J.Computat.Phys. & 1 (1970). 80. J.Killeen and K.D.Marx, in Methods in Computational Physics, 16_, 389 (1976). 81. J.Killeen, A. A.Mirin and NLE.Resnick, Methods in Computational Physics, 9_, 421 (1970). 82. J. Killeen and AJi. Futch, Jr., J. Computet Phys. 2, 236 (1968). 83. N. Corngold, Phys. Rev. A 15_, 2454 (1977). 84. N. Corngold, Phys. Rev. A 2A 656 (1981). 85. A.J.M. Garrett, Phys. Reports 124, (1986). 86. P.W.Atkins, Molecular Quantum Mechanics. (Oxford Univ.Press, 1970) 204-211. 110 87. P.Langhoff, in Theory and Applications of Moment Methods in Many Fermion  Systems, edited by B.J.Dalton, S.M.Grimes, JJP.Vary and C.A.Williams (Plenum Press, New York, 1980). 88. W.P. Reinhardt, in Atomic Scattering Theory, edited by J. Nuttall (Univ. Western Ontario, 1979). 89. D.Montgomery, GJoyce and LTurner, PhysHuids 1L. 2201 (1974). 90. T.E.Sharp, J.R.Eyler and E.Li, IntJ.Mass.Sprectrom.Ion Phys., 9_, 421 (1972). 91. D.Anderson and MJJsak, J.Plasma Physics, 21 107 (1986). 92. T.H.Stix, NucLFusion, i i , 737 (1975). 93. D.Anderson, L.-G.Eriksson and M.Lisak, Plasma Phys.Contr.Fusion, 29, 891 (1987). 94. R.C.Davidson, Methods in Nonlinear Plasma Theory (Academic Press, 1972). 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

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>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0059427/manifest

Comment

Related Items