. 1 . A ..
Rhfismifiwm .5

A

.m.
. firmunmaw m3.
. rwmfiflwm.wum “ha:
“Nita

u n
dim ..
.efifi. Infifiuw.

 

 

 

«.25.. an,

 

 

 

 

 

LIBRARY
Michigan State
University

 

 

 

This is to certify that the

thesis entitled

Numerical Study of Rotating
Core Collapse Supernovae

presented by

Mark Tobias Bollenbach

has been accepted towards fulfillment
of the requirements for

M. S. degree in Physics

 

 

(lit V
5 Q Mg” Wt

Date July 22, 2002

 

0-7639 MS U is an Affirmative Action/Equal Opportunity Institution

 

PLACE IN RETURN BOX to remove this checkout from your record.
To AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

6/01 c:/ClRC/DateDue.p65-p. 15

NUMERICAL STUDY OF ROTATING CORE COLLAPSE SUPERNOVAE

By

Mark Tobias Bollenbach

A THESIS

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

MASTER OF SCIENCE

Department of Physics and Astronomy

2002

ABSTRACT

NUMERICAL STUDY OF ROTATING CORE COLLAPSE SUPERNOVAE

By

Mark Tobias Bollenbach

The explosion mechanism of core collapse supernovae is still far from being un-
derstood. In this work, an overview of the current understanding of core collapse
supernovae and the history of numerical simulations that helped develop it is given.
While it is widely believed that neutrino heating and convection above the neutrino
sphere are the key processes to revive a stalled shock and thus obtain successful explo-
sions, it is still possible that other phenomena are crucial for the explosion mechanism.
For example, recent observations of the polarization of the light emitted by super-
nova explosions indicate that there are large deviations from spherical symmetry in

the very heart of the explosion.

In contrast to most of the previous simulations which were performed in one or
two dimensions, we use the different approach of a three dimensional test particle
based simulation. The underlying microphysics is crudely simplified to make this

computationally possible.

A systematic study of the influence of rotation mainly during the infall phase of
the collapse of a typical iron core is performed. Different equations of state and initial
conditions are used. Indications for significant deviations from spherical symmetry

are found in our simulations.

ACKNOWLEDGMENTS

Most importantly, I would like to thank Professor Wolfgang Bauer for giving me
the opportunity to work at the National Superconducting Cyclotron Laboratory and
for being my advisor for this research project during this fantastic year at Michigan
State University. I enjoyed the relaxed work environment and the beauty of this place

very much.

Further, I would like to thank the German National Scholarship Foundation, the
Studienstiftung des Deutschen Volkes, for the support I obtained during this year and

for an unforgettable meeting in Washington, DC. in April 2002.

Finally, I thank my family for supporting me throughout the years and thus

making all this possible.

iii

Contents

LIST OF TABLES vii
LIST OF FIGURES viii
1 Introduction 1
2 Overview of Type II Supernovae 3
2.1 Presupernova Stellar Evolution and Supernova Progenitor for a 15M®
Star ....................................
2.1.1 Stellar Evolution .........................
2.1.2 Supernova Progenitor ....................... 6
2.2 Explosion Mechanism of Type II Supernovae .............. 10
2.2.1 Collapse .............................. 10
2.2.2 Core Bounce ............................ 14
2.2.3 Prompt Shock Mechanism .................... 15
2.2.4 Delayed Shock Mechanism .................... 17
2.3 Recent Numerical Studies of Core Collapse Supernovae ........ 22
2.3.1 Equation of State ......................... 23
2.3.2 One Dimensional Simulations .................. 26
2.3.3 Two Dimensional Simulations .................. 29
2.3.4 Three Dimensional Simulations ................. 33
3 Our Three Dimensional Simulation 35
3.1 Test Particle Method ........................... 36
3.2 Grid .................................... 37
3.2.1 Cells and Boundaries ....................... 37
3.2.2 Grid Coordinates ......................... 38
3.2.3 Calculation of Densities ..................... 40
3.2.4 Calculation of Derivatives .................... 42

iv

3.3 Equation of State .............................
3.3.1 Cold Nuclear EOS and Polytrope EOS .............
3.3.2 Calculation of Thermodynamic Quantities other than Density
3.3.3 Helmholtz E08 and Lattimer & Swesty EOS ..........
3.3.4 How the EOS Affects the Dynamics ...............

44
44
45
47
49

3.4 Symmetry Assumptions, Boundary Conditions, and Numerical Problems 52

3.4.1 Symmetry Assumptions .....................
3.4.2 Derivatives at Grid Boundaries .................
3.4.3 Background Density .......................
3.4.4 Singularity Treatment ......................
3.5 Time Development ............................
3.6 Calculation of Observables ........................
3.7 Advantages and Weaknesses of this Method ..............
3.8 Computational Requirements ......................
3.9 Output Possibilities ............................
3.9.1 Density Output ..........................
3.9.2 Test Particle Output .......................

Numerical Results and Interpretation

4.1 Angular Momentum Conservation ....................

4.2 Results of Simulation Runs ........................
4.2.1 Cold Nuclear EOS without Electron Contribution .......
4.2.2 Cold Core with Polytrope Electron Contribution ........
4.2.3 Combination of Helmholtz E03 and Lattimer & Swesty EOS .

4.3 Possible Implications of these Results ..................

Summary and Conclusion

Source Code of the Simulation Program

A.1 suno.cpp and suno.h ..........................
A.2 testparticle.cpp and testparticle.h ................
A.3 vector_and-spherical.cpp and vector-and-spherica1.h ......

Source Code of the Output Programs
B.1 output.vbp ................................
B.2 Suno Density 0utput.vbp .......................

C Abbreviations and Symbols

52

54
55
56

57
58
58

59
59
59
60
66
72
77

79

83
83
125
128

133
133
140

148

BIBLIOGRAPHY 151

vi

List of Tables

2.1

4.1

4.2

4.3

4.4

4.5

4.6

Iron core masses of different supernova progenitors (taken from [56],
all masses in units of MG) ........................

Simulation runs using the cold nuclear EOS. The abbreviations are
explained in the text ............................

Key for figure 4.3. to through te are the times (in ms) corresponding
to the density profiles labeled (a) through (e) in the figure. ra through
re are the corresponding radii (in km) of the shown density plots.

Simulation runs using the combination of the cold nuclear EOS and the
polytrope EOS for the electron gas. The abbreviations are explained
in the text. ................................

Key for figure 4.6. ta through te are the times (in ms) corresponding
to the density profiles labeled (a) through (e) in the figure. ra through
re are the corresponding radii (in km) of the shown density plots.

Simulation runs using the combination of the Helmholtz and the LS
EOS. The abbreviations are explained in the text ............

Key for figure 4.8. ta through te are the times (in ms) corresponding
to the density profiles labeled (a) through (e) in the figure. pmmin and
psama:c are the densities (in g/cm3) corresponding to the bottom end
(color: light yellow) and top end (color: black) of the density key.

vii

61

64

72

74

List of Figures

2.1
2.2

2.3

2.4

2.5

3.1
3.2
3.3
3.4
3.5
3.6

4.1
4.2

4.3

Structure of a 15M® star at the onset of core collapse (taken from [56]) 7

Composition of a 15M® star at the onset of core collapse (taken from
[56]) .................................... 8

Infall velocity of the core matter V and local sound velocity A approx-
imately one millisecond before core bounce (taken from [4] who used

the results of [1] to create this plot). . . . .............. 13
Shock revival by neutrino heating after the shock wave has stalled
(taken from [24]). The symbols are explained in the text ........ 18

Angular velocity profile used by Fryer and Heger [18] compared to that
of different models used by Monchmeyer and Miiller (labeled “MM89”)

[35] ..................................... 30
Smearing of test particles for density calculation ............ 41
Cold nuclear EOS and polytrope EOS for electron gas ......... 45
T(p) for the infall phase ......................... 47
Ye(p) for the infall phase ......................... 48
Combination of Helmholtz EOS and LS EOS used in our simulation . 49
Test particle j traveling from one grid cell to another ......... 51

Angular momentum as a function of time for a typical simulation run 60

Time development of the different energies in the late stages of simu-
lation CNEOSO5 ............................. 62

Mass density in a slice in the x-z-plane at five different times for models
CNEOSOB (left), CNEOSOS (center), and CNEOSO7 (right). Events at
the respective times: (a) onset of simulation, (b) clumping and presence
of centrifugal forces become apparent, (c) formation of “vortices” in
CNEOSO7, ((1) core bounce, (e) shortly after core bounce. Note that
the radii of the shown density profiles vary (see table 4.2 for more
data). The black lines below each plot indicate a length of 43.5 km.
The length scale for the plots labeled by (a) is much larger than this.
Images in this thesis are presented in color ................ 65

viii

4.4

4.5
4.6

4.7
4.8

4.9

Initial mass density as a function of the radius calculated by Woosley
and Weaver [56] (solid line) and the variant contracted by a factor of
5 we use in our simulation (dashed line).The location at which the
enclosed mass is 0.7M® (edge of the inner core) is indicated for both
distributions on the abscissa ........................

Time development of the different energies in simulation PCNEOS19

Mass density in a slice in the x-z-plane at five different times for mod-
els PCNEOSO7 (left), PCNEOSI3 (center), and PCNEOSIQ (right).
Events at the respective times: (a) onset of simulation, (b) after 2 ms,
(c) presence of centrifugal forces and vortices become apparent, ((1)
core bounce, (e) shortly after core bounce. Note that the plots have
different radii. The black line below each plot indicates a length of 56.1
km. See table 4.4 for more data. Images in this thesis are presented in
color .....................................

Time development of the different energies in simulation HLSEOSl3.

Mass density in a slice in the x-z-plane at five different times for models
HLSEOSO7 (left), HLSEOS16 (center), and HLSEOS22 (right). Events
at the respective times: (a) onset of simulation, (b) after 2 ms, (c)
presence of centrifugal forces becomes apparent (after 3 ms), ((1) core
bounce, (e) shortly after core bounce. The density scale is only approx-
imately valid: the highest density (indicated by the colors black and
dark red) decreases from left to right. All plots have the same radius
(123.42 km) indicated by the black line in the top left. See table 4.6
for more data. Images in this thesis are presented in color. ......

Density depletion along the z-axis after bounce in model HLSEOS22
at t = 5.21ms. The radius in this plot (measured from the center
horizontally to the right) is 51.43km. Images in this thesis are presented
in color. ..................................

ix

67
69

71
73

75

Chapter 1

Introduction

Supernova explosions are among the most spectacular phenomena that we know of.
It seems that the supernova explosion is also one of the most challenging phenomena

as far as the understanding of the underlying physics is concerned.

The most famous recent supernova was observed in 1987 (therefore labeled SN
1987A). It could be seen with the naked eye and is particularly interesting because it
was the first supernova whose progenitor star had been observed before the explosion.
While such observations of supernovae (that we know of) date back almost 2000 years,
the theory of core collapse supernovae has rapidly developed in the last decades.
Several different explosion mechanisms have been suggested throughout the years.
Most of these turned out to be fundamentally wrong. Apart from being remarkable
optical events, supernovae are believed to play an important role in the synthesis of
heavier elements in the mass range 16 5 A S 60 —— core collapse supernovae are
especially relevant for the synthesis of oxygen. They are possibly also the origin of

the mysterious 7-my bursts [8].

In chapter 2 of this work the state of the art of the knowledge of the physics of
core collapse supernova explosions will be briefly reviewed. The relevant concepts and

processes will be explained. Due to the complexity of the problem numerical studies

play an outstanding role in the understanding of supernova events. An overview of
the history of these studies with a special focus on recent developments will be given.

Special attention will also be paid to the role of rotation in supernovae.

Our new approach to simulate core collapse supernovae will be presented in chap-
ter 3. The numerical techniques used, the necessary approximations and assumptions,
the implementation, and the differences of this method in comparison to previous sim-
ulations will be explained. The strengths and weaknesses of this technique will be
pointed out. At the current state of the art our method can be considered a good
model of reality only during the infall phase, i.e. until core bounce (these terms will

be explained in chapter 2).

In chapter 4 the results of several simulation runs (using different equations of
state for the core matter and a variety of initial conditions) that were performed
using the technique described in chapter 3 will be presented. Conclusions about the
effects rotation may have on supernovae in reality will carefully be drawn. Possible

improvements to our model will be suggested and discussed.
This whole work will be summarized and reflected in chapter 5.

A typographical convention will be to write important terminology in italic type

at the place where it occurs for the first time in this work.

We will further use standard symbols like 0 for the speed of light or .MG for the
solar mass in text and equations. Other abbreviations and symbols will be introduced
and used throughout this work. For convenience a summary of all these is given in

appendix C.

Images in this thesis are presented in color.

Chapter 2

Overview of Type II Supernovae

In this chapter the physics of type II supernovae will be briefly reviewed. Supernova
explosions are labeled either type I or type 11 (yet both with several subdivisions).
This distinction was established due to observable differences: the most important
difference between the two is the absence of hydrogen lines in the spectrum of type Is
while these are present in the case of type 115. The light curves1 of the two phenomena

also differ significantly.

From a theoretical point of view, type Is and 113 are almost utterly distinct phe-
nomena. While the power source of type Is is believed to be thermonuclear burning
initiated by the exceeding of a white dwarf ’32 Chandrasekhar mass3 due to the ac-
cretion of matter from a companion star, type 113 are powered by the gravitational
energy released during the collapse of a star’s iron core (note, however, that there are

also subclasses of type Is that are powered by core collapse).

We will only deal with type 118 in this work and most of the time only with the
typical case of a star in the ZAMS4 mass region z 15MQ. For more information

about type Is one may refer to [56, 11, 4]. Most of what will be said here remains

 

1i.e. the luminosity of the emitted light as a function of time

2a very dense star with a mass of approximately 1M0 but only the size of earth

3the maximum mass a white dwarf can have without collapsing (z 1.4MQ)

4Zero Age Main Sequence: denotes the star’s properties at the beginning of its stellar evolution

valid for stars with a ZAMS mass between 2: 11M@ and a: 40M®. This separation
is necessary because the pr0perties of stars during and at the end of their stellar
evolution are essentially determined by their ZAMS mass (see [11, 4, 55, 20, 52] for
details on stellar evolution). A very light star like our sun for example will never
develop an iron core during its evolution and thus never produce a core collapse
supernova. Also note that for stars in the ZAMS mass region 8M0 < Mata, < 11M®
type II supernovae are believed to occur but the details of these stars’ presupernova
evolution and the supernova event itself depend sensitively on their ZAMS mass and
are somewhat different than what will be dealt with in the rest of this work (see e.g.
[56]). Finally, stars more massive than z 40M® are believed to lose their hydrogen

mantle before the end of their life which disqualifies them as type 113.

For simplicity, we will refer to a star with ZAMS mass ...M® as just a ...M®

star from now on.

2.1 Presupernova Stellar Evolution and Supernova
Progenitor for a 15M® Star

As already mentioned the evolution of a star during its so-called main sequence5 is
mainly determined by its ZAMS mass. It is way beyond the scope of this work to
give a detailed account of the events and processes during the main sequence. Instead
we will briefly follow the evolution of a 15M® star to the point shortly before core

collapse occurs.

 

5the long phase (possibly lasting billions of years) of the star’s life during which hydrogen-burning
is its dominant power source

2.1 .1 Stellar Evolution

Just like our sun all stars “burn” hydrogen from the beginning of their main sequence,
i.e. sequences of nuclear reactions take place which ultimately lead to the fusion of

hydrogen to helium, e.g. by the proton-proton chain which results in the reaction

4]H—+§He+2e++2I/e+27. (2.1)

The intermediate steps have been omitted here. Other reaction chains also con-

tribute to the conversion of hydrogen to helium.

It is well known that during this fusion of hydrogen a lot of energy is released
because (for nuclei lighter than iron) the binding energy per nucleon increases with
increasing mass number of the nucleus. However, for fusion to occur the Coulomb
barrier of the (positively charged) participating nuclei (e.g. two protons) has to
be overcome which requires them to have large kinetic energies. Such energies are
available in the interior of stars where the temperature during the main sequence is

of the order of magnitude 107K.

After significant amounts of helium have been produced and only if the tempera-
ture in the core of the star has become sufficiently high helium itself starts “burning”

by the triple alpha process:

§He+§He +———> fiBe (2.2)

fiBe + 3H8 —> 13,0 + 7.

Similarly, given the needed prerequisites, “burning” of (most importantly) H, He,
C, Ne, O, and Si occurs successively in the center of the star — sometimes even in its
shells. Because of the increasing Coulomb barriers higher and higher temperatures

(e.g. of the order of 4 x 109K for Si-burning) are needed for these fusion reactions

and less and less energy is gained per participating nucleon. After carbon ignition the
energy losses of the star are huge (compared to the previous stages) and dominated
by neutrino emission since then the temperatures are high enough for the occurrence

of (among other processes) the electron capture reaction:

10+ + e- ——> n + V6. (2.3)

The neutrinos created in this reaction hardly interact with the matter of the star and

just radiate away. This way, energy is “carried” out of the star.

Thus the star develops a gigantic power using a much less effective power source
than before. Consequently the time scale of the burning stages becomes smaller with
increasing charge number of the fuel. For example, in a 20M® star the H burning lasts
for approximately 107 years, He burning 106 years, C burning 300 years, O burning

200 days, and Si burning merely 2 days (these numbers were taken from [11])!

The final product of these nuclear burning processes is either 56Ni or 54Fe after

which no more energy can be gained by fusion.

2.1.2 Supernova Progenitor

It shall be mentioned that the knowledge of main sequence stellar evolution is very
sophisticated and in fact much more profound than that of the supernova phenomenon
itself. An incredible effort has been put into numerical simulations of stellar evolution
and there is wide consensus on the development which stars pass through during their
main sequence. These efforts converged to the initial conditions for the core collapse
event of a 15M© star illustrated in figures 2.1 and 2.2 and calculated by the simulations

of Woosley and Weaver [52, 56].

I)” .l 501

 

 

 

 

 

 

 

 

Q
r . i '
O.
c. v
P
l" '7
M
D
a E s
o __ '2 ‘° Q
g " v- v-
: x
m X 8 X
_ 35. .. a
m " N‘
I " n
1: >34
—- °.
:9
at. L“?
v
it
71 "
“o
g 2
g E
:1
_ 3 3
a E
a
a b
u .9
a“. 8
7'; *- e S E
S 2
S In
P
P
i.— P Q
N
D
3 ,
U
JL 1'
ll' 1"— r
; /
s“ ’ -o
q ,7’ r
2 u
8 colt/w:
[' Jsfi -4
I
H
p. I
1 l 1 14’: l 1 L L l 1 1 1 l 1 1 4 o

 

n-3,-“ 510 nor) 1 '9 50|.'(€_m3 a) doc-l

Figure 2.1: Structure of a 15M® star at the onset of core collapse (taken from [56])

 

YITTI T I FTTIT T T I ITIYVTW T l
.—

 

 

 

 

 

 

 

 

 

 

H P-
‘3.
P'- d‘
P
f‘ -r
5
c 2""? o
a
o 0 Z _.
I I I08 0 v If!) 2
'- ¢ PE '- F p
O
O
0
41b In
1 i'v
t' z .1
3

 

 

16
16 o 0
2° No I
ll
interior mm (NI/Mo)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

‘
if I
Q
o
_.2 .5.’ 5 2 °
. o 0 " ' ““5
2 $3 "
o :r: c,
= K 1
,t _ \ //
s. O
r—- o a 0) < 8 CI 0-“ ’
:e s 8 3 .5 ‘L N
.2 N iii
3 x
_ it _.o.
a I—
O
: IL
0
p} S «
11111 I ILLILJLJ 1 [1111 L1 1 I o
O. '7 ‘2‘ 'i’
F
9 S .9

none.» 352“

Figure 2.2: Composition of a 15M® star at the onset of core collapse (taken from
[56])

Note that their progenitor was calculated using a one dimensional model assum-
ing spherical symmetry and neglecting possible effects due to rotation, which is the
case for almost all stellar evolution simulations (an exception including the effects of
rotation is e.g. [20]). Instead of the distance r = [5| from the center of the star, the
interior mass M (r) in units of M0 is used on the abscissa. If p(r) is the mass density

of the star as a function of the distance from its center

1W(r) :2 /|‘|< d3xp(§:') 2 47r fr dr'r'2p(r'). (2.4)

0

Note that the iron core consisting of all the heavy nuclei (48 S A S 65) formed in
the nuclear burning processes mentioned in 2.1.1 (denoted by “Fe” and 54Fe in figure
2.2) whose collapse will get the main focus in this work ends in a relatively abrupt
way at an interior mass of approximately 1.35MG — a typical value as most stellar
evolution simulations for stars of various ZAMS masses greater than 11M® yield iron
core masses between 7%: 1.1MQ and z 2.5MG (see table 2.1). Apart from the heavy
elements (“Fe”) very small amounts of free neutrons, protons and alpha-particles are
present in the iron core. Above the iron core, lighter elements are stratified in an

“onion-skin” structure.

Before collapse the iron core is almost a perfect 7 = g polytrope, i.e. the pressure

p inside the core is only a function of the density p satisfying

17 0< p7 (2-5)

 

ZAMS mass ]] 12 15 20 25 35 50 100
Iron core mass [I 1.31 1.33 1.70 2.05 1.80 2.45 2.3

 

 

 

 

 

 

 

 

 

 

Table 2.12 Iron core masses of different supernova progenitors (taken from [56], all
masses in units of MG)

with 'y = g. This is just due to the fact that the pressure comes mainly from the

gas of relativistic, degenerate electron gas inside the core for which it is known from

statistical mechanics that the pressure is independent of the temperature and follows
4

equation 2.5 with 7 = 3 (y is usually called adiabatic indem, n z: #7 is called

polytropic index).

As pressure comes mainly from the electrons it is important to mention that the
electron fraction6 Y; in the iron core has dropped (mainly by reaction 2.3) to about
0.44 at the onset of collapse. Reaction 2.3 also leads to a drastic decrease of the
entropy per baryon in the iron core region where it is roughly just 1kg (where [is
denotes Boltzmann’s constant) directly before collapse compared to a value of about
23kg at the beginning of the main sequence. Entropy is “carried” from the iron core
to the envelope by neutrinos where an entropy per baryon of a: 40163 is typical directly

before collapse (all these entropy values are taken from [56]).

2.2 Explosion Mechanism of Type II Supernovae

In this section the predominating theories for the explosion mechanism of a typical

core collapse supernova will be described.

2.2.1 Collapse

During Si burning the iron core obviously gains mass (as iron is produced in this
reaction). This essentially goes on until it reaches a mass that results in gravitational
forces which can no longer be supported by the pressure of the present degenerate
electron gas. Also note that the pressure at the edge of the iron core is not zero but

the outer layers (mantle and envelope) of the star help “squeezing” its core.

 

6the number of electrons per baryon where “baryon” is used here as a generic term for protons
and neutrons

10

Electron Capture and Photodisintegration

As soon as collapse begins, two instabilities are of importance. Ongoing electron
capture reduces the electron fraction in the core thus obviously further reducing the
pressure created by the electron gas. The neutrinos created in the electron capture
reactions are radiated away almost freely (at least before densities high enough for the
occurrence of neutrino trapping are reached) which ultimately results in a reduction of
entropy in the core — a phenomenon known as neutrino cooling. This helps the ongoing
collapse even further as a reduction of entropy leads to a reduction of temperature

which on its part implies a pressure decrease.

The second instability is due to a process called photodisintegration that is pos-
sible at the extremely high temperatures now present in the core: heavy nuclei are
fragmented to their constituents by extremely high energetic photons, for example

(and most importantly) by the following reactions:

SgFe + 7 —> 13 gHe + 4 n (2.6)

§He+7 ——> 2p++2n.

These photodisintegration processes require a lot of energy (as they reverse the
nuclear “burning” reactions by which the star was powered during its whole life).
Therefore the temperature increase due to the increase of density in the core during
its collapse is intensely weakened resulting again in a pressure decrement: gravity
can no longer be compensated by pressure. In stars with ZAMS mass greater than
z 20M® photodisintegration is considered to be the dominant cause of collapse while

for lighter stars electron capture dominates.

11

Inner Core, Outer Core, and Neutrino Trapping

Once core collapse has begun, things develop very rapidly. The iron core matter
falls almost (never quite though) at free-fall velocity towards the center of the star.
Hence, the time scale for collapse is merely z 100ms what justifies the assumption of

an approximately adiabatic process.

Two regions in the iron core must be separated: the inner core and the outer core.
In the inner core, the infall velocity of the matter is proportional to the distance from
the center at any given time which obviously causes all the matter in the inner core to
finally arrive at the center of the star simultaneously. The collapse of the inner core
is homologous in the sense that the radial distributions of all important quantities
(like density, temperature, electron fraction etc.) remain similar to themselves during
collapse — just the respective scales change. A typical mass for the inner core is 0.6
to 0.8M®, i.e. the iron core is almost equally split [56]. The inner core ends at
the distance where the infall velocity of the matter exceeds the local sound velocity.
Figure 2.3 illustrates this showing the infall velocity of the matter and the local sound

velocity as a function of the radius approximately one millisecond before core bounce.

Beyond that (obviously time-dependent) distance the outer core falls towards the
center at supersonic velocities. It has decoupled from the inner core and arrives at

the center of the star later.

Contrary to what was believed before 1979 when the large heat capacities of
excitations of heavy nuclei were found which cause the matter to remain relatively
cool, this collapse does not stop before nuclear density is reached. It actually goes
on till the nuclei touch each other and even beyond that: it is believed that roughly

three times the density of isospin symmetric nuclear matter (i.e. 3 x 2.4 x 101435;, 2

12

 

50

IO

 

l l l l l
5 IO 20 50 I00 200 500

Km

 

 

 

Figure 2.3: Infall velocity of the core matter V and local sound velocity A approxi—
mately one millisecond before core bounce (taken from [4] who used the results of [1]
to create this plot).

7 x 10143513) is reached at maximum compression of the core7 [4, 3]. Note that
the matter is not isospin symmetric in the present situation: significant amounts of

protons have been converted to neutrons by electron capture.

However, due to a phenomenon called neutrino trapping the matter is not as
asymmetric as one might expect at first. Neutrino trapping occurs at densities higher
than z 10116—53. As suggested by its name, it means that neutrinos can no longer
escape the core freely but are trapped in there since at these densities elastic scattering
of the neutrinos by the nuclei becomes relevant. The mean free path for the neutrinos
is so small now that they can only diffuse through this high density region. For the
density mentioned above the time scale for the diffusion of neutrinos out of the core
clearly exceeds the collapse time scale meaning that the neutrinos cannot get out of

the core fast enough — they are trapped.

 

7This value is just the result of most numerical simulations.

13

Note that because of the very low neutrino mass the presence of the trapped neu-
trinos hardly affects the nuclei at all in the sense that their pressure contribution is
negligible. However, they are able to heat the matter by neutrino-electron scattering
events. The cross section for these is (under the given conditions) roughly two or-
ders of magnitude smaller than that for neutrino scattering by heavy nuclei but the
electrons are so highly degenerate that they can virtually only gain energy in these

events (hence the matter can only be heated).

As electron capture results in the production of neutrinos and these cannot es-
cape anymore, their chemical potential rises rapidly thus obstructing further electron
capture. In fact reaction 2.3 occurs in both directions until (dynamic) equilibrium is
accomplished:

19+ + e‘ <—> n + V8. (2.7)

Note that there are also positrons present which enable the reaction
n+e+ +—>p++17, (2.8)

also contributing to the neutrino production. After that, the total lepton fraction8 YL
remains essentially constant (till core bounce, the end of collapse) at a value YL z 0.36
(according to [2]). A realistic value for the electron fraction in the core’s center at

that time is Y}, z 0.3 (according to [13]).

2.2.2 Core Bounce

As soon as a density greater than (the isospin dependent) nuclear matter saturation
density is accomplished, the strong interaction between the nuclei becomes repul—
sive as a consequence of the Pauli principle for neutrons (which are fermions). As

mentioned above, at a density of typically z 7 x 10143:; this repulsion becomes so

 

8the number of leptons (here mainly electrons and neutrinos) per baryon

14

strong that the matter suddenly stiffens, the collapse is halted, and the inner core
rebounds - a little bit like a spring that is first compressed and then released. In
doing so, shockwaves are sent out to the infalling outer core. This event is known as

core bounce.

Theoretical predictions of the maximum density reached at bounce and its vigor-
ousness (i.e. the energy of the created shockwave) are certainly dependent on factors
like the inner core mass, temperature, electron fraction and others but most impor-
tant is the equation of state (EOS) of nuclear matter above nuclear density. The
nuclear matter EOS in the density and temperature region present at core bounce is
still unknown since it is extremely difficult to mimic these conditions experimentally.
However, an enormous theoretical effort has been invested in finding the correct EOS
and several theoretical predictions exist (see e.g. [3, 28, 43, 41, 39, 30]). Numerical
simulations have shown that a “softer” nuclear EOS results in the achievement of
higher densities at bounce and more vigorous shock waves in comparison to those

COmputed using a “stiffer” EOS — a somewhat intuitive result.

2.2.3 Prompt Shock Mechanism

During collapse an enormous amount of gravitational energy is released. The main
question is by which mechanism even a small fraction of this energy can be coupled
to the mantle and the envelope of the star in order to eject them. The gravitational
energy released is so immense that the coupling of just z 1% of it would be ample
to eject the outer layers of the star and explain the observed supernova explosion

energies of as 1051erg 2: lfoe.

An appealing coupling mechanism, called the prompt shock mechanism was fa-
vored till the mid 19808: after core bounce the created shockwave moves outward

through the infalling matter of the outer core. Analytical arguments and numerical

15

simulations suggest that the energy of this shockwave is a: 4 to 7 x 1051erg [26] —
completely sufficient to power the explosion. So, after this shockwave has reached the
outer layers of the star, these get enough kinetic energy to escape, are ejected, and

the star explodes.

Unfortunately, things turned out to be not quite this simple: the problem is that
the shock loses enormous amounts of energy while it beats its way through the infalling
matter of the outer core. The densities and temperatures in the shock region are so
high that again electron capture resulting in neutrino losses (at least when these can
diffuse away fast enough) and photodisintegration occur massively. Large amounts of
heavy nuclei are broken up completely to free nucleons. This costs gigantic amounts
of energy (z 1.5 x 1051erg for every 0.1MQ thus converted) that is taken from the
shock energy. Consequently, if the outer core is sufficiently large, the shock eventually
completely stalls (approximately at a distance of 100 to 300km from the center) and
never even reaches the outer layers — the explosion has failed. Note that possible
energy losses of the shock wave in the outer layers of the star are a lot smaller than
those suffered inside the core because densities and temperatures are too low there
for photodisintegration to occur massively. This is why a successful explosion can be

anticipated as soon as the shock wave managed to leave the iron core.

The success of the prompt shock mechanism depends crucially on the iron core
mass: only if the outer core is small enough (that is essentially the case if the total
core mass is small enough) can the prompt explosion work. Numerical simulations
indicate that prompt explosions can only occur in stars with iron core masses smaller
than z 1.25MQ. Stellar structure calculations indicate that there are stars with
such light iron cores. The majority of type II progenitors however has iron cores too

massive for a successful prompt explosion.

16

The prompt mechanism is more likely to succeed if a “softer” nuclear EOS is used

as it results in greater shock energies.

2.2.4 Delayed Shock Mechanism

So what happens once the shock has stalled? There is still an enormous amount
of energy in the core, mainly in the form of thermal excitations and neutrino and
electron chemical potentials. Due to large number density gradients the neutrinos
diffuse outward. As they reach regions with lower densities their mean free path

becomes larger and larger until they can finally radiate away almost freely.

Consequently, the whole matter above the neutrino sphere9 (located at z 40km
from the center) is bathed in a very intense neutrino flux. This can actually help
to revive the stalled shock that, by now, has turned into a so-called accretion shock
(caused by the infalling outer core matter) in the following way (called the delayed
shock mechanism): the shocked matter above the neutrino sphere now contains a
lot of free nucleons which can absorb some of the neutrinos (the cross section for
these processes is small but not negligible). This ultimately leads to the heating
and expansion of the matter (neutrino heating) so that the shock can resume its way
outward. It shall be mentioned that neutrinos can also be absorbed by nucleons inside
nuclei but (for reasons that will not be discussed here in more detail) are re-emitted

shortly afterwards so that these events do not help in reviving the shock.

 

9the distance at which the optical depth for neutrinos is 1. The neutrinos can be considered to
radiate away freely approximately from there. More precisely, all this is dependent on the energy of
the neutrinos.

17

E.
3%.
Q)
Q a
92, an
0o a:
it) To
O
c.
3
" m

R118
PNS
(convective)

gain
radius

 

I
M
shock

Figure 2.4: Shock revival by neutrino heating after the shock wave has stalled (taken
from [24]). The symbols are explained in the text.

18

Figure 2.4 schematically shows the situation after the shock wave has stalled: the
gain radius is the distance from the center at which neutrino cooling is dominated
by neutrino heating (i.e. more neutrinos are captured than emitted), overturn means
essentially the same as convection, M denotes the infalling mass of the outer layers,
R", is the radius of the proto-neutron star10 (“PNS”), and RV denotes the neutrino

sphere.

