NUCLEAR REACTIONS IN THE CRUST OF ACCRETING NEUTRON STARS
By
Kit Yu Lau

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
DOCTOR OF PHILOSOPHY
Physics
2012

ABSTRACT
NUCLEAR REACTIONS IN THE CRUST OF ACCRETING NEUTRON
STARS
By
Kit Yu Lau
There have been many discoveries from observations of accreting neutron stars in X-ray
binaries. Many of the observed phenomena such as superbursts or the cooling of quasipersistent transients during their quiescent state are affected by the thermal properties and
the composition of the crust. To model the nuclear energy release and crust composition,
we build up a first complete network with pycnonuclear fusion. We run a consistent nuclear
reaction network that follows the evolution of an accreted fluid element from the atmosphere
down to the inner crust. We take into account a majority of the most important nuclear
processes including electron capture, neutron capture, neutron emissions, β − decay, and
pycnonuclear fusion reactions. The result of the model shows that there is nuclear reaction
path splitting in the crust of accreting neutron stars due to the usage of finite electron capture
rates. The pycnonuclear fusion reactions can occur at a shallower depth than previously
thought. The composition deep inside the inner crust is mainly 40 Mg, independent of the
initial composition of the ashes of the outer crust. The inner crust is found to be very pure
no matter what the initial abundance of the ashes is in the outer crust. The neutron drip
which is the density that neutrons start dropping out from the nuclei, locates at a higher
density in our model. In general, the nuclear reaction path and the heat energy generation
in the inner crust are significantly different from the previous work.

ACKNOWLEDGMENTS

Without the help and support of many people, getting a PHD is impossible for me. First, I
would like to thank my advisor Hendrik Schatz for his guidelines and help during my PHD
life. Without his guidelines and support, I cannot finish my PhD.
I would like to thank my committee members Edward Brown, Carl Schmidt, Remco
Zegers and Vladimir Zelevinsky for their opinions and help, especially Edward Brown for his
advice on astrophysics. I would also like to thank Betty Tsang, Bill Lynch and Scott Pratt
for their help in my PhD life.
I would also like to thank my group mates Alfredo, Giuseppe, Anna, Zachary for their
help in my work.
Thank you my lab mates Ngoc, Ling-Ying, Brent, Jun, Hui, Liyuan for their help and
emotional support.
Thank you for my friends Mei, Pampa, Shraddha, Kai Chung, Satayaki, Saurabh, Rakhi,
Eric, Cecilia, Stephen for their company and emotional support.
Thank you my friends in Hong Kong for their international emotional support.
Finally, thank you my parents and my brother for their support, especially to my mother
for her continuous support.

iii

TABLE OF CONTENTS

List of Figures

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1 Introduction
2 Background
2.1 Birth of a neutron star . .
2.2 Structure of a neutron star
2.3 Accreting neutron stars . .
2.4 X-ray Bursts . . . . . . . .
2.5 Superbursts . . . . . . . .
2.6 Soft X-ray transients . . .

vi
1

.
.
.
.
.
.

4
5
5
8
9
12
13

.
.
.
.
.
.
.
.
.
.
.
.
.

19
20
21
22
22
23
24
26
26
32
37
37
38
40

.
.
.
.

46
46
69
72
76

5 Astrophysical applications
5.1 Ionic and electronic plasma frequency and temperature . . . . . . . . . . . .
5.2 Coulomb coupling parameter and crust melting . . . . . . . . . . . . . . . .
5.3 Shear modulus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

87
87
89
93

.
.
.
.
.
.

.
.
.
.
.
.

3 Model
3.1 Equation of State . . . . . . .
3.2 Thermal Structure . . . . . .
3.2.1 Neutrino cooling . . .
3.2.2 Heat Transport . . . .
3.3 Nuclear Reactions . . . . . . .
3.3.1 Electron capture rates
3.3.2 Beta decay rates . . .
3.3.3 Neutron capture rates
3.3.4 Pycnonuclear fusion .
3.4 Mass Models . . . . . . . . . .
3.4.1 FRDM model . . . . .
3.4.2 Hilf mass model . . . .
3.5 Numerical Methods . . . . . .

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

4 Results
4.1 56 Fe . . . . . . . . . . . . . . . .
4.2 56 Fe without pycnonuclear fusion
4.3 Other abundances . . . . . . . . .
4.4 Superburst ashes . . . . . . . . .

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.

iv

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.

5.4
5.5
5.6

Superfluidity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
Electrical and thermal conductivities . . . . . . . . . . . . . . . . . . . . . . 97
Shear viscosity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100

6 Conclusion & Discussion

105

Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

109

v

LIST OF FIGURES

2.1

The main stages of evolution of a neutron star. For interpretation of the
references to color in this and all other figures, the reader is referred to the
electronic version of this dissertation. . . . . . . . . . . . . . . . . . . . . . .

6

2.2

Schematic cross section of a neutron star . . . . . . . . . . . . . . . . . . . .

8

2.3

Roche Lobe of a binary system . . . . . . . . . . . . . . . . . . . . . . . . .

10

2.4

Abundances as a function of mass number of X-ray burst ashes. Abundances
are defined as the mass fraction over the mass number [32] . . . . . . . . . .

12

2.5

Light curve of Superburst from 4U1820-30 [35] . . . . . . . . . . . . . . . . .

14

2.6

Comparison of the light curve of the X-ray transient MXB1659-29 during
quiescence with model calculations [58] . . . . . . . . . . . . . . . . . . . . .

18

3.1

The neutron chemical potential versus the neutron density (m−3 ) in log scale

33

3.2

The electron capture Q-value for FRDM and HFB21 mass models for mass
A = 56 as a function of proton number. . . . . . . . . . . . . . . . . . . . .

39

TMass excess calculated with the Hilf mass model and from experiment for
carbon isotopes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

41

The initial composition in the chart of nuclides. Abundances are shown on a
logarithmic color scale from blue (lowest) to green to yellow to orange to red
(highest) with a lower cutoff at 1.0×10−15 . Thick black squares indicate stable
nuclei and thin black squares unstable but particle bound nuclie. Particle
unbound nuclei don’t have an outline but can still be part of the network.
Nuclides that are not part of the network are shaded gray. . . . . . . . . . .

47

Final composition. The lines denote “instantaneous” net reaction flows for the
timestep displayed. Red lines show flows in the direction of increasing neutron
number, blue lines flows in the direction of decreasing neutron number. See.
Fig. 4.1. for additional information. . . . . . . . . . . . . . . . . . . . . . . .

48

3.3

4.1

4.2

vi

4.3

4.4

4.5

4.6

4.7

4.8

4.9

