.. 72:. v 3.9.... _ .. A... 4... c}- at) . .1“ 3.... . ._ . 1.... 2 ei l - _l.’li‘i.', ‘.V émbvr 3“ ii... h. , garmmmr x. rankmrt ‘ 1 .. . id. 1.9.". . ‘ 3.x 9.5.. .: u u.) a: I. \5 4 ~ .1 . - a Vl.!. \v .qu1 s. . “TE 0 LIBRARY N Michigan State University This is to certify that the dissertation entitled A KINETIC THEORY BASED NUMERICAL STUDY OF CORE COLLAPSE SUPERNOVA DYNAMICS presented by TERRANCE T. STROTHER has been accepted towards fulfillment of the requirements for the Ph.D. degree in Physics Lfifim 3L0 Wfichrofessor’ 5 Signature December 11, 2009 Date MSU is an Affirmative Action/Equal Opportunity Employer -u-o----------u- PLACE IN RETURN Box to remove this checkout from your record. TO AVOID FINES return on or before date due. MAY BE RECALLED with earlier due date if requested. DATE DUE DATE DUE DATE DUE 5108 K1IPrCIjIAoc8ProsICIRCIDateDue.hdd A KINETIC THEORY BASED NUMERICAL STUDY OF CORE COLLAPSE SUPERNOVA DYNAMICS By Terrance T. Strother .A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Physics 2009 ABSTRACT A KINETIC THEORY BASED NUMERICAL STUDY OF CORE COLLAPSE SUPERN OVA DYNAMICS By Terrance T. Strother The explosion mechanism of core collapse supernovae remains an unsolved problem in astrophysics after many decades of theoretical and numerical study. The complex nature of this problem forces its consideration to rely heavily upon numerical simula- tions. Current state-of-the—art core collapse supernova simulations typically make use of hydrodynamic codes for the modeling of baryon dynamics coupled to a Boltzmann transport simulation for the neutrinos and other leptons. The results generated by such numerical simulations have given rise to the widely accepted notion that neu- trino heating and convection are crucial for the explosion mechanism. However the precise roles that some factors such as neutrinos production and propagation, rota- tion, three-dimensional effects, the equation of state for asymmetric nuclear matter, general relativity, instabilities, magnetic fields, as well as others play in the explo- sion mechanism remain to be fully determined. In this work, we review some of the current methods used to simulate core collapse supernovae and the various scenarios that have been developed by numerical studies are discussed. Unlike most of the numerical simulations of core collapse supernovae, we employ a kinetic theory based approach that allows us to explicitly model the propagation of neutrinos and a full ensemble of nuclei. Both of these are significant advantages. The ability to explicitly model the propagation of neutrinos puts their treatment on equal footing with the modeling of baryon dynamics. No simplifying assumptions about the nature of neutrino-matter interactions need to be made and consequently our code is capable of producing output about the flow of neutrinos that most other simulations are inherently incapable of. Furthermore, neutrino flavor oscillations are readily incorporated with our approach. The ability to model the propagation of a full ensemble of nuclei is superior to the standard tracking of free baryons, a particles. and a “representative heavy nucleus”. Modeling the weak reactions that free baryons and hundreds of species of nuclei undergo results in a more realistic evolution of the nuclear composition. The explicit knowledge of nuclear composition at all times not only allows us to study its evolution in greater detail than it has before, but it also puts us in the unique position to directly model certain nuclear decay modes and the effects that nuclear structure have on non-weak nuclear reaction that occur in supernovae quite straightforwardly. A systematic study of the influence that electron capture rates and the nuclear equations of state have on the collapse and explosion phase is conducted. The algo- rithmic implementations and motivations for using the various values and expressions for the electron capture rates and nuclear equations of state are explained and the new forms of output that our code is singularly capable of producing are discussed. Dynamics that may prove to be an entirely new neutrino capture driven explosion mechanism were observed in all of our simulations. ACKNOWLEDGMENT Most importantly, I would like to thank Professor Wolfgang Bauer for being my advisor for this research project. Without his willingness to share his expertise and his patients when doing so, I would not have been able to acquire the knowledge and skill set required to be a productive research scientist. I would also like to thank my parents and grandparents for supporting me throughout my journey at Michigan State University, and my wife and two children for surviving it. iv TABLE OF CONTENTS List of Tables ................................. List of Figures ................................ Introduction ................. Review of Core Collapse Supernovae ............... 2.1 Review of Pre-Supernova Evolution ................... 2.2 Review of Core Collapse and Bounce .................. 2.3 The 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 2.2.1 2.2.2 2.2.3 2.2.4 2.2.5 Proposed Explosion Mechanisms ................. The Prompt Shock Mechanism ................. The Neutrino Revival Mechanism ................ The Convection Revival Mechanism ............... The Acoustic Vibration Mechanism ............... Recent Numerical Studies ........................ 2.3.1 2.3.2 2.3.3 2.3.4 2.3.5 Test Particle Approach Equation of State ......................... Weak Interaction Rates ...................... One-Dimensional Simulations .................. Two-Dimensional Simulations .................. Three-Dimensional Simulations ................. Origins of the Test Particle Approach .................. Modeling Macroscopic Systems With Test Particles .......... Modeling Supernovae with Test Particles ................ Matter Test Particles ........................... Matter Test Particle Initial Conditions ................. 3.5.1 3.5.2 Initial Matter Test Particle Momentum Space Coordinates . . Initial Matter Test Particle Nuclear Properties ......... Neutrino Test Particles .......................... Neutrino Test Particle Initial Conditions ................ Data Management ............................ 3.8.1 3.8.2 3.8.3 Spherically Symmetric Spacial Test Particle Grouping ..... Three-Dimensional Spacial Test Particle Grouping ....... Grouping Matter Test Particles By Nuclear Properties Calculating Distributions ......................... 3.9.1 Spherically Symmetric Distributions .............. Radial Derivatives of Spherically Symmetric Distributions ix Linearly Interpolating Spherically Symmetric Distributions . . 3.9.2 T hree-Dimensional Distributions ................ 3.10 Matter Test Particle Motion ....................... 3.10.1 Gravitation ............................ Spherically Symmetric Gravitation ............... Three-Dimensional Gravitation ................. 3.10.2 Nucleonic Force .......................... Symmetric Mean Field N ucleon Potentials ........... Asymmetric Mean Field Nucleon Potentials .......... 3.10.3 Electron Pressure ......................... 3.10.4 Matter Test Particle Scattering ................. 3.11 Neutrino Test Particle Production and Propagation ................................ 3.12 Neutrino Test Particle Production . .- .................. 3.12.1 Testing For Neutrino Test Particle Creation .......... 3.12.2 Creating a Neutrino Test Particle ................ 3.13 Neutrino Test Particle Matter Interactions ............... 3.13.1 Neutrino Oscillations in the Core ................ 3.13.2 Testing for Neutrino-Matter Interactions ............ 3.13.3 Modeling Neutrino Test Particle Captures ........... 3.13.4 Modeling Elastic Neutrino-Matter Interactions ......... 3.13.5 Neutrino Test Particle Movement ................ 3.14 Updating Matter Test Particle Nuclear Properties ................................. 3.15 Updating Matter Test Particle Temperatures .............. 3.16 Weak Interaction Induced Temperature Changes .................................. 3.17 Heat Exhange ............................... 3.18 Fusion ................................... Electron Gas Statistical Mechanics Table ............ 4.1 Arbitrary Temperature Statistical Mechanics .............. 4.1.1 Determining 5 ........................... 4.1.2 Pressure Integral ......................... 4.1.3 Partial Derivatives of Pressure .................. 4.1.4 Kinetic Energy Density and its Partial Derivatives ....... 4.2 Low Temperature Statistical Mechanics ................. 4.2.1 General Sommerfield Expansions ................ 4.2.2 Low Temperature Chemical Potential .............. 4.2.3 Low Temperature Kinetic Energy Density ........... 4.2.4 Low Temperature Pressure .................... 4.3 High Temperature Statistical Mechanics ................ 4.3.1 High Temperature § Parameter ................. 4.3.2 High Temperature Kinetic Energy Density ........... vi 48 48 50 50 53 55 55 57 60 62 64 67 68 69 81 82 83 87 88 94 96 97 98 100 101 106 107 110 112 113 114 4.3.3 High Temperature Pressure ................... 115 4.4 Electron Gas Statistical Mechanics in the Non-Relativistic Limit . . . 117 4.4.1 Non-Relativistic Limit of Arbitrary Electron Gases ...... 117 4.4.2 Non-Relativistic Limit of Low Temperature Electron Gases . . 119 4.4.3 Non-Relativistic Limit of High Temperature Electron Gases . 121 4.5 Comparing Relativistic and Non-Relativistic Electron Gas Statistical Mechanics ................................. 123 4.5.1 Numerical Non-Relativistic Convergence ............ 124 4.5.2 Relativistic and Non-Relativistic Electron Gases in Supernovae 128 Weak and Strong Reaction Tables . . . . . . . . . . . . . . . . 136 5.1 Electron Capture Rate Tables ...................... 136 5.2 Neutrino Production and Interaction Tables .............. 138 5.2.1 Neutrino Production Tables ................... 138 5.2.2 Relativistic Thermal Energy Distributions ........... 139 Relativistic Energy Distributions of Arbitrary Gases ...... 140 Electron Gas Energy Distribution ................ 142 Free Baryon and Nuclear Energy Distributions ......... 142 Maxwell Limit ........................... 144 Comparing Relativistic and Non-Relativistic Energy Distribu- tions ........................... 145 5.2.3 Neutrino Capture Tables ..................... 149 5.2.4 Elastic Neutrino-Baryon/ Nucleus Interaction Table ...... 153 5.2.5 Elastic Neutrino-Electron Interaction Table .......... 155 5.3 Strong Interaction Tables ......................... 157 Computational Requirements . . . . . . . . . . . . . . . . . . 159 6.1 Primary Source Code ........................... 159 6.2 Input Generating Codes ......................... 161 6.3 Output Manipulating Codes ....................... 162 Diagnostic Tests . . ..... . . . . . . . . . . . . . . . . . . 163 7.1 Data Management ............................ 164 7.2 Energy Conservation ........................... 164 7.2.1 Energy Conservation by Elastic Scattering ........... 165 7.2.2 Energy Conservation by Gravitation .............. 165 7.2.3 Energy Conservation by Local Forces .............. 167 Degenerate Electron Gas Pressure ................ 170 Degenerate Electron Gas Forces ................. 173 Degenerate Electron Gas Average Potential Energy ...... 176 Testing the Energy Conservation of Local Forces . . . . . . . . 178 vii 8 Numerical Results ........................ 182 8.1 New Dynamics Observed ......................... 184 8.1.1 Evolution of the Electron Fraction ................ 184 8.1.2 A Possible New Explosion Mechanism ............. 192 8.1.3 Explosion Dynamics ....................... 200 8.2 Genesis of New Dynamics ........................ 210 8.2.1 Qualitative Analysis of Genesis ................. 210 8.2.2 Neutrino Capture Probability Distribution ........... 212 8.2.3 Evolution of the Nuclear Composition ............. 219 8.2.4 Neutrino Energy Distribution .................. 226 8.3 Role of the Nucleon Potential ...................... 229 8.4 Role of Electron Capture Rates ..................... 232 8.5 Robustness of Explosion Mechanism with Respect to Variation of N u- merical Parameters ............................ 234 9 Conclusion ..................... . ...... 238 Bibliography .............. . ..... . . . . . . 242 viii 8.1 8.2 8.3 8.4 8.5 LIST OF TABLES Core characteristics at the time the outward explosion of matter formed as calculated by simulations that used different nucleon potentials. . . Proto-remnant data generated by simulations that used different nu- cleon potentials ............................... Core characteristics at the time the outward explosion of matter formed as calculated by simulations that used different electron capture rates. Proto—remnant data generated by simulations that used different elec- tron capture rates. ............................ Core characteristics at the time the outward explosion of matter formed as calculated by simulations that used different numbers of spherical shells to calculate statistical distributions. ............... ix 230 231 232 233 2.1 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 LIST OF FIGURES Images in this dissertation are presented in color Schematic depiction of the core mass increasing as the layers of nuclear burning above it continue to produce nuclear ash. ........... Table of stable and unstable nuclei currently included in our simulation. 35 Plot of the initial density and temperature distributions of the progen- itor as functions of enclosed solar mass taken calculated by Woosely and Weaver [18]. ............................. Plot of the initial angular velocity distributions calculated by Heger and Langer for their 15 MO rotating progenitor taken from a private communication from Chris Fryer. .................... Initial chemical composition of the progenitor as a function of enclosed solar mass taken from Woosely and Weaver [18] ............. The soft and stiff BKD isoscalar potentials plotted over the density range 0 S p/pO S 3 ............................. The lowest order expansion of U (p, 6) plotted over the density range 0 S p/pO S 3 for the cases 6 = 0, 0.2, and 0.4 .............. Schematic depiction of how two-body elastic scattering is modeled. Neutrino transition probability amplitudes for the P1_,2 transition plotted as a function of [U(k)/lo(ne) over the range 0.1 g lv(k)/lo(ne) g 10. ..................................... 38 39 40 56 59 63 3.9 3.10 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 Plot of the neutrino transition probability amplitudes for the P1_,2 transition for 0.1, 1, 10, and 100 MeV neutrinos 7progagati n§4 through matter with electron densities 1n the range 1037 to 104 Schematic depiction of the neutron drip line nucleus 458 capturing an electron and emitting 3 neutrons to become a 42F nucleus, the nucleus with the largest A S 45 with Z = 15 currently included in our table of nuclei ................................... Plot of the 6 parameter for relativistic and non-relativistic electron gases with 36temperature T=107 K and number densities 1033m 3 S ne < 1036 m ............................. Plots of 8P/ 8T and 811 / 8T for relativistic and non-relativistic electron gases with 36tem§erature T=107 K and number densities 1033m 3 g 718 < 1036 ............................. Plots of 8P/8n and 811/871 for relativistic and non-relativistic electron gases with 36tnelmgerature T=107 K and number densities 1033m 3 S 718 _<_ 1036 ............................. Plot of the g parameter for relativistic and non-relativistic electron gases with 5temperature T=109 K and number densities 1037 m 3 g 716 < 1045 ............................. Plots of 8P/ 8T and 8a / 8T for relativistic and non-relativistic electron gases with 5temperature T=109 K and number densities 1037 m"3 _<_ ne < 104‘5 ............................. Plots of 8P/8n and 811 / 8n for relativistic and non-relativistic electron gases with 5temgerature T=109 K and number densities 1037m "'3 < ne _<_ 104!5 ............................. Plot of the § parameter for relativistic and non-relativistic electron gases with 45temtgerature T— — 1011 K and number densities 1037 m 3 < ne<1045 ............................. Plots of 8P/ 8T and 811/ 8T for relativistic and non-relativistic electron gases with tern erature T = 1011 K and number densities 1037 m_3 S 116 S 1045 m’ . ............................. xi 74 82 4.9 5.1 5.2 5.3 5.4 7.1 7.2 7.3 8.1 8.2 8.3 8.4 8.5 Plots of 8P/8n and 811/871 for relativistic and non-relativistic electron gases with 45temt§erature T— - 1011 K and number densities 1037m “3 < <10 m ............................. Schematic depiction of how two-body inelastic scattering is modeled. Plot of the relativistic energy distribution and non-relativistic 3energy distribution for a gas of electrons with number density 716— — 1033 m and temperature T— - 107 K expressed as functions of velocity. . . . . Plot of the relativistic energy distribution and non-relativistic 3enelrgy distribution for a gas of electrons with number density 116— — 103 and temperature T— — 1010 K expressed as functions of velocity. Plot of the relativistic energy distribution and non-relativistic energy distribution for a gas of non-degenerate protons with temperature T— - 1011 K expressed as functions of velocity. ............. Plot of the matter test particle kinetic energy for toy model with no forces only scattering ............................ Plot of the matter test particle kinetic, gravitational potential, and total energy for a toy model that only includes gravity. ........ Plot of the matter test particle kinetic, local interaction potential, and total energy for a toy model that only includes local forces ....... Plot of the initial electron fraction in the spherical shells ........ Plot of the electron fraction in the spherical shells shortly after the inner region of the core’s rate of deleptonization began to increase. . . Plot of the electron fraction in the spherical shells shortly after the inner region of the core has begun to rapidly deleptonize ........ Plot of the electron fraction in the spherical shells shortly after the spike develops. .............................. Plot of the electron fraction in the spherical shells showing fully relep- tonized matter in the spike having an electron fraction approximately equal to matter in the outermost shells where the effects of electron capture are negligible. .......................... xii 146 147 148 166 179 186 187 188 8.6 8.7 8.8 8.9 8.10 8.11 8.12 8.13 8.14 8.15 Plot of the electron fraction in the spherical shells showing hyper- leptonized matter in the spike having an electron fraction greater than matter had at any point in the core before the collapse began. Plot of the initial electron fraction and 57‘ of matter in the spherical Shells. ................................... Plot of the electron fraction and Br of matter in the spherical shells shortly after the inner region of the core’s rate of deleptonization began to increase. ................................ Plot of the electron fraction and Br of matter in the spherical shells 191 194 shortly after the inner region of the core has begun to rapidly deleptonize.196 Plot of the electron fraction and Br of matter in the spherical shells at the time the outward explosion of matter forms shortly after the spike develops and the matter between the Spike’s peak and the outer edge of the rapidly deleptonized inner region of the core where the electron fraction is depressed has clearly decelerated. .............. Plot of the electron fraction and Br of matter in the spherical shells showing fully releptonized matter in the spike having an electron frac- tion approximately equal to matter in the outermost shells where the effects of electron capture are negligible and a large outward moving density wave is visible. .......................... Plot of the electron fraction and Br of matter in the spherical shells showing hyper-leptonized matter in the spike having an electron frac- tion greater than matter had at any point in the core before the col- lapse began and the density wave has successfully propagated almost completely out of the core ......................... Plot of the initial electron fraction, Br, and density in the spherical shells. ................................... Plot of the electron fraction, ,Br, and density in the spherical shells shortly after the inner region of the core’s rate of deleptonization began to increase. ................................ Plot of the electron fraction, Br, and density in the spherical shells shortly after the inner region of the core has begun to rapidly delep- tonize and the density inversion in the innermost region of the core has formed. .................................. 197 198 199 204 205 8.16 8.17 8.18 8.19 8.20 8.21 8.22 8.23 8.24 Plot of the electron fraction, Br, and density in the spherical shells at the time the outward explosion of matter begins shortly after the spike develops and the matter between the Spike’s peak and the outer edge of the rapidly deleptonized inner region of the core where the electron fraction is depressed has clearly decelerated. .............. Plot of the electron fraction, fir, and density in the spherical shells showing fully releptonized matter in the spike having an electron frac- tion approximately equal to matter in the outermost shells where the effects of electron capture are negligible and a large outward moving density wave is seen beginning its separation from the dense proto- remnant ................................... Plot of the electron fraction, Br, and density in the spherical shells showing hyper-leptonized matter in the spike having an electron frac- tion greater than matter had at any point in the core before the col- lapse began and the density wave has successfully propagated almost completely out of the core and has clearly separated from the dense proto—remnant. .............................. Plot of the neutrino-matter interaction probabilities in the spherical shells before the innermost region of the core has begun to rapidly deleptonize. ................................ Plot of the neutrino-matter interaction probabilities in the spherical shells shortly after the innermost region of the core has begun to rapidly deleptonize. ................................ Plot of the neutrino-matter interaction probabilities in the spherical shells after the rapid deleptonization of the innermost region of the core has led to a large abundance of free neutrons ............ Plot of the neutrino-matter interaction probabilities in the spherical shells shortly at the time the outward explosion of matter begins and the electron gas in the inner portion of the innermost region has become very degenerate ............................... Plot of the initial nuclear composition of the core ............ Plot of the nuclear composition of the core after its inner region has be- gun to rapidly deleptonize, but before the spike in the electron fraction has started to develop. .......................... xiv 207 208 209 215 216 217 218 221 222 8.25 Plot of the nuclear composition of the core at the time the outward explosion of matter begins ......................... 8.26 Plot of the number of free baryons per “heavy” nucleus (A > 1) in the core after its inner region has begun to rapidly deleptonize, but before the spike in the electron fraction has started to develop ......... 8.27 Plot of the number of free baryons per “heavy” nucleus (A > 1) in the core at the time the outward explosion of matter began. ....... 8.28 Plot of the average neutrino energy distribution in the core at the time the outward explosion of matter began .................. 8.29 Plot of the average neutrino energy distribution in the core after a large outward moving density wave has formed and is beginning its separation from the dense proto—remnant ................. XV 223 224 225 227 228 Chapter 1 Introduction Supernova explosions are believed to be one of the main sources for heavy element production (A > 56) in the universe. The origins of the heavy elements and the role that supernova explosions play in their creation remains an area of intense research [1-14]. Studies of supernovae are conducted in two complementary ways, through construction and operation of rare isotope accelerators such as the Facility for Rare Isotope Beams (F RIB), and through astronomical observation of supernova explosions and their numerical modeling. It is the latter subject to which this thesis attempts to contribute. In chapter 2, a brief review of the generally accepted picture of core collapse super- novae is presented. Due to the complex nature of the supernova problem, numerical simulations play a critical role in understanding the specific mechanisms that drive the core collapse and subsequent supernova explosion. We therefore include a sue- cinct overview of recent numerical supernova simulations in this chapter as well. we pay particular attention to those that influenced our work. This is not intended to be a complete overview of recent numerical supernova simulations, rather its purpose is to highlight the different approaches that groups have taken. In chapter 3, we introduce the approach that we used to model core collapse 1 supernovae. Here the algorithmic implementations, numerical techniques used, and approximations and assumptions made are all explained in detail. We also discuss the code’s capabilities and limitations on a single processor as well as what it can do when ported to a massively parallel computer cluster. In chapter 4, we rigorously derive the expressions from statistical mechanics used to characterize the electron gas and describe the numerical evaluation and tabulation of said quantities. All derivations make use of the fully relativistic formalism. In chapter 5, we explain how all of the weak processes currently included in our simulation are modeled, present the formulae used for the associated weak reaction cross sections, and describe the tabulation of all of the quantities needed by the simulation to model the net effects that weak reactions have on neutrino test particle dynamics and the temperature distribution. All particle scatterings are modeled using semi-classical relativistic mechanics. In chapter 7, we describe the ways in which the code has been tested for potential numerical problems. Both the dynamic data management tests as well as the external tests of the numerical implementation of the physical models in our code are explained and their results are discussed. In chapter 8, we discuss the results generated by several simulations with some parameters systematically varied in such a way that we can draw conclusions about the dependencies the collapse and explosion have upon them. These parameters are discussed at length in both in this chapter and chapter 3. The neutrino capture induced explosion mechanism observed in all calculations is discussed in detail and differences between the results are explored as well. In chapter 9, our work is summarized. The phenomenon seen in all calculations that may prove to be the first glimpse at a new explosion mechanism presented in chapter 8 is reflected upon and the validity of these calculations is discussed. The short term and long term plans for our code are given as well. 2 In this thesis, we will make use of standard symbols in the body of the text and the equations displayed, such as c for the speed of light, ME) for the solar mass, G for Newton’s gravitational constant, It for Plank‘s constant, it for or Plank’s constant divided by 27r, k for the Boltzmann constant, and A for the atomic mass number of a nucleus. Other abbreviations and symbols used are explained as they are introduced. The primary source code used to generate the results presented in this thesis is written in C as are all of the supporting codes that create the input files it needs to run and process its output. Due to the primary source code’s length, approximately 25,000 lines, it is not included as part of this thesis. Instead electronic copies are available at http: / / www.pa.msu.edu / ~bauer / code / SuperNova / . Chapter 2 Review of Core Collapse Supernovae In this section, we divide our brief reviews of the physics of core collapse supernovae and the recent efforts to numerically model them in the following way: First we present a short summary of the evolution of a core collapse supernova candidate star. Then we give a synopsis of the accepted picture of the core collapse and bounce followed by discussions of the proposed driving mechanisms for the subsequent explosion. Finally we discuss the different approaches that selected groups have taken, in the last ~ 15 years, taken to simulate core collapse supernovae numerically. Again we stress that the reviews presented here are not complete. It is far beyond the scope of this work to exhaustively discuss all of the relevant physical processes here as well as the efforts made to numerically model them. For in—depth analyses on any of these topics, the reader is referred to the associated citations. In addition, there are other groups that are not mentioned here for brevity’s sake. These groups are undoubtedly making meaningful contributions to the field. The works discussed here were chosen due in part to relevance to our own wok and also because we feel that they represent a broad sampling of the different techniques used to simulate core collapse supernovae. 4 2.1 Review of Pre-Supernova Evolution Stars spend most of their lifetimes in hydrostatic equilibrium. During this time, nu- clear fusion reactions synthesize heavier elements from lighter elements. Each time a fusion reaction occurs, the mass of the fusion product is less than that of the fuel. The mass-defect is changed into thermal energy that counteracts the gravitational force and allows for a hydrostatic equilibrium in the star. When the nuclear burning fuel at the center is exhausted, the thermal pressure decreases and the star expe- riences gravitational contraction. Due to this gravitational compression the central temperature rises until the temperature becomes sufficiently high to ignite the next nuclear burning phase. This sequence of nuclear burning to central fuel exhaustion, contraction, and ignition of a next burning phase repeats itself a number of times that depends on the initial mass of the star. If the star’s initial mass was great enough, this cycle occurs until the fusion of the nuclear fuel in the core is no longer energetically advantageous. At the end of a star’s hydrostatic life, it is left with an onion-skin-like stratification where each layer consists of the ashes of the previous burning phases. As an example, let us consider the pre—supernova evolution of a newly born star that is comprised almost entirely of hydrogen and a has a mass large enough to fuse all of the nuclei present in the core until fusion becomes endothermic. In this case, the star first burns through its 1H fuel by fusing it into 4He via one of the three branches of the proton-proton (p-p) chain once the core temperatures exceed ~ 4 x 106 K [15]. If there are any trace amounts of carbon, nitrogen, and oxygen (CNO) nuclei initially present, it is also possible fuse 1H into 4He, or an a particle, via one of the branches of the CNO bi-cycle at core temperatures greater than 1.5 x 107 K [15]. The threshold temperature for this reaction is higher because the Coulomb repulsion is greater and consequently requires higher thermal energies to penetrate. Once the core has consumed all of its 1H fuel, the star contracts and the core heats 5 up until its temperatures reach 108 K when it can burn its 4He by fusing into 12C via the triple-a process [15]. Once the 4He abundance becomes sufficiently low, the otherwise much less likely reaction 4He-I-12C—>160 also occurs [15]. For 4He burning. the higher temperatures are not only needed to overcome the Coulomb repulsion between the particles involved, but also to ensure that the rate at which a particles fuse together is large enough to offset the very short lifetime of the intermediate 4He+4He—>8Be nucleus. After the core has burned through all of its 4He supply, the star contracts again until the core reaches a temperature of 5 x 108 K when carbon burning, the fusion of 12C nuclei, can occur [15]. Oxygen burning, having a higher Coulomb barrier to contend with, only occurs at temperatures above 109 K [15]. The fusion of 12C and 160 nuclei is negligible since the intermediate temperature at which the 12C- 160 Coulomb barrier becomes penetrable, the 12C nuclei supply is quickly consumed through carbon burning [15]. Carbon and oxygen burning can yield a variety of final state heavy nuclei. For carbon burning, the possible final state heavy nuclei are 24Mg, 23Mg, 23Na, 20Ne, and 160 [15]. For oxygen burning, the possible final state heavy nuclei are 328, 318, 31F, 28Si, and 24Mg [15]. These reactions also produce free neutrons, protons, and (1 particles. These are all immediately captured by the heavy nuclei present since there is no Coulomb repulsion between neutrons and the heavy nuclei, and the temperature is sufficiently high at this point to render the Coulomb barrier between the protons and (1 particles and the heavy nuclei easily penetrable. Thus, in addition to the primary heavy nuclei created by carbon and oxygen burning, secondary reactions can generate a variety of isotopes with non-negligible abundances [15]. Following the consumption of the core’s 12C and 160 supply, the star’s contrac- tion resumes until core temperatures become high enough to burn the nuclei created by carbon and oxygen burning. Since 288i plays a dominant role in the reactions 6 that generate the nuclei present after this phase of hydrostatic nuclear burning is complete, this phase is referred to as silicon burning. The silicon burning phase is more complicated. Another process that disintegrates heavy nuclei becomes impor- tant. At temperatures above the oxygen burning threshold, but well below those required to allow the very large mutual Coulomb barrier between these heavier nuclei to have a non-negligible penetration probabilities, photodisintegration begins to occur [15]. Through this process, free baryons and light nuclei are emitted when energetic photons produced by fusion reactions disintegrate heavy nuclei. These free baryons and light nuclei are then immediately recaptured by the heavier nuclei present. These complementary processes tend to equilibrium, however the resultant state of nuclear statistical equilibrium (NSE) is not perfect. There is a leakage towards the stable iron group nuclei Fe, Co, and Ni generated largely by fusion reactions involving Si nuclei at temperatures above ~ 3 x 109 K [15]. These nuclei are resistant to photo- disintegration up to temperatures of approximately 7 x 109 K because they have the highest binding energy per nucleon. These nuclei are also the end of the exothermic fusion process. Therefore after the silicon burning phase, fusion no longer offers any resistance to gravitational collapse. The contraction of the star recommences. Now the core relies upon the pressure exerted by the electrons present in it to resist the inward crush of gravity. As the cores contraction ensues, the electron gas becomes increasingly degenerate. If the cores mass is not too large, the degenerate electron gas pressure can halt the gravitational contraction. The condition that must be satisfied for this to happen is where Moore is the mass of the core, 77 is its electron fraction, and MC’h is the so- called Chandrasekhar mass [15]. The Chandrasekhar mass limit is the maximum mass 7 of a self-gravitating sphere which can be supported by the pressure of a degenerate electron gas. Small finite temperature corrections can be added to the right hand side of condition (2.1) [16], but they are omitted here. 2.2 Review of Core Collapse and Bounce If the star being considered had an initial mass that was sufficiently large, at least ~ 8 to 11 MO [17], condition (2.1) can only be satisfied temporarily. Electron captures by heavy nuclei lower the electron fraction 7), and thereby the Chandrasekhar mass, and the core’s mass increases with time. The core’s mass increase is a result of the previously mentioned onion-skin-like stratification the star develops as it burns through its nuclear fuel supply. Each layer in which nuclear burning is still occurring loses mass in the process as the nuclear ash it creates becomes part of the mass of the layer below it. This is schematically depicted in figure 2.1 for the innermost two layers. In particular, figure 2.1 shows how the core’s mass increases as the silicon O burning —> Si burning —> Figure 2.1: Schematic depiction of the core mass increasing as the layers of nuclear burning above it continue to produce nuclear ash. burning layer surrounding it continues to produce the iron group nuclei that it is made of. Once condition (2.1) is no longer satisfied, the core begins to collapse. Two other major instabilities then develop. First the rate at which electrons are captured by heavy nuclei increases as the 8 core collapses and densities increase. The core is therefore deprived of its primary pressure source at an accelerating pace and consequently the collapse accelerates [15]. Secondly, since the electron gas is highly degenerate, it is largely insensitive to changes in the core’s temperature. Thus only drastic increases in temperature can significantly increase the electron gas pressure. Two factors prevent such increases in the core temperature. One is that electron captures are a thermal energy sink in the early stages of the collapse since the neutrinos created stream freely out of the core and carry away energy [17]. Second, at temperatures well below those that can significantly increase the electron gas pressure, photodisintegrations take place at a rapid rate [15]. While the temperatures are sufficiently high that nuclear interactions mediated by the strong and electromagnetic interaction are in equilibrium, those mediated by weak interactions are not [16]. This allows some photodisintegrations to result in irreversible losses of a large amounts of thermal energy and accelerate the rate at which electrons are captured [15]. In particular the reactions 56Fe + 7 —> 13 4He + 4 n (2.2) 4He+’7 —> 2p+2n (2.3) which take ~ 100 MeV and ~ 24 MeV of thermal energy to induce respectively can occur with a high frequency [15]. Some of the free protons produced by this reactions quickly capture electrons preventing the inverse fusion process from happening. Since the electron capture rate of free protons is larger than those of the nuclei present in the core at this stage [16], this source of free protons further accelerates the rate of electron capture and hence the collapse continues to accelerate as well. Thus it is the thermal energy lost through the electron capture induced NSE leakage described above and the rapidly increasing rate at which electron captures pull thermal energy out of the core that keep the temperature from drastically increasing and substantially raising the pressure exerted by the electron gas. With the electron gas rendered incapable of halting the collapse, it becomes nearly a free fall. At this point, the dynamics of the collapsing core naturally divide the core into two regions: an inner core and an outer core. The inner core collapses homologically and infall velocities are subsonic [16]. It’s outer boundary is the point at which infall velocities begin to exceed the local sound velocity [16]. Much of the outer core moves at supersonic velocities [16, 18]. Sound signals can only propagate inside the inner core. It is generally believed that the collapse continues until the nuclear matter is sufficiently dense that the neutrons become degenerate and generate a pressure strong enough to stop the collapse and repel the infalling matter [15]. This moment is referred to as bounce. Since the core is not pure neutronium, the density at which bounce occurs is super-nuclear. This super-nuclear density is found by most numerical simulations to be approximately three times normal nuclear matter density [16, 17]. During the collapse a large amount of gravitational energy is released. This energy AEgmv can be approximated for a core with mass MC = 1.5 MG) if we assume that the radius of the core at bounce is Rb z 20 km [15] and that its initial radius at the beginning of the collapse R,- >> Rb. In this case 2 2 1 1 GMC 46 AE mu m —GMC (— — —) z —— z 3 x 10 J (2.4) 9 Hi Rb Rb Coupling as little as ~ 1% of this energy to the rest of the infalling star would supply more than enough energy to eject the outer layers of the star and explain the observed supernova explosion energies of ~ 1 Foe = 1044 J [15, 16, 17]. We turn our attention now to the theories that attempt to explain this phenomenon. 10 2.2.1 Proposed Explosion Mechanisms All of the proposed explosion mechanisms discussed in this section share the generally accepted picture of how the core dynamics unfold immediately following bounce. After the stiffened nuclear matter pressure repels the infall, a mild pressure wave propagates through the inner core [16]. As this pressure wave propagates through the sonic boundary that divides the inner and outer cores, it turns into a shock. After this, the various proposals diverge and we now consider them individually. 2.2.2 The Prompt Shock Mechanism A tempting driving mechanism for supernova explosions is the so—called prompt shock mechanism. In this scenario, the shock simply blasts its way through the outer core and outer layers of the star and results in an explosion. Analytical arguments and numerical calculations initially suggested that the energy the shock has when it is formed is sufficiently large to accomplish this [19]. After extensive examination, this picture was deemed to be overly simplistic [21, 22, 23, 24, 25]. The shock loses a tremendous amount of energy as it propagates through the outer core. The temperatures and densities in the shock are high enough to significantly enhance the electron capture rates that occur in it and result in major losses in thermal energy. Additionally, photodisintegrations result in the loss of a large amount of thermal energy and supply a source of free protons in the way described in section 2.2. These free protons very rapidly capture electrons and result in further thermal energy losses. A successful prompt shock explosion requires the shock to propagate all the way through the outer core where the density profile drops off steeply and the hot silicon burning layer can resupply the shock with thermal energy [16, 17, 18]. Unfortunately calculations have found that this only occurs in the very lightest of stars that can end their lives with supernova explosions [16, 17]. The 11 outer cores of more massive stars are too large for shocks to propagate to their outer boundaries. At some point inside the outer core that is too large, the shock turns into an accretion shock in which additional infalling material accretes to the existing core and the outward motion has stopped. Thus the explosion fails. Therefore, if this picture is correct, it explains only a small fraction the observed supernova explosions and therefore cannot serve as a complete theory. Another driving mechanism for supernova explosions must exist. 2.2.3 The Neutrino Revival Mechanism Numerical studies conducted by Wilson et al. led to the discovery of the possibility that neutrino absorption by nucleons can revive a stalled shock [26, 27]. This model, in which a prompt shock is initiated, then stagnates, and is later revived by neutrinos is called the delayed shock model. A detailed description of how this process occurs is beyond the scope of this brief review and for this purpose we refer the reader to a useful summary given by Janka [19]. For our purposes here, it suffices to note that a non-negligible fraction of the neutrinos produced by electron capture and electron-positron pair annihilations in regions of the core with higher densities and temperatures can transport energy to the shock and that some of this energy is translated into kinetic energy of matter [16, 17]. Ultimately the shock resumes its outward path and results in a supernova explosion. This is still widely believed to be a correct model of a core collapse supernova [16, 17, 19, 20]. However this theory is not free of problems. Not only is this mechanism critically dependent upon the values used for neutrino capture cross sections and production rates as well as other quantities not known with great certainty, but the modeling of neutrino transport between the very short and infinitely long mean free path limits is quite problematic as well [17]. As a result of this, simulations have yielded mixed results for the delayed shock mechanism [16, 17, 28, 29, 30, 31, 32]. These results have been analyzed in 12 great detail [20]. Since we expect the supernova explosion driving mechanism to be quite robust, a more satisfactory explanation is still sought. 2.2.4 The Convection Revival Mechanism Convection is yet another way by which thermal energy can by supplied to a stalled shock to revive it. This mechanism is fundamentally different from the delayed shock mechanism since it heats the shock by mixing it with hot matter from deeper within the core rather than by direct neutrino capture. As in the previous section, we wish to avoid a lengthy rigorous discussion of the hydrodynamics of convection in supernovae and instead we simply note that in the wake of the shock, much of the matter is unsta- ble to convection [17]. In particular it is the region between the neutrino sphere, the hot dense part of the core in which neutrinos are in thermal equilibrium with matter, and the stalled shock where rapid convection quickly develops. Matter just above the neutrino sphere is heated by neutrino capture quite efficiently. It is then moved up the convection cell where it mixes with the cooler material in the stalled shock and warms it. Its place just above the neutrino sphere is taken by cooler matter that has sunk there from the cooler regions above and is receptive to heating by neutrino capture and the cycle repeats itself. In this way, more of the gravitational energy that has been converted to other forms and stored in the core as thermal excitations, electron and neutrino chemical potentials, etc. is made available to the explosion. Additionally, convection inside what will end up being the neutron star can further aid the explosion by intensifying neutrino luminosities. Numerical simulations that model convection, which is most naturally done in at least two spatial dimensions, seem to suggest that conveCtion plays a vital role in generating supernova explosions [20, 33, 34, 35, 36]. With this notion in mind, we note that the two-dimensional sim- ulations by Herant et al. that modeled convection yielded successful explosions while their one-dimensional simulations that did not model convection failed to explode 13 [35, 36], even though these simulations made use of the same microphysics. 2.2.5 The Acoustic Vibration Mechanism A very interesting supernova explosion driving mechanism that relies upon acoustic power generated in the inner core as the driver was recently proposed by Burrows et al. [37, 38]. Two-dimensional axisymmetric simulations over the full 180° of rotating and non-rotating progenitors of various masses found that the acoustic power gener- ated early on in the inner turbulent region stirred by the accretion plumes and most importantly the subsequent excitation and sonic damping of core g-mode oscillations to be the driving mechanism of supernova explosions [37]. The spherically asymetric sound pulses radiated from the core steepen into shock waves that are found to merge as they propagate into the outer mantle and deposit their energy and momentum with high efficiency [37]. One of the most appealing features of this model is that the acoustic power that drives the explosion does not diminish until accretion on to the core subsides. Therefore this power source is available as long as it may be needed to generate an explosion. All Newtonian simulations conducted by Burrows et al. have yielded successful explosions that, when top-bottom asymmetric, are self-collimating [38]. This model is quite encouraging. However this theory requires further investi- gation as other simulations capable of detecting core oscillations have failed to do so while they have yielded successful explosions driven by more traditional mechanisms since the acoustic driving mechanism was suggested [33]. We note that these results may not be mutually exclusive. It is possible that the acoustic driving mechanism is a fallback mechanism of sorts employed when the shock is not revived [38]. 14 2.3 Recent Numerical Studies Our review of the different techniques recently used to numerically model core collapse supernovae is presented in the following way. First the equation of state (E08) and weak interaction rate input needs of numerical calculations are discussed. Then the approaches used by selected groups are divided by the number of spatial dimensions they attempt to model and are briefly described. 2.3.1 Equation of State All simulations of supernovae require some sort of EOS that can characterize the ther- modynamic properties of matter in local thermodynamic equilibrium (LTE). Typical input needed by an EOS are the local temperature, mass density, and electron frac- tion. All other statistical mechanical quantities of interest, such as pressure or internal energy density, are calculated using the EOS. Intensive efforts are have been made to determine a realistic EOS for hot dense matter, in particular the nuclear contribution [39, 40, 41, 42, 43, 44]. The results of these efforts cannot generally be expressed in a simple analytic form. Instead these results must be tabulated, or a code numerically calculates them as needed. We note that some groups run semi-analytic codes that do use a simple formula for the E08 [19, 45, 46, 47, 48]. These formulae are not rigorously derived. They are designed to mimic the key features that a realistic EOS is expected to possess. Since this is not the direction we ultimately want to take our code in, none of the works discussed here make use of such simple approaches, and we discuss them no further. In this section we focus on the efforts made to numerically determine the electron-positron and nuclear gas EOS. A very popular EOS is the Lattimer and Swesty EOS. It has been used by many simulations recently [34, 36, 49, 28, 50, 51, 52, 53]. A formal discussion about the assumptions and techniques that Lattimer and Swesty used when deriving their gen- 15 eralized EOS for hot dense matter can be found in their 1991 publication [39]. We note here that Lattimer and Swesty worked out their EOS in the compressible liquid drop model based on the non-relativistic Skyrme Hartree-Fock framework. It models the presence of the electron-positron and nuclear gases and is applicable only in re- gions where NSE can be assumed. Therefore it is ideal for the progenitor core, but must be replaced by something else in lower density regions. Where the density is low enough, the nuclear contribution to the EOS is negligible and an EOS that models the presence of the electron—position gas is sufficient. There are several equations of state (EsOS) that do this, and their computational speeds and accuracies have been compared [54]. The Helmholtz EOS, developed by Timmes and Swesty, executes the fastest, displays perfect thermodynamic consistency, and has a maximum fractional error of 10’6 when compared to T immes’ very slow exact calculation. It is typically regarded as the EOS of choice when the nuclear contribu- tion to the EOS can be neglected. A methodical discussion about the assumptions and techniques that Timmes and Swesty used when deriving their electron-postiron EOS can be found in [54]. Here we state that it is a sophisticated EOS that, in addition to modeling the presence of a non-interacting electron-positron gas, models other phenomena such as the interaction between electrons and ionized nuclei and the pressure exerted by radiation. Furthermore, it is applicable over the enormous density and temperature ranges 10"9 kg/m3 to 1018 kg/m3 and 103 K to 1013 K respectively. The fact that the nuclear contribution to the EOS does not transition from dom- inant to irrelevant instantaneously means that there can be discontinuities in EOS knowledge in the transition between regions where the Lattimer Swesty EOS is ap- plicable and where the Helmholtz EOS is sufficient [40]. Some groups attempt to circumvent this difficulty with a patchwork approach in which multiple EsOS are used to span the gap between these two limits [34]. Other groups opt to use one E08 16 that can be used in all regions of the star [37]. One such EOS has been developed by Shen et al [40]. Shen et al. constructed their EOS of nuclear matter using the relativistic mean field theory for a wide range of densities, temperatures, and electron fractions. They treat the lepton gases as uniform and non-interacting. Their EOS is calculated specifically for use in supernova simulations and neutron star calculations. A rigorous description of their derivation of this EOS is beyond the scope of this discussion. A very concise description of their approach can be given in the following way: First they use relativistic mean field theory to construct the EOS of homogeneous nuclear matter. Then they use the Thomas-Fermi approximation to describe inhomogeneous matter, where heavy nuclei and free nucleons coexist. This EOS is preferred by some [37] over the Lattimer Swesty EOS since it is relativistic, more realistically models the presence of 0: particles, and again is continuously applicable everywhere it is needed. 2.3.2 Weak Interaction Rates All supernova simulations make use of some weak reaction rates. Since these rates govern the deleptonization and leptonization rates of matter, they singificantly im- pact the dynamics of core collapse supernovae. Consequently a tremendous amount of resources have been invested into calculating them theoretically [55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65]. The theoretical calculation of these rates can be quite involved; so here we simply mention that older calculations used the independent particle model and more recent calculations are made using shell model. These rates are tabulated over a range of temperatures, mass densities, and electron fractions relevant to super- novae for the purpose of interpolation. A slightly more in-depth consideration of the differences between the older and more modern rate calculations is given in section 3.12. We also note that efforts are currently underway to experimentally measure and tabulate electron capture rates relevant to a variety of astrophysical environ- 17 ments as well. Some of these efforts are being made at the National Superconducting Laboratory at Michigan State University (N SCL) [66]. Once FRIB is online, these efforts can be expanded closer to the proton and neutron drip lines. Ultimately an extensive tabulation of experimentally measured astrophysical rates will not only aid in our theoretical understanding of how nuclear structure effects these rates, but will be directly available to those that need them for numerical calculations. 2.3.3 One-Dimensional Simulations One-dimensional simulations that assume spherical symmetry are still useful. Even though most one-dimensional Simulations fail to produce explosions, much of the un- derlying microphysical processes that govern the dynamics of core collapse supernovae can only be modeled in great detail by using them. That is why some groups use data generated by one-dimensional simulations and map it onto their multidimensional simulations [17, 36]. Additionally, they can still be used to accurately model cer- tain isolated portions of the collapse and subsequent explosion [9]. One-dimensional simulations are also useful because their results can be compared to those produced by multidimensional simulations so that the importance of multidimensional effects can be studied [36]. The fact that realistic one-dimensional simulations tend to fail to explode suggests that multidimensional effects are of critical importance to the explosion mechanism. Liebendfirfer, Mezzacappa et a1. [28, 52] Liebendc’irfer, Mezzacappa et al. recently revisited one-dimensional studies of type 11 core collapse supernovae [28, 52]. There goal was to employ an advanced tech- nique to model the transport of neutrinos called multi-group flux-limited diffusion in their one-dimensional simulations in order to determine if previous failures with 18 one-dimensional simulations were the result of approximations made in this area. The multi-group flux-limited diffusion method is a sophisticated algorithm that realisti- cally models the transport of neutrinos of different flavors and energies [67]. More accurate knowledge of neutrino flavor and energy distributions is essential for accu- rately modeling neutrino-matter interactions, as all of these cross sections depend strongly on neutrino energy and many on neutrino flavor as well [68]. Liebendorfer and Mezzacappa’s simulations followed the collapse of a 13 MO star using both New- tonian [52] and general relativistic [28] gravity. They made no attempt to artificially model the effect of convection. Their simulations made use of the aforementioned Lattimer Swesty EOS and an adaptive grid. Their adaptive grid uses cells of vari- able volumes that can contract to resolve rapidly changing distributions and enlarge in region where distributions change slowly to save computation time and memory. None of their one-dimensional simulations resulted in an explosion. The fact that the inclusion of the multi-group flux-limited method did not result in successful ex- plosions lends further credence to the notion that multi—dimensional effects are key ingredients to the supernova explosion. Herant et al. [36] Another recent study of type 11 core collapse supernovae in one dimension was conducted by Herant et a1 [52]. They modeled the collapse of a 15 MO progenitor calculated by Woosley and Weaver [18] using a second order Runge-Kutta hydrody- namic grid based code with adaptive cell size. Their model uses Newtonian gravity, makes no attempt to artificially model the effects of convection, and uses a very sim- ple treatment of neutrino transport. By labeling each grid cell optically “thick” or “thin” based on the density they contain, neutrino transport is modeled in one of two ways in them. 19 In optically “thick” cells, neutrino transport is approximated by diffusion. Dif- fusion is modeled using the so-called flux-limited method which assumes that the neutrino number density n in a cell satisfies dn c .. .. -. ——-=min —Vn,V-DVn) 2.5 .. (.I I < > where D is a diffusion coefficient. This coefficient can place limitations on neutrino diffusion into cells with high neutrino number densities where neutrino degeneracy is non—negligible. Equation (2.5) clearly places an upper limit in the neutrino flux, hence the name flux-limited. In optically “thin” cells, the widely used central lightbulb approximation [69, 70, 71, 72] is employed. This approximation assumes that matter is bathed in a neutrino flux originating from the center of the star. The magnitude of the flux is essentially determined by the neutrino production rates in optically “thin” regions and the rate at which neutrinos can diffuse from optically “thick” regions to optically “thin” ones. The neutrino production rates used where calculated by Takahashi et al [55]. The one-dimensional calculations made by Herant et al. fail to yield explosions. 2.3.4 Two-Dimensional Simulations As previously stated, it is widely believed that convection plays an essential role in the revival of stalled shocks. To model convection naturally, at least two dimensions are required. This is the obvious appeal that two-dimensional simulations have over their one-dimentional counterparts. Presumably it is not merely a matter of pure coincidence that two-dimensional simulations tend to be more successful at yielding explosions than one-dimensional simulations. 20 Herant et al. [36] The initial conditions and algorithms used to model the gravitational force and microphysics for the two-dimensional calculations performed by Herant et al. are identical to their one-dimensional simulations discussed above. The only difference between the two models is the numerical technique and the number of dimensions used. The numerical technique they employ is called smooth particle hydrodynam- ics. This technique is used to model a variety of two-dimensional simulations of astrophysical environments and its defining feature is that it represents continuous mediums with a finite number of discrete particles [35, 36, 73, 74, 75]. While spherical symmetry is assumed for gravitation, cylindrical symmetry is assumed for all other purposes. This means that at most only half of the two-dimensional space needs to be considered. Most of their simulations do not even do that much and only use an opening angle of 90°. Preliminary comparisons of simulations they ran with various opening angles from 90° up to 180° showed little difference. The primary conclusion of Herant et al. is that convection above the neutrino sphere provides supernovae a robust and self-regulating explosion mechanism that is effective under a wide range of physical parameters. They compare the star after its shock has stalled to a powerful convective thermodynamic engine. The region above the neutrino sphere serves as the engine’s hot reservoir and the star’s envelope serves as its cold reservoir. Due to the large temperature difference between these reservoirs, the efficiency of this thermodynamic engine is high and leads to a self regulating ex- plosion mechanism. Once the envelope is heated enough to explode, the cold reservoir is gone and the engine stops. This may be why most measured supernova energies seem to fall in the 1 to 2 Fee range [15, 17]. They do concede however that convection is truly a three-dimensional phenomenon and that the inclusion of rotation into their model may alter their findings. 21 Fryer and Heger [36] Fryer and Heger have also recently modeled type II core collapse supernova in two dimensions [36]. Their model is very similar to the two-dimensional model of Herant et al. discussed above. The only significant differences are that they use spherically symmetric general relativity to model gravitation, a patchwork of ESOS for regions of different densities with the Lattimer Swesty EOS for the high density region, and their initial conditions are different. They use a rotating 15 MO progenitor with an equatorial rotation velocity of ~ 200 km/s that was calculated by Heger and Langer [76]. They enforce the conservation of angular momentum for each of the particles used in the smooth particle hydrodynamic scheme separately. They found that rotation had the effect of limiting convection overall and restrict- ing it to the polar regions. As a result of this, the convective region takes longer to overcome the accretion shock and the explosion occurs at later times and is weak- ened. These explosion are highly asymmetric. The mean velocity of matter in the polar region is roughly twice that of the equatorial region. These asymmetries are expected and may be possible explanations for the observed polarization of emitted light [77, 78, 79]. While some of their results are encouraging, this model is not a complete explanation. Uncertainties in the numerical implementation of such things as neutrino transport, neutrino-matter cross sections, the E308, and the spherical symmetric gravity limit their ability to make quantitative estimates and the specific results they present should be interpreted as “best-estimates”. The trends they ob- serve, like weaker explosions and lower black hole formation limits, are more secure. 22 Burrows et a1. [37, 38] The most recent two-dimensional calculations made by Burrows et al. have led to the possible discovery of the previously mentioned acoustic powered supernova explo- sion. Since their results were already summarized in section 2.2.5, here we just briefly review the numerical techniques they used. They made use of the VULCAN / 2D code developed by Livne [80]. VULCAN / 2D was developed with astrophysical application in mind and is an implicit method for compressible multidimensional flows. It consists of a purely Lagrangian step followed by an explicit remapping step. The remapping limits the time step size by the “particle crossing time” or accuracy considerations, whichever is more restrictive in a given instance. VULCAN / 2D is a Newtonian, two-dimensional, multi-group, multi-angle radiation hydrodynamics code that uses multi-group flux-limited diffusion. Its calculations are axially/azimuthally symmet- ric, extend over the full 180°, makes use of a Cartesian grid for the innermost ~ 20 km, and a spherical grid for regions further out. The fact that the code uses the full 1800 opening angel is an essential requirement for any simulation that hopes to follow vibrational modes of the core. Furthermore the central Cartesian grid allows the best zoning to be maintained in the core. This Cartesian grid can be made to follow the core as it moves as well. They employ the Shen EOS briefly described in section 2.3.1 and all of their neutrino-matter interaction physics is taken from the works Burrows and Thompson, and Burrows, Thompson, and Pinto [68, 81]. 2.3.5 Three-Dimensional Simulations Due to the enormous computational power required by three-dimensional supernova simulations and their tremendous complexity, few have been done. Fryer has been one of the most prominent figures on the forefront in this regard. He has used the smooth 23 particle hydrodynamics technique to model core collapse supernovae in three dimen- sions as well as collapsars and the merging of two white dwarfs [33, 82, 83, 84, 85]. In this section we review the most recent three-dimensional supernova simulation Fryer ran with Young. Fryer and Young [33] The progenitor used by Fryer and Young is a non—rotating 23 MO calculated by Young and Arnett [86]. It was calculated by the state of the art Tycho stellar evolution code. Instead of the classic mixing-length theory technique of modeling convection, Tycho uses a more realistic algorithm based on multidimensional studies of convection in the progenitor star. One of the most significant differences in its output is the lack of kinks in the density and temperature profiles seen in traditional calculations [18] that are artifacts of mixing-length theory which ignores hydrodynamic transport processes at and outside of the convective boundary. They use the SNSPH code [82] which is a parallel, three-dimensional, radiation hydrodynamics code implementing tree code Newtonian gravity, smooth particle hy- drodynamics, and flux-limited diffusion transport schemes. For comparison purposes with other groups, they make use of the Lattimer Swesty EOS down to densities of 1012 kg/m3. This intentional over usage of the Lattimer and Swesty EOS has the effect of significantly altering the entropy profile of the convective region in their model. Their simulation did not launch an explosion until over 600 milliseconds after bounce. This provided them with an ideal opportunity to study the evolution and structure of the convection below the accretion shock to late times. Ultimately con- vection did revive the shock and power the explosion. Convective down-flows buffeted the neutron star, giving it velocities in the 150 — 200 km/s range. Such neutron star 24 velocities are comparable to the low velocity simulations of Scheck et al. [87]. The E = 1 mode was not found to dominate the convection in their simulation and this precluded the neutron star from acquiring velocities > 450 km/s. It is however pos- sible that these larger velocities could be achieved with only minor modifications to the initial conditions and numerical setup [87, 88]. Additionally, their simulation did observe neutron star movement, but this movement did not develop into strong oscil- lations and become the energy source for the supernova explosion like Burrows et al. found. This too might be due to the lack of 8 = 1 mode dominance of the convection. Again we stress that these different findings are not necessarily contradictions. It is possible that numerical viscosity in the particles the smooth particle hydrodynamic approach uses to represent the inner regions of the star damp out these oscillations. As previously mentioned it might also be the case if the convective engine would not have worked, the acoustic driving mechanism would have developed. Further investi- gation into the different results of the Fryer and Young simulation and the Burrows et al. simulation is indeed warranted. 25 Chapter 3 The Test Particle Approach Despite the many advances hydrodynamic based calculations have made in the field of supernova science, a great deal more must be done before these simulations can be regarded as complete. A truly complete hydrodynamic simulation would in principle have to model the dynamics of multiple fluids with strongly time dependent viscosities to simulate the presence of a full ensemble of nuclei and neutrinos. Additionally, the sysem being modeled is three-dimensional, relativistic, has huge magnetic fields and length scales that vary drastically in time, and it must make use of relativistic radiation transport algorithms and Boltzmann transport algorithms for neutrinos. Most state-of-the-art hydrodynamic supernova simulations typically only track the abundances of free baryons, a particles, and an “average heavy nucleus”. A notable exception to this is the work of Hix [89]. Fhrthermore, state-of-the—art hydrodynamic simulations make simplifying assumptions about the flow of neutrinos, and few are done in three dimensions. Failing to model the propagation of a full ensemble of nuclei can possibly average away nuclear structure effects that may influence core dynamics. Between the limits of neutrino trapping and free streaming, many of the simplifying assumptions made by traditional hydrodynamic treatments of the flow of neutrinos can be problematic. Since these state-of-the-art hydrodynamic calculations already 26 strain the capabilities of high performance supercomputers, there may be a long wait before more complex models that simulate the presence of hundreds of different species of nuclei and propagate neutrinos in a general way in a three-dimensional system can be realized. This motivates us to move away from the traditional hydrodynamic approach that we are familiar with and draw from other disciplines of physics in an attempt to circumvent this technological roadblock. It turns out that the field nuclear collision modeling is an ideal candidate for this purpose. Numerical simulations of intermediate and high energy nuclear collisions must be able to model particle production, shock wave formation, collective deflection, as well as the interplay between regular and chaotic collective dynamics. Transport theories based on semi-classical implementation of kinetic theory [90, 91, 92, 93, 94, 95, 96] have been highly successful in meeting these requirements and in doing so reproducing experimental observables and pointing the way to new physical insight into these systems. The simulation of a core collapse supernova poses similar challenges [30, 34, 36, 52, 81, 97, 98]. It is therefore tempting to implement these types of kinetic theory based approaches for the physics and astrophysics of supernova explosions. This is the aim of our work. Our simulation of a core collapse supernova focuses on the stellar core. The motivation for doing this is that the genesis of the supernova explosion is believed to be in the core and we therefore concentrate our efforts on realistically following the core through the explosion phase in the hopes of using the unique advantages that our approach has to gain new insight into the mechanics of the collapse and subsequent explosion. In particular we hope that, through explicitly modeling the propagation of neutrinos and a full ensemble of nuclei, we can probe the dependence of weak reactions on certain nuclear structure effects that other supernova simulations are inherently incapable of and study the impact that these dependencies have on the collapse and explosion mechanisms. 27 3.1 Origins of the Test Particle Approach Efforts to advance the field of nuclear collision simulation beyond the hydrodynamic [99], mean field [100], and cascade [101] approaches required modeling phase space dy- namics by numerically solving a non-relativistic semi-classical transport equation with a two-body correlation source term. This equation, called the Boltzmann-Uehling- Uhlenbeck (BUU) equation, is a semi-classical approximation of the Wigner trans- formed time dependent Hartree-Fock equation that neglects all particle correlations higher than two-body correlations. It is given by [96] a _. _. fl _. _o _p d _. at -O -a—t- (T,p, t) + g ' VrffT,p, t) — VTU ' fo(r7pv t) ‘ 273.7] d ql’d ”d ‘12’ ' 72:2 ‘3'” 1 2 2 2 2 3—. ~ ~ ~ '5 (277141? + (12 — (111 " (121)) '5 (17+ ‘12 ’ 41’ — (12’) ' [fcffaé’ Iat)f(F1§2/at) (1 — f(779fi2t)) (1 _ f(Faq-.2?t)) — f(F'.132 we: «a t) (1 — f(F,ci‘1/,t)) (1 — f0", (72,. 0)] H where f is the Wigner transformation of the one-body reduced density matrix, U is the mean-field potential, and g is an isospin degeneracy factor. The left hand side of the BUU equation is the so-called Vlasov term that describes the temporal change of the phase space density f due to the interactions of the nucleons with the mean field. The right hand side of the BUU equation is the collision integral that represents the effects of the correlations due to the two-body collisions on the phase space density. The collision integral runs over all possible initial incident phase space element momentum vectors (72 and all possible final phase space element momentum vectors (71, and (72,. The delta functions conserve energy and momentum and the exit and entrance channels given by the final two terms respectively are weighted by the phase space availabilities of the states the phase space elements 28 are scattering into and out of respectively in accordance with the Pauli Exclusion Principle. Some initial attempts to numerically solve the BUU equation involved fully dis- critizing the relevant six-dimensional phase space and calculating the phase space densities in each grid cell in every time step. However even the coarsest of grids constructed for this purpose were prohibitively large containing ~ 109 lattice sites. To sidestep this limitation, it was proposed that, instead of tracking the value of the phase space densities in each cell, one could only follow the initially occupied phase space cells in time and represent them by imaginary particles, henceforth referred to as test particles. These imaginary test particles were then propagated in a way that modeled the physical evolution of the phase space. They interact with one an- other via mean field one-body potentials and scatter with realistic cross sections. For the simulation of nuclear collisions, using 100 to 1000 test particles per baryon was sufficient to accommodate the complexity of the phase space dynamics. The test particle method formally approximates the phase space density with a sum over delta functions [105] N: t) = 2: 63(F- 17,-(0)63 (5— are) (3.2) where N is the total number of test particles. The initial coordinates of these delta function point particles, or test particles, have to be determined by some physical input model. For simulations of nuclear collisions, a local Thomas-Fermi approxi- mation, properly Lorentz-boosted, is sufficient. Inserting this solution into the BUU equation generates the following simple first-order linear differential equations that govern the motion of the centroid coordinates of these test particles in the full six 29 dimensional phase space d .. -' _. -' a -' _. 32102 = —VUE05(7‘2') — VUC(qiari) +C(Pi) d _, _ p727 dtrz — mi (3.3) 2 = 1, .,N where again N is the total number of test particles, U E0 S is the mean field one-body nuclear EOS potential, U0 is the effective one-body Coulomb interaction potential between the ith test particle with charge q,- and the rest of the N — 1 test particles, and C(15)) symbolizes the effects that two-body collisions with other test particles have on the 2th test particle’s momentum. Solving the BUU equation with the test particle method has reproduced exper- imental data quite effectively [96, 102, 103]. Proceeding in a completely analogous fashion, the test particle approach can be used to generate simple semi-classical equa- tions of motion for the centroid coordinates of test particles used to model the phase space dynamics of systems with more complicated coupled transport equations. This has been done successfully for relativistic systems in which particle production is im- portant and coupled transport equations had to be simultaneously numerically solved [104]. 30 3.2 Modeling Macroscopic Systems With Test Par- ticles In the previous section, we discussed the origin of the test particle approach and described how it can be used to model a microscopic system. In that case, there were many more test particles than there were physical particles. If we wish to model a macroscopic system with a very large number of physical particles, a mole or possibly much more, it is clearly impossible to have the number of test particles, th, exceed the number of physical particles, N Computational limitations require that in phyS' situations such as this that Nphys/th >> 1. The test particle approach can still be applicable in these cases so long as th is sufficiently large to capture the gross dynamics of the macroscopic system’s phase space. The ratio Nphys NV”, effectively determines a scale cutoff of sorts below which details cannot be resolved. When Nphys/Mp becomes sufficiently large, some truly microscopic phenomena become impossible to directly simulate with test particles. Therefore it must be established that these unresolvable details do not impact the gross phase space dynamics and / or can be taken into account indirectly. This can be accomplished with convergence tests. These types of scale issues are certainly not unique to the test particle approach. Hydrodynamic calculations have to spatially discritize the systems being simulated and when the systems’ total volumes are large enough, the number and size of the cells the volumes are divided into can raise the same scale and resolution concerns. Representative particle models, which are very similar to test particle calculations in the large Nphys/Aftp limit, have to contend with scale and resolution issues as well when the importance sampling of the particles the system is comprised of is made. Failure to represent physical particles with certain characteristics with a sufficient number of representative particles can prevent calculations from resolving details 31 that are essential to the gross system dynamics. Having established that the test particle approach that has been used so success- fully in the study of nuclear collisions can be applied to macroscopic systems as well, we have opened the door to its application to the study of supernovae. This matter is discussed in the next section. 3.3 Modeling Supernovae with Test Particles For supernovae, the one-body transport equation for the baryon phase space density fb(:z:p) for the particular state b of the baryon is given by [106] 8fb(:1:p) Hi mi . Mg . at + Emmvffbfxpl - WVfquftlvjjfbffcpl + WVstngfbffFP) = 15mm + 13,,(mp) (3.4) Here II the phase space element momentum, E; (p) and M; are the effective en- ergy and mass of the baryon in the particular state b, Up(a:) and Us are the mean field nucleon vector and scalar potentials, and 15b and 151/ are the baryon-baryon and baryon-neutrino two-body correlation source terms that take into account how collision between baryons and baryons and neutrinos impact the baryon phase space density. The latter term is what couples the one-body transport equation for the baryon phase space density to that of the neutrino phase space density. The form and evaluation of 15b are discussed in our previous work [107]. 15V can be derived in analogy to work on coupled transport equations for heavy ion collisions [108]. The relativistic quantum nature of these source terms makes them more complicated than the BUU source term discussed in section 3.1, however their structure is the same. For any neutrino species, the transport equation simplifies to an equation of mo- 32 tion that contains only the streaming and baryon-neutrino collision terms since there are no mean field contributions and the effects of neutrino-neutrino collisions are neglected, at Eu(k) fy(a:k) = [gt/(13k) (3.5) It should be noted that employing this formalism implicitly enables us to incorporate matter oscillation of neutrinos into our model. These coupled transport equations are solved as discussed in section 3.1 by ap- proximating the phase space densities with sums over delta functions. Semi-classical equations of motion for the six-dimensional phase space centroid coordinates of test particles that represent baryonic matter and those that represent neutrinos are readily obtained and are tracked at all times. In this way, the dynamics of baryonic matter and neutrinos are treated in an identical fashion. This feature is very different from traditional hydrodynamic calculations. The initial conditions for baryonic or matter test particles are determined by the chosen progenitor. For neutrino test particles, the initial conditions are determined by the local kinematics at the site of their creation. The initial conditions and equa- tions of motion for both types of test particles are discussed in greater detail sections 3.4 and 3.7. Matter test particles must represent a very large number of nucleons. To ease the modeling of weak interactions, each neutrino test particle represents the same number of neutrinos as a matter test particle represents nuclei. This number is taken to be the core’s mass divided by the 56Fe nuclear mass divided by the number of matter test particles used to model the core. Beyond these similarities, matter and neutrino test particles are quite different and are hence described separately. 33 3.4 Matter Test Particles In addition to a tracking the location of matter test particles in six dimensional phase space, the locations of the matter test particles in multiple cubical grids is always known as well. These grids are used to model the scattering of nuclei and free baryons and calculate three dimensional statistical distributions and gravitational force fields. They are discussed in detail in section 3.9.2. Every matter test particle has its own temperature. Initially this is determined by the temperature distribution of the chosen progenitor. A matter test particle’s temperature can change in two ways. Local weak reactions can induce changes in the temperature of the matter represented by a matter test particle. So can exposure to other matter test particles representing matter at different temperatures. Both of these processes are exhaustively discussed in sections 3.11 and 3.17 respectively. Each matter test particle explicitly represents a fixed number of nuclei and can explicitly represent free baryons as well. These free baryons come in multiples of the number of nuclei a matter test particle represents. Initially the type of nuclei and number of free baryons a matter test particle represents is determined by the chemical composition of the chosen progenitor. At later times, these properties are determined by local electron capture rates and the number and type of interactions each matter test particle has with neutrino test particles. For simulations that include the fusion of free baryons and nuclei, the number of free baryons a matter test particle captures can also influence these quantities. These weak and strong processes are described in sections 3.11 and 3.18 respectively. There are a total of 385 nuclei with 2 S A S 60 that a matter test particle can represent. These nuclei extend from the proton drip lines to the neutron drip lines for all nuclei with 2 S A S 60 in the table of nuclear masses we use [109]. This table is displayed in figure 3.1. We note here that any table of nuclear masses can be used so long as we can determine the electron capture 34 a Included Stable Nuclei [j Included Unstable Nuclei Figure 3.1: Table of stable and unstable nuclei currently includcd in our simulation. rates for all of the nuclei it contains. This qualifier shall be. discussed further in section 3.11. The propagation of a full ensemble of nuclei is a significant advantage our approach has over traditional hydrodynamic calculations. Matter test particles implicitly represent electrons as well. \Ve assume that for each proton represented by a matter test particle. free or bound in a nucleus. there is an electron nearby so that macroscopic charge neutrality is preserved. This assump— tion renders the matter test particles charge neutral and advantageously permits us to avoid modeling Coulomb forces between them. Note that the charge neutrality of matter test particles is essential as the Coulomb interaction is N 40 orders of mag- nitude stronger than gravity. Electromagnetic interactions would therefore dominate 35 the the matter test particle dynamics in an unrealistic way. For the purposes of calculating gravitational forces and density distributions, all matter test particles are assigned the same average mass. This average mass is taken to be the mass of the core divided by the number of matter test particles used to model it. After the insertion of the delta function approximation of the baryon phase space density into equation (3.4), it is found that centroid of each matter test particle is subject to three forces: gravitation, a mean field nucleonic force, and a force exerted by the surrounding electron gas on the electrons it implicitly represents. The latter two forces are dependent upon local statistical distributions and are notationally lumped together and denoted by Floc- Matter test particles can also scatter with one another. The equations of motion for the centroid coordinates of the matter test particles are given by the following first-order differential equations 65.. -' -' .. ~_. a-tpj = FG,j+F10c(Tj)+C(Pj) d 173' —i~'- = (3.6) dt3 2 2 2 \/m +PJ/C j = 1,...,N th “ where F’ - is the ravitational force actin on the j matter test particle and C 13' -) Ga g g J symbolizes the effects that two-body collisions with other matter test particles have th matter test particle’s momentum. N is the number of matter test particles on the 3' used to model the core and is constant. These first-order differential equations are numerically solved using the time tested 4th order Runge—Kutta method [110]. For calculations run on a single processor, N cannot greatly exceed 10°, otherwise run times become prohibitively long [111]. This does impose some statistical limita- 36 tions on single processor simulations. These limitations shall be addresses in section 7.2 where conservation of energy is discussed. Computational errors scale almost uni- versally with l/x/N in our approach, so these limitations are expected to vanish in the large N limit. Values of N such as 108 or even larger are currently feasible on parallel computer clusters [112]. 3.5 Matter Test Particle Initial Conditions The evolution of stars on the main sequence is thought to be well understood [15, 113, 114]. Exhaustive efforts to numerically model main sequence stars have led to a multitude of models that successfully reproduce observed behavior [18, 76, 86]. Consequently there are many sets of calculations that could potentially serve as initial conditions for our model. For our preliminary single processor calculations, we chose to work with the spherically symmetric non-rotating 15 MG Woosely and Weaver progenitor [18]. This progenitor was chosen since it is very well known and has been used as the initial conditions for several other supernova simulations, making it particularly useful for comparison purposes [18, 115]. The core of this progenitor is 1.33 MO [18]. Its density and temperature distributions are depicted in figure 3.2. The matter test particles are initially spatially configured in such a way that the density distribution of our core matches the spherically symmetric Woosely and Weaver core. The specific way in which we calculate density and other distributions is discussed in section 3.9. Once the test particles have been been assigned their initial spatial coordinates, the determination of the initial test particle temperatures is straightforward. It is simply read from a fit of the temperature distribution given in figure 3.2. The initial momentum space coordinates and nuclear properties of the matter test particles are slightly more involved and are addressed separately below. 37 Loam) [g cm"’] 0 1:0 Interior Mass [MIMG] Figure 3.2: Plot of the initial density and temperature distributions of the progenitor as functions of enclosed solar mass taken calculated by Woosely and Weaver [18]. 3.5.1 Initial Matter Test Particle Momentum Space Coordi- nates The Woosley and Weaver progenitor is calculated up to the point at which the im- plosion begins. The outer edge of the core is just beginning to collapse at 1000 km/s [18]. This initial condition is suflicient to determine the initial momentum space co- ordinates of the matter test particles in simulations of non-rotating cores. To see this, consider the following. Recall that the dynamics of a collapsing core naturally divide the core into two regions: an inner core and an outer core. In the inner core, collapse velocities are proportional to radial distance 7' [18]. In the outer core, collapse velocities can be approximated as being proportional to 1/\/F [18]. The inner core can be taken to enclose 0.6 MG to 0.8 MG) [18]. Thus having knowledge of the initial collapse velocity at the outer boundary of the core and the dynamics of the collapse discussed above is enough to generate an initial radial collapse profile. For any simulations that include rotation, we simply add to the initial velocity 38 collapse velocity of each matter test particle the rotational velocity at its location determined by the selected initial angular velocity profile. Since the Woosley and Weaver progenitor is assumed to be spherically symmetric, we cannot choose an ini- tial angular velocity distribution that is too large. Otherwise the conservation of angular momentum would result in spherical asymmetry of the progenitor’s core. Thus we should restrict ourselves to more modest profiles like those shown in figure 3.3. The three angular velocity distributions depicted in figure 3.3 are the initial ,I.‘ ..... :n. v v I r r 1' r l r 1' 1 v I r 7 1, 10 L .................... . SN,“ 1 ........ 3N15B ’ - - -s~15c F"! 'T m 0.1 “U .3 0.01 (3 0.001 1 1 J 1 l a 1 1 L I 1 1 ‘ hitcrior Mass [MIMQ] Figure 3.3: Plot of the initial angular velocity distributions calculated by Heger and Langer for their 15 A19 rotating progenitor taken from a private communication from Chris Fryer. angular velocities calculated by Heger and Langer for their 15 MO rotating progen- itor [76, 116]. Our earliest works did discuss more rapidly rotating cores, but these calculations where geared towards probing the role of angular momentum in the early stages of the collapse [107, 112, 117]. These early calculations can only meaningfully follow the very beginning of the collapse since, beyond using initial conditions which indicated that some electron captures had occurred in the core’s center, weak reac- tions were not modeled and the electron gas was assumed to be perfectly degenerate at all times. 39 3.5.2 Initial Matter Test Particle Nuclear Properties The chemical composition of the Woosley and Weaver progenitor is shown in figure 3.4. The chemical composition of the core is rather simple. It is comprised of 54Fe, 56Fe, and “Fe” nuclei with 48 S, A S, 65 with neutron excesses greater than 56Fe. Due to a limited availability of electron capture rates, we restrict of considerations to nuclei with A g 60. We now have ways to extrapolate electron capture rates to nuclei with A > 60, however this technique was not developed until quite recently and all of the results presented in this thesis shall assume that A g 60. The number of 1.0‘ - r E “Fe” ' 56se- .1r .510 5 E : 3 E 10‘25' ‘3: 1 1° 0 1.0 Interior Mass [MIMQJ Figure 3.4: Initial chemical composition of the progenitor as a function of enclosed solar mass taken from Woosely and Weaver [18]. matter test particles that represent 54Fe, 56Fe, and “Fe” nuclei within a given radius 40 is read from a fit of figure 3.4. The particular species of nuclei represented by a matter test particle that represents “Fe” nuclei is selected by a Monte Carlo algorithm that samples the nuclei included in our simulation that satisfy the definition of “Fe” nuclei given above. 3.6 Neutrino Test Particles Our simulations explicitly model the propagation of neutrinos by representing them by test particles as well. In this way, we treat the dynamics of baryons and neutrinos on equal footing. This is another advantage that our approach has over traditional hydrodynamic calculations. Neutrino test particles are assumed to be massless, move at the speed of light, and subject to no mean-field-type forces. The only way they can interact with other test particles is through scattering with or being captured by matter test particles. Thus the propagation of neutrino test particles is quite simple. Merely multiplying the unit momentum vector of a neutrino test particle by the speed of light and propagation time determines its new location. No complicated numerical method of approximating the solutions to differential equations is required. This light speed propagation does put limits in on the time step size. To realistically model the propagation of neutrino test particles, the time step size should be no larger than 10-5 s. For the purpose of updating three dimensional statistical distributions that can be altered by weak reactions, the locations on the neutrino test particles in one of the cubical grids the matter test particles are tracked in is always known. This cubical distribution grid is described in detail in section 3.9.2. Unlike matter test particles, the number of neutrino test particles is not constant. Neutrino test particles can be created and destroyed. The latter process can be induced by a weak interaction or by the neutrino test particle escaping the core. The number of neutrino test particles therefore varies greatly at different times during the 41 simulation. The specific mechanisms by which neutrinos are created and destroyed that are included in our simulation and how we model them are discussed at length in section 3.11. For the remainder of this consideration of neutrino test particles, it suffices to note that the only neutrino production mechanism currently included in our model is electron capture by nuclei and free protons. 3.7 Neutrino Test Particle Initial Conditions We do not model the presence of the neutrinos created before the simulation begins. Some neutrinos created by the electron captures that lower the central electron frac- tion enough to begin the collapse will undoubtedly still be in the present in the core when our simulation starts. However the densities are low enough to make impact that neutrino-matter interactions have on the core dynamics at this stage of the col- lapse negligible. Thus we need not model the propagation of these neutrinos. We model only the presence of neutrinos created after our simulation begins with neutrino test particles. Since all neutrino test particles are created when the nuclei and/ or free protons represented by a matter test particle capture an electron, all neutrino test particles represent electron neutrinos. However neutrinos and anti-neutrinos of other species can be represented by test particles as well. For reasons discussed at length in section 3.11, only the presence of electron neutrinos is currently modeled. The initial spatial coordinates of a neutrino test particle are determined by the location where the electron captures that result in its production are modeled. The initial mo- mentum space coordinates are determined by the species of the nuclei that capture the electrons and the temperature and number density of the gas the electrons are captured from. A rigorous description of precisely how newly created neutrino test particles’ momentum vectors are determined is presented in chapter 5. 42 3.8 Data Management For the purposes of calculating statistical distributions, gravitational forces, modeling weak and strong reactions and two-body collisions, it is necessary to organize test particles according to their spatial locations. In addition to spatially grouping the test particles, it is also necessary to organize the matter test particles in these groupings by their nuclear properties so the effects of weak and strong reactions that occur in the them can be efficiently modeled. The spatial grouping of test particles is done in two ways throughout the simulation. They are addressed separately below. 3.8.1 Spherically Symmetric Spacial Test Particle Grouping To group test particles by their radial distances from the origin, spherical shells defined by matter test particle occupancy are used. To locate the inner and outer radii of these shells, a quicksort algorithm is used that sorts the matter test particles by their radial distance from the origin in a one-dimensional array. To minimize the vulnerability of distributions calculated with these spherical shells to statistical fluctuations, we require the number of matter test particles contained in each shell to be at least 104. For single processor simulations with 106 matter test particles, this limits us to 100 spherical shells in which statistical distributions can be determined. These spherical shells are redetermined after the matter test particles have been moves to their new locations in each time step. 3.8.2 Three-Dimensional Spacial Test Particle Grouping Spacially grouping test particles in three dimensions is more complicated. Some of the cubical grids used for this purpose contain 106 grid cells or more. Therefore systematically searching for the grid cell a given test particle is contained in would be very wasteful. To make efficient use of these grids, one must have a clever way of 43 accessing test particles that are located in a given grid cell. We address this issue by assigning each test particle a number called the grid index, [grid’ Ig’f‘id = 01223; +A/Iiy+iz (3.7) where 2'33, 2'3), and z' z are the 1:, y, and z indices of the grid cell containing the test particle, and 0 3 ix, 2y, 2'; < M, where M 3 is the total number of cells in the relevant grid. Employing a quicksort algorithm after each time step ensures that test particles with the same value of I grz' d are stored next to each other in a one-dimensional array. In this way, test particles inside a given grid cell can be found in a very efficient way, with an algorithm that scales as N log N. Usage of this algorithm only requires the recalculation of the 11:, y, and 2 grid cell indices of a test particle in any grid it must be tracked in after it is moved to a new location. 3.8.3 Grouping Matter Test Particles By Nuclear Properties Once it has been determined that a specific weak or strong test particle reaction is to be modeled in a volume, spherical shell or grid cell, we must be able to efficiently access matter test particles in that volume that represent the nuclei and free baryons involved in the reaction. To address the first of these needs, we proceed in a similar fashion to the way in which the three dimensional spatial grouping needs were met. We assigning each matter test particle a number called the nuclear composition index, [camp where A and Z are the atomic mass number and nuclear charge of the nuclei repre- sented by the matter test particle and Zmag; is the maximum charge of any nucleus included in our simulation. For simulations with Amax = 60, Zmag; = 30. A quick- 44 sort algorithm is employed in each volume in every time step to ensure that the matter test particles with the same value of [comp are stored next to each other in a one-dimensional array. To efficiently access matter test particles in a given volume that represent a cer- tain number of free protons or neutrons for each nucleus they represent, we assign each matter test particle numbers called the free proton and free neutron indices. These numbers are precisely equal to the number of free protons and free neutrons represented by each matter test particle for each nucleus they represent. Again we use a quicksort algorithm in each time step to ensure that the matter test particles with the same free proton and free neutron indices are stored next to each other in one-dimensional arrays and update the arrays every time a matter test particle’s free proton or neutron index is changed by a weak or strong reaction. Each time a matter test particle’s nuclear composition index and or one of its free baryon indices are changed during a weak or strong reaction, it is immediately moved to its new location in the effected index array or arrays for the spherical shell or grid cell containing it. This is done to ensure that weak reactions can still be efficiently modeled in this volume in the current time step. 3.9 Calculating Distributions After having explained the algorithms used to divide the core’s volume in such a way that spherically symmetric and three-dimensional distributions can be calculated in section 3.8, a description of how distributions are calculated in these volumes is warranted. There are four distributions that must be directly calculated with matter test particles to model the dynamics of the core’s collapse: density, electron fraction, temperature, and the average ,5 of matter. We turn our attention to the calculation of these distributions when spherical symmetry is assumed first. 45 3.9.1 Spherically Symmetric Distributions Let N be the number of matter test particles contained in each spherical shell. The density at Fk, the average radius of the kth spherical shell, with average volume Vic is taken to be = N - Mtp Vk (3.9) Pf'Fk) where Mtp is the mass assigned to each matter test particle and fit is the average radius of the kth spherical shell. The electron fraction at the average radius of the kth spherical shell is taken to be 7707,) - M (3.10) — Atotfk) where Ztotfk) is the total sum of the nuclear charges of the nuclei and free protons represented by matter test particles in the kth spherical shell and Atot(k) is the total sum of the atomic mass numbers of the nuclei and free baryons represented by matter test particles in the kth spherical shell. These quantities are readily calculated N ztotcc) = Z(Z+np).- '=1 (3.11) N Atotfk) = Z (A + "P + "Nb 221 km spherical shell. where the above sums run over the N matter test particles in the Z and A are the nuclear charge and atomic mass number of the nuclei represented by a matter test particle and n P and n N are the number of free protons and neutrons a matter test particle represents for every nucleus it represents. The temperature at 46 kth the average radius of the spherical shell is taken to be _ 1 N T(rk) = NET,- (3.12) i=1 where T,- is the temperature of the matter represented by the 2th matter test particle. In the spherically symmetric case, one can only speak sensibly of the radial component of the average 6 of matter fir. We take average fir of the matter at the average radius of the kth spherical shell to be given by — _1_ N . Z gs €17 ll MR 31' 'fiz' (3-13) 1 where B;- is the 6 vector of the it” matter test particle and ii,- is the unit vector h pointing to the centroid of the it matter test particle. These spherical distributions are calculated and stored in each time step. Radial Derivatives of Spherically Symmetric Distributions The radial derivative of any of the above distributions at the average radii of the kth spherical is not the spherical shells are approximated in the following way. If the innermost or outermost shell, the derivative of any distribution at its average radius is approximated by linearly interpolating the slopes of the lines connecting the value of the distribution at fk to its values at Fk+1 and fk—l' Consider some distribution Q in the kth spherical shell. Its radial derivative at 77k is taken to be 362319) 7-‘llr.+1 ‘ 77k X 62071;) - 62919—1) 57" 7Ik+1_ fie—1 Fig - Fk-l 1" -c i" — 7" + _k ’E‘1 x C“ 5+1) f2”) (3.14) 1"k+1 — Tic—1 7‘k+1 - Tie 47 The expression used to approximate the radial derivatives of spherically symmetric distributions at the average radii of the innermost and outermost spherical shells is less complicated since there is only one neighboring data point. It is simply taken to be the slope of the line connecting the value of the distribution at the average radius of the innermost or outermost spherical shell to its value at the average radius of its only neighboring shell. The radial derivatives of the density, electron fraction, and temperature distributions are calculated and stored in each time step. Linearly Interpolating Spherically Symmetric Distributions To calculate the value of a spherically symmetric distribution or its radial derivative at an arbitrary radius that is between the innermost and outermost radii at which it is known, we linearly interpolate it using the closest two data points. For all radii outside the outermost data point, we equate the value of the distribution or its radial derivative with its value at the outermost data point. For all radii inside the innermost data point, the value of a distribution is equated with the its value at the innermost data point. Due to the assumption of spherical symmetry, the radial derivative of all distributions must be zero at the origin. Thus for all radii less than the innermost data point, the radial derivative of a distribution is linearly interpolated between zero and its value at the innermost data point. Note that by the construction of the spherical shell that these inner and outer radii approximation will at most affect one percent of the matter test particles. 3.9.2 Three-Dimensional Distributions For multi-processor simulation, where memory and speed limitations do not confine our considerations to spherically symmetric matter distributions, in addition. to cal- culating spherical distributions, we use two types of cubical grids to calculate and interpolate three-dimensional distributions. Distribution grids, in which statistical 48 distributions are calculated, are static cubical grids centered about the origin that contain 101 x 101 x 101 grid cells of size 1 km3. The density, electron fraction, and temperature at the centers of these grid cells are calculated the same way equations (3.9), (3.10), (3.11), and (3.12) show for the spherically symmetric case. The only difference is that now N is not a fixed number, rather it is the number of matter test particles that happen to be in a given fixed grid cell during a particular time step. The calculation of the average 3 of matter in a distribution grid cell is the three-dimensional analog of equation (3.13) and as such does not involve matter test particle unit position vectors. Since there are more than 106 cells in the distribution grid, some cells will contain a very low number of matter test particles and others will contain none at all. There— fore unphysical fluctuations in the distributions can be caused by these occupancy problems. It is necessary to use this many grid cells not only to resolve fluctuations in the distributions, but also to allow super-nuclear densities to be measured. To circum- vent this pitfall, we impose some minimum value that N must be to directly measure the value of distributions in a given cell with the matter test particles it contains. If N is below this threshold in a given distribution cell, the value of distributions at the center of that cell are linearly interpolated using the spherical shells. Gradients calculated at the center of distribution grid cells that are not on one of the faces of the cubical distribution grid are approximated using a three-dimensional analog of equation (3.14). Since the data points in this case are equidistant, the coef- ficients will always be 1 / 2 allowing some further algebraic simplifications to be made to the formula that are trivial and shall not be commented on further here. Deriva— tives of a distribution in a grid cells on one of the faces of the cubical distribution grid with respect to the coordinate that runs perpendicular to the face are taken to be be given by the slope of the line connecting the value of the distribution at the center of that cell to its value at the center of its only neighboring cell in that coordinate. 49 Once the distribution grids have been calculated, the centers of its cells become the corners of an interpolation lattice with 100 x 100 x 100 1 km3 grid cells inside of which we can linearly interpolate three-dimensional distributions and their gradients. While the core is initially much larger than this cube, we intuitively expect that most deviations from spherical symmetry will occur, at least initially, near the center of the core when it is contracted enough to have good occupation in many of the distribution grid cells. Thus we focus our efforts on determining three-dimensional distributions in the region in which expect all of the phenomenon of interest to begin. It is important to note that the fact that the distribution grid must contain at least 106 cells precludes single processor calculations from making use of it. Since single processor calculations can only use approximately 106 matter test particles, it is clear that the distribution grid is far too fine to use. Thus single processor simulations all assume spherical distributions of matter and exclusively rely upon the spherical shells to calculate one-dimensional distributions. 3.10 Matter Test Particle Motion As explained in section 3.4, the motion of the centroid of each matter test particle is influenced by three forces and scattering with other matter test particles. The three forces are gravitation, a mean field nucleonic force, and a force exerted by the surrounding electron gas on the electrons a matter test particle implicitly represents. These way these forces and collisions are all modeled are discussed separately in the sections below. 3.10. 1 Gravitation Currently gravitation is modeled with Newtonian mechanics. We cannot however use exact Newtonian mechanics as this would be an N 2 algorithm. With N = 106, the run 50 time would be unacceptably long. To avoid these long run times, we have devised two ways to model gravitation. For simulations that assume spherical symmetry, we use a modified Newtonian monopole that makes assumptions about the central density distribution that remove the numerical singularity at the origin. For simulations that makes no assumptions about the distribution of matter, we have an algorithm that is three dimensional, does not violate causality, and is also free of numerical singularities. We discuss the spherically symmetric model of gravitation first. Spherically Symmetric Gravitation The Newtonian monopole approximation is an appealing alternative for simulations that assume a spherical matter distribution. Its application requires only knowledge of the locations of a matter test particles and their radial rankings, both of which are known quantities for all matter test particles. In this picture, the gravitational force acting on a matter test particle with position vector 7'" and radial rank N + 1 is simply given by Fa=—G-N-M(2;Fr3 (3.15) where G is the gravitational constant and MG is the mass assigned to each matter test particle for the purpose of calculating gravitational forces. The problem with this approach is that gravitational forces exerted on the innermost matter test particles can result in unphysical motion. Obviously the innermost test particle never feels a gravitational force since ra- dial rank 1 implies that the N in equation (3.15) is zero. Furthermore, and more importantly, when the central density is sufficiently high, many matter test parti- cles can be very close to the origin. All of these test particles except for the very innermost are subject to very large gravitational forces directed toward the origin. 51 Since our time step size is finite, matter test particles move finite distances between times at which the gravitational forces they feel are computed. Thus any one of these highly accelerated inner test particles may pass through the origin and move to a new location that is much farther away from the origin than the last location where the gravitational force it experiences was computed before it is recalculated. At this more distant location, the gravitational force may be much lower and result in the test particle artificially gaining kinetic energy. This violation of energy conservation can result in many matter test particles being ejected from the core in a given time step once the central densities become high enough. To circumvent these inner test particle problems, we assume the density is constant in the region containing the innermost 50 matter test particles. This assumption is compatible with the way the density distribution is calculated for the following reasons. For spherically symmetric simulations, the innermost point at which the density is known is at the average radius of the spherical shell containing innermost 104 matter test particles which is always larger than R50, the radius defined by the radial distance to the 50th innermost matter test particle. Thus changes in the density distribution cannot be measured for r 3 R50 in the spherically symmetric case. Recall that for three dimensional calculations, to directly calculate the density inside a cubical grid cell instead of linearly interpolating it off the spherical density shells, it must contain at least 250 matter test particles in order to avoid fluctuations due to low occupation. Therefore the sphere defined by R50 will always be much smaller than the central grid cell if it can be used to directly measure the density in it. Thus the density does not change much on the interval 0 S r S R50 in the three dimensional case either. From classical mechanics, it is known that inside a sphere of constant density that the gravitational field generates a linear restoring force [118]. Therefore the assumption of a constant density inside the sphere with radius R50 yields the following 52 simple result for the gravitational force acting in a matter test particle with position vector 7" and radial rank N + 1 F’ —G.N-Mg 77/73, N250 G = (3.16) —G- 50- Mg; F/Rgo, N < 50 The advantages that this approach has over the unmodified monopole model are as follows. Unlike the unmodified monopole model, this model can exert a gravitational force on the innermost matter test particle with radial rank 1. Furthermore, since the magnitude of the gravitational force decreases linearly to zero inside the sphere defined by R50, the violations of energy conservation caused by the radially asymmetric sampling of the gravitational field as test particles pass through the origin discussed above are negligible. As a result, no matter test particles are ejected from the core since the origin can no longer act as a numerical singularity. Taking only fractions of a second to run, this approach is ideal for spherically symmetric calculations where the minimization of run times can justify minor sacrifices in accuracy associated with violations of causality resultant from the implicit assumption that gravitational force act instantaneously in this model. Three-Dimensional Gravitation For fully three dimensional simulations of the collapse, an approach is needed that makes no assumptions about the distribution of matter, does not violate causality, and is also free of numerical singularities. To fulfill these requirements, we generate a gravitational acceleration lattice off of which we linearly interpolate the gravitational acceleration at arbitrary points in the following way. First we divide the volume the matter test particles occupy into identical cubical cells, the corners of which will form the lattice sites of what shall be referred to as the mass lattice. Then we “smear” the mass of each matter test particle over the eight closest lattice sites 53 defined by the corners of the cell it is in using what could be described as an inverse linear interpolation. The fraction of the mass of a given matter test particle that is “smeared” onto one of the eight closest cites to it is linearly proportional to its distance to it. Once the mass lattice is constructed, we exactly calculate the Newtonian gravitational forces that would be experienced by a matter test particle sitting at the center of each mass lattice cell exerted by the masses accumulated at all of its lattice sites. In doing so, we create a gravity lattice off of which the gravitational force felt by a matter test particle at any point inside it can be linearly interpolated. To avoid action-at-a—distance forces that violate causality, when we calculate the gravitational acceleration at a particular gravitational lattice site due to the mass accumulated at mass lattice site at a distance r, we should access the mass that was accumulated at that mass lattice site at the time r/c in the past, where c is the speed of light. Since we will generally not have a mass lattice that was calculated exactly at a time r/c in the past, we linearly interpolate between the two lattices that were calculated closest to that time. In this way we model the finite propagation speed of gravitational waves. Once the gravity lattice has been computed, we can determine the gravitational force acting on a matter test particle at arbitrary points inside it by linearly interpolation using the eight closest gravity lattice sites to the point of interest. In effect, this is a 4 dimensional space-time linear interpolation of the gravitational force. While the calculation of the gravitational force at each site in the lattice is an N 2 algorithm, in this case N is not the number of matter test particles used in the simulation, it is the number of lattice sites. Choosing the number of lattice sites to be approximately 104 allows us to calculate the gravitational forces acting on 106 matter test particles in two to three minutes on a single processor with minimal sacrifices in accuracy. The lattice sites can be chosen to be static or follow the core as it contracts but remain at safe distances from each other at all times, thus we do not encounter any numerical singularities. 54 Simulations that model weak reactions require approximately 105 time steps to follow the core through the explosion phase. Therefore on a single processor more than two weeks would be spent calculating gravitational forces using the model described above. It is clear that this algorithm must be implemented on a parallel computing cluster to keep the run times from becoming prohibitively long. 3. 10.2 Nucleonic Force The mean field nucleonic force acting on the nucleons represented by the matter test particles is generated by the mean field nucleon potential energy. Currently all of the mean field nucleon potential energies are momentum / temperature independent. This is not a result of necessity, but rather one of circumstance. The decision to make the simulation capable of following the core through the explosion phase using fi- nite temperature statistical mechanics was made relatively recently and consequently time restrictions prevented us from incorporating momentum / temperature dependent mean field nucleon potential tables. This would be a straightforward exercise. Ad- ditionally, all of the potentials used by our simulation are symmetric with respect to proton and neutron exchange. For comparison purposes, our code makes use of two mean field nucleonic potentials that neglect nuclear asymmetry and one that it takes into account. These potentials are discussed separately below. Symmetric Mean Field Nucleon Potentials Our simulations currently make use of the soft and stiff momentum—independer1t Bertsch-Kruse—DasGupta (BKD) isoscalar potentials [90]. The soft BKD potential is given by 7 6 U(p) = a1 + b (i) / (3.17) p0 P0 55 where a = —356 MeV, b = 303 MeV, and p0 = 0.16 frn_3 is the normal nuclear matter density [119]. The stiff BDK potential is given by p p 2 U(p) = a— + b (~—-) (3.18) where a = —l24 MeV and b = 70.5 MeV. Both of these potentials are plotted in figure 3.5 over the density range 0 S p/pO S 3. In figure 3.5, it is seen that both of the BKD 300 — Soft BKD / 250 rm -'- Stiff BKD / 200 / 150 / 100 / 50 / / MeV] nucleon '100 T l I l r 0.0 0.5 1.0 1.5 2.0 2.5 3.0 9,911 Figure 3.5: The soft and stiff BKD isoscalar potentials plotted over the density range 0 5 NM) S 3- potentials exhibit realistic features such as a minimum energy per nucleon at some density pmz'n near p0, an attraction for p < pmz’nv and a repulsion for p > pmz’n- For the soft BDK potential, it can easily be confirmed that pmin z p0. For the stiff BDK potential, the lower value pmin z 0.9 p0 is obtained and the repulsion for p > pmz’n is much stronger. 56 Asymmetric Mean Field Nucleon Potentials The mean field nucleon potential of isospin asymmetric nuclear matter can generally be expressed as a power series in the isospin asymmetry 5 = (pn — pp)/(pn + pp) [120]. In this notation, pn and pp are the neutron and proton densities respectively. Since we assume that the potential is symmetric with respect to proton and neutron exchange, there are no odd-order terms in this expansion [120]. Thus us the lowest order expansion in (5 is given by U(p, 6) = Uofp) + time) 62 + 0(64) (3.19) where U0(p) = U (p, 6 = 0) is the mean field nucleon potential in symmetric nuclear matter and Usymfpl = (3.20) 1 32mm 2 as? 520 where Usym(p) is the so—called symmetry energy. Many advances have been made recently in our understanding of this symmetry energy and its role in nuclear collisions and astrophysics [119, 120, 121, 122, 123, 124, 125]. Despite this, the symmetry energy is poorly understood in dense neutron-rich matter [1, 122, 126, 127]. Unfortunately this is precisely the environment in which we expect to see the mean field nucleon potential become dominant. Therefore the best we hope for the inclusion of a potential that has a symmetry term to yield is a “ballpark” estimate of the effects that isospin asymmetry has on core dynamics. Further adding to our woes is the fact that some of our calculations yield sufficiently low electron fractions in the central region of the core that the isospin asymmetry 6 the becomes large enough to require additional terms to be retained in the expansion made in equation (3.19), about which very little is known. As we shall see in chapter 57 8, this issue is primarily a concern only for simulations that make use of electron capture rates that are considered to be quite large. Given the uncertainty in the form the symmetry energy assumes and its possible insufficiency in capturing the potentials total dependence upon isospin asymmetry, the particular parameterizations of U0(p) and Usym(p) that we choose to use is arguably immaterial. So long as U0(p) looks something like the BKD potentials for symmetric nuclear matter, we are safe to use any of the parameterizations from any discipline of nuclear physics. For convenience sake, we chose to take a result from the study of nuclear collisions for which the following parameterizations of U0( p) and Usym(p) were generated [119] =a1 1, a modified Direct Simulation Monte Carlo technique randomly selects k matter test particle scattering pairs from the cell for elastic scattering. Once a matter test particle pair has been selected, we boost into their center of mass (C-O-M) frame and randomize their C-O—M frame momenta. This is schematically depicted in figure 3.7. All scatterings are modeled using relativistic El ti -_ as c p cm frame Figure 3.7: Two-body elastic scattering is modeled by randomly repositioning the C—O—M frame momentum vectors to opposite positions on the surface of the momen- tum sphere defined by their C-O-M frame momentum. kinematics. A similar approach was previously used in the simulation of heavy ion collisions [57]. For simulations that allow free baryons to be captured by nuclei, fusion is modeled during the matter test particle scattering process. The algorithm used to model the fusion of free baryons and nuclei is sufficiently complicated that it is explained separately in section 3.18. Weak reactions are not modeled during matter test particle collisions. They are modeled using different algorithms discussed 63 in section 3.11. 3.11 Neutrino Test Particle Production and Propagation The weak processes of critical importance to the collapse and post-bounce evolution of a core collapse supernova are [29, 59, 129, 130, 131, 132] p+e_ (A,Z)+e— u+e— V+N 1/+(A,Z) ll n+1/e (A,Z—l)+l/e u+e_ V+N u+(A,Z) p+176 (A,Z+1)+z7e N+N+V+z7 u+17 (A,Z)+z/+17 u+(A,Z)* ”WWW (3.26) (3.27) (3.28) (3.29) (3.30) (3.31) (3.32) (3.33) (3.34) (3.35) (3.36) (3.37) where a nucleus is symbolized by its atomic mass number A and nuclear charge Z, N represents a free proton or neutron, and U represents a neutrino or antineutrino of any flavor. Let us now review the roles these reactions play in the generally accepted picture of the collapse [1]. During the infall phase, the elastic scattering of neutrinos and nuclei (3.30) is 64 mainly responsible for the trapping as its channel abundance dominates those of the free baryons and its cross section are significantly larger than the neutrino capture by nuclei cross sections [68]. Shortly after trapping, the neutrinos are thermalized by energy downscattering, experienced mainly through collisions with electrons (3.28) whose phase space restrictions favor this process. Reactions (3.26) and (3.27) are also important during the infall phase and after, as they control the neutronization of the matter and thereby significantly influence the collapse dynamics. In the postbounce phase, when free baryon abundances have increased substan- tially, it is neutrino—nucleon scattering (3.29) that provides the main neutrino opacity. Lepton capture by nucelons (3.28) and (3.31) are now responsible for the dominant creation and absorption of electron flavor neutrinos. The dominant source of of u and 7' production is nucleon-nucleon bremsstrahlung [133, 134] (3.33) and V3175 annihi- lation [135, 136, 137] (3.37). Electron-positron annihilation (3.34) is a subdominant source for p and '7' neutrinos. Reactions (3.32), (3.35), and (3.36) have not yet been widely implemented in numerical models so the role that they play is still not clear. We certainly do not expect all of the nuclei to remain in their ground states through- out the collapse and explosion, so processes (3.35) and (3.36) must occur. Similarly if positrons are present, reaction (3.32) must occur as well. Since our simulation is an attempt to model the dynamics of a core during a core collapse supernova with a totally new approach, it is sensible to proceed incre— mentally. Therefore instead of trying to implement all of the above weak reaction simultaneously, we first incorporate those thought to be important during the infall phase and bounce. We do this with the intent of studying the results generated by calculations using the limited infall and bounce weak reaction network and proceeding with the integration of the rest of the above weak reactions once it has been deter- mined that these results are valid. This is the current status of the code. As such only weak reactions (3.26), (3.27), (3.28), (3.29), and (3.30) are have been included 65 so far. Our code does not yet model the presence of positrons, p or 7' neutrinos, anti-neutrinos, or excited nuclei. Rates for many of the processes involving positrons, p and 7' neutrinos, and anti-neutrinos exist [68], and it is clear how to include them in our model. Additionally, our model is uniquely poised to accurately simulate the presence of excited nuclei since it explicitly propagates an ensemble of nuclei. We look forward to investigating this topic further. Having established what weak reactions we currently model and our motivation for temporarily restricting our weak reaction network, we may now proceed with the descriptions of the specific algorithms used to model these weak reactions. Before delving into the individual considerations of these algorithms, we make note of a few characteristics that they all possess. All weak reaction algorithms make use of tabulated average particle energies and some make use of tabulated average neutrino— matter interactions cross sections as well. This was done for two main reasons. First, it circumvents the need to repetitiously randomly sample the thermal energy distribu— tions of nuclei, free baryons, and electrons during the simulation, that would otherwise slow the code down. Second, it also allows us to model the phase space restrictions of average final electron energies and average cross sections of neutrino-matter inter- actions involving a final state electron that is both faster and more accurate than dynamically testing each potential weak reaction that is sensitive to these restrictions for Pauli blocking. Additionally, all weak reaction are modeled in the frame in which LTE is assumed. The origin of this frame is defined to move with the average 3 of matter at the site of the reaction. Therefore it is necessary to boost into and/ or out of this frame when modeling weak reactions of any kind. 66 3.12 Neutrino Test Particle Production Since we do not explicitly simulate the presence of electrons, we use electron capture rate tables to model neutrino test particle production. Our input needs in this regard exceed those of other simulations since we need electron capture rates for all nuclei from the proton drip line to the neutrino drip line for all A S 60. As previously mentioned, most standard hydrodynamic supernova simulations only track the abun- dances of free baryon, possibly a particles, and an ”average heavy nucleus” [16, 115]. Since this “average heavy nucleus” is generally not a drip line nucleus and its A value does not drop down into the low single digits, most electron capture rate tables are confined to nuclei near the valley of )6 stability and, with the exception of free protons, usually have an A 2 20 [56]. Therefore, until more complete tables become readily available, we are forced to extrapolate the rates from whatever table we use to the drip lines and to low A values. Currently the source for electron capture rates use the Fuller-Fowler-Newman (FFN) table [56]. More recent calculations of weak reaction rates using new shell models of the distribution of Gamow—Teller strength have resulted in an improved and often reduced estimate of its strength compared to those the FFN calculations yielded using extrapolations of the known experimental rates and a simple single- state representation of this resonance [138]. Thus one could make the argument that we should not use the FFN rates [58]. However the rates calculated more recently using shell model that we might use instead of the FFN rates can be off by orders of magnitude compared to the experimentally measured values in some instances [66]. Thus given the uncertainty with the calculated rates, our need to extrapolate rates from any table we use to many nuclei, and the additional uncertainty associated with these extrapolated rates, no table can suit our needs perfectly. Therefore, we take the electron capture rates from the FFN table, extrapolate them to the nuclei we need 67 them at, and compare the results of calculations preformed with the tabulated F FN rates and/or the extrapolated FF N rates altered in various ways. In this way, we can see how the collapse and explosion mechanisms depends on the electron capture rates. The specific ways in which we change the capture rates are explained in section 5.1. For the remainder of this discussions, it suflices. to note that the electron capture rates for free protons and all 385 nuclei included in our simulation are tabulated over a range of densities, electron fractions, and temperatures relevant to] a core collapse supernova. 3.12.1 Testing For Neutrino Test Particle Creation To model neutrino test particle production in a given volume, spherical shell or grid cell, for each species of nucleus (A, Z) present in that volume, we do the following. First we interpolate the volumetric electron capture rate per nucleus 7‘(A, Z) of the nucleus being considered from the closest density, electron fraction, and temperature lattice sites at which it is tabulated. Then we calculate the number of nuclei (A,Z) present in the volume N(A, Z). Then we compute N(A,Z)-r(A,Z)-dt-V=m+n (3.38) where dt is the time step size and V is the spherical shell or grid cell volume, m is an integer, and n is a number in the set [0,1). We interpret n as the probability that m + 1 neutrino test particles are produced as the result of electron captures by the nucleus (A,Z) and 1 — n as the probability that m neutrino test particles are produced. A simple Monte Carlo algorithm decides between the two possibilities. By repeating this process for each species of nucleus present in each volume of the core, we model the capture of electrons by all of the nuclei present in our simulation. We note here in passing that for multi-proccesor simulations that model neutrino 68 in the three-dimensional distribution grid, we model neutrino production here first in all of the cells that contain enough matter test particles. Then when we model neutrino test particle production in the spherical shells that contain distribution grid cells in which this has already been modeled, we subtract the volume of these cells from the shell volumes. This way we avoid double considerations of neutrino test particle production in the distribution grid cells. 3.12.2 Creating a Neutrino Test Particle To create a neutrino test particle in a volume, spherical shell or grid cell, with electron number density n and temperature T produced when nuclei (A, Z) capture electrons, we proceed in the following way. First we interpolate EV(A, Z, n, T), the average LTE frame energy of the neutrino produced when an electron is captured by the nucleus (A, Z) from a gas with number density n and temperature T from table of such energies discussed in section 5.2.1. This energy times the number of neutrinos that are represented by a neutrino test particle determine the magnitude of the neutrino test particle’s LTE frame momentum vector. This LTE frame momentum vector is oriented randomly. Then the neutrino test particle is randomly placed in the volume. The average 6 of matter is interpolated at that location and the neutrino test particle’s momentum vector is boosted into the lab frame. This weak test particle reaction changes the local nuclear properties and temperature of matter. The way that these changes are modeled are discussed in sections 3.14 and 3.16 respectively. Test particle momentum is not conserved in individual neutrino test particle productions. The total test particle momentum however can be expected to be conserved relatively well since such a large number of neutrino test particles are produced with randomly oriented momentum vectors. Therefore the neutrino test particle momentum vectors will tend to add to zero. 69 3.13 Neutrino Test Particle Matter Interactions Since we explicitly model the propagation of neutrino test particles, unlike when we model the production of neutrinos, we do not need to defer to rate tables when we model neutrino-matter interactions. Instead we employ probabilistic algorithms that determine if and when the neutrinos represented by a neutrino test particle interact with matter. The general prescription for modeling the propagation of a neutrino test particle during a given time step is as follows. A beam attenuation argument is used to determine if the neutrino test particle is captured in the volume containing it. If it is, we model the capture of the neutrinos it represents. If it is not captured, we proceed with its propagation. To determine if the neutrino test particle elastically scatters with matter in the volume containing it, a beam attenuation argument is again used. If it does, we model the elastic scattering of the neutrinos the test particle represents. We test for neutrino-matter interactions between the neutrinos represented by the neutrino test particle during the given time step by repeating this process for every volume the neutrino passes through during the given time step. In the sections below, we individually describe how the aforementioned processes are modeled. However before immediately proceeding into the discussion of the algorithmic implementation of these processes, a brief discussion about the impact that neutrino oscillations have on the flow of neutrinos in the core is warranted. 3.13.1 Neutrino Oscillations in the Core The effect of coherent forward scattering must be taken into account when considering the oscillations of neutrinos traveling through matter [139]. In particular one must consider the interaction energy arising from W-exchange—induced electron neutrino forward scattering from ambient electrons when calculating transition probabilities [140]. If the mixing of electron neutrinos and tauon neutrinos is neglected, the prob- 70 ability of an electron neutrino with wavenumber k that has propagated a distance L through matter with electron number density 718 oscillating into the muon neutrino state is given by [139] P1_,2(k, L, 716) = sin2(2612) (%)2 sin2 (mi—Ln?) (3.39) where 1 = 149, 2 = u“, 012 is the so-called mixing angle for electron and muon neutrinos, L is the distance the neutrino has traveled since its creation, and [U(k) and lm(k,ne) are the characteristic neutrino oscillation lengths of a neutrino with wavenumber k in vacuum and matter with electron number density ne respectively. These characteristic oscillation lengths are given by [139] 4ah2k lva) [771% _ mglc2 (3.40) 2 -1/2 lm(k,ne) = 111(k) 1 + (lg/(5:3)) - 2COS(2912) (lg/(5:3)] (3.41) with 10(ne) defined by 100...) = 3;: (3.42) where G F is the Fermi coupling constant. We assume that formula (3.39) is applicable, since all of our neutrino are initially electron neutrinos and 613 is sufficiently small compared to 612 that the mixing of electron and tauon neutrinos can be neglected [141]. This simplifies our consid- erations significantly as it spares us from have to use the much more complicated three-neutrino oscillation formula [140]. Now let us consider the general behavior of the transition probability for different values of the ratio [U(k)/lo(ne). In the limit lv(k)/l0(ne) << 1, equation (3.41) gives lm(k,ne) % lv(k) and if we insert this result into equation (3.39) we recover the 71 vacuum oscillation transition probability [140] L P1__,2(k, L,ne) it: P1_,2(k, L) : sin2(2012) Sin2 ([7110) (3.43) U which has a transition probability amplitude sin2(2612) = 0.8597 [140, 141, 142]. In the limit lv(k)/l0(ne) >> 1, equation (3.41) gives lm(k,ne) z [0(ne) and the transition probability becomes P1_,2(k, L,ne) z sin2(2012) (’0("6))2sin2 ( 107(1)) (344) and therefore it is clear that by our assumption the transition probability is strongly suppressed. Since [U(k) or k and [0(ne) or 1/ne, as k and ne become sufliciently large, we expect the transition probability amplitude to be significantly reduced. This feature is of particular interest as it implies that high energy neutrinos pro- duced in regions of the core near or above nuclear matter density have much smaller probability transition amplitudes than lower energy neutrinos produced in regions of the core with lower densities. The transition probability amplitudes for intermediate values of lv(k) / 10(ne) are plotted in figure 3.8. To gain a quantitative sense of the extent to which transition probability amplitudes suppression effects the oscillations of neutrinos propagating through the core being modeled, we numerically evaluate formulae (3.40) through (3.42) for the case of 0.1, 1, 10, and 100 MeV neutrinos trav- eling through matter with electron number densities ranging from 1037 m—3 to 1044 —3 m , the minimum and maximum values respectively expected to be encountered in the core. To do this, we make use of [m% — mglc4 = 8 x 10"5 eV2 [140, 142]. The results are displayed in figure 3.9. It is clear from an examination of figure 3.9 that transition probability amplitudes are already strongly suppressed for low energy neutrinos ~ 0.1 MeV in the most diffuse 72 H O .0 co P as _o .p. Transition Probability Amplitude / P o l l 4 6 8 10 Iv(k)/ lame) —< O N Figure 3.8: Neutrino transition probability amplitudes for the P1__,2 transition plot- ted as a function of lv(k)/l0(ne) over the range 0.1 S lv(k)/lo(ne) S 10. regions of the core and for larger neutrino energies and electron number densities this suppression is enhanced. Therefore even though our model is can easily model neutrino oscillations, we do not, since we can safely neglect the effects of neutrino oscillations in the core. That is not to say that no neutrinos oscillations occur in the core since a 10_7 or even smaller chance of an oscillation occurring is not irrelevant in a system containing ~ 1057 neutrinos. However there is no need to represent any of these neutrinos with test particles. Their relative abundance is simply too low. Outside of the core, the transition probability amplitude can be much larger and any future simulations that attempt to model portions of the collapsing star outside of the core can easily model neutrino oscillations. We also note that if instead of the two-neutrino oscillation V5 S V“, we should consider quasi-two-neutrino oscillations 118 S 113;, where 113: is some admixture of V” and V7- that is mostly 11,”, is possibly 73 \ —o.1 MeV 8 1.5-09 \\ _1 MeV __ 3 .— E‘ 1.5-11 s 10 MeV __ a \\\ --- 100 MeV E 1.5-13 1.E'27 T r l 1 r I 37 38 39 40 41 42 43 44 |0910(ne) Im-B] Figure 3. 9: Plot of the neutrino transition probability amplitudes for the P1_,2 transition for 0.1, 1, 10, and 100 MeV neutrin4os 4propagating through matter with electron densities 1n the range 1037 m —3 to 1044 immaterial. This is the case since a formula just like (3.39) can be derived with a slightly different mixing angle 611. and mass energy square difference [771% — mg|c4 for such quasi-two-neutrino oscillations [140]. The same amplitude suppression argu- ments would apply to these hybrid state oscillations as well. 3.13.2 Testing for Neutrino-Matter Interactions For the purposes of calculating the capture and elastic scattering probabilities of the co—moving neutrinos represented by a neutrino test particle, simple beam attenuation arguments are used. The significant advantage that this approach has over tradi- tional hydrodynamic treatments of neutrino-matter interactions is that it is equally applicable in all regions of the core. No problematic assumption about the behavior 74 of neutrinos between the limits of trapping and free streaming are required. To see how these beam arguments are constructed, consider the following. Let a beam of co—moving identical neutrinos (a neutrino test particle) pass through a differential slab of matter with finite cross sectional area A and differential thickness dx centered about a point located a distance 2: from the beam’s starting point. Let 5.,- and ”i be an average neutrino-matter cross section and number density for particle 2' in this differential slab respectively. In this scenario, it is clear from elementary geometric probability that the probability P that one of the neutrinos in the beam interacts with the matter in the differential slab is given by P = Z 6,-(Lr)nz-(:r) A Adm (3.45) Now let N be the number of neutrinos in a “slice” of the beam before it enters the differential slab of matter. The number of neutrinos removed from the beam of identical neutrinos as it passes through the slab, those that have interacted with matter in it, dN is given by dN -N-P da: From this we conclude that a: :r’ N(:i:) = N0 exp (—/0 ACE—$0) (3.47) From this we can calculate the probability that a neutrino will interact with matter 75 after traveling a distance at through it Pings) by computing Pintfl‘) = [_VQLI'éV—(“Q :r dz’ = 1—exp (-—/0 WT) (3.48) It is the above formula that is interpreted as the probability that the neutrinos repre- sented by a neutrino test particle interact with the matter it encounters as it moves through a volume, spherical shell or grid cell. The integral is numerically evaluated using the sub-rectangle approximation. For single processor simulations that only make use of the spherical shells, it is permissible to use only one sub—rectange when approximating the integral in equation (3.48). This is so for the following reasons. In the early stages of the collapse, before the core has contracted much, the value of A is well approximated by a constant in a spherical shell. Thus one sub-rectangle is sufficient. At later stages in the collapse, in regions of the core the core is highly contracted, while A may very greatly over a radial interval of c - dt ~ 3 km, because the shells are defined by matter test particle occupancy, their radial thickness will be quite small. Therefore, the value of A can still be taken to be roughly constant in a given spherical shell. No such approximations can be made for multi-processor simulations that make use of the three-dimensional grid cells as they are fixed in space and have 1 km side lengths, over which A may vary significantly during the late stages of the collapse. Once this integral has been evaluated using the relevant average capture or elas- tic scattering cross sections described in sections 3.13.3 and 3.13.4, a simple Monte Carlo algorithm uses it to determine if the neutrinos represented by a neutrino test particle are captured or scatter elastically with matter in the volume containing it. The particle number densities 712' are easily calculated since we know the density, electron fraction, as well as how many matter test particles in a given volume repre- 76 sent what nuclei and free baryons. The average neutrino-matter cross sections 6.,- are interpolated from tables discussed in detail in chapter 5. 3.13.3 Modeling Neutrino Test Particle Captures Once it has been decided that the neutrinos represented by a neutrino test particle are captured by nuclei or free neutrons in the volume containing it, a capture channel must be selected. This is done by generating relative capture channel selection probabilities. These relative channel selection probabilities are calculated by weighting the average capture cross sections that were interpolated from a table for the purpose of evaluating the interaction probability (3.48) by the number of matter test particles in the volume that represent the nuclei and free neutrons that correspond to each capture channel. The probability of selecting a capture with index i is therefore given by N." . P(2') — 202 — —— 3.49 27' N151 ( ’ where here 6k and Nk are the average capture cross section and number of matter test particles in the volume that represent the nuclei or free neutrons that correspond to capture channel It respectively. Once these relative probabilities have been calculated, a Monte Carlo algorithm uses them to select a capture channel. This weak test particle reaction changes the local nuclear properties and temperature of matter. The way that these changes are modeled are discussed in sections 3.14 and 3.16 respectively. After the capture of the neutrinos represented by a neutrino test particle has been modeled, the memory allocated to it for the purpose of storing the information needed to model its propagation is de—allocated so that it can be used by another neutrino test particle created at a later time. As was the case with neutrino test particle production, test particle momentum is not conserved during individual neutrino test particle captures. However we do expect the average test particle momentum to be 77 conserved since many neutrino test particles with different momentum vectors will be captured this way so the net effect will tend to cancel. 3.13.4 Modeling Elastic Neutrino-Matter Interactions If it has been determined that the neutrinos represented by a neutrino test parti- cle elastically scatter with matter in the volume containing it, an elastic scattering channel must be selected. This is done in an identical fashion to the capture case considered in section 3.13.3. The probability of selecting an elastic scatting channel with index 2' is still given by formula (3.49). The only difference is that now, an electron channel is present and 5k and N k are the average elastic scattering cross section and relative elastic scatting channel abundances that correspond to elastic scattering channel It respectively. For free nuclear and free baryon channels, N k is still the number of matter test particles in the volume that represent the nuclei or free baryons with elastic scattering index k. For the electron channel, charge neutrality requires that N k = Zj Z, N j where N,- is the number of matter test particles rep- resenting free protons or nuclei with charge Z in the volume with elastic scattering index 1'. These 5k’s are interpolated from a table, discussed in section 5.2.4, that contain the average cross section of elastic neutrino scattering off free baryons and all the nuclei included in our simulation calculated at selected incident LTE frame neutrino energies and matter temperatures. As before a Monte Carlo algorithm uses these channel selection probabilities to select an elastic scattering channel. Once an elastic scattering channel has been selected, we model the elastic scat- tering of the neutrino test particle in the following way. Here our considerations are divided into two categories: elastic scattering off of nuclei and free baryons and elastic scattering off of electrons. For the former case, we first calculate temperature at the centroid coordinates of the neutrino test particle and then we calculate the energy of the co-moving neutrinos represented by the neutrino test particle in the LTE frame 78 defined by the average 5” of matter at the centroid coordinates of the neutrino test particle. Then we interpolate the average final state LTE frame energy of a neutrino with the given incident LTE frame energy that have elastically scattered by the se- lected nucleus of free baryon in a gas at the given temperature from a table of such energies. This table is rigorously described in section 5.2.4. This final LTE frame neutrino energy gives us the magnitude of the final LTE frame 3-momenutm vector. This vector is then randomly oriented in the LTE fame and then boosted back into the lab frame where it serves as the new momentum vector of the co—moving neutrinos represented by the elastically scattered neutrino test particle. If the selected elastic channel is the electron channel, we proceed in a similar fashion to the nuclear and free baryon case. The only difference is that since there is a final state electron, electron phase space restrictions come into play. Therefore the final state electron energies of neutrinos elastically scattered by electrons are tabulated over a range of electron gas temperatures, number densities, and incident LTE frame neutrino energies. So it is necessary to also compute the electron number density at the centroid coordinates of the neutrino test particle. Elastic neutrino-matter interactions do not change the local nuclear properties of matter. They can however alter the local temperature of matter. A detailed expla- nation of how changes in the temperature distribution induced by elastic neutrino- matter interactions are calculated for some of these interactions and in which cases they are modeled is given in section 3.16. Like the neutrino test particle production and capture processes, test particle momentum is not conserved during individual neutrino test particle elastic scatterings. However we do expect the average test par- ticle momentum to be conserved since many neutrino test particles will have their momentum vector randomly re-oreiented so the net effect will tend to cancel. 79 3.13.5 Neutrino Test Particle Movement Any neutrino test particle that survives the modeling of neutrino-matter interactions described above is moved to its new location determined by the number and type of elastic interactions it has taken part in during the given time step, or lack there of. In the free streaming case, the neutrino test particle’s movement during the given time step is simply determined by A5 = 73. c - dt, where 71 is the unit momentum vector it had at the beginning of the given time step and dt is again the time step size. If the neutrino test particle passes through multiple volumes, spherical shells and/ or grid cells, during a given time step and scatters elastically with matter in N of these volumes, it will have has N + 1 momentum vectors during the given time step. Consequently, its movement during the given time step is determined by N+1 A53: 2 fii C- 61:7; (3.50) i=1 h unit momentum vector the neutrino test particle had between where 71,; is the it its (N — 1)th and N M elastic scatterings and 6t,- is the time that elapsed between these scatterings. If the neutrino test particle’s position vector at the beginning of the given time step was i", then its final position vector at the end of the given time step is given by :73" = f + A27. Whenever it is the case that 3:, > Reore, where the core’s radius Reore is taken to be the outer radius of the outermost spherical shell, the neutrino test particle has escaped the core. Once a neutrino test particle has es- caped the core, the memory allocated to it for the purpose of storing the information needed to model its propagation is de-allocated so that it can be used by another neutrino test particle created at a later time. 80 3.14 Updating Matter Test Particle Nuclear Properties Every time a weak test particle reaction that changes the nuclear properties of matter is modeled in a volume, spherical shell or grid cell, we must do two things. First we must we locate a matter test particle with the appropriate nuclear properties and change them in a way that reflects the weak reaction that was modeled. Then we must rearrange the arrays that group the matter test particles in the volume by the nuclei and/ or free baryons that they represent described in section 3.8.3 so that matter test particles in this volume with specific nuclear properties can still be efficiently located when modeling additional weak reactions during the current time step. The latter task'is a straightforward one and shall not be commented on further. The former task can be complicated and requires further consideration. Updating the nuclear properties a of matter test particles that represent nuclei that have captured electrons or neutrinos is quite simple as long as the nuclei were not originally on the neutron or proton drip line respectively. All that is required in this case de—incrementing or incrementing the nuclear charge of the nuclei the matter test particle represents. Now consider the case of a matter test particle that represents nuclei sitting on the neutron drip line capturing electrons. In this case simply de—incrementing the nuclear charge of the nuclei would result in matter test particles that represent highly unstable nuclei existing only briefly over the neutron drip line. Instead we model the emission of the number of free neutrons required to make the final state nuclei one of the nuclei included in our mass table [109]. This is depicted schematically in figure 3.10. In figure 3.10, we see the neutron drip line nucleus 458 capturing an electron and emitting 3 neutrons to become the nucleus 42F, the nucleus with the largest A g 45 with Z = 15 included in our table of nuclei. It is through processes like this that matter test particle can acquire free neutron 81 — Neutron Drip Line Figure 3.10: Schematic depiction of the neutron drip line nucleus 458 capturing an electron and emitting 3 neutrons to become a 42F nucleus, the nucleus with the largest A S 45 with Z = 15 currently included in our table of nuclei. components and why they come in multiples of the number of nuclei represented by matter test particles. The capture of a neutrino test particle by nuclei represented by a matter test particle is handled in a completely analogous way when the nuclei sit on the proton drip line. In this case, free protons are emitted to obtain one of the nuclei included in our table and a matter test particle acquires one or more free proton components in the process that also come in multiples of the number of nuclei represented by matter test particles. The free baryon components of a matter test particle stream with the nuclei it represents. All matter represented by a matter test particle has the same temperature. 3.15 Updating Matter Test Particle Temperatures The ability to realistically model the changes in the temperature distribution is of critically importance for any simulation that follows the a core collapse supernova 82 through the explosion phase. Not only do the forces generated by the electron gas pressure and nuclear EOS generally depend upon the local temperature [113, 119, 121, 122], so do the weak reactions that can change the properties of the local electron and nuclear gases that determine the magnitude and direction of the forces they exert [56, 58, 62, 68]. Matter test particle’s temperature can be changed by two processes. One is if the electrons it implicitly represents are in a part of the electron gas where the tempera- ture is locally changed by weak interactions [58, 62]. The other possibility is that its temperature is changed due to exposure to nearby matter test particles that represent matter at different temperatures. Let us explore the former possibility first. 3.16 Weak Interaction Induced Temperature Changes Recall from section 3.11 that the weak processes currently modeled in our simulation are (A, Z) + e_ S (A, Z — 1) + V3 (3.51) 1/ + e_ S: u + e" (3.52) 1x+ (A,Z) ‘3 11+ (A,Z) (3.53) where a free baryon or nucleus is symbolized by its atomic mass number A and nuclear charge Z and 11 represents a electron or muon neutrino. Note that unlike in section 3.11 where separate notation was used to symbolize nuclei and free baryons, (A, Z) and N respectively, here it is much more convenient to notationally lump them together and symbolize them both as (A, Z). In principle, all of these interactions influence the matter temperature, however not all of the temperature changes they 83 induce are explicitly calculated. The processes depicted in equations (3.51) and (3.53) do in fact result in changes in the kinetic energy density of the local gas of (A, Z) and (A, Z —1) particles that can induce changes in its temperature. However, with minimal sacrifices in accuracy, the effects of the recoil of these particles can be neglected during these interactions. In this approximation, there is no change in the kinetic energy density of the gas of (A, Z) and (A, Z — 1) particles and therefore no change in its temperature either. Therefore we only explicitly calculate the changes in the temperature of the local electron gas induced during the processes given in equations (3.51) and (3.52) and neglect the changes in the matter temperature that result from elastic interaction given in equation (3.53). So for the purpose of calculating weak reaction induces temperature changes, we consider only the changes in the electron gas temperature that occur as a result of processes (3.51) and (3.52). Turing our attention to the former processes first, consider the following arguments. Let U V(n, T) be the kinetic energy density of an element of an electron gas with number density n and temperature T. Each time an electron is removed from or added to the gas element via electron or neutrino capture, UV is changed by some small amount 6UV. This process induces a small change in the number density 671 and, depending on the energy of the electron that was removed or added, a small change in the temperature 6T. Thus, to first-order, we may write 8U a 5Uv = (3%)T6n + (‘aU_f/‘)n5T (3.54) From this equation, it is clear that if one has knowledge of 6UV, 6n, and the partial derivatives (BUV /8n)T and (EUV /8T)n, the temperature changes 6T can be solved for. The relevancy of the above discussion to our test particle model can be elucidated in the following way. For the time being, let us assume that the partial derivatives of UV with respect 84 to n and T can be determined at all points in the core. An in depth description of the specific way in which the partial derivatives of UV with respect to n and T are calculated is presented in section 4 where we address the calculation and tabulation of statistical quantities of interest. Thus, to solve equation (3.54) for 6T, all that we must directly determine with test particles are the quantities 6UV and 6n. If a neutrino test particle is to be created or captured in a given volume, spherical shell or grid cell, the electron number density of the gas in that volume is changed by 6n = IFNV / V respectively, where NV is the number of neutrinos represented by a neutrino test particle and V is the spherical shell or grid cell volume. Thus the 671 in equation (3.54) is readily calculated. If a neutrino test particle is created in a spherical shell or grid cell with electron number density n and temperature T when the nuclei or free protons (A, Z) capture electrons, the electron kinetic energy density changes by 6UV = Eave(A, Z, n, T) 671., where Eave(A, Z, n, T) is the average kinetic energy of an electron captured by the nucleus or free baryon (A, Z) from an electron gas with number density n and temperature T and 671 < 0. This average energy is interpolated from a table that is discussed in detail in chapter 5 where the weak interaction tables used in our simulation are described. Similarly if a neutrino test particle representing neutrinos with energy EV is cap- tured in a spherical shell or grid cell, with electron number density n and temper- ature T, by nuclei or free neutrons, the electron kinetic energy density changes by 6UV = Eave(A,Z,n,T,Ez/) - 671, where Eave(A,Z,n,T,E1/) is the average kinetic energy of an electron produced when the nucleus or free neutron (A, Z), surrounded by an electron gas with number density n and temperature T, captures a neutrino with energy E; and 671 > 0. This average energy is interpolated from a table as well that is thoroughly described in chapter 5. Therefore, since 6UV is also readily calcu- lated when a neutrino test particle is created or captured during one of the processes given in equation (3.51) in a spherical shell or grid cell, equation (3.54) is easily used 85 to solve for the temperature change 6T in that volume induced by this test particle process. Now we address how the temperature changes that result from the elastic neutrino electron interactions depicted in equation (3.52) are calculated. If a neutrino passing through an electron gas element elastically scatters with one of the electrons in it, the kinetic energy density will generally change by some amount 6UV. Unlike the capture processes considered above, this process does not change 72.. Therefore 6 U V must correspond to some small change in the gas element temperature 6T. Again working in the first-order limit, we write 6U (SUV = (3%)” 6T (3.55) When a neutrino test particle representing neutrinos with energy EV elastically scat- ters with electrons in a volume, spherical shell or grid cell, with electron num- ber density n and temperature T, the electron kinetic energy density changes by 6UV = Eave(n,T,E1/) ~NV/V, where Eave(n,T,E1/) is the average change in the energy of an electron elastically scattered by a neutrino with an incident energy E, in an electron gas with the given 71 and T and NV and V are again the number of neutrinos represented by a neutrino test particle and the spherical shell or grid cell volume. This average electron energy change is also interpolated off of a table that is rigorously discussed in chapter 5. The ability to calculate 6 UV in the case of neutrino test particle electron elastic scattering in a spherical shell or grid cell allows equation (3.55) to be used to solve for the temperature change 6T in that volume induced by this weak test particle process. Thus all temperature changes due to weak reactions currently included in our model are accounted for. Each time a weak test particle reaction induces a tem- perature change in a spherical shell or grid cell, the temperature in that volume is updated. However the temperatures of the matter test particles contained in that 86 volume are only updated once in a given time step after all of the weak reactions have been modeled in it. It is assumed that the weak reaction induced temperature change of the electron gas in a given volume is spread homogeneously over all of the matter test particles contained in it. 3. 17 Heat Exhange Intuitively one would expect that if a “hot” matter test particle found itself sur- rounded by “cool” matter test particles, that the hotter test particle would warm up the cooler ones as the cooler ones cool it off. In other words, the subsystem ought to approach thermodynamic equilibrium in accordance with the zeroth law of ther- modynamics. Let the average matter test particle temperature in a given volume, spherical shell or grid cell, be Tave and let tmz'a: be the time it takes for the gas in this volume to mix and achieve a state of thermodynamic equilibrium. Assume that the 2th matter test particle is one of the matter test particles inside this volume h and has temperature Ti- To calculate the change in the it matter test particle’s temperature due to exposure to the other matter test particles in this volume during a time step At, we write At mix ATz' mia: = (Tave " T') ' , (3.56) where it is assumed that At / tmz'a: < 1. If this is somehow not the case, the ratio can be replaced by unity. Notice that this simple expression meets qualitative expecta- tions such as, relative to the local average temperature, comparatively hot matter test particles cool off, comparatively cool matter test particles heat up, and this heating or cooling towards equilibrium is directly proportional to the ratio of the exposure time At and the mixing time tmicr' To implement this approach into our test particle 87 model, we must calculate tmix- We take tmz’x to be given by average time it takes a matter test particle to “cross” the volume containing it, spherical shell or grid cell, as it moves along a straight line. In a spherical shell, “crossing” the volume means moving a distance equal its radial thickness Ar. Since the average radial matter test particle velocity 'Dr is'easily calculated in each spherical shell, in the spherical shells we may define _A: tmz’a: [177" (3.57) In a three dimensional grid cell, “crossing” the volume means moving a distance equal to its side length 6. Since the average matter test particle velocity vector 17am? is also easily calculated in each grid cell, in the grid cells we may define tmia: E I“ | (3.58) Thus tmz’x is readily determined everywhere in the core. Once per time step, after the core’s temperature distribution has been determined, equation (3.56) is used to update each matter test particle’s temperature in a way that models heat exchange with the other matter test particles contained in the same volume. In this way, mat- ter test particles retain a “memory” of their exposure to test particles representing matter at different temperatures. 3.18 Fusion Like our weak reaction network, our nuclear interaction network is very much a work in progress. Recall from section 2.1 that at temperatures above 7 x 109 K, interac- 88 tions between all of the nuclei present, including the highly stable iron group nuclei, that are intermediated by the strong and electromagnetic interactions are nearly in equilibrium. Thus a full nuclear interaction network would have to model fusion, photodisintegration, and possibly fragmentation simultaneously. We have developed some techniques to model these processes, however their simultaneous implementa— tion is not a trivial exercise. Instead, like the weak reaction network, we choose to proceed incrementally and first implement a simple fusion algorithm. Unlike the weak reaction network where we were able to first implemented the processes expected to be dominant in the infall and bounce phases, all of the aforementioned nuclear reac- tions are important at once. Thus no physically realistic data can be extracted from simulations that only model one or two of them. Therefore simulations that make use of the simple fusion algorithm described below can at present only be regarded as preliminary test calculations. For the purpose of conducting simple preliminary studies of the role that the fusion of free baryons and “light” nuclei (A < 60) plays in a core collapse supernovae, we have devised the following prescription. When we model the elastic scatting of matter test particle pairs as described in section 3.10.4, we also allow the matter test particles involved in each collision to capture some or all of the free baryons the other matter test particle might represent so long as the final state nuclei are among the 385 nuclei with A S 60 included in our simulation. For all nuclei (A, Z) such that (A + 1, Z) is a nucleus we propagate, the cross section of neutron capture is typically quite large (barns) [143]. Thus we always allow the nuclei represented by a matter test particle to capture free neutrons represented by another test particle that has been selected as its elastic scattering partner as long as the final state nuclei are among the 385 nuclei with A S 60 included in our simulation. Currently all that is done to model this process is the incrementing of the atomic mass number of the nuclei represented by the neutron capturing matter test particle and the de-incrementing of the free 89 neutron number of the other matter test particle. Due to the Coulomb repulsion between protons and nuclei, the free proton capture process is more complicated. In our model, Coulomb barrier penetration is allowed when the nuclear rest frame energy of the proton is greater than the so—called touching spheres potential given by 1.44 - Z Utouch 5 MeV (3.59) 0.8+ 1.2-A1/3 where we took the proton radius to be 0.8 fm [144] and the radius of the nucleus (A, Z) to be 1.2 - A1/3 fm [145]. To test for Coulomb barrier penetration in our test particle picture, the following is done. For this purpose, we assume that the matter represented by both matter test particles has the same temperature. This temperature may be taken to be the temperature at the center of the scattering grid cell containing the two test particles. With knowledge of the temperature of the gases of free protons and nuclei, their LTE thermal energy distributions can be determined, assuming that the gases are non-degenerate. However it is imperative that the relative motion of these test particles is taken into account when we test for Coulomb barrier penetration. We do this be defining the origins of LTE frames of the gases of free protons and nuclei to be the centroid coordinates of the matter test particles that represent them. Armed with the knowledge of the energy distributions of the free protons and nuclei in their respective LTE frames, one could sample these distributions, calculate the relative 6 of the matter test particles, and with a few Lorentz boosts measure the proton nuclear rest frame energy and test for Coulomb barrier penetration. However to avoid the repetitious sampling of the energy distributions and to obtain better statistics, we interpolate the probability P(A, Z, T, (3) that a proton will penetrate the Coulomb barrier of the nucleus (A, Z) in a gas at temperature T with the given relative LTE frame )8 from a table. This table is discussed thoroughly in section 5.3. 90 A simple Monte Carlo algorithm determines if the Coulomb barrier is penetrated. If it is, all that we currently do to model this process is the incrementing of the nuclear charge and atomic mass number of the nuclei represented by the proton capturing matter test particle and the de—incrementing of the free proton number of the other matter test particle. We note that this is a rather crude algorithm. In addition to not conserving energy, it does not allow matter test particles that represent nuclei with A = 60 to capture free baryons. Neither of these situations are indicative of inherent shortcomings with our model. We know precisely how to incorporate energy conservation into our fusion algorithm and allowing fusion to produce nuclei with A > 60 is also something that we know how to do, only time restrictions prevented us from accomplished these tasks before the results presented in this thesis were generated. Again we emphasis this algorithm must work in concert with photodisintegration and possibly fragmentation algorithms for the results to be realistic. This is discussed further in section 8. 91 Chapter 4 Electron Gas Statistical Mechanics Table Here we describe the process of calculating and tabulating the partial derivatives of the electron gas pressure and kinetic energy density with respect to electron number density and temperature needed for the reasons explained in sections 3.10.3 and 3.16 respectively. We also explain the calculation and tabulation of the electron gas 5 parameter, where 5 E p/kT. This quantity is needed to generate the weak reaction tables discussed in chapter 5. All of our considerations make use of the fully relativistic formalism. When doing this, we rigorously derive any expressions not taken from another source. We do not contend that all expression derived here are new results, however since most treatment of statistical mechanics are confined to the simple non-relativistic limit, a collection of all of the relativistic expressions we need to use are not readily found. Thus we elected to derive such a collection of relativistic expressions and present them as part of this thesis. All relativistic expressions derived are demonstrated to converge with their well known non-relativisitc limits. The statistical mechanical properties of the electron gas are tabulated logarithmi- cally over a range of electron gas number densities and temperatures relevant to the 92 collapsing iron core of our chosen progenitor [18]. We tabulate over a temperature range of 107 K to 1013 K. The electron number density over which we tabulate is determined by the minimum and maximum values obtained from densities in the set [10-10p0, 10 p0], where p0 is nuclear matter density, and electron fractions from the set [0.1, 1]. These are very safe ranges that we fully expect to exceed our needs. We divide our considerations into three categories: low temperature, arbitrary temperature, and high temperature. This is done because the low and high tempera- ture limits are much easier to work in, so we use them where ever they are applicable. We tabulate the exact arbitrary temperature expressions until they converge to within one percent of the high and low temperature limits. After this point, we tabulate the simple high or low temperature limit expressions. 4.1 Arbitrary Temperature Statistical Mechanics From elementary statistical mechanics [113, 146], it is known that exact expressions for many statistical properties of an ideal Fermi gas at an arbitrary temperature can be written as indefinite integrals of functions proportional to the Fermi-Dirac function given by 1 = exp(£ + E/kT) +1 f(E) (4.1) where T is the temperature of the gas, 5 is related to the chemical potential p via E —,u/kT, and E is the fermion kinetic energy. While the such integrals cannot generally be analytically evaluated for arbitrary temperatures, if 5 and T are known, they can be computed numerically [113, 146, 147]. The temperature of the non-degenerate matter and the density of the degenerate matter in all phases of a core collapse supernova are such that the momentum dis- 93 tribution of the electron gas cannot generally be assumed to be non-relativistic or ultra—relativistic [18, 113]. We must therefore work with the fully relativistic formal- ism when deriving the expressions used to determine the electron gas 5, pressure, and kinetic energy density distributions. Since most considerations of this and related topics found elsewhere are confined to the simpler non-relativistic limit, we rigor- ously derive all of these expressions and demonstrate their convergence with their well known non-relativistic limits. Since the pressure and kinetic energy density dis- tributions depend upon 5, we consider it first. 4.1.1 Determining 5 The 5 parameter of an ideal Fermi gas is known to be implicitly related to its number density n and temperature T through the condition [113] n _ 81/00 p2dp (4 2) _ 53 0 exp(€ + E(p)/kT) +1 ' where E (p) is the fermion kinetic energy expressed as a function of momentum. By numerically determining the value of 5 that best satisfies the above integral equation for a given number density and temperature, 5(n,T) can be determined for any n and T pair at which it is needed. To ease the numerical evaluation of the integral in equation (4.2), we re—express it in terms of kinetic energy E instead of momentum p. Working with the fully relativistic formalism, we write E00) 2 (1/1 + 01/ch — 1) me2 _, fl (p/mc)2 = M (4.3) m02 94 l —+ p2 = —2 (E(p)2 + 2mc2E(p)) (4.4) From which we conclude that (E2 + 2m62E) _1/2 (E + meg) dE (4.5) OIH dp = If we combine equations (4.4) and (4.5), we get pzdp = (E2 + 2mc2E)1/2 (E + 77162) dE (4.6) c3 After inserting the above into equation (4.2), we find that the number density is given by (4.7) 87r /00 (E2 + 2mc2E)1/2(E + m02)dE (he)3 0 exp(5 + E/kT) + 1 To further simplify the number density integral, we introduce the dimensionless vari- able u = E/kT kT 3 00 (112 + 2au)1/2(u + a)du —) n = 87’ (776) f0 exp(5 + u) + 1 (48) where a E ch/kT. 5(71, T) is determined by a program that systematically searches for the value of 5 for which the equality /oo (112 + 2au)1/2(u + a)du = n7r2 (3)3 (4.9) 0 exp(5 + 11) +1 [CT is most accurately satisfied for the given 71 and T. This process is fairly tedious and it would therefore be very inefficient to dynamically calculate 5 each time it is needed during the simulation of the collapse. Instead an external code is used 95 to calculate and tabulate 5 and any statistical quantity of interest that depends on it at selected electron number density and temperature pairs. These tables are then used by the supernova simulation as two dimensional lattices that are used to interpolate statistical quantities of interest at arbitrary electron number densities and temperatures. Armed with the ability to calculate 5(n, T), general expressions for the pressure and kinetic energy density of an electron gas with a given number density and temperature can be evaluated. We turn our attention to the first of these matters IlOW. 4.1.2 Pressure Integral From elementary statistical mechanics [113], it is known that the pressure of the electron gas is given by _ :81 0° 11311066119 P(n’T) _ 3113 f0 exp(€ + E(p)/kT) + 1 (4'10) where the electron velocity U(p) is most generally given by _ ii ”(’0’ ‘ me) = _1_ P m 1+(p/mc)2 . 87r 00 124dp(1+(p/mc)2)-1/2 7””) = 33.13/11 exp<§+E

/k1">+1 “'1” Again we re—express the integral in terms of relativistic kinetic energy E (p) Squaring both sides of equation (4.4) and multiplying the result by equation (4.5) gives us 3 2 p4dp = is (E2 + 2mc2E) / (E + mc2) dE (4.12) c 96 Substituting this result and equation (4.3) into the pressure integral yields 3 2 8,, 00 (E2 + 2mc2E) / dE P(n,T) - f0 — 3(hc)3 exp(5 + E/kT) + 1 After again introducing the dimensionless variable 11 = E / kT, we arrive at 3/2 3 00 112 + 2cm d'u P(n,T) = 810“”) kT/O ( ) 3 712 exp(5+u)+1 (4.13) (4.14) where a _=_ mcz/kT as before. As stated in section 3.10.3, for the purposes of cal- culating the average force acting on an electron in an element of an electron gas, we must tabulate (8P/8n)T and (BF/8T)” We address the evaluation of the partial derivatives of the pressure below. 4.1.3 Partial Derivatives of Pressure Equation (4.14) tells us that the electron gas pressure depends explicitly upon the temperature and implicitly upon both the electron number density and temperature through 5(71, T). Therefore we have 8P ‘ 95(5) (ff—")7’ — 85 371 T a: 5e) c) (57'),,‘ (95 BTn 8T6 (4.15) The partial derivatives of the pressure 8P/85 and (UP/Mk are readily obtained by first differentiating equation (4.14) and then carrying out the numerical integration. The partial derivatives of 5 however require further consideration. An examination of equation (4.8) used to numerically determine 5 for a given 71. 97 and T, tells us that n depends implicitly upon itself through 5. Therefore if we take the partial derivative of both sides of (4.8) with respect to n and constant T, we derive 1-512(5) _ 35 571T a a ‘1 +2 +2 871 T 85 where we made use of the obvious fact that (071/ 8n)T = 1. Since 871/85 is easily determined by first differentiating equation (4.8) and then numerically evaluating the resultant integral, (65/6n)T can be solved for and therefore so can (BF/870T. To calculate (35/3T)n, we proceed in an identical fashion and take the partial derivative of both sides of equation (4.8) with respect to T at constant 71. Doing so yields o—a—“c—se) ‘agarn 6T5 —1 _, (fi) = _ (€13) (£92) (4.17) 8T 71 65 8T 5 where in this case we made use of (871/ 8T)n = 0. Like (971/65, (8n/3T)€ is sim- ply calculated by first differentiating equation (4.8) and then numerically evaluating the resultant integral. This allows (85/6T),n to be solved for, which in turn allows (BF/8T)” to be calculated. 4.1.4 Kinetic Energy Density and its Partial Derivatives In addition to tabulating the partial derivatives of the electron gas pressure with respect to number density and temperature, the same partial derivatives of the kinetic energy density are tabulated as well. Recall from section 3.16 that these quantities 98 are needed to calculate the changes in the electron gas temperature induced by weak reactions. Calculating these derivatives requires us to first generate an expression for the kinetic energy density. To do this, we use the density of states approach [147]. In this picture, the number density n and kinetic energy density UV of the electron gas are given by _ 1 00 D(E)dE " ‘ v [0 exp(5 + E/kT) + 1 (4'18) _ 1 00 E - D(E)dE UV — V [0 exp(5 + E/k‘T) +1 (419) where E is still the relativistic kinetic energy of the electrons and D(E) is the so-called density of states. Comparing equations (4.7) and (4.18), we conclude that 87r (1..)3 (E2 + 2mc2E)1/2 (E + mc2) (4.20) 85 foo E (E2 + 2mc2E)1/2 (E + mc2) (IE 0 —> U = 4.21 V (hc)3 exp(5 + E/kT) + 1 ( ) Again after letting it = E / kT, we find that l 2 kT 3 oo 11 (112 + 2au) / (u + a) du UV = 8flkT —— / (4.22) he 0 exp(5 + u) + 1 with a E mc2/kT. It is the above forms of UV that is differentiated and then numerically evaluated and used to determine the partial derivatives (BU V / 8n)T and (3U V / 8T )7). needed to model the efl'ects that weak reactions involving electrons have on the temperature of the electron gas. 99 Proceeding just as before when we calculated the partial derivatives of the pres- sure, we note that equation (4.22) tells us that UV depends explicitly upon the temperature and implicitly upon both the electron number density and temperature through 5 (n, T). Therefore we have BUV _ BUV 55 . (811 )T — 35(55)T (4.23) BUV _ _UEZ _5_5 aUv) A (aTln ’ 1% (wlfifar g “'2‘” The partial derivatives of the kinetic energy density BUV/85 and (OUV/BT)€ are readily obtained by first differentiating equation (4.22) and then carrying out the numerical integration. The partial derivatives of 5 are taken from equations (4.16) and (4.17) to yield. aUV _ Buy an ‘1 ((9—an — afar) “'25) —1 (WV = _WV (2%) (3’3) + (2211) (4.26) 3T n (35 35 8T 6 3T 6 The inclusion of the partial derivatives of (BUV/Bn)T and (6UV/8T)n to our list of tabulated statistical quantities completes the table. 4.2 Low Temperature Statistical Mechanics When an ideal gas of fermions is said to be in the low temperature limit, its chemical potential is positive and sufficiently large that p/kT >> 1. In this case, the Fermi- Dirac distribution approaches a step function centered about the local Fermi energy and its sensitivity to the temperature is greatly reduced. Since exact expressions for 100 many statistical properties of an ideal Fermi gas can be written as indefinite integrals of functions proportional to the Fermi—Dirac function, this limiting behavior makes it is possible to derive analytic expressions for the chemical potential, pressure, and kinetic energy density of an electron gas to any order of accuracy in temperature corrections [113, 146, 147]. Here we generate these expressions with the lowest order temperature corrections using fully relativistic kinematics. Since discussions of the Sommerfield expansion technique [113, 146, 147] usually employed to derive statistical mechanics in the low temperature limit are typically confined to the simpler non-relativistic limit, we rig- orously derive this treatment of indefinite integrals of a specific class of functions in a completely general way and then apply it to the fully relativistic integral expressions for the statistical quantities of interest. Once these expressions have been obtained, we demonstrate their convergence with their well known non-relativistic limits. Recall that the Fermi energy in an increasing function of the number density of particles in the fermion gas [146]. Thus when the number density of a degenerate Fermi gas becomes sufficiently large, particles near the Fermi surface can become relativistic. The conditions in the core are such that the electron gas present in it be both degenerate and relativistic, so it necessary for our treatment of degenerate electron gases to be fully relativistic as well. 4.2.1 General Sommerfield Expansions Let us consider integrals of the form I = [00° g(E)f(E)dE (4.27) 101 "-5... where g(E) is any smooth function of fermion kinetic energy E and f (E) is the Fermi-Dirac function given by 1 f(E) = m (4.28) where [3 E l/kT and T and p are the fermion gas temperature and chemical potential respectively. Integrating (4.27) by part, we define u = f(E) —4 du = %dE (4.29) do = g(E)dE ——» v = C(E) (4.30) where G’ (E) = g(E) and conclude that I = E5300 G(E)f(E) — G(0)f(0) — [000 mag—£1113 (4.31) If we restrict our considerations to functions g(E) such that G (E) f (E) —> 0 as E —> co, the above expression for I becomes _ 43(0) _ °° BL I—1+e_fl“ /0 G(E)aEdE (4.32) If we assume that fly > O, the first term in the above can be expanded into a geometric series in —e-fl“. Doing so yields 00 _. n 00 8f I=—G(0) 1+ —-e 31‘ — G(E),, dE (4.33) z ( ) f, If we further restrict our considerations to cases in which 1311 >> 1 in our approximate treatment of integrals I , we may neglect the additive contributions to terms of the 102 order of unity or larger from terms of the order e_fifm with n 2 1. _. :r = —G(0) + [000 C(E) (€45) dE (4.34) Additionally, the assumption that 611 >> 1 implies that —8 f / 8E is narrowly peaked about E = ,u [113, 146, 147]. This opens the door to approximating C(E) with a Taylor expansion about E = p. _Carrying out this expansion yields 0° 1 am 00 8f) I=—00 — — E— n —— dE 4.35 (’+,§,n1(aE")E=,./o ( 10 ( 6E ( 1 Now let us turn our attention to the evaluation of —0 f / (9E. From the definition of the Fermi—Dirac function in equation (4.28), it is clear that __8_f_ = 7 65(E-fll 8E ’ (exam—10,1)? (4.30) Substituting this result into the above yields 00 n 00 813(E— ) I = —G(0) + 71::0% (3%) E=#/0 (E _ Hlnfi (BME—H): 1)2dE (4.37) If we introduce a dimensionless variable t E )3( E — p), the above becomes 00 n 00 n t I=—G(0)+Z-1-71— LG / ——t8 dt _ n14" 6E" E:.. 4101.11)? 00 1 1 BnG t — t [/°° ___._.__ _ j W _4] (.3... —00 (et + 1) —00 (et + 1) 103 Note that the second integral in each term of the above sum is of the order of e—fi” and can they can therefore be neglected in our approximate treatment of integrals I if it can be demonstrated that their inclusion is equivalent to adding terms terms of the order e—fi“ to terms of the order of unity or larger. To prove that this is indeed the case, consider the following. Since et/(et + 1)2 is even, 00 tnet In E / ————2dt 75 0 V even 72. (4.39) —00 (et + 1) by parity. However, since it is easy to demonstrate that 10 = 1 and it is clear that In is a strictly increasing function of even n’s, it must be that In is of the order of unity or larger for all even 72 Z 0. Therefore the inclusion of the second integral in each term of the sum in equation (4.38) is equivalent to adding terms of the order 6’5” to terms of the of unity or larger as for each In of the order of unity or larger there are at most two terms of the order e‘fi“ added to it. Therefore without further sacrificing accuracy in our approximate treatment of integrals 1', we may write 001 1 6W6 ‘” 1 1 6W0 = G01) — C(O) + Z FEE (w) E_ In n=2 _ 00 ll \ 91 .9 9.. tr: + EIH ml” 3 /"\ as :63 \__/ ta 5“ MAM Where in the second line in the above we made use of the fact that 10 = 1 and 11 = 0 and in the third line we used the fact that G(E) is by definition the anti-derivative 0fg(13) 104 We know that the more degenerate a gas of fermions becomes, the more closely the Fermi-Dirac distributions approaches a step function centered about the local Fermi energy 6 F [113, 146, 147]. This leads us to expect that when a fermion gas is sufficiently highly degenerate that (3;; >> 1 and u —-> 617, that for any smooth function g(E) that 00 617 p. /0 g(E)f(E)dE—+/0 g(E)dE~/0 g(E)dE (4.41) Thus the contributions made by the terms in the sum in equation (4.40) ought to be small corrections. Therefore we sensibly seek to derive an expression for our approximate treatment of integrals I with the lowest order correction in temperature. Thus we equivalently retain only the term with the lowest order in fi’l in equation (4.40) and in doing so derive _ fl 1 1 89 where we again made use of the fact that G" (E) = g(E). We can further refine the above by noting that the assumption that fly >> 1 —> )u = eF(1 -6) where 0 < 6 << 1. Therefore the integral in the above may be written in the following way. p EF 6F [0 g(E)dE— f0 g(E)dE— f“ g(E)dE (4.43) To lowest order in 6, the second integral in the right hand side of the above may be approximated by g(e F)(e F — ,u) giving M 6F / g(E)dE a: /0 g(E)dE — g(eF) (6F — p) (4.44) 0 Inserting this and the readily confirmed result [2 = 7r2 / 3 into equation (4.42) and 105 expressing fl in terms of T yields 6 7r2 fooo g(EmEldE “/0 F g(E)dE - g(CF) (6F — )u) + g(k-Tflg’m) (4.45) Approximation (4.45) is the general Sommerfield expansion for a smooth function g(E). It is applicable for any smooth function g(E) with anti-derivative G (E) that satisfies G(E) f (E) —> 0 as E —> oo whenever it can be assumed that )3/1 >> 1. 4.2.2 Low Temperature Chemical Potential To derive a fully relativistic expression for the electron gas chemical potential in the low temperature limit, it is necessary to consider the number density integral given in equation (4.18) by 1 (X) n— [0 D(E)f(E)dE (4.46) _ V where D(E) is again the fully relativistic expression for the so—called density of states previously shown in equation (4.20) to satisfy as v = C (E2 + 2mc2E)1/2 (E + mc2) (4.47) where C = 87r/ (hc)3. D(E ) / V can easily be demonstrated to satisfy the requirements that g(E) did when approximation (4.45) was derived. Therefore if we assume that flit >> 1, we are free to use approximation (4.45) when evaluating the number density integral simply by replacing g(E) by D(E) / V in the formula. Doing so gives 6 D e 7r2 I n a: A; F ng — $2 (61: — p) + -—6—(kT)22‘-£—'ul (4.48) 106 By the definition of D(E), the first integral is 72.. Therefore we arrive at 2 I N W 21904) " " CIT—Fair) D(eF) 7r2 D'(,u) kT 2 = .F 1‘F6F0(ep)(¥) (4'49) Furthermore, if ,u = eF(1 — (5) with 0 < 6 << 1, any attempt to expand D, (E) about 61: and retain anything beyond the zeroth order term would introduce terms at least second order in 6. Therefore to lowest order in temperature corrections, we have 2 I 2 7r D (6F) kT z 1 — — — 4. 0 “ 6F 66F D(EF) (F) ( 5) Using equation (4.47), one can simply demonstrate that D’ 262 + 4mc26 + m2c4 (6F) — F F (4.51) 5 _ F D(EF) 6%. + 3m626F + 2m2c4 which is a strictly increasing function of c F and it is easily seen that in the zero and infinite limits of 6131, it approaches 1 / 2 and 2 respectively. Therefore eFD' (6F) /D(€F) is always of the order of unity. Thus if 0 < 6 << 1, it must be that kT/ e F < 1 and it serves as the low temperature expansion parameter just as it does in the non-relativistic limit [113, 146, 147]. For the purpose of generating our statistical mechanics table in the low temperature limit, we calculate p. at the desired electron gas number density and temperature pairs and tabulate the parameter 5 E —p./kT. 4.2.3 Low Temperature Kinetic Energy Density To derive the a fully relativistic expression for the electron gas kinetic energy density UV and its derivatives, we start with the integral expression given in equation (4.22) 107 to be 7 — 1 DO 1 4 LV _ V/o E-D(E)f(E)dE (4.52) E . D(E) / V can, too, be demonstrated to satisfy the requirements that g(E) did when approximation (4.45) was derived. Therefore if we assume that By >> 1, we are free to use approximation (4.45) when evaluating the electron gas kinetic energy integral simply by replacing g(E) by E - D(E) / V in the formula. Doing so gives 6 6 W2 2 UV z /0 F E- Eng — GFLVFZ (EF — p) + F9?— (D(p) + uD'(u))4.53) The first integral in the above is clearly the ground state kinetic energy density UVgs.‘ ' c 7r2 2 —» UV z UVg.s. — eFD(VF) (6F — ,u) + Fig/1)— (D(p) + [JD/(m) (4.54) Inserting ,u = e F(1 — 6) into the above gives ~ 2 D(€Fl UV NUVgs. " 6F—V—‘S 7r2 2 + Fig/1)— (D(6F(1 — 5)) + spa — 6)D’(eF(1 — 6») (4.55) If we expand D and D’ into a Taylor series about 6 F and retain only terms up to first-order in 6, we find D(e ) 7r2D e ) 7r2 D, e ) UVzUVgfi. — a} VF 6+F (VF (kT)2—FeF—%—L6(kT)2 2 2 7r kT + —eF(1—5)(D’(6F)—6FD”(4F)6)LV—l (4.56) Carrying out the above products and again retaining only terms up to first-order in 6 gives 2 D(EF) 7r2 D(EF) 2 7r2 D’(€F) 2 2 l 2 7r D (6F) 2 7r D (61?) 2 — — — 6 + 6 V (kT) 6 (F V (kT) 2 D” — 7% F Sflmr)? (4.57) Inserting the expression for 6 implied by equation (4.50) and regrouping results in 2 2 ~ 7r 2 D(ep) k‘T UV ~ UVg.S. + FéF V (E; X 7r2 D’(6F) 2 kT 2 «2219”(6F) kT 2 .- 1-0—3—(CF D(EFl) (a?) _FEF M61?) (5) (408) It has already been established that the assumption /3,u >> 1 implies that kT/ c F << 1. It is therefore logical to suspect that the second two terms in the parenthesis above are negligible higher order temperature corrections. This can be confirmed by demonstrating that their coefficients are at most of the order of unity. Equation (4.51) was already used to prove that eFD’(eF)/D(6F) E (1/2,2), so the second term’s coefficient is of the order of unity and is therefore negligible. To conclude that same about the third term, we again use equation (4.47) to show 2 D”(€F) _ 2€%+4m0261:~ —m (5 _ F D(EF) 6%. + 4m026F + m2c4 2C4 (4.59) WhiCh is a strictly increasing function of e F and it is easily seen that in the zero and in- finite limits of 617, it approaches —1 and 2 respectively. Therefore 6%.D”(6F)/D(€F) 109 is always of the order of unity or less and the third term in the parenthesis is negligible as well. Therefore, to lowest in order in temperature correction, the electron kinetic energy density is given by Uv z UVg.s. 1,3636%?) (If—3:)2 (4.60) For the purpose of generating our statistical mechanics table in the low temperature limit, we use equation (4.60) to calculate and tabulate (8UV/6n)T and (BUV/BT)n at the desired electron gas number density and temperature pairs. The temperature derivative is trivial and the number density derivative is most conveniently calculated using the chain rule BUV = avak—F (4 61) an T (961: an ' 4.2.4 Low Temperature Pressure We start the derivation of the fully relativistic expression for the electron gas pressure by considering its general integral expression given in equation (4.13) by p: /0 h(E)f(E)dE (4.62) where = 87r 2 mcz 3/2 h(E) _ 3016)?) (E +2 E) (4.63) h(E) and its anti-derivative satisfy all the requirements that g(E) did when approx- imation (4.45) was derived. Therefore, if we assume that fly >> 1, we are free to use approximation (4.45) when evaluating the electron gas pressure by replacing g(E) by 110 h( E) in the formula. Doing so gives 6 2 P z [0 F h(E)dE —— th(cF) (CF — ,u) + %—(kT)2hl(p) (4.64) The first integral in the above is clearly the ground state pressure Pg_3,. Additionally. after a quick comparison of equations (4.47) and (4.63), one can see that h’(E) = D(E ) / V. Thus the above can be written 7r2 P z 109.3. — EFh(€F) (EF — ,u) + ?(kT)2l—)T(/Hl (4.65) Inserting ,u = 6 F(1 — 6) into the above, Taylor expanding D(E) about 61:, and retaining only terms up to first-order in 6 gives 7r2 . 2 P z Pg, — th(eF)5 + g(k—‘T/l— [D(EF) — eFD’(eF)6] (4.66) Inserting the expression for 6 implied by equation (4.50) and regroupings gives 2 2 7r 2 D(CF) kT P m P — — g.s.+ 6 {F V (6F X l—V MEF) 6FD’(‘5F) ”2 (M)2(E)2 (4.67) eFD(eF) D(eF) —6 D(eF) 617 It has been shown previously that eFD’(eF)/D(6F) E (1/2,2), therefore the last term is the parentheses above is negligible relative to the first. To determine the Order of the second term, consider the following. Using equations (4.47) and (4.63), We find that h(eF) _ 145151+2mc2 _ _ 4.68 €FD(€F) 3 5F+mc2 ( ) 111 which is a strictly decreasing function of e F that has a maximum value in the e F ——+ 0 limit of 2 / 3 and is minimized when 6 F —+ 00 as it approaches 1 / 3. Thus the second term in the parentheses above is for all e F of the order of unity and must therefore be kept. Therefore, to lowest in order in temperature correction, the electron gas pressure is given by h(eF) eFD’(€F) eFD(eF) D(eF) 1 — V (4.69) (W)? PgP9'3‘+_€F V 6F Proceeding in an identical fashion as before, for the purpose of generating our statis- tical mechanics table in the low temperature limit, we use equation (4.69) to calculate and tabulate (8P/8n)T and (BF/8T») at the desired electron gas number density and temperature pairs. Again we find that the temperature derivative is trivial and the number density derivative is most conveniently calculated using the chain rule 8 (all?) = 313.2 (4.70) (977. T (96F an 4.3 High Temperature Statistical Mechanics When an ideal gas of fermions in is said to be in the high temperature limit, the Chemical potential is negative and sufficiently large that 5 E —)u/kT >> 1 [113, 146, 147]. This means that for any fermion kinetic energy E, exp(€ + E / kT) >> 1 and the F ermi—Dirac distribution simplifies to f(E) -’ eXIX—6 — E/kT) (4-71) Since exact expressions for many statistical properties of an ideal Fermi gas can be WI‘it.ten as indefinite integrals of functions proportional to the Fermi-Dirac function, the above simplification allows the g parameter dependence to be factored out of 112 these integrals. In addition to this factorization simplifying the evaluation of these integrals, it also allows equations (4.2) or (4.7) that relate 5 to the number density and temperature of the electron gas to be solved algebraically instead of the resorting to the laborious numerical techniques used in the arbitrary temperature case to solve for g. As was done for the arbitrary and low temperature cases, the fully relativistic high temperature expressions for the electron gas 5 parameter, kinetic energy density, and pressure are rigorously derived and their convergence with their well known non- relativistic limits are shown. 4.3.1 High Temperature .5 Parameter To derive an expression for the 6 parameter of an ideal Fermi gas in the high tem- perature limit, we start with the general equation (4.7) that implicitly relates the 5 parameter to the fermion gas number density n and temperature T through the condition _ 00% n— /0 v f(E)dE (4.72) Where again D(E) is the fully relativistic expression for the density of states shown in equation (4.20) to be given by ESL) = 0:733 (E2 -+— 2m62E)1/2 (E + 777.02) (4.73) 113 If the high temperature limit of the Fermi-Dirac distribution given in equation (4.71) is substituted into equation (4.72), g is easily solved for. Doing so yields 5 = 111 (l [000 DE/E) exp(—E/kT)dE) In (87” (1%)} 15m) (4.74) with [g(T) defined by [€(T) E /()00 (142 + 2au)1/2(u + a) exp(—u)du (4.75) where the dimensionless variable 11. E E / kT was introduced to obtain a dimensionless integral and a E mcz/kT. While the above integral must be evaluated numerically, since it depends only upon the temperature, the tabulation process is quite fast and the high temperature end of our 5 parameter table is easily generated. 4.3.2 High Temperature Kinetic Energy Density It is clear from the general integral expression for the electron gas kinetic energy density given in equation (4.19) that in the high temperature limit UV —> exp(—.§)/OOO§.—€(—E;)exp(—E/kT)dE _ kT 3 E exp(—€) 87r - kT (E) [U(T) (4.76) with I U(T ) defined by [U(T) E ADO 24(142 + 2au)1/2 (u + a.) exp(—u)du (4.77) 114 where again the dimensionless variable 24 E E / [CT was introduced to obtain a dimen- sionless integral and a E mc2/kT. This integral must be evaluated numerically as well. If the expression for 5 from equation (4.74) is inserted into the above equation for UV, we find that in the high temperature limit - nkT (4.78) For the purpose of generating our statistical mechanics table in the high temperature limit, we use equation (4.78) to calculate and tabulate (8UV/0n)T and (aw/air)“ at the desired electron gas number density and temperature pairs. The number density derivative is trivial and the temperature derivative is calculated by first dif- ferentiating (4.78) using the chain rule and then performing the resultant numerical integrations. 4.3.3 High Temperature Pressure Recall from equation (4.13) that the general integral expression for the electron gas pressure is 00 P =/0 h(E)f(E)dE (4.79) where h(E) is the anti-derivative of D(E ) / V given by 77 3 2 h(E) = 353—3 (E2 + 2m02E) / (4.80) 115 In the high temperature limit, equation (4.71) implies 00 P —> exp(—€)/0 h(E)exp(—E/kT)dE 7r 3 E exp(—§)-8§kT (Eh—70:) [P(T) (4.81) with I P(T) defined by Ip(T) E [000 (142 + 2au) 3/2 exp(—u)du (4.82) where again the dimensionless variable u E E / [CT was introduced to obtain a dimen- sionless integral and a E ch/kT. This integral must be evaluated numerically as well. If the expression for g from equation (4.74) is inserted into the above equation for P, we find that in the high temperature limit ‘ 1 For the purpose of generating our statistical mechanics table in the high temperature limit, we use equation (4.83) to calculate and tabulate (BF/870T and (8P/3T)n at the desired electron gas number density and temperature pairs. The number density derivative is trivial and the temperature derivative is calculated by first differentiating (4.83) using the chain rule and then performing the resultant numerical integrations. 116 4.4 Electron Gas Statistical Mechanics in the Non- Relativistic Limit While it is only a necessary condition for the fully relativistic electron gas statistical mechanics derived in sections 4.1, 4.2, and 4.3 to converge with non-relativistic elec- tron statistical mechanics and not a sufficient case for its validity, it is still instructive to demonstrate this convergence. In this section, the non-relativistic limits of the expressions derived for the relativistic electron gas statistical mechanical quantities of interest with arbitrary, low, and high temperatures are all determined separately. We begin with the non-relativistic limit of electron gases with arbitrary temperatures. 4.4.1 Non-Relativistic Limit of Arbitrary Electron Gases For arbitrary electron gases, we consider the non-relativistic limits of the electron gas density of states, number density, pressure, and kinetic energy density. Considering first the density of states given in equation (4.20), we notice that the non-relativistic limit E << mc2 implies D(E) = V 0:733 (E2 + 2mc2E)1/2 (E + 77102) 1 2 -—> V 87f (2mc2E) / mc2 (hc)3 = V313 (2m)3/2 131/2 (4.84) which is known to be the expression for the density of states of a non-relativistic electron gas [147]. In terms of the dimensionless variable 14 = E //6T and parameter a. E meg/kl" defined in the integral equations (4.8), (4.14), and (4.22) for n, 1’, and 117 UV the condition E << mc2 implies u << 0.. Thus for these integrals we have (kT)3 /00 (212 + 2au)1/2(u + a)du 877 — hc 0 exp(€ + 14) +1 3 00 1/2 87r (E) 21/2a3/2/ u (121 h 0 exp c (g + 11) +1 43 3/2 00 ul/zdu h3 (2ka) /0 exp“ + 10+ 1 (4.85) 3/2 3 oo 62+2au du 871' (hr) kT/o ( ) _3_ E exp({ + u) + 1 87rkT kl)3 (2a)3/2 foo 113/2cm 3 he 0 exp(£ + u) + 1 87rkT 3/2 foo u3/2du 2 [CT 4.86 3h3 ( m ) 0 exp(€ +14) +1 ( ) 3 2 1 2 87rkT (Ell) [00 u(u + 2au) / (u + a)du c 0 exp(€ + 11) +1 3 3 2 87rkT (LT) 21/2a3/2 foo “ / d“ he 0 exp(€ + u) + 1 47rkT 3/2 00 113/2dr]. . 53 (2ka) f0 exp(€+u)+1 (4.87) where equations (4.85), (4.86), and (4.87) are known to be the non-relativistic ex— pressions for the number density, pressure, and kinetic energy density of a Fermi gas [146]. 118 4.4.2 Non-Relativistic Limit of Low Temperature Electron Gases For low temperature electron gases, we consider the non-relativistic limits of the electron gas chemical potential, kinetic energy density, and pressure. Starting with the chemical potential, we note that from equation (4.51), it is clear that the condition E << mc2 implies D’(EF) 1 5F D(EF) —* 2 (4.88) Substituting this result into the fully relativistic expression for the low temperature electron gas chemical potential given in equation (4.50) yields 1_ 71;: (g) 2] (4.89) which known to be the non-relativistic expression of the chemical potential in the low #361: temperature limit [146]. It is also known that in the non-relativistic limit that the density of states can be written as [147] v 7’ E 3/2 6F E 3 ” E1/2 (4.90) Inserting this into the fully relativistic expression for the low temperature electron gas kinetic energy density given in equation (4.60) gives 2 2 77 [CT 2 _ — 4. UV Uvg_3_+ 46pm (5F) ( 91) 119 Since the ground state kinetic energy density in the non-relativistic limit is given by UVg.s. = 3/5neF [146], the above can be rewritten as 2 2 Sr kT (4.92) which known to be the non-relativistic expression of the kinetic energy density in the low temperature limit [146]. Finally, we consider the non-relativistic expression for the low temperature electron gas pressure. From equation (4.68), we see that in the non-relativistic limit that V h(€F) _) 2 €FD(€F) g (4.93) Inserting this result and equations (4.88) and (4.90) into the fully relativistic expres- sion of the low temperature pressure of the electron gas given in equation (4.69) leads to 2 H 2 EF Since the ground state electron gas pressure in the non-relativistic limit is given by P95, = 2/ 5ne F [146], the above can be rewritten as 2 2 57r kT wa. [1+17(—) (4.95) 6F which known to be the non-relativistic expression of the electron gas pressure in the low temperature limit [146]. 120 4.4.3 Non-Relativistic Limit of High Temperature Electron Gases For high temperature electron. gases, we consider the non-relativistic limits of the electron gas 5 parameter, kinetic energy density, and pressure. Beginning with the 5 parameter, we observe that for the I €(T) integral defined in equation 4.75, the non-relativistic condition 11 << 0. implies that [€(T) = /000 (112 + 2au)1/2(u + a) exp(——u)d’u 00 __, %(2(1)3/2~/O 211/2 exp(—u)du (4-96) The above integral is recognized as the generalized factorial (1 / 2)! = (1 / 2)( -1 / 2)! : (1/2)\/F. Therefore we have 1 I€(T) —+ 21——(27m)3/ 2 (4.97) 71' Plugging this result into the high temperature 5 equation (4.74) gives 3 n h _. _1 _ _— E “[2( 27rka) — 16033) + ln(2) (4.98) where A is the so-called called the thermal wavelength [146]. Apart from the additional ln(2), the above is the well known non-relativistic Boltzmann expression for the 5 parameter [146, 147]. Two things should be noted about this extra ln(2) term. One is its origin in spin degeneracy. It is possible to build a spin correction into Boltzmann theory and reproduce this extra ln(2) term by simply taking into account fermion spin 121 degeneracy when constructing the partition functions needed to derive the expression for the 5 parameter. The other is that from the standpoint of numerical convergence with the Boltzmann result, it is inconsequential. Recall that in the high temperature limit, 5 >> 1 and since ln(2) a: 0.693, the dominant contribution must come from the expected term. As the system moves farther and farther toward the high temperature limit, this discrepancy will become completely negligible. For the I U(T ) integral defined in equation 4.77, the u < a condition implies [U(T) = 400 11(112 + 2au)1/2(u + a) exp(—u)du —> $(2a)3/2 ADO 113/2 exp(—u)du = é—(26)3/2(3/2)! = §;(27ra)3/2 (4.99) Inserting this result and the non-relativistic limit of the integral I g(T) from equation (4.97) into equation (4.78). for the electron gas kinetic energy density gives UV _. gnkT (4.100) which is identified as the non-relativistivic Boltzmann expression for the kinetic energy density of a gas [146]. Finally, for the [P(T) integral defined in equation 4.82, the 122 u << a condition implies 1P(:r) = I/Ooo(u2+2au)3/2exp(—u)du -—* 03/2 00U3/28X —-u U (2) [0 p( )d = (2a)3/2(3/2)! = 4%(2m3/2 (4.101) Inserting this result and the non-relativistic limit of the integral I {(T) from equation (4.97) into equation (4.83). for the electron gas pressure gives P —> nkT (4.102) which is identified as the non—relativistivic Boltzmann expression for the pressure of a gas [146]. 4.5 Comparing Relativistic and Non-Relativistic Electron Gas Statistical Mechanics Having established the convergence of the relativistic electron gas statistical mechan- ics theory derived in sections 4.1, 4.2, and 4.3 with non-relativistic electron gas sta- tistical mechanics, we now numerically compare the tabulated electron gas statistical mechanical quantities computed with the relativistic and non-relativistic formalism at selected electron gas number densities and temperatures. This is done to gain a quantitative sense of what temperatures and electron number densities at which 123 the relativistic and non-relativistic electron gas statistical mechanics converge and what how they differ over the range of temperatures and electron number densities expected to be relevant to collapse of the core being modeled. 4.5.1 Numerical Non-Relativistic Convergence To demonstrate the numerical convergence of the relativistic and non-relativistic elec- tron gas statistical mechanics, we access the low temperature low electron number density corner of our tables. Specifically we consider electron gases with tempera— ture T = 107 K and number densities 1033 rn—3 S ne 3 1036 m—3. We choose to consider an electron gases with the lowest tabulated temperature 107 K for the obvious reason that they will have the lowest thermal velocities and hence will be the least likely to require relativistic corrections to their thermal velocity distributions. We consider electron gases with the four lowest number densities tabulated since, at T = 107 K, the effects of degeneracy are not negligible for 716 Z 1033 m‘3 and as number densities become sufficiently large, the degeneracy becomes extreme and relativistic. In figure 4.1, the 5 parameter is plotted for relativistic (blue) and non-relativistic (red) electron gases with temperatures and number densities specified above. The 5 values are seen to be in excellent agreement at T = 107 K up to approximately 3 x 1034 m_3. For all greater number densities, the non-relativistic 5 parameter be- comes more negative. Recall that 5 E ,u/kT, therefore a more negative 5 parameter means a more degenerate gas. This tendency of non-relativistic electron gases to become more degenerate than relativistic gases with the same number densities and temperatures is reflected in figures 4.2 and 4.3. In figure 4.2 it is seen that as the number densities increase, the non-relativistic gases become increasingly less sensitive to temperature compared to their relativistic counterparts which is a direct result of their greater degeneracy. In figure 4.2 it is seen that as the number densities increase, 124 — Relativistic ‘-- Non-Relativistic r -100 -150 \ -200 \\ ”.2... \\ -300 \\ \\ \ '450 1 I 1 r T 33 33.5 34 34.5 _3 35 35.5 36 '0910("e) [m ] Figure 4.1: Plot of the 5 parameter for relativistic (blue) and non-relativistic (red) electron 3:888 with temperature T = 10 K and number densities 1033 m’3 g 716 g 1036 m— . the non-relativistic gases become increasingly more sensitive to number density com- pared to their relativistic counterparts which is also a direct result of their greater degeneracy. At sufficiently low number densities, the partial derivatives of the electron gas pressures and kinetic energy densities converge in figures 4.2 and 4.3. The divergence of the temperature derivatives of the relativistic and non-relativistic kinetic energy densities occur at lower number densities in figure 4.2. Conversely the divergence of the number density derivatives of the relativistic and non-relativistic pressures occur at lower number densities in figure 4.3. These asymmetries in the relativistic and non-relativistic divergencies are attributable to subtle differences in the temperature 125 and number density dependancies between the relativistic and the non-relativistic pressure and kinetic energy densities discussed in previous sections in this chapter. 115 — Relativistic log(aP/6T) [Pa/K] T /1 11.3 — — Relativistic log(au/aT) [J/(mA3K)] - Non-Relativistic log(aP/aT) [Pa/K] / 11'1 ‘ wwwNon-Relativistic log(au/8T) [J/(mABV ‘ 10.9 10.7 10.5 // 10.3 ,_,/ 10.1 ”if / 9.9 9.7 9.5 I 1 i . I 33 33.5 34 34.5 35 35.5 36 -3 '0910(“e) [m ] Figure 4.2: Plots of BP/BT and Bu/ET for relativistic (red / blue) and non-relativistic (green/ turquoise) electron gases with temperature T = 107 K and number densities 1033 m—3 g 71.5 g 1036 66-3. 126 -13 — Relativistic log(aP/an) [Nm] — Relativistic log(au/an) [J] -— Non-Relativistic Iog(aP/an) [Nm] - Non-Relativistic log(au/an) [J] -13.5 ~- -14 -14.5 -15 ‘15.5 i 1 1 i l 33 33.5 34 34.5 35 35.5 36 '0910(“e) [m'3] Figure 4.3: Plots of BP/Bn and 311/877. for relativistic (red/ blue) and non-relativistic (green/ turquoise) electron gases with temperature T = 107 K and number densities 1033 m_3 S 116 S 1036 m-3. 127 4.5.2 Relativistic and Non-Relativistic Electron Gases in Su- pernovae In all of the finite temperature calculations our code has performed, not just those presented in this thesis, the lowest temperatures and electron number densities en- countered were on the order of 107 K and 1037 m”3 respectively and the highest were on the order of 1011 K and 1045 m"3 respectively. In this section we compare the tabulated statistical mechanical quantities computed with the relativistic and non-relativistic formalism of electron gases at temperatures 109 K and 1011 K over the number density range 1037 m-3 5 715 _<_ 1045 m—3 to obtain a measure of how the results differ in the environment similar to the system we model. We begin these considerations with a plot of the 5 parameter of relativistic (blue) and non-relativistc (red) electron gases at 109 K presented in figure 4.4. Clearly the non-relativistic gases are much more degenerate than the relativistic gases. This greater degeneracy of non-relativistic electron gases is merely a continuation of the trends seen emerging in the last section and is also seen in the plots of the temperature and number density derivatives of the pressure and kinetic energy density of electron gases at 109 K shown in figures 4.5 and 4.6. Again it is found that the non-relativistic gases are less sensitive to temperature changes and more sensitive to number density changes than the relativistic gases due to their higher level of degeneracy. Clearly relativistic effects are not negligible. 128 -1.E+04 - Relativistic — Non-Relativistic -2.E+04 w! -3.E+04 '4.E+04 ( ‘5.E+04 t i 1 i 37 39 41 43 45 '0910("e) [m8] Figure 4. 4: Plot of the 5 parameter for relativistic (blue) and non-relativistic (red) electron In__g3ases with temperature T— - 109 K and number densities 1037 m 3 S 716 S 1045 129 20 — Relativistic log(aP/aT) [Pa/K] 19 __ — Relativistic log(au/aT) [J/(mABK)] / -- Non-Relativistic log(aP/aT) [Pa/K] 18 ........ Non-Relativistic Iog(au/aT) [J/(mA3K)] / 17 / l I 41 43 45 logiome) [in-3] Figure 4.5: Plots of BP/aT and Bu/BT for relativistic (red/ blue) and non-relativistic (green/ turquoise) electron gases with temperature T = 109 K and number densities 1037 m-3 g 77.6 g 1045 111-3. 130 — Relativistic Iog(aP/an) [Nm] -— Relativistic log(au/an) [J] / ‘8 i" —Non-Relativistic Iog(aP/an) [Nm] / _9 Non-Relativistic Iog(au/an) [J] -10 . ‘////////f -11 1///’//,/' 4‘,,.——“" -12 II22/,zc’I/rzI:::T’7’5:::::::::::”’F’Ffidr -13 2”:::::::::::1——————'r I -.._.l -14 . . . 37 39 41 43 45 -3 '0910("e) [m ] Figure 4.6: Plots of BP/an and 014/071 for relativistic (red/ blue) and non-relativistic (green/ turquoise) electron gases with temperature T = 109 K and number densities 1037 m’3 S 126 S 1045 m—3. 131 Finally we consider electron gases 1011 K. In figure 4.7, the 5 parameter is plot- ted for relativistic (blue) and non-relativistic (red) electron gases with temperature T =1011 K and number densities 1037 m—3 S ne 3 1045 m‘3. The non- 20 1 — Relativistic -- Non-Relativistic 0- -20 “100 T i 1 37 39 41 43 45 '0910(“e) [m'3] Figure 4.7: Plot of the 5 parameter for relativistic (blue) and non-relativistic (red) electron g‘aases with temperature T = 1011 K and number densities 1037 m-3 _<_ 713 S 1045 m' . relativistic electron gases are once again found to be more degenerate than the rela- tivistic gases, however at T = 1011 K both relativistic and non-relativistic electron gases are seen to be partially non-degenerate (5 Z 0) at sufficiently low number densities. As in our previous considerations of the temperature and number density deriva- tives of the pressure and kinetic energy density of electron gases, the plots displayed in figures 4.8 and 4.9 indicate that the more degenerate nature of the non-relativistic 132 22 — —- Relativistic log(aP/aT) [Pa/K] 21 __ —Re|ativistic log(au/aT) [J/(mA3K)] / ._ Non-Relativistic Iog(aP/aT) [Pa/K] //' 20 -3 WW Non-Relativistic Iog(au/aT) [J/(mA3K)] ; 19 1’1”:::::;”’r 14 . . 1 37 39 41 _3 43 4s I0910(ne) [m ] Figure 4.8: Plots of BP/BT and (911 / ET for relativistic (red/ blue) and non-relativistic (green/ turquoise) electron gases with temperature T = 1011 K and number densities 1037 m_3 S 716 S 1045 m—3. electron gases renders them less sensitive to temperature changes and heightens their sensitivity to increases in number densities relative to relativistic gases. Additionally the relativistic and non-relativistic temperature and number density derivatives of the pressure converge on this range of number densities. 133 _7 5 +5 —Relativistic log(aP/an) [Nm] 2 ' —Re|ativistic Iog(au/an) [J] /t[ -8 -_. —- Non-Relativistic log(aP/an) [Nm] ,1 Non-Relativistic log(au/an) [J] / [ -10.5 // / - 1 1 ' -12 , , . 37 39 41 43 45 “3910019) [m8] \ Figure 4.9: Plots of 8P/ 8n and 811/871 for relativistic (red/ blue) and non-relativistic (green/ turquoise) electron gases with temperature T = 1011 K and number densities 1037 m’3 g ne 5 1045 m“3. 134 Unlike the pressure derivatives, the temperature and number density derivatives of the kinetic energy density do not converge on this range of number densities. This asymmetry can be attributed to the differences between relativistic corrections to the temperature dependence of the pressure and kinetic energy density and the relevancy of these corrections at 1011 K. At 1011 K, the difference between the relativistic and non-relativistic results are greater than they were at 109 K. Since the difference was already significant at 109 K, it is clear that using the relativistic formalism to characterize the electron gas is imperative. Not only does the relativistic nature of the electron gas profoundly impact core dynamics directly through the pressure it exerts, but its relativistic statistical properties also do so in a more indirect fashion by governing weak reactions that involve initial and/ or final state electrons and the impact that they have on matter. 135 Chapter 5 Weak and Strong Reaction Tables In this chapter, the tables used to model weak and strong reactions are discussed in detail. These tables are divided into three groups: those that govern the rate of elec- tron capture, those that govern the impact the neutrino production and interactions have on matter and the flow of neutrinos, and those that govern fusion. We describe these tables separately in the sections below. 5.1 Electron Capture Rate Tables A total of four electron capture rate tables are used to study the role that electron capture plays in the collapse and explosion mechanics. As stated in section 3.12, the FFN electron capture table serves as the foundation for our treatment of electron capture. The FFN table contains the base 10 logarithm of the electron capture rates for free protons and 187 nuclei with atomic mass numbers 21 S A _<_ 60, most of which are in or near the valley of beta stability. Since we model the propagation of 385 nuclei that exist between the proton and neutron drip lines with A S 60, it is necessary to extrapolate the FFN electron capture rates along each isobar out to the proton and neutron drip line as well down to smaller nuclei with 2 g A g 20. The 136 isobaric extrapolations are done in one of three ways described below. The first technique is to simply linearly extrapolate the base 10 logs of the electron capture rates using the two closest tabulated values to the proton-rich and neutron- rich boundaries of the FF N table for each A. The second way is to use an enhanced linear extrapolation that doubles the rate at which the linearly extrapolated base 10 logs of the rates increase as one moves toward the proton drip line and halves the rate at which they decrease as one moves towards the neutron drip line. The third way is to use a reduced linear extrapolation that conversely halves the rate at which the linearly extrapolated base 10 logs of the rates increase as one moves toward the proton drip line and doubles the rate at which they decrease as one moves towards the neutron drip line. The extrapolation to smaller nuclei with A < 21 is currently only done one way. To extrapolate the electron capture rate to a nucleus (Z, A) with A S 20, we simply equate it with the smallest nucleus with the same proton to neutron ratio for which an electron capture rate is either given by FFN table or has been extrapolated. Other extrapolation techniques that take into account the parity of proton and neutron numbers of the nuclei are possible, but have not been implemented yet. Thus far, we have described the methods used to generate three modified version of the F FN rate table. The results produced by simulations that use these different tables can be used to study the sensitivity of the core dynamics to the values used for the extrapolated F FN electron capture rates. To study the core dynamic’s overall dependence on electron capture rates, another table is needed that uses different rates from FFN. Motivated by the fact that more modern shell model calculations typically yield electron capture rates a full order of magnitude or more lower that the FFN rates [16], we satisfy this need by simply subtracting one from each of the base 10 logarithms of the electron capture rates from the linearly extrapolated F FN table. This order of magnitude reduced version of the linearly extrapolated FFN table serves 137 as the fourth electron capture rate table used to probe the role of electron capture rates in core collapse supernova mechanics. 5.2 Neutrino Production and Interaction Tables Recall from sections 3.12.2, 3.13.3, 3.13.4, and 3.16, that the algorithms used to model the weak reactions included in our simulation and their effects on the temperature distribution all make use of tables of average interactions cross sections and / or average particle energies. We tabulate these all of these average cross sections and average particle energies logarithmically over the same temperature range as we did for the electron gas statistical mechanics table in section 4. For weak reactions involving electrons, these cross sections and energies are tabulated over the same range of electron number densities as we did for the electron gas statistical mechanics table. For weak reactions that involve an initial state neutrino, these cross sections and energies are tabulated logarithmically over incident LTE frame neutrino energies from 10"4 MeV to 104 MeV. In the sections below, we describe how we tabulate the average cross sections and/ or average particle energies needed to model each weak reaction. 5.2.1 Neutrino Production Tables To model the production of neutrino test particles via electron captures by the nucleus or free proton (A, Z) and the effects that these reactions have on the temperature distribution, it was explained in sections 3.12.2 and 3.16 that it is necessary to tabulate the average energy of an electron captured by the nucleus or free proton (A, Z) from an electron gas with selected number densities n and temperatures T as well as the average energy of the neutrino produced during the process for the purpose of interpolation. This is done in the following way. 138 At each (n, T) pair in the interpolation lattice, we systematically sample the ther- mal energy distributions of the electrons 5000 times and model the potential capture of each of these electrons by 5000 nuclei or free protons (A, Z) systematically sampled from the thermal energy distributions of (A, Z) particles. The way that these thermal energy distributions are calculated are explained in section 5.2.2. These potential cap- tures are modeled using relativistic classical mechanics. Once the LTE frame kinetic energies of the electron and (A, Z) particle have been generated, the magnitude of their momentum vectors are determined and then the vectors are randomly oriented. Then we simply boost into the C-O-M frame of the selected electron-(A, Z) pair and check to see if the C-O-M frame energy is large enough to overcome the Q-value of the electron capture. If it is not, we reject the capture and move on to the next pair. If it is, we model the capture by changing randomly re-orienting the C-O-M frame momentum vectors redefining the C-O-M frame to be that of the new neutrino and final state (A, Z — 1) nucleus or free neutron, and then boosting the neutrino’s mo- mentum into the LTE frame. This process is similar to the two—body elastic scattering model described in section 3.13.4 and is schematically depicted in figure 5.1. The final LTE frame energy of the neutrino is calculated and then its contribution, along with that of the initial state electron’s LTE frame kinetic energy, are added to the average energies tabulated for the particular (71, T) pair in the interpolation lattice. Since 25 million potential electron captures are modeled at each (n, T) pair in the interpolation lattice, we are confident that our statistics are good. This procedure is followed for all nuclei and free protons (A, Z) included in our simulation. 5.2.2 Relativistic Thermal Energy Distributions To generate fully relativistic expressions for the energy distributions for the ideal gases of electrons, free baryons, and nuclei involved in weak reactions, we divide our efforts in the following way. First we derive a fully relativistic energy distribution valid for 139 Inelastic cm frame Figure 5.1: Two-body inelastic scattering is modeled by randomly repositioning the C-O—M frame momentum vectors to opposite positions on the surface of the momen- tum sphere defined by their C-O-M frame momentum and then changing the mag- nitude of the C-O~M frame momentum vectors in the way that energy conservation dictates. any ideal gas of indistinguishable particles, fermion, bosons, or Boltzmann particles. Once this general expression has been obtained, we tailor it to describe the energy distributions of an arbitrary electron gas and the gases of free baryons and nuclei which are assumed to be non-degenerate. After these specific expressions have been derived, we demonstrate their convergence with the non-relativistic non-degenerate Maxwell velocity distribution. Relativistic Energy Distributions of Arbitrary Gases Let (n(E)) be the average number of particles in the state with kinetic energy E. By the definition of the density of states D(E), the number of particles in a gas is N is 140 given by [146] N =/0 D(E)(n(E))dE 1 oo =>1 = W/o D(E)(n(E))dE (5.1) From the above is is clear that if the range of integration is restricted to some finite energy interval, the integral would yield the fraction of particles in that energy inter- val. This is interpreted as the probability of finding a particle in this energy range. Thus we identify the probability of a particle being found within some infinitesimal energy range dE centered about kinetic energy E as being given by 1 P(E)dE = ND(E)(71(E))dE (5.2) Recall that up to a constant independent of kinetic energy E, the fully relativistic expression for the density of states satisfies D(E) oc (E2 + 2mc2E)1/2(E + ch). Therefore the above can be written as P(E)dE = $132 + 2mC2E)1/ 2(13 + mc2)(n(E))dE (5.3) where we have defined N = [000 (E2 + 2mc2E)1/ 2(13 + mc2)(n(E))dE (5.4) Therefore it is clear that the probability of finding a particle in the energy interval [E1, E2] is given by 1 E2 2 2 1/2 2 P(El g E 3 E2): W /E (E + 2mc E) (E+ mc )(n(E))dE (5.5) I 141 Armed with this general expression for the energy distribution of a gas of indistin- guishable particles, we can now proceed with the derivations of the energy distribu- tions of specific gases. Electron Gas Energy Distribution During a core collapse supernova, the electron gas can initially be assumed to be approximately degenerate [18, 62], but as the collapse ensues the degeneracy condition is eventually lifted as electrons are “up-scattered” by increasingly energetic neutrinos [62]. Thus, for the electron gas, the most general expression for (n(E)) must be used as no simplifying assumptions can be made about the degeneracy of the energy distribution that would remain valid at all times. Since electrons are fermions, the most general expression for the average number of particles in the state with certain kinetic energy E is given by the Fermi-Dirac distribution 1 = exp(5+ E/kT) +1 => MB) (5.6) where again 5 is related to the chemical potential 11 via 5 E ——;1/kT. The electron gas 5 parameter is tabulated at all the number density temperature pairs that the electron gas energy distribution must be determined and sampled. Thus we can divide the electron energy distribution into bins for the purpose of systematic sampling by numerically evaluating the integrals in equations (5.4) and (5.5) using the techniques extensively described in section 4.3 where similar integrals were encountered when electron gas statistical mechanics at arbitrary temperatures was discussed. Free Baryon and Nuclear Energy Distributions In our simulation, the gases of free baryons and nuclei are assumed to be non- degenerate or equivalently in the high temperature limit. If we prescribe to the 142 generally excepted picture of core collapse supernovae though bounce described in section 2.2, this assumption is safe to make for the free proton and nuclear gases. Due to their lower number densities and there larger particulate masses relative to the electron gas, in all regions of the core, their average thermal kinetic energies will be much larger than their Fermi energies rendering the effects of degeneracy negligible. However this picture, the assumption of non-degeneracy is not valid for the neutron gas which is believed to become degenerate and halt the collapse. This neutron gas degeneracy would then result in the phase space blocking of some electron captures by free protons and neutron drip line nuclei that would otherwise produce more free neutrons. The effect of this would be to reduce the number of neutrinos produced in regions of the core with sufficiently high free neutron number densities. In these regions, general statistical mechanics must be used to describe the neutron gas. It is clear how to integrate this into our approach, but currently it is not implemented. For gases of fermions and bosons, in the high temperature limit where exp(5) >> 1, the average number of particles in the state with a kinetic energy E approaches the Boltzmann result 1 (n(E)) = exp(5 + E/kT) i1 —-> exp(—5 — E/kT) (5.7) This simplification allows the 5 ’s to be factored out of both integrals in equations (5.4) and (5.5) and cancel out when the energy distributions are determined by their ratio. This is quite convenient as it removes the implicit dependence the energy distributions had on particle number density through the 5 parameter and makes them functions of temperature alone. The resultant integrals of functions proportional to Boltzmann factors can also be numerically evaluated using the techniques thoroughly discussed in section 4.3 allowing us to divide these energy distributions into bins for the purpose of systematic sampling. 143 Maxwell Limit While it is only a necessary condition for the fully relativistic energy distributions de- rived above to converge with their well known non-relativistic non-degenerate Maxwell limits, and not a sufficient case for their validity, it is still instructive to demonstrate this convergence. In the non-degenerate limit, as demonstrated above, in equations (5.4) and (5.5), (n(E)) can be replaced by exp(—E/kT). Additionally, it was shown in section 4.4.1 that in the non-relativistic limit the condition E << mc2 implies that D(E) —2 CEl/Z, where C is a constant independent of energy. Therefore in the non-degenerate non-relativistic limit, the N integral defined in equation (5.4) becomes N = C/OOOEl/2exp(—E/kT)dE CH = C(kT)3/2 [000 u1/2exp(—u)du ( .8) Again the above integral is recognized as the generalized factorial (1 / 2)! = (1 / 2)( —1 / 2)! where (—1/2)! = fl, and the above becomes N = C(kT)3/Qi2f- (5.9) The energy distribution given in equation (5.3) therefore becomes P(E)dE = %E1/ 2 exp(-E/kT)dE = WE”? exp(—E/kT)dE (5.10) 144 Furthermore, in the non-relativistic limit, the kinetic energy is given by E = mv2 / 2 and the energy differential is given by dE = mv dv. After inserting these into the above, we arrive at 3 2 P(E)dE = (27:21,) / 47w2 exp(—m'v2/2kT) dv (5.11) The right hand side of the above is recognized as the probability of finding a particle in a Maxwellian gas of indistinguishable particles at temperature T within a differential neighborhood dv of a velocity 0. Thus the non-relativistic non—degenerate limit of the fully relativistic energy distributions derived in the pervious sections are consistent with the Maxwell picture as they must be. Comparing Relativistic and Non-Relativistic Energy Distributions Having established the convergence of the relativistic energy distributions with the non-degenerate non-relativistic Maxwell limit, we now numerically compare selected energy distributions computed with the relativistic and non-relativistic formalism. First we consider electrons. In figure 5.2, the relativistic energy distribution (blue) and non-relativistic energy distribution (red) for a gas of electrons with number density 713 = 1033 m_3 and temperature T = 107 K are expressed as functions of velocity and plotted over the range 0 S v < 0. Note that these electron gas number densities and temperature are not found in the hot dense core. Such conditions exist outside the core in the cooler more diffuse regions of the progenitor. Two observations are made about these energy distributions. One is that they agree very well with one another. This is expected as relativistic corrections to finite temperature thermal distributions of non-degenerate electron gases are not significant for T << 109 K [113]. Secondly, they both exhibit the classic Maxwellian curve shape and in velocity space predict a most probable velocity ”prob = ‘/2kT/me. This too is due to the lower 145 0.25 — Relativistic Energy Distribution A — Non-Relativistic Energy Distribution 0.2 1 0.15 0.1 0.05 l 0 l l l l 0 0.2 0.4 0.6 0.8 1 VIC Figure 5.2: Plot of the relativistic energy distribution (blue) and non-relativistic energy distribution (red) for a gas of electrons with number density 713 = 1033 m—3 and temperature T = 10 K expressed as functions of velocity. temperature and non-degenerate nature of the electron gas. Next we consider an electron gas with number density 713 = 1039 m‘3 and tem- perature T = 1010 K. The relativistic energy distribution (blue) and non-relativistic energy distribution (red) are expressed as functions of velocity and plotted for this gas are plotted over the range 0 S v < c in figure 5.3. Electron gases with these properties can be found in the core at intermediates stages of the collapse outside the hottest and densest region of the core. From figure 5.3, it is clear that relativistic and effects are quite pronounced for this electron gas. The electron gas treated with the fully rel- ativistic formalism is found to be partially degenerate and its velocity distribution is narrowly peaked about U x 0.9945 c and predicts no electron having velocities greater than c. On the other hand, the electron gas treated with the non-relativistic formal- 146 0.25 - Relativistic Energy Distribution — Non-Relativistic Energy Distribution 0.2 0.15 0.1 0.05 0 . . 0 0.2 0.4 0.6 0.8 1 VIC Figure 5.3: Plot of the relativistic energy distribution (blue) and non-relativistic energy distribution (red) for a gas of electrons with number density 713 = 1039 In”3 and temperature T = 1010 K expressed as functions of velocity. ism is found to be highly degenerate with an unphysical Fermi velocity UF 3 12 c and incorrectly predicts that virtually all of the electrons have velocities greater than c. Since the density of velocity states in this case is well approximated as being homo- geneously distributed over the velocity sphere v S v F1 the fraction of electrons with v S c z 1/12 ”F is approximately given by (1/12)3 x 5.79 x 1041. This is why the area under the non-relativistic energy distribution curve is so small over the velocity range 0 S v < 0. Finally we consider the energy distributions of baryonic matter. Since we assume that baryonic matter is non-degenerate, here we compare the non-degenerate rela- tivistic energy distribution of baryonic matter expressed as a function of velocity to the Maxwell distribution. We suspect that relativistic effects will not be as significant 147 in this case due to the larger particulate mass. Thus we begin our considerations of baryonic energy distributions with the baryonic gas that is most likely to become relativistic: a high temperature free proton gas. Since free protons are the lightest type of baryonic matter included in our simulation, they are the most likely to ac- quire large thermal velocities. The highest temperatures expected to be relevant to supernova dynamics are on the order of T = 1011 K [56]. In figure 5.4, the relativistic energy distribution (blue) and non-relativistic energy distribution (red) for a gas of non-degenerate protons with temperature T = 1011 K are expressed as functions of velocity and plotted over the range 0 S v < 0. Even at this extreme temperature, 0.5 — Non-Degenerate Relativistic Energy Distribution — Maxwell Energy Distribution 0.4 0.3 I 0.2 0.1 0 l l 1 l T 0 0 2 0.4 0 6 0 8 1 VIC Figure 5.4: Plot of the relativistic energy distribution (blue) and non-relativistic energy distribution (red) for a gas of non-degenerate protons with temperature T = 1011 K expressed as functions of velocity. the free proton gas is seen in figure 5.4 to be well approximated as non-relativistic and this approximation will become more accurate at lower temperatures. This result 148 will hold for gases of slightly more massive neutrons and will certainly be applica- ble for gases of the 385 more massive nuclei included in our simulation. Thus only the electron gas needs to be treated with the fully relativistic formalism, however for symmetry’s sake, relativistic energy distributions are used for baryonic matter as well. 5.2.3 Neutrino Capture Tables To model the capture of neutrino test particles by the nuclei or free neutrons (A, Z) represented by matter test particles and the effects that these reactions have on the temperature distribution, it was explained in sections 3.13.3 and 3.16 that it is necessary to tabulate the average energy of an electron produced by the capture of a neutrino with selected incident LTE frame energies EV by the nucleus or free neutron (A, Z) surrounded by an electron gas with selected number densities n and temperatures T as well as the average cross section of this capture for the purpose of interpolation. This is done in the following way. At each (n, T, EV) site in the interpolation lattice, we systematically sample the thermal energy distributions of the (A, Z) particles 10 million times and model their potential captures of neutrinos with energy EV. Once the LTE frame kinetic energy of the (A, Z) particle has been generated, the magnitude of the momentum vectors of the neutrino and (A, Z) particle are determined and then the vectors are randomly oriented. These LTE frame momentum vectors are then used to calculate the energy of the electron that might be produced during the capture of the neutrino by the (A,Z) particle and the cross section of this particular capture. We consider the former task first. Just like the way the electron captures were modeled, we model neutrino captures using relativistic classical mechanics and simply boost into the C-O-M frame of the selected neutrino-(A, Z) pair and check to see if the C-O-M frame energy is large 149 enough to overcome the Q-value of the neutrino capture. If it is not, we reject the capture and move on to the next pair. If it is, we model the capture by changing the C-O-M frame energy in the way that energy conservation dictates, redefining the C-O-M frame to be that of the new electron and final state (A, Z + 1) nucleus or free proton, randomly re-orienting the new C-O-M frame momentum vectors, and then boosting the electron’s momentum into the LTE frame. The LTE frame kinetic energy of the new electron is calculated and then its contribution is added to the average value tabulated for the particular (71, T) pair in the interpolation lattice. To calculate the cross section of the capture of the neutrino being considered by the nucleus or free neutron (A, Z) being considered, a different prescription is followed. The neutrino—matter cross section that we use are taken from a very convenient com- pilation of such cross sections and rates presented by Burrows and Thompson [68]. The effective cross section of neutrino capture by free neutrons given there can be written as 2 1+39 E +A . Ugap=(1—fe(Ee))0'0 X (————4 A) X (_V 2111)) 772,36 2 2 1/2 meC X l— —— W, 5.12 (El/+Anp) A! ( ) where E, is the incident nuclear rest frame neutrino energy Anp = mnc2 — mpc2, E3 is the final state electron nuclear rest frame energy, and fe(Ee) is the electron gas Fermi-Dirac function evaulated at E5. The scale cross section 00 is given by 402 (mecz) 2 7r (66)4 00 = 2 1.705 x 10_20 barns (5.13) 150 the axial-vector coupling constant 9 A is given by 9A 2 —1.26 (5.14) and WM is the weak magnetic recoil correction and is approximately given by WM = 1+ 1.1- Ey/mnc2 [148]. To calculate the cross section of neutrino capture by the nucleus (A, Z), detailed balance is used and it is assumed that electron captures are dominated by the 1 f7/2 ——> 1f5/2 Gammow—Teller resonance [149]. Without the nuclear phase space blocking term, the effective capture cross section is given by Burrows and Thompson to be 2 2 1/2 cap- _ 902 M _ _Tlei2_ 0.4 —(1 fe(Ee))C 49A( .6662 ) 1 (El/+6?) (5.15) where E; is the nuclear rest frame neutrino energy, 62’ = MECQ — M A’C2 E [VI/102 — M A,c2 +A, with A’ symbolizing the final state nucleus and A symbolizing the energy of the neutron 1f5/2 state, and C is an estimate of the relevant nuclear spin sums. For “ballpark” figures for A and C, we use those generated for neutron-rich nuclei with Z > 20 and N < 40. In this case, A ~ 3 MeV [150] and using zero-order shell model, C is roughly approximated by a simple piece-wise function of N and Z [68, 129]. c z ng(Z)Nh(N) (5.16) where Np(Z) is an the estimate of the number of protons in the single-particle 1f7/2 level and N h(N ) is an estimate of the number of neutron holes in the single—particle 1f5/2 level. The zero-order shell model predictions for Np(Z) and N 71“" ) are given 151 by [129] Z — 20, 20 < Z < 28 8, Z > 28 _ 6, N < 34 Nh(N) = (5.18) 40—N, 34 8T (7 6) The above requirement that the electron gas pressure be independent of tempera- ture was compatible with our previous works that focused on the early stages of the collapse phase when the electron gas can still be reasonably well approximated as de- generate [107, 111, 112, 117, 154]. When working with finite temperature electron gas statistical mechanics, as we do now, it is not possible to construct an average electron potential from which one can derive the average force acting on an electron in an electron gas element by the electrons surrounding it. Therefore one cannot sensibly speak of an average matter test particle electron potential energy unless the electron gas is assumed to be perfectly degenerate so that condition (7.6) is automatically satisfied. That is not to say that there is no way of monitoring the conservation of energy by our latest version of the code which works with finite temperature statistical mechanics even if we wish to avoid considerations of matter test particle shapes and sizes. We could switch to a volumetric formula for the matter test particle nucleonic 169 potential and then divide the core into volumes and use volumetric hydrodynamic-like conservation laws to test the energy conservation of the local forces. However we are reluctant to do this as it is a step away from our test particle approach and have not implemented anything like it yet although it appears to be inevitable. Instead we defer to test of the conservation of energy by local forces that assume the electron gas is perfectly degenerate to gain insight into the finite temperature code’s ability to conserve energy. We feel that this is reasonable temporary fix since the accuracy with which the local forces conserve energy boils down to our ability to determine the local value of statistical distributions and their derivatives. Since the gradient of the local degenerate electron gas pressure will depend on the local values of the electron number density and its gradient, we will obtain a measure of how well the code can determine the local value of distributions and their derivatives and we can do so using matter test particle potentials. Before delving into such considerations, it is necessary to derive a convenient analytic formula for the pressure exerted by a degenerate electron gas. Degenerate Electron Gas Pressure In the degenerate limit, the pressure integral given in (4.11) simplifies to the analytic form [113] :9 Pdeg= 3 8,13/0 F p4 “p (7.7) m \/1+ ()p/mc2 where p F is the Fermi momentum in the degenerate electron gas element and m is the electron mass. Introducing the dimensionless variable a: = p/mc yields P _ m4c5 /$F x4dx deg 352113 0 \/1+:r2 170 The above dimensionless integral is easily solved with a hyperbolic substitution. How- ever, in order to gain further physical insight into the system that will aid in our attempt determine the average potential interaction energy uave(n), we proceed in the following less than intuitive way. Via integration by part, we may write :1: 4 x ,/ f F ___:1: d5“ =x%‘/1+x%—3/ F232 1+a:2d:r (7.9) 0 \/l+:1:2 0 After plugging this result into equation (7.8), we arrive at 4 5 m4 C5 :1: _ m c 3 2_ F 2:2 /——— Now let us physically interpret the two terms above. Turning out attention to the first of these, we note that since the electron number density is given by n = pe/m, where pg is the mass density of the degenerate electron gas element, the definition of :1: F can be written as . = ___’i(3“2>1/3(P_e)1/3 F mc m 3 2 3 _ h 37r _p_e Inserting the above result into the first term on the right hand side of equation (7.10) and recognizing 1 + 5%. as y F gives us m4c5 3 / 2 2 gn—7r—c—2h3xF 1+:rF: prec (7.12) which is simply the total energy density the electron gas element would have if all of the electrons in it moved with the local Fermi momentum. To shed light on the physical meaning of the second term on the right hand side of equation (7. 10), consider 171 the following. If the electron gas element is perfectly degenerate, none of the electrons it contains have a momentum greater then the local Fermi momentum. Thus if the Pauli Exclu- sion Principle is taken into account, the total energy of the electrons in the element, U total, can be written as UtOtal = 2 Z (/m2c4+p2c2 (7.13) lflSpF where the above sum runs over all allowed momentum states with magnitudes no larger than the local Fermi momentum p F' In the statistical limit [146], the above summation becomes V/h3 times an integral over a continuous distribution. In spher- ical coordinates, after integrating over all angles, the above expression becomes P , UtOtal = 8772-1) F dpp2(/m2c4 + pgc2 (7.14) Again defining a: = p/mc and substituting it into the above integral gives 4 5 a: U10tal =Vm—9- F dxmz 1+3:2 (7.15) 7r2h3 0 A comparison of the above expression to the second term in the right hand side of (7.10) leads us to the conclusion that this term is none other than the total energy density of the degenerate electron gas element U‘t/Otal. Thus, according to the above conclusion and equation ( 7.12), the degenerate elec- tron gas pressure given by equation (7.10) reduces to Pdeg = 7Fp302 — Uf/Otal (7.16) which is merely the energy per unit volume required to promote all of the electrons 172 in the element of degenerate electron gas to the local Fermi surface. With this simple analytic expression for Pdeg’ we can now proceed with the evaluation of its gradient to determine the average force, Fave, acting on an electron in an element of degenerate electron gas by those surrounding it. Degenerate Electron Gas Forces Using the chain rule, the gradient of the pressure is given by BPde —-o _ g —o VPdeg — 3.731: VICF _ 8 . 2 total " For organizational purposes, we evaluate the two terms in the parentheses separately first and then we address the gradient of :1: F' Doing so, we write the first term in the parentheses of equation (7.17) as 5 2 _ 2.6 2 E(vppec) — mc 3$F( 1+a:F n) = mc2——$—F—n + mc2 1+ :52 an 2 F.6‘x ‘/1+a:F F Again recalling the definition of 231:, we have (7.18) CL‘F = —h—(3772n)1/3 _,,, = (7393:3351 (7.19) Therefore the partial derivative of n with respect to :r F is simply 0n _ 3(fl9)3 33% 8:13p _ h 3W2 = 31 (7.20) xF Plugging equation (7.20) into equation (7.18) and replacing the square roots with 'y F gives 5 2 2W 2 n _ = _ 3— amp (7Fpec ) mc 7Fn+mc 7F 53F 2 $17 3 = 7ch n —2— + -— (7.21) 717 33F Now consider the second term in the parentheses of equation (7.17). aUtotal 4 5 a: V = me 3 Fd222V1+$2 3331:) ”2&3 83F 0 4 5 _ TL 2 2 Regrouping and replacing the square root with ”y F leads to aUtotal 3 3 V _ A 2 m C 2 With the right hand side of equation (7.20) in mind, the term in the parentheses above is identified as 3n /:1:F. Thus the above expression simplifies to aUtOtal 2 3 = — .24 32F 127ch 17F (7 ) 174 Substituting equations (7.21) and (7.24) into equation (7.17) gives ~ _ 5 2 total ~, VP — m (VFpeC _UV )VIF = nvpmc2$§ {7.2T (7.25) 7F Finally, we turn our attention to the gradient of 23F. 1 3 ~ h(3"2) / ~ 1 3 Vsz = ————V[n /] mc 1/3 h 37r2 = ( ) 1n‘2/3vn me 1 3 1h(322n) / _‘ = — —Vn 3 me n — 435% (7 26) _ 3 n . Combining this result with equation (7.25) gives .. 1 2137671 7F Note that 2 2 2 2 2 2 (I? ’U 717 C 175 Plugging this into the expression for fipdeg given in equation (7.27) yields _. 1 s vpdeg = gvpmv%Vn 1 _. = EvaFVn (7.29) Finally, with this result, we are poised to generate an expression for the average force acting on an electron in an element of the degenerate electron gas by the surrounding electrons. 1 Fave : —; VPdeg 1 6n = —— — 7. 329va n ( 30) Degenerate Electron Gas Average Potential Energy While it is assumed that the degenerate electron gas is ideal and therefore non- interacting, the fact that it obeys the laws of quantum statistics means that there is a quantum statistical interaction between the electrons that allows them to influ- ence each others’ kinetic energies. To determine this energy, consider the following arguments. If .. 1 .. 2 Fave = —;VPdeg E —VUa/()e (7.31) it must be that for any i = 1,2,3 B-uave : lapdeg (7 32) (9.2:,- n 8:5,; ' 176 Therefore, by equation (7.30), we arrive at auave _ 119171111727: 8:13,,- 3 n 83:,- (733) In the degenerate limit, Uave is a function of n only. Therefore the chain rule tells us that auave : EPFUF 8n 3 12 TI. (9115122601,) / 1 [TI pF(n,)’UF(Tl,) I _, -—————d = _ d 7.34 [0 8n’ n 3 0 n, n ( l where the upper limit of integration, n, is the local electron number density. Carrying out the simple integration on the left hand side of the above and expressing vF(n’) in terms of p F (n, ) in the right hand side yields 2 —1/2 1+ (M) dn’ (7.35) 193:0“) ’ me 1 n Uavefn) = 3772/0 n where we took uave(0) = 0 since 77. = 0 implies that there is no electron gas. Eir- thermore pF(n’) = h(37r2n’)1/3 I 1 3 —>’n : _— 3772h3pF —->an = ——1———p%dpF (7.36) 22h3 177 The insertion of the the above expressions for n, and dn’ into equation (7.35) allows that integral to be expressed solely in terms of Fermi momentum. 1 PF 19' 2 _ Uave(n) = ‘TE/O Pl? 1+('m_1:) dP’F = {(1+ ((7%)?)1/2 —1}mc2 = (,F _ 1) mc2 (7.37) Thus the average quantum statistical interaction potential energy that satisfies F ave = —6uave, with Fave = —1/n- 6Pdeg is simply the kinetic energy of an electron that moves with the local Fermi momentum. Since we are free add to this potential energy 2 an arbitrary constant, for convenience, we add me and define average interaction po- tential energy of an electron in an element of degenerate electron gas uave(n) to be the total energy of an electron moving with the local Fermi momentum uave(n) E ’7ch2 (7.38) Thus we can define an average matter test particle electron potential energy Uave simply by evaluating Have at each matter test particle’s centroid coordinates and multiplying the result by the number of electrons it implicitly represent. Testing the Energy Conservation of Local Forces Armed with a general expression for Have in the degenerate electron gas limit, we can now test the ability of local forces to conserve energy. The toy model constructed to do this only simulates the local forces generated by the local mean field nucleonic potential and electron gas pressure. We choose to we use static initial conditions so 178 the matter test particles initially only have local interaction potential energy and ac- quire kinetic energy as they are accelerated by the effective interaction potential field. Since this the matter test particles will repel one another as they interact through local forces in this toy model, at some point the roughly 100 data points we have to interpolate distributions and their radial derivatives off of will start to become widely separated. Since the data points must be sufficiently close to one another to accurately resolve radial changes in distributions, the fact that the test particles are moving away from one anther will eventually lead to problems. Thus we only follow this toy model for 0.05 5, during which the core significantly expands, and do not expect energy to be conserved perfectly. In figure 7.3, we present a plot of the matter test particle kinetic energy (red), local interaction potential energy (green), and total energy (blue) as functions of time. It is clear from figure 7.3 that this algorithm con- 1.E+45 1.E+44 TE L'- 1.E+42 >- E / 2 1.E+41 Ill 1.E+40 —Total Energy 1.E+39 -— Kinetic Energy -- Local Interaction Potential Energy 1.E+38 . 1 ' ‘ o 0.01 0.02 0.03 0.04 0.05 ‘ Time [s] Figure 7.3: Plot of the matter test particle kinetic, local interaction potential, and total energy for a toy model that only includes local forces. 179 serves energy reasonably well. For this toy model, violations of energy conservation 4th order uncertainties with the Runge—Kutta integrator are not just generated by the we use to solve the matter test particle equations of motion and computer roundoff errors. Limitations on our ability to calculate statistical distributions and their gra- dients come into play as well. As the core expands, the total energy decreases by approximately 7%. This is attributable to the expected statistical issues associated with large radial separation between data points mentioned above. This does increase our confidence that we are capable of calculating the values of distributions and their radial derivatives with fair accuracy at arbitrary points. Although this is admittedly not the a very tough test. In full calculations, as the core collapses, distributions will change much more rapidly radially and it will become more important to accurately calculate the local values of distributions and their radial derivatives. One could propose using a toy model that includes gravitation which would cause the core to contract and keep the data points from getting too far apart from one another. However in the absence of weak reactions, which cannot sensibly be included in a model that assumes the electron gas remains perfectly degenerate at all times, other statistical problems are expected to arise. To see that this is the case, consider the following. Single processor simulations usually have about 100 data points spread over the radius of the core at which statistical distributions are known. This means that if any of the statistical distributions begin to change too rapidly radially before the core is contracted enough so that the data points at the average radii of the spherical shells are sufficiently close enough together to resolve these radial changes, we cannot safely interpolate between the data points. In full calculations where weak reactions are modeled, electrons captures reduce the electron gas’s ability to resist the inward crush of gravity and the collapse rapidly accelerates. Consequently none of the statistical distributions radially change too much before the core is contracted enough to resolve 180 their radial behavior. Such is not the case for the proposed toy model. The lack of electron captures would allow the electron gas to repel the inward crush of gravity after the core has only slightly contracted. This would then generates ripples in the density distribution that would occur when the data points are not close enough together to resolve them. This was found to be the case in some preliminary calculations. In these studies, many matter test particles were subject to unphysical forces and energy was not conserved. In essence, this toy model would be less of a test of energy conservation as it would be a test of how far we must push the code in a direction in which we expect it to perform poorly before energy conservation completely breaks down. Thus we opt to save this proposed toy model for calculations performed after the code has been parallelized and converted to 64-bit so 107 or even 108 matter test particles could be used to construct many more more spherical shells all containing substantially more matter test particles than we currently do. These calculations would be capable of easily resolving any density ripples that would be formed as the electron degeneracy pressure begins to halt the collapse in toy models with no weak reactions. 181 Chapter 8 Numerical Results In this chapter, we discuss the results generated by several calculations. Simulations were performed using all three nucleon potentials described in section 3.10.2, the four different modified version of the FF N electron capture rate tables discussed in section 5.1, and various numbers of spherical shells that divided the core into volumes for the purpose of calculating distributions. The motivation for systematically varying the nucleon potentials and electron capture rates is to study any physical dependencies the collapse and explosion mechanisms may have upon them. The purpose of sys- tematically varying the number of spherical shells is to test for unwanted numerical dependence upon the number of spherical shells used to calculate statistical distribu- tions. Additionally some preliminary calculations that model fusion were performed in each series. The inclusion of the fusion calculations is merely confirm that the subroutines designed to model this are functioning properly since the nuclear inter- action network is incomplete at this time and no realistic data can be generated by simulations that make use of it and they are discussed no further. Other than these differences, all of the simulations discussed here had identical inputs. All simulations were performed on a single processor, made use of the 106 matter test particles, and the same initial conditions. All of the results presented here were generated by sim- 182 ulations of the collapse of a non-rotating 1.33 ME) core with the initial density and temperature distributions shown in figure 3.2, initial chemical composition given in figure 3.4, and the initial collapse profile described in section 3.5.1. Before we proceed into a discussion about the distinctions between results obtained by the different simulations, we first address what these simulations have in common. This is done because irrespective of which nucleon potential or electron capture rate table was utilized or how many spherical shells the simulations made use of, a similar phenomenon that profoundly impacted core dynamics was seen emerging in all of them. The resultant dynamics may prove to be the first glimpse of a completely new explosion mechanism that is entirely different from the accepted picture. Hence an in depth description of this phenomenon, what precipitates it, and how it effects the core’s evolution is warranted. Once this description is complete, we will compare the results of different simulations. To aid in our description of the new dynamics observed in all of our simulations, while making general comments applicable to all simulations presented in this thesis, we specifically reference the output of one particular calculation. This calculation made use of the soft BKD nucleon potential and the linearly extrapolated FFN elec- tron capture rate table with all entries reduced by an order of magnitude discussed in sections 3.10.2 and 5.1 respectively. We chose to reference this calculation for two main reasons. One is that isospin independent nucleon potentials are generally better understood than isospin dependent nucleon potentials for the reasons explained in section 3.10.2. It turns out that of the two isospin independent nucleon potentials we use, the phenomenon we have repeatedly mentioned is slightly more pronounced in the soft potential case and thus we chose to make an example out of calculation that used it. The second reason we elected to reference this calculation is that the electron capture rate table it used is probably the most realistic of the four different electron capture rate tables we generated for the reasons explained in section 5.1. 183 8.1 New Dynamics Observed In all of the calculations we performed, the first deviation form the accepted picture observed in statistical distributions was found in the electron fraction. Therefore we commence our description of the new dynamics seen in our simulations with an account of how the electron fraction distribution evolves in time. After we explore this topic and the impact that these electron fraction configurations have on the overall collapse and explosion mechanics, we examine the genesis of the resultant dynamics in neutrino—matter interactions and nuclear structure. 8.1.1 Evolution of the Electron Fraction Initially the electron fraction distribution was always found to change in time in such a manner that is consistent with the accepted picture of a core collapse supernova. Electron capture rates initially slowly increased as the collapse progressed and at later times rapidly increased and eventually led to rapid deleptonization of the inner core. This can be seen in the electron fraction plots from our reference calculation in figures 8.1, 8.2, and 8.3. In these plots each red dot symbolizes the electron fraction in a spherical shell. Not long after the inner region of the core begins to rapidly deleptonize, an unex- pected phenomenon is observed in each calculation. In the accepted picture of core collapse and bounce, this rapid deleptonization continues throughout the entire inner region of the core further reducing its primary pressure source thus allowing the rest of the core to easily collapse into it and generate a bounce through the canonical mechanism described in section 2.2. Contrary to this notion, all of our calculations found that neutrino captures at intermediate radii inside the inner region of the core deposited a large amount of electrons in a narrow spherical shell centered about a radius of approximately 50—60 km. The electron fraction in this region became in- 184 creasingly elevated as did the pressure of the electron gas. Eventually a well defined narrow spike in the electron fraction always formed. This electron fraction spike is seen in various stages of development in figures 8.4, 8.5, and 8.6. The significance of the spike in the electron fraction distribution can be understood in the following way. It is seen in figures 8.4, 8.5, and 8.6 that the spike is sufficiently narrow that the magnitude of the radial derivative of the electron fraction almost everywhere in it becomes large compared to other locations in the core. Thus the radial derivative of the electron fraction in the spike at radii greater than the radius of its peak can become extremely large and negative. Conversely the radial derivative of the electron fraction in the spike at radii less than the radius of its peak can become extremely large and positive. This coupled with the a presumably large negative radial derivative of the density distribution in this region, demonstrated to be true in section 8.1.3, gives rise to two important features in the electron number density distribution. One is an much steeper drop off in electron number density at radii just outside the radius of the spike’s peak than would otherwise occur if the spike was absent. The other is a less steep drop off in electron number density at radii just inside the radius of the Spike’s peak than would otherwise occur if the spike was absent. The former results in a large increase in the outward pressure exerted by the electron gas present at radii just outside the electron fraction Spike’s peak radius that works against the inward collapse of the outlying infalling matter. The latter results in a reduction in the pressure exerted by the electron gas that would otherwise resist the tendency of the matter inside the electron fraction spike’s peak radius to collapse under its own weight. To see the extent to which these unexpected changes in the pressure profile of the electron gas impact core dynamics, we next study the radial component of the average 3 of the matter in the spherical shells. 185 1 J . shell electron fraction .3- .6- .4_ ’ °'°._________._......--',..uuuu. .2. 00 250 5'00 7'50 1000 r [km] t = 0.00000 3 Figure 8.1: Plot of the initial electron fraction in the spherical shells. 186 . shell electron fraction 0 2'50 500 7'50 1600 r [km] t = 0.09700 3 Figure 8.2: Plot of the electron fraction in the spherical shells shortly after the inner region of the core’s rate of deleptonization began to increase. 187 14 . shell electron fraction .8u .6- 4 ./ .24 0 r . . . 0 250 500 750 1000 r [km] t = 0.10860 3 Figure 8.3: Plot of the electron fraction in the spherical shells shortly after the inner region of the core has begun to rapidly deleptonize. 188 1i . shell electron fraction .8- .6- .41,” .2- 0 . . . . 0 250 500 750 1000 r[km] t=0.112653 Figure 8.4: Plot of the electron fraction in the spherical shells shortly after the spike develops. 189 l 1 . shell electron fraction 0 2'50 500 750 1000 r[km] t=0.119058 Figure 8.5: Plot of the electron fraction in the spherical shells showing fully relep- tonized matter in the spike having an electron fraction approximately equal to matter in the outermost shells where the effects of electron capture are negligible. 190 .shell electron fraction 0 250 500 750 1000 r [km] t = 0.12400 3 Figure 8.6: Plot of the electron fraction in the spherical shells showing hyper- leptonized matter in the spike having an electron fraction greater than matter had at any point in the core before the collapse began. 191 8.1.2 A Possible New Explosion Mechanism In this section we present and analyze a series of plots that show both the electron fractions and the radial component of the average 5 of matter, denoted by 57-, in the spherical shells. This is done to highlight the dependence the radial motion of matter has on the electron fraction configuration. To ensure that the details of the Br distribution can be resolved, it is plotted on a base 10 logarithmic scale. This precludes us from displaying negative values. Therefore color is used to differentiate between negative Er’s, which indicate inward radial motion and positive Br’s, which indicate outward radial motion. The former is symbolized by a red dot and the latter by a green dot. In these plots each magenta dot symbolizes the electron fraction in a shell. These plots are taken from our reference calculation at the same times the previous series of plots were in the last section. As such the first three plots presented in figures 8.7, 8.8, and 8.9 depict intuitively expected Br distributions as no unexpected features have yet developed in the electron fraction distribution. Initially the 37‘ distribution is determined by the initial collapse profile. At later times, all of our simulations found that inward acceleration was seen everywhere in the core and eventually became the largest for matter at intermediate radii. Matter at this distance from the origin is close enough to experience a larger gravitational acceleration than matter on the outer region of the core, but it is far enough from the origin that its inward motion is not impeded by particle collisions that occur in the high density inner region of the core. The last three plots of the electron fraction and Br distributions require more in depth considerations as the presence of the electron fraction spike complicates the 31' distribution in every simulation. In figure 8.10 shows the electron fraction and Br distributions shortly after the spike develops. Two important statements must be made about this plot. One is that matter is clearly decelerated in the region between 192 the electron fraction Spike’s peak and the outer edge of the rapidly deleptonized inner region of the core where the the electron fraction distribution is depressed. The other is that matter is seen beginning to move outward in this plot as indicate by the green dot. In our reference calculation, after this time step an increasing number of spherical shells were found to have positive Br’s that increased in magnitude and a well defined outward moving density wave was formed. Therefore the significance of the green dot in figure 8.10 is that it is the radius at which the outward explosion of matter begins to occur. In our reference calculations the radius at which the explosion began was approximately 58 km. Similar results were obtained by the rest of the simulations. This is a significant deviation from the accepted picture. Instead of the outward explosion of matter being generated by the pressure exerted by matter at supernuclear densities accumulated in a small volume near the origin, all of our calculations found that the an outward explosion of matter formed occurs at radii greater than 50 km and is the result of a build up of electrons produced during neutrino captures that significantly increases the outward electron gas pressure. The propagation of the density wave formed by the neutrino capture driven explosion can be seen in figures 8.11 and 8.12. In these figures, it is seen that the ejecta acquires velocities in the neighborhood of one tenth the speed of light. The proto—remnant is more difficult to discern from these plots and we therefore postpone in depth discussions about it until the next section where density plots are considered. It suffices to say here that while the proto—remnant is tightly gravitationally bound, violent collisions between energetic matter test particles can cause the value of Br in a some of the shells that comprise the proto—remnant to become momentarily positive. Thus examining a single plot of the 37' distribution can make the identification of the proto—remnant perplexing. Only when all of the 37‘ distribution plots are considered, particularly in the form of a movie, can the outer edge of the proto-remnant be distinguished from the ejecta. 193 oshell electron fraction 4 . ... I”..oooooooo o l l l I 0 250 500 750 1 000 r [km] 1 0 - . , - posmve shell beta 1 onegative shell beta 10 -1 10-2 1 O -3" o o ‘ 1o -4 10 -5- 10 -6. 10-7 . r . . 0 250 500 750 1 000 r [km] 1 = 0.00000 s Figure 8.7: Plot. of the initial electron fraction and Br of matter in the spherical shells. 194 oshell electron fraction .4 - CM... .2. O l T 0 250 560 750 1000 r [km] 10 . . , vposmve shell beta 1 onegative shell beta 10 -1- 10-2 K . 10-34 10-44 1 0 -5. 1 0-6. 10 -7 . . . . 0 250 500 750 1000 r [km] t = 0.09700 8 Figure 8.8: Plot of the electron fraction and Br of matter in the spherical shells shortly after the inner region of the core’s rate of deleptonization began to increase. 195 oshell electron fraction 10-6- 0 250 500 7550 1000 r [km] opositive shell beta onegative shell beta 10-7 0 2'50 €00 7E0 1000 r [km] t = 0.10860 3 Figure 8.9: Plot of the electron fraction and 37‘ of matter in the spherical shells shortly after the inner region of the core has begun to rapidly deleptonize. 196 oshell electron fraction .2 l 0 l 1 r l 0 250 500 750 1000 r [km] 10. . . oposmve shell beta 1 onegative shell beta 10'1- A... . , . 10- E .1. 10'3~ 2 10-4- ’ 10-5 10-6- 10-7 . . . . 0 250 500 750 1000 r[km] t=0.11265S Figure 8.10: Plot of the electron fraction and Br of matter in the spherical shells at the time the outward explosion of matter forms shortly after the spike develops and the matter between the Spike’s peak and the outer edge of the rapidly deleptonized inner region of the core where the electron fraction is depressed has clearly decelerated. 197 ,5 ~shell electron fraction 0 250 500 750 1000 r [km] ~positive shell beta onegative shell beta lo-t 10 '1 10-3 10*4 1042' 104i 10=7 . .4 . 0 250 500 750 1000 r[km] t=0.11905 s Figure 8.11: Plot of the electron fraction and Br of matter in the spherical shells showing fully releptonized matter in the spike having an electron fraction approxi- mately equal to matter in the outermost shells where the effects of electron capture are negligible and a large outward moving density wave is visible. 198 6 shall electron fraction ; O O 0 O O .4 .2 O l I I F T 0 250 500 750 1 000 r [km] 10 . , , ~posulve shell beta 1 onegative shell beta 1 0 -1- Afigwarma' ' . alga? 1 0 -2- fiaflfif O .. 10-3 10- 10- * 10- 10- 0 Figure 8.12: Plot of the electron fraction and Br of matter in the spherical shells showing hyper-leptonized matter in the spike having an electron fraction greater than matter had at any point in the core before the collapse began and the density wave has successfully propagated almost completely out of the core. 250 500 7'50 1000 r [km] t = 0.12400 3 199 This neutrino-mater interaction driven explosion mechanism appears to be com- pletely consistent with the physics that have been included in our model so far. Not only is it clear what. powers it, its origins in neutrino-matter interactions, nuclear structure, and quantum statistical mechanics discussed in section 8.2 are well under- stood, and it has survived every numerical convergence test it has been subject to. More convergence tests must conducted before we can be completely confident that the mechanism is entirely physically consistent with the physical processes currently being modeled. This subject is discussed further in section 8.5. Additionally, more physical processes must be built into our model before we can safely conclude that this explosion mechanism is completely physically realistic. This subject is discussed further in chapter 9. 8. 1.3 Explosion Dynamics To gain further insight into the dynamics of the new explosion mechanism our calcu- lations are consistently yielding thus far, we now study a series of plots that show the density distribution as well as the electron fraction and Br distributions. This shall allow us to examine the role that density plays in this scenario and how it quanti- tatively differs from the accepted picture of the bounce and explosion. The density is plotted on the same base ten logarithmic scale as the 31" distribution is in units of nuclear density p0. Other than the appearance of the black dots that symbolize the density in the shells, these plots are identical to those that appeared in the last section. The first three plots are taken at times before the electron fraction spike develops and as such no phenomena that are completely at odds with the canonical story of this stage of the collapse are observed. The first plot seen in figure 8.13 again simply shows the initial conditions of the Woosley and Weaver progenitor. The analyses of the density distributions in the next two figures are more complicated and are 200 presented separately below. The second plot given in figure 8.14 shows an expected increase in the density in most of the core leveling off at the origin. However it also shows some fluctuations that are due to the low densities in the outer edge of the core. These low density fluctuations were found in all calculations. Since none of the important processes that govern the core’s evolution take place in the outermost region of the core, we need not concern ourselves with these fluctuations as the high density regions of the core where our statistics are much better truly drive the dynamics. We do note, however, that the magnitude of these fluctuations will be reduced when the number of matter test particles is increased. They appear in all plots to varying extents and we shall comment on them no further. The third plot displayed in figure 8.15 an expected increase in the density of shells with intermediate radii and a slight density inversion in the innermost region of the core. This density inversion is seen at this stage in the collapse in all calculations. Thus far its origins appear to be physically consistent with the density and electron fraction distributions and the resultant forces generated by the gravitational field and electron pressure. It turns out that the interplay between these two forces gives rise to an inward acceleration profile that is slightly larger at radii in a small range outside the outer radius of the density inversion zone. This result should not come as to great of a shock since the gravitational acceleration profile of a spherically symmetric star peaks at some distance from the origin that depends on its density profile and the core is no longer in hydrostatic equilibrium. Not only has this result been observed by simulations with all of the physical parameters systematically varied, it has survived a battery of preliminary numerical convergence test. It has been observed in simulations that use 106 matter test particles and 50, 75, and 100 spherical shells to calculate statistical distributions. Additional convergence test using substantially more than 106 matter test particles and more than 100 spherical shells must be run as well, 201 however these test will have to wait until the code has been converted to 64—bit and parallelized. The fourth plot displayed in figure 8.16 corresponds to the time the neutrino- matter interaction driven explosion begins. Several things must be noted about the density distribution at this time. The density is now steeply increasing for radii less than 125 km. Over this range, the density increases by roughly three orders of magnitude. The interplay between the forces generated by the gravitational field and electron pressure no longer result in an density inversion at small radii. A close examination of the innermost density data points reveals that the innermost two shells have a lower density than the other shells. This is due to the fact that these shells contain half as many matter test particles as the rest and as such are twice as susceptible to statistical fluctuations. These two shells are really one shell that was cut in half in an attempt to generate a data point as close as possible to the origin in order to minimize substantial numerical artifacts caused by a lack of data over too large of a distance from the origin. When the physical density inversion vanishes, these two shells lag slightly behind. Their densities are quickly driven up to the intuitively expected values by the dynamics in the shells with better statistics. Furthermore since this explosion mechanism is not directly sensitive to the value of the density or its gradient near the origin, we have no reason to suspect that these momentary fluctuations impact the explosion dynamics in an appreciable way. Since this is the plot of the density distribution at the time of explosion, it is im- portant to point out that nowhere in the core have we achieved supernuclear densities. This is yet another way in which this neutrino capture driven explosion mechanism significantly differs from the accepted bounce and explosion picture. In fact the high- est density at this point in our reference calculation only slightly more than 0.36 p0 and more importantly at the radius at which the outward explosion of matter forms the density is on the order of 10‘3 p0. All of or calculations yield results similar 202 to this. The maximum density in the core when the explosion begins is always well below nuclear matter density and the density at the radius where it is always on the order of 10-3 p0. The fourth and fifth plots presented in figures 8.17 and 8.18 are very similar. Except for the outermost region of the core where low density fluctuations occur, the density in both plots is now a strictly decreasing function of radius and the central densities are both marginally greater than nuclear matter density. The division between the proto-remnant and the is somewhat clear in the fifth plot and very clear in the sixth plot. Since the shells used to calculate statistical distributions are defined by matter test particle occupancy, as the ejecta continues to separate from the proto- remnant, the shell “between” them ends up with a very low density compared to its two neighboring shells. This depression in the density in this region is physical. however having only the one data point to describe it is insufficient. Only when many more matter test particles and spherical shells are used can this part of the density distribution be resolved. Thus we must terminate calculations before this occurs. This is the statistical barrier 32-bit single processor simulations encounter that requires the simulation to be terminated that was alluded to in section 6.1. When all of the plots presented in this section are considered, it is clear that the presence of the spike in the electron fraction has the effect dividing into a well defined high density proto-remnant and outward moving ejecta in the fashion described in section 8.1.1. At the time our reference calculation was stopped, 0.124 s after the collapse began, the mass of the proto-remnant was approximately 0.253 MG) and it’s radius was roughly 7.4 km. Its central density, electron fraction and temperature were 1.653 p0, 0.285, and 2.01 x 1011 K respectively. Results similar to these were obtained by all simulation and are discussed and compared in sections 8.3 and 8.4. 203 oshell electron fraction . ... ”..ooocoooo .4 ~ W""""'" .2 . 0 . . . . 0 250 500 750 1000 r [km] 10 - . . oshell denSIty [nuclear densnty] 1 - negative shell beta 1 positive shell beta 10 -1- 10-4 10-3“ .0‘ 10 -4- 10-5- 10-64 0...... 10-7 . . . . 0 250 500 750 1000 r [km] t = 0.00000 3 Figure 8.13: Plot of the initial electron fraction, Br, and density in the spherical shells. 204 6; oshell electron fraction 14“.”... .2- 0 I T I l 0 250 500 750 1000 r [km] 10 . . . oshell denSlty [nuclear denSIty] 1 onegative shell beta ‘ ~positive shell beta 10-1- 10-2 K 0 10-31 10-4. 10-5- - 0 10'6. ° 0 o 10-7 . . I r 0 250 500 750 1000 r [km] t = 0.09700 3 Figure 8.14: Plot of the electron fraction, ET, and density in the spherical shells shortly after the inner region of the core’s rate of deleptonization began to increase. 205 6. she" electron fraction Hf, . . . . .2 ~ 0 . . . . 0 250 500 750 1000 r [km] 10 . . . -she|l denSIty [nuclear denSlty] 1 . negative shell beta ' . positive shell beta 10 -1- 10- ' 10-3- 10-4- 10-5- 10'6. . O o 10 '7 . . . . 0 250 500 750 1000 r [km] t= 0.10860 3 Figure 8.15: Plot of the electron fraction, Br, and density in the spherical shells shortly after the inner region of the core has begun to rapidly deleptonize and the density inversion in the innermost region of the core has formed. 206 6; -shell electron fraction .4 4/0 0 .2 ~ 1 0 l 1 T l 0 250 500 750 1000 r [km] 10 . oshell density [nuclear density] onegative shell beta opositive shell beta A... . O 0 105- ° . 10-7 . . . . 0 250 500 750 1000 r[km] t=0.112655 Figure 8.16: Plot of the electron fraction, Br, and density in the spherical shells at the time the outward explosion of matter begins shortly after the spike develops and the matter between the Spike’s peak and the outer edge of the rapidly deleptonized inner region of the core where the electron fraction is depressed has clearly decelerated. 207 .6: ~shel| electron fraction AM" . . . . .2- 00 250 500 750 1600 r [km] oshell density [nuclear density] onegative shell beta opositive shell beta -51! °. 10 1 . . 10-6- - . 1 10-7 . . . 0 250 500 730 1000 r[km] t=0.11905s Figure 8.17: Plot of the electron fraction, Br, and density in the spherical shells show- ing fully releptonized matter in the spike having an electron fraction approximately equal to matter in the outermost shells where the effects of electron capture are neg- ligible and a large outward moving density wave is seen beginning its separation from the dense proto—remnant. 208 .6 ~shel| electron fraction ; O. O O O .4 ' .2- 0 . . . . 0 250 500 750 1000 r [km] 10 . . oshell denSIty [nuclear denSIty] 1 onegative shell beta opositive shell beta 10-1 1 10- . 10-3 10-4 10-5 10-6 - . . 10-7 . . . 0 250 500 750 1000 r [km] t = 0.12400 3 Figure 8.18: Plot of the electron fraction, Br, and density in the spherical shells showing hyper-leptonized matter in the spike having an electron fraction greater than matter had at any point in the core before the collapse began and the density wave has successfully propagated almost completely out of the core and has clearly separated from the dense proto—remnant. 209 8.2 Genesis of New Dynamics Thus far we have described how the neutrino capture driven explosion mechanism effects core dynamics, but we have not discussed the origin of this mechanism. We have not given any reason for why the build up of electrons produced by neutrino capture seems to be so concentrated in the region that it is always found to be and results in the spike in the electron fraction that drives the dynamics in the fashion described in the previous two sections. This topic is explored in this section. For the spike in the electron fraction to be created, it must be that in the narrow shell in which the spike forms neutrino capture strongly dominates electron capture. Recall that in all calculations the spike always formed in a region containing neutron- rich matter with a density on the order of 10—3 p0. The importance of this is fourfold. Obviously the farther away a nucleus is from the valley of beta stability and the closer it is to the neutron drip line, the lower its electron capture rate is. Furthermore, since it is energetically favorable for very neutron-rich nuclei to capture neutrinos and move towards the valley of stability, the neutrino capture Q-value becomes positive and large. This increases the bare capture cross section given in equation (5.15) and therefore the neutrino capture probability increases in the region. Additionally at densities on the order of 10—3 p0, none of the electron capture rates are insurmountably large and neutrino captures are not likely to be forbidden by the Pauli Exclusion Principle. This sets the stage for neutrino capture dominance. 8.2.1 Qualitative Analysis of Genesis Before we begin an analysis of the code’s output to describe the genesis of the neutrino capture explosion mechanism, let us qualitatively walk through how this scenario unfolds in our simulaitons. Recall that the region of the core comprised of neutron- rich matter with a density approximately equal to 10—3 p0 where the spike in the 210 electron fraction distribution forms is centered about a radius roughly 50 km from the origin. To its interior is a region of hotter denser matter in which the degeneracy of the electron gas suppresses neutrino capture rates and electron capture rates are strongly dominant. The energy of the neutrinos escaping this hotter denser region increases as it continues to contract as does the rate at which it emanates neutrinos. Thus the region where the electron fraction spike forms is bathed in an increasingly intense flux of increasingly energetic neutrinos. For the reasons stated in the previous paragraph, neutrino capture rates in this region quickly outpace electron capture rates and the local electron fraction rises. Not only can neutrino captures push neutron- rich nuclei back along isobars towards the valley of stability, as the energies of the neutrinos propagating through the region increase, they can push some nuclei past the valley of stability into proton-rich territory. In addition to this leptonization mechanism, it turns out free neutrons are highly abundant in this region and they capture a large amount of neutrinos as well. The heavy presence of free neutrons is the result of numerous electron captures by neutron drip line nuclei in the region. Through these two mechanisms, the region can become hyper-leptonized. As it does, it absorbs nearly all of the high energy neutrinos from the flux bathing it. This deprives outlying regions of the opportunity to become hyper-leptonized and is the reason that only the base of the outer edge of the electron fraction spike thickens after it is formed. The tallest part of the electron fraction spike remains narrow and radically alters the radial derivative of the electron number density distribution. This in turn divides the core into a proto—remnant and ejecta in the fashion described in section 8.1.1. We now undertake the task of confirming that the picture painted above is indeed an accurate description of our results. The density and the level of leptonization of matter in the core has already been established by the plots presented in the previous sections. However, the the fact that neutrino capture probabilities are significantly 211 higher in the region where the spike in the electron fraction distribution forms than they are in the region immediately to its interior due to the effects of electron gas degeneracy must be substantiated. We must also verify that the changes in the nuclear composition of the core described above regarding nuclei and free baryon distributions is accurate. Furthermore we must confirm that the description of the neutrino energy distribution given in the qualitative overview, in particular the fact that the matter in the electron fraction spike filters out all of the high energy neutrinos, are correct. Various forms of output generated by the code are used to validate the statements made in the qualitative overview. 8.2.2 Neutrino Capture Probability Distribution To confirm that neutrino capture probabilities in the region where the electron frac- tion spike forms are comparatively much higher than they are in the region interior to it where statistical constraints imposed by electron gas degeneracy suppresses neu- trino captures, we consider four plots of the neutrino-matter interaction probability distributions. These four plots show the elastic scattering and capture probability distributions in the core at times chosen to highlight the emergence of important features. Before we consider the plots of the neutrino-matter interaction probability distributions just before and after the spike in the electron electron fraction distribu- tion froms, we must consider two plots of these distributions from early times in the collapse to make sense of one of the features seen in these plots. Figure 8.19 shows the neutrino—matter interaction probability distributions before the innermost region of the core has begun to rapidly deleptonize. Not surprisingly both the elastic scattering and capture probabilities increase as the temperatures and densities increase. Figure 8.20 shows the neutrino-matter interaction probability distributions shortly after the innermost region of the core has begun to rapidly deleptonize. Here a steep drop off is seen in the capture probability at smaller radii 212 that corresponds exactly to the region that is now neutron rich. This turns out to be a nuclear structure effect. While the bare neutrino capture cross sections are larger for neutron rich nuclei because such captures are energetically advantageous, recall from equations (5.15) and (5.16) that this cross section is weighted by a zero-order shell model estimation of the relevant nuclear spin sums that approximates the number of neutron holes in the single particle 1f5/2 level. It was shown in section 5.2.3 that for all nuclei that this weighting factor is appli- cable for, it is a linearly decreasing function of neutron number N for N 2 34. Thus for neutron-rich nuclei there is a linear reduction in the capture cross section as its neutron number increases. The fact that the capture probability plot shows a linear decrease on a logarithmic scale is simply a consequence of the fact that cross sections are exponentiated in the beam attenuation arguments used to calculate interaction probabilities. Thus plotting the interaction probabilities on a logarithmic scale makes these plots linearly dependent on the cross section weighting factors. Clearly sensi- tivity to this nuclear structure dependence of the capture probability distribution can only be realistically modeled by a simulation that propagates a full ensemble of nuclei. In figure 8.21 we present a plot of the neutrino-matter interaction probability dis- tributions shortly before the electron fraction spike forms. Here it is seen that there is a rapid increase in the capture probabilities in the very innermost region of the core. This increase is due to the now large abundance of free neutrons in this region of the core produced by c0pious amounts of electron captures by neutron drip line nuclei. This is explicitly demonstrated to be the case in section 8.2.3. Finally in figure 8.22 a plot of the neutrino—matter interaction probability distributions at the time the outward explosion of matter begins is shown. This plot is rather similar to the previous plot except for a slight zigzag in the capture probability distribution seen at very small radii. The region in which the capture very steeply drops off cor- responds to a region in which the temperature increases by roughly a factor 2.5 and 213 the electron number density increases by slightly more than an order of magnitude. This is readily confirmed by examining plots density, electron fraction, and temper- ature distributions at this time. It is also easily confirmed by examining the of the table of electron gas 5 parameters (5 E —/.L/kT) generated by the code the produces the electron gas statistical mechanics tables that the electron gas goes from being somewhat degenerate to very degenerate over this radial interval. Thus the rapid decrease seen here is due to the fact that many neutrino captures are forbidden by the Pauli Exclusion Principle. Inside this region the temperature increases enough to ease these statistical constraints. Thus we have confirmed that the neutrino capture probability is found to be significantly higher in the region where the spike in the electron fraction distribution forms than it is in the region immediately to its interior due to the effects of electron gas degeneracy as stated in the qualitative overview of the neutrino capture explosion mechanism given in section 8.2.1. Here in passing we note that in figures 8.21 and 8.22 it is seen that the neutrinos in the innermost shells have an elastic scattering probability of 1. The neutrinos in this region have come into thermodynamic matter. This is the criterion we use to determine if neutrinos are in thermodynamic equilibrium with matter that was alluded to in section 5.2.4. 214 1‘ o elastic scattering probability a capture probability 10'“ 104* 103‘ 10"“ 105- . 10'6 ' . ' ‘ . 0 2.50 500 750 1000 r [km] t = 0.06000 8 Figure 8.19: Plot of the neutrino-matter interaction probabilities in the spherical shells before the innermost region of the core has begun to rapidly deleptonize. 215 1 ~ 0 elastic scattering probability 9 capture probability 10"1 ‘ 10'2“ 1 \ 10"3 I 1.251111 to. 3% .. .1 O. Q “.8? “ _ 1 4 ~ 6 a; 31:92.1 . . . 10. 6% Task . 0 10'5‘ . . . 10-6 250 I 3r T 0 500 750 1000 r [km] 1 = 0.09000 8 Figure 8.20: Plot of the neutrino-matter interaction probabilities in the spherical shells shortly after the innermost region of the core has begun to rapidly deleptonize. 216 1 1 o elastic scattering probability - capture probability 10+ 10'” ”I . .. 0 0 103‘ :51). o . .s“. o . .1, O :: o 0 104‘ '1’ ”. I“; 0 10'5‘ . . . 10‘6 ' ' ' ' 0 250 500 750 1000 r [km] t= 0.10860 8 Figure 8.21: Plot of the neutrino-matter interaction probabilities in the spherical shells after the rapid deleptonization of the innermost region of the core has led to a large abundance of free neutrons. 217 1 \ e elastic scattering probability a capture probability 0 e 10.1 1, .e o t... . 1 ) 10'24. ° 10-3 . 10'4‘ 10'5“ .6 , , 1° 0 250 500 r [km] t = 0.11265 8 Figure 8.22: Plot of the neutrino-matter interaction probabilities in the spherical shells shortly at the time the outward explosion of matter begins and the electron gas in the inner portion of the innermost region has become very degenerate. 218 8.2.3 Evolution of the Nuclear Composition To confirm that neutrino captures are creating proton-rich nuclei in the region where the spike in the electron fraction forms, we consider three plots of the core’s nuclear composition. In figures 8.23, 8.24, and 8.25, the initial nuclear composition of the core, its composition after its inner region has begun to rapidly deleptonize, but before the spike in the electron fraction has started to develop, and its composition at the time the outward explosion of matter began are presented. Before proceeding directly into the comparison of these plots, we briefly explain their meaning. In these figures, each cell that corresponds to a free baryon or nucleus with A > 1 whose presence is being modeled in the core is assigned a color based on their “test particle count”. For nuclei with A > 1 the test particle count is simply the number of matter test particles that represent a given nucleus. For free protons and neutrons, the test particle count refers to the total number of free protons and neutrons all matter test particles represent for each nucleus they represent. Therefore if a large number of matter test particles participate in drip line weak reactions, they may represent many free baryons for each nucleus they represent. Thus it is possible for the test particle count as it is defined here to greatly exceed 106 and it is for this reason that the test particle count legend assigns colors up to values of 108. Having established how to correctly interpret figures 8.23, 8.24, and 8.25, we proceed with their comparison. Figure 8.23 shows the initial composition of the Woosley and Weaver core comprised exclusively of Fe group nuclei. Figure 8.24 shows the nuclear composition of the core after its inner region has begun to rapidly deleptonize, but before the spike in the electron fraction has started to develop. Notice that in addition to many more neutrino-rich nuclei, free baryons are present in the core, as are some non-neutron-rich and proton-rich nuclei. Printouts of the text files used to store the nuclear composition of matter in the spherical shells confirm 219 that this build up of non-neutron-rich and proton-rich nuclei is concentrated in the region were the spike in the electron fraction distribution will form. This is significant because the electron fraction in this region has not changed yet. Only the dispersion of nuclei it contains has. This is another example of why it is of critical importance that the propagation of a full ensemble of nuclei is modeled. Figure 8.25 shows the nuclear composition of the core at the time the outward explosion of matter began. At this time there has been approximately an order of magnitude increase in the number of free protons is observed'as well as a noticeable increase-in the number of proton-rich nuclei. Again printouts of the text files used to store the nuclear composition of matter in the spherical shells confirm that this build up of proton-rich nuclei is concentrated in the region were the spike in the electron fraction is. Thus the changes in the distribution of nuclei described in the qualitative overview of the neutrino capture driven explosion mechanism are confirmed. Next we confirm that the region where the electron fraction spike develops contains a large abundance of free neutrons and that an appreciable fraction of these free neutrons capture neutrinos and become free protons. To do this, we consider two plots of the number of free baryons per “heavy” nucleus (A > 1). In figure 8.26, we display a plot of the number of free baryons per heavy nucleus in the core after its inner region has begun to rapidly deleptonize, but before the spike in the electron fraction has started to develop. There it is seen that there is indeed a large abundance of free neutrons in the region of the core that is rapidly deleptonizing. Notice that despite the fact that the spike in the electron fraction distribution has yet to deve10p, there is a spike in the free proton abundance at the radius where it will form. At this point, the only source of free protons in neutrino capture by free neutons. In figure 8.27, we present a plot of the number of free baryons per heavy nucleus in the core at the time the outward explosion of matter began. In this figure it is seen that the spike in the free proton abundance that precisely corresponds to the electron fraction 220 spike has increased substantially and that these gains in free protons have come at the cost of lowering the local neutron abundance. Thus the changes in the distribution of free baryons described in the qualitative overview of the neutrino capture driven explosion mechanism given in section 8.2.1 are confirmed. Test Particle Count ’ 108 Figure 8.23: Plot of the initial nuclear composition of the core. Test Particle Count " 10 8 106 ,104 Figure 8.24: Plot of the nuclear composition of the core after its inner region has begun to rapidly deleptonize, but before the spike in the electron fraction has started to develop. Test Particle Count ' 10 8 1303112658 Figure 8.25: Plot of the nuclear composition of the core at the time the outward explosion of matter begins. [‘0 [\3 CO etree protons per heavy nucleus siree neutrons per heavy nucleus “2.. 2 10"1 1 O C 2 e 10" ‘ . .e 10'3 ‘ 0000...: 10"“ 10-5 ~ .1 - 1. ~=~ . . "' P ' ' 0 250 500 750 1000 t: 0.10860 8 Figure 8.26: Plot of the number of free baryons per “heavy” nucleus (A > 1) in the core after its inner region has begun to rapidly deleptonize, but before the spike in the electron fraction has started to deve10p. 224 ofree protons per heavy nucleus .free neutrons per heavy nucleus C“ 0 10'2‘ :0 : 103 3" 10“i -~ 0 313 e 500 750 1600 r [km] t = 0.11265 3 Figure 8.27: Plot of the number of free baryons per “heavy” nucleus (A > 1) in the core at the time the outward explosion of matter began. 225 8.2.4 Neutrino Energy Distribution Finally we confirm that energetic neutrinos produced in the region interior to the electron fraction spike and are filtered out of the flux of neutrinos that passes through the spike. This is done by considering two plots of the average neutrino energy distribution propagating through the spherical shells. In figure 8.28, a plot if the average neutrino energy distribution in the core at the time the outward explosion of matter began is given. The core is obviously divided into a region containing energetic neutrinos and a region containing noticeably much less energetic neutrinos. The location of this division corresponds exactly to the radial location of the peak of the spike in the electron fraction distribution. Therefore it is clear that the energetic neutrinos are being captured as they pass through the spike. The slight bump in the average energy distribution seen exterior to the location of the electron fraction spike is due to the fact that electron captures dominate neutrino captures in this region. Since the captured electron energies are much lower in this region, the energy of the neutrinos produced is much lower too. To see that the trend of energetic neutrinos to be captured in the spike is sustained, we consider another plot of the average neutrino energy distribution in the core after a large outward moving density wave has formed and is beginning its separation from the dense proto—remnant is presented in figure 8.29. In this figure, it is seen that the average neutrino energy drops off slightly less abruptly only because the base of the outer edge of the electron fraction spike at this time has broadened in the fashion described above. The small kink in near the origin corresponds to the radius at which the ejecta is separating from the proto—remnant. Thus the energetics of the neutrinos everywhere in the core are accounted for and are found to be consistent with the qualitative overview of the neutrino capture explosion mechanism given in section 8.2.1. 226 104‘ . average neutrino energy [MeV] 10'2“ 10'4 . . 4 . 7* o 250 500 750 1000 r [km] 1: 0.11265 3 Figure 8.28: Plot of the average neutrino energy distribution in the core at the time the outward explosion of matter began. 227 104‘ 0 average neutrino energy [MeV] 10'2- 10. 4 I I I r 0 250 500 750 1000 r[km] t=0.11905s Figure 8.29: Plot of the average neutrino energy distribution in the core after a large outward moving density wave has formed and is beginning its separation from the dense proto-remnant. 228 8.3 Role of the Nucleon Potential In this section we explore how the neutrino capture explosion mechanism introduced in the previous section depends on the nucleon potential. This is done by comparing the results generated by simulations that used the three nucleon potentials described in section 3.10.2. The three simulations compared here all made use of the linearly extrapolated FFN electron capture rate table with all entries reduced by an order of magnitude described in section 5.1. Again we chose to work with simulations that made use of this electron capture rate table since it is probably the most realistic of the four different electron capture rate tables we generated for the reasons explained in section 5.1. Since the neutrino capture driven explosion mechanism does not require the pres- sure exerted by matter at supernuclear densities to generate the bounce and it in fact begins well below the densities where the nucleon potential becomes important, the nucleon potential cannot directly have a significant impact of the explosion dynam- ics. It can do so indirectly however by regulating the dynamics of the proto—remnant that begins to form shortly after the explosion. Recall that it is the flux of neutrinos emanating from the proto-remnant that supplies the explosion mechanism with the source electrons it needs drive up the spike in the electron fraction and generate the pressure profile that continues to accelerate the ejecta out of the core. The nucleon potentials used in this work determine the energy and intensity of the neutrino flux emanating from the proto—remnant by determining its density profile. Since the tab- ulated electron capture rates and average neutrino energies produced during electron captures both increase as the density increases, the denser a nucleon potential allows the proto—remnant to become, the more neutrinos it will radiate. In table 8.1 some basic information describing the core’s dynamics at .the time the outward explosion of matter began are presented for simulations that use the isospin 229 independent soft and stiff BKD potentials and the isospin dependent potential with a symmetry energy term. Here the time of explosion, the radius at it begins Temp, and the density at rap and the origin at the time the outward explosion of matter began are tabulated. I—EOS temp [S] . Temp [km] P(Texp) [Pol 10(0) [100] soft 0.1127 58.35 1.297x10—3 0.3627 stiff 0.1147 56.75 1.589x10—3 0.1466 symm 0.1189 51.74 2.529x10—3 0.1230 Table 8.1: Core characteristics at the time the outward explosion of matter formed as calculated by simulations that used different nucleon potentials. In table 8.1 it is seen that the soft BDK simulation the central density increased the fastest, it generated an explosion the fastest, and was able to generate an outward explosion of matter at the farthest distance from the origin where the density was the lowest. The stiff BKD potential and the isospin dependent potential, denoted by symm in the table, came in second and third in these regards respectively. A careful examination of figures 3.5 and 3.6 reveals that for densities in the range 0 S p/pO S 0.35 the soft BKD potential has the steepest radial derivative and is therefore the most attractive followed closely by the stiff BKD potential and less closely by the more shallow isospin dependent potential. Obviously the potential that is the most attractive at densities below 0.35 p0 will initially allow densities to accumulate the fastest in the central region of the core and radiate an intense flux of energetic neutrinos capable generating an outward explosion of matter at farther radii where the matter is more diffuse. Thus the results given in table 8.1 are to be expected. In table 8.2 some basic information describing the proto—remnants calculated by simulations using the three different nucleon potentials at the time the simulations were terminated is presented. Here the time at which the simulation was terminated and the central density, total mass, and radius of the proto—remnant at the time of 230 termination are tabulated. The central densities in all simulations assumed their final values well before the termination time. The differences between these central densities is consistent with the different densities of minimum energy discussed in section 3.10.2 and seen in figures 3.5 and 3.6. The differences between the proto- remnant masses is related to the time taken generate an outward explosion of matter. Even though explosions that took longer to form started at smaller radii, the extra time taken to form always allowed more mass to fall into region that would become part of the proto—remnant. The difference in their radii is again related to the different densities of minimum energy as well as the general differences in the shapes of the potentials. EOS tstop [3] 10(0) [100] M [fl/[O] R [km] . soft 0.1240 1.653 0.253 7.348 stiff 0.1240 1.290 0.279 7.689 symm 0.1350 0.8908 0.372 11.32 Table 8.2: Proto—remnant data generated by simulations that used different nucleon potentials. Similar comparison analyses were conducted for simulations that used the other three versions of the F FN electron capture rate table described in section 5.1. Both for brevity’s sake and to avoid redundancy, we do not present them here. It suffices to say that the results are very similar to those given above. Simulations that used different electron capture rates did in some ways yield some substantially different results, however the role the nucleon potential played in all cases was identical and the relative differences between the results produced by simulations that used the same electron capture rate table were the same as those seen in the discussion above. 231 8.4 Role of Electron Capture Rates In this section we probe the sensitivity of the neutrino capture explosion mechanism to the value of electron capture rates. This is done by comparing the results generated by simulations that used the four versions of the FFN electron capture rate table described in section 5.1. Again since the differences between the accepted picture of bounce and explosion and the neutrino capture explosion mechanism are slightly more pronounced in simulation that use the soft BKD potential, all of the simulations compared here used it. In this section we refer to the previously considered simulation that used the lin- early extrapolated FFN table with all of its entries reduced by an order of magnitude simulation 1. The simulation that used the so—called reduced linearly extrapolated FFN rate table described in section 5.1 is called simulation 2. The simulation that used the linearly extrapolated FFN rate table is called simulation 3. The simulation that used the so-called enhanced linearly extrapolated FFN rate table described in section 5.1 is called simulation 4. In table 8.3 the same basic information describing the core’s dynamics at the time the outward explosion of matter began are presented for simulations 1, 2, 3, and 4 is presented as was given in table 8.1. Again the time of explosion, the radius at it begins ramp, and the density at Temp and the origin at the time the outward explosion of matter began are tabulated. # temp [8] Temp [km] P(Texp) [P0] {9(0) [P0] 1 0.1127 58.35 1.297x10“3 0.3627 2 3 0.1100 51.33 2.124x10—3 0.2005 0.1092 53.32 1.801x10—3 0.2830 4 0.1091 55.66 1.521x10—3 0.3026 Table 8.3: Core characteristics at the time the outward explosion of matter formed as calculated by simulations that used different electron capture rates. 232 Despite the fact that simulations 2, 3, and 4 use tables that are now thought to be less realistic overestimates of electron capture rates, the results they generate are quite similar to simulation 1. The larger electron capture rates initially allowed the core to collapse more quickly more quickly in these simulations and consequently the core contracts more in all cases before the explosion is generated. At later stages in the collapse, the extrapolated electron capture rates were used with greater frequency and some divergences between these three simulation were observed. Simulation 2 had lower rates simulation 1 and consequently took longest of the three to produce enough energetic neutrinos to generate the explosion. That is also why the outward explosion of matter was formed at the smallest radius of the three. Simulation 4 had the highest extrapolated capture rates and therefore was able to generate an explosion the fastest at the greatest radius of the three. In table 8.2 some basic information describing the proto—remnants calculated by simulations 1 and 2 at the time of termination is presented. It turns out that the high values of the extrapolated electron capture rates used by simulations 3 and 4 lead to sufficiently violent separations between the proto-remnant and ejecta that the simulations had to be terminated in the earliest stages of the separation for the reasons explained in section 8.1.3. Therefore no meaningful data about the proto- remnant fit for comparison data from simulations 1 or 2 could be extracted. Thus for simulations 1 and 2 only do we tabulate the time at which the simulation was terminated and the central density, total mass, and radius of the proto—remnant at the time of termination. # tstop [s] 19(0) [100] 1” [MO] R [km] 1 0.1240 1.653 0.253 7.348 2 0.1224 1.565 0.253 7.436 Table 8.4: Proto-remnant data generated by simulations that used different electron capture rates. 233 It is seen here that despite the subtle differences in the way the cores evolved in simulations 1 and 2, the end result is quite similar. The two proto—remnants are iden- tical in mass and have only minor differences in density and volume. Comparisons like the one above for sets of simulations that used the other two nucleon potentials were also conducted. The results obtained about the dynamics at the time the outward explosion of matter began were very similar to those presented in table 8.3. The study of the proto—remnant proved to be much more difficult in these cases. Due to the more shallow but also more forceful explosions generated the stiff BKD potential, the ad- dition of the larger capture rates led to problems for all three simulations in resolving the even the earliest stages of the separation of the proto—remnant and ejecta, just like simulations 3 and 4 did in the above consideration. For the series of simulations that used the isospin dependent potential, the situation unfolded just like it did in the case above, but for different reasons. In this case proto-remnant calculated by isospin dependent potential simulation 2 was nearly identical to the proto—remnant calcu- lated by isospin dependent potential simulation 1, but the other two proto—remnants were lost due to the expected unphysical behavior of the isospin dependent nucleon potential for highly asymmetric nuclear matter described in section 3.10.2. 8.5 Robustness of Explosion Mechanism with Re- spect to Variation of Numerical Parameters Before we attempt to advance our model through the inclusion of additional physical processes, it is imperative that we demonstrate that the new dynamics observed thus far are independent of the code’s numerical parameters. In this section, we discuss the ways the neutrino capture driven explosion mechanism has been tested for depen- dence on the code’s numerical parameters has been subject to as well as those it shall be subject to in the future. Obvious candidates for numerical parameters the results 234 might depend on are the number of matter test particles used to model the core and the number of spherical shells used to calculate statistical distributions. The number of tests for dependancies on these parameters that single processor 32-bit simulations can be subject to is rather small. It has been previously mentioned that simulations that use 106 matter test particles and 100 spherical shells are considered low reso- lution. It has also been stated that we cannot increase the number of matter test particles significantly beyond 106 for these simulations and without doing that we cannot substantially increase the number of spherical shells either. Previous calcula- tions have shown that 106 is a rather firm lower limit for the number of matter test particles used irrespective of the number of spherical shells the core is divided into. We can however follow the core through the beginning of the explosion phase using less spherical shells with reasonable accuracy. The agreement will not be perfect since reducing the number of shells reduces the accuracy of our calculations, but “ballpark” similarities between them can be expected. The same statistical difficulties that 100 shell simulations encounter when the attempting to calculate the density distribution when the ejecta separates from the proto—remnant; they will just manifest themselves much more quickly when substantially less than 100 shells are used. Thus we compare the output of our reference calculation when 100, 75, and 50 spherical shells are used to calculate statistical distributions. In table 8.5 the same basic information describing the core’s dynamics at the time the outward explosion of matter began are presented for simulations that use 100, 75, and 50 spherical shells as was given in table 8.1 and 8.3. Again the time of explosion, the radius at it begins ramp, and the density at Temp and the origin at the time the outward explosion of matter began are tabulated. The agreement between the calculations is not perfect, but the neutrino capture driven explosion mechanism still worked and its physical description given in section 8.2 is valid in all cases. Again, since we know that 50 and 75 shell simulations 235 shell # temp [s] Temp [km] p(rea:p) [pol 10(0) [100] 100 0.1127 58.35 1.297x 10—3 0.3627 75 0.1147 53.67 1.501x10"3 0.2037 50 0.1189 61.12 1.407x10’3 0.2805 Table 8.5: Core characteristics at the time the outward explosion of matter formed as calculated by simulations that used different numbers of spherical shells to calculate statistical distributions. are less accurate than 100 shell simulations, we did not expect to see an actual numerical convergence of the tabulated quantities. To see this, calculations using more spherical shells each containing more matter test particles than they currently do must be performed and this cannot happen until the code is converted to 64—bit and parallelized. The fact that the radius at which the outward explosion of matter formed in the 50 and 75 shell simulations is within i5 km of the 100 shell result is encouraging as is the fact that the density at these radii at the time of explosion were all quite similar. In addition to running convergence tests using more matter test particles and spherical shells, convergence tests can be performed in conjunction with the works of other groups. To do this, several things must be noted. First the main difference between the neutrino capture driven explosion mechanism introduced here and the traditional picture of bounce is the capture of neutrinos that exist between the limits of trapping and free streaming. Traditional hydrodynamic treatments of neutrino transport are not as reliable in this realm while our beam attenuation arguments are equally applicable everywhere in the core at all times. Thus we want to avoid comparisons that involve modeling neutrino-matter interactions. Therefore two ob- vious avenues present themselves. One is simply to compare our simulation to that of hydrodynamics groups with all neutrino-matter interaction algorithms deactivated in both codes. In this case good agreement can be expected. Another possible test would be to run our code through the time the outward explosion of matter began 236 and then use the data generated as initial conditions for a hydrodynamic code and see if its core explodes in a similar fashion. All of these test will be performed once the code is converted to 64-bit and parallelized. We note that we can make some comparisons between our neutrinos modeled in our simulations and those modeled in traditional hydrodynamic calculations. As long we restrict these comparisons to regions in which the neutrinos are either trapped of free streaming where traditional hydrodynamic treatments of neutrino transport are sufficient, convergence can be expected. This is seen to be the case when figure 8.29 is considered. The energies in the hottest densest region of the core have the energies of 100—150 MeV and the escaping neutrinos have energies in the 5-10 MeV range. These are the results predicted by hydrodynamic calculations [16]. This is an encouraging result. The only region in which our neutrinos behave unexpectedly is in the intermediate range in which only our formalism can realistically model their propagation. 237 Chapter 9 Conclusion Clearly the main achievement of this work is the possible discovery of an entirely new supernova explosion mechanism. The confluence of events that lead to the rapid and highly localized accumulation of electrons produced by neutrino captures at radii on the order of 50 km that alters the gradient of the in such a way that it cuts the core into two parts, a well defined proto-remnant and ejecta, is a completely novel concept that are code is uniquely poised to model. We are encouraged to see that this neutrino capture driven mechanism is robust enough to launch explosions in all simulations using three different nucleon potentials, four different electron cap- ture rate tables, and three different numbers of spherical shells. However we realize that in addition to converting the code to 64—bit and parallelizing the code, activating the three-dimensional subroutines, and conducting more numerical convergence tests, there is still more physics that needs to be added to our model. Fusion and photodis- integration must be built in, so do several weak reactions mentioned in section 3.11, the effects of neutrino and neutron degeneracy, additional nuclear decay modes, and a temperature dependent nucleon potential. It is clear how to accomplish all of this with the test particle approach. We do not yet regard this work as complete, but we are confident that we are 238 working in the right direction with our approach. The ability of the of the test particle method to treat the dynamics of baryons and neutrinos on equal footing, explicitly model the propagation of neutrinos in a general way that is applicable in all regions of the core at all times, and explicitly model the propagation of a full ensemble of nuclei are all significant steps forward. As explained in section 8.2, without all of these assets that only our code possesses, our code would not be capable of generating a neutrino capture driven explosion. They each play a critical role in its realization. We look forward to continuing to improve out model in the ways mentioned above and seeing if the neutrino capture driven explosion mechanism is still observed. Some physics related to weak reactions yet to be built in to our model is expected to have the effect of reducing the magnitude of the localized accumulation of electrons that powers the explosion. The extent to which they do this remains to be seen. Specif— ically the effects of neutrino degeneracy and neutrino pair production will tend to stabilize the core against developing hyper-leptonizationed regions through neutrino capture. Clearly the inclusion of neutrino degeneracy can have the effect of forbidding some electron captures from occurring in regions where neutrino number densities be- come sufficiently large. Furthermore neutrino pair production generates antineutrinos that can be captured by free protons and nuclei thereby creating positrons that can deleptonize the core without necessarily producing neutrinos. Since neutrino degeneracy and the highly temperature dependent neutrino pair production rates do not become relevant until high densities and temperatures are achieved in the innermost region of the core [68], the neutrino capture driven explo- sion is well underway before we need to concern ourselves with their effects. However continuous neutrino capture in the hyper-leptonized region is responsible for increas- ing the rate at which the ejecta is accelerated out of the core. Thus the impact that neutrino degeneracy and pair production have on the core dynamics my still prove to be problematic for the neutrino capture driven explosion mechanism. At this point we 239 do not believe that the effects of neutrino degeneracy and pair production will prove to be fatal for the neutrino capture driven explosion mechanism, however we cannot be certain of this until all of the physics discussed above is built into our model. Once full three-dimensional calculations have been run, we expect to see devia- tions from spherical symmetry in the density distribution that are both intuitively anticipated and have been seen in previous preliminary calculations [107]. We have long since hypothesized that the resultant density depletion along the axis of rotation will lead to focussing of neutrino emission along the poles which will amplify the par- ity violation induced recoil kick scenario proposed for the neutron star remnant by Horowitz et al. [155, 156]. It will be interesting to see how this phenomenon might effect the neutrino capture driven explosion mechanism assuming that it proves to still be viable. Furthermore, three-dimensional simulations will explicitly enforce the conservation of momentum in each neutrino test particle matter test particle inter- action that takes place inside the three-dimensional distribution grid described in section 3.9.2. This is expected to enhance the explosion mechanism as, in addition to the ejecta being subject to the outward force exerted by the electron gas, it will acquire net outward momentum through the enforcement of momentum conservation in neutrino-matter interactions. However, before such observations can be made, we must implement the aforementioned changes. This thesis has laid out the path for a multi—year investigation into the new ways to understand the physics of supernova explosions and has taken important first steps. The future of this work looks extremely bright, and it promises to provide a new intellectual bridge between nuclear physics and astrophysics. In particular, it has the capability to serve as the nexus for efforts made at nuclear physics labo- ratories such as N SCL and FRIB and laboratories that study neutrino physics like DUSEL, the Deep Underground Science and Engineering Laboratory. The explicit modeling the propagation of neutrinos and a full ensemble of nuclei, coupled with 240 the completely general treatment of neutrino-matter interactions the kinetic theory formalism naturally employs, gives our model the ability to use output generated by the aforementioned facilities to simulate the evolution of the nuclear composition of the core in ways that no other supernova simulation is currently capable of. As this model advances, we hope to see its kinetic theory based approach establish itself as the new gold standard for supernova calculations. 241 [ll [2] [3] [4] [5] [6] [7] l8] [9] BIBLIOGRAPHY H.A. Bethe. Supernova mechanisms. Rev. Mod. Phys, 62:801-866, October 1990. H.A. Bethe. Supernova theory. Nucl. Phys. A, 606:95-117, 1996. H.A. Bethe. Nuclear physics needed for the theory of supernovae. Ann. Rev. Nucl. Part. Sci, 38:1-28, 1988. D. Arnett, J .W. Truran, and SE. Woosley. Nucleosynthesis in supernova models II. the 12C detonation model. Astrophys. J., 165:87—103, April 1971. D. Arnett. Supernovae and Nucleosynthesis: An investigation of the history of matter, from the big bang to the present. Princeton Univ. Press, Princeton, 1996. D. Arnett. Origin and Evolution of the Elements. in Carnegie Observatories Cen- tennial Symposia. eds. A. McWilliam and M. Rauch, (Cambridge Univ. Press, Cambridge, 2004) p. 12. T. Nakamura, H. Umeda, K. Nomot, F. Thielemann, and A. Burrows. Nucleosyn— thesis in type II supernovae and the abundances in metal-poor stars. Astrophys. J., 517:193—208, May 1999. C. Fryer, P. Young, M. Bennett, S. Diehl, F. Herwig, R. Hirschi, A. Hunger- ford, M. Pignatari, G. Magkotsios, G. Rockefeller, and F.X. Timmes. Nucle- osynthesis calculations from core-collapse supernovae. to appear in 10th Sym- posium on Nuclei in the Cosmos (NIC X) Conf. Proc., http://pos.sissa.it/cgi- bin / reader / conf. cgi?confid=53. P. Young and C. Fryer. Uncertainties in supernova yields. I. one-dimensional ex- plosions. Astrophys. J., 664:1033-1044, August 2007. [10] C. Fryer, F. Herwig, A. Hungerford, and F.X. Timmes. Supernova fallback: a possible site for the r-process. Astrophys. J. Lett. 646 L131-L134, July 2006. [11] M. Limongi and A. Chieffi. Evolution, explosion, and nucleosynthesis of core- collapse supernovae. Astrophys. J., 592:404—433, July 2003. 242 [12] A. Chieffi and M. Limongi. Explosive yields of massive stars from Z = 0 to Z = Z9. Astrophys. J., 608:405-410, June 2004 [13] G. Martinez. Nuclear physics aspects of supernovae evolution and nucleosynthe- sis. J. Phys. C: Nucl. Part. Phys, 352014057, January 2009. [14] C. Fréihlich, C. Martinez, M. Liebendfirfer, F.-K. Thielemann, E. Bravo, W. R. Hix, K. Langanke, and N .T. Zinner. Neutrino-induced nucleosynthesis of A>64 nuclei: The Vp-process. Phys. Rev. Lett., 96:142502, 2006. [15] D. Prialnik. An Introduction to the Theory of Stellar Structure and Evolution. Cambridge Univ. Press, Cambridge, 2000. [16] K. Langanke. Lectures in Nuclear Astrophysics. TUD Darmstadt, summer semester 2006, http: / / theory. gsi.de / ~langanke / . [17] MT. Bollenbach. Numerical Study of Rotating Core Collapse Supernovae. MS thesis, Michigan State Univ., 2002. [18] SE. Woosley and TA. Weaver. The physics of supernova explosions. Ann. Rev. Astron. Astrophys., 24:205-253, 1986. [19] H.-Th. Janka. Conditions for shock revival by neutrino heating in core-collapse supernovae. Astron. Astrophys., 268:360—368, 1993. [20] J. Murphy and A. Burrows. Criteria for core-collapse supernova explosions by the neutrino mechanism. submitted to Astrophys. J., July 2008. [21] J. Wilson. A Numerical Study of Gravitational Stellar Collapse. Astrophys. J., 163:209—219 January 1971. [22] W. Hillebrandt. An exploding 10 solar mass star - A model for the Crab super- nova. Astron. Astrophys., 110:1:L3—L6, June 1982. [23] W.D. Arnett. Neutrino escape, nuclear dissociation, and core collapse and/ or explosion. Astrophys. J., 263:L55-L57, December 1982. [24] W. Hillebrandt, K. Nomoto, K., and RC. Wolff. Supernova explosions of massive stars - the mass range 8 to 10 solar masses. Astron. Astrophys., 133:1:175—184, April 1984. - [25] E. Baron, J. Cooperstein, and S. Kahana. Type-II supernovae in 12-solar-mass and 15-solar—mass stars the equation of state and general relativity. Phys. Rev. Lett., 55:126-129, July 1985. [26] J .R. Wilson, R. Mayle, S.E. Woosley, and T. Weaver. Stellar core collapse and supernova. in Twelfth Texas Symposium 0n Relativistic Astrophysics. eds. M Livio and G. Shavi, 470:267—293, December 1986. 243 [27] H. Bethe and JR. Wilson. Revival of a stalled supernova shock by neutrino heating. Astrophys. J., 295214-26, August 1985. [28] M. Liebendérfer, A. Mezzacappa, F.-K. Thielemann, OE. Messer, R. Hix, and S. Bruenn. Probing the gravitational well: No supernova explosion in spherical symmetry with general relativistic boltzmann neutrino transport. Phys. Rev. D., 63, 10, 103004, May 2001. [29] M. Rampp and H.-Th. Janka. Radiation hydrodynamics with neutrinos. Astron. and Astophys., 396:361-392, September 2002. [30] R.Buras, M. Rampp, H.-Th. Janka, and K. Kifonidis. Improved models of stellar core collapse and still no explosions: What is missing?. Phys. Rev. Lett., 90, 241101, 2003. [31] M. Liebendrfer, M. Rampp, H.-Th. Janka, and A. Mezzacappa. Supernova sim- ulations with boltzmann neutrino transport: a comparison of methods. Astrophys J., 620:840—860, January 2005. [32] A. Marek and H.-Th. Janka. Delayed neutrino-driven supernova explosions aided by the standing accretion-shock instability. Astrophys. J., 694:664—696, 2009. [33] C. FTyer and P. Young, Late-time convection in the collapse of a 23 MG) star. Astrophys. J., 659:1438—1448, April 2007. [34] C. Fryer and A. Heger. Core-collapse simulations of rotating stars. Astrophys. J., 541:1033-1050, October 2000. [35] M. Herant, W. Benz, and S. Colgate. Postcollapse hydrodynamics of SN 1987A: Two-dimensional simulations of the early evolution. Astrophys. J., 395:642-653, August 1992. [36] M. Herant, W. Benz, W. Hix, C. Fryer and S. Colgate. Inside the supernova: A powerful convective engine. Astrophys. J., 435:339—361, November 1994. [37] A. Burrows, E. Livne, L. Dessart, C.D. Ott, and J. Murphy. A new mechanism for core-collapse supernovae explosions. Astrophys. J., 640:878-890, 2006. [38] A. Burrows, E. Livne, L. Dessart, C.D. Ott, and J. Murphy. Features of the acoustic mechanism of core-collapse supernova explosions. Astrophys. J., 655:416- 433, 2007. [39] J.M. Lattimer and FD. Swesty. A generalized equation of state for hot, dense matter. Nucl. Phys. A, 535:331-367, 1991. [40] H. Shen, H. Toki, K. Oyamatsu, and K. Sumiyoshi. Relaistic equation of state of nuclear matter for supernova and neutron star. Nucl. Phys. A, 637(3):435—450. May 1998. 244 [41] K. Sumiyoshi, H. Shen, K. Oyamatsu, M. Terasawa, H. Suzuki, S. yamada, and H. Toki. Unstable nuclei and EOS table for supernova explosion and r—process in relativistic many body approach. RIKEN Review, 26, January 2000. [42] K. Sumiyoshi, M. Terasawa, H. Suzuki, S. Yamada, H. Toki, G.J. Mathews, and T. Kajino. Relativistic simulations of supernovae and the r-process: a new relativistic EOS and nuclear reaction network. Nucl. Phys. A, 688:4780—4800, 2001. [43] K. Sumiyoshi, and H. Toki. Relativistic equation of state of nuclear matter for the supernova explosion and the birth of neutron stars. Astrophys. J., 422:700—718, February 1994. [44] F.D. Swesty, J .M. Lattimer, and ES. Myra. The role of the equation of state in the “prompt phase” of type II supernovae. Astrophys. J., 425:195—204, April 1994. [45] A. Burrows and J. Goshy. A theory of supernova explosions. Astrophys. J., 416:L75-L78, October 1993. [46] T. Zwerger and E. Miiller. Dynamics and gravitational wave signature of axisym- metric core collapse. Astron. Astrophys., 320:209-227, 1997. [47] S. Yamada and K. Sato. Gravitational radiation from rotational collapse of a supernova core. Astrophys. J., 450:245—252, September 1995. [48] S. Bonazzola and J .A. Marek. Efficiency of gravitational radiation from axisym- metric and 3D stellar collpase. Astron. Astrophys., 267:623—633, 1993. [49] S. Bruenn, K.R. De Nisco, and A. Mezzacappa. General relativistic effects in the core collapse supernova mechanism. Astrophys. J., 560:326—338, October 2001. [50] O.E.B. Messer, A. Mezzacappa, S. Bruenn, and M. Guidry. A comparison of boltzmann and multigroup flux-limited diffusion neutrino transport during the postbounce shock reheating phase in the core-collapse supernovae. Astrophys. J., 507:353—360, November 1998. [51] A. Mezzacappa, A.C. Calder, S. Bruenn, J.M Blondin, M.W. Guidry, M.R. Strayer, and AS. Umar. The interplay between proto—neutron star convection and neutrino transport in core-collapse supernovae. Astrophys. J., 493:848-862, February 1998. [52] A. Mezzacappa, M. Liebendorfer, O.E.B. Messer, W. Hix, F.-K. Thielemann, and S. Bruenn. Simulation of the spherically symmetric stellar core collapse, bounce, and postbounce evolution of a star of 13 solar masses with boltzmann neutrino transport, and its implications for the supernova mechanism. Phys. Rev. Lett., 86(10):1935, March 2001. [53] M. Rampp and H.-Th. Janka. Spherically symmetric simulation with boltzmann neutrino transport of core collapse and postbounce evolution of a 15 MO star. Astrophys. J., 539:L33-L36, August 2000. 245 [54] F.X. Timmes and D. Arnett. The accuracy, consistancy, and speed of five equa- tions of state for stellar hydrodynamics. Astrophys. J. Suppl. S., 125:277—294, November 1999. [55] K. Takahashi, M. F. El Bid, and W. Hillebrandt. Beta transition rates in hot dense matter. Astron. Astrophys., 267:623-633, 1978. [56] G. Fuller, W. Fowler, and M. Newman. Stellar weak interaction rates for intermediate-mass nuclei IV. Interpolation procedures for rapidly varying lepton capture rates using effective log( ft) values. Astrophys. J., 293:1-16, June 1985. [57] G. Kortemeyer, F. Daflin, and W. Bauer, Nuclear flow in consistent boltzmann algorithm models. Phys. Lett., B 374225-30, 1996. [58] G. Martinez, K. Langanke, and DJ. Dean. Competition of electron capture and beta-decay rates in supernova collapse. Astrophys. J. Suppl. S., 126:493-499, February 2000. [59] G. Martinez, M. Liebendoerfer, and D. Frekers. Nuclear input for core-collapse models. Nucl.Phys., A 777:395-423. 2006. [60] G. Martinez and K. Langanke. Nuclear weak-interaction processes in stars. Rev. Mod. Phys, 75:819-862, 2003. [61] A. Juodagalvis, K. Langanke, W. Hix, G. Martinez, and J .M. Sampaio. Improved estimate of stellar electron capture rates on nuclei. submitted to Phys. Rev. C, September 2009. [62] G. Martinez. Shell-model applications in supernova. physics. Eur. Phys. J., 25:s01:659—664, June 2005. [63] K. Langangke and G. Martinez. Nuclear weak-interaction processes in stars. Rev. Mod. Phys, 75:819-862, June 2003. [64] K. Langanke. Nuclear physics and core collapse supernovae. Nucl. Phys. A, 690229—40, July 2001. [65] J. Pruet and G. Fuller. Estimates of stellar weak interaction rates for nuclei in the mass range A = 65—80. Astrophys. J. Suppl. 3., 149:189—203, November 2003. [66] G.W. Hitt, assistant professor of physics at Khalifa University, private commu- nication, August 2009. [67] O. E. B. Messer, A. Mezzacappa, S. Bruenn, and M.W. Guidry. A compari- son of boltzmann and multigroup flux-limited diffusion neutrino transport during the postbounce shock reheating phase in core-collapse supernovae. Astrophys. J., 507:353—360, November 1998. 246 [68] A. Burrows and T. Thompson. Neutrino-matter interaction rates in supernovae: The essential microphysics of core collapse. to be published in Core Collapse of Massive Stars, ed. C. Fryer, Kluwer Academic Publishers. [69] L. Scheck, T. Plewa, H.-Th. Janka, K. Kifonidis, and E. Mueller. Pulsar recoil by large-scale anisotropies in supernova explosions. Phys.Rev.Lett., 92:11-103, 2004. [70] K. Kotake and S. Yamada. Neutrino— and gravitational-wave astronomy of massive-star collapse. AAPPS Bulletin, 19:3:30—35, June 2009. [71] A. Mezzacappa. Toward a standard model of core collapse supernovae. Nucl. Phys. A, 688:158-167, June 2000. [72] K. Kotake, N. Ohnishi, and S. Yamada. Gravitational radiation from standing accretion shock instability in core-collapse supernovae. Astrophys. J., 655:406—415, January 2007. [73] C. Fryer. Mass limits for black hole formation. Astrophys. J., 522:413—418, September 1999. [74] RA. Gingold and J .J . Monaghan. Smooth particle hydrodynaimcs: theory and application to non-spherical stars. Mon. Not. Roy. Astro. Soc., 181:375—389, 1997. [75] M. Herant and W. Benz. Postexplosion hydrodynamics of SN 1987A. Astrophys. J., 387:294-308, March 1992. [76] A. Heger and N. Langer. Presupernova evolution of rotating massive stars. II. Evolution of the surface properties. Astrophys. J., 544:1016—1035, December 2000. [77] L. Wang, D.A. Howell, P. Héiflich, and J.C. Wheeler. Bipolar supernova explo- sions. Astrophys. J., 550:1030-1035, April 2001. [78] L. Wang and J.C. Wheeler. Supernovae are not round. Sky and Telescope, Jan- uary 2002. [79] L. Wang, J.C. Wheeler, Z. Li, and A. Clocchiatti. Broadband polarimetry of supernovae: SN 1991D, SN 1994Y, SN 1994AE, SN 1995D, and SN 1995H. As- trophys. J., 467:435—445, August 1996. [80] E.Livne. An implicit method for two-dimensional hydrodynamics. Astrophys. J., 412:634—647, August 1993. [81] T. Thompson, A. Burrows, and P. Pinto. Shock breakout in core-collapse super- novae and its neutrino signature. Astrophys. J., 592:434-456, 2003. [82] C. Fryer, G. Rockefeller, and M. Warren. SNSPH: A parallel three-dimensional smoothed particle radiation hydrodynamics code. Astrophys. J., 643:292-305, May 2006. 247 [83] C. Fryer and M. Warren. Modeling core-collapse supernovae in three dimensions. Astrophys. J., 574:L65-L68, July 2002. [84] G. Rockefeller, C. Fryer, P. Young, M. Bennett, S. Diehl, F. Herwig, R. Hirschi, A. Hungerford, M. Pignatari, G. Magkotsios, and F .X. Timmes. Nucleosynthetic yields from collapsars. to appear in 10th Symposium on Nuclei in the Cosmos (NIC X) Conf. Proc., http://arxiv.org/pdf/0811.4650v1. [85] S. Diehl, C. Fryer, A. Hungerford, G. Rockefeller, M. Bennett, F.Herwig, R. Hirschi, M. Pignatari, G. Magkotsios, F .X. Timmes, P. Young, G. Clayton, P. Motl, and J. Tohline. NuGrid: Toward high precision double-degenerate merger simulations with SPH in 3D. to appear in 10th Symposium on Nuclei in the Cos- mos (NIC X) Conf. Proc., http://arxiv.org/pdf/0811.4646v1. [86] P. Young, C. Meakin, D. Arnett, and C. Fryer. The impact of hydrodynamic mixing on supernova progenitors. Astrophys. J., 629:L101-L104, August 2005. [87] L. Scheck, T. Plewa, H.-Th. Janka, K. Kifonidis, and E. Miiller. Pulsar recoil by large-scale anisotropies in supernova explosions. Phys. Rev. Lett., 92:011103, January 2004. [88] C. Fryer. Neutron star kicks from asymmetric collapse. Astrophys. J., 601:L175- L178, February 2004. [89] W. Hix, O. E. B. Messer, A. Mezzacappa, M. Liebendorfer, J. Sampaio, K. Langanke, D. Dean, and G. Martinez. Consequences of nuclear electron capture in core collapse supernovae. Phys. Rev. Lett., 91: 201102, 2003. [90] G.F. Bertsch, H. Kruse, and S. Das Gupta. Boltzmann equation for heavy ion collisions. Phys. Rev. C, 29:673-675, 1984. [91] H. Kruse, B.V. Jacak, and H. Stocker. Microscopic theory of pion production and sidewards flow in heavy-ion collisions. Phys. Rev. Lett., 4:289—292, 1985. [92] W. Bauer, G.F. Bertsch, W. Cassing, and U. Mosel. Energetic photons from inter- mediate energy proton- and heavy-ion-induced reactions. Phys. Rev. C, 34:2127- 2133, December 1986. [93] H. Stocker, and W. Greiner. High energy heavy ion collisionsprobing the equation of state of highly excited hardronic matter. Phys. Rep, 137:277-392, May 1986. [94] G.F. Bertsch and S. Das Gupta. A guide to microscopic models for intermediate- energy heavy ion collisions. Phys, Rep., 160:189—233, 1988. [95] P. Schuck, R.W. Hasse, J. Jaenicke, C. Grégoire, B. Rémaud, F. Sébille, and E. Suraud. Semiclassical and phase space approaches to dynamic and collisional problems of nuclei. Prog. Part. Nucl. Phys. 22:181-278, 1989. 248 [96] W.G. Gong, W. Bauer, C.K. Gelbke, and S. Pratt. Space-time evolution of nuclear reactions probed by two-proton intensity interferometry. Phys. Rev. C, 432781-800, February 1991. [97] H.-Th. Janka, R. Buras, and M. Rampp. The mechanism of core-collapse su- pernovae and the ejection of heavy elements. Nucl. Phys. A, 718:269—276, May 2003. [98] JR. Wilson, Numerical Astrophysics. Jones and Bartlett, Boston, 1985. [99] M.A. Cirit. Hydrodynamic model of the proton-nucleus collision. Phys. Rev. D, 21(7):1854—1860, April 1980. [100] J.J. Bai, R.Y. Cusson, J. Wu, P.-G. Reinhard, H. Stoecker, W. Greiner, and MR. Strayer. Mean field model for relativistic heavy ion collisions. Z. Phys. A- Atomic Nuclei, 326:269-277, 1987. [101] Y. Kitazoe. Cascade-model analysis of pion production in relativistic nuclear collisions. Prog. Theor. Phys, 73(5):1191-1197, May 1985. [102] W.G. Gong, W. Bauer, C.K. Gelbke, N. Carlin, R.T. de Souza, Y.D. Kim, W.G. Lynch, T. Murakami, G. Poggi, D.P. Sanderson, M.B. Tsang, H.M. Xu, S. Pratt, D.E. Fields, K. Kwiatkowski, R. Planeta, V.E. Viola, Jr., and SJ. Yennello. Intensity-interferometric test of nuclear collision geometries obtained from the boltzmann-uehling—uhlenbeck equation. Phys. Rev. Lett., 65:2114—2117, October 1990. [103] W.G. Gong, W. Bauer, C.K. Gelbke, N. Carlin, R.T. de Souza, Y.D. Kim, W.G. Lynch, T. Murakami, G. Poggi, D.P. Sanderson, M.B. Tsang, H.M. Xu, S. Pratt, D.E. Fields, K. Kwiatkowski, R. Planeta, V.E. Viola, Jr., and SJ. Yennello. Space-time evolution of the reactions 14N+27A1, 197Au at E /A = 75 MeV and 129Xe+27Al, 122Sn at E /A = 31 MeV probed by two-proton intensity interferometry. Phys. Rev. C, 43(4):1804—1820, April 1991. [104] BA. Li, W. Bauer, and G.F. Bertsch. Preferential emission of pions in asym- metric nucleus-nucleus collisions. Phys. Rev. C, 44(5):2095-2099, November 1991. [105] C.—Y. Wong. Dynamics of nuclear fluid. VIII. Time—dependent Hartree-Fock ap- proximation from a classical point of View. Phys. Rev. C, 25(3):1460—1475, March 1982. [106] W. Bauer. Dynamical simulations of supernovae collapse and nuclear collisions via the test particle method - similarities and differences. Acta Phys. Hung. A, 21:371-376, 2004. [107] W. Bauer and T. Strother. Collective Motion in Nuclear Collisions and Super- nova Explosions. Int. J. Phys. E, 14:129-137, 2005. 249 [108] S.J. Wang, B.A. Li, W. Bauer, and J. Randrup. Relativistic transport theory for hadronic matter. Ann. Phys, 209:251-305, August 1991. [109] W. Bauer and G. Westfall. University Physics. McGraw Hill, (2010) [110] W.H. Press, S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery. Numerical Recipes in C. Press syndicate of the University of Cambridge, Cambridge, UK, second edition, 1988. [111] T. Strother and W. Bauer. Modeling weak interaction rates during the super- nova collapse phase. Prog. Part. Nucl. Phys, 622468-472, April 2009. [112] W. Bauer and T. Strother. Nuclear Reactions and Stellar Evolution: Unified Dynamics. in AIP Conf. Proc., eds. R. Alarcon, P. L. Cole, C. Djalali, and F. Umeres, 947:333—340, Springer, Cusco, 2007. [113] D. Clayton. Principles of Stellar Evolution and Nucleosynthesis. McGraw Hill, New York, 1968. [114] F. Shu. The Physical Universe An Introduction to Astronomy. University Sci- ence Books, Sausalito, 1982. [115] C. Fryer. private communication, April 2009. [116] C. Fryer. private communication, May 2004. [117] T. Strother and W. Bauer. The nuclear and many-body physics of supernova explosions. Int. J. Mod. Phys. E, 16:1073—1081, 2007. [118] J. Marion and S. Thorton. Classical Dynamics of Particles and Systems. Brooks/ Cole, Belmont, 2004. [119] L. Shi and P. Danielewicz. Effect of cluster formation on isospin asymmetry in the liquid-gas phase transition region. Europhys. Lett., 51:34-40, July 2000. [120] L.-W. Chen, B.-J. Cai, C.M. Ko, B.-A. Li, C. Shen, and J. Xu. Higher-order effects on the incompressibility of isospin asymmetric nuclear matter. Phys. Rev. C, 80:014322, July 2009. [121] B.-A. Li, L.—W. Chen, C.M. Ko. Recent progress and new challenges in isospin physics with heavy-ion reactions. Phys. Rep., 464:113-281, August 2008. [122] B.—A. Li, C. Das, S. Das Gupta, and C. Galeb. Effects of momentum-dependent symmetry potential on heavy-ion collisions induced by neutron-rich nuclei. Nucl. Phys. A, 735:563—584, May 2004. [123] J. Xu, L.-W. Chen, B.-A. Li, and H.-R. Ma, Nuclear constraints on properties of neutron star crusts. Astrophys. J., 697:1549—1568, June 2009. 250 [124] C. Das,1 S. Das Gupta, C. Gale, and B.-A. Li. Momentum dependence of sym- metry potential in asymmetric nuclear matter for transport model calculations. Phys. Rev. C, 67: 034611, March 2003. [125] P. Danielewicz and J. Lee. Symmetry energy I: Semi-infinite matter. Nucl. Phys. A, 818236-96, 2009 [126] J. M. Lattimer and M. Prakash. Neutron star structure and the equation of state. Astrophys. J., 550:426—442, March 2001. [127] B.-A. Li, C.M. Ko, W. Bauer. Isospin physics in heavy ion collisions at intermediate-energies. Int. J. Mod. Phys. E, 72147-230, 1998. [128] Need to cite L. Landau and E. Lifshitz. Fluid Mechanics. Pergamon Press, Oxford, 1987. [129] S. Bruenn. Stellar core collapse: Numerical model and infall epoch. Astrophys. J. Supp, 582771-841, August 1985 [130] M. Rampp and H.-T. Janka. Radiation hydrodynamics with neutrinos. Vari- able Eddington factor method for core-collapse supernova simulations. Astron. Astrophys., 396:361-392, December 2002. [131] K. Langanke and C. Martinez. Nuclear weak-interaction processes in stars. Rev. Mod. Phys, 752819-863, June 2003. [132] A. Burrows, S. Reddy, and T. Thompson. Neutrino opacities in nuclear matter. Nucl. Phys. A, 777:356-394, October 2006. [133] S. Hannestad and G. Raffelt. Supernova neutrino opacity from nucleon-nucleon bremsstrahlung and related processes. Astrophys. J., 507:339—352, November 1998. [134] T. Thompson, A. Burrows, J. Horvath. u and r neutrino thermalization and production in supernovae: Processes and time scales. Phys. Rev. C, 62:035802, September 2000. [135] J. R. Bond, Ph.D. thesis, Caltech, 1979. [136] R. Buras, H.-T. Janka, M. T. Keil, G. Raffelt, and M. Rampp. Electron neu- trino pair annihilation: A new source for muon and tau neutrinos in supernovae. Astrophys. J., 578:329—326, April 2003. [137] M.T. Keil, G. Raffelt, and H.-T. Janka. Monte carlo study of supernova neutrino spectra formation. Astrophys. J., 590:971-991, June 2003. [138] A. Heger and S. E. Woosley. Presupernova evolution with improved rates for weak interactions. Astrophys. J., 560:307—325, October 2001. [139] L. Wolfenstein. Neutrino oscillations in matter. Phys. Rev. D, 17:2369—2374, May 1978. 251 [140] B. Kayser. 13. Neutrino mass mixing and flavor change. Particle Data Group Rev. , http:/ / pdg.lb1. gov / 2009 / reviews / rpp2009—rev—neutrino—mixing.pdf, March 2008. [141] L. Lu, W.-Y. Wang, and Z.-H. Xiong. Parameterization for neutrino mixing matrix with deviated unitarity. Chin. Phys. Lett., 262081401, August 2008. [142] N. Lia and B.-Q. Ma. A new parametrization of the neutrino mixing matrix. Phys. Lett. B, 600:248—254, October 2004. [143] W. Bauer. private communication, August, 2009. [144] A. Watson. The Quantum Quark. Cambridge Univ. Press, Cambridge, 2004. [145] K. Krane. Introductory Nuclear Physics. New York, Wiley, 1987. [146] K. Huang. Statistical Mechanics. New York, John Wiley & Sons, 1987. [147] D. Schroeder. An Introduction to Thermal Physics. Addison Wesley Longman, San Francisco, 2000. [148] P. Vogel. Analysis of the antineutrino capture on protons. Phys. Rev. D, 29:1918— 1922,1984. [149] H. Bethe, G. Brown, J. Applegate, and J.M. Lattimer. Equation of state in the gravitational collapse of stars. Nucl. Phys. A, 324:487—533, 1979. [150] G. Fuller. Neutron shell blocking of electron capture during gravitational col- lapse. Astrophys. J., 252:741-764, January 1982. [151] S. Bruenn and A. Mezzacappa. Ion screening effects and stellar collapse. Phys. Rev. D, 56:7529-7547, 1997. [152] R. Bowers and JR. Wilson. A numerical model for stellar core collapse calcu- lations. Astrophys. J. Supp, 502115-160, October 1982. [153] D.Tubbs and D. Schramm. Neutrino Opacities at high temperatures and den- sities. Astrophys. J., 201:467—488, October 1975. [154] T. Strother and W. Bauer. Nuclear physics and supernova explosions: Unified dynamics. to appear in Int. J. Mod. Phys. D, 2009. [155] C.J. Horowitz and G. Li. Cumulative parity violation in supernovae. Phys. Rev. Lett., 80:3694—3699, 1998. [156] C.J. Horowitz and J. Piekarewicz. Macroscopic parity violation and supernova asymmetries. Nucl. Phys. A, 640:281-300, 1998. 252