A 1-D SPECIAL RELATIVISTIC SHOCK TUBE WHERE STELLAR FLUID UNDERGOES NEUTRINO HEATING by Gregory Mohammed B.Sc., University of Northern British Columbia, 2005 THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE IN MATHEMATICAL, COMPUTER, AND PHYSICAL SCIENCES (PHYSICS) UNIVERSITY OF NORTHERN BRITISH COLUMBIA December 2012 (c) Gregory Mohammed, 2012 1+1 Library and Archives Canada Bibliotheque et Archives Canada Published Heritage Branch Direction du Patrimoine de I'edition 395 Wellington Street Ottawa ON K1A0N4 Canada 395, rue Wellington Ottawa ON K1A 0N4 Canada Your file Votre reference ISBN: 978-0-494-94127-0 Our file Notre reference ISBN: 978-0-494-94127-0 NOTICE: AVIS: The author has granted a non­ exclusive license allowing Library and Archives Canada to reproduce, publish, archive, preserve, conserve, communicate to the public by telecommunication or on the Internet, loan, distrbute and sell theses worldwide, for commercial or non­ commercial purposes, in microform, paper, electronic and/or any other formats. L'auteur a accorde une licence non exclusive permettant a la Bibliotheque et Archives Canada de reproduire, publier, archiver, sauvegarder, conserver, transmettre au public par telecommunication ou par I'lnternet, preter, distribuer et vendre des theses partout dans le monde, a des fins commerciales ou autres, sur support microforme, papier, electronique et/ou autres formats. The author retains copyright ownership and moral rights in this thesis. Neither the thesis nor substantial extracts from it may be printed or otherwise reproduced without the author's permission. L'auteur conserve la propriete du droit d'auteur et des droits moraux qui protege cette these. Ni la these ni des extraits substantiels de celle-ci ne doivent etre imprimes ou autrement reproduits sans son autorisation. In compliance with the Canadian Privacy Act some supporting forms may have been removed from this thesis. Conform em ent a la loi canadienne sur la protection de la vie privee, quelques formulaires secondaires ont ete enleves de cette these. W hile these forms may be included in the document page count, their removal does not represent any loss of content from the thesis. Bien que ces formulaires aient inclus dans la pagination, il n'y aura aucun contenu manquant. Canada Table of Contents Table of Contents ii List of Figures v List of Tables vii Acknowledgements ix 1 Introduction 1 1.1 Stars and Their W ay s............................................................................................................ 2 1.2 Notation and C onventions.................................................................................................... 3 1.3 Introduction to the Godunov M e th o d ................................................................................ 4 2 3 Literature Review 7 2.1 The Supernova P a ra d ig m .................................................................................................... 7 2.2 Background on the Godunov M eth o d ................................................................................ 9 2.3 Two Exact Riemann S o l v e r s ............................................................................................. 10 2.4 High Resolution Shock C a p tu r in g ................................................................................... 11 2.5 Concerning N eu trin o s.......................................................................................................... 12 The 1-D Special Relativistic Hydrodynamics Equations for the Stellar Fluid including the Neutrino Model 14 3.1 1-D Special Relativistic H ydrodynam ics........................................................................... 15 3.1.1 18 Fluid Dynamics Equation of S t a t e ...................................................................... TABLE OF CO NTENTS 4 3.2 1-D Special Relativistic Neutrino H ydrodynam ics........................................................ 18 3.3 Prelude to the Neutrino M o d e l ......................................................................................... 21 3.4 The Neutrino Model .......................................................................................................... 24 3.5 Summary of E q u a tio n s ....................................................................................................... 30 The Godunov Method 31 4 The Finite Volume M e th o d ................................................................................................. 32 4.1.1 Methods of Calculating the F l u x .......................................................................... 34 The Godunov M eth o d .......................................................................................................... 36 4.2.1 Godunov A lg o r ith m .............................................................................................. 37 High Resolution M eth o d s................................................................................................... 38 4.2 4.3 4.4 5 4.3.1 The CFL C o n d itio n......................................................................................................40 4.3.2 Total Variation D im in is h in g ..................................................................................... 41 4.3.3 Slope-Limiter M eth o d s............................................................................................... 42 S u m m a r y .................................................................................................................................. 43 Simulation Results and Analysis 44 5.1 Sod Shock Tube Tests and R e s u lt s ....................................................................................... 46 5.2 Simulation Results using KKT data and Neutrino Fluxes Z e r o e d .............................. 67 5.3 Simulation Results with Neutrino Fluxes A p p lie d .......................................................... 85 5.3.1 Only the Neutrino Energy Flux A p p lie d ............................................................. 85 5.3.2 Only the Neutrino Momentum Flux Applied 5.3.3 Both Neutrino Fluxes A p p lie d .................................................................................103 5.4 6 iii S u m m a r y .................................................................................................................................112 Conclusions and Future Directions 6.1 113 Future W o rk ............................................................................................................................. 115 A Derivations and Notes A .l ........................................................94 118 Static Einstein Tensor C om ponents......................................................................................119 A.2 Static Einstein Field Equation S o lu tio n s ............................................................................152 TABLE OF CO NTENTS iv A.3 Static GR Equations of M o tio n ...........................................................................................153 A.4 Non-Static Einstein Tensor C o m p o n e n ts.......................................................................... 154 A.5 Non-Static Einstein Field Equation Solutions A .6 Geometrizations of Constants and Crucial Terms A.7 Development E n v iro n m e n t..................................................................................................169 A .8 Comments and Notes on the Programming P r o c e s s .......................................................170 A.9 The Parameter F i l e ...............................................................................................................171 ................................................................ 165 ..........................................................167 Bibliography 180 Index 183 List of Figures 1.1 The Grid Representation of the Finite Volume Method.................................................... 5 3.1 A volume through which the fluid flows.............................................................................. 15 3.2 Neutrino Model - Location of Shock T u b e ....................................................................... 25 3.3 Neutrino Model - Velocity Profile ..................................................................................... 27 3.4 Neutrino Model - Density at the N eu trin o -S p h ere.......................................................... 28 3.5 Neutrino Model - Neutrino E n e r g y .................................................................................... 29 4.1 Finite Volume Grid showing Fluxes..................................................................................... 32 4.2 Upwind Method for Advection............................................................................................. 35 4.3 The Godunov method and Riemann problems.................................................................... 37 5.1 Sod result for density................................................................................................................... 48 5.2 Sod result for pressure.................................................................................................................49 5.3 Sod result for internal energy................................................................................................. 50 5.4 Sod result for velocity............................................................................................................. 51 5.5 Sod result for density............................................................................................................... 52 5.6 Sod result for pressure............................................................................................................. 53 5.7 Sod result for internal energy................................................................................................. 54 5.8 Sod result for velocity............................................................................................................. 55 5.9 The standard deviations obtained for resolutions of 100, 200, 400, 800, 1600 and 10000 cells............................................................................................................................... 57 5.10 The order of error in the method........................................................................................... 57 v L IS T OF FIGURES vi 5.11 Reverse Sod rest density...................................................................................................... 59 5.12 Reverse Sod pressure............................................................................................................ 60 5.13 Reverse Sod result for internal energy................................................................................ 61 5.14 Reverse Sod result for velocity............................................................................................. 62 5.15 Percent Difference between Sod data in both directions.................................................. 62 5.16 The breakdown of rest density for v = 0.65............................................................................63 5.17 The breakdown of pressure for v = 0.65............................................................................. 64 5.18 The breakdown of internal energy for v = 0.65..................................................................... 65 5.19 The breakdown of velocity for v = 0.65.................................................................................. 6 6 5.20 Initial rest density data........................................................................................................... 68 5.21 Initial pressure data................................................................................................................ 69 5.22 Initial internal energy data.................................................................................................... 70 5.23 Initial velocity data................................................................................................................. 71 5.24 Kuroda rest density data with neutrinos off at 400 cells................................................. 72 5.25 Kuroda pressure data with neutrinos off at 400 cells...................................................... 73 5.26 Kuroda internal energy data with neutrinos off and at 400 cells................................... 74 5.27 Kuroda velocity data with neutrinos off and at 400 cells....................................................75 5.28 Kuroda rest density data with neutrinos off at 1600 cells...................................................76 5.29 Kuroda pressure data with neutrinos off at 1600 cells........................................................ 77 5.30 Kuroda internal energy data with neutrinos off and at 1600 cells..................................... 78 5.31 Kuroda velocity data with neutrinos off and at 1600 cells............................................. 79 5.32 Overlay of the Kuroda results for density at 400 and 1600 cells................................... 81 5.33 Overlay o f the Kuroda results for pressure at 400 and 1600 cells................................. 82 5.34 Overlay of the Kuroda results for internal energy at 400 and 1600 cells......................... 83 5.35 Overlay of the Kuroda results for velocity at 400 and 1600 cells................................. 84 5.36 Kuroda rest density data with neutrinos on at 400 and 1600 cells................................ 86 5.37 Kuroda pressure data with neutrinos on at 400 and 1600 cells..........................................87 5.38 Kuroda internal energy data with neutrinos on at 400 and 1600 cells......................... 88 5.39 Kuroda velocity data with neutrinos on at 400 and 1600 cells...................................... 89 L IS T OF FIGURES vii 5.40 Rest density data with neutrino energy flux on at times of 2.5s and 3.33s.................. 90 5.41 91 Pressure data with neutrino energy flux on at times of 2.5s and 3.33s......................... 5.42 Internal energy data with neutrino energy flux on at times of 2.5s and 3.33s..................92 5.43 Velocity data with neutrino energy flux on at times of 2.5s and 3.33s......................... 93 5.44 Kuroda rest density data with neutrinos on at 400 and 95 1600 cells............................. 5.45 Kuroda pressure data with neutrinos on at 400 and 1600 cells.......................................... 96 5.46 Kuroda internal energy data with neutrinos on at 400 and 1600 cells.............................. 97 5.47 Kuroda velocity data with neutrinos on at 400 and 1600 cells....................................... 98 5.48 Rest density data with neutrino momentum flux on at times of 2.5s and 3.33s. . . . 99 5.49 Pressure data with neutrino momentum flux on at times of 2.5s and 3.33s................... 100 5.50 Internal energy data with neutrino momentum flux on at times of 2.5s and 3.33s. . . 101 5.51 Velocity data with neutrino momentum flux on at times of 2.5s and 3.33s....................102 5.52 Kuroda rest density data with neutrinos on at 400 and 1600 cells................................104 5.53 Kuroda pressure data with neutrinos on at 400 and 1600 cells.........................................105 5.54 Kuroda internal energy data with neutrinos on at 400 and 1600 cells.............................106 5.55 Kuroda velocity data with neutrinos on at 400 and 1600 cells......................................... 107 5.56 Rest density data with both neutrino fluxes on at times of 2.5s and 3.33s......................108 5.57 Pressure data with both neutrino fluxes on at times of 2.5s and 3.33s............................ 109 5.58 Internal energy data with both neutrino fluxes on at times of 2.5s and 3.33s................ 110 5.59 Velocity data with both neutrino fluxes on at times of 2.5s and 3.33s.............................I l l List of Tables A .l The Christoffel symbols for a Static, Spherically Symmetric Spacetime........................130 A.2 The Ricci tensor components for a Static, Spherically Symmetric Spacetime A.3 Einstein components - Time Independent Spacetime........................................................ 151 A.4 Christoffel components - Time Dependent Spacetime....................................................159 A.5 Ricci components - Time Dependent Spacetime................................................................ 162 A .6 Einstein components - Time Dependent Spacetime...........................................................165 A.7 A listing of all the Geometrized Quantities......................................................................... 168 Vlll 146 Acknowledgements The author at this juncture would like to extend acknowledgments to the role of various people in the realization of this thesis, and also to extend the greatest of thanks and compliments to others who played a pivotal role in the completion of the thesis. First to the Physics department and its professors, I extend the greatest of thanks for their support and advice during the later parts of the thesis document and in the learning process in terms of the programming involved in the numerical realization of the special relativistic hydrodynamics code, written in Java. Second to the Counselors and Psychiatrists at the University of Northern British Columbia, for their incredibly important role in my mental well being and health. Also to their understanding of the need for the costs of my medications to be supplemented, due to my living in the poverty level. These people’s work is usually received with arrogance and ignorance, and I am here to say that I would not be here today without the acknowledgment and treatment of my mental health disorder. Third I would like to extend thanks and appreciation to the native community of one o f my best friend’s sweat family. These people accommodated me in times of need when academia could not and would not help. Without the support of my friends in the sweat family, I fear I would have turned down dark paths. Fourth I would like to extend greatest thanks to my best friend Jayme Fougere, who has nagged me to work on the thesis during periods of depression. Even when he could not convince me to do so, he got me out of the house to go hiking and just hang out. His influence on my physical well-being is paramount to my mental health. Fifth I would like to thank my mother, Ruby Mohammed, for her beating me over the head to learn and dedicate myself to education. She instilled a belief in education as a path to success and ix Acknowledgements x fulfillment, and this lead to the production of this thesis (albeit after much suffering and poverty!). Finally, but most assuredly not lastly or unimportantly, I would like to extend my thanks and eternal gratitude to my thesis supervisor and mentor, Dr. Patrick Mann, who taught me everything I know in the physics of relativity and of the technique of numerical work, and who also taught me how to function in Canadian society and not get into trouble. He has also supported me financially and now I owe him an arm and a leg. I consider him a most excellent and admirable hobbit, and the greatest of men. 1 Introduction Stellar models attempt to reproduce supernovae via core collapse processes which are created in numerical simulations. Recent advances in the understanding of supernovae processes as well as in computer hardware technologies are producing increasingly sophisticated multi-dimensional supernova simulations. Some o f the processes involved are neutron capture, radioactive decay, the different fusion branches and the neutrinos they produce, neutrino high-energy physics, neutrinomatter interaction and neutrino scattering. This thesis is a simulation of the behavior of stellar fluid undergoing neutrino heating in a shock tube. The model is spherically symmetric, using the usual special relativistic hydrodynamics equations. Within the shock tube Cartesian coordinates are used, and the shock tube is long enough that this system of coordinates is a good approximation. This work focuses on the electron-neutrino (where “electron-neutrino” is shortened to “neutrino” in this work, see Section 1.1 for the justification) flux and energy terms as source terms in these equations. The numerical method being used is the Godunov Finite Volume Method, utilizing an exact Riemann solver on the nodes (boundaries of cells). The exact Riemann solver is a general relativistic one, originally written in Fortran 77, by Pons, Marti and Muller ([1]). This thesis will show that this simple model, where the neutrino physics employed by Kuroda et 1 Sec. 1.1 Stars and Their Ways ... 2 al. ([2 ]) is reduced to simple expressions for the neutrino energy flux and neutrino momentum flux, produces results which are physically reasonable. One such result is that the net result of either a neutrino energy flux or a neutrino momentum flux, or both, acting on the stellar fluid (hereafter referred to as “fluid”) is a flow outward from the neutrino-sphere. The neutrino-sphere is the region within the star where the optical depth, r > | . This region extends out from the core of the star to the boundary where r = | . Another result is a confirmation of the approach to the source terms (neutrino fluxes). The approach is to sum the resultant fluxes due to the neutrino energy flux and neutrino momentum flux (see Shibata et al. ([3])), in order to produce a source term which adds to the fluid’s energy and momentum evolution equations. In the presence of both fluxes, it is expected that the fluid would be subject to the sum of those fluxes. This is observed in this thesis’ simulations. An exciting result is that the fluid energy generated by the inclusion of both neutrino fluxes acting together is on the same order of magnitude as that observed. This is significant as it shows that this thesis’ simple model has brought out at least some of the critical factors needed in order to produce an explosion. It has also brought to light that the active area of the star is in the region just outside the neutrino-sphere. This is where neutrinos emerge with the temperature they attained in the neutrino bath within the neutrino-sphere. Once outside the neutrino-sphere, these neutrinos are decoupled from the fluid, but have a large cross-section comparable to the elements of the fluid at the densities found there. So neutrino heating occurs, and neutrino cooling is minimal (an assumption in this thesis’ model). This increases the efficiency of the neutrino heating, and leads to explosion energies being attained. 1.1 Stars and Their Ways ... It is generally accepted that most of the elements of nature are created in stars and released to the universe upon a supernova. Supernova are therefore important, ultimately, to the existence of life. These explosions are incredibly powerful and can outshine a galaxy for weeks. The mechanism of a supernova is thought to depend on the transfer of energy just outside the neutrino-sphere. Neutrinos appear to be the mediators of this transfer due to the vast numbers produced by nuclear interactions during core collapse. In the neutrino-sphere, the neutrinos are coupled with the stellar Sec. 1.2 Notation and Conventions 3 fluid, and adopt their temperature. Outside of the neutrino-sphere, the neutrinos decouple from the stellar fluid and carry away that energy outwards. In the “gray” zone between the surface of the neutrino-sphere and the atmospheric regions, neutrinos may encounter the stellar fluid and transfer some of their energy to the fluid, thus heating the fluid. Neutrino cross-section is the property critical to neutrino heating. The species of neutrino with the largest cross-section is the electron-neutrino, ue. It interacts with the stellar fluid in different ways. The one of concern in this thesis is i/e scattering. This is represented in general as ur + x -» vc + x, where x is a proton, neutron, electron or nuclei. The opacity due to ve + e —> ue+ e scattering is large compared to all other processes (not just scattering), at low neutrino energies (e„ < 5 MeV) and high matter temperatures (Burrows and Thompson ([4])), and so is the major process to consider. Shock systems set up by the collapse may become energetic enough to drive the fluid outwards. This may produce explosions on the order of those observed. Modeling these scenarios depends on in-depth knowledge of subatomic physics, thermodynamics, relativity and quantum mechanics. 1.2 Notation and Conventions The notations and conventions used in relativity are applied in this thesis. Greek indices take values fj = 0 ,1 ,2 .3 . The components of the index start with time followed by the spatial components. In the example, fi = t , x , y , z for Cartesian coordinates. Latin indices (i,j,k, . . . ) denote spatial components. The signature applied is {—, + , + , +}. Partial differentiation along the time coordinate use the “dot” notation : 4-Vectors are denoted by an arrow over a symbol, such that x is a 4-vector. The usual threevectors are denoted by bold-face, such as x. A “,” denotes partial differentiation: = T a$ Ox** ~ 3 where T ni3 is a tensor. Also, (i 2) ( } Sec. 1.3 Introduction to the Godunov Method 4 a ( v v ) ; = ( V j? ) (1.3) which denotes the covariant derivative by the semi-colon in the lower index. Finally, there is the Einstein summation convention which is the expression of repeated indices. 3 = (1.4) - ,= o The Lorentz factor is defined as -?= z , where geometrized units are used such that c = G = 1. Appendix A .6 provides details of the geometrizations used in this thesis. 1.3 Introduction to the Godunov Method The terminology used in this thesis follows the conventions applied by Leveque ([5]). Lower case “q” is used to refer to primitive quantities, for instance a density on a node. A lower case “f” is used to refer to a flux term as a function of “q” . So, a conservation law can be written as, q{x, t) 2. Since the simulations produce shock reheating and re-acceleration Lorentz factors of this magnitude might occur. Nevertheless, this thesis is focused on a Godunov method with the employment of an exact Riemann solver as described in Chapter 1, Section 1.3. Sec. 2.3 2.3 Two Exact Riemann Solvers 10 Two Exact Riemann Solvers This thesis is concerned with solving special relativistic hydrodynamics using the Finite Volume Godunov method, which employs an exact Riemann solver put forward by Pons, Marti and Muller ([1]). There are two exact Riemann solvers which were explored in this thesis. These are the codes written by Pons, Marti and Muller ([1]) and by Rezzolla and Zanotti ([19]). These codes are both in Fortran, and had to be converted to Java in order to be used in the Godunov evolution code, which is also in Java. Both of these codes are excellent, with the Rezolla and Zanotti code being the superior, as it handles cases where the states on both sides of the jump in the Riemann problem are equal, and also the cases where the tangential velocity is not equal to zero. The code written by Pons, Marti and Muller does not handle state equality at all, and works only for tangential velocities equal to zero. However, due to Rezzolla and Zanotti’s code not being very well documented, and the fact that the Pons et al. code was converted to C++ prior to this thesis by the supervisor, the latter was used (since converting C++ to Java is easy). Sod did the early work on a set of initial conditions now referred to as the Sod Shock Tube, applying different numerical techniques to the solution of the problem in Newtonian non-linear hydrodynamics equations for an ideal gas. His results are the benchmark against which any code must be compared. The problem is in fact just the Riemann shock tube so that the known exact (Riemann) solution can be used for comparison. A similar test example for the relativistic shock tube was given by Pons et al. and is the benchmark for relativistic solvers. The equations of relativistic hydrodynamics can be be written in conservative form as (for the details see Chapter 3): <9fU + d t F ^ — 0 ( 2 .1 ) where U are the vectors of the conserved variables and F (l) are the fluxes. U = { D , S \ t) T ( 2 .2 ) Sec. 2.4 High Resolution Shock Capturing F w = ( D v \ S Jv l 4 pSJ\ S l - D v l) r - 11 (2.3) The Pons et al. ([1]) code uses arbitrary tangential velocities in its solution of the Riemann problem, and focuses only on the modulus of v* and not the direction of the tangential velocity. This is where the Rezzolla and Zanotti ([19]) code differs from the Pons et al. code. In the former, the different possible directions of the tangential velocities are taken into account: 1 . two shock waves, one moving towards the left initial state and the other to the right initial state, 2 . one shock wave and one rarefaction wave, the shock moving to the right and the rarefaction to the left, and 3. two rarefaction waves, one to the left and the other to the right. The details are given in Rezzolla and Zanotti ([19]). 2.4 High Resolution Shock Capturing The application of high resolution shock capturing methods caused a revolution in the field of numerical SRHD. Prior to this, methods used by Wilson ([6 , 7]) and others ([ 8 , 9]) applied artificial viscosity to make their first-order Eulerian methods stable. However, this smoothed the shock systems sometimes to such a degree that the information of the original system was lost. The artificial viscosity approach also limited the simulations to those with Lorentz factors < 2.0. The Godunov method itself also yields smoothed results. A high resolution method isabsolutely necessary to extract the details of the solutions. Due to this limitation, it was desirable to find second-order or higher order methods which could preserve the shock systems and so develop a relativistic model of the the solution to the SRHD equations. The basic properties of any such method are: • high order of accuracy, • stable and sharp description of discontinuities, and Sec. 2.5 Concerning Neutrinos 12 • convergence to the physically correct solution The fluxes at the zone interfaces where the HRSC is applied is usually found by either an exact Riemann solver or by methods which approximate them based on the solutions on either side of the interface. Due to the computational expense of the exact Riemann solver, approximate solvers are usually used. Approximate solvers are faster than exact Riemann solvers and are not as computationally expensive. The exact Riemann solving approach has the advantage that no further work is necessary once the fluxes have been calculated using the exact solutions from the exact Riemann solver. The approximate approaches need to apply corrections and methods to ensure that their solutions converge. This adds extra computations. Because o f this, it was decided to use the exact Riemann solver, in this case the one by Pons et al. ([1]). 2.5 Concerning Neutrinos The mechanism of core-collapse supernovae is thought to depend upon the transfer of energy from the core to the mantle o f the inner regions of a massive star (8 — 12M0 ) after it becomes unstable to collapse. It appears that neutrinos are the primary carriers of this energy transfer. One approach to understanding the role o f the neutrinos and the particulars o f the mechanisms of the energy transfers is to consider neutrino cross-sections for interaction with the stellar fluid. Tubbs and Schramm ([20]) and Bowers and Wilson ([17]) are excellent resources for the neutrino crosssections of the three flavors of neutrino type, and for the development of the particulars of the neutrino interactions. However, much work has been done since Tubbs and Schramm, and Bowers and Wilson ([20, 17]). Burrows and Thompson ([4]) provide a thorough summary of work done with neutrino cross-sections and opacity. This thesis was initially going to follow these works and model the neutrino interactions using cross-sections and special equations of state for the various types of interactions. However, at a late stage a switch in the model has been made, with primary reference to Kuroda et al. ([2]). Kuroda et al. ([2]) use a M l closure (this is just an equation which closes the set of neutrino equations and is of a form which depends on the Lorentz factor (Levermore ([21]))) to solve the Sec. 2.5 Concerning Neutrinos 13 energy-independent set of radiation energy and momentum based on Thorne’s ([22]) momentum formalism. In that paper, two moment formalisms were presented, with the pertinent one being that which applies to spherical symmetry. Shibata et al. [3] extends this formalism by providing a truncated moment formalism for the radiation hydrodynamics. Kuroda et al. make a nice separation of the fluid and radiation tensors, like this, + V „ T “» (2.4) Since V nT l',i3 = 0 this can be written as, V -T S L , = -Q '1 = Q" (2.5) (2.6) where Q !i are the source terms that describe the exchange of energy and momentum between the fluid and radiation. This thesis will beconcerned with only the absorption heating). of neutrinos by the fluid (neutrino The derivation of thespecial relativistic neutrino fluid equations and the simplifications made in this thesis to the general relativistic work of Shibata et al. are shown in detail in Chapter 3. 3 The 1-D Special Relativistic Hydrodynamics Equations for the Stellar Fluid including the Neutrino Model This chapter presents the thesis model in its entirety, from the development of the special relativistic fluid hydrodynamics to the development of the special relativistic neutrino hydrodynamics and source terms. The fluid equations are obtained by working in general relativity first, with the relevant simplifications applied to reduce the equations to flat space. Cartesian coordinates are used in this case. The model is 1 -dimensional, which further simplify the equations. This neutrino transport model is based upon a truncated moment formalism for radiation (neutrino flux) hydrodynamics. This is a formalism where a set of covariant equations are defined from the distribution function of radiation, and an approximation is made where the higher order moments can be neglected. Together with a closure relation, the causal relation can be preserved, and a solution of the radiation transfer in the optically thick and optically thin regions can be derived. 14 Sec. 3.1 3.1 I-D Special Relativistic Hydrodynamics 15 1-D Special Relativistic Hydrodynamics Consider a cube of some volume immersed in a fluid flow (see Figure 3.1). The amount of mass entering the volume must equal the amount of mass leaving the volume. The rest density, p0, is referred to as the rest mass, and ua is the velocity field. u„ is the four-velocity, and satisfies the relation unua = —1. The three-velocity is then defined as \ l = p-. Not ctiang>e of mass across iti'ovolume is zero Fluid m a ss flcw n a in FUd ftvass ci 4- V This implies that Equation (3.26) can be decomposed as, = 0 (3.26) Sec. 3.2 1-D Special Relativistic Neutrino Hydrodynamics 19 (3.27) (3.28) where Q 0 represents the source terms describing the exchange of energy and momentum between the fluid and the neutrino radiation. Conservation of mass also needs to be addressed. The contribution o f neutrino mass to the fluid is ignored. Similarly the fluid does not create neutrinos and therefore the fluid does not lose mass. These approximations are reasonable in the outer regions of a supernova where neutrinos are free-streaming (see below for more details). Thorne’s work ([22]) uses the ADM formalism ([23]) (the 3 + 1 split, where the metric is ds 2 = ( - a 2 + 3, ft1) dt2 + 2Qidtdxi + 7 ijd xldx3). The 3 + 1 split allows for numerical work in general relativity to be carried out using time slices. Each time slice is a surface, and the transition from one time slice to the other time slice corresponding to the next time step involves a lapse function and a shift vector to take into account the general relativistic effects in such a transition. The 4-vector normal to the surface (time slice) is n a, and is defined by the lapse function a and the shift vector f3\ written as n a = . Since the 4-Dimensional surface is split into a 3-surface at some particular time, then a three metric is needed to describe the shape of the 3-surface. This is defined as yag = gai3 + nQng. The following uses a truncated moment formalism, which was introduced in the work of Thorne ([22]). Start by decomposing the neutrino energy-momentum as, (3.29) where E („) is the radiation internal energy, F („) is the radiation flux and _P(„) is the radiation internal pressure. (3.30) (3.31) (3.32) Sec. 3.2 1-D Special Relativistic Neutrino Hydrodynamics 20 Kuroda et al. ([2]) develop a system of energy-momentum evolution equations for radiation energy and radiation flux in a general spacetime. Here, for simplicity, a flat space is used to derive the equivalent equations. Recalling Equation 3.28, V a7 $ = Q& In flat space, na — { 1 ,0 ,0 ,0 } and (3.33) is the projection of the energy-momentum tensor into the spatial slice, so the t component vanishes, that is, Ffo — {0, F ^ } . The same reasoning holds for the radiation pressure, P^fy P ^ is composed of Pfo = 0 , P ^ = /(F = 0 and /{{I is non-zero. The covariant derivative V , tT (aJj is now just a partial derivative, 0 „ T ^ . Using Equation (3.33), and setting ft = t, then, T $ ,Q = Q* (3-34) [F ((,}n V + F ^ n 1 + F{u)n a + P g ] a = Q l (3.35) Applying the flat space restrictions gives, [£ M » " + % ] ^ = Q< (3.36) Now expand the a sum to get, (3.37) => E,( u ) X (3.38) For (5 — i, '-pai (3.39) l {u) [E(U), P E + F^n> + F(U),P + / $ ] Applying the flat space restrictions gives, = Q' (3.40) Sec. 3.3 Prelude to the Neutrino Model 21 (3.41) Now expand the a sum to get, (3.42) (3.43) To summarize, Equations (3.38, 3.43) are the evolution equations for the energy and momentum of the neutrino fluid. However, there are only two equations and three unknown quantities (F, F \ P ij). Another equation is needed to uniquely define the system. With most fluids, such an equation is an equation of state for the fluid. Levermore ([21]) investigated approximate methods to study transport phenomena, and in so doing developed a neutrino equation of state which is applicable to this work: pij ^ (tA (u) _ 3\ 2si 1 pij * ihirt thin '' 3 (1 r\ 2 x) p ifh:] thick (3.44) The development of this equation is not discussed here; it is only presented as a demonstration of how Equations (3.38, 3.43) can be closed. The evolution model for E v and Fv is not employed in the model used in this thesis. Instead, the E u and Fv required are calculated at the surface of the neutrino-sphere (discussed in Section 3.4). 3.3 Prelude to the Neutrino Model This section uses work described in Kuroda et al.’s paper ([2]). Here the source terms, and 7^ need to be defined. The electron-neutrinos are thought of as being in two kinds of regions: trapped and free-streaming. The trapped region occurs inside the neutrino-sphere, and is of no concern in this thesis, as the shock tube begins at the surface o f the neutrino-sphere. The free-streaming neutrinos are in the optically thin region, and provide the neutrino heating. The neutrino source term, Q11 is composed of the neutrino cooling term, Q 11c and the neutrino heating term, Q ^ H. Kuroda et al. ([2]) develop a system of evolution equations (which are not Sec. 3.3 Prelude to the Neutrino Model 22 implemented in this thesis) where source terms appear composed of the neutrino cooling and neutrino heating terms. - Q T l = - (Q I,C - < T " ) l,n (3.45) Q "n l l = { Q ll-c - Q l,-, l) n ll (3.46) This thesis ignores the cooling term, so the source terms are really the neutrino heating terms. These are developed by Kuroda et al. to give, -Q ^i = (esJ 2 KVe { - W F VrJ + P V kctu k) Q W = c r ^ ( r , , J 2 s Vc ( W E Vr - F * u k) (3.47) (3.48) Here, t Svr is the neutrino energy (remember the thesis only deals with electron-neutrinos). The expression for the opacity, k 8c = (e„) 2 K, is used from the work of Janka ([24]). The subscript sc represents “scattering”, since the only process of concern here is thetransfer of energy and momentum by electron-neutrino scattering by interaction with free nucleons. Janka develops the scattering term to be, 24 where, rv = - 1 .2 6 , Yn + Yp « 3.69 x 10 57 km. m cc2 (m ec2)2 m u <349) 1, vc + ve (Tubbs and Schramm ([20])). This thesis uses the approaches developed by Matteo et al., Liu et al. and Zhang and Dai. The neutrino pressure is given by, P<,) The neutrino energy density , (3-64) = U(„), is given by (Matteo et al. ([25])): Tw | i 8a7f- [ 2 ^ V3 Hv) ~ IM , J_ , _J_ 2 t (3.65) l / 3 1 ' 3 r(l/) With the optical depth T(„) = |2 at the surface of the neutrino-sphere: E oT (4 \ [7 + 7v/3] ------ — 1 = — w 20 + 8 \/3 (3.66) where a is the radiation constant , which is a = 7.5657e-15 erg cm - 3 K-4 . o is the StefanBoltzmann constant and a — 5.6704e^5 erg cm - 2 s~x K "4. The neutrino flux, F(v) is given by (Matteo et al. ([25])): Sec. 3.4 ^ The Neutrino Model 7aTM = 5 ^ 25 <3'68' The temperature of the neutrinos, T^), is defined by, T(u) = E{v)/ks (3.69) since the neutrinos are decoupled (this thesis assumes that there is total decoupling immediately at r- < I) from the stellar fluid outside the neutrino-sphere the temperature is constant. Using the values found earlier, then T(„) = 1.74 x 10n K. At this point it is good to show a visual representation of the supernova model. SHOCK TUBE rSTELLAR ATMOSPHERE NEUTRINOSPHERE OKM 1000 KM Figure 3.2: The Position of the Shock Tube in relation to the Neutrino-Sphere. Sec. 3.4 The Neutrino Model 26 The neutrino-sphere is the region where T(„) > §. Using the graphs in the work done by Kuroda et al., it was determined that the best position to place the Riemann shock tube for the purposes of this thesis was on the surface of the neutrino-sphere. This is the only location which, according to Kuroda’s data shown in Figure 3.3, would capture a shock in the arbitrarily chosen time of 40 ms (see Figure 3.5). The time of 40 ms was chosen as it was the longest interval for which good data from Kuroda’s graphs could be extracted. The neutrino-sphere’s surface corresponds to .r = 800 km. The length of the shock tube is 1.8 x 102 km. Based on the previous discussion, T(„) has been calculated. In this thesis, using Equations (3.66, 3.68), the E(„) and F (,/) can be computed and taken as constant along the shock tube. The data of Kuroda et al. can be used to obtain the values for the velocity profile, the density profile and the neutrino specific energy at time = 10 ms to time = 40 ms. Figure 3.3 was used to determine the length of the shock tube and its placement by using the time span obtained from Figure 3.5. In Figure 3.5 the time of 40 ms was chosen because it intersects the energy profile for a 1dimensional special relativistic simulation at an energy o f 15 MeV which, according to Burrows and Thompson ([4]) is in the range of energies for electron-neutrinos in the core collapse. Also, the velocity profile, Figure 3.3, only extends to 37 ms so no greater time data could be utilized. Knowing the time constraint and the location of the neutrino-sphere as detailed in Kuroda et al. ([2 ]) then using the velocity profile allowed the determination of the length of the shock tube and the location of the initial shock (as midway between the neutrino-sphere and the location at a time of 37 ms). Sec. 3.4 o> O X o X .................................. ho4 ■ I io5 io6 I. I l l I ■I HI x (cm) Figure 3.3: The velocity profile o f the stellar fluid in the shock tube (Fig. 16 [2\). Sec. 3.4 1 ■i m m T I » I M ill I I I I III K) i m i mi i ii n i Neutrino Sphere 10 14 1 0 »' The Neutrino Model 13 s1 0 1 E10 n- o 10 o >10 o. 1 0 ' 8 -*> r ^ 2 x l 0 Al l g c m A-3 r SR r GR 1 0 1 10 ' 1 0 * 1 0 IQ - 1 0 ' x 1 0 7 1 0 ® (cm) Figure 3.4: The density at the surface o f the neutrino-sphere (Fig. 15 [2]). Sec. 3.4 The Neutrino Model > a> tfi E __ — A V ° 3DGR 1 DGR 3DSR 80 100 Figure 3.5: The neutrino energies at different locations in the shock tube (Fig. 13 [2]). 29 Sec. 3.5 3.5 Summary of Equations 30 Summary of Equations The relevant equations which need to be evolved and used in the calculations at each time step are summarized below. The fluid equations are given with the associated source terms. D,t + {D vl),t = 0 (3.70) E,t + ((E + p y ) , = (3.71) S \ + (.S V ) (3.72) = -Q * The source terms are, in geometrized units, and W = y/l i—v 2 ’ Q l = 1.69 x 10Vo [ W E (V) - F(l/)v) (3.73) ~ Q X = 1.69 x 104po { - W F {v) + PH v) (3.74) The quantities P ^ , E ^ and EtlJ] are given by the equations, E = 1 } aTfv) [7 + 7n/3] J 20 + 8 \/3 P{u) = (3.75) (3.76) 7aTU F«">- I r i k (3 '77> T (u) = ^ (3.78) This is the complete set of fluid and neutrino equations implemented in the thesis code in order to simulate a 1-D special relativistic core collapse supernova. 4 The Godunov Method In general, gas dynamics equations are of a type similar to, Qt + g{q)x = o (4 -i) This general case is non-linear and discontinuous solutions may exist. In particular, Equations (3.70, 3.71, 3.72) in Chapter 3 are in the form, Qt + g(q)x =< source > (4.2) which are generally non-linear and have discontinuous solutions. In this chapter the Finite Volume Method is introduced. Finite volume methods are closely related to finite difference methods; however, finite volume methods are derived from the integral form of the conservation law. The upwind method is a finite volume method and is in fact the linear version of Godunov’s method. The discussions presented next follow this sequence. The Finite Volume Method is presented, then the Upwind method followed by the Godunov method which is used in this thesis. 31 Sec. 4.1 4.1 The Finite Volume Method 32 The Finite Volume Method Differential equations assume continuity, which means that discontinuities are explicitly excluded. A numerical method which approximates the differential form of the evolution equations is expected to break down in the presence of discontinuities. If the integral form of the equations are used, then shocks are expected in the system, and a numerical method developed based on the integral form would be able to cope with those discontinuities. Such a method is the Finite Volume Method. 0 FLUX O SPACE TIME t(n+2) t(n+1) F (n,i-l/2) t(n) Q(n,i-1) F |n,i+l/2) Q(n,i) Q(n,i+1) Figure 4.1: A representation of the grid-based nature of the Finite Volume Method. Sec. 4.1 The Finite Volume Method 33 Figure 4.1 is a representation of the grid-based nature of the Finite Volume Method (FVM). The vertical axis is time and the horizontal axis is the spatial x dimension. The following discussion is a summary of the discussion by Leveque ([5]). An /th grid cell is defined as, (4.3) C i = (.Ti_ i , x i+i ) Let Q f, be an approximation of q on a cell C,. Then with A x = x l + 1 - x,_ i write, rx+i Qi ~ f A * J x- \ q(x. tn) dx = -^~ f q(x, t„) dx A x Jc . (4.4) Integrate over a cell C, to get: (4.5) ^1(" ddx+zb ( f (" ith I i t i Q, ) ‘ “ 1 ( « (x< - *))) - =0 (4.6) (4.7) 0 Now integrate over t: I * 7 t { Q ,) ‘ i t + C g ) + = 0 (4-8> <4 .9 ) o?“ = or - £ ( / («(*»{.«)) - / {i (*<-}•*))) rf‘ <4,10> Then define the flux approximation by: <4-" > which means that Equation (4.10) is now written as: o r ’= o ? - ^ ( c i - ^ i ) <4.i2) Sec. 4.1 The Finite Volume Method 34 The FVM approximation of the conservation equation in 1-dimension has now been derived. The remaining task is to consider approaches to calculating the integrated flux, F. F is the approximation to the flux on the boundary of a cell, which is the flux on the nodes. This explicitly satisfies conservation as the flux flowing out of one cell is the flux flowing into the adjacent cell. 4.1.1 Methods of Calculating the Flux There are various ways of calculating the flux in Equation (4.12). A first attempt would involve ) + /« ? ;■ ) ] <4.i3) Using this in Equation (4.12) gives, A f Q7+l = Q 7 - ^ [ / (Q?+1) - / (Q?-0 ] (4.14) This method is the simple, centered in space, forward differenced Euler scheme , and is well known to be unstable. For hyperbolic problems (where a linear system of the form qt + Aqx = 0 is referred to as hyperbolic if the n x n matrix A is diagonalizable with real eigenvalues (Leveque ([5]))), information propagates as waves moving along characteristics, where a characteristic is a curve along which a partial differential equation reduces to an ordinary differential equation (Haberman ([28])). For a system of equations there are several waves propagating at different speeds and possibly different directions. “Upwind” methods are those whereby the information for each characteristic is obtained by looking in the direction from which the information is arriving. As an example (Leveque ([5])), consider the advection (where advection refers to transport of a substance by a fluid) equation qt + uqx = 0, where u is a constant, and let u > 0. The u < 0 case is similar. The upwind method for advection is illustrated in Figure 4.2. Here it can be seen that if Q" represents a variable value at a grid point, then the characteristic can be traced back and an interpolation be done to yield flux as, - TiAt. This suggests defining the numerical Sec. 4.1 The Finite Volume Method 35 Q(n+1, i) t (n+1) t ( n) Q (n, i-1) Q (n, i) x, - u At Figure 4.2: Upwind Method for Advection. F?_x= uQ U 1 (4-15) 2 This leads to the standard First-Order Upwind Method for the advection equation, Q7+1 = Q i - ^ ( Q ? ~ Q i - 1 ) (4-16) This exercise has shown that the advection equation has only one characteristic; the velocity u, and in only one direction, left being negative and right being positive (in 1-Dimension of course). Using this characteristic and dividing the grid into cells, a numerically implementable equation can Sec. 4.2 The Godunov Method 36 be written where knowledge o f the characteristic leads to knowledge in the upwind direction of the next value of Q. In order to maintain stability in this method, the Courant condition (a criteria for stability, defined as ^ (Leveque ([5]))) must lie between 0 and 1, that is, 0 < ^ < 1. For a system of conservation equations there will be many waves and characteristic velocities. Each wave needs to be suitably upwinded to guarantee stability. The Godunov method is a way of identifying each wave and upwinding it. 4.2 The Godunov Method A Riemann problem is a hyperbolic equation together with piece-wise constant data with a single jum p discontinuity at some point, say x = 0 (Leveque ([5])). This data is expressed as: (4.17) Given a continuous function q(x, t), the Godunov method imposes local Riemann problems at the boundaries between cells and the intersection of q(x. t) with the cell boundary. These are the boundaries between the horizontal lines in Figure 4.3. The horizontal lines represent the averaged quantity of q(x, t), Q (x, t) in the “cell” which is the region captured by the boundaries. The flux across the boundaries is a function of the Riemann solution on the boundary, as seen in Equation (4.18). Using an exact Riemann solver means that an exact value for the flux F n x is obtained. Let J+ 2 q1 (Qj n , Qj) be the exact q from the Riemann solution. It is dependent on the jump between Q 's from the left and right cells. / {Qj+i, Q j)) is the flux as a function of the Riemann solution, ql , using the values of the state on the left, Qj, and the state on the right, Qj+i, of the cell boundary (the “node”). (4.18) The implication holds true because ql is constant along rays in the Riemann solution, and therefore Sec. 4.2 r a l e l .i c u o i i shock ;fv ; [ 37 d iic o iitin u ily ; jj/f/'-'Y' : jli/f \ Vvj\ ' / Wi 1 / '& 1 / c o n ta c e The Godunov Method ' Hit/ \ I {!$ ; \ 1 Ii : \ \ \ \ I ! <•/// ' !,:// [!$ I X ) .| X j . |/ J X) X jt j .2 Xjf 2 co n tin u u m solution / _.- diarrelc so&tlioci f* ' \ xyi X jf | X| r................. *«♦, x j+ | X j+3 Figure 4.3: The Godunov method and Riemann problems. the integral can be evaluated. So the Godunov method can be summarized for a general system of conservation laws (Leveque ([5])): • Solve the Riemann problem at x t_ i to obtain ( Q f -i, Q?) • • Define the flux FP_ j_ = T (Q"_i, Ql) by Equation (4.18). • Apply the flux differencing formula. 4.2.1 Godunov Algorithm At this point it is useful to present an algorithm which will lead to a numerical implementation of the Godunov method. • create grid • set initial data • for each time step do - for each node do * compute Riemann solution between left and right cells to get q4 Sec. 4.3 High Resolution Methods 38 * compute flux from q - end do - for each cell do * evolve cell averages from the flux on the nodes (explicitly conservative) - end do • end do 4.3 High Resolution Methods The Godunov method for the advection equation is stable, but has poor accuracy due to large amounts of smearing, that is, smoothing of any discontinuities. As such, some method of sharpening the solution is needed. This leads to the need for High Resolution Methods , which are generally in the form, The F ’s are calculated from the Riemann solutions. The other term is the diffusion term, where are the left- and right-going waves. If the diffusion is negative (an anti-diffusion), then the numerical solution is sharpened. The amount of sharpening depends on the method used to derive the diffusion term. It turns out that if the Courant condition is satisfied, then |u | is positive, and there is anti-diffusion, resulting in a solution which approximates the exact solution to a high degree of accuracy (Leveque ([5])). The anti-diffusion term is more appropriately written as a flux which corrects the Godunov method. One requirement for this method is a methodology to control the amount of anti-diffusion. This is achieved using a limiter that changes the amount of correction used, depending on how the solution is behaving. A general form of how the flux can be limited is, F 1n 2i — (Qi-i.QO + 2 [Ph ( Q i - u Q i ) - F l ( Q i - u Q i ) ] (4.20) Sec. 4.3 High Resolution Methods 39 where the subscript II refers to a high-order (sharp) method and the subscript L refers to a loworder (smooth) method. If 4>" 2 = 0 then a low-order method is obtained and if <£>"_ t = 1, 1 2 1 2 a high-order method is obtained. So adjusting 4>n_i between 0 and 1 determines the degree of 1 2 accuracy. An interesting approach to the computation of high-order anti-diffusion terms is to use highorder fitting to the cell values to estimate the jum p between cells. The correction can be implemented as a linear manipulation of the solution at boundaries where there is a significant change in the solution on the left and the right of the boundary. Such a manipulation is a piecewise linear reconstruction , o f the form, qntn) = Ql + °1 Q - ■*',) • for x , < x < x i+}_ (4.21) where, = \ (•'' j + - ' m ) = x i-h + l Ax (42 2) - m « - <■) <4-»> Godunov’s method with anti-diffusion is written, « ? " - or - $ (or - o r . , ) - ^ It can be readily observed that choosing slopes erf = 0 gives Godunov’smethod. To obtain second-order accurate methods the slopes should be non-zero in such a way that erf approximates the derivative qx over the ith grid cell. Three possibilities are, = 0!±1 - _ g L . ( 4 .2 4 ) 2 A.c <4.25, Various combinations of Equations (4.24,4.25,4.26) can be used as long as the slope is limited to smooth out oscillations and instabilities. Sec. 4.3 4.3.1 High Resolution Methods 40 The CFL Condition The CFL condition (which is the acronym for its orginators, Courant, Friedrichs and Lewy) is a necessary condition that must be satisfied by any FVM in order for it to be stable and convergent (see Leveque ([5])). It is a statement that the method must be used in such a way as to allow information to have a chance to propagate at the correct physical speeds. The CFL condition is usually expressed as a ratio called the Courant number, which is a necessary but not sufficient condition for stability of the scheme. The Courant number is, ^ ( 4 27) Ax where smax denotes the maximum wave speed in that grid cell with intervals A t and Ax. Irrespective of the numerical techniques employed, an important condition which must be satisfied is that of convergence. Key to the issue of convergence is the Lax-Wendroff Theorem T heorem 1. Consider a sequence o f grids indexed by j = 1 , 2 . . . . , with mesh parameters AtSJ\ A x (j> -» 0 as j —¥ oo. Let ( f j ) {x, t) denote the numerical approximations computed with a consistent and conservative method on the jth grid. Suppose that ( f J) converges to a function q as j —> oo, in a manner made precise in [5, pages 240-243]. Then q(x, t) is a weak solution o f the conservation law. Theorem 1 applies equally well to conservative methods for non-linear systems of conservation laws as well as to scalar equations. If a sequence of numerical solutions converge to a function q(x. t) as the grid is refined, then that function q(x, t ) must be a weak solution of the conservation law. Thus, because of Theorem 1, one can have confidence that the method is converging to a valid solution. Since it is not desirable to work with nonconservative methods in the first place, then one may be tempted to say that the Lax-Wendroff theorem is sufficient for convergence. For this hydrodynamic problem, there are physical constraints which need to be obeyed. One of these is that of entropy and the part it plays in the physics of the problem. This condition is actually quite simple to understand; due to the second law of thermodynamics, the total entropy can never decrease in a system. Thus, in each cell of the grid, the total entropy of that system Sec. 4.3 High Resolution Methods 41 cannot decrease. The use of the word “total” is intended to mean that the average entropy in a cell is constant or has increased in the next time-step. This condition implies that the entropy can be represented by a conserved function, say //( 0 b if \b\ < |a| and ab > 0 0 if ab < 0 (4.32) The minmod compares the two slopes and chooses the one that is smaller in magnitude. If the slopeshave different sign, then thevalue o f Q f must be a local minimum or maximum, and so erf = 0 in order to satisfy TVD (Leveque ([5])). A better choice of limiter which gives second-order accuracy as well is the superbee limiter, defined below: erf = maxmod (a-, of) (4.33) where, Here each one-sided slope is twice the other. The maxmod function in Equation (4.33) selects the larger modulus. In regions where the solution is smooth, this will return the larger of the one-sided slopes. Sec. 4.4 4.4 Summary 43 Summary This chapter presented the Finite Volume Method, in particular the upwind method and Godunov’s method, which is the numerical method used in this thesis. The reconstruction is a linear one employing a minmod limiter. A few steps of the algorithm are: • Calculate the slope sigmas using minmod. • Calculate the left and right primitive variable values using linear reconstruction. • Solve the Riemann problems across the node. • Use the exact Riemann solutions to compute high-order fluxes on the node. 5 Simulation Results and Analysis The code is tested using the classic test of implementing the Sod shock tube (Sod ([29])). This test does not include neutrino source terms. If the Sod shock tube simulation yields results which conform to Sod’s results within acceptable limits, then there is confidence that the simulation code is correct. In the present work an exact Riemann solver, created by Pons, Marti and Muller ([1]), is used to compare the approximate solution from the evolution code against the exact solution obtained from the exact Riemann solver. A quantitative test using the exact Riemann solution is implemented where standard deviations are calculated and plotted against resolution. Once the Sod shock tube results are verified, then it is known that the Riemann solving code is working properly. At that point the full evolution code with neutrino source terms can be run. Implementing the Riemann solver at all times resulted in very long runtimes. Instead a selecting mechanism is employed which selects the Riemann solver only when a specified minimum jump is detected. Otherwise an upwinding method is used to calculate the fluxes. This greatly reduces the runtimes to manageable values on the order of seconds. Once these issues were addressed, several runs were made. 1. Classic Sod Runs 44 45 2. 400 Cells : • The hard-coded Sod data was used to make runs o f the exact solution against the approximate solution for 400 cells. • A plot of the standard deviations obtained for this run is also shown in the graph. • Standard Deviations of (1-3)% were obtained. 3. 1600 C ells: • Standard deviations were reduced to (% 1)% 4. Runs were done left to right and right to left to ensure the Riemann solver was consistent in both directions. 5. Also, tests were done using the Sod shock tube for velocities approaching the speed of light. As the runs were carried out it was found that the Riemann solver becomes increasingly inaccurate as the velocity becomes a significant fraction of the speed of light. A velocity o f 0.65 c was chosen as the break-point for acceptable error. This is correct behavior, as velocities should never be allowed to exceed the speed of light. According to Chapter 2, this would mean that this thesis’ evolution code is no better than codes using artificial viscosity. In fact, there is debate about this finding. Some argue that there is no need to use exact Riemann solvers at all. One of the things this thesis will show is that such an approach will miss good results. 6 . The Kuroda et al. data, with the Neutrino terms off. This was done only with 400 and 1600 cells at a runtime of 750000 km. 7. The Kuroda et al. data, with the Neutrino terms on. This was also done only with 400 and 1600 cells at a runtime of 750000 km. 8 . There were two cases. • The parameter file contains a parameter called the energy tuner. This is a variable which can be used to set the value of the neutrino energy flux used in the evolution of Sec. 5.1 Sod Shock Tube Tests and Results 46 the fluid dynamics. This is simply a dimensionless number. If it gets beyond a zeroth order of magnitude, something is wrong! The energy tuner was set to 1.0 x 10 ~ 2 with the flux tuner zeroed. • The parameter file contains a parameter called the flux tuner. This is a variable which can be used to set the value of the neutrino momentum flux used in the evolution of the fluid dynamics. This is simply a dimensionless number. If it gets beyond a zeroth order of magnitude, something is again wrong! The flux tuner was set to 1.0 x 10" 2 with the energy tuner zeroed. 5.1 Sod Shock Tube Tests and Results The test consists of a one-dimensional Riemann problem for an ideal gas: 1 .0 Pi Pi = -I'l 0 .0 0.125 Pr Pr VT 1 .0 = 0 .1 0 .0 where, pi is the fluid density to the left of the location of the discontinuity and pr is the fluid density to the right of the discontinuity. Likewise, pt is the fluid pressure to the left and pr is the fluid pressure on the right, and vt is the fluid velocity on the right and ly is the fluid velocity on the right. Sod used initial conditions in a “tube” of some length, x. The jum p was located at a position midway along x. Then the simulation was run out to the total time, T = 0.25 s. The classic results are shown in the following graphs of p(), p, u and v, where the green graphs are the exact solutions and the red graphs are the approximate (evolution) solutions from this thesis code. The exact solutions are known because there is an analytic solution to the special relativistic Riemann problem. The Sod results using the Sod data are also known due to Sod’s work ([29]). Sec. 5.1 Sod Shock Tube Tests and Results 47 The blue graphs are the standard deviations for each primitive variable. Runs were made of 100, 200, 400, 800, 1600 and 10000 cells for a total runtime of 0.25 s. Only the 400 and 1600 resolution results are shown through Figure 5.1 - Figure 5.8. Sec. 5.1 Sod Shock Tube Tests and Results 48 tm im tm m m Approximate Solution Exact Solution Difference The Fluid Rest Density (rhoO) 0.8 0.6 llltlillllllllltllllllitllllllllllillM lttl 0.4 0.2 -0.3 -0 .2 -0.1 0 0.1 0.2 Length of Shock Tube (km) Figure 5.1: The plot of fluid density when applying the Sod data with a resolution of 400 cells. 0.3 Sec. 5.1 Sod Shock Tube Tests and Results 49 1 jmHMHUIHUHI Approximate Solution Exact Solution Difference The Fluid Pressure 0.8 0.6 0.4 IliH liH filttilliH lilllH IIH IilU ililliU lillH IIM IIIH H IIlIiH lllttlH IIII 0.2 iiiililllllllllll 0 *i -0.3 - 0.2 - 0.1 0 0.1 0.2 Length of Shock Tube (km) Figure 5.2: The plot of fluid pressure when applying the Sod data with a resolution of 400 cells. 0.3 Sec. 5.1 Sod Shock Tube Tests and Results 50 1.6 Approximate Solution + Exact Solution Difference — *■ 1.4 The Fluid Internal Energy 1.2 1 IIHIinUltUlllllll! 0.8 0.6 lilH H H I H K II I tlltillilM I I I H im if llH U I 0.4 0.2 0 -0.3 -0.2 -0.1 0 0.1 0.2 Length of Shock Tube (km) Figure 5.3: The plot of fluid internal energy when applying the Sod data with a resolution of 400 cells. 0.3 Sec. 5.1 Sod Shock Tube Tests and Results 51 0.45 ItllllilliH ltlllllllH H IIIN IIH IttllllllH IM M IH H H IIIIH IIlH IIH llI 0.4 The Fluid Velocity 0.35 0.3 0.25 Approximate Solution Exact Solution Difference 0.2 0.15 0.1 0.05 0 -0.3 -0.2 -0.1 0 0.1 0.2 Length of Shock Tube (km) Figure 5.4: The plot of fluid velocity when applying the Sod data with a resolution of 400 cells. 0.3 Sec. 5.1 Sod Shock Tube Tests and Results 52 1 Approximate Solution Exact Solution Difference The Fluid Rest Density (rhoO) 0.8 0.6 0.4 0.2 - 0.2 - 0.1 0 0.1 0.2 Length of Shock Tube (km) Figure 5.5: The plot of fluid density when applying the Sod data with a resolution of 1600 cells. 0.3 Sec. 5.1 Sod Shock Tube Tests and Results 53 1 Approximate Solution Exact Solution Difference 0.8 0.6 0.4 0.2 -0.3 - 0.2 - 0.1 0 0.1 0.2 Length of Shock Tube (km) Figure 5.6: The plot of fluid pressure when applying the Sod data with a resolution of 1600 cells. 0.3 Sec. 5.1 Sod Shock Tube Tests and Results 54 1.6 Approximate Solution + The Fluid Internal Energy 1.2 1 0.8 0.6 -0.3 - 0.2 - 0.1 0 0.1 0.2 Length of Shock Tube (km) Figure 5.7: The plot of fluid internal energy when applying the Sod data with a resolution of 1600 cells. 0.3 Sec. 5.1 Sod Shock Tube Tests and Results 55 0.45 0.4 The Fluid Velocity 0.35 0.3 0.25 Approximate Solution Exact Solution Difference 0.2 0.15 0.1 0.05 0.3 -0.2 - 0.1 0 0.1 0.2 Length of Shock Tube (km) F igure 5.8: The plot of fluid velocity when applying the Sod data with a resolution of 1600 cells. 0.3 Sec. 5.1 Sod Shock Tube Tests and Results 56 A graph of the error is a much better analytical tool, and this is obtained by finding the standard deviation across the whole grid to estimate the average error for a particular resolution. Six resolutions were used to generate six data points with which to plot Figure 5.10. The order of the error is given by, A q « a A x 11 (5.3) where n is the order of the error , q is the variable and A.r is the step size. Taking the logs of this equation will give the order of the error as, log A q = log a + n log A x (5.4) This plot is shown in Figure 5.10, and the order can be seen to decrease as the resolution increases. Sec. 5.1 Sod Shock Tube Tests and Results 57 Resolution^} sdv(2) sdrho0(3) sdu(4) sdg(5) 100 1 52E-3 1.11E-3 ~ 3 08E-3 8 323175080969736E4 200 6 736801383463994E4 5 259081611732966E4 1 41E-3 4.119972844994704E4 400 2 8170552508135935E4 2 741770745680064E4 6 305408088619393E4 2 2536437449718904E-4 800 1 186603519756628E4 1 5295908656314994E4 3.053930398538809E4 1 3909234596653616E-4 1600 6 032028873662157E-5 9 239357015461792E-6 1 5687832872800372E4 9 293153940622746E-5 10000 1 1658732942278758E-5 2 25516189404611E-5 2.863559G239531658E-5 2 634178622789456E-5 100 200 400 300 1600 10000 1 00E+5 1 00E+5 1.00E+5 1.00E+5 1.52E+002 O.OOE+OOO 1 52E+002 O.OOE+OOO 1.52E+002 0.00E+00Q 1 11E+002 O.OOE+OOO 1 11E+002 0 Q0E+000 1 11E+002 0 Q0E+000 3.08E+002 141E-Q03 3 08E+002 141E-003 3.08E+002 141E-003 1 OOE+O05 O.OOE+OOO 1.00E+Q05 0 OOE+OOO 1 OOE+005 0 00E+000 Figure 5.9: The standard deviations obtained for resolutions of 100, 200, 400, 800, 1600 and 1 0 0 0 0 cells. 1.2QE+005 = -13.24 ln{x)* 141.74 RJ = 0.13 1.0CE-005 A A A .......................................... tx'i = -11976.03 Irux) 128173.62 FF= 0.13 Standard Deviation &00E*004 tx}=-13.24ln^)v 19525 FP = 0.13 600E*004 \ tx) = *36.38Irtfx) + 39463 R*= 0.13 4.00E»0(H 200E+OW O.OOE-OOO » • 0 ■------ ■ 2000 4000 6000 £800 10000 Resdutian i n Number of Cel Is ■ s&.Q) v LogarithmicRegression fersA C ) sAXA) sdrho0(3) LogarithmcRegression for sdui+) * sdpt5) Icgarithmc Regression for sdrho0( 3} Logarithrric Regression tbrsdprS) Figure 5.10: The order of error in the method. 12000 Sec. 5.1 Sod Shock Tube Tests and Results 58 This simulation was also performed in the other direction, that is, with a net flow from right to left. The purpose of this was to test the Riemann solver to ensure it performed the same no matter the direction of velocity. It was found that similar standard deviations were obtained. The results are shown in Figure 5 .11- Figure 5.14. Figure 5.15 shows the percent differences between the Sod data in the left to right and right to left directions. They range from 9.5 x 10“ 3% to 5.9 x 10_3%, which shows very close agreement. Sec. 5.1 Sod Shock Tube Tests and Results 59 WHHIHtlttHHHtU Approximate Solution + Exact Solution Difference — *■ The Fluid Rest Density (rhoO) 0.8 0.6 itnviiiiiiiiiitiiiiiiiiitiiiiiiiiiiiiiiiiii 0.4 0 .2 -0.3 iiiiiiiiniiitiifliiiiKi -0 .2 - 0.1 0 0.1 0 .2 Length of Shock Tube (km) Figure 5.11: The plot of fluid density when applying the Sod data with a resolution of 400 cells in the right to left direction. 0.3 Sec. 5.1 Sod Shock Tube Tests and Results The Fluid Pressure Approximate Solution + Exact Solution Difference — #■ Length of Shock Tube (km) Figure 5.12: The plot of fluid pressure when applying the Sod data with a resolution of 400 cells in the right to left direction. 60 Sec. 5.1 Sod Shock Tube Tests and Results 61 | _i The Fluid Internal Energy Approximate Solution + Exact Solution Difference — * - IIIMIIHIIHHIlllll iifiiiiififini# ^MttltlltUtlllllllltflltlftllllllillllllltllll 0^ -0.3 -0.2 -0.1 0 0.1 0 .2 Length of Shock Tube (km) Figure 5.13: The plot of fluid internal energy when applying the Sod data with a resolution of 400 cells in the right to left direction. 0.3 Sec. 5.1 Sod Shock Tube Tests and Results 62 0.4 Approximate Solution Exact Solution Difference 0.3 The Fluid Velocity 0.2 0.1 0 - 0.1 - 0.2 -0.3 -0.4 -0.5 -0.3 - 0.2 - 0.1 0.1 0 0.2 Length of Shock Tube (km ) Figure 5.14: The plot o f fluid velocity when applying the Sod data with a resolution of 400 cells in the right to left direction. Resolution(l) s*/(2) 400 2.817E-4 400 1.867E4 Difference 9.50E-5 sdrho0(3) 2.741E-4 3.331E-4 sdu<4) 6 305E-4 6.906E-4 sd£(5) 2 253E-4 3 030E-4 5.9QE-5 6 01E-5 7 77E-5 % Difference 9.50E-003 5.90E-003 6.01E-003 7 77E-003 Figure 5.15: Percent Difference between Sod data in both directions. 0.3 Sec. 5.1 Sod Shock Tube Tests and Results 63 Figures 5.16 - 5.19 show that the Riemann solver begins to breakdown as the velocity approaches the speed of light, v = 1. This data shows results for v = 0.65 c. 1.2 Approximate Solution Exact Solution Difference iiiiiiiiiiiiiiiiii:iiiiiiiii)i!iiiiiii:::iiiiii{i:ii! iiimiiiiiimiiiHHiiiiiiiitiiiiiiiiHiiiiiiiiui o o .ck_ 0.8 £• w IIrt > "5 _2 0.4 LL Approximate Solution + Exact Solution Difference — «- 0.2 0.1 0 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 Length of Shock Tube (km) Figure 5.19: The increase in the difference in velocity for v = 0.65. Observation of Figure 5.4 and Figure 5.19 reveals that the standard deviation increases from a; 0.34 to w 0.65 as the velocity increases. This is an indication that the solution becomes unstable as the velocity approaches the speed of light. Figure 5.3 shows a clear breakdown in the accuracy of the approximate solution as opposed to Figure 5.18. A value of v = 0.65 was chosen as the limit to which the results obtained from this thesis code can be trusted. Sec. 5.2 5.2 Simulation Results using KKT data and Neutrino Fluxes Zeroed 67 Simulation Results using KKT data and Neutrino Fluxes Zeroed Now that confidence in the exact Riemann solver had been established, the code was switched to perform simulations using the Kuroda et al. data with the neutrino fluxes zeroed. This would set the benchmark results against which the later simulations with neutrino fluxes set to particular values would be tested. The initial data which is referred to as the KKT data are shown in Figures 5.20 5.23. The results at 400 and 1600 cells, for a runtime of 750000 km (2.5 s) are shown through Figures 5.24 - 5.31. Sec. 5.2 Simulation Results using KKT data and Neutrino Fluxes Zeroed 68 The Initial Fluid Rest Density. Rest Density in kg/kmA3 Rest Mass 1.4e+036 1.2e+036 le + 0 3 6 8e+035 6e+035 2e+035 0 100 200 300 400 500 600 Length of Shock Tube (km) Figure 5.20: Initial rest density data. 700 800 900 1000 Sec. 5.2 Simulation Results using KKT data and Neutrino Fluxes Zeroed 69 The Initial Fluid Pressure. Pressure Pressure in Pascals le+043 8e+042 6e+042 4e+042 2e+042 0 100 200 300 400 500 600 Length of Shock Tube (km) Figure 5.21: Initial pressure data. 700 800 900 1000 Sec. 5.2 Simulation Results using KKT data and Neutrino Fluxes Zeroed 70 The Initial Fluid Internal Energy. 8e+029 Internal Energy Internal Energy In Joules 7e+029 6e+029 4e+029 3e+029 2e+029 le + 0 2 9 0 100 200 300 400 500 600 Length of Shock Tube (km) Figure 5.22: Initial internal energy data. 700 800 900 1000 Sec. 5.2 Simulation Results using KKT data and Neutrino Fluxes Zeroed 71 The Initial Fluid Velocity. 100000 Velocity + 80000 60000 Velocity in km/s 40000 20000 -20000 -40000 -60000 -80000 100000 0 100 200 300 400 500 600 Length of Shock Tube (km) Figure 5.23: Initial velocity data. 700 800 900 1000 Sec. 5.2 Simulation Results using KKT data and Neutrino Fluxes Zeroed 72 Approximate Fluid Rest Density. 3e+036 Rest Density (kg/km Rest Density + m le+036 0 100 200 300 400 500 600 700 800 900 Length of Shock Tube (km) Figure 5.24: The results under Kuroda data and the neutrino fluxes zeroed, for the rest density and at 400 cells. 1000 Sec. 5.2 Simulation Results using KKT data and Neutrino Fluxes Zeroed 73 Approximate Fluid Pressure. 4.5e+043 Approximate Fluid Pressure 4e+043 Pressure in Pascals 3.5e+043 3e+043 2.5e+043 2e+043 + 1.5e+043 le + 0 4 3 5e+042 0 0 100 i i 200 300 400 500 600 700 800 900 Length of Shock Tube (km ) Figure 5.25: The results under Kuroda data and the neutrino fluxes zeroed, for the pressure and at 400 cells. 1000 Sec. 5.2 Simulation Results using KKT data and Neutrino Fluxes Zeroed 74 Approximate Fluid Internal Energy. 3.5e+030 ■........ "i— ... i' Internal Enerqy r ... r + Internal Energy In Joules 3e+030 2.5e+030 + 2e+030 ★ W* t 1.5e+030 *** le+030 - + 5e+029 0 0 100 200 300 400 500 600 700 800 900 Length of Shock Tube (km) F igure 5.26: The results under Kuroda data and the neutrino fluxes zeroed, for the internal energy and at 400 cells. 1000 Sec. 5.2 Simulation Results using KKT data and Neutrino Fluxes Zeroed 75 Approximate Fluid Velocity. 60000 Velocity 50000 Velocity (km/s) 40000 30000 20000 10000 -10000 -20000 -30000 0 100 200 300 400 500 600 700 800 900 Length of Shock Tube (km) Figure 5.27: The results under Kuroda data and the neutrino fluxes zeroed, for the velocity and at 400 cells. 1000 Sec. 5.2 Simulation Results using KKT data and Neutrino Fluxes Zeroed 76 Approximate Fluid Rest Density. 3e+036 Rest Density (kg/km A3) Rest Density 2e+036 le+036 5e+035 - 0 100 200 300 400 500 600 700 800 900 Length of Shock Tube (km) Figure 5.28: The results under Kuroda data and the neutrino fluxes zeroed, for the rest density and at 1600 cells. 1000 Sec. 5.2 Simulation Results using KKT data and Neutrino Fluxes Zeroed 77 Approximate Fluid Pressure. 4.5e+043 Approximate Fluid Pressure + 4e+043 Pressure in Pascals 3.5e+Q43 2.5e+043 + 1.5e+043 le+043 5e+042 i 0 100 200 i 300 400 500 600 700 800 900 Length of Shock Tube (km) Figure 5.29: The results under Kuroda data and the neutrino fluxes zeroed, for the pressure and at 1600 cells. 1000 Sec. 5.2 Simulation Results using KKT data and Neutrino Fluxes Zeroed 78 Approximate Fluid Internal Energy. 3.5e+030 Internal Energy + Internal Energy In Joules 3e+030 2e+030 1.5e+030 le+030 0 100 200 300 400 500 600 700 800 900 Length of Shock Tube (km) Figure 5.30: The results under Kuroda data and the neutrino fluxes zeroed, for the internal energy and at 1600 cells. 1000 Sec. 5.2 Simulation Results using KKT data and Neutrino Fluxes Zeroed 79 Approximate Fluid Velocity. 60000 Velocity 50000 40000 30000 20000 10000 -10000 -20000 -30000 0 100 200 300 400 500 600 700 800 900 Length of Shock Tube (km) Figure 5.31: The results under Kuroda data and the neutrino fluxes zeroed, for the velocity and at 1600 cells. 1000 Sec. 5.2 Simulation Results using KKT data and Neutrino Fluxes Zeroed 80 A comparison between the results in Figures 5.24 - 5.31 is shown in Figures 5.32 - 5.35. In the overlay of the rest density, Figure 5.32, it can be observed that there is an oscillation at the left end of the solution, that is, at the front of the outgoing shock. This is a physical instability due to the interaction of fluid going out and incoming fluid; it is not numerical as the magnitude of the spike does not change with increasing resolution. Looking at Figure 5.32 for the density overlay and Figure 5.34 for the internal energy, it can be seen that there is an instance of the waveform which is in the middle. This is a numerical instability because its magnitude of oscillation decreases with increasing resolution. This interesting effect is clearly visible between 500 - 600 km, where there is a spike in the rest density (rest mass). Since the net flow is outward in that region, with no such spike in the velocity (the velocity is actually a perfect ramp with positive slope), this can be viewed as the manifestation of mass being transported outward. This spike at the top of the density profile indicates that the net shock (result of two shocks moving in on each other) moves outward, thus carrying mass outward. This conclusion is confirmed in Figure 5.34, where it is observed that the internal energy of the fluid increases outward, with a positive energy behind that mass spike observed in Figure 5.32. The pressure in Figure 5.33 also confirms this conclusion as it shows a ramp with a high pressure behind the mass, which therefore pushes material outward to the “surface” of the shock tube. The velocity ramp is of special importance. It shows a very well-defined ramp between the two shocks with a positive slope, attaining a maximum velocity of 6.0 x 104 km s-1 , or 0.2c. This is roughly 1 2 times the velocity outside the outer shock and 2 times the velocity behind the inner shock. The net effect is a strong push outward. The strong linearity of the velocity suggests that an exact solution to this problem may exist. Sec. 5.2 Simulation Results using KKT data and Neutrino Fluxes Zeroed 81 Overlay of Density for resolutions of 400 and 1600 cells. 3e+036 t------------- 1------------- 1 r ' I I Resolution = 400 Resolution = 1600 f + 2.5e+036 2e+036 1.5e+036 le + 0 3 6 5e+035 + 0 0 100 200 300 400 500 600 700 800 Length of Shock Tube (km ) Figure 5.32: Overlay o f the Kuroda results for rest density at 400 and 1600 cells. 900 1000 Sec. 5.2 Simulation Results using KKT data and Neutrino Fluxes Zeroed 82 Overlay of Pressure for resolutions of 400 and 1600 cells. Resolution = 400 Resolution = 1600 4e+043 Fluid Pressure in Pascals 3.5e+043 2e+043 le + 0 4 3 5e+042 *#• 0 100 200 300 400 500 600 700 800 Length of Shock Tube (km) F igure 5.33: Overlay of the Kuroda results for pressure at 400 and 1600 cells. 900 1000 Sec. 5.2 Simulation Results using KKT data and Neutrino Fluxes Zeroed 83 Overlay of Internal Energy for resolutions of 400 and 1600 cells. 3.5e+030 i Resolution = 400 Resolution = 1600 Internal Energy In Joules 3e+030 — r i i i i i r 700 800 900 + ft v 2.5e+030 2e+030 & 1.5e+030 le+030 5e+029 I___________ I___________ I___________ I___________ I___________ L 0 0 100 200 300 400 500 600 Length of Shock Tube (km) Figure 5.34: Overlay of the Kuroda results for internal energy at 400 and 1600 cells. 1000 Sec. 5.2 Simulation Results using KKT data and Neutrino Fluxes Zeroed 84 Overlay of Velocity for resolutions of 400 and 1600 cells. 60000 Resolution = 400 Resolution = 1600 50000 Velocity (km/s) 40000 30000 20000 10000 -10000 -20000 -30000 0 100 200 300 400 500 600 Length of Shock Tube (km) 700 800 ' Figure 5.35: Overlay of the Kuroda results for velocity at 400 and 1600 cells. 900 1000 Sec. 5.3 5.3 Simulation Results with Neutrino Fluxes Applied 85 Simulation Results with Neutrino Fluxes Applied The thesis code was then run with the Kuroda et al. data and neutrino fluxes applied. Since there are two fluxes supplied by the neutrino radiation, one being the energy flux and the other being the momentum flux, then three test cases were used: only the neutrino energy flux, only the neutrino momentum flux and both fluxes applied. The scenarios are discussed in this order. 5.3.1 Only the Neutrino Energy Flux Applied The thesis code was run using the Kuroda et al. data and the neutrinos switched on, with the specification that there is only a neutrino energy flux present. The value of this neutrino energy flux was given as 1.0 x 10~2, a factor multiplying the value of 2.92 x 10~n km - 1 in geometrized units (referring to Appendix A Section A. 6 ). The results are shown in Figures 5.36 - 5.43. In Figure 5.36 it can be observed that there is a large amount of mass pushed against the mass to the right of « 495 km. It is seen that the mass to the right of 500 km is reacting to the large influx of mass provided by the neutrino energy flux in the time this data was generated (750000 km or 2.5 s). This reaction is indicated by the bump just before 700 km. The bump at around 565 km is the net push of mass from left to right. The fluid’s internal energy (see Figure 5.38) is curious, as it shows a rapid decrease in energy for the right moving shock but a huge amount of energy in front of the left moving shock. The conclusion to be garnered from this graph is that the shock system is driving the outer shock in the outward direction. Figures 5.40 - 5.43 show the time overlays of each fluid variable, for times of 2.5 s and 3.33 s. It can clearly be seen from these overlays that there is a net movement of mass outward (to the right as seen in Figure 5.40). There is a movement of mass inward to the neutrino-sphere, but this is expected as the cooling effects, which are ignored in this thesis but physically are still present, cause some matter to fall onto the nascent neutron star. However, what is being observed here is that there is a net explosion. The velocity overlay, Figure 5.43, shows this clearly, as there is a powerful shockwave moving outward. Sec. 5.3 Simulation Results with Neutrino Fluxes Applied 86 Fluid Rest Density for Neutrino Energy Flux (kgAm A3) Overlay of Rest Density for resolutions of 400 and 1600 cells. 6e+036 Resolution = 400 Resolution = 1600 5e+036 4e+036 3e+036 2e+036 le+036 0 x x 100 200 300 400 500 600 700 800 900 Length of Shock Tube (km) Figure 5.36: The results under Kuroda data and only a neutrino energy flux of 1.0e-2 geometrized units, for the rest density at 400 and 1600 cells. 1000 Sec. 5.3 Simulation Results with Neutrino Fluxes Applied 87 Overlay of Pressure for resolutions of 400 and 1600 cells. Fluid Pressure for Neutrino Energy Flux (Pa) 4e+043 Resolution = 400 Resolution = 1600 3.5e+043 + 3e+Q43 2.5e+043 +4 2e+043 + 1.5e+043 le+043 5e+042 0 0 100 200 300 400 500 600 700 800 900 Length of Shock Tube (km) Figure 5.37: The results under Kuroda data and only a neutrino energy flux of 1,0e-2 geometrized units, for the pressure at 400 and 1600 cells. 1000 Sec. 5.3 Simulation Results with Neutrino Fluxes Applied 88 Fluid Internal Energy for Neutrino Energy Flux (J) Overlay of Internal Energy for resolutions of 400 and 1600 cells. Resolution = 400 Resolution = 1600 % + 2e+030 1.5e+030 le+030 5e+029 t#> 0 100 200 300 400 500 600 700 800 900 Length of Shock Tube (km) Figure 5.38: The results under Kuroda data and only a neutrino energy flux of 1,0e-2 geometrized units, for the internal energy at 400 and 1600 cells. 1000 Sec. 5.3 Simulation Results with Neutrino Fluxes Applied 89 Overlay of Velocity for resolutions of 400 and 1600 cells. Fluid Velocity for Neutrino Energy Flux (km/s) 40000 Resolution = 400 Resolution = 1600 30000 20000 10000 -H I -10000 -20000 -30000 0 100 200 300 400 500 600 700 800 900 Length of Shock Tube (km) Figure 5.39: The results under Kuroda data and only a neutrino energy flux of 1.Oe-2 geometrized units, for the velocity at 400 and 1600 cells. 1000 Sec. 5.3 Simulation Results with Neutrino Fluxes Applied 90 Fluid Rest Density for Neutrino Energy Flux (kg/km•"'3) Overlay of Rest Density for times of 7.5e5 km and 1.0e6 km @ 1600 cells. 7e+036 Time = 7.5e5 km Time = 1.0e6 km 4e+036 3e+036 2e+036 le+036 i 0 100 200 l 300 400 500 600 700 800 900 Length of Shock Tube (km) Figure 5.40: Rest density data with neutrino energy flux on at times of 2.5s and 3.33s. 1000 Sec. 5.3 Simulation Results with Neutrino Fluxes Applied 91 Fluid Pressure for Neutrino Energy Flux (Pa) Overlay of Pressure for times of 7.5e5 km and 1.0e6 km @ 1600 cells. Time = 7.5e5 km Time = 1.0e6 km 2.5e+043 2e+043 1.5e+043 le+043 5e+042 0 100 200 300 400 500 600 700 800 Length of Shock Tube (km) Figure 5.41: Pressure data with neutrino energy flux on at times of 2.5s and 3.33s. 900 1000 Sec. 5.3 Simulation Results with Neutrino Fluxes Applied 92 Fluid Internal Energy for Nautrino Energy Flux (J) Overlay of Internal Energy for times of 7,5e5 km and 1.0e6 km @ 1600 cells. Time = 7.5e5 km Time = 1.0e6 km 2e+030 1.5e+030 le+030 5e+029 0 100 200 300 400 500 600 700 800 900 Length of Shock Tube (km) F igure 5.42: Internal energy data with neutrino energy flux on at times of 2.5s and 3.33s. 1000 Sec. 5.3 Simulation Results with Neutrino Fluxes Applied 93 Overlay of Velocity for times of 7.5e5 km and 1.0e6 km @ 1600 cells. Fluid Velocity for Neutrino Energy Flux (km/s) 50000 Time = 7.5e5 km Time = 1.0e6 km 40000 30000 20000 10000 -10000 -20000 -30000 0 100 200 300 400 500 600 700 800 Length of Shock T ube (km) F igure 5.43: Velocity data with neutrino energy flux on at times of 2.5s and 3.33s. 900 1000 Sec. 5.3 5.3.2 Simulation Results with Neutrino Fluxes Applied 94 Only the Neutrino Momentum Flux Applied The thesis code was also run with the neutrino energy flux set to 0.0 and the neutrino momentum flux set to 1.0 x 10-2 which multiplies a factor of 2.19 x 10~8 km - 3 in geometrized units. The results are shown in Figures 5.44 - 5.51. In Figure 5.44 it can be observed that there is a much greater effect as opposed to Figure 5.36, where the neutrino energy injected into the fluid has a less dramatic effect on the density. The jum p in fluid energy is consistent with the assumption that the type of neutrino under consideration, that is, electron-neutrinos, is a good assumption, since their cross-sections are the largest amongst the three species of neutrinos. Also, their population being the highest lends credence to the large effect observed in this thesis. The conclusion is that the neutrino energy flux has to be much greater than the neutrino momentum flux in order to have a similar effect on the fluid density. It is most likely due to a discrepancy in the choice of neutrino energy flux relative to the choice for the neutrino momentum flux. This indicates a limitation of this thesis’ model. This conclusion is further bolstered by the graph of the fluid’s internal energy (see Figure 5.46). Here it can be seen that there is one energy spike which corresponds to the one mass spike observed at « 600 km. There is a small spike at the inner shock, but this is greater than the spike in the case where the Kuroda et al. data was used with the neutrino fluxes zeroed (see Figure 5.34). It is reasonable to suggest that the neutrino flux has a dispersion effect on the fluid properties of density and energy. This is reflected in the graphs o f pressure and velocity (see Figures 5.45 - 5.47) Figures 5.48 - 5.51 show the time overlays for all the fluid variables. Times of 2.5 s and 3.33 s are shown. In the rest density graph, Figure 5.48, it can bee seen that the neutrino momentum flux drives matter inward and outward, with a large rarefaction between these. This leads to the conclusion that the neutrino momentum flux involves a large amount of interaction between the neutrinos and fluid. This effect is expected as this thesis chose to use electron-neutrinos due to their large cross-section and large numbers. Fluid Rest Density for Neutrino Momentum Flux (kg/km A3) Sec. 5.3 Simulation Results with Neutrino Fluxes Applied 95 Overlay of Rest Density for resolutions of 400 and 1600 cells. 3e+036 T-------------- 1-------------- 1-------------- 1---------------1-------------- 1.............. I......... ...... I I------- Resolution = 400 Resolution = 1600 + 2.5e+036 2e+036 1.5e+036 le+036 t + 5e+035 0 0 100 200 300 400 500 600 700 800 900 Length of Shock Tube (km) Figure 5.44: The results under Kuroda data and only a neutrino flux of 1.0e-2 geometrized units, for the rest density at 400 and 1600 cells. 1000 Sec. 5.3 Simulation Results with Neutrino Fluxes Applied 96 Fluid Pressure for Neutrino Momentum Flux (Pa) Overlay of Pressure for resolutions of 400 and 1600 cells. Resolution = 400 Resolution = 1600 4e+043 2.5e+043 1.5e+043 le+043 5e+042 UK 0 100 200 300 400 500 600 700 800 900 Length of Shock Tube (km) Figure 5.45: The results under Kuroda data and only a neutrino flux of 1.0e-2 geometrized units, for the pressure at 400 and 1600 cells. 1000 Sec. 5.3 Simulation Results with Neutrino Fluxes Applied 97 Fluid Interna I Energy for Neutrino Momentum Flux (J) Overlay of Internal Energy for resolutions of 400 and 1600 cells. 4e+030 3.5e+030 . 1 1 Resolution = 400 Resolution = 1600 1 " I ....... I ' I I I I 800 900 + 3e+030 , 2.5e+030 I f 2e+030 + y j : I 1.5e+030 le+030 + 5e+029 0 0 100 200 300 400 500 600 700 Length of Shock Tube (km) Figure 5.46: The results under Kuroda data and only a neutrino flux of 1,0e-2 geometrized units, for the internal energy at 400 and 1600 cells. 1000 Sec. 5.3 Simulation Results with Neutrino Fluxes Applied 98 Flux (km/'s) Overlay of Velocity for resolutions of 400 and 1600 cells. 50000 40000 _ Resolution = 400 Resolution = 1600 Fluid Velocity for Neutrino Momentum 30000 20000 10000 -10000 -20000 -30000 0 100 200 300 400 500 600 700 800 900 Length of Shock Tube (km) Figure 5.47: The results under Kuroda data and only a neutrino flux of 1.0e-2 geometrized units, for the velocity at 400 and 1600 cells. 1000 Sec. 5.3 Simulation Results with Neutrino Fluxes Applied 99 Fluid Rest Density for Neutrino Momentum Flux (kg/km Overlay of Rest Density for times of 7.5e5 km and 1.0e6 km @ 1600 cells. 3e+036 Time = 7.5e5 km Time = 1.0e6 km + 2.5e+036 2e+036 1.5e+036 le+036 0 100 200 300 400 500 600 700 800 900 Length of Shock Tube (km) F igure 5.48: Rest density data with neutrino momentum flux on at times of 2.5s and 3.33s. 1000 Sec. 5.3 Simulation Results with Neutrino Fluxes Applied 100 Fluid Pressure for Neutrino Momentum Flux (Pa) Overlay of Pressure for times of 7.5e5 km and 1.0e6 km @ 1600 cells. 5e+043 Time = 7.5e5 km Time = 1.0e6 km 4e+043 3.5e+043 3e+043 2.5e+043 2e+043 le+043 5e+042 0 100 200 300 400 500 600 700 800 900 Length of Shock Tube (km) Figure 5.49: Pressure data with neutrino momentum flux on at times o f 2.5s and 3.33s. 1000 Sec. 5.3 Simulation Results with Neutrino Fluxes Applied 101 Fluid Internal Energy for Neutrino Momentum Flux (J) Overlay of Internal Energy for times of 7.5e5 km and 1.0e6 km @ 1600 cells. 4e+030 Time = 7.5e5 km Time = 1.0e6 km 3e+030 2e+030 le+030 5e+029 0 100 200 300 400 500 600 700 800 900 Length of Shock Tube (km) Figure 5.50: Internal energy data with neutrino momentum flux on at times of 2.5s and 3.33s. 1000 Sec. 5.3 Simulation Results with Neutrino Fluxes Applied 102 Fluid Velocity for Neutrino Momentum Flux (km /s) Overlay of Velocity for times of 7.5e5 km and 1.0e6 km <§> 1600 cells. 60000 Time = 7.5e5 km Time = 1.0e6 km 50000 40000 30000 20000 10000 -10000 -20000 -30000 0 100 200 300 400 500 600 700 800 900 Length of Shock Tube (km) Figure 5.51: Velocity data with neutrino momentum flux on at times of 2.5s and 3.33s. 1000 Sec. 5.3 5.3.3 Simulation Results with Neutrino Fluxes Applied 103 Both Neutrino Fluxes Applied The thesis code was run using the Kuroda et al. data and the neutrinos switched on, with the specification that there is a neutrino energy flux and a neutrino momentum flux present. The value of this neutrino energy flux was given as 1 .0 x 1 0 ” 2 and the value of the neutrino momentum flux is 1.0 x 10 ” 2 in geometrized units. The results for a time of 2.5 s are shown in Figures 5.52 - 5.55. The overlays for times of 2.5 s and 3.33 s are shown in Figures 5.56 - 5.59. In Figure 5.52 it can be observed that it is the summation of the density profiles of Figures 5.36 5.44. At a distance of « 495 km there is now a clearly defined mass spike, along with the one in the region between 500 - 600 km. The spike in the density at the inner shock is now the sum of the previous spikes in Figures 5.36 - 5.44. This is to be due to the form of the source term being the sum of the neutrino energy term and neutrino momentum term. This conclusion is further bolstered by the graph of the fluid’s internal energy (see Figure 5.54). Here it can be seen that there are two energy spikes which correspond to two mass spikes observed at « 495 km and between 500 —600 km as in Figure 5.38. The spikes in the internal energy drive the inner and outer shocks, producing the spikes found there. This instability is physical in nature, as it is preserved across different resolutions and different physical conditions. Figure 5.57 is the most interesting graph. It reveals a reversal of the fluid pressure at 3.33 s. This only occurs when both the neutrino fluxes are applied. As such, this is a more realistic phenomenon and shows both the supernova implosion and explosion. Material is being deposited onto the compact object on the left side (at the inner shock) where a spike in pressure is observed. The inner shock does not move very much inward between 2.5 s and 3.33 s. The outer shock does keep accelerating outward and there is a high pressure behind it, suggesting that if the simulation were run for a longer than 3.33 s time, the outer shock would approach the right end of the shock tube. Sec. 5.3 Simulation Results with Neutrino Fluxes Applied 104 Overlay of Rest Density for resolutions of 400 and 1600 cells. Fluid Rest Density for Both Neutrino Fluxes (kg/km 4.5e+036 n------------ 1------------ 1------------r r i Resolution = 400 Resolution = 1600 4e+036 + 3.5e+036 3e+036 i 2.5e+036 2e+036 1.5e+036 le+036 5e+035 J _________ I_________ I_________ I_________ I_________ I_________ I_________ L 0 0 100 200 300 400 500 600 700 800 900 Length of Shock Tube (km) Figure 5.52: The results under Kuroda data and both a neutrino energy flux of 1,0e-2 and a neutrino flux of 1.Oe-2 in geometrized units, for the rest density at 400 and 1600 cells. 1000 Sec. 5.3 Simulation Results with Neutrino Fluxes Applied 105 Overlay of Pressure for resolutions of 400 and 1600 cells Fluid Pressure for Both Neutrino Fluxes (Pa) 4e+043 T — i---------- 1----------- 1--Resolution = 400 Resolution = 1600 3.5e+043 + 3e+043 2.5e+043 -H- 2e+043 '+ 1.5e+043 le+043 5e+042 -Hfcf 0 0 100 200 300 •kM L 400 500 600 700 800 900 Length of Shock Tube (km) Figure 5.53: The results under Kuroda data and both a neutrino energy flux o f 1,0e-2 and a neutrino flux o f 1,0e-2 in geometrized units, for the pressure at 400 and 1600 cells. 1000 Sec. 5.3 Simulation Results with Neutrino Fluxes Applied 106 Fluid Internal Energy for Both Neutrino Fluxes (J) Overlay of Internal Energy for resolutions of 400 and 1600 cells. 2.5e+030 Resolution = 400 Resolution = 1600 * + -H- + 1.5e+030 le+030 0 100 200 300 400 500 600 700 800 900 Length of Shock Tube (km) Figure 5.54: The results under Kuroda data and both a neutrino energy flux of 1,0e-2 and a neutrino flux of 1.Oe-2 in geometrized units, for the internal energy at 400 and 1600 cells. 1000 Sec. 5.3 Simulation Results with Neutrino Fluxes Applied 107 Overlay of Velocity for resolutions of 400 and 1600 cells. Fluid Velocity for Both Neutrino Fluxes (km/s) 40000 Resolution = 400 Resolution = 1600 30000 20000 10000 -10000 -20000 -30000 0 100 200 300 400 500 600 700 800 900 Length of Shock Tube (km) Figure 5.55: The results under Kuroda data and both a neutrino energy flux of 1.Oe-2 and a neutrino flux of 1,0e-2 in geometrized units, for the velocity at 400 and 1600 cells. 1000 Sec. 5.3 Simulation Results with Neutrino Fluxes Applied 108 Fluid Rest Density for Both Neutrino Fluxes (kg/km A3) Overlay of Rest Density for times of 7.5e5 km and 1.0e6 km @ 1600 cells. Time = 7.5e5 km Time = 1.0e6 km 4e+036 3e+036 2.5e+036 2e+036 le+036 5e+035 0 100 200 300 400 500 600 700 800 900 Length of Shock Tube (km) Figure 5.56: Rest density data with both neutrino fluxes on at times of 2.5s and 3.33s. 1000 Sec. 5.3 Simulation Results with Neutrino Fluxes Applied 109 Overlay of Pressure for times of 7.5e5 km and 1.0e6 km @ 1600 cells. Fluid Pressure for Both Neutrino Fluxes (Pa) 4e+043 Time = 7.5e5 km Time = 1.0e6 km 3e+043 1.5e+043 le+043 0 100 200 300 400 500 600 700 800 Length of Shock Tube (km) Figure 5.57: Pressure data with both neutrino fluxes on at times of 2.5s and 3.33s. 900 1000 Sec. 5.3 Simulation Results with Neutrino Fluxes Applied 110 Fluid Internal Energy for Both Neutrino Fluxes (J) Overlay of Internal Energy for times of 7.5e5 km and 1.0e6 km @ 1600 cells. 2.5e+030 Time = 7.5e5 km Time = 1.0e6 km + 2e+030 1.5e+030 le+030 5e+029 0 100 200 300 400 500 600 700 800 900 Length of Shock Tube (km) Figure 5.58: Internal energy data with both neutrino fluxes on at times of 2.5s and 3.33s. 1000 Sec. 5.3 Simulation Results with Neutrino Fluxes Applied 111 Overlay of Velocity for times of 7.5e5 km and 1.0e6 km 1600 cells. Fluid Velocity for Both Neutrino Fluxes (km/s) 40000 Time = 7.5e5 km Time = 1.0e6 km 30000 20000 10000 -10000 -20000 -30000 0 100 200 300 400 500 600 700 800 Length of Shock Tube (km) Figure 5.59: Velocity data with both neutrino fluxes on at times of 2.5s and 3.33s. 900 1000 Sec. 5.4 5.4 Summary 112 Summary The results obtained agree with expectations. Electron-neutrinos were considered in this thesis for two reasons: they have the largest cross-sections amongst the three neutrino species, and they are produced in the largest amounts by nuclear fusion reactions. It is expected that this is a recipe for the greatest interaction with the stellar fluid. When the neutrino model, which uses the results of Kuroda et al., was activated, and the neutrino fluxes were disengaged, results were obtained which served as the benchmark for the results obtained when the fluxes were activated. In the case where only the neutrino momentum flux was active, it was found to interact with the stellar fluid to a greater degree than the case where only the neutrino energy flux was activated. When both fluxes were used, the results showed that the effects from the neutrino momentum flux summed with the effects of neutrino energy flux, resulting in greater fluid pressure and density. This is expected because the source term in the energy-momentum tensor for the neutrinos was split as the sum of the source provided by the neutrino momentum and the source provided by the neutrino energy. The very exciting result of a mixture of implosion and explosion when both the neutrino fluxes were applied and the simulation run for 3.33 s is a major triumph for this work. It is also what is to be intuitively expected. Having both neutrino fluxes active is the realistic approach, and the observation that the results are stable and TVD is great. In this scenario, the fluid pressure reverses but the outer shock is still accelerating outward. The inner shock moves inward very little, and a pressure spike is observed. This spike is not numerical, but is physical as it shows the deposition of material onto the compact object. It is conceivable that neutrinos within the neutrino-sphere can be reheated and another burst would be produced into the shock tube, thereby bolstering the outgoing shock and powering the supernova explosion. The fluid energy and pressure in all three cases have been found to be the same order of magnitude as observed in 3-dimensional core collapse simulations. This is excellent as it lends credibility to the model employed in this thesis. It is also impressive because this thesis’ model is 1 -dimensional and special relativistic, with simplified neutrino physics. 6 Conclusions and Future Directions Neutrinos dominate the process behind core collapse supernovae. Only ab o u t« 1 % or !=s 1 0 44 J of the gravitational binding energy released in the formation process of the compact remnant end up as kinetic energy of the expanding material, and 99 % of this energy is radiated away in neutrinos. Colgate & White ([13]) were the first to suggest that neutrinos may play a crucial role for the explosion by taking up gravitational binding energy from the core and depositing it in the rest of the star. Improvements in models since then and more realistic equations of state have changed the perception of the collapse processes compared to Colgate & W hite’s pioneering work. The discovery of weak neutral currents and the corresponding importance of neutrino scattering off nucleons lead to the realization that the forming neutron star is highly opaque to neutrinos. Thus, the neutrino luminosities were too low and the energy transfer rate was not large enough to invert the infall of the surrounding gas into an explosion (see Janka ([24])). For a number of years efforts were concentrated on the prompt bounce-shock mechanism, which is the process whereby the energy given to the hydrodynamical shockwave in the moment of core bounce was thought to lead directly to the ejection of the stellar mantle and envelope (see Janka ([24])). More realistic models showed that, due to severe energy loss experienced by the 113 114 shock, the shock’s outward propagation stops well inside the iron core (see Bruenn et al. ([14])). Later work showed that neutrinos can produce an explosion on a timescale much longer than previously thought. More than 0.1s after core bounce the conditions for neutrino energy deposition have significantly improved. However, these simulations produced lower than observed explosion energy. This thesis’ simulations fall into this latter category of work. The thesis model is a stripped down version of primarily the work done by Kuroda et al. ([2]), and also using stripped down equations provided by Matteo et al. ([25]), Liu et al. ([26]) and Zhang & Dai ([27]). The physical model is that of stellar fluid undergoing neutrino heating within a 1-dimensional shock tube placed against the surface of the neutrino-sphere. The tube extends 1000 km out from this location. The coordinate system is Cartesian. The isolation of the system from its surroundings make it an excellent test tube. The results of Kuroda et al. ([2]) were used to set up initial conditions for the simulation when the neutrino model was activated and runs done for each of three cases: the fluid under heating provided by only a neutrino energy flux, the fluid under heating provided by only a neutrino momentum flux and the fluid under heating provided by both neutrino energy flux and neutrino momentum flux. Two timescales were used, 2.5 s and 3.33 s. This is well after the 0.1 s quoted by Janka ([24]). The objective was to determine what would happen to the shock system over longer timescales. It was found that the shock not only continues to move outward, but is also driven by large energy and pressure behind the shock. The density profiles show that fluid mass is pushed outward, and that the fluid velocity actually increases from one timescale (2.5 s) to the next (3.33 s). According to Bruenn et al. ([14]), shocks stall at 100 — 200 km. In this thesis, such a stall was not observed; in fact, net mass movement occurs at 500 —600 km and on the longer timescale o f 3.33 s mass has moved out to 800 km. Bruenn et al. ([14]) also mention that the work in supernova simulation makes use of fluxlimited diffusion techniques. This estimates the correct flux and limits any deviation to that flux value. However, this method fails when the transition from large to small is abrupt. This thesis uses an exact Riemann relativistic solver provided by Pons, Marti and Muller ([ 1]), which is used to calculate the exact flux at cell boundaries. This solver is only implemented when a certain threshold jum p is detected. Otherwise an upwinding method is used to approximate the flux. Having the Sec. 6.1 Future Work 115 simulation run in this way greatly reduced runtimes as opposed to running the Riemann solver all the time. The Godunov method was tested thoroughly against this solver using shocks moving left to right and right to left, and also at velocities close to the speed of light. It was found that the Godunov method performs very well using Sod data, with standard deviations of between 0.22 - 0.04 %. Thus, there is high confidence that the exact Riemann solver is working correctly and that the Godunov method is also reproducing the exact solutions to the Sod data with high confidence. The conclusion is that the code is trustworthy. It is worthwhile to note that this trim model and code has reproduced the findings of much more complicated 3-dimensional codes running much more realistic neutrino equations of state and detailed microphysics. This thesis has also confirmed that on longer timescales there is no stalling of the shock system. Instead there is an explosion which reach energies in the 1030 J range. Suggestions would be to use an exact Riemann solver in the more detailed codes in order to exactly calculate the fluxes, perform simulations on longer timescales, and concentrate on electronneutrino scattering (which is what this thesis did). 6.1 Future Work This thesis employs a number of ad-hoc terms. The model is a toy model, being in 1-dimension and not taking into account a term even though the shock tube is very long, extending from the surface of the neutrino-sphere at 8000km to a region in the “atmosphere” at 100000km. The other assumption is that outside of the neutrino-sphere, it can be assumed that the optical depth drops to 0 and remains at 0 out to the actual atmosphere of the dying star. This eliminates a e r (where t is the optical depth) term in the source term, where the source is the neutrino stream produced from a hot spot in the left side of the shock tube. That source term involves two assumptions. One is that the neutrino momentum flux produced at the hot spot is that value throughout the length of the shock tube, that is, a constant. The second is that the neutrino energy flux is also a constant along the shock tube, due to decoupling of the neutrino stream from the stellar fluid outside of the neutrino-sphere. This made that source term very easy to implement, and easy to “tune” if necessary. These assumptions were done in special relativity. Sec. 6.1 Future Work 116 The optical depth needs to be considered carefully, with a full integral of d r = npo d.r carried out to find r along the shock tube. This can then be used to determine e~r which would scale the neutrino energy and neutrino flux properly along the shock tube. Finally, the tube itself needs to be changed to a frustum, to better model the spherical symmetry in a 1 -dimensional model, which would introduce a -p- term. The neutrino momentum and energy fluxes need to be calculated at each time slice using Shibata et al.’s ([3]) evolution equations for each of these fluxes. This would accurately model the neutrino fluxes along the shock tube, instead of assuming them to be constant once produced at the neutrinosphere hot spot. The work of Kuroda et al. ([2]) implement these equations. It would be interesting to determine how closely this thesis’ results agree with evolved neutrino fluxes. One effect that is known to be resolved in this case is that of the spike at the inner shock (see, for example, Figure 5.36). This occurs because of the inclusion in the fluid evolution of a constant neutrino flux. The neutrino evolution equations will correctly control this value, and lead to more stable results. Realistic neutrino and fluid equations of state need to be used. This thesis used the ideal gas equation of state, both for the fluid and neutrino gas. The work of Janka ([24]) makes it clear that more realistic equations of state produce better neutrino evolution results. This may lead to higher explosion energies; it is not clear if this may actually be the case, but integrating the more realistic equations of state into this thesis’ code will be a good test ground. This thesis has shown that delayed start shock systems produce explosion energies comparable to those observed. The literature reveals that there is opposition to the delayed start approach. Performing more simulations using this thesis’ code at different timescales will yield results which will support the delayed start mechanism. However, more realistic neutrino evolution along the shock tube may produce less interaction with the fluid, and so counter the delayed start mechanism. It is not clear which it will be, so future work should investigate this. Since it has been shown that employing an exact Riemann solver at cell boundaries in order to get the exact fluxes across them is much better than using flux-limiters, then developing this thesis’ code using more realistic fluid and neutrino equations of state, and using the neutrino flux evolution equations, will be a good test to determine runtimes. This thesis did run the code with the exact Riemann solver employed at all times, which yielded very large runtimes. It is thus known that Sec. 6.1 Future Work 117 such a code as just proposed will be slow. The results may be worth it. With increased computing power and reduced financial cost, the computational expense may balance. A Derivations and Notes 118 Sec. A .l A .l Static Einstein Tensor Components 119 Derivation of the Einstein Tensor components for a Static spacetime Starting with the calculations for combinations of a = t and ft from t - » . o = I, rf = I: l S ( d 9/it 1 d gfh dgt l \ ( d x y + dz‘ dx'J- iA) r , _ 1 „(0g„ On,., I,1 _ 2 9 { s ^ + 8 ^ \ d x'J- ) , =* r* = (A. 6 ) o = l, 0 = r: r' r '-’ “ 29 8 = t since the off diagonals are zero: + !i‘!' ' _ I I ar-> + d x ' 0x> 1 ■ (A 7i and A are functions of r only, and the off diagonals for the metric are zeros, then, for 7 = 6, r*rt = \ H ~ 2*) ( - 2 $ 'e M ) , = (A. 12) ft = /, ft — 9: 1 ^ = 29 m / dgse dgs-y + dg0~ (A' 13) d = i since the off diagonals are zero: Since \ I*=0, = (A. 19) 0, (A.20) r % = o, (A.21) = (A.22) 0. Carrying on with the calculations for combinations of a = r and /3 from t —» 0, a = ?•, ,<3 = t: pr _ 1 7—2 7(5 { dgs-f ax* dgty\ ax'J- $ !)i-7 ^0t~; \ 0 ( } r since the off diagonals are zero: 1 rr / ^ Q rt r- = 2?rT ^ + # " # J ' Since $ and A are functions of r < A ' 2 4 ) only, and the off diagonals for the metric are zeros, then, for 7 = A 0, & I7r = 0, (A.25) 17, = 0, (A.26) IT* = 0. (A.27) Sec. A. 1 Static Einstein Tensor Components 122 For 7 = /: 1 _ 1 rr ( d g H dgrt o fl 1 n ..t ' 2 y d x 1 Dx* 0gu U — dxr = I (e~2A) d l (_e2*) . 2 17, = v ' dxr y ' i ( e '“ ) (2 4 > 'e » ). =» r [ t = $ 'e 2_2A. a = (A.28) /3 = r: 1 r<5 f , d 9Sj dgn <5 = r since the off diagonals are zero: r n1 = W 2 ' r (\ dI xT71 + td tx r2 - T dT x r1 )) • (A -3°) Since $ and A are functions of r only, and the off diagonals for the metric are zeros, then, for 7 = ^ 0, & Frt = 0 , (A.31) Frff = 0 , (A.32) I* = 0. (A.33) For 7 = r: T rr r = -cyo qr A dg,rr ^ d)(Jr dx c)xr l d = - (e 2A) 2 v df/r dxr (e2A) < rr - (e_2A'' v) (2A'e2A) , 2 1 jr. rr = a '. (A.34) Sec. A .l Static Einstein Tensor Components 123 a = r, /i = 0: T-r r"T _ dgs 7 l „ r s ( d 9s e d(/0 l \ + Tyr - o J ) • 29 ( ' 5 = r since the off diagonals are zero: rr _ i, / r _ 'W } [ d x - f + dx<> dxr)- fA36) (A -36> Since <1>and A are functions of r only, and the off diagonals for the metric are zeros, then, for 7 = l',r, 4>, n;, = 0, (A.37) T 0rr = 0, (A.38) - 0. (A.39) For 7 = 8: _ 1 rr / dgr9 1 00 — o 9 2' \d xd t Oyr0 dx° Ogoe dxr = ^ (e“ 2A) f - ^ r 2 =}► = —r e ” 2A. (A.40) Likewise: a = ?’, 6 = 0 : r r - 1 ar,)( i ^ (JSl * \ _ ^ rr f dgr£ 29 \ d x6 dgr± _ dx* dxr 1 =» TL = - r s i n 2 0 e -2A. (A.46) Carrying on with the calculations for combinations of a = 6 and (3 from t —> 4>, « = 0 , 0 = /,: ^0 _ I „es ( d 9*i , d 9ti7 ^ - 5 *2 " ^ + f l ~ dflA (A '47) 6 = 0 since the off diagonals are zero: Since $ and A are functions of r only, and the off diagonals for the metric are zeros, then, for = /, r, 0, , a = 9, (3 = i': (A.49) II o I tt = 0 , l-J 7 (A.50) r1 0at — — nu) (A.51) 1 td> -~ (A.52) nu- Sec. A.l r» r '-' - 2® I an + Static Einstein Tensor Components an a.,:‘ 1 ' dxr dx9 125 (A 53 ) ( ’ <) = 0 since the off diagonals are zero: n 2 \ dx'1 Since <&and A are functions of r only, and the off diagonals for the metric are zeros, then, for 7 = i. r, d>, T% = 0, (A.55) T9rr = 0, (A.56) r "rd = 0. (A.57) For 7 = 9: r 6 = - n " ( - t L _l ^E2!L _ d 11,11\ r0 29 \ d x 9 dx? 8x°) ’ => r% = r - 1. (A.58) a = 9, 0 — 9: r 0 1 es ( ^ e dgSl dg(h ] 6 since the off diagonals are zero: 1 oe ( dgoo , dg6l dge 0, - o i r Iy 7737 + 773T - 773T ) ■ ().r '! dx9 dx6 (A.60) Since and A are functions of r only, and the off diagonals for the metric are zeros, then, for 7 = ^ 0, Sec. A .l Static Einstein Tensor Components 126 Ti = 0, (A.61) Vie = 0, (A.62) rL - o. (A.63) For 7 = r: dg0r dg0r \ . r " — 1 oo ( dgoo j ______________ er 2 \ dxr dxe dxe ) ! 52 0 - ’ ) \ d x ' ■r r Or L = ~ r - 1. (A.64) er = 0, (i = , dgoy _ ^ 9 & A ^ - 2 9 V 8Xr 8x* dx* ) ■ (A 66x and A are functions of r only, and the off diagonals for the metric are zeros, then, for 7 = C r , 9, For 7 = 0: It* = 0, (A.67) T% = 0, (A.6 8 ) T% = 0. (A.69) Sec. A. 1 Static Einstein Tensor Components ro 00 oe ( 29 \ d x * dge$ dx* 1 = 5 (r~7 ( ^ 127 dg^A 0xB) ’ rhia2ff) ’ =>T } 0 = -s in flc o s fl. (A.70) Finishing with the calculations for combinations of a = 0 and (3 from t a — . 13 = /: r* - ( a 7 i) = 0 since the off diagonals are zero: F 0 = - r / 0 ( ^ 9(t>l 4 - ^ g<:n - ^ E l l \ ty 2 \ dx7 dx1 dx* ) (A 72) Since $ and A are functions of r only, and theoff diagonals for the metric are zeros, then, for 7 = r, 0, 0, Tft = 0, (A.73) rfr = 0 , (A.74) Tie = 0, (A.75) r f , = 0. (A.76) = (]). fj = r: r f 7 = \g*6 + 2^ \ 0 .x7 ' <9.rr Ax'5 . (A.77) 5 — 0 since the off diagonals are zero: ■'0 _ ^ „4>Q ( dg T^ - y ° [ l ^ J- dgjr, _ dgn + 7 ^ - l * i ) ' a — 0, j3 = 6\ r 4 5 = (p since the off diagonals are zero: p0 _ * 4> (^96 dfjcn d( \dx* Static Einstein Tensor Components 129 . d _ ^9(>d> dxe dx* = I (r~ 2 sin-2 0) ( ~ r 2 sin2 9 2 v ' \ dxs T 6t) ■ /* om 29 I, & T + 1 S = 0 since the off diagonals are zero: Since $ and A are functions o f r only, and the off diagonals for the metric are zeros, then, for 7 = r< tt>t = 0, (A.91) F ii> I4>= 0. (A.92) For 7 = r: _ }^n4> , d g ^ ^ dg r* , = cot». (A.94) To keep all these results manageable, a table is used, II r J, = o r ‘r = *' r lo = 0 rL to = 0 fc- C , = *' Ko = 0 r ‘* = 0 a = f, 3 — 9 r« = ° r ‘,. = 0 r* —nU 1Or r ofl = 0 n .=0 17, = «'e2* - 2A r*r = 0 r;;. = 0 II “55. IT* = 0 •Tr = A' f> = r, 3 — 1) •« = 0 •« = 0 •’Sr = 0 •T. = 0 •Te = - f e " 2A >Tr = 0 1 *8 = 0 11 II II c II ft. II 7= r II 1 = t -> II ■e- Table A.l: The Christoffel symbols for a Static, Spherically Symmetric Spacetime. tv = /■, 8 — t II II II II e II cr^ a = 9 , ,3 = r a = 0 . 11 = 0 II rr^ II G a = . 3 = t r^ = ° 1T« = 0 •To = 0 •A r = 0 •To = 0 = - r sin2 0e~ 2A rU = 0 r?, = o r?r = o I'?» = 0 rf^ = o r1 rr 9 =0 r r8 = r ” 1 r?o = 0 r 0t = 0 rSr = r - ‘ r 88= 0 r *t = ° •1r = ° =0 r, 3 = 4> *1. = ° 1 t r -- rr" 1 il 0 a = , 3 = 9 >1 a = . 8 = r =° r fo = ° rM=° r?o = *’" 1 C = c o te r to = cot 9 r 00 f, = 0 Now the Ricci curvature elements need to be calculated. RaS = - # " l F L+ - r^ r and the previous table of results for the Christoffel symbols. Carrying on then: For a — t, 8 = t\ (A.95) Sec. A .l DT1 R" = ~ 131 BP1 “7# + = Static Einstein Tensor Components “ r"fIA (A'%) « ) - ± (I7 J + rjr*, - rjrf, = £9.xr ; (* '« “ ■'“ ) - tL ( r « + r * + P , + it* + fr!,r!, + r;,rt, + r X + rf,r‘ + r 'x + r;,r;r + r» r;r+ r?,r . pt p(? . rr pff i p(? p(? I p<£p0 ' ft / + + 1 t r 1 t ^ 1 tOl t(j>' tL tj = (24>'2 - 2 $ 'A' + $ " ) e2^ - 2A - ($ ') 'e2$~2A + ( $ ' + A' + r ” 1 + r " 1 - $ J\ ') $A'g 2^ - 2A fltt = ( $ ' 2 - $ 'A ' + ^ + $ " ) e2* - 2A “3 For a = r, fj — r: BP'1 BP1 Bx-i Dx5 f># u (A.97) ^2 I P7 p<5 a3 16 _ p7 p5 aS 01 (A.98) Sec. A .l Static Einstein Tensor Components 132 fir 7 r> Q7, p7 ^ ~ dx1 dx0 p<5 ' p 7 p'S R" = i f f (r- ) “ I ? (r^ (A.99) (>s ,iy T<* + r;'r 's ‘ r " r - = ^ (A,)^ ( r‘‘ + r" + r"«+ r- ) i ( 'p t p t i p r pt p# jpt T* r r i- ff ~r I r r i r( T I r r l Qi i i- r r Afit i r f r r _L p pr _i_ p e p r _j_ P0 p r tr ' rr rr ' rr Or ~ rr _i_ F * F ^ _i_ F r F ^ i_ p O p O ■ +■ 1 r r i t 0 ~r I r r l r 0 “r i r r l qq i pi ' p4> i p r p4> * r r L f f i "1“ L r r I r fi L j-r 4>r , ptj) pO t i- r r l fiQ , p 9 p(f> , p<£ p4> ] “T 1 r r A ( ) f i *F -L r r l f i f i J - (p.r', + r;rp, + rf8r», + rf/* i pr pt . pr pr I pr p9 rr ’ rr ' \pB pt ' r t rO , ' p# p r , rr r0 ' i p4> p t _i_ p p r ‘ rt ' rt*- r