This delayed mechanism was first suggested by Wilson in 1985 [53]. Today it is
the commonly accepted theory for core collapse supernova explosions. Still, there is
a lot of trouble with it: it depends crucially on the cross sections for neutrino capture
on nucleons, the neutrino production rates, the details of the neutrino transport, and
other factors many of which are not known with great certainty. In particular, the
simulation of neutrino transport in the “gray” region, i.e. at densities where neutrinos
are neither trapped nor able to freely radiate away is most problematic. Depending on
these input parameters, simulations of the delayed mechanism (performed by several
different groups) yielded successful explosions even for stars with very massive iron
cores while in other cases the explosions “fizzled” again. A more detailed overview of
these simulations will be given in 2.3. A very good summary of the delayed mechanism

and the skepticism about its correctness was recently done by Janka [24].
Convection

As just mentioned, the delayed mechanism is most likely not a completely satisfactory
explanation for supernova explosions as it has lead to failures in several (mostly one
dimensional) simulations that were performed since Wilson suggested it. That is why

the quest for a reliable explosion mechanism went on.

Given the more powerful computers that became available in the 19903, it was

 

10the innermost part of the core that is believed to form a neutron star later

19

an obvious approach to go to higher dimensions: several two dimensional simulations
have been carried out (see for example [23, 18, 22]) in which especially the impact
of convection has been studied. Convection is going to occur for sure in the region
between the shock and the neutrino sphere where the shock leaves behind an entropy
profile that is instable to convection. This convection takes place very rapidly: large
amounts of matter move up and down between the hot proto—neutron star and the
colder outer regions —- a little bit similar to the convection occurring when cold water
in a pot is put on a hot plate, yet on entirely different scales. A detailed discussion of
convection is clearly beyond the scope of this work. It shall be mentioned, however,
that it is widely believed today that convection helps the success of the explosion
by transporting energy from above the neutrino sphere (where the matter can be
heated by neutrino absorption in a relatively easy way) to the outer layers of the star.
Thus a larger fraction of the gravitational energy that was released during collapse is
made available for the explosion. Two dimensional simulations done by Herant et a1.
[23, 22] for example resulted in successful explosions while their one dimensional ones,

that did not include convection but used the same microphysics otherwise, failed.

Apart from the convection between the neutrino sphere and the stalled shock there
is a possibility that convection in the proto—neutron star may help the explosion by

intensifying the neutrino luminosities [24].

Still, this is not necessarily the final solution to the supernova problem. In partic-
ular, it is conceivable that phenomena that can only occur in three dimensions play

a role.

Rotation

It is well known that stars carry angular momentum. The impact of this rotation

on the core collapse and the explosion mechanism has been studied relatively little

20

(note however that studies including possible effects of rotation date back to 1970).
On the other hand, in the beginning there were ideas (now considered wrong) that
the explosion might actually be powered by the conversion of rotational energy via

magnetic fields [29].

As simplifications are almost always a necessity to model real phenomena, it was a
most reasonable approach to neglect rotation for the beginning. Another reason why
relatively little attention has been given to rotation so far is a lack of computer power
that is needed because it can only be simulated sensibly in more than one dimen-
sion. Simulating convection and neutrino transport is computationally so expensive
that simulations including these effects are still restricted to one and two dimensions
nowadays. There is also a lack of rotating progenitor models, so the precise initial
conditions for simulations of rotating core collapse are relatively uncertain. Some

studies were made, however, that will be presented in 2.3.

Indications have been found by several groups that rotation does not affect the ex-
plosion mechanism dramatically (see for example [18, 57]). It appears that explosions

are slightly delayed and a little bit weaker if rotation is included.

It must be mentioned, however, that often rotation is mimicked in questionable
ways necessary because most simulations are performed in two dimensions. Well
known phenomena from earth such as vortices in flowing liquids or tornadoes in
the atmosphere could hardly occur in less than three dimensions. We would also
like to point out that rotation can become extremely rapid in the late stages of core
collapse because the inner core contracts more and more. This enforces higher angular
velocities as its angular momentum is essentially conserved which can at smaller radii

and equal mass only be assured by faster rotation.

Another argument for studying the impact of rotation are recent observations

21

of the polarization of the light emitted by supernovae made by Wang and Wheeler
[50, 49, 51]. Using a method in which they observe the polarization and the wavelength
of the light simultaneously (“spectropolarimetry”), they found that the light emitted
by type II supernovae (and actually all core collapse supernovae which also include e. g.
types Ib and Ic) is significantly polarized while this is not (or not nearly as strong) the
case for supernovae not powered by core collapse. They draw the conclusion that there
must be large asymmetries in the explosion mechanism to cause this polarization.
They also found that the polarization of the light gets greater the closer to the center
of the explosion they measure it. This gives even more support to the conclusion that

the asymmetry originates in the very heart of the explosion.

Rotation is a very good candidate (and the only obvious one) to explain such
deviations from spherical symmetry as it defines a direction — that of the rotation

axis.

2.3 Recent Numerical Studies of Core Collapse Su-
pernovae

The separate treatment of previous numerical studies of supernovae is a bit artificial as
it was hopefully made clear in the previous sections that the currently predominating
theories for the explosion mechanism were chiefly developed by the excessive use of

numerical simulations.

This section will mainly be focused on the techniques, assumptions, and simplifi-
cations used in specific recent simulations (performed in the last ten years) while the
previous sections primarily dealt with their results. Here, this can hardly be done in
great detail for which one may refer to the respective citations. This is not meant to

be a complete overview of the work that was done in the last decade, we just want to

22

exemplify different approaches.

It shall be mentioned that there are also some attempts to deal with the supernova
problem on a more analytical level (not completely relying on numeric). The works

of Burrows and Goshy [9] or the very recent one of Janka [24] are examples for this.

2.3.1 Equation of State

An EOS describing the thermodynamic properties of the matter is a necessary input
for all supernova simulations. The purpose of an EOS is to make use of the fact that
the star is in (local) thermodynamic equilibrium. Thus only three thermodynamic
quantities (in core collapse simulations usually the temperature T, density p, and
electron fraction Y6) have to be calculated locally, all other quantities (the most
important among these the pressure and the specific internal energy) are given by the

EOS.

A great theoretical effort has been invested into finding a correct EOS, especially
the correct nuclear part of it (see for example [28, 42, 41, 43, 44]). The theory for this
is virtually a whole separate science (like several aforementioned aspects of supernova

theory). So we will again just give a brief account of the state of the art in this field.
Electron-Positron EOS

The electron-positron part of the EOS is less problematic than the nuclear part.
Its contribution is often just modeled by the EOS for a relativistic degenerate ideal
fermion gas, neglecting the Coulomb interaction, for which statistical mechanics pre-

dicts the temperature-independent pressure-density relation

pox.) = §(§)§hc(7:)§p%. (2.9)

 

23

where m3 is the baryon mass [59]. If an adiabatic change (662 = 0) is assumed (which
is reasonable at least for the infall phase) and the conversion of particle species is also
neglected (dN, = 0, where N, denotes the number of particles of species i), the first

law of thermodynamics reduces to

dU = 6Q — pdV + Z p,dN,~ .—_ —pdV (2.10)

and an expression for the internal energy density as? = Q5“ can easily be derived

from equation 2.9 by integration:

u:,:;><p,m=—j—(i)inc(”e )ipa. (2...)

 

However, this is a crude approximation. There is a variety of more sophisticated
EsOS for the electron-positron gas which cannot be expressed by such simple analytic
expressions as equations 2.9 and 2.11. They are available in the form of data tables or
computer programs. Five different ones of these are compared in [47]. Sophistication
is achieved by adding effects due to such phenomena like the presence of antiparticles
(i.e. positrons), the interaction of the electrons with the ionized nuclei present in the
matter, and radiation pressure prad. The latter is usually done using the well-known

blackbody approximation that yields for the radiation pressure:

4 0
prad : 5ZT4, (2.12)

where 0 denotes the Stefan-Boltzmann constant. The consideration of the former

effects is more complicated and will not be dealt with here.

In recent simulations, EsOS can be called over 109 times during a simulation run.
Hence, it is important for the applicability of these more sophisticated EsOS (more
precisely: the programs that calculate their output values or interpolate the tabulated

data) that they work reasonably fast. This is true for both the electron EOS and the

24

nuclear EOS. As far as the electron EOS is concerned the so—called Helmholtz EOS
developed by Timmes and Swesty [48] can be considered a great compromise between

execution speed, thermodynamic consistency, and accuracy.
Nuclear EOS

As mentioned in 2.2.3, the nuclear EOS (in particular at densities above the nuclear
matter saturation density) is of crucial importance for key features of the explosion

mechanism — especially for the vigorousness of core bounce.

Mainly for the purpose of saving CPU time, various extremely simplified analytical
EsOS have been proposed by different groups. All these are not rigorously derived
from theory but simply mimic essential key features such as the drastic rise in pressure
and specific internal energy above nuclear saturation density. For example, Zwerger
and Miiller [59], Yamada and Sato [57], and Bonazzola and Marck [5] use a simple
polytrope EOS (see equation 2.5) imitating the stiffening above nuclear matter density

with a large adiabatic index like e.g. 7 = 2.5.

A very similar example is the high density EOS introduced by Baron et al. [3]:

p(p) = [20:0 [(p—i)7 — 1], (2.13)

 

where K0 is the compression modulus, p0 the nuclear matter density, and 7 a param-

eter. All these parameters are isospin dependent in their model.

Just as in the case of the electron EOS much more sophisticated EsOS are avail-
able in the form of computer programs (e.g. [27, 40]). In their derivation one has to
deal with a lot of difliculties: at lower densities during core collapse the matter forms
essentially a two phase system — a “gas” phase of free nucleons and alpha particles
and a “condensed” phase consisting of heavier nuclei. At densities around one-half

p0 a so called “Swiss cheese” phase is formed: the nuclei touch each other uniformly

25

filling large space regions except for certain holes. Finally, at very high densities, the
nuclei are “squeezed” beyond their saturation density. Thus, the nuclear matter goes
through phase transitions which does not make things easier. These more sophis-
ticated EsOS also rely on certain parameters (such as the nuclear incompressibility
modulus or the nuclear surface tension) for which no exact values (for the densities

and temperatures of relevance here) could be obtained from experiments so far.

The grand masters of these EsOS are Lattimer and Swesty whose EOS (from now
on referred to as “LS EOS”) [27] has been used in many recent simulations (see e.g.
[18, 31, 34, 32, 23, 33, 7, 37] just to name a few) and can be considered state of the art.
It is valid only when the matter is in nuclear statistical equilibrium which is achieved
at the high temperatures in a supernova core. The assumptions and techniques used

in the derivation of the LS EOS are described in great detail in [28].

2.3.2 One Dimensional Simulations

One dimensional simulations in which spherical symmetry is assumed still play an
important role in the process of better understanding the supernova explosion mech-
anism. They are still the only way to simulate the underlying microphysical processes
in great detail. However, while one dimensional simulations became more and more
SOphisticated over the years, most of them failed to produce explosions. It thus
became clear that phenomena that can only appear in higher dimensions, such as
convection or rotation, are very likely to play an important role in the explosion

mechanism.

Mezzacappa, Liebendorfer et al. [34, 31]

A very recent example of such a “fizzle” are the state of the art one dimensional

simulations performed by Mezzacappa, Liebendbrfer et al. [34, 31]. They follow the

26

whole explosion process of a 13M® star using Newtonian [34] and general relativistic
[31] gravity, explicitly not including convection (that is sometimes effectively included
in one dimensional models) and paying special attention to the details of neutrino
transport in order to make sure that the failure of previous simulations is not due
to approximations in this area. They use a method called multigroup flux-limited
diffusion including effects due to different neutrino flavors and energies which can be
important because the cross sections for the interactions of neutrinos with the matter
strongly depend on these. The LS EOS (see 2.3.1) is used in their hydrodynamic
simulation that uses an adaptive grid. The latter means that the size of the cells for
which they calculate quantities like mass density, temperature, or neutrino density is
variable. This is used to obtain a preferably higher resolution in regions where one

or more of these quantities change rapidly by going to smaller cell sizes there.

As mentioned above, their main result is that they do not obtain an explosion

thus giving further support to the importance of higher dimensional effects.
Herant et al. [23]

Most of the recent work has been done in two dimensions. Quite often, however,
results from one dimensional simulations are used to simplify the two dimensional ones
by mapping data obtained from the former ones on the latter ones. One dimensional
simulations are also performed in order to compare their results to those achieved

with two dimensional ones.

An example for this is the simulation by Herant et al. [23]. They use a second
order Runge-Kutta hydrodynamic grid-based code with adaptive cell size without
including eflects of convection in their one dimensional model. The way they model
neutrino transport is strongly simplified: cells in their grid are either labeled optically

“thick” or “thin” depending on the density currently available. In the “thick” regions

27

neutrino transport is approximated by diffusion (see below) and in the “thin” regions
a so—called “central light bulb approximation” is used. This means that the matter
in these regions is bathed in a neutrino flux originating from the center of the star,
the magnitude of which is essentially determined by neutrino production rates in the
“thin” regions and the rate at which neutrinos diffuse from “thick” to “thin” regions.

They use the neutrino production rates calculated by Takahashi et al. [45].

For the diffusion they use a so—called “flux-limited” method. This means that

dn . c
a — m1n(§[Vn[, V - DVn), (2.14)

where n denotes the neutrino concentration in a cell and D is a diffusion coefficient.
This obviously sets an upper bound for the flux thus limiting it. If neutrino concen-
trations are very high, diffusion may be obstructed by their fermi properties which
disallow them to share identical quantum states. This effect is also included in their

model of neutrino transport.

Their one dimensional model does not explode.
Burrows et a1. [10]

Finally, the work of Burrows et al. [10] shall be mentioned. They also ran both one
and two dimensional simulations using monopole Newtonian gravity in a hydrody-
namic grid-based code in which (contrary to the previously mentioned works where
the grid is fixed to the mass elements) the cells are fixed in space, which forces them
to use a higher cell resolution the closer they get to the center in order to get an
appropriate resolution of core bounce and shockwave formation. Neutrino transport
is approximated by assuming diffusion. Three difl'erent neutrino species (Ve, 17., and
12,, where 11,, represents all remaining neutrino species) are treated separately. The

neutrino absorption and emission rates in the “gray” (semi translucent) regions are

28

approximated by relatively simple analytic expressions. The EOS they use is a slightly

modified version of the LS EOS.

The results of their simulations will be discussed in 2.3.3.

2.3.3 Two Dimensional Simulations

Fryer and Heger [18]

The work of Fryer and Heger [18] is one of the most interesting with respect to our
own simulation as they study the core collapse of the rotating progenitor calculated
by Heger [20] and try to find possible explanations for the polarization measurements

by Wang and Wheeler [50, 49, 51].

The angular velocity profile they use is shown in fig. 2.5. The discontinuities
of the angular velocity are located at the boundaries of well-defined regions in the
progenitor model in which angular velocity is equilibrated by convection. They use
a numerical technique called smooth particle hydrodynamics [19, 21, 22, 23, 12] an
essential feature of which is that the star’s mass is represented by discrete particles.
Neutrino transport is treated in the same way as by Herant et al. (see 2.3.2). A
patchwork of different EsOS for the different density regions is used, among these the
LS EOS for high densities. Gravity is treated using general relativity and assuming
spherical symmetry. Cylindrical symmetry around the rotation axis but no equatorial
symmetry is assumed for the rest (the former is a requirement to reduce the number
of dimensions from three to two). Angular momentum conservation is enforced for

each particle separately.

They find that rotation limits the convection overall and restricts it to the polar
regions. This ultimately delays the explosion and also weakens it. Further, asym-

metries due to rotation are apparent: the mean velocity of the matter in the polar

29

 

 
   

 

 

 

 

 

 

'5 I l I l I c; I I] I I I r I r I I l I fl
h 2 29 e -
«2 o z
10 .5” o .—i o. -:
: - - - - " u N :
- I‘\ “L _
- \ ~~ :
S .-
‘\
1 :— \\\ _:
5‘ : '0...\ j
.m : a»... “\ :
v 0.. \
c: . --------- MM89 Model A \ \ -
- «0.. \‘ \ _
--- MMBQ Model a s‘ \\
0.1 .— ...'. \\ \ —:
: —- MM89 Model C \ -
_ .‘0, \ q
.- '. \ -I
_ '0. \ -I
- o ‘ u
o. \
.0 \

.. '. \ '

0" ‘

.5. ‘
0.01 :- ‘i

: 1 L l l l 1 1 11 1 l l l 1 l l I l 1 .5.

107 10° 10°

Stellar Radius (cm)

Figure 2.5: Angular velocity profile used by Fryer and Heger [18] compared to that
of different models used by Monchmeyer and Miiller (labeled “MM89”) [35]

regions is larger by a factor of 2 compared to the equatorial region. They speculate

that these asymmetries might explain the polarization of the emitted light.
Herant et al. [23]

The two dimensional model by Herant et al. [23] is very similar to the one of Fryer

and Heger just described (the latter can actually be considered an advancement of the

30

former). Both use essentially the same numerical technique. Herant et al., however,
approximate gravity using the Newtonian limit assuming spherical symmetry. More
importantly, the effects of rotation are not studied in their work. They use a (non-

rotating) progenitor calculated by Woosley and Weaver.

Their main result is that convection above the proto—neutron star is the key to
obtaining an efficient and robust explosion mechanism. Their model turned out to be
insensitive to microphysical details like neutrino cross sections (for which they used
a variety of values). They compare the star in the phase after the shock has stalled
to a thermodynamic engine in which the convective matter transports energy from
the hot reservoir (the proto-neutron star) to the cold reservoir (the star’s envelope)
thus powering the explosion. Due to the large temperature difference between these
two reservoirs the efficiency of this “Carnot engine” is high. This leads to a self-
regulated explosion mechanism: the envelope is just heated enough to explode, which
is a possible explanation for the fact that all observed supernova explosion energies
are approximately between 1 and 2 x 1051erg. However, they admit that there are
certain subtleties in simulating convection in two dimensions and call the possibility

that rotation may be important “a very real one”.
Burrows et al. [10]

Most of what was said (e.g. the way neutrino transport is treated) in 2.3.2 about the
one dimensional model of Burrows et al. remains valid for their two dimensional one.
Actually, they also use the one dimensional simulation to simplify the two dimensional
one by mapping data from the former on the latter. For example, they just use the
information for the inner core obtained from the one dimensional model in the two

dimensional one.

Just like the two groups just mentioned, they also find that convection is helpful

31

to obtain a successful explosion. However, unlike Herant et al. they draw the conclu-
sion that the mechanism still significantly depends on the way neutrino transport is

modeled.
Yamada and Sato [57]

Yamada and Sato [57] follow a completely different approach than all the aforemen-
tioned works. They do not deal with neutrino transport, a realistic EOS, convection,
or other microphysical details but study the impact of rotation in a simple two di-
mensional model. Essentially, a polytrope EOS (equation 2.5) is used in which they
make the adiabatic index 7 a function of density: at low densities (p < 4.0 x 109*)
and at densities at which neutrino trapping occurs but the nuclear forces can still be
neglected (10123373 < ,0 < 2.8 x 10143313), 7 = g for the relativistic degenerate elec-
tron gas is used. At intermediate densities (4 X 10935;; < p < 10123,;fi) the pressure

deficit due to electron capture and photodisintegration is mimicked by a drop of the

5

adiabatic index below 3

(e.g. 7 = 1.3). Finally, at densities greater than nuclear
(p > 2.8 x 101431;?) the stiffening of the nuclear matter is modeled by 7 = 2.5. Their

EOS also contains an ideal-gas-like temperature dependence.

The main purpose of their model is to find deviations from spherical symmetry in
a two dimensional hydrodynamic simulation . They study the impact of rotation on
the collapse of the iron core by artificially adding angular velocity to a (non-rotating)

15M® progenitor calculated by Woosley and Weaver.

According to their results, the maximum density at core bounce, the explosion
energy, and the ejected mass all decrease monotonically with increasing angular mo-

mentum of the core.

32

Zwerger and Miiller [59]

The approach of Zwerger and Miiller [59] is very similar to that of Yamada and Sato as
far as the EOS and the treatment (or rather “non-treatment”) of neutrino transport is
concerned. However, as initial models Zwerger and Miiller use 7 = g - polytropes11 in
rotational equilibrium whose collapse is initiated by suddenly decreasing the adiabatic

index below 3‘:- (to values between 1.325 and 1.28).

Their main focus is on the gravitational wave signal of the supernova event, a topic
also studied by several other groups (e.g. [5, 58, 38]) that will not be discussed here
any further. An interesting result is that in some of their rapidly rotating models core

bounce occurs due to huge centrifugal forces even before nuclear density is reached.

This work must be seen in the context of a series of similar studies (among these
[38, 15, 25]) performed by different people at the Max Planck Institut fiir Astrophysik

in Garching, Germany.

2.3.4 Three Dimensional Simulations

Very few core collapse simulations have been done in three dimensions. The ones
we know of mainly deal with the gravitational wave signal [5, 38] emitted by the
explosion and use the simplified polytropic EOS described above (mimicking key fea-
tures by manipulating parameters in one way or another). The effects of rotation
are included in these models. To our best knowledge, there is currently no three
dimensional simulation of a core collapse supernova in which a realistic EOS, a real-
istic progenitor model, or even neutrino transport (to say nothing of convection) is
included. This situation is naturally due to the immense computational requirements

of such a simulation.

 

11equilibrium solutions of the polytrope EOS 2.5

33

It shall be mentioned, however, that Chris Fryer is currently working on such a

model [16, 17] but has not published any results yet.

34

Chapter 3

Our Three Dimensional Simulation

As pointed out in 2.3.4 there were hardly any core collapse simulations in three
dimensions so far. The basic idea of our own simulation is to simplify the underlying
microphysics, mainly by using input from one and two dimensional simulations, in
order to be able to follow the core’s dynamics in three dimensions during collapse
and bounce. The main focus is put on the impact of rotation during collapse. Our
hope was to thus find deviations from spherical symmetry that are so significant that
they may deliver alternatives to the currently favored complicated convection-driven
explosion mechanism. After all, it is well known that much weaker pressure gradients
(than those present during core collapse) and slower rotation can lead to very major
phenomena like tornadoes in the atmosphere or vortices in flowing liquids (yet on a

much larger time scale).

We think that the effects of rotation can truly only be studied in three dimensions
because adding rotation in a two dimensional model is always somewhat artificial. A
three dimensional approach is further motivated by the aforementioned polarization
observations by Wang and Wheeler [50, 49, 51] for which no satisfactory explanation
could be obtained from two dimensional simulations so far. The dependence of many

one and two dimensional models on such (uncertain) parameters as the viscosity of

35

the matter and neutrino cross sections also raised our hope that there might be a

solution to the supernova problem that is largely independent of these details.

3.1 Test Particle Method

The method we use is similar to the so called test particle method or pseudoparticle
method that has been used extensively in nuclear physics [54]. Its original purpose
is to approximate the time-dependent Hartree-Fock equations in a way that allows

classical interpretation.

While in nuclear physics usually several thousand test particles represent one
nucleon, things are vice versa in our model for the collapse of an iron core with
a mass around 1.3MQ: (assuming that we use 106 test particles) one test particle
represents a mass of 1.3 x 10’6MQ z 2.5857 x 1024kg — just a little less than the mass

of the entire earth!

In our model all th particles have the same mass MIC/Nu, 2: mtp, where M [C
denotes the mass of the whole iron core. For each individual particle, position 7‘}
and momentum p‘, are tracked (as classical three vectors). The equations of motion
for the test particles are the relativistic versions of the Newtonian ones known from

classical mechanics:

 

. d". -'.
Ti; : -d—T£J— : p] .. (3'1)
mt. + (’35)2
:4 dfi. -o -o -o _.
Pj = % " F0001, MN...) + FEOSJ(TJ) (3 2)

j : 13"‘3th’

where Eag- denotes the force on particle j due to gravity and fieosg the force due

to the equation of state. Gravity is modeled assuming spherical symmetry by using

36

the Newtonian monopole approximation. This means that FIGU- is the gravitational
force test particle j would experience if the mass of all test particles located closer
to the center of the iron core (which is identical to the origin of the coordinates) was

combined in a point mass located at the center:

- mg, #{ie {1,...,N,,,} : [m < [73]}

Gt:— Im3

 

77,-. (3.3)

This approximation is only good as long as the deviations from spherical symmetry
(of the mass distribution) are sufficiently small. The calculation of FEOSJ is a bit

more complicated and will be explained in 3.3.4.

3.2 Grid

In order to be able to locally define thermodynamic quantities, most notably the local

mass density, we introduced a spherical coordinate grid.

3.2.1 Cells and Boundaries

The boundaries between the cells of this grid are defined by the surfaces of constant
r, constant d), and constant 0 in spherical coordinates (r, (19,0) using the standard
notation (r = [is], (b = azimuth angle, 6 = polar angle). For the 45 coordinate the

boundaries are located at

27r 27r 27r
0 1 — 2 —,..., N —1 x ——

where N¢ denotes the total number of boundaries (for the a) coordinate). For the

0 coordinate the boundaries are chosen so that the difference between cosd of two

arbitrary neighboring boundaries is constant, i.e.

[cos 01 — cos 02[ = const.

37

if 61, 62 are neighboring boundaries (for the 0 coordinate). This is made so to make
the volume of a grid cell independent of 9. Finally, for the r coordinate the boundaries

are usually chosen equidistant so there are boundaries at

R R
— 2 —... —=
ler, er, ,erNr R,

where R denotes the radial location of the edge of the spherical grid and N, the
number of boundaries for the r coordinate. Note that this way the volume of the cells
depends only on r (not on d or 0) and increases with increasing distance from the

center.
Radial Cell Width Feature

However, an option which enables us to use non-equidistant r-boundaries is included
in our code: a monotonic growing function f (y) : [0 : 1] —> [0 : 1], y I——) f (y) satisfying
f (0) = 0 and f (1) = 1 can be defined. The boundaries are then located wherever
(f (fi) x N.) E N. For example, f (y) = 3109 yields a higher grid resolution near the

center. f (y) = id(y) = y gives equidistant r-boundaries as before.

3.2.2 Grid Coordinates

Recapitulating, we have a spherical grid consisting of Nee”, :2 N, x N4, x Ncos0 cells.
These are conveniently labeled by three integers n, E {1, . . . , N,}, n¢ E {1, . . . , N¢},
ncosg E {1,...,Ncosg} where increasing nx corresponds to increasing X for X E
{r, (1), cos 6}1. These three integers (nr, n¢, ncosg) will be referred to as grid coordinates

from now on.

Given a point inside the grid in spherical coordinates (r, (15, 9), we can find its grid

 

1In the C++ code, these integers run from 0 to N x — 1.

38

coordinates using the following trivial relations:

 

FT

r z —N,.] 1 3.4

n .R + ( )
"(b

= —N] 1 .
71¢ 2” ¢ + (3 5)
r 0 1
”cost? COS 2+ NeosO] + 1: (36)

where the squared brackets [Y] denote truncation of the quantity Y. Note that

increasing ncosg corresponds to increasing cos 6 but to decreasing 6.
Cell Volume

Using this, we can now explicitly give the volume of a grid cell as a function of its

grid coordinates:

Vol(nr,n¢,ncosg) = V0102.) = fl {f-1(fl)3 — f-1(3’;1)3}, (3.7)

 