The relationship between time (s), log pressure (dyn cm−2 ), log column depth
(g cm−2 ), log density (g cm−3 and electron Fermi energy(MeV) . . . . . . .

49

The two step electron capture process from 56 Fe to 56 Cr (see Fig.4.1 and 4.2
for details) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

51

Abundance as a function of column depth for 56 Fe, 56 Mn, 56 Cr, 56 V, 56 Ti,
56 Sc, 56 Ca . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

53

The fluxplot shows the splitting of the path at time 9.638 × 1010 s. (See Fig.
4.1. and 4.2. for details) . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

54

The fluxplot shows the small amount of 56 Ar does electron capture down to
18 C and the first pycnonuclear fusion. The red lines from low Z to high Z
show pycnonuclear fusion. (see Fig.4.1 and 4.2 for details) . . . . . . . . . .

57

The abundance change of isotopes 56 Ar, 62 Ar, 61 Cl, 42 Si and 40 Mg versus the
column depth. The decrease of 40 Mg at high column depth is because the
continusous depletion by the pycnonuclear fusion cycles converting 40 Mg into
neutrons. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

60

Time integrated reaction flows for the complete calculation with initial 56 Fe
composition. See Fig 4.1. and 4.2. for details. The reaction flows towards
lower neutron number are here green instead of blue. . . . . . . . . . . . . .

61

4.10 < Z > (This work: blue line, Haensel & Zdunik: black circles) and < A >
(This work: green line, Haensel and Zdunik: red squares) versus column depth. 63
4.11 The energy generation (MeV per nucleon) versus the column depth (g cm−2 ).
The black line shows the heat deposited into the crust by our model and the
red circles shows the table from P.Hasensl & Zdunik 1990 paper. . . . . . . .

65

4.12 The level scheme of electron capture [19]. . . . . . . . . . . . . . . . . . . . .

66

4.13 The impurity parameter versus column depth (g cm−2 ) for an initial 56 Fe
composition. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

67

4.14 The electron Fermi energy versus column depth (g cm−2 ) for an initial composition of 56 Fe. The line shows the electron Fermi energy of our model while
the red circles show the electron Fermi energy calculated by Haensel & Zdunik
[83]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

68

vii

4.15 The integrated heat energy deposited into the crust with (red) and without
(black) pycnonuclear fusion versus column depth (g cm−2 ). . . . . . . . . . .

70

4.16 The electron Fermi energy with and without pycnonuclear fusion versus column depth (g cm−2 ). The red line shows the calculation with pycnonuclear
fusion while the black line shows the calculation without pycnonuclear fusion

71

4.17 < Z > versus column depth (g cm−2 ) for an initial abundance 38 Ar, 42 Si,
41 K, 47 Sc, 46 Ti, 50 Cr, 61 Ni and 70 Ge . . . . . . . . . . . . . . . . . . . . . .

73

4.18 < N > versus column depth (g cm−2 ) for and initial abundance 38 Ar, 42 Si,
41 K, 47 Sc, 46 Ti, 50 Cr, 61 Ni and 70 Ge . . . . . . . . . . . . . . . . . . . . . .

74

4.19 < A > versus column depth (g cm−2 ) for and initial abundance 38 Ar, 42 Si,
41 K, 47 Sc, 46 Ti, 50 Cr, 61 Ni and 70 Ge . . . . . . . . . . . . . . . . . . . . . .

75

4.20 Reaction flows for an initial composition of 38 Ar (see Fig. 4.1 and 4.2 for
details). Here unbound nuclei are not distingusihed from bound nuclei . . . .

76

4.21 The initial composition of the calculation with superburst ashes . . . . . . .

77

4.22 Composition at the end of the calculation with initial superburst ashes(see
Fig. 4.1 and 4.2 for details) . . . . . . . . . . . . . . . . . . . . . . . . . . .

78

4.23 The first fusion of the superburst ashes . . . . . . . . . . . . . . . . . . . . .

80

4.24 Time integrated reaction flows for an initial composition of superburst ashes

81

4.25 The abundance of 72 Ti,62 Ar,70 Ca,46 Si,40 Mg versus column depth. . . . . . .

82

4.26 < Z >, < N >, < A > versus column depth. . . . . . . . . . . . . . . . . . .

83

4.27 The integrated energy deposition versus column depth for initial superburst
ashes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

84

4.28 Impurity parameter versus column depth.

. . . . . . . . . . . . . . . . . . .

86

The plasma temperature in log10 (K) versus column depth in log10 (g cm−2 )
from this work for an initial composition of 56 Fe (black) and superburst ashes
(red) as well as from Haensel and Zdunik (green) . . . . . . . . . . . . . . .

90

5.1

viii

The electron temperature log10(K) versus column depth in log10 (g cm−2 )
from this work for an initial composition of 56 Fe (black) and superburst ashes
(read) as well as from Haensel and Zdunik (green) . . . . . . . . . . . . . . .

91

5.3

The Γ versus pressure log10 (dyn cm−2 ). . . . . . . . . . . . . . . . . . . . .

94

5.4

The shear modulus in log10 (erg cm−3) versus column depth in log10 (g
cm−2 ) from this work for an initial composition 56 Fe (black) and superburst
ashes (red) as well as from Haensel and Zdunik (green) . . . . . . . . . . . .

96

The neutron abundance in log10 versus column depth in log10 (g cm−2 ) from
this work for an initial composition of 56 Fe (red) and superburst ashes (black)
as well as from Haensel & Zdunik (green) . . . . . . . . . . . . . . . . . . . .

98

5.2

5.5

5.6

The electron-impurity scattering frequency and electron-phonton scattering
frequency log10 (s−1 ) of superburst . . . . . . . . . . . . . . . . . . . . . . . 101

5.7

The electrical conductivity in log10(s−1 )versus column depth log10 (g cm−2 )
for an initial composition of 56 Fe (black), Superburst ashes (red) and Haensel
& Zdunik (green) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102

5.8

The shear viscosity in log10( g cm−1 c−1 ) versus column depth in log10 (g
cm−2 ) for initial abundance of 56 Fe,(black), superburst (red) and Hasenel
&Zdunik (green) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104

ix

Chapter 1
Introduction
Neutron stars consist of the densest matter in the Universe [1, 2]. They are the remnants
of the gravitational collapse of massive stars during supernovae events. Neutron stars are
mostly composed of neutrons. There are ions and electrons in the outer crust and electrons,
neutrons and nuclei in the inner curst. The core consists of superdense matter composed
of neutrons, protons, muons, and possibly hyperons, quarks, pions, and kaons [3]. Neutron
stars are supported by quantum degeneracy pressure due to the Pauli exclusion Principle,
which does not allow two neutrons to occupy the same quantum state simultaneously. X-ray
observations can help us understand these objects. Neutron stars are either observed in
isolation [4] or in binary systems where the neutron star orbits a companion star - a regular
star, a white dwarf, a black hole, or another neutron star. Some highly magnetized neutron
stars emit strong radio radiation and are radio pulsars. Radio pulsars in binary systems can
be used to determine accurately the mass of the neutron star through relativistic orbital
effects that affect the observed radio pulsations [5]. Neutron star masses between 1.35 [6] to
1.97 solar masses[7] have been observed. If the companion star is a white dwarf or a main
1

sequence star then mass transfer onto the neutron star can occur, leading to X-ray emissions
that can be highly variable. Examples for extreme variations in X-ray flux are recurrent
X-ray bursts, soft X-ray transients and Superbursts [8]. X-ray bursts are periodic and rapid
increases in luminosity peaked in the X-ray regime of the electromagnetic spectrum. Soft
X-ray transients emit X-rays for some time before they turn off for extended periods of
time, possibly due to instabilities in the accretion process. Superbursts are particularly long
energetic X-ray bursts. These surface phenomena probe the neutron star surface and can
also provide information on compactness [9]. Many of these observations are influenced by
nuclear processes. Nuclear reactions heat up the surface as well as the crust. When the
accretion stops and the crust starts cooling down, it continues to emit X-rays due to the
residual nuclear heat. The goal of our work is to understand the nuclear reactions that heat
the neutron star crust during accretion. These nuclear reactions also set the composition
as a function of depth. The crust composition in turn determines crust properties such as
thermal and electrical conductivities [10], elastic properties [10], viscosity [10], superfluidity
and superconductivity [10], neutrino emissions [10], g-mode oscillations [11] and evolution of
magnetic fields [12]. The reactions in the crust might also deform the crust leading to the
possibility of gravitational wave emission if the neutron star is spinning rapidly [13].

Detailed studies of the composition of accreting neutron star crusts were carried out
first by [14] and by [15]. The energy released by nuclear reactions in the accreting neutron
star crust was first calculated by Vartanyan and Ovakimova [16]. Haensel & Zdunik [17]
further developed the work by using different S factors for the pycnonuclear fusion rates and
different mass models. They calculated the compositional changes with depth by minimizing
the Gibbs Free energy of a Wigner-Seitz cell of accreted matter assuming only electron
2

captures, neutron emissions and absorption and pycnonculear fusion occur. They assume
the nuclear reactions happen at zero temperature, a single nuclide composition and infinitely
fast reaction rates [18]. Gupta et al [19, 20] went the next step by using a reaction network
which included realistic electron capture rates and neutron emissions. In our work, we use
a complete reaction network calculation with a self-consistent thermal structure and finite
electron and neutron capture rates. A new formalism for pycnonuclear reaction rates in a
multicomponent plasma is also used [21]. The results differ significantly from previous work
due to the more comprehensive and detailed modeling of the nuclear processes. Examples
for such differences include the splitting of the nuclear reaction path due to the usage of
finite electron capture rates. The pycnonuclear fusion reactions occur at a shallower depth.
The composition in the inner crust is substantially different and consists of mainly 40 Mg,
independent of the initial composition of the ashes of the outer crust. The inner crust is
found to be very pure no matter what the initial abundance of the ashes is in the outer crust.
In general, the nuclear reaction path and the heat energy generation in the inner crust are
significantly different from the previous work.
The thesis will be organized as follows: In chapter 2 some background information on
neutron stars and related phenomena is given. In chapter 3, the detail of our thermal and
nuclear model for accreting neutron star crusts is discussed. In chapter 4 follows a detailed
discussion of the nuclear reactions and heat generation. In chapter 5 a brief discussion of
the potential impact on various crust properties and observables follows.

3

Chapter 2

Background

Since the purpose of the thesis is to investigate the nuclear reactions in the crust of accreting
neutron stars, an introduction to formation and principal properties of neutron stars is
given in the following subsections. The idea of neutron stars was first proposed by Baade
and Zwicky in 1934 [22]. They suggested that a neutron star can be formed in a supernova
explosion. The first calculation of the neutron star structure was performed by Oppenhemier
& Volkoff in 1939 [23]. They assumed neutron star matter is composed of an ideal gas of
free neutrons at high density. However, this simple model failed to predict observed neutron
star masses. Later, Cameron calculated the structure of neutron stars using a Skyrme force
model in the 1950’s [24]. Neutron stars in the form of pulsars were first discovered in 1967
[25]. About a year later the discovery of a similar source in the Crab Nebula [26] further
supported the idea that neutron stars are the remnants of supernova explosions.
4

2.1

Birth of a neutron star

When a massive star(≥ 8M⊙ ) reaches the end of its life and undergoes gravitational collapse,
the core collapses into a neutron star or a black hole and a supernova explosion occurs.
The explosion mechanism is not well understood, but the most promising mechanism is
the so called delayed neutrino driven mechanism: the core collapse stops when the star’s
interior density reaches nuclear saturation density and a shock wave is produced at the
boundary between the collapsing inner core and the infalling outer core. The outgoing
shock wave loses energy due to neutrino emissions and nuclear dissociation of the infalling
material and therefore stalls at radii of about 100 to 200 km. The neutrinos from the core
resuscitate the shock, which then accelerates the stellar envelope outwards leading to mass
ejection. The proto-neutron star left behind shrinks rapidly because of pressure losses from
neutrino emissions. At the resulting high densities the high electron Fermi energy leads
to neutralization by fusing electrons and protons into neutrons. As a result the core heats
to temperatures around 6.0 × 1011 K. After 10 to 20 s, neutrino emission starts to cool
the neutron star [1].Figure 2.1 is a diagram showing the main stages of the evolution of a
neutron star. The radius R and central temperatures T for the neutron star are indicated
as it evolves in time t [1].

2.2

Structure of a neutron star

Neutron star masses in the range of 1.35 [6] to 1.97 [7] solar masses have been measured,
with radii to be less than 20 km. There are five major regions in a neutron star: atmosphere,
ocean, crust, outer core and inner core. The atmosphere and ocean are made of normal atoms.
5

Figure 2.1: The main stages of evolution of a neutron star [1]. For interpretation of the
references to color in this and all other figures, the reader is referred to the electronic version
of this dissertation.

6

They are thin shells with negligible mass compared to the total stellar mass. Although
the atmosphere is not hot and dense enough for nuclear reactions to happen, it plays an
important role in shaping the emergent photon spectrum. Beneath the atmosphere is the
ocean, which influences the transport and release of thermal energy from the star’s surface.
The crust which is about 1 to 2 km thick contains mainly of ionized nuclei and neutrons.
The dominant nuclei in the crust vary with density, and it is the main focus of this thesis to
investigate the composition of the crust of an accreting neutron star.
The crust is further divided into outer crust and inner crust. The outer crust consists of
ordinary nuclei and electrons and starts at a density around 1.0×106 g cm−3 . It reaches down
to neutron drip density, which is the density where neutrons drip out from nuclei, around
4.0 − 6.0 × 1011 g cm−3 . The inner crust begins at neutron drip and consists mostly of a
mix of neutrons, nuclei and electrons, and extends to the bottom of the crust at a density
around 1014 g cm−3 . Just above the core-crust interface, theory [1] predicts a transition
from three-dimensional (3D) nuclei, to 2D cylindrical nuclei, to 1D slabs of nuclei interlaid
with planar voids, to 2D cylindrical voids, and to 3D voids before a transition to uniform
nucleonic matter in the core. This series of transitions is known as nuclear pasta. At the
core-crust interface, nuclei are so close together that they are nearly touching each other.
Beneath the crust is the core which constitutes most of the mass of the star. The
composition of the core is less well known due to uncertainties in the equation of state for
matter exceeding nuclear matter saturation density (∼ 0.17f m3 ). The outer core is made
of nucleons, electrons, and muons. The inner core region possibly consists of more exotic
particles such as hyperons, pions or kaons. A Quark condensate may be found in the center
of the neutron star [1].
7

Figure 2.2: Schematic cross section of a neutron star

2.3

Accreting neutron stars

Accreting neutron stars accrete matter from a companion star and are bright X-ray sources,
so called X-ray binaries. They can be further divided into low mass X-ray binaries, intermediate X-ray binaries and high mass X-ray binaries depending on the mass of the companion
stars. Low mass X-ray binaries have donor stars’ masses ≤ 1M⊙ . High mass X-ray binaries
have companion stars with masses larger than about 10M⊙ while intermediate X-ray binaries have companion stars with masses in between [27]. In this thesis, our interest is in low
mass X-ray binaries (LMXB). These systems are usually found in the galactic bulge and in
globular clusters because the companion star tends to be old. Matter is typically transferred
from the companion star (primary star) to the compact star (secondary star) by Roche lobe
overflow. The Roche lobes is the equipotential suface that includes the inner Lagrange point
8

(see Fig. 2.3). Interior evolutionary changes of a star cause it to expand and fill up its
Roche lobe at which point matter can flow slowly over the gravitational potential saddle
point (Lagrange point) between the two stars. Because of angular momentum conservation,
the overflowing matter forms an accretion disk around the neutron star. Total accretion
rates are typically 1015 − 1018 g s−1 . When matter approaches the neutron star surface, it
is decelerated heating up the gas and leading to the emission of X-rays [2]. The Eddington
accretion rate represents an upper limit on the accretion rate as long as matter is distributed
evenly across the neutron star surface. The Eddington mass accretion rate is defined as
the rate at which the gravitational energy released by the in-falling matter generates the
Eddington luminosity, which is the maximum luminosity a star can emit while remaining
in hydrostatic equilibrium. Assuming Thomson scattering is the only opacity source, the
Eddington limit mEdd can be calculated as [28].
˙

mEdd =
˙

2mp c
10km
)gcm−2 s−1
= 8.8 × 104 (1.71/(1 + XH ))(
(1 + XH )RσT
R

(2.1)

where XH is the hydrogen mass fraction of the accreted material, mp is the proton mass, σT
is the Thomson cross section, R is the radius of the star and c is the speed of light.

2.4

X-ray Bursts

A new class of bursts, known as type I X-ray bursts, was discovered in 1975 [29]. Type I
X-ray bursts are thought to be powered by thermonuclear explosions on the surface of weakly
magnetized neutron stars [30]. These X-ray bursts last for seconds to minutes and have a
recurrence time of the order of hours or days. Their luminosity ranges from 1036 to 1038
9

Figure 2.3: Roche Lobe of a binary system

10

erg s−1 . Typical rise times of X-ray bursts are 2 s but can be up to 10 s, while the roughly
exponential decay time ranges from 10 s to several minutes. X-ray bursts are produced
when ionized hydrogen and helium are accreted from the envelope of the companion star,
accumulated on the neutron star, and ignited. During accretion the gravitational energy
release is about 200 MeV/u and powers the continuous X-ray emission of the X-ray binary.
Some of the X-ray bursts are powered by the rp-process (rapid proton capture process) which
releases less than 6 MeV/u, which is small compared to the gravitational energy. However,
X-ray bursts are observable as the nuclear energy is released in very short period of time [31].
Many bursts are powered by mixed hydrogen and helium burning. In this case the burst is
triggered by helium burning via the 3-α reaction. Cooling through expansion is prevented
because the matter is degenerate. When matter is degenerate, pressure is independent of
temperature, and therefore a temperature rise does not induce an expansion. The result is
a thermonuclear runaway as the rise in temperatures accelerates the thermonuclear reaction
rates. The acceleration in reaction rates further increases the temperature, which increases
the thermonuclear reactions and so on. Thus the temperature becomes quickly high enough
to allow 15 O and 18 Ne to undergo a (α, p) reactions leading to breakout from the hot CNO
cycle marking the onset of the (α, p) process. After the (α, p) process (a sequence of (α, p)
and (p, γ) reactions) burning of hydrogen via the rapid proton capture process (rp-process)
follows. At the drip line, the (γ, p) process slows (p, γ) down. In those cases, β + decays can
compete with proton capture. Peak temperatures as high as 1-2 GK have been predicted to
be reached in X-ray bursts [32]. The furthest point that the rp-process can reach are the
isotopes where a SnSbTe cycle is formed by (γ, α) reactions on α unbound proton-rich Te
isotopes [32]. Figure 2.4 shows the predicted composition of the X-ray burst ashes. As the

11

Abundance
10^(-2)
10^(-3)
10^(-4)
10^(-5)
10^(-6)
0

20

40

60

80

Mass Number

100 120

Figure 2.4: Abundances as a function of mass number of X-ray burst ashes. Abundances are
defined as the mass fraction over the mass number [32]
burst cools, the burning layer shrinks and forms a layer on the neutron star surface. The
fuel for the next burst is then accreted on top of this layer of ashes.

2.5

Superbursts

Superbursts are rare, extremely powerful X-ray bursts (see Fig. 2.5). They are 1000 times
more energetic than X-ray bursts, last for hours and have recurrence times of years. Superbursts are expected to occur at a density around 109 g cm−3 , which is 3 orders of magnitude
occurs deeper in density than regular X-ray bursts. The much higher ignition depth explains
the longer recurrence time and the higher total energy as more matter needs to be accumulated to trigger the burst. They have a rapid rise time and a long decay time of typically
a few hours, with similar spectra as X-ray bursts, indicating their thermonuclear origin.
˙
Superbursts are typically observed in a system with mass accretion rates of 0.1 − 0.3MEdd
12

but have also been found in more rapidly accreting sources [33]. Superbursts are believed
to be produced by the ignition of carbon in the ashes ocean on the neutron star surface.
However, the origin of significant amounts of carbon is still an open question. Superbursts
are an important probe for the properties of neutron stars since they are more sensitive to
the thermal properties of the crust and core cooling models as they occur in a deeper region
compared to ordinary X-ray bursts. Superbursts can help to constrain the thermal conductivity of the neutron star crust and the rate of neutrino cooling in the core [34]. Gupta
et al [19] found that there is roughly 4 times more heat deposited in the outer crust than
predicted before due to the electron captures to excited states. The higher heat release in the
outer crust leads to a higher temperature reducing the ignition depth of carbon burning and
shortening the recurrence time of superbursts. This helps to reduce the discrepancy between
the observation and theory, though it cannot completely account for the discrepancy.

2.6

Soft X-ray transients

Transient X-ray sources are quiescent for most of the time but are bright X-ray sources
during their active state, which lasts typically for 10-100 days [36]. They were initially
classified based on their spectral hardness because their nature was unclear [37]. Hard
X-ray transients are characterized by an equivalent temperature of > 15keV. Equivalent
temperature is the temperature of an air parcel from which all the water vapor has been
extracted by adiabatic process. They are believed to contain a young, pulsating neutron
star orbiting a Be star companion [38]. Soft X-ray transients (SXRTs) are characterized by
an equivalent temperatures of < 15keV and usually show type I X-ray bursts during their
active state. They are believed to belong to the class of low mass X-ray binaries (LMXRBs)
13

Figure 2.5: Light curve of Superburst from 4U1820-30 [35]

14

containing an old neutron star. There are also ultrasoft X-ray transients, which do not
contain any bursts or pulsations, and are believed to contain a black hole [39].
Soft X-ray transients are binary systems where the accretion onto the neutron star turns
on and off for periods ranging from months to decades, possibly due to the thermal accretion
disc instabilites [36]. During the accretion outburst, the neutron star emits X-rays similar to
persistent X-ray binaries. At the beginning of an outburst the flux increases over a few days
and reaches X-ray luminosities of ∼ 1037 − 1038 erg s−1 . It then decays by a slow, nearly
exponential decay with a timescale of weeks to months. When the accretion turns off, there
is a residual quiescent luminosity. For example, Aq1 X-1 emits ∼ 1032 − 1033 erg s−1 in the
quiescent state.
The spectral evolution can be described as follows: First there is the beginning of the
outburst, which shows relatively soft spectra when they are close to the outburst peak, (kT
∼ 5 keV). The outburst phase has intermediate luminosities (∼ 1035 − 1037 erg s−1 ) during
the rise and decay of the outburst. It declines with a high energy tail extending to at least
∼ 100 keV.
Most models of SXRT involve accretion disk instabilities. A compact star is continuously
accreting from a disk fed by a companion star filling its Roche lobe and transferring matter.
In disk instability models, the accretion disk can become thermally and viscously unstable
due to strong opacity variations caused by partial ionization of hydrogen [40]. [41] modified
the model by pointing out the importance of X-ray heating of the accretion disk which leads
to recurrence times.
The nature of the quiescent luminosity is not known with certainty. Four types of models
have been discussed: residual accretion onto the neutron star surface [36], residual accretion
15

to the magnetosphere radius which defines as the radiation pressure dominated region of
an accretion disc surrounding a neutron star [36], non-thermal processes powered by the
rotational energy loss of a rapidly spinning neutron star [42] and thermal emission from the
hot neutron star [36]. The latter is an important motivation for this thesis. When matter
accretes on the neutron star, it may burn steadily or have X-ray bursts in the atmosphere
of the neutron star. The ashes from the atmosphere are continuously pushed down to the
crust and undergo further nuclear reactions such as electron capture, neutron emission and
absorption and pycnonculear fusions [17]. These nuclear reactions release heat energy which
further heats the crust of the neutron star. While the core is isothermal, the crust is not
and can be at a higher temperature than the core due to the crust heating processes. When
accretion halts, we can see the luminosity from the crust.
There are observational evidences for the thermal emission model. The quiescent spectra
measured with Chandra ACIS-S of the transient neutron star Aq1 X-1 after the accretion
outburst is consistent with that predicted by Brown, Bildsten and Rutledge [43] from deep
crustal heating [44]. However, the observation of hard tails and flucutations in Aq1 X-1
and Cen X-4 [45] can hardly be reconciled with quiescent emission being powered by this
mechanism alone.
Observations of soft X-ray transients are important for constraining properties of the
neutron star, as probing the thermal state of the crust helps to determine the cooling mechanism of the neutron star. In that context quasi-persistent X-ray transients (long period
transients) provide the unique opportunity to measure the crust cooling behaviour of an
accreting neutron star as a function of time.
KS1731-260 and MXB1659-29 [46] are two such sources, where the quiescent emission is
16

believed to be powered by the thermal emissions from the cooling crust [47]. KS 1731-260
was first observed in August 1989. It was in its active state [48, 49] until early 2001 [50].
It then went into the quiescent phase within two months. It had a luminosity ∼ 1 × 1033
erg s−1 at that time. In 2001 September the source was found to have a luminosity of
approximately (2-5) ×1032 erg s−1 [51]. According to neutron star cooling models, this
decrease in luminosity can be interpreted as the rapid cooling of the neutron star crust.
MXB 1659-29 was first detected [52] in October 1979. In July 1979, it was in its quiescence
state [53]. It continued to exist in the state for 21 years until in April 1999 (Zand et al.
1999) [54] it switched to the active state. After 2.5 years in the active state, MXB 1659-29
went back to quiescence in September 2001 [55] with a residual luminosity of (3.2-4.3) ×1033
erg s−1 assuming a distance of 10kpc [56]. When MXB 1659-29 was in quiescence, it was
observed that the luminosity was decaying exponentially, decreasing by a factor of 7-9 over
1.5 years [57].
The cooling light curves of these quiescent transients can be used to probe the thermal
structure of the neutron star crust. An example for this is given by Brown et al who fitted
the light curve with a broken power law (see Fig. 2.6) [58]. The break in the power law is
explained by a change in heat capacity at the corresponding depth in the neutron star crust.
Such a change in heat capacity is predicted by the model because of the onset of a significant
fraction of free neutrons. The quantitative interpretation of the cooling light curves depends
on the understanding of the thermal structure of the neutron star crust, in particular the
distribution and strength of various nuclear heat sources. This is the topic of the thesis.

17

kTeff(eV)
140

Obs
Numerical
Model

120

100

80

60

10

100

1000

time(days)
Figure 2.6: Comparison of the light curve of the X-ray transient MXB1659-29 during quiescence with model calculations [58]

18

Chapter 3

Model

In order to understand the nuclear reactions and heat generation in the accreting neutron star
crust, we developed a dynamic model that includes electron capture rates, neutron capture
rates, and pycnonuclear fusion rates coupled to a thermal neutron star crust model. We use
a large nuclear reaction network of ∼ 3000 isotopes covering the mass range A=1-181 and
extending from stable nuclei to beyond the neutron-drip line in order to follow the nuclear
reaction sequences in detail. This approach has been used before by [19] up to neutron drip.
[20] then extendend their model to just beyond neutron drip line. We extend the [19] model
to include neutron emission and absorption processes, and for the first time, pycnonuclear
fusion rates. Our model follows the dominant nuclear processes in detail. It is very different
from most earlier approaches (except Gupta et al [19, 20]) where the Gibbs free energy of a
Wigner-Seitz cell at constant pressure is minimized, assuming zero temperature and where
the composition was limited to a single species at a given depth as in ref [17]. I will discuss
the equation of state, thermal structure, mass model and nuclear reactions used in the model.
19

3.1

Equation of State

The neutron star structure is taken from Brown et al [59, 34]. We assume a neutron star
radius of R = 10.8 km and a mass M of 1.6M⊙ . The resulting surface gravity GM/R2 (1 + z)
is 2.43 × 1014 cm s−2 where G is the gravitational constant. The gravitational red shift z
is [1 − 2(GM )Rc2 ]−1/2 − 1 = 0.33]. The crust mass is 0.01M⊙ and the crust thickness is
˙
about 0.3 km. A mass accretion rate in the rest frame on the surface of M = 3.0 × 1017 g
s−1 is chosen based on typically observed rates. Assuming spherical, isotropic accretion
this corresponds to a local accretion rate per unit area of 2.1 × 104 g s−1 cm−2 . For solar
composition the resulting luminosity generated by the gravitational energy release is 30%
of the Eddinton luminosity. The sensitivity of the calculations on neutron stars parameters
and accretion rates should be explored in future work.