p 0 , rr 7‘ p(f> \ rB r

'2 - A ' 2 - — => R rr = - |V 2- $ 'A ' + $ " - ^ A '^ i> 0,5 For a = 9, f3 = 9: (A. 100) oj al . p7 p5 fop' r).r J 0(3 75 p 7 ptf 01 (A .101) Sec. A. 1 Static Einstein Tensor Components p a/3 aj3 ~ Ox1 ft7 | p 7 p 5 dx@ p 7 p <5 133 (A. 102) (V<* ^ 7 Itm= xl . ( r y - XL. ( r j j + rj#r 4„, - r^ rj, = ^ ( - re"A) ^ ( r“ + r- + r«»+ r«) + (r‘„r!, + r„ r‘H+ r»„r‘, + r*0 ptr;t + l y i +r;„r;r + r»,r;r +r* a , pt p0 1 eev to i pr p0 . p0 p0 . p0 p0 1 0 0 * re ^ 1 0 0 1 00 + 1 ee1 0 0 . pt p0 . pr p0 . p0 p0 . p0 p0 1 00At0 "r 1 00Ar0 ~r 1 001 00 ^ V00A00 - (■>», + r'rr;, + r‘#r«, + r^r* 1 pr pt 1pr pr . pr p0 , pr p0 , p0 pr , p0 p0 . p0 p0 , p0 pt , p0 pr | p0 p0 , p0 t^0 1 + 1 Od 00 1 0 r V 6 + 1 001 00 ^ T-1 0 tL Or . p0 pt + 1 011 00 1 1 Or Or ' 00l Or ' 1 6 1 Or + 1 Or1 00~r i g g l gg -T I g^L (w 1 001 00J = 3? + ( $ ' + A' + r - 1 + r _1) • ( - r e _2A) - ^ - - e _2A - - e “ 2A 4- cot 2 9^ = 2rA 'e~2A - e " 2A + esc 2 9 - r $ 'e _2A - rA 'e_2A - 2e~2A + 2e“ 2A - cot' 0 Rge = 1 + (rA' - 1 - Ra3 For tv = 4>. ft = (p: r $ ')>\ e- —2A (A. 103) d T l t3 d T a7 l dx1 d x :i 4 - r7 r r ' r° 1 a *1 07 Sec. A .l ^ orla or7 n" /> a1 ~ ()X1 Static Einstein Tensor Components , A i p7 p 7 po Qx a + 7* 134 (A. 104) <‘<5 ^7 k„ = ^ (ry - ^ (rj,) + r y ^ - r y j , = (_,. r f 9 e - “ ) + A (- sin ^cos (?) ~ — (r^, + r;. + C + r >«PJ^ i' If r 1(f>4)rt .4t ' 4- r r4r>rt4 4r 4 a- r 4 ' (f>*}> Of ' (p t + r y +r y ; r + r"y;r +r y ; r , p i p0 ip r 1 1 tO "T 1 p6 , , p i p I p0 p<£ * A t |pr r0 ' p0 p# i p pO <£ 06 "*■1 M >1 r . p<£ p<£ A 4>(f> 6 4>4> (M>J — lf r 1c£tr4>4t 7^4- r*4 r*4 , pt i pr pt i p r p r , p r pW , p r p<£ * <£t 4>r ' (pr pr ' p 6 (pr PP (pr ■p { p t . p# p r , p0 p0 T-*- (M1 po “r 1 (jir1- pt) 7" i . p0 p pf . p<£ t'P \ , p p r , p 0 p0 + 1 7" 1 0 0 1 Hfyj, = sin 2 9 (rA 'e“ 2A - r4>'e~2A - e~2A+ l) Andnow the cross terms are needed ... ;7i'7 D For o = /. 8 = r: _ < *0 Q7 <9tc7 dx3 I p7 P<5 Q'^ 7,5 p7 p it Ql5 4,7 (A. 106) Sec. A .l Static Einstein Tensor Components 135 (A. 107) + (r[rrj, + pr;, + r X + r X + r x + r x +r X + F X i p£ p0 + 1 ,p r p + 1 tr r+ i p<£ p0 1t r 1 $0 . p ^ p i 1 t r 1 9+ 1 t r 1 4>4>) - (i'X + r X + r‘„r* + r ; x I pr pt Ip r p r ' tt r r * ~ tr rr , p r ptt t6 rr . p r p t rr I p0 p t , p0 p r . p(? p(9 . p0 p<£ t t 1 re r 1 t r L r e r 1 t e L r e + a r0 p 0 pt i ppr + i t t 1 r 4> + 1 t r 1 r 4> ' , p p 0 . p<^> prf> \ t 0 L r ^ 1 tJ Rtr = 0 r? — n3 aa dx~< F o r « = t., 8 = 9: (A. 108) dx3 + tf ^ ^ _ r"1, 0(5 ^ (A. 109) Sec. A .l Static Einstein Tensor Components 136 o r7 = + F A - F A (A.110) R i> = o f ( r 'y ~ i ? ( r ? >) + r ‘»r ? i “ r 'A = ^ (0 )- ^ ( r « + r '' + r “ + r - + (r!,rj, + r;,rt, + r ' X + r ^ , + i > , ; + r?#i - + r J X + rjA i r t rS , r8 r* i i _L F ' f ’ F ^ i pr p0 F0 + 1 teL te + 1 teL re + 1 te1 ee ~r L teL e T t T& i i T r T& + 1 teL t4>+ 1 te1 r(t>* F ^ F^J I te1 G4>+ 1 te1 4>) - (rj.rj, + r*rr;t + r*,r« + rj*r*, +r;,r‘r + + r y l + ry £ +rf,r‘„+ rfrr;„ + rJX + r*rt i r 'f'r t i p<£ p r , p<£ p 0 , p<£ p<£ \ ^ tt.* Q(j> ' * fr* Q ' * £0 06 ’ £<£ 06J > R(o = 0 (A .l 11) dri< I'o i = R"a = For a = t, 8 = : an _ ■ # + r “A - rA ( A ' I12) Sec. A .l nri Static Einstein Tensor Components 137 nri = ^ ( r w ) - g | (r?7) + = ^<°»“ ^ ( r« + r- + r"»+ r+ ( i y + ry ;, + ry * + r y , I p t p r I p r p r . pf? p r , p P p r ' tp tr ' tip rr ‘ ^ t(f> Or ‘ tp pr i p t p , p r tip t p ' p | pt? p<*> I p

p0 1 f0l , p0 p 0 (pH ' t(p (pH . p0 p«t> \ ~r 1 ftf,1 <£0 J R,* = 0 (A. 114) 9 rL Raff = = ^ For tv = r, (3 = t: an + r “A _ r " i p p 1 + 1 r t 1 r0 1 r t1 00 ^ 1 r t 1 00 j _ fI p1 rt P tt4.' P r rrLrtt -j'- P1 rOP1 t t 4' - P Tr<($t>>r tt 1 , p r p r , p r p0 , p r p0 p0 rQ*~ tr ' 1 r< rch*rrt*t L tr ""T* ^1 rr* r r L tr ' 1 rO 1 tr & tr rr pt $ -pt rt. 10 _i_p 0 p f ' * rt* t 0 1 p0 p r , p0 p0 ' * rr^ 10 ' * rO^ 10 1 ' , p# p<£ ' A r<£A tO p p r 1 p<^> p# , p<£ p \ rr t(f> ' rO t' r0 10J =>Rrt = 0 (A .l 17) am 0(10 — an Raa = Zd ^ - - W For « = r, 0 = 9 + ^ ~ iAAm Sec. A. 1 Static Einstein Tensor Components ,o p n .0F 7 n:^ _____ 22 _ 139 Ox1 dv3 F7 F^ _ "V R ' ° = a ? ( r ; w " 4 ? (I7 l) + ^ F 7 F'-1 l3~i (A.l 19) “ r ; < r »’ = 4 ? (r"> “ i b ( r ; ' + F ' - +r:'»+r'*) + (rX!. + r;»r;, + r^rj, + r* i*, + i x +r;„r;,+r;„r;r + r* r* i Pt p 0 + 1 rOL tO i p r pt? . p 1 r 9L ( 0 + . p r p 0 | p0 p . p /?,.«= o on* Raft = <9x7 For o- P = o^ nx c*7 + i y t , — F a7<5F*5 fti <9xd (A .121) Sec. A .l n /? r 7 nv'i o;^ ____ 22. 4- r ’1' dx1 O x _ ^ - J r p r , p# p0 ' r(p (pr , p0 p0 -h 1 7.0i ^ -h I r0 l , p r -pfj> _j_ pf? p 0 p f p ' , p0 ' 1 r

p<^ ) i r0 i 0 0 l r 0 i 0 0 I t (A. 123) (14 d Qfl For a = 6, (3 = t: &P 22. i p7 dxs p<5 _ p 7 p<5 7 " (A' 125) « ) - ± (ry + rj,r‘f - r^rf, = J ; (D - ( n . + n ,. + C + i* ) + (r',r;, + r;,r‘t + r»,r', + r *r\, i p t p r i P r p r i P® P r i l ^ r r + 1 6/t1 tr + 1 0 t l rr + 1 0 t L 6 r + 1 flt1 p0 1 e t 1 (6» + 1 e t L re + 1 e t l oe + 1 e t l e i pt p^ i pr i p® p ^ 1 _i_ p ^ p ^ i + 1 e t 1 t + 1 0 tA r«4 + 1 e t 1 e + 1 w 1 4>) - (r‘,r;, + r‘rr;, + r*„rf, + r^r* i P r pt i p f pr , p fl p t i _i p r p(? i p r p^> .p0 p0 i P0 P 0 + 1 e t L tr + 1 0 r L tr + l 60l tr + 1 0 1 tr p0 p r + 1 01L 1.0 +■ 1 Or W 1 00L 10 ^ 1 01 It) . p p# . p 0 p r . p p0 i p ^ p $ \ + i 0 t l t + 1 0 r l t + 1 001 t4 + 1 001 t 4>J >Ro, = 0 (A. 126) (9 F 7 ^ For o = 0. si = r: = a?r - (Ori apr + r« A - r« A (a -127) Sec. A. I o _ (9F7 BV'1 ____ 2 2 u- 3x2 Bx3 !i ~ «» = Static Einstein Tensor Components ^ (A. 128) ,,rt ^ n ) - i p K ) + r s .il, - W 142 , = ^ ( ’^ - ^ ( r»‘ + r»'+rS»+ r«) + (rlrj, + r;rr;, + r«rr‘, + r tC , p t p r | p r p r , p6) p r . prf> p r + 1 Or1 tr + 1 Or1 rr + 1 Or1 Or + 1 0 rV r , p t riff ' dr tff I p r pO , riff p£* Or i ff ' OrA06 i xi riff ffr 6 , p t p i p r p<;6 . p ' Or t o~' Or r r p ' Oi Oip p 0 p \ Or rpifiJ - (F,r‘, + rjrr;, + n,r* + r‘,rf, ,r r rf i r r r r _l v r v () _l r r i . pO p r ' 8t rr pt Or rr ' 00 rr ’ , p0 pO 0 (j> rr i p0 p4> - t i 0(L r 6 -+ -1 0ri r 0 -t- 1 0 0 1 r& -r i 0 ^ 1 r 0 i p4> p t i p p r i p4> p0 , p4> p \ + 1 0 t [ r$ + 1 Or1 r

+ 1 001 rJ (A. 129) => R()r = 0 _ p rt/s For a = 8, ft = ft: 3T 7 BT1 dx'< B x3 2!2 Qf3 — p"* 04 ^ (A. 130) Sec. A .l /ir7 c> Static Einstein Tensor Components /9T7 rt7 , p 7 p t9 .p0 p 0 ' 00 00 i pr p 0 | p0 p 0 ' 0 r 0 ' 0 0 0 0 , pt p 0 ] pr p 0 , p 0 p 0 ' 0 0 t ' 0 0 r 0 ’ 00 00 ' i p0 p0 0 0 0 0 _ I f r 4^r 44- r*9r (fit ‘ r##r4- r (4-’ r* &4> 4>t ( y +n.rtr + r;,r;r + p„r»r + r y t . p 0 pt i p 0 pr , p0 p 0 , p0 p0 t * m.L 00 "T" 1 Or 00 1 00 x 00 ■ +‘ 1 00 00 . p 0 pi 0 t 0 0 % 0 . p 0 pr 0 r 0 0 . p0 p0 ' 00 0 0 , p0 p0 \ ' 0 0 0 0 y (A. 132) = o r9T7 /? — ,>>j _ a/5 /9r 7 f>7 , p7 p<5 dx-y d x 3 ad lS _ p 7 P<5 °5 131 (A .133) Sec. A .l Static Einstein Tensor Components d n i _ iE ° * * = dx'i dxt + T< ^ - r " A 0 ,^ -y \ d ri<"= s ? ( r «> ~ i b ( r « } + r « ^ s ~ r « r f -> 144 (A-1 3 4) = ^ ( ° ) - | ? ( r *'+ r * + r - + r - ) p t + (r^rj, + r;,n, + r»,r‘, + r*r; r r i ^ thl ■*• rr ■*■rht ^ fir ^< AJti>tL *■4>r p t p0 . p r p0 , p0 p0 , p?> p0 1 $ tL (0 + 1 ^ t1 r 6 + 1 ^ t1 00 + 1 ^ t1 ^ pt p0 1 (?tL t | p r p(4 | p L tL r 6 '+' 1 tl 0 1 ^ t1 ^ - ( n , i 'n + r v n , + i y ? , + i y f , . p r p t i p r p r , p r pO , p r p '+' i tL tr "+■ 1 4>rl tr "r L (j>0 1 tr ^ 1 H >1 tr i p0 p t . pO p r , p0 p0 , p0 p + i 4>tL 10 ~T~ 1 ijjr1 W ~l~ 1 0 L 10L 10 i ['<> p( j_ p<£ p r . p<£ p0 . p<^ p \ + 1 0 t l t

rl t + 1 (fOl t + 1 (jxj>1 td>J Ra = 0 (A. 135) Sec. A. 1 Static Einstein Tensor Components p tia/3 — nr1 t (iji 7/> p7 p i (A. 137) (l(5 ^ (rj,) + r j A - rlsr;., = £ > (’■' "> - ^ ( r;» + r * + r « + r « ) = R * ' tW1 a 7 I p7 p i ( r 5r) - 145 + (r*.r;f + r;..r‘f + r»rr;, + r*.rj, , p f p r I p r p r , p0 p r , p 0 p r ' 0/ fr ' r r r ' 0r* 9r "T" ^ 0r 4>r . p f p0 . p r p0 , p0 p0 I p 0 p0 '^ 0r^ td r rO r r 09 ' 6 I p f p 0 i p r p . p(? p 0 . p 0 p 0 \ ' 0 1 f0 0r rrf> r* 9cj> rfrr 4>J - (r‘„r;(+ r^p, + rj.it, + r'Mr?, , p r p f i p r p r , p r p0 , p r p 0 ' 0f rr ' 0r rr ' 00 rr * ^ <50 rr . p0 p f , p0 p r i p0 p0 i p0 p 0 * •*'0/ r0 ’ 0r r0 * 00^ r() ' * 00 rO , p 0 p t , p 0 p r . p 0 p0 l p 0 p 0 \ * 0f r0 ’ 0r r0 00 r0 * 00 r0^ (A .138) /?0 r = 0 p _ 013 ;) r 7 " J _____ '22 I p 7 p i dx't dx3 0(3 75 _ p7 p i ^ (A. 139) Sec. A .l Static Einstein Tensor Components R "s = Tpr - + r -A ~ r '.»r ^ = a ? - 5 ? (r^ = ^ (cot#)- ^ + r^ 146 (A-l40) ' “ r« r^ ( I*« + r " + I * + r « ) \(p pt + (r^rj, + r^r‘rt + rj,r>, + r*,r d>t + r ^ ,r ;r + i pt + 1 r y ;r + r X pf? i p r te ^ pf? i pf? pf? + i1 , r ; r . p 0 pf? 1 r&^ 1 eL ee ^ L 4>oL e , pt p . p p "T" 1 ciidL 0 0 ^ 1 W 1 J Ip l p( i pt pr , p t p0 Or ' 00 Oi i 1 rf-t1 et + 1 4>rl ot + i eL ot + L tl 6 1 4>r 0 1 4>6 l 64> 1 /? R,p, = o a = r -S !3 II o Rrr - ($'2 _ $'A' + " - 2A') =0 at — 0 Ret = 0 Rgr = 0 /?„„ = 1 + (rA' - 1 - r')e 2A AV_> = o rt = (j) S3 & II © ('2 - 'A' + 2±1 +4>") e2 2A 0 C II Rtt = a —t j3-r 0? 0 = t S3 II o Table A.2: The Ricci tensor components for a Static, Spherically Symmetric Spacetime. A>r = 0 Ae = 0 = s'" 2 0 e~2A (rA' —r' - 1 + e2A) R ,g — 0 Finally, the Einstein tensor components can be determined, using: A =o Sec. A .l Static Einstein Tensor Components Crnf} = R q0 - \ g as R where, R = ga0Ra$- 147 (A. 142) (A. 143) Recall that the metric gives, 9tt , 2 (A. 144) (A. 145) !'h r = e'2\ (A. 146) 9i>0 = sin2 Ogee- (A. 147) So using the results in Table A.2, the fact that the cross terms of the metric are zero (for this type of star), and the above metric components, then the derivation of the Einstein tensor components proceeds as follows. R = ga0Rafi (A. 148) = gn R tt + grrR rr + ge0R ee + R = { - * - * * ) ■ ( [ $'2 _ + 2 $ ' + $ " ) e24,- 2A - (e-2A) • | V 2 - $ 'A ' + " + ( r ” 2) - ( l + ( r A ' - l - v $ /) e - 2A) + ( r “ 2 sin ~ 2 9) ■(sin 2 0e~2A (rA ' — r' + — + " )1 e - [ (j) '2 _ 148 2A — i ( - e 2*) • ( ( - e - 2$) • ( ( V 2 - 3>'A' + ^ + $ ti \ g2"2A - (e“ 2A) • ^" - - A '^ + (r~ 2) • ( l + (rA ' — 1 - r $ ') e _2A) + ( r -2 sin - 2 9) ■(sin 2 9e~2A (rA ' - r $ ' - 1 + e2A) ) ) - ( $ '2 - $ 'A ' + — + $ " ) e2*“ 2A r , 1 / V 2 _ $ 'A ' _ ^ + $ " ) e2*- 2A 2 <& + t ( l + ( r A ' - l - r*')e-“ ) ^ + ^ e 24> ( r ~ 2 sin - 2 9) (sin 2 0e~2A (rA ' — r $ ' - 1 + e2A)) 1 / 2-2 A + 4> l ea2 2 2 ^ p2-2A + - (rA ' - 1 - r $ ') —— b - (rA ' - 1 - r $ ' + e2A) — 1 2 \ r r 2 rz 1 1 e2* j_e '24>-2A _j_ A 7,-2 2A' 2 r2 «2$ 1 24>-2A + -Tr I e (A. 150) Sec. A .l '2 _ + " Static Einstein Tensor Components 149 2 2 \' i (e2A) - ((-e-2*) • ((V2- $ 'A ' + —r + V - (e~2A) • (V2- $ 'A ' + $ " - ^ 2' _ n \ „2') e~2A) + ( r “ 2 sin-2 9) • (sin2 0e -2A (rA ' - r $ ' - 1 4- e2A) ) ) fV2- $ 'A ' + $ " 1e ^ f $ ' 2 - $ 'A ' + 4- $ " 2A 2r 2 $' *2A r r 2$' | 1 2 r2 „2A A' a2A 1 (rA ' - r $ ' - 1) 2 r2 (A .151) Sec. A .l Static Einstein Tensor Components 150 ( ’00 — It00 — - gooR = 1 + (rA' - 1 — r') e 2A - i - 2 ^ ( - e - 2*) • ^ ( V 2 - $ 'A ' + ^ + e2 - (e~2A) • ^ $ '2 - $ 'A ' + $ " - -A '^j + (r~ 2) ■( l + (rA ' — 1 —r') e^ 2A) + (r~ 2 sin~2 6) ■(sin2 0e~2A (rA ' — r' — 1 + e2A) ) ) = 1 + (rA ' - 1 - r $ ') e~2A + ~2e~2A ^ 4- y e ^2A ^ - e 2 2 _ ^ A' + ~~~ + - - ( l + (rA ' - 1 - r $ ') e~2A) (rA ' —r<&' — 1 + e2A) = ~ ( l + (rA' - 1 - r $ ') e _2A) + r 2e~2A - 'A' + $"') Z Z Z Z v / + r e “ 2A$ ' - re~ 2AA' - ^ r e - 2AA' + ^ re ~ 2A$ ' + (U~2A - \ £ & z z Gee = r 2e ” 2A f $ ' 2 - $ 'A ' + $ " + — - — r r For a = 4>, 0 = 0: (A. 152) Sec. A .l ^00 ^00 2^ Static Einstein Tensor Components 151 ^ = sin2 0e " 2A (rA ' —r<&' — 1 + e2A) - $ 'A ' + — + <£\\ff >' \ e 24>—2A - ~ r 2 sin2 8 ^ ( - e “ 2$) • ^ - (e " 2A) • - $ 'A ' + $ " - ^A'^j + ( r " 2) • ( l + (rA ' — 1 —r (T>') e " 2A) + ( r " 2 sin "2 0) • (sin2 0e"2A (rA ' —r<3?' — 1 + e2A))) = sin2 0e " 2A (rA ' —r<&' — 1 + e2A) + ^ r 2 sin2 0e " 2A + - r 2 sin 2 0e"2A ( $ ' 2 - $ 'A ' - — + 2 \ r - ^ sin 2 0 e " 2A (rA ' - r $ ' - 2 — $ 'A ' + + "^ i sin 2 0 (1 + (rA ' - 1 - r $ ') e " 2A) 1 + e2A) = ^ r sin 2 0 e 2 AA' — ^ r sin 2 0 e 2A$ ' — ^ sin 2 0 e 2A + ^ sin 2 8 Z Li Zj Cd + r 2 sin 2 0e~2A ( V 2 - 3>'A' + <&") + r sin 2 0e"2A$ ' - r sin 2 0e"2AA 1 2 sin 2 8 — - r sin 2 0e 2AA' + - sin 2 0e 2A + - r sin 2 8e 2A<5' 2 2 = sin 2 8 r 2e 2A ( $ ' 2 - = ^ + <5" + ^ - A ) o II £ It = 0 G0r = 0 G0e = 0 G.pt, = Sill2 0G»« Sec. A.2 A.2 Static Einstein Field Equation Solutions 152 Derivation of the Einstein Field Equation solutions for a Static spacetime Using the it components of T and G , then (with the details shown in Appendix A.2), G ft = 8 nT it 2 A' 1 \ 2(„ 2. e2<1> n ---------e + — = 8npe r r1 ) r1 -^-e2$ - f ( 2 ?n(r)) = 87rpe2'I, rl dr ^ dm (r) 2 4n G p dr One more equation can be found usng the r r components of T and (7. Thus Grr = 8nTrr 2V ( + 1 e2A 2\ 5- = 8rrpe 24>V -f 1 - e 2A = 87r r 2pe2A 2 $>V = 8irr2pe2X + e2A - 1 87t r 2pe2A + e2A — 1 dr 2r _ r2 e 2A 8n r2pe2A + e2A ^ e - 2A 2r 4 ttr 3p + § - ^e r 2e - 2 A m( r ) + 47T7'3p r 2e -2A (A. 154) Sec. A.3 Static GR Equations of Motion (i m (r) + A u r p dr r 2 — r 2 + r 2e 2A ///(/•) + 47r rAp 153 r (r —r + re~2A) ??i(r) -b 47r?’3p = r (r - r (1 - e~2A)) A.3 G? m (r) + 47r r 3p dr r ( r — 2m{r)) (A. 155) Derivation of the General Relativistic Equations of Motion for a Static spacetime The equations of motion of the fluid are equivalent to the vanishing of the divergence of Tag. Proceeding in this line, then, V a T “* = (p + p) v aV au0 + V a [(p + p) ua}a0 + V apga0 (A. 156) and u g V 0 = (p + p) UQV aUgU0 + V Q {(p + = ~ Va [{P + P) p) ?iQ] UgU0 + UgVap ( f 0 + U ^ a P 9a0 = -unv np - unVap - {p + p) + UgV^pg”0 = - (Pu");a - (Pua)-a ~ \ip + P) M“]:a + i p f ^ ) .a = - [(A) + <0 K°];„ - 0™“ );,, - [(A) + + P) + {PP^'»s) ;a = - (A)«°);a - (eu°)-a - (?«“ );* “ {Po^l.a - (eM“ );a ~ (Pua)-a + {PgaSu&),a = - (2V=~gP0>n,a - ( 2 y / ^ r u % - (2v ^ p u “ ) a + (pgo8u g );n => 0 = - (2 v c r g;Pona) a so, 0 = ( s f ^ g p o i f ) (> (A. 157) (/"/“% ) . a = (2v/ z g(/) + e )va) a . (A. 158) Sec. A.4 Non-Static Einstein Tensor Components 154 The equations (A. 157), and (A. 158) can now be developed as follows. Recognize that v l = 0, since the spacetime is static. Then, for (A. 157), let D = y/^g p o u 1, so that, o = { V ^ g p o u 1) t + { V ^ opqu1) = + (\ 4u l = D t + (D j)'. ^ > 0 = £> + V - ( D v ) , (A. 159) but v = 0, so the right-hand side is zero. In a similar way, for (A. 158), let S = 2 y /^ g ( p +p)ii?, so that, with v = 0 , then, tf + V - ( S v ) = (/«,“% ) ;a 0 = ( p < T % ) ;Q = P9a0u?a + PU09?n = vtf*« + I t y ) + pu0 (< ,° i + K s!, sa + r a5r ) =» - p g ^ u ^ - pu0g ^ a = p g ^ T i ^ + p u g ( r “sg5a + T*sgaS) , where, A.4 = \ g a6 (ggs,y + tM /3 - 9^,s)- Derivation of the Einstein Tensor components for Time-Dependent spacetime In the proceeding,' = Jj; and ' = ■§;• For a = I : (A.160) Sec. A.4 Non-Static Einstein Tensor Components d OL drd{ f ) d 0 d 0 \2 T rlU £ ) T r ] +T -H ~ dt +r =► - e - a d rfr 2 dr d { 1 - '5 I ./ 7TT - r dt 2 - Iat- 2- dr i e, d a dr 2 d r a ( f ) 0t 2‘ 0 ~r = 0 dt 1 dr d d r dr 1 ^ d t dt + — .-V -= 0 d t dr dr 2 d r dr ■2 i i e A( $ dr dr dr, . „ dt dt d dr d r \ v dt dt dr dr dr dr dt 2 dr dr * V + I^ -e - f - V dr ) 2 \ dr ) - d^ v ( dt = a i + a i o . 2 /i d(p + — - r sin 0 . , dt 2 \ d'.T r , d2t d2t d 1 , (dr e dt 2 dr , d2t -si 'd 9 \2 2 . j / d\2 —e — + r 2 sin 0 - . dr ) \d r I a i dr0 (£ ) 2 2 . 3 / d(f)\2 + r sm 0 1 — \d r J - S - o at i/ dr dt 155 Sec. A.4 Non-Static Einstein Tensor Components d d 0 -r- _o_ dr d 3 DL DL d rO (^) dr 'd c b V + r sin 9 + r + r 2 sin2 9 —e d 1 a e dr 2 o i j /dr 1 dr d ( ^ ) 2 \d r = t „ (d C 2 + r - - e 156 dr dr = C . + a~r 2r <9 1 2 . 2 ( d d ' 2 _ 0_1 „ ( d t ' 2 = C ~d?2S + S 2 r s m S (37 d d(dr) 1 A / d r \ 2 A/dr V d d 1 A d r \2 :eA I — + — rrrre dr d (jd) 2 \ d r ) ' d r <9 ( ^ ) 2 V^r / ^ dr d (dr) 2 \dr/ d d 1 iA 'e A 2 + r sin + r d A dr 2 0 — d dr x dr dr dr dr = ( d dt d r 1 x dr dr d r dl dr 2 d r dr . d9 \ 2 . 2 n (d(f>\2 1 , J d t \ 2 + r ( — ) + r sm d ( — ) — - v e I — 2 \ dr / Xdr d t^ 2 dr d,r J = ( X' x / dr \ 2 2 P { Ih ) 2 .M Y . 2 0 ( d 4> y v' „ ■r ( ~— 1 - r sm 9 | — ) 4- —e \d r J => e d2r dr2 X' / d r V 2 \ dr J X dr dt 2 dr dr • 2 n ( d4>\ V V ( dt\ -r sm 2 9 - + —e" — \dr J 2 \dr) For « = 9: = C = t (A. 162) Sec. A.4 Non-Static Einstein Tensor Components 01 0 0 1 dr 0 ( ^ ) d dr 0 Tt ] + r T rW (£ ) _0_ ~ 09 —e 09 dr ^ V + r 2l sm - 22 9n 1( d4>' + r l2 I — dr e" d . + r 1 sin 2 9 . dj ) \ dr 0 1 2 ( d9' 0_ 1 07 09 2( dr 0 ( f )' ~2 T 1 • =>• r — 10 2 . 2 (d4>' r sin u I — 2 50 Vd r . ,020 + r — d r (I9 0r 0r 0r , , 1 0 dr 1 2 09 dr 2 0r 00 + 2 r - — :------ r sin 0 cos 0 0t 0r 0r2 ,0^0 1 + 0 9 2r ’ - jLlel 09 2 3 1 ■— - r 2 sm 2 9a — 50 2 I dr, , 0 20 0 dr 0 1 , / 00 i -r-z + — - —. 7 - r d r 2 d r dr 0 [jOj 2 V^r 7 = 0 5 ’ 157 — 00 0r 0r dr d,9 ’T t T t = ° 2 ■ n n ( (^ } r s m 0 co s 0 = 0 (A .163) Sec. A.4 Non-Static Einstein Tensor Components 8L 8L dr 8 ( f ) 80 0 d dr 8 _a_ dr 00 dr 00 + r dr d 8 1 2 ■ 2 n f d(t>\2 T-TTT-—7 sm 9 I I d r d (7 h )2 \ (lT ' + r 2 sin 2 0 1 1 8 1. r 2 sin 2 0 dr 8 ( f ) 2 dr J d(j>\ dt d r J1 ^ ' T t 8ct>2 \d r J 002 8 1 dt 002 \0r d dr 8 1 00 4+ - r sin 0t y 0 t dr 8 ( ^ ) 2 \ dd° ° l r2 zm 2 o ( d +, — — 90 2 \ 0r d d tV + r 2 ( M \ 2 + r 2 sin2 0 ( r r dr \d r J dr 158 00 \ 0 dr J - - — - r 2 sin2 9 ^ ^ 00 0r 2 0r 0r 0 001 00 00 - r sin 0 - — — = 0 0000 2 0r 0r 00 00 =» r 2 sin 2 0 ^ - ^ + 2 r sin 2 0 ^ - ^ + 2 r 2 s in 0 co s 0 - ^ 0r2 0r 0r 0r 0r 2 / i z i ^00 00 . 2 /i 00 -r s m 0 cos 0 — r sin 0 — — = 0 0r 0r 0r 0r 2 • 2 /i^2^ • 2/1 dr d(p 2 . 00 0 0 =>• r sin 0—— + r sm 0 - — — + r sm 0 cos 0— — = 0. 0r2 0r 0r 0 r dr The equations (A. 161), (A. 164) (A. 162), (A. 163), (A. 164) can now be compared to the geodesic equation and the Christoffel tensor components read off, resulting in the table below. From the Riemann tensor: ^ the following is obtained, - c a - Lft 2 pi _ v' L rt ~~ 2 pi _ ApA-i/ 1 rr 2 r** = o if = o II 0 St, II 33. n f= o r^r = 0 II II op r* = o a = o rt = r, ft = t pr _ v ' v - \ 1 tt — 2 a = r, ft = r tv = r, 8 = 9 pr _ A ri 2 pr _ n 1 et ~ u pr _ A fr 2 pr _ Y rr 2 cv = r, ft — ft if = o rV = ° r le = 0 if = o r re = o r;4 = o F0r = 0 r ee = - re_A if, = o if, = o rjr = 0 r^ = o r « = - r s i n s 9e_ i a = 9, ft = t rf, = o rfr = o r?„ = o ri = o rfr = o r1 6te* ~— uo p0 _ 1 o = 9, 8 = 9 if = o L Or ~ r fV = 0, ft = ft if = o rv = (A, ft = t if . = o if = o rt = ft. ft = r if , = o a = ft, ft = 9 if , = o a = ft, ft = ft if, = 0 rr = o pf> _ Q 1 dr ~ U p _ ^1 P II = o 'Ct II tv = /, 8 = t T0 pfl _ 1 = 0 J 1 4>r r rf = 0 T ro n 1 00 — ~ u r?« = o v1 d4,0 — —u n r *j = - s i n » cos# pf> _ Q 1 te u p If, = 0 p<^ __ 1 1 r^> ” r r» = o i f , = co t# i f , = co t# r

U = o 14„ = rjw - 1'*,. + rj.iv + r*„rv - n„r%, (A.i«) = IV , - r v , + ij,r v + r y v - rv rv , (a . i «7) = r;„,„ - r"„,„ + r j . r v - (A.i68> V !>%,,, = - ' V + r J ^ V + I m .iV - r j r v < a.m ») Considering Table A .l, it can be seen that unless ft = u or (ft, u) — (r, f) these components vanish. Also, Rr$^t = 0- So it is now possible to derive the relevant Ricci tensor components. Sec. A.4 pr ~ L tt,r pr Non-Static Einstein Tensor Components i pr pr , p r pt p r per L t.r,t "r 1 r r L tt "T 1 r t L tt pu/ . py p r yte,t ,0 p^/ t "t" 160 o t L tr pw 1 roL tt T f f i1 to 0 pr I p0 i tt ,0 t v —\ (A. 170) - K tr + pt K ( ) r + K r(j)r pt 1 rr l | pt pr pfl i pfl p r r6,r ' pfl p a 1 o r 1 rd rff rr i p0 pr j p^ .p** p0 pc^ r<^,r 1 err rt 0 p<£ r<£ rr 1 ~7 pt pa , p t pf 1 r^r I J- r^l rr I 1 £fl rr , / 2r 1 rr 1 X/ ^ ar*- r$ 1 1 \ \ - u ~ 7 + 2 ? A _ 7 + 2 Ae ( a - v ) eA- + p " + p 'A ' + 7*>AeA- " - 7 t/'2 + \ k 2ex~L' 4 4 4 // /2 i/ ^ A t/A' A' 2 A2 _ At> T + _T + 7 + 2 + T ~ 4 (A. 171) Sec. A.4 Non-Static Einstein Tensor Components Roe — Rho + Ro,-o + ~ jyGO.t v 't s jytit.9 I pr + i + Re n p r _i_ T't p * 1 Tt1 00 o p r ^ y r Y_ P r P r > f P 00 Qt ^ _i_ P r Jp P 1 00,r ” > 0 p 0 "T" 1 r r 1 00 P r P^ 00 ~~ 1 r+ R°4>i + *° *0 o — r * / - r 4/ -i- r 4 r r 4- r f j p ^ — r* y,t y ip t.fi rt 0 0 ^ ^ M r 0 0 - _o tf>t + i v . - rv ^ + n ,r v + i^ jT - i W r + r'00,0 _ *•y/Cdfh 4'- r■*■0rf}* r a00 ^p a00 r0lr- .00 = - - i / s i n 2 0e A sin 2 6e A + rA' sin 2 9e A — -A ' sin 2 0e A + .s itrM e ^ + sin 2 0 —-eos2"# —s i n ^ e ^ + ^es2! ! - s i n 2 0 e A 1 + 2 K - A') sin 2 0 /? 6t6i sin 0 (A .173) Sec. A.4 Non-Static Einstein Tensor Components 162 R ,, = R frtt + Ifrrt + R%t + r* r r + r/ i p<*> p r , p<£ x*v y ^ r . t "+■ 1 r4> rt rt T p' ,0 (A. 174) R rt To make it easy to refer to these values, a table is used. Table A.5: The Ricci tensor components for a Time-Dependent, Spherically Symmetric Spacetime. 3 = 7- 0=6 Rl r = i n tH = o O II 11 95. 0 = t R )t> = 0 Rr* = 0 R tt = « = t i f . * ? - - * X2 , Xi. “ T + T (» = V Krl = £ lirr = - V V - + 3^ T i + , r A , A2 At/1 -A 1/ + [ 7 + *1-------- r j e /?«,. = 0 <\ ~ rt,:,, = 0 Ft./i, - 0 Run e ~ A [l + J (* ' S) « 11 O 11 C R«i = 0 = A' )] + 1 R(i* ~ o = s in 2 At this stage, all the information is present to proceed with the derivation of the Einstein tensor components. Sec. A.4 Non-Static Einstein Tensor Components 163 The Ricci scalar, R, has to be evaluated as follows (where gtr = 0). R = < f l i f t + g rr R r r + g 0d R oe + g H R u + & * £ '2 A A2 XV \ —— 4 4— " -j- --—A — —— b ----- r - + -- cv 4 4 2 4 r 2 —e - 4- e + -A '2 u" u'X! A' ’a ---- — --- + —r - + -- + — 4 4 2 2 V 1 Xu — --4 4 A2 1+ -(I/-V) + 1 1 sin 2 9c r2 sin 2 9 (A. 176) + sin 2 9 Opting to leave this expression as it is, then the Einstein tensor components can now be resolved. Gtt — Rtt — 2 gttR u" u AV u' , 4---4 r 2 + 4 4* ----- _L. ------- A A2 Af> u" u 2 u'X' X' e~A ---------------- 1-----_ _|------- L — "b cx~u 2 4 ~T CNI 4" Xu 2 ~ T +T AV u'] , A A — --A2 Xu 4 -------- H------ e + — - tt 4 2 4 4 r 4 i/" 2 —e A_ ^v—A i + i 4" i 4" sin 9c r2 sin 9 1 + - ( u ' - A') + sin 2 9 i/—A Gtt = x e ‘ =► Gt, = e" - + J This simplifies to: r + e" ' 1 - e - A l + - (l/-A ') r2 L [—e_A [1 - rA'] + 1] Proceeding to the C rr term, then: + 1 (A. 177) Sec. A.4 (t, r —Rrr Non-Static Einstein Tensor Components 164 ^ 9rr R A2 Xu 2 + T ~T A u" lJ2 u 'X X ------------- y -------- h — + 2 4 4 r u" v '2 Xu' u' v" v2 u 'X X 4 4 — ------ -J- 2 -A 1 sin 2 Be' r 2 sin B ;/ + A' A 1 •e 2r A -|~ A A2 Xu A2 Xu ^A—v 2 + T ~ T r j V ~ A') + 1 1 + 2 _|_ A 1+ - V - l + -(u'-X) -e V) + sin 2 6 + 1 - i [ - e ' A[1 + ru'\ + 1] => G rr ~ — (A. 178) The Gee term becomes: Goo — Roe ~ £ BeoR + 1 = —e ^2 Xu' 2 + 4 4 i/' 2 1 . 21 2 => Goo = The G 00 term is: r —e 1 r 2 sin 6 \ r u" e A -----2 - sin Be 4 A _ A^ Az> _ - _ + _ A A2 + r X u f22 u 'X ~4T ^----4i— ^----r 1 + 2 V ^ja-A Xu aA—v 2 + T ~ T - A') + 1 1+ 2 V - X ) + sin B 1 u tt to -ir 2 2 + ...— H 4 r (A. 179) Sec. A.5 .sin2 6c •2 sin 2 0 \u" v' 2 ——— I" 2 4 AV u" i/X ' r -A ur + u'2 \ e"~A X' +7 + —e A 2 r. 4 T ~ T +~ Xu ^— 4 4 A2 A A^ _ Xu .A~ v —+ — — + 1 1 1 - sm ■ 2 9a -r~ 2 r2 sin 2 9 => G ^ = - r 2 sin 2 9 165 1 + - ( 1/ - V) + sin 2 9 r 2 sin 2 6 —e r 2 sin 2 9 Non-Static Einstein Field Equation Solutions sin 9c~ u"v 2 2 + sin 9 i/X! 4 + X'-u' + 4 r (A. 180) = sin 2 9Gee Finally, the G tr term is, G tr = R tr - ^ R (A. 181) = ^ Table A.6 : The Einstein tensor components for a Time-Dependent, Spherically Symmetric Spacetime. ft = * o = t G t t ~ c1' [ A = r 0 = 0 fi — & /"• _ A G t r ~ 7: G ,„ = 0 —0 G,„ = 0 GrO = 0 - ' “ U l-rA '1 + l]] l> = V a= e G ,t = $ G „t = 0 G rl. = [ 4 j [ - e - x [1 + A.5 G '„ = 0 r.6r = o —0 Derivation of the Einstein Field Equation solutions for a Time-Dependent spacetime Using the tt components of T and G, then, 0 G y (t Sec. A.5 Non-Static Einstein Field Equation Solutions 166 8ttT u = G u 87 TPcv = e" 8 tt pe" = - [ - e - A[1 - rV] 4 - 1] A' - - e " “ A + -r e ‘'" A+ -r ,* 87tpev = fr ( l - e “ A)l 8 ttps.u = l e'7^ -(2 m (r.O ) H ar dm 47rr2p. ^ r)r (A. 182) The metric component, yrr, can be obtained as: 9rr C — 1 2m (A. 183) The second constraint equation can be found using the r r components of T and G. Here, the reader is reminded of the fact that v = 2 $ , and as such, Sec. A.6 Geometrizations of Constants and Crucial Terms 167 Stt 1 •pi* — ('j j'f' [1 4- rv'] + l] 87rpeA = —eA 8Trpe 2' 1 + eA n A 7 = 87rpe r rz r2V + 1 - eA = 8 n r2pex + 1 2 'r = 87tr 2peA + eA (94> dr 1 87r r 2peA + eA 1 2r |e ~ A87r?;2peA+ eA 47ir3p + r _-A 2 r 2e A 7Tl(r, t) + 47T7’3p ?.2g-A m (r, i) + Atti'^p r2 —?'2 4- r 2e~~A m (r, /.) + 47r r 3p r (?■ —?• 4- ?-e “A) m (r, f ) 4- 47rr3p r (r —r (1 —e~A)) dr m (r, /) + 4irr3p r (r - 2m (r, t)) ’ (A. 184) Using this, then <3> can be found by numerical integration and it would follow that g u = e". The other components o f the metric, gee and are, by definition, g6e = r 2 and — r 2 sin 1 9. Now the general relativistic hydrodynamics evolution equations can be developed. A.6 Geometrizations of Constants and Crucial Terms A large amount of work was done with geometrizing constants and certain variables used in the thesis code. The geometrizations were carried out by conforming to the system shown in Table A.7. The internal energy per electron-neutrino, e t /, was geometrized as follows. From Kuroda Sec. A.6 Geometrizations of Constants and Crucial Terms kmKs Units Geometrization Length xl Time xc Mass c Velocity x lC Energy V G Pressure X-T C4 Density v G Temperature K 168 x “T cl Table A.7: A listing of al the Geometrized Quantities. et al.’s data, ev = 15 MeV, where 1 MeV = 1.602 x 10“ 13 kgm V "2. This yields ev 2.403 x 10~18 k g k m V 2. £„ = 2.403 x 10“ 18 x ~ km2. cl c = 1.78 x 10"61 x 1.0 x 109 km2, =►£■„ = 1.78 x 10-52 km -2 . (A. 185) (A.186) (A. 187) The temperature at the surface of the neutrino-sphere, and therefore of the neutrinos, was found by, T = (A. 188) kB 2.403 x 10~18 k g k m V 2 1.381 x 10~29 k g k m V 2K - r (A. 189) 1.74 x 1011 K. (A. 190) It was necessary to geometrize some constants. For the Stefan-Boltzmann constant, rr = T 5.6704 x 10“ 8 kgs 3K '4, this works out to a&eo = 1.26 x 10"52 km~3K “ 4. The Radiation Constant, a = 7.5657 x 10~16 Jm -3K -4 became ageo = 1.68 x 10-55 km _1K ~ l. Sec. A.7 Development Environment 169 The Boltzmann constant, kB = 1.381 x 10 23 JK~4 resulted in k BgC0 = 3.07 x 10 72 km2K " 1. The value of the neutrino luminosity as obtained from Kuroda et al.’s paper is = 0.8 x 1053 erg/s = 8 x 1039 kgkm2s-3 . This yields L v%gt0 = 1.78 x 104 km -2 . The atomic mass unit, m u = 1.66 x 10~27 kg, became m t).ge0 = 3.69 x 10~51 km. The electron’s mass is m ec2 = 0.511 MeV = 8.19 x 1CT20 kgkm2s-2 . This geometrizes to (m cc2)gc0 = 1.82 x 10~53 km, and (m ec2)2eo = 3.31 x 10-106 km2. Now k , the opacity of the neutrino fluid, can be calculated. From Janka, ([24]), this is given by, 24 \m eczy rnu Also by Janka are the values: cv = -1 .2 6 , Yn + Yp « 1 and the electron-neutrino cross-section, (T() = 1.76 x 10~°4 km2. The other values have already been presented. pQ \ V3.C9 x 1 0 -51) ' ( —1.262) + 1 \ f 1.76 x 10-54 x (1.78 x 10“52)2\ f 24 = k = ) ^ 3.31 x 1 0 -106 ) { ' } 0.37 x 1.69 x 10“ 52 x 2.71 x 105Vo, (A. 193) 1.69 x 1 0 "V o. (A. 194) Finally, the right fluid pressure is taken as pr = 0.0, since the thesis model uses a pressure ramp with zero on the right. The left pressure is calculated from p = ( 7 — 1)pi£v, where eu is the the fluid internal energy at the surface of the neutrino-sphere, where the shock tube sits. This works out to pi = 1.6 x 108 kgkms~2 . A.7 Development Environment Rather than use a text editor to write Java, and then compile on the command line, an Integrated Development Environment (IDE) is used. Both Netbeans and Eclipse are good choices when it comes to a Java IDE, they both offer beginners a host of features and plugins that will make learning Java a lot easier. Also, they provide auto-completion, which is essential as Java’s list of libraries are massive, and it is not possible to know all of the methods and classes. An IDE is indispensible for this main reason. In the end, the popular choice of Netbeans was made. The Sec. A.8 Comments and Notes on the Programming Process 170 large community backing and providing support for Netbeans, and its ease of use, made the choice trivial. A.8 Comments and Notes on the Programming Process In the course of writing this code, many different systems were used before settling to the current setup. Initially, a Linux environment was used, where code was written in C++ using Vim, and then late KATE, the KDE Advanced Text Editor. KATE has many great settings for code sequence identification by colors and indents. KATE does this all automatically once it has been initially set up. Its GUI is also very advanced and easy to use. Linux itself was very painful to use. In the end, so much time was being spent getting Linux to work all the time and with different software when such was needed, that it was deemed useless to work with, since the actual work time was minimized. The next step was to go to a Mac system, using MacOS Snow Leopard. The Mac turned out to be excellent in terms of installing new programs and having everything work right out of the box. LTEXran flawlessly, and XCode was a dream. However, this was when Leopard was being used, which is 32 bit. Snow Leopard is 64 bit, and when it was installed, everything crashed (everything non-Mac, that is, free). It was then back to Linux-like frustration with trying to get things to work, until finally the Mac was abandoned. The current setup is a Windows XP system, where NetBeans is used. Windows is hated in the business, but in this case it is found to be the best thing to work with. Everything works! Windows XP is stable and is probably the best Windows OS ever. NetBeans does all the right things, with no crashes to date. Java is seamless, and it is easy to generate a .jar file in NetBeans, which can be run on the command line on the university’s supercomputer. In fact, this is what is being done. Notice that early on C++ was mentioned as the language for the thesis code. Now Java is being used. This switch was only due to frustrating problems with segmentation faults in C++. After years of being unsuccessful in dealing with these, the switch to Java was made, since Java handles memory leaks and garbage collection. Actually, it is possible to generate an unhandled memory leak in Java, but this only happens when something really silly has been done, which should not be done in the first place. Sec. A.9 A.9 The Parameter File 171 The Parameter File The evolution code and the Riemann solver use information passed by reading the contents of a “param.dat” file. This file is extensive, and utilizes a number of switches which tell the code what to do and how to handle certain groups of data. It is important to understand this file before carrying out various evolutions. Here is an example of it: # ============================================================== # NOTES # : Written by Gregory Mohammed, 31 August, Masters Thesis, 2012 # # There are a lot of switches and data in here. Read through # carefully so that you understand the comments and what each # switch and set of data do and are for. Do some test runs if # you are not clear. # ============================================================= # Indicate choices: # Initial conditions: 0 = smooth shock, 1 = shock # ============================================================= initShock = 1 # ============================================================= # Show initial slice: Show = 1, Don't show = 0 # This also controls the plotting of an exact Riemann solution # on the evolution graph. If 0 the exact solution is not shown. # This does not work in version 6.0.0. Just ignore the "exact" # solution, # if you set showInitialSlice to 1. = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = showInitialSlice = 1 showExact = 0 # ============================================================= # Method: Upwind = 0, Godunov = 1 Sec. A .9 T he P aram eter File # ============================================================ method = 1 # ========================================================== # Neutrinos: Do Not Include = 0, Include = 1 # = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = == = = = = = = neutrinos = 1 # = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = # Tunes the value of the neutrino flux and energy term. # Make both 0.0 when neutrinos = 0. # ========================================================== #fluxTuner = 1.0e-2 fluxTuner = 0.0e-0 #energyTuner = 0.0e-0 energyTuner = 1.0e-2 # ========================================================== # High Resolution Method: ZeroSigma = 0, MinMod = 1 , MC = 2 # ========================================================== highResType = 1 # ========================================================== # Constant = 0, Linear = 1, CubicHermite = 2 # ========================================================== reconstructionType = 1 # = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = # Resolution, measured in the number of cells. # ========================================================== resolution = 400 # =========================================================== # Test file with parameters for the program Godunov-v.2.0.0, # which is written in Java 17 2 Sec. A.9 # = = = = = = = = = = = = = = = = m The Parameter File 173 = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = ==== = = = = = = = = gamma = 2.0 # ============================================================= # for gamma=2, k=epsilon/rhoO # ============================================================= k = 0.125 # =======================™ ==============================™ == # number of evolutions # ============================================================= totalSteps = 10000 # totalSteps = 1 # ============================================================= # Total time for the Sod tube. # # The Sod values are hard-coded for the Sod shock tube. # This allows the param.dat file to be flexible for "playing" # with data to observe results for different physical situations # (or unphysical ones t o o ) . # # THE ONLY VARIABLE WHICH NEEDS TO BE EDITED HERE IS THE # totalTimeSod (below). # Set to 0.25 when dataSwitch = 0. These times are already # geometrized. # Use 250.0 for resolution = 10000. # = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = totalTimeSod = 0.25 # = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = m = = ™ = = = = = = = = = = # Total time for the Kuroda et a l . data. This gives a good # evolution which shows good data. Sec. A.9 The Parameter File 174 # 500000 <= totalTimeKKT <= 100000000 is a good range. # # Set showInitialSlice = 0. Otherwise the KKT data doesn't show # details in the graph. # vprofile can be what you want to see. # = = = = = = = = = = = = = = = = = « = = = = = = = = = M = = = = = = = = = = = = = = = = S = = = = M = = = = M = = = ItotalTimeKKT = 750000.0 totalTimeKKT = 1000000.0 # ============================================================= # dump control: allowed change in evolution variables # ============================================================= epsdmp = 1.0 # = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = # interval between dumps # = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = dmpinterval = 5000 # = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = # number of correctors # = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = corrector = 1 # =================================«==========*=============== # artificial viscosity # ====================================r========m=============== artvis.kl = 0.0 artvis.k2 = 0.0 # = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = # courant # initial delta = pl*freefall # = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = Sec. A.9 The Parameter File 175 #courant = 0.4 courant = 0.4e-0 # pi small for KKT, and Sod with total time > 0.25. pi = 0.4e-0 # pi = 0.4 for classic Sod. #pl = 0.4 # = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = # dataSwitch: # Sod = 0, Experimental Sod = 1, Neutrino Model = 2 = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = dataSwitch = 2 # = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = # If plotGeo = 1, then plot geometrized units, otherwise plot # the ungeometrized units. # = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = plotGeo = 1 # = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = # Data for the Sod shock tube. Feel free to experiment here, # the classic Sod data are hard-coded, # here. as and unaffected by changes If dataSwitch = 0, then the classic Sod tube will be # plotted using the hard-coded data. # your data will be plotted. If dataSwitch = 1, then If dataSwitch = 2, then the KKT # data will be plotted. # = = = = = = = = = = = = = = = = = = = = = = = = = = == = = = == = = = == = = = == = = == ;= ;= = = == = = = == = = = x_left = -0.3 x_right = 0.3 # = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = v_left = 0.0 v_right = 0.0 # = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = Sec. A.9 The Parameter File rhoO_left = 1.0 rho0_right = 0.125 # ============================================================ p_left = 1.0 p_right = 0.1 # = = = = = = = = = = = = = „ = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = „ = # This data below is ONLY for the Kuroda et al. data against # the data used in the exact Riemann solver. # = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = # Neutrino data from Kuroda et a l . dataSwitch = 1 # = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = # Here use cgs units. # = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = xLeft = 8.0e6 xRight = 1.0e8 # = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = vLeft = 1.0e7 vRight = -1.0e7 # = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = lrhoOnu = 2.0el4 rrhoOnu = 1.0e9 # ========= = = « = = = = = = = = = = = = = M = = = = = = = = = ™ = = M = = = = = = = = = = = « = lpnu = 1.07e8 rpnu = 0.0 # = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = # Make different velocity profiles. Other variable profiles # are hard-coded as noted below. # vprofile = 1 => Discontinuous, velocity # rhoO and p are steps. is zero, and This is compared with the 176 Sec. A.9 # 177 the original Sod data. # vprofile = 2 => parabolic, left increasing, right decreasing. # Other variables are negative ramps. This is # compared with the Sod tube for similar input. # The Parameter File = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = vprofile = 2 # = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = # Mass of the neutronized core, or black hole. # This is specified by the user, and may be 0 if the user does # not want to consider it. However, doing so zeros out the # ad-hoc gravity term, which leads to a numerical explosion, # which is non-physical. # # Mass of the Sun = 1.98892e30 kg # Geometrized Mass of the Sun = 4.5 km # = = = = = = = = = = = = = _ ================================================== #massCore = 4.5e-14 massCore = 0.0 # = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = As can be seen, there are numerous comments in the file, which are intended to help the user understand the purpose of each value and switch, init Shock cab take two values, one o f either 0 or 1. A “0” indicates that the user wants to evolve a smoothed shock, while a “ 1” indicates the evolution of a discontinuous initial condition, ie, a Riemann problem. This thesis is only concerned with case 1. showInitialSlice is either a 0 or 1. “0” indicates that the initial time slice data is not printed to a file, or anything else. A “ 1” shows the initial slice. In the case where the user wants to plot the exact Riemann solution and the evolution data on the same axes, then this switch must be set to “ 1” . However, the user would only want to do this for the Sod data. The exact solution for the data of interest in this thesis is not well defined, and so this switch should be “0”. Sec. A.9 The Parameter File 178 method can be set to “0” which uses an upwinding numerical scheme, or “ 1” which implements the Godunov numerical method. For this thesis, method = 1. neutrinos is a switch which, if “0”, turns off the neutrino equations and the code runs with no neutrino terms. If it is “ 1”, then the neutrino terms are implemented. f luxTuner and energyTuner are values which are type double, and are used to to “tune” the neutrino flux term and the neutrino energy term. These combinations add versatility to the code for the purposes of testing and experimentation. highResType is a value which determines the type of high resolution method to be implemented. “0” turns off high resolution, “ 1” switches to a minmod method and “2” uses a MC method. reconst ructionType is a switch which determines the accuracy of the reconstruction of the Godunov solution. “0” employs no reconstruction, “ 1” employs a linear fit and a “2” is a cubic hermite reconstruction, resolution is an integer which gives the number of cells in the evolution grid, gamma and k are values which may have specific values, gamma is | for non-relativistic fluids and | for relativistic fluids. If gamma is 2, then k is just totalSteps is a value which gives the total number of time steps for the evolution. This value is ad-hoc, and some experimentation is needed to find a value which provides a good compromise between the evolution results obtained and the time to run the program. It should not be so short that no evolution is observed in the graphs, nor too long that the evolution progresses beyond informative data. totalTimeSod and totalTimeKKT are values which set the total runtime of the evolution in geometrized units. The Sod time is 0.25 for the Sod data. The totalTimeKKT is so named because it pertains to the data obtained from the work of Kuroda et al. ([2]). Its value is unknown, and can only be obtained by experimentation. epsdmp = 1. 0 is a value which is not used, so leave this set to “ 1” . The dmpinterval is an integer which determines the time step interval on which data is dumped to a file. In this case, every 1000th time step is output to a file which has its name coded with the step value, corrector determines the number of correctors used in the evolution. This is primarily used in the upwinding scheme, which is not used in this thesis, artvis .kl and artvis .k2 are parameters for the artifical viscosity implemented in the evolution. This is never used, so just leave Sec. A.9 The Parameter File 179 them set to “0” . The courant and pi fields are used to set the courant number and a term (pi) which together allow for huge total times to be used. This is necessary for the Kuroda et al. data. dataSwitch can be one of three values. “0” uses the hard-coded Sod data to plot the Sod solution against the exact Riemann solution. For this to work, showInitialSlice = 1, totalTimeSod = 0.25 and dataSwitch = 0, and vprofile = 1. “ 1” switches to an experimental Sod, where the user can experiment with the values of x_left, p_left, x_right, rho0__left, rhoO_right and p_right with totalTimeSod = 0 .25 and showInitialSlice = 1 to plot their values against the exact Riemann solution to the hard-coded Sod data. “2” switches to the Kuroda et al. data, set in xLeft, xRight, vLeft, vRight, lrhoOnu, lpnu, rrhoOnu and rpnu. In this case use only the totalTimeKKT to set the runtime. When using the Kuroda et al. data, the vprofile is usually “2”, but can be “ 1”. The mass of the neutronized core, or black hole, mas sCore, is a normalized mass to the usual mass of the Sun. This value is necessary to provide the ad-hoc gravitational potential which serves to provide a gravitational field in this toy model. If massCore = 0, then this potential is zero, and does not affect the code. The purpose is to investigate the effect of a gravitational field on the neutrino flux. Bibliography [1] J. A. Pons, J. M. Marti, and E. Muller, “An exact riemann solver for multidimensional special relativistic hydrodynamics,” 1999, arXiv:astro-ph/9910462vl. [2] T. Kuroda, K. Kotake, and T. Takiwaki, “Fully general relativistic simulations of core-collapse supernovae with an approximate neutrino transport,” 2012, arXiv:astroph/1202.2487v 1. [3] M. Shibata, K. Kiuchi, Y. ichiro Sekiguchi, and Y. Suwa Progress o f Theoretical Physics, vol. 125, pp. 1255-1287, 2011. [4] A. Burrows and T. A. Thompson, “Neutrino-matter interaction rates in supernovae,” 2002, arXiv:astro-ph/0211404vl. [5] R. J. Leveque, Finite Volume Methods fo r Hyperbolic Problems. New York: Cambridge University Press, 2002. [6] J. R. Wilson Astrophysical Journal, vol. 173, p. 431, 1972. [7] J. R. Wilson, Sources o f Gravitational Radiation. Cambridge, United Kingdom: Cambridge University Press, 1979. [8] J. Centrella and J. R. Wilson Astrophysical Journal, vol. 54, p. 229, 1984. [9] J. Hawley, L. Smarr, and J. R. Wilson Astrophysical Journal Supplement Series, vol. 55, p. 211, 1984. [10] P. J. Mann Computational Physics Communication, vol. 67, p. 245, 1991. 180 Bibliography 181 [11] P. J. Mann Computational Physics Communication, vol. 107, p. 188, 1993. [12] P. J. Mann Computational Physics Communication, vol. 75, p. 10, 1993. [13] S. A. Colgate and R. H. White The Astrophysical Journal, vol. 143, p. 626, 1966. [14] S. W. Bruenn, K. R. de Nisco, and A. Mezzacappa The Astrophysical Journal, vol. 560, pp. 326-338, 2001. [15] E. Tadmor, “Approximate solutions of nonlinear conservation laws and related equations,” 1997. [16] M. May and R. H. White Methods Comput. Phys., vol. 73, p. 219, 1967. [17] R. L. Bowers and J. R. Wilson Astrophysical Journal Supplement Series, vol. 50, pp. 115— 159, 1982. [18] P. Anninos and P. C. Fragile The Astrophysical Journal Supplement Series, vol. 144, pp. 243257, 2002. [19] L. Rezzolla and O. Zanotti, “An improved exact riemann solver for relativistic hydrodynamics,” 2001, arXiv:gr-qc/0103005v2. [20] D. L. Tubbs and D. N. Schramm The Astrophysical Journal, vol. 201, pp. 467-488, 1975. [21] C. Levermore Journal Quantum Spectroscopy Radiation Transfer, vol. 31, pp. 149-160, 1984. [22] K. S. Thorne Mon. Not. R. astr. Soc, vol. 194, pp. 439-473, 1981. [23] R. Arnowitt, S. Deser, and C. W. Misner Physical Review, vol. 116, 1959. [24] H.-T. Janka Astronomy and Astrophysics, vol. 368, pp. 527-560, 2001. [25] T. D. Matteo, R. Perna, and R. Narayan The Astrophysical Journal, vol. 579, pp. 706-715, 2002 . Bibliography 182 [261 T. Liu, W.-M. Gu, L. Xue, and J.-F. Lu The Astrophysical Journal, vol. 661, pp. 1025-1033, 2007. [27] D. Zhang and Z. G. Dai The Astrophysical Journal, vol. 703, p. 461, 2009. [28] R. Haberman, Applied Partial Differential Equations. United States of America: Pearson Education, 2004. [29] G. A. Sod, “A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws,” Journal o f Computational Physics, vol. 27, pp. 1-31, 1978. Index 3 + 1 split, 19 gain radius, 8 3-vector notation, 3 Godunov’s method with anti-diffusion, 39 4-vector notation, 3 accretion shock, 8 ADM formalism, 19 advection, 34 High Resolution Methods, 38 hyperbolic, 34 lapse function, 19 Lax-Wendroff Theorem, 40 advection equation, 34 maxmod, 42 bounce shock, 8 cell definition, 4 characteristic, 34 minmod, 42 neutrino energy density, 24 neutrino flux, 24 Conservation of energy-momentum, 16 neutrino pressure, 24 Conservation of rest mass, 16 neutrino-sphere, 2 conservation of total energy-momentum, 18 notations and conventions, 3 Courant condition, 36 Courant number, 40 covariant derivative, 4 order of the error, 56 piecewise linear reconstruction, 39 dot notation, 3 prompt bounce-shock mechanism, 113 Einstein summation convention, 4 radiation constant, 24 entropy conservation, 41 Riemann problem defined, 36 exact Riemann solver, 10 shift vector, 19 First-Order Upwind Method, 35 special relativistic fluid hydrodynamics, 17 forward differenced Euler scheme, 34 Stefan-Boltzmann constant, 24 183 Index superbee, 42 Total-Variation Diminishing, 41 184