3N¢Ncosg N, N,
or for the most important case f(y) = f‘1(y) = id(y) = y:
Vol(nr) = Egg; {(1%)3 — ("N— 1)3}. (3.8)
Innermost Cell
The volume of the cells with n, = 1 is usually very small. Numerical problems

also arise from the “wedge shape” of these cells. These problems will become more
transparent in 3.2.3 and 3.2.4. Anyway, all these cells at the center are treated as one

(spherical) cell.
Grid Scaling

During a core collapse simulation an enormous change in the length scale occurs: the

inner iron core contracts roughly by a factor of 102 to 103. Even if N, is chosen very

39

large2 (e.g. N, = 500) the inner core at core bounce would still be almost completely
in the innermost cell if the grid were fixed in space (as long as f (y) = id(y) is chosen).
This way an appropriate resolution of core bounce and other phenomena would be
destroyed. That is why we chose to scale down the grid simultaneously with the core.
More precisely the size of the grid (which is given by R, the distance of its outer
edge from the center) is modified during collapse as follows: in every time step, the
distance of the %gth outermost test particle is determined and multiplied by a factor
slightly larger than 1 (e.g. 1.12). The advantage of this apparently odd procedure
is that it proved to be very functional: almost all test particles remain in the grid
during collapse and fill a reasonable portion of it, and if a few particles are “lost”
(what this means more precisely and how it can happen will be illuminated in 3.4)

the simulation is not ruined (which can be the case if just the outermost test particle

is used to define the size of the grid).

3.2.3 Calculation of Densities

A mass density p(n,, n¢, ncosg) can be calculated for each cell of the grid using the
following evident way. The number of test particles in the cell th(n,,n¢,, ncosg) is
counted, multiplied by the mass of a test particle mtp, and divided by the cell volume
Vol (n,):

N (nryn 7nC080) m
p(nran¢ancos(9) : tp V:l(n ) tp. (39)

 

In order to minimize errors due to fluctuations and to smooth the density distribution,
we decided to use a slightly more sophisticated method in which the test particles’
mass is smeared over the cell it is in and the seven (well-defined) neighboring cells
which are located nearest to the test particle. Figure 3.1 illustrates this method for

the two dimensional case: first, the cells which are closest to the test particle are

 

2Note that the total number of cells is clearly limited by the condition Nee", << NW.

40

 

Figure 3.1: Smearing of test particles for density calculation

determined by detecting in which of the sub cells I.a through I.d it is located. Then

the number fractions X of the test particle these cells “get” are calculated as follows:

 

 

 

x: = (1-%.;Z’f.,)(1-%f;¥;%) (3.10)
xu —-— (l-éfifléfrfil (311)
= <;.;:T:.><:.4:Is;) «12>
= 6.17;)(1-éif’ti-ile «m»

where r and 43 are the coordinates of the test particle and r3, rc, (b3 and dc denote the
coordinates of the boundaries of cell I.a (with r3 > re and d3 > ¢c)- 21:, x.- = 1 is

ensured. It is obvious how this generalizes to three dimensional spherical coordinates.

The method yields for example that a test particle that is located exactly in

41

the corner of a cell is uniformly shared among the eight adjacent cells (in the three
dimensional case), and belongs completely to the cell it is located in if it is exactly in

its center.

The procedure may seem artificial and unjustified because if you visualize the
smearing as being caused by giving the test particles a finite size, the technique just
described implies that the size of the test particles depends on the distance from the
center (as the cell size depends on it). But the use of the grid is only sensible if the
grid resolution can be chosen high enough for the variations of the density (and other
quantities defined for each cell) in neighboring cells to be small. Thus the smearing
is just a technique to smooth possible fluctuations the details of which should not

significantly affect our results.

3.2.4 Calculation of Derivatives

To follow the dynamics of core collapse, the calculation of (spacial) derivatives of
thermodynamic quantities is necessary. More precisely, this is needed to calculate

fieosy from equation 3.2 which will be explained in 3.3.4.

Let Q be a thermodynamic quantity defined on the grid (meaning that for each
cell (n,, n¢, ncosg) a real number Q(n,, n¢, ncosg) E R is defined). To approximate the
gradient VQ(r, (b, 0), we use a modification of the standard technique for calculating
numerical derivatives described in [36] and the well known expression for the gradient

in spherical coordinates [6]:

8f)
VQ(r, 45,0) = 5; (Mme, +

100 1652 _,

— " — ' 6 —— .14
rsin06¢ (r,¢,0)e¢+( sm ) r8(cos 6) (r,¢,9)69’ (3 )

 

where 5,, é}, é}; denote the orthonormal basis vectors for spherical coordinates. These

are dependent on the coordinates (r, (b, 0). What remains to be done is the numerical

42

definition of
(99
, an .
(we) 0(608 0) (me)

if}
Br

69
may 8(1)

 

 

 

Linear Interpolation of Derivatives

We will exemplarily describe our technique for this for 96%. Let (r, (b, 6) be located in

the cell with grid coordinates (n,, n45, ncosg). Then two obvious approximations for

 

 

 

 

(Ii—(1? are
80 Q( r ’ ’ C05 )—9( r: I cos )
(5; (r,¢,6))right : ;x+{lf:¢(:%l:) _ f_7:(:):;r;) }9 (3.15)
89 Q( T: acos)_Q(r" a 3 cos)
(I)? (r,¢,6))left = gxnff71(:rN——ri) ff_11(:%%n)}0 a (3.16)

where the denominators are just the distances in the r coordinate between the centers
of the neighboring cells (2 NA: for f (y) = id(y) ). We decided to linearly interpolate

the derivative between these two values:

 

 

09 _ r—rs X(6f2

51: rB—rx(8§2

E

 

(3.17)

 

87" (rm?) — TB — 7‘s (r,¢,9))ri9ht TB — 7’3 (man) left,

where rs and r3 denote the r-coordinates of the cell boundaries of cell (n,, n¢, ncosg)

with r5 3 r 3 r3 (just as in figure 3.1). 3::- and

80
6 cos 6

 

are in principle interpolated

 

the same way. For 6320 note that

an _ —1 99
6cos6 (r,¢>,6) ‘ sin6 66 (gas)

 

which was already used in equation 3.14.

Problems at Grid Boundaries

At the boundaries of the grid and the z-axis, boundary conditions for the derivatives

have to be fixed. To a certain extent these are always arbitrary and physical reasoning

43

is needed. Other problems may arise from the fact that test particles can be located
outside the grid. We will explain in 3.4 (when the physical meaning of the spatial

derivatives will have become clear) how these problems are dealt with.

3.3 Equation of State

Once we have defined the needed thermodynamic quantities on the grid, an appro-

priate EOS can be used to obtain remaining ones.

3.3.1 Cold Nuclear EOS and Polytrope EOS

As a first approximation to a realistic EOS for supernova matter we used the following
well-known parameterization for the energy per baryon of cold (isospin symmetric)

nuclear matter [30] which only depends on the density:

p p "
unuc = a. — + b — 3.18
(p) m (p0) ( >

with pa = 2.4 x 10143517 2 nuclear matter density, a = —218MeV, b = 164MeV, and
0 = g. This EOS yields such realistic features as a minimum of the energy per baryon
at p0, a slight attraction for the matter at densities below p0, and a strong repulsion

at densities higher than p0.

Once an electron fraction is defined (see 3.3.2), the presence of electrons can be
modeled by adding equation 2.11 (multiplied by mp2 to get the internal energy per
baryon) to equation 3.18:

p p 0 9 7r § 1 p %
umt(p, Y6) = a Ed + b (3;) + 1(3) the3 (577;) . (3.19)
Note that an electron fraction other than 0.5 implies that an asymmetry energy term
should be added to the nuclear part. As the correct shape of this term does not stand

firm yet [30] and the isospin asymmetry of the matter remains relatively small (as

44

58-11 . . -....., . .. ..., . . ......, 4 -
TotalEOS—
NuclearPart -------
ElectronPart --------

 

4e-11

36-1 1

29-11

Energy per Baryon [J]

1e-11

 

 

 

 

 

Mass Density in Units of Nuclear Matter Density

Figure 3.2: Cold nuclear EOS and polytrope EOS for electron gas

described in 2.2.1), it is neglected here. Figure 3.2 shows the different contributions

to this EOS (using an electron fraction as defined in 3.3.2).

Equation 3.19 is obviously a crude approximation because (among other things)

temperature does not even appear.

3.3.2 Calculation of Thermodynamic Quantities other than
Density

More sophisticated EsOS generally need more input than just the density of the
matter. The ones that will be used in this work will be described in 3.3.3 and need
the temperature T, electron fraction Y8, and the composition of the matter (the latter
is given by the average charge and mass numbers of the present nuclei). A realistic

calculation of these quantities is not nearly as simple as that of the density: for the

45

electron fraction, for example, one has to use proper electron capture rates (such as
the ones calculated in [45]) and once the produced neutrinos are trapped, neutrino

diffusion has to be dealt with.

Knowing that in a three dimensional simulation including all these detailed calcu-
lations is not yet feasible, we decided to circumvent these difficulties by mapping data
from one dimensional simulations on our model. More precisely, the data calculated
by Cooperstein and Wambach [14, 13] was used. Their tabulated data enables us to
define the temperature T(p) and electron fraction Y8(p) as a function of density and

should yield a good approximation for the state of the matter until core bounce3.

Temperature

The tabulated data obtained from Cooperstein’s work was interpolated using appro-
priate fit functions. The relation T(p) that was finally used in our simulations is

shown in figure 3.3 together with Cooperstein’s values.
Electron Fraction

The electron fraction as a function of density Y8(p) was obtained from Cooperstein’s
data and interpolated just like the temperature. The resulting relation for our sim-
ulation is shown in figure 3.4. Note that Ye(p) is continued relatively arbitrarily for
p > 10125353, yet the convergence of Ye to a value Ye z 0.31 for very high densities is

realistic as described in 2.2.1.

Further Thermodynamic Quantities

More thermodynamic quantities such as the entropy per baryon or the electron chem-

ical potential could be defined using Cooperstein’s data but this is not necessary

 

3Our simulation can currently only be considered realistic until core bounce as will be further
explained in 3.7.

46

 

 

 

 

 

6 I- — —-4
5 .— -_
5'
O
a. 4 » -
2
3
‘5
8
E 3 u. m ..
0
'—
2 ~ ~ ~-
. Cooperstein Data1 +
Cooperstein Data 2 X
Used in Simulation
0 . 1 1 1 . . . 1 L . . l . . . l . . A
0.1 1 10 100 1000 10000

Density [10“ g/cm3]

Figure 3.3: T(p) for the infall phase

because the E808 we use do not need more input. Au contraire — virtually any

imaginable quantity of interest can be returned by the E508.

3.3.3 Helmholtz EOS and Lattimer & Swesty EOS

As a more realistic EOS (than the one presented in 3.3.1) we used a combination of the
nuclear EOS by Lattimer & Swesty and the Helmholtz EOS by T immes (mentioned
in 2.3.1). The former is used for p Z 1011:53’ the latter for p < 10110—55 where the
nuclear contribution to the pressure is negligible. Note that the LS EOS also includes
an electron gas contribution. Moreover, there are different parameter sets for the LS
EOS which mainly affect the behavior of the matter at densities above nuclear. We

used the parameter set in which the nuclear compressibility is K = 180 MeV.

Both 13508 are available as FORTRAN-programs [27, 46]. The input for the LS EOS

47

 

0.44 V V V I —r V V r Y f 1 I V —V_ ‘ Y r. I .
Used In Simulation
Cooperstein Data +

 

0.42 r

0.4 --

0.38 -

0.36

Electron Fraction

0.34

0.32

 

 

 

 

0.3 A l A l A A 4 l A I A l A L A l A_ L

0.1 1 10 100 ~ 1000 10000
Density [10" g/cm3]

 

Figure 3.4: Y6(p) for the infall phase

is p, T, and Y... The Helmholtz EOS gets the composition of the matter instead of
Y... This composition is defined by the mass and charge numbers of an arbitrary
number of different ions present in the matter and their respective number fractions.
We used a small constant hydrogen fraction of 2.69 x 10'4 (taken from [14]). The
rest of the matter got a charge number Z = 26 and a mass number A(p) = Z / Y6(p)
mimicking the change of the matter composition due to ongoing electron capture
while taking into account that most of the nuclei present are iron isotopes. Note that
this approximation does not lead to absurd isotOpes (with extreme neutron excess)
because the Helmholtz EOS is only used for densities at which Y; z 0.41 which yields
A g 63.

The internal energy per baryon uint(p) = um, (p, T (p), Y8(p)) this ultimately yields

using T(p) and Ye(p) from 3.3.2 is shown in figure 3.5. There is no longer a density

48

 

160 V I 7" V U I] V I 'l V V 7" V V II V Y " V V 7'

140 l-- _,__, - --—~—

120 -

100 - ~ --

80 -. ...._._

so - - -1-

lntemal Energy per Baryon[MeV]

20b

 

 

 

L Al A A Al I A l A A A1 A I ‘1

0.0001 0.001 0.01 0.1 1
Density [Pol

 

0 A - 1 x 1 1
1 9-07 1 9-06 1 9-05

Figure 3.5: Combination of Helmholtz EOS and LS EOS used in our simulation

region in which the EOS is attractive.

3.3.4 How the EOS Affects the Dynamics

One major question still needs to be answered: how does the EOS affect the dynamics
of our core collapse simulation? This happens through F303,], the force due to the

EOS on test particle j as introduced in equation 3.2.

The EOS is used to calculate the internal energy per baryon um for every grid cell.
This is easily converted into an internal energy per test particle by multiplying with
mtp/mB =: V, the number of baryons per test particle. A gradient of this quantity
at the location of test particle j is calculated using the technique described in 3.2.4.
Then

F - = -l/ Vuin 3.20
E03,] t (T,¢,9)j ( )

49

where (r, (19,0)j are the spherical coordinates of test particle j’s position.

This technique is justified by energy conservation: during the infall phase it is a
good approximation to neglect energy losses due to neutrinos (and photons) radiating
away from the core because the magnitude of these losses is relatively small. Thus the
sum of the core’s kinetic energy Ekm, gravitational energy E0, and thermal (internal)
energy Em

Etot = Ekin + Ea + Em: (3.21)

should be approximately constant during collapse. After bounce this is no longer

valid as neutrino losses cannot be neglected anymore.

We will now explain how our method of calculating F1303 ,3- implies the (approxi-
mate) conservation of E“. For simplicity, assume the situation shown in figure 3.6:
two neighboring cells (n,, 7145.71.20.29) and (n, + 1, 72¢, ncosg) with the respective internal
energies per test particle ul and 712. Moreover let uz > ul, a test particle j located in
cell 1, and the internal energy per test particle in all remaining neighboring cells of
cell 1 be equal to ul so that FEOSJ only gets a contribution from the gradient between
u} and 11.2. Let gravity be “turned off” for the moment. Then, the force on j due to

the EOS is approximately

U2*U1..

6.-
Ar

 

FEOS,j '2 —Vu z —

where Ar is the radial cell width. Now assume that j travels from cell 1 to cell 2 as
indicated by the dashed arrow. In doing so, due to the action of 13” 303,33 its kinetic

energy is reduced by

112—111

Ar

AT 2 11.1—11.2.

 

2
AEkz'n =/ FEOS,j ° dl‘ z -
1

Thus, the internal energy j gains (it is now located in a cell with a higher internal

energy per test particle) is (approximately) compensated by its loss of kinetic energy.

50

 

Figure 3.6: Test particle j traveling from one grid cell to another

In this derivation, the way gradients are calculated was simplified and the fact that
the appearance of j in cell 2 (and its disappearance in cell 1) affects ul and U2 was

neglected.

One might consider using the pressure as returned by the EOS used (instead of
the internal energy) to calculate FEOSJ. However, there is no evident way of doing
this in our model without making further assumptions about the size and shape of
the test particles. We wanted to avoid this because the test particles are after all

completely imaginary objects.

In our simulations, energy (and angular momentum) is kept track of to assure that
the method just described works properly and the conservation of these quantities
is not destroyed by numerical limitations, a too large time step size, a too low grid

resolution, or other effects.

51

3.4 Symmetry Assumptions, Boundary Conditions,
and Numerical Problems

Like in almost all numerical simulations there are certain subtleties in ours that
require special attention. It is also desirable to make simplifications wherever possible

to invest the available CPU time where it does most good.

3.4.1 Symmetry Assumptions

The total number of grid cells News is limited by the condition Nceu, << th because
otherwise there are inevitably many cells which contain very few, just one, or no test
particles. This leads to unphysical density fluctuations on a length scale essentially
given by the cell size. Most of the time we used th = 106 and found that this kind of
trouble occurs when News 3 2 x 104 is used. To get a good resolution of interesting
features, it is thus desirable to make use of symmetries. Two symmetries can be

anticipated almost a priori:

o equatorial symmetry: there is virtually no imaginable reason (other than fluc-
tuations or numerical errors) why equatorial symmetry should be broken in our
model. Nevertheless, we ran several simulations in which it was not enforced
without finding significant deviations therefrom. Therefore equatorial symme-
try is enforced by averaging the thermodynamic quantities in grid cells mirror

symmetric about the equatorial plane.

0 cylindrical symmetry: exactly the same way, it was found that there are no
significant deviations from cylindrical symmetry (about the rotation axis) in

our model. Hence this symmetry is enforced by setting N¢ = 1.

52

Note that no symmetry whatsoever was enforced for the positions and momenta of
the test particles (i.e. these are not located in positions mirror symmetric about the

equatorial plane or something of the kind).

3.4.2 Derivatives at Grid Boundaries

z-Axis

The %2— derivative at the z-axis is simply set to zero. This causes the test particles
initially located near the z-axis to tend to “stick” to it. This is not good and only
done this way because we could not think of a better way to deal with it. Note,
however, that this affects just a small fraction of the test particles. Those particles
near the z-axis are initially rotating very slowly around this axis so that there is at

least no obvious reason why they would move away from it. This treatment can also

be seen as a consequence of assuming cylindrical symmetry.
Grid Surface

At the radial edge of the grid, % is calculated by assuming that the density of
the matter outside the grid has a value pm,” slightly lower than that present at the
outermost layer of the iron core at the onset of the simulation. Thus the presence of

the star’s mantle and envelope beyond the surface of the iron core is imitated in a

simple way.
Center

The derivative 9%! at the center is also set to zero. This is reasonable because matter

should not traverse the core’s center in large amounts.

53

3.4.3 Background Density

As described in 3.2.2 the grid is usually chosen a bit larger than the space region
in which the vast majority of the test particles is located. Hence, empty cells may
occur. To assume p = 0 in these would be completely unrealistic as even outside the
iron core (in the star’s mantle) densities are high. Thus a minimum density pmz-n (the
magnitude of which is determined as described in 3.4.2) is enforced throughout the

grid by setting p to pm,-n wherever p < pmm.
3.4.4 Singularity Treatment

The force on particle j due to gravity 130,]- calculated as in equation 3.3 has a (nu-
merical) singularity at IFJI = 0. This can cause numerical trouble for very small |1"'j|
which yield extremely large forces. Thus these particles can be vigorously accelerated
and more or less “shot” out of the core’s center. This is unrealistic and not desirable.
Therefore we decided to “switch off” gravity for the 11551 or 110‘; innermost particles.

Note that IFGJ| is relatively small for the innermost particles anyway because the

mass they enclose is comparatively small.

3.5 Time Development

The time development of the collapsing core is essentially obtained by numerically
integrating the system of 2 x 3 x th coupled first-order ordinary differential equations
3.1 and 3.2. This is done using a fourth-order Runge-Kutta algorithm, a standard
method precisely described e.g. in [36]. The whole simulation program (except for

the E808 mentioned in 3.3.3 which are written in FORTRAN) is written in OH.

This algorithm as such uses a constant time step size. As it is clear that the core’s

matter will be continuously accelerated during collapse we implemented a very simple

54

way to adapt the step size nonetheless: once the core has contracted by a certain
factor (e.g. 10) the step size is suddenly made sufficiently smaller (e.g. divided by
10) to guarantee a good resolution of the further time development. Depending on
the initial conditions used in the simulation this step size modification can be done

several times.

3.6 Calculation of Observables

Several physical observables are kept track of during the time development of the
collapsing core. Among these are the total kinetic energy Ekin, (Newtonian) gravi-
tational energy EC, internal energy Em, and the total angular momentum E of the
core. It is important to do this in order to make sure that the core’s total energy and
angular momentum are (at least approximately) conserved. It is also instructive to
follow the conversion of the different energy types into one another. The observables

are calculated as follows:

0 E0 is calculated assuming spherical symmetry:

MP ' 1,.. .,N ,
E0 2 -ZG#{Z 6{ ..}]Ir1<lr.l}m.. (3.22)

lrjl

 

i=1

0 For Ekm the well-known expression for the relativistic energy is used:

 

 

IV”, an9
Ek,n = Z(\/m,2pc4 + p'3-2c2 — mtpcz) = Z \/m,2pc4 + 153-2c2 — thmtp02 (3.23)
i=1

i=1

0 Em is evidently calculated by summing up the internal energies of all particles:

Ncosa

Nr N¢
Eva-1: Z 2 ”E15“ (nr7n¢3nC080) X th(nr1n¢)nC089)) (3-24)

nr=l 71¢: l "€080: —1

where u] m) now denotes the internal energy per test particle not per baryon.

55

o The calculation of f: is also obvious:

3
'3

11' = (3.25)

H.
II
p—o
5.31
X
k l

3.7 Advantages and Weaknesses of this Method

A big advantage of the method described in the preceding sections is its simplicity: it
enables us to simulate the collapse of a rotating iron core in three space dimensions
until core bounce on a quite ordinary home computer. Our code also conserves an-
gular momentum (apart from tiny numerical errors) which was verified in countless
different simulation runs using different initial conditions and E808. Angular momen-
tum conservation is even good if cylindrical symmetry is not assumed. If the step size
is chosen sufficiently small, energy conservation is given at least approximately (this

will be deepened in 4.2.1, 4.2.2, and 4.2.3).

Yet, one has to keep in mind that it is a model. There is no way whatsoever to
guarantee that the results our simulation yields are completely realistic (this is more or
less the case for all core collapse simulations). Possible new phenomena found in such
models may however guide our intuition in understanding the supernova explosion

mechanism.

A weakness of our model at this time is that there are no collisions between the
test particles. This is the main reason why our simulation is currently limited to
the infall phase: after bounce the test particles that rebound from the core just fly
through the infalling outer core matter without resistance. No shockwave is created
making this a completely unrealistic scenario. Even during the infall phase some

trouble arises from the lack of collisions (that will be described in chapter 4).

The simplicity of the model mentioned above can also be considered a weakness:

56

the use of data obtained from lower dimensional simulations certainly restricts the

possibilities for the occurrence of new phenomena.

3.8 Computational Requirements

The computationally most expensive operations will be depicted now. In order to
calculate FGJ the particles essentially have to be sorted by their distance from the
center. Moreover, to calculate FEOSJ the EOS must be evaluated for every cell of the
grid. Both things need to be done four times in every time step (due to the way the

fourth-order Runge-Kutta algorithm works).

A typical simulation run uses 106 test particles, approximately 12 x 103 grid cells,
and about 103 time steps. Therefore the quicksort algorithm as described in [36]
(which is used to sort the particles) is called 4000 times. Most expensive, however,
are the EOS calls: 4 x 12 x 103 x 103 z 5 x 107 of these occur in a typical run. The
time needed to calculate one time step4 can become up to z 4min when the LS EOS

is used, leading to a total time of roughly 4000min z 3days.

3.9 Output Possibilities

Two different output programs were written to visualize the data created by the
simulation program. For simplicity, Microsoft Visual Basic© was used for these
programs which involve a lot of graphical output. The source code for both programs

is reproduced in appendix B.

 

4All these times are valid for a 1GHz Intel Pentium© III with 512 MB RAM that was used
throughout this project.

57

3.9.1 Density Output

The mass density in a slice through the core that includes the rotational axis (the
z-axis) can be plotted for every time step. The different densities are indicated by

colors. Due to the large density changes a logarithmic density scale had to be used.

3.9.2 Test Particle Output

Another program shows (up to) 2000 test particles in a (pseudo) three dimensional

5. The motion of these individual particles can be followed during collapse.

picture
The perspective can be changed at any time. Apart from giving a good impression of
the collapse dynamics, this program is very useful to identify errors like those arising

from a too large number of grid cells.

 

5Showing more particles is not sensible because separate test particles could than hardly be
discerned.

58

Chapter 4

Numerical Results and
Interpretation

We will now describe some results obtained using the code described in chapter 3.

4.1 Angular Momentum Conservation

In all simulation runs angular momentum conservation was almost perfect. This
should not come as a surprise because cylindrical symmetry about the rotation axis
was assumed. The plots of the total angular momentum IE] as a function of time are
often just a horizontal line. A typical angular momentum evolution (for the initial
conditions and EOS that will be described in 4.2.1) is shown in figure 4.1. Note that

during the fluctuations at the end of the simulation run IEI varies by less than 0.01%.

4.2 Results of Simulation Runs

We performed three series of simulations in which different EsOS were used.

59

 

2.604088+42 j I I I I I I r
2.604066-1-42 - 1 -- .. .- . Q , -4.
2504046442 - —- » ~77 . - . .. _ _ , _,_

2.604026+42 r- ~~ - - L- 7 - .- , .,, . ”-7 -_. - - -...

Total Angular Momentum [Js]

 

 

 

 

2.604e+42 . , . , - - . . . . . -
2.60398e+42 ~ . . . . . _
+
2.60396e+42 ~ f -- -- . .. _- . .1.
it
2.60394e+42 ‘ ‘1 , I .
l l l J L l l l
0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18
Time[s]

Figure 4.1: Angular momentum as a function of time for a typical simulation run

4.2.1 Cold Nuclear EOS without Electron Contribution

The first series is more of a toy model we created as a first test for the code: we only
used the cold nuclear contribution (as in equation 3.18) to the EOS. Nevertheless this
model exhibits some realistic features. The iron core mass was chosen to be 1.5MQ
and the background density set to pmin = 0. To achieve a homologous collapse, a
spherical homogeneous mass distribution had to be chosen: the test particles were
randomly distributed in a sphere with radius 1.676307x106m which is about 100 times
the radius this mass would have at nuclear matter density (if also homogeneously
distributed in a sphere). This yields an (approximately) homologous collapse because
the cold nuclear EOS hardly affects the matter at all at densities p < p0 and the

collapse of a homogeneous self-gravitating sphere is homologous.

60

 

 

 

Name of Run wo [Lg—d] IZOI [1042.18] giggling tbounce [ms] pm“ [p0]
CNEOSOO 0.0 0.0 0.0% 173.8 9.5
CNEOSO3 0.3 0.977 0.073% 174.1 3.7
CNEOSO4 0.4 1.30 0.13% 174.2 5.1
CNEOSO5 0.5 1.63 0.20% 174.5 3.5
CNEOSOG 0.6 1.95 0.29% 174.8 3.6
CNEOSO7 0.7 2.28 0.40% 175.0 3.0
CNEOSO8 0.8 2.60 0.52% 175.4 2.8

 

 

 

 

 

 

 

 

 

Table 4.1: Simulation runs using the cold nuclear EOS. The abbreviations are ex-
plained in the text.

In this series the cores were started rotating like a rigid body around the z-axis.
Different initial angular velocities 020 were used. The grid parameters N, = 120,
N4, = 1, and Ncosfl = 100 were applied. 106 test particles were employed. All runs of

the series are summarized in table 4.1.
Energy Conservation

Energy conservation is only fulfilled approximately. It is very good before the EOS
starts to strongly affect the dynamics at high densities. Figure 4.2 shows a typical
energy development. The time on the abscissa is measured from the onset of the
simulation. At earlier times (than shown in figure 4.2), only Eh," and EG contribute
significantly to Eat. Although the total energy is obviously not conserved at all
times it is good to see that the energy deficit that arose shortly before core bounce is

ultimately corrected.

Time of Core Bounce and Maximum Density

Note that a time them“ at which core bounce occurs can be identified in the energy
graph as a minimum of E0. Another common observable in supernova simulations

is the maximum density pm” achieved at core bounce. We decided to calculate pm“

61

 

3e+46

2e+46

1e+46

Energy [J]
O

-1e+46

-2e+46

 

 

 

 

-3e+46
01734 01736 01738 0174 01742 01744 01746 01748 0175

Time [s]

Figure 4.2: Time development of the different energies in the late stages of simulation
CNEOSO5

by averaging over a certain number of highest densities encountered in the grid cells

during the simulation to avoid possible errors due to density fluctuationsl.

The values for tbounce and pmax obtained from the different simulation runs of this
series are shown in table 4.1. Obviously, this data confirms that a greater initial
angular momentum leads to a later core bounce. There is also a clear tendency that
rotation causes lower densities at core bounce (although model CNEOSO4 unfortu-

nately breaks ranks without identifiable reason).

[EELglma is the ratio of the initial rotational energy to the initial gravitational
energy. It is common to use this ratio as a measure for how rapidly a supernova core

is rotating.

 

1Usually the roughly 200 highest densities encountered during the whole time development were
used for this.

62

Deviations from Spherical Symmetry

Figure 4.3 shows the development of the density distribution during the infall phase till
shortly after bounce for models CNEOSO3 (slowly rotating), CNEOSOS (moderately

rotating), and CNEOSO7 (rapidly rotating).

The clumping of the matter that becomes obvious in the plots labeled with (b) is
due to the wide attractive density region of the EOS: once there is a small density
fluctuation somewhere, it is amplified until nuclear density is reached because this
goes along with a decrease of the total energy. This leads to clumps of matter at
nuclear density whose size is essentially determined by the size of the grid cells. This

is obviously not a realistic feature and should be ignored.

A trivial feature of all three models is the oblate shape of the mass distribution
during the late stages of the infall phase (b). In the slowly rotating model a prolate
shape is apparent after bounce (e) - a feature that appeared in other simulations
before (e.g. [57]). In the rapidly rotating model two other interesting deviations from
spherical symmetry can be seen: during the late stages of the infall phase (c) the infall
velocity of the matter along the rotation axis becomes very large which cannot just
be due to the lack of centrifugal forces there because it leads to a slight “doughnut”
shape of the mass distribution. This higher infall velocity also leads to a significantly
larger density along the z-axis (the vertical in the plots) close to the center. One
might envision that along the rotation axis vortices are formed which help the matter

swirl to the center faster.

The assumption of spherically symmetric gravity is no longer completely justified
when such large deviations from spherical symmetry as in model CNEOSO7 occur.
Note, however, that this assumption tends to stabilize the core against large scale

deviations from spherical symmetry.

63

 

 

 

 

 

 

 

[ CNE0803 ] CNE0805 [ CNEOSO7]

1., 0 0 0

1,, 173.352 173.745 174.305
1. 173.907 174.287 174.645
t3 174.102 174.5 175.000
1,. 174.417 174.812 175.205
7., 1844 1844 1844
7,, 87.81 83.06 75.57
7. 58.72 53.60 53.29
rd 115.2 57.99 43.49
7. 148.7 95.05 54.23

 

 

 

 

Table 4.2: Key for figure 4.3. ta through te are the times (in ms) corresponding
to the density profiles labeled (a) through (e) in the figure. ra through re are the
corresponding radii (in km) of the shown density plots.

In model CNEOSO7 shortly after core bounce (e) the mass density in the equatorial
plane is much higher than that along the z-axis. A slight tendency towards this feature

can already be seen in the moderately rotating model CNEOSO5.

In the three core bounce plots (d) and after bounce (e) it can be seen how rotation
leads to a more diffuse mass distribution. Keep in mind, however, that the radius in

plot (d) for model CNEOSO3 is roughly two times larger than in the other two plots.

64

(a) .10ls g/cm"3
’ [OM

1013

 

 

1012

1011

108

   

  

107

       

H H O—.
CNEOSO3 CNEOSOS CNEOSO7
Figure 4.3: Mass density in a slice in the x-z-plane at five different times for models
CNEOSO3 (left), CNEOSOS (center), and CNEOSO7 (right). Events at the respective
times: (a) onset of simulation, (b) clumping and presence of centrifugal forces become
apparent, (c) formation of “vortices” in CNEOSO7, ((1) core bounce, (e) shortly after
core bounce. Note that the radii of the shown density profiles vary (see table 4.2 for
more data). The black lines below each plot indicate a length of 43.5 km. The length
scale for the plots labeled by (a) is much larger than this. Images in this thesis are

presented in color.

65

4.2.2 Cold Core with Polytrope Electron Contribution

It is desirable to use a realistic (inhomogeneous) density profile of a presupernova
core like that shown in figure 2.1 as initial condition. Unfortunately, this turned out
to be impossible in our model. It is quite important that the inner core matter (or
in our simulation the test particles representing it) arrive at the center of the core
simultaneously. Great accuracy is needed to achieve this given a change in the length

scale by a factor of about 102 during collapse.

Many other groups have tripped over this problem [5, 59]. One remedy (followed
e.g. by [59]) is to use (artificial) initial conditions that are an equilibrium solution of
the EOS used and then initiate collapse by suddenly reducing the electron fraction

(or the adiabatic index if a polytropic EOS is used) a little hit.

As the eflects of rotation should mainly play a role during the late stages of collapse
(when angular momentum conservation can lead to very high angular velocities) we
decided to use the following approach: the simulation is started when collapse is
already going on and the (inner) core has already contracted by a factor of 5. Thus
the problems arising from the scale change are simplified. The fact that the matter is
already falling inward is mimicked by imprinting an initial velocity profile 17(7') = kvr
(in addition to the rotation) on the test particles, where a realistic value 19,, = 89.1 s‘1
(estimated using the results of previous simulations [4]) is chosen. This is a good

model only for the inner core, so the mass is chosen to be 0.7MQ.

Since centrifugal forces are relatively weak during the early collapse stages the
initial mass distribution is still assumed to be spherically symmetric. The initial
mass density as a function of the radius calculated by Woosley and Weaver2 [56] and

the variant contracted by a factor of 5 we used in our simulation are shown in figure

 

2Strictly speaking, Woosley and Weaver’s data was approximated (extremely accurately) by the
analytic expression p(r) = 3:51-43; with appropriate parameters 161 and kg for figure 4.4.

66

 

 

 

 

 

l I I F l
L) I I I I I
”“57 <——-r ————— r ————— T ————— 1 ————— 1—-———i
'_ 1e+14
"a
a.
as
i“
g
c 1e+13
o .
Q
0’
3
5
1e+12
- I I I I
I I I I I
I I I I I
I I I I I ‘
l9-7Me I I I I 07”.
1e,“ 11 I I I I I
0 100000 200000 800000 400000 600000 600000

Distance from Center r [m]

Figure 4.4: Initial mass density as a function of the radius calculated by Woosley
and Weaver [56] (solid line) and the variant contracted by a factor of 5 we use in our
simulation (dashed line).The location at which the enclosed mass is 0.7MQ (edge of
the inner core) is indicated for both distributions on the abscissa.

4.4.

If collisions between the test particles were included in our code the use of fully
realistic initial conditions might become possible because these collisions should lead
to a local averaging of the infall velocity thus facilitating the simultaneous arrival of

most test particles at the center.

The angular velocity profile of Heger’s (rotating) progenitor (figure 2.5) indicates
that it is an excellent approximation to assume that the inner core (initially) rotates
like a rigid body. Therefore we used a constant angular velocity 010 as initial rotation.

However, some simulation runs were performed using a parameterized r-dependent

67

 

 

 

Name of Run (.00 [£23] [£0] [1041.18] :g‘glinit tbounce [ms] pm” [p0]
PCNEOSOl 10 0.467 0.027% 3.01 5.68
PCNEOSO7 70 3.27 1.3% 3.20 2.93
PCNEOSI3 130 6.06 4.5% 3.53 2.09
PCNEOSIQ 190 8.86 9.6% 3.86 1.66

 

 

 

 

 

 

 

 

 

Table 4.3: Simulation runs using the combination of the cold nuclear EOS and the
polytrope EOS for the electron gas. The abbreviations are explained in the text.

initial angular velocity profile [57]:

2
7‘0

—T2 + 73' (4.1)

020(7') = 02.;

This yields a very good approximation of the angular velocity of Heger’s progenitor
if the parameters wC (the central angular velocity) and r0 (which defines the degree
of differential rotation) are chosen suitably. As anticipated no significant differences

to the rigidly rotating models occurred.

In this series, the simple polytrope electron contribution was added to the nuclear
part of the EOS (see equation 3.19). The parameters th = 106, N, = 120, Ncosg =

100, and pm," 2 1.3 x 101155“;1 were applied.

Table 4.3 shows the data for all runs of this series.
Energy Conservation

Again, energy conservation is approximately fulfilled. Figure 4.5 shows a typical
energy development. We will not conceal, however, that in some of the other runs of
this series (especially PCNEOSOI and PCNEOSO7), energy conservation was not as

good as in figure 4.5, most likely because of a bad choice for the initial time step size.

68

3e+45

 

2e+45

1e+45

Energy [J]

-1645 ‘

 

 

 

 

|

-2e+45 1'-
l |
l |
-3e+45 [ :
| |
l l
-4e+45 1 l

0 0.0005 0.001 0.0015 0.002 0.0025 0.003 0.0035 0.004 0.0045
Time [s]

Figure 4.5: Time development of the different energies in simulation PCN E0819

Deviations from Spherical Symmetry

Figure 4.6 shows the development of the mass density for the three simulation runs
PCNEOSO7 (slowly rotating), PCNEOSIB (moderately rotating), and PCNEOS19

(rapidly rotating).

The peculiar clumping of the matter that was apparent in the previous series does
not occur in these models as can be seen in the plots labeled by (b) and (c). The
oblate shape of the density distributions is again visible in plots (b) and (c) and
obviously more pronounced in the more rapidly rotating models. The slowly rotating
model PCNEOSO7 shows a slight prolate shape shortly after core bounce in plots ((1)
and (e).

The formation of “vortices” along the rotation axis (the vertical in the plots)

69

 

I

[PCNEOSO7] PCNEOSlB [ PCNEOSl9 ]

 

 

 

 

 

 

 

 

ta 0 0 0

tb 2.00 2.00 2.00
to 3.38 2.82 2.53
td 3.86 3.53 3.19
t,3 4.14 3.93 3.64
ra 123.4 123.4 123.4
Tb 102.9 102.9 102.9
rC 56.10 56.10 56.10
rd 56.10 56.10 56.10
re 68.57 68.57 68.57

 

 

Table 4.4: Key for figure 4.6. to through t. are the times (in ms) corresponding
to the density profiles labeled (a) through (e) in the figure. ra through re are the
corresponding radii (in km) of the shown density plots.

during the infall phase can be seen in the two more rapidly rotating models in plots
(0) and (d). It is also quite obvious in plots (e) of these models that most of the mass

rebounds from the core in the equatorial plane.

The core bounce plots (d) show how rotation leads to a drastic deformation of the

core at bounce.

A curiosity occurs in model PCN E0807: in plot (c) a slight prolate shape can be
seen which is quite odd during the infall phase. This is most likely a freak due to the

(aforementioned) bad time step size used in this simulation run.

It is interesting to note that almost all of the features that appeared in the toy

model using the cold nuclear EOS reappeared in this much more realistic one.

70

(b)

 

 

 

_ l 01 1
“saw;-

o——o

PCNEOSO7 PCNEOS13 PCNEOS19

   

Figure 4.6: Mass density in a slice in the x-z—plane at five different times for mod-
els PCNEOSO7 (left), PCNEOSl3 (center), and PCNEOS19 (right). Events at the
respective times: (a) onset of simulation, (b) after 2 ms, (c) presence of centrifugal
forces and vortices become apparent, (d) core bounce, (e) shortly after core bounce.
Note that the plots have different radii. The black line below each plot indicates a
length of 56.1 km. See table 4.4 for more data. Images in this thesis are presented in
color.

71

 

 

 

 

Name of Run wo [L23] lfol [10“JS] 53%|...“ tbounce [1118] pm... [pol
HLSEOSOl 10 0.466 0.027% 3.44 2.99
HLSEOSO7 70 3.26 1.3% 3.61 2.35
HLSEOSIO 100 4.67 2.7% 3.76 1.63
HLSEOS13 130 6.06 4.5% 3.89 1.03
HLSEOSl6 160 7.46 6.8% 4.03 0.56
HLSE0819 190 8.86 9.6% 4.17 0.34
HLSEOS22 220 10.3 13% 4.31 0.21

 

 

 

 

 

 

 

 

 

Table 4.5: Simulation runs using the combination of the Helmholtz and the LS EOS.
The abbreviations are explained in the text.

4.2.3 Combination of Helmholtz EOS and Lattimer & Swesty
EOS

In this series the same initial conditions as in 4.2.2 were used. But as. EOS the most
realistic approach to actual core collapse conditions in this work was implemented —
the combination of the Helmholtz EOS for the e"/e+ gas at lower densities and the
LS EOS at higher densities. Both EsOS were described in 3.3.3. The parameters

M, = 106, N, = 110, Ncoso = 100, and pm," 2 1.3 x 1011;53- were applied.

Table 4.5 shows the data for all runs of the series. With this realistic EOS the
time of core bounce tbmmce and the maximum density pm, obviously depend on the

initial angular momentum in a strictly monotonous way.
Energy Conservation

Figure 4.7 shows a typical energy development for the models of this series. Again,
energy conservation is only fulfilled approximately. Temporary violations appear
especially near core bounce. However, the conversion of gravitational energy into

kinetic and internal energy during collapse is rendered quite well.

72

 

4e+~45

3e+45

2e+45

1e+45

  

 

Energy [J]
O

-1e+45

   

 

 

 

I
l I
I l
-uws ———r——— T
I I I
I I I
-uws————e-——e——— + +
I I I I
I I I I
«+45 I I l L I I l I
0 0.0005 0.001 0.0015 0.002 0.0025 0.003 0.0035 0.004 00045

Figure 4.7: Time development of the different energies in simulation HLSEOSl3.

Deviations from Spherical Symmetry

Figure 4.8 shows the development of the mass density distribution during the infall
phase till shortly after bounce for models HLSEOSO7 (slowly rotating), HLSEOSlfi

(moderately rotating), and HLSEOS22 (rapidly rotating).

The apparently different sizes of the density profiles are due to the grid scaling
described in 3.2.2 (the same radius is used in all plots in this figure): the abrupt
change between the light yellow and the white regions is located at the radial edge
of the grid. Again, the oblate shape of the more rapidly rotating models is obvious
in (b) and (c). A “vortex” formation along the z-axis during the infall phase is not
identifiable in any of the models. It is obvious that core bounce gets more diffuse with

increasing rotation (d). A slightly prolate shape shortly after bounce is apparent in

73

 

lab
the
den
end

1110
310.

thi:

 

 

 

 

 

 

F ]] HLSEOSO7 HLSEOSlG HLSEOS22 ]
ta 0 0 0
1,, 2.00 2.00 2.00
1.. 3.00 3.00 3.00
to: 3.61 4.03 4.31
t. 3.80 4.51 5.33
pm”... 5.13 x 1010 5.30 x 1010 5.30 x 1010
pm”... 6.15 x 1014 1.48 x 1014 5.93 x 1013

 

 

 

 

 

Table 4.6: Key for figure 4.8. ta through te are the times (in ms) corresponding to
the density profiles labeled (a) through (e) in the figure. pscmin and psama, are the
densities (in g/cm3) corresponding to the bottom end (color: light yellow) and top

end (color: black) of the

model HLSEOSO7. In the very rapidly rotating model HLSEOS22 a density depletion

along the z-axis (vertical in the plot) can be seen shortly after bounce (figure 4.9 shows

this in a magnified way).

density key.

74

I mm?!

(a) H g/cmA3

 

 

1014
(b)
1013
(C)
E 1012
(d) E
(e)
‘ -*‘-:~._‘ “A"- ‘6 1011
HLSEOSO7 HLSEOS16 HLSEOS22

Figure 4.8: Mass density in a slice in the x-z-plane at five different times for mod-
els HLSEOSO7 (left), HLSEOSlG (center), and HLSEOS22 (right). Events at the
respective times: (a) onset of simulation, (1)) after 2 ms, (c) presence of centrifugal
forces becomes apparent (after 3 ms), ((1) core bounce, (e) shortly after core bounce.
The density scale is only approximately valid: the highest density (indicated by the
colors black and dark red) decreases from left to right. All plots have the same radius
(123.42 km) indicated by the black line in the top left. See table 4.6 for more data.
Images in this thesis are presented in color.

75

g/cm"3

 

loll

Figure 4.9: Density depletion along the z-axis after bounce in model HLSEOS22 at
t = 5.21ms. The radius in this plot (measured from the center horizontally to the
right) is 51.43km. Images in this thesis are presented in color.

76

 

4.3 Possible Implications of these Results

Our simulation succeeds in reproducing several results of previous ones which used
different techniques of modeling rotating core collapse: we can confirm that core
bounce occurs later the faster the core is rotating initially. Also the ballpark for the
maximum density at bounce of pm, z 3p0 is substantiated by our results (obtained
from the slowly rotating models). More interesting is the validation that increasing
rotation leads to lower densities at bounce. Further, we can acknowledge that in
very rapidly rotating models (like HLSEOS220, HLSEOSl90, and HLSEOSI60) core
bounce occurs before nuclear matter density is reached (due to extreme centrifugal

forces) — a feature observed before (e.g. in [59]).

As “vortices” along the rotation axis appeared during the infall phase in the
rapidly rotating models in which simple artificial EsOS were used (sections 4.2.1 and
4.2.2) but not in the most realistic model (section 4.2.3), one may assume that this
feature is due to the (unrealistic) attractive density region of the E303 used in the

former models which is not existent in the EOS of the latter one.

A very interesting feature of all rapidly rotating models is the density depletion
along the rotation axis shortly after core bounce. In the rapidly rotating models,
most of the mass is shot out in (or near) the equatorial plane after bounce3. A highly
speculative conclusion out of this is that the shockwave created after core bounce
might be much more vigorous in the equatorial plane. Thus a prompt shock might
succeed there: the star might explode in the equatorial plane first ejecting a ring of

matter in this plane.

The density depletion along the rotation axis is extremely relevant for neutrino

 

3It is nearly impossible to see this in figures 4.3 and 4.8. But in the movies showing the full
motion of the density distributions it is a quite obvious feature.

77

transport. Even though neutrinos are not dealt with in our simulation at all, it
is evident that they would be able to escape from the central core much easier by
diffusing through the low density region along the rotation axis. The peculiar concave
shape of this low density region might actually “focus” a beam of neutrinos along the
rotation axis creating such a large neutrino flux there that neutrino heating could

lead to a jet-like explosion along this axis.

However, a lot more work needs to be done to underpin these speculations. Of par-
ticular importance is the implementation of collisions between individual test particles
in our model or some other way that actually yields the formation of a shockwave
after bounce. The inclusion of neutrino transport (at least in some simplified pa-
rameterized way) in our model is also most desirable if serious conclusions about the
eflects of neutrinos shall be drawn. After bounce our simple definitions of the tem-
perature and electron fraction as functions of the mass density break down — a more
sophisticated way of mapping data obtained from lower dimensional simulations on

our model is needed to remedy this.

Both phenomena we just suggested (the prompt explosion in the equatorial plane
and the delayed jet-like explosion along the rotation axis) are possible explanations for
the polarization observations by Wang and Wheeler [50, 49, 51] because it is known
that the light scattered from or through asymmetric surfaces is likely to be polarized

(see [50] for more detail).

78

T._—~_-
. I.. .

Chapter 5

Summary and Conclusion

First of all, the currently favored explosion mechanism of core collapse supernovae
was described. It was made clear that this mechanism is complicated and could not
be verified with great certainty so far since several different numerical simulations of
supernovae performed by different groups all over the world yielded both successful
explosions and failures. These deviating results are naturally due to uncertainties
in the input physics (e.g. neutrino transport, convection) and to different numerical

techniques used in the respective works.

The possible importance of rotation for the explosion mechanism was pointed
out. The most relevant indication for the breakdown of spherical symmetry is the
strong polarization of the light emitted by core collapse supernovae. Rotation is a
very appealing candidate to explain the large deviations from spherical symmetry
necessary to create such polarizations. This possible importance of rotation clearly

motivates a three dimensional simulation.

Thus, a simplified three dimensional simulation we introduced to simulate the
collapse of a supernova core during the infall phase until shortly after core bounce
was presented in detail. In this model the core’s mass is represented by discrete test

particles. A spherical coordinate grid is used to locally define thermodynamic quanti-

79

ties. Input from other simulations is used to reduce the computational requirements
of this simulation program to a minimum. The equations of motion for the test par-
ticles were derived making use of the assumption of approximate energy conservation
during the infall phase. This assumption is justified by the short time scale of col-
lapse which disables the neutrinos to carry energy out of the core in large amounts.
A weakness of our technique (at the current state of the art) is the lack of collisions
between the test particles which makes the formation of a shock wave after bounce

impossible.

Finally, our simulation program was used with different EsOS and initial condi-
tions to study the impact of rotation. Several (mostly predictable) effects of rotation
found in previous works could be validated. Two interesting new features appeared

in some of our rapidly rotating models:

1. A slight density depletion along the z-axis shortly after core bounce occurred
in all rapidly rotating models independent of the EOS used. At the same time
significantly more mass rebounded from the core in the equatorial plane than

along the z-axis.

2. During the late stages of the infall phase the formation of “vortices” along the
rotation axis was observed. However, this feature was extremely weak (hardly

visible) in the model using the most realistic EOS.

Careful conclusions about the possible effects these features may have for the ex-
plosion mechanism were drawn. We suggested that the explosion might be initiated
in the equatorial plane and along the rotation axis. This would deliver a possible ex-
planation for the light polarization. However, it was noted that several improvements
(such as the implementation of neutrino transport) should be made in our model to

assert the validity of these new features.

80

To conclude, let us state that the supernova problem remains unsolved. Great
theoretical challenges still lie ahead. It seems highly unlikely that computers powerful
enough to make a three dimensional simulation including realistic neutrino transport
and convection possible will become available in the near future. But probably not
before then will there be certainty about the explosion mechanism. If nothing else,
the present work has shown that the possible impact of rotation may have been

underestimated for a long time.

81

APPENDIX

82

Appendix A

Source Code of the Simulation
Program

In this appendix, the source code used for the simulation program which is written

in C++ will be reproduced.

A.1 suno.cpp and suno.h

These two files contain all main parts of the simulation C++ code. Most important is
the large C++ class Star which contains the memory structures and functions that de-
liver all information about the supernova core. The access functions for the FORTRAN
routines for the E803 are members of Star. Such trivial things as physical constants
and functions for unit conversions are also declared and defined in these files. The

whole code is extensively commented.

suno .h:

//////////////////////////////////////////////////////////////////
// suno.h

// Header file for suno.cpp
//////////////////////////////////////////////////////////////////

tirade! SUNO-H
#define SUNO_H

83

tinclude "fortran.h"
tinclude <stdio.h>
tinclude <iostream.h>
tinclude <time.h>
linclude <stdlib.h>
Cinclude <nath.h>
tinclude <fstream.h>
tinclude (vector)
tinclude "vector_and_epherical.h"
Cinclude "testparticle.h"
Sinclude <string>
tinclude <strstream>

const double c I 299792458; // global constant: speed of light [m/s]

const double 6 = 6.67259e-11; // gl. const.: gravitation constant [m‘3/kg/s‘2]
const double K-B = 1.3807e-23; // g1. const.: Boltzmann constant [J/K]

const double HASS_NEUTRON 8 1.6749286e-27; Ilkg

const double HBAR = 6.6262e-34 / 2.0 / PI; // global constant: h bar [Js]
//const double 6 = 0; // "switch off" gravity

// declarations of global functions
int compare_tp_dist(const voidt a, const void* b);

// functions for unit conversions

double nev_joule(double in_nev);

double joule_mev(double in-joule);

double mev_kelvin(double in_mev);

double kelvin_mev(double in-ke1vin);

double baryonperfm3_kgperm3(double in-bpf3);
double kgperm3-baryonperfm3(double in-kpm3);

// FORTRAN-function for electron gas equation of state

extern “C" void FORTRAN_NAME(helmholtzeos)(const doubles xmass, const
doublet aion, const doublet zion, const integer! ionmax, const doublet
temp, const doublet den. const doublet energs. const doublet
pressure);

// FORTRAN-function: initialization stuff for Lattimer t Swesty EOS
extern "C“ void FORTRAN_NAHE(loadmx)();

// FORTRAN-function for LtS EOS

extern "C” void FORTRAN_NAME(lseos)(const doublet ipvar, const doublet t_old,
const doublet y-e, const doublet brydns. const integer. iflag.
const integers eosflg. const integert forflg. const integert sf,
const doublet xprev, const doublet pprev, const doublet ls-out);

// forward declarations
class Star;

class TPHap;

class TPDoubleMap;
Cendif

suno . cpp:

tinclude "suno.h"
// global functions for unit conversions
double nev-joule(double in_nev){

return in-mev I 1.602e-13;}

double joule-nev(double in-joule){
return in-joule / 1.602e-13;}

double nev_kelvin(double in_nev){
return nev_joule(in_nev)/K_B;}

84

double kelvin_mev(double in_kelvin){
return joule-mev(in,kelvin*K-B);}

double baryonperfm3_kgperm3(double in,bpf3){
return HASS_NEUTRON*1e45#in_bpf3;}

double kgperm3_baryonperfm3(double in_kpm3){
return in_kpm3/MASS,NEUTRON*1e-45;}

//////////////////////////////////////////////////////////////////////
// class from which class TPDoubleMap is derived
//////////////////////////////////////////////////////////////////////
class TPHap

{
protected:

int i-r-max, i_phi-nax, i_cos_theta-max;

double r-max; // largest r coordinate of all partilcles to be saved
public:

bool SetIRMax(int nmax)

{

if(nmax>0)
{

i-r_max=nmax;
return false;

}

else

{
i-r_max=0;

cout << "Error in TPHap::SetIRMax(): i_r-max not positive.\n";

return true;

}
}
bool SetIPhiHax(int nmax)
{
if(nmax>0)
{
i_phi_max=nmax;
return false;
}
else
{
i_phi-max=0;
cout << “Error in TPHap::SetIPbiHax(): i_phi_max not positive.\n";
return true;
}
}
bool SetICosThetaHax(int nmax)
{
if(nmax>0)
{
i_cos-theta_max=nmax;
return false;
}
else
{

i_cos-theta_nar-O;

cout << "Error in TPMap::SetICosThetaHax():
i_cos_theta_mar not positive.\n";

return true;

}
}
bool SetRHax(double nr_mar)
{

if(nr_nar>-O)

{
r_nax I nr_mar;
return false;

85

else

r_max=1;
cout << "Error in TPMap::SetRMax(): r_max

not positive.\n";

};

return true;
}
}
// access functions (contd.)
int GetIRMax()

{
return i-r_max;
}
int GetIPhiMax()
{
return i_phi_max;
}
int GetICosThetaMax()
{
return i_cos_theta_max;
}
double GetRMax()
{
return r_max;
}

//virtual void SetAll2Zero()=0;

// calculates the volume of one box
double VolOfBox(int i_r);

// given a point as a Vector object, calculates
// i_r index of the box the point is in
int GetIR(Spherical);

// given a point as a Vector object, calculates
// i_phi index of the box the point is in
int GetIPhi(Spherical);

// given a point as a Vector object, calculates
// i_cos_theta index of the box the point is in
int GetICosTheta(Spherical);

// function defining box length in r direction
double F-R(doub1e rorm);

// inverse function of F-R()
double F_R,Inv(doub1e iroirm);

/##¢##*###$t*#¢i*fifiltt¢tlfifiittttfitifitttltttifitttttttt\

#

given a point as a Spherical object, calculates

* i-r index of the box the point is in.

0
j
m

i_r may be greater than i_r_max.
Gets: - Spherical spos
Returns: - i_r index of box as int

\sseeaeeeeaeseesseseasesssseeeasaesesstseeaeassattaee/
int TPHap::GetIR(Spherical spas)

{

int i-r;

i_r I (int)((i_r-nax) * F_R(spos.GetR()/r-nax));
// makes sure that i-r<i_r-max (appears obsolete)
//i_r 8 i_r < i_r_nax ? i-r : i-r_max - 1;

// changed in order to get zero potetial gradients
II for distant particles

return i-r;

86

}

Ititttttttttttitttttfitttitttttttltttttitttttttttttttt\
* given a point as a Spherical object, calculates

# i_phi index of the box the point is in

4 Gets: - Spherical spos

* Returns: - i_phi index of box as int
\fittififittltttttttttttttttttitfi*tttttfiltttttt##tittlti/
int TPMap::GetIPhi(Spherical spos)

{
int i,phi;
i_phi 8 (int)((i_phi-max) # spos.GetPhi()/2.0/PI);
// makes sure that i_phi<i_phi_max (appears obsolete)
i_phi 8 i_phi < i_phi_max ? i_phi : i_phi_max - 1;
return i_phi;

}

[seeaeaeeaaeaseaaessaeeeeaeasaeeateeaaseaeaaaaeesases\
t given a point as a Spherical object, calculates

* i_cos-theta index of the box the point is in

4 Gets: - Spherical spos

a Returns: - i_cos_theta index of box as int
\asaatssasaeecassettes:atsseeeastsceaeaeaeaeseesassse/
int TPMap::GetICosTheta(Spherical spos)

{
int i-c0s_theta;
i_cos_theta I (int)((i-cos-theta_max) # (cos(spos.GetTheta())+1)/2.0);
i-cos_theta 8 i_cos-theta < i-cos_theta_max ?
i_cos_theta : i_cos-theta_max-1;
return i_cos,theta;
}

/fit$#$##tifififitttttitififitiltitifittttfiiltittfifitittit‘fitttttt¥\
# Function defining box length in r direction by the

t relation F_R(r/r_max) I i_r/i_r_max

* for future experiments

4 Gets: - double r/r-max

4 Returns: - double F-R(r/r_max)
\eseesassaseeaeeceasaeeseseeeseaasestseaaaeaeetaeeesaeseeat/
double TPHap::F-R(double rorm)

{
return rorm;
//return pov(rorm,0.7);
//return pov(rorm,0.3);
}

lifititfifiittttttitlltttifitit¥$¥¥100ttfifittttttttttltt¥¥t¥¥#t*\

t Inverse Function of F_R
\aseeaeeeaesasasseeaseseasesaaesesaeaseaeeeaaeaeeaeaeeeseaa/

double TPMap::F-R_Inv(double iroirm)

{
return iroirm;
l/return pov(iroirm,1.0/0.7);
//return pov(iroirm,1.0/0.3);
}

[estatescease-assess:sesaeaeeeeeaaaeaeeeeeat03440334030aeaa\
t Calculates the volume of one box.
t This obviously depends on r respectively i_r.
t Gets: - i-r
t Returns: - Volume of all boxes on the shell labeled

t by i_r:

t Vol-Vol(shell)/(i_phi.maxti-eos-theta_max)
\eaasasaeasesaseeeeeeesaseaaaeaeaaaeeeeeeaeaesaaaeeeeeeeses/

87

double TPMap::VolOfBox(int i-r)

{

}

double vol;
double r0, r1; // radius of inner/outer shell boundary

// calculate inner/outer radii
r0 = r_max * F_R_Inv((double)i_r/(double)i_r-max);

r1 8 r-max e F-R_Inv((double)(i-r+1)/(double)i_r_max);

// volume of shell devided by no. of boxes
vol = 4.0*PI/3.0*(pou(r1,3)-pou(r0.3)) / i_phi_max / i_cos_theta_max;

return vol;

/////////////////////////////l///////////////////////////////////////
// class for storing mass density of test particles or (smeared)

// particle numbers in certain boxes
/////////////////////////////////////////////////////////////////////
class TPDoubleHapzpublic TPHap

{

private:

public:

double*** valueinbox; // ' the mass density in the boxes

// constructor 1

TPDoubleMap()

{
SetIRHax(1):
SetIPhiH81(1);
SetICosThetaMax(1):

}

// constructor 2
TPDoubleMap(int ni_r_max, int ni-phi_max. int ni-cos-theta-max)
{

ReserveMemory(ni_r-max, ni_phi_max, ni_cos_theta_max);

}
// destructor
‘TPDoubleMap()
{
int i=0, j=0;
// free memory
for(i=0; i<i-r-max; i++)
for(j=0; j<i_phi,max; j++)
delete[] valueinboinJEj];
for(i=0; i<i_r_max; i++)
delete[] valueinboxfi];
deleteE] valueinbox;
}

// reserves memory for the array containing the densities
void ReserveMemory(int ni_r_max, int ni_phi_max. int ni_cos_theta-max);

// sets test particle mass density in box (i_r. i-phi, i_cos_theta)
void SetValueInBox(int i-r. int i_phi, int i_cos,theta. double dens)
{

valueinbox[i_r][i_phi][i_cos_theta] I dens;

return;

}

// returns mass density in box (i_r. i_phi, i,cos_theta)
double GetValueInBox(int i_r. int i_phi. int i_cos_theta)
{

return valueinbox[i_r][i_phi][i-cos_theta];

}

// sets all entries in tpinboxflflfl to O
void SetAll22ero()

8E!

int i=0,j=0,k=0;
for(i=0; i<i_r_max; i++)
for(j=0; j<i_phi_max; j++)
for(k=0; k<i_cos_theta_max; k++)
valueinbox[1][j][k]=0;
}

// adds 1.0 to valueinbox in the box in which the
// TestParticle tp is located
void AddTP(Spherical spos);

// adds a total of 1.0 to the box in which the test particle is located
// and the neighboring boxes depending on the distance between the

// test particle’s position and the boundaries of the box

void AddTPSmeared(Spherical spos);

// function needed by AddTPSmeared
double Smear(double x. double x-0. double x,bound);

// sets values in boxes (i-r, i_phi. i_cos_theta) and
// ('.', i_cos_theta_max-1 - i-cos_theta) to the common average
void AverageCosThetaBoxes();

):

leeeeeeseeeeseessaeeeeeseseaeasseseeeseaeeeeeeeseeaeeea\
I Function of x: 0.5 if x==x_bound. 0 if
I xIIx_c, linearly interpolated in between.
I Gets: - doubles x, x_c. x_bound
I Returns: - value as described above
\eeeeeeeeeeeeeeeeeeeeeeeeeeassseeseeaassaeaeeeeseseeses/
double TPDoubleHap::Smear(double x, double x-c. double x_bound)
{
return O.5I(x-x-c)/(x-bound-x_c):

}

[it‘ll##tt‘ttttttttttttt‘t.‘0‘fittilittttttttttttttttttt\
I Smears a test particle that is located in a certain
box of the TPDoubleHap over the eight closest
neighboring boxes. The part of the test paticle (i.e.
a fraction of 1.0) distributed to each box depends
linearly on the distance between the particle and

the boundary of the box (in spherical coord’s).

For example: the box in which the test particle is
located gets at least 1/8 if the particle is on a
corner point of the box and 1.0 if the

particle is in its very center.

Gets: - Spherical spos (location of test particle)
Returns: nothing
\eeeeeesssaesa4006340033040assesesssaaeeeasaessseeeseee/
void TPDoubleMap::AddTPSmeared(Spherical spos)

{

Ilfiiiifiill

// indices of the box in which tp is located

int i-r. i-phi, i_cos_theta;

// indices of box whose value is to be modified (see below)
int i_r-mod. i_phi_mod, i_cos,theta-mod;

int k,m.n; // loop counters

double r.phi,theta; // tp’s position in spherical coord’s
double r-c. phi-c. theta-c; // center of tp’s box in s/c
double r_s. phi-s. theta_s; // lower boundaries of box
double r-b, phi,b, theta_b; // upper bounds of box

// don’t add particle if out of bounds of map
if(spos.GetR()>r_max)return;

// shorter names for variables...

r I spos.GetR();
phi I spos.GetPhi();

89

}
/

§§I§§

theta I spos.GetTheta();

// calculate the indices of the box in which
// the test particle is located

i_r I GetIR(spos);

i_phi I GetIPhi(spos);

i_cos-theta I GetICosTheta(spos);

// calculates boundaries of box in s/c
r_b I r_max I F-R_Inv((double)(i_r+1) / i_r-max);
r-s I r_max I F_R-Inv((doub1e)i-r / i-r_max);

phi_b I 2IPII (double)(i,phi+1)/i-phi_max;
phi_s I 2IPII (double)i_phi/i_phi_max;

theta_b I acos(2.0*(double)i_cos_theta/i-c0s_theta_max-1.0);
theta_s I acos(2.0I(double)(i_cos_theta+1)/i_cos_theta_max-1.0);

// calculates center of box in s/c
r_c I 0.5*(r-b+r_s);

phi_c I 0.5I(phi_b+phi-s);

theta-c I O.5I(theta_b+theta_s):

// actually smears test particle over 8 closest neighboring boxes
for(kI0; k<2; k*+)

{

i_r_mod I i_r+k*(r>r_c?1:-1);

// makes sure that i-r_mod >I0 and <i_r_max

i_r-mod I i_r_mod<0?0:(i_r_mod<i_r_max?i_r-mod:i_r_max-1);

r_b I r>r c7r b:r s;

for(m=0; m<2; m++)

{
i_phi_mod I i_phi+mI(phi>phi_c?1:-1);
// makes sure that i_phi-mod >=O and <i_phi_max
i-phi_mod I i_phi_mod<0?i-phi_max-1:
(i_phi_mod<i_phi_max?i_phi_mod:0);
phi,b I phi>phi,c?phi_b:phi_s;
for(nI0; n<2; n++)
{
i_cos-theta_mod I i_cos_theta+n*(theta>theta-c?1:-1);
// makes sure that i_cos_theta_mod >=0
// and <i-cos_theta_max
i_cos_theta_mod I i_cos_theta_mod<0?0:
(i,cos_theta_mod<i_cos_theta-max?
i-c0s_theta_mod:i_cos_theta_max-1);
theta-b I theta>theta_c?theta_b:theta_s;
valueinboin.r.mod][i_phi.mod][i_cos.theta_mod]+I
(1.O-k+(kIIO?-1:1)ISmear(r,r-c.r-b))I
(1.0-m+(m==0?-1:1)ISmear(phi.phi_c,phi-b))I
(1.0-n+(nII0?-1:1)*Smear(theta.theta_c.theta_b));
}
}
}
return;

 

::t:t:t:::::;:::t:::::::::::::::::::::::::::::essee¢eee\
Adds 1.0 to valueinbox in the box in which

the test particle at position spec is located.

(valueinbox I 8test particles)

Gets: - Spherical spos (location of test particle)

Returns: - nothing

\teaeeeeesseseeeaeseeeeeeeeseeaseeteasesseesasasaseseetasset]
void TPDoubleHap::AddTP(Spherical spos)

{

// indices of the box in which tp is located

90

int i_r, i_phi. i-cos_theta;

// don’t add particle if out of bounds of map
if(spos.GetR()>r_max)return;

// calculate the indices of the box in
// which the test particle is located
i_r I GetIR(spos):

i_phi I GetIPhi(spos);

i_cos_theta I GetICosTheta(spos);

// actually adds the test particle
valueinbox[i_r][i_phi][i_cos-theta]++;
return;

}

[IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII\

I Reserves memory for the requested number

I of boxes and sets all entries in valueinbox[][][]

I to zero.

I Gets: - three ints indicating the number of

I separations in the three coordinates

I Returns: nothing

\IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII/

void TPDoubleMap::ReserveHemory(int ni_r_max, int ni,phi,max,
int ni_cos_theta_max)

{
int i=0. j=0;
// memory reservation
if(SetIRHax(ni_r_max) II SetIPhiMax(ni_phi_max) ll
SetICosThetaMax(ni_cos_theta_max))
{
return;
}
else
{
valueinbox I new doubleII[i_r_max];
for(i=0; i<i_r_max; i++)
valueinboin] I new doubleI[i_phi-max];
for(iIO; i(i_r_max; i++)
for(jIO; j<i_phi_max; j++)
valueinbox[i][j] I new double[i_cos_theta_max];
SetAll2Zero();
}
}

IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII\
I Sets values in boxes (i-r, i_phi. i_cos_theta) and

I (',', i_cos_theta_max-1 - i_cos_theta) to their

I common average. Sensible if mirror symmetry about

I the equatorial plane is assumed.

I Gets: nothing

I Returns: nothing
\IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII/
void TPDoubleHap::AverageCosThetaBoxes()

{

double av-va1I0; // average value for current pair of boxes

// average all the boxes...
for(int i=0; i<i_r_max; i++)
for(int jIO; j<i_phi_max; jII)
for(int kIO; k<i_cos_theta_max/2; k++)
{
av-val I 0.5 I (valueinbox[i][j][k] +
valueinbox[i1[j][i_cos_theta-max-1-k]);
valueinboxfi][j][k] I av_val;

Ell

}

valueinboinJEj][i_cos_theta_max-1-k] I av_val;

return;

//I////////////////////////////////////////////////////////

// main

class of program

//////////////////////////////////////////////////////////I

class Star

{

private:

public:

TestParticleI tp;

int number_tp;

double radius; // radius of star [m]

double m_star; // mass of star [kg]

double m_tp; // mass of test particle [kg] (IIm_star/number_tp)
double A; // parameter in nuclear eqn. of state (EOS): -3.49236e-11 J
double SIGMA; // parameter in nuclear EOS: 4.0/3.0

double 8; // parameter in nuclear EOS: 2.62728e-11 J

double DENS_HIN; // minimum density enforced

// tp_idEi] is the ID of test particle i

intI tp-id;

// id_tp[i] is the test particle number of the t/p with ID i

intI id-tp;

// test particle mass density map (using a spherical coordinate grid)
TPDoubleMap dmap;

// potential map, potential due to £05

TPDoubleHap potmap;

l/flag indicating if thermalization is activated

bool flag-thermalize;

// flag indicating averaging of mass density in cos(theta)-direction
bool flag_average_theta;

// map containing indices of particles located in each box
vector<int>eee box_tp;

// input variables for L8 EOS. saved for faster convergence of £08
doubleIIII lseos_ipvar;

// another input variable for LS EOS

doubleIII lseos_pprev;

// flag indicating if "energy method" is used instead of "force method"
bool flag-energy_method;

double time_passed; // time in seconds

// constuctor
Star(int number, double radius. double mass, int i-r_max,
int i_phi_max. int i-c0s_theta_max. bool thermalize,
bool average-theta, double newdens_min, bool energy-method)

number-tp I number;

// reserve memory for test particles
tp I new TestParticle [number_tp];
tp_id I new int [number_tp];

id_tp I new int [number_tp];

// set test particle IDs

for(int iIO; i<number_tp; i++)

{

tpti].SetID(i);
}
SetRadius(radius);

SetHass(mass);
RefreshTP();

dmap.ReserveHemory(i-r-max, i_phi-max, i,cos_theta-max);
potmap.ReserveMemory(i_r_max, i_phi-max. i-cos_theta_max);

ll reserves mom for box-tp only if thermalization is activated

92

 

if((flag-thermalizeIthermalize) II true)
ReserveMemory4BoxTP(i_r_max, i-phi_max, i_cos_theta_max);

// reserves mom for lseos_ipvar and lseos,pprev
ReserveMemory4LSEOSInput(i_r_max. i_phi-max, i-eos_theta_max);

flag-average_theta I average_theta;

// sets parameters for E08 (old)
AI-3.49236e-11;

SIGMA=4.0/3.0;

B=2.62728e-1l;

DENS_MIN I newdens-min;

// sets flag showing if energy or force method is used
f1ag_energy-method I energy_method;

time_passed I 0;

}

// destructor
'Star()
{
// free memory reserved for test particles
deleteE] tp;
deleteE] tp-id;
deleteE] id_tp;

// delete box_tp array if thermalization is activated
if(flag_thermalize II true)
DeleteBoxTP();

// delete arrays lseos_ipvar and lseos_pprev
DeleteLSEOSInput();
}

// refreshs tp_id and id_tp
void RefreshTP();

// sets radius of star
void SetRadius(double newr)
{
radius I newr;
return;

}

// sets the star’s mass
void SetMass(doub1e newmass)

{
m,star I newmass;
m-tp I m-star/number_tp;
return;

}

// returns the mass of a test particle
double GetMTP()
{

return m_tp;

}

// returns radius of star
double GetRadius()
{

return radius;

}

// returns the number of test particles
int GetNumberTP()

93

{

return number_tp;

}

// returns test particle i
TestParticle GetTP(int i)
{
if(i>I0 it i<number-tp)
return tp[i];
else{
cout << "Error in class Star function GetTP():
Bad test particle number.\n";
return tp[01;} // return garbage
}

// sets test particle i
void SetTP(int i, TestParticle newtp)
{
tpEi] I newtp;
}

// returns the test particle mass density map
TPDoubleMap GetTPDensityMap()
{

return dmap;

}

// distributes test particles in position space
void DistributeTPPos(int no_shells);

// puts test particle i to a random position in a sphere of radius r
bool PutTPInSphere(int i, double r);

// mass density as a function of distance from
// the center (initial condition)
double initial_mass-density(doub1e r, double parameter_d);

// distributes test particles in momentum space imposing rigid
// body rotation
void DistributeTPMom(Vector omega, double r_O, bool diff-rot);

// determines no. of test particles with smaller distance
// to the origin than current t/p
int CountTPInside(int current_tp);

// calculates & returns the derivatives of position
// and momentum of particle i
TPChange Deriv(int i);

// calculates the next step in the star’s time-development
// using the Euler method
void NextStepEuler(double stepsize);

// calculates the next step in the star’s time-development
// using a 4th order Runge Kutta algorithm
double NextStepRK4(double h. bool modify_stepsize_avvel);

ll calculates the current TPChange for all test particles
void CalculateK(TPChangeI k. double dt);

// moves all test particles to Iinitial I kfactor I Ik
void HoveTPs(TestParticleI initial, TPChangeI k, double kfactor);

// performs a certain number of steps

void Stepper(double stepsize. int number_of_steps.
int points_per_step_save. bool modify_stepsize,
bool modify_stepsize-avve1);

94

// modifies stepsize used by Star::Stepper()

// (based on density distribution)

double ModifyStepsize(double dt, boolt stepsize-modified,
boolt stepsize_modified_tvice);

// modifies stepsize (based on averag velocity of t/ps)
double ModifyStepsizeAvVel(double dt, double av_disp);

// sorts the test particles by their distance from the origin
void SortTP():

// saves the cordinates of every (number_tp/100)th test particle
// in a file readable by SuNoOutput.exe created with Visual Basic
void SaveCoordinatesVB(int points,per_file);

// saves the test particle mass density in a slice through the star
void SaveTPMDensityVB();

// saves the whole data (positions and momenta) for one time step
void SaveAllData(doub1e max-radius, int timestep);

// reads the whole data (positions and momenta) for use as
// initial cconditions
void ReadAllData();

// calculates the TPSmearedNumberHap for the current config of the star
void Star::CalculateTPSmearedNumberHap(TPDoubleMapt nmap);

// calculates the mass density map
void CalculateTPMDensityHap();

// calculates the potential map, potential due to 805
void CalculateEOSPotentialMap();

// calculates the temperature map
void CalculateTemperatureMap();

// calculates the temperature in HeV as a function of the density
double CalculateTemperature(double rho);

// calculates the gradient of the EOS potential \nabla V
// at the location of test particle i
Vector GetGradEOSPotential(int i);

// calculates the derivative of the EOS potential in the
// r-direction in a more sophisticated vay
double CalculateInterpolatedDVDR(int i_r, int i_phi.

int i_cos_theta. double r);

// calculates the derivative of the EOS potential in the

// phi direction in a more sophisticated way

double CalculateInterpolatedDVDPhi(int i_r. int i_phi,
int i_cos_theta. double phi);

// calculates the derivative of the EOS potential in the

// ”cos(theta)-direction" in a more sophisticated vay

double CalculateInterpolatedDVDCosTheta(int i_r, int i-phi.
int i_cos-theta, double costheta);

// returns the force on test particle i
Vector GetForceOnTP(int i);

// equation of state. returns potential in box i_r, i_phi, .
double CalculateEOSPotential(int i_r. int i-phi. int i_cos_theta);

// returns the flag indicating if thermalization is activated

bool GetFlagThermalize(){
return flag-thermalize;}

95

// randomizes momenta of test particles in box (i_r.i_phi,i-cos_theta)
void ThermalizeBox(int i_r. int i_phi, int i_cos_theta);

// reserves memory for the box_tp
void ReserveMemory4BoxTP(int i-r-max. int i_phi_max, int i_cos_theta_max)
{
int i=0, j=0;
box-tp = new vector<int>tt [i_r_max];
for(i=0; i<i_r_max; i++)
box_tp[i] = new vector<int>t [i_phi_max];
for(i=0; i<i_r_max; i++)
for(j=0; j<i_phi_max; j++)
box_tp[i][j] = new vector<int> [i_cos_theta_max];
return;

}

// reserves men for lseos-ipvar and lseos_pprev, sets initial guesses
void ReserveMemory4LSEOSInput(int i,r_max, int i_phi_max.
int i_cos_theta,max)

 

{
lseos_pprev 8 new doublet! [i_r-max];
lseos_ipvar = new double*** [i_r_max];
for(int i=0; i<i,r_max; i++){
lseos-pprev[i] s new double. [i-phi_max]; g
lseos_ipvar[i] I new doublett [i-phi_max];}
for(int i=0; i<i_r-max; i++)
for(int jao; j<i_phi_max; j++){
lseos_pprev[i][j] 8 new double [i_cos_theta,max];
lseos_ipvarEiJEj] 8 new doublet [i-cos_theta_max];}
for(int i=0; i<i-r-max; i++)
for(int j=0; j<i-phi-max; j++)
for(int k=0; k<i_cos_theta_max; k++)
lseos-ipvar[i][j][k] I new double [4];
double pprev 8 DENS_HIN/MASS-NEUTRON*1e-45 t 0.3;
// set values for initial guesses
for(int i=0; i<i-r_max; i++)
for(int j-O; j<i_phi_max; j++)
for(int k-O; k<i-cos,theta_max; k++){
lseos_ipvar[i][j][k][1]=0.155;
lseos,ipvar[i][j][k][2]=-15.0;
lseos_ipvarEiJEj][k][3]=-10.0;
lseos_pprev[i][j][k]= pprev;}
return;
}

// deletes box_tp array
void DeleteBoxTP()
{
int i=0, j=0;
for(i=0; i<dmap.GetIRMax(); i++)
for(j=0; j<dmap.GetIPhiMax(); j++)
delete[] box_tp[i][j];
for(i=0; i<dmap.GetIRMax(); i++)
delete[] box_tp[i];
delete box_tp;
return;

}

// deletes arrays lseos_ipvar and lseos_pprev
void DeleteLSEOSInput()
{

for(int i=0; i<dmap.GetIRHax(); i++)

for(int j-o; j<dmap.GetIPhiMax(); jH)
for(int k80; k<dmap.GetICosThetaHax(); k++)
delete[] lseos_ipvar[i][j][k];
for(int i=0; i<dmap.GetIRHax(); i++)

96

};

for(int j=0; j<dmap.GetIPhiHax(); j++){
delete[] lseos,pprev[i][j];
delete[] lseos-ipvar[i][j];}
for(int i=0; i<dmap.GetIRMax(); i++){
delete[] lseos-pprev[i];
delete[] lseos_ipvarEi];}
delete[] lseos_pprev;
delete[] lseos_ipvar;
return;

}