Electrons, ions and neutrons contribute to the equation of state of the crust. Electrons
are relativistic and degenerate and the table of the Helmholtz free energy (a thermodynamic
potential that measures the ”useful” work obtainable from a closed thermodynamic system
at a constant temperature and volume) from [60] Timmes & Swesty is used. The pressure
from ions is calculated from the derivatives of the Helmholtz free energy with respect to the
ionic free energy. The equation of state for the ions lattice comes from fits by [61] Farouki&
Hamaguchi (1993), who modeled the molecular dynamics of a one component, stronglycoupled plasma. The neutron pressure is calculated from the zero temperature compressible
liquid drop model of Ref [62].
20

3.2

Thermal Structure

The thermal structure of the neutron star is calculated as in Brown (2000) [59]. The model
assumes hydrostatic equilibrium. The pressure is approximated to be nearly independent
of temperature. The gravitational potential Φ ≈ 1/2c2 ln(1 − 2GM/rc2 ) is nearly constant
in the crust. Under these approximations, the equations for temperature, luminosity, and
nuclear heating are solved as follows [59],

−L(1 + z)
dT
=
dr
4πr2 K

(3.1)

dL
dLnuc
=
− 4πr2 (1 + z)ǫν
dr
dr

(3.2)

dLnuc
˙
= M Q(r)m−1
amu
dr

(3.3)

where z = (1 − 2GM/rc2 )−1/2 − 1 is the gravitational redshift, Q(r) is the heat deposited
per accreted nucleon, and mamu = 1.66 × 10−24 g = 1 amu, K is the thermal conductivity,
ǫν is the neutrino energy loss rate. Since the potential of the crust Φ is nearly constant,
2
equations (1) and (2) are simplified by holding eΦ/c constant. The physics of interest is

contained in the neutrino cooling ǫν , the conductivity K and the nuclear heating Lnuc . We
will discuss the first two in this section. Nuclear heating is one of the main themes of this
thesis and will be discussed in detail in section 4. The equations are solved by a relaxation
method [63].
We use the thermal profile obtained by Gupta et al [19] using an iterative process, where
an initial profile is used to calculate the heat sources in the crust, which is then used to
recalculate the profile. The nuclear processes are not very sensitive to the thermal profile.
21

This approach is therefore sufficient for the goal of this thesis, which is to identify the critical
nuclear reaction sequences. Further iteration processes to obtain a self consistent thermal
profile should be done in future work.

3.2.1

Neutrino cooling

We include two neutrino cooling processes in addition to the energy lost by neutrinos emitted
in nuclear electron capture and beta decay processes discussed in section 4. Neutrino pair
bremsstrahlung is generated when electrons scatter off the Coulomb field of nuclei in the
crust [18] (Haensel et al. 1996):

e− + (A, Z) → e− + (A, Z) + ν ν
¯

(3.4)

The rate calculation of Haensel et al is used. We also consider Cooper pairing of neutrons
in the 1 S0 state using an analytical equation from [64].

n + n −→ n + n + ν + ν
¯

(3.5)

We neglect neutrino emission in the outer crust such as Plasmon neutrino emissivity since
these processes are not significant for the temperatures of interest here.

3.2.2

Heat Transport

Heat is transported by relativistic and degenerate electrons. The conductivity is given by
the Wiedemann-Franz law:
2
K = (π 2 /3)kB T ne /m∗ τ
e

22

(3.6)

where m∗ = EF /c2 (EF is the electron Fermi energy) is the effective electron mass and τ is
e
the relaxation time, kB is the Boltzmann constant and ne is the electron number density.

1/τ = 1/τee + 1/τeQ + 1/τep

(3.7)

where τee , τeQ , τep represents relaxation times for electron-electron, electron-impurity, and
electron-phonon scattering. Impurity scattering occurs when there is more than one species
of nuclei. Nuclei other than the dominat species can then be considered as an impurity in the
lattice structure of the crust, hampering the heat transport. Electron-impurity scattering
[65] often dominates the heat transport and is calculated according to Itoh & Kohyama 1993

K = 2.363 × 1017 T8 ρ6

Xi
×
Ai
(1 + 1.018[
i

1
2/3 ρ2/3 )
i (Zi /Ai )Xi ]
6

<S>×

<Z>
Q

(3.8)

where T8 is temperature T /108 , ρ6 is density ρ/106 , X is the mass fraction, A is the
mass number, < Z > is the average proton number, < S > is the structure factor which
determines how a material scatters electrons, Q is the impurity parameter.
Electron-electron scattering can be neglected as the strongly degenerate electrons restrict
the available phase space. For the electron-phonon relaxation time we use Potekhin et al.[66].

3.3

Nuclear Reactions

The main focus of the thesis is to study the composition and heat generation by nuclear
reactions as a function of depth in the neutron star crust. Rates for electron captures with
emissions up to 20 neutrons, beta decay, neutron capture, photodisintegration, and pycnonu23

clear fusions are used in our model. Because of the extreme environment with high density
(ρ ∼ 109 g cm−3 to 1013 g cm−3 ) and temperature (∼ 0.5GK ) it is impossible to obtain
the rates of these processes in the laboratory. Instead, rates are calculated theoretically.

3.3.1

Electron capture rates

Electron capture rates (EC) are calculated theoretically as in Gupta et al [19]. The electron
capture rate Rij between a parent nucleus in state i and a daughter nucleus in state j is
given by:
Rij = ln(2)fij (T, ρ, Ef )/(f t1/2 )ij

(3.9)

where fij (T, ρ, Ef ) is the phase space factor and (f t1/2 )ij is the comparative half-life. fij is
defined:
fij =

wl

w2 (q + w)2 G(Z, w)S− (w, µe )(1 − Sτ )dw

(3.10)

where the electron total energy, w, is in electron rest mass units (me c2 ) and q is the
nuclear mass difference between parent and daughter ground states in (me c2 ) units corrected by the excitation energies of the initial and final states to yield the channel Q-value:
q = qgs + (Ei − Ef )/me c2 . The capture threshold, ωl, is unity if q ≥ −1 and |q| otherwise.
S− (w, µe ) is the Fermi-Dirac distribution for relativistic degenerate electrons. It is very
sharply peaked for low temperatures. Tabulation of the phase space integral is therefore
difficult. Numerical integration, while challenging, is possible but computationally too expensive as it has to be carried for each step of the model because of the changes in electron
density. A fast analytic phase space approximation [67] for fij that is valid for low temperature is used to solve the problem. The integral is separated into two separate terms at
24

the chemical potential, thus eliminating the need to find the peak of the integrand. The
Coulomb correction G(Z,w), is assumed to be constant for the energy range considered here,
and is taken outside the integral. Using this approximation, the electron capture phase space
can be calculated for each timestep in a reaction network simulation.
The comparative half life (f t1/2 )ij is related to the nuclear transition matrix elements
for transition probability Bij by

(f t1/2 )ij = D/Bij , D = 2π 7 /g 2 m5 c4
e

(3.11)

where me is the mass of the electron, c is the speed of light, and g is the β decay strength
constant. Calculated electron capture strenghth functions were taken from the results of
[68]. These strength functions were calculated using a QRPA model. In the code, a table
of electron capture transition strength functions as a function of excitation energy in the
daughter nucleus is used. Temperatures are low and the parent nuclei are assumed to be in
their ground states. The neutrino losses for each transition are calculated. This is different
from earlier work where fixed fraction of (µe − Ethr ) lost to neutrinos was assumed [17].
Electron capture rates and up to 20 induced neutron emissions are included in the network. When the excitation of a daughter nucleus is greater than the neutron separation
energy, electron capture induced neutron emission is allowed. Here it is assumed that the
largest number of neutrons allowed to be emitted is emitted with a 100 % branching.
The above formalism is not developed for nuclei with Z < 8 and for some nuclei that
are beyond the neutron drip line. These nuclei are however encountered in some cases in
our network. For those nuclei, we assume a ground state to ground state transition, with a
log((f t1/2 )ij ) value of 6 and maintain the same fast phase space approximation. This log
25

(f t1/2 )ij value is a typical value for the dominant transitions in this mass region.

3.3.2

Beta decay rates

For beta decay rates, the same formulism as for electron capture rates is used. However,
instead of using the fast analytic phase space approximation, we use a table of rates as a
function of temperature and density. In addition, we set the rate to zero if the Fermi energy
is higher than the beta decay threshold.

3.3.3

Neutron capture rates

For two nuclear components i and j in an astrophysical plasma with number densities ni and
nj and a Maxwell-Boltzmann distribution of relative velocities the number of reactions, r,
per (cm3 s) is given by the following equation:

r = ni nj σν

(3.12)

This form assumes neutrons and ions have a Maxwell-Boltzmann velocity distribution.
This might not be the case in the deeper crust. Below a correction for the calculation of
reverse rates under neutron degenerate conditions is discussed, which has been used in this
work. However, we neglect the effect of neutron degeneracy on the forward neutron capture
rates. This is justified as these rates have very high theoretical uncertainties.

σν = (

1
πµ1/2 (kT )2
8

)

σ(E)E exp(

−E
)dE
kT

(3.13)

where E is the energy and µ is the reduced mass. The neutron capture rates are calculated
26

assuming a Maxwell-Boltzmann velocity distribution. Target states µ in an astrophysical
plasma of temperature T are thermally populated and the astrophysical cross section σ is
given by
σ(Eij ) =

µ
µ
µν
µ (2Ji + 1)exp(−Ei /kT ) ν σ (Eij )
µ
µ
µ (2Ji + 1)exp(−Ei /kT )

(3.14)

Photodisintegration reactions are the reverse reactions of neutron capture. They are
derived using detailed balance [70]. The commonly used equation is valid only if the neutrons
are non-degenerate:

λγ = (

Ai Aj 3/2 (2Ji + 1)(2Jj + 1) Gi (T )
)
× (T )3/2 F e−Q/kT < σi ν >∗
Am
(2Jm + 1)
Gm (T )

(3.15)

where λγ in s−1 is the photodisintegration rate, temperature T9 = T /109 K, A is the mass
number and NA < σν >∗ is the neutron capture rate in cm3 s−1 mol−1 . F is a constant.

3/2

(T )3/2 F = T9

9.8685 × 109 molcm−3 .

(3.16)

G(T ) is the temperature-dependent partition function given by

µ

(2Ji0 + 1)G(T ) =
µ

(2Ji + 1)e

µ
−Ei /kT

+

(2J µ + 1)e−ǫ/kT ρ(ǫ, J µ , π µ )dǫ

(3.17)

with ρ being the level density and µm the last included experimentally known state and Ei
are the energy of excited states.

The reaction rates are represented in our code as seven parameter ’reaclib’ fits using the
27

expression [71]

−1/3
1/3
5/3
−1
NA σν ∗ = exp(a0 + a1 T9 + a2 T9
+ a3 T9 + a4 T9 + a5 T9 + a6 lnT9 )

(3.18)

This parameterization is flexible enough to accommodate the different temperature dependencies of the various reaction types across the fitted temperature range 0.01 ≤ T9 ≤ 10.

We use radiative direct capture for the neutron capture rates since the nuclei we are
interested in our calculation have low neutron separation energy. The level density may
be so low that there are no or only very few compound nuclear states availabe for the
neutron to be captured to. Thus the Hauser-Feshbach formalism is not a good approach
here. Neutron capture occurs predominantly via a direct electromagnetic transition, direct
radiactive capture (DRC) [69], to a bound final state. Following Mathews [69], the DRC
cross section σn,γ can be approximated as:

σn,γ =

