A VLASOV PLASMA SIMULATION CODE By PETER KARL LOEWENHARDT B.Sc. Hons., The University of Saskatchewan, 1987 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in THE FACULTY OF GRADUATE STUDIES (DEPARTMENT OF PHYSICS) We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA February 1989 © Peter Karl Loewenhardt, 1989 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 Physics The University of British Columbia 1956 Main Mall Vancouver, Canada V6T 1Y3 7 DF.fin/ft-n ii ABSTRACT A 1-1/2 dimensional, electromagnetic Vlasov plasma simulation code, relativistically correct, was constructed and tested. The code can deal with one or two species, each with a Maxwellian distribution and a possible drift velocity. The code also allows external fields, such as a laser, to be included in the simulations. Simulations of Landau damping, the two-stream instability and stimulated Raman scattering were carried out and compared with theory and with the results of the particle code EMI. When dealing with electrostatic problems, the Vlasov code gave results agreeing very closely with theory in both non-relativistic and relativistic regimes. Here the Vlasov code preformed better than EMI which had problems with background noise and longer run times. When a laser field was included within the simulations, however, the Vlasov code produced some spurious features, unlike EMI. These anomalous features may be caused by aliasing, recurrence or by some other unknown effect. Since realistic results are produced as well, it is believed that this problem can be overcome. When a very high intensity laser was included in the simulations the Vlasov code produced much better results but was plagued by very long run times. iii TABLE OF CONTENTS Abstract ii Table of Contents iii List of Tables vi List of Figures vii Acknowledgement xi Introduction 1 Chapter 1 The Plasma Dispersion Relation 1.1 Theory 4 1.2 Dispersion Relation for a Maxwellian Distribution 9 Chapter 2 The Two-Stream Instability 2.1 Counter Streaming Cold Beams 12 2.2 Counter Streaming Warm Beams 14 2.3 Modification of the Warm Beam Treatment 19 2.4 The Corrected Growth Rates 22 Chapter 3 Stimulated Raman Scattering 3.1 Coupled Equations Describing Raman Scattering 24 3.2 Dispersion Relation 30 Chapter 4 Particle Code Theory 4.1 Computational Cycle 34 4.2 Integration of Equations of Motion 35 iv 4.3 Loading the Particle Distribution 41 Chapter 5 The Program EMI 5.1 Initialization 44 5.2 The Simulation Loop • 48 Chapter 6 Scheme for the Vlasov Algorithm 6.1 Introduction 51 6.2 Non-Relativistic Landau Damping 54 6.3 Relativistic Landau Damping 57 6.4 How to Calculate the Longitudinal Field 59 6.5 Including External fields 60 Chapter 7 The Vlasov Code 7.1 Introduction 65 7.2 Initialization Procedure 69 7.3 The Simulation Loop 72 Chapter 8 Landau Damping Simulations 8.1 Non-relativistic Landau Damping 78 8.2 Relativistic Landau Damping 87 8.3 Phase Space Waves 92 8.4 Particle Trapping 98 Chapter 9 Warm Two-Stream Instability 9.1 Non-Relativistic Case 103 9.2 Relativistic Case 116 Chapter 10 Stimulated Raman Scattering Simulations 10.1 Low Laser Intensity: Vlasov Code 123 V 10.2 Low Laser Intensity: EMI 136 10.3 High Laser Intensity 141 Chapter 11 Conclusions 146 References 150 Appendix A Program Listings for Theoretical Results A.1 Non-Relativistic Landau Damping 153 A.2 Fluid Theory Growth Rates for the Two Stream Instability =..157 A.3 Necessary Two Stream Instability Velocity Spread , 159 A. 4 Determination of the Two Stream Instability Growth Rates 160 Appendix B Listing of the Vlasov Code B. l Initialization Procedure 171 B.2 Simulation Loop 182 B.3 Cubic Spline Routines 200 B. 4 Fast Fourier Transform Routines 204 Appendix C User's Guide C l Compiling and Running the Vlasov Code 207 C. 2 Plotting Routines 209 vi List of Tables 1 Radix-2 digit reversal 43 2 Radix-3 digit reversal 43 3 EMI input parameters 45 4 Vlasov input parameters 70 5 Vlasov Landau damping input parameters 79 6 EMI Landau damping input parameters 83 7 Non-relativistic Landau damping results 86 8 Vlasov relativistic Landau damping parameters 88 9 EMI relativistic Landau damping parameters 90 10 Vlasov relativistic Landau damping results 92 11 EMI relativistic Landau damping results 92 12 EMI two stream instability input parameters 104 13 Vlasov two stream instability input parameters Ill 14 Two stream instability results 116 15 EMI relativistic two stream instability input parameters 117 16 Vlasov relativistic two stream instability input parameters 120 17 Relativistic two stream instability results 122 18 Vlasov Raman scattering input parameters 125 19 EMI Raman scattering input parameters 136 20 Vlasov diagnostic routines 211 21 EMI diagnostic routines 225 vii Lis t of F igures 1 The Landau Contour 8 2 Vlasov treatment (solid line) and linear approximation (dotted line) of the real frequency u>r 10 3 Vlasov treatment (solid line) and linear approximation (dotted line) of the damping rate — u>j 11 4 Warm beam growth rates for T = lOOeV (fluid treatment) 17 5 Maximum Mode excited versus beam velocity spread 22 6 Warm beam growth rates for T = lOOeV (Vlasov treatment) 23 7 Vlasov (solid line) and fluid treatment (dotted line) warm beam growth rates 23 8 Computational cycle in EMlt 2 0! 34 9 Leap-Frog method of integration!20' 36 10 Initial variable placement in time 66 11 A Vlasov code time step advancement 68 12 Vlasov electrostatic energy mode 1 time evolution 81 13 Vlasov total electrostatic energy time evolution 82 14 EMI electrostatic energy mode 1 time evolution 85 15 EMI total electrostatic energy time evolution 85 16 Comparison with theoretical real frequency 86 17 Comparison with theoretical damping rate 87 18 Vlasov electrostatic energy mode 1 time evolution 88 19 Vlasov total electrostatic energy time evolution 89 20 EMI electrostatic energy mode 1 time evolution 90 21 EMI total electrostatic energy time evolution 91 22 Vlasov distribution function at t = 0 94 viii 23 Vlasov distribution function at t — Su;"1 94 24 Vlasov distribution function at t = I2u~x 95 25 Vlasov distribution function at t = 40CJ~1 95 26 EMI phase space at t = 0 96 27 EMI phase space at t = SLO'1 96 28 EMI phase space at t = 12a;-1 97 29 EMI phase space at t = 40a;-1 97 30 EMI velocity distribution function at t = 0 98 31 EMI velocity distribution function at t = 40a;-1 98 32 Vlasov electrostatic energy mode 1 time evolution 99 33 Vlasov electrostatic energy mode 2 time evolution 100 34 Vlasov electrostatic energy mode 3 time evolution 100 35 EMI electrostatic energy mode 1 time evolution 101 36 EMI electrostatic energy mode 2 time evolution 10L 37 EMI electrostatic energy mode 3 time evolution 102 38 EMI mode 1 time evolution 106 39 EMI mode 2 time evolution 106 40 EMI mode 3 time evolution 107 41 EMI total electrostatic energy time evolution 107 42 EMI phase space at t = 0 108 43 EMI phase space at t = 12a;-1 108 44 EMI phase space at t = 24a;"1 109 45 EMI particle velocity distribution at t = 0 109 46 EMI particle velocity distribution at t = 24a)"1 110 47 Vlasov mode 1 time evolution 112 48 Vlasov mode 2 time evolution 112 49 Vlasov mode 3 time evolution 113 ix 50 Vlasov total electrostatic energy time evolution 113 51 Vlasov distribution function at t = 0 114 52 Vlasov distribution function at t = 15W"1 115 53 Vlasov distribution function at t = 32a;"1 115 54 Comparison of measured and theoretical growth rates 116 55 EMI mode 1 time evolution 118 56 EMI total electrostatic energy time evolution 119 57 Vlasov mode 1 time evolution 120 58 Vlasov total electrostatic energy time evolution 121 59 Comparison of measured and theoretical growth rates 122 60 Vlasov electrostatic energy mode 5 time evolution 127 61 Vlasov electrostatic field Ex at t = 65CJ"1 127 62 Vlasov electrostatic field Ex at t = l"iQu>~1 128 63 Vlasov transverse electric field mode 4 time evolution : 128 64 Vlasov transverse electric field mode 1 time evolution 129 65 Vlasov average mode distribution 129 66 Vlasov electron distribution function at t = 65a?"1 130 67 Vlasov electron transverse velocity at t = 65u;~1 131 68 Vlasov electron distribution function at t = 104a;"1 131 69 Vlasov electron distribution function at t = 130W"1 132 70 Vlasov electrostatic energy mode 5 time evolution 133 71 Vlasov transverse electric field mode 4 time evolution 133 72 Vlasov transverse electric field mode 1 time evolution 134 73 Vlasov electron distribution function at t = 120W"1 134 74 Vlasov ion distribution function at t = 120W"1 135 75 Vlasov total electrostatic energy time evolution 135 76 EMI electrostatic energy mode 5 time evolution 138 X 77 EMI transverse electric field mode 4 138 78 EMI transverse electric field mode 1 139 79 EMI electron phase space at t = 220a--1 139 80 EMI electrostatic energy mode 5 time evolution 140 81 EMI electron phase space at t = 350a;-1 141 82 EMI ion phase space at t = 350a)-1 " 141 83 Vlasov electrostatic energy mode 5 time evolution 143 84 Vlasov transverse electric field mode 1 time evolution 143 85 Vlasov transverse electric field mode 4 time evolution 144 86 EMI electrostatic energy mode 5 time evolution 144 87 EMI transverse electric field mode 1 time evolution 145 88 EMI transverse electric field mode 4 time evolution 145 xi A C K N O W L E D G E M E N T Special thanks must go to my supervisor, Dr. A. J. Barnard, for all his patience and help with this project. I would also like to thank Marc Majka for supplying the Superplot routines and C. K. Birdsall and Bruce I. Cohen for the listing of EMI. The financial support of the National Sciences and Engineering Research Council was greatly appreciated. I would also like to thank my wife for her support and putting up with all my long nights. INTRODUCTION Computer simulation in physics bridges the gap between theory and experiment. These simulations can present a physical concept which narrows possible scenarios in the-ory and experiment. Simulations can look into regions where analytic theory fails and can help an experimentalist explain some physics behind an observable process. Computer simulation of plasmas usually deals with solving the fluid equations, moving representa-tive particles self-consistently (particle codes) or with kinetic theory (solving the Vlasov equation). The construction of a code solving the Vlasov equation is presented here and comparison is made with theory and the results of a well known particle code EMI. A 1-1/2 dimensional electromagnetic Vlasov code was constructed for many reasons. Particle codes dominate plasma simulations, so a fully viable Vlasov code would provide a different approach to plasma problems, perhaps bringing new insights into view. It would certainly provide a check on the particle codes' results. The restriction to 1-1/2 dimen-sions still allows many phenomena to be simulated and reduces the computational power required to a reasonable level. The capability of the Vlasov code to include transverse radiation fields allows a much larger number of processes to be modelled than a simple electrostatic code. Processes of interest to inertial confinement fusion occurring when an intense laser exists within a plasma, such as stimulated Raman and Brillouin scattering, can be modelled. Even without the fusion application these and other electromagnetic processes are of great physical interest. Periodic boundary conditions only allow the mod-elling of large homogeneous plasmas, leaving out features such as density gradients and thin foils. A process should be understood first within these restrictions before looking at a non-homogeneous, higher dimensional system. The basic single distribution Vlasov code written by Dr. A. J . Barnard was a starting point in this endeavor. The code was transferred from a main-frame computer system at U.B.C. to a Sun 3/60 workstation where all the alterations and simulations 1 Introduction 2 took place. The particle code EMI (and its counterpart EM1B) were also transferred from the main-frame to the Sun. The FORTRAN compiler on the Sun did not support as many features as FORTRANVS and thus modifications had to be made to each code. An example of this is that namelist was not supported, requiring a reconstruction of input techniques. As well, variable declarations and include files had to be re-written for EMI to run correctly on the Sun. Diagnostics were added to the Vlasov code to facilitate observation of the simulations. As well, extensive graphics routines were written for both the Vlasov code and EMI to view these diagnostics, utilizing the Superplot graphics package. Some diagnostics in EMI also had to be altered since EMI was originally written using another graphics routine. The algorithm used in the Vlasov code was checked and found to be correct rela-tivistically for electrostatic problems. However, it was found that when external radiation fields were present in the plasma the basic algorithm had to be modified. Modifications to the velocity and spatial shifts were calculated and these second order effects were im-plemented in the Vlasov code. The transverse velocity in the code was also found to be a function of velocity as well as space, requiring major modifications to the code, especially in shifting the distribution function and in calculating the transverse current. The Vlasov code was also expanded to include the distribution functions for two species. This again required extensive alterations to the simulation code. The modi-fications were carried out so that the second species' distribution could be dealt with either non-relativistically or relativistically, depending on the mass ratio. As well, the addition of a second distribution required a change in the way the longitudinal field and the transverse current were calculated. Input and diagnostic routines were modified to accommodate this second distribution function. As well, the normalization of the initial Maxwellian distribution for each species, originally calculated by integration, was altered to be carried out using the modified Bessel function. The Vlasov code was then tested against theory and the results of EMI. Landau damping simulations were done and damping rates and frequencies compared with the Introduction 3 solution of the plasma dispersion function. Relativistic modifications to these values pre-dicted by theory also provided a test for the Vlasov code. The two-stream instability was tested in non-relativistic and relativistic regions, the growth rates from the simulations being again compared with the solution of the plasma dispersion function. Lastly, stim-ulated Raman scattering simulations were carried out for a lower intensity and a higher intensity radiation field. The organization of the thesis is as follows: Chapter 1 discusses basic theory of the plasma dispersion relation and the resulting numerical solution is compared with linear theory. Chapter 2 discusses the theory and physics behind the two-stream instability. It is shown that fluid theory does not produce the correct growth rates for the instability with warm Maxwellian beams. Instead, the plasma dispersion relation had to be solved to obtain theoretical growth rates. Chap-ter 3 discusses stimulated Raman scattering and the dispersion relation describing the instability is derived. Chapters 4 and 5 deal with the construction of the particle code EMI and the method in which it operates. Chapter 6 shows that the Vlasov algorithm is correct electrostatically but with external radiation fields present modifications must be made. Chapter 7 explains the layout and operation of the Vlasov code. Chapters 8, 9 and 10 give the results of the Landau damping, two-stream instability and stimulated Raman scattering simulations respectively. Appendix A gives an explanation and listing of the programs written to derive the theoretical results the simulations were compared with. Appendix B gives a complete listing of the Vlasov code. Appendix C explains how to use the Vlasov code and gives a listing of the plotting routines developed for the Vlasov code and EMI. CHAPTER 1 THE PLASMA DISPERSION RELATION 1.1 Theory The plasma dispersion relation is an important mathematical tool used to describe phenomena in plasmas. It relates basic properties of a plasma to the function describing the particle distribution. One of the initial endeavors to reveal the properties of the plasma dispersion relation was carried out by J. D Jackson in 1960^. The Vlasov equation in one dimension with no applied external fields can be written as at ox m av which, when coupled with Poisson's equation — = innee(l- I f dv) (1-2) form a set of fundamental equations describing a simple collisionless plasma. Here / is the distribution function f(v,x,t), normalized in an unperturbed state (/ = /o) so that / fQdv = l . (1-3) We have assumed, for the moment, that any behaviour of the plasma is of high enough frequency that the electrons can be seen as moving in a background of fixed neutralizing ions. Chapter 1 The Plasma Dispersion Relation 5 One can linearize these equations by assuming that there is a small perturbation to the system, i.e. / = /o + / i , where f\ is the small perturbation to the unperturbed distribution function. Assuming that there are no external fields so that E in the Vlasov equation is of first order and only keeping terms up to and including first order, one finds £ 4 + t,££._f^£ft = 0 dt dx m dv (1-4) and ^ = -47rn e e / dv Ox J (1-5) Multiplying the Vlasov equation by exp (—ikx + iut) and integrating over time and space, one obtains + oo oo + oo oo kx + iut -oo 0 -oo 0 + 00 oo - -1 / dx fdtE^ = 0 m J J av -oo o (1-6) The time integration in the first term can be done by parts, oo Jdt-bWl fie tut oo : J dt fiueiut (1-7) Similarly, the spatial part of the second term of the Vlasov equation can be integrated by parts, +oo J dzv?£l-e-ik*=vf1t-ik* dt + i / dxkvfxut™1 , (1-8) — oo J Chapter 1 The Plasma Dispersion Relation 6 where the first term vanishes for a realistic f\. Thus we have +°° t=oo +°° 0 0 j dx f i e - i k x + i u t + J dx J dte-ikx+i^{(-iw+ikv)fi-^E^} = 0 . (1-9) -00 0 Assuming that u is in the upper half of the complex w-plane (then tut oc —U{t), the first term of equation (1-9) becomes +00 J cixe-'*x/i(^z,0) = *(",*) • (1-10) From the properties of the Fourier transform, it follows that equation (1-9) can be written as •$(v,k) + (-iu + ikv)f!{v,k,uj) - — E{k,u)~ = 0 ( I - H ) m uv giving h(vtk,u>) = 1 U{Vlk) + -E(klU)^S.\ . (1-12) i[kv — LO) < m av ) The Fourier transform of Poisson's equation can be solved to find the electric field; ikE = -4?rnee J fx dv (1 - 13) +00 +00 \k2 — / J 1 ,, dv }E(Lo,k) = iimee / v /, dv . (1-14) I. m J v—LO/k J J v — LO/k One can then define a function +°° k LO J V — LO/k —00 Chapter 1 The Plasma Dispersion Relation 7 where dfo \jr\V) = and, as always, G(v) = ^ (1 - 16) Ane2ne F m The expression for the electric field then becomes ,-,/, s '47rnee t §(v,k) , . „. E(k,u) = 2 / ^ / dv . 1-18 e(k,u>/k) J v — uj/k To obtain the behaviour of a mode with respect to time, this can be Fourier transformed to E(k, t): B ( M )« 3 2 i / ^ . - . " _ > / i ^ L . ( 1_1 9 ) Jc e(k,uj/k)J v-u/k The behaviour of the system is governed by the zeros of e(k,w/k), that is, from the equation 4=/"^* • -20) since for any realistic initial conditions no other terms have zeros. Equation (1-20) is the usual dispersion relation obtained from the Vlasov equation. Due to the assumptions made above in time-integrating the Vlasov equation, e is only defined in the upper half of the complex-u; plane. Therefore, to find the solution to equation (1-20), e must be analytically continued to the entire complex-u; plane. To do this, one makes use of the relation Chapter 1 The Plasma Dispersion Relation 8 where p is vanishingly small. Note that where e is discontinuous across the real-w axis, a jump is found from equation (1-21), - - P I ( G(v) wj J V — u — + p \ r G(v) J V — u = -2iriG(u) . The analytic continuation of e then becomes £(w) = { _ t _GiiL d v J v—ufk ' Im(u>) > 0 (1 - 22) % ~ I d v ~ 2lTiG(w/*) ' J f f l M < 0 • This is equivalent to evaluating the integral along a Landau contour as in figure 1 (1 - 23) Vi Vr V=w/k Figure 1: The Landau Contour. Chapter 1 The Plasma Dispersion Relation 9 1.2 Dispersion Relation for a Maxwellian Distribution The usual form of the unperturbed velocity distribution function is that of a non-relativistic Maxwellian, /o = -j^-e-^l^ (1 - 24) V27T Vth where vth is the thermal velocity in one-dimension. This gives G(v) = 7=L— ve-vi/2v>» . (1-25) th With an interest in Landau damping, assume that there are no growing oscillations so that the roots of e lie in the lower half of the complex-w plane. The dispersion relation for the plasma is then derived from e(u,u>/k) = 0, giving k2v2h 1 +Fve~v2/2 1 j v±Z—dv + V^izt-*2'2 =0 (1-26) J v-z where the substitution z = u>/kvth was made for simplicity. The analytic continuation of e was used since it was assumed that Im(u>) < 0. The complex error function'2! W(z) = e-z2(l + -^= j* ef d?) (1-27) can be used to write this dispersion relation as 0 = ^ k + l + i.fl2W(^) . (1-28) u2 p 2 v/^2 A comparison of the numerical solution of this equation with usual linear theory is pre-sented in figures 2 and 3. The approximations are given by the dotted curves. Figure 2 Chapter 1 The Plasma Dispersion Relation 10 compares the real frequency with the Bohm and Gross dispersion relation (1 - 29) Figure 3 compares the resulting damping rate with the linear theory result (1 - 30) This natural damping is Landau damping, which arises from the pole in the plasma dispersion function and is a collisionless process. As can be seen from these graphs, linear theory is only valid for small damping rates. Figure 2: Vlasov treatment (solid line) and linear approximation (dotted line) of the real frequency ur. Chapter 1 The Plasma Dispersion Relation '-1.0 -0.S 0.0 0.6 1.0 1.9 2-0 LOG(K*Vth/Wp) Figure 3: Vlasov treatment (solid line) and linear approximation (dotted line) df the damping rate — a>j. CHAPTER 2 THE TWO-STREAM INSTABILITY The two-stream instability occurs when two species of particles in a plasma are moving relative to each other. This instability can arise from driving a current through a plasma, inducing a relative drift velocity between different species'3!. It can also occur for two counter-streaming beams of particles. The drift energy of the particles excites density perturbations, the drift energy being converted to oscillation energy'4!. Initially a small density perturbation causes modulation in the velocity of a beam'5!. This causes the space charge to bunch and large potentials are set up. Thus the perturbation supports itself and begins to influence the other beam, resulting in an amplified feed-back loop. The perturbations grow exponentially in time'6' until saturation, since the forces between the streams reinforce each other. When viewed from phase space, the two-stream instability forms a characteristic vortex shape at and beyond saturation'7,8!. This vortex shape is a result of an exchange of energy between particles and the wave created by the instability. The particles rotate in phase space from their loss and gain of energy from the wave, the system oscillating around some mean energy level'9!. This behaviour actually reveals particle trapping'10!, since the particles which become trapped in potential troughs arising from the instability follow closed orbits in phase space, thus forming the vortex pattern'4!. Examples of this behaviour in phase space will be presented later with the simulation results. 12 Chapter 2 The Two-Stream Instability 13 2.1 Counter Streaming Cold Beams One can consider the simple case of two electron beams counter-streaming in a background of neutralizing and fixed ions'4'5). Suppose velocities of the beams are equal and opposite, ±VQ , and for simplicity assume a cold plasma and no applied magnetic field (T = 0, B = 0). Since this problem deals with streams of particles, the plasma fluid equations cou-pled with some of Maxwell's equations will most likely applyt4'1^. The fluid equations are derived from taking moments of the Vlasov equation'1 J' % + Z-Vf + ^-(E + vxB)-% = 0 (2-1) at m av giving the continuity and momentum equations dn3 dt + V • K ) = 0 (2-2) l + ^ . V ^ i i ^ ^ x i ) - - ^ . (2-3) at m3 c msn3 Here the subscript s denotes the species of particles, ns the density of species s and vs the fluid velocity of species s. Let s = a , j3 for the individual beams, corresponding to the drift velocities +VQ, —VQ respectively. To obtain a dispersion relation for the linear phase of this instability, one linearizes and Fourier transforms the fluid equations, including the above assumptions. From this, the following dispersion relation is obtained'11,4J, + 7—^-72 • (2-4) up (to — kvo) (to + kvo) This is the dispersion relation describing the two-stream instability for two cold counter-streaming electron beams each of density no in a background of fixed ions. Chapter 2 The Two-Stream Instability 14 Conditions for the growth of this two-stream instability can found and the theoret-ical maximum growth rate can be calculated from the above dispersion relation'4!. The requirement ^ < V 2 . (2-5) Up must be satisfied to have an onset of the instability. The maximum growth rate of this instability turns out to be 2.2 Counter Streaming Warm Beams One can now generalize to a more realistic situation of warm (T ^ 0) electron beams streaming with equal and opposite velocities in a background of fixed neutralizing ions. This means that the pressure term dropped in the cold beam analysis must be retained. Assuming that to zeroth order there are no significant temperature gradients during the initial growth of the instability, one can write the pressure term in a more convenient form by noting that p = nT: Vp TVna mens men3 (2-7) Let T = T/me for simplicity and linearize the fluid equations with 5 = 0. Therefore, after Fourier transforming the momentum equation, one obtains • - i u ~ ~ e ^ i -Tknsi , . —iuvs\ + tkv30vsi = i . (2 — a) Chapter 2 The Two-Stream Instability 15 The continuity equation is unchanged. Dropping the tilde for the Fourier coefficients and solving for v3i, one obtains = -ieEJrrie _ 3 1 (u-kv3Q)-k2T/(uj-kvao) " ^ ' This then gives an expression for nsl: -ikns0eEi/me n3i = 7y . (2 — iu) (w - kva0)2 - k2T Poisson's equation is used to link the density and the electric field perturbations. Noting that each beam has a density of no and that va,po = ±urj> the following dispersion relation for the case of warm counter-streaming beams is obtained: \ = K + K • (2-H) u2 (u-kvQ)2-k*T (w + kvQ)2 - k2T If T = 0 the cold beam dispersion relation is recovered. As before, the above dispersion relation can be used to estimate the maximum growth rate for this instability. Combining the denominators in the the dispersion relation, one obtains ' _1_ = (u + kv0)2 -2k2T + (uJ-kv0)2 <^P (w2 -k2v2)2 + k*T2 -{uj-kvQ)2k2T-{u; + kv0)2k2T and then, ^-2u2(k2(v2 + T)+u2p) + k2v2(k2(v2-2T) + ^ - + - ^ . (2-13) Chapter 2 The Two-Stream Instability 16 This can again be solved using the quadratic formula as in the cold beam case, giving 4 = ^ + * £ + 1 A + « a i + « a . ( 2 _ 1 4 ) Since we are searching for an imaginary u we chose the negative sign in the quadratic formula (the only way u>2 can be negative). To simplify this expression let • *-£ (2-15) P y2 = - 2 (2-16) uj and p = vl . (2-17) Also set LO — i-y as before to search for a pure imaginary root. Then equation (2-14) can be written as 7 2 = y/l + \x2pT + ixp - (xp + xT + 1) . (2-18) Fix the temperature T. Also fix the mode excited so that comparisons with the cold beam case can be made. One can then differentiate equation (2-18) with respect to p = VQ to maximize the growth rate with respect to the beams' relative velocities. Differentiating equation (2-18), y + x = 0 (2-19) dp 2y/l +Ax2pT + ixp gives after some simplification Chapter 2 The Two-Stream Instability 17 The maximum value of the growth rate can then be found by substituting this result back into equation (2-18) and simplifying: 7 L * = 2\A + *T + x2T2 - 1 - 2xT -4(xT + l) ' (2-21) Noting that Im(u>) = 7wp from the definition of y2, one finds Im(u>) = Lop n / k2T k4T2 2k2T 2jl + —r- + —r--l-tot io2 4(k*T/uj + 1) (2 - 22) for the maximum growth of warm counter-streaming electron beams as derived from the fluid equations. A plot of the theoretical growth rates for a range of relative beam velocities is given in figure 4, where equation (2-14) was solved for T = lOOeV (see appendix A). o a: o Beam Velocity (VO/c) Figure 4: Warm beam growth rates for T = lOOeV (Fluid treatment). Chapter 2 The Two-Stream Instability 18 An estimate can now be made of the necessary velocity difference between the two warm beams for the instability to occur. Looking back to equation (2-14), we find that the condition for growth of the instability is only satisfied if yjio* + 4kH2T + 4k2v2u2 > k2{vl + T)+Wp . After collecting terms of T together, this can be written as 0 > T2 + 2T(^-v2) + v*--j^ . (2-23) The roots of the quadratic equation on the right hand side of this inequality are T = v2 (2 - 24) and (2 - 25) In order to find the region of T that gives the correct inequality for growth, let T = VQ — u2/k2. This value of T lies between the two roots given above. Substituting this value back into the inequality one obtains 0 > —u>p/k2, which is always true. Therefore, since equation (2-23) is quadratic in T, the instability will occur only if T lies between the two roots equation (2-24) and equation (2-25). The lower bound on T is a generalization of the instability condition of the cold beam case, since if T = 0 the inequality for growth of the cold two-stream instability, kvo/up < y/2, is obtained. The upper bound of T < v2, can be seen as vth < VQ, where vth is the thermal velocity of the particles in the beam. This upper bound implies that the relative velocity spread of the beams should at least be greater than the thermal spread of the beams for the instability to occur. In light of the assumptions made in the fluid treatment, this result is somewhat surprising and will Chapter 2 The Two-Stream Instability 19 be shown later to be close to the correct value of 1.307Vth for a Maxwellian distribution of particles. 2.3 Modification of the W a r m Beam Treatment The fluid theory derivation just presented for warm counter-streaming beams has one major assumption within it, namely neglecting the velocity spread in the warm beams. Indeed, all the parameters dealt with were implicitly assumed to be single-valued functions of space and time, or average values'1). The physical assumption made is that the relative velocity between the two beams is much greater than the thermal spread of the beams. In general, the thermal spread of the beams will lessen the growth rate, since particles of different velocity will react differently to the perturbation, bunching at different places and thus reducing the amplification causing the instability'5!. Other methods can be used in dealing with the necessary relative velocity of two counter-streaming warm beams. Indeed, it can be shown'1! that a single-humped distri-bution function (i.e. a Maxwellian) will not give rise to an instability. Nyquist analysis can be used to show, however, that a significantly double-humped distribution will give rise to an instability'11!. The requirement presented in the Nyquist analysis for the insta-bility to occur is, as before, that the relative velocity of the beams must be greater than their thermal spread. The only method which derives an exact answer, however, is the one which utilizes the fundamental Vlasov equation and the resulting plasma dispersion relation derived in chapter 1. From the Vlasov approach, one can find the minimum relative velocity of the warm beams that will allow growing oscillations'1'. Imagine starting with a system of two counter streaming warm beams with only a small relative velocity, thereby having a velocity distribution close to a single-humped distribution, and thus e(u>) only having roots in the lower half of the complex-u; plane. Increasing the relative speed of the beams Chapter 2 The Two-Stream Instability 20 makes the velocity distribution deviate farther from a single-humped distribution, and the roots of e(u>) will begin to rise towards the real-w axis. If the relative speed of the beams is increased further, the roots will eventually cross the real-w axis and a transition from damped to growing oscillations will be made'1'. The critical state of the system is, then, when the roots of e(ui) lie on the real-w axis. We look for solutions of the dispersion relation with u/k = x + iy,y = 0, to find this critical point. The real and imaginary parts of e(u) = 0 then become Re(e) = K;-P f dv = 0 (2 - 26) L>* J V — X and Im(e) = -nG(x) = 0 , (2-27) where the term in the imaginary part comes from the analytic continuation of e(a>) to the entire complex-^ plane (see chapter 1). Note that equation (2-27) implies that the phase velocity x = tur/k of the wave equals the position of zero slope in the distribution function. Also note that G(x) = 0 turns the principal value integral of equation (2-26) into an ordinary integral. Now we have a relation to estimate the minimum relative velocity necessary for an instability to occur. Assuming that each streaming electron beam has a non-relativistic Maxwellian distribution, then fo is and (2 - 28) (2 - 29) Letting and G(v) becomes Chapter 2 The Two-Stream Instability 21 m = (2 - 30) m = . (2-31) Vth G = - ^ L e - ^ - m e - ^ \ . (2-32) We want to find the minimum velocity difference V = vo - {-vo) = 2v0 = vth(m + (2-33) for which the instability occurs. Note that the condition G(x) = 0 is satisfied trivially for r]i = rj2, other solutions not being important for the problem at hand'1!. The critical condition (equation (2-26) and equation (2-27)) becomes with this change in variables —oo —oo where as well the variable change of r\ —> t was made and TJ(X) is t] to be evaluated at v = x = ur/k. Utilizing the complex error function W(z) and defining rj = 771 =772, one finds that the above relation can be written as k2v] = 2 ( n y | / m { W { ± ) } - 1) . (2-35) This equation was solved numerically (see Appendix B) and a graph of the function is given in figure 5. As can be seen from this graph, the minimum velocity spread for any growth is 2.614ut/>, giving the requirement of VQ > 1.307UM for an instability to exist. Chapter 2 The Two-Stream Instability 22 o °1 -o 10 11 12 13 1« 15 Velocity spread / Vth Figure 5: Maximum mode excited versus beam velocity spread. 2.4 The Corrected Growth Rates One can directly obtain the growth rates of the warm two-stream instability by numerically integrating the real part of the dispersion relation to2 with a relativistically correct G(v) (listing in Appendix B). An example of this is given in figure 6 and a comparison with the less correct fluid equation treatment in figure 7. The parameter values chosen in all the figures were k = 0.5, cop = 0.1, c = 3.927 and T = lOOeV, obtained from runs of the particle code EMI. Indeed, from these graphs and the simulations to be presented, it seems that neglecting the temperature spread of the beams can lead to quite incorrect results. vGjv) v2 + y2 dv (2 - 36) Chapter 2 The Two-Stream Instability 5 ? -o -O or o o - ,— i—•—i—i— i—'— i—>—i— '— i— '— i— ' i—• i 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 o.oa 0.00 B e a m V e l o c i t y ( V O / c ) Figure 6: Warm beam growth rates for T = lOOeV (Vlasov treatment). B e a m V e l o c i t y ( V O / c ) Figure 7: Vlasov (solid line) and fluid treatment (dotted line) warm beam growth rates, T = lOOeV. CHAPTER 3 STIMULATED RAMAN SCATTERING The process of stimulated scattering arises when an intense electromagnetic wave interacts with a plasma, coupling into a scattered electromagnetic wave and an excited plasma mode'12'. Stimulated Raman scattering occurs if a plasma wave is excited, while stimulated Brillouin scattering occurs if an ion acoustic wave is excited. The incident radiation must be of large enough amplitude to have the plasma perturbation overcome Landau damping and for the electromagnetic waves to overcome collisional damping'13!. Stimulated scattering is of considerable interest, since the process impedes the transmis-sion of laser light to the critical plasma density at which enhanced absorption occurs'14'. This then can prevent the correct absorption of energy in inertial confinement targets and can interfere with laser heating of magnetically confined plasmas'15'. This discussion will concentrate on stimulated Raman scattering. 3.1 Coupled Equations Describing Raman Scattering Stimulated Raman scattering occurs when a large electromagnetic wave is incident upon a plasma, resonantly decaying into a scattered electromagnetic wave and an elec-tron plasma wave'16'. As already stated, this process can possibly prevent laser light from reaching the critical surface in laser fusion applications. It can also prevent efficient com-pression of the fusion target by preheating the target core with energetic electrons. The electron plasma wave produced in the instability grows until its associated electrostatic fields are large enough to accelerate individual electrons, saturating the growth of the plasma wave and producing the energetic electron pre-heat problem'14,13,17'. 24 Chapter 3 Stimulated Raman Scattering 25 The process of stimulated Raman scattering can be seen as a feedback loop'16!. Imagine a plasma with a density fluctuation 8n due to a plasma wave. The incident radi-ation electric field E accelerates electrons in the perturbed plasma creating a transverse velocity V?L. This induces a transverse current jt = —eui^n which, if the wave numbers and frequencies of the system are properly matched, generates a scattered electromagnetic wave with an electric field of amplitude SE. This scattered radiation interferes with the incident radiation, producing variations in the wave pressure V(E2/8ir) = V(E • 6E)/4it which accentuate the density perturbations. Thus a feedback loop is possible. The necessary frequency and wave number matching conditions are wo = W j+ue (3-1) and k0 = ka + ke (3 - 2) where u>, k denote frequency and wavenumber respectively, with the subscripts 0, 5 re-ferring to the incident and scattered electromagnetic waves and the subscript e to the electron plasma wave. Recalling that u2 = u2p + c2k2, u2a=u2p + c2k2 (3-3) and u2e=u2p + 3v2hk2e , (3-4) we see that equation (3-1) becomes U0 = y/uj+<?k* + y/u,j+3v*kkl , (3-5) implying that LOQ > 2u>p . (3 — 6) Chapter 3 Stimulated Raman Scattering 26 This means that stimulated Raman scattering can only occur at densities n < ncr/4, simplicity, imagine that a large electromagnetic wave is incident upon a homogeneous plasma of uniform temperature and density. Assume that the density fluctuations asso-ciated with the plasma waves of the process are rapid enough in nature that there is no coupling with ion motion and therefore the ions can be treated as a fixed, neutralizing background. Also assume that the plasma perturbations are one-dimensional in nature, that is, Sn is along the direction of propagation of the incident light wave. One begins with Maxwell's equations After combining the curl of equation (3-7) and the time derivative of equation (3-8) and using the vector identity ncr being the critical density of the plasma at which U)Q = Anncre2/me. One can derive coupled equations describing the Raman instability'12'15'16!. For (3-7) „ - IDE 4TT-V x B = -— + —j c at c (3-8) V x V x £ = V ( V . | ) - V2l (3-9) one obtains /I d2 r-,o.^ An dj _ , _ ( ? ^ - v 2 ) £ = - ^ - v ( v . £ ) (3 - 10) Now separate terms into components: E = E~T + Ei, ]=JT + j~L (3-11) Chapter 3 Stimulated Raman Scattering 27 where the subscript T denotes the transverse component which has a curl but no diver-gence and the subscript L denotes the longitudinal component that has a divergence but no curl. Combining Poisson's equation V • E = iirp ( 3 - 1 2 ) and the conservation of charge equation ^ = - V T T J ( 3 - 1 3 ) one obtains 7^ . ( dt Now V • = 0 and V • ET = 0 , implying that that 1 d 2EL _ dfL 4TT dt2 dt Noting that V • E = V • EL, equation ( 3 - 1 0 ) can be written as V - ( ^ + 47rj) = 0 . (3-14) (3-15) (_L|_ _ ^ )ET + ^ = V(V • EL) + V2EL . (3 - 16) Using the vector identity of equation (3-9), the right hand side of this equation can be written as —V x V x EL which is zero since EL has no curl. As well, the transverse current JT can be written as j? = —neeur with dux/dt = — eEfjm^ to zeroth order. Zeroth order is only required here since this term will later be multiplied by a first order term. Then equation (3-16) becomes < ™ - * ) & - - ^ & • (3-") cz dt2 c2me Chapter 3 Stimulated Raman Scattering 28 For the scattering of the large incident light ray off a small plasma wave, set E~r = Et + E and ne = no + rie, where E and he are first order quantities representing the scattered light wave electric field amplitude and the electron perturbation, respectively. After this linearization and again using dux/dt = —eET/me, equation (3-17) can be written as - c2V2 + ul)u3 = -—neut (3 - 18) where ut refers to the velocity of the electrons due to the incident light wave and us refers to the electron velocity due to the scattered light wave. This equation then represents the transverse current (tike right hand side) producing a scattered electromagnetic wave (the left hand side)!16'18*.' The connection between the scattered electromagnetic wave and amplification of the density perturbation® can also be shown. Since we have assumed that the ions are fixed, the electrons can h& fcueated as a warm fluid. Then the plasma fluid equations ^ + V • (netfe) = 0 (3 - 19) ^ + ue • Vue = -—(E + iT. x -) - ^ (3 - 20) at me c neme can be used, ne being (the electron density, ue the electron fluid velocity and pe the pressure. Separate the velocity into transverse and longitudinal parts, v7e = UL + U~T. The vector identity u2 ue • Vue = V(^) - ue x (V x uT) (3-21) (V x n't = 0) will assist in rewriting equation (3-20). Also note that Chapter 3 Stimulated Raman Scattering 29 where A is the vector potential, so that the term V x u~r on right hand side of equation (3-21) equals eB/mec. Finally, after substitution of all of this, equation (3-20) can be written as - W j - I B L . (3-23) at me 2 mene The V(u2) term represents the ponderomotive force'16'19', proportional to the intensity of the longitudinal and transverse components of E. To deal with the pressure term, assume the adiabatic equation of state, pe/"e = constant, so that equation (3-23) can be written as where vth = yjTjmt is the electron thermal velocity. Linearize this equation and the continuity equation by letting U"L = u"i, ne = no + ne, EL = E\ and u\j< = ut + u3 and keep terms of first order only. Equation (3-24) becomes ^ = _ V( « V tf.) - 2 & V n « • (3-25) at m no The continuity equation can be linearized to give + n0V • = 0 . (3-26) As well, the linearized Poisson's equation can be used to show that ^ V - ^ = -W?ne . (3-27) me v Equations (3-25), (3-26) and (3-27) then give (^+u2p-3v2hV2)ne = n0V2(ut-v:3) . (3-28) Chapter 3 Stimulated Raman Scattering 30 This equation shows how a fluctuation in electron density (left hand side) is produced by variations in the wave pressure (right hand side). The coupled equations (3-28) and (3-18) describe the feedback loop causing stimulated Raman scattering. 3.2 Dispersion Relation A dispersion relation can be derived for the Raman instability from these coupled equations'12,16'. Taking the incident light wave to be oscillatory, u~r = UQ COS(&O • x — u>ot), equations (3-18) and (3-28) can be Fourier transformed to give —* 4r7T G ^ (u>2 - k2c2 - u2)u3(k,ui) = u0{ne(k - k0,u; - u0) + ne(k + k0,u> + too)} (3 - 29) and (to2 - u>l)rie(k,u) = ——UQ • {us(k - k0,u - w0) + ua(k + k0,u) + ui0)} . (3-30) Solving equation (3-29) for us gives r f2ne2\uo{ne(k-ko,Lo-io0) + ne(k + ko,uj + ujo)} = ( ^ — J {u2_c2k2_u2p) • (3-31) Substituting this back into equation (3-30) and neglecting the terms with k±2ko, w±2u>o as off-resonant and therefore having little effect on the instability, one obtains the following relation tf-vD.Zpi L _ _ + L _ _ \ ( 3 _ 3 2) 4 lD(w-WD,*-ifco) D{u+w0,k + k0)) Chapter 3 Stimulated Raman Scattering 31 after the cancellation of the remaining ne terms. The function D, introduced for ease of notation, is defined as D(u, k) = u 2 - k 2 c 2 - w2 . (3 - 33) This is a dispersion relation describing stimulated Raman scattering for fixed ions. The maximum growth rate of the instability in this linear theory can be determined from equation (3-32). For the specific case of stimulated backscatter, which usually is stronger than forward scatter'12', one can neglect the upshifted light "wave as non-resonant. This is only contradicted for very long wavelength electrostatic perturbations'12'. The above dispersion relation then gives (u2-c,2){(W-W 0)2-(fc-r0)2c2-^} = ^ ^ . (3-34) Now assume that the frequency of the plasma perturbation is the natural plasma wave frequency with a smalli growth term, i.e. u> = u>e + if where 7 <C u>e. Substituting this into the above equation and keeping only the real part one finds u2k2u2 412ue(uJQ-uJe)= \ 0 , (3-35) or k u ^ l _ l _ ( 3 _ 3 6 ) 4 y U>e(U0 - We) for the maximum growth rate. Here we have also assumed that the incident and backscat-tered light waves are polarized in the same direction, as is the case in the one-dimensional simulations carried out, since the growth rate then reaches its maximum value!12'. The wavenumber of the plasma perturbation which gives a maximum growth rate can be found by noting that for maximum growth the direct backscattered light wave is Chapter 3 Stimulated Raman Scattering 32 resonant. This implies that, with w « we , (w2-Wo)-(*-*o)V-u;J = 0 . (3-37) One then finds jb=sJb0 + ^ / l - ^ + 4 _ 4 . (3-38) c y w0 u>5 wj If one can neglect 3v2hkl/u>o due to the size of wo (and for not too high a temperature), equation (3-38) gives k = k o + ^.Jl-^L (3-39) c y w0 for the natural plasma wave vector giving the maximum growth rate. Note that for densities n « ncr the wave number is k = ko and goes to k = 2ko for n <C ncr. CHAPTER 4 PARTICLE CODE THEORY The basic concept of a particle code is to move a number of representative particles self-consistently with respect to internally and externally produced forces. The simulation system is initially set up to reflect the problem at hand and then left to develop through time, hopefully producing realistic results. Discussion of the theory and operation of particle codes has been carried out elsewhere in detail'20'. An overview of the material directly applicable to the particle code EMI is given here. There are substantial, basic approximations used within particle codes. Indeed, by the very nature of time and size constraints, a particle code utilizes a much smaller number of particles than the actual plasma it is trying to represent. Typically, hundreds to thousands of particles are used to simulate a possible laboratory plasma of realistic size (~ 10 1 4 to 10 2 4 particles). Also, in the simulations presented here, only one-dimensional movement of the plasma is allowed (an extra "half-dimension" is added if external fields are present) and thus some natural two- and three-dimensional behaviour of a plasma is left out. The restriction of the simulation plasma to one-dimension, however, more easily allows the behaviour of a plasma to be modelled by a small number of simulation particles. This is due to the smaller number of particles within a Debye sphere in one-dimension'20'. Other problems of non-physical instabilities and aliasing also arise from approximations and methods used within particle codes. It is therefore necessary to realize that simulation parameters must be chosen to avoid problems, and most importantly, that simulation models are intended to produce the essence of plasma behaviour without every specific detail being correct'20'. 33 Chapter 4 Particle Code Theory 34 4.1 Computational Cycle The basic endeavor of a particle code, as already stated, is to advance particles with respect to the forces between them and applied to them. This requires the integration of the equations of motion of the particles. The integration, however, must somehow be done discretely due to the nature of computers. This is carried out by advancing time in small steps and utilizing a spatial grid upon which fields and forces are stored. The basic, self-consistent computational cycle used by EMI is as follows. The position x and the velocity v are stored for each particle. These values are then used to assign charge and current densities to the spatial grid. This requires a form of weighting, since there usually are many more particles than grid points. The charge and current densities are then used to calculate electric and magnetic fields at each spatial grid point. These fields are interpolated from the grid to a force on each particle, again using the same form of weighting. These forces are then used to advance the particles in time and the cycle repeats itself. This can be diagrammatically seen in figure 8. Integrate equations of motion, move particles. F{ —• Hi —> x,-Weighting Grid —> Particles ($,8); - Fi Integrate field equations on grid. {E,B)j «- (pj)j Weighting Particles —* Grid {x,v)i -» {p,])j Figure 8: Computational cycle in EMl'20'. Chapter 4 Particle Code Theory 35 First-order weighting is used in EMI to relate particle quantities to grid quantities and vice versa. As an example, the point charge of a particle is assigned to the nearest grid point by linear interpolation. This is equivalent to an assignment to the spatial grid of a particle of charge q, lying between the spatial grid points j and j + 1, of the form Xj+i — X{ . . qi =q Ax • ( 4 _ 1 ) and Here qj is the charge assigned to the jth grid point, Ax is the spacing between the grid points, Xj is the jth. grid point and x,- is the position of the ith particle. This weighting produces a triangular particle shape of width 2Ax '20'. 4.2 Integration of Equations of Motion The advancement cycle used in EMI requires the determination of particles' veloc-ities and positions from forces acting upon them. In particular, the following equations of motion must be solved, m f = F (4-3) where m is the mass of a particle and F is the force acting upon a particle. To integrate these equations over finite time steps, they are replaced with finite-difference equations, vt+r - vt = p T (4-5) Chapter 4 Particle Code Theory 36 = vt+1 (4-6) where r is the time step size. This method is correct to second order in the time step size'20!, and obviously becomes exact as r —• 0. To solve the above equations, it is required to stagger in time the velocity and position stored for the particles. For example, in E M I , x is stored at the full integral time step, while v is stored at the half-integral time step. In effect, the terms on the right hand side of equations (4-5) and (4-6) are centered in time with respect to the terms on the left hand side. This procedure is called the leap-frog method and is graphically represented in figure 9. As a result of this method, the initial values with which the simulation is set up with must be altered to obey this time staggering. This is done by shifting v at t = 0 to v at t = — r/2 using forces calculated at t = 0. Note also that energies calculated from v and x must then be altered to appear at the same time'20'. Time t+At/2 t + A t / 2 t+At Figure 9: Leap-Frog method of integration'20'. Chapter 4 Particle Code Theory 37 In the computational cycle just given, it is necessary to calculate the electric and magnetic fields from the charge and current densities assigned to the spatial grid points. The electrostatic field is found by utilizing a combination of the fast Fourier transform and finite-differencing schemes. Poisson's equation ox1 e0 or dEx p dx e0 (4 - 8.) are used to find the electrostatic potential <f> and electric field Ex from the charge density p. The finite-difference version of the above equations are (f>j-i - 2(f>j + <j>j+\ _ p (4-9) Ax 2 e0 { and where j again denotes the spatial grid points of separation Ax. Combining these equations and a Fourier series representation of p, <f> and Ex over a finite grid, equation (4-7) becomes'20' WO = (4-11) with M 2 K2 = k2 sin(fcAx/2) JfcAx/2 (4 - 12) The electrostatic field is then found as follows. The charge density p(x) is initially found on the spatial grid. This is then Fourier transformed to p(k) and divided by K 2 to obtain 4>(k). Units are chosen in EMI so that 6Q = 1- Finally 4>(k) is inverse transformed to 4>(x), Chapter 4 Particle Code Theory 38 which is used in the finite-difference equation (4-10) to find Ex(x). Unfortunately, the use of the spatial grid in determining the force acting on the particles causes perturbations at a wavenumber k to be coupled to perturbations at wavenumbers 27rfcn/Ax, where n is an integer. This aliasing, however, is not important if the Debye length of the simulation system is much greater the the grid spacing Ax t2 0l EMI is an electromagnetic code, so transverse fields must be calculated from the spatial grid as well. This is done by assigning a transverse current density to the spatial grid from the particles. For a polarization of external fields as given in EMI, Maxwell's equations are in Heaviside-Lorentz units. Here Ey and Bz are the transverse electric and magnetic fields, and jy is the transverse current. Adding and subtracting these equations gives where is defined as the left and right going fields in vacuum, T± = \(Ey±Bz) . (4-15) The terms with the brackets of equation (4-14) form a total derivative dT'/dt for an observer moving at the speed of light ±c. This means that equation (4-14) can be inte-grated along x = ±ct + constant (vacuum characteristics) for a grid spacing constraint Ax = cr. This constraint is used in EMI and thus equation (4-14) can be written in a finite-difference form, rHt + r t I ^ ) - T ^ ) ^ J v t ( t ^ x ± c r _ ) ( 4 _ l g ) Chapter 4 Particle Code Theory 39 The transverse current on the right hand side must be space and time centered with respect to the terms on the left hand side. This correct centering gives rise to ^ - + 1 = ^ ; - l 0'-i T 1+j+i) ( 4 - 1 7 ) as the finite-difference scheme used in E M I , where j denotes a spatial grid point and n a time step. The current j~ is defined as to be calculated from the ith particle's transverse velocity at the half integral time step, qvy(t + r/2), linearly weighted to the spatial grid from the particle's position xl(t). The transverse current jy is defined similarly, except that it is linearly weighted from the particle's position xt(t + r). These values are defined in this manner so that equation (4-17) is correct to second order in r . The stability of this method depends mainly on the way in which the transverse current is found so that it is correctly time and space centered'20!. With the inclusion of external fields, the relativistically correct equations of motion to be integrated are and dx -=v (4-19) where u is the momentum per unit mass and a transverse velocity vy is allowed. Here u = with 7 = \/l + u2/c2. As mentioned before, these equations are solved using the leap-frog method. In dealing with the v x B force, one must note that this is really just a rotation of the vector u by an angle 6 = 2 t&n-1 (qBr/2jmc). This rotation is carried out in EMI correct to second order in r. Firstly define u~ and u + as S- = t ? t _ r / 2 + ^-ET- (4 - 20) ' m 2 Chapter 4 Particle Code Theory 40 and m 2 The rotation carried out is then u + ut+T/2 = u+ - 2-ET- . (4-21) — = —^— (it + + u~) x B . (4-22) / 2fmc The equations of motion are integrated by initially finding u~ from ut_r/2 and equation (4-20). The rotation of equation (4-22) is then executed, producing u + . Lastly, ut-LT/2 is found from adding the last half of the electric field impulse to u-+ in equation (4-21). The actual scheme used in the 1 ^-dimensional code EMI is a more efficient version of this procedure developed by Buneman'21'. With the abbreviations r = - tan( | ) , (4-23) * = -™d = rh*- ' (4_24) and c = cos 6 = 1 + r the sequence for the rotation becomes c s = \—4 , (4-25) ux = ux + uy r (4 - 26) u+ = u~ - u'x s (4 - 27) and u+ = u'x + u+r . (4-28) Chapter 4 Particle Code Theory 41 The position of a particle is advanced as in equation (4-19), x t + r = xt + vXtt+z T = xt + '—T . (4-29) 7t+r/2 In this way, the equations of motion of the particles are solved using the leap-frog method. 4.3 Loading the Particle Distribution EMI is a code with periodic boundary conditions modelling an infinite plasma. Initial conditions, especially with respect to the establishment of the initial distribution f(x,v,t — 0), are very important in such a code since they determine the later behaviour of the system. EMI uses two methods to load particles into an initial distribution. For the case of non-relativistic temperatures, the distribution / in EMI is a Maxwellian, normalized so that J fdv = N (4-30) where N is the number of particles of the species being loaded. To load the particles in the simulation, the integral JQV f dv is evaluated. Every time the integral's value increases by one, that value of v is assigned to a particle. This integral is carried out from 0 to v until an upper maximum of vmax = 4vt/, is reached (essentially oo). Since this integral's lower limit is 0 and not —oo, only half of the particles are assigned a velocity in this manner. The remaining particles are then given the negative of these velocities. The vx and vy distributions are determined independently and the velocities of randomly selected particles are exchanged to discourage correlations. Chapter 4 Particle Code Theory 42 For relativistic temperatures, the distribution function in EMI is loaded using a relativistic two-dimensional Maxwellian of the form J [ P ) 2TT T (1 + T/mc2) K ' where p = f\v\/c. This distribution is normalized so that 27T +00 J d9 J fpdp = N , (4-32) 0 0 where the angle 6 is defined by px = pcosO, py = psinfl. To load the distribution, the integral J0P f pdp is done numerically, the result being multiplied by 7r/2. The maximum value of this integral is then jiV, the number of particles in the first quadrant. Every time 7r/2 J0P f pdp increases by one, four particles are assigned a velocity from p and an angle 9,9 + 7r/2, 0 + 7r, 6 + 3TT/2, respectively. Here 9 is a random angle between 0 and 7r/2. The actual particle velocity is found from - - p C O S < ? v9 = J^L . (4-33) V T + P2 v / T T p 1 Other modifications are necessary to load a satisfactory Maxwellian. For either non-relativistic or relativistic temperatures it may be desired to add a possibly relativistic drift velocity to the initial distribution. This is done using'22'20' V*-l + vxv0/c> l + ^ 0 / c2 V y ( 4 3 4 ) where vQ is the drift velocity to be added to the distribution. Lastly, after the initial distribution is loaded, whether relativistic or not, a random exchange of particle positions is carried out. This is done to reduce correlations between particles, since velocities are assigned to the particles in an ordered sense. Chapter 4 Particle Code Theory 43 The above distributions may be loaded using either a random or a quiet start. The random start utilizes a uniform random number set. These random numbers are used to determine the random angle 9 and to randomly exchange particle positions to discourage correlations. Random numbers, however, produce bunching and correlations in phase space'20!. These fluctuations can be much larger than regular experimental levels, and this can alter the behaviour of sensitive systems to be simulated. Indeed, on a rough scale the variations in density may be correct but finer details can be left out. The quiet start attempts to cure this problem by loading particles into phase space as smoothly as possible. This method assigns a 9 and a position x to the particles using mixed-radix digit reversal. For radix-two digit reversal, used to scramble the positions of the particles, a fraction between 0 and 1 is found by mirroring base two numbers and converting them to a decimal fraction as shown in table 1. Radix-Two Digit Reversal Base Two Number Mirrored Fraction Decimal Equivalent 10 101 0.01 0.101 1/4 = 0.25 1/2 + 1/8 = 5/8 = 0.625 Table 1: Radix-2 digit reversal. The velocity angles 9 are scrambled using radix-three digit reversal, as shown in table 2. Radix-Three Digit Reversal Base Three Number Mirrored Fraction Decimal Equivalent * 1 12 0.1 0.21 1/3 = 0.33... 1/3 + 1/9 = 7/9 = 0.77... Table 2: Radix-3 digit reversal. This method allows for loading an initial distribution function with noise levels much less than those of a random start. CHAPTER 5 THE PROGRAM EMI As already stated, E M I is a fully relativistic, l|-dimensional electromagnetic par-ticle code with periodic boundary conditions. The code allows for up to two species of particles. Both the initialization of the simulation parameters and the computational cy-cle are combined into one program. Major regions of calculation are, however, subdivided into a number of subroutines. 5.1 Initialization At the beginning of E M I , parameters to be used in the simulation are calculated. Initial basic parameters are input from a file EMI.INPUT. A table of these input param-eters along with their default values is given on the next page. The subroutine HISTRY, to be described later, is called at the beginning of the program to zero matrices used in the computational cycle. Fundamental values used in the simulation are then calculated from the input parameters. The system's length L is found from L = 2TT/KO, KO being the system's fundamental wave number. From this and the number of spatial grid points, NG, the grid point spacing is calculated, A x = L/NG. The cyclotron period is calculated from the cyclotron frequency in the case of an applied static magnetic field. The charge per particle of each species, Q, is found from the supplied plasma frequency WI and charge/mass ratio, QM. That is, Q — u>p/Ns QM, where Na is the number of particles of that species and units are chosen so that eo = 1. As well, the density perturbation amplitude for each species is normalized by multiplying by L/NP, NP being the total number of particles. 44 Chapter 5 The Program EMI 45 EMI INPUT PARAMETERS Variable Description Default Value Nl Number of particles of species 1 128 N2 Number of particles of species 2 128 NT Number of time steps 400 NG Number of grid points 32 MODE Mode to perturb 1 IMODE Mode of initial transverse field 0 Xll Amplitude of perturbation in species 1 0.0 X21 Amplitude of perturbation in species 2 0.0 V10 Drift velocity of species 1 / c 0.0 V20 Drift velocity of species 2 / c 0.0 WI Plasma frequency of species 1 1.0 W2 Plasma frequency of species 2 1.0 WC1 Cyclotron frequency of species 1 0.0 QM1 Charge to mass ratio of species 1 -1.0 QM2 Charge to mass ratio of species 2 -1.0 KO System fundamental mode 0.5 DT Time step size 0.1 TEMPI Temperature of species 1 0.0 TEMP2 Temperature of species 2 0.0 IREAD Start from the IREAD'th save 0 ISAVE Save every ISAVE steps 0 TSTART Starting time 0.0 NKAYS Number of electrostatic modes to save 8 KAY Matrix containing modes to save 2,3,4,5,6,7,8,9 Al Compensating factor 0.0 A2 Smoothing factor 0.0 ISD Random or quiet start 1 (Random) EO Amplitude of radiation fields in the plasma 0.0 Table 3: EMI input parameters. Chapter 5 The Program EMI 46 EMI INPUT PARAMETERS Variable Description Default Value IRHO Save RHO every IRHO time steps 128 IRHOS Save smoothed RHO every IRHOS time steps 128 IPHI Save PHI every IPHI time steps 128 IXVX Save X and VX data every IXVX time steps 128 IVXVY Save VX and VY data every IVXVY time steps 0 IEX Save EX every IEX time steps 128 IEY Save EY every IEY time steps 128 IBZ Save BZ every IBZ time steps 128 IJY Save JY every IJY time steps 128 IEYL Save left going fields every IEYL time steps 128 IEYR Save right going fields every IEYR time steps 128 ITHERM Save temperature every ITHERM time steps 128 IFT Save FT. of radiation fields every IFT time steps 128 NTH Call HISTRY. every NTH time steps 128 IH Save Mode energies every IH time steps 128 Table 3 (continued) After the determination of the basic simulation parameters, the loading of the par-ticles is then carried out. For a quiet start, the particles' positions are found from mixed-radix digit reversal. The fraction from this digit reversal is generated by the function RADREV. For a random start, the positions are initially set to equal divisions along the length of the plasma. Perturbations of the particle density are carried out here by alter-ing particle positions. The following subroutines are then called to finish establishing the simulation plasma at time zero. Subroutine CREATR The subroutine CREATR loads the initial velocity distribution of the particles using either a random or a quiet start. This is done, as described in Chapter 5, by integrating Chapter 5 The Program EMI 47 f dv. The quantity assigned to each particle is actually v/c, the reason for which will be explained later in the discussion of the subroutine SETRHO. For relativistic temperatures, a velocity angle is assigned to each particle using either a random number set or mixed-radix digit reversal. Any drift velocity as input from EMI.INPUT is also added to each species. Subroutine SMEAR This subroutine is only called in the case of a random start. It de-correlates the positions of the particles by randomly exchanging their positions. At this point, any initial radiation fields within the plasma are loaded into the simulation. Subroutine SETRHO SETRHO calculates the initial accumulated charge density at the spatial grid points and enforces periodic boundary conditions on the particle positions. As well, this sub-routine converts the particle positions x to a;/Ax. This normalization accompanies the assignment of v/c = VT/AX (AX = cr) to each particle in CREATR. These basic vari-ables are normalized in this fashion to reduce the number of calculation steps used in advancing particle in time. Then, (5-1) is the new advancement of the particle velocities, a being the acceleration, and (5-2) is the new advancement calculation for the particle positions. In this way, the two multi-plications by r necessary per particle per time step are avoided. Chapter 5 The Program EMI 48 Subroutine SETV This subroutine shifts v backwards half, a time step so that x and v are staggered in time. This allows leap-frog integration of the equations of motion, as explained previously. The fields acting on the particles at t = 0 are found from a call to the subroutine FIELDS (explained later). These fields are then used to find the initial force on the particles, from which the backwards shift in v is done. This backwards shift is a combination of a backwards half step rotation due to Bz and a reverse half step acceleration due to E. This shift in v is carried out by a call to the subroutine ACCEL, to be explained later. Restart It is important to note that EMI has a restart option, that is, the simulation can begin where a previous simulation run ended. In this case a file of data from the previous run is read, giving all the input parameters and the particles' positions and distribution function. All of the above subroutines are then skipped in such a restart. 5.2 The Simulation Loop After initialization of the simulation plasma, particles are advanced in the com-putational cycle. This cycle is iterated through the number of times requested in the initialization procedure. The cycle consists of the subroutines ACCEL, HISTRY, MOVE and FIELDS. Subroutine ACCEL This subroutine's main purpose is to advance the particles' velocities by one time step. The acceleration is normalized here as described before to avoid unnecessary multi-plications by r. The quick rotation scheme, as described in Chapter 5, is used to deal with Chapter 5 The Program EMI 49 any magnetic field existing in the plasma. As well, particle kinetic energies are calculated here using TXX Kinetic Energy = — v0idVnew . (5 — 3) Note that only species of particles with a large enough charge to mass ratio are treated relativistically in ACCEL, heavier particles such as ions being treated non-relativistically. Subroutine MOVE The other half of the leap-frog integration, that of advancing the particles' positions, is completed within this subroutine. Periodic boundary conditions are enforced, so that positions xnew < 0 or xnew > L are replaced with xnew + L or xnew — L respectively. As well, the charge and current densities are accumulated on the spatial grid using linear weighting. Subroutine HISTRY HISTRY saves diagnostic results from the simulation. It is called in the initialization procedure at rj = 0 to zero the diagnostic matrices used within the code. Data such as particle kinetic energies and field energies are then stored in these matrices during a simulation run. The matrices are then saved and emptied every time HISTRY is called, ready to collect the next sequence of diagnostic data. In this way, most output is done only every few hundred time steps. Subroutine FIELDS Calculation of the field quantities from the charge and current densities stored on the spatial grid is carried out by this subroutine. This includes finding the electrostatic potential <j> and electric field Ex. As well, Maxwell's equations are integrated utilizing the transverse current stored on the spatial grid. The field energies per mode are also Chapter 5 The Program EMI 50 calculated. The electrostatic energy is found from 2 kmax Electrostatic Energy = — ^ Pk4>k (5 — 4) k=K0 utilizing the Parseval equality'20'. This overestimates the energy at short wavelengths, while the other possible calculation Ylk-El/L would underestimate the energy at short wavelengths. Expression (5-4) is used to calculate the electrostatic energy since Coulomb's law (p —> Ex) is modified by the use of the spatial grid'20'. After completion of the number of iterations requested of the computational cycle, HISTRY is called again to store the final diagnostic data. In some of the above subrou-tines, a fast Fourier transform is required. This is done by a set of efficient subroutines. CPFT transforms a(x) + ib(x) —> A(k) + iB(k). To gain speed in transforming real quan-tities used in EMI, a pair of real sequences a(x) and b(x) are set up and transformed together. This is done by RPFT2 and the inverse by RPFTI2. The succession of calls RPFT2, CPFT then constitute a fast Fourier transform and CPFT, RPFTI2 the inverse transform. A listing of the particle code EMI can be obtained from C. K. Birdsall or Bruce I. Cohen (University of California, Berkley). Since EMI is fairly well known and available its listing was not included here. CHAPTER 6 SCHEME FOR THE VLASOV ALGORITHM 6.1 Introduction In deriving an algorithm for simulating a plasma using the Vlasov equation, a dif-ficulty arises from the fact that computers are discrete in nature. This means that the continuous behaviour of a plasma as described by the Vlasov equation must somehow be simulated discretely. The method of fractional steps developed by N. N. Yanenko'23'was utilized to accomplish this. The plasma in question is considered to be one-dimensional along x and infinite in extent so that periodic boundary conditions can be used. However, since electromagnetic radiation will be treated in the simulations, other dimensions will be included in the algorithm. The plasma will then be allowed to acquire a transverse velocity in the y direction. We also assume that, for the moment, any phenomena in the plasma is of high frequency so that the ions form a fixed, uniform positive background. The plasma has an unperturbed electron density ne. With respect to the Vlasov equation, we are only concerned with the electron distribution function fe(x,ux,t), where the subscript e will be dropped from now on for simplicity. For convenience in relativistic cases, ux will be taken as the momentum per unit mass. Then the velocities vx and vy can be found from the fact that'22! u Uy = (6-1) 51 Chapter 6 Scheme for the Vlasov Algorithm 52 giving vx= , U x vy= , U y (6-2) where u2=ul + ul . (6-3) In the relativistic case, fo, the unperturbed distribution function, is given by the relativistic Maxwellian / o ( u « ) = 7 e x p + «»/<:*) . (6-4) Here mc2 * = ^ (6-5) and A = 2eqK1{q) (6-6) where K\ is the modified Bessel function. The factor eq is kept in the numerator and denominator to keep values finite for calculation. In the non-relativistic case vx = ux and fo is given by the usual non-relativistic Maxwellian distribution. The normalization of fo in either case is chosen so that J fQdux = 1 . (6-7) The Vlasov equation describes the evolution of the distribution function / and is given by df df df ( ( { — -vr- a r - — . ( b - 8 ) dt dx du X The assumption in using the Vlasov equation is that collisions in the plasma are negligible. If no external fields are incident on the plasma the acceleration ax (here the change in Chapter 6 Scheme for the Vlasov Algorithm 53 momentum per unit mass in time) arises from longitudinal electric forces obeying ^ - 4 x , (6-9) in c.g.s. units. Here Ex is the longitudinal electric field and p is given by p = nec(l - J f dux) (6-10) for a neutralizing positive background charge which we are assuming. If externally applied fields are present, the acceleration ax also contains the Lorentz term, giving ax = ~(Ex+vyBt/c) (6-11) and one must include the velocity in the y direction, with an equation of motion for uy given by ^l = -±(Ey-vxBz/c) . (6-12) The evolution of the electromagnetic fields are given by Maxwell's equations dBz 3Ey dEy dBz ~W = " C & f ' -to = ~C-dx- ~ 4*'» ( 6 " 1 3 ) where the transverse current jy is given by jy(x,t) = -nee J duxf(x,ux,t)vy(x,ux,t) . (6-14) It is important to note the significant difference between the coordinates x, ux, t involved in / and the rest of the variables described here. The coordinates of / are simply coordinates (not variables) and are independent of each other. That is, they are chosen in the problem and thus are not functions of each other. A result of this is that any derivative of these coordinates with respect to one another is zero. Their choice, however, Chapter 6 Scheme for the Vlasov Algorithm 54 influences other variables that depend on them. Visually one can imagine that x, ux and t form a three-dimensional grid in which / is advanced. 6.2 Non-Relativistic Landau Damping Now consider a plasma with no external fields in the non-relativistic limit. A per-turbation in the electron distribution function will produce Landau damping. This phe-nomenon should be seen in the simulations under these conditions and the results can be compared with theory and the particle code EMI. The algorithm developed from the method of fraction steps by Cheng and Knorr'24' is used here in the non-relativistic case to advance / in time. As will be shown, this method is correct to second order in the time step size. We assume that the distribution function is known after n time steps, each time step of length r, fn(x,ux) = f(x,ux,nr). (6-15) To see how this method works we will look at the advancement of / by one time step, that is, how fn+1(x,ux) = f(x,ux,nr + r) (6-16) is obtained. The first step in the algorithm shifts / along x at the one half time step, f'(x,ux) - fn(x - vx^,ux) . (6-17) Then / is shifted along ux at the full time step, f"(x,ux) = f'(x,ux - CLXT) (6-18) Chapter 6 Scheme for the Vlasov Algorithm 55 but ax is the acceleration evaluated at the half time step, a'x = ax(x,tn + ^) = ax + ( 6 - 1 9 ) where a Taylor series expansion to first order in r has been carried out. Lastly, the distribution function for the next time step is found by completing the shift in x, II T fn+1(x,ux) = f (x - vx-,ux) . (6-20) Recalling that we are dealing with a non-relativistic situation where ux = vx, we can check that this procedure is correct to second order in the time step size r. Substituting equation (6-18) into equation (6-20) gives fn+1(x,ux) = / (x - (vx - axr)^,ux - axr) . (6-21) Now here the acceleration ax must be calculated at x — VXT/2 due to the shift in x. To first order in r this is 1 i T \ daxr daxr ax(x - vx-) = ax + — - - V X - Q - - (6 - 22) so that equation (6-21) becomes / T i T2 fn+i(x,ux) = / (x - vx- + ax — ,ux - 8ux) (6 - 23) where 8ux is given by dax dax r2 8ux = axr -| 5 — vx——— . (6 — 24) 3x 2 v ; Chapter 6 Scheme for the Vlasov Algorithm 56 Finally, using equation (6-17), one arrives at r2 fn+i(x,ux) = fn(x - VXT - ax — ,ux + 6ux) . (6-25) This is a statement of how the algorithm advances /. We can check how correct the algorithm is by Taylor expanding equation (6-25), keeping terms up to second order in r r2 fn+i - /« = TSI + —S2 • (6 - 26) Here S\ is the coefficient of r in the expansion, Si =-vx- ax-— . (6-27) In evaluating the coefficient of r2/2, one must remember that x, ux (here = vx) and t are coordinates and not dependent variables. As well, for the case of Landau damping where no external fields are applied, ax = —e/mEx which is independent of ux [ux is integrated out in finding Ex). One then obtains _v2&J_ 2v d2/ q 2a2/ ^ daLdf__da±df_ q df_ _ 2 x dx2 x xduxdx xdu2 x dx dux dt dux z dx Now compare these coefficients with those obtained from a Taylor expansion to second order in r of the Vlasov equation itself, df d2fr2 r2 /•«-- f - i T + #T- T' r + I i T • ( 6 " 2 9 ) From the non-relativistic Vlasov equation one obtains Chapter 6 Scheme for the Vlasov Algorithm 57 so that the first order coefficients are equal. Using some of the same reasoning as for finding S2, one obtains the second order terms: 2 x dxdt xduxdt dt dux Utilizing the Vlasov equation again one finds id2f n d2f 2d2f dax df dax df . df ID o n , T 2 = v 2 x ^ + 2 v x a x ^ + a 2 x ^ + v x ^ - ^ + ax/x . (6-32) This agrees with the second order coefficient Si and thus the algorithm is correct to second order in the time step size. 6.3 Relativistic Landau Damping When including relativistic effects, one no longer has ux = vx, but rather the relation given in equation (6-2) . To concentrate on the changes included in the relativistic treatment, we will assume no external fields are present. Landau damping then results from a perturbation in /. Again we assume we know the distribution function after n time steps fn(x,ux) = f{x,ux,nr) (6-33) and we want the distribution function at the next time step, fn+i{x,ux) = f(x,ux,nT + r) . (6-34) The algorithm chosen for advancing / in the relativistic case is the same as for the non-relativistic case, and will now be shown to be correct to second order in r. Substituting backwards, we obtain Chapter 6 Scheme for the Vlasov Algorithm 58 » r2 fn+1(x,ux) = fn(x - VXT - ax — ,ux - Su) (6 - 35) where Su is the same as in the non-relativistic case. Here the acceleration ax must compensate for the fact that the acceleration used in the ux shift is different from the one used to shift vx in the x shift, due to the relativistic effects. Thus ax is given by ' dvx dvx . . to zeroth order in r (only zeroth order is necessary since ax will be multiplied by r2). We can then Taylor expand equation (6-35) up to second order in r, / n + l - / n = 5 1T + y 52 (6-37) where S\ is Vxlh~axdu~ (6-38) and the coefficient of r2/2, is given by S<-v2d2f 1 2V , J , ^ I a 2 9 ' 1 I vxdQx df dQx df I a x d V x d f . (6-39) 2 x dx2 1 xduxdx x du2 x dx dux dt dux xduxdx We must now compare this with the result from the Vlasov equation itself. The time dependence of the Vlasov equation as expanded in a Taylor series to second order in T is df d2fr2 r2 From the Vlasov equation, rp df df = ~Vx~dx ~ axdux~ ' ( 6 " 4 1 ) Chapter 6 Scheme for the Vlasov Algorithm 59 and the first order coefficients agree. Using the Vlasov equation twice and remembering that ux is not equal to vx, the second order term becomes = + 2v^dxlu-x +a*dvJx+ V^du~x ~ -bTduZ + axdv7xTx • (6~42) Thus, since the second order terms agree, the algorithm is correct to second order in time step size for relativistic cases as well. 6.4 How to Calculate the Longitudinal Field The acceleration used in this previous electrostatic algorithm is evaluated at the half integral time step. This leaves the problem of how to calculate the longitudinal field at the half time step a«(Mn + = -—£*(*»» + £) • (6-43) Z rn Z This longitudinal electric field, E'x, is calculated by integrating the charge density p at the half time step p = ne|l - J duxf(x,ux,tn + ^ ) | . (6-44) The distribution function evaluated at the half time step can be written, to first order in r, as T T T f(x,ux,tn + -) = fn(x - vx-,ux - ax-) . (6-45) The charge density can then be written as (6 - 46) Chapter 6 Scheme for the Vlasov Algorithm 60 Since this integral is taken over the entire ux axis, a shift can be made in ux. This is equivalent to a shift in the ux origin, which leaves the integral unchanged. Shifting the ux origin by +axr/2, the above integral becomes p = ne|l - Jduxf'(x,ux)^ (6-47) where /' is the distribution function given by equation (6-17) in the time evolution algo-rithm. The longitudinal electric field should therefore be computed after the first shift in x is carried out in the algorithm. This then presents a way of shortening the algorithm in practice. We can store the shifted distribution function / and shift x by vxr, combining the initial and final shifts in x into one shift. The stored distribution function / can then be used to calculate the longitudinal electric field for the.next ux shift. 6.5 Including External fields Applying external fields to the simulation plasma requires a modification of the Vlasov algorithm. Assume that the incident electromagnetic radiation is polarized so that we have an Ey and a Bz. The plasma should be treated relativistically to be able to simulate large incident fields such as lasers. Then, as stated before, the acceleration in the one-dimensional Vlasov equation is given by ax = --{Ex+vy^) . (6-48) m c The electric field component Ey in the applied field causes the plasma to experience a transverse velocity vy which is given by (6-2). Note that this velocity and the momentum Chapter 6 Scheme for the Vlasov Algorithm 61 per unit mass uy depend on the coordinates and t. Now, vx = Ux (6 - 49) - (ul + u})/c* implies that vx depends partly on uy, so that unlike before, vx is a function of the coordinates XL % j 3? 3-11 d t. This dependence gives rise to extra terms in the Taylor expansion of the Vlasov equation. As before, we start with the Vlasov equation d* df df (fi ^ Tt=-Vx^-axdu-x (6"50) and Taylor expand the time dependence to second order in the time step size, a/ r2a2/ / n + 1 - / n = r - + Y — . (6-51) Using the Vlasov equation, the first order term remains correct, so that the applied fields only cause a second order correction to be necessary in the algorithm. Therefore look at the second order term, the coefficient of T2/2, a 2 / 9 dl df dt2 ~~dt[Vxdx +axduJ (6 - 52 = _dv±dl_v ^L_dax_df__a d2f dt dx x dxdt dt dux x duxdt Re-substituting the Vlasov equation in two of the terms allows their expansion, &f_=d_,_v dL _a df_\ dxdt dx x dx x dux 6 - 53 dvx df a2/ dax df d2f = vx-z—-r a r -dx dx dx2 dx dux dxdux Chapter 6 Scheme for the Vlasov Algorithm 6 2 and a 2 / = d d£ _a df_ duxdt dux 1 dx 1 dux dvxdf d2f _ 3ax df d2f duxdx 1 duxdx duxdux 2' du\ (6 - 54) Therefore, we obtain dt2 2 dx dt dx xduxdux (6 - 55) where T2 is given by equation (6-42) and represents the expansion in absence of external fields. The modifying terms of equation (6-55) can be represented more easily by using equation (6-1), giving dvx _ uxuy duy (6-56) dx c2F dx and dvx uxuy duy ~dt = ~ c2F dt (6-57) and using equation (6-11), dax du2 e Bz 1 m c F ul, du (l + 7 f ) c2 du, (6 - 58) where u 2\ 3/2 (6 - 59) Equation (6-57) can be expressed more basically by utilizing the equation of motion for uy (equation (6-12)), giving dv euxuy Bz at m cLb c (6 - 60) Chapter 6 Scheme for the Vlasov Algorithm 63 We can now consider how to modify the Vlasov algorithm. Recall that the correction terms are coefficients of r2/2. Thus they need only be correct to zeroth order in r for the algorithm to be correct to second order in the time step size. First shift x by 6x/2 with dvx dvx T2 8 x = -V x T + K_ _ _ _ ) _ . (6-61) This gives a modified distribution function /'. In dealing with the next shift in ux, one must look at how the longitudinal field will be calculated, which utilizes the integral I = J duxfn(x - vx^,ux - ax^) . (6-62) The problem involved with this is the fact that ax = --(Ex + u„ — ) . (6-63) m c Indeed, since vy (and therefore ax) now depends on ux, one cannot shift the ux origin by axr/2 without affecting the integral /. Using / as before is thus no longer correct to second order. Note that Ex does not depend on ux (ux is integrated out in finding Ex), so that the only modification comes from vy. We need to advance ux in / to the half time step used in I. Therefore, first partially shift ux by 6\ux, where 8\ux = — vy — ^ . (6-64) m c 2 This gives an intermediate distribution function/* = f(x + 6x/2,ux+6iux) which can be used to obtain a correct value to the above integral. This value of i" is correct to second order since the vy part of the acceleration has been accounted for up to the time of the Chapter 6 Scheme for the Vlasov Algorithm 64 integral, r/2. As before, the shift by Ex is just a shift in the ux origin and does not change the value of the integral, so that 1 = J duxf*(x,ux) . (6-65) This integral can then be used to find the longitudinal field Ex. After this, shift ux by $2UX = ^{Ex + vy^)r + ax^Y (6-66) m which completes the shift in ux. Finally shift x by 6x/2 again to complete advancing / by one time step. In practice this last shift in x is combined with the first shift in a;, as alluded to before. As well, if the distribution function /* is stored instead of / , then SiUx and ^ i i i can be combined into one shift in ux. C H A P T E R 7 T H E V L A S O V C O D E 7.1 Introduction The Vlasov code presented here is based upon the algorithm of the previous chapter. It is a code with periodic spatial boundary conditions which allows external fields Ey and Bz and is correct to second order in the time step size. The restriction of mobile electrons in a neutralizing background of fixed ions is dropped with the introduction of a distribution function for the ions themselves. This second distribution is advanced using the algorithm in the previous chapter, but it introduces some modifications to calculation of the charge and current densities. One major difference to note between the Vlasov code and particle codes is that there are no positions of particles to keep. In fact, there are no particles. The distribution function / is simply defined over an x-grid (/periodic over x) and a ux-grid (/ assumed to be zero at its endpoints). The Vlasov code then simply shifts the arguments of f(x, ux) as alluded to in the previous chapter. The algorithm for integrating Maxwell's equations by Birdsall and Langdon'20' was attempted but was found to be unstable for the code. All the requirements of the algorithm were met, including the correct averaging of the transverse current jy, except that jy could not be linearly weighted over any "particle" positions. The stability of the advective differencing scheme used in EMI and E M 1 B is attributed mostly to the unique and correct determination of the of the transverse current, time and space centered'20', so this 65 Chapter 7 The Vlasov Code 6 6 deviation is the most likely reason for its failure. Another algorithm was then used, to be explained below. One of the first considerations in developing an algorithm for advancing the variables in a code is to determine at what time step each variable should be stored at. Initially, the variables can be set up as in figure 10, with ux at t = 0, uy, Ex, Bz and /* at the half time step t = 1/2, and Ey at the full time step t = 1. Then ux can be advanced (the component of f*(x,ux), not the ux-grid points), correct to second order in time step size, using since ux is defined at the full time step and Ex, vy (calculated from t t „ ) , and Bx are defined at the half time step. Bz is then advanced using Maxwell's equations a ,x = -±(Ex+vy^) m c (7-1) dBz dt-(7-2) with Ey at the full time step. To advance uy we need to use Bz old Ux Ey -1/2 0 1/2 1 TIME Figure 10: Initial variable placement in time. Chapter 7 The Vlasov Code 67 where B* is the magnetic field at the full time step. A quadratic fit in time is used to determine a value of B* correct to second order (see equation (7-23)). Lastly Ey is advanced using dEy dBz dt dx - 4 * j f • ( 7 - 4 ) with Bz and jy (calculated from /*, uy) at the half time step. A time step advancement is graphically displayed in figure 11. It is important to note that the derivatives dEy/dx and dBz/dx are found at the same x-grid spacing as Ey and Bz, respectively, since cubic splines are used to find the derivatives of the functions. This is more convenient than in particle codes where derivatives are found at x + Ax/2, using the last and present positions of the function in question. As was already noted, the Vlasov code allows up to two distribution functions (i.e. electrons and ions). This requires saving each distribution on its own Uj-grid and saving the transverse velocity for each species. The distribution function fe is always dealt with relativistically, while fi is only dealt with relativistically if the mass ratio mi/me between the species is less than 100. In a non-relativistic treatment, the ux dependence of vx and vy is dropped. When both distribution functions are used in the code, the method of calculating the charge density and the transverse current is modified. In this case, p = nQe J(fi-fe)dux ( 7 - 5 ) and jy =n0e J(fiVVli- fevy,e)dux , ( 7 - 6 ) a generalization of the method used in the previous chapter. Chapter 7 The Vlasov Code t t+T/2 t+T t+3772 t+27" t t + T / ? t+T X+3T/2 t + 2 T Fig 11: A Vlasov code time step advancement. Chapter 7 The Vlasov Code 69 A listing of the Vlasov code can be found in appendix B. The variables used in the code are defined in the listing of COMM (comments on Vlas). VINC, MINC and CINC are include files which define the variables for the compiler. The code itself is separated into two distinct parts. Firstly, an initialization procedure is run to calculate the parameters of the algorithm and to set up the variables at the correct time positionings. These values are then stored. The main body of the program then reads these stored initial values and iterates through the requested number of time steps of the algorithm. 7.2 Initialization Procedure The main body of the initialization section of the Vlasov code is VINI. The initial default values in electrostatic units of the parameters are defined at the beginning of VINI. The subroutine INPT is then called, which inputs the initial values of some of the parameters. These parameters define the system to be studied and are explained in table 4. Any entry of zero substitutes the default value of the corresponding variable. Basic plasma parameters are then calculated. The length of the plasma slab XMX is calculated from a wavelength WVL= 1 x IO-4, the number of wavelengths in the slab XWL, and the ratio of the plasma density to the critical density DEN; The x-grid spacing DLX is then found from DLX=XMX/NDX, where NDX is the number of x-grid points. The time step size DLT is found from XMX = WVL * XWL (7-7) VI - DEN DLT = DLX CTX (7-8) VOL Chapter 7 The Vlasov Code 70 where VOL is the speed of light and CTX is the Courant factor, usually set equal to 0.5 for convergence. The minimum UMN and maximum UMX velocities for integration purposes are found (usually ±5x the thermal velocity), from which each i^ -grid's spacing is found, D U I . E M ™ ( 7 _ 9 ) where NDU is the number of ux-gvid points and DUI is the grid spacing for the distribution function fi. The frequency and amplitude of a laser field are also calculated, utilizing the fact that B = y/l — ne/nc E. At this time, coefficients for the cubic spline interpolation procedure are determined by calling the subroutine SPQR. INPUT PARAMETERS Variable Description Default Value NDX Number of x-grid points 16 NDU Number of ux-grid points 51 MIE Ion/electron mass ratio -1.0 ISPEC Number of species (1 or 2) 1 UT1 Thermal velocity (species 1) / c 0.20 UT2 Thermal velocity (species 2)/ c 0.20 DEN Density of electrons / critical density 0.360 XWL Slab thickness in laser wavelengths 0.533333 VOC Electron quiver velocity / c lO"16 NPA Number of perturbations in the plasma slab 1 PAME Amplitude of the plasma perturbations (electrons) 0.010 PAMI Amplitude of the plasma perturbations (ions) 0.010 V01 Drift velocity (electrons) 0.0 V02 Drift velocity (ions) 0.0 NFX Number of Fourier components to save in EX 2 NFY Number of Fourier components to save in EY 2 ITF Store distribution function(s) every ITF steps 1000 IDUMP Store diagnostic data every IDUMP steps 100 Table 4: Vlasov input parameters. Chapter 7 The Vlasov Code 71 Once these initial values are determined, fe(x,ux) and fi(x,ux) are set up on the x-and ux-grids by calling the subroutines INIFE and INIFI respectively. If only one species is requested, only INIFE is called. These subroutines calculate the values of a relativistic Maxwellian at the ux grid points. These distribution functions are then normalized using the modified Bessel function so that J feo dux = 1 , J fi0 dux = 1 . (7 - 10) The modified Bessel function K\(q) is calculated using the subroutine BEK1. The nor-malization fo = o Q \ f x e xP(g /i / 9 , TfT=f) ( 7 - u ) 2e«.Ki(g) V 1 - (vx±v0)2/c2 as given in the previous chapter is used to calculate the Maxwellian distributions, with q = mc2 /ksT and ±VQ being a drift velocity. Any requested sinusoidal perturbation in the x component of the distributions is then added before returning to VINI. After establishing fe and fi on a grid system, each distribution's x component is shifted a half time step forward. This is done to prepare the distribution functions for the combination of shifts in the efficient version of the vlasov algorithm. At this point, any radiation fields are loaded into the simulation. Ey is loaded at a full time step and Bz at the half time step. To use a quadratic fit in integrating Maxwell's equations, Bz is also stored at a half time step in reverse. A transverse velocity is also loaded consistently with respect to Ey. A shift in ux is then carried out, 8\UX, by calling the subroutine VLSVU1. This is done along with the previous x shift so that the distributions functions are stored at the correct points in time (/* in the previous chapter). This allows the necessary integrations to be carried out and the efficient version of the Vlasov algorithm to be used. After this, the subroutine LONG or LONGONE (to be described later) is Chapter 7 The Vlasov Code 72 called, which calculates the longitudinal electric field at the half time step. Lastly, STOR is called to store the description of the initial simulation state created by VINI. 7.3 The Simulation Loop The main program containing the simulation loop is VRUN. This program calls a number of subroutines to advance the distribution functions the number of time steps requested. Initially the subroutine RECL is called to establish the initial state of the system as set up by VINI. As well, initial parameter values pertainent to each diagnostic are stored in respective data files. The iteration loop then begins, consisting of the subroutines VLSV, VLSVI, VLSVION, LONG, LONGONE, MXWL and SAFT. Inside of this time step loop, the distribution functions, field values and transverse velocities are stored every certain number of time steps. After this, STOR is called to save the final state of the system. The storage of this final state enables the program to be restarted from where a previous simulation ended. Subroutine VLSV This subroutine relativistically shifts the electron distribution function in both its x and ux components. Since /* is stored instead of /, these shifts are the combination of the two x shifts and the two ux shifts in the algorithm of the previous chapter. The shifts are carried out using cubic spline routines. The shift in x is done using SHIP, which shifts the x component of the distribution function assuming periodic boundary conditions. Similarly, the shift in ux is done using the subroutine SHIF, which shifts the ux component of the distribution function. This shift assumes that this ux component is zero at each end of the distribution function. Subroutine VLSVI carries out the same shift in the x and ux components of the sec-ond species' distribution function if the mass ratio between this species and the electrons Chapter 7 The Vlasov Code 73 is less than 100. If the second species is more than 100 times as heavy as the electrons (and then can suitably be called ions), the subroutine VLSVION is called instead of VLSVI. VLSVION shifts the x and ux components of the ion distribution function non-relativistically. This treatment neglects the higher order corrections to the shifting as given in the previous chapter. Subroutine L O N G This subroutine calculates the longitudinal electric field and the transverse current by integration. As stated before, the charge density is calculated from (7-12) The longitudinal electric field Ex can then be found from Ex=4n f pdx . (7 - 13) Jo As well, the transverse current is found from jy = nee J(fivy>i - feVy>e)dux . (7 - 14) The integrals (7-12) and (7-14) are carried out using Simpson's rule with an end correction, assuming that the integrands are zero at each limit of integration. Depending upon the ion to electron mass ratio, the relativistic or non-relativistic version of uy is used in the transverse current integration. A precautionary correction is carried out in the final calculation of the charge den-sity. This correction ensures that there is no imbalance of charge over the periodic length Chapter 7 The Vlasov Code 74 of the simulation plasma slab. For the case of fixed ions, the ions are seen as having a total fixed charge normalized to one. Then p = l - J fedux . (7-15) Suppose that the spatial integral of the originally calculated charge density p is Jpdx = R (7-16) instead of zero as hoped for. Let g — J fedux. Then Jgdx = l-R (7-17) which should be 1 for there to be no imbalance of charge. Therefore, let g = g/(l — R). This implies that the corrected charge density is In this way, the electron charge is kept to the correct normalization with respect to the ions. As well, the spatial average of the charge density is zero. This must be the case in a periodic system'20', since Ex is periodic [X + L an / dx = Ex(x + L) - Ex(x) = 0 = Aw <p> . (7 - 19) J x If both species are introduced as being electrons (as in the two-stream instability), then the correction has the form Chapter 7 The Vlasov Code 75 The factor of 2 comes from the need to balance each of the electron distributions. For both an electron and an ion distribution, the previous correction was tested and used. Define a normalization function g such that g = l-P • (7-21) Then for < p >= 0, < g >= 1, giving the same correction factor as before, where R = j pdx. For any of these cases, the correction is small but over many time steps it can become significant. This is especially true if the distribution functions are severly deformed, as when an intense laser field is incident upon the plasma. The longitudinal electric field Ex is then calculated from the charge density p. The integral (7-13) is carried out using the cubic spline routine SINT, which assumes periodic boundary conditions in x. If only one species exists, LONGONE is called instead of LONG. LONGONE dis-regards the second distribution function in the calculation of Ex and jy. This reduces the number of computations necessary. Subroutine M X W L MXWL is where the external fields Bz and Ey and the transverse velocity uy are advanced one time step. Bz is advanced using BTW ~ B°zld dEy T ~ ° dx where the derivative of Ey with respect to x is done utilizing the cubic spline routine SLIP. Note that here B* at the full time step is calculated with a qaudratic fit, using the value of Bz from the previous time step: Bprev. W o l d Wnew B*2 = + ^ + . (7-23) Chapter 7 The Vlasov Code 76 Then uy is advanced in time utilizing this value of 5*, unew _ old y— = (E« - —B*) . 7 - 24 r m * c Again the transverse velocity of the ions is treated relativistically only if the ion to electron mass ratio is small. The non-relativistic version of (7-24) neglects the vxB term as small. If the ions are fixed, only the electron transverse velocity uy>e is advanced. Lastly, the relation Enyew - Ef dBz . ( . is used to advance Ey in time. Here the derivative of Bz is carried out by cubic splines with a call to the subroutine SLIP. The transverse current jy was calculated by integration in LONG or LONGONE. Subroutine SAFT SAFT calculates the Fourier transform of the transverse and longitudinal electric fields for diagnostic purposes. This is done using the subroutines FOUR and FOUC. As well, the total electrostatic energy is found from k Since Ex was found by direct integration in the Vlasov code, it avoids the problems of using a Fourier transform to find the longitudinal electric field. This version of calculating the electrostatic energy was then used rather than the method used in EMI. In calculating the above sum individual electrostatic mode energies were also stored. The mode energies are divided by 1 x I O 2 0 so that the computer can deal with their large values. The transverse electric field energy per mode is found as well, using Ey(k)/Sw. Chapter 7 The Vlasov Code 77 When a pre-deterrnined number of time steps have gone by, the subroutine DUMPIT is called from SAFT. This subroutine stores the above diagnostics, freeing memory storage space for new diagnotic data. CHAPTER 8 LANDAU DAMPING SIMULATIONS As was shown in Chapter 1, Landau damping can be derived from a purely math-ematical analysis. In fact, the process was originally found in this manner and was later shown to correspond to physical reality'4'. The physical explanation of Landau damping involves the resonant behaviour of particles with an electrostatic wave'11'. The particles with speeds very close to the phase velocity = u/k of the wave interact strongly with the wave. Particles moving slightly slower than u^, gain energy to speed up to a velocity v^. Those particles moving slightly faster than v$ lose energy to the wave and and slow to a speed v^. In a Maxwellian distribution (as assumed in the simulations) there are more particles moving slower than than faster, resulting in a net gain of energy in the par-ticles from the wave. Therefore the wave undergoes Landau damping. At large enough amplitudes, particles with velocities close to v$ are trapped within the large potential troughs of the electrostatic wave. Non-relativistic and relativistic damping simulations were carried out. All simula-tions had fixed ions with mobile electrons. Amplitudes of the initial perturbation were kept small (< 1%) when damping rates were determined so that linear theory assumptions were not violated. Larger amplitude simulation results are also given, revealing Case-Van Kampen mode coupling and particle trapping. 8.1 Non-relativistic Landau Damping Non-relativistic Landau damping simulations were carried out at plasma tempera-78 Chapter 8 Landau Damping Simulations 79 tures where Vth < 0.04c. Simulations with different values of kvth/wp were run by chang-ing the temperature of the plasma. This then gave a range of damping values to compare with the theory already presented. As an example, the following input parameters were used in a run of the Vlasov code. V L A S O V I N P U T P A R A M E T E R S Variable Description Value NDX Number of x-grid points 32 NDU Number of ux-grid points 51 ISPEC Number of species (1 or 2) 1 UT1 Thermal velocity (species 1) / c 0.04 DEN Density of electrons / critical density 0.000225 XWL Slab thickness in laser wavelengths 0.533333 NPA Number of perturbations in the plasma slab 1 PAME Amplitude of the plasma perturbations (electrons) 0.010 Table 5: Vlasov Landau damping input parameters. All the other input values have default values. With these parameters, the length of the simulation plasma slab was „ „v WVL.XWL 1 x 10-^ 0.533333) . . . , „ , „ _ VI - DEN VI - 0.000225 v ' giving as the x grid spacing and XMX DLX = — — = 1.667 x 10~6 cm (8-2) N DX. DLT = ^ C T X = i-66?*10"?m(0.5) = 2.778 x 10"17 sec (8 - 3) VOL 3 x 1010cm/secv ; v ; Chapter 8 Landau Damping Simulations 80 as the time step size. The temperature of vth = UT1 * c = 0.04c gives rise to a spacing on the ux grid of TyT-\ 10 x 0.04c-10 x-0.04c A A ,n 8 . ' ,A A\ DUE = = 4.8 x 108 cm/sec (8 - 4) 51 — 1 where UMX, UMN = ±10 x vth- Since mode 1 was excited (NPA= 1), the wavevector of the initial perturbation had a value k = — — = 1.178 x lO'cm-1 . (8-5) XMX As well, with the given value of DEN, the plasma frequency was up = VDEN — ^ - = 2.827 x 1013 sec-1 . (8-6) W V Li A value of kvth = 5.0 (8-7) is obtained from these parameters. Since kvth/up is a rather large in this case, very strong damping should occur. Figure 12 shows a semi-log plot of the electrostatic energy's mode 1 versus time. The damping is so strong that the perturbation damps away before one plasma cycle is over. Recall that the energy (in e.s.u.) was divided by 1 x 1020 before the logarithm was taken. This of course does not affect the damping rate calculation. Chapter 8 Landau Damping Simulations 81 Figure 12: Vlasov electrostatic energy mode 1 time evolution. Another interesting effect is also visible in figure 12. The mode 1 energy repeats its previous history after a certain amount of time. This is the "recurrence" effect which occurs in a simulation plamsa with periodic boundary conditions and a velocity grid'24' and E — 0. Only a small perturbation was used in this run, so that E = 0 is approximately correct. Fourier transforming the spatial part of the Vlasov equation after linearization gives ^ + ivkt = 0 (8-8) at which after integration gives rise to f1=f1(t = 0)e~ivht . (8-9) This result also refers to free-streaming, which will be dealt with later in more detail. With the ux grid of the Vlasov code, it can then be seen that the system repeats itself Chapter 8 Landau Damping Simulations 82 with a period of 2n/Avk. The recurrence time of this run should then be close to 2TT UpTreciir — A , LOp — 3.14 Avk (8 - 10) which it is, as seen in figure 12. Figure 13 shows the progression of the total electrostatic energy over time. A large damping rate is evident. Figure 13: Vlasov total electrostatic energy time evolution. To make a comparison of the Vlasov code simulation runs and the particle code EMI, the input parameters of EMI must be chosen to set up an equivalent initial state. The non-default input parameters of EMI are shown in table 6. Chapter 8 Landau Damping Simulations 83 EMI INPUT PARAMETERS Variable Description Value Nl Number of particles of species 1 32768 N2 Number of particles of species 2 0 NG Number of grid points 32 MODE Mode to perturb 1 Xll Amplitude of perturbation in species 1 26.08 WI Plasma frequency of species 1 0.02 WC1 Cyclotron frequency of species 1 0.0 QM1 Charge to mass ratio of species 1 -1.0 KO System fundamental mode 0.5 DT Time step size 0.07854 TEMPI Temperature of species 1 0.0016 ISD Random or quiet start 0 (quiet start) Table 6: EMI Landau damping input parameters. The number of grid points NG= 32 were chosen to be the same as the Vlasov run. This meant that since the plasma simulation slab has a set length of L= 47r, the grid spacing was Ax = ^ = 0.3927 . (8 - 11) Mode one was excited giving k — 2TT/L = 0.5. The time step size was chosen to be Ai = 0.07854 (8 - 12) so that the speed of light c = Ax/At.= 5. All these values are in EMI simulation units. The temperature of the plasma in EMI is normalized with respect to the electron rest mass, meaning that TE M 1 =-f = 0.0016 (8 - 13) Chapter 8 Landau Damping Simulations 84 or a plasma temperature of T = 800 eV. A quiet start was used since Landau damping requires a good representation of the distribution function and a low noise level. No damping was observed for a random start with only a 1% perturbation, even with 216 particles. This agrees with previous findings'9,20]. 215 particles were used in this quiet start, the 1 % perturbation then being NP Xll = — x 0.01 = 26.08 . (8 - 14) 1J The choice of the plasma frequency Wl= 0.02 kept kvth/wp = 5.0. Figure 14 shows the time evolution of mode 1. Again strong damping is evident, but in this case background noise of the particle simulation clouds later behaviour. The recurrence effect has been seen in other particle codes where the initial distribution was loaded as a large number of individual beams'20!. The total electrostatic energy is plotted in figure 15, showing the strong initial damping and then the growth of background noise. This noise could have been overcome somewhat if a larger initial perturbation was used, but linearity would be violated. Energy values in EMI are given in terms of the rest mass energy of all the electrons in EMI simulation units. The results from non-relativistic, small amplitude Landau damping runs are given in table 7. The resulting frequencies obtained are compared with the solution of the plasma dispersion relation (see appendix A) in figures 16 and 17. As can be seen, the Vlasov code is more accurate than EMI for Landau damping, most likely due to the background noise difference. Small deviations do occur at larger damping rates even for the Vlasov code due to relativistic and non-linear effects. The Vlasov code was also more efficient, taking approximately 2.3 hours to run this short simulation on a Sun 3/60 workstation, while it took EMI 16.4 hours. Chapter 8 Landau Damping Simulations t o « Ener > \ A / A "T3 O f o 3 V ' 1 1 1 1 ' ' ' 1 1 ' ' ' ' 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 3 T i m e ( 1 / W p ) 4 ' • ' • 1 5 Figure 14: EMI electrostatic energy mode 1 time evolution. s " 0.0 1.0 2.0 3.0 4 .0 S.O T i m e ( 1 / W p ) Figure 15: EMI total electrostatic energy time evolution. Chapter 8 Landau Damping Simulations 86 NON-RELATIVISTIC LANDAU DAMPING VLASOV 5.000 3.270 -9.2286 1.581 2.701 -1.7704 0.631 1.555 -0.3287 0.312 1.227 -0.0348 EMI 5.000 3.140 -3.3530 1.581 2.654 -1.7620 0.631 1.543 -0.3211 0.312 1.220 -0.0340 Table 7: Non-relativistic Landau damping results. 8-«r o-VLASOV A-EM1 8-Ji 2.00 o 1 XI 1.00 o i l , -7*^ > i i i . | . i i , | i , i i | , i , , | , , . 0 - 0 . 5 0 . 0 0 . 5 1 .0 1 .5 LOG(K*Vth/Wp) • i i i i < 2 . 0 Figure 16: Comparison with theoretical real frequency. Chapter 8 Landau Damping Simulations 87 LOG(K*Vth/Wp) Figure 17: Comparison with theoretical damping rate. 8.2 Relativistic Landau Damping Table 8 shows the values used in the Vlasov code for a relativistic Landau damping run, along with the resultant simulation parameters. Note the high temperature of the plasma (vth > 0.04c). These parameters give kvth/oJP = 0.5 for vth = 0.2c. The results of this run are given in figures 18 and 19. Here the value of kvth/u)p is smaller and correspondingly the damping rate is much lower. Figures 18 and 19 show that the damping extends over many plasma periods. Chapter 8 Landau Damping Simulations 88 VLASOV PARAMETERS Variable Description Value NDX Number of x-grid points 32 NDU Number of ux-grid points 51 ISPEC Number of species (1 or 2) 1 UT1 Thermal velocity (electrons) / c 0.2 DEN Density of electrons / critical density 0.36 XWL Slab thickness in laser wavelengths 0.533333 NPA Number of perturbations in the plasma slab 1 PAME Amplitude of perturbation (electrons) 0.010 XMX Plasma slab length 6.667 x 10-5 cm DLX x grid spacing 2.084 x 10"6 cm DUE ux grid spacing 4.800 x 109 cm/sec DLT Time step size 3.473 x 10~17 sec K Wavenumber of perturbation 9.420 x 104 cm"1 WP Plasma frequency 1.131 x 1015 sec"1 Table 8: Vlasov relativistic Landau damping parameters. o o I Q " O _ <N I 1 0.0 5.0 10.0 15.0 20.0 25.0 T ime(1/Wp) Figure 18: Vlasov electrostatic energy mode 1 time evolution. Chapter 8 Landau Damping Simulations 89 0 . 0 5 . 0 1 0 . 0 1 5 . 0 2 0 . 0 2 5 . 0 T i m e ( 1 / W p ) Figure 19: Vlasov total electrostatic energy time evolution. An equivalent run with EMI utilized the parameters in table 9. These values had kvth/u!p = 0.5 with vth = 0.2c, as in the Vlasov case. The results from this run are given in figures 20 and 21. Background noise again interrupts the time evolution of the energy but not until many plasma periods have passed. Simulations with identical parameters were run with a range of particle numbers with the expected result that as the number of particles increased, the background noise decreased. Chapter 8 Landau Damping Simulations E M I P A R A M E T E R S Variable Description Value Nl Number of particles of species 1 65536 N2 Number of particles of species 2 0 NG Number of grid points 32 MODE Mode to perturb 1 XII Amplitude of perturbation in species 1 52.15 WI Plasma frequency of species 1 2.0 WC1 Cyclotron frequency of species 1 0.0 QM1 Charge to mass ratio of species 1 -1.0 KO System fundamental mode 0.5 DT Time step size 0.07854 TEMPI Temperature of species 1 0.04 ISD Random or quiet start 0 (quiet start) Ax Grid spacing 0.3927 At Time step size 0.07854 c Speed of light 5 Table 9: EMI relativistic Landau damping parameters. Figure 20: EMI electrostatic energy mode 1 time evolution. Chapter 8 Landau Damping Simulations 91 T i m e ( 1 / W p ) Figure 21: EMI total electrostatic energy time evolution. A number of relativistic runs were carried out keeping kvth/wp fixed at 0.5 but varying ut/j. The comparison of the simulation runs and the theoretical values (calculated by Dr. A. J. Barnard) is given in tables 10 and 11. VLASOV RELATIVISTIC LANDAU DAMPING (KVTH/WP = 0.5) Vth (L0r/Up)th (iOr/u!p)mea3 {Ui/0Jp)th (Ui/u>p)meas 0.200c 1.36394 1.3659 -0.08645 -0.0869 0.150c 1.38846 1.3812 -0.11373 -0.1144 0.100c 1.40424 1.3963 -0.13518 -0.1323 0.050c 1.41291 1.4120 -0.14874 -0.1475 0.025c 1.41500 1.4280 -0.15522 -0.1540 Table 10: Vlasov relativistic Landau damping results. Chapter 8 Landau Damping Simulations 92 EMI RELATIVISTIC LANDAU DAMPING (KVTH /WP = 0.5) Vth (ur/up)th (ur/up)meas (Ui/Up)th ( u i / u p ) m e a s 0.200c 1.36394 1.2823 -0.08645 -0.0651 0.150c 1.38846 1.3659 -0.11373 -0.0998 0.100c 1.40424 1.3963 -0.13518 -0.1336 0.050c 1.41291 1.3963 -0.14874 -0.1567 0.025c 1.41500 1.4120 -0.15522 -0.1796 Table 11: EMI relativistic Landau damping results. 8.3 Phase Space Waves When larger perturbations are applied to a plasma, significant non-linear effects arise. Figure 22 shows the initial Maxwellian distribution function of a plasma with a 10 % spatial perturbation. This and the following results are from the Vlasov code, the distribution's x and ux component axes being labelled by grid points. In this run, vth = 0.2c, DEN= 0.36 and kvth,/u>p = 0.5. As well, more ux grid points were used so that the three-dimensional graphics routine would display a smooth function. Figure 23 shows the distribution function at a time t = 8u>~1, where the excited mode has become contorted and waves have begun to form on the distribution. Figure 24 at t = 12a*-1 shows that as time progresses, more and more of these undulations arise. Figure 25 shows the distribution at a much later time t = 20a;-1. These phase space waves result from the coupling of the initial perturbation to the free-streaming mode of the Vlasov equation'24!. The free-streaming solution to the Vlasov equation with E — 0 was derived before, fx oc eikxeikvt (8 - 15) where exp(ikx) represents the initial perturbation. This result plus additional factors can be derived from the Vlasov equation with E ^ 0 (as for large perturbations)'11'. These modes are in fact Case-Van Kampen modes, derived from the spatial Fourier transform of the Vlasov equation'11'. Equation (8-15) describes the propagation of a wave, Chapter 8 Landau Damping Simulations 93 / i oc c''(**-*«*> (8 - 16) with a velocity "frequency" of kt. A velocity period v can then be defined as As an example, take the distribution function at t = 12a;-1. In this simulation k = 9.424 x 104 cm, up = 1.131 x 1015 sec-1 and DUE= 3 x 108 cm/sec. This gives 27T TTIT = 6.4 x 109 cm/sec . (8 - 18) (9.424 x: 104 cm-1) (12wf') The right half of the distribution function in figure 24 goes to zero at roughly ux grid point 30, which implies thai the right half of the distribution covers about 70 ux grid points or 2.1 x 1010 cm/ssc This then implies that "X'°rCm,/SeC*3 6.4 x 10y cm/sec wavelengths should be ^isaMe on each side of the distribution at t = 12a;-1. This is indeed the case in figure 24. The particle code EMI revealed that this behaviour has a more physical meaning than just a mathematical formalism. The particles in EMI are moved self-consistently without referring to the Vlasov equation, giving a good check on this behaviour. At a 10% perturbation background noise prevented these waves from being seen, but at a 50% perturbation the free-streaming coupling was evident. Figures 26 through 29 show the evolution of these waves as viewed from directly above the distribution function (phase space). Figure 30 and 31 show the velocity distribution of the particles at t = 0 and t = 40a;-1, again revealing the free-streaming convolutions. Chapter 8 Landau Damping Simulations 94 Figure 23: Vlasov distribution function at t = 8o>_1. Chapter 8 Landau Damping Simulations Figure 25: Vlasov distribution function at t = AOLO' 1. Chapter 8 Landau Damping Simulations 96 Figure 26: E M I phase space at t = 0. Figure 27: E M I phase space at t — 8a ; - 1 . Chapter 8 Landau Damping Simulations Figure 29: EMI phase space at t = 40a;-1. Chapter 8 Landau Damping Simulations 98 VX (v/c) Figure 30: EMI velocity distribution function at t = 0. . vx (v/c) Figure 31: EMI velocity distribution function at t = 40wp 1. 8.4 Particle Trapping Chapter 8 Landau Damping Simulations 99 Particle trapping occurs in cases where the amplitude of the perturbation is so large that particles are trapped in the potential troughs of the waves. This would then lead to an initial damping period followed by a short period of growth when particles are trapped. Figures 32 through 34 show the mode energies for a Vlasov code run where a 50% per-turbation was introduced into the initial distribution function. The trapping behaviour is evident, the higher modes having stronger damping rates and trapping growth rates. Figures 35 to 37 show an equivalent run from EMI. The damping rates are much larger than predicted by linear theory as expected in the strongly non-linear regime. Figure 32: Vlasov electrostatic energy mode 1 time evolution. Chapter 8 Landau Damping Simulations 100 Figure 33: Vlasov electrostatic energy mode 2 time evolution. Figure 34: Vlasov electrostatic energy mode 3 time evolution. Chapter 8 Landau Damping Simulations 101 Figure 35: EMI electrostatic energy mode 1 time evolution. Figure 36: EMI electrostatic energy mode 2 time evolution. Chapter 8 Landau Damping Simulations 102 Figure 37: EMI electrostatic energy mode 3 time evolution. CHAPTER 9 WARM TWO-STREAM INSTABILITY Simulating the two-stream instability provided a chance to test the use of two dis-tribution functions in the Vlasov Code. As well, both relativistic and non-relativistic sim-ulations were carried out, providing another check on the relativistic correctness of both codes. The growth rates of the instability measured from the simulations were compared with the numerical solution of the plasma dispersion relation and phase space behaviour was observed. Al l the simulations presented in this chapter represent counter-streaming electron beams with Maxwellian distributions. 9.1 Non-Relativistic Case In the non-relativistic version of the simulations, the plasma temperature was set so that vth = 0.014c, or equivalently T = 0.0002 in E M I . This is equivalent to T = lOOeV since the E M I temperature is normalized to the electron rest mass. For the two-stream instability to occur, it is required that kvo/up < \f2 and VQ > 1.307vth as already shown. As an example, table 12 shows parameters used in an E M I run. For these values, Air A x = — = 0.3927 (9 - 1) and c = 4^ = 3-927 (9 - 2) At K ' 103 Chapter 9 Warm Two-Stream Instability 104 EMI INPUT PARAMETERS Variable Description Value Nl Number of particles of species 1 32768 N2 Number of particles of species 2 32768 NG Number of grid points 32 WI Plasma frequency of species 1 0.1 W2 Plasma frequency of species 2 0.1 V01 Drift velocity of species 1 / c 0.05 V02 Drift velocity of species 2 / c -0.05 WC1 Cyclotron frequency of species 1 0.0 QM1 Charge to mass ratio of species 1 -1.0 QM2 Charge to mass ratio of species 1 -1.0 KO System fundamental mode 0.5 DT Time step size 0.1 TEMPI Temperature of species 1 0.0002 TEMP2 Temperature of species 1 0.0002 ISD Random or quiet start 0 (quiet start) Table 12: EMI two stream instability input parameters. in EMI simulation units. The choice of v0 satisfies v0 > 1.307ut/,. As well, since mode 1 is excited in the instability, k — 0.5 and = 0.5 (0.051(3.927) = 0 9 8 <^ _ u>p 0.1 The instability should therefore occur. A number of EMI runs were carried out using the above parameters but varying the particle number and starting techniques. It was found that as the particle number decreased the growth rate of the instability was not greatly affected but the system reached saturation more quickly. The rapid saturation was probably due to the noisier background associated with fewer particles. If the noise level is high, the electrostatic energy has a higher initial energy and therefore does not have to increase as much to reach saturation. Chapter 9 Warm Two-Stream Instability 105 As expected, it was found that a quiet start with a large number of particles gave the best results. Figure 38 shows the time evolution of the electrostatic energy of mode 1 in a EMI simulation using the above parameters. The exponential growth of the instability (7/wp = 0.4105) is evident, the system saturating around t = 22a;-1. After saturation the energy of mode 1 displays a wavy structure correponding to the bounce frequency of the trapped particles. This is caused by an exhange of energy between the particles and the electostatic wave they are trapped in. Roughly, the period of this phenomena was measured to be Tt, « 5a;-1 implying that the bounce frequency was u>6 = 27r/Ti = 0.126 (in EMI simulation units). From the energy level in figure 38 after saturation, the bounce frequency should be = y/(1.0)v/5~x IO"4 (0.5) = 0.106, (9-4) in fair agreement with the previous value. Figure 39 shows the time evolution of mode 2, figure 40 the evolution of mode 3. Non-linear coupling causes these modes to grow at twice and three times the growth rate of mode l ' 2 0 ' . Indeed, the growth rate of mode 2 was "f/u>p — 0.7542 and of mode 3 7/wp = 1.0527. Figure 41 shows the evolution of the total electrostatic energy over time, the exchange of energy between the particles and the electrostatic wave being clearly seen after saturation. When particles are trapped, vortex patterns form in phase space since the particles oscillate back and forth in a potential trough of an electrostatic disturbance. Figure 42 shows the phase space of an EMI simulation at t = 0. During the exponential growth of the instability, the beams bend towards each other in phase space, as in figure 43 (t — 12a;-1). Figure 44 shows the characteristic vortex pattern of trapping at t = 24a;"1, just after saturation. The thermal energy gained by the electrons is evident in figures 45 and 46 showing the velocity distribution function of the particles at t = 0 and t = 24a;"1. Chapter 9 Warm Two-Stream Instability 106 T i m e (1/Wp) Figure 38: EMI mode 1 time evolution. Figure 39: EMI mode 2 time evolution. Chapter 9 Warm Two-Stream Instability 107 • I • " i I • i i i i •• i • i • i i i T' 1 1 1 I i . i | i . . i | i i i . | 0 5 10 16 20 25 30 35 40 45 50 T i m e ( 1 / W p ) Figure 40: EMI mode 3 time evolution. o 0.0 5.0 10.0 15.0 20.0 25.0 30.0 35.0 40.0 45.0 50.0 T i m e ( 1 / W p ) Figure 41: EMI total electrostatic energy time evolution. Chapter 9 Warm Two-Stream Instability 108 X ( g r i d p o i n t s ) Figure 42: EMI phase space at t = 0. Figure 43: EMI phase space at t = \2u>p 1. Chapter 9 Warm Two-Stream Instability 109 Figure 45: E M I particle velocity distribution at t = 0. Chapter 9 Warm Two-Stream Instability 110 Figure 46: EMI particle velocity distribution at t = 24up 1 . An equivalent initial system was set up in the Vlasov code, utilizing the input param-eters shown in table 13. These parameters were chosen so that the relative.drift velocity of the beams were the same as in EMI and so kvo/u>p = 0.98. The value of DEN was found from the fact that k = 2TT _ 2TT\/1 - DEN XMX " ( l x lO"4) (0.53333) = 1.1781 x l O V l - DEN , (9-5) 2nc v WVL VMN = 1.885 x 1015v/DEN (9-6) and using the requirement of kv0 = 0.98u>„ gave Vl - DEN = 10.4536VTJEN (9-7) Chapter 9 Warm Two-Stream Instability 111 VLASOV INPUT PARAMETERS Variable Description Value NDX Number of x-grid points 32 NDU Number of ux-grid points 101 ISPEC Number of species (1 or 2) 2 MIE Mass ratio between species -1.0 UT1 Thermal velocity (species 1) / c 0.014 UT2 Thermal velocity (species 2) / c 0.014 DEN Density of electrons / critical density 0.009 XWL Slab thickness in laser wavelengths 0.533333 V01 Drift velocity of species 1 / c 0.05 V02 Drift velocity of species 1 / c -0.05 Table 13: Vlasov two stream instability input parameters. or DEN= 0.009. Figure 47 shows the electrostatic energy of mode 1 over time in a Vlasov run using these parameters. Once again the exponential growth of the instability is clearly seen (7/u>p = 0.4329) and the system saturates at t a 30a;-1. The later saturation time is due to the low level of noise in the Vlasov simulations. The measured bounce period Tb RS 2.5a;"1 after saturation implies that u>b — 4.494 x 1014sec_1. This is self-consistent with the average energy level of mode 1 after saturation in figure 47, giving [eEk V m (5.2728 x 1017statcoul./gm)(<y87r(8 x 10~9 1 x 1020) ergs)(1.1728 x 105 cm) (9-8) or u>b — 5.266 x 1014sec-1. All Vlasov values are in electrostatic units. Figures 48 and 49 show mode 2 and 3 respectively. The growth rates j/u!p = 0.8659 and 7/wp = 1.2944 were closer to twice and three times the growth rate of mode 1 than produced by EMI. Figure 50 shows the time evolution of the total electrostatic energy, again showing the exchange of energy between particles and the electrostatic wave of the instability. Chapter 9 Warm Two-Stream Instability 112 Figure 47: Vlasov mode 1 time evolution. Figure 48: Vlasov mode 2 time evolution. Chapter 9 Warm Two-Stream Instability 113 T i m e ( 1 /Wp) Figure 49: Vlasov mode 3 time evolution. 0 . 0 5 .0 10.0 15.0 2 0 . 0 2 5 . 0 .30.0 3 5 . 0 4 0 . 0 4 5 . 0 5 0 . 0 T i m e (1/Wp) Figure 50: Vlasov total electrostatic energy time evolution. Chapter 9 Warm Two-Stream Instability 114 The phase space behaviour is also visible in the Vlasov simulations. Figure 51 shows the distribution functions at t = 0. Figure 52 show the distributions at t = 15o;-^ during the growth phase of the instability where the two beams are bending towards each other in phase space. Figure 53 shows the vortex formation due to particle trapping after saturation (t = 32a;-1). A number of simulations were carried out with differing relative beam velocities, the results being given in table 14. The growth rates obtained are compared with the numer-ical solution of the plasma dispersion relation (see appendix A) in figure 54. As seen in this figure, EMI and the Vlasov code results agree with theory except at small growth rates. In this region the Vlasov code agrees with theory but EMI deviates probably due to background noise interference. The Vlasov code was again more efficient, the above run taking approximately 12.6 hours on a Sun 3/60 workstation, while EMI took 54.6 hours. Runs with smaller growth rates took much longer (up to the order of days) due to later saturation times. Figure 51: Vlasov distribution function at t = 0. Chapter 9 Warm Two-Stream Instability 115 Figure 53: Vlasov distribution function at t = 32a;-1. Chapter 9 Warm Two-Stream Instability 1 1 6 Figure 54: Comparison of measured and theoretical growth rates. NON-RELATIVISTIC TWO STREAM INSTABILITY EMI ±v0 Sat-OO 7/u;p Mode 1 j/up Mode 2 ~f/wp Mode 3 Tb 0.030 45 0.1756 . 0.1685 0.9955 21 0.035 37 0.2169 0.3819 0.8907 11 0.040 28 0.3831 0.5762 1.0021 8 0.050 22 0.4105 0.7542 1.0527 •5 0.060 25 0.3738 0.6367 1.0497 7 0.070 32 0.2819 0.4042 1.0431 8 0.075 62 0.1365 0.2264 0.7673 10 VLASOV 0.030 51 0.2576 0.5246 0.7865 24 0.035 44 0.2950 0.6180 0.8772 13 0.040 35 0.3811 0.7587 1.0819 6 0.050 30 0.4329 0.8659 1.2944 2.5 0.060 36 0.4001 0.8073 1.1998 5 0.070 45 0.2930 0.5821 0.8909 8 0.075 71 0.1365 0.1365 0.3950 12 Table 14: Two stream instability results. Chapter 9 Warm Two-Stream Instability 117 9.2 Relativistic Case To provide a check on how correct both EMI and the Vlasov code are in a relativistic region, a two-stream instability was simulated with a plasma temperature Vth = 0.142c. Table 15 shows the input parameters used in the EMI simulations. Using these values, the speed of light had a value of c = 3.927 as before. As well, the choice of VQ = 0.33c satisfies VQ > 1.307ft/, and = (0.5) (0.33) (3.927) = 1 2 % <^ Up 0.5 so the instability should arise. E M I I N P U T P A R A M E T E R S Variable Description Value Nl Number of particles of species 1 32768 N2. Number of particles of species 2 32768 NG Number of grid points 32 WI Plasma frequency of species 1 0.5 W2 Plasma frequency of species 2 0.5 V01 Drift velocity of species 1 / c 0.33 V02 Drift velocity of species 2 / c -0.33 WC1 Cyclotron frequency of species 1 0.0 QM1 Charge to mass ratio of species 1 -1.0 QM2 Charge to mass ratio of species 1 -1.0 K0 System fundamental mode 0.5 DT Time step size 0.1 TEMPI Temperature of species 1 0.02 TEMP2 Temperature of species 1 0.02 ISD Random or quiet start 0 (quiet start) Table 15: EMI relativistic two stream instability input parameters. Figure 55 shows the growth and saturation of the instability in EMI. A wavy struc-ture appears quite clear in this graph after saturation. Figure 56 shows the time evolution Chapter 9 Warm Two-Stream Instability 118 of the total electrostatic energy. The oscillations after saturation in both graphs are due to the energy exchange between the particles and the electrostatic wave of the instability as in the non-relativistic case. The amplitude of these oscillations die out in time as the particles which do not have the exact bounce frequency of the potential trough escape'20'. The higher modes in EMI sometimes showed non-linear coupling as before, but not for very small growth rates. As well, when this coupling existed, the higher modes did not have growth rates that were even multiples of the growth rate of mode 1. Figure 55: EMI mode 1 time evolution. Chapter 9 Warm Two-Stream Instability 119 § 1 A oo o A a — p /1 1 Enei § ~o ros S_ a Elecl 300 0.0002 1 i | i r*i i . | i | i | i | i | i i . | i i i i i | i | i i 3 10 20 30 40 50 60 70 60 90 100 110 120 130 140 150 T i m e ( 1 / W p ) Figure 56: EMI total electrostatic energy time evolution. The equivalent Vlasov run utilized the parameters in table 16. With these parameters, k = 1.063 x 105 cm_1, uip = 8.1296 x 1014 sec_1 giving kv0/ojp = 1.296 as in the EMI simulations. Figure 57 shows the progression of the electrostatic energy's mode 1. The exponential growth can be seen with oscillations after saturation. The amplitude of these oscillations decreases with time as in the EMI simulation. Figure 58 shows the time evolution of the total electrostatic energy with a smooth decrease in the amplitude of the oscillations after saturation. Modes 2 and 3 were also not non-linearly coupled when the growth rate was very small. When this coupling existed, however, the growth rates of the higher modes were closer multiples of the growth rate of mode 1 than EMI produced. Chapter 9 Warm Two-Stream Instability 120 VLASOV INPUT PARAMETERS Variable Description Value NDX Number of x-grid points 32 NDU Number of Uz-grid points 101 ISPEC Number of species (1 or 2) 2 MIE Mass ratio between species -1.0 UT1 Thermal velocity (species 1) / c 0.142 UT2 Thermal velocity (species 2) / c 0.142 DEN Density of electrons / critical density 0.186 XWL Slab thickness in laser wavelengths 0.533333 V01 Drift velocity of species 1 / c 0.33 V02 Drift velocity of species 1 / c -0.33 Table 16: Vlasov relativistic two stream instability input parameters. -5.0 0.0 1 ^ § U J 1 • )de -— o 5- -~ cT 1 ~ o i i Q ' o — O *ri CN ' 0 , i i , | i i i i | i , i i | i , i i | , , , , | , , , , | , i , i | i . i . | . . . i ! i i i i | 0 2 5 . 0 5 0 . 0 7 5 . 0 1 0 0 . 0 1 2 5 . 0 1 5 0 . 0 175.0 2 0 0 . 0 2 2 5 . 0 2 5 0 . 0 T i m e ( 1 / W p ) Figure 57: Vlasov mode 1 time evolution. Chapter 9 Warm Two-Stream Instability 121 T i m e (1/Wp) Figure 58: Vlasov total electrostatic energy time evolution. A number of simulations were carried out with differing relative beam velocities, the results being given in table 17. Figure 59 compares the simulation growth rates with the theoretical growth rates obtained from numerically solving the dispersion relation (rela-tivistically correct — see appendix A). EMI again deviates at low growth rates probably due to background noise. Both EMI and the Vlasov code underestimate the maximum growth rate. Chapter 9 Warm Two-Stream Instability 122 RELATIVISTIC TWO STREAM INSTABILITY EMI ±VQ Sat.(u,-i) 7/wp Mode 1 7/wp Mode 2 j/tjjp Mode 3 Tb 0.230 100 0.0135 0.0000 0.0000 80 0.255 50 0.1083 0.1734 0.2641 70 0.285 90 0.1372 0.2238 0.3168 60 0.330 45 0.1671 0.1811 0.2631 25 0.370 50 0.1302 0.1794 0.2898 75 0.390 55 0.0412 0.0888 0.0235 80 VLASOV 0.230 110 0.0221 0.0000 0.0000 140 0.255 120 0.0899 0.1739 0.2529 60 0.285 90 0.1534 0.2891 0.4262 40 0.330 75 0.1784 0.3537 0.4909 30 0.370 75 0.1451 0.2664 0.4158 35 0.390 73 0.0917 0.2262 0.3552 40 Table 17: Relativistic two stream instability results. © - o-VLASOV o • A-EM1 /Wn '"P, ' o . <L> D or : o -CM —, -C o -"5 o o© -o -8 : O 0. • ' i ' ' ' ' i ' ' ' ' i 1& 0.20 0.25 0.30 0.35 0.40 B e a m V e l o c i t y ( V O / c ) Figure 59: Comparison of measured and theoretical growth rates. CHAPTER 10 STIMULATED RAMAN SCATTERING SIMULATIONS Raman scattering simulations were carried out with both codes. The inclusion of intense external fields deformed the distribution function of the plasma significantly, requiring that the number of grid points in both codes be increased. The Vlasov code required this increase mainly in the x grid to keep periodicity exact and to deal with the contorting distribution function. This increase in the number of x grid points in the Vlasov code unfortunately lengthened its run times. The increase in the number of grid points necessary for EMI was not as great as for the Vlasov code. Generally the Vlasov code took longer to run these simulations than EMI. 10.1 Low Laser Intensity: Vlasov Code The initial parameters were set up in the Vlasov code to have a simple Raman scattering system adhering to periodic boundary conditions. The incident pump wave was chosen to have four wavelengths Ao inside the simulation plasma slab. The backscatter wavelength As was set to be a simple multiple of Ao (As = TIXQ), which implied that ks = ko/n and therefore ws = to2 + c2k2s = to2 + c2k\jn2. The system must keep the Raman matching conditions for backscatter, so that 123 Chapter 10 Stimulated Raman Scattering Simulations 1 2 4 where the approximation toe « u>p was made. Squaring both sides gives ^(i-4)=i+2./i+44 ooi n2, V loin* Letting p = c2&o/w2 and simplifying the above expression, one obtains the quadratic equation P2( l - ir) - 2 p ( l + - ir)-3 = 0 with the solution (1 + 1/n2) + y/(l + 1/n2)2 + 3(1 - 1/n2)2 P (1-1/n2)2 where the positive sign was chosen so that p is positive (as it must be). Also note that col i i DEN = - 4 LO 2 1 + cHl/uj l+p • Choosing n = 4, we found p = 3.4665 and DEN= 0.2239. This value of density then should create a system where four wavelengths of incident light and one wavelength of backscatter light fits into the plasma slab. Since ke = ka + ks for Raman backscatter, the produced plasma perturbation would have five wavelengths fitting into the simulation plasma slab. These values coincide with the natural plasma wavevector giving the maximum growth rate (equation (3-39)), ke = k0 + —4/1 - — - • c V w0 Chapter 10 Stimulated Raman Scattering Simulations 125 Since u>e « u>p, one obtains DEN k e = k o + ^L c V u>0 ^ 1/ (ke-k0)2C2\ "2o A ul J DEN = I f i - ^ i £ y = 0.22 With the wavevector values chosen, &o = 4, k$ = 1 and fce = 5, 4~( 1 _ 42 in agreement with the value of DEN previously obtained. A Vlasov code simulation was run with the parameters in table 18. VLASOV INPUT PARAMETERS Variable Description Value NDX Number of x-grid points 256 NDU Number of u^ -grid points 51 ISPEC Number of species (1 or 2) 1 UT1 Thermal velocity (species 1) / c 0.1 VOC Electron quiver velocity / c 0.04 DEN Density of electrons / critical density 0.2239 XWL Slab thickness in laser wavelengths 4 NPA Number of perturbations in the plasma slab 5 PAME Amplitude of the plasma perturbations (electrons) 0.01 Table 18: Vlasov Raman scattering input parameters. The value of the oscillatory velocity of the electrons in the incident field VQ/C = 0.04 is large enough for Raman scattering to occur but is below the plasma electron thermal 1 Chapter 10 Stimulated Raman Scattering Simulations 126 velocity vth = 0.1c. With these values, XMX= 4.54 x 10-4 cm, implying that DLX= 1.773 x 10-6 cm and DLT = 2.956 x 10-17 sec . As well, u>p = vTJM-^- = v / 0^239 2?3* 1Q*° = 8.919 x 1 0 1 4 sec"1 . v WVL 1 x IO-4 With these values a run to t = 200a;-1 takes about 7500 time steps. A small initial perturbation in mode 5 of the x distribution was introduced to shorten the saturation time of the instability, therby shortening run times. As well, recall that the electrons are given an initial transverse velocity in an attempt to shorten run times. Figure 60 shows the growth of mode 5 of the electrostatic energy, the natural mode to be excited. Initially Landau damping of the excited wave and the growth of the in-stability interact until about t = lOOo;"1 where the instability begins to dominate. The growth rate is roughly 7 = 0.065, in agreement with numerical results obtained by Drake et al'12'. The recurrence time for this simulation was Trecur « 68a;-1. Figure 60 reveals that the recurrence effect does not occur in this system. Figure 61 displays the electro-static field over the x grid at t = 65a;-1, and it is clear that the expected mode is excited (five wavelengths in the plasma slab). Figure 62 is a plot of the electrostatic field around saturation (i = 130a;-1). The increase in amplitude of the electrostatic field is due to the growth of the instability. Figure 63 shows the time evolution of mode 4 of the trans-verse electric field (the incident pump field). Pump depletion occurs around saturation as expected. However, the initial growth of this mode is anomalous and may be caused by aliasing. Figure 64 shows the time evolution of mode 1 of the transverse electric field Chapter 10 Stimulated Raman Scattering Simulations 127 (backscatter). Its increase is evident and growth seems to continue at a lesser rate after saturation. Figure 65 shows the average mode amplitude distribution for the transverse electric field. Note the the significant backscatter at mode 1. T i m e ( 1 / W p ) Figure 60: Vlasov electrostatic energy mode 5 time evolution. 5-2 ° -?; o T • i i i i i i i i i i i i i i i i i i i . . . > i 3 5 0 1 0 0 1 5 0 2 0 0 2 5 0 X - G r i d P o i n t Figure 61: Vlasov electrostatic field Ex at t = 65wp 1. Chapter 10 Stimulated Raman Scattering Simulations 128 Figure 63: Vlasov transverse electric field mode 4 time evolution. Chapter 10 Stimulated Raman Scattering Simulations 129 Figure 65: Vlasov average mode distribution. Chapter 10 Stimulated Raman Scattering Simulations 130 Figure 66 shows the distribution function at t = 65a;-1, again revealing that mode 5 of the plasma perturbation is growing. Figure 67 shows the transverse velocity of the electrons over the x-ux grid. The electrons follow the large electric field of the incident pump wave. Four wavelengths of the incident light wave fit into the simulation plasma slab as shown in this figure. Figure 68 depicts the distribution function at a later time, t — 104a;-1, closer to saturation. The increase in the amplitude of mode 5 is evident in comparison with figure 66. Figure 69 shows the distribution function at t — 130a;-1, around saturation. The distribution function is highly contorted with the natural mode excitation at its fullest. The vortex structures visible reveal particle trapping in the large potential trough of the instability. This feature agrees with the phase space view of trapping in the Raman instability seen by Klein et al'14!. This trapping is part of the saturation mechanism in Raman scattering'14]. Figure 66: Vlasov electron distribution function at t = 65a;"1. Chapter 10 Stimulated Raman Scattering Simulations 131 Chapter 10 Stimulated Raman Scattering Simulations 132 Figure 69: Vlasov electron distribution function at t = 130o;p 1. A similar run including ions was then carried out, with mi/me — 400 and Te/T,- = 40. Figure 70 shows the growth of the instability in mode 5 of the electrostatic energy. The measured growth rate was f/u>p = 0.060. This, was lower than for the case when m,- —• oo, in agreement with other findings'25,261. Saturation again occurs around t = 130a;-1. The addition of a second distribution function nearly doubles the run time required of the code, and runs beyond t « 130a;-1 become difficult (run time > 3 days). Figure 71 shows a semi-log plot of the mode 4 energy of the transverse electric field (incident pump wave), again showing the anomalous initial growth and pump depletion around saturation time. Figure 72 shows a similar graph of the growth of the backscattered light. Figure 73 shows the electron distribution function near saturation at t = 120a;-1. The vortex pattern does not seem to form at this or any time to the same extent as when the ions were kept fixed. Figure 74 shows the distribution function of the ions at t = 120a;-1. A perturbation in the ion distribution function can be seen. This reaction may slow the vortex formation Chapter 10 Stimulated Raman Scattering Simulations 133 of particle trapping. Figure 75 shows the evolution of the total electrostatic energy. Figure 70: Vlasov electrostatic energy mode 5 time evolution. Figure 71: Vlasov transverse electric field mode 4 time evolution. Chapter 10 Stimulated Raman Scattering Simulations 134 8-f _ cn (D c: U J OJ S -"D O 1 o o 1 m — c 1 I ' 1 ' i ' I 1 1 ' i—'—I—'—i—i—i—<—i—i—i—•—i—i—i 10 20 30 40 60 ©0 70 80 00 100 1 10 120 130 T i m e ( 1 / W p ) Figure 72: Vlasov transverse electric field mode 1 time evolution. Figure 73: Vlasov electron distribution function at t = 120u>_1. Chapter 10 Stimulated Raman Scattering Simulations 135 Figure 75: Vlasov total electrostatic energy time evolution. Chapter 10 Stimulated Raman Scattering Simulations 136 10.2 Low Laser Intensity: EMI An equivalent simulation was run with the particle code EMI using the parameters in table 19. E M I INPUT P A R A M E T E R S Variable Description Value Nl Number of particles of species 1 32768 N2 Number of particles of species 2 0 NG Number of grid points 256 MODE Mode to perturb 5 IMODE Mode of initial transverse field 4 EO Amplitude of radiation fields in plasma 0.091 XII - Amplitude of perturbation in species 1 0.01 WI Plasma frequency of species 1 1.074 WC1 Cyclotron frequency of species 1 0.0 QM1 Charge to mass ratio of species 1 -1.0 KO System fundamental mode 0.5 DT Time step size 0.049 TEMPI Temperature of species 1 0.01 ISD Random or quiet start 0 (quiet start) Table 19: EMI Raman scattering input parameters. With these parameters, c = Ax/ At = 1. The plasma frequency and the amplitude of the input field were chosen to make the system correspond to the Vlasov simulation. Since the incident light wave has four wavelengths in the plasma slab of length 47r, AO = n, implying that kQ = 2. As well, the EMI plasma frequency WI must be chosen from the fact that Wl= \/ne/ncrwo = \/0.2239 UQ. Knowing that 2 2 , 2 ; 2 LOQ =LOp+C kQ one obtains UJO = 2.27 Chapter 10 Stimulated Raman Scattering Simulations 137 Also, the transverse electric field amplitude can be found from requiring VQ/C = 0.04 and from the relation EQ e m W0 Since |e/m| = 1.0 and c = 1, the amplitude of the initial transverse electric field is EQ = w0— = (2.27) (^) = 0.091 c 1 in EMI simulation units. Figure 76 shows the time evolution of the mode 5 electrostatic energy. The growth of the instability is evident, saturating at t = 230W-1. The saturation time is longer than for the Vlasov code because of a much smaller initial perturbation. A smaller perturbation was given because the run times for EMI were reasonable and the small growth rate would be hidden in any larger perturbation (especially with the higher noise levels of EMI). The measured growth rate was "f/u>p = 0.026. This value is lower than that of the Vlasov code and the value obtained by Drake et al'12!, but does agree more closely to the simple linear theory growth rate of Chapter 3: 7 kuo / 1 wp 4 y wp(wo - wp) = 0.022 with the assumption that we « wp. Figure 77 shows the time evolution of the transverse electric field mode 4 (pump). Pump depletion is evident at saturation and the anoma-lous growth before saturation does not occur. Figure 78 pictures the time evolution of mode 1 of the transverse electric field (backscatter). The backscatter amplitude grows until saturation where it levels off for a short time and then begins a slower growth after saturation. Figure 79 shows a phase space plot of the electrons at saturation, with the same features as the Vlasov distribution function at saturation. Chapter 10 Stimulated Raman Scattering Simulations 138 o —. in 7 "1 ' ' ' • ' " * ' ' I 1 1 1 ' I ' ' ' 1 I ' ' 1 1 I ' 1 1 1 I ' 0 SO lOO 160 200 250 300 T i m e ( 1 / W p ) Figure 76: EMI electrostatic energy mode 5 time evolution. Time (1/Wp) Figure 77: EMI transverse electric field mode 4. Chapter 10 Stimulated Raman Scattering Simulations 139 >-UJ 3 Time (1/Wp) Figure 78: EMI transverse electric field mode 1. > X • • i • • 1 • i • • 1 1 1 • • • • i • • 1 1 1 • • 1 1 1 1 • 1 1 1 1 • 1 1 1 1 1 1 1 1 1 • • • i • 0 . 0 2 5 0 5 0 . 0 7 5 . 0 1 0 0 . 0 1 2 5 . 0 1 5 0 . 0 1 7 5 . 0 2 0 0 . 0 2 2 5 . 0 2 5 0 . 0 X ( g r i d p o i n t s ) Figure 79: EMI electron phase space at t = 220up 1 Chapter 10 Stimulated Raman Scattering Simulations 140 An EMI simulation was also carried out with ions of finite mass, m,/m e = 400, Te/Ti — 40. Figure 80 shows the growth of mode 5 of the electrostatic energy. The growth rate ~f/up = 0.021 is lower than when m; oo, in agreement with the Vlasov code. As well, the system saturates at a later time. Figure 81 shows an electron phase space plot and figure 82 shows an ion phase space plot, both close to saturation at t = 350a;-1. The ions are slightly perturbed as in the Vlasov case. 430 Time (1/Wp) Figure 80: EMI electrostatic energy mode 5 time evolution. Chapter 10 Stimulated Raman Scattering Simulations 141 > 1 1 1 ' 1 • • i • • • • i •' • 1 1 1 • • 1 1 1 • 1 1 1 1 1 1 1 1 1 1 1 1 1 • 1 1 • i • 1 • 1 1 • 0 2 5 SO 75 1 0 0 125 I S O 1 7 5 2 0 0 2 2 S 2 5 0 X ( g r i d p o i n t s ) Figure 81: EMI electron phase space at t = 350u> - l Q I 1 1 1 1 I 1 ' 1 1 I 1 1 1 ' I 1 1 1 1 I 1 1 1 1 I 1 1 1 ' I 1 1 1 1 I 1 1 1 1 ! 1 1 1 ' I 1 1 1 1 I 1 1 0 . 0 2 5 . 0 5 0 . 0 7 5 . 0 1 0 0 . 0 1 2 5 . 0 1 5 0 . 0 1 7 5 . 0 2 0 0 . 0 2 2 5 . 0 2 5 0 . 0 X ( g r i d p o i n t s ) Figure 82: EMI ion phase space at t = 350^ 1 . Chapter 10 Stimulated Raman Scattering Simulations 142 10.3 High Laser Intensity Similar Raman scattering simulations were repeated with VQ/C = 0.4. This intense field severly deformed the plasma distribution function and necessitated 1024 x grid points in the Vlasov code. Only one run was carried out in this region due to very long run times (6.5 days !). A larger perturbation at mode 5 was introduced into the plasma in both codes in an attempt to quicken saturation times. Figure 83 shows the rapid growth of mode 5 of the electrostatic energy, saturating at t = 15a;-1. The measured growth rate was -j/up — 0.35, significantly higher than the low intensity expected. After saturation the energy level decreased somewhat and then began increasing, indicating an oscillation was starting. This oscillation may be due to particle trapping and energy exchange between the plasma and the transverse fields^27). The erratic behaviour at the far end of the run is numerical in origin, the distribution function being displaced off the grid system due to extreme deformation. Figure 84 shows the backscatter (mode 1, transverse field) which grows until saturation. After saturation the amplitude level decreases slightly. Figure 85 shows the time evolution of the initial pump wave. This time no growth occurs before saturation and depletion occurs after saturation, in agreement with EMI. Figure 86 shows the time evolution of mode 5 of the electrostatic energy in an EMI equivalent simulation. Growth is again seen until saturation with evidence of an oscilla-tion starting afterwards as in the Vlasov simulation. Figure 87 shows the increase in the backscatter, again reducing in intensity somewhat after saturation. Figure 88 shows the pump depletion in mode 4 of the transverse electric field. Evidently the Vlasov code and EMI agree much better at this higher level of initial pump intensity. Chapter 10 Stimulated Raman Scattering Simulations 143 7 I ' • ' — I - " — • — ' — [ — ' — 1 ' 1 I 1 1 — 1 ' I 0 . 0 5.0 10.0 15.0 2 0 . 0 2 5 . 0 3 0 . 0 T i m e ( 1 / W p ) Figure 83: Vlasov electrostatic energy mode 5 time evolution. T i m e ( 1 / W p ) Figure 84: Vlasov transverse electric field mode 1 time evolution. Chapter 10 Stimulated Raman Scattering Simulations 144 T i m e ( 1/Wp) Figure 85: Ylasov transverse electric field mode 4 time evolution. T i m e (1/Wp) Figure 86: E M I electrostatic energy mode 5 time evolution. Chapter 10 Stimulated Raman Scattering Simulations 145 400 500 1 >^ L U ) 100 200 I I > I , | i | i i'"-!^ 1 ' 1 ' 1 1 1 10 2 0 3 0 4 0 5 0 6 0 7 0 BO 9 0 1 0 0 110 T i m e Figure 87: E M I transverse electric field mode 1 time evolution. i—'— i—•—i—«—i—<—i—>—i—>— i—>— i— 1 —i—<—i—'—i— '—i 0 10 2 0 3 0 4 0 S O 6 0 7 0 BO 9 0 1 0 0 110 T i m e Figure 88: E M I transverse electric field mode 4 time evolution. CHAPTER 11 CONCLUSIONS The object of the work presented here was to construct and test the performance of a 1-1/2 dimensional, electromagnetic, relativistic Vlasov code capable of dealing with a two-component plasma. Simulations were carried out for Landau damping, the two-stream instability and stimulated Raman scattering. The comparison of these simulation results, both with theory and with simulations by the particle code EMI, revealed that in electrostatic problems the Vlasov code preforms very well. It was found to be relativisti-cally correct even with two distribution functions present. The Vlasov code also required much less computation time to run an electrostatic simulation due to the low noise level in the code. Landau damping simulations revealed that the Vlasov code is very accurate and fast in simulating electrostatic problems. The damping rates and frequencies obtained agreed very closely with theory, closer than EMI results. The observation of damping was obscured in EMI by noise, especially after a number of plasma periods. Relativistic results were also tested, the Vlasov code again being the most accurate and efficient code. Case-Van Kampen modes were observed with both codes. EMI required a larger perturbation than the Vlasov code to observe these modes, even with a quiet start. This was again due to background noise in EMI. The energy evolution characteristic of particle trapping was observed for a large initial perturbation with both codes. The good performance of the Vlasov code in Landau damping simulations is attributed to a lower background noise than EMI and a much better representation of a Maxwellian distribution function. Two stream instability simulations revealed that even with two species the Vlasov code was faster and more accurate than EMI in electrostatic problems. Both codes produced growth rates in non-relativistic and relativistic regimes agreeing with theory. 146 Chapter 11 Conclusions 147 EMI, however, deviated at low growth rates. Both codes showed non-linear coupling of higher electrostatic modes except for small growth rates in the relativistic case. Vortex formation in phase space characteristic of particle trapping was also observed with both codes after saturation. Saturation times for the Vlasov code were generally longer than for EMI simulations, again due to background noise in the particle code. The Vlasov code was found to have anomalous features not in agreement with EMI when external radiation fields were present. In stimulated Raman scattering simulations the initial pump field grew until saturation, unlike EMI results. Both codes, however, showed pump depletion at saturation and backscatter growth . As well, both showed growth in the expected plasma mode, although the Vlasov code produced a higher growth rate than EMI. Vortex formation was seen in both simulations, revealing electron trapping as a saturation mechanism for stimulated Raman scattering. When very intense fields were present the Vlasov code became very inefficient due to the necessity of having a large number of grid points to describe the distribution function. Therefore the Vlasov code is not yet a viable choice in dealing with very intense radiation. The results produced from these runs, however, were in much better agreement with EMI. The incident pump wave did not have an initial growth sequence and all the other necessary features, such as pump depletion, were evident. It is believed that the anomalous features produced by the Vlasov code can be overcome. The better performance of the Vlasov code for higher laser intensities indicates that the code can produce realistic results with external radiation fields present in the simulation plasma. Aliasing may cause the undesirable features. The high laser intensity simulation that produced good results used a large number of x grid points. This increase in the number of grid points would reduce aliasing significantly by sampling functions well below the Nyquist frequency. This indicates that increasing the number of grid points in future simulations would improve results. As well, the application of a digital filter to certain parameters may also improve results by reducing aliasing effects. The recurrence Chapter 11 Conclusions 148 effect does not seem to be a cause of the anomalous behaviour of the Vlasov code, since no noticeable effect was seen in Raman scattering simulations after a number of recurrence periods. In the course of debugging the Vlasov code it was noticed that the integration to find the transverse current was especially sensitive to minor alterations. A change in the transverse current altered the anomalous growth of the incident pump wave in Raman scattering simulations. This sensitivity may be due to a multiplication of aliasing effects, since two functions are multiplied and then integrated to find the transverse current. Corrections in this area of the code should be concentrated upon in further endeavors. The Vlasov code provides a method of simulation that can be very accurate. Vlasov simulations also have a very low noise level and a good representation of the velocity distribution function. With these features, interesting effects could be seen in simulations that a particle code would have difficulty producing. Such an effect would be the plasma echo, where information is stored in the convolutions of the distribution function. This requires a very good description of the distribution function. Once the anomalies in the Vlasov code are corrected, more reliable results can be obtained from simulations of Raman and Brillouin scattering. Interesting studies of the dependence of the these scattering processes on laser intensity and plasma temperature could be carried out and compared with theoretical results derived by Drake et al'12'. Other phenomena that require external fields could also be modelled, such as laser beat heating of a plasma and the observation of a precursor. The excellent observation of the distribution function provided by the Vlasov code is unique and can lead to physical insight. Particle trapping phenomena and plasma perturbations can therefore be seen readily. Development of the Vlasov code to include other dimensions should not be as difficult as for a particle code, finite differencing being necessary only in time since cubic splines are used to find spatial derivatives. Including other dimensions would allow a greater range of phenomena to be modelled, such as the two-plasmon decay. All these features encourage further development of the Vlasov code. Chapter 11 Conclusions 149 The majority of this work was carried out on a Sun 3/60 workstation. This enabled program development and the debugging process to be carried out easily and conveniently. The main-frame system at U.B.C. was plagued by slow reaction time due to the large number of users, frequent system shutdowns and by the cost for computing time. This is characteristic of any main-frame system. The workstation used in this work had a quick response time since only one user was operating it at any moment. Periodic system shutdowns were unnecessary for the workstation, a data backup being done individually. Most attractively, there was no cost for computing time (only the initial cost of the system and a service contract) which saves a large amount of money during the development of any program. These features and ease of use make a workstation attractive even when compared to a supercomputer such as the Cray. The computational speed of a main-frame computer system is at the moment greater than that of a workstation. Simulation runs of any great length, however, would be put on low priority in a main-frame system, increasing the time between requesting a simulation run and receiving the resulting data. This delay makes even a heavy task's run times comparable for a main-frame system and a workstation. Very demanding runs, such as a high intensity radiation simulation, take too long to run feasibly on a workstation. More powerful workstations are presently being marketed which should decrease the computational speed difference between workstations and main-frames. References 150 REFERENCES [1] Jackson, J. D., Longitudinal plasma oscillations. Nucl. Energy, Part C: Plasma Physics 1:171-189, 1960. • [2] Abramowitz, Milton, and Irene A. Stegun, eds. Handbook of Mathematical Functions With Formulas, Graphs, and Mathematical Tables, National Bureau of Standards, Department of Commerce, U.S.A., November 1970. [3] Stringer, T. E., Electrostatic instabilities in current-carrying and counterstreaming plasmas. Nucl. Energy, Part C: Plasma Physics, 6:267-279, 1964. [4] Chen, F. F. Introduction to Plasma Physics and Controlled Fusion: Volume 1, Plenum Press, 1985. [5] Bohm D., and E. P. Gross, Theory of plasma oscillations. Phys. Rev. 75:1851-1864, 1949. [6] Jackson, E. Atlee, Drift instabilities in a maxwellian plasma. Phys. Fluids 3(5):786-792, 1960. [7] Roberts, K. V., and H. L. Berk, Nonlinear evolution of a two-stream instability. Phys. Rev. Letters, 19(6):297-300, 1967. [8] McKee, Christopher F., One-dimensional simulation of relativistic streaming instabilities. Phys. Fluids, 14(10):2164-2176, 1971. [9] Parachoniak, Ron. Numerical experiments using a 1-D electrostatic relativistic plasma simulation model. Master's thesis, University of British Columbia, September 1986. [10] Lampe, Martin, and P. Sprangle, Saturation of the relativistic two-stream instability by electron trapping. Phys. Fluids, 18(4):475-481, 1975. [11] Nicholson, D. R. Introduction to Plasma Theory, John Wiley and Sons, New York, 1983. [12] Drake, J. F., P. K. Kaw, Y. C. Lee, G. Schmidt, C. S. Liu and MarshaU N. Rosenbluth, Parametric instabilities of electromagnetic waves in plasmas. Phys. Fluids, 17(4):778-785, 1974. References 151 [13] Hiob, Eric, and A. J. Barnard, Vlasov simulation of stimulated Raman scattering in one dimension. Phys. Fluids, 26(10):3119-3126, 1983. [14] Klein, H. H., Edward Ott, and Wallace M. Manheimer, Effect of trapping on the saturation of the Raman backscatter instability. Phys. Fluids, 18(8):1031-1033, 1975. [15] Forslund, D. W., J. M. Kindel, and E. L. Lindman, Theory of stimulated scattering processes in laser-irradiated plasmas. Phys. Fluids, 18(8):1002-1016, 1975. [16] Kruer, William L. The Physics of Laser Plasma Interactions. Frontiers in Physics, Addison-Wesley Publishing Company, 1988. [17] Aithal, Shridhar, Pierre Lavigne, Henri Pepin, and Tudor Wyatt Johnston, Kent Estabrook, Superthermal electron production from hot underdense plasmas. Phys. Fluids, 30(12):3825-3831, 1987. [18] Dr. T. J. Boyd, University of Wales, Bangor, personal correspondence. [19] Hora, Heinrich, Nonlinear confining and deconfining forces associated with the interaction of laser radiation with plasma. Phys. Fluids, 12(1):182-191, 1969. [20] Birdsall, Charles. K., and A. Bruce Langdon. Plasma Physics via Computer Simulation, McGraw-Hill, New York, 1985. [21] Buneman, 0., Time reversible difference procedures. J. Comput. Phys., 1:1211-1219, 1967 [22] Jackson, J. D. Classical Electrodynamics, John Wiley & Sons, New York, 1975. [23] Yanenko, N. N. The Method of Fractional Steps, Springer-Verlag, New York, 1970. [24] Cheng, C. Z. and Georg Knorr, The integration of the Vlasov equation in configuration space. J. Comput. Phys, 22:330-337, 1976. [25] Eastabrook, Kent, and W. L. Kruer, Theory and simulation of one-dimensional Raman backward and forward scattering. Phys. Fluids, 26(7):1892-1903, 1983. [26] Kruer, W. L., Kent Estabrook, B. F. Lasinski, and A. B. Langdon, Raman backscatter in high temperature, inhomogenous plasmas. Phys. Fluids, 23(7):1326-1329, 1980. References 152 [27] Forslund, D. W., J. M. Kindel, and E. L. Lindman, Plasma simulation studies of stimulated scattering processes in laser-irradiated plasmas. Phys. Fluids, 26(10): 3119-3126, 1983. [28] Burden, Richard L., and J. Douglas Faires. Numerical Analysis. Third edition. Prindle, Weber and Schmidt, 1985 [29] Estabrook, Kent, W. L. Kruer, and B. F. Lasinski, Heating by Raman backscatter and forward scatter. Phys. Rev. Letters, 45(17):1399-1403, 1980 [30] Liu, C. S., Marshall N. Rosenbluth, and Roscoe B. White, Raman and Brillouin scattering of electromagnetic waves in inhomogeneous plasmas. Phys. Fluids, 17(6):1211-1219, 1974 APPENDIX A PROGRAM LISTINGS FOR THEORETICAL RESULTS A.l Non-Relativistic Landau Damping The dispersion relation obtained in Chapter 1, is solved by estimating the damping rate and finding a self-consistent real frequency from the imaginary part of the above expression. Once this has been done the real part of the dispersion relation can be solved for kvt^/ujp. This procedure utilizes the subroutine from the U.B.C. Computing Center which calculates the complex error function. complex*16 z,wz real*8 x1yJtry,lhs(10000)1rhs(10000),xst(10000) real*8 k,pi)spi21de,co,si1re,im)s2,step,wr,wi1kvp,trs real*8 a,b,stork(10000))storx(10000),story(10000) integer iij open(l,file='data') write(6,*)"Real frequency (1)" write(6,*)" Damping rate (2)" read(5,*)iij ( A - l ) C-C c-c c c c c c This programs solves the non-relativistic dispersion relation for both the real frequency and the damping rate. Non-relativistic dispersion relation - Landau damping 153 Appendix A: Program Listings for Theoretical Results 154 write(6,*)iij pi=3.141592654 trs=0.434294481903d0 spi2=dsqrt(pi/2.d0) s2=dsqrt(2.d0) C C ***Loop through estimates for the damping rate. C do 1000 y=0.00000001d0,4.d0,0.01d0 step=0.1d0 a=0.1d0 b=100.d0 C C ***Loop through real frequency values to find a C ***self-consistent value. C 50 do 10 x=a,b,step de=2.d0*dexp((y*y-x*x)/2.d0) co=dcos(x*y) si=dsin(x*y) j=j+l lhs(j)=de*(x*co+y*si) z=dcmplx(x/s2,y/s2) call w(z,wz) re=dreal(wz) im=dimag(wz) rhs(j)=x*re-y*im xst(j)=x C C ***Lhs and rhs are the left and right sides of the C ***imaginary part of the dispersion function C try=abs(lhs(j)-rhs(j)) if(lhs(l) .gt. rhs(l)) then if(lhs(j).lt.rhs(j)) then step=step/2.d0 a=xst(j-l) b=xst(j) 1=0 go to 50 endif endif if(lhs(l) .It. rhs(l)) then if(lhs(j).gt.rhs(j)) then step=step/2.d0 a=xst(j-l) b=xst(j) j=0 go to 50 endif Appendix A: Program Listings for Theoretical Results 155 end if if(try .It. 0.000001) then C C ***If the imaginary part of the dispersion function is C ***satisfied then calculate K*Vth/Wp from the real C ***part of the dispersion function. C i=i+l k=-l.d0+spi2*de*(x*si-y*co)+spi2*(x*im-|-y*re) . k=dsqrt(k) stork(i)=log(k)*trs storx(i)=log(x*k)*trs story(i)=log(y*k)*trs go to 1000 endif 10 continue step=step/2.d0 go to 50 1000 continue C C C ***Write out the LOG of the results. C if(iij.gt.l) go to 1500 do 1200 jj=i,l,-l write(l,100)stork(jj),storx(jj) 1200 continue if(iij.eq.l) go to 1600 1500 continue do 1300 jj=i, 1,-1 write(l,100)stork(jj),story(jj) 1300 continue 1600 continue C C ***Calculate the linear approximation to C ***the real frequency and the damping rate. C kvp=dsqrt(pi/8.d0) do 2000 y=0.1d0,180.d0,0.01d0 wr=dsqrt(l+3.d0*y*y) if(iij.eq.l) then write(l,100)log(y)*trs,log(wr)*trs go to 2000 endif if (y.lt.0.5d0) go to 2000 x=l.d0/y wi=kvp*x**3*dexp(-wr*wr*x*x/2.d0) write(l,100)log(y)*trs,log(wi)*trs 2000 continue 100 format(3fl5.10) end Appendix A: Program Listings for Theoretical Results 156 The Complex Error Function Subroutine subroutine w(z,erfz) complex*16 z,erfz implicit real*8 (a-h,o-z) x=dreal(z) y=dimag(z) call wz(x,y,wx,wy) erfz=dcmplx(wx,wy) return end subroutine wz(x,y,wx,wy) implicit real*8 (a-h,o-z) xx=dabs(x) yy=dabs(y) call wwl(xx,yy,wx,wy) if (x. ge. O.dO .and. y .ge. O.dO) return if (x .It. O.dO .and. y. ge. O.dO) go to 10 dcst=2.d0*dexp(yy*yy-xx*xx) wx=dcst*dcos(2.d0*xx*yy)-wx wy=dcst*dsin(2.d0*xx*yy)+wy if (x .gt. O.dO) return wy=-wy return 10 wy=-wy return end subroutine wwl(x,y,wx,wy) real*8 h,h2,lambda,rl,r2,s)sl,s2,tl,t2,c,x,y,wx,wy integer capn,nu,n,npl,ndo logical b if (y .It. 4.29 .and. x .It. 5.33) go to 10 h=0.d0 capn=0 nu=8 h2=0.d0 lambda=0.d0 s=0.d0 go to 20 10 s=(l.d0-y/(4.29d0))*dsqrt(l.d0-x*x/(28.41d0)) h=1.6d0*s h2=h+h capn=6+23*s lambda=h2**capn nu=9+21*s 20 b=. false. if (h .eq. O.dO .or. lambda .eq. O.dO) b=.true. Appendix A: Program Listings for Theoretical Results 157 rl=0.dO r2=0.d0 sl=0.dO s2=0.d0 ndo=nu+l do 40 i=l,ndo n=ndo-i npl=n+l tl=y+h+npl*rl • t2=x-npl*r2 c=0.5d0/(tl*tl+t2*t2) rl=c*tl r2=c*t2 if(h .gt. O.dO. and. n .le. capn) go to 30 go to 40 30 tl=lambda+sl sl=rl*tl-r2*s2 s2=r2*tl+rl*s2 lambda=lambda/h2 40 continue wy=1.12837916709551d0 if'(y .eq. O.dO) go to 50 wx=1.12837916709551d0 if (b) wx=wx*rl if(.not. b) wx=wx*sl go to 60 50 wx=dexp(-x*x) 60 if (b) wy=wy*r2 if(.not. b) wy=wy*s2 return end A.2 Fluid Theory Growth Rates for the Two Stream Instability The growth rates as derived from fluid theory are found for cold and warm coun-terstreaming beams. C : C Plots the theoretical growth rates C for the two-stream instability. C Cold versus warm (not including the C temperature spread). C : C real*8 v0,wp,k,storv(10000),storw(10000) real*8 t,vl1v2,vi)fzero,ftemp,c C C ***Defining the theoretical growth functions Appendix A: Program Listings for Theoretical Results 158 fzero(v0)=dsqrt(wp**4+4*k**2*v0**2*wp**2)-(k**2*v0**2-|-wp**2) ftemp(v0)=dsqrt(wp**4+4*k**2*v0**2*wp**2+4*k**4*v0**2*t) 1 -(k**2*(v0**2+t)+wp**2) C open(l,file='data') wp=0.1dO k=0.5d0 c=3.927 write(6,*)" Minimum Speed (Simulation Units) ?" read(5,*)vl write(6,*)" Maximum Speed (Simulation Units) ?" read(5,*)v2 write(6,*)" Speed iteration (Simulation Units) ?" read(5,*)vi \vrite(6,*)"do you want: cold beams (1)" write(6,*)" : warm beams (2)" write(6,*)" : both (3)" read(5,*)iij ii'(iij.eq.l) go to 10 write(6,*)"Temperature of the warm beams ?" read(5,*)t 10 v2=v2+vi if(iij.eq.2) go to 13 do 50 v0=vl,v2,vi i=i+l storw(i)=sqrt(fzero(v0)) storv(i)=v0/c 50 continue do 60 j=i,l,-l write( 1,100)stor v(j) ,storw(j) 60 continue if(iij.eq.l) stop 13 continue ii=0 do 70 v0=vl,v2,vi ii=ii+l storw(ii)=sqrt(ftemp(v0)) storv(ii)=v0/c 70 continue do 80 jj=l,ii write(l,100)storv(jj),storw(jj) 80 continue 100 format(2fl5.10) end Appendix A: Program Listings for Theoretical Results 159 A.3 Necessary Two-Stream Instability Velocity Spread The critical parameter values for the two-stream instability are found from the dis-persion relation . ^ = 2 ( , | M i y ( l ) } - l ) (A - 2) as shown in Chapter 3. The algorithm obtained from the U.B.C. Computing Center was again used to find the necessary values of the complex error function (listing on previous pages). The program iterates through values of velocity spread / thermal velocity = 2VQ I vth, the value of k2v2h/u>2 being found from the above dispersion relation. C — C Find the k values which allow the C two-stream instability (non-relativistic). C C C C The allowed range of modes for the two-stream C instability are found* by solving the real part C of the non-relativistic dispersion relation. C This utilizes the aamplex error function routine. C C C complex*16 z.erfz real*8 a,b,c,k2,eta;#ta0,etal,etai,spi,pi C open(l,file='data'jl C write(6,*)"Initia8 and final values of 2*Vo/vth ?" read(5,*)eta0,etal write(6,*)" Iteration, size ?" read(5,*)etai C pi=3.141592654d0) spi=dsqrt(pi) a=dsqrt(2.d0) write(6,*)eta0,etal,etai C C Appendix A: Program Listings for Theoretical Results 160 C ***determine the k values*** C C do 10 eta=etaO,etal,etai b=eta/a z=dcmplx(b,0.d0) C the complex error function call w(z,erfz) c=dimag(erfz) C the dispersion relation k2=2*(eta*(spi/a)*c - 1) write(l,20)(2*eta),dsqrt(k2) 20 format(2fl5.10) 10 continue end A . 4 D e t e r m i n a t i o n o f the T w o - S t r e a m I n s t a b i l i t y G r o w t h R a t e s The theoretical growth rates of the warm two-stream instability were found by solving the plasma dispersion relation 4=7 m* (A - 3) J v - iy — oo where a pure growth rate u> — iy has been assumed. This integral can be done along the real u-axis since the pole in the integrand is in the upper complex plane. The actual integral carried out was the real part of the above dispersion relation, J vl+yl UJ 2 p — oo Here G(v) = dfo/dp with the normalization J fodp =1 (A - 5) / Appendix A: Program Listings for Theoretical Results 161 and p is the momentum per unit mass. A relativistic Maxwellian distribution f ° = 2. .g,(«) { A ~ 6 ) was used for each beam, where q = mc2/kgT, K\ is the modified Bessel function and eq is placed in both the numerator and denominator of fo to keep evaluated values within the range of a computer. The relativistic relation between p and v from Chapter 5 can be used to show that since dv v2 3^2 sr'1-?' {A-7) the integral in question can be written as 4=T^* (A - 8) ui J v2 + yl and e?-9/v/l-t,2/<:2 '•=-si*<ir • ( A- 9 ) Now, since we are dealing with the two-stream instability, a drift velocity must be added relativistically to the distribution for each beam. Assume that each beam has an equal and opposite velocity. Using the fact that4 1 + VVQ/C1 where VQ is the drift velocity, the distribution function for the system of counter-streaming beams becomes Appendix A: Program Listings for Theoretical Results 162 Using this distribution function dfo/dv then becomes 2e<*Kxc2 (A-12) The plasma dispersion relation equation (A-4) was then solved iteratively, finding the growth rate for possible beam velocities. Parameter values used in solving this were determined from EMI particle code runs, and then converted to dimensionless quantities so that comparisons with the Vlasov code could be made. The program utilized, an algorithm from the U.B.C. Computing Center for determining the values of the modified Bessel function, converted to run on the Sun 3/60 by Dr. A. J. Barnard. The upper and lower limits, theoretically iinfinity, had to be chosen so that the computer could deal with the size of the numbers in question. This restriction on the size of the limits of the integration would not affect the resulting answer significantly due to the small size of the integrand at large values. The smallest produced numbers would come from the exponential in dfo/dv. Letting the smallest reasonable number be approximately 1 x 10-50, then the restriction becomes 1 ) = /n(l x 10"50) = -115 (A - 13) •y/l — V2/6 This gives Appendix A: Program Listings for Theoretical Results 163 as the maximum value of the upper limit on the integration. The lower limit would then be the negative of this. C C Two-Stream Instability growth rates C using the plasma dispersion relation. C C C C ***This program solves the real part of the dispersion C ***function using numerical integration (adaptive quadrature C ***procedure for Simpson's rule [ "Numerical Analysis, C ***Burden, R. L. and Faires, J. D. ] ) l 2 8 l The algorithm C ***loops through all the required relative speeds of the C ***beams. For each velocity it finds an estimate for the C ***growth rate of the instability which gives k=0.5, C ***as chosen in the particle simulations. From these C ***estimates the approximate growth rate is calculated C ***for the two-stream instability. All units are as in C ***particle code EMI. C c implicit real*8 (a-h,o-z) real*8 a,b,tol,app,toli(10000),ai(10000),hi( 1 0 0°0) real*8 fai(10000),fci(10000),fbi(10000),si(10000),y,pi,spi real*8 Ii(10000),fd,fe,sl)s2,vl,v21v3)v4)v5,v6,v7,v8,t,x real*8 try,kay(10000),v0,wp,rk)c)bek,tmp,fp)fm real*8 bekl real*8 g integer n,i jj C C C ***The distribution function of the beams C g(x) (Relativistic Maxwellian). The C functions fp and fm are for simplifying C writing the function. Note that the C drift velocities are added C relativistically to the distribution. C Due to the units used (not e.s.u.), an extra speed C of light was put into the denominator. *** C C C C C fp(x)=l.d0-(x+v0)**2/(c*c*(l.d0+x*v0/(c*c))**2) fm(x) = l.dO-(x-vO)**2/(c*c*(l.dO-x*vO/(c*c))**2) g(x)=-tmp/(2.d0*bek*c*c*c)*x/(x*x+y*y) . • * ( dexp(tmp-tmp/dsqrt(fp(x)))/(dsqrt(fp(x))*fp(x)) *(x+v0)*(l.d0-v0*v0/(c*c))/(l+v*v0/(c*c))**3 + dexp(tmp-tmp/dsqrt(fm(x)))/(dsqrt(fm(x))*fm(x)) Appendix A: Program Listings for Theoretical Results 164 *(x-v0)*(l.d0-v0*v0/(c*c))/(l-v*v0/(c*c))**3 ) C open(l,file='../data/data') C C C ***Set initial parameters in plasma*** C C C C speed of light c=3.927d0 t=0.02d0 wp=0.5d0 rk=0.5d0 temperature(T/mc*c) plasma frequency Mode C C C ***Get Bessel function value for normalization*** C ( Function subroutine from U.B.C. Computing C center and modified to work on the sun 3/60 C by Dr. A. J . Barnard ). C tmp=l.d0/t bek=bekl(tmp,ierr) if(ierr.eq.O) then ' write(6,*)" Error in Modified Bessel Function" stop endif tmp=l.d0/t C C C ***Set integration parameters*** C C C upper limit C b=c*dsqrt(l.d0-l.d0/(l.d0+115.d0*t)**2) write(6,*)b C lower limit a=-b C tolerance tol=0.0000001d0 C maximum number of steps n=1000 C C C pi=3.141592654d0 spi=dsqrt(2.d0*pi) Appendix A: Program Listings for Theoretical Results 165 C C ***loop through required speeds C write(6,*)"bekr\bek do 40 v0=0.7d0,3.2d0,0.1d0 C write(6,*)v0 C C ***loop through estimates for Wi/k C do 50 y=0.01d0,2.0d0,0.001d0' C C ***Integration of the real part of the dispersion function C app=0.d0 i=l toli(i) = 10.d0*tol ai(i)=a hi(i)=(b-a)/2.d0 x=a fai(i)=g(x) x=a+hi(i) fci(i)=g(x) x=b fbi(i)=g(x) si(i)=hi(i)*(fai(i)+4*fci(i)+fbi(i))/3.d0 li(i)=l.d0 10 if(i.lt.l) go to 1000 x=ai(i)+hi(i)/2.d0 fd=g(x) x=ai(i)+3.d0*hi(i)/2.d0 fe=g(x) sl=hi(i)*(fai(i)+4*fd+fci(i))/6.d0 s2=hi(i)*(fci(i)+4*fe+fbi(i))/6.d0 vl=ai(i) v2=fai(i) v3=fci(i) v4=fbi(i) v5=hi(i) v6=toli(i) v7=si(i) v8=li(i) i=i-l if((abs(sl+s2-v7)) .It. v6) then app=app+sl+s2 else if(v8 .ge. n) then write(6/)" Procedure Failed after" ,n,"steps" stop else i=i+l ai(i)=vl-|-v5 fai(i)=v3 Appendix A: Program Listings for Theoretical Results 166 fci(i)=fe fbi(i)=v4 hi(i)=v5/2.d0 toli(i)=v6/2.d0 si(i)=s2 . Ii(i)=v8+1 i=i+l ai(i)=vl fai(i)=v2 fci(i)=fd fbi(i)=v3 hi(i)=hi(i-l) toli(i)=toli(i-l) si(i)=sl li(i)=li(i-l) endif endif go to 10 1000 continue C C ***Find which estimate of the growth rate is most correct C jj=jj+l kay(jj)=dsqrt(wp*wp*app) if(jj.lt.2) try=abs(kay(l)-rk) if(abs(kay(jj)-rk).lt.try) then try=abs(kay(jj)-rk) endif if(try.lt.0.005) then write(l,100)v0/c,y*kay(jj)/wp C write(6,*)v0/c,y*kay(jj)/wp go to 2000 endif 100 format(2fl5.10) 50 continue 2000 jj=0 try=0.d0 40 continue end The Modified Bessel Function Subroutine c real*8 function bskl(x,int) c c 3 May 88 c. routine from cntr via Tom Nicol c. modified for the SUN by ajb c. bskl = kl(x) = DBESK1 in amdahl library c. bekl = kl*exp(x) = DBSEK1 Appendix A: Program Listings for Theoretical Results 167 implicit real*8 (a-h,o-z) dimension p(5),q(4),pp(10)lqq(8)lf(5),g(2) data xinf,argmin,xmax/ 1 7.237005577e+75, 1.381786969e-76, 1.778556746e+02/ data p/ 1 1.203176761e-01, 2.499784339e+01, 1.797134565e+03, 2 4.433331009e+04, 1.798473002e+05/ data q/ 1 2.500000000e-01, -7.035978939e+01, 9.316074668e+03, 2 -5.537343720e+05/ data f/ 1 2.107618740e-03, 4.076347371e-01, 3.092138339e+01, 2 9.517918567e+02, 8.501049475e+03/ data g/ -2.216786554e+02, 1.700209895e+04/ data pp/ 1 6.982646013e-02, 6.943429490e+00, 1.016245135e+02, 2 5.190226277e+02, 1.207119909e+03, 1.430398225e+03, 3 9.005826148e+02, 2.986452401e+02, 4.833007694e+01, 4 2.958171782e+00/ data qq/ 1 3.017103765e+01, 2.281150284e+02, 6.730124200e+02, 2 9.246057876e+02, 6.385720888e+02, 2.244322525e+02, 3 3.767671737e+01, 2.360279593e+00/ int=l if(x.lt.argmin) go to 200 a l.OdO 1 if (x .gt. xmax) go to 220 if (x .le. l.OOdO) go to 4 a = dexp(-x) go to 20 c entry for bekl entry bekl(x,int) int - 2 if(x.lt.argmin) go to 200 a = l.dO 3 if (x .gt. l.OdO) go to 20 a = dexp(x) c 0 .It. arg .le. 1.00 4 if (x .It. 3.0d-39) go to 5 XX = X * X sump = ((((p(l)*xx+p(2))*xx+p(3))*xx+p(4))*xx+p(5))*xx+q(4) sumq = ((q(l) * xx + q(2)) * xx + q(3)) * xx + q(4) sumf = (((f(l)*xx + f(2))*xx + f(3))*xx + f(4))*xx + f(5) sumg = (xx + g(l)) * xx + g(2) bskl = ((xx*dlog(x)*sumf/sumg + sump/sumq)/x)*a return c return for small arg 5 bekl = l.OdO / x return c 1.00 .It. arg 20 xx = l.OdO / x sump = pp(l) Appendix A: Program Listings for Theoretical Results 168 do 30 i = 2, 10 sump = sump * xx + pp(i) 30 continue c sumq = xx c do 40 i = 1, 7 sumq = (sumq + qq(i)) * xx 40 continue c sumq = sumq + qq(8) bskl = a * sump / sumq /dsqrt(x) return c error return for arg .It. argmin 200 int=0 resl = xinf go to 260 c error return for arg .gt. 177.855-220 int=0 resl = O.OdO 260 bskl = resl return end APPENDIX B LISTING OF THE VLASOV CODE Comments on the Vlasov Code Parameters ctx Courant factor for timestep - usually 0.5 dlt = dlx/vol*ctx den ne/nc = ratio of density to critical density dlt t interval in sec dix x interval in cm mie ion to electron mass ratio erne electron charge/mass ratio esu units emi ion charge/mass ratio esu units emx laser field amplitude in esu ndx number of cells - x-intervals ndu number of u grid points nfx number of Fourier components of Ex saved nfy Ey nts current time step npa Number of perturbations in the slab omg laser angular frequency pame Perturbation amplitude in electron distribution pami Perturbation amplitude in ion distribution rkd k/kd where k = 2*pi*npa/length kd = sqrt(den)*2pi*c/uth/lambda ubc max velocity in integration / uth umn min velocity in integration umx max velocity in integration use umn/uth utl uth/c = sqrt(kB*T/mc**2) (electrons) ut2 uth/c (ions) ute thermal velocity of electrons in cm/s uti thermal velocity of ions in cm/s voc vosc/c ...defines laser irradiance vol c = velocity of light wvl laser wavelength in cm xmx slab thickness in cm xwl xmx/wvl = slab thickness in laser wavelengths Arrays fe(x,u) Electron distribution function fi(x.u) Ion distribution function uye(x,u) Transverse velocities of electrons uyi(x.u) Transverse velocities of ions 169 Appendix B: Listing of the Vlasov Code 170 bz(x) Magnetic field az(x) Magnetic field stored one time step in reverse ex(x) Longitudinal electic field ey(x) Transverse electric field ge(x) Transverse current sp(x) Spline interpolation sq(x) ...coefficients -sr(x) ... calculated once ss(x) Used for temporary st(x) ... storage in the su(x) ... spline routines sv(x) tx(x) FT of EX ty(x) FT of EY uxe(x) Velocity grid for fe uxi(x) Velocity grid for fi xv(x) Temporary storage, yv(x) ... i.e. for FFT ek(mode,time) storage of electrostatic energy per mode tye(mode,time) storage of transverse electric field energy per mode es(time) storage of total electrostatic energy Include file V I N C C v i n e real*8 variables IMPLICIT REAL*8 (A-H,0-Z) COMMON/PARAR/ VOL.MIE.EME.EMI.PIE^OC.WVL.EMX.OMG.XMX.DEN, 1 UT1,UT2,UTE)UTI,DLX,DLT,DUE,DUI,CDT,CTX, 2 PAME,PAMI,RKD,XWL,V01,V02 REAL*8 VOL.MIE.EME.EMI.PIE.VOC.WVL^MX.OMG.XMX.DEN, 1 UT1,UT2,UTE,UTI,DLX,DLT,DUE)DUI,CDT,CTX, 2 PAME,PAMI,RKD,XWL,V01,V02 Include file C I N C C C c i n e integer variables COMMON/PARAI/ NPA,NDX,NDU,NX2)NTS,NFX,NFY) 1 rTF,IDUMP,NT,ISGN,ISPEC INTEGER NPA,NDX,NDU,NX2) 2 NTS,NFX,NFY,ITF,IDUMP,NT,ISGN,ISPEC Include file M I N C C C m i n e matricies C . COMMON/VASP/FE(512,201),FI(512,201),UYE(512,201), 2 UYI(512,201),EK(20,500), 3 TYE(20,500),GE(512),EX(512), 4 EY(512),BZ(512),AZ(512), 5 UXE(201),UXI(201),XV(512)1YV(512), Appendix B: Listing of the Vlasov Code 171 6 SH(512),SP(512),SQ(512))SR(512), 7 SS(512),ST(512),SU(512),SV(512)> 8 ES(IOOO), 9 TX(16384),TY(16384) REAL*8 FE(512,201),FI(512,201),UYE(512,201), 2 UYI(512,201),EK(20,500), 3 TYE(20,500),GE(512),EX(512), 4 EY(512),BZ(512)IAZ(512), 5 UXE(201),UXI(201),XV(512),YV(512), 6 SH(512),SP(512),SQ(512),SR(512), 7 SS(512),ST(512),SU(512),SV(512)1 8 ES(IOOO), 9 TX(16384),TY(16384) B.l Initialization Procedure C v i n i C C C Initialize data for Vlasov shift program C Set up distribution f(x,v) at half step C #include "vine" #include "cine" C C C, E / M (GAUSSIAN) OPEN(4,file=,../datav/DUMP',form='unformatted') OPEN(10,file='VLAS.IN') VOL = 2.99792456D10 EOM = 5.2727638D17 C PI,1/SQRT(2*PI) PIE = 3.14159265358979312 C C C *****DEFAULT VALUES ***** C C c UTH/C=SQRT(KB*T/ME) UTC -.- 0.20d0" C NE/NC DEN = 0.360d0 c THICKNESS/PLASMA WL XWL = 0.533333333d0 C CAN READ IN K/KD INSTEAD OF XWL RKD = 0.0 c CONST. IN MAXWELL INTEGR. CTX = 0.5 c MIN, MAX. V E L / U T H u s e = -5.00 Appendix B: Listing of the Vlasov Code 172 C C C C C C c c c c c c UBC = 5.00 PAM = O.OlOOOdO NPA = 1 VOC = 1.0D-16 WVL = 1.0D-04 NDX = 16 NDU = 51 NFY - 2 NFX = 2 IRL = 1 IDUMP = 100 CALLINPT PERTURBATION AMPLITUDE NO. OF PERT. IN PLASMA VOSC/C (LASER INTENS.) VAC. WAVELENGTH GRID SIZE NDX * NDU NO. COMPONENTS SAVED DUMP DATA EVERY IDUMP STEPS ***INPUT INTIAL PARAMETER VALUES *** AUXILIARY CONSTANTS *** IF(RKD.EQ. 0.) GO TO 30 THICKNESS/WAVELENGTH XWL = DFLOAT(NPA)*UTC/RKD + *DSQRT(1.DO/DEN - 1.D0) SLAB THICKNESS 30 XMX = WVL*XWL/DSQRT(1.DO-DEN) C C C C C C DLX = XMX/DFLOAT(NDX) DLT = DLX/VOL*CTX DX, DT, DU/UTH CDT = VOL*DLT UTH = UTC*VOL THERMAL VELOCITY VELOCITY RANGE IN F(X.V) UMN = USC*UTH UMX = UBC*UTH P/M INTERVAL DLU = (UMX-UMN)/DFLOAT(NDU-l) ***SET UP UX GRID M = NDU/2 N = M + 1 U = DLU UX(N)= O.DO DO 40 I = 1, M UX(N+I) = U UX(N-I) = -U Appendix B: Listing of the Vlasov Code 173 40 U = U + DLU C LASER FIELD AMPLITUDE EMX = 10709.8d0*VOC/WVL C LASER ANGULAR FREQUENCY OMG = 2.0D0*PIE*VOL/WVL C C ***SPLINE INTERPOLATION COEFFICIENTS N = MAX0(NDU,NDX) CALL SPQR(SP,SQ,SR,N) C ***SET UP F AT TIMESTEP 0 CALL INIF C ***X SHIFT : HALF STEP DO 90 J = 1, NDU C H = -UX(J)*DLT/DLX/2.D0 DO 60 I = 1, NDX C SH(I)= H 60 XV(I)= FE(I,J) C N = NDX CALL SHIP(N) DO 80 I = 1, NDX 80 FE(I,J) = YV(I) 90 CONTINUE C C *** DELTA Ul SHIFT C CALL VLSVU1 C C ***NOW F AT HALF TIMESTEP C C remove for now, put after uy is determined so to C find transverse current correctly. C CALL LONG C C *** BZ FIELD: DLT/2 C P =-OMG*DLT/2.D0 D = DLX*DSQRT(1.D0 - DEN)*2.D0*PIE/WVL F = DSQRT(1.D0 - DEN)*EMX DO 140 I = 1, NDX BZ(I)=F*DCOS(P) 140 P = P + D C C *** AZ (BZ stored one time C step in reverse) C P = OMG*DLT/2.D0 DO 145 I = 1, NDX AZ(I)=F*DCOS(P) 145 P = P + D Appendix B: Listing of the Vlasov Code 174 C C 150 C C C C C *** EY FIELD: DLT P = -OMG*DLT DO 150 I=1,NDX EY(I)=EMX*DCOS(P) P =P + D SET UP PARTICLES WITH A SMALL TRANSVERSE VELOCITY **** VOSC (UY): DLT/2 C 170 160 C C C C F = VOC*VOL P = - OMG*DLT/2.D0 DO 160 I = 1, NDX PP=DSIN(P) PP=0.D0 DO 170 J = 1,NDU UY(I,J)= F*PP CONTINUE P = P + D CALL LONG EX FIELD: DLT/2 ***STORE DATA CALL STOR STOP END SUBROUTINE INPT C c C Enter selected input values C C #include "vine" #include "mine" REAL*8 M,K #include "cine" INTEGER NX,IA,IX,IY,IT C C READ(10,*) NX,NU,IS,M READ(10,*) UE,UI,D READ(10,*) X,V READ(10 *.) IA.PE.PI READ(10,*) V1.V2 READ(10 *) IX.IY READ(10,*) IT READ(10,*) ID Appendix B: Listing of the Vlasov Code 175 IF(NX.EQ.O) GO TO 10 NDX=NX 10 IF(NU.EQ.O) GO TO 15 NDU=NU 15 IF(M.EQ.O) GO TO 17 MIE=M 17 IF(IS.EQ.O) GO TO 20 ISPEC=IS 20 IF(UE.EQ.0.) GO TO 30 UT1=UE 30 IF(UI.EQ.0.) GO TO 40 UT2=UI 40 IF(D.EQ.0.) GO TO 50 DEN=D 50 IF(X.EQ.0.) GO TO 60 XWL=X 60 IF(V.EQ.0.) GO TO 70 VOC=V 70 IF(IA.EQ.O) GO TO 80 NPA=IA 80 IF(PE.EQ.0.) GO TO 82 PAME=PE 82 IF(PI.EQ.0.) GO TO 84 PAMI=PI 84 IF(V1.EQ.0.) GO TO 86 V01=V1 86 IF(V2.EQ.O.) GO TO 90 V02=V2 90 IF(IX.EQ.O) GO TO 100 NFX=IX 100 IF(IY.EQ.O) GO TO 120 NFY=IY 120 IF(IT.EQ.O) GO TO 130 ITF=IT 130 IF(ID.EQ.O) GO TO 140 IDUMP=ID 140 CONTINUE RETURN END C S U B R O U T I N E SPQR(SP,SQ,SR,NS) C C 23 Jun 83 C Spline interpolation initialization C Calculate coefficients P, Q, R C REAL*8 SP(1),SQ(1), SR(1) Integer NS Appendix B: Listing of the Vlasov Code 176 SP(1) = 4.0D0 SQ(1) = -0.25D0 SR(1) = -0.25D0 C DO 10 I = 2, NS SP(I) = SQ(I - 1) + 4.0D0 SQ(I) = -1.0D0 / SP(I) 10 SR(I) = -SR(I-1)/SP(I) RETURN END C S U B R O U T I N E I N I F E C c C -Initialize electron distribution function FE(X,V) C -Start with relativistic Maxwellian C -Density (1+PAME*SIN(KX)) C #include "vine" #include "mine" REAL*8 BEK,TMP,BEK1 #include "cine" C SCALED U INTERVAL W=DUE/VOL UTC = UTE/VOL USQ = UTC*UTC TMP = 1.D0/USQ write(6,*)"TMP",TMP BEK=BEK1(TMP,IERR) IF(IERR.EQ.O) THEN WRITE(6,*)"ERROR IN MODIFIED BESSEL FUNCTION" STOP ENDIF BEK=0.5D0/BEK UC = UXE(l)/VOL C NORMALIZATION DONE WITH C BESSEL FUNCTION S = O.DO F = O.DO C UNPERTURBED FE(V) DO 20 J = 1, NDU C C ADD DRIFT VELOCITY RELATIVISTICALLY C V = (UC - V01)/(1.D0+UC*V01) X = (DSQRT(1.D0+V**2)-1.D0)/USQ XV(J)= BEK*W*DEXP(-X) Appendix B: Listing of the Vlasov Code 177 20 UC = UC + W C NORMALIZE FE(V) DO30J = l ,NDU XV(J)= XV(J)*UTE/DUE C C INTEGRATE USING SIMPSON'S RULE AS A CHECK C S = S + (1.D0 + F)*XV(J) 30 F =1 .D0-F S = S/1.5D0*DUE/UTE write(6,*)'TNTEGRAL = ",S C C PERTURBED FE(X,V) C NPA PEAKS OF AMPLITUDE PAME DO 50 I = 1, NDX T = DFLOAT(2*(I-l))*PIE/DFLOAT(NDX) T = T*DFLOAT(NPA) F = 1.0D0 + PAME*DSIN(T) C DO 40 J = 1, NDU 40 FE(I,J)= F*XV(J) 50 CONTINUE RETURN END C C S U B R O U T I N E I N I F I C C C -Initialize ion distribution function FI(X,V) C -Start with relativistic Maxwellian C -Density (1+PAMI*SIN(KX)) C #include "vine" #include "mine" REAL*8 BEK.TMP.BEKl #include "cine" C SCALED U INTERVAL W=DUI/VOL UTC = UTI/VOL USQ = UTC**2 TMP = 1.D0/USQ BEK=BEK1(TMP,IERR) IF(IERR.EQ.O) THEN WRITE(6,*)"ERROR IN MODIFIED BESSEL FUNCTION1 STOP ENDIF BEK=0.5D0/BEK Appendix B: Listing of the Vlasov Code 178 UC = UXI(1)/V0L -C NORMALIZATION USING C BESSEL.EUNCTION S = O.DO F = O.DO C UNPERTURBED FI(V) DO 20 J = 1, NDU C C ADD DRIFT VELOCITY RELATIVISTICALLY C V = (UC - V02)/(1.D0+UC*V02) X = (DSQRT(1.D0+V**2)-1.D0)/USQ XV(J)= BEK*W*DEXP(-X) 20 UC = UC + W C NORMALIZE. EI(V) DO 30 J = 1, NDU XV(J)= XV(J)*UTI/DUI C INTEGRATE USING SIMPSON'S RULE S = S + (1.D0 + F)*XV(J) 30 F =1 .D0-F S = S/1.5D0*DUI/UTI write(6,*)"INTEGRAL = ",S C C PERTURBED FI(X,V) C NPA PEAKS OF AMPLITUDE PAMI DO 50 I = 1, NDX T = DFLOAT(2*(I-l))*PIE/DFLOAT(NDX) T = T*DFLOAT(NPA) F = 1.0D0 + PAMI*DSIN(T) C DO 40 J = 1, NDU 40 FI(I,J)= F*XV(J) 50 CONTINUE RETURN END C : S U B R O U T I N E B E K 1 c : real*8 function bskl(x,int) c c 3 May 88 c. routine from cntr via torn nicol c. modified for the SUN by ajb c. bskl - kl(x) = DBESK1 in amdahl library c. bekl = kl*exp(x) = DBSEK1 c implicit real*8 (a-h,o-z) dimension p(5),q(4),pp(10)Jqq(8),f(5),g(2) data xinf,argmin,xmax/ Appendix B: Listing of the Vlasov Code 179 1 7.237005577e+75, 1.381786969e-76, 1.778556746e+02/ data p/ 1 1.203176761e-01, 2.499784339e+01, 1.797134565e+03, 2 4.433331009e+04, 1.798473002e+05/ data q/ 1 2.500000000e-01, -7.035978939e+01, 9.316074668e+03, 2 -5.537343720e+05/ data f/ 1 2.107618740e-03, 4.076347371e-01, 3.092138339e+01, 2 9.517918567e+02, 8.501049475e+03/ data g/ -2.216786554e+02, 1.700209895e+04/ data pp/ 1 6.982646013e-02, 6.943429490e+00, 1.016245135e+02, 2 5.190226277e+02, 1.207119909e+03, 1.430398225e+03, 3 9.005826148e+02, 2.986452401e+02, 4.833007694e+01, 4 2.958171782e+00/ data qq/ 1 3.017103765e+01, 2.281150284e+02, 6.730124200e+02, 2 9.246057876e+02, 6.385720888e+02, 2.244322525e+02, 3 3.767671737e+01, 2.360279593e+00/ int=l if(x.lt.argmin) go to 200 a = l.OdO 1 if (x .gt. xmax) go to 220 if (x .le. l.OOdO) go to 4 a = dexp(-x) go to 20 c entry for bekl entry bekl(x,int) int = 2 if(x.lt.argmin) go to 200 a = l.dO 3 if (x .gt. l.OdO) go to 20 a = dexp(x) c 0 .It. arg .le. 1.00 4 if (x .It. 3.0d-39) go to 5 XX = X * X sump = ((((p(l)*xx+p(2))*xx+p(3))*xx+p(4))*xx+p(5))*xx+q(4) sumq = ((q(l) * xx + q(2)) * xx + q(3)) * xx + q(4) sumf = (((f(l)*xx + f(2))*xx + f(3))*xx + f(4))*xx + f(5) sumg = (xx + g(l)) * xx + g(2) bskl = ((xx*dlog(x)*sumf/sumg + sump/sumq)/x)*a return c return for small arg 5-bekl = l.OdO / x return c 1.00 .It. arg 20 xx = l.OdO / x sump = pp(l) c do 30 i = 2, 10 sump = sump * xx + pp(i) Appendix B: Listing of the Vlasov Code 180 30 continue c sumq = xx c do 40 i = 1, 7 sumq = (sumq + qq(i)) * xx 40 continue c sumq = sumq + qq(8) bskl = a * sump / sumq /dsqrt(x) return c error return for arg .It. argmin 200 int=0 resl = xinf go to 260 c error return for arg .gt. 177.855 220 int=0 resl = O.OdO 260 bskl = resl return end C SUBROUTINE VLSVU1 c c C Shift U by DELTA Ul C ,C #include "vine" #include "mine" #include "cine" INTEGER N C C E=EME*0.5D0*DLT/DUE C C C DO 190 I=1,NDX C DO 160 J=1,NDU C C Ul=(UXE(J)/VOL)**2 U2=(UYE(I,J)/VOL)**2 R=DSQRT(1.D0+U1+U2) SH(J)=E*(UYE(I,J)/R)*(BZ(I)/VOL) 160 XV(J)=FE(I,J) N = NDU Appendix B: Listing of the Vlasov Code 181 CALL SHIF(N) DO 170 J=1,NDU 170 FE(I,J)=YV(J) C 190 CONTINUE C C IF TWO SPECIES SHIFT SECOND C DISTRIBUTION AS WELL C IF(ISPEC.GT.l) THEN C E=(EMI/2.D0)*DLT/DUI • C DO 290 I=1,NDX DO 260 J = 1,NDU C Ul=(UXI(J)/VOL)**2 U2=(UYI(I,J)/VOL)**2 R=DSQRT(1.D0+U1+U2) SH(J)=E*(UYl(I,J)/R)*(BZ(I)/VOL) 260 XV(J)=FI(I,J) N=NDU CALL SHIF(N) DO 270 J = 1,NDU 270 FI(I,J)=YV(J) 290 CONTINUE C ENDIF C C RETURN END S U B R O U T I N E S T O R C C C Recall data from from Unit 4 C #include "vine" #include "mine" #include "cine" REWIND (4) WRITE(4) VOL,MIE,EME,EMI,PIE, 1 VOC,WVL,EMX,OMG, 2 XMX,DEN,UTE,UTI, 3 DLX,DUE,DUI,DLT, 4 CDT,CTX,PAME,PAMI, Appendix B: Listing of the Vlasov Code 182 5 NTS,NDX,NDU,ISPEC, 6 NPA,NFX,NFY, 7 ITF,IDUMP,ISGN C 20 WRITE(4) (EX(I), I = 1, NDX) WRITE(4) (EY(I), I = 1, NDX) WRITE(4) (BZ(I), I = 1, NDX) WRITE(4) (AZ(I), I = 1, NDX) WRITE(4) (GE(I), I = 1, NDX) WRITE(4) ((UYE(I.J), I = 1, NDX),J = 1,NDU) IF(ISPEC.GT.l) WRITE(4) ((UYI(I,J), I = 1, NDX),J = 1,NDU) C WRITE(4) ((FE(I,J), I=1,NDX), J = 1,NDU) IF(ISPEC.GT.l) WRITE(4) ((FI(I,J), I=1,NDX), J = 1,NDU) C M = MAX0(NDX,NDU) WRITE(4) (SP(I), 1=1,M) VVRITE(4) (SQ(I), 1=1,M) WRITE(4) (SR(I), I=1,M) VVRITE(4) (UXE(J), J = 1, NDU) IF(ISPEC.GT.l) WRITE(4) (UXI(J), J = 1, NDU) C RETURN END B.2 Simulation Loop C* v r u n C c C— Electromagnetic Vlasov code C— Periodic boundary conditions C C #include "vine" #include "mine" ^include "cine" C OPEN DATA FILES C open(4,nle='DUMP',form='unformatted') open(ll,file=,../datav/FE',form='unformatted') open(12,nle='../datav/ESE',form='unformatted') open(13,nle='../datav/LRMODES',form='unformatted') open( 14,nle='../datav/EK',form='unformatted') open(21,file='../datav/UY',form='unformatted') open(22,nle=,../datav/EY',form='unformatted') open(23,file='../datav/BZ',form=,unformatted') open(24,file='../datav/GE',form='unformatted') Appendix B: Listing of the Vlasov Code 183 open(25)file='../datav/EX',form= 'unformatted') C C READ DATA FOR RESTART C CALL RECL C C INTIALIZE DATA FILES C VVRITE(6,6)NTS WRITE(6,7) READ(5*)NTM IF(NTM.EQ.O) GO TO 60 IF(IDUMP.GT.NTM) THEN WRITE(6,*)" IDUMP must be less than or equal to NTM" STOP ENDIF IF(IDUMP.GT.NTM) THEN WRITE(6,*)" IDUMP must be. less than NTM" STOP ENDIF IF(IDUMP.GT.500) THEN WRITE(6,*)"IDUMP must be less than or equal to 500" STOP ENDIF NX2=NDX/2+l WRITE(ll) NDX.NDU.NTS.NTM.ITF.ISPEC \VRITE(12)IDUMP,NTS,NTM,DLT,DEN,OMG WRITE(13)IDUMP,NTS,NTM,NFX,NFY,DLT,DEN,OMG WRITE(14)IDUMP,NTS,NTM,NFX,NFY,DLT,DEN,OMG WRITE(21) NDX,NDU,NTS,NTM,ITF,ISPEC WRITE(22)NTS,NTM,ITF,NDX VVRITE(23)NTS,NTM,ITF,NDX WRITE(24)NTS,NTM,ITF,NDX VVRITE(25)NTS,NTM,ITF,NDX C WRITE(6,200) WRITE(6,210) VVRITE(6,220)NDX,NDU,NTS,NPA,PAME,NFX,NFY)ITF WRITE(6,200) WRITE(6,230) VVRITE(6,240) XMX.VOC WRITE(6,200) " C C C *****TIMESTEP LOOP***** C C DO 20 1=1, NTM NTS=NTS+1 NT=NT+1 C C SAVE THE DISTRIBUTION FUNCTION(S) Appendix B: Listing of the Vlasov Code 184 C C C C C AND FIELDS EVERY ITF STEPS IF((NTS/ITF)*ITF.EQ.NTS) THEN WRITE(11)((FE(II,JJ)1II=1,NDX),JJ=1INDU) WRITE(21)((UYE(II,JJ))II=1,NDX),JJ=1)NDU) IF(ISPEC.GT.l) THEN WRITE(11)((FI(III,JJJ),III=1,NDX),JJJ=1)NDU) WRITE(21)((UYI(II)JJ)III=1,NDX)1JJ=1,NDU) ENDIF WRITE(22)(EY(II),II=1,NDX) WRITE(23)(BZ(II),II=1,NDX) WRITE(24)(GE(II) ,11= 1 ,NDX) WRITE(25)(EX(II),II=1,NDX) ENDIF SHIFT FE BY DELTA U AND DELTA X C C C C C C C C C C C CALL VLSV IF(ISPEC.GT.l) THEN IF(MIE.LT. 100) THEN CALL VLSVI ELSE CALL VLSVION ENDIF ENDIF EX AT NTS+1/2 IF(ISPEC.GT.l) THEN CALL LONG ELSE CALL LONGONE ENDIF IF MORE THAN ONE SPECIES SHIFT FI AS WELL CALL MXWL UPDATE UY.BZ TO NTS+1/2 EY TO NTS + 1 CALCULATE AND SAVE FOURIER TRANSFORMS CALL SAFT 20 CONTINUE CALL DUMPIT 60 CALL STOR STOP STORE DATA Appendix B: Listing of the Vlasov Code 185 6 FORMAT(' Now at step ',15) 7 FORMAT('Enter # of steps - 0 stops:') 200 FORMAT(4X,' ') 210 FORMAT^X.'NDX'^X.'NDU'^X.'NTS'.SX/NPA'^X/PAME', 1 5X,'NFX',4X,'NFY',4X,'ITF') 220 FORMAT(3X,2I6,2I8,3X,F6.4,I6,2I8) 230 FORMAT^X/XMX'^X. 'VOC') 240 FORMAT(3X,F8.6,3X,F6.3) END C C SUBROUTINE RECL C c C Recall data from from Unit 4 #include "vine" #include "mine" #include "cine" INTEGER M,N C REWIND (4) READ(4) VOL, MIE, EME, EMI, PIE, 1 VOC, WVL, EMX, OMG, 2 XMX, DEN, UTE, UTI, 3 DLX, DUE, DUI, DLT, 4 CDT, CTX, PAME, PAMI, 5 NTS, NDX, NDU, ISPEC, 6 NPA, NFX, NFY, 7 ITF, IDUMP, ISGN C 20 READ(4) (EX(I), I = 1, NDX) READ(4) (EY(I), I = 1, NDX) READ(4) (BZ(I), I = 1, NDX) READ(4) (AZ(I), 1=1, NDX) READ(4) (GE(I), I = 1, NDX) READ(4) ((UYE(I,J), I = 1, NDX),J=1,NDU) IF(ISPEC.GT.1).READ(4) ((UYI(I,J), I = 1, NDX),J=1,NDU) C READ(4) ((FE(I,J), I=1,NDX), J=1,NDU) IF(ISPEC.GT.l) READ(4) ((FI(I,J), I=1,NDX), J=1,NDU) C M = MAXO(NDX.NDU) READ(4) (SP(I), I=1,M) READ(4) (SQ(I), I=1,M) READ(4) (SR(I), I=1,M) READ(4) (UXE(J), J = 1, NDU) IF(ISPEC.GT.l) READ(4) (UXI(J), J = 1, NDU) Appendix B: Listing of the Vlasov Code 186 RETURN END C SUBROUTINE VLSV C c C Shift FE* by DELTA U = DELTA Ul + DELTA U2 C and DELTA x C #include "vine" #include "mine" REAL*8 E,E2 #include "cine" INTEGER N.NDl C C *** SHIFT FE *** C ND1=NDU-1 E=EME*DLT/DUE E2=(EME*DLT)**2*0.5D0/DUE C DO 190 I=1,NDX C C DERIVATIVE OF UY W.R.T. UX. C DO 175 J=2,ND1 175 XV(J)=(UYE(I,J-1)-UYE(I,J+1))*0.5D0/DUE XV(1)=XV(2) XV(NDU)=XV(ND1) C C DO 160 J = 1,NDU C C Ul=(UXE(J)/VOL)**2 U2=(UYE(I,J)/VOL)**2 RR=1.D0+U1+U2 R=DSQRT(RR) F=R*RR C C CALCULATE THE SECOND ORDER CORRECTION C TO THE SHIFT C BEF=E2*BZ(I)/(VOL*F) AUX=((l+Ul)*XV(J)-(UXE(I)/VOL)*(UYE(I,J)/VOL)) C C C C SH(J)=E*(EX(I)+(UYE(I,J)/R)*BZ(I)/(2.D0*VOL)) Appendix B: Listing of the Vlasov Code +BEF*AUX*(EX(I)+(UYE(I,J)/R)*BZ(I)/VOL) CARRY OUT THE SHIFT SH(J)=(EX(I)+(UYE(I,J)/R)*BZ(I)/VOL)*(E+BEF*AUX) SH(J)=EX(I)*E XV(J)=FE(I,J) N=NDU CALL SHIF(N) DO 170 J = 1,NDU FE(I,J)=YV(J) CONTINUE SHIFT X BY DELTA X E=DLT/DLX E2=DLT*DLT*0.5D0/DLX DO 290 J = 1,NDU Ul=(UXE(J)/VOL)**2 DO 270 I=1,NDX XV(I)=UYE(I,J) N=NDX DERIVATIVE OF UY W.R.T. X CALL SLIP(N) DO 240 I=1,NDX U2=(UYE(I,J)/VOL)**2 RR=1.D0+U1+U2 R=DSQRT(RR) F=R*RR S=(UXE(J)/VOL)*(UYE(I,J)/VOL)/F DERIVATIVE OVER X-GRID, SO DIVIDE BY DLX (REAL DISTANCE). DVX=-S*SS(I)/DLX DVT=S*EME*(EY(I)-(UXE(J)/R)*BZ(I)/VOL) THE SHIFT IN X Appendix B: Listing of the Vlasov Code 188 SH(I)=-E*UXE(J)/R 1 +E2*((UXE(J)/R)*DVX-DVT) 240 XV(I)=FE(I,J) C N=NDX CALL SHIP(N) C DO 280 I=1,NDX 280 FE(I,J)=YV(I) C 290 CONTINUE RETURN END C S U B R O U T I N E V L S V I c . c C Shift FI* by DELTA U = DELTA Ul + DELTA U2 C and DELTA x (relativistic) C #include "vine" #include "mine" #include "cine" INTEGER N.NDl C C ND1=NDU-1 E=EMI*DLT/DUI E2=(EMI*DLT)**2*0.5D0/DUI C DO 190 I=1,NDX C C DERIVATIVE OF UY W.R.T. UX. C DO 175 J=2,ND1 175 XV(J)=(UYI(I,J-1)-UYI(I,J+1))*0.5D0/DUI XV(1)=XV(2) XV(NDU)=XV(ND1) C C DO 160 J = 1,NDU C C Ul=(UXI(J)/VOL)**2 U2=(UYI(I,J)/VOL)**2 R=DSQRT(1.D0+U1+U2) F=R**P C C CALCULATE THE SECOND ORDER CORRECTION Appendix B: Listing of the Vlasov Code 189 C TO THE SHIFT C BEF=E2*BZ(I)/(VOL*F) AUX=((l+Ul)*XV(J)-(UXI(I)/VOL)*(UYI(I,J)/VOL)) C C C C SH(J)=E*(EX(I)+(UYI(I,J)/R)*BZ(I)/(2.D0*VOL)) C 1 +BEF*AUX*(EX(I)+(UYI(I,J)/R)*BZ(I)/VOL) C C CARRY OUT THE SHIFT C SH(J)=(EX(I)+(UYI(I,J)/R)*BZ(I)/VOL)*(E+BEF*AUX) C 160 XV(J)=FI(I,J) N=NDU CALL SHIF(N) C DO170J = l,NDU 170 FI(I,J)=YV(J) C 190 CONTINUE C C SHIFT X BY DELTA X C E=DLT/DLX E2=DLT*DLT*0.5D0/DLX C C DO 290 J = 1,NDU Ul=(UXI(J)/VOL)**2 C DO 270 I=1,NDX 270 XV(I) = UYI(I,J) N=NDX C C DERIVATIVE OF UY W.R.T. X C CALL SLIP(N) C C DO 240 I=1,NDX C U2=(UYI(I,J)/VOL)**2 RR=1.D0+U1+U2 R=DSQRT(RR) F=R*RR S=(UXI(J)/VOL)*(UYI(I,J)/VOL)/F C C DERIVATIVE OVER X-GRID, SO DIVIDE BY C DLX (REAL DISTANCE). C Appendix B: Listing of the Vlasov Code 190 DVX=-S*SS(I)/DLX DVT=S*EMI*(EY(I)-(UXI(J)/R)*BZ(I)/VOL) C C THE SHIFT IN X C SH(I)=-E*UXI(J)/R 1 +E2*((UXI(J)/R)*DVX-DVT) 240 XV(I)=FI(I,J) C N=NDX CALL SHIP(N) C DO 280 I=1,NDX 280 FI(I,J)=YV(I) C 290 CONTINUE RETURN END SUBROUTINE VLSVION C C 16 Jun 88 C— Shift FI* by DELTA U = DELTA Ul. + DELTA U2 C— and DELTA X (non-relativistic) C #include "vine" #include "mine" #include "cine" C (DLT/DUI FACTOR FROM SHIF) F = EMI*DLT/DUI C U - SHIFT DO 190 I = 1, NDX DO 160 J = 1, NDU C V = U/R 150 SH(J)= F*(EX(I) + UYI(I,l)*BZ(I)/VOL) 160 XV(J)= FI(I,J) N = NDU CALL SHIF(N) DO 170 J = 1, NDU 170 FI(I,J) = YV(J) 190 CONTINUE C X - SHIFT DO 290 J = 1, NDU DO 240 I = 1, NDX C V = U/R 230 SH(I)= - UXI(J)*DLT/DLX 240 XV(I)= FI(I,J) N = NDX CALL SHIP(N) Appendix B: Listing of the Vlasov Code 191 DO 260 1 = 1, NDX 260 FI(I,J) = YV(I) 290 CONTINUE C RETURN END C S U B R O U T I N E L O N G C C Two species. C Longitudinal electric field EX (esu) from fe(x,v) and fi(x,v). C Also calculate transverse current GE for use in MXWL. C #include "vine" #include "mine" #include "cine" INTEGER NUl,NU2,ONE,MONE,N C ONE=l MONE=-l NU1 = NDU - 1 NU2 = NDU - 2 SI = DFLOAT(ISGN) C CONST. FOR U INTEGR. C A has the same value as for dui/uti A = DUE/UTE/7.5D0 C F = DEN* OMG**2/EME c *** Integrations *** c XV = INT. (Fi(X,V)-Fe(X,V))DUx c ge = INT. (Fi*Vyi-Fe*Vye) DUx SI = 0.D0 S2 = 0.D0 Q = 0.D0 Rl = 0.D0 R2 = 0.D0 X = 1.D0 DO 30 I = 1, NDX c FOR V-INTEGRATION, USE c SIMPSON WITH END CORRECTION c S = 2.75D0*((SI*FI(I,1)-FE(I,1))+(SI*FI(I,NDU)-FE(I,NDU))) 2 + ((SI*FI(I,2)-FE(I,2))+(SI*FI(I,NU1)-FE(I,NU1))) 3 - 0.25D0*((SI*FI(I,3)-FE(I,3))+(SI*FI(I,NU2)-FE(I,NU2))) C C Ul=(UXE(l)/VOL)**2 U2=(UYE(I,l)/VOL)**2 R1=DSQRT(1+U1+U2) Appendix B: Listing of the Vlasov Code 192 C Ul=(UXE(2)/VOL)**2 U2=(UYE(I,2)/VOL)**2 R2=DSQRT(1+U1+U2) C Ul=(UXE(3)/VOL)**2 U2=(UYE(I,3)/VOL)**2 R3=DSQRT(1+U1+U2) C U1=(UXE(NDU)/V0L)**2 U2=(UYE(I,NDU)/VOL)**2 RDU=DSQRT(1+U1+U2) C U1 = (UXE(NU1)/V0L)**2 U2=(UYE(I,NU1)/V0L)**2 RU1=DSQRT(1+U1+U2) C Ul=(UXE(NU2)/VOL)**2 U2=(UYE(I,NU2)/VOL)**2 RU2=DSQRT(1+U1+U2) C C IF(ABS(MIE).LT.100.D0) THEN C U1=(UXI(1)/V0L)**2 U2=(UYI(I,l)/VOL)**2 RI1=DSQRT(1+U1+U2) C Ul=(UXI(2)/VOL)**2 U2=(UYI(I,2)/VOL)**2 RI2=DSQRT(1+U1+U2) C Ul=(UXI(3)/VOL)**2 U2=(UYI(I,3)/VOL)**2 RI3=DSQRT(1+U1+U2) C U1=(UXI(NDU)/V0L)**2 U2=(UYI(I,NDU)/VOL)**2 RIDU=DSQRT(1+U1+U2) C U1=(UXI(NU1)/V0L)**2 U2=(UYI(I,NU1)/V0L)**2 RIU1=DSQRT(1+U1+U2) C Ul=(UXI(NU2)/VOL)**2 U2=(UY1(I,NU2)/V0L)**2 RIU2=DSQRT(1+U1+U2) C Q = 2.75D0*(FE(I,1)*UYE(I,1)/R1-SI*FI(I,1)*UY1(I,1)/RI1 2 + FE(I,NDU)*UYE(I,NDU)/RDU -SI*FI(I,NDU)*UYI(I,NDU)/RIDU) 3 + (FE(I,2)*UYE(I,2)/R2-SI*FI(I12)*UYI(I,2)/RI2 4 + FE(I,NU1)*UYE(I,NU1)/RUI-SI*FI(I)NU1)*UYI(I)NU1)/RIU1) Appendix B: Listing of the Vlasov Code 193 5 - 0.25D0*(FE(I,3)*UYE(I13)/R3-SI*FI(I,3)*UYI(II3)/RI3 6 + FE(I,NU2)*UYE(I,NU2)/RU2-SI*FI(I,NU2)*UYI(I,NU2)/RJU2) C Y - l.DO C DO 20 J = 2, NU1 Ul = (UXE(J)/VOL)**2 U2 = (UYE(I,J)/VOL)**2 RE = DSQRT(1+U1+U2) Ul = (UXI(J)/VOL)**2 U2 = (UYI(I,J)/VOL)**2 RI = DSQRT(1+U1+U2) S = S + (7.DO + Y)*(SI*FI(I,J)-FE(I,J)) 20 Y = l.DO - Y ELSE C C C Q = 2.75D0*(FE(I)1)*UYE(I,1)/R1-FI(I,1)*UYI(I,1) 2 + FE(I,NDU)*UYE(I,NDU)/RDU - FI(I,NDU)*UYI(I,1)) 3 + (FE(I,2)*UYE(I,2)/R2-FI(I,2)*UYI(I,1) 4 + FE(I1NU1)*UYE(I,NU1)/RU1-FI(I1NU1)*UYI(I,1)) 5 - 0.25D0*(FE(I,3)*UYE(I,3)/R3-FI(II3)*UYI(I,1) 6 + FE(I)NU2)*UYE(I)NU2)/RU2-FI(I,NU2)*UYI(I,1)) C C Y = 1.0D0 C DO 25 J = 2, NU1 Ul = (UXE(J)/VOL)**2 U2 = (UYE(I,J)/VOL)**2 RR = DSQRT(1+U1+U2) C C INTEGRATE FV*UY TO FIND TRANSVERSE CURRENT C (IONS NON-RELATIVISTIC) C C C C INTEGRATE FV*VY TO FIND TRANSVERSE CURRENT (IONS RELATIVISTIC) Q = Q + (7.DO + Y)*(FE(I,J)*UYE(I,J)/RE-SI*FI(I,J)*UYI(I)J)/RI) C C C INTEGRATE FV TO FIND CHARGE DENSITY C Q = Q + (7.D0 + Y)*(FE(I,J)*UYE(I,J)/RR-FI(I,J)*UYI(I,1)) C C C INTEGRATE FV TO FIND CHARGE DENSITY S 25 Y = S + (7.DO + Y)*(FI(I,J)-FE(I,J)) = 1.D0-Y C Appendix B: Listing of the Vlasov Code 194 ENDIF C GE(I)= A*Q IF(ISGN.EQ.-l) THEN XV(I)= 2+A*S ELSE XV(I)=A*S ENDIF C C C R = R + (l.ODO + X)*XV(I) 30 X = l.DO - X C C MAKE NET CHARGE ZERO C CHARGE DENSITY CORRECTION C R = R/1.5D0/DFLOAT(NDX) C IF(ISIGN.EQ.-l) THEN C DO 40 1=1, NDX 40 XV(I)=2.D0*(XV(I)-R)/(2.DO-R) C ELSE C DO 45 I = 1, NDX XV(I)=(XV(I)-R)/(1.D0-R) 45 CONTINUE C ENDIF C 2.*** EX *** C C INTEGRATE: PER. SPLINE *** N = NDX CALL SINT(N) C ' AVERAGE FIELD E = 0.D0 DO 50 I = 1, NDX 50 E = E + YV(I) E = E/DFLOAT(NDX) C MAKE AVERAGE EX = 0 C CONVERT TO E.S.U. DO 60 I = 1, NDX 60 EX(I)= (YV(I) - E)*F*DLX C RETURN END C Appendix B: Listing of the Vlasov Code 195 S U B R O U T I N E L O N G O N E c C One species. C Longitudinal electric field EX (esu) from fe(x,v) C Also calculate transverse current GE for use in MXWL. C ' #include "vine" #include "mine" #include "cine" INTEGER NUl,NU2,ONE,MONE C ONE=l MONE=-l NU1 =NDU - 1 NU2 =NDU - 2 C CONST. FOR U INTEGR. A = DUE/UTE/7.5D0 F = DEN*OMG**2/EME c *** integrate *** c ge = INT. Fe(X,V)*VY(X,V) DU c XV = INT. Fe(X,V) DU pp = O.DO R = O.DO S = O.DO Q = O.DO X - l.DO DO 30 I = 1, NDX c FOR V-INTEGRATION, USE c SIMPSON WITH END CORRECTION c S = 2.75D0*(FE(I,1)+FE(I,NDU)) • 2 + (FE(I,2)+FE(I,NU1)) 3 - 0.25D0*(FE(I,3)+FE(I,NU2)) C C Ul=(UXE(l)/VOL)**2 U2=(UYE(I,l)/VOL)**2 R1=DSQRT(1+U1+U2) C Ul=(UXE(2)/VOL)**2 U2=(UYE(I,2)/VOL)**2 R2=DSQRT(1+U1+U2) C Ul=(UXE(3)/VOL)**2 U2=(UYE(I,3)/VOL)**2 R3=DSQRT(1+U1+U2) C Ul=(UXE(NDU)/VOL)**2 U2=(UYE(I,NDU)/VOL)**2 RDU=DSQRT(1+U1+U2) C U1=(UXE(NU1)/V0L)**2 Appendix B: Listing of the Vlasov Code 196 U2=(UYE(I,NU1)/V0L)**2 RU1=DSQRT(1+U1+U2) C Ul = (UXE(NU2)/VOL)**2 U2=(UYE(I,NU2)/VOL)**2 RU2=DSQRT(1+U1+U2) C C C Q = 2.75D0*(FE(I,1)*UYE(I,1)/R1+FE(I,NDU)*UYE(I,NDU)/RDU) 4 + (FE(I,2)*UYE(I,2)/R2+FE(I,NU1)*UYE(I,NU1)/RU1) 5 - 0.25D0*(FE(I,3)*UYE(I13)/R3+FE(I,NU2)*UYE(I,NU2)/RU2) C C Y = l.ODO C DO 20 J = 2, NU1 Ul = (UXE(J)/VOL)**2 U2 = (UYE(I,J)/VOL)**2 RR = DSQRT(1+U1+U2) C C INTEGRATE FE*VY TO OBTAIN TRANSVERSE CURRENT C Q = Q + (7.D0 + Y)*FE(I,J)*UYE(I,J)/RR C C INTEGRATE FE TO OBTAIN CHARGE DENSITY C S = S + (7.DO + Y)*FE(I,J) 20 Y = l.DO - Y C C GE(I)= A*Q XV(I) = l.ODO - A*S X=1.D0-X 30 R=R+(1.D0+X)*XV(I) C C MAKE NET CHARGE ZERO C CHARGE DENSITY CORRECTION C R=R/1.5D0/DFLOAT(NDX) DO 40 1=1,NDX 40 XV(I)=(XV(I)-R)/(1.D0-R) C C 2 *** EX *** C C INTEGRATE: PER. SPLINE *** N=NDX CALL SINT(N) C AVERAGE FIELD E = O.DO DO 50 I = 1, NDX 50 E = E + YV(I) Appendix B: Listing of the Vlasov Code 197 E = E/DFLOAT(NDX) C MAKE AVERAGE EX = 0 C CONVERT TO E.S.U. DO 60 I = 1, NDX 60 EX(I)= (YV(I) - E)*F*DLX RETURN END C C SUBROUTINE MXWL C C Sept 88 C Integrate Maxwell equations C BZ, UY at half time step C EY at next integer time step C #include "vine" #include "mine" REAL*8 F.FIO #include "cine" INTEGER N.M.ONE C F = EME*DLT FIO = EMI*DLT C C UPDATE BZ TO NTS+1/2 C YV = D EY/DLX C DO 10 I = 1, NDX 10 XV(I) = EY(I) N = NDX CALL SLIP(N) C C ALSO FIND BZ(DLT) USING C QUADRATIC FIT. C AZ(X) PREVIOUS STEP'S BZ(X) C DO 20 I = 1, NDX P = BZ(I) - CTX*SS(I) YV(I) = -AZ(I)/8.D0+3.D0*BZ(I)/4.D0+3.D0*P/8.D0 AZ(I) = BZ(I) 20 BZ(I) = P C C C UPDATE UY TO NTS+1/2 DO 30 I = 1, NDX DO 40 J = 1, NDU Ul=(UXE(J)/VOL)**2 U2=(UYE(I,J)/VOL)**2 Appendix B: Listing of the Vlasov Code 198 R=DSQRT(1+U1+U2) UYE(I,J)= UYE(I.J) - F*(EY(I)-(UXE(J)/R)+(YV(I)/VOL)) IF(ISPEC.GT.l) THEN C UPDATE UY1 REL. OR NON-REL C DEPENDING UPON MASS RATIO IF(MIE.LT.IOO.DO) THEN Ul=(UXI(J)/VOL)**2 U2=(UYI(I,J)/VOL)**2 R=DSQRT(1+U1+U2) UYI(I,J)= UYI(I,J) - FIO*(EY(I)-(UXI(J)/R)*(YV(I)/VOL)) ELSE UYI(I,1)= UYI(I,1) - FIO*EY(I) ENDIF ENDIF 40 CONTINUE 30 CONTINUE C UPDATE EY TO NTS + 1 C YV = D BZ/DLX DO 50 1= 1, NDX 50 XV(I) = BZ(I) N = NDX CALL SLIP(N) F=DEN*OMG**2/EME*DLT DO 60 I = 1, NDX 60 EY(I)= EY(I) - CTX*SS(I) 2 + F*GE(I) RETURN END C S U B R O U T I N E S A F T C C 16 Jun 88 C Save Fourier transforms C Calculate energies C #include "vine" #include "mine" REAL*8 ESE,EKEN,RHOK,EYEN #include "cine" INTEGER ONE.N ESE = 0.D0 RHOK = 0.D0 C PIL=8.D0*PIE C C USE YV : FOUR OVERWRITES C F O U R : F F T ONE=l N=NDX Appendix B: Listing of the Vlasov Code 199 DO 30 I = 1, NDX 30 YV(I)= EY(I) CALL FOUR(YV,N,ONE) ONE=l N=NDX C STORE IN T Y J = 1 M = (NT-1)*NFY DO 40 I = 1, NFY J = J + 2 EYEN = DSQRT(YV(J)**2 + YV(J+1)**2) TY(M+I)= EYEN TYE(I,NT)=EYEN*EYEN/(PIL) 40 CONTINUE C C FIND ELECTROSTATIC ENERGIES C DO 90 1=1, NDX 90 YV(I) = EX(I) C CALL FOUR(YV,NDX,ONE) C J = l M = (NT-1)*NFX DO 100 1=1, NX2 3=J+2 RHOK=DSQRT(YV(J)**2 + YV(J+1)**2) EKEN = (RHOK)*(RHOK)/(PIL*1.0d20) lF(I.LE.NFX) THEN TX(M+I) = RHOK EK(I,NT) = EKEN ENDIF ESE = ESE + EKEN RHOK = O.DO EKEN = O.DO 100 CONTINUE ES(NT) = ESE/XMX ESE = O.DO C C DUMP DATA EVERY IDUMP STEPS C IF(NT.EQ.IDUMP) THEN CALL DUMPIT NT = 0 ENDIF C 200 FORMAT(I5,3X,F40.15) RETURN END Appendix B: Listing of the Vlasov Code 200 C SUBROUTINE DUMPIT c c #include "vine" #include "mine" #include "cine" Integer N C C N = NFX'IDUMP VVRITE(13) (TX(I), I = 1, N) N = NFY*IDUMP WRITE(13) (TY(I), I = 1, N) C DO 10 1=1,NT 10 VVRITE(12) ES(I) C WRITE(14) ((TYE(I,J), I = 1, NFY),J = 1,NT) . WRITE(14) ((EK(I,J),I=1,NFX),J = 1,NT) RETURN END B.3 Cubic Spline Routines C SUBROUTINE SHIP(N) C C 16Jun88 C Given periodic function at unit intervals, C XV(N+1) = XV(1), and spline coeff. SP, SQ C interpolate to find YV(X) = XV(X + SH) C #include "vine" #include "mine" #include "cine" Integer M,N C M = N - 1 C AUX. COEFF. SU SU(1)= 3.0D0*(XV(2)-XV(N))/SP(1) C DO 10 I = 2, M SU(I)= (3.0D0*(XV(I+1)-XV(I-1)) 2 - SU(I-1))/SP(I) 10 CONTINUE C AUX. COEFF. ST, SV ST(N)= 1.D0 SV(N)= 0.D0 Appendix B: Listing of the Vlasov Code 201 C D O 2 0 I = l , M J = N -1 ST(J)= SQ(J)*ST(J+1)+SR(J) 20 SV(J)= SQ(J)*SV(J+1)+SU(J) C CALCULATE SLOPES SS SS(N)= (3.0D0*(XV(1)-XV(M))-SV(1) 2 -SV(M))/(ST(1)+ST(M)+4.D0) C DO 30 I = 1, M J = N - I 30 SS(J)= ST(J)*SS(N) + SV(J) C DO 40 I = 1, N D = SH(I) C SHF=INTEGER(K)+ FRAC.(D) K = IDINT(D) IF(D.LT.ODO) K = K - 1 D = D - K DL1 = D *(1.D0-D)**2 DL2 = (1.D0-D) *D**2 DL3 = (l.D0-D)**2 *(1.D0+2.D0*D) DL4 = D**2 *(1.D0+2.D0*(1.D0-D)) C J,K IN RANGE 1 TO N J = MOD(N+I+K-l, N)+ 1 K = MOD(N+I+K, N)+ 1 40 YV(I)= SS(J)*DL1 - SS(K)*DL2 1 + XV(J)*DL3 + XV(K)*DL4 RETURN END C S U B R O U T I N E S H I F ( N ) C C 16 Jun 83 C Given fn. F at unit intervals, 0 at ends C (XV(1)=XV(N)=0 ), and spline coeff. SP.SQ C Shift : YV(X) = XV(X + SH) C #include "vine" #include "mine" #include "cine" INTEGER N1,N,K C Nl = N - 1 C AUX. COEFF. SU SU(1)= 3.0D0*XV(2)/SP(1) C DO 10 I = 2, Nl 10 SU(I)= (3.0D0*(XV(I+1)-XV(I-1))-SU(I-1))/SP(I) Appendix B: Listing of the Vlasov Code 202 C CALCULATE SLOPES -SV(N)= -(3.0D0*XV(N1)+SU(N1))/SP(N) C DO 30 I = 1, Nl K = N - I 30 SV(K)= SQ(K)*SV(K+1) + SU(K) C DO SHIFT DO 80 I = 1, N C SH = J + POS. FRAC. D S = SH(I) J = IDINT(X) IF(S .LT. 0.D0) J = J-l D = S - J Dl = D*(1.D0-D)**2 D2 = (1.0D0 - D)*D**2 D3 = (1.0D0 - D)**2 *(1.0D0 +2.0D0*D) D4 = D**2*(1.D0 + 2.D0*(1.D0 - D)) K = I+J + l IF(K .GT. N .OR. K .LT. 1) GO TO 40 G2 = XV(K) S2 = SV(K) GO TO 50 C OUTSIDE GRID, FNC.SLOPE = 0 40 G2 - O.DO S2 = O.DO 50 K = I + J IF(K .GT. N .OR. K .LT. 1) GO TO 60 GI = XV(K) SI = SV(K) GO TO 70 60 GI = O.DO SI = O.DO 70 CONTINUE 80 YV(I)= S1*D1 - S2*D2 + G1*D3 + G2*D4 RETURN END C S U B R O U T I N E S L I P ( N ) C C 24 Jun 83 C Given periodic function at UNIT intervals, C XV(N+1) = XV(1), and spline coeff. SP,SQ,SR C find SS(N) = df/dx by spline interpolation C #include "vine" #include "mine" #include "cine" INTEGER M,N C Appendix B: Listing of the Vlasov Code 203 M = N - 1 SU(1)= 3.D0*(XV(2)-XV(N))/SP(1) C DO 10 I = 2, M 10 SU(I)= (3.D0*(XV(I+1)-XV(I-1)) 1 - SU(I-l)) /SP(I) C ST(N)= l.DO SV(N) = O.DO C DO 20 I = 1, M 3 = N - I ST(J)= SQ(J)*ST(J+1) + SR(J) 20 SV(3)= SQ(J)*SV(J + 1) + SU(J) C SS(N)= (3.D0*(XV(1)-XV(M))-SV(1)-SV(M)) 1 /(ST(1) + ST(M) + 4.D0) C DO 30 I = 1, M J = N - I 30 SS(J)= ST(J)*SS(N) + SV(J) RETURN END C S U B R O U T I N E S I N T ( N ) C C 16 Jun 88 C Given periodic function at unit intervals, C XV(N + 1) = XV(1), and spline coeff. SP, SQ, SR C Integrate : YV(N) = Int. XV up to N C #include "vine" #include "mine" #include "cine" Integer N,M C M = N - 1 C COEFF. SU SU(1)= 3.0D0*(XV(2) - XV(N))/SP(1) C DO 10 I = 2, M 10 SU(I)= (3.0D0*(XV(I+1)-XV(I-1)) 2 -SU(I-1))/SP(I) C COEFF. ST, SV ST(N)= l.DO SV(N)= O.DO C DO 20 I = 1, M J = N - I Appendix B: Listing of the Vlasov Code 204 ST(J)= SQ(J)*ST(J+1) + SR(J) 20 SV(J)= SQ(J)*SV(J+1) + SU(J) C - SLOPES SS SS(N)= (3.0D0*(XV(1)-XV(M))-SV(1)-SV(M)) 1 /(ST(1) + ST(M) + 4.D0) C 1. DO 30 I = 1, M J = N - I 30 SS(J)= ST(J)*SS(N)'-fSV(J) C *** INTEGRATION YV(1)= (XV(N) + XV(1))/2.D0 1 + (SS(N) - SS(1))/12.D0 DO 40 I = 2, N YV(I)= YV(I-l) + (XV(I-l) + XV(I))/2.D0 1 + (SS(I-l) - SS(I))/12.D0 40 CONTINUE RETURN END B.4 Fast Fourier Transform Routines subroutine four(data,na,is) c 20 May 88 c is= 1 : transform nn reals to nn/2 complex c -1 : transform nn/2 conjugate complex to real c (both for positive frequency half only C-NR400 implicit real*8 (a-h,o-y) dimension data(l) n - nn/2 phi=3.141592653589793d0/dfloat(n) cl - 0.5d0 if(is.eq.l) then c2 = -0.5d0 call fouc(data,nn,l) else c2= 0.5d0 phi = -phi endif yr = -2.0d0*(dsin(0.5d0*phi))**2 yi = dsin(phi) wr = l.OdO + yr wi = yi do 11 i=2, n/2+1 11 = 2*i- 1 12 = il + 1 Appendix B: Listing of the Vlasov Code 205 13 = nn + 3 - i2 14 = i3 + 1 gr = cl*(data(il) + data(i3)) gi = cl*(data(i2) - data(i4)) hr = -c2*(data(i2) + data(i4)) hi = c2*(data(il) - data(i3)) data(il)= gr + wr*hr - wi*hi data(i2)= gi + wr*hi + wi*hr data(i3)= gr - wr*hr + wi*hi data(i4)=-gi + wr*hi + wi*hr wt = wr wr = wr*yr - wi*yi + wr wi = wi*yr + wt*yi + wi 11 continue if(is.eq.l) then gr = data(l) data(l)= gr + data(2) data(2)= gr - data(2) else gr = data(l) data(l) = cl*(gr + data(2)) data(2) = cl*(gr - data(2)) call fouc(data,nn,-l) endif return end subroutine fouc(data,nn,is) c c 19 May 88 c is = 1 : transform of nn/2 complex number pairs c -1 : inverse transform c-NR c implicit real*8 (a-h, o-y) dimension data(l) c j = 1 do 11 i = 1, nn, 2 if (j .gt. i) then tr' = data(j) ti = data(j+l) data(j) = data(i) data(j+l)= data(i+l) data(i) = tr data(i+l)= ti endif m — nn/2 1 if((m.ge."2) .and. (j.gt.m)) then Appendix B: Listing of the Vlasov Code 206 j = j - m m = m/2 goto 1 endif j = j + m 11 continue mmax — 2 2 if (nn. gt. mmax) then ix = 2*mmax phi = 6.28318530717959d0/dfloat(is*mmax) yr = -2.d0*(dsin(0.5d0*phi))**2 yi = dsin(phi) wr = l.OdO wi - O.dO do 13 m = 1, mmax, 2 do 12 i = m,nn,ix j = i + mmax tr = wr*data(j) - wi*data(j+l) ti = wr*data(j+l)+wi*data(j) data(j) = data(i) - tr data(j+l) = data(i+l)-ti data(i) — data(i) + tr data(i+l)= data(i+l)+ti 12 continue wt = wr wr = wr*yr - wi*yi + wr wi = wi*yr + wt*yi + wi 13 continue mmax = ix goto 2 endif return end A P P E N D I X C USER'S G U I D E C l Compiling and Running the Vlasov Code The Vlasov code simulations presented here were run on a Sun 3/60 Workstation. The subroutines of the Vlasov code were kept as separate entities. These were compiled separately and then linked together using a "makefile". The features of this makefile made this very advantageous, since only subroutines that were changed during the de-bugging process were recompiled, after which all the subroutines were linked back together. An example of the makefile for the Vlasov Code is shown below. All that is required after this file (called "makefile") has been produced is to type "make" and the code compiles any changed components and links all the components back together. The resulting object codes were VINI and VRUN. One executes VINI, where the initial system is set up using the values in VLAS.IN. After this, VRUN is executed, iterating through the Vlasov algorithm the number of times requested. m a k e f i l e prog: vini.o spqr.o inpt.o inife.o inifi.o stor.o vrun.o recl.o vlsv.o vlsvi.o vlsvion.o vlsvul.o long.o mxwl.o shif.o ship.o sint.o slip.o four.o fouc.o saft.o dumpit.o longone.o bekl.o f77 vini.o inpt.o spqr.o inife.o inifi.o stor.o ship.o sint.o long.o longone.o vlsvul.o shif.o bekl.o four.o fouc.o -o vini f77 -0 vrun.o recl.o vlsv.o vlsvul.o vlsvi.o vlsvion.o long.o longone.o mxwl.o dumpit.o stor.o shif.o ship.o sint.o slip.o four.o fouc.o saft.o -o vrun vini.o: vini.F 177 -c vini.F spqr.o: spqr.f 207 Appendix C: User's Guide 208 f77 -0 -c spqr.f inpt.o: inpt.F f77 -c inpt.F bekl.o: ../func/bekl.f f77 -c ../func/bekl.f inife.o: inife.F f77 -0 -c inife.F inifi.o: inifi.F f77 -0 -c inifi.F stor.o: stor.F f77 -0 -c stor.F vrun.o: vrun.F f77 -0 -c vrun.F recl.o: recl.F f77 -0 -c recl.F vlsv.o: vlsv.F f77 -0 -c vlsv.F vlsvul.o: vlsvul.F f77 -0 -c vlsvul.F vlsvi.o: vlsvi.F f77 -0 -c vlsvi.F vlsvion.o: vlsvion.F f77 -0 -c vlsvion.F long.o: long.F f77 -0 -c long.F longone.o: longone.F 177 -0 -c longone.F mxwl.o: mxwl.F f77 -0 -c mxwl.F shif.o: shif.F 177 -0 -c shif.F ship.o: ship.F f77 -O -c ship.F sint.o: sint.F f77 -0 -c sint.F slip.o: slip.F f77 -0 -c slip.F four.o: four.f f77 -0-c four.f fouc.o: fouc.f f77 -0 -c fouc.f saft.o: saft.F f77 -0 -c saft.F dumpit.o: dumpit.F f77 -0 -c dumpit.F Appendix C: User's Guide 209 An unfortunate feature of Fortran on the Sun 3/60 was that it did not support the namelist read capability. This meant that input parameters had to be stored unlabelled in a file. The following file is an example of this input file, VLAS.IN, with its corresponding entries. VLAS.IN 32,51,1,0. DO 0.04D0,0.D0,0.000225D0 0.D0,0.D0 1,0.01D0,0.D0 0.D0.0.D0 5,1 500 100 Desciption NDX,NDU, No. of species ISPEC, Ion to electron mass ratio MIE Electron thermal velocity / c UT1, Ion thermal velocity UT2, DEN Plasma slab size in wavelengths XVL, Vosc /c VOC NPA.PAME.PAMI Drift velocities / c : V01,V02 Number of modes to store: NFX.NFY Store energy diagnostics every IDUMP steps Store disribution function(s) and fields every ITF steps C.2 Plotting Routines Plotting routines were written using the graphics package Superplot. This package, which can be obtained from Marc Majka at the U.B.C laboratory for Computational Vision (for free!), gives rudimentary 2-D plotting capabilities that can be integrated into T E X . The 2-D graphics routine simply reads a data file of (x,y) pairs creating a graph of these values with instructions given in an instruction file. In this respect the program is "stand alone" and required a different approach than other types of graphics routines. Because a data file of (x,y) pairs had to be created for every graph made, the plotting routines written for the Vlasov code and EMI were separate for each diagnostic. In this way the computer disk storage would not be overcome with huge amounts of data. In Appendix C: User's Guide 210 fact, it would not be possible to fit both the diagnostic data from a simulation and all the required (x,y) files into the disk storage area in most applications. Superplot has two main commands. Firstly, txyplot takes (x,y) data from a specified file and writes a superplot language version of the resulting plot in a file named plot.file. The instruction set as to how the graph should be made is in file plot.format. A number of these graphing formats are stored, one for each diagnostic. An example of a plot.format file is given below. plot.format x 0 140 y -25.0 0.0 ylabel "LOG(Mode Energy)" xlabel "Time(l/Wp)" xtick 0.0 140 5 ytick -25.0 0.0 1.0 xnumber 0.0 140 10 ynumber -25.0 0.0 5.0 xformat %.lf yformat %g join output plot.file After this sunplot reads the plot.file and produces a graph on any of the Sun win-dows. If a hardcopy plot is required, pips converts the superplot plot.file to a postscript file. The following is a macro (called "plot") to carry out the process of plotting a graph on the screen of the Sun. If the (x,y) pairs are in a file called "data", then one simply carries out the command "plot data". plot #! /bin/sh if (test $# -It 1) then echo Name of data file ? read filename set Sfilename echo $1 fi txyplot plot.format $1 sunplot plot.file Appendix C: User's Guide 211 Table 20 names and describes the plotting routines developed for the Vlasov code. This table is followed by a listing of the routines. After this, table 21 lists the plotting routines developed for EMI after which follows the listing of the named routines. These routines read and analyze the unformatted diagnostic data from simulations and write the resulting products in (x,y) pairs in a file called "data". VLASOV DIAGNOSTIC PLOTTING ROUTINES Routine Description ese Electrostatic energy over time exf Time average of the electrostatic energy modes ext Electrostatic energy per mode versus time fdist Distribution function(s) at a time interval fmts Writes the distribution function(s) to MTS for 3-D plotting field x-grid values of EX,EY,GE,BZ at a time interval modeav Time average of the transverse electric field modes plotmode Energies per mode versus time (electrostatic and EY) slope Slope of a graph using a least-squares fit timef Transverse electric field modes versus time uydist Electron/ion transverse velocity at a time interval fmts Writes the electron transverse velocity to MTS for 3-D plotting Table 20: Vlasov diagnostic routines. ESE C C Outputs electrostatic energy versus time C c real*8 es,dlt,wp,time,omg,den integer n open(12,me='ESE,,form='unformatted') open(10,file='data') C C read in initial parameters C read(12)idump,nts,ntm,dlt,den,omg wp=omg*dsqrt(den) write(6,*)idump,nts,ntm,dlt,wp Appendix C: User's Guide 212 ntm=ntm+nts nts=nts+l C C output electrostatic energy C versus time C do 10 i=nts,ntm time=dfioat(i)*dlt*wp read(12)es write(10,*)float(time),float(es) 10 continue end E X F C : C Output average of each Ex mode C C real*8 x(200000),y,dlt,omg,den real*4 stl(20) integer idumpjndx.ndu.nts.ntm.itfe.nfxjnfy open(10,file=:7usr/Pe,;er/datav/data') open(13,file=7usr/peter/datav/LRMODES',form=,unformatted') C C read in initial parameters C read(13)idump,nts,ntm,nfx,nfy,dlt,den,omg write(6,*)idump,nts,ntm,nfx,nfy,dlt,den,omg C C read in Fourier transform C data of Ex and Ey C j=ntm/idump ml=l m2=l nl=0 n2=0 do 5 i=l j nl=nl+nfx*idump read(13)(x(ii),ii=ml,nl) ml=ml+nfx*idump n2=n2+nfy*idump read(13)(y,ii=m2,n2) m2=m2+nfy*idump 5 continue C ntm=ntm+nts nts=nts+l Appendix C: User's Guide 213 C C find the average o f each Ex C mode and write result t o data C do 10 j = l,nfx jj=j do 20 i=nts,ntm stl(j)=stl(j)+x(jj) jj=jj+nfx 20 continue 10 continue write(10,*)0,0.0 do 25 j=l,nfx stl(j)=stl(j)/dfioat(ntm-nts) 25 write(10,*)j,stl(j) end EXT c C Output an electrostatic mode versus time C C real*8 x(200000),y,dlt,wp,time,den,omg integer nts,nfx,nfy,ntm open(13,file=7usr/peter/datav/LRMODES,,form=,unformatted') open(10,file=7usr/peter/datav/data') C C read initial parameters C -read(13)idump,nts,ntm,nfx,nfy,dlt,den,omg wp=omg*dsqrt(den) j=ntm/idump write(6,*)idump,nts,ntm,nfx,nfyj C C read in FFT of Ex, Ey C ml=l m2=l nl=0 n2=0 do 5 i= 1 j nl=nl+nfx*idump read(13)(x(ii),ii=ml,nl) ml=ml+nfx*idump n2=n2+nfy*idump read(13)(y,ii=m2,n2) m2=m2+nfy*idump 5 continue C Appendix C: User's Guide 214 ntm=ntm+nts nts=nts+l write(6,*)ntm C C output requested mode versus time C write(6,*)"Fourier Component ?" write(6,*)"betweea" ,l,"and" ,nfx. read(5,*)ifc j=ifc do 10 i=nts,ntm time=dfloat(i)*dlt*wp write(10,*)float(time)1float(x(j)) " 10 j=j+nfx F D I S T C C Output distribution function(s) over x-grid or C Ux-grid at requested time interval C-C real*8 fe(1024,301),fi(1024,301) integer ndx.ndu.ntSjntirijitfjt^te^pec C C open( 11, file='. ./datav/ FE ',form= 'unformatted') open(10)file='../datav/data') read in initial parameters read(ll)ndx,ndu,nts,ntm,itf,ispec write(6,*)ndx,ndu,nts,ntm,itf,ispec vvrite(6,*)"Number of species ",ispec write(6,*)"NDX =",ndx write(6,*)"NDU =",ndu write(6,*)"Time Step Range =",nts,"-?",ntm.. write(6*)"ITFE =",itf time interval request te=(nts+ntm)/itf ti=nts/itf+l write(6,*)"Which time interval?" write(6,*)"must be between",ti,"and",te read(5,*)itt do 10 i=ti,itt read(ll)((fe(iijj),ii=l,ndx),jj=l,ndu) if(ispec.gt.l) read(ll)((fi(iijj)Jii=l,ndx),jj=-l,n'du) write(6,*)fe(l,l) 10 write(6,*)" Fetching data of time internal-!"ji Appendix C: User's Guide 215 C if(ispec.gt.l) then write(6*)" " write(6,*)"==================== write(6,*)" " write(6,*)" Do you want FE (1)" write(6,*)" FI (2)" write(6,*)" or FE + FI (3) ?" write(6,*)" " write(6,*)" ========== ========== write(6,*)" " read(5,*)ifunc else ifunc=l endif C C 15 write(6,*)" Which component of the distribution?;" write(6,*)"Velocity distribution (1)" write(6,*)" Position distribution (2)" read(5,*)iij if(iij.gt.2.or.iij.lt.l) then write(6,*)"Incorrect Answer, try again:" -write(6,*)" " write(6,*)" =-=-=-=-=-=-=-=-=-=-=-=-=-=-= write(6,*)" " go to 15 endif C C C if(iij.eq.l) then 17 if(ifunc.eq.l) write(6,*)"Position Component of fe(x,u) ?" if(ifunc.eq.2) write(6,*)"Position Component of fl(x,u) ?" if(ifunc.eq.3) write(6,*)"Position Component of fe(x,u)+fi(x,u) ?" write(6,*)"Must be between 1 and",ndx read(5,*)iii if(iii.gt.ndx.or.iii.lt.1) then write(6,*)"Incorrect Answer, tryvagain:" write(6 *)" " write(6,*)" =-=-=-=-=-=-=-=-=-=-=-=-=-=-= write(6,*)"' " go to 17 endif do 20 i=l,ndu if(ifunc.eq.l) write(10,*)i,float(fe(iii,i)) if(ifunc.eq.2) write(10,*)i,float(fi(iii,i)) 20 if(ifunc.eq.3) write(10,*)i,fioat(fe(iii,i)+fi(iii,i)) endif C C C Appendix C: User's Guide 216 if(iij.eq.'2) then 45 if(ifunc.eq.l) write(6,*)" Which velocity component of fe(x,u) ?" if(ifunc.eq.2) write(6,*)" Which velocity component of fi(x,u) ?" if(ifunc.eq.3) write(6,*)" Which velocity component of fe(x,u)+fi(x,u) ?" write(6,*)"Must be between 1 and",ndu read(5,*)iv if(iv.gt.ndu.or.iv.lt.l) then write(6,*)"Incorrect Answer, try again:" write(6,*)" " write(6,*)" =-=-=-=-=-=-=-=-=-=-=-=-=-=-=" write(6,*)" " go to 45 endif do 50 i=l,ndx if(ifunc.eq.l) write(10,*)i,float(fe(i,iv)) if(ifunc.eq.2) write(10,*)i,float(fi(i,iv)) 50 if(ifunc.eq.3) write(10,*)i,float(fe(i,iv)+fi(i,iv)) endif C end F M T S c : C Write distribution function(s) at a requested time interval C to an MTS account for 3-D plotting C C output in file "sunstuffl" C C real*8 fe(10241 01),fi(1024,101),ux(101) integer ndx,ndu,nts,ntm,itfe,ti,te open(ll,nle='../datav/FE',form='unformatted') open(20,file='../datav/sunstuffl') C read(ll)ndx,ndu,nts,ntm,itfe,ispec write(6,*)"NDX =",ndx write(6,*)"NDU =",ndu write(6,*)"Time Step Range =" ,nts,"-" ,ntm+nts write(6,*)"ITFE =",itfe write(6,*)"Time interval to be transferred?" read(5,*)itime te=(nts+ntm)/itfe ti=nts/itfe write(6, *)" Dumps" ,ti," to", te C ti=l do 30 ij—ti,itime Appendix C: User's Guide 217 read(M)|^uy.i),vii=l,ndx),jj=l,ndu) if(ispec.gtU!)i uead(ll)((fi(ii jj),ii=l,ndx) jj=l,ndu) 30 continue if(ispec.glr.li):; write(20,100)((fe(iijj)+fi(ii,jj),ii=l,ndx),ij=l)ndu) write(2aVMQ$(i'(lfe(:ii jj) ,ii= 1 ,ndx) jj=l,ndu) C 100 format(USf:13») end FIELD C-c c c-c c c c c c c 10 20 c Output fields Ey, Bz, Ex or the transvenae current GE over x-grids real*8 sttwfip)24) integer. nits;,ptrEn,itf,ndx,ti,time charactaar*1!0 fil open^r^e^'skta') wri;te(6*)riMfir file to view, EY, BZ, GE or EX ?" read(5-,*);iil open(ll1fife=ffii,form='unformatted') re ad (li)teta,Bteii,i t f, ndx t i=(Etelr!Efis)) / i f wriite^^yWhich time interval," write( VjTiSetween", 1 " and" , t i , " ? " read(5,*))ltanifi: do 10: ii=l,j!lranie r e adi(i 1) )'(j's feufE(l i ) , i i = 1 ,ndx) do 20 i=l,ndx write(4L*).fl:Qat(i),fioat(stuff(ii)) read in initial parameters output matrix at time interval requested MODEAV/ C C OutpjJt the; average of. easdis Ey mode C C reat*S:xry(,150000),dlt,den,omg real'S atui(:>29) Appendix C: User's Guide 218 integer n,nts,ntm,nfx,nfy open(ll,file=7usr/peter/datav/LRMODES',form='unformatted') open(10,file=7usr/peter/datav/data') C C read initial parameters C read(ll)idump,nts,ntm,nfx,nfy,dlt,den,omg write(6,*)idump,nts,ntm,nfx,nfy,dlt,den,omg j=ntm/idump write(6,*)j ntm=ntm+nts C read FFT of Ex, Ey nts=nts+l nl=0 n2=0 ml=l m2=l do 5 i=l j nl=nl+nfx*idump read(ll)(x,ii=ml,nl) ml=ml+nfx*idump n2=n2+nfy*idump read(ll)(y(ii),ii=m2,n2) m2=m2+nfy*idump 5 continue C Calculate and output average do 65 j = l,nfy jj=j do 60 i=nts,ntm str(j)=str(j)+y(jj) jj=jj+nfy 60 continue 65 continue write(10,*)0,0.0 do 70 j=l,nfy str(j)=str(j)/dfioat(ntm-nts) 70 write(10,*)j,float(str(j)) C PLOTMODE C C Output electrostatic or transverse electric field C mode energies versus time. LOG, LN or regular. C c real*8 ek(20,14000),ey(20,14000),convert,lnten,dlt,wp,time,omg,den integer n,idump,nts,ntm,nfx,nfy,ipos,nm,ns open(14,file='EK',form='unformatted') Appendix C: User's Guide 219 open(10,file='data') C read initial parameters read(14)idump,nts,ntm,nfx,nfy,dlt,den,omg wp=omg*dsqrt(den) n=ntm/idump write(6,*)idump,nts,ntm,nfx,nfy,n,dlt,wp convert=0.434294481903d0 is=l id=idump C read in Ey, Ex energies C • per mode do 10 i=l,n : read(14)((ey(iijj),ii=l,nfy)jj=is,id) read(14)((ek(ii,jj),ii=l,nfx)jj=is,id) is=is+idump id=;id+idump 10 continue nm=ntm+nts ns=nts+l write(6,*)"EXE (1) or EYE (2) ?" read(5,*)imat if(imat.eq.l) then "write(6,*)"Which Fourier component of the EXE matrix ?" write(6,*)"between" ,l,"and" ,nfx read(5,*)ipos write(6,*)"LOG (1), LN (2) or plain(3) graph ?" read(5,*)iform j = l do 20 i=ns,nm time=dfl.oat(i)*dlt*wp C C lnten=log(ek(ipos j)) if(iform.eq>I) write(10,*)float(time),noat(lnten*convert) if'(iform.eq.2) write(10,*)float(time),float(lnten) if(iform.eq.3) write(10,*)float(time),float(ek(iposj)) j=j.+ l 20 continue endif if(imat.eq.2) then write(6,*)"Which Fourier component of the EYE matrix ?" write(6 *)"between",l,"and",nfy read(5,*)ipos write(6,T LOG (!). L N (2) or plain(3) graph ?" read(5,*)iform j=l do 30 i=ns,nm time=dfloat(i)*dlt*wp C lnten=log(ey(ipos j)) Appendix C: User's Guide 220 30 C endif end if(iform.eq.l) write(10 *)float(time),float(lnten*convert) if(iform.eq.2) write(10,*)float(time),fioat(lnten) if(iform.eq,3) write(10,*)float(time),float(ey(iposj)) j=j+l continue S L O P E C-c c -c c c c 5 C C Least-squares slope determination (Ex mode energy) real*8 ek(20,14000),ey(20,14000))tt,lten,x(14000),y(14000),den,omg real*8 a,b,c,d,xsq,top,bot,dlt,wp,ti,tf,time integer n,idump,nts,ntm,nfx,nfy,ipos,nm,ns open(14,file='EK',form='unformatted') read initial parameters read(14)idump,nts,ntrn,nfx,nfy,dlt,den,omg n=ntm/idurnp wp=omg*dsqrt(den) write(6,*)idump,nts,ntm,nfx,n,dlt,wp is—1 id=idump nm=ntm+nts ns=nts+l read Ey, Ex energies per mode do 5 i=l,n read(14)((ey(ii,jj),ii=l,nfy)jj=is,id) read(14)((ek(ii,jj),ii=l,nfx) jj=is,id) is=is+idump id=id+idump continue request mode and time range for slope calculation write(6,*)ek(l,2) write(6,*)"position of mode in array ?" read(5,*)ip write(6,*)"Start and end timesteps for slope evaluation" write(6,*)"( In units of W~l )" read(5,*)ti,tf jj=0 kk=0 Least-squares do 10 i=ns,nm time=dfloat(i)*dlt*wp if(time.gt.tf) go to 35 Appendix C: User's Guide 221 if(time.lt.ti) then jj=jj+l go to 10 endif C C ***Need to take the natural log in this case C kk=kk+l x(i-jj)=time y(i-jj)=log(ek(ip,i)) 10 continue 35 continue write(6,*)"kk",kk do 50 j=l,kk b=b+x(j)*y(J) c=c+x(j) d=d+y(j) xsq=xsq+x(j)**2 50 continue top=dfloat(kk)*b-c*d bot=dfloat(kk)*xsq-c**2 C a=top/bot write(6,*)"The approximate slope is = " write(6,*)a write(6',*)"The approximate imaginary frequency is=" a=a/2.D0 write(6,*)a C end T I M E F C ; • C Outputs transverse field modes versus time C C real*8 x,y(200000),dlt,omg,den,wp,t integer nts,ntm,nfx,nfy open(ll,nle=7usr/peter/datav/LRMODES',form='unformatted') open(10,file=7usr/Peter/clatav/data') C read initial parameters read(ll)idump,nts,ntm,nfx,nfy,dlt,den,omg write(6,*)idump)nts,ntm,nfx,nfy,dlt,den,omg wp=omg*dsqrt(den) j=ntm/idump ntm=ntm+nts nts=nts+l C read FFT of Ex, Ey ml=l Appendix C: User's Guide 222 m2=l nl=0 n2=0 do 5 i=l j nl=nl+nfx*idump read(ll)(x,ii=ml,nl) ml=ml+nfx*idump n2=n2+nfy*idump read(ll)(y(ii),ii=m2,n2) m2=m2+nfy*idump 5 continue C output mode versus time C of Ey write(6,*)"Fourier Component ?" write(6,*)"between",l,"and",nfy read(5,*)ifc j=ifc do 100 i=nts,ntm t=dfloat(i)*dlt*wp write(10,*)float(t),float(y(j)) 100 j=j+nfy C end U Y D I S T C C Output electron and/or ion transverse velocity over C x-grid or Ux-grid at requested time interval C C . real*8 uye(1024J301),uyi(1024,301) integer ndx,ndu,nts,ntm,itf,ti,te,ispec open(ll,file='../datav/UY,,form='unformatted') open(10,file='../datav/data') C read initial parameters read(ll)ndx)ndu)nts,ntm,itf,ispec write(6,*)ndx,ndu)nts)ntm,itf1ispec write(6,*)"Number of species ",ispec write(6,*)"NDX =",ndx write(6,*)"NDU =",ndu write(6,*)"Time Step Range =",nts,"-",ntm write(6,*)"ITFE =",itf te=(nts+ntm)/itf ti=nts/itf+l write(6,*)"Which time interval?" write(6,*)"must be between",ti,"and",te read(5,*)itt C read Uy,e; Uy,i Appendix C: User's Guide 223 do 10 i=ti,itt read(ll)((uye(ii,jj),ii=l,ndx)jj = l,ndu) if(ispec.gt.l) read(ll)((uyi(ii,jj),ii=l,ndx),jj=l,ndu) 10 write(6,*)" Fetching data of time interval" ,i C if(ispec.gt.l) then write(6,*)" " Write(6,*)"============ =============" write(6,*)" " write(6,*)" Do you want UYE (1)" write(6 *)" UYI (2)" write(6,*)" or UYE + UYI (3) ?" write(6,*)" " write(6,*)" === ======= = ==============" write(6,*)" " read(5,*)ifunc else ifunc=l endif C C 15 write(6,*)" Which component of the distribution?;" write(6,*)" Velocity distribution (1)" write(6,*)"Position distribution (2)" read(5,*)iij if(iij.gt.2.or.iij.lt.l) then write(6,*)"Incorrect Answer, try again:" write(6,*)" " write(6,*)"=-=-=-=-=-=-=-=-=-=-=-=-=-=-=" write(6,*)" " go to 15 endif C C C if(iij.eq.l) then 17 if(ifunc.eq.l) write(6,*)"Position Component of uye(x,u) ?" if(ifunc.eq.2) write(6,*)" Position Component of uyi(x,u) ?" if(ifunc.eq.3) write(6,*)" Position Component of uye(x,u)+uyi(x,u) ?" write(6,*)"Must be between 1 and",ndx read(5,*)iii if(iii.gt.ndx.or.iii.lt.l) then write(6,*)"Incorrect Answer, try again:" write(6,*)" " write^,*)" =-=-=-=-=-=-=-=-=-=-=-=-=-=-=" write(6,*)" " go to 17 endif do 20 i=l,ndu if(ifunc.eq.l) write(10,*)i,float(uye(iii,i)) if(ifunc.eq.2) write( 10,*)i,float(uyi(iii,i)) 20 if(ifunc.eq.3) write(10,*)i,float(uye(iii,i)+uyi(iii,i)) Appendix C: User's Guide 224 endif C C C if(iij.eq.2) then 45 if(ifunc.eq.l) write(6,*)"Which velocity component of uye(x,u) ?" if(ifunc.eq.2) write(6,*)" Which velocity component of uyi(x,u) ?" if(ifunc.eq.3) write(6,*)" Which velocity component of uye(x,u)+uyi(x,u) ?" write(6,*)"Must be between 1 and",ndu read(5,*)iv if(iv.gt.ndu.or.iv.lt.l) then write(6,*)"Incorrect Answer, try again:" write(6,*)" " write(6,*)" =-=-=-=-=-=-=-=-=-=-=-=-=-=-=" write(6,*)" " go to 45 endif do 50 i=l,ndx if(ifunc.eq.l) write(10,*)i,float(uye(i,iv)) if(ifunc.eq.2) write(10,*)i,float(uyi(i,iv)) 50 if(ifunc.eq.3) write(10,*)i,float(uye(i,iv)+uyi(i,iv)) endif C end UYMTS C C Write electron transverse velocity to C MTS for 3-D plotting C real*8 uye(1024,101),uyi(1024,101),ux(101) integer ndx^du.nts.ntrrijitfe.tijte open(ll,file='../datav/UY,,form= 'unformatted') open(20,file='../datav/sunstuffl') C read initial parameters read(ll)ndx,ndu,nts,ntm,itfe,ispec write(6,*)"NDX =",ndx write(6,*)"NDU =",ndu write(6,*)"Time Step Range =",nts,"-",ntm+nts write(6,*)"ITFE =",itfe te=(nts+ntm)/itfe ti=nts/itfe write(6,*)" Dumps" ,ti," to" ,te C write(6,*)"Time interval to write?" read(5,*)itime C output Uye at time interval ti=l do 30 ij=ti,itime Appendix C: User's Guide 225 k=k+l read(ll)((uye(ii,jj),ii=l,ndx)jj=l,ndu) if(ispec.gt.l) read(ll)((uyi(iijj),ii=l,ndx),jj=l,ndu) 30 continue write(20,*)((float(uye(ii,jj)),ii=l,ndx)jj=l,ndu) C end EMI DIAGNOSTIC PLOTTING ROUTINES Routine Description phase Phase space (x,v) at specific time interval slope Slope determination by least squares plotmode Electrostatic mode energy versus time plotek Electrostatic mode distribution at a time interval ploteka Time average electrostatic mode distribution plote Electrostatic or transverse field energy versus time fv Particle velocity distribution fx Particle spatial distribution field Grid point values of charge density, potential and Ex timef Transverse field mode versus time modey Transverse fields' mode distribution at a time interval modeyav Transverse fields' average mode distribution Table 21: EMI diagnostic routines. PHASE C-C This program reads phase space data from XVX1 and XVX2 C (phase space data from EMI and EM1B). It deposits C the data in separate files for each time interval C for use with Superplot. C real*8 x,vx integer a,n,m,k,bb,ng,nt,ixvx,nl,n2 character*10 b open( 1 ,file='X VX 1 ',form= 'unformatted') open(2,file='XVX2',form='unformatted') open(7,file='datar) open(8,file='data2') Appendix C: User's Guide 226 open open open open open open open open open open open open open open open open open open open C (9,file= (10,file= (ll,file= (12,file= (13,file= (14,file= (15,file= (16,file= (17,file= (18,file= (19,filer (20,file= (21,file= (22,file= (23,file= (24,file= (25,file= (26,file= (27,file= 'data3') = 'data4') = 'data5') ='data6') ='data7') ='data8') ='data9') ='datalO') ='datall') ='datal2') = 'datal3') = 'datal4') = 'datal5') ='datal6') = 'datal7') ='datal8') ='datal9') = 'data20') ='data21') read(l)ng,nt,ixvx,nl,n2 bb=0 write(6,*)"Nuniber of species to plot?" read(5,*)k if(k.eq.l) then write(6,*)"Which file; XVX1 or XVX2 ?" read(5,*)b if(b.eq.."XVXl") then bb=l n=nl endif if(b.eq."XVX2") then bb=2 n=n2 endif if(k.eq.O) then write(6,*)"File must be XVX1 or XVX2" go to 5 endif endif C C Skip a requested number of time intervals write(6,*)"Skip how many initial dumps ?" write(6,*)"between",0,"and",nt/ixvx read(5,*)m write(6,*)"*** Writing data ***" write(6,*)" " write(6,*)" " a=6 if(k.eq.2) then n=nl do 30 j = l,m Appendix C: User's Guide 227 do 10 i=l,n read(l)x,vx read(2)x,vx 10 continue 30 continue endif if(k.eq.l) then 32 do34j=l,m do 36 i=l,n read(b)x,vx 36 continue 34 continue endif kount=0 if(k.lt.2) go to 55 do 50 j=l,21 a=a+l write(6,*)"Writing to data file",a-6," . . ." do 40 i=l,n read(l)x,vx write(a,60)x,vx read(2)x,vx write(a,60)x,vx 40 continue 50 continue 55 a=6 do 56 j=l,21 a=a+l write(6,*)" Writing to data file" ,a-6," . . ." do 58 i=l,n read(bb)x,vx vvrite(a,60)x,vx 58 continue 56 continue 60 format(2fl0.5) write(6,*)" " write(6,*)" " write(6,*)"Finished !" C end SLOPE C C Determines slope by least-squares. C electrostatic energy C C real*8 h(8,10000),dt,tt,lten,x(10000),y(10000) real*8 a,b,c,d,xsq,tf,ti,top,bot,wp Appendix C: User's Guide 228 integer ih,nks,it,ip open(l,file='MODEN') C read initial parameters read(l,*)nks,dt,nt,ih,wp C requested mode of electrostatic C energy C write(6,*)"position of mode in array ?" read(5,*)ip C time range for slope C determination C write(6,*)"Start and end times for slope evaluation" write(6,*)"in terms of (1/Wp)" read(5,*)ti,tf it=nt/ih do 10 i=l,it tt=i*dt*ih*wp read(l,*)(h(j,i)j=l,nks) if(tt.gt.tf) go to 35 if(tt.lt.ti) then jj=jj+l go to 10 endif C C ***Need to take the natural log in this case C kk=kk+l x(i-jj)=tt y(i-jj)=log(h(ip,i)) 10 continue 35 continue C least squares do 50 j = l,kk b=b+x(j)*y(j) c=c+x(j) d=d+y(j) xsq=xsq+x(j)**2 50 continue top=kk*b-c*d bot=kk*xsq-c**2 a=top/bot write(6,70)a 70 format("The approximate slope is = ",f20.15) a=a/2 write(6,80)a 80 format("The approximate imaginary frequency is = " ,f20.15) C Appendix C: User's Guide 229 P L O T M O D E C C Output electrostatic energy modes versus time C LOG, LN or regular.. C C real*8 h(8,1400a),dt,tt,lnh,wp integer ih,nks,it,nt,ip,iform open(l,nle=,MOE5EW) open(2,file='dafca') C read initial parameters re ad (1, *) n ks, dfc,pit ,jh, wp write(6,*)nks,dfe,rct,fh,wp C requested mode to view C write(6,*)"positiaEi' of mode in array ?" read(5,*)ip it=nt/ih write(6*)"LOGC% » ( 2 ) or plain graph(3) ?" read(5,*)iform do 10 i=l,it tt=ii*ih;*dt*wp readCl*Mh(j,i),j=l,nks) C M;=fog|hXip,i)) iif(iifanTKr.eq;.l) write(2,20)tt,lnh*(0.434294481903d0) if(ifer.xn..eq.2) write(2,20)tt,lnh if(faiR*q,3) write(2,20)tt,h(ip,i) 10 continue 20 forniat(2fL5'.Wj C end P L O T E K C C Outprat electesafcitk energy over all modes C for a certain iSrae; interval C C real*8 dc(rr),xflJj,wl,dt integer j,nt,ng,ifii opea("5 JBife= 'PLOTIlKJMP',form= 'unformatted') openi(2,fiiIe='data') C write(6,*)''Whaeh!tirme:Step ?" read(l)nt,ng,i:ft. write(6,*)"betwW'„l,"and",nt/ift read(5,,*)i Appendix C: User's Guide 230 ng=ng/2+l C do 10 i=l,ng 10 x(i)=(dfloat(i)-l) C do 20 i= l j 20 read(l)(ek(jj),jj=l,ng) G do 30 i=l,ng write(2,*)float(x(i)),float(ek(i)) 30 continue C end P L O T E K A C-c c c-C C C c 10 C 25 20 C 30 C Output time average of electrostatic energy modes real*8 ek(129),x(129),ekk(129) integer j,nt,ng open(l,file=,/usr/peter/data/PLOTDUMP',form=,unformatted') open(2,file=7usr/peter/data/data') read(l)nt,ng,ift ng=ng/2+l nt=nt/ift do 10 i=l,ng x(i)=dfloat(i)-l do 20 i=l,nt read(l)(ekk(jj) jj=l,ng) do 25 j=l,ng ek(j)=ek(j)+ekk(j) continue continue do 30 i=l,ng ek(i)=ek(i)/dfloat(nt) write(2,*)float(x(i)),float(ek(i)) continue end read initial parameters find average of each mode Appendix C: User's Guide 231 P L O T E C C Output total electrostatic or transverse electric field energy C C real*8 tim(9600),eyle(9600),eyre(9600) real*8 eye(9600),exe(9600),f(9600),g(9600) real*8 l(9600),m(9600) real*8 h(8,9600),wp,dt integer nkays,nth,nt,en open(4,file='PLOTHIST,)form='unformatted') open(l,file='data') C read initial parameters read (4)nkays,nth,nt,wp,dt write(6,*)nkays,nth,nt,wp,dt write(6,*)"ESE (1) or EYE (2) ?" read(5,*)imat do 10 i=l,nt read(4)tim(i),eyle(i),eyre(i),eye(i),exe(i),f(i),g(i),l(i),m(i),(h(jJi)j=l)nkays) if(imat.eq.l) write(l,*)float(dfloat(i)*wp*dt),float(exe(i)) if(imat.eq.2) write(l,*)float(dftoat(i)*wp*dt),float(eye(i)) 10 continue C end FV c C Find the particle velocity distribution function C at a certain time interval C C real*8 vx(131072),vmin,vmax,bin(256),vxbin(256) real*8 x(131072),xxx,xi,xf integer nbins.nl,n2,n,m,ng,nt,ixvx,m,nn open(l,file=7 u s r/P e t e r/ (l a t ; a/data') open(10,file=,XVXl,,form=,unformatted') open(ll)file=,XVX2',form=,unformatted,) C read initial parameters read(10)ng,nt,ixvx,nl,n2 write(6,*)ng,nt,ixvx,nl,n2 write(6,*)"Skip how many data dumps?" write(6,*)"between" ,0,"and" ,nt/ixvx read(5,*)m C number of bins in C distribution calculation write(6,*)"Number of bins?" read(5,*)nbins write(6,*)"fe only (1)" write(6,*)"fi only (2)" Appendix C: User's Guide 232 10 20 25 C C C 30 35 40 50 100 C write(6,*)"fi + fe (3)" read(5,*)inumb n=nl+l nn=nl+n2 do 25 j=0,m j=nl if(inumb.eq.l.or.inumb.eq.3) then do 10 i=l,nl read(10)x(i),vx(i) continue endif if(inumb.gt.l) then do 20 i=n,nn read(ll)x(i),vx(i) continue endif continue call minmax(nn,vmin,vmax,vx) vvrite(6,*)vmin,vmax read in (x,v) find minimum, maximum v find number of particles in velocity bin do 30 j = l,nbins bin(j)=0.D0 vxbin(j)=vmin+(vmax-vmin)*(dfloat(j)-0.5D0)/dfioat(nbins) continue do 40 i=l,nn xxx=nbins*(vx(i)-vmin)/(vmax-vmin) jj=xxx+l if(jj.gt.nbins.or.jj.lt.l) go to 40 bin(jj)=bin(jj)+l continue do 50 i=l,nbins write(l,100)vxbin(i),bin(i) continue format(fl0.8,3x,f20.1) end subroutine minmax(n,vmin,vmax,vx) real*8 vx(131072),vmin,vmax C integer n vmin=vx(l) vmax=vx(l) do 10 i=2,n vmin=dminl(vmin,vx(i)) vmax=dmaxl(vmax,vx(i)) 10 continue return c end Appendix C: User's Guide 233 F X C . C Find particle x-distribution C . C real*8 x(131072),denx(257) real*8 vx(131072) integer nl,n2,n,m,ng,nt,ixvx open(l,file='/usr/peter/dat a/data') open(10,file=,XVXl',form='unformatted') open(ll,file=:,XVX2,,form=,unformatted,) C read(10)ng,nt,ixvx,nl)n2 write(6,*)ng,nt)ixvx,nl,n2 write(6,*)"Skip how many data dumps?" write(6,*)"between" ,0,"and" ,nt/ixvx read(5,*)m write(6,*)"Number of bins?" read(5,*)nbins write(6,*)"xe (1)" ' write(6,*)"xi (2)" read(5,*)inumb n=nl+l nn=nl+n2 do 25 j=0,m j=nl C read in (x,v) if(inumb.eq.l) then do 10 i=l,nl read(10)x(i),vx(i) 10 continue endif if(inumb.eq.2) then do 20 i=n,nn read(ll)x(i),vx(i) 20 continue endif 25 continue C find number of particles C within spatial bin do 50 i=l,nn do 60 j = l,ng if(x(i).gt.j.and.x(i).lt.(j+l)) denx(j)=denx(j)+l if(x(i).lt.0) denx(ng+l)=denx(ng+l)+l 60 continue Appendix C: User's Guide 234 50 continue C ng=ng+l write(l,*)0.0,float(denx(ng+l)) do 100 i=l,ng 100 write(l *)(float(i)-0.5),float(denx(i)) write(l,*)float(ng),fioat(denx(ng+l)) C end F I E L D C-C c c-c c c c 10 c c c 20 30 C Output charge density RHO, smoothed charge density RHOS, potential PHI or electrostatic field EX real*4 xj(256),f(256) integer a,n,ng,nt,idump character*10 fil open(l,file='data') field to plot write(6,*)"File to read data from (RHO, RHOS, PHI, EX) ?" read(5,*)fil open(ll,file=fil) read(ll)ng,nt,idump n=nt/idump do 10 j=l,ng xj(j)=j-l write(6,*)"IntervaL to plot ?",l,"to",n read(5,*)a if(a.lt.0) stop do 20 i=l,a read(ll)f . do 30 i=l,ng write(l,*)float(xj(i)),float(f(i)) end read initial parameters set up grid points plot field at requested interval T I M E F C C Output right and left going field Appendix C: User's Guide 235 C energy for a mode over time C C real*8 time(10000),hl(8)10000),hr(8,10000),wp real*4 t,h integer nt,nkays,kay(8) open(8,file=7usr/peter/data/TEMPDUMP'>form='unformatted') open(9,file='/usr/peter/data/data') C read(8)nt,nkays,wp write(6,*)"NT\nt,"NKAYS" ,nkays read(8)(kay(j),j=l,nkays) write(6,*)"KAY",(kay(j)j=l,nkays) C C input mode energies do 10 i=l,nt read(8) time(i),(hl(k,i),k=l,nkays),(hr(k,i),k=l,nkays) time(i)=time(i)*wp 10 continue C write(6,*)" Which component of KAY do you wish to look at ?" write(6,*)"KAY:",(kay(j),j=l,nkays) write(6,*)"(choose between",1,"and",nkays,")" read(5,*)im write(6,*)"Right going (1)" write(6,*)"or left going (2)" write(6,*)"fields?" read(5,*)id C if(id.gt.l) go to 100 C do 20 i=l,nt t=i h=hr(im,i) 20 write(9,*)t,h C if(id.eq.l) stop C 100 continue do 110 i=l,nt t=time(i) h=hl(im,i) 110 write(9,*)t,h C end C C MODEY Output mode distribution of right and left moving fields Appendix C: User's Guide 236 C c real*8 a(10000),b(10000),x(256),dt,wp integer ij^ij.npts.nf.nt.ift.ng open(10,file=7usr/peter/data/LRMODES',form=,unformatted') open(2,file='/usr/peter/data/data') read(10)ift,ng,nt,dt,wp npts=ng/2+l ij=int(nt/ift) write(6,*)ij,npts,nt 5 write(6,*)"Time Interval ? " read(5,*)i if'(i.gt.ij) tlien write(6,*)"Interval must be less than",ij go to 5 endif do 7 j=l,npts 7 x(j)=dfloat(j-l) do 10 j = l,i read(10)(a(k),k=l,npts) read(10)(b(k),k=l,npts) 10 continue do20j = l,npts write(2,100)x(j),a(j) 20 continue do 30 j=npts,l,-l write(2,100)x(j),b(j) 30 continue 100 format (2f 15.10) C MODEYAV C C Output average mode distibution of C right and left moving fields C C real*8 a(100000),b(100000),x(300),dt>wp integer ij.k.ij.npts.nt^mp.ift.ng open(10,file='LRMODES,,form='unformatted') open(2,file='data',form='formatted') C input initial vlaues read(10)ift,ng,nt,dt,wp write(6,*)ift,ng,nt,dt,wp ij=int(nt/ift) npts=ng/2+l write(6,*)ij,npts,nt C set up grid points Appendix C: User's Guide 237 do 5 i=l,npts 5 x(i)=dfloat(i-l) C k=l i=l do 10 i=l,ij read(10)(a(j)j=k,npts+k-l) C read(10)(b(j),j=k,npts+k-l) k=k+npts 10 continue k=l i=0 j=0 tmp=ij*npts do 20 k=l,npts do 30 j=npts,tmp,npts a(k)=a(k)+a(k+j) b(k)=b(k)+b(k+j) 30 continue 20 continue write(6,*)"Output leftgoing(l);" write(6,*)"Output rightgoing (2);" write(6,*)"Output both (3);" read(5,*)ii if(ii.lt.2)go to 55 do 100 k=l,npts 100 write(2,130)x(k),(a(k)/dfloat(tmp/npts)) if(ii.eq.2) stop 55 do 120 k=npts,l,-l 120 write(2,130)x(k),(b(k)/dfloat(tmp/npts)) 130 format(2fl5.10) C end
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- A Vlasov plasma simulation code
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
A Vlasov plasma simulation code Loewenhardt, Peter Karl 1989
pdf
Page Metadata
Item Metadata
Title | A Vlasov plasma simulation code |
Creator |
Loewenhardt, Peter Karl |
Publisher | University of British Columbia |
Date Issued | 1989 |
Description | A 1-1/2 dimensional, electromagnetic Vlasov plasma simulation code, relativistically correct, was constructed and tested. The code can deal with one or two species, each with a Maxwellian distribution and a possible drift velocity. The code also allows external fields, such as a laser, to be included in the simulations. Simulations of Landau damping, the two-stream instability and stimulated Raman scattering were carried out and compared with theory and with the results of the particle code EM1. When dealing with electrostatic problems, the Vlasov code gave results agreeing very closely with theory in both non-relativistic and relativistic regimes. Here the Vlasov code preformed better than EM1 which had problems with background noise and longer run times. When a laser field was included within the simulations, however, the Vlasov code produced some spurious features, unlike EM1. These anomalous features may be caused by aliasing, recurrence or by some other unknown effect. Since realistic results are produced as well, it is believed that this problem can be overcome. When a very high intensity laser was included in the simulations the Vlasov code produced much better results but was plagued by very long run times. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-08-21 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0085071 |
URI | http://hdl.handle.net/2429/27586 |
Degree |
Master of Science - MSc |
Program |
Physics |
Affiliation |
Science, Faculty of Physics and Astronomy, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1989_A6_7 L63.pdf [ 9.91MB ]
- Metadata
- JSON: 831-1.0085071.json
- JSON-LD: 831-1.0085071-ld.json
- RDF/XML (Pretty): 831-1.0085071-rdf.xml
- RDF/JSON: 831-1.0085071-rdf.json
- Turtle: 831-1.0085071-turtle.txt
- N-Triples: 831-1.0085071-rdf-ntriples.txt
- Original Record: 831-1.0085071-source.json
- Full Text
- 831-1.0085071-fulltext.txt
- Citation
- 831-1.0085071.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0085071/manifest