// deletes all entries in box_tp array
void C1earBoxTP(){

int i,j,k;

for(i=0; i<dmap.GetIRMax(): i++)

for(j=0; j<dmap.GetIPhiMax(); j++)
for(k=0; k<dmap.GetICosThetaHax(); k++)
box_tpEi][j][k].c1ear();
return;}

// measures the total angular momentum of the star
Vector HeasureTotalAngularMomentum();

// measures the kinetic energy of the star
double HeasureKineticEnergy(bool relativistic):

// measures the gravitational energy of the star
double MeasureGravitationalEnergy();

// measures the internal energy of the star
double MeasureInternalEnergy();

// internal energy per baryon via Helmholtz EOS from Frank Timmes

double CalculateInternalEnerngelmholtz(const doublet xmass.
const doublet aion, const doublet zion, integer ionmax,
double temp, double den, const doublet pressure);

// internal energy (and other thermodynamic qutts.)

// via Lattimer & Swesty EOS

double CalculateLSEOS(doublet ipvar. double y_e, double density,
double! pprev, const doublet ls_out);

// modifies the momenta of all t/ps conserving the total energy
void HodifyMomentaConservingEnergy();

// sets the initial internal energies for the t/ps
void InitializeInternalEnergies();

// gives t/ps a momentum towards the center of the core
void KickTPInside(double max_velocity);