(2J + 1)
Re2
Y +3 2
8π Z 2
ξ
×(
) Sj (Qj + E)
3 E)1/2 c3 (2l + 1)(2s + 1)(2I + 1)
3 A (2µ
Y +1

(3.19)

where E is the incident neutron energy on a nucleus with mass A, and I is the spin of the
target nucleus, J is the spin of the final state, R = 1.35A1/3 is the effective sharp nuclear
radius, and µ is the reduced mass in the entrance channel. The factor ξ is the sum of incident
channel spins which lead to the same final spin J .SJ is the spectroscopic factor for the final
state which we set to 1 to obtain a crude approximation, QJ is the binding energy of the
final state, and the dimensionless parameter Y is

Y =

[2µ(QJ + E)]1/2

R = 0.295[µ(QJ + E)]1/2 A1/3
28

(3.20)

Typically, for (n, γ) and (γ, n) reactions, the energies ω of electromagnetic transitions
can be lower than ωp where ωp is the plasma frequency. It has been suggested that (γ, n)
rates may be suppressed because the γ rays will be screened by the electron plasma [19, 20].
For example, the plasma energy is ω ∼ 1.3 MeV at ρYe ∼ 1011 g/cm3 and 4.1 MeV at
1012 g/cm3 . This is larger than the typical neutron separation energy Sn which is ∼ 1.4
MeV. However, recent theoretical work demonstrated that this is not the case. [72]. Even
though the emission of real transverse and longitudinal plasmons will be prohibited when
ω < ωp , the processes with virtual plasmons (especially the longitudinal virtual plasmon)
strongly increase the radioactive decay [72]. Thus the de-excitation due to electron capture
will not be prohibited even when the plasma frequency is high. Consequently plasma effects
cannot suppress neutron capture reactions (n, γ) at ωp > ω. At the temperature range that
our model concerns, there will be a high level of fluctuating plasma micro fields (associated
with virtual plasmon) to power the inverse reaction (γ, n) and the associated rates will be
modified. However, the rates are very uncertain due to lack of nuclear structure information
on the exotic nuclei in question. Thus we decided to ignore the problem for now and factor
it into the (n,γ) rate uncertainties. Further work should be done in the future to improve
the model to include the plasma enhancements properly.
However, the velocity distribution of the neutrons or ions are not a Maxwell-Boltzmann
distribution anymore when the neutron density is too high or the temperature is too low.
The neutron gas becomes degenerate when the thermal energy is less than the Fermi energy
2
3
kB T <
(3π 2 n)2/3
2
2m

(3.21)

where m refers to the mass, and n to be the number density of neutrons. Assume a temper29

ature of around 0.5 GK, the neutron number density should be larger then ∼ 1034 cm−3 .
This regime is reached quickly at a column depth (pressure over gravity) 8.0 × 1015 g cm−2 .

In order to address this problem, we derive the detailed balance equation for degenerate
neutrons. For a Boltzmann gas, the relationship between number density n and the chemical
potential is [73]
n=

(2J + 1)eµ/kB T (2πmkB T )3/2
h3

(3.22)

where n is the number density and µ is the chemical potential, J is the angular momentum,
kB is the Boltzmann constant, T is the temperature, m is the mass, h is the Planck constant.
Consider a nuclear reaction involving a photon with photon energy Eγ ,

X A + n −→ γ + X A+1

(3.23)

The relationship between number density n and chemical potential in the approximation of
an ideal Boltzmann gas is:
(2πmA kB T )3/2 µ /k T
nA = (2JA + 1)
e A B
h3

nA+1 = (2JA+1 + 1)

(2πmA+1 kB T )3/2 µ
e A+1 /kB T
3
h

(3.24)

(3.25)

Furthermore, the chemical potentials is in equilibrium and have the relationship:

µA + mA c2 + µn + mn c2 = µA+1 + mA+1 c2

(3.26)

where for totally degenerate neutrons µn is equal to the Fermi energy EFn . The rest masses
30

are related to the reaction Q-value Eγ as mA c2 + mn c2 − mA+1 c2 = Eγ . Thus combining
the 3 equations above one obtains

A 3/2 −Eγ −µn
2(2JA + 1)
nA nn
nn
e
(
) e
=
nA+1
(2JA+1 + 1) A + 1

(3.27)

Comparing with the standard Saha equation, which assumes a Maxwell-Boltzmann distribution for the velocities of all particles:
nA nn
A 3/2 −Eγ 2π 2 −3/2
2(2JA + 1)
(
(
) e
)
=
nA+1
(2JA+1 + 1) A + 1
mn kB T

(3.28)

one obtains a correction factor:

f=

nn
2(mn kT )(2π 2 )−3/2

exp(−µn /kT )

(3.29)

where nn is the neutron density, and mn is the neutron mass. In equilibrium the ratio of
the (n,γ) rate NA < σν > to the (γ,n) rate λγ ratio is then for fully degenerate neutrons:
NA < σν >
n nn
= A f
λγ
nA+1

(3.30)

Thus, in our model, we correct (γ, n) rates by the factor f given by Eq.3.29. For the neutron
chemical potential, we use the classical chemical potential for low neutron density and the
neutron Fermi energy for high neutron density. For the neutron chemical potential in the
intermediate region between low neutron density and high neutron density, we use a fitting
curve. Figure 3.1 shows the neutron chemical potential versus the neutron density (m−3 ) in
log scale for a temperature of 0.5 GK. After the fitting, the code is rerun and calculates the
31

new thermal structure.
Here we just take into account neutron degeneracy to calculate the ratio of the neutron
capture rates to (γ,n) rates. Further work on the calculation of neutron capture rates under neutron degenerate conditions is still required for more accurate results. However, the
neutron capture rates are quite uncertain. The nuclei involved are far from stability and
nothing is known about them. In addition, the neutrons can form Cooper pairs resulting in
superfluidity. Below a critical temperature, these Cooper pairs condense into the lowest energy single particle state which form a superfluid. The difference between Cooper pairs and
bosons is that Cooper pairs only form below the critical temperature. Microscopic calculations in uniform nuclear matter predict a critical temperature of the order of T ∼ 109 − 1010
K which is above the temperature in our model. This might also alter neutron capture rates
and could be taken into account in future work.

3.3.4

Pycnonuclear fusion

Pycnonuclear fusion is fusion of two nuclei due to high density. Due to the uncertainty
principle, nuclei still have a zero point energy even at zero temperature. Thus, they have
oscillations even at zero temperature. At high density, the nuclei are so closely packed
that they vibrate and overlap with each other and fuse. Pycnonuclear fusion occurs in the
crusts of accreting neutron stars. In this work we employ the pycnonuclear fusion rates
for a multicomponent plasma derived by Yakovlev et al [74]. In this approach, one of the
difficulties in the calculation of these rates is the need for reliable astrophysical S factors
at astrophysical energies .These energies are low compared to the presently accessible range
of laboratory fusion experiments. The S factor is therefore calculated theoretically. At low
32

neutron chemcial potential (MeV)

5

0

-5

-10

-60

-40

-20

0

20

40

neutron density log(m^(-3))
Figure 3.1: The neutron chemical potential versus the neutron density (m−3 ) in log scale

33

energies, σ(E) is dominated by the l=0 (s-wave) channel. The relationship between the
fusion cross section σ(E) at low energie, and the astrophysical S factor is then

S(E) = σ(E)Ee(2πη)

where η = (Z1 Z2 e2 / )

(3.31)

µ/(2E) is the Gamow parameter, Z1 and Z2 are charge numbers of

the nuclei. This parameterization removes from the fusion cross section the strong nonnuclear
energy dependence associated with Coulomb barrier penetration.

Another difficulty is related to the plasma physics. One needs to consider the Coulomb
fields of the surrounding plasma particles to calculate the Coulomb barrier penetration.
These Coulomb fields modify the reaction rates and lead to a screening effect that tends
to enhance the reaction rate. One distinguishes five nuclear burning regimes (two thermonuclear regimes, with weak and strong plasma screening respectively; two pycnonuclear
regimes for zero-temperature and thermally enhanced burning; and an intermediate regime).
Previous model work of the inner crust was limited to a single component plasma where
pycnonuclear fusion only occurs among the same nuclides [14, 18]. In the most recent study
[18] pycnonuclear fusion rates were estimated using the formulae of [75]. In order to incorporate the reaction rates into the crust model, the phenomenological expression of [21]. is
used to describe the pycnonuclear fusion processes in these five regimes. In this approach
the fusion cross section is calculated using a one-dimensional barrier penetration formalism
and the S˜o Paulo potential to describe the real part of the nuclear interaction VN (r, E):
a

VN (r, E) = VSP (r, E) = VF (r)exp(−4v 2 /c2 ).
34

(3.32)

VF (r)is the density-dependent double-folding potential, c is the speed of light, E is the
particle collision energy (in the center of mass reference frame), v is the local relative velocity
of two nuclei 1 and 2,

v 2 (r, E) = (2/µ)[E − VC (r) − VN (r, E)]

(3.33)

VC (r) is the Coulomb potential, and µ is the reduced mass. The barrier penetration model
calculates the reaction cross section σ(E) using the standard partial wave (l=0,1,...) decomposition of the effective potential Vef f (r, E). Coulomb(VC ), nuclear(VSP ) and centrifugal
potentials contribute to the potential:

Vef f (r, E) = VC (r) + VSP (r, E) +

2 l(l + 1)

2µr2

.

(3.34)

S(E) can be extrapolated to lower energies relevant to stellar burning. Following [21], the
energy dependence of S(E) is parameterized by an universal 9-parameter analytic expression:

S(E) = exp(B1 + B2 E + B3 E 2 +

C1 + C2 E + C3 E 2 + C4 E 3
)
1 + exp[(Ec − E)/D]

(3.35)

where E is the center-of mass energy of reacting nuclei in MeV, and EC , D; B1 , B2 , B3 ,
C1 , C2 , C3 and C4 are nine fit parameters. The parameter EC is approximately equal to
the height of the Coulomb barrier. It divides the energy range into that below the barrier
(E ≤ EC ) and that above the barrier (E ≥ EC ), where S(E) behaves distinctly different.
The parameter D characterizes a narrow energy width of the transition region between
the energy ranges below and above the Coulomb barrier. The parameter B1 determines
35

S(0). The parameters B2 and B3 specify the energy dependence of S(E) at E ≤ EC .
The parameters C1 , C2 , C3 , C4 ,with B1 , B2 and B3 contribute to the S(E) dependence at
E ≥ EC .
A unified phenomenological expression for the temperature and density-dependent reaction rate, which combines all five burning regimes and assumes a uniformly mixed multicomponent plasma at low temperatures is expressed as [74]:

pycn

Rij (ρ, T ) = Rij

(ρ) + ∆Rij (ρ, T )

(3.36)

pyc

where Rij (ρ) is the density-dependent pycnonuclear reaction rate at zero temperature and
∆Rij (ρ, T ) is a correction taking into account the thermonuclear burning regime. The S
factor enters both of the above two terms. The reaction rates in the thermonuclear regimes
(with weak and strong plasma screening) can be calculated reasonably well. The reaction
rates in other regimes (zero-temperature and thermally enhanced pycnonuclear regimes;
intermediate thermo-pycnonuclear regime) are much less certain. They are very sensitive
to currently unknown microphysical correlation properties in a multi-component plasma (a
uniform mix, a regular crystalline lattice, a phase-separated matter, a matter with impurities
and defects) [74]. The expression above assumes a uniform mix. Other microstructures can
strongly change the reaction rates. For example, reactions in a regular multi-component
plasma can be strongly suppressed due to the absence of suitable nearby nuclei [74]. For
future work it is important to know the actual microstructure of the multi-component plasma
at low temperatures. For example, the availability of closest neighbors, local separations,
and oscillation frequencies of neighboring nuclei, especially in the presence of impurities and
lattice defects.
36

We neglect all these effects and just use Eq 3.35 and Eq 3.36 in our code. We use the
946 fusion reactions involving stable and neutron-rich isotopes of C, O, Ne and Mg which
have been parameterized according to Eq3.35 [21]. For the other pycnonuclear reactions, we
extrapolate these fusion rates.

3.4

Mass Models

Nuclear masses are important input parameters for the network calculation. They are used
to calculate the phase space for the electron capture and beta decay rates, and for the
calculation of (γ,n) rates from (n,γ) rates. They are also used to calculate the (n,γ) rates.
They affect energy generation and determine the position of neutron drip. Masses also
determine the Q-value of the electron capture and beta decay process which in turn strongly
affects the phase space (Eq 3.13) and the depth at which an electron capture reaction occurs
[76]. In our nuclear reaction network of ∼ 3000 isotopes, we use theoretical mass model
when experimental masses [77] are not available. For the mass model and strength functions
as well as deformations we use the Finite Range Droplet Mass Model (FRDM) [78], which
also provides the input for the calculation of electron capture and beta strength functions
providing some self-consistency. For Z ≤ 7 and very neutron rich nuclei that the FRDM
does not cover, we use the analytical Hilf [79] mass formula.

3.4.1

FRDM model

The Finite Range Droplet Mass Model (FRDM) [78] is based on a macroscopic nuclear
droplet model with added microscopic corrections for shell and pairing effects. The pairing and shell correction with a folded-Yukawa single-particle microscopic model. The total
37

nuclear potential energy can be written as

Epot (Z, N, shape) = Emac (Z, N, shape) + Es+p (Z, N, shape)

(3.37)

The values of nine constants are determined from a least squares adjustment to the ground
state masses of 1654 nuclei ranging from 16 O to 263 106 and to 28 fission-barrier heights.
The macroscopic part of the model is mainly a droplet model. The Yukawa-plus-exponential
model for the surface tension was incorporated into the droplet model to account for the
finite range of the nuclear force and an empirical exponential term was added to account for
compressibility effects for light nuclei and for other higher order effects.
The microscopic term is the shell-plus-pairing correction Es+p (Z, N, shape) which is the
sum of the proton shell-plus-pairing correction and the neutron shell-plus-pairing correction.
The shell correction is calculated by use of Strutinsky’s method. For the pairing correction
the Lipkin-Nogami version of the BCS method is used.
The rms error of the mass model is 0.669 MeV for the known masses, and it extrapolates
very well towards neutron rich nuclei, with an average error of 0.448 MeV for the region
N >65. However, there are not huge differences to other mass models such as HFB21 [80].
As an example Fig 3.2 shows the predicted Q values for the electron capture for the mass
chain A =56.

3.4.2

Hilf mass model

The Hilf mass formula is a semi-empirical atomic mass formula [79]. The mass M is separated
¯
into an average part M descried by the droplet model (DM) of Myers and Swiatecki [81] and
38

0

EC Q values (MeV)

FRDM
HFB
-10

-20

-30

-40
15

20

25

Proton number Z

30

Figure 3.2: The electron capture Q-value for FRDM and HFB21 mass models for mass
A = 56 as a function of proton number.

39

˜
a microscopic term M such that:
¯
˜
M =M +M

(3.38)

¯
˜
M (N, Z) = M (N, Z, θ) + Mpair (N, Z) + MW (N, Z) + M (N, Z, θ)

(3.39)

˜
where Mpair is the pairing term, MW is the Wigner term, and M is the shell correction. The
˜
shell correction term M is derived by bunching the average single-particle spectrum from
a Fermi gas model. Inner shells’ effects, changing magicities along neutron magic isotone
chains, damping of the off-Fermi-energy shell contributions, and nuclear deformation are all
included. The mass formula is an analytic expression with 50 parameters fit, and a rms error
of 670 keV.
In this thesis the mass formula is mainly needed for extremely neutron rich light nuclei.
Figure 3 shows the Masses for the Hilf mass model and the experiment mass with Z = 6.

3.5

Numerical Methods

The crust model code calculates composition and energy generation as a function of depth
in a steady state approximation. Following the fate of an accreted fluid element increasing
time means increasing depth and pressure. For each time step pressure is evolved according
to the accretion rate:
δP (t) = g mδt
˙

(3.40)

where P is the pressure, g is the surface gravity, and m is the mass accretion rate per
˙
unit area. For a temperature profile provided as input, the EOS routine (as mentioned in
chapter3.1) calculates the density taking into account electron, neutron and ion abundances.
40

Mass Excess (MeV)
70

Hilf mass
Experiment mass

50

30

10

10

15

20

25

30

Mass Number
Figure 3.3: TMass excess calculated with the Hilf mass model and from experiment for
carbon isotopes.

41

For the resulting temperature and density conditions the nuclear abundance changes are
calculated using a reaction network solver for the abundances Yi (t) provided by Hix and
Thielemann [82]. The energy released is calculated using:

dQ = dEnuc + dEe − dEν

where dEnuc =

(3.41)

i dYi BEi is the nuclear energy generation calculated from the abundance

change dYi and binding energy BEi of each isotope i. dEe = µe dYe where µe is the Fermi
energy of the electrons is the energy change due to changes in electron abundance. dEν is
the energy released as neutrinos from electron capture. Once the calculation is completed
and all heat sources are calculated, a new thermal profile is calculated using the thermal
model from Ed Brown described in section 3.1.

The nuclear reaction network solver by Hix and Thielemann [82] solves the coupled first
order differential equations for the abundance of each isotope as a function of time. The
rate equations can be divided into 3 main categories. One category is reactions involving a
single nucleus, such as electron captures, beta decay and photodisintegrations. For reaction
involving two nuclei, such as neutron capture and pycnonuclear fusion, the reaction rate
depends on the number densities of both target and projectile nuclei. There are also a few
three particle processes but they are not important in our network. The reaction network is
described by the following set of differential equations [49]:

dYi
=
dt

i
Nj λj Yj +

i
Nj,k ρNA < σν >j,k Yj Yk

(3.42)

where the three sums are over reactions which produce or destroy a nucleus of species i with
42

1, 2 reactant nuclei, respectively. λ represents the rate for a one body reaction and < σν >
represents the energy averaged astrophysical reaction rate of the two body or three body
reactions. The N i ’s can be positive or negative numbers that specify how many particles
of species i are created or destroyed in a reaction. The nuclear abundance depends on the
number density ni as Yi = ni /ρNA , where NA is Avogadro’s number. For a nucleus with
atomic weight Ai , Ai Yi is the mass fraction of this nucleus, and

Ai Yi = 1.

Because our network includes strong, weak and electromagnetic reactions, we have a
wide range of timescales. Systems whose solutions depend on a wide range of timescales are
termed stiff. Thus an implicit approach for solving the reaction network is essential. There
are explicit differencing schemes and implicit differencing schemes. An explicit differencing
scheme is
1
[Y (t + ∆t) − Yi (t)] = ri (Y (t))
∆t i

(3.43)

where ri defines the rate of change of species i. Explicit schemes require too small timesteps
for a practical calculation. An implicit differencing scheme is

1
[Y (t + ∆t) − Yi (t)] = ri (Y (t + ∆t))
∆t i

(3.44)

where the abundance change is calculated using rates evaluated with the new abundances.
The code employed here use an implicit differencing scheme which allows larger timesteps.
The equation above can be solved by the Backward Euler Method. ri can be expanded as:

ri (Y (t + ∆t)) = ri (Y (t)) +
j

where fi = dYi /dt
43

δfi
∆Yj + ...
δYj

(3.45)

Replacing ri (Y (t + δt)) by equation 3.44, we get
1
[Y (t + ∆t) − Yi (t)] = ri (Y (t)) +
∆t i

j

δfi
∆Yj + ...
δYj

(3.46)

δf

Writing the Jacobian matrix Jij = δYi and [Yi (t + ∆t) − Yi (t)] = ∆Yi :
j
1
∆Y = ri (Y (t)) +
∆t

Jij ∆Yj + ...

(3.47)

j

Rearranging and writing in vector/matri form we get

ri (Y (t)) =

∆Y
˜
− J ∆Y
∆t

(3.48)

Solving for ∆Y , we get:
1
∆Y = ri (Y (t))(
− J)−1
∆t

(3.49)

Here we use the subroutine DEGSV in Lapack to invert the Jacobian matrix. In our code,
we first determine a trial timestep ∆t using the following equation:

∆t =

f × ∆max × Y
˙
Y

(3.50)

where f is 1.00 × 10−2 and ∆max is set by the user to be 1.00 × 10−1 . In order to minimize
noise effects in the network, abundances lower than 1.0 × 10−30 are set to zero. In order to
avoid unnecessary tracking of very small abundances, only nuclei with abundances greater
than 1.0 × 10−6 are used for timestep determination. The smallest timestep allowed by the
program is time t × 1.0e × 10−15 s because of numerical precision. The timestep is also

44

further restricted to be at most 2 times the previous timestep.
The network is evolved to time t+tdel and mass convergence is tested. The mass fraction
Ai × dYi is compared with the previous timestep to check if the change is larger than
1.0 × 10−6 . If yes, the timestep is cut by half and the evolution is redone. If the mass
convergence test is passed, the next step is a test to ensure the code has not stepped over
the onset of an electron capture reaction. The electron capture rates have sharp thresholds
with the rates prior to threshold being very small (resulting in the code to choose very long
time steps) and with the rates beyond the threshold dramatically increasing. In order to
make sure the code did not step over the onset of an electron capture rate, we compute dy/dt
at the new time and ensure:
y(t + ∆t)
< ∆t × 10
y
˙

(3.51)

If the code fails in the above test, tdel will be cut in half, and the abundance evolution
will be redone again. The process above will be continued until all the tests are passed.

45

Chapter 4
Results
The initial composition for crust model calculations varies from system to system depending
on the composition of the accreted matter and the characteristics of surface burning processes
such as X-ray bursts and superbursts. As a starting point a calculation using an initial
abundance of 56 Fe is presented. This is not necessarily a realistic case as typically the initial
composition will be a mix of different isotopes. But it is a good study example to simplify
the situation and to understand the general picture of the fate of the ashes. In addition, the
results can be directly compared to several previous studies such as [17], which also used an
initial 56 Fe composition. Once 56 Fe is discussed, we investigate the crust processes for other
initial single species compositions. Finally, realistic superburst ashes are used to understand
the crust processes for more realistic multi species compositions.

4.1

56

Fe

The neutron star crust simulation was performed assuming an initial 56 Fe abundance with
a mass fraction equal to 1 (see Fig.4.1.). The evolution starts at time 1.00 × 106 s, a
46

Time:
Temp:
Density:
n-Abundan:
EF_e:
EF_n:
Max Flow:

1.000e+06 s
0.36 GK
3.77e+07 g/cm^3
6.59e+00
0.92 MeV
0.00 MeV
0.00e+00

20
50
8
28
20
Figure 4.1: The initial composition in the chart of nuclides. Abundances are shown on a
logarithmic color scale from blue (lowest) to green to yellow to orange to red (highest) with a
lower cutoff at 1.0 × 10−15 . Thick black squares indicate stable nuclei and thin black squares
unstable but particle bound nuclie. Particle unbound nuclei don’t have an outline but can
still be part of the network. Nuclides that are not part of the network are shaded gray.
mass density of 3.77 ×107 g cm−3 , and an electron Fermi energy of 0.919 MeV. The initial
temperature is 0.36 GK. The initial mass density, electron Fermi energy and temperature
are determined by the choice of initial time. The initial time is directly related to the initial
depth at which the calculation begins. This initial depth was chosen such that it is well above
the point where 56 Fe begins to capture electrons.The simulation ends at time 1.00 × 1013 s
with a final temperature of 0.46 GK and a final mass density of 2.63 × 1013 g cm−3 near the
bottom of the crust (see Fig.4.2.). The electron Fermi energy is 35.42 MeV at this density.
For readers’ convenience, Figure 4.3 shows the relationship of different physical quantities
(time, pressure, column depth, density, electron Fermi energy).
At time 5.4 ×108 s, density 3.94 × 109 g cm−3 and EF ∼ 6.0 MeV, 56 Fe which has
47

Time:
Temp:
Density:
n-Abundan:
EF_e:
EF_n:
Max Flow:

5.806e+12 s
0.47 GK
1.79e+13 g/cm^3
9.33e-01
36.20 MeV
4.80 MeV
1.84e-05

20
50
8
28
Figure 4.2: Final composition. The lines denote “instantaneous” net reaction flows for the
timestep displayed. Red lines show flows in the direction of increasing neutron number, blue
lines flows in the direction of decreasing neutron number. See. Fig. 4.1. for additional
information.

48

40
pressure
column depth
density
electron fermi energy

30

20

10

0
6

8

10

12

log time (s)

14

16

Figure 4.3: The relationship between time (s), log pressure (dyn cm−2 ), log column depth
(g cm−2 ), log density (g cm−3 and electron Fermi energy(MeV)

49

an electron capture threshold of 6.2156 MeV (Q values plus excitation energy) undergoes
electron capture to 56 Mn. Experimentally it was found that there is lower lying GT strength
in 56 Mn but as experimental information is not avaliable for the vast majority of transitions
this work uses only theoretical strength functions. 56 Mn has an electron capture threshold of
4.8096 MeV, which is lower than the 56 Fe electron capture threshold and further undergoes
electron capture to 56 Cr. Thus a two step electron capture occurs due to the odd-even mass
staggering of the nuclei. The electron capture occurs for a Fermi energy below the threshold
because the thermal energy ∼ kT excites part of the electrons to higher energies. Figure 4.4
shows the reaction flow of the electron capture from 56 Fe to 56 Cr. The flow F is calculated
dY

by time integration of the flux dti by the following equation:

Fi→j =

(

dYi
dY
− i
)dt
dt i→j
dt j→i

(4.1)

dY

where flux dti is calculated by the following equations:
For 1 reactant
dYi
= ρ < σν >1,i yi
dt i→j

(4.2)

where < σν > is the product of cross section and relative velocity averaged over the MaxwellBoltzmann velocity distribution. Y is the abundance and ρ is the density.
For 2 reactants,
dYi
= ρ < σν >2,i,j yi yj /δi,j
dt i→j

(4.3)

The flows shown in Figures 4.2. and 4.4. are instantaneous flows. That is, they are
integrated over one timestep and the timestep is determined by the code as explained in
section 3.4. The figure 4.2 also shows the main flows only while the weak flows are suppressed.
50

Time:
Temp:
Density:
n-Abundan:
EF_e:
EF_n:
Max Flow:

4.783e+08 s
0.50 GK
3.57e+09 g/cm^3
1.74e-40
5.62 MeV
0.00 MeV
1.23e-06

20
50
8
28
Figure 4.4: The two step electron capture process from 56 Fe to 56 Cr (see Fig.4.1 and 4.2 for
details)

At time 3.4 ×109 s, density 1.67 × 1010 g cm−3 , and EF ∼9.00 MeV, 56 Cr, which has an
electron capture threshold of 9.25 MeV (Q value plus excitation energy), undergoes electron
capture to 56 V. 56 V has an electron capture threshold of 7.143 MeV and undergoes electron
capture to 56 Ti. 56 V has a β − decay Q-value of 9.20 MeV. The electron Fermi energy is not
high enough therefore to block the β − decay and some 56 V β decay back to 56 Cr. This leads
to some URCA type cycling between 56 Cr and 56 V and the associated neutron emission
leads to some cooling at this location (the possibility of such an effect had been pointed out
by Gupta et al.)
At time 2.0 ×1010 s, density 7.02 × 1010 g cm−3 , neutron abundance 1.78 ×10−22 , and
EF ∼ 15.09 MeV, 56 Ti, which has an electron capture threshold of 15.75 MeV undergoes
electron capture to 56 Sc. 56 Sc, which has an electron capture threshold of 14.56 MeV under51

goes electron capture to 56 Ca. Figure 4.6 shows the abundance versus time for the electron
captures so far. Clearly the transition regions are rather narrow owing to the sharp turnon
of the electron captures, though the width is not negligible.
At time 8.6 × 1010 s, density 2.26 × 1011 g cm−3 , neutron abundance 9.81 ×10−20 , and
EF ∼ 21.79 MeV, 56 Ca undergoes electron capture to 56 K. 56 K then electron captures and
emits neutrons to 54 Ar while at the same time the released neutrons are captured by 56 Ca
producing 57 Ca, which then electron capture to 57 K and undergoes further electron capture
with neutron emission to 56 Ar. The simultaneous occurrence of neutron capture and electron
capture on 56 Ca happens because the immediately following electron capture on 56 K leads to
significant neutron emission. This leads to a buildup of a neutron abundance of 9.81 ×10−20 ,
leading to a significant neutron capture branch before all 56 Ca has had a chance to capture
electrons. The neutron captures are effective as 56 Ca lies on the neutron deficient side of
the (n,γ)-(γ,n) equilibrium point. The result is a split of the reaction path with one branch
moving to neutron rich Ca isotopes via neutron capture and the other branch moving to
lower Z via electron captures (and of course neutron captures as well). This is very different
from Haensel & Zdunik [17] where 56 Ca is directly converted to 56 Ar without a splitting of
the reaction path. The difference comes from our consideration of realistic finite electron
capture rates, as opposed to infinite fast electron capture in Haensel & Zdunik.
At time 1.4 ×1011 s, density 3.69×1011 g cm−3 , EF ∼ 23.61 MeV, and neutron abundance
3.0 ×10−17 , 52 Ar electron captures and neutrons emissions to 51 Cl and then electron capture
and neutron emissions to 48 S. At the same time 58 Ca further neutron captures to 60 Ca and
62 Ca. 58 Ca also has an electron capture branch ane neutron emissions to 57 K and then to
54 Ar, 56 Ar.

52

Abundance
0.1
0.01
0.001
0.0001
1e-05
1e-06
1e-07
1e-08
1e-09
1e-10
1e-11
1e-12
1e-13
1e-14
1e-15
1e+12

Fe56
mn56
cr56
sc56
ti56
sc56
ca56

1e+13

1e+15

1e+14

1e+16

1e+17

Column depth (g cm^(-2))

Figure 4.5: Abundance as a function of column depth for 56 Fe, 56 Mn, 56 Cr, 56 V, 56 Ti, 56 Sc,
56 Ca

53

Time:
9.638e+10 s
Temp:
0.51 GK
Density:
2.73e+11 g/cm^3
n-Abundan: 2.12e-18
EF_e:
22.44 MeV
EF_n:
0.00 MeV
Max Flow: 3.65e-07
20
50
8
28
Figure 4.6: The fluxplot shows the splitting of the path at time 9.638 × 1010 s. (See Fig.
4.1. and 4.2. for details)
At time 2.0 ×1011 s, density 4.70 × 1011 g cm−3 , neutron abundance 2.75 ×10−17 , and
EF ∼ 26.99 MeV, 48 S is partially converted by a sequence of electron captures with neutron
emission into 42 Si, with the remainder neutron capturing to 50 S.
At time 2.2 ×1011 s, density 5.12×1011 g cm−3 , EF ∼27.78 MeV, and neutron abundance
1.10 × 10−10 , 50 S also undergoes electron capture with neutron emissions ending up in 42 Si.
62 Ca finally also starts electron capturing into 56 Ar.

At time 3.0 ×1011 s, density 6.43 × 1011 g cm−3 , EF ∼ 30 MeV, and neutron abundance
2.08 ×10−08 , a small amount of 42 Si (EC threshold 30.39 MeV) does electron capture to
41 Al (EC threshold 29.72 MeV) and then to 36 Mg. Most 36 Mg is then converted by neutron

capture into 40 Mg because of the high neutron capture rate. However, there is again a
splitting of the path as electron capture can compete with neutron capture. Therefore some
36 Mg (EC threshold 26.85 MeV) does further electron capture to 34 Na and 35 Na. Within

54

the time 3.0 × 1011 − 3.1 × 1011 s, the net flux for the neutron capture of 36 Mg is 5.25 × 10−4 .
The flux for going to 34 Na, 35 Na are 1.33 × 10−4 and 1.479 × 10−4 respectively. The flux
towards lighter nuclei below Mg is therefore significant (or the order of 50%). 34 Na, 35 Na (EC
threshold 24.11 MeV, 25.72 MeV) do further electron capture to 31 Ne and neutron capture to
32 Ne. 32 Ne (EC threshold 30.35 MeV) undergoes electron capture to 31 F. 31 F (EC threshold

29.02 MeV) coverts to 24 O (EC threshold 27.93 MeV) and then into 23 N. 23 N (EC threshold
28.63 MeV) then undergoes electron capture to 18 C. 18 C (EC threshold 28.47MeV) does
electron capture to 17 B. Some of the 18 C neutron captures to 20 C. As we can see, some of
the heavier isotopes are converted all the way down to very light isotopes (17 B) within very
short range of electron Fermi energy (∼ 1MeV). This is the super electron capture cascade
also found in [20]. When nuclei reach low neutron separation energies electron capture
tends to occur into excited states of the daughter nuclei that can lie significantly above the
thresholds for multiple neutron emission. Thus, the nuclei undergo electron capture and
emit neutrons. We find neutron emission up to 7 neutrons to be significant at this early
stage. Because the neutron emission leads to the formation of a daughter nucleus with lower
electron capture threshold, an electron capture, again, accompanied by neutron emission,
immediately follows. This happens because the electron capture rate for such low threshold
nuclei is higher than the neutron capture rates. The process then repeats until neutron
capture wins (as is partially the case for 36 Mg, where a neutron capture branch leads to
40 Mg with a high electron capture threshold ending the cascade) or pycnonuclear fusion sets

in. As a result high Z nuclei can be converted to low Z nuclei and a large amount of neutrons
within a short range of electron Fermi energy. Thus, the nuclei now have Z values that are
low enough for pycnonuclear fusion to occur. 17 B, 18 C, 23 N, 24 O,32 Ne fuse to form heavier

55

elements. This is a major difference to Haensel and Zdunik, where fusion of neon at much
greater depth (4.379 × 1012 g cm−3 as opposed to 6.40 × 1011 g cm−3 here) are the first
pycnonuclear reactions to occur. In our calculation, the super electron capture cascade effect
[Gupta] converts some fraction of the material into very low Z isotopes at relatively shallow
depths. Because of the lower Coulomb repulsion those nuclei can readily undergo fusion.
Figure 4.7 shows the reaction flow from 56 Ar down to 18 C and 18 C neutron captures to 20 C
and it fuses with itself to form 40 Mg. Some 32 Ne also fuses with 20 C to form 52 S. These are
the first pycnonuclear fusions.
However, the super electron capture cascade does not covert all nuclei into lighter species.
Along the way towards lighter elements there are branchings where neutron capture can
compete with electron capture, moving some abundance into more neutron rich species
with high electron capture thresholds, where the reaction flow has to wait for the electron
Fermi energy to rise to include electron capture, or the neutron density to raise to shift the
equilibrium towards more neutron rich isotopes. These waiting points are 56 Ar, 62 Ar, 50 S,
42 Si, 40 Mg, and 32 Ne (see Fig 4.7).

One of the nuclei produced in the initial super electron cascade is 32 Ne. It is subsequently
destroyed by fusion with lighter isotopes. So not only does one have early fusion, one also
destroys neon causing the next major fusion step to be Mg fusion instead of Ne fusion as in
Haensel and Zdunik [17]. Because of this in our model Ne fusion is essentially skipped and
after the first initial fusion events, the major next pycnonuclear fusion is delayed compared
to Haensel and Zdunik.
At time 3.4×1011 s, density 7.69×1011 g cm−3 , EF ∼30.84 MeV, and neutron abundance
5.43 ×10−06 , 50 S neutron captures to 56 S while 42 Si neutron captures to 46 Si.
56

Time:
Temp:
Density:
n-Abundan:
EF_e:
EF_n:

3.055e+11s
0.50 GK
6.49e+11 g/cm^3
8.70e-10
30.11 MeV
0.00 MeV

50
8
28
Figure 4.7: The fluxplot shows the small amount of 56 Ar does electron capture down to 18 C
and the first pycnonuclear fusion. The red lines from low Z to high Z show pycnonuclear
fusion. (see Fig.4.1 and 4.2 for details)
At time 4.046 × 1011 s, density 9.35 × 1011 g cm−3 , neutron abundance 1.68 × 10−02 ,
and EF ∼ 32.32 MeV, some 62 Ar neutron captures to 64 Ar and some electron captures to
57 Cl followed by neutron captures to 61 Cl. 56 S neutron captures to 58 S.

At time 4.071×1011 s, density 9.40×1011 g cm−3 , and EF ∼ 33.69 MeV, 62 Ar has a large
electron capture rate which acts as a pump to pump away the equilibrium Ar isotopes. As a
consequence 64 Ar (γ, n) to 62 Ar which undergoes electron captures with neutron emissions
as well as neutron capture to 61 Cl. 58 S (γ, n) to 56 S and undergoes electron capture and 13
neutron emissions followed by neutron capture to 45 P. 45 P undergoes electron capture and 5
neutron emissions to 42 Si and neutron capture to 46 Si. At this point in time the composition
consists mainly of 46 Si, 61 Cl, and 40 Mg. 40 Mg is still unchanged because of its very large
threshold for electron capture and the very low neutron capture Q-value of 4.69 MeV, which
leads to an extremely small neutron capture rate of 2.06 × 10−45 . In fact, (n, γ)-(γ, n)
57

equilibrium is located towards more neutron rich isotopes, which however cannot be reached
by the system because of the small neutron capture rate.
At time 4.801 × 1011 s, density 1.06 × 1012 g cm−3 and EF ∼ 33.75 MeV, 40 Mg fuses with
40 Mg to form 80 Cr. This is the first major fusion event where the isotope that constitutes

the main composition is fusing with itself. 80 Cr is immediately converted to light nuclei
via a super electron capture cascade, with some fraction branching back to 40 Mg and some
fraction further being converted to lighter nuclei.
At time 5.512 × 1011 s, density 1.56 × 1012 g cm−3 , neutron abundance 3.35 × 10−1 ,
and EF ∼ 34.42 MeV, electron capture with neutron emissions converts 61 Cl into 44 S,
which neutron captures into 50 S. 50 S does electron capture to 43 P then neutron captures to
47 P.47 P then electron captures to 40 Si followed by neutron captures to 42 Si. At 42 Si another

branching occurs. Some material undergoes electron capture down to 41 Al and then to 36 Mg
which then neutron captures to 40 Mg. Some 42 Si however undergoes further neutron capture
to 46 Si.
At time 6.5 × 1011 s, density 1.75 × 1012 g cm−3 , EF ∼ 35.97 MeV, neutron abundance
3.43 × 10−1 , 46 Si then electron captures to 41 Al and then to 36 Mg, which neutron captures
to 40 Mg. At this point all the major waiting points have been converted into 40 Mg and free
neutrons. 40 Mg is a special waiting point because with N = 28 it has a closed neutron shell.
Its neutron capture Q-value is 4.69 MeV resulting in a essentially non-existent neutron
capture rate, and a very high electron capture threshold of 37.03 MeV. β decay is Fermi
blocked. Therefore it is essentially stable at this stage, and all other nuclei are eventually
converted to 40 Mg. This is very different from the Haensel and Zdunik [17] calculation,
which did not include shell effects. 40 Mg continues to fuse with itself to form 80 Cr which is
58

immediately converted into lighter nuclei by super electron capture cascades. The flow to
lighter nuclei branches in part into 40 Mg essentially forming a cycle, and in part continues
into the nitrogen and carbon region, where these isotopes then fuse mostly with 40 Mg forming
subcycles. During the process, neutrons are kept being emitted. Essentially one loop in the
40 Mg-80 Cr cycle converts one 40 Mg nucleus into neutrons. Figure 4.2 shows the abundance

distribution and reaction flows during this time. As the fusion reaction with 40 Mg is the
slowest reaction in the cycle, 40 Mg is the only isotope that builds up a significant abundance.
Figure 4.8 shows the most dominant abundance change of the isotopes in the deeper crust.
The branchings into different subcycles at 40 Mg can be analyzed quantitatively. During
1.0×1012 −1.0×1013 s, there is a flux for 4.528×10−3 down to 33 N and a flux of 2.447×10−3
fluxesfor 40 Mg fusion to form 80 Cr. The flux for 23 N to fuse with 40 Mg is 2.905 × 10−3 and
the flux for 20 C to fuse with 40 Mg is 1.00 × 10−3 . The 3 most dominant fusions overall are:

23 N +40 Mg →63 K

(4.4)

40 Mg +40 Mg →80 Cr

(4.5)

20 C +40 Mg →60 Ar

(4.6)

The pycnonuclear fusion paths are very different from P.Haensel & Zdunik(1990)’s paper
[17]. The 3 pycnonculear fusions that P.Haensel & Zdunik propose are 34 Ne + 34 Ne →68 Ca,
36 Ne +36 Ne →72 Ca and 48 Mg + 48 Mg →96 Cr. In addition, they propose that after each

fusion there is a successive step-wise electron capture sequence that transforms the nuclei
back to lower Z isotopes while we find an essentially instantaneous conversion back to the
initially fusing isotope (with some branchings to other isotopes). There are two main reasons
59

0.02
ar56
ar62
cl61
si46
mg40

Abundance

0.015

0.01

0.005

1e+15

1e+16

Column depth (g cm^(-2))

1e+17

Figure 4.8: The abundance change of isotopes 56 Ar, 62 Ar, 61 Cl, 42 Si and 40 Mg versus the
column depth. The decrease of 40 Mg at high column depth is because the continusous
depletion by the pycnonuclear fusion cycles converting 40 Mg into neutrons.

60

Fe (26)
Ca (20)
50
O (8)
28
20
Figure 4.9: Time integrated reaction flows for the complete calculation with initial 56 Fe
composition. See Fig 4.1. and 4.2. for details. The reaction flows towards lower neutron
number are here green instead of blue.

for the difference. One is that the path splitting leads to formation of more than one nuclear
species resulting in fusion reactions of unlike species. The other reason are shell effects
resulting in the extremely low neutron capture rate and extraordinary high electron capture
threshold of 40 Mg making it a major waiting point. Finally, the super electron capture
cascade, a consequence of taking into account non-equilibrium neutron emission processes
during electron capture, and the competition between electron capture and neutron capture,
lead to a much more rapid conversion of pycnonuclear fusion products into light nuclei.
The final abundance distribution at the end of our calculation (EF =36.20 MeV) is shown
in the Fig 4.2. The composition is dominated by neutrons with an abundance of 0.956.
The next dominant nuclide is 40 Mg with an abundance of 1.100 × 10−3 while 46 Si has an
abundance of only 4.689 × 10−9 .
Figure 4.8 shows the total reaction flows integrated over the entire calculation of the
crust.
Figure 4.10 shows the average charge number < Z >, and the average mass number
61

< A > (which excludes neutrons) of the composition versus column depth (g cm−2 ). At
the beginning, in our model, < A > remains unchanged while < Z > stepwise decreases
because of the electron capture processes. This is similar to Haensel & Zdunik [17] though
their transitions occur at different locations because they assume ground state-ground state
electron capture and because they use different nuclear masses. At a column depth of
8.0 × 1015 g cm−2 , < A > increases because of a significant neutron capture branch in the
reaction path. At column depth of 1.0 × 1016 g cm−2 , a sharp decrease in < A >, < N >
and < Z > occurs due to the super electron capture cascade effect which is similar to Gupta
et al for a 106 Pd case. Beyond a column depth of 1.0 × 1016 g cm−2 , < A >, < Z > and
< N > remain constant at < A >= 40, < Z >= 20 and < N >= 28 since the composition
remains mainly 40 Mg. This is a major difference to Haensel & Zdunik [17] where < A > and
< Z > continue to increase by a factor of 2 due to pycnonculear fusion and then decrease
due to electron capture and neutron emission. This pattern repeats for two cycles. However,
we find that due to shell effect 40 Mg dominates the composition for a long time. In addition,
after a pycnonuclear fusion reaction electron capture and neutron emission quickly process
material back to 40 Mg forming rapid cycles due to the super electron capture cascade effect.

Figure 4.11 shows the integrated heat energy deposited into the crust of 56 Fe versus
column depth (g cm−2 ). The heat energy per nucleon deposited into the crust is calculated
by equation 3.41. One heat source is electron capture reactions. To illustrate the heating
process Figure 4.12 shows the nuclear level scheme and transitions involved. An electron
is captured by the parent nucleus in its ground state (we neglect thermal excitations in
our approach as temperatures are relatively low) to some excited state in the daughter
nucleus. About 3/4 of the transition energy is carried away by the neutrino that is emitted
62

100

<Z>, <A>

80
60
40
20
0

1e+12

1e+13

1e+14

1e+15

1e+16

1e+17

Column depth (g cm^(-2))

Figure 4.10: < Z > (This work: blue line, Haensel & Zdunik: black circles) and < A >
(This work: green line, Haensel and Zdunik: red squares) versus column depth.

63

in this process. The excited daughter nucleus then deexcites to the ground state depositing
its excitation energy as heat into the crust. Taking into account excited daughter states
therefore reduces the neutrino energy loss significantly (Gupta et al) [19]. This is different
from the Haensel&Zdunik [17] who assume electron capture proceeds from ground state
to ground state which release much less energy because of the higher neutrino losses. We
therefore obtain more heating. At a depth of 2.0 × 1016 g cm−2 the integrated deposited
energy continues to increase smoothly and is not governed by individual electron capture
transitions anymore. This is caused by energy released in pycnonuclear fusion reactions and
the rapid electron capture sequences they induce. Our code shows a much higher energy
deposition into the crust because of different pycnonuclear fusion paths.
Fig 4.13 shows another quantity characterizing the composition of the neutron star crust,
the so called impurity parameter, which measures the diversity of nuclide charge numbers.It
is an important parameter setting the conductivity of the crust [12]

Qimp = n−1
ion

ni (Zi − Z )2

(4.7)

where nion is the total number density of ions, ni is the number density of nucleus i and Zi
is the charge number of nucleus i. Figure 4.13 shows the impurity parameter versus column
depth (g cm−2 ). In the figure, five peaks appear at around 1.0 × 1013 g cm−2 , 8.0 × 1013 g
cm−2 , 6.0 × 1014 g cm−2 , 3.04 × 1015 g cm−2 and 1.0 × 1016 g cm−2 . The first four peaks
are due to the electron capture process which leads to coexistence of two species of nuclei at
the same time. The fifth peak is due to the splitting of the reaction path. Beyond the fifth
peak, the composition is dominated by 40 Mg which leads to a very pure crust with a very
small impurity parameter.
64

Integrated heat deposited into the crust (MeV)

3
2.5
2
1.5
1
0.5
0

1e+13

1e+14

1e+15

1e+16

1e+17

Column depth (g cm^(-2))

Figure 4.11: The energy generation (MeV per nucleon) versus the column depth (g cm−2 ).
The black line shows the heat deposited into the crust by our model and the red circles shows
the table from P.Hasensl & Zdunik 1990 paper.

65

Electron capture for which the first transition “1” proceeds to the ground state of nucleus Z-1
but the second transition instead goes to an excited state “2” followed by a radiative de-excitation “3”
Figure 4.12: The level scheme of electron capture [19].

Another interesting difference between our model and Haensel & Zdunik [83] is the behavior of the electron Fermi energy. Figure 4.14 shows the electron Fermi energy versus
column depth. Up to neutron drip around 1 × 1016 g cm −2 electrons supply the pressure.
Therefore the electron Fermi energy is simply determined by the pressure (with the mass
density adjusting itself for compositional changes) and identical in both calculations. Beyond neutron drip Haensel & Zdunik [17] ’s electron Fermi energy continues to increase,
albeit at a decreasing rate, while in our model the electron Fermi energy remains more or
less constant at a value of around 36.5 MeV. This is because the density is defined by the
total pressure of the crust. The outer curst (lower density), is mainly supported by electron
degeneracy pressure. But the inner crust, is mainly supported by neutron degeneracy pressure. In our model, the neutron pressure (as well as the neutron abundance) constrains the
electron pressure in the inner crust. If the electron Fermi energy kept increasing beyond the
electron capture threshold of 40 Mg (which has an electron capture threshold at around 37
66

7

Impurity parameter

6
5
4
3
2
1
0

1e+13

1e+14

1e+15

1e+16

1e+17

Column depth (g cm^(-2))

Figure 4.13: The impurity parameter versus column depth (g cm−2 ) for an initial 56 Fe
composition.

67

electron Fermi energy (MeV)

50
40
30
20
10
0

1e+13

1e+14

1e+15

1e+16

1e+17

Column depth (g cm^(-2))
Figure 4.14: The electron Fermi energy versus column depth (g cm−2 ) for an initial composition of 56 Fe. The line shows the electron Fermi energy of our model while the red circles
show the electron Fermi energy calculated by Haensel & Zdunik [83].

68

MeV), 40 Mg would further electron capture. The subsequent electron capture and neutron
emissions would quickly proceed to carbon, which would then trigger a pycnonuclear fusion
cycle with remaining 40 Mg nuclei ending in 40 Mg, effectively converting a 40 Mg nucleus into
neutrons. This would lead to an increase in neutron abundance and a decrease in density
and therefore electron Fermi energy turning electron capture on 40 Mg off again. The electron
Fermi energy is therefore stablilized near the electron capture threshold of 40 Mg.

4.2

56

Fe without pycnonuclear fusion

In order to study the effect of pycnonuclear fusion on the model, we ran a calculation with
an initial composition of 56 Fe but turned off the pycnonuclear fusion rates. As expected,
instead of getting stuck as 40 Mg the reaction flow proceeds all the way down to the lightest
nuclei in the network.
Figure 4.15 shows the heat energy generation with and without pycnonuclear fusion versus
column depth (g cm−2 ). Clearly without pycnonculear fusion there is much less heating in
the inner crust.
By comparing the heat energy generation with and without pycnonuclear fusion, we can
see that the pycnonuclear fusion first sets in at around 8.0 × 1015 g cm−2 due to the fusion
of light isotopes such as carbon, nitrogen oxygen and neon isotopes. After 2.0 × 1016 g cm−2
there is a significant amount of heat deposition due to the pycnonuclear fusion of 40 Mg with
other light isotopes of carbon and nitrogen. We discussed the 3 main fusion cycles in the
previous subsection and they are the main heat sources in that region.
Figure 4.16 shows the electron Fermi energy versus column depth (g cm−2 ) for a calculation with and without pycnonuclear fusion. Without pycnonuclear fusion the stablilization
69

Integrated heat deposited into the crust (MeV)

3

2.5
2

1.5
1

0.5
0

1e+13

1e+14

1e+15

1e+16

1e+17

Column depth (g cm^(-2))

Figure 4.15: The integrated heat energy deposited into the crust with (red) and without
(black) pycnonuclear fusion versus column depth (g cm−2 ).

70

electron Fermi energy (MeV)

60
without pycnonculear fusion
with pycnonuclear fusion

50
40
30
20
10
0

1e+13

1e+14

1e+15

1e+16

1e+17

Column depth (g cm^(-2))
Figure 4.16: The electron Fermi energy with and without pycnonuclear fusion versus column
depth (g cm−2 ). The red line shows the calculation with pycnonuclear fusion while the black
line shows the calculation without pycnonuclear fusion

mechanism described above is disabled. Rather, nuclei at the lower end of the network are
essentially stable and can neither fuse or electron capture. Therefore, the electron abundances is fixed and the electron Fermi energy simply continues to increase with increasing
density.
71

4.3

Other abundances

In order to investigate the dependence of the final abundance on initial composition, the
model was run for a number of single species initial compositions. In each case, the initial
mass fraction of the initial species chosen is 1. The simulation conditions are the same as for
the 56 Fe case. Figures 4.17, 4.18 and 4.19 shows < Z >, < N > and < A > versus column
depth for calculations with 38 Ar,42 Si, 41 K, 47 Sc, 46 Ti, 50 Cr, 61 Ni ,70 Ge as initial abundance.

As for an initial composition of 56 Fe, initially during the electron capture phase < A >
remains constant while < Z > decreases and < N > increases stepwise. Overall the behavior for different initial isotopes is qualitatively similar with transitions simply occurring
at different locations owing to different electron capture thresholds. Neutron drip occurs at
around column depth 1.0 × 1016 g cm−2 with only a very small dependence on the initial
composition. And independent of the initial nucleus the composition is quickly converted
into 40 Mg. Therefore regardless of initial composition, beyond a column depth 2.0 × 1016
the crust is composed of 40 Mg.

38 Ar behaves differently as it reaches neutron drip at a lower Z (oxygen and neon isotopes)

than 40 Mg. Oxygen and Neon isotopes fuse to form heavier isotopes and they electron
capture and neutron emission down to 40 Mg. Figure 4.20 shows the time integrated reaction
flows for an initial composition of 38 Ar. As we can see from the graph, oxygen and neon
isotopes fuse to form heavier isotopes.
72

40
ar38
ge70
ni61
cr50
ti46
sc47
k41

35

<Z>

30
25

20
15
10

1e+12

1e+14

1e+16

Column depth (g cm^(-2))
Figure 4.17: < Z > versus column depth (g cm−2 ) for an initial abundance 38 Ar, 42 Si, 41 K,
47 Sc, 46 Ti, 50 Cr, 61 Ni and 70 Ge

73

50
ar38
ge60
ni61
cr50
ti46
sc47
k41
si42

45

<N>

40
35
30
25
20

1e+13

1e+14

1e+15

1e+16

1e+17

Column depth (g cm^(-2))
Figure 4.18: < N > versus column depth (g cm−2 ) for and initial abundance 38 Ar, 42 Si,
41 K, 47 Sc, 46 Ti, 50 Cr, 61 Ni and 70 Ge

74

80
ar38
ge70
ni61
cr50
ti46
sc47
k41
si42

70

<A>

60

50

40

30

1e+13

1e+14

1e+15

1e+16

1e+17

Column depth (g cm^(-2))
33
Figure 4.19: < A > versus column depth (g cm−2 ) for and initial abundance 38 Ar, 42 Si,
41 K, 47 Sc, 46 Ti, 50 Cr, 61 Ni and 70 Ge

75

Ca (20)
50
O (8)

28

Figure 4.20: Reaction flows for an initial composition of 38 Ar (see Fig. 4.1 and 4.2 for
details). Here unbound nuclei are not distingusihed from bound nuclei .

4.4

Superburst ashes

After studying the reaction path of individual initial nuclides, we now go to a more realistic
composition, the superburst ashes. This would reflect a system where all the accreted matter
is processed by superbursts. The superburst ashes from [84] are used. It is dominated by
66 Ni with a mass fraction of 0.348, 64 Ni with a mass fraction of 0.151, 60 Fe with a mass

fraction of 0.140, and 54 Cr with a mass fraction of 0.058. The evolution starts at time
1.0 × 108 s, temperature 0.48 GK, density 1.20 × 109 g cm−3 , and electron Fermi energy
3.66 MeV. The evolution ends at time 1.0 × 1013 s, with temperature 0.46 GK, and electron
Fermi energy 35.45 MeV. Figure 4.21 shows the initial composition and Fig 4.22 shows the
final composition at the end of the calculation. As expected from the results discussed in
4.3, the initially diverse ashes are again converted to 40 Mg at the end of the evolution.
At the beginning of the evolution, some nuclei do β decay back to more proton-rich nuclei
because of a mismatch of initial electron Fermi energy in our calculation, and the electron
Fermi energy corresponding to the initial composition. However, the composition quickly
settles into the correct equilibrium composition for the initial depth. The nuclei then undergo
electron capture and in some cases the inverse β − decay can also occur (green lines for some
76

Time:
Temp:
Density:
n-Abundan:
EF_e:
EF_n:
Max Flow:

1.000e+08 s
0.48 GK
1.20e+09 g/cm^3
2.82e-11
3.66 MeV
0.00 MeV
0.00e+00

50

28
20
50
8
20

28

Figure 4.21: The initial composition of the calculation with superburst ashes

77

2.992e+12 s
Time:
0.48 GK
Temp:
1.10e+13 g/cm^3
Density:
n-Abundan: 8.90e-01
36.24 MeV
EF_e:
3.71 MeV
EF_n:
Max Flow: 1.20e-04

50

28
20

8
20

28

Figure 4.22: Composition at the end of the calculation with initial superburst ashes(see Fig.
4.1 and 4.2 for details)

78

nuclei in Fig 4.23). This leads to a cycling behavior. The associated neutrino emission
leads to cooling. Beyond an electron Fermi energy of 12 MeV, neutron emission leads to
an increase in neutron abundance. The neutrons emitted are absorbed by nuclei (especially
even Z nuclei). The nuclei continue to be transformed into more neutron rich species by
electron capture and neutron capture while the neutron abundance keeps increasing.
At time 2.964 × 1011 s, density 6.56 × 1011 g cm−3 , neutron abundance 4.94 × 10−10 ,
and electron Fermi energy around 29.92 MeV, a super electron capture cascade as described
in Gupta et al [19] occurs. As a consequence high Z nuclei are converted into low Z nuclei
and a large number of neutrons are emitted within a small range of electron Fermi energy.
The situation is similar to the case of initial 56 Fe composition. One difference are the higher
abundances of 32 Ne, 42 Si and 70 Ca due to the different initial abundances in the calculation
with superburst ashes. Because of the higher abundance of 32 Ne, the dominant initial fusion
is 32 Ne fuses with 32 Ne to form 64 Ca and 23 N fuses with 23 N forming 46 Si as well as the
fusion of 23 N with 32 Ne to form 55 Cl. These first fusion reactions are different from the case
with an initial 56 Fe composition. Figure 4.23 shows the first fusion of the superburst ashes.
At time 3.957 × 1011 s, density 8.85 × 1011 g cm−3 , neutron abundance 4.84 × 10−3 , and
electron Fermi energy around 32.16 MeV, the dominant nuclei are 70 Ca, 62 Ar, 61 Cl, 46 Si and
40 Mg. The difference from the 56 Fe case is the appearance of significant amounts of 70 Ca

which is produced partly by electron capture and neutron emission from heavier isotopes
and partly by the fusion of neon isotopes as discussed above.
At time 5.153 × 1011 s, density 1.56 × 1012 g cm−3 , neutron abundance 3.78 × 10−1 ,
and electron Fermi energy around 34.00 MeV, 62 Ar and 61 Cl electron capture and neutron
emission partly to 46 Si and partly to 40 Mg similar to the 56 Fe calculation.
79

Time:
Temp:
Density:
n-Abundan:
EF_e:
EF_n:
Max Flow:

2.964e+11 s
0.50 GK
6.56e+11 g/cm^3
4.94e-10
29.89 MeV
0.00 MeV
1.39e-05

50

28
20
50
8
20

28

Figure 4.23: The first fusion of the superburst ashes

80

Ca (20)
50
O (8)
28
Figure 4.24: Time integrated reaction flows for an initial composition of superburst ashes
At time 6.689 × 1011 s, density 1.90 × 1012 g cm−3 , neutron abundance 3.83 × 10−1 , and
electron Fermi energy around 35.95 MeV, all isotopes convert to 40 Mg. Beyond this point,
the three most dominant fusion cycles are the same as in the 56 Fe calculation. Figure 4.24
shows the reaction flows.
Figure 4.25 shows the abundances of 72 Ti, 62 Ar, 70 Ca, 46 Si, and 40 Mg versus time. 40 Mg
starts building up at 3.0 × 1011 g cm−2 and starts decaying due to the fusion of 40 Mg at
1.0 × 1012 g cm−2 . The decay of 40 Mg occurs somewhat later than in the 56 Fe case.
Figure 4.26 shows < Z >, < N >, < A > versus column depth (g cm−2 ). At the
beginning, < N > increases and < Z > decreases while < A > remains the same. Then at
a column density of around 8.0 × 1015 g cm−2 , < Z > decreases while < N > and < A >
increase because of neutron capture. At around 1.0 × 1016 g cm−2 , < Z > and < N > and
thus < A > decrease because the super electron capture cascade. Beyond a column depth of
1.0 × 1016 g cm−2 the composition besides of neutrons is dominated by 40 Mg. An interesting
point is that the final abundance is the same as the 56 Fe case, this implies that the final
81

0.01

Abudance

ti72
ar62
ca70
si46
mg40

0.001

0.0001

1e+15

1e+16

1e+17

Column depth ( g cm^(-2))
Figure 4.25: The abundance of 72 Ti,62 Ar,70 Ca,46 Si,40 Mg versus column depth.

82

70
<Z>
<N>
<A>

60

<Z>, <N>, <A>

50
40
30
20
10
0

1e+15

1e+16

1e+17

Column depth ( g cm^(-2))

Figure 4.26: < Z >, < N >, < A > versus column depth.
composition does not depend on the initial abundance.
Figure 4.27 shows the integrated energy deposited into the crust as a function of column
depth. At the beginning, the energy is negative is due to the cooling via the β decay/electron
capture cycles. At column depth 5.0 × 1014 g cm −2 , heat generated by the electron captures
leads to an increase of the integrated energy deposition. After 1.0 × 1016 g cm−2 , there is a
jump in the energy deposition due to the onset of the super electron capture cascade. The
pycnonuclear fusion reactions further increase the heat energy deposition beyond that depth.
Figure 4.28 shows the impurity parameter as a function of column depth. After some
initial fluctuations the main effect is a drastic decrease in impurity around a column depth
83

Integrated heat energy deposited into the crust (MeV)

3
2.5
2
1.5
1
0.5
0
1e+12

1e+13

1e+14

1e+15

1e+16

1e+17

Column depth ( g cm ^(-2))

Figure 4.27: The integrated energy deposition versus column depth for initial superburst
ashes

84

of 1.0 × 1016 g cm−2 . This implies that the ashes transition from a diverse distribution to a
more uniform distribution. At the end of the evolution, the impurity parameter drops much
below 1 implying that the crust is very pure. The small increase in the impurity parameter
is because the 40 Mg keeps decreasing and converts to neutrons while the other abundances
(which are small) remain more or less constant.

85

Impurity parameter

10

1

0.1

0.01

0.001
1e+12

1e+13

1e+14

1e+15

1e+16

Column depth (g cm^(-2))
Figure 4.28: Impurity parameter versus column depth.

86

1e+17

Chapter 5
Astrophysical applications
The composition and heat generation of the crust of accreting neutron stars affect various
crust properties which are linked to observables. In this section, we use the results to briefly
discuss on the electrical and thermal conductivities, ionic coupling parameter, ionic and
electronic plasma frequency, plasma temperature, shear modulus and viscosity.

5.1

Ionic and electronic plasma frequency and temperature

Plasma properties affect the charge screening properties of nuclear reactions and the transmission of electromagnetic fields. In this context plasma frequency and plasma temperature
are two important parameters. The plasma frequency is defined as a resonant frequency
where the whole plasma behaves collectively as an oscillating system. For a solid state the
plasma frequency is an isotropic measure of the refractive index for transmission and reflection of radiation inside the plasma. The plasma in the magnetosphere of the neutron
87

stars can produce coherent radio emissions which are pulsars. In our calculation of the pycnonuclear reaction rates, we use the plasma frequency to take into account its effects on
the pycnonuclear reaction rates. The ion plasma frequency for a one component plasma is
defined as [10]
ωpi = (

4πe2 nN Z 2 1/2
)
Amu

(5.1)

where e is the electron charge, nN is the number density of ions, Z is the proton number of
the nuclei, A is the mass number of nuclei, and mu is the atomic mass unit. The plasma
temperature Tpi is defined as
Tpi = ωpi /kB

(5.2)

For multiple species of nuclei, the ion plasma frequency is defined as [74]
2
4πZj e2 nj

2
ωpi =

Aj mu

j

(5.3)

If the electrons are cold, (the electron thermal speed is negligible), the electron plasma
frequency is defined as [10]
ωpe = (

4πe2 ne 1/2
)
m∗
e

(5.4)

where ne is the electron density, m∗ = me (1 + x2 )1/2 and xr = pFe /(me c) and pFe =
e
r
(3π 2 ne )1/3
Figure 5.1 shows the ion plasma temperature (K) versus column depth (g cm−2 ) for an
initial composition of 56 Fe compared to previous work. The ion plasma temperature for our
model with initial superburst ashes is also shown. We can see that the plasma temperature
in general increases with the density. This is expected as the ions density increases with mass
88

density. The three cases are very close to each other in the outer crust. Beyond a column
depth ∼ 1016 g cm−2 , the plasma temperature increases much less for all three cases. This
is because neutrons drip out from the nuclei decreasing the density of ions. There are some
smaller discrepancies in the inner crust between our model and previous work because the
final abundance and therefore Z 2 /A are different. In addition, the electron density in our
model tends to remain constant in the inner crust while Haensel and Zdunik [17] predict a
continued increase (see Fig. 4.14) leading to different slopes in the plasma frequency.
Figure 5.2 shows the electron plasma temperature in log10 (K) versus column depth in
log10 (g cm−2 ) for 56 Fe compared to previous work. The ion plasma temperature for our
model with initial superburst ashes is also shown. The electron plasma temperature shows
similar behavior as the ion plasma temperature.

5.2

Coulomb coupling parameter and crust melting

Using the composition predicted by Haensel and Zdunik, it has been shown that for neutron
star crusts with low thermal conductivity where electron-ion scattering dominates, alternating layers of liquid and solid material can form [59]. Following Farouki & Hamaguchi [61],
the crust is liquid when the ionic coupling parameter Γ ∼< 170, otherwise, the crust forms a
crystallized solid. Electron capture reactions lower the Z values which leads to a decrease in
Γ below 173. Pycnonuclear fusions of identical species in the model of Haensel and Zdunik
double the Z values and increase the Γ value above 173. That is why the crust forms alternating layers. The strain in the crust disappears in the melted layers because the liquid
does not have shear stress. This limits the mass quadrupole of the neutron star induced by
thermal perturbations to the electron capture rate (Bildsten 1998) [11].
89

plasma temeprature log(K )

9.5
Fe56
Superburst
Haensel & Zdunik

9

8.5

8

7.5
10

12

14

16

18

20

Column depth log( g cm^(-2))
Figure 5.1: The plasma temperature in log10 (K) versus column depth in log10 (g cm−2 )
from this work for an initial composition of 56 Fe (black) and superburst ashes (red) as well
as from Haensel and Zdunik (green)

90

electron temperature log(K)

11
Fe56
Superburst
Haensel & Zdunik

10.5

10

9.5

9
10

12

14

16

18

20

Column depth ( g cm^(-2))
Figure 5.2: The electron temperature log10(K) versus column depth in log10 (g cm−2 ) from
this work for an initial composition of 56 Fe (black) and superburst ashes (read) as well as
from Haensel and Zdunik (green)

91

The ionic coupling parameter for a crust composed of a single species of nuclei is defined
as:
Γ=

Z 2 e2 4π
( n )1/3 ;
kB T 3 N

(5.5)

where Z is the charge number, e is the electron charge, kB is the Boltzmann constant, T
is the temperature, and nN is the number density of nuclei. The exact value of Γ for the
transition from solid to liquid is not known, but it is predicted to be approximately ∼ 173.
For a multicomponent plasma, the ionic coupling parameter can be defined as [74]:

Γ=

Γj =
j

j

2
Zj e2 4π
( n )1/3 =
kB T 3 N

5/3 2
e

Zj
j

kB T

(

4π
ne )1/3
3

(5.6)

Figure 5.3 shows Γ versus pressure for an initial composition of 56 Fe. Γ increases with
the mass density as the electron density increases until a pressure of around 1030 dyn cm−2
at which point the super electron capture cascade leads to a significant decrease in electron
density and Z value. After this sharp drop, differences from the composition of Haensel &
Zdunik appear. The main difference is that Γ remains fairly constant in the pycnonculear
regime because the abundance is locked in 40 Mg due to shell effects, and because the super
electron cascade effect leads to rapid cycling of material back to 40 Mg once a fusion reaction
has occurred. Therefore there is no temporary high Z phase following a fusion as in Haensel
& Zdunik and consequently Γ stays low.
However, Γ does not go as low as in Brown [59] because shell effects favour 40 Mg as
opposed to lighter nuclei. Therefore, when taking Γ = 173 as transition the crust in my
model is solid beyond a pressure of ∼ 1028.5 dyn cm−2 .
The predicted values of Γ vary only slightly between 183 and 216 in the inner crust,
92

rather close to 173. Given the uncertainty in where exactly the liquid to solid transition
takes place, there could still be liquid to solid transitions.

5.3

Shear modulus

The shear modulus of the neutron star crust can affect crust oscillations [86]. Specific
frequencies are observed for these oscillations. The shear modulus is also an important
parameter for the calculation of star quakes [87]. Magnetic mountains are irregularities on
the neutron star surface [87]. They have been proposed to be created by matter that is
accreted onto the star and is channeled to the poles by a strong magnetic field. It has been
suggested that if the field is strong enough it can maintain a moutain at the poles that is
not flattened by the strong gravitational force. The typical height of magnetic moutians is a
few meters. Magnetic mountains on rapidly rotationg neutron stars emit gravitational waves
because of their large masses and high accelerations and the height of the mountains depends
on breaking strain of the neutron star crust [87]. An effective shear modulus µ assuming a
body-centered cubic polycrystal is calculated by Ogata & luchimaru [88]. They performed
directional averages over rotations of the Cartesian axes. At T=0, they obtain:

1
µ = (2b11 + 3c44 )
5

(5.7)

where
b11 = 0.0245nN (Ze2 )(

4π
n )1/3
3 N

(5.8)

c44 = 0.1827nN (Ze2 )(

4π
n )1/3
3 N

(5.9)

93

400
350
300

Gamma

250
200
150
100
50
0
24

26

28

30

Pressure log ( dyn cm^(-2))

Figure 5.3: The Γ versus pressure log10 (dyn cm−2 ).

94

32

Figure 5.4 shows the effective shear modulus for an initial composition of 56 Fe from
Haensel & Zdunik’s and our model. The effective shear modulus for an initial composition
of superburst ashes is also shown. For a mix of nuclei, we assume Z is the average Z in
equations 5.8 and 5.9. The shear modulus in general increases with mass density as the
density of ions increases. Beyond neutron drip at a density around ∼ 1012 g cm−3 free
neutrons reduce the ion number density leading to a weaker increase of the shear modulus
with depth. In all 3 cases the calculated shear moduli are similar in the outer crust but they
differ in the inner crust. In the model of Haesnel & Zdunik the shear modulus is a factor
of 2 to 3 higher than in our model beyond a column depth of 1.0 × 1016 g cm−2 . This is
because the higher Z and higher ion density predicted by Haensel and Zdunik.

5.4

Superfluidity

Neutrons in the neutron star can form a superfluid. Neutrons of opposite spins form pairs
with zero total angular momentum and behave like bosons. These Cooper pairs condense
into a lowest energy single particle state below a critical temperature.
Pulsar glitches and neutron star oscillations are believed to be the observational evidence
of superfluidity in neutron stars. Besides, quiescent luminosity also provides clues for neutron
superfluidity [58]. Brown et al (2009) models the cooling curves of KS 1731-260 and MXB
1659-29. They found that the light curve of a cooling crust can be described as a broken
power law transitioning a constant at later times. From the model, the fit suggests that the
heat capacity is low, implying that the inner crust is a superfluid. Thus, the location of
neutron drip is critical.
Figure 5.5 shows the neutron abundance in log10 versus column depth for an initial com95

shear modulus log ( erg cm^(-3))

30
29
fe56
superburst
Haensel & Zdunik

28
27
26
25
24
23
10

12

14

16

18

20

Column depth ( g cm^(-2))
Figure 5.4: The shear modulus in log10 (erg cm−3) versus column depth in log10 (g cm−2 )
from this work for an initial composition 56 Fe (black) and superburst ashes (red) as well as
from Haensel and Zdunik (green)

96

position of 56 Fe from Haensel & Zdunik and from our crust model. The neutron abundance
of the initial composition of superburst ashes is also shown. In our model, neutron drip
occurs later than in previous work. This is because some of the nuclei emit neutrons and
immediately some nuclei capture the neutrons emitted which postpone neutron drip. The
differences in the mass models also lead to different neutron drip. After 2.0 × 1016 g cm−2 ,
our models emit neutrons more efficiently than Haensel & Zdunik because of the more rapid
pycnonuclear fusion cycles. In Haensel & Zdunik nuclei fuse together to form heavier nuclei
and electron capture then down to heavier nuclei (heavier than 40 Mg) which leads to less
neutrons emission.

5.5

Electrical and thermal conductivities

The main carriers in the thermal and electrical transport processes in the crust are electrons. Electron-electron scattering is negligible because the electrons are degenerate [10].
Thus, electric and thermal conductivities in the crust are mainly set by electron-phonon and
electron-impurity scattering. The electronic conductivity affects ohmic diffusion timescales
and therefore magnetic evolution [12]. The thermal conductivity affects the cooling in transients during the off state as well as the thermal profile. In the relaxation -time approximation, the electrical (σ) and thermal (κ) conductivities are (Yakovlev & Urpin 1980) [89]
e 2 ne
m∗ ν σ
e

(5.10)

2
π 2 kB T ne
2m∗ νκ
e

(5.11)

σ=

κ=

97

Abundance (log)

1
Superburst
Fe56
Haensel & Zdunik
0.1

0.01

0.001

1e+17

1e+16

Column depth log ( g cm^(-2))
Figure 5.5: The neutron abundance in log10 versus column depth in log10 (g cm−2 ) from
this work for an initial composition of 56 Fe (red) and superburst ashes (black) as well as
from Haensel & Zdunik (green)

98

Assuming electron-electron scattering is negligible:

ν = νκ = νσ

(5.12)

where ν is the sum of the electron-phonon (νph ) and electron-impurity collision frequencies
(νimp ).
ν = νph + νimp

(5.13)

The electron-phonon and electron-impurity collision frequencies can be calculated as Urpin
& Yakovlev (1980) [68]:
2c
∆ 2 1/2
1
[1 + (
) ] ;
=
2k T
νph
3.5kT
13e B

(5.14)

with the ∆ is the Debye temperature:

∆ = 1.76 × 108 K(

ρ
2Z
)( 10
)1/2
A 10 gcm−3

(5.15)

For mixed nuclei, the Debye temperature is calculated approximately by averaging the Z
and A values.
1
νimp

=

13π 3 Z me c
8me e4 Q pFe

(5.16)

where Q is the impurity parameter and pF e is the Fermi momentum of the electrons. Figure
5.6 shows the electron-impurity scattering frequency for the most impure crust composition
studied here, the superburst ashes, and compares it to the electron-phonon frequency for the
same composition for a temperature of 0.5 GK. As we can see, the electron-phonon frequency
is around an order of magnitude higher than the electron-impurity scattering frequency. This
implies that electron-impurity scattering is negligibile in the crust of accreting neutron stars
99

that process all their accreted material through superbursts. But if we compare it to the
electron-phonon scattering at 0.01 GK, we can see that the electron-impurity scattering is
comparable with the electron-phonon scattering.
Figure 5.7 shows the electrical conductivity calculated from the composition and electron
density of Haensel & Zdunik and from this work for both, initial 56 Fe and initial superburst
composition. They are calculated for a temperature of 0.5 GK. As we can see, the electrical
conductivity increases with the mass density as the electron density increases. Beyond a
density of ∼ 1012 g cm−3 the conductivities diverge. This is because the Z/A values are
different in our model from Haensel & Zdunik. The electron densities are also different. In
our model, the electron density remains constant while in Hasenel and Zdunik the electron
density increases.
The thermal conductivities show a similar behavior as the electrical conductivities.

5.6

Shear viscosity

The shear viscosity of the crust might have a damping effect on the amplitude of r-modes
oscillations in rotating neutron stars [90]. Such damping has implications for the emission
of gravitational waves. Here we assume the crust is a solid plasma where the transport is
dominated by electrons. The electron shear viscosity can be expressed as [10]:

η=

ne pFe vFe
5νph

(5.17)

Figure 5.9 shows the electron shear viscosity for an initial composition of 56 Fe calculated from
the results of Haensel & Zdunik and from this work. The electron-shear viscosity calculated
100

19
v_ph-0.01GK
v-ph-0.5GK
v_imp

log (s^(-1))

18

17

16

15
14

12

14

16

18

20

Column depth log( g cm^(-2))
Figure 5.6: The electron-impurity scattering frequency and electron-phonton scattering frequency log10 (s−1 ) of superburst

101

Electrical conductivity log( s^(-1))
1e+23
Fe56
Superburst
Haensel & Zdunik
1e+22

1e+21

10

12

14

16

18

20

column depth log(g cm^(-2))

Figure 5.7: The electrical conductivity in log10(s−1 )versus column depth log10 (g cm−2 ) for
an initial composition of 56 Fe (black), Superburst ashes (red) and Haensel & Zdunik (green)

102

for an initial composition of superburst ashes is also shown. The shear viscosity increases
with mass density. This is because the electron density increases with mass density. There
are discrepancies between our models and Haensel & Zdunik’s results beyond column depth
of 1.0 × 1016 g cm−2 because Hasenel & Zdunik find higher electron densities.

103

shear viscosity log ( g cm ^(-1) c^(-1))

Fe56
Superburst
Haensel & Zdunik

1e+16

1e+15

1e+14

14

16

18

20

Column depth log (g cm^(-2))
Figure 5.8: The shear viscosity in log10( g cm−1 c−1 ) versus column depth in log10 (g cm−2 )
for initial abundance of 56 Fe,(black), superburst (red) and Hasenel &Zdunik (green)

104

Chapter 6
Conclusion & Discussion
The goal of the thesis is to understand the nuclear reactions and associated heat generation in
the crust of accreting neutron stars. This is the first network reaction calculation of neutron
star crust processes that uses a complete set of nuclear reactions, including electron capture
with neutron emissions, β decay, neutron capture, photodisintegration, and pycnonuclear
fusion. The main differences to the results from previous work are:
(1) Even for single species initial composition, the reaction path splits and multiple
nuclear species are present beyond neutron drip.
This is because we use finite and realistic electron capture rates. When nuclei undergo
electron capture and neutron emissions, the high neutron capture rates capture the neutrons
emitted and lead to path splitting.
(2) Pycnonucler fusion sets in earlier than predicted previously and during this earlier
phase involves lighter species.
This is because the super electron cascade effect which was found by Gupta et al [20] leads
to electron capture to production of lighter nuclei at a low density which allows pycnonuclear
105

fusions to set in at a lower depth.
(3) 40 Mg plays a special role. Regardless of initial composition all nuclei are converted
into 40 Mg beyond a column density of 1.0 × 1016 g cm−2 . Therefore 40 Mg is the main fuel
for the beginning of the main phase of pycnonuclear fusion
This is because we use a mass model considering shell effects. Because 40 Mg has a
neutron magic number N=28 and therefore a high electron capture threshold and a high
neutron separation energy. Neutron capture rates and electron capture rates are low.
(4) Pycnonculear fusion proceeds in cycles where the fusion reaction is the slowest reaction. Therefore, fusion is not followed by a layer of high Z nuclei but the composition in the
inner crust remains predominantly 40 Mg throughout.
This is because of the super electron capture cascade effect that leads to nuclei emit all
the way down to light nuclei expect 40 Mg which has a high electron capture threshold as
explained in (3).
(5) Neutron drip happens at a higher density than in previous work.
This is because some nuclei emit neutrons and other nuclei which have high neutron
capture rates immediately capture the neutrons emitted. This postpones the neutron drip
to a higher density. The difference is also caused by a different mass model we use.
As a result one finds significant differences to the predictions from previous work.
In this thesis we use an initial set of nuclear physics input. Although the nuclear physics
input is uncertain because the nuclei are very neutron rich and experimental data are unavailable, significant improvements are still possible in the nuclear physics. First, the neutron
capture rates we use here are calculated using a rather crude approximate analytical formula.
More sophisticated predictions of the dominant direct capture processes should be possible.
106

The neutron degeneracy and neutron superfluidity in the calculation of the neutron capture
rates are neglected. This is justified for an initial study as neutron capture rate calculations have large uncertainties, but it will be interesting to explore the effects of the neutron
degeneracy and superfluidity in the future. Although a correction factor for the photodisintegration rates is included to account for fully degenerate neutrons, a descriptions allowing
for arbitrary degrees of degeneracy could be implemented. Furthermore, the pycnonculear
reaction rates for odd-even nuclei are calculated by extrapolations from the even-even nuclei.
In principle, the S-factor of the odd-even nuclei should be calculated to give more accurate
pycnonuclear fusion rates. However, the pycnonuclear fusion rates themselves are also very
uncertain. The structure of the neutron star crust (whether it is body-centered cubic or face
centered cubic, lattice imperfections, distribution of nuclei etc) are expected to affect the
pycnonculear fusion rates greatly [74]. The uncertainties from a lack understanding of the
masses and structure of neutron rich exotic nuclei also affect the calculation. In our model,
we use the FRDM model, complemented with the Hilf mass for nuclei lighter than oxygen.
In the future, it would be preferable to use a consistent set of mass which covers the neutron
rich and light nuclei, though the accurate prediction of the masses of light nuclei remains a
major challenge in nuclear theory. There are also large uncertainties in the electron capture
rates, in particular because the first forbidden transitions are neglected. This could also be
improved in the future, should global predictive models of first forbidden strength become
avaliable.

107

BIBLIOGRAPHY

108

BIBLIOGRAPHY

[1] J.M. Lattimer. & M. Prakash., Science., 304:536, 2004.
[2] S.L. Shapiro & S A. Teukolsky., Black Holes, White Dwarfs and Neutron Stars: The
Physics of Compact Objects., 1983.
[3] P. Haensel et al., Neutron Stars., Springer, 2007.
[4] Kaspi et al., Compact stellar X-ray sources., Cambridge University Press, 2006.
[5] S.E. Thorsett & D. Chakrabarty., ApJ., 512:288, 2010.
[6] Kiziltan et al., arXiv:1011.4291., 2010
[7] Demorest et al., Nature. 467:1081-1083, 2010
[8] Cornelisse et al., A & A.,382:174, 2000
[9] Steiner et al., ApJ, 722:33, 2010
[10] N. Chamel & P. Haensel, Living Reviews in relativity, volume 11, 2008
[11] L. Bildsten., ApJ, 506:842, 1998
[12] E.F. Brown & L. Bildsten., ApJ, 496:915, 1998
[13] Ushomirsky et al., MNRAS, 319:902, 2000
109

[14] K. Sato., Prog. Theor. Phys., 62:957, 1979
[15] G. Bisnovatyi-Kogan & V. Chechetkin., Sov. Phys. Usp., 22:89, 1979
[16] Y. Vartanyan & N.Ovakimova., Soob. Byur. Obs., 49:87, 1976
[17] P. Haensel & J. Zdunik., A & A., 227:431, 1990
[18] P. Haensel et al., A & A., 314:328, 1996
[19] S.S. Gupta et al., ApJ, 662:1188, 2007
[20] S.S. Gupta et al., Phys. Rev. Lett. 101:23101, 2008
[21] M. Beard et al., Atomic and Nuclear Data Tables, 96:541, 2010.
[22] W. Baade, & F. Zwicky., Proc. Nat. Ac. Sci., 20 :254, 1934
[23] J.R. Oppenheimer, & G.M. Volkoff., Phys. Rev., 55:374, 1939
[24] A.G.W. Cameron., ApJ., 130:884, 1959
[25] A.Hewish et al., Nat., 217:709, 1968
[26] E. Schreier et al., ApJ., 172:L79, 1972
[27] H. J. Grimm et al., MNRAS., 339:793, 2003
[28] H. Schatz et al., ApJ., 524:1014, 1999
[29] Grindlay et al., ApJ., 205:L127, 1976
[30] L. Maraschi & A. Cavaliere., X-ray binaries and compact objects., 1977
[31] H. Schatz & Rehm., Nuclear physics A 777:601, 2006
[32] H. Schatz et al., Phys. Rev. Lett., 86:3471, 2001
110

[33] J.J.M. in’t Zand et al., A & A., 426:257, 2004
[34] E.F. Brown., ApJ., 614:L57-L60, 2004
[35] E. Kuulkers., Nuclear Physics B-Proceedings Supplements., 132:466, 2004
[36] S. Campana et al., The Astronomy and Astrophysics review., 8:279, 1998
[37] L. Cominsky et al., ApJ., 224:46, 1978
[38] L. Maraschi et al., Nat., 259:292, 1976
[39] N. E. White & F. E. Marshall., ApJ. 281:354, 1984
[40] J. K. Cannizzo et al., in Accretion Disks in Compact Stellar Systems., World Scientific
Publishing, 1993
[41] J. van Paradijs, ApJ., 464:L139, 1996
[42] M. Tavani & J. Arons.,ApJ., 477:439, 1997
[43] E. Brown et al., ApJ., 504:L95, 1998
[44] R. E. Rutledge et al., ApJ., 559:1054-1059, 2001
[45] S. Campana et al., A & A. 324:941, 1997
[46] E. M. Cackett et al., MNRAS., 372:479 2006
[47] E. M. Cackett et al, ApJ. 687:L87, 2008
[48] Sunyaev R., IAU Circ., 4839, 1, 1989
[49] Sunyaev R., IAU Circ., 5104, 1, 1990
[50] R. Wijnands et al., ApJ., 560:L159, 2001
[51] R. Wijnands et al., ApJ., 573:L45, 2002
111

[52] W. H. G. Lewin et al., IAU Circ., 2994, 1976
[53] L.Cominsky et al., ApJ., 270:226, 1983
[54] in’t Zand J. et al., IAU Circ., 7138, 1999
[55] R. Wijnands.et al., ApJ., 566:1060, 2002
[56] R. Wijnands et al., ApJ., 594:952, 2003
[57] R. Wijnands et al., Rev. Mex. Astron, Astrofis. Conf. Ser., Vol. 20, Compact Binaries
in the Galaxy and Beyound. UNAM, Mexico, 2004
[58] E. F. Brown & A. Cumming., ApJ., 698:1020, 2009
[59] E. F. Brown., ApJ., 531:988, 2000
[60] F. X. Timmes & F. D. Swesty., ApJ, 126:501, 2000
[61] R. Farouki,& S. Hamaguchi., Phys. Rev. E., 47:4330, 1993
[62] F. D. Mackie & G. Baym., Nucl. Phys. A., 285:332, 1977
[63] W.H. Press et al., Numerical Recipes in FORTRAN, Cambridge Univ. Press, 1992
[64] D. G. Yakovlev et al., A & A, 343:650, 1999
[65] N. Ioth & Y. Kohyama, ApJ, 404:268, 1993
[66] A.Y. Potekhin et al., A & A. 346:345, 1999
[67] A. D. Becerril Reyes et al., Symposium on Nuclear Astrophysics -Nuclei in the Cosmos
-IX, 2006
[68] P. M oller & J. Randrup., Nuclear Physics A., 514:1, 1990
¨
[69] G.J. Mathews et al., ApJ., 270:740, 1983
[70] Christian iliadis, Nuclear Physics of Stars, 2007
112

[71] T. Rauscher & F. -K. Thielemann., Atomic Data and Nuclear Data Tables., 75:1, 2000
[72] P. S. Shternin & D. G. Yakovlev., Phys. Rev. D., 79:123004, 2009
[73] C. J. Hansen et al., Stellar Interiors - Physical Principles, Strcuture, and Evolution,
1995
[74] D. G. Yakovlev et al, Phys. Rev. C., 74:035803, 2006
[75] E. E Salpeter & H. M. Van Horn, ApJ., 155,:183, 1969
[76] A.Estrade, Phys. Rev Lett., 107:172503, 2011
[77] Private Communication April 2011 by Georges Audi and Wang Meng
[78] P. M¨ller et al., Atomic Data and Nuclear Data Tables., 59:185, 1995
o
[79] E. R. Hilf et al., Proceedings of the 3rd International Conference on Nuclei Far from
stability., pages 142-148, 1976
[80] S. Goriely et al., Phys. Rev. 82:035804, 2010
[81] W. D. Myers & W. J. Swiatecki, Nuclear Physics. 81:1, 1966
[82] W. R. Hix and F.-K. Thielemann., Journal of Computational and Applied Mathematics,
109:321,1999
[83] P. Haensel. & J. L. Zdunik, A & A, 229:1,1990
[84] H.Schatz et al., Nuclear Physics A, 718:247, 2003
[85] B. M. Mirza, International Journal of theoretical physics, 48:723
[86] A.I. Chugunov & D.G. Yakovlev, Astronomy Reports,49:724
[87] C. J. Horowitz & Kai Kadau, Phys. Rev. Lett., 102:191102, 2009
[88] S. Ogata, & S. lchimaru, Phys. Rev. A,42:4867, 1990
113

[89] D. G. Yakovlev & V. A. Urpin, AZh, 57:526, 1980
[90] O.L.Caballero et al., Phys. Rev. C, 78:045805, 2008

114