Combining Hydrodynamic Modeling with Nonthermal Test Particle Tracking to Improve Flare Simulations by Henry deGraffenried Winter III A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy in Physics MONTANA STATE UNIVERSITY Bozeman, Montana April 2009 c©COPYRIGHT by Henry deGraffenried Winter III 2009 All Rights Reserved ii APPROVAL of a dissertation submitted by Henry deGraffenried Winter III This dissertation has been read by each member of the dissertation committee and has been found to be satisfactory regarding content, English usage, format, citations, biblio- graphic style, and consistency, and is ready for submission to the Division of Graduate Education. Dr. Petrus C.H. Martens Approved for the Department of Physics Dr. Richard J. Smith Approved for the Division of Graduate Education Dr. Carl A. Fox iii STATEMENT OF PERMISSION TO USE In presenting this dissertation in partial fulfillment of the requirements for a doctoral degree at Montana State University, I agree that the Library shall make it available to bor- rowers under rules of the Library. I further agree that copying of this dissertation is al- lowable only for scholarly purpose, consistent with “fair use” as prescribed in the U.S. Copyright Law. Requests for extensive copying or reproduction of this dissertation should be referred to Bell & Howell Information and Learning, 300 North Zeeb Road, Ann Arbor, Michigan 48106, to whom I have granted “the exclusive right to reproduce and distribute my dissertation in and from microform along with the non-exclusive right to reproduce and distribute my abstract in any format in whole or in part.” Henry deGraffenried Winter III April 2009 iv DEDICATION To my wife Denlyn, Thank you for tremendous amount of support. I know that you have sacrificed much and for many years so that I could achieve this dream and many others. Yet, through it all, you stayed with me and we worked together to build a life that makes me proud. That is amazing to me, but I’ve come to expect amazing from you. To my mother, Dr. Betty Winter, You showed me that someone could change their life to become the person they wanted to be, at any point in their life. You are a shining example of someone who can switch career paths with grace and dignity. I am sorry I that I did not emulate the grace and dignity part, but that is part of me being the person I want to be. I could not have done this without you. To you, I dedicate this dissertation. v ACKNOWLEDGMENTS This work would not have been possible without the support of my advisor, Dr. Petrus Martens. Patience does have some rewards. The author would also like to thank Dr. David McKenzie. Dr. McKenzie was in- strumental in providing the initial computing resources that made the testing of the code possible. Dr. McKenzie also provided time and support that aided the author in the comple- tion of this project. Jeremy Gay and Dr. Ladean McKittrick generously provided additional computing resources that made the completion of this work possible. A special thanks to Margaret Jarrett. Any student that matriculates through the Physics Graduate Program at Montana State University realizes what a special position she holds. Every physics student has been the recipient of her aid at one time or another, whether they realize it or not. Thank you. Thanks to Dr. Angela C. Des Jardins. Good friends are hard to find, good office mates are even harder. This research made extensive use of NASA’s Astrophysics Data System and was sup- ported by NASA grant NAG5-12820. vi TABLE OF CONTENTS 1. INTRODUCTION TO SOLAR FLARES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 General Properties. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .1 Emission . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Energy Release . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 Looptop X-ray Sources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 Research Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2. THEORY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 Equations Describing the Thermal Distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 Thermal Distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 MHD Equations . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 Grids. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 Spatial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 Temporal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 Radiative Losses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .32 Corona. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 Chromosphere and Transition Region . . . . . . . . . . . . . . . . . . . . . . . . . . 34 Radiative Loss Time Scale. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 Notes on Radiative Losses. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 Conduction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .37 Saturated Flux Limitation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .38 Conduction Time Scale. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 Tests. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 Nonthermal Particle Transport. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 Stochastic Simulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 Coulomb Scattering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 Close Collisions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 Far Collisions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 Collision Time Scales. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 Effects of a Nonuniform Magnetic Field. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 PaTC Time Scales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 Tests. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 Emission Calculations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 3. STATIC VERSUS DYNAMIC ATMOSPHERES . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 Experimental Setup. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 Loop Geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 vii TABLE OF CONTENTS-CONTINUED Nonthermal Particles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 Beam Species. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 Beam Time and Energy Distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 Pitch-Angle Cosine Distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 Error Tests. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 Grid Spacing Test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 Monte Carlo Tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 The Effect of a Hydrodynamic Plasma on a Nonthermal Particle Beam . . . . . . . . . . . 81 Effects on HXR emission. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 Apparent HXR Source Motions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 4. THE EFFECTS OF DENSITY GRADIENTS ON NONTHERMAL PARTICLES 109 Experimental Setup. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 Results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 Discussion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 5. THE EFFECTS OF DIFFERENT PITCH ANGLE DISTRIBUTIONS . . . . . . . . . . 118 Experimental Setup. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 Thermal Evolution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 SXR Signatures. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 HXR Signatures. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 6. CONCLUSION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137 Future Studies. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139 Alfvén Wave Acceleration of Particles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139 Code Improvements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141 Improved Test Particle Sampling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142 CUDA Programming. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144 APPENDICES. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146 APPENDIX A: List of Symbols. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146 APPENDIX B: Numerical Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149 REFERENCES. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160
 viii LIST OF TABLES Table Page 1.1 Mount Wilson sunspot group classification system. Sunspot groups can also be described as a mixture of the above classifications. A complete description is available at http://www.spaceweather.com/glossary/magneticclasses.html. . . . . . . . . 3 1.2 Solar flare classifications based upon the 1− 8 Å channel of the GOES satellites. Each classification is subdivided into 10 subsections (1-10), each with increasing energy. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 3.1 Initial properties of the post flare loop. . . . . . . . . . . . . . . . . . . . . 71 4.1 A table showing the properties of the constant cross-section tube. . . . . . . 110 ix LIST OF FIGURES Figure Page 1.1 Sample X-ray flare spectrum as observed by the RHESSI satellite. The SXR emission is fitted with a thermal component (green) having a temperature of 16.7 MK. The HXR emission is fitted by power laws, which appear as straight lines on a log-log plot. The power law has breaks, or “knees”, at 12 keV and 50 keV. The origin of these breaks is still a matter of scientific debate. The black dashes are observed data points and associated errors (Benz, 2008). . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.2 A diagram of a flare from Martens and Kuin (1989). Note the lack of a concentrated looptop emission source. The frequency at which such sources occur was a surprise to theorists and remains largely unexplained. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 1.3 The famous Masuda flare (Masuda et al., 1994) of 1992 Jan. 13. Hard X-rays are shown in contours and the results from soft X-rays are in color. From left to right, the first two columns show images from the SXT thick aluminum and beryllium filters respectively and the temperature and emission measure derived from their ratio in the last two columns. The rows represent images obtained form HXT in three hard X- ray energy bands: 13.9-22.7 keV, 22.7-32.7 keV, and 32.7-52.7 keV respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.1 A visual representation of the staggered grid system. The energy and electron number densities, ε , ne are defined in N volumes. Other quantities, such as velocity, v, magnetic field strength, B, area, A, and parallel component of acceleration due to gravity, g, are defined on N − 1 surfaces. The first and last volume serve as static boundary conditions for the grids and are not altered by the simulation. . . . . . . . . . . 29 2.2 The radiative loss function as used by the MSULOOP code for the current set of experiments. The diamonds represent data points and the line shows interpolated values. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 2.3 A plot showing the normalized temperature profile for a uniformly heated loop. The red line shows the analytic solution based on (Martens, 2008) and the green dashed line shows the result from MSULOOP. In this simulation saturated flux effects are not taken into account. The percent differences between the numerical and are less than 2%. . . . . . . . . . . . . . . 42 x LIST OF FIGURES-CONTINUED Figure Page 2.4 A diagram showing the interaction of two particles. The impact parameter, b, is defined as the distance of closest approach between the test and field particle. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 2.5 A simple diagram showing the relationship between the collision frame and the loop frame. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 2.6 A plot showing the percent differences in PATC’s calculation of the change in a particle’s energy in a 0.0001 second time step and an analytic approximation to the energy change. The differences are interpreted as being caused by the numerical treatment’s increased accuracy, since it uses fewer approximations than the analytic expression. . . . . . . . . . . . . . . 58 2.7 A plot showing the ratio of the analytic change in nonthermal particle en- ergy to PATC’s calculation of the energy change in a particle’s energy in a 0.0001 second time step. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 2.8 A plot showing the percent differences between an analytical calculation of the magnetic field strength where a test particle should mirror and the numerical calculation. The red, green, and blue lines represents a particle with a pitch-angle cosine of .4, .2, and one respectively. . . . . . . . . . . . . . . 61 2.9 A diagram representing the calculation of instrument responses on an arbitrarily sized pixel grid. The SHOW_ LOOP software performs a path length integration through the volume emission of a segment of the loop observed by a pixel. This combined with the area of the loop observed by the pixel provides a total signal from a loop of arbitrary inclination angle observed by pixel. Once the signal for the entire grid is computed, the grid can then by convolved with a function to take effects such as point spread functions, or defocus into account. . . . . . . . . . . . . . . . . . . . . . . . . . 64 2.10 An image built up using multiple, simulated loops. The image on the far left shows a system of loop strands as would be viewed by the XRT instrument on Hinode in the Al-Poly channel with a 2′′ resolution. The center image shows the same system, in the same bandpass but with a 0.5′′ resolution. The last image shows the system in a theoretical bandpass for the AIA mission but with an increased resolution of 0.1′′ . . . 65 xi LIST OF FIGURES-CONTINUED Figure Page 3.1 A contour map of the magnetic field generated by a Green (1965) style current sheet using the formalisms of Bungey and Priest (1995). Notice the two null points which are regions of zero magnetic field. The characteristic field line for the flare loop geometry is shown at bottom. . . 69 3.2 The field lines calculated by the Green current sheet model. The blue line represents the current sheet in the half plane. The red line illustrates the field line chosen for the basis of the simulated loop’s geometry. From this figure it is easy to see why this is called Y-type. The current sheet connects with the separatrix field lines, just above the red line, to form an inverted Y. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 3.3 A log-log plot of the number of particles as a function of energy. The slope of this corresponds to the electron spectral index δ . A value of δ = 3 was input into the beam generator. The fitted value is shown as δ = 2.8. 75 3.4 A plot of various pitch-angle cosine distribution functions, f (µPA), corresponding to γPA = −5, −3, 0, 3, 5. A pitch-angle cosine of 1 (−1) corresponds to the nonthermal particle’s momentum being directed parallel (anti-parallel) to the loop axis. A pitch-angle cosine of zero corresponds to the nonthermal particle’s momentum being directed perpendicular to the loop axis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 3.5 A plot of the pitch-angle cosine distribution to be used in the current experiment with γPA = −4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 3.6 A graph showing the percent difference between the test variables, temperature and density, at a given number of grid cells in comparison to a simulation run with 1400 grid cells. This still was taken at the beginning of the simulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 3.7 Same as 3.6 but taken after 5 second s of simulated time. . . . . . . . . . . . . . 80 3.8 Same as figures 3.17 and 3.7 but at 300 second s. By analyzing plots such as this for every time step it becomes apparent that increasing the number of grid cells past 700 only yields minor improvements in accuracy. . . . . . 80 xii LIST OF FIGURES-CONTINUED Figure Page 3.9 This plot shows how the average of the apex density and temperature and the associated standard deviation changes as a function of the number of runs five second s into the simulation, . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 3.10 The same as 3.9 but 10 second s into the simulation. . . . . . . . . . . . . . . . . . 82 3.11 The same as 3.9 and 3.10 but 300 second s into the simulation. . . . . . . . . 83 3.12 A plot showing how the percent standard deviation for temperature and density changes as a function of time for five runs. . . . . . . . . . . . . . . . . . . 83 3.13 A plot showing how the percent standard deviation for temperature and density changes as a function of time for eight runs. . . . . . . . . . . . . . 84 3.14 A plot showing how the percent standard deviation for temperature and density changes as a function of time for 12 runs. . . . . . . . . . . . . . . . . . . 84 3.15 A plot showing how the relative standard deviation for temperature changes as a function of time. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 3.16 A plot showing how the percent standard deviation for 3 − 6 keV and 12 − 25 kev emission as a function of time. . . . . . . . . . . . . . . . . . . . . . . . . . . 85 3.17 Hard X-ray emission in the 3 − 6 keV passband for a dynamic atmosphere simulation. The top row shows emission simulated in an instrument with 7′′ pixels and a Gaussian point spread function with a FWHM of 3 pixels. The first column shows the total emission, the second shows the nonthermal component only, and the last column shows the thermal component only. The lower row shows the light curve of the total (black), nonthermal (green) and thermal (red) emission in the box defining the loop apex. Images are binned in 4 second increments. Contours enclose the 40%, 60%, and 80% levels. . . . . . . . . . . . . . . . . . . . 88 3.18 Same as figure 3.17 but from 4 − 8 second time bin. . . . . . . . . . . . . . . . . . 89 3.19 Same as figure 3.17 but from 296 − 300 second time bin. . . . . . . . . . . . . . 90 xiii LIST OF FIGURES-CONTINUED Figure Page 
 3.20 Hard X-ray emission in the 3 − 6 keV passband for a static atmosphere simulation. This plot shows emission in the 0 − 4 second time bin. . . . . . 91 3.21 Hard X-ray emission in the 3 − 6 keV passband for a static atmosphere simulation in the 4 − 8 second time bin. . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 3.22 Hard X-ray emission in the 3 − 6 keV passband for a static atmosphere simulation in the 0 − 4 second time bin. . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 3.23 A comparison of the Monte Carlo estimator value of the apex hard X-ray emission in the 3 − 6 keV passband in a dynamic (red) and static (blue) atmosphere. Error bars show the +/ − sN-1 value. . . . . . . . . . . . . . . . . . . . . 94 3.24 Same as figure 3.23 but plotted on a log scale. Error bars have been removed for clarity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 3.25 Distributions of sampled data with the dynamic case in red and the static case in blue. The means of the distribution, which comprise the value for the Monte Carlo estimator, are over plotted as a dashed line. Flux units are arbitrary. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 3.26 Kolmogorov-Smirnov test of the differences hard X-ray emission in the 3 − 6 keV bandpass in the dynamic and static atmosphere simulations. The red line maps the evolution of the D statistic as a function of time. The blue lines shows the confidence level with which the null hypothesis can be rejected. This provides strong, quantitative evidence that the evolution of the nonthermal beam is significantly statistically different in the dynamic atmosphere than in the static atmosphere. . . . . . . . . . . . . . 97 3.27 Plots of the pressure, density, average velocity, and temperature as a function of loop coordinate. This plot is at the beginning of the simulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 3.28 Same as figure 3.27 but after four seconds of simulated time. . . . . . . . . . . 100 3.29 State variables plotted after 125 seconds of simulated time. . . . . . . . . . . . 100 3.30 Plots of the state variables at 300 seconds. . . . . . . . . . . . . . . . . . . . . . . . . . 101 xiv LIST OF FIGURES-CONTINUED Figure Page 
 3.31 A plot of the density at four seconds of simulated time. Red represents the dynamic atmosphere and blue represents the static atmosphere, which is the same as the dynamic atmosphere at t = 0 seconds. . . . . . . . . . . . . . . 102 3.32 A plot of the density at 15 seconds of simulated time. . . . . . . . . . . . . . . . . 103 3.33 A plot comparing the static and dynamic atmosphere density at 125 seconds. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 3.34 Density plot at 300 seconds. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 3.35 The initial pitch-angle cosine distribution for the dynamic atmosphere simulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 3.36 The pitch-angle cosine distribution for the dynamic atmosphere after four seconds of simulated time. Note that at this time the beam pitch-angle cosine distribution is closer to a Gaussian centered around zero. Particles with large absolute values of pitch-angle cosine have already lost all of their energy in the denser loop legs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 3.37 The pitch-angle distribution for the dynamic atmosphere after 30 seconds of simulated time. Note that the distribution of pitch-angle cosines is steeply peaked around zero. Particles that are at the loop apex, with such a small pitch-angle cosine will scatter less due to the low density of the loop apex. With less scattering, and a pitch-angles near zero, any particle near the apex will now be effectively trapped there. . . . . . . . . . . . . . . . . . 106 4.1 A plot showing the density profiles of the tubes used in this experiment. . 110 4.2 A plot showing the initial pitch-angle cosine distributions of the nonthermal particle beams injected into the center of the tube. The quantity IPA is defined as the number of particles with a pitch-angle cosine between the two blue dotted lines, in the central portion of the distribution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 4.3 A plot showing the median values of the inner pitch-angle ratio, IPA/N . Error bars were taken using a bootstrap with replacement method and represent the 95% confidence level. A completely uniform distribution would have a {I pos/N } med = 0.5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 xv LIST OF FIGURES-CONTINUED Figure Page 
 4.4 This plot shows the median values inner position ratio, I pos /N, for each value of the density gradient factor as a function of time. . . . . . . . . . . . . . 113 4.5 This plot is the same as figure 4.4 but zoomed in at 2 ≤ t ≤ 6 . . . . . . . . . . 114 4.6 A plot showing the nonthermal beam alive times as a function of gradient factor averaged over the ten runs. . . . . . . . . . . . . . . . . . . . . . . 115 4.7 A plot showing the skew of the distribution of alive times for all runs. . . . 116 5.1 A plot showing pitch-angle cosine distributions for the three simulations to be conducted. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 5.2 This plot shows the mean apex temperature, in kelvins, of a flaring loop for each simulation. The γPA = 4 case became almost three times as hot as the γPA = −4 case. The error bars denote one standard deviation. . . . . . 120 5.3 This plot shows the mean radiative loss rates, as a function of time, for each of the simulated cases. The error bars denote one standard deviation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 5.4 This plot shows the mean apex thermal electron density for each simula- tion. The maximum density enhancement of the γPA = 4 came nearly a minute before the γPA = −4 case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 5.5 This plot shows the emission in XRT filter channels for the γPA = −4 pitch-angle cosine distribution case. The plot at the bottom of the figure the apex to footpoint emission ratio for each case as a function of time. . . 124 5.6 This plot shows the emission in XRT filter channels for the γPA = 0 pitch- angle cosine distribution case. The plot at the bottom of the figure the apex to footpoint emission ratio for each case as a function of time. . . . . . 125 5.7 This plot shows the emission in XRT filter channels for the γPA = 4 pitch- angle cosine distribution case. The plot at the bottom of the figure the apex to footpoint emission ratio for each case as a function of time. NB This was the only case that showed an initial peak in apex to footpoint emission ratio. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 xvi LIST OF FIGURES-CONTINUED Figure Page 
 5.8 This plot shows the apex to footpoint emission ratio for each filter channels. While the ratios differed slightly in magnitude from filter to filter, the obvious signal was the time profile of the ratio. . . . . . . . . . . . . . 127 5.9 This plot shows the HXR emission for the γPA = 0 case in the t = 0 − 4 time bin. The bottom plot shows the total apex emission as defined by the bounding box as a function of time with the nonthermal signal in green, the thermal signal in red, and the total signal in black. . . . . . . . . . . . . . . . . 128 5.10 This plot shows the HXR emission for the γPA = 0 case in the t = 4 − 8 time bin. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 5.11 This plot shows the HXR emission for the γPA = 0 case in the t = 52 - 56 time bin. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 5.12 This plot shows the HXR emission for the γPA = 0 case in the t = 296 − 300 time bin. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131 5.13 This plot shows the HXR emission for the γPA = 4 case in the t = 0 − 4 time bin. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 5.14 This plot shows the HXR emission for the γPA = 4 case in the t = 4 − 8 time bin. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 5.15 This plot shows the HXR emission for the γPA = 4 case in the t = 40 − 44 time bin. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 5.16 This plot shows the HXR emission for the γPA = 4 case in the t = 296 − 300 time bin. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 6.1 This plot illustrates how the use of Sobol numbers can increase the fractional accuracy the Monte Carlo calculation of pi while using fewer pairs of random numbers. The upper dashed line represents N−1/2 while the lower dashed line represents N −1 . The blue line charts the fractional accuracy as a function of number of uniformly distributed random pairs, which matches the N−1/2 line well. The black solid line is the fractional accuracy as a function of number of Sobol pairs. While the Sobol pairs do not exactly track the N −1 trend they do show an orders of magnitude improvement over uniform pairs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143 xvii CONVENTIONS In this work modified Gaussian units (mgu) are used. This system of units may not be familiar to scientists in other fields who commonly use the International System of Units (SI). The mgu unit system uses centimeters, grams and seconds for non-electromagnetic quantities, electrostatic units (esu) for electrical quantities and electromagnetic units for magnetic quantities (Sturrock, 1994). Despite the adoption of SI in most fields of science, mgu remains the most widely used system in plasma physics (Benz, 1993b). The mgu unit of energy is the erg. In all formulae dealing with the rate of change of the energy, this unit is assumed. However, due to historical reasons, the unit of energy most often used by high energy astronomers is the electron-volt (eV ) and particle kinetic energies and photon energies are expressed in kilo-electron-volts (keV ) or mega-electron-volts (meV ). Temperatures are always expressed in kelvins in this work, though many other plasma physics and kinetics texts express temperatures in terms of eV as well. There is much confusion regarding the precise definition of the kurtosis of a distribution in current literature. Properties of distributions are normally described by taking successive moments. In older and some current texts, the kurtosis is defined as the fourth moment of the distribution. However, some newer texts now label the kurtosis as being the fourth mo- ment minus three, the kurtosis for a normal distribution. This quantity used to be known as the kurtosis excess. However, this work uses fourth moment minus three for the definition of the kurtosis, since this seems to be the modern, if not more confusing, standard. This work deals heavily in the construction and use of software routines to simulate physical process. Individual routine names are italicized, E.g. get_loop_mfp.pro, and soft- ware suites containing many subroutines are written in an all capital notation, E.g. PATC. xviii ABSTRACT Solar flares remain a subject of intense study in the solar physics community. These huge releases of energy on the Sun have direct consequences for humans on Earth and in space. The processes that impart tremendous amounts of energy are not well understood. In order to test theoretical models of flare formation and evolution, state of the art, numerical codes must be created that can accurately simulate the wide range of electromagnetic radi- ation emitted by flares. A direct comparison of simulated radiation to increasingly detailed observations will allow scientists to test the validity of theoretical models. To accomplish this task, numerical codes were developed that can simulate both the thermal and nonther- mal components of a flaring plasma, their interactions, and their emissions. The HYLOOP code combines a hydrodynamic equation solver with a nonthermal particle tracking code in order to simulate the thermal and nonthermal aspects of a flare. A solar flare was simulated using this new code with a static atmosphere and with a dynamic atmosphere, to illustrate the importance of considering hydrodynamic effects on nonthermal beam evolution. The importance of density gradients in the evolution of nonthermal electron beams was investi- gated by studying their effects in isolation. The importance of the initial pitch-angle cosine distribution to flare dynamics was investigated. Emission in XRT filters were calculated and analyzed to see if there were soft X-ray signatures that could give clues to the nonther- mal particle distributions. Finally the HXR source motions that appeared in the simulations were compared to real observations of this phenomena. 1 CHAPTER 1 INTRODUCTION TO SOLAR FLARES The study of the Sun has a long and rich history, from ancient notions that it was a hot ball of iron, to utilizing it as a research laboratory in order to solve basic questions of particle physics. The interested reader is directed to Phillips (1995) and Zirker (2002) for in-depth and fascinating accounts of solar studies. Despite centuries of study and over 40 years of intense scrutiny by satellite observatories, there is still much about the Sun that is not understood. In particular, solar flares, the most powerful explosions in our solar system, remain a mysterious phenomenon. A basic overview of solar flares is provided to aid in the understanding the work as a whole. General Properties On September 1, 1859, Richard Carrington and R. Hodgson observed an intense bright- ening in the visible light emission of the Sun. This was the first recorded observation of a solar flare (Phillips, 1995). This was a serendipitous observation since such white light flares are actually quite rare with only 60 being recorded from Carrington’s initial obser- vation to 1992 (Phillips, 1995). Even though the Sun emits primarily in the visible wave- lengths, flares typically do not. Only a fraction of the emitted energy of a flare, about 25% or so, is emitted in the visible spectrum and usually is insignificant when compared to the photospheric background (Phillips, 1995). Little did Carrington and Hodgson know that they had witnessed one of the largest releases of energy in the solar system. Later measurements showed that flares could release energy up to 1033 ergs in time scales of tens of seconds to minutes. While this is a small amount compared to the Sun’s 2 total luminous output, L% & 3.83×1033 ergsec−1, it still easily dwarfs that of thermonu- clear weapons1 or even that of the Shoemaker-Levi 9 comet’s impact with Jupiter (Ivlev et al., 1995)!2 It is easy to see why the understanding of solar flares is not only of impor- tance to solar theorists but has practical applications as well. The field of Space Weather studies how flares and other solar phenomena affect the Earth’s local environment. Some of these effects can influence life on Earth directly. Flares are often associated with massive clouds of plasma, with their own magnetic fields, called Coronal Mass Ejections (CMEs). Massive rearrangements of the Earth’s magnetic field can occur if the magnetic field of a CME is aligned opposite to the Earth’s magnetic field (Zirker, 2002). This rearrangement generates electrical currents in sensitive electronics or in long wires and pipelines. In 1989 a magnetic disturbance caused a power blackout in Quebec that affected millions (Oden- wald, 1999). Communication blackouts are also common during these “magnetic storms”. In 1984 Air Force One lost communications with the United States as President Reagan flew to China (Hilts, 1988). Understanding flares is even more important now that low Earth orbit is “permanently” inhabited. During solar flares astronauts can be exposed to high levels of ionizing radiation that can cause long-term illnesses, such as an increased risk for cancer, or even death. A solar storm in 1990 gave Mir cosmonauts a year’s radi- ation dosage in just a few days (Odenwald, 1999). Due to the relativistic speed of these particles there is very little warning time before the damage is done. Flares can occur over many different spatial scales including, theoretically, those too small to detect with current instrumentation. A very large flare can cover an area up to 3×109 square kilometers or 1% of the visible solar surface (Phillips, 1995). The location of solar flares, however, is not uniformly distributed across the solar disk. Ninety percent 1A Megatonne equivalent of TNT has an energy output of∼ 4.2×1023 ergs (Thompson and Taylor, 2008). The largest fusion warhead developed by the United States of America, the B-41, had an estimated yield of 25 Megatonnes (Sublette, 2007). If true, then it would take 107 weapons of this type to release the same energy of an X class solar flare. 2The largest fragment of the Shoemaker-Levy 9 comet was estimated to have a kinetic energy of 1030 erg. 3 Table 1.1: Mount Wilson sunspot group classification system. Sunspot groups can also be described as a mixture of the above classifications. A complete description is available at http://www.spaceweather.com/glossary/magneticclasses.html. Type Description α Unipolar sunspot group β A sunspot group having distinct positive and negative magnetic polarities γ A complex active region with irregular grouping of positive and negative polarities δ A qualifier to the other sunspot groups that indicates that the umbrae of sunspots of different polarities are separated by less than 2 degrees within one penumbra have opposite polarity of flares occur in active regions (Strong et al., 1999). In the corona, the hot, tenuous, outer solar atmosphere, active regions are areas that contain large loop structures emitting in the extreme ultraviolet and X-ray portions of the electromagnetic spectrum. In the solar photosphere, the visible surface of the Sun, these active regions are associated with dark sunspot groups (Zirker, 2002). These sunspot regions are areas of intense magnetic fields with strengths of about 4× 103 Gauss on average (Phillips, 1995). Sunspot groups can have many different configurations. The configurations most associated with flares are βγ sunspot regions as defined by the Mount Wilson Classification system (Phillips, 1992). In this classification scheme, a βγ region is defined as a region of bipolar positive and negative magnetic field with irregular groupings of flux that give rise to complex magnetic field configurations (see table 1). This association of flares with regions of strong and complex magnetic fields gives strong clues about the origin of these events. The Sun has shown a surprisingly regular, 11 year cycle of magnetic activity (Phillips, 1995) which also correlates well to flare activity. At the peak of solar magnetic activity, the lower solar latitudes have high occurrences of sunspots and flares are more frequent (Zirker, 2002). 4 Emission The exact definition of hard X-rays (HXR) tends to vary slightly depending upon the bandpasses of a particular instrument but are usually defined as photons having energies of 10keV to 120keV (1.03× 10−1 ≤ λ ≤ 1.24 Å). HXR emission is strongly associated with the impulsive onset of solar flares. HXR events are usually impulsive, short-lived events with fluctuations on the order of 0.02 seconds and lifetimes on the order of minutes. Some, very rare, features have started as soft X-ray sources and increased steadily in en- ergy to become HXR features. These events can have lifetimes on the order of hours. If the emissions from the more commonly seen impulsive HXR flare events are from thermal electrons that would imply a plasma temperature on the order of 109 K (Phillips, 1995). However, high-energy, nonthermal electrons being deflected by slow protons in a free-free, Bremsstrahlung radiation mechanism can also explain the emission and many theories of flare evolution predict the emission of nonthermal electron beams (Zirker, 2002). While newer instruments have aided in determining whether HXR emission is thermal or nonther- mal, many events remain hard to interpret. The HXR spectra of flares tend to be dominated by power law distributions in photon energy (Benz, 2008) expressed as I ( εγ ) = aε−γp γ , (1.1) where I ( εγ ) is the HXR Bremsstrahlung flux ( photonscm−2 s−1 keV−1), a is a constant, εγ is the photon energy and γp is the photon spectral index. Such distributions are described as having a hard spectrum if γp is small (γp ! 4) and having a soft spectrum if γp is large (γp " 5) . Brown (1971) showed that one could infer a power law distribution of the injected electron energy spectrum from the HXR photon spectrum. The power law distribution of 5 the electron energies is expressed by F = F0E−δE , (1.2) where F is in units of electronscm−2 s−1 keV−1 and δE is the electron energy spectral in- dex. In the simplified case this relationship takes two simple forms; γp = δE +1 when the observation time is much shorter than the energy loss time of the electrons or a thin-target model, and γp = δE−1 for a thick-target model where the nonthermal particle loses enough energy to become thermalized (Brown, 1971; Tandberg-Hanssen and Emslie, 1988). The observations of a power law distribution in emitted HXR photon energy in flares strongly suggest the presence of a nonthermal distribution of electrons. However, it should be noted that this simplistic relationship has recently come under increased scrutiny and may not hold to be true under more rigorous analysis (Brown et al., 2007). Determining the true distributions of the nonthermal particles and the mechanisms responsible for their acceler- ation is a matter of intense research. Flare regions emit strongly in soft X-rays (SXR), photons having energies between 0.1 and 10keV (1.24≤ λ ≤ 124 Å), after their initial impulsive onset. Solar soft X-rays consist of both continuum and spectral soft X-ray emission produced by thermal plasma and possi- bly a difficult to detect nonthermal component. Measurements of SXR spectral lines by the Bent Crystal Spectrometer (BCS) on board the Solar Maximum Mission (SMM) show flare plasma can have temperatures up to 2× 107 K (Strong et al., 1999). An interesting fact discovered by SXR observations is that cooler plasma with temperatures of ∼ 2×106 K can be observed in a flare region co-existing with the hotter plasma. The only way that this can be possible is if the two plasmas are contained within separate magnetic field flux tubes. This means that some flaring regions consist of several loop features at different characteristic temperatures (Phillips, 1995). 6 Figure 1.1: Sample X-ray flare spectrum as observed by the RHESSI satellite. The SXR emission is fitted with a thermal component (green) having a temperature of 16.7 MK. The HXR emission is fitted by power laws, which appear as straight lines on a log-log plot. The power law has breaks, or “knees”, at 12 keV and 50 keV. The origin of these breaks is still a matter of scientific debate. The black dashes are observed data points and associated errors (Benz, 2008). 7 Table 1.2: Solar flare classifications based upon the 1−8 Å channel of the GOES satellites. Each classification is subdivided into 10 subsections (1-10), each with increasing energy. Class Flux range measured in W m−2 Description B F ≤ 10−6 Common during times of high magnetic activity C 10−6 ≤ F ≤ 10−5 Few noticeable effects at Earth M 10−5 ≤ F ≤ 10−4 Can cause radio blackouts in the polar regions and aurora X F ≥ 10−4 Can cause world-wide blackouts, major magnetic storms and aurorae, major threat to astronauts The strength of a flare’s SXR emission provides an important measure of the flare strength and the potential impact it can have on the Earth-space environment. For this reason the current flare classification system is based on the intensity of the flare between 1− 8 Å. This system was established by the National Oceanic and Atmospheric Admin- istration (NOAA), which is charged with monitoring the Earth-space environment by the U.S. government. X-ray flux is measured by a series of Geostationary Operational Envi- ronmental Satellites (GOES) and flares are given a rating of B, C, M, or X with X flares producing the most powerful emission (see table 1 for details). The data from the GOES satellites are freely available on the Internet3. Very intense flares can produce gamma rays whose energies are measured in mega- electron volts (MeV) . These high energy photons were first detected by balloon borne instrumentation and later studied more in-depth by the Gamma Ray Spectrometer (GRS) on board the SMM satellite (Strong et al., 1999). Many different gamma ray spectral lines were detected, each with different causes. A line at 2.2 MeV is produced when high- energy free neutrons, formed from proton collisions, are slowed down by further collisional 3Data can be obtained from NOAA’s Space Weather Prediction Center (http://www.swpc.noaa.gov/) or from N3KL.ORG (http://www.n3kl.org/sun/noaa.html). 8 process until they are captured by a proton to produce an alpha particle. Particles in flares are raised to such high energies that proton collisions with a nucleus can raise the nucleus’ energy state. A strong line at 4.4 MeV is produced by the de-excitation of carbon nuclei that were excited by the above process. Another strong line is observed in flares at 0.511 MeV. This line is caused by the annihilation of electrons and their positron anti-particles. The anti-particles are produced by the decay of radioactive isotopes, which are produced by the high-speed proton collisions with atomic nuclei (Strong et al., 1999). Energy Release Science is lacking a complete model for flare formation. It is generally agreed that the energy for flares comes from stressed magnetic fields. In a near perfectly conducting plasma, magnetic field lines are tied to plasma motions in regions where the gas pressure exceeds the magnetic pressure. The turbulent motions of the photospheric plasma can twist magnetic fields in the upper corona and stress them into highly non-potential states. Energy is stored in these non-potential states until a critical point is reached. In a process called reconnection, the magnetic topology of the region changes to a more relaxed, or potential, state. The energy liberated in this process is available for thermal heating of the plasma and for the acceleration of particles out of the thermal distribution. One model for reconnection consists of magnetic field lines, that are oppositely aligned, being pushed together in a reconnection region where they recombine to form two new loops, one moving downward and one moving upward. The problem lies in the time scales for reconnection. In order for field lines to reconnect they must diffuse through the plasma to form new connections. The change in a magnetic field in a plasma is governed by Faraday’s law, which in a plasma can 9 be combined with Ampére’s and Ohm’s laws to write the induction equation ∂B ∂ t = ∇× (v×B)+ηB∇2B. (1.3) The first term on the right is how the magnetic field is advected with the plasma and the second term defines the rate that the magnetic field can diffuse through a plasma and is dependent upon the magnetic diffusivity as defined by ηB = c2 4πσ , (1.4) where c is the speed of light and σ is the electrical conductivity of the plasma. For most applications in the solar corona the plasma is considered to be in the ideal MHD limit with σ = ∞, making nB = 0 and reconnection impossible. In order to study non-ideal processes a ratio of the first right hand side term of equation 1.3 to the second term is defined to give the relative importance of each in a given situation. This Magnetic Reynolds Number, also known as the Lundquist Number, RM, is defined by replacing derivatives with l−1 0 where l0 is a characteristic length scale of the system and replacing the velocity with a characteristic velocity scale of the system, v0, to give the following expression RM = l0v0 ηB . (1.5) Typical coronal values for the corona give RM ≈ 108− 1012 (Priest and Forbes, 2000a), which is an indicator of how difficult it is for magnetic fields to diffuse and reconnect in the corona. Via some mechanism, the nature of which is still a matter of intense debate and study, the above limitations are overcome and reconnection occurs. The change in magnetic field implies the generation of an electric field via Faraday’s law. This electric field is a potential 10 acceleration mechanism for the nonthermal electron populations that cause the aforemen- tioned HXR power law spectrum. Numerous works have investigated the direct acceler- ation of electrons by electric fields caused by reconnection. Wood and Neukirch (2005) investigated the acceleration of test particles in a current sheet around a localized recon- nection region. In the presence of a strong guide field in the ignorable direction, it was found that an initially thermal population of particles could be accelerated to a nonthermal beam, with energies comparable to those seen in flares and with a power law energy spec- trum. The hardness of this spectrum was strongly determined by the maximum value of the electric field parallel to the magnetic field. This “simple” model for particle accelera- tion is attractive though there are numerous unresolved issues (See Fletcher et al., 2007a; Fletcher and Hudson, 2008). However, the work described in this thesis is based on the assumption that nonthermal electrons are directly accelerated by electric fields in the re- gion of the reconnection site and not secondary electric fields that may arise from waves, turbulence, etc. While the exact method of particle acceleration may not be fully understood, it is known that the tenuous solar corona is not dense enough to produce the observed emission either in soft or hard X-rays. Material must be coming from some source and then heated to the point of emitting high energy radiation. The most popular theory for this transport of material is chromospheric evaporation. In this model a flare event in the corona heats the underlying layer in the solar atmosphere called the chromosphere either via a thermal con- duction front or by superthermal electron beams colliding with the dense target provided by the chromosphere . The chromospheric plasma is heated to coronal temperatures and then expands and moves upward along magnetic field lines into the corona. There it emits the observed soft X-rays and drains back down into the chromosphere as it cools. This the- ory matches well with the observations of the Russian Intercosmos-4 satellite that observed upward and downward Doppler shifts of flare plasma with velocities of up to 100 kilome- 11 ters per second. Years later Skylab detected similar motions during flares with velocities of 50kms−1 (Strong et al., 1999). The instruments aboard Skylab allowed researchers to make crude electron density estimates of flaring regions. The densities observed were on the or- der of 1012 particles per cubic centimeter (Strong et al., 1999). These densities were far too high for the corona but were similar to chromospheric densities, which added weight to the chromospheric evaporation theory. However, this theory does not perfectly match obser- vations either. The improved instruments on SMM allowed the measurements of Doppler shifts and SXR emission. It was found that high velocity material started moving long before the HXR signature for high-speed electrons encountering the denser chromospheric plasma and continued after the HXR signal ended. If this was the case then the heating of the chromospheric material could not be a result of a flare event in the corona and the chromosphere was somehow independently heated (Strong et al., 1999). Fisher et al. (1985c) modeled two regimes in which chromospheric evaporation could occur. These regions are differentiated by the nonthermal energy flux by electrons with energies over 20keV, F20. Explosive evaporation occurred when energy fluxes F20 ≥ 5× 109ergscm−2 sec−1. In this regime chromospheric material is pushed downward into the lower chromosphere due to the overpressure relative to the denser underlying chromosphere as well as upward. This form of chromospheric evaporation was observed by Milligan et al. (2006a) with the Coronal Diagnostic Spectrometer (CDS; Harrison et al., 1995) onboard the SOHO satellite with measured upward velocities of ∼ 2.70× 107 cms−1 and downflows of 3.5× 106 − 4.5× 106 cms−1. In gentle evaporation, which occurs when F20 ≤ 109ergscm−2 sec−1, a quasi-steady equilibrium between radiative losses and the nonthermal energy input develops allowing the heated material to expand upward in re- sponse to temperature increases. Milligan et al. (2006b) observed an example of gentle chromospheric evaporation, which showed no downflows and with upward velocities of ∼ 1.10×107 cms−1. Both events were also observed by the Reuven Ramaty High-Energy 12 Figure 1.2: A diagram of a flare from Martens and Kuin (1989). Note the lack of a concen- trated looptop emission source. The frequency at which such sources occur was a surprise to theorists and remains largely unexplained. Solar Spectroscopic Imager (RHESSI; Lin et al., 2002) with its unprecedented spectral range and resolution in both the SXR and HXR ranges. Using RHESSI it was possible to determine nonthermal electron flux penetrating the chromosphere for the first time for evaporative events. It was found that the explosive case had an order of magnitude larger electron flux than the gentle case, which was predicted twenty-one years prior by Fisher et al. (1985a; 1985b; 1985c). 13 Looptop X-ray Sources The first evidence of an emission source at the apex of flaring loops was provided by observations of the 1973 June 15 flare by the NRL slit-less objective grating spectrograph onboard Skylab. These observations showed a 1.7×107 K region at the loop apex with the temperature decreasing down the loop length (Cheng, 1977). Unfortunately, Skylab was launched during a time of low solar activity and high resolution images of flares were rare (Acton et al., 1992). The fact that such an emission source existed was a surprise due to the low densities in the corona. Even more surprising was that a majority of flares had looptop emission sources. The Hinotori (Tanaka, 1983) and Solar Maximum Mission (SMM) (Strong et al., 1999) satellites observed numerous flares in concert. Hinotori was able to make images of flares in the 6− 12keV and the 15− 40keV bandpasses. Hinotori observations showed that ex- tended HXR sources are a common feature of two types of flares: impulsive, showing an HXR burst with rapidly varying spikes; and gradual-hard with long (> 30 min) HXR emis- sion with a broad, slowly varying peak and a hard photon spectrum (γ = 5−8) (Tanaka, 1987). In impulsive flares small kernels of HXR emission were seen between the flare ribbons. Observations of limb flares of this type reinforced the interpretation that these ker- nels where actually located at the loop apex. Looptop sources were found to be the major contributor of emission in the 20− 30keV bandpass in the majority of examples (Tanaka, 1987). This seems to contradict the thick-target model of HXR emission which would indi- cate that higher plasma densities are needed for significant HXR emission (Brown, 1971). Whether the emission was thermal or nonthermal in nature remained an open question (Tanaka, 1987). Gradual-hard flares showed an extended source at the top of flare loop arcades throughout the flare. In later phases of the flare, co-spatial SXR sources appeared, with only a few showing a small displacement (≤ 10′′) from the centroid of the HXR emis- 14 sion (Tanaka, 1987). The HXR emission was interpreted as Bremsstrahlung radiation from a nonthermal population of electrons trapped in the coronal portion of the post-flare loops. HXR spectral information for the 1981 May 13 flare was available from SMM’s Hard X-ray Burst Spectrometer (HXRBS; Orwig et al., 1980) and showed a hard HXR photon spectral index of γp ∼ 3.6. Attempts to model the HXR photon emission using a perfect magnetic trap model (Bai and Dennis, 1985) yielded a target plasma density that was two orders of magnitude smaller than the density needed to explain the SXR emission seen by Hinotori (Tanaka, 1987). At the higher density, nonthermal electrons with an energy of < 30keV would be thermalized quickly, necessitating a continuous replenishment of nonthermal par- ticles at the loop apex over the life of the flare (Tanaka, 1987). Resolving this discrepancy remains a central topic in solar flare research. Acton et al. (1992) compiled a survey of ten flares observed with Soft X-ray Telescope (SXT; Tsuneta et al., 1991) on board the Yohkoh satellite (Ogawara et al., 1991). These flares all showed a compact (often only seen in one 2′′ square pixel) emission source at the flare loop apex with a temperature of ∼ 2× 107 K. These sources appeared almost at the onset of the flare and continued well into the decay phase. The sources also remained compact throughout the flare and did not spread down the loop legs. The compactness of the source, its long life and the early onset, made it difficult to explain the increased emission as a density enhancement caused by chromospheric evaporation (Acton et al., 1992). The best known observation of a looptop emission source is the famous Masuda flare of 1992 January 13 (Masuda et al., 1994). This flare was observed simultaneously by SXT and the Hard X-ray Telescope (HXT; Kosugi et al., 1991). The observations clearly showed a SXR (thermal) loop in the SXT thick aluminum and beryllium channels. It also showed three HXR sources, two at the footpoints and one above the SXR loop apex, in the HXT L (14−23 keV ), M(23−33 keV ), and M1(33−53 keV ) bands. This was the first observation 15 of a source above the SXR loop and was thought to be an indication of the reconnection site. The original interpretation of the emission was that the source was a thermal plasma whose temperature was raised to 2×108 K by an unknown mechanism, most likely shock heating (Masuda et al., 1994). This high temperature was an order of magnitude larger than previously seen in flares. The thermal interpretation was initially questioned on the grounds that the source varied on a time scale shorter than the local thermalization time scale (Hudson and Ryan, 1995). It was later determined that a thermal description of the emission was not consistent with the observations when ratios of multiple HXT filters were taken into account (Alexander and Metcalf , 1997). Alexander and Metcalf (1997) also found that HXR emission at both the footpoints and the loop apex could be explained by a single source of nonthermal particles. In this scenario, the loop apex acted as a thick target for electrons with an energy ≤ 20keV and a thin-target to electrons with a higher energy. The HXR emission source then became entirely thermal after the impulsive phase of the flare with a temperature of 2×107 K. Further observations seemingly invalidated the theory that looptop sources are solely the product of chromospheric evaporation. The thermal emission scales as the square of the density increase caused by chormospheric evaporation (Ith ∝ n2 eV , where V is the vol- ume of the emitting region). The increased density also provides a thicker target for non- thermal Bremsstrahlung, which linearly increases the nonthermal emission. The expected morphology of a flare loop would then be a strong brightening of footpoint sources, in both HXR and SXR which would then be followed by a somewhat weaker emission at the loop apex. However, a survey of flares that were observed by HXT and SXT found that this morphology was a rarity, occurring only eight times in the thirty-six flares studied (Nitta et al., 2001). Fletcher and Martens (1998a) numerically modeled the HXR emission of nonthermal electron beams in a converging magnetic trap model. It was found that in a magnetic field geometry of a Syrovatskii current sheet (Syrovatskii, 1971) the apex HXR 16 Figure 1.3: The famous Masuda flare (Masuda et al., 1994) of 1992 Jan. 13. Hard X-rays are shown in contours and the results from soft X-rays are in color. From left to right, the first two columns show images from the SXT thick aluminum and beryllium filters respectively and the temperature and emission measure derived from their ratio in the last two columns. The rows represent images obtained form HXT in three hard X-ray energy bands: 13.9-22.7 keV, 22.7-32.7 keV, and 32.7-52.7 keV respectively. 17 emission was larger and more compact than in previous models which used semi-circular geometries for flare loops. However, this model only applied to the HXR emission. The source of the SXR emission at the loop apex remained unexplained. Sui et al. (2006) made a serendipitous RHESSI observation of X-ray sources moving along the legs of a flare loop. A survey was performed of early impulsive flares, where impulsive hard X-ray emission is seen before a significant rise in soft X-ray emission. This would allow for a spectral analysis of the nonthermal spectrum down to energies of < 10 keV without contamination from thermal spectra. This is important for determining the low energy cutoffs in the nonthermal electron energy spectrum. RHESSI made obser- vations of a C1.1 early impulsive flare on 2002 November 28, near the southwest solar limb, starting at 04:35:30 UT. This was a rare observation since none of the attenuators were in front of the germanium detectors. This allowed RHESSI to make observations at energies as low as 3keV . These low energy observations showed a 3−6keV X-ray source that started at the loop apex, bifurcated, moved down the flare legs at speeds of ∼ 500 and 700kms−1, then moved back upward at an average speed of∼ 340kms−1 (Sui et al., 2006). The initial interpretation of this newly observed motion was that it was an effect of the soft-hard-soft (SHS) spectral evolution of the nonthermal particle beam (Sui et al., 2006). A large number of flares have been observed to have a power law photon spectrum that goes from a large γ (soft), to a small γ (hard), then back to a small γ (Grigis and Benz, 2004; Benz, 2008). If the relationship between γ and δ is assumed to be linear, as the thick and thin target models suggest, then the SHS pattern is suggestive of a nonthermal particle accelerator that shows a similar evolution in time. Sui et al. (2006) proposed that as the nonthermal electron energy spectrum hardens in time, there is a larger population of high energy particles, which can then penetrate more deeply into the loop plasma. This gives rise to an apparent downward motion. Likewise, as the spectrum softens again, the lack of higher energy particles in the distribution restricts the depth that the nonthermal particles 18 can penetrate, giving rise to an apparent upward motion (Sui et al., 2006). However, spec- tral observations of this flare have yet to be published and the ratio of thermal to nonthermal energy contained within these sources remains unknown. Research Motivation While observations from ever more sophisticated instruments have revealed many new aspects of solar flares, these new observations have the annoying habit of raising as many new questions as they answer. Clearly, a true understanding of solar flare processes and features will only be gained by combining observations with models of flares based on theoretical considerations. Many attempts have been made to simulate solar flares numer- ically. For the most part these can be separated into two groups: Those that model the hydrodynamic aspects of a flare (See Mariska et al., 1989; Reeves et al., 2007) while us- ing an analytic expression for heating due to nonthermal electrons, and those that model the evolution of the nonthermal particle beam (See Leach and Petrosian, 1981; Fletcher and Martens, 1998a; McClements and Alexander, 2005) while using a static thermal atmo- sphere. In this work a hybrid loop model is devised that combines models of a hydrodynamic plasma contained within a solar loop and models of nonthermal particle evolution in a self- consistent way. There have been other attempts to achieve this (Miller and Mariska, 2005; Liu, 2006). This work is unique in the way that it treats the nonthermal particles. Previous works treat nonthermal particle evolution via a numerical scheme to solve the Fokker- Planck equation. In this work, nonthermal particles are directly modeled as test particles which undergo stochastic processes as they evolve. This allows the use of measurements of statistical error in the simulations that are not present in previous works. In this work the basic equations involved in the evolution of the thermal plasma and the 19 nonthermal particles are outlined, and the ways in which they are treated numerically are explained. An experiment is then conducted where the effects of including a hydrodynamic plasma model and a nonthermal particle model are compared to a nonthermal model with a static atmosphere in the generation of HXR looptop emission sources. 20 CHAPTER 2 THEORY When dealing with plasmas it is often useful to describe the distribution of a species of particles in a 7 dimensional phase space, f (r,v, t) where r is the position vector, v is the velocity vector. and t is time. The evolution in time of this distribution function is given by taking the total time derivative of f to yield d f dt = ∂ f ∂ t + ∂ f ∂x dx dt + ∂ f ∂y dy dt + ∂ f ∂ z dz dt + ∂ f ∂vx dvx dt + ∂ f ∂vy dvy dt + ∂ f ∂vz dvz dt . (2.1) When d f dt represents the effects of two-body, large angle collisions, Eq. (2.1) is known as Boltzmann’s Equation first derived by Ludwig Boltzmann in 1872 (Lifshitz and Pitaevskii, 1981; Parks, 2004; Sturrock, 1994). By recognizing velocity and acceleration terms in Eq. (2.1) it can be simplified to d f dt = ∂ f ∂ t + v · ∂ f ∂r +a · ∂ f ∂v . (2.2) In the study of astrophysical plasmas, often only electromagnetic forces are considered. Applying Newton’s Second Law and the Lorentz Force Law to Eq. (2.2) yields d f dt = ∂ f ∂ t + v · ∂ f ∂r + q m ( E+ v c ×B ) · ∂ f ∂v , (2.3) where m and q are the particles’ mass and charge, respectively and E and B are the electric and magnetic fields. In the absence of close or binary collisions, Eq. (2.3) becomes the Vlasov Equation first studied by A.A. Vlasov in 1945 (Parks, 2004). An accurate picture 21 of space plasmas can be gained by coupling Eq. (2.3) with Maxwell’s Equations. How- ever, the non-linearity of these coupled equations makes finding solutions a challenging endeavor. One of the difficulties in modeling solar flares is the distribution of particles often con- sists of two distinct components, one thermal and one nonthermal. These two components interact, exchanging energy and momentum, but have drastically different characteristic scale lengths that govern their evolution. Not surprisingly, there are also two broadly differ- ent statistical approaches in plasma physics. The magnetohydrodynamic (MHD) approach assumes that the plasma can be treated as a thermal gas and represented by properties av- eraged over the velocity distribution. The kinematic approach makes no such assumptions and expresses the evolution of the distribution of particles directly instead of by averages. In this work the MHD approach is applied to the thermal component of the distribution plasma while the kinematic approach is applied to the nonthermal component. This chapter is devoted to the description of the new hybrid loop modeling software suite, HYLOOP. HYLOOP is composed of two main components that treat the thermal component and the nonthermal component as separate but interacting distributions of par- ticles. MSULOOP uses the MHD approach in one dimension to calculate the evolution of a thermal plasma. It is the successor to the code used to model X-ray bright points by Kankelborg and Longcope (1999), with modifications made to aid in the study of highly dynamic phenomena such as flares. The Particle Tracking Codes, PATC, is a completely new set of programs designed to use the kinematic approach and track the evolution of nonthermal particle beams. Equations Describing the Thermal Distribution Many numerical modeling codes have been written for coronal loop studies (Dmitruk 22 and Gómez, 1999 & Aschwanden and Schrijver, 2002 are a couple of examples). A recent survey of loop modeling codes (Mulu-Moore et al., 2007) showed that not every code could reproduce the RTV scaling laws that have been a cornerstone in the study non-flaring loops for decades (Rosner et al., 1978). These differences emphasize the fact that numerical loop models are not to be treated as black boxes. The underlying equations and assumptions that comprise any code must be understood in order to determine the regimes where the code’s results have validity. With this fact in mind, the equations and assumptions that go into the hydrodynamic equation solver MSULOOP, will be examined. Thermal Distribution First, a characteristic scale length is derived that will aid in the discussion of a plasma. Consider the electrostatic potential, φ , of a charge surrounded by freely moving particles with an equal number of either charge. The particles will align themselves in a way to minimize the energy of the potential, qφ . Thus the effective potential will be screened and have the form φ = q r exp [ −r λDebye ] , (2.4) where the Debye length, λDebye, has been defined as the e-folding distance, i.e. distance required to reduce the potential by e−1 (Parks, 2004). If a point charge is placed in a neutral distribution of particles, its potential is effectively screened out at distances larger than the Debye length. For a plasma of temperature T and particle density of n the Debye length is given by λDebye = ( kBT 8πne2 )1/2 [cm] (2.5) where kB is the ubiquitous Boltzmann’s constant in cgs units and e is the basic unit of charge in esu units. Assuming a flare plasma with a temperature of T = 107 K and an 23 electron density of ne = 109 cm−3 which is assumed to be half of the total number density, the Debye length is λDebye & 0.2 cm. This is an important length scale that will be used in numerous other definitions and only for scales larger than this length can the plasma be considered to be a neutral fluid. Another useful quantity when determining which equations are applicable to a given distribution of particles is the plasma parameter, gp, which is defined as gp ≡ 1 ND , (2.6) where ND≡ 4πnλ 3 Debye 3 and is the number of particles contained within a sphere with a radius equal to the Debye length. If gp , 1 then there are a large enough number of particles to apply statistical arguments and treat the system as a fluid plasma instead of individual particles. A small value of gp also ensures that the number of binary collisions will be low and an individual particle will interact with the plasma via long range collective effects (Parks, 2004; Choudhuri, 1998). After just a few collisions, an initial distribution of particles will eventually evolve into a Maxwell-Boltzmann distribution which has the form in velocity space of fMB (v)≡ ( m 2πkBT )3/2 exp [ −mv2 2kBT ] , (2.7) The definition of the thermal plasma then becomes a system of particles where deviations from the thermal, Maxwell-Boltzmann distribution are considered to be small and therefore ignored. This occurs in time scales defined by the thermal electron and proton relaxation times. These are the time scales under which slightly non-Maxwellian distributions of elec- trons and protons exchange energy and momentum and relax to a Maxwellian distribution of velocities. The thermal electron and proton relaxation times are estimated numerically 24 by Sturrock (1994) as τMe ≈ 10−1.49 T 3/2 2ne [s] , (2.8) τMp ≈ 10−0.12 T 3/2 2ne [s] , (2.9) where an assumption of a fully ionized hydrogen plasma has been made and the Coulomb logarithm is assumed to be lnΛ ≈ 20 in the solar corona. Again, using T = 107 K and ne = 109cm−3 gives a value of the thermal electron relaxation time of τMe≈ 0.5seconds and a thermal proton relaxation time a τMp ≈ 12seconds. The use of a Maxwellian distribution function is a good approximation for effects that have characteristic time sclaes larger than these relaxation times. A mean free path for an electron can be defined as λm f p ≡ τMevrmse [cm] (2.10) where vrmse is the thermal electron speed given by vrmse = ( 3kBT me )1/2 [ cm s−1] , where me is the electron mass. Again, using the typical flare values defined previously vrmse ≈ 2.13×109cms−1 . Using this yields a mean free path of λm f p ≈ 8.6×108 cm. As long as the characteristic scale lengths of the density, Lne ≡ ( dne dx )−1 and the temperature, LT ≡ (dT dx )−1 , are larger than λm f p then an assumption of local thermal equilibrium is a good one and the approximation that the distribution is Maxwellian can be used (Choud- huri, 1998). MHD Equations With a suitably small value of gp the macroscopic properties of plasma are adequate 25 to describe the state of the system. The macroscopic properties can be found by taking averages over the distribution. These averages are defined by 〈Q〉 ≡ ∫ ∞ −∞ Q f (v)d3v, (2.11) where Q is the variable of interest. The state of a thermal plasma can be determined by three velocity averaged variables; density, n, the mean or bulk velocity, V, and internal energy density, ε as defined by n(x, t)≡ ∫ f (x,v, t)d3v (2.12) V(x, t)≡ 1 n ∫ v f (x,v, t)d3v (2.13) ε (x, t)≡ m 2n ∫ (v−V)2 f (x,v, t)d3v. (2.14) Just as the state variables were defined by taking moments of the distribution func- tion, the equations that govern their evolution are derived by taking successive moments of Boltzmann’s equation to yield, ∂ρ ∂ t +∇ · (ρV) = 0 (2.15) ∂V ∂ t +(V ·∇)V = F− 1 ρ ∇P+ 1 4πρ (∇×B)×B+ 1 ρ ∇ · (µν∇ ·V) (2.16) ∂ε ∂ t +∇ · (εV) =−P∇ ·V+∇ · (κ∇T )−n2 eΛrad(T )+EH + µν (∇ ·V)2 . (2.17) These equations are the continuity equation, Eq. (2.15), the momentum equation, Eq. (2.16), and the energy equation, Eq. 2.17. The evolution of the magnetic field must also be 26 taken into account and is given by the induction equation, ∂B ∂ t = ∇× (V×B)+λ∇2B. (2.18) In order to make the MHD equations more tractable, a variety of simplifications are used in the MSULOOP code. Firstly, the effects of viscosity are ignored in the experi- ments described here and the kinematic viscosity, ν = µν/ρ , has been set to zero. It is acknowledged that viscosity can play an important role in the evolution of flares (Peres and Reale, 1993) and this role will be investigated in subsequent work. Also, the ideal MHD limit is taken which states that the conductivity of the coronal plasma is so high that the diffusivity, λ , can be essentially set to zero. This along with the additional assumption that the magnetic field is reasonably static within the time frame of the simulation, allows for the removal of the induction equation, Eq. (2.18) (McMullen, 2002). The magnetic field, however, is not ignored as the plasma is constrained to flow along the magnetic field lines and defines the coordinate s along the loop length. The MHD equations have now been reduced to 1D hydrodynamic equations with the magnetic field providing the geometry of the coronal loop. In this work any cross-field diffusion effects are ignored and particles are constrained to a characteristic field line for the duration of the simulation. With the magnetic field providing a convenient axis of symmetry, the three dimensional properties of a coronal loop are mapped onto a one dimensional grid. The coronal loop is specified by the coordinate along the loop length, s, the component of acceleration due to gravity along the magnetic field, g‖ (s), the non-uniform cross-sectional area of the loop, A(s) and the three plasma state variables: the electron density ne (s, t), bulk velocity V(s, t), and the internal energy 27 density ε (s, t), without any loss of generality. The equations used by MSULOOP are ∂ne ∂ t =− 1 A ∂ ∂ s (AneV ) (2.19) ∂V ∂ t = −1 nemp ∂P ∂ s +g‖ −V ∂V (2.20) ∂ε ∂ t =− 1 A ∂ ∂ s (AεV )− P A ∂ ∂ s (AV )−ERad + 1 A ∂ ∂ s (AFc)+EH (2.21) In the above equations the divergence operator (∇·) has been replaced with 1 A ∂ ∂ s A in order to ensure that flux conservation laws are strictly obeyed on a nonuniform grid. The conductive flux (FC) and the radiative losses (ERad) are discussed in detail in later sections. The heating term (EH) is not well known in the corona and is usually a free parameter to be tested under simulation. Another assumption being made in the simulations is that the ratio of gas pressure to magnetic pressure is small, β,1, throughout the loop length. Along with an equation of state, P = 2 3ε = 2nekBT , where P denotes the gas pressure, the evolution of the thermal plasma in a loop is completely described with the three predefined state variables. These equations are then solved on a staggered spatial grid (defined below) by using a second order Runge-Kutta Modified Midpoint Method (McMullen, 2002; Press, 2002). This method is not the most computationally efficient. However, Runge-Kutta methods are self-starting, which means that only the initial conditions for the equations need to be known in order to perform the calculation ad infinitum and small errors are not amplified as the calculations progress (Arfken and Weber, 1995). Considering the current trend towards inexpensive computational time and power, the ease of use and stability of the Runge-Kutta method seems to far outweigh its computational efficiency shortcomings. 28 Grids On the Sun, the plasma contained within coronal loops evolves in a continuum of space and time. However, numerical methods cannot deal with continua. The volume that the plasma occupies and the time over which it evolves have to be broken down to discrete chunks to form grids. The definition and size of the grids is a non-trivial matter. Spatial MSULOOP solves Eqs. (2.19)-(2.21) on a one-dimensional, staggered grid that consists of a series of volumes and cross-sectional surfaces that bound those volumes on either side. Since any processes that would allow for particles to diffuse across field lines are not currently included in the simulations, the surfaces perpendicular to the loop axis, as defined by a magnetic field line, are ignored. The volumetric state variables ε and ne are defined within volumes. The bulk velocity, V , the cross-sectional areas, A, parallel components of the gravitational acceleration, g‖, and the magnetic field strength, B are defined on surface grids parallel to the magnetic field. To understand why the effort is spent on defining a staggered grid, which makes book keeping more difficult, one has to look at the numerical methods chosen to deal with the advective portions of the Eqs. (2.19)-(2.21). The simplest stable method for solving an advection type equation is an upwind differencing scheme, which uses different expressions of the differential equation depending upon the sign of the velocity (Morton and Mayers, 1994). When all variables are defined the same grid this is only a first order accurate scheme in space and time. See Appendix B for a discussion on order accuracy. The stability criterion for upwind differencing is the Courant, Friedricks, Lewy (CFL) condition, ∆t ≤ ∆s v , (2.22) which is the condition for convergence for a surprising number of numerical schemes (Mor- 29 Figure 2.1: A visual representation of the staggered grid system. The energy and electron number densities, ε, ne are defined in N volumes. Other quantities, such as velocity, v, magnetic field strength, B, area, A, and parallel component of acceleration due to gravity, g, are defined on N−1 surfaces. The first and last volume serve as static boundary conditions for the grids and are not altered by the simulation. ton and Mayers, 1994). Making smaller step sizes in time, which increases the stability, decreases the spatial accuracy of the result. The accuracy of this calculation can be im- proved by using values interpolated between the volume grids (Press, 2002). By having certain quantities defined on the surface grid and others defined in the volumes bounded by that grid, the upwind differencing is second order accurate with a minimum amount of additional interpolations. It is not possible to have the spatial grids be arbitrarily small and have a simulation finish within a graduate student’s lifetime. Decisions must be made on how to distribute grid spacings in order to yield a physically meaningful simulation within the constraints of the current computing resources. Thin transition region like layers may form in simulations and move with time. These layers might not be resolved with grids set for the corona and lead to numerical inconsistencies (Antiochos et al., 1999). How to define the resolution of a grid is not straightforward and currently there has not been a definitive criterion to base grid size upon. Antiochos et al. (1999) found that defining grids based on the local density 30 worked well due to the density’s strong spatial dependence in both the transition region and the chromosphere. In accordance with this finding, grid spacings are distributed as a Gaussian function ensuring that the smallest grids are at the base of the loop where the density is the largest and the largest grid steps are at the loop apex where the density is at a minimum. The choice of grid spacing can alter the outcome of a simulation. If grids are chosen too coarsely then the local truncation error increases to the point where the errors overwhelm the state variables and the simulation is invalid. If grid spaces are chosen to be too fine then computation time unnecessarily increases. Adaptive mesh refinement techniques exist to calculate grid spacings that minimize both truncation error and run time (See khattri, 2006; Dreher and Grauer, 2006, for examples). Since plasma parameters rapidly vary in solar flares adaptive mesh refinement (AMR) schemes for changing the spatial grid spacing during a simulation were studied. However, it was decided to use a simpler method of determining grid size for the initial use of this new code. Each simulation is initially run on a coarse, static grid with the spacing distributed as a Gaussian function with the maximum percent difference between the step size of one grid point to the next being set to 10%. Successive simulations were run, each having a greater number of cells than the last. State variables are monitored for each run and percent differences are taken of these quantities, comparing each run to a reference grid of arbitrarily small grid spacing. This procedure is repeated until the percent difference fell below a tolerance level defined for each run. The results for these tests are reported in Chapter 3. Temporal It is easy to see how changing a characteristic time scale can alter the fidelity of a simulation. However, just as it was the case with spatial grids, temporal grids cannot be arbitrarily small or the real time of simulation will be longer than is tolerable. The HYLOOP codes have four characteristic time scales, each of which must be properly calculated for 31 a simulation to finish and have some bearing on reality. The time step that governs the nonthermal particles’ evolution, ∆Pt is defined in the section describing nonthermal particle transport. The CFL condition, Eq. (2.22), requires that shocks do not propagate through a grid to another cell before the calculation in the next grid finishes. If this condition is not met then one side of a grid would not be “aware” of what is happening on the other side, and the solution would quickly became unstable (Press, 2002). HYLOOP uses a modified form of the CFL condition and forces the hydrodynamic equations to be iterated forward in time by a step proportional to the following formula ∆gt ≤ min [ ∆gs(i) cs(i)+ |v(i)| ] (2.23) where i is the cell index, ∆gs(i) is the length of the ith spatial grid cell, cs (i) is the local sound speed, and |v(i)| is the absolute value of the velocity of the thermal plasma. Equation (2.23) is an upper bound for a stable simulation. The time step derived in Eq. (2.23) is divided by a user defined safety factor, S0, usually set to 5 in order to increase the stability and fidelity of the upwind differenced quantities (Press, 2002). It is easy to see from Eq. (2.23) how refining the spatial grid can quickly lead to an inordinate increase in run time. When the spatial grid scale is reduced, either locally or globally, not only must the code solve the hydrodynamic equations on more spatial grids but the time step is also reduced, leading to even more calculations per reporting time step. The next time step defines the rate at which information is passed between MSULOOP and PATC. This communication time step, ∆Ct, is user defined. Ideally this time step would be the same as ∆gt. However, due to the large number of test particles used in each simula- tion and the number of calculations that have to be performed per particle, the run times for ∆Ct = ∆gt became prohibitively large on the current computing resources even with mod- 32 est improvements in parallel processing capability and efficient code writing techniques. In general practice ∆Ct ∼ 10−3 sec. Particles are allowed to propagate along a static plasma during this time. The energy and momentum loss of the particles then become input terms for the hydrodynamic plasma. The plasma is allowed to evolve until it catches up in time with the particle simulation. The particles then have an updated target plasma and the process repeats. It should be noted that the value of ∆Ct used in HyLoop represents an improvement of three orders of magnitude from the previous state of the art codes (Liu, 2008). The next time step to be defined is the reporting time step, ∆Rt. This is the largest of the time scales and is user defined. It would strain the capacity of the available storage disk infrastructure to maintain records of all state and particle variables on the smallest time scales used to solve the myriad of equations contained within HYLOOP. Instead, all of the information calculated by the codes is written to disk every ∆Rt of simulated time. While having this time scale be larger than the other characteristic time scales does entail a loss of information, it is deemed a necessary loss. Not only would many gigabytes of data be compiled by the codes if the reporting time step is not orders of magnitude larger than the other time scales, but since the writing of information to disk still involves the moving of physical parts for most systems, the run time would also increase significantly. While this time scale does not affect the fidelity of the simulation, it can affect the estimation of thermal emission, which is calculated post simulation. Radiative Losses Corona In the corona, ions are excited primarily by collisions and de-excited by via dielectronic recombination and the emission of photons (Aschwanden, 2005). For the most part, the simulations are dealing with an optically thin plasma and the photons are free to 33 radiate into free space without re-absorption. This results in a net loss of energy from the corona. The description of the radiative losses of an optically thin plasma is given by ERad = n2 eΛRad (T ) . (2.24) To determine radiative losses a radiative loss curve, ΛRad (T ), has to be computed in a man- ner that is both numerically efficient and yields a good approximation to the true radiative losses in the solar corona. A piecewise-continuous function expressing the radiative losses of an optically thin plasma is calculated using the CHIANTI spectral database package (Dere et al., 1997) in the SOLARSOFT software distribution. CHIANTI calculates spectral and continuum photon emission with a given set of elemental abundances. It is well known that the abundances of elements in the corona differ significantly than their photospheric values and that this variation is strongly dependent upon the first ion- ization potential (FIP) of each element. Since hydrogen is fully ionized in the corona, it is usually not possible to determine the absolute abundances, i.e., abundances with respect to hydrogen (see White et al., 2000 for a counter example). Instead, only measurements of relative elemental abundances are possible. Relative abundances are calculated as a ratio between the fraction of a given element and another element, whose abundance is also in doubt. This leads to much debate on whether the abundances of high FIP elements are be- ing depleted by a factor of ∼ 4, low FIP elements are being enhanced by a factor of ∼ 4, or a combination of enhancement and depletion (see Feldman and Widing, 2003 for a review). Recently, a mechanism that explains the observed coronal abundance variations has been proposed that has been very successful in matching observed values (Laming, 2004). This model supports the view that, on average, high FIP elements are depleted by a factor of∼ 4 in the corona as suggested by Meyer (1985). The current set of simulations use a radiative loss curve calculated with Meyer abundances under an assumption of constant pressure as 34 described by Martens et al. (2000). Even though abundances vary in time and from place to place (Feldman and Widing, 2003) it is asserted that this average treatment is suitable for modeling the radiative loss curve and it should be noted that this treatment is consistent with other state of the art coronal loop models (Klimchuk et al., 2007). Chromosphere and Transition Region It was shown by McClymont and Canfield (1983) that optical depth effects in the transition region play an important role in deter- mining the stability of coronal loops. This complicates things greatly from the optically thin case since the nonlocal nature of the opacity means that a full treatment of radiative losses in the lower transition region and chromosphere requires a full 3D radiative hydro- dynamics code (Allred et al., 2005). There are other effects that conspire to make radiation loss calculations more difficult in lower temperature regions. The assumption of a fully ionized plasma, which is implicit in the calculation of Eq. (2.24), is not valid for solar plasmas with temperatures ≤ 20,000K, which can be composed of 10% neutral hydrogen and thus drastically increase the radiative losses (Fontenla et al., 1991). In the ambipolar diffusion process, electric fields are set up in weakly ionized plasma that counteract charge separation and ensure that ions and electrons diffuse from regions of high density to low density at the same rate (Choudhuri, 1998). In gravitationally stratified atmospheres this has the curious effect of setting up an electric field that “buoys up” heavier ions and “weighs down” electrons. One effect of this is that the plasma scale height is twice that expected by a non-ionized gas made up of the same ions (Kivelson and Russell, 1995). Also, neutral species can be found at a temperature higher than could be predicted by local ionization rates (Fontenla et al., 1991). A proper treatment of the radiative loss curve requires delicate energy balance calcula- tions that include the effects of radiative transfer that include opacity, ambipolar diffusion and mass motions. Fortunately, this daunting task is tackled by Fontela, Avrett, and Loeser 35 Figure 2.2: The radiative loss function as used by the MSULOOP code for the current set of experiments. The diamonds represent data points and the line shows interpolated values. (1991, hereafter FAL) who published the radiative loss rates and conductive fluxes for a so- lar plasma that includes the previously mentioned effects. The first use of these improved radiative loss rates in a numerical loop model was by HYLOOP’s predecessor as described in Kankelborg and Longcope (1999). The results from Eq. (2.24) and the corrections from FAL are tabulated as a piecewise continuous function of temperature and saved as a system variable that all programs can access and then interpolated for each temperature calculated on the volume grid. The total radiative loss curve for the experiments conducted in this work is shown in figure 2.2. 36 Radiative Loss Time Scale In order to choose the most effective method of dealing with equations numerically, it is often instructive to examine their characteristic time scales. To do this, the energy equation, Eq. (2.21) is rewritten using the equation of state to put the right hand side in terms of temperature instead of internal energy. All left hand side terms are set to zero except for the radiative loss term, in the form of Eq. (2.24) yielding ∂T ∂ t ≈ neΛRad (T ) 3kB . (2.25) Rearranging Eq. (2.25) produces the following estimation of the radiative cooling time, in seconds at a given temperature, T0, τRad ≈ T0 3kB neΛRad (T0) . (2.26) A flare plasma with ne = 109 cm−3, and T0 = 107 K yields a radiative cooling time of τRad ≈ 1.1×105 s. Since this time scale is orders of magnitude larger that any of the temporal grid scales, ∆Pt, ∆Rt, ∆gt, ∆Ct, it is determined that this term could be handled by a single step Euler’s method which is only first order accurate but computationally inexpensive (Press, 2002). Notes on Radiative Losses For more complete treatments on the calculation Eq. (2.24) the interested reader is encouraged to review Jefferies et al. (1972), Rybicki and Lightman (1986), and Landi and Landini (1999). It should also be mentioned that all of the radiative loss equations for the corona are calculated under the assumption of local thermal equi- librium (LTE), where matter is in thermal equilibrium with itself if not adjacent cells or the radiation it emits and in a state of ionization equilibrium. For the corona, ionization equilibrium means that collisional ionization rates have come to equilibrium with radiative 37 and dielectric recombination rates Aschwanden (2005). However, it has been shown that in cases of rapid heating, such as flares, or rapid cooling, the ionization equilibrium assump- tion is not a good one (Bradshaw and Mason, 2003). Since these non-equilibrium effects have a minimal impact on spectrum averaged radiative loss rates and therefore the evolu- tion of the thermal plasma (Klimchuk, 2006), they are ignored for now. This does limit the current study to only using spectrally averaged quantities as observable outputs. Conduction Many loop s