// calculates and returns electron fraction as function of density
double CalculateYE(double rho);

[lifittllttfifit‘ittllittititiltfilititttltltlitifittttt.tt¢¢fitt‘¢¥t\

.iG‘IOIOi‘I

t

Calculates and returns electron fraction as function of
density using the values from cooperstein.

This is a good approximation for

(and only for) the infall phase till core-bounce. The
data is taken from Cooperstein I Hambach, Nucl. Phys. A
420, p.591 . 1984 and Cooperstein, Nucl. Phys. A 438.
p.722, 1985 and was approximated by analytic expressions
using the gnuplot fit function.

Gets: - mass density rho in kg/m“3 as double

Returns: - electron fraction as double

\eteeoeeeeeeeseeeeteeeseeeeeesseeeeetsseeeteetaeesseerettttwete/
double Star::CalculateYE(double rho)

97

}

double dens 8 rho/1e14;
return O.921883/(8.02162+dens)+0.307112;

Iiii...tittit###*t###tfittfifittttttttfit##tttfit#¢#¥##$##$#¢t¢*tttt\

§§I§§§§§§§CGGOOI§§§§§§

.

Calculates and returns the internal energy using the
EOS by Lattimer and Swesty which includes a nuclear
contributions and is thus good at high densities.
FORTRAN_NAHE(loadmx)() must be called before this
function can be used. Parameters from bound.atb and
maxwe1.atb are used (see EOS documentation for details),
so these files have to be in the directory of suno.exe.
Gets: - inputvariables ipvar as double array.

ipvar[O] =8 temperature in HeV
ipvar[1] 8* guess at nuclear density in fm‘-3
ipvar[2] =8 guess at proton eta
ipvar[3] =8 guess at neutron eta

(save old values of these box for every new call)

- electron fraction y-e as double

- baryon mass density (in kg/m‘3) dens as double

- pprev as double reference. this is the guess at
the exterior proton fraction

- ls-out[4] as double array. Here pressure[MeV/fm‘3],
internal energy/baryon[MeV], entropy/baryon [k-B]
and free energy/baryon [7] are stored (in this order).

Returns: - internal energy per baryon [J] as double
Note: The L8 EOS is able to calculate a LOT more thermodynamic

quantities than just these.

\tetsweetest:assesses:twas:attettttsaettettstttttestostetteeteeeet/
double Star::CalculateLSEOS(double* ipvar, double y_e. double density,

{

doublet pprev, const doublet ls-out)

static bool firstcallstrue;

// initial guess for temperature (only used if

//temperature is not input variable)

static double t_old;

// needed by LS EOS. can be ignored

static double xprev;

// indicates ipvar[O]: 1=temperature, 2=interna1 energy, 3=entropy
static integer iflagai;

// indicates which EOS was used

static integer eosflg;

// forces use of £08 indicated by eosflg if set to 1 (don’t...)
static integer forf1g=0;

// sf=1 means successful EOS call, everything else not

static integer sf;

// conversion kg/m‘3 --> fm‘-3

double dens 8 kgperm3_baryonperfm3(density);

// initialization if this is first call of function
if(firstcall)
{

FORTRAN_NAHE(loadmx)():

firstcall I false;

}

// check if density is ok, set initial guesses for input
// variables if not
if(density > 1e17){
ipvar[1]-0.155;
ipvar[2]=-15.0;
ipvar[3]=-10.0;
tpprev=DENS_HIN/HASS_NEUTRON¥1e-45 # 0.3;
}

FORTRAN_NAHE(lseos)(ipvar, tt_old, &y_e. tdens, tiflag, leosflg,
tforflg. tsf. txprev, pprev, ls_out);

98

return mev_jou1e(ls_out[1]): // conversion MeV --> J

}

Itit.ttttttttttttttttttttt!ttttttttttttttttttitttttttttttt\
* Calculates and returns the internal energy using the
* Helmholtz EOS from Frank Timmes which does NOT include
t a nuclear contribution.
* Gets: - ionmax as integer. This is the number of
different ion types in the matter
- xmass as double array containing the mass
fractions of the different ions
- aion as double array containing the mass numbers
of the ions (=number of nucleons)
- zion as double array containing the charge
numbers of the ions (=number of protons)
- temperature temp in Kelvin as double
- mass density den in kg/m‘3 as double
Returns: - specific internal energy [J/kg]
Note: - the limits of the EOS are:
1e4<= temp (8 1e11, 1e-7 <= y_e*den (3 1e14
- the Helmholtz EOS is able to calculate a LOT
t more thermodynamic quantities than just this
\teasesotetttemseetttttttssat:eteeettsseteeeaeeseestsstttt/
double Star::CalculateInternalEnerngelmholtz(const doublet xmass,
const doublet aion, const doublet zion, integer ionmax.
double temp, double den, const doublet pressure)

I§§QII§§I§§I§

double energ=0;

// check if temperature and density are ok
if(den < 1e-7 ll den > 1e14 ll temp < 1e4 ll temp > 1e11)
{
cout << "Error in Star::CalculateInternalEnerngelmholtz():
bad value for density or temperature.\n";

}

// convert den to g/cm“3
den *8 1e-3;

// calls the Helmholtz EOS. saves specific internal energy in energ
FORTRAN_NAHE(helmholtzeos)(xmass.aion,zion,tionmax,&temp.tden,
lenerg.pressure);

// convert erg/g to J/kg
energt-le-4;

return energ;

}

[titlttttittififilttti$.01!!!it!tit!tifitttttitlttttltlttt‘tt\
Calculates and returns the temperature as a function

of the mass density. This is a good approximation for
(and only for) the infall phase till core-bounce. The
data is taken from Cooperstein t Hambach, Nucl. Phys. A
420. p.591 , 1984 and Cooperstein, Nucl. Phys. A 438,
p.722. 1985 and was approximated by analytic expressions
using the gnuplot fit function.

Gets: - mass density rho in kg/m‘3 as double

t Returns: - temperature in K as double
\tttttttttttttttttltltttttt*tttetttttttttttttittttttetteve]
double Star::CalculateTemperature(double rho)

{

Di‘lil’fili

double tIO;

// converts rho from kg/m‘3 to 1e11 g/cm‘3
rhottle-14;

99

if(rho<-1e1)
t=1.48962tpow(rho,0.182321)-0.280647tpow(rho,0.446682);
else
if(rho>1e1 at rho<=4e1) // interpolates between the two data sets
t=(4e1-rho)/3e1t(1.48962tpow(rho.0.182321)
-O.280647tpow(rho,0.446682))
+ (rho-1e1)/3e1¢(0.148522tpow(rho,0.485484)+0.728965);
else
if(rho>4e1)
t80.148522tpow(rho,0.485484)+0.728965;

return mev_kelvin(t);

}

[tettttttttttttttttttttttttttttttt*ttttttttttttttttttttttt\
a Calculates and returns the total internal energy (i.e.

t the energy due to the EOS) of the star.

* Gets: - nothing

I Returns: - total internal energy as double
\tittittttttttttttttttttttttttttttttfitttttttttitttttitties]
double Star::MeasureInternalEnergy()

{
double e-int=0;
for(int i=0; i<potmap.GetIRMax(); i++)
for(int j=0; j<potmap.GetIPhiHax(); j++)
for(int k=0; k<potmap.GetICosThetaHaxC); k++)
{
e_int+=potmap.GetValueInBox(i,j.k)‘
dmap.GetValueInBox(i,j.k)/m-tp*dmap.VolOfBox(i);
// could be simplified if numbermap was saved
}
It // alternative method
for(int i=0; i<number_tp; i++)
e_int += tp[i].GetInternalEnergy();
1"/
return e_int;
}

Ittettttttttttittttetttttttttttit.Itttttttttttttttttlttttt\
* Calculates and returns the gravitational energy of the

a star (using a crude approximation).

* Gets: nothing

i Returns: - gravitational energy as double
\tttttl*tittttttitfitltlttttttttttttitttttttttttttfitfitttttt/
double Star::HeasureGravitationalEnergy()

{
double e_grav=0;
SortTP();
for(int i=0; i<number-tp; i++)
e_grav += -G*m_tp*itm_tp/tp[i].GetPos().GetNorm();
return e_grav;
}

[testtestttteststeeeeoseevaassesseaaesesttestosteeeetestes\
* Calculates and returns the total kinetic energy of the

t star. The flag relativistic indicates if calculation is

t to be performed relativistically or not.

i Gets: - relativistic as bool

1 Returns: - total kinetic energy as double
\etttttfitttttttttttittttittttttttettettstlttttttttetiltits]
double Star::HeasureKineticEnergy(bool relativistic)

{

double e-kin=0;

100

if(relativistic)

{
for(int i=0; i<number-tp; i++)
e_kin +8 sqrt(pow(m_tp,2)tpow(c,4)+
tp[i].GetHom().GetNorquuared()#pow(c,2))-m-tptpow(c,2);
}
else
{
for(int i=0; i<number_tp; i++)
e_kin+=tp[i].GetMom().GetNorquuared() / 2.0 / m_tp;
}

return e-kin;

}

[eteesseeeaseeaseteesesteeeseeeeeeeaeeeeteeeeeeeeeeeetetet\
* Calculates and returns the total angular momentum

t vector of the star.

t Gets: nothing

i Returns: - angular momentum as Vector
\eeeeeeeeeeeeeeeeeeeeeeesteeeeeeeeeeeeseeeetattesteeeteeet/
Vector Star::HeasureTotalAngularHomentum()

{
Vector angmom;
for(int i=0; i<number_tp; i++)
angmom 8 angmom + (tpEiJ.GetPos()8tp[i].GetMom());
return angmom;
}

[attesttsetestimatesesteem.teeteeeseeeeeeeteeeeteteeete#e#\
Averages the motion of the test particles in box

(i_r, i_phi. i-eos-theta) not (E) conserving kinetic
energy, but consering angular momentum.

Only the projections of the momenta of the t/ps in the
radial direction are modified.

Gets: - subscripts of box

Returns: nothing
\eeeeeeeeeeeeetteteeetttsetse:esteem!teetetteeeeeeeeeeeset/
void Star::ThermalizeBox(int i_r, int i_phi, int i_cos_theta)

{

{O‘Ifi'i‘l’

Vector new_mom; // new momentum vector
int n80;
int no_tp_box=0; // number of test particles in current box

no_tp-box 8 box-tp[i-r][i_phi][i_cos_theta].size();
if(no_tp_box)
{

Vector e-r; // unit vector in spherical coordinates

double p_r80; // projection of moentum on e_r

for(n=0; n<no_tp_box; n++)

{
// calculates e_r for current test particle
e_r 8 tp[box_tp[i_r][i-phi][i_cos_theta][n]].GetPos()
* (-1/tp[box_tp[i_r][i_phi][i_cos_theta][n]].GetPos().GetNorm());
// calculates and adds p_r of the particles
p_r+=tp[box_tp[i-r][i-phi][i_cos_theta][n]].GetHom()‘e_r;
}

// calculates average p_r
p_r /8 no_tp_box;

for(n80; n<no_tp_box; n++)

{
e_r 8 tp[box_tp[i_r][i_phi][i_cos_theta][n]].GetPos()

101

}

8 (-1/tp[box_tp[i_r][i-phi][i_cos_theta][nJJ.GetPos().GetNorm());
new_mom 8 tp[box-tp[i_r][i_phi][i-cos_theta][n]].GetHom();
// subtracts old radial part
new,mom 8 new_mom - (new_mom‘e_r)8e-r;
// adds new radial part
new_mom 8 new_mom + p_r8e_r;
tp[box_tp[i_r][i_phi][i-cos_theta][n]].SetMom(new_mom);

return;

}

[fit#tfitfifitltttttit#itttttififiltittittittfittittttfl¢#¢#$\

{09%}!

Calculates and returns the gradient of the potential
due to the EOS at the position of test

particle i as a Vector.

TPDoubleHap potmap must be calculated before this works
in a reasonable way.

Gets: - int i identifying the test particle

8 Returns: - gradient Vector
\ttttfittttttfitttttttttittttetetetttttttttttttettttttt/

Vector Star::GetGradEOSPotential(int i)

{

l/ for the resulting gradient

Vector grad;

Spherical spos;

// indices of the box test particle i is in
int i_r80, i_phi80. i_cos-theta80;
//derivatives

double dV_dr80, dV_dphi8O. dV_dcostheta80;
// unit vectors for spherical coordinates
Vector e_r, e_phi. e_theta;

// converts the test particle’s position vector to spherical coord’s
spos 8 tp[i].GetPos().GetSpherical();

// locates test particle in TPHap

i_r 8 potmap.GetIR(spos);

i-phi 8 potmap.GetIPhi(spos);
i_cos_theta 8 potmap.GetICosTheta(spos);

// calculates gradient only if test particle is in map.
// i.e. i_r<i-r_max
if(i_r<potmap.GetIRMax())
{
ll calculates derivatives in spherical coordinates
dV_dr 8 CalculateInterpolatedDVDR(i_r, i-phi.
i_cos_theta, spos.GetR());
//dV_dr 8 Ca1culateDVDR(i-r. i_phi, i_cos_theta);
dV_dphi 8 CalculateInterpolatedDVDPhi(i_r, i_phi,
i-eos_theta, spos.GetPhi());
//dV_dphi 8 Ca1culateDVDPhi(i-r, i-phi, i_cos_theta);
dV_dcostheta 8 CalculateInterpolatedDVDCosTheta(i_r;
i-phi, i_cos_theta, cos(spos.GetTheta()));
//dV_dcostheta 8 CalculateDVDCosTheta(i_r, i_phi, i_cos_theta);

// sets cartesian basis vectors for spherical coordinates
e_r 8 tp[i].GetPos() 8 (1/tp[i].GetPos().GetNorm());

e_phi.SetCoords(-sin(spos.GetTheta())*sin(spos.GetPhi()).
sin(spos.GetTheta())#cos(spos.GetPhi()). 0);

e_theta.SetCoords(cos(spos.GetTheta())tcos(spos.GetPhi()),
cos(spos.GetTheta())tsin(spos.GetPhi()).
-sin(spos.GetTheta())):

// calculates gradient in cartesian coordinates

102

}

grad 8 dV,dr 8 e_r + (1/sin(spos.GetTheta())/spos.GetR())
8 dV_dphi 8 e-phi + (-sin(spos.GetTheta()))8
(dV-dcostheta/spos.GetR()) 8 e_theta;

return grad;

else
return grad; // grad is null vector in this control path

/¥#¥*#¥fitttfil##t*$flfi*¥t#¥**#*ttitttttttfitittttt\

e
e
e
e
e
e

e

Calculates the derivative of the potential

due to the eqn. of state in the r-direction at
a radial position r located in the box

i_r, i-phi, i-cos_theta. The derivative is
linearly interpolated depending on r.

Gets: - ints i_r. i-phi, i_cos_theta, double r
Returns: - derivative

\tltttttttttfiti.ttittttittttfifiittitfiitttttitttti/
double Star::CalculateInterpolatedDVDR(int i_r, int i_phi.

{

int i_cos_theta, double r)

// derivative

double dV_dr;

// potential in the neigboring box with greater r
double v_r_b;

// potential in current box

double v_r_c;

// potential in neighboring box with smaller r
double v_r_s;

// radial distance between center of current box
// and that of box with greater r

double dr_b;

ll radial distance between center of current box
// and that of box with smaller r

double dr_s;

// r value of smaller boundary of current box
double r-s;

// r value of greater boundary of current box
double r_b;

// calculates/sets values as described above
v_r,c 8 potmap.GetValueInBox(i_r. i_phi, i_cos_theta);
r_b 8 potmap.GetRMax() 8 potmap.P-R_Inv((double)(i_r+1)/
(potmap.GetIRMax())):
dr_b 8 potmap.GetRHax() 8 (potmap.F_R_Inv((double)(i-r+1.5)/
(potmap.GetIRMax())) - potmap.F_R_Inv((double)(i_r+0.5)/
(potmap.GetIRHa1())));
dr_s 8 potmap.GetRHax() 8 (potmap.F-R_Inv((double)(i_r+0.5)/
(potmap.GetIRHax()))
- potmap.F_R_Inv((double)(i_r-0.6>0?i_r-0.5:0)/
(potmap.GetIRHax()))):
r_s 8 potmap.GetRMax() 8 potmap.P-R_Inv((double)i_r/
(potmap.GetIRHax())):

// sets v_r_b equal to DENS_HIN if current box is outermost one
// (assuming mass density is low there)
v-r_b 8 i_r<potmap.GetIRHax()-1?

potmap.GetValueInBox(i-r+1, i_phi. i_cos_theta):DENS_MIH;

// sets v_r_s equal to v-r-c if current box is innermost one
II this sets derivative at center to zero
v_r_s 8 i_r>0? potmap.GetValueInBox(i-r-1, i_phi, i-cos_theta):v_r-c;

// linear interploation between the two derivatives

dV_dr 8 (r-r-s)/(r-b-r-s)8(v_r_b-v_r_c)/dr_b + (r_b-r)/(r_b-r-s)8
(v_r_c-v_r_s)/dr-s;

103

“,..!

I“. H
L.

}

return dV_dr;

/88#8888888888888888888888888888888888888888888\
Calculates the derivative of the potential

due to the eqn. of state in the phi-direction at
an angular position phi located in the box

i_r. i_phi, i_cos,theta. The derivative is
linearly interpolated depending on phi.

Gets: - ints i_r, i_phi, i,cos-theta, double phi

‘I‘I'I'ii‘l

*

Returns:

- derivative

\88888888888888888888888888888888888888888888888]
double Star::CalculateInterpolatedDVDPhi(int i_r, int i_phi,
int i_cos_theta. double phi)

{

}

double
double
double
double
double
double
double

dV_dphi;// derivative

v-phi-b;// potential in the neigboring box with greater phi
v,phi_c;// potential in current box

v-phi_s;// potential in neighboring box with smaller phi
dphi;// angular distance between neighboring boxes

phi-s;// phi value of smaller boundary of current box
phi_b;// phi value if greater boundary of current box

// calculates values as described above

v-phi_c 8 potmap.GetValueInBox(i_r, i_phi. i_cos_theta);
phi-s = (double)i_phi/potmap.GetIPhiMax()828PI;

phi_b 8 (double)(i_phi+1)/potmap.GetIPhiMax()828PI;

dphi 8 2.08PI/(double)potmap.GetIPhiHax();

// makes sure i-phi is in correct boundaries
v_phi_b 8 i_phi<potmap.GetIPhiMax()-1?

potmap.GetValueInBox(i_r, i-phi+1, i_cos_theta):
potmap.GetValueInBox(i_r, O , i_cos_theta);

v-phi_s 8 i_phi>0?potmap.GetValueInBox(i_r. i_phi-l, i_cos_theta):

potmap.GetValueInBox(i-r, potmap.GetIPhiMax()-1. i_cos_theta);

// linearly interpolates between the two derivatives at the boundaries
dV_dphi 8 (phi-phi_s)/(phi_b-phi_s)8(v-phi-b-v-phi_c)/dphi +

(phi-b-phi)/(phi_b-phi-s)8(v,phi_c-v_phi-s)/dphi;

return dV-dphi;

letteeeeeeeteeteeseeeeeeteeeseeeeeeeeeeeeeeeeee\

Calculates the derivative of the potential

due to the eqn. of state in the "cos(theta)-direction" at
a position costheta located in the box

i_r, i_phi. i_cos_theta. The derivative is

linearly interpolated depending on cos(theta).

.GC’OI'.

*

Gets:
Returns:

- ints

i_r. i_phi, i_cos_theta. double costheta

- derivative

\eeeeseeteeeeeeeoeeeeeeateeeeeseeeeeeeeeeeeeeete/
double Star::CalculateInterpolatedDVDCosTheta(int i_r. int i_phi,
int i-cos_theta, double costheta)

{

double
double
double
double
double
double
double

dV_dct;// derivative

v-ct_b;// potential in the neigboring box with box with greater r
v,ct_c;// potential in current box

v_ct_s;// potential in neighboring box with smaller r

dct;// difference of cos(theta) for two neighbouring boxes
ct-s;// cos(theta) value of smaller boundary of current box
ct_b:// cos(theta) value of greater boundary of current box

// calculates derivative only if there are separations at all
if(potmap.GetICosThetaHax()>1)

{

// sets values as described above
v_ct-c 8 potmap.GetValueInBox(i-r, i_phi. i_cos-theta);

104

dct 8 2.0/(double)potmap.GetICosThetaMax();
ct_s 8 2.0 8 (double)i_cos_theta/potmap.GetICosThetaMax()-1.0;
ct-b 8 2.0 8 (double)(i_cos_theta+1)/potmap.GetICosThetaHax()-1.0;

// makes sure i_cos_theta is in correct boundaries.

// if current box is at a boundary the potential

// for the "missing" next box is set to the value

// of the potential in the current box

// HIGHLY ARBITRARY!!

v_ct_b 8 i_cos_theta<potmap.GetICosThetaMax()-1?
potmap.GetValueInBox(i_r, i_phi, i_cos_theta+1):
v_ct_c;

v_ct_s 8 i_cos-theta>0?
potmap.GetValueInBox(i_r, i,phi, i_cos_theta-l):
v-ct_c;

// linearly interpolates between the two derivatives

// at the boundaries

dV_dct 8 (costheta-ct-s)/(ct_b-ct-s)8(v_ct_b-v_ct-c)/dct +
(ct_b-costheta)/(ct_b-ct_s)8(v_ct_c-v_ct_s)/dct;

else
dV_dct80;

return dV_dct;

}

It8.88888888888888888888888888888888888888888888888888¢\
8 Refreshs the information stored in tp_id[] and

8 id-tp[].

8 Gets: nothing

8 Returns: nothing
\888888888888888888888888888888888888888888888888888888/
void Star::RefreshTP()

{
int i=0;
for(i=0; i<number_tp; 18+)
{
tp_id[i] a tp[i].GetID();
id_tp[tpEi].GetID()] 8 i:
}
return; .
}

[eeeeeeeeeeeeeeeeeeeeeeeteeeteeeeeoeeeeeeeeeeeeeeeeteeeee\
8 saves the coordinates of every
8 (number_tp/points_per_step)th test particle in a file
8 ("suno_all.dat") readable by output.exe
8 Needs: - function Star::RefreshTP()
8 Gets: - int points-per-step. i.e. no. of points to be
8 saved per time step
8 Returns: nothing
\88888888888888888888888888888888888888888888888.88888888/
void Star::SaveCoordinatesVB(int points_per_step)
{
int j.k; // loop counters
int skip; // number of points to be skipped

// Update arrays tp_id and id,tp
RefreshTP();

// open output file
ofstream sunoall ("suno_all.dat", ios::app);

// calculates skip

105

skip 8 number-tp/points_per-step ? number_tp/points_per_step : 1;

// save coordinates in output file
if (sunoall.is_open())
for(j80; j<number-tp; j+8skip)

{
for(k80; k<3; k++)
{
sunoall << tp[id_tp[j]].GetPos().GetCoord(k) << endl;
}
}

// set end-of-timestep mark
sunoall << "end of timestep\n";

// close output file
sunoall.close();

// the following is just for error diagnosis
// begins here///////////////////////////////////////
ofstream denstp ("dens-tp.dat", ios::app);
if (denstp.is_open())
for(j80; j<38skip; j+8skip)

{
Spherical spos;
spos8tp[id_tp[j]].GetPos().GetSpherical();
int i_r8dmap.GetIR(spos);
if(i_r<dmap.GetIRMax())
{
denstp << "Particle 8" << j << ": Density: " <<
dmap.GetValueInBox(i_r.dmap.GetIPhi(spos).
dmap.GetICosTheta(spos)) <<
", Force due to 808:" <<
GetGradEOSPotential(id_tp[j]).GetNorm() <<
", Force due to gravity:" <<
G8m_tp8m_tp8id_tp[j]/
tp[id-tp[j]].GetPos().GetNorquuared()
<< endl;
}
}

denstp.close();
// ends here/////////////////////////////////////////

return;

[sessseeeeeeeseeeseeesesseeeesessssssessetseesstsssssesss\
Saves all the positions and momenta of all test
particles inside a sphere of radius max_radius to
a binary data file. Other parameters such as the
number of test particles saved, the total mass of the
saved star etc. are alos saved.
Gets: - max_radius as double

- timestep as int (used for filename)
Returns: nothing
\88888888888888888888888888888888888888888888888888888888]
void Star::SaveAllData(double max_radius. int timestep)

{

{ififiifl‘li

double tempBO; // for temporary storage

int max-i80; // number of last particle saved

strstream s;

s << timestep << ends;

string timestepstring 8 s.str(); // convert timestep to string
string filenamestr 8 "suno-full-data";

string filenameext 8 ".dat“;

filenamestr 8 filenamestr 8 timestepstring + filenameext;
const char8 filename 8 filenamestr.c_str();

106

}

ofstream out_bin (filename, ios::binary);

// I’m not sure if this is necessary but t/ps are sorted
SortTP();

// finds last test particle within sphere of radius max_radius

for(max_i80; tp[max-i].GetPos().GetNorm() < max_radius It
max_i<number_tp; max_i++);

cout << "Number of t/ps saved: “ << max-i << endl;

// calculates the total mass of the saved particles
double new_m_star 8 m_tp 8 max_i;

// save all the (more or less...) relevant parameters
out-bin.write(tmax_i, sizeof(int));
out_bin.write(tmax_radius, sizeof(double));
out_bin.write(&new_m_star, sizeof(double));
out_bin.write(&m_tp, sizeof(double));
/8out_bin.write(&GAMMA, sizeof(double));
out_bin.write(tY_E, sizeof(double));
out-bin.write(&flag_thermalize, sizeof(bool));
out_bin.write(tflag_average-theta, sizeof(bool));*/

// save all the positions and momenta
for(int i=0; i<number_tp; i++)

{
for(int j80; j<3; j+8){
temp 8 tp[i].GetPos().GetCoord(j);
out_bin.write(&temp, sizeof(double));
}
for(int j80; j<3; j++){
temp 8 tp[i].GetHom().GetCoord(j);
out_bin.write(&temp, sizeof(double));
}
}

// close filestream
out_bin.close();

return;

/tttt#¥tfifittfitttttfitttt‘tfi.it.iii!!!ttttttltfitttfififitfitttt\

.
*
$
.
Q

Reads the whole data (positions, momenta and parameters)
saved with Star::SaveAllData() for use as initial
conditions of a new simulation.

Gets: nothing

Returns: nothing

\.*.$..#‘..‘$¥#¥.I¥¥ti¥¥fit.fififififittfifitfifi##fitfifitttifiittfittt/
void Star::ReadAllData()

{

double temp80; // temporary
Vector vec; // also just temporary
ifstream in_bin (“suno-full_data.dat", ios::binary);

// read all the (more or less...) relevant parameters
in_bin.read(tnumber_tp, sizeof(int));
in_bin.read(kradius, sizeof(double));
in_bin.read(&m-star, sizeof(double));

in_bin.read(tm_tp, sizeof(double));

/8 in_bin.read(&GAHHA. sizeof(double));
in-bin.read(tY-B, sizeof(double));
in_bin.read(&f1ag_thermalize, sizeof(bool));
in_bin.read(tflag_average-theta. sizeof(bool)):*/

// read all the positions and moments

for(int i80; i<number_tp; 1+8)
{

107

for(int j80; j<3; j++)

{
in_bin.read(&temp, sizeof(double));
vec.SetCoord(j. temp);

}

tp[i].SetPos(vec);

for(int j80; j<3; j++)

{
in_bin.read(&temp, sizeof(double));
vec.SetCoord(j. temp);

}

tp[i].SetMom(vec);

}

// close file stream
in_bin.close();

return;

}

/888e8888888.88888888888888888888888888888888888888888888\
8 Saves a slice through the density map in the x-z-plane

8 to "suno_dens.dat". This file is readable by "SuNo

8 Density 0utput.exe".

8 Gets: nothing

8 Returns: nothing
\88888888888888888888888888888888888888888888888888888888/
void Star::SaveTPHDensityVB()

{

int i-r80, i_cos_theta80, i=0; // loop counters
ofstream tpmdensity("suno_dens.dat", ios::app);

tpmdensity << dmap.GetRHax() << endl;
tpmdensity << time_passed << endl;
for(i_r80; i,r<dmap.GetIRMax(); i_r++)
for(i=0; i<2; i++)
for(i_cos-theta80; i_cos-theta<dmap.GetICosThetaMax();
i_cos-theta++)
tpmdensity << dmap.GetValueInBox(i_r.
i8dmap.GetIPhiHax()/2,i_cos_theta)
<< endl;

tpmdensity.close();

return;

}

[eeeeeeeeeseeeeeeeeeeeeeeeeteeeteeeseeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeseeeeeeee\
8 Compares the distance from the origin of two test particles. This function
8 is needed to use qsort() from stdlib.h.
8 Gets: - pointers to the two test particles to be compared
8 Returns: - 1 if a is farther away from origin than b, -1 otherwise
\8888888888888888888888888888888888888888888888888888888888888888888888888888888]
int compare_tp_dist(const void8 a. const void8 b)
{

return ((TestParticle8)a)->GetPos().GetNorquuared() >

((TestParticle8)b)->GetPos().GetNorquuared()? 1 : -1;

[8888888888888888888888888888888888888888888888888888888888888888888888888888888\
8 Sorts the test particles by their distance from the origin.

8 Needs: stdlib.h. global function compare_tp_dist

8 Gets: nothing

8 Returns: nothing
\‘ttttfitlttttttttttttttttfittitttttttttttttttfitttttttttlttttttttttttltttfitittitti/
void Star::SortTP()

108

qsort(tp, number_tp. sizeof(TestParticle),compare_tp_dist);

/#¥#t*#*$***fittitfi‘lfitt*i‘tfiittittitfitittifit##fittttttttt*#\

iii...

Calculates & returns the time derivatives of test paricle
i’s position and momentum, those are its velocity

and the force on it.

TEST PARTICLES MUST BE SORTED!

Gets: - number of test particle i

Returns: - derivatives as a TestParticle object

ttfitttfi*t#¥#*$##*‘¥#¥#¥fitfiti##tfitlfifiiiififiilfifittiitttttfitill
TPChange Star::Deriv(int i)

{

TPChange derivs;

// stores derivatives of position and momentum of particle i in derivs
derivs.SetPosChg(tp[i].GetHom()8(1/
sqrt(m_tp8m_tp+tp[i].GetHom().GetNorquuared()/(c8c))));
// the long expression above is the relativistic velocity
// of test particle i
derivs.SetHomChg(GetForce0nTP(i));

return derivs;

Itit...ttttttttttttttttitttttitttt0...!!!it.it!!!##0##.tlfittlttttitfittlttifiit¢t¢\
8 Calculates the next step in the star’s time development using the Euler

8 method.

8 Uses function Star::Deriv().

8 Gets: - time stepsize dt for the step
\ttt‘t‘it‘tit‘ti‘##tfit##‘ttfitttfiittfitttttttt##8##.titit##*¥¥#*##tfittttttttfittttt/
void Star::NextStepEuler(double dt)

{

}

int i=0;
TestParticle8 initial;

initial 8 new TestParticle [number_tp];

// sorts the test particles
SortTP();
RefreshTP();

for(i=0; i<number_tp; i++)
initialEi] 8 tpEi]:

for(i80; i<number_tp; i++)
{
tp[i] 8 initialEi] 8 dt8Deriv(i);
/8 tp[i].SetMom(initial[i].GetHom() 8 Deriv(i).GetHom() 8 dt);
tp[i].SetPos(initial[i].GetPos() 8 Deriv(i).GetPos() 8 dt);8/
}

delete[] initial;

return;

ltfittitfifitfitfi$tfitttli.ttitttltfittl##t‘ltfittt##ttttttfittfilttttifitfittttttt\

.

Calculates number_of-steps steps with stepsize dt. Density of test

8 particles is written into a file ("suno_out.txt").

*

5....

Coordinates of (some) test particles are saved using SaveCoordinates().

Gets: - time stepsize dt for the step

- number of steps as int

- number of points to be saved to the output file by function

Star::SaveCoordinatesVB() per time step as int

109

points_per_step_save
- bool modify_stepsize which indicates if the stepsize may be
modified during the calculation
Returns: nothing

*I'filfi

soease888888888888888888888eeeeeeeseseeeeseeeeeeeeeeeeeeeeeeeeeeeeseesee/

void Star::Stepper(double dt, int number_of_steps, int points-per_step_save,
bool modify_stepsize. bool modify_stepsize_avvel)

{
int i=0;
bool stepsize_modified8false; // indicates if stepsize already modified
bool stepsize_modified_twice8false;
double av,disp=0; // average position change of t/ps in last step

// deletes old file created with Star::SaveCoordinatesVB() and

// writes the number of time steps to the 1st line of the output file
// open output file

ofstream sunoall ("suno-all.dat");

// write number of stime steps to file
//(+1 because the initial configuration is also saved)
sunoall << number_of_steps + 1 << endl;

// close output stream
sunoa11.close();

// deletes old file created with Star::SaveTPMDensityVB() and
// writes the no. of time steps, i-r-max, i_cos_theta-max and
// F-R_Inv(i) to the new output file

ofstream tpmdensity ("suno_dens.dat");

tpmdensity << number_of_steps + 1 << endl;
tpmdensity << dmap.GetIRMax() << endl;
tpmdensity << dmap.GetICosThetaHax() << endl;
for(i=0; i<dmap.GetIRNax(); i++)
tpmdensity << dmap.F_R-Inv((double)i/dmap.GetIRHax()) << endl;

tpmdensity.close();

// files to save energies and angular momentum
ofstream sunoenergy ("suno-energy.dat");
ofstream sunoangmom ("suno_angmom.dat");

// the following is just for dignosis

// begins here//////////////////////////////////
// deletes old dens_tp.dat - file

ofstream denstp ("dens_tp.dat");

denstp.close();

// ends here///////////////////////////////////

SortTP():
CalculateTPHDensityMap():
SaveTPHDensityVB();

// saves initial configuration
// SaveCoordinates(i);
SaveCoordinatesVB(points_per_step-save):

// set internal energies for test particles if "energy method" is used
if(flag-energy_method)
InitializeInternalEnergies();

for(i=0; i<number-of-steps; i++)
{
// sets new time
time-passed +8 dt;

// unnecessary output

110

}

cout << "Calculating time step " << i+1 << ”.\n";
if(modify_stepsize it !stepsize_modified-twice)
dt 8 HodifyStepsize(dt, stepsize-modified,
stepsize-modified_twice);
av_disp 8 NextStepRK4(dt. modify-stepsize_avvel);
// enforce energy conservation by modifying momentum
// so that the change in kinetic energy compensates
// the change in internal energy
if(flag-energy_method)
HodifyMomentaConservingEnergy();

if(i>5 it modify_stepsize_avvel)
dt 8 HodifyStepsizeAvVel(dt, av,disp);

SaveCoordinatesVB(points_per_step_save);
SaveTPHDensityVB();

// save energies

sunoenergy << time_passed << "\t" << HeasureKineticEnergy(true)
+ MeasureGravitationalEnergy() + MeasureInternalEnergy() <<
"\t " << MeasureKineticEnergy(true) << “\t" <<
HeasureGravitationalEnergy() <<
"\t" << MeasureInternalEnergy() << endl;

// save angular momentum
sunoangmom << time_passed << "\t" <<
MeasureTotalAngularHomentum().GetNorm() << endl;
}

// close file streams
sunoenergy.close();

sunoangmom.close();

return;

/*#$fi#iittttttfiitttttfittiitifiiitttfiiifilifitttttififiltttti*¥*#l*#\

.
i
I
t
t

t

Modifies the stepsize used by Star::Stepper(). Checks the
average displacement of all t/ps in last time step. If too
big, stepsize is divided by 2.

Gets: - current stepsize dt as double

- average displacement av_disp as double

Returns: - new stepsize as double

\8888888888888888888888888888888888888888888888888888888888888/
double Star::HodifyStepsizeAvVel(double dt, double av_disp)

{

}

if(av_disp>radius/6e2)

{
dt /8 2.0;
cout << "Step size divided by 2.\n";
cout << "Average displacement of t/ps was: " << av_disp << "m.\n";
}
else
if(av_disp<radius/6e3)
{
dt 88 2.0;
cout << "Step size multiplied by 2.\n";
cout << "Average displacement of tips was: " <<
av_disp << "m.\n";
}

return dt;

leeeeeeeeeeeeseeeeteeee88888888888888eeeeeeeeeeeeeeeeeeeeeeeee\
8 Modifies the stepsize used by Star::Stepper(). Stepsize is

8 divided by 10 if 70% of the test particles are in a shell

8 with 1/4 the radius of the initial star.

111

8 Gets: - currently used stepsize dt as double
8 Returns: - new stepsize as double
\eee8888888888888888888888888888888888888888888888888888888888/

double Star::HodifyStepsize(double dt, boolt stepsize_modified.

{

}

boolk stepsize_modified_twice)

double new_dt80;

// check if 70% of t/ps are in the 1/4 - radius - shell
if(!stepsize_modified)
new_dt8 (stepsize_modified 8
(tp[(int)(0.7 8 number_tp)].GetPos().GetNorquuared()
< radius8radius/16)) ? dt/10 : dt;
/e
// check if 70% of t/ps are in the 1/5 - radius - shell
if(Estepsize_modified)
new_dt8 (stepsize_modified 8
(tp[(int)(0.7 8 number-tp)].GetPos().GetNorquuared()
< radius8radius/25)) ? dt/lO : dt;
8/
else
// check if 70% of t/ps are in a 1/25 - radius - shell
if(!stepsize_modified_twice)
new_dt8 (stepsize,modified-twice 8
(tp[(int)(0.7 8 number_tp)].GetPos().GetNorquuared()
< radius8radius/625)) ? dt/20 : dt;

return new_dt;

[##tttfitififitttttfifiitfit*ltttfittifitfitlttil‘tt‘t‘ttfififittfifitlt*##fi#\

#
a
e
t
0

mass density as a function of distance from the center
(initial condition)

- distance r from center as double
- parameter_d as double

Returns: - density at that distance

\*##‘ll.##ttittfififillfittitfiififilittfi##lttttttltfitltfitt.##**#*¥###/

{

}

double Star::initial_mass_density(double r, double parameter_d)
double rho; // mass density

// realistic
//rho 8 1.0/(4.08PI/3.08parameter_d8pow(r810.3)+1.005e-13);

// realistic inner core 0.7 H_\odot contracted by 5
rho 8 pow(5.0,3)/(8.2436e-308pow(r86,3)+1.005e-13);

// homogeneous
//rho 8 m_star/4.0/pow(radius,3);

// 'bonazzola fig. 1 c
l/rho 8 2.4e17 8 0.08 8 exp(-8e-108pow(r,2));

return rho;

[tilttittitfiilfilltiltfifitfifitfitfililittt*fitfitfittfiifitttttllflitfittt\

§ C'il i Q l

Sets test particle tp[i] to a random position in a sphere
of radius r around the star’s center.

Random number generator should be initialized before
calling this function.

Gets: - int i identifying the t/p

- radius r of sphere as double

Returns: - true if an error occurs. false otherwise

\8888888888888888888888888888888888888888888888888888888888888/
bool Star::PutTPInSphere(int 1, double r)

112

Vector vec; // for temporary use

// check if argument values are ok
if(i<0)
{
cout << "Error: argument i<0 in function PutTPInSphere().\n";
return true;
}
else
if(i>8number_tp)
{
cout << "Error: argument i>=number_tp in function
PutTPInSphere().\n";
return true;

}
else
if(r <8 0)
{
cout << "Error: argument r<80 in function PutTPInSphere().\n";
return true;
}

// finds random vector
do
for(int j80; j<3; j++)
vec.SetCoord(j, 2.08rand()/RAND_HAX-1.0);
while(vec.GetNorquuared()>1);//rejection method(see Numerical Recipies)

// sets test particle i to the random position
tp[i].SetPos(vec8r);

return false;

}

lttttttttttittfittttttittt‘fifitittfit##titfi¥¢##$#ttttfittttttt#####\
Distributes test particles in position space. Uses function
Star::initial,mass-density(). normalizes it (so that it

gives the correct total mass of the star) and distributes

test particles accordingly.

Needs: - number of shells with different density as int
Returns: nothing
\88888888888888888888888888888888888888888888888888888888888888/
void Star::DistributeTPPos(int no_shells)

{

§§§§§§

// shell thickness

double delta_r80;

// array containing radii of spheres to be filled

double r[no-shells];

// mass density for current shell

double rh080;

// array containig the numbers of test particles distributed
int no_tps_sphere[no_shells];

int no_tps_current_sphere80;

// total number of t/ps already distributed

int tps-dist80;

// bulk mass of core resulting from density distribution
double bulk_mass80;

// parameter for initial density distribution

double par-d81.5e-30;

// unnecessary output
cout << "Positionz";

// initializes random generator
srand(time(NULL)):

// calculates shell thickness

113

delta_r 8 radius / no_shells;

// calculate the bulk mass resulting from initial_mass_density()
// in order to normalize density distribution so that m_star comes
// out correctly. also vary parameter par_d to get correct mass
//do{

bulk_mass 8 0;

for(int k80; k<no_shells; k++)

{
r[k] 8 radius-delta-r8k;
rho 8 initial_mass_density(r[k] - delta-r/2, par_d);
bulk_mass +8 rh084.08PI/3.08(pow(r[k],3)-pow(r[k]-
delta_r,3));
}

//par_d 880.0001e-30;
//}while(bulk,mass>m_star):

cout << "bulk mass: " << bulk_mass << endl;

// assuming that density increases when going inward, distribute

// test particles in shells (starting outside) in a way that

// creates the correct densities

for(int k80; k<no_shells; k++)

{

// get desired mass density for current shell

rho 8 m_star/bulk-mass 8 initial_mass_density(r[k] -
delta_r/2, par-d);

// calculate how many test particles are needed in

// current sphere (not shell!!)

// always too small due to round-off

no_tps_current-sphere 8 (int)(rho84.0/3.08PI8pow(r[k],3)/m_tp);

no_tps_current_sphere 8 no_tps_current-sphere<number_tp?
no_tps-current_sphere:number_tp;

// subtract the (approximate) number of particles that

// is already in the sphere

for(int n80; k-n>0; n+8)

{
no_tps_current_sphere -8 (int)(no-tps_sphere[n]8

pow(r[k]/r[n].3)+1); llalways subtract to much

}

no-tps_current_sphere 8 no-tps_current-sphere>0?
no_tps_current_sphere:0;

no_tps_sphere[k] 8 no_tps_current_sphere;

// finally put the corrected number of test particles
// in the sphere
for(int j80; j<no-tps_current-sphere; j++)
{
if(tps_dist+j>8number-tp)
{
cout << "Ran out of test particles
while creating initial position space configuration.\n";
break;
}
PutTPInSphere(tps_dist+j, r[k]);

// output of progress (unnecessary)
if((tps_dist+j)X(number_tp/10)880)
cout << ".";
}
// count the distributed test particles
tps_dist+8no-tps_current_sphere;

}
cout << "\nResidual test particles: “ << number_tp-tps_dist << endl;

// distribute residual test particles
// completely arbitrary but shouldn’t matter

114

for(int j8tps_dist; j<number_tp; j++)
PutTPInSphere(j. radius/3);

return;

}

/*#tt##t*#*#t¥tttfitttfiti'tittititfitttttlttttttttttttfitttlfitttfifittifitttttt*tttfit*$\

 

8 Distributes test particles in momentum space imposing a "rigid body rotation“
8 at angular velocity omega if diff_rot 8 false.

8 Distributes test particles in momentum space using a radial

8 dependence of the angular velocity

8 omega_0 8 r_0‘2

8 omega 8

8 r‘2 + r_0‘2

8 if diff_rot 8 true.

i

8 Gets: - Vector omega (corresponds to omega_0 in the above formula)
8 - r_0 as double

8 - diff_rot as bool

8 Returns: nothing

\eseseeeeesseessssssssssetssets-seeeeeeeeeeeeseeeesseeesseeeeeeeseeestestesstses/
void Star::DistributeTPMom(Vector omega, double r_0, bool diff_rot)
{

Vector v; // for temporary storage / local velocity

double max80; // maximum velocity squared

double vc80;

// unnecessary output of rotation type used
switch(diff_rot)
{
case true:
cout << "Using differential rotation with r_0 8 "
<< r_O << "m.\n";
break;
default:
cout << "Using rigid body rotation.\n”;
}

// unnecessary output
cout << "Momentumz";

for(int i=0; i<number-tp; i++)
{
switch(diff_rot)
{
case true:
v 8 r_08r_0/(tp[i].GetPos().GetHorquuared() +
r_08r_0) 8 (omega 8 tp[i].GetPos());
break;
default:
v 8 omega 8 tp[i].GetPos();
}
//v 8 v + (-1) 8 tp[i].GetPos() 85;// infall
if((vc8v.GetNorquuared())>max)max8vc;
// set relativistic momentum
tp[i].SetHom(m_tp 8 v 8 (1/sqrt(1-(v‘v)/c/c)));
// output of progress (unnecessary)
if(iX(number_tp/10)=80)
cout << ".";

}

// Warning if velocities are close to c (or greater than c!)
if(max>0.818c8c at max < c8c) cout

<< "Warning: initial velocities greater than 0.9c.\n”;
if(max > C8c) cout

<< "Warning: initial velocity greater than c found!\n";

cout << endl;

115

return;

}

[essessseeessseseeseeesassttestes:assessesssseetesseesesestsesesssssesse\
8 Gives the t/ps a moomentum towards the center of the core. This

8 momentum increases linearly with the t/p’s distance from the center.

8 Gets: - double max_velocity indicating the maximum velocity (for the

8 outermost t/p

8 Returns: nothing
\ssssseeeeesseesseeesassessassesseseesesssesstsmsssssesssessesetesssssee/
void Star::KickTPInside(double max-velocity)

{
Vector newmom;
for(int i=0; i<number_tp; i++)
{
newmom 8 tp[i].GetMom() + tp[i].GetPos() 8
(-1/radius) 8 m_tp 8 max_velocity;
tp[i].SetMom(newmom);
}
return;
}

[88888888888888888888888888888888888888888888888888888888888888888888888\
8 Calculates the next step in the star’s time development using a
8 fourth order Runge-Kutta-Algorithm

e
8 Gets: - time stepsize dt for the step

8 - flag modify_stepsize_avvel as bool indicating if if the
8 average velocity of all t/ps shall be calculated

\88888888888888eeeeeeeseeeeeee888888888888888888888esteeeeeeeeeeeeeee888/
double Star::NextStepRK4(double dt, bool modify_stepsize_avvel)
{
int i;
// for saving the initial configuration
TestParticle 8initial;
// changes of test particle momentum and positions
TPChange 8k1, 8k2, 8k3, 8k4. tota1-k;
// average position change of all test particles in this step
double av_pos_chg80;

// memory reservation

initial 8 new TestParticleEnumber_tp];
k1 8 new TPChangeEnumber_tp];

k2 8 new TPChangeEnumber_tp];

k3 8 new TPChangeEnumber_tp];

k4 8 new TPChangeEnumber_tp];

// 1st step
CalculateK(k1,dt);
// copies initial configuration
// (after 1st step because new test particles are sorted)
for (i=0 ;i<number-tp ;i++)
initialEtp_id[i]] 8 tp[i];
// moves system to 1st intermediate position
HoveTPs(initia1. k1, 0.5);

// 2nd step

CalculateK(k2,dt);

// moves system to new intermediate position
HoveTPs(initial, k2, 0.6):

// 3rd step

CalculateK(k3 . dt) ;

// moves system to final position
HoveTPs(initial, k3, 1);

116

// 4th step and performing overall step

CalculateK(k4,dt);
for (i=0; i<number_tp ;i++)
{

total-k 8 k1[i]‘(1.0/6.0) 8 k2[i]*(1.0/3.0)
+ k3[i]8(1.0/3.0) + k4[i]8(1.0/6.0);
tp[id_tp[i]] 8 initialEi] + tota1_k;
// calculate average displacement only if needed
if(modify_stepsize_avvel)
av,pos-chg +8 total_k.GetPosChg().GetNorm();
}

// free memory
delete[] initial;
delete[] k1;
delete[] k2;
delete[] k3;
delete[] k4;

// normalize av_pos_chg
av-pos_chg /8 number_tp;

return av_pos-chg;

}

[888888888888888888888888888888888888888888888888888888888888888\
8 Calculates the current TPChange for all test particles and

8 saves it in the TPChange array k. dt is the step size.

8 Gets: - a pointer to an array of number-tp TPChange objects

8 - desired step size

8 Returns: nothing. but changes 8k
\ttttttttttttfitt‘#0.!tit...‘ttttittittltilfittltttittittttttttttt/
void Star::CalculateK(TPChange8 k, double dt)

{
int i=0;
SortTP():
RefreshTP():
CalculateEDSPotentialMap();
for (i=0 ;i<number-tp ;i++)
k[tp-id[i]] 8 Deriv(i) 8 dt; // calculates k
return;
}

ltfitififiltttltttttltfitttttfifittfifi‘l‘ltfitttttfifitittfiit#tttttfi¥¥###‘#\

8 Moves all of the star’s test particles to 8initial+8k 8 kfactor
8 Gets: - a pointer to an array of number_tp TestParticles initial
8 - a pointer to an array of number_tp TPChanges k

8 - a double kfactor

8 Returns: nothing, but changes 8tp

\8888888888888888888888888888888888888888888888888888888888888888/
void Star::HoveTPs(TestParticle8 initial. TPChange8 k, double kfactor)
{

int i;

for (i=0 ;i<number_tp ;i++)
tp[id-tp[i]] 8 initialEi] + k[i]8kfactor;

return;

}

/ssessessessseeeeeeeeeeestseseeeeeseeeeeeeesseeeeeeeeeeees¢\
8 Modifies the momenta of all test particles to enforce

8 energy conservation. Should be called after every time

8 step. Previous internal energy of the test particle is

8 compared to its current internal energy. The difference

117

is subtracted from its kinetic energy so that also
the total angular momentum is conserved.

e
t
8 Gets: nothing
t

Returns: nothing
\eeeeeeeeeeseeeeseseestateseeeseseeeeeeeseeeeeeeeeeeeeeseet/

void Star::HodifyMomentaConservingEnergy()

{

Spherical spos;
int i_r80, i_phi80. i_cos_theta80;
double e_kin_change80, momentum-change80,

p80. p_r80, newenergy80, temp80. real_e_kin_change80;

Vector newmom, e-r;

CalculateEOSPotentialMap(); // might be obsolete
for (int i=0 ;i<number_tp ;i++)

{

}

return;

}

// converts position to spherical coords
spos 8 tp[i].GetPos().GetSpherical();
// r-direction unit vector
e-r 8 tp[i].GetPos() 8 (-1/spos.GetR());
// indices of box in which t/p is located
i_r 8 potmap.GetIR(spos);
i-phi 8 potmap.GetIPhi(spos);
i_cos,theta 8 potmap.GetICosTheta(spos);
if(i_r<potmap.GetIRHax()) // t/p inside map?

newenergy 8 potmap.GetValueInBox(i-r,i,phi,i-cos-theta);
else

newenergy 8 tp[i].GetInternalEnergy();//new energy88old energy
e_kin_change 8 newenergy-tpEi].GetInternalEnergy();
p8tp[i].GetMom().GetNorm(); // momentum of t/p

p_r8tp[i].GetHom()‘e_r; // radial proj. of momentum

temp8p_r8p,r-p8p+pow(e_kin-change/c+
sqrt(pow(m_tp8c,2)+p8p).2)-pow(m_tp8c,2);
if(temp>80)
{
momentum_change 8 p_r - sqrt(temp);
if(momentum_change < -p_r)
momentum_change 8 -p-r; // stop test particle

else

// arbitrary. doesn’t matter very much if p is very small
momentum-change 8 0;

1

// calculates new momentum

newmom 8 tp[i].GetHom() + momentum_change 8 e_r;
// sets new momentum

tp[i].SetHom(newmom);

real_e,kin_change 8 sqrt(pow(m_tp8c8c,2)+
newmom.GetNorquuared()8c8c)-
sqrt(pow(m_tp8c8c.2)+pow(p,2)8c8c);

// sets new internal energy for t/p, allowing "energy debt"

tp[i].SetInternalEnergy(newenergy-e-kin_change-real-e_kin_change);

//tp[i].SetInternalEnergy(newenergy); // no "energy debt"

[seasonsesseeseseseetosses:eeeeeeeeseeeeeeeeeesseeeeeesssse\
8 Sets the initial values for the internal energies of

8 all test particles. Should be called before the beginning
8 of the simulation.

118

8 Gets: nothing
8 Returns: nothing
\ttttttttttttttttttttttttttttttttttttfittt888*88888888888888/

void Star::InitializeInternalEnergies()

{

Spherical spos;

int i_r80. i-phi80, i_cos_theta80;

CalculateEOSPotentialMap();

for (int i=0 ;i<number_tp ;i++)

{
// get t/p’s position in spherical coordinates
spos 8 tp[i].GetPos().GetSpherical();
// indices of box in which t/p is located
i_r 8 potmap.GetIR(spos);
// make sure t/p is in map
i_r 8 i_r < potmap.GetIRHax() ? i-r : potmap.GetIRMax()-1;
i_phi 8 potmap.GetIPhi(spos);
i_cos-theta 8 potmap.GetICosTheta(spos);
// set internal energy value
tpEi].SetInternalEnergy(potmap.GetValueInBox(i_r,

i_phi,i_cos,theta));
}
return;
}

[seeseseseeeeeeseeseesseeesseeassess:eseeeseeeeeeeseseeeeee\
8 Calculates the TPSmearedNumberHap for the current configuration
8 of the star.
8 Gets: - a reference to a TPDoubleHap object
8 Returns: nothing
\eeseeseeeesseseeeeeseeestatesmesses:eeeeeeeeseeeeeseemsate]
void Star::CalculateTPSmearedNumberHap(TPDoubleHapt nmap)
{
int i=0;
Spherical spos;
double new_r-max; // new value for r_max, calculated below

//SortTP(); // might become obsolete

// sets r-max to 10.1/9.0 times the distance form the

// origin of the (number-tp/100)th most distant test

// particle (in order to include all "important" particles

// that do not escape due to fluctuations)

// r-max is conjectured to be smaller than 1.18radius

new_r-max 8 tp[(int)(0.998(number-tp-1))].GetPos().GetNorm()810.1/9.0;
nmap.SetRHax(new-r_max < 1.18radius ? new_r_max : 1.18radius);

// deletes entries in box-tp array if thermalization
// is activated
if(GetFlagThermalize())

ClearBoxTP();

// adds all test particles to the number map
for(i=0; i<number_tp; i++)
{
// convert tp’s position vector to spherical coordinates
spos 8 tp[i].GetPos().GetSpherical();
nmap.AddTPSmeared(spos): // nmap.AddTP() for unsmeared
// add test particle’s index to box-tp if
// thermalization is activated
if(GetFlagThermalize())
{
// add index to box-tp if test particle in map
int i_r8nmap.GetIR(spos);
if(i_r<nmap.GetIRMax())
box_tp[i_r][nmap.GetIPhi(spos)]
[nmap.GetICosTheta(spos)].push-back(i);

119

}

return;

}

[888888eeeeeeeeeeeee88888888888888easeeeeteeeeseeeeeeeeeee\
8 Calculates the mass density map.
8 Gets: - nothing
8 returns: nothing
\esseesseeeteessssssssesesessseeseeessteesteeeeseseeeesese/
void Star::CalculateTPHDensityMap()
{

TPDoubleMap nmap(dmap.GetIRMax(), dmap.GetIPhiHax(),

dmap.GetICosThetaMax());

double vol, dens;

int i_r80, i_phi80, i_cos_theta80;

double tpincsph80;

static const double RH0_180.0; //2.4e17 , density for thermalization

// calculates number map first
CalculateTPSmearedNumberHap(nmap);

// sets r-max to the value calculated by CalculateTPNumberHap()
dmap.SetRHax(nmap.GetRMax());

// set the density in the central sphere to an averaged value
// (independent of i_phi and i-theta)
vol 8 dmap.VolDfBox(0) 8 dmap.GetIPhiMax() 8
dmap.GetICosThetaMax(); // volume of central sphere
// counts all particles in central sphere
for(i_phi=0; i_phi<dmap.GetIPhiMax(); i_phi++)
for(i-cos,theta80; i_cos_theta<dmap.GetICosThetaMax(); i-cos_theta++)
tpincsph +8 nmap.GetValueInBox(0,i_phi,i_cos_theta);
// sets the averaged value
dens 8 m-tp 8 tpincsph/vol;
dens 8 dens > DENS_MIN ? dens : DENS-MIN;
for(i_phi80; i_phi<dmap.GetIPhiMax(); i-phi++)
for(i_cos_theta80; i_cos_theta<dmap.GetICosThetaHax(); i_cos_theta++)
dmap.SetValueInBox(0.i_phi,i-cos-theta, dens);

// calculates mass density map by dividing number map by volume of
// corresponding boxes and multiplying with mass of a test particle
// (i-r80 if seperations shall be made for the central sphere)
for(i_r81; i-r<dmap.GetIRHax(); i_r++)

{
vol 8 dmap.VolOfBox(i-r);
for(i_phi80; i_phi<dmap.GetIPhiHax(); i_phi++)
for(i-cos,theta80; i_cos-theta<dmap.GetICosThetaMax();
i_cos-theta++)
{
dens 8 m_tp 8 (double)nmap.GetValueInBox(i_r,
i_phi.i-cos_theta)/vol;
dens 8 dens > DENS_HIN ? dens : DENS_HIN;
dmap.SetValueInBox(i-r,i_phi.i-cos_theta, dens);
// thermalize test particles in box
// if density is greater than RHO_1
if(dmap.GetValueInBox(i-r. i-phi. i_cos_theta)
>RHO_1 It GetFlagThermalize())
ThermalizeBox(i_r.i_phi,i-cos-theta);
}
}

// averages the density values of the cos_theta-boxes (if flag is true)
I] which are mirror-symmetric about the equatorial plane
if(flag_average_theta)

dmap.AverageCosThetaBoxes();

120

return;

}

/8888eeeeeeeee88888888888888888888888888888888888888888888\
8 Calculates the EOS potential (=specific internal energy)
8 map.

8 Gets: - nothing

8 returns: nothing
\888888888888888888888888888888888888888888888888888888888/
void Star::CalculateEDSPotentialMap()

{

int i_r80, i_phi80, i,cos_theta80;

// calculates density map
CalculateTPHDensityMap();

// sets r_max to the value calculated by Ca1cu1ateTPNumberMap()
potmap.SetRHax(dmap.GetRHax());

// calculates EDS potential map
for(i-r80; i_r<potmap.GetIRHax(); i_r++)
for(i_phi80; i_phi<potmap.GetIPhiHaxC): i_phi++)
for(i_cos_theta80; i-c0s_theta<potmap.GetICosThetaMax();
i_cos,theta++)
{
potmap.SetValueInBox(i_r, i_phi. i_cos_theta.
CalculateEOSPotential(i_r, i_phi. i_cos_theta));

return;

}

[seeeeeeeesees88888eeeeeeeeeeeeeeeeeeeeeeeeseeeseeee\
8 Equation of state for the star. Calculates and
8 returns the internal energy (due to the EOS)
8 of a test particle (not a nucleon!) in a certain
8 box of the map.
8 Gets: - integers i_r, i-phi, i_cos-theta which
8 specify the box for which the potential
8 energy is calculated and returned
8 Returns: - potential energy in the box as double
\eassesseseeeeeeeeeeeeeasesseeeeeeeeeeeeeeeseeeeeeee/
double Star::CalculateEDSPotential(int i-r. int i_phi. int i_cos_theta)
{
// nucleons per test particle. neutron mass is used as nucleon mass
const static double nucleons_per_tp 8 m_tp/HASS-NEUTRDN;
// density in kg/m'3
double density 8 dmap.GetValueInBox(i_r, i_phi, i_cos-theta);
double potential80;
double y_e8Ca1culateYE(density);

// flag indicating if this is the first function call
static bool firstcall8true;

// transition density between the two EOS

static double density-at-gap 8 1e14;

static double potentia1_gap80;

double pressureEl];

// temperature in K

double temperature 8 CalculateTemperature(density);

// stuff for LS EOS
double 1s,out[4];

// stuff for Helmholtz EOS

integer ionmax82; // fortran int, number of different ions
double xmass[2], aion[2]. zion[2]; // composition of matter
xmass[0]82.69e-4; // low density value from cooperstein
xmass[1]81-xmass[o]; // should be made a function of density

121

aion[0]81;

aion[1]826/y_e;

zion[0]=1; // use hydrogen and
zion[1]826; // "iron" only

// calculate the gap at the transition density between the
// two EOS in first call
if(firstcall)
{
double temperature_at-gap 8 CalculateTemperature(density_at_gap);
lseos_ipvar[i_r][i_phi][i_cos-theta][0] 8
kelvin_mev(temperature_at_gap);

potential_gap 8 m_tp 8 CalculateInternalEnerngelmholtz(xmass.
aion, zion, ionmax, temperature_at_gap, density_at_gap, pressure)
- nucleons_per_tp 8
CalculateLSEOS(lseos-ipvar[i_r][i-phi][i-cos_theta],
y_e, density-at-gap. &(lseos_pprev[i_r][i_phi][i_cos_thetaJ),
ls_out);
firstcall8false;

}

if(density<density_at_gap)
{// use Helmholtz EDS
// contribution only from electron-positron-gas for
// densities below 1e14
potential 8 m_tp 8 CalculateInternalEnerngelmholtz(xmass.
aion, zion, ionmax, temperature, density. pressure);

}
else
{
// use Lattimer and Swesty EDS
// set temperature input variable
// temperature in HeV
lseos_ipvar[i-r][i-phi][i_cos_theta][0] 8 kelvin_mev(temperature);
// nuclear and electron-positron contribution at densities above 1e14
potential 8 nucleons_per_tp 8
CalculateLSEOS(lseos_ipvar[i_r][i-phi][i_cos-theta], y_e.
density. k(lseos_pprev[i-r][i-phi][i-cos_theta]),ls_out)
+ potential_gap;
}
[e

// nuclear contribution (old way) and polytrope EOS for electron gas
//static const double GAHHA81.33;

//static double K80.75 8 pow(PI/3.0,2.0/3.0) 8 HBAR 8 c

e pow(y_e/HASS_NEUTRON,4.0/3.0);

//static double E81.0/(GAHMA-1.0)8K;

const static double RH0_082.4e17; // nuclear density [kg/m‘3] 2.4e17

potential 8 nucleons_per_tp 8 (A 8 density / RHO_0
+ B 8 pow(density/RHO_0. SIGHA));

//potential +8 E8pow(dmap.GetValueInBox(i_r,
i-phi, i-c0s_theta).GAHHA-1.0) 8 m_tp;
8/

return potential;

}

[eeeeeeeeeeeeeeeeeeeeee88888888888888888888888888888\
Returns the force on test particle 1. Test
particles must be sorted. TPHDensityHap dmap must
be calculated for proper results!

Gets: - int i identifying test particle

Returns: - force Vector
\‘ttttttt‘ttttltttttttfitittttttttt¢tt¢ttlttt88888888]

{ii}!

122

Vector Star::GetForceOnTP(int i)

double r; // distance of test particle i from origin

// if test particles are unsorted. i must be replaced by CountTPInside(i)

if(i>number_tp/10000) // singularity treatment (very, very simple...)
f 8 tp[i].GetPos()8(-G8m_tp8m_tp8i/pow(r,3));

if(!flag_energy_method) // add only if "force method" is used
f 8 f - GetGradEDSPotential(i):

{
Vector f(0.0,0); // null vector
// force due to gravity
r 8 tpEi].GetPos().GetNorm();
// force due to equation of state
return f;

}

[seeeesseeeeesesseeeeeeeeesseeseesteases:sssessseseeeseeee\
Reads the parameters for the simulation from a

datafile (filename).

- a whole bunch of references to the parameters

to be read

t
l
8 Gets: - name of the datafile as char[]
0
a
e

Returns: - true, if an error occurs

8 - false, if all goes fine

\888888888888888888888888888888888888888888888888888888888/

bool read_parameters(char8 filename. intt number-tp, intt i-r_max,
intt i_phi_max, intt i-cos_theta_max, doublet radius, doublet mass,
doublet omega_norm, doublet stepsize. intt numb_part_save,
boolt modify_stepsize, boolt thermalize, boolt diff-rotation.
doublet r_0, boolt average_theta. boolt modify_stepsize_avvel,
doublet dens_min. boolt energy_method. doublet kick)

ifstream parameters (filename);

if(parameters)
{
parameters
parameters
parameters
parameters
parameters
parameters
parameters
parameters
parameters
parameters
parameters
parameters
parameters
parameters
parameters
parameters
parameters
parameters

else

)>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>

number_tp;
i_r_max;
i_phi-max;
i_cos_theta_max;
radius;
mass;
omega_norm;
stepsize;
numb-part-save;
modify_stepsize;
modify_stepsize-avvel;
thermalize;
diff-rotation;
r_O;
average-theta;
dens-min;
energy_method;
kick;

cout << "Error opening parameter file.\n";
return true;

}
parameters.close():

return false;

123

int main()

{

int number_tp81000, i_r-max840, i_phi_max81, i_cos-theta_max840;

double radius81.256e6, mass82e30, omega-norm82.1, stepsize=5e-4,
r_081.256e5;

double dens,min=1e1, kick80;

bool modify_stepsize8true, thermalize8false, diff_rotation8false.
average_theta8false;

bool modify_stepsize_avvel, new_simulation8true,
energy-method8false;

int steps=0; // number of steps

int numb_part_save82000; // number of particles saved in every time step

char8 parameter_file8"suno_parameters_new.dat“;

 

cout << "SuNo Version 0.3\n \n";
cout << "Use saved data (80) or start new simulation (81)? ";
cin >> new_simulation;

if(!new-simu1ation)
parameter_file8"suno_parameters_ctd.dat";

// reads the parameters from a data file

if(read_parameters(parameter_file, number-tp. i_r_max, i_phi,max,
i_cos_theta_max, radius. mass. omega_norm, stepsize, numb_part-save,
modify_stepsize, thermalize, diff_rotation, r_0. average_theta.
modify_stepsize-avvel. dens-min. energy_method. kick))
return 1; // exit with error code 1 if error in

//read-parameters() occurs

// creates a Star object (number of test particles. radius. mass, i-r_max,
// i_phi_max. i_cos_theta_max. thermalization flag)
Star star(number_tp. radius. mass. i_r-max, i-phi-max. i_cos_theta-max,
thermalize, average_theta, dens,min, energy_method);
Vector omega(0.0,omega_norm);//init. angular velocity of star

if(new_simulation)

{
// distribute test particles in position space
// (argument8number of shells)
star.DistributeTPPos(i_r_max);

// distribute test particles in momentum space
star.DistributeTPMom(omega, r-0, diff_rotation);

cout << "Total angular momentum: " <<
star.MeasureTotalAngularHomentum().GetNorm()
<< " Js" <<endl;
cout << "Ratio lE_rot/E-potl: " <<
fabs(star.HeasureKineticEnergy(true)/
star.MeasureGravitationalEnergy()) << endl;

// give t/ps a kick towards the center of the core
star.RickTPInside(kick);

else
star.ReadAllData(); // changes some parameters read from the file

cout << "Enter number of time steps: ";
if(cin >> steps)

{
// calculate time development: step size. number of steps
// , ..., flags indicating if stepsize
// may be modified during the calculation
star.Stepper(stepsize, steps. numb_part_save.
modify_stepsize, modify-stepsize-avvel);
}

// saves final configuration if this is a new simulation

124

if(new-simulation)
star.SaveAllData(radius/1, steps);

cout << "Bye.\n";
getchar();

return 0;

A.2 testparticle.cpp and testparticle.h

These files contain the C++ classes TestParticle and TPChange. All functions and
data directly related to individual test particles (like their position and momentum

vectors) are included in these.

testparticle.h:

////////////////////////////////////////////////////////l/////////
// testparticle.h

// Header file for testparticle.cpp
//////////////////////////////////////////////////////////////////

tifndef TESTPARTICLE_H
Odefine TESTPARTICLE_H

Sinclude "vector_and_spherical.h"

// forward decleration
class TPChange;

class TestParticle
{
private:

Vector pos; // position vector

Vector mom; // momentum vector

int id; // identification no. of test particle

double internal_energy; // internal energy of this test particle
public:

TestParticle() // constructor

{
Vector nullvec(0.0.0);

pos 8 nullvec;
mom 8 nullvec;

id 8 0;
internal_energy80;

}
double GetInternalEnergy()
{
return internal-energy;
}

125

void SetInternalEnergy(double newinterg)

{
internal_energy 8 newinterg;
return;
}
int GetID() // access function for ID
{
return id;
}
void SetID(int newid) // set new ID
{
id 8 newid;
return;
}
Vector GetPos() // access function for position
{
return pos;
}
Vector GetHom() // access function for momentum
{
return mom;
}
void SetPos(Vector newpos) // set new position
{
pos 8 newpos;
return;
}
void SetMom(Vector newmom) // set new momentum
{
mom 8 newmom;
return;
}

// binary operator for addition of a TPChange object
friend TestParticle operator 8 (TestParticle tp, TPChange tpc);

};
////////////////////////////////////////////////l/l/////////////////

// this class is simply used to perform a position and a momentum
// change simultaneously
////////////////////////////////////////////////////////////////////
class TPChange

{
private:
Vector pos_chg;
Vector mom-chg;
public:
void SetPosChg(Vector newpc)
{
pos_chg 8 newpc;
return;
}
Vector GetPosChg()
{
return pos_chg;
}
void SetHomChg(Vector newmc)
{
mom-chg 8 newmc;
return;
}
Vector GetHomChg()
{
return mom_chg;
}

// binary operators for ...

friend TestParticle operator + (TestParticle tp, TPChange tpc);
friend TPChange operator 8 (TPChange tpc, double x);

friend TPChange operator 8 (double x, TPChange tpc);

126

friend TestParticle operator + (TPChange tpc, TestParticle tp);
friend TPChange operator + (TPChange tpcl, TPChange tpc2);

tendif

testparticle.cpp:

///////////////////////////////////////
// test particle class

///l//////////////////////////////////
// operators to multiply both the momentum and the position of a TPChange
// object by a constant
TPChange operator 8 (TPChange tpc, double x)
{
TPChange temp;

temp.SetPosChg(tpc.GetPosChg() 8 x);
temp.SetMomChg(tpc.GetMomChg() 8 x);

return temp;

}

TPChange operator 8 (double x, TPChange tpc)

{ TPChange temp;
temp.SetPosChg(tpc.GetPosCth) 8 x):
temp.SetMomChg(tpc.GetMomChg() 8 x):
return temp;

}

// binary operators for addition of a TPChange object to a TestParticle obj
TestParticle operator + (TestParticle tp, TPChange tpc)
{
TestParticle temp;
temp.SetPos(tp.GetPos()+tpc.GetPosChg());
temp.SetHom(tp.GetHom()+tpc.GetMomChg());
temp.SetID(tp.GetID()):
temp.SetInternalEnergy(tp.GetInternalEnergy()):

return temp;

}

TestParticle operator + (TPChange tpc, TestParticle tp)
{
TestParticle temp;
temp.SetPos(tp.GetPos()+tpc.GetPosChg());
temp.SetHom(tp.GetHom()+tpc.GetHomChg());
temp.SetID(tp.GetID()):
temp.SetInternalEnergy(tp.GetInternalEnergy());

return temp;

}

// binary operator for adding two TPChange objects
TPChange operator 8 (TPChange tpc1, TPChange tpc2)

TPChange temp;
temp.SetPosChg(tpc1.GetPosChg()+tpc2.GetPosChg());
temp.SetHomChg(tpc1.GetHomChg()+tpc2.GetHomChg());

return temp;

127

A.3 vector_and_spherical . cpp and
vector_and_spherica1 . h

These files contain the C++ classes Vector and Spherical. These are needed to store
vectors in cartesian or spherical coordinates. Functions for operations like vector
addition, cross and dot products, or the conversion from cartesian to spherical coor-

dinates (i.e. Vector to Spherical objects) are also implemented here.

vector_and-spherica1.h:

71"
l

tifndef VECTOR-AND_SPHERICAL_H
tifndef IOSTREAM-H

tinclude <iostream.h>

8define IDSTREAH-H

Sendif

Oifndef HATH-H

Sinclude <math.h>

tdefine HATH_H

Sendif

// forward declaration
class Spherical;

III////////////////////////////////
// Vector with three coordinates
///////////////////////////////////
class Vector
{
private:
double coord[3];
public:
// constructor
Vector (double x, double y, double 2)

 

{

SetCoords(x,y,z):
}
Vector ()
{

SetCoords(0,0,0);
}

// returns the coordinates of the vector
double GetCoord(int i)

{
if(i>80 tt i<3)
return coordEi];
else
cout << "Error in Vector::GetCoord: invalid index." << endl;
return 0;
}

128

};

// sets the coordinates
void SetCoord(int i, double val)

{
if(i>80 tt i<3)
{
coord[i] 8 val;
return;
}
else
cout << "Error in Vector::SetCoord" << endl;
}

// sets all coordinates at once

void SetCoords(doub1e x, double y, double 2)

{
coord[0] 8 x, coord[1] 8 y, coord[2] 8 2;
return;

}

// returns the squared norm of the vector
double GetNorquuared()

{
double normsq80.0; // squared norm of vector
for(int i=0; i<3;i++)
normsq+8coord[i]8coord[i];
return normsq;
}

// returns the norm of the vector
double GetNorm()
{

return sqrt(GetNorquuared());
}

// converts to a Spherical object
Spherical GetSpherical();

friend Vector operator 8 (double, Vector);

friend double operator ‘ (Vector avec, Vector bvec);
friend Vector operator 8 (Vector avec, double lambda);
friend Vector operator 8 (Vector bvec, Vector avec);
friend Vector operator + (Vector avec, Vector bvec);
friend Vector operator - (Vector avec, Vector bvec);

////////////////////////////////////////////////

// 3 component Vector in spherical coordinates

////////////////////////////////////////////////

class Spherical

{

private:

public:

double r, phi, theta;

// constructors
Spherical()
{
r80. phi80, theta80;
}
Spherical(double nr, double nphi, double ntheta)
{
SetR(nr);
SetPhi(nphi);
SetTheta(ntheta);
}
// access functions
double GetR()

129

a...

{
return r;
}
double GetPhi()
{
return phi;
}
double GetTheta()
{
return theta;
}

// sets r. r must be >80, returns true if error occurs
bool SetR(double nr)

 

{
if(nr>=0)
{
r8nr;
return false;
} 0
else
{ I
cout << "Error in Spherical::SetR(): Invalid value for r: "
<< nr << endl;
r80;
return true; _
}
}

// sets phi, nphi is adjusted so that 0<8phi<=2PI
void SetPhi(double nphi)

{
phi8nphi>=0?nphi-(int)(nphi/(28PI))828PI:
nphi-(int)(nphi/(28PI)-1)828PI;
return;
}

// sets theta. theta must be >80 and <8 PI, returns true if error occurs
bool SetTheta(double ntheta)

{
if(ntheta>80 tt ntheta <8 PI)
{
theta 8 ntheta;
return false;
}
else
{
cout << "Error in Spherical::SetTheta(): invalid value
for theta.\n";
theta 8 0;
return true;
}
}

// returns vector in cartesian coordinates
Vector GetCartesian();

 

};

Cdefine VECTOR-AND_SPHERICAL_H
tendif

vector-and-spherica1 . cpp:
tinclude "vector_and_spherical.h"

// binary operator for scalar product
double operator “ (Vector avec, Vector bvec)

{
double temp 8 0;

130

for(int i=0; i<3; i++)
{

temp+=bvec.GetCoord(i)8avec.GetCoord(i);

}

return temp;

// s-multiplication (right)
Vector operator 8 (Vector avec, double lambda)

{
Vector temp;
for(int i=0; i<3; i++)
temp.SetCoord(i,1ambda8avec.GetCoord(i));
return temp;
} I

// s-multiplication (left)
Vector operator 8 (double lambda, Vector avec)

{
Vector temp; L .
for(int i=0; i<3; i++)
temp.SetCoord(i,lambda8avec.GetCoord(i)):
return temp;
}

// binary operator for vector product ("cross")
Vector operator 8 (Vector bvec, Vector avec)

{
Vector temp;
for(int i=0; i<3; i++)
temp.SetCoord(i,bvec.GetCoord((i+1)%3)8avec.GetCoord((i+2)%3) -
bvgc,GQtCoord((i+2)23)8avoc.GetCOOrd((i*1)z3))i
return temp;

}

// binary operator for vector addition
Vector operator + (Vector avec, Vector bvec)

 

{
Vector temp;
for(int i=0; i<3; i++)
temp.SetCoord(i,avec.GetCoord(i)+bvec.GetCoord(i));
return temp;

}

// binary operator for vector subtraction
Vector operator - (Vector avec, Vector bvec)

{
Vector temp;
for(int i=0; i<3; i++)
temp.SetCoord(i,avec.GetCoord(i)-bvec.GetCoord(i));
return temp;

}

losesseeeeeeeeesesassesseseessseeeeeeseeeeeseeeteseeeeesssse\
8 Converts a Spherical object (r,phi,theta) to a Vector

131

8 object (cartesian coordinates).

8 Needs: - math.h

8 Gets: - nothing

8 Returns: - Vector
\eeeeeeeeeeeeeeeeeeeeeeseeeeeeeeseeeeeeeseeeeeeseseeeeeeeeee/

Vector Spherical::GetCartesian()

{
Vector cart;
cart.SetCoord(0, r8sin(theta)8cos(phi)):
cart.SetCoord(1, r8sin(theta)8sin(phi));
cart.SetCoord(2, r8cos(theta));
return cart;

}

[88888888888888888888888888888888888888888888888888888888888\
8 Coverts a Vector onject to a Spherical object.

8 Needs: - math.h

8 Gets: - nothing

8 Returns: - Spherical
\eeessenseeeeeeseessseesesteesseeeeesesseeeeessseessseeeeeee/

Spherical Vector::GetSpherical()
{

Spherical spher;

double temp;

temp 8 atan2(this->GetCoord(1),this->GetCoord(0));
spher.SetR(this->GetNorm()):
spher.SetPhi(temp>=0?temp:temp+28PI);
spher.SetTheta(acos(this->GetCoord(2)/spher.GetR()));

return spher;

132

 

Appendix B

Source Code of the Output
Programs

The source code of the two output programs written in Microsoft Visual Basic©

will be reproduced in this part of the appendix.

B. 1 output . vbp

output.vbp is a Microsoft Windows© based program that reads and displays the
whole time development of the positions of up to 2000 test particles during a sim-
ulation run from a data file created by the simulation program. The particles are
shown in a pseudo three dimensional plot for each time step. One can zoom in and
out and rotate the particles around three axes. To better visualize the dynamics a

line indicating the particle’s velocity can be shown in its stead.

output.vbp:

Dim pointsperstep As Integer ’ number of points per time step
Dim pointshow As Integer ’ number of points actually shown
Dim step() As Double ’ array in which all points will be saved
Dim xax(2) As Double ’ screen coordinates for x axis

Dim yax(2) As Double ’ screen coordinates for y axis

Dim zax(2) As Double ’ screen coordinates for z axis

133

Public stepno As Integer ’ number of steps

Public flag1 As Boolean ’ flag indicating if files point have been loaded yet
Public nofs As Integer ’Number of time steps

Public lambda As Double ’ scale factor for points

Dim gamma As Double ’scale factor for velocity lines

Public longestaxis As Double ’ longest axis (on screen)

Public maxi As Double ’ largest coordinate of all points’ coord.s

Public mfx As Double ’ middle of picture1 (horizontal)

Public mfy As Double ’ middle of picture1 (vertical)

Const pi 8 3.14159265358979 ’guess what...

’ calculates and returns minimum of two doubles
Function minimum(a As Double, b As Double) As Double
If a < b Then

minimum 8 a

Else

minimum 8 b

End If

End Function

’calculates and returns maximun of three doubles
Function maximum(a As Double, b As Double, c As Double) As Double
If a > b Then

temp 8 a

Else

temp 8 b

End If

If temp > c Then

maximum 8 temp

Else

maximum 8 c

End If

End Function

’ show / don’t show axis

Private Sub CheckAxis-Click()

Call draw_main

End Sub

’ show / don’t show velocity lines
Private Sub Checleines_Click()
Call draw_main

End Sub

’ draw all pictures

Private Sub Command1_Click()
Call draw-main

End Sub

Private Sub Form_Load()

’ auto redraw form and pictures when form is scaled...
Form1.AutoRedraw 8 True

Picture1.AutoRedraw 8 True

’ use full form for picture1
Picture1.Height 8 Form1.Height - 620
Picture1.Hidth 8 Form1.Hidth - 3420

pointsperstep 8 find-pointsperstep ’finds points per step
pointshow 8 pointsperstep ’show all points by default

nofs 8 find-number,of-steps ’finds number of steps
gamma 8 38 ’set scale factor for velocity lines to 3.0

’ set flag indicating that data was not loaded yet
flag1 8 False

’ set maximum value for time step scrollbar

134

 

 

ScrollStep.Hax 8 nofs - 1

’ dimensionates array for saving points
ReDim step(nofs, pointsperstep - 1, 2)

End Sub

’resize picture1 if form is resized
Private Sub Form_Resize()

If Form1.Height - 620 > 0 Then
Picture1.Height 8 Form1.Height - 620
End If

If Form1.Hidth - 3420 > 0 Then
Picture1.Hidth 8 Form1.Hidth - 3420
End If

Call draw-main

End Sub

Private Sub HScrollgamma_Change()

gamma 8 168 8 HScrollgamma.Value / HScrollgamma.Max
Call draw_main

End Sub

Private Sub HScrollgamma_Scroll()
Call HScrollgamma_Change
End Sub

Private Sub ScrollPhi-Change()
Call draw-main
End Sub

Private Sub ScrollPhi,Scroll()
Call draw_main
End Sub

Private Sub ScrollPsi_Change()
Call draw_main
End Sub

Private Sub ScrollPsi_Scroll()
Call draw_main
End Sub

Private Sub ScrollShowPoints_Change()
Call draw_main
End Sub

Private Sub ScrollShowPoints_Scroll()
Call draw-main
End Sub

Private Sub ScrollStep-Change()
stepno 8 ScrollStep.Value

Call draw_main

End Sub

Private Sub ScrollStep_Scroll()
stepno 8 ScrollStep.Value

Call draw_main

End Sub

Private Sub ScrollTheta_Change()
Call draw_main
End Sub

Private Sub ScrollTheta_Scroll()

Call draw_main

End Sub

’ loads all points from "suno_all.dat" and saves ’em to array step(,,)

135

.1 —"_d nm- a"

Function load_points()

’On Error Resume Next

Dim fs, a ’filesystem object, filestream

Dim ret As String ’for temprary storage of read lines
Dim i As Integer ’point index

Dim s As Integer ’step index

Dim c As Integer ’coordinate index

’ opens file stream

Set fs 8 CreateObject("Scripting.FileSystemObject")
Set a 8 fs.0penTextFile("suno_all.dat", 1, False)
If Err Then

Exit Function

End If

’ skip 1st line (containing the number of time steps)
a.skip1ine

’ reads and saves the points
s 8 0

Do while s < nofs
For i 8 0 To pointsperstep - 1
For c 8 0 To 2
ret 8 a.readline
step(s, i, c) 8 CDbl(ret) ’stores read components to array step
’ check if current maxi is exceeded
If Abs(step(s, i, c)) > maxi Then
maxi 8 Abs(step(s, i, c))
End If
Next c

Next i

s 8 s + 1

a.skip1ine ’skips line containing the "end of time step" string
Loop

a.Close ’close file stream

End Function

’ calculates mfx and mfy

Function calculate-mfx_mfy()

mfx 8 Picture1.ScaleHidth / 2

mfy 8 Picture1.ScaleHeight / 2

End Function

’ calculates and draws axes, necessary before calling function draw-main()
Function draw-axis()

theta 8 pi 8 (-0.5 + ScrollTheta.Value / ScrollTheta.Hax) ’1st rotation around 2
phi 8 pi 8 (-0.5 + ScrollPhi.Value / ScrollPhi.Hax) ’1st rotation around x
psi 8 pi 8 ScrollPsi.Value / ScrollPsi.Hax ’2nd rotation around 2

’ rotation z -> x -> 2 by theta,phi,psi
xax(O) 8 Cos(psi) 8 Cos(theta) - Sin(psi) 8 Cos(phi) 8 Sin(theta)

xax(1) 8 Sin(psi) 8 Cos(theta) + Cos(psi) 8 Cos(phi) 8 Sin(theta)
xax(2) 8 Sin(phi) 8 Sin(theta)

yax(O) 8 -Cos(psi) 8 Sin(theta) - Sin(psi) 8 Cos(phi) 8 Cos(theta)
yax(1) 8 -Sin(psi) 8 Sin(theta) + Cos(psi) 8 Cos(phi) 8 Cos(theta)
yax(2) 8 Sin(phi) 8 Cos(theta)

zax(0) 8 Sin(psi) 8 Sin(phi)

zax(1) 8 -Cos(psi) 8 Sin(phi)

zax(2) Cos(phi)

’ set longestaxis to the length of the longest axis on screen display
longestaxis 8 Sqr(maximum(xax(0) ‘ 2 + xax(1) “ 2, yax(O) ‘ 2 +
yax(1) ‘ 2. zax(0) ‘ 2 + zax(1) ‘ 2))

’ scale factor for axes, makes sure axes are completely visible
axfac 8 minimum(Picture1.ScaleNidth, Picture1.ScaleHeight) / 2.1 / longestaxis

136

 

l_

’ draws and labels coordinate axes
If CheckAxis.Value 8 1 Then

Picture1.Line (mfx - xax(0) 8 axfac, mfy

axfac, mfy - xax(1) 8 axfac)

Picture1.PSet (mfx + xax(0) 8 axfac, mfy

Picture1.Print "x"

Picture1.Line (mfx - yax(O) 8 axfac, mfy

axfac, mfy - yax(1) 8 axfac)

Picture1.PSet (mfx + yax(0) 8 axfac, mfy

Picture1.Print "y"

Picture1.Line (mfx - zax(0) 8 axfac, mfy

axfac, mfy - zax(1) 8 axfac)

Picture1.PSet (mfx + zax(0) 8 axfac, mfy

Picture1.Print "2"

’displays units, scale and tics for axis

’positive x-axis

+ xax(1)
- xax(1)
+ yax(1)
- yax(1)
+ zax(1)

- zax(1)

axfac)-(mfx + xax(O) 8
axfaC). RGB(255, O, 0)
axfac)-(mfx + yax(O) 8
axfac). RGB(255, 0, 0)
axfac)-(mfx + zax(0) 8

axfac), RGB(255, 0, 0)

lambda 8 minimum(Picture1.ScaleVidth, Picture1.ScaleHeight) / 2 / maxi
/ longestaxis 8 ScrollZoom / ScrollZoom.Hax 8 200

plotx 8 mfx + lambda 8 maxi 8 xax(O)
ploty 8 mfy - lambda 8 maxi 8 xax(1)
plotxend 8 mfx + lambda 8 maxi 8 xax(O)
plotyend 8 mfy - lambda 8 maxi 8 xax(1)
Picture1.Line (plotx, ploty)-(plotxend,
Picture1.PSet (plotxend, plotyend)
Picture1.Print ”+" t maxi t "Oil"

’ negative x-axis

plotx 8 mfx - lambda 8 maxi 8 xax(O)
ploty 8 mfy + lambda 8 maxi 8 xax(1)
plotxend 8 mfx — lambda 8 maxi 8 zax(0)
plotyend 8 mfy + lambda 8 maxi 8 xax(1)
Picture1.Line (plotx, ploty)-(plotxend,
Picture1.PSet (plotxend, plotyend)
Picture1.Print “-" t maxi t "[m]"
’positive y-axis

plotx 8 mfx + lambda 8 maxi 8 yax(O)
ploty 8 mfy - lambda 8 maxi 8 yax(1)
plotxend 8 mfx + lambda 8 maxi 8 yax(0)
plotyend 8 mfy - lambda 8 maxi 8 yax(1)
Picture1.Line (plotx, ploty)-(plotxend,
Picture1.PSet (plotxend, plotyend)
Picture1.Print "+" t maxi t "[m]"
’negative y-axis

plotx 8 mfx - lambda 8 maxi 8 yax(0)
ploty 8 mfy + lambda 8 maxi 8 yax(1)
plotxend 8 mfx - lambda 8 maxi 8 yax(O)
plotyend 8 mfy + lambda 8 maxi 8 yax(1)
Picture1.Line (plotx, ploty)-(plotxend,
Picture1.PSet (plotxend, plotyend)
Picture1.Print "-" t maxi t " [m] "

’ positive z-axis

plotx 8 mfx + lambda 8 maxi 8 zax(0)
ploty 8 mfy - lambda 8 maxi 8 zax(1)
plotxend 8 mfx + lambda 8 maxi 8 zax(0)
plotyend 8 mfy - lambda 8 maxi 8 zax(1)
Picture1.Line (plotx, ploty)-(plotxend,
Picture1.PSet (plotxend, plotyend)
Picture1.Print "+" t maxi t "[m]"

’ negative z-axis

plotx 8 mfx - lambda 8 maxi 8 zax(0)
ploty 8 mfy + lambda 8 maxi 8 zax(1)
plotxend 8 mfx - lambda 8 maxi 8 zax(0)
plotyend 8 mfy + lambda 8 maxi 8 zax(1)
Picture1.Line (plotx, ploty)-(plotxend,
Picture1.PSet (plotxend, plotyend)
Picture1.Print "-" t maxi t "[mfl"

’displays tics on half axes

+ axfac /
- axfac /
plotyend)

+ axfac /
- axfac /
plotyend)

+ axfac /
- axfac /
plotyend)

+ axfac /
- axfac /
plotyend)

+ axfac /
- axfac /
plotyend)

+ axfac /
- axfac /
plotyend)

20 8
20 8

20
20

20 8

20

20
20

20
20

20
20

yax(O)
yax(1)

8 yax(O)
8 yax(1)

xax(O)
8 xax(1)

8 xax(O)
8 xax(1)

8 xax(O)
8 xax(1)

8 xax(O)
8 xax(1)

137

’positive x-axis

lambda 8 minimum(PictureI.ScaleHidth, Picture1.ScaleHeight) / 2

I maxi I longestaxis 8 ScrollZoom I 50

plotx 8 mix 8 0.5 8 lambda 8 maxi 8 xax(O)

ploty 8 mfy - 0.5 8 lambda 8 maxi 8 xax(1)

plotxend 8 mix 8 0.5 8 lambda 8 maxi 8 xax(O) 8 0.8 8 axfac I 20 8 yax(O)
plotyend 8 mfy - 0.5 8 lambda 8 maxi 8 xax(1) - 0.8 8 axfac / 20 8 yax(1)
Picture1.Line (plotx, ploty)-(plotxend, plotyend)

’ negative x-axis

plotx 8 m1: - 0.6 8 lambda 8 maxi 8 zax(0)

ploty 8 mfy 8 0.5 8 lambda 8 maxi 8 xax(1)

plotxend 8 mix - 0.5 8 lambda 8 maxi 8 xax(O) 8 0.8 8 axfac I 20 8 yax(O)
plotyend 8 m1y 8 0.5 8 lambda 8 maxi 8 xax(1) - 0.8 8 axfac I 20 8 yax(1)
Picture1.Line (plotx, ploty)-(plotxend, plotyend)

'positive y-axis

plotx 8 mix 8 0.5 8 lambda 8 maxi 8 yax(O)

ploty 8 mfy - 0.6 8 lambda 8 maxi 8 yax(1)

plotxend 8 mix 8 0.5 8 lambda 8 maxi 8 yax(O) 8 0.8 8 axfac / 2O 8 zax(0)
plotyend 8 mfy - 0.5 8 lambda 8 maxi 8 yax(1) - 0.8 8 axfac / 2O 8 xax(1)
Picture1.Line (plotx, ploty)-(plotxend, plotyend)

’negative y-axis

plotx 8 mix - 0.5 8 lambda 8 maxi 8 yax(O)

ploty 8 mfy 8 0.5 8 lambda 8 maxi 8 yax(1)

plotxend 8 mix - 0.6 8 lambda 8 maxi 8 yax(O) 8 0.8 8 axfac / 2O 8 xax(O)
plotyend 8 mfy 8 0.6 8 lambda 8 maxi 8 yax(1) - 0.8 8 axfac I 20 8 xax(1)
Picture1.Line (plotx. ploty)-(plotxend. plotyend)

’ positive z-axis

plotx 8 mix 8 0.5 8 lambda 8 maxi 8 zax(0)

ploty 8 mfy - 0.6 8 lambda 8 maxi 8 zax(1)

plotxend 8 m1: 8 0.5 8 lambda 8 maxi 8 zax(0) 8 0.8 8 axfac / 20 8 xax(O)
plotyend 8 mfy - 0.5 8 lambda 8 maxi 8 zax(1) - 0.8 8 axfac I 20 8 xax(1)
Picture1.Line (plotx, ploty)-(plotxend. plotyend)

’ negative z-axis

plotx 8 mix - 0.5 8 lambda 8 maxi 8 zax(0)

ploty 8 m1y 8 0.5 8 lambda 8 maxi 8 zax(1)

plotxend 8 mix - 0.5 8 lambda 8 maxi 8 zax(0) 8 0.8 8 axfac / 20 8 zax(0)
plotyend 8 m1y 8 0.5 8 lambda 8 maxi 8 zax(1) - 0.8 8 axfac I 20 8 xax(1)
Picture1.Line (plotx, ploty)-(plotxend, plotyend)

End If

End Function
Function drau_points()

’scale factor for points
lambda 8 minimum(PictureI.Scalewidth, Picture1.ScaleHeight) I 2 / maxi
I longestaxis 8 ScrollZoom I 50

’scale factor for distz
sc-distz 8 127 I maxi I longestaxis / 2

For i 8 0 To pointshow — 1

’ calculates actual screen locations for plotting points

plotx 8 mfx 8 lambda 8 (step(stepno, i, 0) 8 xax(O) 8 step(stepno. i,1)
8 yax(O) 8 step(stepno. i, 2) 8 zax(0))

ploty 8 mfy - lambda 8 (step(stepno, i, 0) 8 xax(1) 8 step(stepno, i,1)
8 yax(1) 8 step(stepno, i, 2) 8 zax(1))

’ calculates distance orthogonal to the screen from point to screen
distzl 8 127 8 sc_distz 8 (step(stepno, i. 0) 8 xax(2) 8 step(stepno,i, 1)
8 yax(2) 8 step(stepno, i, 2) 8 zax(2))

’ displays velocity lines it corresponding box is checked
If Checleines.Value 8 1 And stepno > 0 Then

138

 

plotxold 8 mfx 8 lambda 8 (step(stepno - 1, i, O) 8 zax(0) 8

step(stepno - 1, i. 1) 8 yax(O) 8 step(stepno - 1, i, 2) 8 zax(0))
plotyold 8 mfy - lambda 8 (step(stepno - 1, i, O) 8 xax(1) 8

step(stepno - 1, i, 1) 8 yax(1) 8 step(stepno - 1. i. 2) 8 zax(1))
Picture1.Line (plotx, ploty)-(plotx 8 gamma 8 (plotxold - plotx). ploty

8 gamma 8 (plotyold - ploty)), RGB(255 - distzZ, 255 - distzx. 255 - distzZ)
Else

’ displays dots otherwise, dots are bigger and darker

’ the closer they are to the screen

Picture1.DravHidth 8 disth / 80 8 1

Picture1.PSet (plotx, ploty). RGB(255 - distzZ, 255 - disth, 255 - disth)
'Picture1.Circle (plotx, ploty), distzZ / 5, RGB(255, 0, 0)

’Picture1.Line (plotx. ploty)-(plotx 8 20, ploty 8 20). , BF

End If

Next 1

Picture1.DravWidth 8 1

End Function

Private Sub ScrollZoom_Change()
Call draw_main
End Sub

Private Sub ScrollZoom_Scroll()

Call draw_main

End Sub

’ finds and returns the number of points per time step
Function find_pointsperstep() As Integer

Dim fs, a 'filesystem object. file stream
Dim ret As String
Dim i As Integer ’counts number of lines

’ open file stream

Set fs 8 CreateDbject("Scripting.FileSystemObject")
Set a 8 fs.0penTextFile("suno_all.dat", 1, False)
If Err Then

Exit Function

End If

’ skip 1st line (containing the total number of steps)
a.skip1ine

i = '1
ret 8 MI

’ counts number of lines until end of timestep is reached
Do While ret <> "end of timestep"

rat 8 a.readline

i 8 i 8 1

Loop

’ 3 lines 8 one point
find_pointsperstep 8 i / 3
a.Close

End Function

Function drau_main()

’ loads data if this is the first call of the function

If flag1 8 False Then

schrott 8 HsgBox("All data will be loaded now. This may take a while...“,
vbOKDnly)

’ save points for all time steps in array step(,,)
Call load_points

flag1 8 True

End If

139

’ calculate number of points to be shown
pointshow 8 pointsperstep 8 (ScrollShowPoints.Value / ScrollShowPoints.Max)

’ display how many points are shown
Label7.Caption 8 "Showing " & CStr(pointshow) & " pts."

’ cleans old picture, draws new one
Picture1.Cls

Call calculate-mfx_mfy

Call draw-axis

Call draw_points

’ shows number of time step
Form1.Label4 8 stepno

End Function

’ reads the total number of time steps from the let line of suno_a11.dat
Function find_number-of_steps() As Integer

Dim fs, a

Dim x As Integer

Set fs 8 CreateObject("Scripting.FileSystemObject")
Set a 8 fs.0penTextFile("suno_all.dat", 1. False)
If Err Then

Exit Function

End If

’ read 1st line of file, convert to integer
x 8 a.readline

find_number-of_steps 8 x

a.Close
End Function

B.2 Suno Density Output.vbp

Suno Density Output .vbp is the second Microsoft Windows© based program that

reads and displays the whole time development of the mass density in a slice through

the x-z-plane of the core from a data file created by the simulation program. The

densities are indicated by different colors. A logarithmic density scale is implemented

to cover the full density range occuring in our simulations. A “cutoff” density below

which the color corresponding to the lowest density is always taken can be adjusted.

One can also zoom in and out.

Suno Density Output.vbp:

140

Dim big-r_max As Double ’greatest r_max value

Dim r-max() As Double ’array containing r_max-values for all time steps
Dim time() As Double ’array containing time for all time steps

Dim dens() As Double ’array containing densities for all boxes and time steps
Dim max_dens As Double ’maximum density in all time steps and boxes

Dim r_sep() As Double ’actual r-positions of seperations in r-direction
Dim theta_sep() As Double ’theta positions of separators in theta-direction
Dim i_cos_theta_max As Integer ’number of seperations in theta-direction
Dim i_r_max As Integer ’number of seperations in r-direction

Dim step As Integer ’ current time step

Dim PI As Double

Dim centerx As Integer

Dim centery As Integer

Dim max_radius As Integer

Dim nofs As Integer ’number of time steps

Function load_data()

Dim fs, a ’filesystem object, filestream

Dim ret As String ’for temprary storage of read lines
Dim i_r As Integer

Dim i_cos_theta As Integer

Dim i As Integer

Dim s As Integer ’step index

’ opens file stream

Set fs 8 CreateObject("Scripting.FileSystemObject")
Set a 8 fs.OpenTextFile("suno_dens.dat". 1, False)
If Err Then

Exit Function

End If

’ reads number of time steps from file
rat 8 a.readline
nofs 8 CInt(ret)

’reads i-r_max and i_cos_theta_max from file
ret 8 a.readline

i-r-max 8 CInt(ret)

ret 8 a.readline

i-cos_theta_max 8 CInt(ret)

’ reads r-positions of seperations in r-direction from file
ReDim r-sep(i_r_max)

For i_r 8 0 To i-r-max - 1

rat 8 a.readline

r_sep(i_r) 8 CDbl(ret)

Next i_r

r_sep(i-r-max) 8 1

’dimensionates array for theta-positions of separators in theta-direction
ReDim theta_sep(i_cos-theta_max)

’reads density maps and r_max for all steps from file
ReDim r_max(nofs - 1)

ReDim time(nofs - 1)

ReDim dens(nofs - 1, i_r_max - 1, 1, i_cos_theta_max - 1)

For s 8 0 To nofs - 1

ret 8 a.readline

r_max(s) 8 CDbl(ret)

ret 8 a.readline

time(s) 8 CDbl(ret)

If r_max(s) > big_r_max Then
big_r_max 8 r-max(s)

End If

For i-r 8 0 To i_r_max - 1

141

 

For i 8 0 To 1

For i_cos_theta 8 0 To i_cos_theta_max - 1
ret 8 a.readline

dens(s, i-r, i, i_cos_theta) 8 CDbl(ret)
If dens(s, i-r, i, i_cos_theta) > max_dens Then
max,dens 8 dens(s, i_r, i, i_cos-theta)
End If

Next i_cos_theta

Next i

Next i-r

Next s

a.Close
End Function

Private Sub Command1_Click()

If Check1.Value 8 1 Then

VScrollZoom.Value 8 5 ’set standard zoom factor
End If

For ix 8 0 To nofs - 1

HScrollTimeStep.Value 8 i1

’ modify zoomfactor if necessary and check1 checked

If Check1.Value 8 1 Then

If r_max(step) I big_r_max 8 VScrollZoom.Value / VScrollZoom.Max 8 200 < 0.5 Then
VScrollZoom.Value 8 VScrollZoom.Value 8 2

Else

If r-max(step) / big_r_max 8 VScrollZoom.Value / VScrollZoom.Hax 8 200 > 1.1 Then
VScrollZoom.Value 8 VScrollZoom.Value / 2

End If

End If

End If

Call HScrollTimeStep-Change

SavePicture Picture1.Image, "SuNo" t get4digitnumber(iZ) t ".bmp"
Next i1

End Sub

Private Sub Command2_Click()

SavePicture Picture1.Image. "SuNo" I get4digitnumber(HScrollTimeStep.Value)
& ll . bmpll

End Sub

Private Sub Form_Load()

’ auto redraw form and pictures when form is scaled...
PI 8 3.141592

Form1.AutoRedraw 8 True

Picture1.AutoRedraw 8 True

Picture1.FillStyle 8 1

big_r-max 8 08

step 8 O

maxdens 8 08

’ use full form for picture1
Picture1.Neight 8 Form1.Height - 620
Picture1.Hidth 8 Form1.Hidth - 3420

HsgBox ("All data will be loaded now. This may take a while...“)
Call load_data

HScrollTimeStep.Hax 8 note - 1
HScrollTimeStep.Value 8 0

Call draw_density,key

Label2.Caption 8 CStr(max_dens) I " kg/m‘a"

142

Label3.Caption 8 CStr(Exp(HScrollDens.Value I HScrollDens.Max 8
(Log(max_dens) - Ill) 1 " kg/m‘3"

Label4.Caption 8 "Radius in plot: " t CStr(big-r_max) & “ m"
End Sub

Private Sub Form-Resize()
’resize picture1 if form is resized
Picture1.Cls

If Form1.Height - 620 > 0 Then
Picture1.Height 8 Form1.Height - 620
End If

If Form1.Hidth - 3420 > 0 Then
Picture1.Hidth 8 Form1.Hidth - 3420
End If

Call calculate_drawseps

Call draw-densities

If Form1.Neight - 685 > 0 Then
Labe13.Top 8 Form1.Height - 685
End If

Call draw-density_key

End Sub

Function calculate_drawseps()
Dim i_cos-theta As Integer
Dim costheta As Double

’ calculates center of picture1
centerx 8 Picture1.ScaleHidth / 2
centery 8 Picture1.ScaleHeight / 2

’ set maximum radius for const.-r-circles
max_radius 8 minimum(CInt(centerx), CInt(centery))

’ calculates theta-separators

For i_cos_theta 8 0 To i_cos-theta_max

costheta 8 (28 8 CDbl(i_cos_theta) I CDbl(i_cos_theta_max) - 18)
theta_sep(i-cos-theta) 8 arccos(costheta)

Next i_cos_theta

End Function

Function minimum(a As Integer, b As Integer) As Integer
If a > h Then

minimum 8 b

Else

minimum 8 a

End If

End Function

Function draw_densities()

Picture1.Cls

Dim i_r. i, i-eos-theta As Integer

Dim nextrdraw, rdraw As Long

Dim color As Long

Dim costheta. sintheta. nextcostheta. nextsintheta As Double
Dim rfactor As Double

Dim temp As Double

Picture1.DrawHidth 8 VScrollRes.Value 8 1

rfactor 8 r_max(step) 8 max_radius I big-r_max 8 VScrollZoom.Value

143

n
O .

/ VScrollZoom.Max 8 200

For i 8 0 To 1

For i_r 8 0 To i_r-max - 1

nextrdraw 8 r_sep(i_r 8 1) 8 rfactor

For i-cos-theta 8 0 To i_cos_theta_max - 1

rdraw 8 r_sep(i-r) 8 rfactor

’color 8 CInt(dens(step, i_r, i, i_cos-theta) I (max_dens 8
(1 - HScrollDens.Value / (HScrollDens.Max 8 1))) 8 255)
’If color > 255 Then

’color 8 255

’End If

If i 8 0 Then

theta_b 8 mod2pi(-theta_sep(i_cos_theta 8 1) 8 5 8 PI I 2)
theta-s 8 mod2pi(-theta-sep(i_cos_theta) 8 5 8 PI I 2)
Else

theta_s 8 theta_sep(i_cos_theta 8 1) 8 PI I 2

theta-b 8 theta-sep(i_cos-theta) 8 PI / 2

End If

’color 8 CLng(dens(step, i_r, i. i_cos_theta) / max_dens I
(1 - HScrollDens.Value I (HScrollDens.Hax 8 1)) 8 1020)

If dens(step, i-r, i, i-cos_theta) <> 0 Then

temp 8 HScrollDens.Value I HScrollDens.Max 8 (Log(max_dens) - 1)

color 8 CLng((Log(dens(step, i-r, i, i_cos_theta)) — temp) I
(Log(max-dens) - temp) 8 1275)

Else

color 8 0

End If

 

If color > 1275 Then
color 8 1275

Else

If color < 0 Then
color 8 0

End If

End If

Do While rdraw < nextrdraw
If color <8 255 Then
Picture1.Circle (centerx, centery), rdraw, RGB(255, 255, 255 - color),
theta_s. theta_b
Else
If color > 255 And color <8 511 Then
Picture1.Circle (centerx, centery), rdraw, RGB(511 - color, 255, 0),
theta_s, theta_b
Else
If color > 511 And color <8 767 Then
Picture1.Circle (centerx, centery), rdraw, RGB(0, 767 - color. color -
511).
theta_s, theta_b
Else
If color > 767 And color <8 1023 Then
Picture1.Circle (centerx. centery), rdraw, RGB(color - 768, 0, 1023 -
color).
theta_s, theta_b
Else
If color > 1023 Then
Picture1.Circle (centerx, centery), rdraw, RGB((1 - ((color - 1024) /
255) ‘ 4) 8 255, 0. 0),
theta_s, theta_b
End If
End If
End If
End If
End If

rdraw 8 rdraw 8 VScrollRes.Value

144

Loop

Next i_cos,theta
Next i_r

Next i

Call draw-density_key
End Function

Private Sub NScrollDens_Change()

’Labe12.Caption 8 CStr(max_dens I (10 ‘ (HScrollDens.Value / 10))) & "[kg/m‘SJ"
Labe13.Caption 8 CStr(Exp(HScrollDens.Value / HScrollDens.Hax 8

(Log(max_dens) - 1)))

t I. kS/m‘afl

Call draw_densities

End Sub

Private Sub HScrollDens_Scroll()

’Labe12.Caption 8 CStr(max_dens I (10 ‘ (HScrollDens.Value / 10))) & "[kg/m‘3J"
Labe13.Caption 8 CStr(Exp(HScrollDens.Value I HScrollDens.Hax 8

(Log(max_dens) - 1)))

t n kg/m-au
Call draw_densities
End Sub

Private Sub HScrollTimeStep_Change()

Picture1.DrawNidth 8 1

step 8 HScrollTimeStep.Value

Labe11.Caption 8 “Time step: " & CStr(step)
Call draw-densities

End Sub

Function draw,density_key()

’draws density key

’For i2 8 0 To 255

’Picture2.PSet (Picture2.Hidth I 2, Picture2.Height / 255 8 i2). RGB(255 - i2. i1. 0)
’Next ix

Picture1.DrawHidth 8 20

’ show radius in plot

Picture1.CurrentX 8 Picture1.ScaleNidth - 180
Picture1.CurrentY 8 1

Picture1.Print Label4.Caption

Picture1.CurrentX 8 25
Picture1.CurrentY 8 1
Picture1.Print Label2.Caption

’ show time in plot
Picture1.CurrentX 8 Picture1.ScaleHidth - 120
Picture1.CurrentY 8 Picture1.ScaleHeight - 15

Picture1.Print CStr(time(step)) & "s"

’ labels for scale

Picture1.CurrentX 8 25

Picture1.CurrentY 8 Picture1.ScaleHeight - 15
Picture1.Print Label3.Caption

Picture1.CurrentX 8 25
Picture1.CurrentY 8 1
Picture1.Print Label2.Caption

temp 8 HScrollDens.Value / HScrollDens.Hax 8 (Log(max_dens) - 1)

For denslabel 8 1 To 17
If 10 “ denslabel < max-dens Then

145

labelpos 8 Picture1.ScaleNeight 8 (1 - (denslabel 8 Log(10) - temp)

I (Log(max_dens) - temp))
Picture1.CurrentX 8 20
Picture1.CurrentY 8 labelpos
Picture1.Print "1e" & CStr(denslabel)
End If
Next denslabel

For ix 8 0 To 255

’Picture2.PSet (Picture2.Nidth / 2, Picture2.Neight - Picture2
I 1278 8 1%). RGB(255. 255, 255 - iZ)

Picture1.PSet (10, Picture1.ScaleHeight - Picture1.ScaleHeight
8 i1). RGB(255, 255, 255 - 1%)

Next 1%

For i1 8 256 To 511

'Picture2.PSet (Picture2.Nidth I 2, Picture2.Neight - Picture2
/ 1278 8 1%). RGB(511 - i1, 255, 0)

Picture1.PSet (10. Picture1.ScaleHeight - Picture1.ScaleHeight
/ 1278 8 1%). RGB(511 - i1. 255, 0)

Next i1

For i1 8 512 To 767

’Picture2.PSet (Picture2.Hidth I 2. Picture2.Neight - Picture2
/ 1278 8 i1). RGB(O. 767 - i2. i1 - 511)

Picture1.PSet (10, Picture1.ScaleHeight - Picture1.ScaleHeight
/ 1278 8 i1). RGB(0, 767 - i1. i1 - 511)

Next ii

For i% 8 768 To 1023

’Picture2.PSet (Picture2.Hidth I 2, Picture2.Height - Picture2
/ 1278 8 1%). RGB(iX - 768, 0. 1023 - 1%)

Picture1.PSet (10, Picture1.ScaleHeight - Picture1.ScaleHeight
I 1278 8 ix). RGB(iZ - 768. 0. 1023 - i2)

Next ii

For ix 8 1023 To 1278

’Picture2.PSet (Picture2.Hidth I 2, Picture2.Height - Picture2
I 1278 8 iI), RGB((1 - ((iZ - 1023) I 255) ‘ 4) 8 255, 0, O)

Picture1.PSet (10. Picture1.ScaleHeight - Picture1.ScaleHeight
/ 1278 8 ix). RGB((1 - ((11 ' 1023) / 255) ‘ 4) 8 255. 0, 0)

Next i1

End Function

Function arcsin(x As Double) As Double
If x <> 1 And x <> -1 Then

arcsin 8 Atn(x / Sqr(-x 8 x 8 1))
ElseIf x 8 1 Then

arcsin 8 PI I 2
ElseIf x 8 -1 Then
arcsin 8 -PI / 2
End If

End Function

Function arccos(x As Double) As Double

If x <> 1 And x <> -1 Then

arccos 8 Atn(-x I Sqr(-x 8 x 8 1)) 8 2 8 Atn(1)
ElseIf x 8 1 Then

arccos 8 0

ElseIf x 8 -1 Then
arccos 8 PI

End If

End Function

Private Sub VScrollRes_Change()
Picture1.Cls

146

.Height

/ 1278

.Height

.Height

.Height

.Neight

Call draw_densities
End Sub

Private Sub VScrollRes_Scroll()
Picture1.Cls

Call draw_densities

End Sub

Private Sub VScrollZoom_Change()
Label4.Caption 8 "Radius in plot: " t CStr(big_r_max I
(VScrollZoom.Value
I VScrollZoom.Max 8 200)) t " m"
Picture1.Cls
Call draw_densities
End Sub

Private Sub VScrollZoom_Scroll()
Label4.Caption 8 "Radius in plot: " t CStr(big_r,max /
(VScrollZoom.Value
I VScrollZoom.Max 8 200)) I " m"
Picture1.Cls
Call draw_densities
End Sub
'converts an arbitrary angle to a number between 0 and 2pi
Function mod2pi(x As Double) As Double
Dim n As Integer
n 8 CInt(x I (2 8 PI) - 0.5)
mod2pi 8 x - n 8 2 8 PI
End Function
’ converts as integer > 0 into a 4 digit string by adding zeros in front
Function get4digitnumber(i As Integer) As String
If i < 10 Then
get4digitnumber 8 ”000" t CStr(i)
Else
If i < 100 Then
get4digitnumber 8 "00" & CStr(i)
Else
If i < 1000 Then
get4digitnumber 8 "0" & CStr(i)
Else
If i < 10000 Then
get4digitnumber 8 CStr(i)
Else
Form1 . Pr int "Error”
Exit Function
End If
End If
End If
End If
End Function

147

 

Appendix C

Abbreviations and Symbols

M star
Nceus

(nr, ”4): ”cos 9)

mass number of a nucleus

speed of light z 2.99792458 x 108%
gravitational energy of core

(total) internal energy of core
(total) kinetic energy of core

total energy of core

equation of state

equations of state

function defining the locations of the grid boundaries

for the r coordinate

force on test particle 3' due to the EOS
force on test particle j due to gravity
gravitation constant 2 6.67259 x 10811;“;
Planck’s constant = % z %2- x 10‘34Js
Boltzmann’s constant z 1.3807 x 1043;“:-
Lattimer & Swesty EOS

(total) angular momentum of core
solar mass z 1.989 x 103°kg
baryon mass 8* 1.6749286 x 10‘27kg

mass of the iron core of the star

mass of a test particle

mass of the entire star

number of grid cells

grid coordinates (integers defining a grid cell)

 

148

 

pmin

pmax

T

tbounce

uint

Vol(n,., n¢, Tlcosg), Vol(n,)

number of grid boundaries for the ab coordinate
number of grid boundaries for the 0 coordinate
number of grid boundaries for the r coordinate
number of test particles

pressure

momentum vector of test particle j

distance from the center of the star (or coordinate
system)

radial location of the edge of the grid

position vector of test particle j

(mass) density

saturation density of (isospin symmetric) nuclear
matter

minimum density enforced

maximum density in a simulation run
temperature

time of core bounce in a simulation run
internal energy per baryon

volume of the grid cell (n,, 72.4,, ncosg)

electron fraction

lepton fraction

charge number of a nucleus

zero age main sequence

149

BIBLIOGRAPHY

150

Bibliography

[1] W.D. Arnett. Neutrino trapping during gravitational collapse of stars. Astrophys.
J., 218:815—833, December 1977.

[2] E. Baron and J. Cooperstein. The effect of iron core structure on supernovae.
Astrophys. J., 353:597--611, April 1990.

[3] E. Baron, J. Cooperstein, and S. Kahana. Supernovae and the nuclear equation
of state at high densities. Nucl. Phys. A, 440:744—754, 1985.

[4] H.A. Bethe. Supernova mechanisms. Rev. Mod. Phys, 62(4):801—866, October
1990.

I.

[5] S. Bonazzola and J .A. Marck. Efficiency of gravitational radiation from axisym-
metric and 3d stellar collapse. Astron. Astrophys., 267:623-633, 1993.

[6] IN. Bronstein, K.A. Semendjajew, G. Musiol, and H. Miihlig. Taschenbuch der
Mathematik. Verlag Harri Deutsch, Frankfurt am Main, third edition, 1997.

[7] SW. Bruenn, K.R. De Nisco, and A. Mezzacappa. General relativistic effects
in the core collapse supernova mechanism. Astrophys. J., 560:326—338, October
2001.

[8] A. Burrows. Supernova explosions in the universe. Nature, 403:727—733, Febru-
ary 2000.

[9] A. Burrows and J. Goshy. A theory of supernova explosions. Astrophys. J.,
416:L75—L78, October 1993.

 

[10] A. Burrows, J. Hayes, and B. Fryxell. On the nature of core-collapse supernova
explosions. Astrophys. J., 450:830—850, September 1995.

[11] B.W. Carrol and DA. Ostlie. Modem Astrophysics. Addison-Wesley, 1996.

[12] C.L.Fryer. Mass limits for black hole formation. Astrophys. J., 522:413—418,
September 1999.

[13] J. Cooperstein. The equation of state in supernovae. Nucl. Phys. A, 438:722—739,
1985.

[14] J. Cooperstein and J. Wambach. Electron capture in stellar collapse. Nucl. Phys.
A, 420:591—620, 1984.

151

[15] H. Dimmelmeier, J .A. Font, and E. Mueller. Gravitational waves from relativistic
rotational core collapse. Astrophys. J., 560:L163—L166, October 2001.

[16] CL. Fryer. 3d sph core collapse simulations. http://qso.lanl.gov/ clf/main.html,
January 2002.

[17] CL. Fryer. private communication, January 2002.

[18] CL. Fryer and A. Heger. Core-collapse simulations of rotating stars. Astrophys.
J., 541:1033-1050, October 2000.

[19] RA. Gingold and J.J. Monaghan. Smoothed particle hydrodynamics: theory
and application to non-spherical stars. Mon. Not. R. astr. Soc., 181:375-389,
1977.

[20] A. Heger, N. Langer, and SE. Woosley. Presupernova evolution of rotating
massive stars. i. numerical method and evolution of the internal stellar structure.
Astrophys. J., 528:368-396, January 2000.

[21] M. Herant and W. Benz. Postexplosion hydrodynamics of sn 1987a. Astrophys.
J., 387:294-308, March 1992.

[22] M. Herant, W. Benz, and S. Colgate. Postcollapse hydrodynamics of sn1987a:
T wo—dimensional simulations of the early evolution. Astrophys. J., 395:6428653,
August 1992.

[23] M. Herant, W. Benz, W.R. Hix, C.L. Fryer, and SA. Colgate. Inside the su-
pernova: A powerful convective engine. Astrophys. J., 435:339—361, November
1994.

[24] H.-Th. Janka. Conditions for shock revival by neutrino heating in core-collapse
supernovae. Astron. Astrophys, 368:527—560, 2001.

[25] H.-T h. Janka, Th. Zwerger, and R. Monchmeyer. Does artificial viscosity destroy
prompt type-ii supernova explosions ? Astron. Astrophys, 268:360—368, 1993.

[26] J. Lattimer, A. Burrows, and A. Yahil. Type ii supernova energetics. Astrophys.
J., 288:644—652, January 1985.

[27] J.M. Lattimer and FD. Swesty. Equation of state version 2.7 (ls eos v2.7).
http: //www.ess.sunysb.edu/dswesty/lseos.html.

[28] J .M. Lattimer and FD. Swesty. A generalized equation of state for hot, dense
matter. Nucl. Phys. A, 535:331—367, 1991.

[29] J .M. LeBlanc and J .R. Wilson. A numerical example of the collapse of a rotating
magnetized star. Astrophys. J., 161:541—551, August 1970.

[30] B.-A. Li, C.M. Ko, and W. Bauer. Isospin physics in heavy-ion collisions at
intermediate energies. Int. J. Mod. Phys, 7(2):147—229, April 1998.

[31] M. Liebendorfer, A. Mezzacappa, F.-K. Thielemann, O.E.B. Messer, W.R. Hix,
and SW. Bruenn. Probing the gravitational well: No supernova explosion in

152

 

spherical symmetry with general relativistic boltzmann neutrino transport. Phys.
Rev. D, 63(103004), May 2001.

[32] O.E.B. Messer, A. Mezzacappa, S.W. Bruenn, and M.W. Guidry. A comparison
of boltzmann and multigroup flux-limited diffusion neutrino transport during

the postbounce shock reheating phase in core-collapse supernovae. Astrophys.
J., 507:353—360, November 1998.

[33] A. Mezzacappa, A.C. Calder, S.W. Bruenn, J.M. Blondin, M.W. Guidry, M.R.
Strayer, and AS. Umar. The interplay between proto—neutron star convection
and neutrino transport in core-collapse supernovae. Astrophys. J., 493:848-862,
February 1998.

[34] A. Mezzacappa, M. Liebendorfer, O.E.B. Messer, W.R. Hix, F.-K. Thielemann,
and SW. Bruenn. Simulation of the spherically symmetric stellar core collapse,
bounce, and postbounce evolution of a star of 13 solar masses with boltzmann
neutrino transport, and its implications for the supernova mechanism. Phys.

Rev. Lett., 86(10):1935, March 2001.

[35] R. Monchmeyer and E. Miiller. Timing Neutron Stars, volume 262 of NATO
ASI Ser. C. ASI, New York, 1989.

[36] W.H. Press, S.A. Teukolsky, W.T. Vetterling, and BB F lannery. NUMERICAL
RECIPES in C. Press Syndicate of the University of Cambridge, Cambridge,
UK, second edition, 1988.

[37] M. Rampp and H.-T. Janka. Spherically symmetric simulation with boltzmann
neutrino transport of core collapse and postbounce evolution of a 15mg star.
Astrophys. J., 539:L33—L36, August 2000.

[38] M. Rampp, E. Miiller, and M. Ruffert. Simulations of non-axisymmetric rota-
tional core collapse. Astron. Astrophys, 332:969—983, 1998.

[39] H. Shen, H. Toki, K. Oyamatsu, and K. Sumiyoshi. Relativistic equation of state
of nuclear matter for supernova and neutron star. Nucl. Phys. A, 637(3):435—450,
May 1998.

[40] H. Shen, H. Toki, K. Oyamatsu, and K. Sumiyoshi. Relativistic eos table.
http://physics.senkou.numazu-ct.ac.jp/sumi/eos/, March 2001.

[41] K. Sumiyoshi, H. Shen, K. Oyamatsu, M. Terasawa, H. Suzuki, S. Yamada, and
H. Toki. Unstable nuclei and eos table for supernova explosion and r-process in
relativistic many body approach. RIKEN Review, 26, January 2000.

[42] K. Sumiyoshi, M. Terasawa, H. Suzuki, S. Yamada, H. Toki, G.J. Mathews,
and T. Kajino. Relativistic simulations of supernovae and the r-process; a new

relativistic eos and nuclear reaction network. Nucl. Phys. A, 688:4780—4800, 2001.

[43] K. Sumiyoshi and H. Toki. Relativistic equation of state of nuclear matter for the
supernova explosion and the birth of neutron stars. Astrophys. J., 422:700—718,
February 1994.

153

[44] ED. Swesty, J .M. Lattimer, and ES. Myra. The role of the equation of state
in the ”prompt phase” of type ii supernovae. Astrophys. J., 425:195—204, April
1994.

[45] K. Takahashi, M.F. El Bid, and W. Hillebrandt. Beta transition rates in hot and
dense matter. Astron. Astrophys, 67:185—197, 1978.

[46] F.X. Timmes. The helmholtz eos. http://flash.uchicago.edu/ fxt/eos.shtml.

[47] F.X. Timmes and D. Arnett. The accuracy, consistency, and speed of five equa-
tions of state for stellar hydrodynamics. Astrophys. J. Suppl. 3., 125:277—294,
November 1999.

[48] F.X. Timmes and FD. Swesty. The accuracy, consistency, and speed of an
electron-positron equation of state based on table interpolation of the helmholtz
free energy. Astrophys. J. Suppl. 5., 126:501—516, February 2000.

[49] L. Wang, D.A. Howell, P. Hoflich, and J .C. Wheeler. Bipolar supernova explo-
sions. Astrophys. J., 550:1030—1035, April 2001.

[50] L. Wang and JG Wheeler. Supernovae are not round. Sky and Telescope,
January 2002.

[51] L. Wang, J.C. Wheeler, Z. Li, and A. Clocchiatti. Broadband polarimetry of
supernovae: Sn 1994d, sn 1994y, sn 1994ae, sn 1995d, and sn 1995 h. Astophys.
J., 467:435—445, August 1996.

[52] T.A. Weaver, G.B. Zimmermann, and SE. Woosley. Presupernova evolution of
massive stars. Astrophys. J., 225:1021-1029, November 1978.

[53] JR. Wilson. Numerical Astrophysics, page 422. Jones and Bartlett, Boston,
1985.

[54] C.-Y. Wong. Dynamics of nuclear fluid. vii. time-dependent hartree-fock approx-
imation from a classical point of view. Phys. Rev. C, 25(3):1460—1475, March
1982.

[55] SE. Woosley, N. Langer, and TA. Weaver. The evolution of massive stars
including mass loss: Presupernova models and explosion. Astrophys. J., 411:823—
839, July 1993.

[56] SE. Woosley and T .A. Weaver. The physics of supernova explosions. Ann. Rev.
Astron. Astrophys, 24:205—253, 1986.

[57] S. Yamada and K. Sato. Numerical study of rotating core collapse in supernova
explosions. Astrophys. J., 434:268-276, October 1994.

[58] S. Yamada and K. Sato. Gravitational radiation from rotational collapse of a
supernova core. Astrophys. J., 450:245—252, September 1995.

[59] T. Zwerger and E. Miiller. Dynamics and gravitational wave signature of ax-
isymmetric rotational core collapse. Astron. Astrophys, 320:209-227, 1997.

154

 

 

1818222”