STARQUAKES, HEATING ANOMALIES, AND NUCLEAR REACTIONS IN THE
NEUTRON STAR CRUST
By
Alex Thomas Deibel

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
Astrophysics and Astronomy — Doctor of Philosophy
2017

ABSTRACT
STARQUAKES, HEATING ANOMALIES, AND NUCLEAR REACTIONS IN THE
NEUTRON STAR CRUST
By
Alex Thomas Deibel
When the most massive stars perish, their cores may remain intact in the form of extremely
dense and compact stars. These stellar remnants, called neutron stars, are on the cusp of becoming
black holes and reach mass densities greater than an atomic nucleus in their centers. Although the
interiors of neutron stars were difficult to investigate at the time of their discovery, the advent of
modern space-based telescopes (e.g., Chandra X-ray Observatory) has pushed our understanding
of the neutron star interior into exciting new realms. It has been shown that the neutron star interior
spans an enormous range of densities and contains many phases of matter, and further theoretical
progress must rely on numerical calculations of neutron star phenomena built with detailed nuclear
physics input.
To further investigate the properties of the neutron star interior, this dissertation constructs
numerical models of neutron stars, applies models to various observations of neutron star highenergy phenomena, and draws new conclusions about the neutron star interior from these analyses.
In particular, we model the neutron star’s outermost ≈ 1 km that encompasses the neutron star’s envelope, ocean, and crust. The model must implement detailed nuclear physics to properly simulate
the hydrostatic and thermal structure of the neutron star. We then apply our model to phenomena
that occur in these layers, such as: thermonuclear bursts in the envelope, g-modes in the ocean,
torsional oscillations of the crust, and crust cooling of neutron star transients.
A comparison of models to observations provides new insights on the properties of dense matter that are often difficult to probe through terrestrial experiments. For example, models of the

quiescent cooling of neutron stars, such as the accreting transient MAXI J0556-332, at late times
into quiescence probe the thermal transport properties of the deep neutron star crust. This modeling
provides independent data from astronomical observations on the nature of neutron superfluidity
and the thermal conductivity of nuclear pasta.
Our neutron star modeling efforts also pose new questions. For instance, reaction networks
find that neutrino emission from cycling nuclear reactions is present in the neutron star ocean and
crust, and potentially cools an accreting neutron star. This is a theory we attempt to verify using
observations of neutron star transients and thermonuclear bursts, although it remains unclear if this
cooling occurs. Furthermore, on some accreting neutron stars, more heat than supplied by nuclear
reactions is needed to explain their high temperatures at the outset of quiescence. Although the
presence of heating anomalies seems common, the source of extra heating is difficult to determine.

To Carol

iv

ACKNOWLEDGMENTS

Throughout my journey to become an astrophysicist many people have inspired me and have generously given me their time. I would like to thank my undergraduate research advisor Monica
Valluri, who patiently guided me through my first research project and sparked my interest in numerical models. It is a pleasure to thank Andrew W. Steiner for being my research mentor early in
my graduate career as we navigated the complicated landscape of the magnetar interior. I am especially grateful for the guidance and encouragement of my doctoral supervisor Edward F. Brown,
which has been vital to my success.
There are many others I would like to thank that over the course of my graduate career have
offered their guidance, support, and collaboration. I thank Andrew Cumming for his collaboration
on neutron star thermal evolution and superburst ignition models. I thank Hendrik Schatz and
Zach Meisel for expanding my knowledge of nuclear physics. I thank Sanjay Reddy and Bob
Rutledge for useful advice. I thank Dany Page and Duncan Galloway for their invitations to to the
International Space Science Institute teams on neutron stars. I also thank the International Space
Science Institute for supporting me as a young scientist at their international meetings. I also
thank the Joint Institute for Nuclear Astrophysics Center for the Evolution of the Elements whose
support allowed me to present my research around the world. I would like to thank members of
my guidance committee: Laura Chomiuk, Sean Couch, and Jim Linnemann, for their perspective
and guidance.
During the course of my research for this dissertation I received support from the National
Aeronautics and Space Administration through Chandra Award Number TM5-16003X issued by
the Chandra X-ray Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of the National Aeronautics and Space Administration under contract
v

NAS8-03060. Much of this work was supported by the National Science Foundation under grant
No. PHY-1430152 (Joint Institute for Nuclear Astrophysics - Center for the Evolution of the Elements) and grant No. AST-1516969. I also received support from the Michigan State University
College of Nature Science Dissertation Completion Fellowship.

vi

TABLE OF CONTENTS

LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

ix

LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

x

Chapter 1

1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Chapter 2
The Neutron Star Ocean and Crust
2.1 Structure . . . . . . . . . . . . . . . . . . .
2.2 Low-density Equation of State . . . . . . .
2.2.1 Nuclei and the Ion-lattice . . . . . .
2.2.2 Quasi-free Neutrons . . . . . . . . .
2.3 Composition . . . . . . . . . . . . . . . . .
2.3.1 Equilibrium . . . . . . . . . . . . .
2.3.2 Accreted . . . . . . . . . . . . . . .
2.4 Thermal Transport . . . . . . . . . . . . . .
2.4.1 Outer Crust . . . . . . . . . . . . .
2.4.2 Inner Crust . . . . . . . . . . . . .
2.4.3 Heat Sources and Sinks . . . . . . .

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

7
7
11
12
15
16
16
17
20
21
24
30

Chapter 3
Magnetar Starquakes . . .
3.1 A Magnetized Crust . . . . . . . . .
3.1.1 Electron Quantization . . . .
3.1.2 Composition . . . . . . . .
3.2 Crust Torsional Modes . . . . . . .
3.3 Constraints from Observations . . .
3.3.1 Nuclear Symmetry Energy .
3.3.2 Neutron Entrainment . . . .
3.3.3 Crust-core Transition Depth
3.3.4 Crust Magnetic Field . . . .

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

32
33
33
35
37
42
42
45
49
49

Chapter 4
Neutron Star Transient Thermal Evolution
4.1 A Cooling Slab . . . . . . . . . . . . . . . . . . . .
4.2 Realistic Cooling Models . . . . . . . . . . . . . . .
4.2.1 Heating Anomalies in MAXI J0556-332 . . .
4.2.2 The Reflare in MAXI J0556-332 . . . . . . .
4.3 Observational Tests for Ocean g-modes . . . . . . .
4.3.1 Neutron Star Ocean Thermal Evolution . . .
4.3.2 g−mode Spectrum of the Neutron Star Ocean
4.4 Late Time Cooling and Nuclear Pasta . . . . . . . .

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

52
53
55
58
64
66
69
72
78

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

vii

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

Chapter 5
Urca Cooling Nuclei Pairs . . . . .
5.1 Urca Pair Formation . . . . . . . . . . . . .
5.2 Crust Urca Pairs . . . . . . . . . . . . . . .
5.2.1 Impact on Neutron Star Transients .
5.2.2 Impact on Neutron Star Superbursts
5.3 Ocean Urca Pairs . . . . . . . . . . . . . .
5.3.1 Observational Signatures . . . . . .
Chapter 6

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

85
86
89
90
92
99
103

Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110

APPENDICES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Appendix A: Local Stability of Non-equilibrium Nuclear Reactions .
Appendix B: Thermal Diffusion . . . . . . . . . . . . . . . . . . .
Appendix C: Axial Perturbations . . . . . . . . . . . . . . . . . . .

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

114
115
117
120

REFERENCES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122

viii

LIST OF TABLES

Table 2.1:

Parameters of the mass model. . . . . . . . . . . . . . . . . . . . . . . . 15

Table 3.1:

Equilibrium nuclei below the crust-core transition in a magnetized crust. . 36

Table 5.1:

Urca pairs in the neutron star crust (Schatz et al., 2014). . . . . . . . . . . 90

Table 5.2:

Properties of the strongest Urca e− -capture parents identified in Meisel
& Deibel (2017) and in Schatz et al. (2014), absent even-A nuclides, excluded by Meisel et al. (2015); Deibel et al. (2016). . . . . . . . . . . . . 96

Table 5.3:

Urca pairs active in the neutron star ocean. . . . . . . . . . . . . . . . . . 102

ix

LIST OF FIGURES

Figure 2.1:

Schematic of a neutron star’s structure. Degenerate electrons provide the
primary pressure support in the outer crust and overlying layers. In the
inner crust, the dripped neutrons begin providing the primary pressure
support. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

Figure 2.2:

Ocean-crust transition density as a function of temperature for a neutron
star with g14 = 1.85 and using the crust equation of state from Haensel &
Zdunik (1990). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

Figure 2.3:

Equilibrium crust composition in the neutron star for the SLy4 crust equation of state. The panels are (a) the crust composition without shell corrections and (b) the crust composition with shell corrections. The shell
correction to the binding energy is shown in Equation (2.14). This figure
is reproduced from Deibel et al. (2014). . . . . . . . . . . . . . . . . . . 18

Figure 2.4:

Nuclear landscape in the inner crust at ρ = 6×1012 g cm−3 . Black squares
indicate locally stable nuclei that satisfy the stability condition in Equation 2.24. Nuclei are colored according to the binding energy per nucleon.

20

Nuclear landscape in the inner crust at ρ = 2×1013 g cm−3 . Black squares
indicate locally stable nuclei that satisfy the stability condition in Equation 2.24. Nuclei are colored according to the binding energy per nucleon.

21

Figure 2.5:

Figure 2.6:

Superfluid critical temperature models for the neutron star inner crust as
function of mass density. The three solid gray curves are from Gandolfi
et al. (2008). The solid black curve is from Schwenk et al. (2003). The
dashed black curve is from Gezerlis & Carlson (2011). The red curve is
from Chen et al. (1993). The blue curve is from Wambach et al. (1993). . 25

Figure 3.1:

Equilibrium composition of the neutron star crust in a strong magnetic
field B = 2×1015 G. (Left panel): Equilibrium crust composition without
a magnetic field. (Right panel): The black arrow indicates the first transition between equilibrium nuclei, which has moved to a greater depth as a
consequence of the strong magnetic field when compared to the left panel. 37

Figure 3.2:

Schematic of a torsional oscillation mode in the magnetar crust. Crust
oscillation frequencies are found using Equation 3.14, using the boundary
condition ξ = 0 (loss of traction) at the base and top of the crust. . . . . . 38

x

Figure 3.3:

Alfv´en velocity (blue curve) and shear velocity (black curve) in the crust
as a function of mass density. The composition is that of a 1.4 M neutron
star using the SLy4 crust EOS. . . . . . . . . . . . . . . . . . . . . . . . 40

Figure 3.4:

Frequency of the fundamental l = 2 mode as a function of the magnetar
mass for the core EOS probability distribution (centroid and ±2σ) from
Steiner et al. (2013) and an SLy4 crust EOS. The dashed black line indicates the observed 29 Hz QPO of SGR 1806-20. The frequencies are
evaluated for a crust-core transition density of 0.12 fm−3 with B = 0 G. . . 41

Figure 3.5:

Magnetar mass as a function of radius for the core EOS probability distribution from Steiner et al. (2013). Frequencies are evaluated using the
SLy4 crust EOS (L = 46 MeV) for nt = 0.12 fm−3 . The thick red solid
line indicates masses and radii for which the fundamental mode has a
frequency of 29 Hz in the case of SLy4 and 18 Hz in the case of Rs.
The black short-dashed line indicates masses and radii for a 626 Hz harmonic mode and B = 0 G. Masses and radii from 626 Hz harmonic modes
with magnetized crusts are labeled accordingly. Arrows indicate masses
and radii that match both the fundamental and the harmonic modes for
the field-free case and the case with the magnetic field of SGR 1806-20
(B = 2.4 × 1015 G). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

Figure 3.6:

Magnetar mass as a function of radius for the core EOS probability distribution from Steiner et al. (2013). Frequencies are evaluated using the
Rs crust EOS (L = 86 MeV) for nt = 0.12 fm−3 . The thick red solid
line indicates masses and radii for which the fundamental mode has a
frequency of 29 Hz in the case of SLy4 and 18 Hz in the case of Rs.
The black short-dashed line indicates masses and radii for a 626 Hz harmonic mode and B = 0 G. Masses and radii from 626 Hz harmonic modes
with magnetized crusts are labeled accordingly. Arrows indicate masses
and radii that match both the fundamental and the harmonic modes for
the field-free case and the case with the magnetic field of SGR 1806-20
(B = 2.4 × 1015 G). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

Figure 3.7:

Magnetar mass as a function of radius for EOS probability distribution
from Steiner et al. (2013). Frequencies are evaluated using the SLy4 crust
EOS with B = 0 G and nt = 0.12 fm−3 . The red dot-dashed, blue dotted,
and black dashed lines indicate masses and radii from fundamental modes
of frequency 29 Hz for different free neutron entrainment fractions fent .
The shaded band indicates masses and radii from 626 Hz harmonic modes
as fent is varied from 0.50 to 1.0. Arrows indicate the masses and radii
that match both the fundamental and the harmonic modes for fent = 1.0,
0.75, and 0.50. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

xi

Figure 3.8:

The same as Fig. 3.7, but for the Rs crust EOS with B = 0 G. Here the
free neutron entrainment fraction fent is varied from 0.20 to 0.30, with
fent labeled next to the corresponding curves. The red dot-dashed, blue
dotted, and black dashed lines indicate masses and radii from the 29 Hz
fundamental mode. The shaded band indicates masses and radii from
the 626 Hz harmonic mode. Arrows indicate the masses and radii for
fent = 0.30, 0.25, and 0.20 that match both the fundamental and harmonic
modes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

Figure 3.9:

The same as Fig. 3.5, but for the SLy4 crust EOS with nt = 0.08 fm−3 .
The thick red solid line indicates masses and radii determined from a
fundamental mode of 29 Hz. Masses and radii from harmonic modes
with magnetized crusts are labeled accordingly. Arrows indicate masses
and radii that match both the fundamental and the harmonic modes for
the field-free case and the case with the magnetic field of SGR 1806-20
(B = 2.4 × 1015 G). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

Figure 4.1:

Quiescent light curve of the neutron star MXB 1659-29. The model uses
a neutron star mass of M = 1.6 M , a neutron star radius of R = 11.2 km,
a shallow heat source of Qshallow = 1.0 MeV per accreted nucleon, and
a core temperature of T core = 4 × 107 K. The model also uses an impurity parameter of Qimp = 2.5 for the entire crust as done in Brown &
Cumming (2009). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

Figure 4.2:

Location of the ocean-crust transition as a function of temperature at the
base of the ocean. The black solid curve is for the one-component accreted composition (Haensel & Zdunik, 1990) and the black dotted curve
is for a multi-component accreted composition (Steiner, 2012). The crystallization density is given by Equation (2). . . . . . . . . . . . . . . . . . 57

Figure 4.3:

Quiescent lightcurve of the neutron star MAXI J0556-332. The solid
black curve corresponds to a model with M = 1.5 M , R = 11 km,
Qshallow = 6.0 MeV, and T core = 108 K; the dashed black curve is for
the same model with T core = 3 × 107 K. The black dotted curves are light
curves with a reheating event ≈ 170 days into quiescence for Qshallow =
6.0 MeV (upper curve) and Qshallow = 3.0 MeV (lower curve). The blue
dashed curve is for a M = 2.1 M , R = 12 km neutron star fit to the
observations by changing the shallow heating depth and strength. the data
above the light curve are contamination from residual accretion. Note that
∞ ∝ g1/4 /(1 + z) which leads to different observed core temperatures
T eff
for different gravities. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

Figure 4.4:

Schematic of the neutron star’s outer layers and the location of the anomalous shallow heating in MAXI J0556-332. . . . . . . . . . . . . . . . . . 61

xii

Figure 4.5:

Quiescent temperature evolution of the neutron star MAXI J0556-332.
Solid black curves indicate the evolution of the crust temperature during
quiescence for the M = 1.5 M and R = 11 km model, shown in Figure 4.3. The red dotted curve is the melting line of the crust (Γ = 175)
for the crust composition in Haensel & Zdunik (1990), the black dotted
curve is the transition from an electron-dominated heat capacity to an
ion-dominated heat capacity (CVe = CVion ), and the blue dotted curve is
where the local neutrino cooling time is equal to the thermal diffusion
time (τν = τtherm ). The gray dashed curve shows the lattice Debye temperature ΘD ; when T
ΘD electron-impurity scattering influences the
thermal conductivity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

Figure 4.6:

Markov chain Monte Carlo fits to the quiescent light curve of MAXI J0556332 using the crust relaxation code crustcool. The contours
√ show the
isodensity surfaces of the likelihood L, corresponding to −2lnL =
0.5, 1, 1.5, 2, for the neutron star mass M, radius R, pressure at the shallow heating depth Ph , and the shallow heating strength Qshallow . . . . . . 65

Figure 4.7:

Schematic of the neutron star’s outer layers and the location of g-modes. . 68

Figure 4.8:

Thermal evolution as a function of time during outburst/quiescence for
MAXI J0556-332 (solid curves), MXB 1659-29 (dashed curves), and
XTE J1701-462 (dot-dashed curves). Panel (a): Fundamental n = 1,
l = 2 g-mode frequencies in the ocean (Equation 4.10). Panel (b): The
ocean-crust transition density (Equation 2.5). Panel (c): Temperature at
the ocean-crust transition. . . . . . . . . . . . . . . . . . . . . . . . . . . 71

Figure 4.9:

Thermal evolution as a function of time during outburst/quiescence for
˙ = 0.1–1.0 M
˙ Edd . Panel (a): Fundamental
neutron star transients with M
n = 1, l = 2 g-mode frequencies in the ocean (Equation 4). Panel (b):
The ocean-crust transition density (Equation 2). Panel (c): Temperature
at the ocean-crust transition. . . . . . . . . . . . . . . . . . . . . . . . . 73

Figure 4.10:

Characteristic timescales in the neutron star ocean. The upper curves are
the thermal time for an equilibrium composition (solid curve) and accreted composition (dotted curve). The lower curves are for the timescale
for changes in ρt for an equilibrium composition (solid curve) and an
accreted composition (dotted curve). . . . . . . . . . . . . . . . . . . . . 77

xiii

Figure 4.11:

Crust cooling models of the late time cooling of MXB 1659-29. The solid
gray curve is a model that uses Qimp = 2.5 throughout the entire crust and
T core = 4×107 K. The solid blue curve is a model with Qimp = 20 for ρ >
8×1013 g cm−3 , Qimp = 1 for ρ < 8×1013 g cm−3 , T core = 3.25×107 K,
and using the G08 pairing gap. The dashed red curve uses the same Qimp
as the solid blue curve, but with the S03 pairing gap. The dotted blue
curve is a model with the G08 pairing gap and Qimp = 1 throughout the
crust, but without a low thermal conductivity pasta layer. . . . . . . . . . 79

Figure 4.12:

Thermal transport in the inner crust of MXB 1659-29 at the start of
quiescence. The gray vertical lines indicate the neutron drip density
and the transition to nuclear pasta. Panel (a): The temperature profile
(solid curve) corresponding to the cooling model in Figure 4.11 with a
Qimp = 20 pasta layer and the G08 pairing gap. The dashed curves show
two choices for T c (ρ); the blue dashed curve corresponds to G08 and the
red dotted curve is S03. Panel (b): The heat capacity profiles for the same
models as Figure 4.11. Solid black curve: Qimp = 3.7 throughout the inner crust. Dashed blue curve: Qimp = 20 for ρ > 8 × 1013 g cm−3 and
Qimp = 1 for ρ < 8 × 1013 g cm−3 using the G08 pairing gap that closes
in the crust. Dotted red curve: same as dashed curve, but with a different
choice for T c (ρ) from the S03 pairing gap that closes in the core. Panel
(c): Thermal conductivity profiles for the same models. . . . . . . . . . . 81

Figure 4.13:

Scattering frequencies for electrons and neutrons in the inner crust at
the beginning of quiescence for the model with Qimp = 20 at ρ > 8 ×
1013 g cm−3 , Qimp = 1 at ρ < 8 × 1013 g cm−3 , and the pairing gap that
closes in the crust (Gandolfi et al., 2008). Subplot: Thermal conductivity
K from electron scattering (dotted red curve), neutron scattering (dashed
blue curve), and from both electrons and neutrons (solid black curve).
The mass density ρ is given in units of 1014 g cm−3 . The region containing nuclear pasta is to the right of the vertical black dotted line. . . . . . . 83

Figure 5.1:

Schematic of an Urca reaction shell. The red colored region is composed
of the parent nucleus, which electron captures into the daughter nucleus
represented by the blue colored region. The Urca reaction layer has a
thermal width QEC ± kB T , where µe ≈ QEC is the center of the Urca shell. 88

Figure 5.2:

Luminosity sources in the neutron star crust as a function of crust temperature. The red curves indicate the neutrino luminosities from Urca pairs.
Blue curves indicate neutrino cooling from other processes in the neutron
star crust. The gray shaded region indicates the heating luminosity from
crustal heating. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91

xiv

Figure 5.3:

Quiescent lightcurve of the neutron star MAXI J0556-332. The crust
model is for a M = 1.5 M , R = 11 km neutron star with Qshallow =
6.0 MeV. The light curve without Urca cooling is shown as a red dashed
curve. Model light curves with Urca pairs, with Lν = 1036 ergs s−1 ,
are shown as black curves. From left to right, shell depths are y/yh =
3.4, 6.7, 17, 43, 140, that correspond to ρ10 = 3, 5, 10, 20, 50, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93

Figure 5.4:

Quiescent lightcurve of the neutron star MAXI J0556-332 with an 33 Al
Urca shell. The crust model is for a M = 1.5 M , R = 11 km neutron star
with Qshallow = 6.0 MeV. . . . . . . . . . . . . . . . . . . . . . . . . . . 94

Figure 5.5:

Schematic of the neutron star’s outer layers and the location of superburst
ignition and Urca cooling pairs. . . . . . . . . . . . . . . . . . . . . . . . 95

Figure 5.6:

Ocean temperature profiles during steady state with Urca pairs in the crust
for various values of Qb . Panel (a): Abundances of Urca pairs calculated
from X-ray burst ashes. Panel (b): Abundances of Urca pairs calculated
from superburst ashes. Panel (c): Ocean temperature profiles without
Urca cooling in the crust. . . . . . . . . . . . . . . . . . . . . . . . . . . 98

Figure 5.7:

Ocean temperature profiles during steady state with Urca pairs in the crust
for uncertainties in f t-values. We test the values Qb = 0.1, 0.25, 1, 2 MeV
per accreted nucleon. Panels (a)–(d): Abundances of crust Urca pairs
calculated from X-ray burst ashes. Panels (e)–(h): Abundances of crust
Urca pairs calculated from superburst ashes. The upper red-dashed curves
are with f t-values enhanced by a factor of 10 and the lower blue-dashed
curves are for f t-values reduced by a factor of 10. . . . . . . . . . . . . . 100

Figure 5.8:

Depth of Urca cooling pairs of a given mass number A. The size of data
points corresponds to the coefficient L34 given in Equation 5.8, and can
be found in Table 1. Urca pairs with L34 ≥ 10 are colored in red. The
gray band indicates empirical constraints on the superburst ignition depth
between yign ≈ 0.5–3 × 1012 g cm−2 (Cumming et al., 2006). . . . . . . . 102

Figure 5.9:

Temperature in the ocean as a function of column depth for a neutron
˙ Edd , M = 1.4 M , R = 10 km, and core
star model with: m
˙ = 0.3 m
temperature T core = 3 × 107 K. The dashed-black curves indicate the
ignition of unstable carbon burning in a mixed iron-carbon ocean with
XFe = 0.8 and XC = 0.2. The solid-black curves indicate ocean’s without
Urca cooling layers. In each panel, the lower dotted-blue curve is for
L34 = 100 and the upper dotted-blue curve is for L34 = 10, for Urca
pairs located at: Panel (a): y = 1.2 × 1013 g cm−2 , Panel (b): y = 1.2 ×
1012 g cm−2 , and Panel (c): y = 1.2 × 1010 g cm−2 . . . . . . . . . . . . . 105

xv

Figure 5.10:

Cooling light curves for superbursts with E18 = 0.2, α = 0.25, g14 = 1.6,
and ignition depths yign = 3.2×1011 g cm−2 and yign = 3.2×1012 g cm−2 .
The red dotted curves are models without Urca cooling. The solid curves
are for a cooling source at y = 0.4 yign with X · L34 = 100, 1000. . . . . . 106

Figure 5.11:

Steady state ocean temperature as a function of column depth for a neu˙ Edd , M = 1.4 M , R = 10 km, and core
tron star model with: m
˙ = 0.3 m
7
temperature T core = 3 × 10 K. Solid black curves are steady state temperature profiles for the heating strengths Qshallow = 1, 2, 4. The dashed
black curve indicates the ignition of unstable carbon burning in a mixed
iron-carbon ocean with XFe = 0.8 and XC = 0.2. . . . . . . . . . . . . . 108

xvi

Chapter 1
Introduction
The theoretical prediction (Rutherford, 1920) and discovery (Chadwick, 1932) of a neutral atomic
particle, the neutron, had far reaching consequences for astronomy. It was soon proposed that
supernovae mark the transition of ordinary stars into compact stars made primarily of neutrons
(Baade & Zwicky, 1934). The properties of the neutron in concert with the equations of stellar
structure, form a theoretical foundation for neutron stars (Tolman, 1934; Oppenheimer & Volkoff,
1939). On the cusp of becoming black holes, neutron stars reach densities greater than an atomic
nucleus in their cores and epitomize the fascinating nature of matter under the most extreme conditions. Decades after their conception, neutron stars were discovered at X-ray wavelengths (Giacconi et al., 1962) and radio wavelengths (Hewish et al., 1968), and with the advent of modern
space-based telescopes (e.g., Chandra X-ray Observatory) it was discovered that neutron stars
produce a diverse array of high-energy phenomena.
Of particular interest to the study of the neutron star interior is the X-ray emission from lowmass X-ray binaries, where a neutron star primary orbits with a low-mass companion. An accretion
disk forms around the neutron star as a result of mass transfer from the low-mass companion. Neutron stars with accretion disks undergo sporadic accretion outburst episodes that may last for weeks
to years; for example, accretion outbursts are observed on the accretion neutron star KS 1731-260
(Sunyaev & Kwant Team, 1989; Wijnands et al., 2001). The accumulating material during an accretion outburst triggers non-equilibrium reactions in the neutron star’s outer layers (Sato, 1979;
Bisnovaty˘i-Kogan & Chechetkin, 1979) which may heat the crust out of thermal equilibrium with
1

the core (Brown & Bildsten, 1998; Brown, 2000). When the accretion outburst ends, stored heat
from accretion-driven reactions diffuses toward the neutron star’s surface where it powers the quiescent light curve (Eichler & Cheng, 1989; Ushomirsky & Rutledge, 2001; Rutledge et al., 2002).
Reproducing the observed quiescent light curve with numerical models reveals new clues about
the physics of the interior.
For instance, the cooling behavior of neutron star transients reveals the thermal properties of
nuclear pasta. Deep in the neutron star’s crust, nuclei are deformed into various complex shapes in
the high density environment (Ravenhall et al., 1983; Hashimoto et al., 1984). Nuclear pasta has
been studied extensively using molecular dynamics simulations (e.g., Horowitz et al. 2004a), but
the thermal and electrical conductivities remain difficult to determine (Horowitz et al., 2015). The
late time cooling of quiescent neutron star transients, however, is sensitive to the thermal transport
properties of the deep inner crust. As a result, late time cooling light curves may reveal the thermal
conductivity of pasta (Horowitz et al., 2015; Deibel et al., 2017). Furthermore, the late time cooling
is also sensitive to the thermal properties of normal and superfluid neutrons that coexist with the
pasta at high densities.
Models of cooling neutron stars also pose new questions. For example, the large quiescent luminosities of some neutron star transients, such as the hottest transient MAXI J0556-332 (Homan
et al., 2014), can not be explained by accretion-driven heating alone (Deibel et al., 2015). Interestingly, reconciling neutron star superburst observations with superburst ignition models also
requires extra heating in the neutron star’s outer layers (Cumming et al., 2006). The source of extra heating is unknown, but may be connected to the deposition of rotational energy in the neutron
star ocean by a spreading boundary layer of accreted material (Inogamov & Sunyaev, 1999, 2010;
Philippov et al., 2016).
Thermonuclear bursts from accretion neutron stars also provide an avenue to examine the
2

physics of the neutron star interior. For example, X-ray bursts observed from neutron star transients
are well fit by unstable ignition of hydrogen and helium in the neutron stars envelope (Woosley
et al., 2004). Long type-I X-ray bursts, or superbursts, are even more energetic explosions in the
neutron star envelope (Woosley & Taam, 1976; Taam & Picklum, 1978). Superbursts are thought
to be triggered by the unstable ignition of the 12 C + 12 C fusion reaction (Cumming & Bildsten,
2001; Strohmayer & Brown, 2002a). Much like the crust cooling models of neutron star transients,
superburst ignition models require anomalous extra heating to ignite carbon at the depths inferred
from superburst cooling light curves (Cumming et al., 2006). The source of extra heating in these
models is unknown, but is required to match observed ignition depths for all superbursts.
The physics of the neutron star interior may be altered by bygone nuclear burning in the neutron
star envelope. Nuclear burning during X-ray bursts and superbursts produces an array of protonrich nuclei (Wallace & Woosley, 1981; Woosley et al., 2004; Keek & Heger, 2011). The ashes of
nuclear burning are compressed to greater depths in the neutron star’s ocean and crust by subsequent accretion outbursts. Not only does compression deposit heat via nuclear reactions (Haensel
& Zdunik, 1990, 2008), but advances in accretion-driven nuclear reaction studies (Schatz et al.,
2014) have revealed that neutrino cooling may occur in nuclear burning ashes. Some nuclei form
e− −capture / β− −decay cycles, known as Urca cycles, that emit neutrinos and effectively cool
the neutron star’s outer layers (Gamow & Schoenberg, 1941; Bahcall & Wolf, 1965). Neutrino
cooling will in turn impact thermal evolution models of neutron star transients (Meisel & Deibel,
2017) and carbon ignition depths in superburst ignition models (Deibel et al., 2016).
Isolated and highly-magnetized neutron stars, or magnetars, produce energetic γ-ray flares
(Barat et al., 1983) which may offer a glimpse of the neutron star’s interior. These giant flares
may be triggered by a reconfiguration of the magnetic field which fractures the crystalline crust
and induces oscillations in the crust, or starquakes (Duncan, 1998). Quasi-periodic oscillations
3

in the γ-ray flare emission (Israel et al., 2005; Strohmayer & Watts, 2005; Watts & Strohmayer,
2006) may be modes of the oscillating crust. Because crust oscillation modes are uniquely determined by neutron star properties, a study of these oscillations may lead to constraints on the
interior physics (Piro, 2005; Steiner & Watts, 2009; Deibel et al., 2014). Although quasi-periodic
oscillations originate from accretion disks around accreting black holes and neutron stars (e.g.,
Lewin & van der Klis 2006), surface oscillations are a compelling explanation for the observed
quasi-periodic oscillations in magnetars; magnetars are isolated and do not have an accretion disk
due to mass transfer from a companion. Furthermore, matching observed oscillation frequencies
to crust oscillation frequencies can produce constraints relevant for nuclear theory because crust
oscillation frequencies are sensitive to the highest density regions of the neutron star crust.
We begin in Chapter 2 by constructing a model of the neutron star’s outer layers. This model
draws from nuclear theory and uses energetic arguments to calculate a ground-state composition
for the neutron star’s outer layers. A model of the equilibrium composition is a necessary starting
point when studying the dynamical, compositional, and thermal evolution of the neutron star’s
outer layers. We also discuss a calculation for the accreted composition by following nuclear
reactions in accreted material using energetic arguments. The thermal transport properties of the
interior, such as the specific heat and thermal conductivity of various layers, are also discussed.
In Chapter 3, the crust model is applied to calculating crust oscillation frequencies for starquaking magnetars. The large magnetic fields of magnetars

1014 G alter the crust composition

by quantizing electrons. Furthermore, magnetic driven Alfv´en waves in the crust must be incorporated into solutions for crust oscillations. We match crust oscillation frequencies to the observed
giant flare quasi-periodic oscillations in the magnetars SGR 1806-20 and SGR 1900+14 to draw
constraints on the structure of these stars and nuclear physics at high mass densities.
Moving away from isolated neutron stars, we examine the observable phenomena of neutron
4

stars in accreting systems. Chapter 4 demonstrates that not only can crust models reproduce observations, but the can probe the neutron star interior in an intuitive manner. By implementing the
thermal transport properties of dense matter into the crust model, we follow the thermal evolution
of cooling neutron star transients (Chapter 4). Such a calculation is applicable to low-mass X-ray
binaries where a neutron star primary orbits with a low-mass companion. During an accretion
outburst onto the neutron star primary, the continual compression of the neutron star’s outer layers induces non-equilibrium reactions in the crust; reactions that heat the neutron star crust out of
thermal equilibrium with the core. We model the thermal relaxation of the crust when accretion
halts and accurately reproduce quiescent light curves of neutron star transients. These fits reveal,
for example, that powerful heating during outburst is required to fit the quiescent observations of
the neutron star MAXI J0556-332 and the source of this extra heating remains unknown.
Being confident in our knowledge of thermal transport in the crust, we outline an observational
test to answer an open question in neutron star transient observations. Observed quasi-periodic
oscillations in the transient neutron star Z-sources, such as Sco X-1, are of unknown origin. Curiously, some observed frequencies are close to the neutron star ocean’s thermal g-mode frequencies.
Reproducing thermal g-modes in our model, we investigate if observed oscillations originate in the
ocean. In particular, we show that as the ocean heats up during accretion, thermal g-mode frequencies increase, and we predict that observed oscillation frequencies should be larger in hotter neutron
star transients.
We then incorporate nuclear reactions into the crust model and demonstrate that light curve
fits to neutron star transients constrain nuclear processes (Chapter 5). Specifically, we use the
quiescent light curve of MAXI J0556-332 to determine if Urca cooling reaction layers, rapid backand-forth cycles of e− -capture and β− -decay, are active in the accreting neutron star crust. We
demonstrate that a crust thermal relaxation model with Urca cooling layers does not produce a
5

light curve consistent with the quiescent observations of MAXI J0556-332. This work motivates
subgrid models of Urca cooling reactions to further resolve the reaction’s properties and delineate
its impact on X-ray bursts, superbursts, and crust thermal relaxation.
Lastly, we discuss the prospects of future research. A high priority is determining the origin of
the anomalous extra heating in neutron star transients, such as the hottest transient MAXI J0556332. Extra heating, from an unknown origin, is also required in superburst ignition models. Preliminary studies suggest that the extra heating may be from friction at the base of the neutron star
ocean, as the ocean is spun up to higher rotation velocities during accretion.

6

Chapter 2
The Neutron Star Ocean and Crust
In this chapter we outline the physics input for numerical models of a neutron star’s outer layers.
In particular, the numerical model includes a calculation of the composition from the low-density
equation of state and the thermal transport properties of the neutron star’s outer layers. The equilibrium composition will then be used in the context of highly-magnetized neutron stars in Chapter 3.
In Chapter 4, the thermal transport properties of the outer layers are used to investigate the thermal
evolution of neutron star transients.

2.1

Structure

An accreted fluid element will first enter the neutron star atmosphere which is rarified hydrogen.
A schematic of neutron star structure can be seen in Figure 2.1. Nuclei heavier than hydrogen
gravitationally settle to the bottom of the atmosphere and form a boundary layer between the
atmosphere and the ocean beneath. The ocean is a plasma of ions and electrons, and contains
nuclear burning shells that alter the ocean’s composition. For example, hydrogen-rich and heliumrich material burn to iron-group elements (Taam et al., 1996). The electrons are partially degenerate
and relativistic, and provide the majority of the pressure support. Moving deeper into the ocean,
the electron gas becomes degenerate and highly relativistic near ρ

7

106 g cm−3 and the pressure

is (see Chapter 2, Shapiro & Teukolsky 1983)

P = Pe ≈

p4F,e c

,
12π2 3

(2.1)

where pF,e = (3π2 ne )1/3 is the electron Fermi momentum and ne is the number density of electrons. When scaled to typical values of the ocean the pressure is

P ≈ 3.6 × 1026 ergs cm−3 ρ4/3
(Ye /0.4)4/3 ,
9

(2.2)

where Ye is the electron fraction and ρ9 ≡ ρ/(109 g cm−3 ). Now we can define column depth
(y ≈ P/g) as
y ≈ 3.6 × 1012 g cm−2 ρ4/3
(Ye /0.4)4/3 g−1
14 ,
9

(2.3)

where the surface gravity is g14 = g/(1014 cm s−2 ).
A phase transition occurs at the bottom of the ocean where the ion density becomes high, and
the ocean nuclei crystallize into a solid crust (Baym et al., 1971b; Horowitz & Berry, 2009) due to
Coulomb repulsion of the constituent ions (Wigner, 1934). The point at which the ion-liquid ocean
crystallizes into an ion-lattice is determined by the plasma coupling parameter of a composition
with average proton number Z ,
Γ=

Z 2 e2
,
akB T

(2.4)

where a = (3 Z /4πne )1/3 is the radius of the Wigner-Seitz cell (i.e., the inter-ionic spacing).
The ocean’s crystallization point (or the crust’s melting point) occurs when Γ = 175 (Farouki
& Hamaguchi, 1993; Potekhin & Chabrier, 2000; Horowitz et al., 2010). Therefore, rearranging

8

10-9 envelope
103 ocean (ion liquid)

e

100

outer crust (ion lattice solid)

400

1011 inner crust (ion lattice + neutrons)

µ

1012
1013

depth [m]

mass density [g cm-3]

106

1

µ

n

nuclear pasta

1000

1014 core

Figure 2.1: Schematic of a neutron star’s structure. Degenerate electrons provide the primary
pressure support in the outer crust and overlying layers. In the inner crust, the dripped neutrons
begin providing the primary pressure support.
Equation 2.4, the ocean-crust transition density can be expressed as

ρt ≈

2.2 × 106 g

cm−3

3 26 6 A
T
Z
56
3 × 107 K

.

(2.5)

where A is the average mass number of the composition. The ocean-crust transition density as
a function of crust temperature can be seen in Figure 2.2. At mass densities ρ > ρt the liquid
ion-ocean crystallizes into a body-centered cubic nuclear lattice known as the outer crust. The
pressure in the outer crust is primarily supported by the degenerate and relativistic electron gas
with a small pressure contribution from the nuclear lattice (Baym et al., 1971a,b).
The crust is composed of a nuclear crystal lattice embedded in an degenerate electron gas
at lower densities (ρt
higher densities ρ

ρ

ρdrip ≈ 4 × 1011 g cm−3 ). A quasi-free neutron gas appears at

ρdrip ; The neutron-drip point, the point at which it becomes energetically

favorable for neutrons to drip out of nuclei, marks the beginning of the inner crust. The inner
9

400
300

1010

200

[ ]

]

1011

[

109
108
107
106 7
10

120
108

[ ]

109

Figure 2.2: Ocean-crust transition density as a function of temperature for a neutron star with
g14 = 1.85 and using the crust equation of state from Haensel & Zdunik (1990).
crust can then be described as a nuclear lattice embedded in both a degenerate electron gas and
a quasi-free neutron gas (Baym et al., 1971a), so the total pressure has three contributions, from
the electrons, lattice, and quasi-free neutrons. The neutron pressure becomes the primary pressure
support at mass densities near ρ

4 × 1012 g cm−3 in the inner crust (Shapiro & Teukolsky,

1983), though the electron gas still has a non-negligible pressure contribution. In the deep inner
crust at mass densities near ρ

6 × 1013 g cm−3 the nuclei deform in an attempt to minimize

their energy at such extreme pressures. These deformed nuclei are called nuclear pasta (Ravenhall
et al., 1983; Pethick & Ravenhall, 1995; Pethick & Potekhin, 1998) and have been studied through
semi-classical molecular dynamics simulations (Horowitz et al., 2004b; Horowitz & Berry, 2008;
Schneider et al., 2013). The thermal transport properties of the nuclear pasta phase, however,
remains unknown (Horowitz et al., 2015). We will discuss the possibility of constraining the
thermal conductivity of pasta using late time observations of cooling neutron stars in Section 4.4.
The crust dissolves into a nucleon liquid near a mass density ρ ≈ 2.8×1014 g cm−3 (Baym et al.,
1971a) via a first order phase transition which marks the beginning of the outer core. The crust-

10

core transition occurs slightly above the nuclear saturation density ρ0 ≈ 2 × 1014 g cm−3 (Baym
et al., 1971a; Pethick et al., 1995) — the density of an atomic nucleus. The outer core contains
µ, e, n, p matter because the electron chemical potential µe approaches the rest-mass energy of the
muon (mµ c2 ≈ 105 MeV). The pressure support comes from the repulsive core of the nuclear force
because nucleons are packed to densities ρ

2.2

ρ0 .

Low-density Equation of State

The crust composition at a given mass density is calculated by examining the total energy density
of the system. The crust is imagined as a series of “drops” of nuclear matter, each inside an
identical Wigner-Seitz cell. The energetics of a single cell, how a nucleus interacts with the electron
and neutron gases sharing the cell, is assumed to represent the entire crust at the mass density
considered. The nucleus occupies a volume fraction χ of the cell and the average baryon density
in the nucleus is nl = nn + n p , where nn is the density of neutrons and n p the density of protons
inside the nucleus. The number density of baryons per unit volume is n = χ(nn + n p ) + (1 − χ)ndrip ,
where ndrip represents the density of quasi-free dripped neutrons. The total energy density of the
crust is

E
(Z, A)
+ (1 − χ) (nn = ndrip , n p = 0) + we (ne ) ,
w(Z, A, n) = χ nn mn + n p m p + nl bind
A

(2.6)

where we and ne are the energy density and number density of electrons, respectively. Note that
we = 3Pe where Pe is the degenerate electron pressure given in Equation 2.1. We use a liquid drop
mass model (Baym et al., 1971a,b; Ravenhall et al., 1983; Steiner, 2008) for the nucleus binding
energy Ebind that contains the energy density contributions from the nucleus and the lattice: the

11

Coulomb energy (wCoul. ), the surface energy (Esurf ), the shell energy (Eshell ), and pairing energy
(Epair ) corrections to the homogeneous bulk matter Hamiltonian ( ) of bulk n, p, e matter.

2.2.1

Nuclei and the Ion-lattice

In both the inner and outer crust the nuclear contribution to the energy density will be that of the
equilibrium nucleus at the given baryon density. The lattice contribution will be that of a crystal
lattice composed of equilibrium nuclei. For these contributions to the energy density we use a
liquid-drop mass model which includes the lattice contribution (Baym et al., 1971a,b; Ravenhall
et al., 1983; Steiner, 2008). We also add shell corrections to the mass model as described in
Dieperink & van Isacker (2009) and an updated neutron-drip line to better predict the neutron
rich nuclei beyond the neutron-drip point. The first term on the right-hand side of Equation 2.6
represents the energy density of the nucleus and the lattice;

E
(Z, A)
n
= (nn , n p ) + l Esurf + Eshell + Epair + wCoul.
nn mn + n p m p + nl bind
A
A

(2.7)

where nn , n p , and nl are respectively the average neutron, proton, and baryon number densities
inside the nucleus with a given proton number Z and baryon number A. The parameters mn and
m p are the neutron and proton masses. The above parameters will be determined from

nl = nn + n p = n0 + n1 I 2 ,

(2.8)

where I = 1 − 2Z/A is the isospin asymmetry. The quantity n0 is the nuclear saturation density
and n1 is a negative quantity that allows for the decrease in saturation density based on the isospin
asymmetry and the increase of saturation density due to the Coulomb interaction (Steiner, 2008).

12

The average neutron and proton densities within the nucleus are

nn = nl (1 + δ)/2 ,

and
n p = nl (1 − δ)/2 ,

(2.9)

where δ = 1 − 2n p /(nn + n p ) is the density asymmetry and η = δ/I = 0.92 is a constant of our
model that will determine the thickness of a neutron skin (i.e., the difference between neutron and
proton radii; Steiner 2008). The bulk matter binding energy is

Ebulk =

A
[ (nn , n p ) − nn mn − n p m p ] ,
nl

(2.10)

where (nn , n p ) is the energy density of homogeneous bulk matter at a given neutron and proton
number density, for which we use the bulk matter Hamiltonian in the Skyrme model (Skyrme,
1959) with SLy4 coefficients (Chabanat et al., 1995). The surface energy goes as

Esurf = σ (nn , n p )

where the unitless function

36π 1/3
,
nl 2 A

(2.11)

that is proportional to the surface tension σ is

(nn , n p ) = 1 − σδ δ2 ,

(2.12)

where σδ represents the surface energy density (Myers & Swiatecki, 1969; Steiner et al., 2005).

13

The Coulomb energy density is

wCoul. =

2π 2 2 2
n p e R p 2 − 3χ1/3 + χ ,
5

(2.13)

where e2 is the Coulomb coupling and R p is the proton radius (3Z = 4πn p R p 3 ). The respective χ
terms in parenthesis correspond to the Coulomb contribution, the lattice contribution, and a high
density correction as χ approaches unity. The shell corrections to the binding energy per baryon
(Sato, 1979; Dieperink & van Isacker, 2009) are

Eshell (Z, N) = a1 S 2 + a2 (S 2 )2 + a3 S 3 + anp S np ,

(2.14)

where
nv n¯v zv z¯v
+
,
Dn
Dz

(2.15)

nv n¯v (nv − n¯v ) zv z¯v (zv − z¯v )
+
,
Dn
Dz

(2.16)

nv n¯v zv z¯v
,
Dn Dz

(2.17)

S2 =
S3 =

S np =
and

n¯v ≡ Dn − nv ,
z¯v ≡ Dz − zv .

(2.18)

The parameters Dn and Dz correspond to the degeneracy of the neutron and proton shells, i.e.,
the difference between the magic numbers enclosing the current amount of neutrons or protons.
The quantities nv and zv are the number of valence neutrons and protons, i.e., the difference between the current number of protons or neutrons and the preceding magic number. The effect of the

14

Table 2.1: Parameters of the mass model.
Parameter
SLy4
n0
0.1740 fm−3
n2
−0.0157 fm−3
η
0.9208
σδ
1.964
σ
1.164 MeV
a1
−1.217 MeV
a2
0.0256 MeV
a3
0.00387 MeV
anp
0.0357 MeV
ap
5.277 MeV

Rs
0.1597 fm−3
0.0244 fm−3
0.9043
1.465
1.041 MeV
−1.298 MeV
0.0311 MeV
0.00349 MeV
0.0287 MeV
5.265 MeV

shell corrections on the composition of the crust can be seen in Figure 2.3. The pairing contribution
to the energy density taken from Brehm (1989) with updated coefficients is





−a p A−1/3 , even-even,







Epair = 
+a p A−1/3 ,
odd-odd,










0,
even-odd,

(2.19)

where a p is a constant of our model. The combined contributions of nuclear matter and the nuclei
crystal lattice result in the expression,

wnuc + wl = (nn + n p )χB(Z, A)/A .

(2.20)

The coefficients used in the mass model above are listed in Table 2.1.

2.2.2

Quasi-free Neutrons

The energy density of dripped neutrons is that of bulk neutron matter,

n,drip

= (1 − χ) (nn,drip , n p = 0) ,
15

(2.21)

where we use the SLy4 bulk matter energy density (Skyrme, 1959) with SLy4 coefficients (Chabanat et al., 1995). The number density of dripped neutrons is determined from (Baym et al.,
1971a)
n = χ(nn + n p ) + nn,drip (1 − χ) ,

(2.22)

where n is the number density of baryons. The approach accounts for the filling of space with
dripped neutrons, as opposed to nucleic matter, at depths greater than the neutron-drip point ρ
ρdrip ≈ 4 × 1011 g cm−3 . As the density approaches nuclear saturation the fraction of space filled
by the neutron gas approaches unity, that is, the filling fraction of nuclei χ goes to zero.

2.3

Composition

The crust equation of state discussed in Section 2.2 can be used to calculate the composition of
the neutron star crust. In a cold-catalyzed crust, a crust made of nuclei frozen in from the neutron
star’s birth, the composition is in its most energetically favorable configuration. A cold-catalyzed
composition is typically adequate when modeling neutron star observables. Accreting neutron
stars, however, have crusts polluted by compressed accreted material and thermonuclear burning
ashes. A calculation of the accreted neutron star crust must account for the nuclear reactions that
take place in an accreted fluid element as it is compressed through the electron gas of the outer
crust and the quasi-free neutron gas of the inner crust.

2.3.1

Equilibrium

The ocean and crust of an isolated neutron star are composed of cold-catalyzed matter remaining
from when the neutron star was first formed. During the hot post-supernova conditions all nuclear

16

reactions are assumed possible and the composition will find the most energetically favorable configuration. The most energetically favorable configuration is found by minimizing the total energy
density (Equation 2.6) at a given mass density. The nucleus that minimizes the total energy density
is called the equilibrium nucleus.
We now use the formulae outlined in the previous section to calculate the equilibrium composition of the ocean and crust of an isolated neutron star. The equilibrium composition for the
SLy4 crust equation of state is shown in Figure 2.3. In the neutron star ocean, at mass densities
ρ < 106 g cm−3 , the equilibrium nucleus is 56 Fe as it has the minimum binding energy per baryon
in our mass model. At mass densities ρ

106 g cm−3 electrons are degnerate and relativistic, and

the electron chemical potential µe increases monotonically with mass density. As a consequence,
the equilibrium nucleus must be increasingly neutron-rich, otherwise it would energetically favorable for the nucleus to undergo e− −capture .

2.3.2

Accreted

On accreting neutron stars, the accumulation of accreted material in the envelope compresses
deeper layers and replaces the existing composition. In fact, the entirety of the ocean and crust
are replaced by accreted material over a timescale of ∼ 106 years. Although an accretion outburst
initially deposits hydrogen-rich and helium-rich material into the envelope, the accreted material eventually burns to heavier elements via the 3α, αp, and rp-processes (Wallace & Woosley,
1981). For example, unstable burning during type-I X-ray bursts produces range of nuclei with
A ∼ 60–100 (Schatz et al., 1998, 2001; Brown et al., 2002; Schatz et al., 2003; Woosley et al.,
2004; Cyburt et al., 2010), which are pushed deeper into the neutron star by subsequent accretion outbursts. A neutron star crust composed of cold-catalyzed matter, that is, a crust of entirely
equilibrium nuclei, is ultimately replaced by ashes compressed by accretion.
17

Z

Z

50

N

0

N

proton, neutron number

200
150
100

200
50

100

150

with shell
corrections

no shell
corrections

0

proton, neutron number

(b)

(a)

6.0

8.0

10.0 12.0 14.0
−3

log10 (ρ) [g cm

]

6.0

8.0

10.0 12.0 14.0
−3

log10 (ρ) [g cm

]

Figure 2.3: Equilibrium crust composition in the neutron star for the SLy4 crust equation of state.
The panels are (a) the crust composition without shell corrections and (b) the crust composition
with shell corrections. The shell correction to the binding energy is shown in Equation (2.14). This
figure is reproduced from Deibel et al. (2014).
In addition to replacing the pristine crust, accreted material undergoes non-equilibrium processes that deposit heat (Bisnovaty˘i-Kogan & Chechetkin, 1979; Sato, 1979; Haensel & Zdunik,
1990). More recent studies of non-equilibrium processes, the excited states of nuclei, and pycnonuclear reactions, have refined our understanding of the amount of heat deposited in the crust
during accretion as well as the final composition of the accreted crust (Haensel & Zdunik, 2003,
2008; Gupta et al., 2007; Gupta et al., 2008; Steiner, 2012; Lau, 2012). The accreted crust composition is thought to have a low impurity parameter (Equation 2.32), a result that is further supported
by observations of neutron star crust cooling in quasi-persistent transients (Brown & Cumming,
2009).
The energetics of the electron and neutron gases in the crust determine the allowed reactions
in an accreted element. The accreted material at a given depth in the crust must be composed of

18

nuclei that stably co-exist with the electron and neutron gas at that depth. The stability criterion
for nuclei in the crust is outlined in Appendix 6.
An accreting neutron star will have a crust composition set by the composition of initially
accreted material and how that material is altered by the subsequent non-equilibrium reactions
(Bisnovaty˘i-Kogan & Chechetkin, 1979; Sato, 1979) that take place as it is compressed deeper in
the star by ongoing accretion. An example of a non-equilibrium reaction is an electron capture,
which occurs in the crust when

µe − B(Z, A) > ∆n − B(Z − 1, A)

(2.23)

where B(Z, A) is the binding energy of the nucleus and ∆n ≡ mn − m p − me is the binding energy
of the neutron. As µe increases with depth in the outer crust, e− -capture reactions occur that
push the composition to neutron-rich isotopes. In addition to changing the crust composition,
non-equilibrium reactions deposit ≈ 1–2 MeV per accreted nucleon of heat into crust (Haensel &
Zdunik, 1990, 2003; Gupta et al., 2007; Gupta et al., 2008; Haensel & Zdunik, 2008; Lau, 2012).
The condition for local stability against non-equilibrium reactions is (derived in Appendix 6)
∂2 B ∂2 B
+
<0,
∂Z 2 ∂A2

(2.24)

which holds for all combinations of e− -capture, β− -decay, n-capture, and n-emission. The chemical potential of electrons µe and the chemical potential of neutrons µn determine the non-equilibrium
nuclear reactions that will take place in accreted material during compression.
In the outer crust electron captures make the composition more neutron-rich until the neutrondrip point is reached. In the inner crust neutron emissions make nuclei lighter and contribute neu-

19

Figure 2.4: Nuclear landscape in the inner crust at ρ = 6 × 1012 g cm−3 . Black squares indicate locally stable nuclei that satisfy the stability condition in Equation 2.24. Nuclei are colored
according to the binding energy per nucleon.
trons to the quasi-free neutron gas until the crust-core transition is reached. Using only energetic
arguments, the nuclei in the accreted element are forced to “land” on stable nuclei at each depth.
In other words, the stable nuclei that make up the composition at a given density are nuclei which
can stably co-exist with the electron and neutron gases at their respective chemical potentials. The
accreted crust is composed of stable isotopes that may or may not be equilibrium nuclei.

2.4

Thermal Transport

For studies of the thermal evolution of accreting neutron stars (discussed in Chapter 4) it is necessary to formulate expressions for the thermal transport properties of the neutron star’s outer layers.
Heat diffuses throughout the entire crust through electron conduction and the evolution of temperature with time follows the thermal diffusion equation
∂T
1
∂
∂T
nuc − ν
=
−
−4πr2 K
2
∂t
CV
∂r
4πr ρCV ∂r
20

,

(2.25)

Figure 2.5: Nuclear landscape in the inner crust at ρ = 2 × 1013 g cm−3 . Black squares indicate locally stable nuclei that satisfy the stability condition in Equation 2.24. Nuclei are colored
according to the binding energy per nucleon.
where CV is the specific heat, ν is the specific neutrino emissivity, nuc is the specific nuclear heating emissivity, K is the thermal conductivity and we have ignored general relativistic corrections
over a thin shell. The derivation of the expression in Equation 2.25 can be found in Appendix 6.

2.4.1

Outer Crust

In the outer crust, the degenerate electron gas mediates the thermal transport. The thermal conductivity of the electrons is
Ke =
where ν =

j 1/τ j

2 2
π2 ne c kB T
,
3 EF,e ν

(2.26)

is the total electron scattering frequency from j sources, τ is the relaxation

time of electron scattering off of scatterer j, and EF,e is the electron Fermi energy. The dominant contribution to the electron scattering frequency depends on the electron temperature relative to the plasma temperature, T P = ( /kB )ωP and the electron plasma frequency is ωP =
( /kB )(4πZ 2 e2 n/Amu )1/2 .
21

In the ocean, the contributions to the total electron scattering frequency are 1/τ = 1/τee +1/τei ,
where νee = 1/τee is the electron-electron scattering frequency and νei = 1/τei is the electron-ion
scattering frequency. Because the electrons are degenerate at densities of interest

106 g cm−3 in

the ocean, only electron-ion scattering is important. If the crust temperature is T > ΘD , where the
Debye temperature is ΘD ≈ 0.45T P , the specific heat is set by the electrons

CVe = π2

2T
ZkB

Am p EF,e

,

(2.27)

and the thermal conductivity is set by electron-ion scattering (Yakovlev & Urpin, 1980; Potekhin
et al., 1997, 1999),
4EF,e Zα2 Λei
νei =
,
3π

(2.28)

where α = e2 / c is the fine structure constant and Λei is the Coulomb logarithm for electron-ion
scattering. This treatment of electron-ion scattering is a good approximation for one-component
compositions, but accreted crusts with multi-component compositions have a lower thermal conductivity than found using this method (Roggero & Reddy, 2016). At mass densities above the
ocean-crust transition (Equation 2.5), the specific heat is set by the ions

CVion ≈

3kB
Amu

(2.29)

and the thermal conductivity is also set by electron-ion scattering. When the crust temperature
is T

ΘD , however, the thermal conductivity is set by electron-phonon scattering (Yakovlev &

Urpin, 1980; Baiko & Yakovlev, 1995, 1996),

k T
νep ≈ 13α B ,

22

(2.30)

where α = e2 / c ≈ 1/137 is the fine structure constant. In this same temperature regime (T

ΘD ),

if impurities exist in the crust then electron-impurity scattering becomes important, with scattering
frequency (Flowers & Itoh, 1976; Brown & Cumming, 2009),

νeQ =

4πQimp e4 nnuc Λimp
p2F,e vF,e

,

(2.31)

where nnuc is the number density of nuclei, pF,e is the electron Fermi momentum, and vF,e is the
electron Fermi velocity. The impurity parameter is given by

Qimp ≡

1
n

ni (Zi − Z )2 ,

(2.32)

i

where ni is the number density of the nuclear species with Zi number of protons and Z is the
average proton number of the composition. A value of Qimp = 0 implies a pure crust and Qimp > 0
contains impurities. The impurity parameter has been constrained from crust cooling models of
quiescent neutron star transients (Brown & Cumming, 2009). For example, Brown & Cumming
(2009) find that the quiescent light curve in MXB 1659-29 is best fit with Qimp ≈ 4 and the light
curve of KS 1731-260 is best fit with Qimp ≈ 1.5. The similar Qimp parameters in these sources,
despite their different accretion histories, would suggest that accreted material is relatively free of
impurities once compressed into the inner crust (as suggested by calculations of nuclear reactions
in accreted material done in Section 2.3.2). By contrast, more recently accreted material in the
outer crust that has yet to be processed completely by non-equilibrium nuclear reactions has a large
impurity parameter. Models of superburst ignition (to be discussed in more detail in Section 5.2.2)
need Qimp ≈ 100 to have a sufficiently large heat flux into the ocean from crustal heating to ignite
superbursts (Cumming et al., 2006).

23

2.4.2

Inner Crust

In the inner crust the electron thermal conductivity is primarily set by electron-impurity scattering,
with scattering frequency (Itoh & Kohyama, 1993; Potekhin et al., 1999)

νeQ =

4πe4 ne Qimp
ΛeQ
p2 vF,e Z

(2.33)

F,e

≈ 3 × 1018 s−1

ρ14 Ye 1/3 Qimp ΛeQ
0.05
Z

.

(2.34)

scaled to typical values of the inner crust. Here pF,e and vF,e are the Fermi momentum and velocity
of the electrons, respectively, ΛeQ is the Coulomb logarithm, Ye is the electron fraction, and ρ14 ≡
ρ/(1014 g cm−3 ). The quantity Qimp ΛeQ / Z is of order unity in the inner crust. The resulting
thermal conductivity is
2

Ke =

π EF,e kB T c
Z
12
Qimp ΛeQ
e4

ρ Ye 1/3
Z
≈ 4 × 1019 erg s−1 cm−1 K−1 T 8 14
,
0.05
Qimp ΛeQ

(2.35)

where T 8 ≡ T/(108 K).
In the inner crust the temperature is T

109 K and typically below the critical temperature for

neutron superfluidity and the free neutron gas is superfluid. The presence of a neutron superfluid
in the neutron star crust is an idea predating the discovery of neutron stars (Migdal, 1959), and
there is both theoretical (see, e.g., Gezerlis & Carlson, 2010) and observational evidence (e.g.,
Shternin et al. 2007) that the neutrons are below their superfluid critical temperature in the deep
crust. The critical temperature T c of the 1 S0 neutron singlet pairing gap is expected to increase
from zero near the neutron drip density ρ ≈ 4 × 1011 g cm−3 to a maximum value near T c

24

109 K

]
[

16
14
12
10
8
6
4
2
1013

core

[

1014

]

1015

Figure 2.6: Superfluid critical temperature models for the neutron star inner crust as function of
mass density. The three solid gray curves are from Gandolfi et al. (2008). The solid black curve is
from Schwenk et al. (2003). The dashed black curve is from Gezerlis & Carlson (2011). The red
curve is from Chen et al. (1993). The blue curve is from Wambach et al. (1993).
before decreasing again at high mass densities near ρ ≈ 1014 g cm−3 where the repulsive core of
the neutron interaction removes the tendency to form pairs. Theoretical models for the superfluid
critical temperature T c (ρ) are shown in Figure 2.6.
The density range spanned by the neutron superfluid is an important input for pulsar glitch
models; for example, a recent study of pulsar glitches suggests that the neutron superfluid extends
from the crust into the core continuously in order to supply adequate inertia for pulsar glitches
(Andersson et al., 2012). Furthermore, the neutron superfluid plays a role in relaxation times of
pulsar glitches (e.g., Baym et al. 1969; Link 2012) and the dissipation of rotational energy through
the coupling of the crust to the neutron superfluid is important for pulsar timing models (e.g., Pines
& Alpar 1984).
If the 1 S0 neutron singlet pairing gap closes in the crust, a layer of normal neutrons forms
before the crust-core transition. The thermal conductivity of the normal neutrons may become
25

comparable to the electron thermal conductivity which will impact the cooling of neutron star
transients (to be discussed in Chapter 4). The neutron thermal conductivity is

Kn =

2T
π2 nn kB

3mn νn

(2.36)

where mn is the neutron effective mass and νn is the scattering frequency from neutron-neutron
and neutron-cluster scattering. We can derive the neutron scattering frequencies in the relaxation
time approximation where the scattering frequency can be expressed as (Flowers & Itoh, 1976;
Potekhin et al., 1999)
mn nion 2kF,n
νn =
dq q3 |V(q)|2 S κ (q) ,
3
3
nn 0
12π

(2.37)

where q is the momentum transfer, pF,n = (3π2 nn )1/3 ≡ kF,n is the neutron Fermi momentum,
and V(q) is the the Fourier transform of the scattering potential. The scattering medium is described
by the structure function



2  k2


S (q, ω)
ω
1
dω ω
 F,n
1 +
 ,
−
S κ (q) =
3
 2


2π
k
T
1
−
exp(−
ω/k
T
)
k
T
2
q
−∞
B
B
B
∞

(2.38)

which is written in terms of the dynamical structure factor S (q, ω).
To describe neutron-nucleus scattering in the inner crust, we assume that the nuclei are spherical and that the surface thickness is negligible compared to the size of the nucleus. Although the
nuclei in the pasta phase are certainly non-spherical, a description of scattering in non-spherical
geometries is beyond the scope of this work. Under these assumptions, the potential seen by the
neutrons can be modeled as a square well with V(r < RA ) = V0 , where RA is the radius of the
scattering center. The depth of the potential V0 ≈ Vin − Vout , where Vin and Vout are the neu26

tron single particle potentials inside and outside the scattering structures, respectively. In the pasta
phase, the density contrast between the scattering structure and the background rapidly decreases
with increasing density, implying a correspondingly rapid decrease in V0 and reduced neutron
scattering.
With the spherical assumption, the effective neutron-nucleus potential in momentum space is
4πR3A

F A (qRA ) ,

(2.39)

3[sin(x) − x cos(x)]
.
x3

(2.40)

Vn,A (q) = V0

3

with a form factor (Flowers & Itoh, 1976)

F A (x) =

The form factor F A → 1 in the limit that momentum transfers are small (x = qRA

1) and is

suppressed when momentum transfers are large. Inserting Equation 2.39 into Equation 2.37, we
find the neutron-phonon scattering frequency

νn,phn =

4 mn c2 nion V0 RA 2
Λn,phn ,
27π
nn
c

(2.41)

where the Coulomb logarithm is given by

Λn,phn =

2kF,n RA
0

phn

dx x3 F A2 (x) S κ (q = x/RA ) .

(2.42)

We evaluate the integral in Equation 2.42 using a Runge-Kutta scheme of order 8(5,3) (Hairer
phn

et al., 1993) and fitting formulae for S κ (q) (Potekhin et al. 1999; Equations 21 and 22) that were
developed in the context of electron-phonon scattering.

27

We find the frequency of neutron-impurity scattering using a similar approach. We assume that
the impurities are uncorrelated elastic scatterers, and write the scattering frequency as a sum over
all impurity species. The neutron-impurity potential for an impurity of species j with radius R j is
 3

 R j F A (qR j )

4πR¯ 3
 ,
¯ 
F A (qR)
−
1
Vn, j = V0

 R¯ 3 F A (qR)
¯
3

(2.43)

where R¯ is the radius of the average ion in the lattice. With this assumption, the dynamical structure
factor for impurity scattering is (Flowers & Itoh, 1976)

S imp (q, ω) =

1
nion

2πn j δ(ω) ,

(2.44)

j

where j is the sum over impurity species and n j is the number density of impurities.
imp

Upon using Equations (2.44) and (2.38) to obtain S κ (q) =

imp
j n j /nion , and inserting S κ (q)

and Vn, j (Equation 2.43) into Equation (2.37), we find the neutron-impurity scattering frequency
4 mn c2 nion V0 R¯ 2
νnQ =
ΛnQ Q.
27π
nn
c

(2.45)

Here we define the Coulomb logarithm for neutron-impurity scattering,

ΛnQ =

2kF,n R¯
0

dx x3 F A2 (x),

(2.46)

and the impurity parameter for neutron scattering,

Q=

1
ΛnQ 0

2kF,n R¯

 3
2


R
¯

n
F
(xR
/
R)
j  j A
j
3
2

dx x F A (x)
− 1 .
 3
¯
n
F A (x)
R
j ion

28

(2.47)

For scattering involving momentum transfers q

¯ ≈ 1. Taking
1/R j the ratio F A (qR j )/F A (qR)

R3j ∝ Z j and R¯ 3 ∝ Z then gives Q ≈ Qimp / Z 2 where Qimp (Equation 2.32) is the impurity
parameter for electron scattering. The neutron-impurity scattering frequency is therefore
4 mn c2 nion V0 R¯ 2 Qimp
ΛnQ .
νnQ ≈
27π
nn
c
Z 2

(2.48)

Since the neutron chemical potentials inside and outside the nucleus are required to be equal
in Gibbs equilibrium, we can estimate V0 as the difference in the single particle kinetic energies
inside and outside the nucleus,

V0 ≈

2 (3π2 n )2/3
in

2mn



2/3 

n

n
1 −
 ,
nin

(2.49)

where nin is the neutron number density inside the nucleus. We take RA to be the proton radius of
the nucleus given by (4π/3)R3A nin = Z, where Z is the proton number of the nucleus. We therefore
expect that V0 RA / c ∼ O(1) in the inner crust. The total scattering frequency is
4 mn c2 nion V0 RA
νn = νn,phn + νnQ ≈
27π
nn
c
mn
= 6.7 × 1020 s−1
mn

2

nion /nn
0.01

Qimp

ΛnQ
Z 2
Qimp
V0 R A 2
ΛnQ (.2.50)
Λn,phn +
c
Z 2

Λn,phn +

The ratio of the phonon and impurity scattering frequencies is
νn,phn
νnQ
and is typically of order unity for Qimp

=

Z 2 Λn,phn
Qimp ΛnQ

10.

29

(2.51)

Scaling the thermal conductivity of normal neutron to typical values in the inner crust,
2
−1
Qimp
nn
c 2
9π3 nn kB T
ΛnQ
Kn =
Λn,phn +
4 mn mn c2 nion V0 RA
Z 2
−2
m
nn /nion
≈ 3 × 1017 erg s−1 cm−1 K−1 T 8 Yn ρ14 n
mn
100
−1
Qimp
c 2
×
Λ
.
Λn,phn +
nQ
V0 RA
Z 2

(2.52)

In this expression mn = pF,n ∂ε(p)/∂p −1
p=pF,n is the Landau effective mass and ε(p) is the neutron
single particle energy including the rest mass (see, e.g., Baym & Chin, 1976). The dimensionless
quantity V0 RA / c is of order unity and measures the strength of the neutron–nucleus interaction;
RA is the typical size of the scattering structure and energy V0 is the magnitude of the scattering
potential.
The specific heat capacity of normal neutrons is

CVn

2
π2 nn kB T
=
ρ pF,n vF,n

≈ 3 × 104 erg g−1 K−1 Yn1/3 ρ−2/3
14

T7
3

,

(2.53)

scaled to typical values of the inner crust, where nn is the number density of normal neutrons, pF,n
is the neutron Fermi momentum, and vF,n is the neutron Fermi velocity.

2.4.3

Heat Sources and Sinks

The emissivity terms ν and nuc in Equation 2.25 represent thermal energy loss via neutrinos
and thermal energy deposition via nuclear reactions, respectively. During active accretion, heat
is deposited by non-equilibrium nuclear reactions (Sato, 1979; Bisnovaty˘i-Kogan & Chechetkin,

30

1979) which deposit ≈ 0.2 MeV per accreted nucleon in the outer crust (Gupta et al., 2007) and
≈ 1.5 MeV per accreted nucleon in the inner crust (Haensel & Zdunik, 1990, 2003, 2008). When
we model the thermal evolution of the crust in neutron star transients in Chapter 4, we parametrize
nuclear heating during active accretion rather than resolve individual reaction layers, following the
approach of Brown (2000).
We parametrize the heating expected from electron captures in the outer crust by distributing
heat over the pressure P = 7.2 × 1026 erg cm−3 to P = 1.3 × 1030 erg cm−3 such that the total
heat deposited is ≈ 0.2 MeV per accreted nucleon. To parametrize pycnonuclear fusion heating
in the inner crust we distribute heat over the pressure interval P = 2.6 × 1030 erg cm−3 to P =
1.5 × 1031 erg cm−3 such that the total heat deposited is ≈ 1.5 MeV per accreted nucleon. In
neutron stars that require extra heating to match quiescent temperatures, we include the additional
nuclear emissivity of nuc = mQ
˙ shallow proportional to the accretion rate m,
˙ where Qshallow is
in MeV per accreted nucleon. For example, the neutron star transient MXB 1659-29 requires
Qshallow ≈ 1 MeV per accreted nucleon to reproduce quiescent temperatures.
Neutrino cooling in the crust is primarily due to neutrino-pair production. For neutrino losses
due to neutrino-pair production we use the neutrino emissivity from Yakovlev et al. (2001) (their
Equation 22). At ocean and crust temperatures T

2×109 K neutrino losses from begin to balance

crust heating. For instance, MAXI J0556-332 reaches a nearly maximally hot state where neutrino
losses offset the heating from extra shallow heating in the crust (Deibel et al., 2015). In addition,
neutrino emission from Urca cycling nuclei pairs (discussed in Chapter 5) limits the ocean and
crust temperature to T

8 × 108 K near the depths inferred for superburst ignition (Deibel et al.,

2016).

31

Chapter 3
Magnetar Starquakes
Magnetars are isolated neutron stars with large magnetic fields between Bdip ∼ 1014 –1015 G,
as determined from their spin periods and spin-down rates assuming a dipole magnetic field1 .
Magnetars were initially classified as either Anomalous X-ray Pulsars (AXPs) or Soft Gamma-ray
Repeaters (SGRs), as determined by emission in either the soft X-ray or hard X-ray/soft γ-ray,
respectively (see Mereghetti (2008) for a review).
Powerful γ-ray flares from some magnetars (Barat et al., 1983) are thought to occur after a
reconfiguration of the strong magnetic field fractures the crystalline crust, triggering a powerful
burst of energy (Thompson & Duncan, 1995; Schwartz et al., 2005). In this model, crust fracturing
not only transports energy into the magnetosphere, where it is reemitted as the giant flare, but
energy is also deposited into the rigid crust causing a starquake (Duncan, 1998). If this is indeed the
case, then quasi-periodic oscillations (hereafter QPOs) observed in the tail of giant flare emission
(Israel et al., 2005; Strohmayer & Watts, 2005; Watts & Strohmayer, 2006) may be interpreted
as torsional oscillation modes of the crust during a starquake (Piro, 2005; Strohmayer & Watts,
2006; Samuelsson & Andersson, 2007). Crust oscillations, however, do not reproduce the lowest
frequency QPOs (Samuelsson & Andersson, 2007). Alternative theories explain QPOs as amplified
core magnetohydrodynamic modes (Glampedakis et al., 2006; Levin, 2006), but these models do
not produce the highest frequency QPOs (van Hoven & Levin, 2012).
1 McGill SGR/AXP Online Catalog,
http://www.physics.mcgill.ca/˜pulsar/magnetar/main.html

32

Within the pure crust mode paradigm for the observed QPOs, studies of starquakes often assume that superfluid neutrons in the neutron star inner crust (Migdal, 1959; Gezerlis & Carlson,
2010; Shternin et al., 2007; Brown & Cumming, 2009) near the crust-core phase transition eliminate viscous drag (Ruderman, 1968). In this case, QPOs arise from pure crust modes and can
therefore constrain crust physics (Steiner & Watts, 2009; Deibel et al., 2014). In particular, as
we discuss in this chapter, crust modes constrain the nuclear symmetry energy at mass densities
near nuclear saturation. Furthermore, because torsional mode frequencies are sensitive to the crust
thickness — which in turn is set by the neutron star’s mass and radius — modeling QPO frequencies constrains the magnetar’s mass and radius.

3.1

A Magnetized Crust

Magnetars are isolated non-accreting neutron stars and it is assumed that their crust composition
has not been polluted from the pristine composition frozen in post-supernova. Therefore, the equilibrium crust composition outlined and computed in Chapter 2 will be a sufficient representation of
the magnetar’s crust. However, the high magnetic fields of magnetars alter the crust composition
at low densities because the electrons become quantized into energy levels above the critical magnetic field (see Equation 3.1). These modifications will alter the energy density contribution from
electrons and consequently the nucleus that minimizes the total energy density in the calculation
for the equilibrium crust.

3.1.1

Electron Quantization

At the magnetic field strengths inferred for magnetars the ocean and crust composition is altered
at low densities by the magnetic field. The electron contribution to the energy density becomes
33

comparable to the other energy density contributions at low densities when the electron cyclotron
energy equals the rest mass energy of the electron, ωc = me c2 , giving a field

Bcrit =

me 2 c3
≈ 4.4 × 1013 G ,
e

(3.1)

and we define the quantity,
B∗ = B/Bcrit ,

(3.2)

which is a convenient quantity when discussing the large magnetic fields of magnetars. In magnetars where the magnetic field is B

Bcrit we see the composition is altered at low densities, as

seen in Table 3.1. When the magnetic field is B
greater than the neutron-drip point (ρ

Bcrit the composition can be altered at densities

ρdrip ≈ 4 × 1011 g cm−3 ). The electron contribution to

the energy density is that of an electron gas embedded in a perpendicular magnetic field B = Bˆz.
Following the work of Broderick et al. (2000), the electrons acquire an effective mass, me , in the
presence of the magnetic field

1 1
2
me = m2e + 2 x + − ν eB ,
2 2

(3.3)

where me , x, and ν are respectively the electron mass, principal quantum number, and electron
spin (Rabi, 1928; Ventura & Potekhin, 2001). The additional term in the electron mass designates
the Landau level which is set almost entirely by the strength of the magnetic field. The maximum
Landau level is found using the fermi momentum of the electron

2

2 −m
kF 2 = EF,e
e ,

34

(3.4)

where kF is the electron Fermi momentum and EF,e is the electron Fermi energy. In the neutron star
crust the electrons are degenerate and relativistic, and EF,e ≈ µe

me c2 , where µe is the chemical

potential of electrons. The maximum Landau level is set as the principal quantum number before
kF 2 becomes negative. The fermi momentum can be found by using the known electron number
density, ne = Zn/A, and solving for the chemical potential after summing over the total number of
Landau levels within the equation

ne =

eB
kF .
2π2 σz n

(3.5)

Knowing the electron fermi momentum we can easily obtain the electronic contribution to the
energy density per baryon,

e (Z, A, n; ne ) =

kF
eB
k2
2
4π 0

k2 + me 2 dk .

(3.6)

Note that we use the relativistic treatment of electrons because magnetar magnetic field strengths
near B

Bcrit imply that the electrons are relativistic ( kF,e

3.1.2

Composition

Electron quantization in a strong magnetic field (B

me c2 ).

Bcrit ) alters the energy density of electrons.

A change in the electron contribution to the total energy density changes the equilibrium nucleus at
a given mass density in the crust. A revised crust composition in various magnetic fields B

Bcrit

are shown in Table 3.1.
The crust composition sets the shear modulus µ, which encodes the crust’s response to elastic

35

Table 3.1: Equilibrium nuclei below the crust-core transition in a magnetized crust.
ρmax (g cm−3 )
Nuclei1

B∗ = 0

64 Ni
28
66 Ni
28
84 Se
34
82 Ge
32
80 Zn
30
78 Ni
28
76 Fe
26
122 Zr
40
120 Sr
38
118 Kr
36
116 Se
34
114 Ge
32
112 Zn
30
110 Ni
28
166 Zr
40

B∗ = 1

B∗ = 10

B∗ = 102

B∗ = 103

2.23 × 108

2.33 × 108

1.63 × 109

1.75 × 1010

1.37 × 109

1.40 × 109

2.92 × 109

2.71 × 1010

5.66 × 109

4.87 × 109

5.29 × 1010

1.73 × 1010

1.69 × 1010

7.62 × 1010

3.99 × 1010

3.94 × 1010

1.01 × 1011

1.56 × 1011

1.57 × 1011

1.61 × 1011

1.86 × 1011

1.85 × 1011

1.76 × 1011

2.51 × 1011

2.52 × 1011

1.98 × 1011

3.54 × 1011

3.54 × 1011

4.04 × 1011

5.17 × 1011

5.15 × 1011

5.77 × 1011

8.11 × 1011

8.13 × 1011

8.56 × 1011

2.35 × 1012

2.25 × 1012

3.94 × 1012

4.02 × 1012

8.64 × 1012

8.65 × 1012

1.07 × 1013

1.08 × 1013

36

Figure 3.1: Equilibrium composition of the neutron star crust in a strong magnetic field B =
2 × 1015 G. (Left panel): Equilibrium crust composition without a magnetic field. (Right panel):
The black arrow indicates the first transition between equilibrium nuclei, which has moved to a
greater depth as a consequence of the strong magnetic field when compared to the left panel.
stresses. We use the shear modulus of a body-centered cubic lattice (Strohmayer et al., 1991),

µ=

0.1194Γ
ni kB T ,
1 + 0.595(Γ0 /Γ)2

(3.7)

where Γ is the plasma coupling parameter of the crust (Equation 2.4) and Γ0 ≡ 175 is the melting
point of the crust. The speed of shear waves in the crust v s =

3.2

µ/ρ is set by the shear modulus.

Crust Torsional Modes

Magnetar giant flares have luminosities Lpeak ∼ 1044 −1047 ergs s−1 (Barat et al., 1983). Operating
under the assumption that a fraction of the energy is transported to the crust, the crust will receive
a strong impulse that will drive oscillation modes, shown schematically in Figure 3.2. The intense

37

Figure 3.2: Schematic of a torsional oscillation mode in the magnetar crust. Crust oscillation
frequencies are found using Equation 3.14, using the boundary condition ξ = 0 (loss of traction)
at the base and top of the crust.
gravity and vertical stratification confine the modes to the axial direction and the mode spectrum
can be well described by a differential equation of motion for axial perturbations. What arises is an
eigenvalue problem, where the solution of the differential equation are the eigenmodes of the crust.
To produce the most realistic solutions we add general relativistic effects and magnetic effects to
our differential equation for torsional oscillations.
The axial perturbation equation, which describes axial torsional oscillation modes, is derived by
Schumaker & Thorne (1983) and applied to modern neutron star models in Samuelsson & Andersson (2007). Here we follow the formalism of Samuelsson & Andersson (2007) when describing
the axial perturbation equation in the absence of a magnetic field, which takes the form

ξ + F ξ + Gξ = 0 ,

(3.8)

where ξ is an axial perturbation and primes indicate derivatives with respect to the radial coordinate. The functions F and G depend on the elastic properties of the crust and the space-time

38

metric,
F = ln r4 eν−λ (ε + pt ) v2r

,

(3.9)

and


e2λ  −2ν 2 v2t ( − 1) ( + 2) 
G = 2 e ω −
 ,
r2
vr
where r is the radial coordinate,

(3.10)

is the energy density, p is the pressure, ω is the oscillation

angular frequency, and l is the angular wave number. Here we use a static, spherically symmetric
space-time metric,
ds2 = −e2ν dt2 + e2λ dr2 + r2 dθ2 + sin2 θ dφ2 ,

(3.11)

and the perturbation equation may be expressed as (Samuelsson & Andersson, 2007),

ξ + ln r4 eν−λ (ε + pt ) v2r



e2λ  −2ν 2 v2t ( − 1) ( + 2) 
ξ + 2 e ω −
 ξ = 0 ,
r2
vr

(3.12)

where e2λ ≈ e−2ν = (1 − 2M(r)/r)−1 in the crust. The large magnetic fields in magnetars drive
magnetic waves, or Alfv´en waves, which can significantly alter crust dynamics where vA > vs
(Piro, 2005). This is the case for the magnetic field of SGR 1806-20 (B = 2.4 × 1015 G) where
vA

vs in the outer crust, as seen in Figure 3.3. Therefore, we add corrections for a finite Alfv´en

velocity to the perturbation equation. The Newtonian expression for axial perturbations with a
magnetic field are (Piro, 2005; Steiner & Watts, 2009),

 
 
(µξ )
2
+ vA ξ + ω2 1 +
ρ



v2A  (l2 + l − 2)µ 
 −
 ξ = 0 ,
c2
ρr2

(3.13)

where A = B/ 4πρi . Combining the general relativistic expression and the Newtonian magnetic

39

0.010

L = 46 MeV
B = 2.4 ⇥ 1015 G

0.006

vs

0.002

0.004

v/c

0.008

vA

1010

1012

1011

log10 (⇢) [g cm

1013
3

1014

]

Figure 3.3: Alfv´en velocity (blue curve) and shear velocity (black curve) in the crust as a function
of mass density. The composition is that of a 1.4 M neutron star using the SLy4 crust EOS.

40

36
●

L = 46 MeV
nt = 0.12 fm

●

3

●
●
●

●

32

●
●
●

●

-2

●

●

30

●
●

●

●
●

●

28

SGR 1806

●

●

20

●

●
●

●

●
●
●

26

Fundamental Frequency [Hz]

34

+2

●
●

●

●
●

●

24

●

●
●
●

0.8

1.0

1.2

1.4

1.6

1.8

2.0

Magnetar Mass [M ]
Figure 3.4: Frequency of the fundamental l = 2 mode as a function of the magnetar mass for the
core EOS probability distribution (centroid and ±2σ) from Steiner et al. (2013) and an SLy4 crust
EOS. The dashed black line indicates the observed 29 Hz QPO of SGR 1806-20. The frequencies
are evaluated for a crust-core transition density of 0.12 fm−3 with B = 0 G.
expression,




d
ln r4 eν−λ (ε + p) 2s ξ + e2λ e−2ν ω2 1 +
( 2s + 2A )ξ + 2s
dr

2
A 
−
c2 




l2 + l − 2 2s 
 ξ = 0.

2
r
(3.14)

In the limit that B → 0, Equation 3.12 is recovered. In the limit that e2λ → 1, e−2ν → 1, the
Newtonian expression (Equation 3.13) is recovered. A full derivation of the perturbation equations
can be found in Appendix 6.
In this study, we calculate crust modes for serveral crust equations of state. The modes will
41

differ between crust equations of state due to a changing shear modulus and crust thickness. We
integrate Equation 3.14 over the crust, that is, from r = Rcore the core radius to the top of the
crust r = R175 where the crust melts where Γ = 175, given by the ocean-crust transition density in
2.5. The traction ξ vanishes at the top of the crust where pressure vanishes, and we also assume
neutron superfluidity at the crust-core interface causes of a loss of traction there by eliminating
viscous drag (Ruderman, 1968).

3.3

Constraints from Observations

Each neutron star mass and radius give a unique crust thickness and therefore a unique set of
torsional crust modes. Because of the orthogonal behavior of fundamental and overtone crust
modes, matching the observed fundamental and harmonic modes from a magnetar can only be
achieved by one mass and radius combination for a given equation of state, an example is shown
in Figure 3.4. Observations of photosphere radius expansion bursts and low-mass X-ray binaries
give empirical constraints on the neutron star core equation of state (Steiner et al., 2010; Steiner
et al., 2013). We investigate which core equation of state contains crusts that host torsional modes
at similar frequencies to those observed in giant flare QPOs. In particular, we examine the giant
flares from the magnetars SGR 1806-20 (Strohmayer & Watts, 2006; Watts & Strohmayer, 2006)
and SGR 1900+14 (Strohmayer & Watts, 2005).

3.3.1

Nuclear Symmetry Energy

The different core equations of state from Steiner et al. (2013) are characterized by different degrees of stiffness of the nuclear symmetry energy. The nuclear symmetry energy, S (n), is defined
as the energy cost from departing from neutron-proton symmetric matter. Equations of state with
42

2.5

(a)

+2

L = 46 MeV
nt = 0.12 fm 3

-1
1.5

-2

●
●
●

➤
➤

0.0 G

●

SGR 1806-20

●
●

●
●

●

●
●

●
●

2.4⇥1015 G

●
●
●

●
●

1.0

0.0 G
1.0⇥1015 G

●

●

5.0⇥1015 G

●
●
●

●

0.5

Magnetar Mass [M ]

2.0

+1

0.0 G SGR 1806-20
11

➤

➤

10

12

13

14

Magnetar Radius [km]
Figure 3.5: Magnetar mass as a function of radius for the core EOS probability distribution from
Steiner et al. (2013). Frequencies are evaluated using the SLy4 crust EOS (L = 46 MeV) for
nt = 0.12 fm−3 . The thick red solid line indicates masses and radii for which the fundamental
mode has a frequency of 29 Hz in the case of SLy4 and 18 Hz in the case of Rs. The black shortdashed line indicates masses and radii for a 626 Hz harmonic mode and B = 0 G. Masses and radii
from 626 Hz harmonic modes with magnetized crusts are labeled accordingly. Arrows indicate
masses and radii that match both the fundamental and the harmonic modes for the field-free case
and the case with the magnetic field of SGR 1806-20 (B = 2.4 × 1015 G).
.

43

2.5

(b)

+2

L = 86 MeV
nt = 0.12 fm 3

-1
-2

●

1.5

●
●

➤

0.0 G

1.0

●

●

●

●

●

●
●

5.0⇥1014 G

●

1.0⇥1015 G

●

2.4⇥1015 G

●

●

●
●
●

●

●

●

0.5

●

0.0 G

●

●

●

SGR 1806-20 ●

➤

Magnetar Mass [M ]

2.0

+1

0.0 G

SGR 1806-20

11

12

➤

➤

10

13

14

Magnetar Radius [km]
Figure 3.6: Magnetar mass as a function of radius for the core EOS probability distribution from
Steiner et al. (2013). Frequencies are evaluated using the Rs crust EOS (L = 86 MeV) for nt =
0.12 fm−3 . The thick red solid line indicates masses and radii for which the fundamental mode
has a frequency of 29 Hz in the case of SLy4 and 18 Hz in the case of Rs. The black short-dashed
line indicates masses and radii for a 626 Hz harmonic mode and B = 0 G. Masses and radii from
626 Hz harmonic modes with magnetized crusts are labeled accordingly. Arrows indicate masses
and radii that match both the fundamental and the harmonic modes for the field-free case and the
case with the magnetic field of SGR 1806-20 (B = 2.4 × 1015 G).

44

different values of the symmetry energy will have varying degrees of stiffness/softness near the nuclear saturation point. The degree of stiffness is determined by the first derivative of the symmetry
energy at nuclear saturation, L ≡ 3n0 (∂S /∂n)n=n0 , where n0 = 0.16 fm−3 . Different values for L
lead to different mass-radius relations for the neutron star. As a result, studies that model magnetar
QPOs as torsional oscillations of the crust can constrain the nuclear symmetry energy (Steiner &
Watts, 2009; Sotani et al., 2012, 2013b; Deibel et al., 2014; Iida & Oyamatsu, 2014; Sotani et al.,
2015).
Magnetar masses and radii that support crust frequencies consistent with observed QPOs can
be seen in Figures 3.5 and 3.6 for the SLy4 crust equation of state with L = 46 MeV and the Rs
equation of state with L = 86 MeV, respectively. Furthermore, two magnetic fields, B = 0 and
2.4 × 1015 G — the inferred dipole magnetic field of SGR 1806-20 — are tested to illustrate the
impact of the magnetic field on determining masses and radii with this method. Both values of L
give magnetar masses and radii that support crust modes consistent with observed QPOs. A softer
equation of state with L = 46 MeV, however, finds magnetar masses and radii more consistent with
the best fits from other empirical constraints from Steiner et al. (2013), when a realistic 2.4×1015 G
magnetic field is included in the crust oscillation calculation. Although an equation of state with
L = 86 MeV may reproduce the observed oscillations with a realistic magnetic field, it requires
magnetar masses and radii well outside the 2σ boundaries from empirical constraints.

3.3.2

Neutron Entrainment

A fraction of the quasi-free neutron gas in the inner crust may be entrained by the nuclear lattice
during a crust oscillation (Chamel et al., 2013a; Chamel, 2013). If this is the case, an axial perturbation will move more mass and the shear velocity of crust oscillations will decrease. This will
in turn alter the oscillation frequencies of the crust (Sotani et al., 2013a). The impact of neutron
45

entrainment on fundamental and overtone crust frequencies can be seen in Figures 3.7 and 3.8 for
the SLy4 crust equation of state with L = 46 MeV and the Rs equation of state, respectively. Here
fent indicates the fraction of quasi-free neutrons entrained by the crust lattice during an oscillation.
We find that the fraction of entrained neutrons has a negligible effect on the harmonic modes of
the crust because harmonic mode energy is uniformly distributed throughout the entire crust (Piro,
2005). The fundamental mode frequency, however, is significantly altered by neutron entrainment,
likely because the fundamental mode’s energy is concentrated in the deepest part of the crust.
Using the SLy4 equation of state, entrainment fractions fent

0.75 are needed to have predicted

crust modes match observed QPOs in SGR 1806-20, as can be seen in Figure 3.7. Entrainment
fractions fent

0.75 require larger magnetar masses and radii that are outside of the best fits from

empirical constraints from observations of PREs and LMXBs (Steiner et al., 2013). For the Rs
equation of state, low entrainment fractions fent

0.3 are preferred, as can be seen in Figure 3.8.

Note that values of fent ≈ 0.35 − 0.90 have been suggested by a study of Bragg scattering with the
crust’s crystal lattice (Chamel et al., 2013b).
Generally, the degree of neutron entrainment alters fundamental and harmonic crust modes
by changing the shear velocity in the crust. Crust equations of state with larger values of L have
smaller fundamental mode frequencies; a similar result was found by Sotani et al. (2013a). Furthermore, these authors also find that a higher entrainment fraction decreases fundamental mode
frequencies (note that N s /Nd = 1 − fent in the notation of Sotani et al. 2010). Our results are only
for equations of state fit to observations of neutron stars (Steiner et al., 2013) and we can not fully
determine the correlation of L and fent using our equations of state with discrete L values.

46

2.5

+2
+1

●

2.0

●
●

➤

0.75

0.50

●

-2

●
●
●

●
●

0.75

●

1.0

●

1.0
0.50

●

1.0

0.5

1.0

1.5

-1

●

0.50

➤
➤

Magnetar Mass [M ]

L = 46 MeV
nt = 0.12 fm 3

1.0
12

➤

11

0.50

➤

➤

10

0.75
13

14

Magnetar Radius [km]
Figure 3.7: Magnetar mass as a function of radius for EOS probability distribution from Steiner
et al. (2013). Frequencies are evaluated using the SLy4 crust EOS with B = 0 G and nt =
0.12 fm−3 . The red dot-dashed, blue dotted, and black dashed lines indicate masses and radii from
fundamental modes of frequency 29 Hz for different free neutron entrainment fractions fent . The
shaded band indicates masses and radii from 626 Hz harmonic modes as fent is varied from 0.50
to 1.0. Arrows indicate the masses and radii that match both the fundamental and the harmonic
modes for fent = 1.0, 0.75, and 0.50.

47

2.5

+2

L = 86 MeV
nt = 0.12 fm 3

-1
-2

●

1.5

●
●

1.0

●

➤
➤
➤

0.20

●

0.20

●

0.25

●

0.25

●
●

●

0.30

●

0.30

●
●

●

0.20

●

0.30

0.5

Magnetar Mass [M ]

2.0

+1

0.30

0.25

11

12

➤

➤

➤

10

0.20
13

14

Magnetar Radius [km]
Figure 3.8: The same as Fig. 3.7, but for the Rs crust EOS with B = 0 G. Here the free neutron
entrainment fraction fent is varied from 0.20 to 0.30, with fent labeled next to the corresponding
curves. The red dot-dashed, blue dotted, and black dashed lines indicate masses and radii from the
29 Hz fundamental mode. The shaded band indicates masses and radii from the 626 Hz harmonic
mode. Arrows indicate the masses and radii for fent = 0.30, 0.25, and 0.20 that match both the
fundamental and harmonic modes.

48

3.3.3

Crust-core Transition Depth

The crust’s harmonic frequencies are sensitive to the location of the crust-core transition. Harmonic
crust modes scale with the crust thickness and therefore a change in the depth of the crust-core
transition will change the crust thickness and harmonic mode frequencies (Deibel et al., 2014).
The pressure of the crust-core transition is unknown (Newton et al., 2013).
To investigate the impact of the crust-core transition density on mode frequencies, we test
the crust-core transition density values nt = 0.08 and 0.12 fm−3 to represent the range of values
obtained in Oyamatsu & Iida (2007). The impact of the crust-core transition density on crust
fundamental and overtone modes can be seen by comparing Figures 3.5 and 3.9.
For a crust-core transition density of nt = 0.12 fm−3 , crust oscillation frequencies match observed QPOs for SGR 1806-20 for a magnetar mass M = 1.25 M and radius R = 12.4 km. A
shallow crust-core transition at nt = 0.08 fm−3 does not give a mass and radius consistent with
constraints from observations of low-mass X-ray binaries (Steiner et al., 2013), requiring a magnetar mass M = 0.96 M and radius R = 13.5 km, as can be seen in Figure 3.9. A small magnetar mass and large magnetar radius is also needed in the case of the Rs equation of state with
nt = 0.12 fm−3 . For this equation of state, predicted crust oscillations match observed QPOs for a
magnetar mass M = 1.10 M and radius R = 13.8 km, as can be seen in Figure 3.6.

3.3.4

Crust Magnetic Field

In crust mode calculations at B = 0 G we assume that the crust thickness remains constant (i.e.,
Rcrust − Rcore = ∆R is constant). In large magnetic fields B

Bcrit the crust thickness decreases

with increasing magnetic field because the melting point (Equation 2.4) moves to greater depths.
Furthermore, the magnetic field alters the composition, as seen in Table 3.1, which alters the shear

49

2.5

+2

L = 46 MeV
nt = 0.08 fm 3

-1
1.5

-2
●
●

0.0 G
15
● 1.0⇥10
G

●

➤

1.0

➤

●
●
●

0.0 G
SGR 1806-20

●
●

●
●

●

●

●

●

●

●

●

●

2.4⇥1015 G

●

0.5

Magnetar Mass [M ]

2.0

+1

0.0 G SGR 1806-20
11

12

➤

➤

10

13

14

Magnetar Radius [km]
Figure 3.9: The same as Fig. 3.5, but for the SLy4 crust EOS with nt = 0.08 fm−3 . The thick
red solid line indicates masses and radii determined from a fundamental mode of 29 Hz. Masses
and radii from harmonic modes with magnetized crusts are labeled accordingly. Arrows indicate
masses and radii that match both the fundamental and the harmonic modes for the field-free case
and the case with the magnetic field of SGR 1806-20 (B = 2.4 × 1015 G).

50

modulus. We find that large magnetic fields, and the resulting decrease in crust thickness and
change in composition, have a negligible effect on crust fundamental modes, but may significantly
impact crust harmonic modes.
The difference between crust frequencies in a crust with a magnetic field of B = 0 G and
B = 2.4 × 1015 G can be seen in Figures 3.5 and 3.6. The fundamental crust mode for both crust
magnetic fields are identical — similar to the results found by Nandi et al. (2012). The fundamental
mode energy is concentrated at the base of the crust (Piro, 2005) where the Alfv´en velocity is small
compared to the shear velocity and the fundamental mode frequency is mostly determined by the
radius of the neutron star (i.e., our choice of equation of state). Harmonic modes have their energy
spread over the crust and are sensitive to the Alfv´en velocity at lower densities, which is large
compared to the shear velocity for the magnetic fields examined here, as can be seen in Figure 3.3.
Therefore, harmonic mode frequencies are greatly impacted by the choice of magnetic field; the
SLy4 equation of state requires B

5 × 1015 G and the Rs equation of state requires B

1015 G to

have magnetar masses and radii consistent with predictions from low-mass X-ray binaries (Steiner
et al., 2013).
The sensitivity of the harmonic mode to the crust magnetic field is worthy of note. If observed
QPOs are indeed crust oscillation modes, crust magnetic fields must not be much greater than
inferred dipole magnetic fields, otherwise QPOs can only be matched by magnetar masses and
radii that are inconsistent with observations of PREs and LMXBs (Steiner et al., 2013). In order
to match QPOs with smaller magnetar radii, for instance R < 11.1 km as found by Guillot et al.
(2013), requires larger values of L and fent .

51

Chapter 4
Neutron Star Transient Thermal Evolution
In this chapter, we explore the thermal evolution of a neutron star transient’s outer layers during
and after an accretion outburst. Accretion outbursts have been observed for many neutron stars
in low-mass X-ray binaries (see Degenaar et al. (2015) for a summary). Accumulated accreted
material compresses existing envelope material deeper in the neutron star. Material compressed
into the neutron star crust induces non-equilibrium reactions (Bisnovaty˘i-Kogan & Chechetkin,
1979; Sato, 1979) that deposit ≈ 1–2 MeV per accreted nucleon of heat (Haensel & Zdunik, 1990,
2003; Gupta et al., 2007; Haensel & Zdunik, 2008). The neutron star enters quiescence when
accretion halts and the neutron star’s outer layers begin to cool back to thermal equilibrium with
the core. The quiescent emission of the neutron star is powered by the thermal emission from the
neutron star surface (Rutledge et al., 2002) as a portion of heat deposited during outburst diffuses
toward the surface. This thermal evolution model can be understood qualitatively as a cooling slab
(Section 4.1).
The cooling slab model can successfully explain cooling neutron star transients in quiescence
once realistic microphysics for the neutron star crust are included (Section 4.2). We apply the
thermal evolution model to examine the thermal evolution of the hottest neutron star transient
observed to date, MAXI J0556-332 (Matsumura et al., 2011; Homan et al., 2011; Sugizaki et al.,
2013; Homan et al., 2014), in Section 4.2.1. The high temperature reached in the outer layers of
MAXI J0556-332 make it an interesting test case for thermal evolution models and reveals heating
anomalies in the outer crust. In Section 4.3 we apply the thermal evolution model to the Z-source
52

transients — which trace out a distinctive “Z” shape in the X-ray color-color diagram (Hasinger &
van der Klis, 1989) — to create an observational test for the presence of g-modes in the neutron
star’s ocean. Finally, in Section 4.4, we use the late-time cooling of transients to examine the
thermal properties of the deepest regions of the inner crust near the crust-core transition.

4.1

A Cooling Slab

During an accretion outburst, heat is deposited in the neutron star’s outer layers by nuclear reactions
which brings the neutron star’s outer layers out of thermal equilibrium with the core. Given the
length of an observed outburst we can predict the temperature of the neutron star at the outset
of quiescence by following the time-dependent thermal diffusion through the neutron star’s outer
layers. We model the thermal evolution of the crust during outburst and quiescence by constructing
a cooling slab that follows the heat flux across the entire crust. The heat flux can be described by
Fick’s Law
F = −K∇T ,

(4.1)

where F is the heat flux, K is the thermal conductivity, and T is the temperature. The divergence
of the heat flux describes the change in local temperature as a function of time, that is, the timedependent thermal diffusion. The thermal diffusion equation can be written as

ρC

dT
= ∇ · (K · ∇T ) ,
dt

(4.2)

where ρ is the mass density, C is the specific heat per unit mass, and T is the local temperature. In

53

one dimension, the thermal evolution is described by

ρC

∂
∂T
dT
=
· K·
,
dt
∂x
∂x

(4.3)

where x is the spatial coordinate. This expression is generically the diffusion equation in one
dimension,
du
∂
∂u
∂D(u(x), x) ∂u 2
∂2 u
=
D(u(x), x)
=
+ D(u(x), x) 2 ,
dt
∂x
∂x
∂u
∂x
∂x

(4.4)

which can be finite differenced depending on the diffusion coefficients dependence on u. Advancing the solution for u with time can be achieved using the method of lines, which can be found in
Press et al. (2007). This method involves integrating the system of ordinary differential equations
formed by finite differencing the right hand side of Equation 4.4. Integrating the resulting equation
with respect to time for each zone’s set of ordinary differential equations gives the solution for u at
the next time step for that zone.
In the case of our cooling slab model of the neutron star crust, u represents the crust temperature
T and for our spatial coordinate x we choose the mass density in the crust ρ. Then, the diffusivity
D is
D(ρ) =

K
,
ρCV

(4.5)

where the thermal conductivity K ≈ Ke throughout the ocean and crust because electron conduction is the primary means of thermal transport.
Scaling to typical values of the inner crust, where electron-impurity scattering sets the thermal
conductivity of electrons, gives

D ≈ 4 cm2 s−1

Ye 1/3
Z
,
0.05Yn
Qimp ΛeQ

54

(4.6)

where Yn is the free neutron fraction, Z is the average proton number of the composition, Qimp
is the impurity parameter, and ΛeQ is the Coulomb logarithm for electron-impurity scattering
(Potekhin et al., 1999).

4.2

Realistic Cooling Models

Neutron star thermal evolution models that include accretion-driven heating during outburst naturally reproduce quiescent light curves of neutron star transients, for example in XTE J1702-462
(Fridriksson et al., 2010, 2011; Page & Reddy, 2013; Turlione et al., 2015) and EXO 0748-676
(Degenaar et al., 2009, 2014). For example, in Figure 4.1 we show a model fit to the quiescent
cooling curve of MXB 1659-29 using the neutron star thermal evolution code dStar (Brown, 2015)
where the cooling slab model reproduces well the broken power-law behavior of the quiescent light
curve. The code solves the fully general relativistic heat diffusion equation using a method of lines
algorithm in the MESA numerical library (Paxton et al., 2011, 2013, 2015) and the microphysics
of the crust follows Brown & Cumming (2009).
We modeled the ≈ 2.5 year outburst in MXB 1659-29 (Wijnands et al., 2003, 2004) using a
local mass accretion rate m
˙ = 0.1 m
˙ Edd , where m
˙ Edd = 8.8 × 104 g cm−2 s−1 is the local Eddington
mass accretion rate. The model uses a neutron star mass of M = 1.6 M and radius of R = 11.2 km
that are consistent with the quiescent light curve fits from Brown & Cumming (2009). Whereas
Brown & Cumming (2009) held T b fixed during accretion to simulate the effect of a shallow heat
source, we instead include the heat source directly and allow T b to evolve as accretion proceeds.
The model includes a Qshallow = 1 MeV per accreted nucleon shallow heat source spread between
y = 2 × 1013 g cm−2 and y = 2 × 1014 g cm−2 (which will be discussed further in Section 4.2.1),
where y ≈ P/g is the column depth (Equation 2.3. For the crust composition we use the accreted

55

[ ]

140
120
100
80
60
40
101

102

[

103
]

104

Figure 4.1: Quiescent lightcurve of the neutron star MXB 1659-29. The model uses a neutron star
mass of M = 1.6 M , a neutron star radius of R = 11.2 km, a shallow heat source of Qshallow =
1.0 MeV per accreted nucleon, and a core temperature of T core = 4 × 107 K. The model also uses
an impurity parameter of Qimp = 2.5 for the entire crust as done in Brown & Cumming (2009).
composition from Haensel & Zdunik (2008) that assumes an initial composition of pure 56 Fe (see
their Table A3). The temperature at the top of the crust T b is mapped to the neutron star’s effective
temperature by a T eff − T b relation (Gudmundsson et al., 1982, 1983) updated for a mixed helium
and iron envelope (Brown et al., 2002).
The model of crust cooling in MXB 1659-29 uses a crust impurity parameter of Qimp = 2.5
and core temperature T core = 4 × 107 K, consistent with the fit from Brown & Cumming (2009).
The impurity parameter is assumed constant throughout the entire crust and was constrained to
Qimp < 10 in MXB 1659-29 (Brown & Cumming, 2009). In this model, the crust reaches thermal
equilibrium with the core by ≈ 1000 days into quiescence, and so predicts a constant temperature

56

]

1011
1010

[

109
108
107
106 7
10

108

[ ]

109

Figure 4.2: Location of the ocean-crust transition as a function of temperature at the base of the
ocean. The black solid curve is for the one-component accreted composition (Haensel & Zdunik,
1990) and the black dotted curve is for a multi-component accreted composition (Steiner, 2012).
The crystallization density is given by Equation (2).
at later times. The predicted quiescent light curve (Figure 4.1) using our thermal evolution model
agrees well with the predicted cooling curve from Brown & Cumming (2009).
Although heat diffuses toward the neutron surface during quiescence, the quiescent light curve
can be simply understood as a cooling front moving inward from the surface toward the crust-core
interface. Therefore, as the neutron star’s outer layers cool the light curve probes successively
deeper layers (Brown & Cumming, 2009) with increasingly longer thermal times. A cooling front
moving through the neutron star’s outer layers will move on a timescale corresponding to the
thermal time of the layer its in. The thermal time to reach a mass density ρ (Henyey & L’Ecuyer,
1969) is
1
τtherm =
4

2
ρCV 1/2
dz
,
z K

(4.7)

where the integration is carried out over the depth of the layer z, CV is the specific heat, and K is

57

the thermal conductivity. The thermal time for the ocean is then

τ∞
therm,liquid ≈ 1.2 days ρ9

g14 −2 Ye 3 Z
2
0.4
34

1+z
,
1.24

(4.8)

where Ye is the electron fraction, g14 is the neutron star gravity in units of 1014 cm s−2 , and
1 + z = (1 − 2GM/(Rc2 ))−1/2 is the redshift to the observer frame at infinity. At higher densities,
where ρ

109 g cm−3 , the heat capacity is set by the ions with CV ≈ CVion = 3kB /Amu (Equa-

tion 2.29). The thermal conductivity is set by lattice phonons with an electron-phonon scattering
ν = νep ≈ 13αkB T/ (Baiko & Yakovlev, 1995, 1996), where α = e2 / c ≈ 1/137 is the fine
structure constant. The resulting thermal time in the solid crust is

τ∞
therm,solid ≈ 2.1 days ρ9

g14 −2 Ye 2 56
2
0.4
A

1+z
,
1.31

(4.9)

as was previously shown in Brown & Cumming (2009) (their Equation 9). As can be seen from the
density dependence of the thermal time, the early quiescent light curve probes low densities in the
neutron star’s ocean and the late time light curve reveals the physics of the inner crust. As we will
show in the following sections, the early light curve in MAXI J0556-332 reveals a strong heating
anomaly in the shallow crust and the late-time light curve of MXB 1659-29 constrains the thermal
properties of the deep inner crust.

4.2.1

Heating Anomalies in MAXI J0556-332

Accretion-driven nuclear reactions deposit ≈ 1 MeV per accreted nucleon of crustal heating during
an accretion outburst; however, some sources require additional heating in the outer crust during
outburst to fit quiescent observations. For example KS 1731-260 (Wijnands et al., 2001, 2002) and

58

MXB 1659-29 (Wijnands et al., 2003, 2004; Cackett et al., 2008) require ≈ 1 MeV per accreted
nucleon (Brown & Cumming, 2009) during outburst to fit quiescent light curve observations. Although the source of extra heating is unknown, some heat may be supplied by compositionallydriven convection in the neutron star ocean (Medin & Cumming, 2011, 2014, 2015) and stabilized
helium burning in the neutron star envelope (Brown & Bildsten, 1998).
Here we discuss the hottest neutron star transient MAXI J0556-332 that is an exceptional and
interesting test case for our thermal evolution model. MAXI J0556-332 was discovered in outburst
in January 2011 (Matsumura et al., 2011). The source entered quiescence 16 months later at a high
temperature (Homan et al., 2014). In following the quiescent thermal evolution of MAXI J0556332 during outburst/quiescence, we demonstrate that a strong heating anomaly is present in its deep
ocean and shallow crust. In fact, the

6 MeV per accreted nucleon of heating is more than can be

supplied by compositionally-driven convection or nuclear reactions. Furthermore, the anomalous
heating required in this source is greater than any other source observed to date.
To fit the cooling light curve, we assume that a M = 1.5 M and R = 11 km neutron star
accreted at the local Eddington rate m
˙ =m
˙ Edd ≡ 8.8 × 104 g cm−2 s−1 for 16 months, matching
the duration of the MAXI J0556-332 outburst (Homan et al., 2014), before cooling began. At
the temperatures observed for MAXI J0556-332, the T eff –T b relation is insensitive to the helium
mass in the envelope. The shallow heat source is uniformly distributed in log y centered on a value
yh = 6.5×1013 g cm−2 (ρ ≈ 1.2×1010 g cm−3 ) and ranging from yh /3 to yh ×3; where the column
depth y is given in Equation 2.3. The strength of shallow heating is assumed to vary proportionally
with the accretion rate.
As can be seen in Figure 4.3, the moderate gravity from a M = 1.5 M and R = 11 km
neutron star best reproduces the quiescent cooling observations of MAXI J0556-332. Late time
∞ ∝ g1/4 /(1 + z) scaling, but also due
observations will change with gravity, not only due to the T eff

59

[ ]

350
300
250
200
150
100
50
100

= .

101

102
[ ]

103

.

104

Figure 4.3: Quiescent lightcurve of the neutron star MAXI J0556-332. The solid black curve
corresponds to a model with M = 1.5 M , R = 11 km, Qshallow = 6.0 MeV, and T core = 108 K; the
dashed black curve is for the same model with T core = 3 × 107 K. The black dotted curves are light
curves with a reheating event ≈ 170 days into quiescence for Qshallow = 6.0 MeV (upper curve)
and Qshallow = 3.0 MeV (lower curve). The blue dashed curve is for a M = 2.1 M , R = 12 km
neutron star fit to the observations by changing the shallow heating depth and strength. the data
∞ ∝ g1/4 /(1 + z)
above the light curve are contamination from residual accretion. Note that T eff
which leads to different observed core temperatures for different gravities.

60

accretion flow
-9

10

envelope
accumulation and compression

106

1

ocean (ion liquid)

100

crust (ion lattice)

1010

anomalous extra heating

depth [m]

mass density [g cm-3]

103

400
e--capture heating
pycnonuclear fusion heating

1014 core

103

1015

105

Figure 4.4: Schematic of the neutron star’s outer layers and the location of the anomalous shallow
heating in MAXI J0556-332.

61

to the thickness of the crust. As shown by Deibel et al. (2015), a fit to the quiescent light curve
requires extra heating of Qshallow = 6.0 MeV per accreted nucleon. The heat must be deposited
over the range y ≈ 2 × 1013 –2 × 1014 g cm−2 in order to fit the break in the light curve near
≈ 20 days, as shown schematically in Figure 4.4. Although the depth of extra heating is well
constrained by the break in the light curve, other neutron star masses and radii require a different
amount of extra heating, as can be seen in Figure 4.6. Similar shallow heating is required in thermal
evolution models of MXB 1659-29 and KS 1731-260 (Brown & Cumming, 2009), each requiring
≈ 1 MeV per accreted nucleon extra heating during outburst near y

2 × 1014 g cm−2 to fit the

quiescent light curve.
Shallow heating raises the temperature of the crust to T b ∼ 109 K by the end of outburst. The
ocean-crust boundary (the crust melting point, Equation 2.4) is pushed significantly deeper as a
result. By the end of outburst the ocean-crust boundary is located at ρt ≈ 2 × 1011 g cm−3 , as
can be seen in Figure 4.5. The crust is in a near maximally hot state, due to the large amount
of shallow heating, with a temperature “ceiling” set by neutrino emission at mass densities ρ
1010 g cm−3 . Because T b > ΘD throughout the outer crust, the thermal conductivity is set by
electron-ion scattering and the thermal time from Equation 4.8 applies.
Note that MAXI J0556-332 requires more shallow heating than other sources; the ≈ 4–10 MeV
per accreted nucleon of shallow heating is larger than required in KS 1731-260 and MXB 1659-29,
each requiring ≈ 1 MeV per accreted nucleon (Brown & Cumming, 2009). This hints at an energy
source much larger than the ≈ 0.2 MeV per accreted nucleon supplied by compositionally driven
convection in the ocean (Medin & Cumming, 2011, 2014, 2015) or the ≈ 2 MeV per accreted
nucleon additional deep crustal heating possible given the uncertainties on the nuclear symmetry
energy (Steiner, 2012). The Keplerian energy of the accretion flow is ∼ 80 MeV per accreted
nucleon (at the inner-most stable circular orbit) and may plausibly provide the shallow heating. As
62

=

[ ]

=

109

108
109

1010

[

]

1011

=

1012

Figure 4.5: Quiescent temperature evolution of the neutron star MAXI J0556-332. Solid black
curves indicate the evolution of the crust temperature during quiescence for the M = 1.5 M and
R = 11 km model, shown in Figure 4.3. The red dotted curve is the melting line of the crust (Γ =
175) for the crust composition in Haensel & Zdunik (1990), the black dotted curve is the transition
from an electron-dominated heat capacity to an ion-dominated heat capacity (CVe = CVion ), and the
blue dotted curve is where the local neutrino cooling time is equal to the thermal diffusion time
(τν = τtherm ). The gray dashed curve shows the lattice Debye temperature ΘD ; when T
ΘD
electron-impurity scattering influences the thermal conductivity.

63

suggested by Inogamov & Sunyaev (1999, 2010), gravitational modes excited in a differentially
rotating envelope may dissipate energy deeper in the star. The mode energies are of the order
required and the dissipation of these modes in the shallow crust is worthy of future study with
realistic ocean and crust models.
The high accretion rate during outburst, when combined with the large amount of shallow
heating, brings the crust into a regime of stable helium burning for helium layers at y ≈ 2 ×
108 g cm−2 (Bildsten & Brown, 1997; Zamfir et al., 2014). During outburst, the crust also enters a
regime of stable carbon burning for carbon layers at y

1010 g cm−2 (Cumming & Bildsten, 2001).

For this reason, an appreciable layer of carbon can not accumulate at the superburst ignition depth
around y ∼ 1012 g cm−2 . There have been no Type I X-ray bursts or superbursts observed from
MAXI to date, consistent with stable burning of helium and carbon. We predict that MAXI J0556332 is unlikely to have either type-I X-ray bursts or superbursts if strong shallow heating occurs
during subsequent accretion outbursts.

4.2.2

The Reflare in MAXI J0556-332

The luminosity of MAXI J0556-332 increased ≈ 170 days into quiescence to a level similar to that
observed at the end of the outburst, and the reflare lasted for ≈ 60 days (Homan et al., 2014). We
˙ Edd for 60 days after a 170 day initial cooling phase to
run a model that accretes at m
˙ ≈ 0.5 m
model the reflare event. We run the model for two values of the shallow heating source, Qshallow =
6.0 MeV and Qshallow = 3.0 MeV, which can be seen in Figure 4.3.
The light curves that include the reflare overshoot the observations hundreds of days into quiescence and the light curve deviation lasts ≈ 500 days before returning to the cooling behavior
seen prior to the reflare. Clearly, extra heating must not be operating during the reflare, or it must
be significantly weaker than during the main accretion outburst. This result may imply that the
64

1 1
.5 2.0 3.5

1.2 1.5 1.8 2.1

log10 Pressure [cgs]
28 28 28 29
.0 .4 .8 .2

M [M ]

log10 Pressure [cgs]

M [M ]

.0
13
.5

12

.5
10

2.1

1.8

1.5

.4
28
.8
29
.2
1.2

28

.0

28

20

16

8
12

4

10

R [km]

Qshallow [MeV]

R [km]

Figure 4.6: Markov chain Monte Carlo fits to the quiescent light curve of MAXI J0556-332 using
the crust relaxation√code crustcool. The contours show the isodensity surfaces of the likelihood L,
corresponding to −2lnL = 0.5, 1, 1.5, 2, for the neutron star mass M, radius R, pressure at the
shallow heating depth Ph , and the shallow heating strength Qshallow .

65

extra heating mechanism is not proportional to the mass accretion rate, as is assumed in our thermal relaxation model. The lack of extra heating during the reflare is similar to the lack of extra
heating in the ≈ 2 month outburst seen in Swift J174805.3-244637 (Degenaar et al., 2015). Extra
heating is needed, however, after a similar ≈ 2.5 month outburst in IGR J17480-2446 (Degenaar
& Wijnands, 2011; Degenaar et al., 2011, 2013). These sources may have fundamental differences
in their structure that make extra heating active in one and not the others, or perhaps the growth
timescale of the extra heating mechanism is ∼ 2 months in all sources.

4.3

Observational Tests for Ocean g-modes

The Z sources are neutron stars the accrete intermittingly near the Eddington rate (Hasinger & van
der Klis, 1989) as inferred from their X-ray luminosities between LX ∼ 0.5–1.0 LX,Edd , see Lewin
& van der Klis (2006) for a review. These sources trace out a distinctive “Z” shape in the X-ray
color-color diagram: the top of the Z is the horizontal branch, the diagonal is the normal branch,
and the bottom of the Z is the flaring branch. The spectral state of the source (i.e., its location on the
Z) was traditionally thought to be correlated with changes in the mass accretion rate (Priedhorsky
et al., 1986; Hasinger et al., 1990; Homan et al., 2007), but more likely cause secular evolution of
a source between a Z-source state and a lower luminosity atoll state (LX ∼ 0.1–0.5 LX,Edd ) (Lin
et al., 2009). Instead, accretion instabilities at a constant accretion rate are likely the determining
factor for changes in the spectral state (Lin et al., 2009; Homan et al., 2010; Fridriksson et al.,
2015).
Quasi-periodic oscillations (hereafter QPOs) have been observed in Z-source power spectra
during their evolution along the Z-track and each branch of the Z-track has a characteristic set of
QPO frequencies, see Lewin & van der Klis (2006) for a review. These QPOs typically appear at

66

low-frequencies: the horizontal branch oscillations (hereafter HBOs) are QPOs with frequencies
between ≈ 10−30 Hz. The normal branch oscillations (hereafter NBOs) are QPOs with frequencies
∼ 6 Hz. The NBOs increase as the source moves toward the flaring branch and the NBOs eventually
blend into flaring branch oscillations (hereafter FBOs) around ≈ 10−50 Hz. For instance, continous
blending have been observed in Sco X-1 (van der Klis, 1989; Casella et al., 2006).
The HBOs are thought to be of a geometric origin, resulting from Lense-Thirring precession of
the inner-most region of the accretion disk (Lense & Thirring, 1918; Bardeen & Petterson, 1975).
Furthermore, the HBOs are consistent with similar frequency QPOs from accreting black holes in
the same spectral state. The NBOs may originate from oscillations on the neutron star or accretion
disk. For example, the neutron star ocean may host a variety of non-radial oscillations; such as,
thermal g−modes (Bildsten & Cutler, 1995; Bildsten et al., 1996; Bildsten & Cumming, 1998).
The FBOs are consistent with shallow surface waves in the ocean, known as crustal interface
modes (Piro & Bildsten, 2005). Although the FBOs typically have large frequencies

100 Hz, the

deformability of the crust may pin these modes deep in the crust and result in lower frequencies
20 Hz (Piro & Bildsten, 2005).
The cores of accreting neutron stars may also contain oscillations, although these oscillations
occur at much high frequencies. Specifically, rotating relativistic stars are unstable to the production of Rossby-modes (hereafter r-modes) in their cores, that are excited by gravitational radiation
(Andersson, 1998) — analogous to hurricane-like r-modes in a terrestrial ocean. Recently, limits
were placed on neutron star core r-mode amplitudes from spin-down rates and thermal emission
from quiescent low-mass X-ray binaries (Mahmoodifar & Strohmayer, 2013). Interestingly, the
core r-modes may propagate into the neutron star’s ocean during an accretion outburst, where
their amplitudes increase (Lee, 2014). A candidate r-mode was observed in XTE J1751-305 as a
≈ 250 Hz oscillation in the X-ray light curve (Strohmayer & Mahmoodifar, 2014). This frequency
67

accretion flow

10-9 envelope
accumulation and compression

1

ocean (ion liquid)
g-modes

106

100

crust (ion lattice)

1010

anomalous extra heating

depth [m]

mass density [g cm-3]

103

400
e--capture heating
pycnonuclear fusion heating

103

1014 core
r-modes

105

1015

Figure 4.7: Schematic of the neutron star’s outer layers and the location of g-modes.
is approximately 2/3 of the star’s rotation frequency (ν ≈ 435 Hz), which is consistent with a core
r-mode (Andersson et al., 2014).
In the following section, we outline an observational test for the origin of NBOs using the
thermal evolution of the neutron star ocean during an accretion outburst, following the approach
of Deibel (2016). In particular, the larger ocean temperatures during an accretion outburst will
raise the ocean’s fundamental g-mode frequencies and therefore lead to larger NBO frequencies
if ocean g-modes are indeed the origin of the NBOs. Moreover, the hottest neutron star transients

68

that require shallow heating in their oceans have the largest g-mode frequencies during outburst
and thereby provide excellent observational tests for the presence of ocean g-modes. In particular,
we will examine the ocean’s thermal evolution in three neutron star transient Z-sources, all with
observed accretion outbursts: MXB 1659-29, XTE J1701-462, and MAXI J0556-332.

4.3.1

Neutron Star Ocean Thermal Evolution

During an accretion outburst, the ocean temperature rises as accretion-driven nuclear reactions in
the crust supply a heat flux into the ocean (Brown, 2000). Although a majority of crustal heating
(≈ 90%) is conducted toward the core, the heat flux entering the ocean is sufficient to raise the
temperature to T b

108 K during outburst. The increase in temperature will be especially pro-

nounced in neutron stars that have active shallow heating during their accretion outbursts, such as
KS 1731-260 and MXB 1659-29 which require ≈ 1 MeV per accreted nucleon (Brown & Cumming, 2009), and MAXI J0556-332 which requires ≈ 6–16 MeV per accreted nucleon (Deibel
et al., 2015).
As can be seen in Figure 2.2, the ocean crystallization point (given by Equation 2.5) shifts
to higher mass densities as the crust temperature increases. Initially the neutron star ocean is in
thermal equilibrium with the neutron star core near T b ≈ T core = 3 × 107 K and the ocean-crust
transition density is ρ ≈ 2.2 × 106 g cm−3 . During an accretion outburst, the ocean temperature
increases, and the ocean-crust transition density increases, as unavoidable accretion-driven nuclear
reactions heat the ocean out of thermal equilibrium with the core. We now examine the ocean’s
thermal evolution in three neutron star transient Z-sources: MXB 1659-29, XTE J1701-462, and
MAXI J0556-332.
We run a model of the MAXI J0556-332 ≈ 460 day accretion outburst (Matsumura et al., 2011;
Homan et al., 2011; Sugizaki et al., 2013; Homan et al., 2014), with the neutron star parameters fit
69

from quiescent observations: M = 1.5 M , R = 11 km, Qshallow = 6 MeV per accreted nucleon,
and Qimp = 1. During outburst, the ocean-crust transition (Equation 2.5) moves from ρ ≈ 2.2 ×
106 g cm−3 at the beginning of outburst to ≈ 2×1011 g cm−3 by the end of outburst when the ocean
temperature is T b ≈ 7.5 × 108 K. The ocean properties as a function of time during the accretion
outburst can be seen in Figure 4.8. The ocean has T b

7 × 108 K for ≈ 500 days into quiescence

and takes ≈ 2900 days to reestablish thermal equilibrium with the core.
We also run thermal evolution models for the observed accretion outbursts in MXB 1659-29
and XTE J1701-462. For MXB 1659-29 we run a model of its previous 2.5 year outburst (Wijnands
˙ Edd . Brown &
˙ ≈ 0.1 M
et al., 2003, 2004) at the accretion rate inferred for this source near M
Cumming (2009) fit the quiescent light curve of this source and found the best-fit parameters:
M = 1.5 M , Qimp ≈ 4, and Qshallow ≈ 1 MeV per accreted nucleon. The properties of the
ocean in MXB 1659-29 during outburst and quiescence can be seen in Figure 4.8. The ocean
has T b

108 K for ≈ 700 days into quiescence and takes ≈ 1400 days to reestablish thermal

equilibrium with the core.
For XTE J1701-462, we run a thermal evolution model for its ≈ 2 year outburst (Fridriksson
et al., 2010, 2011). We use the neutron star parameters fit from quiescent light curve fits: M =
1.6 M , R = 11.6 km, and Qimp ≈ 7. We also add Qshallow = 0.17 MeV per accreted nucleon
from a recent study of the quiescent light curve (Turlione et al., 2015). The properties of the
ocean in XTE J1701-462 during outburst and quiescence can be seen in Figure 4.8. The ocean
has T b

108 K for ≈ 1700 days into quiescence and takes ≈ 3100 days to reestablish thermal

equilibrium with the core.

70

16
14
12
10
8
6
4
2
0
11

,

[ ]

(a)

(b)

[

]

1010
10 9
108
107
106
109
10

[ ]

(c)

108
1071000 500 0

500 1000 1500 2000 2500 3000 3500

[

]

Figure 4.8: Thermal evolution as a function of time during outburst/quiescence for MAXI J0556332 (solid curves), MXB 1659-29 (dashed curves), and XTE J1701-462 (dot-dashed curves). Panel
(a): Fundamental n = 1, l = 2 g-mode frequencies in the ocean (Equation 4.10). Panel (b): The
ocean-crust transition density (Equation 2.5). Panel (c): Temperature at the ocean-crust transition.

71

4.3.2

g−mode Spectrum of the Neutron Star Ocean

The g−modes may be excited in the neutron star ocean as angular momentum is transported into
the ocean by a spreading layer of accreted material (Inogamov & Sunyaev, 1999, 2010) — the
same mechanism that excites global acoustic modes in the neutron star envelope during accretion
(Philippov et al., 2016). The g−mode spectrum has been derived for the ocean (McDermott et al.,
1983; Bildsten & Cutler, 1995). For a non-rotating neutron star with a magnetic field B

1011 G,

the g-mode frequencies can be analytically approximated within ≈ 10 % (Bildsten & Cutler, 1995),
l(l + 1)
Tb
56 1/2 10 km
6 3 × 107 K A
R


−2 −1/2

3nπ 2
ρt


× 1 +
ln
,

−3
2
1 g cm

fn,l ≈ 2.8 Hz

(4.10)

where n is the number of nodes in the ocean, l is the angular wavenumber, and R is the radius of
the neutron star. The upper boundary of the ocean is taken at the base of the envelope near a mass
density ρ ≈ 1 g cm−3 , however, the frequency spectrum is largely insensitive to the location of
this upper boundary. Even though Equation 4.10 is derived for an isothermal ocean, the oscillation
spectrum is primarily set by the temperature at the base of the ocean in non-isothermal models
(Bildsten & Cumming, 1998), and Equation 4.10 gives frequencies accurate within ≈ 10 % (Bildsten & Cutler, 1995). Note that this expression is not valid for rapidly rotating neutron stars with
fspin

300 Hz, where the g-mode frequencies become highly modified (Bildsten et al., 1996).

The thermal time in the ocean (Equation 4.8) determines the time to raise the ocean’s temperature during the accretion outburst. The thermal time in the liquid ocean is shown in Figure 4.2.
Normal-branch oscillations have been observed in several sources, for example Cyg X-2 (Hasinger
& van der Klis, 1989; Wijnands et al., 1997; Dubus et al., 2004), SCO X-1 (Hertz et al., 1992; van

72

10
8
6
4
2
0

,

[ ]

(a)

]

(b)

108

[

107
106

(c)

[ ]

108
107500

.
=.

0

500

[

]

1000

1500

Figure 4.9: Thermal evolution as a function of time during outburst/quiescence for neutron star
˙ = 0.1–1.0 M
˙ Edd . Panel (a): Fundamental n = 1, l = 2 g-mode frequencies
transients with M
in the ocean (Equation 4). Panel (b): The ocean-crust transition density (Equation 2). Panel (c):
Temperature at the ocean-crust transition.
der Klis et al., 1996; Titarchuk et al., 2014), and GX 5-1 (Kuulkers et al., 1994; Jonker et al.,
2002). NBO frequencies are typically between ≈ 5–7 Hz, but larger frequencies are observed; for
example, XTE J1701-462 has NBO frequencies between ≈ 7–9 Hz (Homan et al., 2007; Fridriksson et al., 2010). As we show in this section, the temperature of the neutron star’s ocean during
outburst supports fundamental g-modes of similar frequencies. In particular, we run a model of the
last outburst in XTE J1701-462 and find g-modes consistent with the observed NBO frequencies
in this source.

73

Before accretion begins, the predicted fundamental g-mode is initially ≈ 3 Hz in the cold ocean
(T b ≈ T core ∼ 107 K). During active accretion, the predicted fundamental g-mode frequency increases as unavoidable accretion-driven nuclear heating raises the ocean’s temperature. For exam˙ = 0.1–1.0 M
˙ Edd have predicted fundamental g-modesbetween
ple, sources accreting between M
≈ 3–7 Hz, as can be seen in Figure 4.9. Predicted fundamental g-mode frequencies
quire active shallow heating during outburst to reach ocean temperatures
support these frequencies. Therefore, observing NBOs

7 Hz re-

108 K required to

7 Hz in transients that require shal-

low heating, such as MAXI J0556-332, MXB 1659-29, and XTE J1701-462, provides an observational test to link NBOs and the ocean g-modes. In particular, during steady state accretion
MAXI J0556-332 has predicted fundamental g-mode frequencies between ≈ 8–16 Hz, MXB 165929 and XTE J1701-462 have predicted g-mode frequencies between ≈ 8–10 Hz — frequencies only
possible in the hotter oceans found in sources with shallow heating. The evolution of the predicted
fundamental g-mode in these sources can be seen in Figure 4.8. We predict that MAXI J0556-332
had a fundamental g-mode frequency near ≈ 11 Hz at the time of its recent activity (Negoro et al.,
2016; Russell & Lewis, 2016; Jin & Kong, 2016). If shallow heating is active during the current
outburst in MXB 1659-29 (Negoro et al., 2015), we predict this source will have a fundamental
g-mode frequency near ≈ 9 Hz at the time of the most recent observation (Bahramian et al., 2016)
after being in outburst for ≈ 160 days.
NBO frequencies observed in the last outburst from XTE J1701-462 (Fridriksson et al., 2010)
are consistent with ocean g-modes. NBOs were observed within the first ≈ 10 weeks near ≈
7 Hz (Homan et al., 2007), and our model predicts fundamental g-mode frequencies near ≈ 7 Hz
after 10 weeks of active accretion. Furthermore, once the ocean reaches steady state, our model
of the XTE J1701-462 outburst contains predicted ocean g-modes between ≈ 7–10 Hz, which is
consistent with the observed ≈ 7–9 Hz NBOs observed in this source during outburst (Homan
74

et al., 2010).
Although current instrumentation may not allow observations of NBOs in MXB 1659-29 and
MAXI J0556-332 during quiescence, renewed accretion outbursts in both objects may allow observations of NBOs. We predict that observed NBOs during these sources current accretion outbursts
should be

5 Hz in MXB 1659-29 and

10 Hz in MAXI J0556-332. The renewed activity in

MAXI J0556-332(Negoro et al., 2016; Russell & Lewis, 2016) perhaps holds the best prospects
for catching large NBO frequencies because the ocean in MAXI J0556-332 had not cooled to
the core temperature before the most recent accretion outburst, and the predicted ocean g-mode
frequency is already

7 Hz.

Note that the above observational tests for the presence of ocean g-modes does not require a
complete picture of how X-ray emission is modulated by ocean oscillations, which is outside the
scope of this work. Uncertainties remain in the nature of the accretion flow near the neutron star
surface and how the accretion flow may interact with the neutron star’s outer layers. For example,
a spreading boundary layer of accreted material may extend into the ocean (Inogamov & Sunyaev,
1999, 2010), and may transport angular momentum therein. A coupling of the boundary layer and
the ocean in this way may only modulate X-ray emission near the equator, however, where modes
may even interfere with accretion disk emission (Bildsten et al., 1996). If observations support
a g-mode origin for the NBOs, this will motivate future work to model the ocean’s coupling to a
spreading boundary layer.
If NBOs are indeed ocean g-modes, observed NBO frequencies

7 Hz can then be used as

a diagnostic of the shallow heating strength. For example, as shown in this work, ocean g-mode
frequencies near ≈ 10 Hz (≈ 16 Hz) are found in sources in steady state with Qshallow ≈ 1 MeV
(Qshallow ≈ 6 MeV). It is worth noting, however, that during outburst the ocean in MAXI J0556332 approaches a maximum temperature set by neutrino emission around T b ≈ 2 × 109 K at the
75

depth of the shallow heat source (Deibel et al., 2015). Because the ocean temperature is near
maximum, a g-mode frequency of ≈ 16 Hz is likely the largest n = 1 l = 2 g-mode that can be
supported in Z-source oceans. Therefore, NBO frequencies of ≈ 16 Hz only give the lower limit
Qshallow

6 MeV for the shallow heating strength.

The depth of the ocean-crust transition during outburst is relevant to the study of the shallow
heating mechanism. For example, the ocean-crust interface in MAXI J0556-332 moves to densities
between ρt ∼ 109 –1011 g cm−3 during outburst. This density range is characteristic of the location
of the shallow heating inferred from quiescent light curves of quasi-persistent transients (Brown &
Cumming, 2009), which is needed at mass densities near ρshallow

3 × 1010 g cm−3 . For example,

shallow heating is required between ρshallow ∼ 5 × 109 –3 × 1010 g cm−3 in MAXI J0556-332,
as determined from its quiescent light curve (Deibel et al., 2015). This suggests that the shallow
heating mechanism is connected to the ocean-crust phase transition, and work in this direction is
ongoing.
Often the observed ≈ 5–7 Hz NBOs blend into ≈ 10 − 20 Hz flaring branch oscillations (hereafter FBOs); for example, a continuous blending was observed in SCO X-1 (Casella et al., 2006).
The timescale of observed NBO-FBO blending disfavors a g-mode origin for the FBOs. The
thermal time at the depth of the ocean-crust transition (Figure 4.10) determines the timescale for
variations in the g-mode frequencies. For g-mode frequencies
timescales of hours; variations in frequencies

7 Hz, variations may occur on

7 Hz occur on timescales of days — much longer

than the timescale of the observed blending. The NBO-FBO blending phenomena, however, is
consistent with a superposition of ocean g-modes and ocean interface waves. The interface modes
were determined to have frequencies ∼ 200 Hz when calculated with a rigid crust boundary condition (McDermott et al., 1988), but may have lower frequencies near ∼ 20 Hz when the flexibility of
the crust is taken into account (Piro & Bildsten, 2005). It is unclear why oscillation energy might
76

]
[

103
102
101
100
10 1
10 2
10 3
106

107

108

[

109

]

1010

1011

Figure 4.10: Characteristic timescales in the neutron star ocean. The upper curves are the thermal
time for an equilibrium composition (solid curve) and accreted composition (dotted curve). The
lower curves are for the timescale for changes in ρt for an equilibrium composition (solid curve)
and an accreted composition (dotted curve).
transition between the different modes, but the coupling of the g-modes to the interface modes is
worthy of future study.
Quasi-periodic oscillations between ≈ 5–7 Hz in the atoll sources (LX ∼ 0.1–0.5 LX,Edd ) are
also consistent with ocean g-modes. For example, the ≈ 7 Hz oscillation observed in 4U 182030 (Wijnands et al., 1999; Belloni et al., 2004) is consistent with a n = 1 l = 2 g-mode in a
lighter ocean. In a lighter ocean, such as those considered in Bildsten & Cutler (1995), g-mode
frequencies are larger by a factor of ∼ (56/16)1/2 compared to those studied here. This may
explain the observed oscillations in 4U 1820-30, where accretion from the helium companion star
(Wijnands et al., 1999) would result in a different ocean composition than the one considered here.

77

4.4

Late Time Cooling and Nuclear Pasta

As discussed in Section 4.2, as a neutron star cools in quiescence the X-ray light curve reveals
layers with increasingly longer thermal times. For the cooling front to reach the deepest layers
of the inner crust, close to the crust-core transition, the neutron star must cool for

1000 days.

Though many quiescent neutron stars experience a new accretion outburst before cooling entirely,
one cooling neutron star transient MXB 1659-29 cooled for

2500 days and allows a unique view

of the inner crust’s thermal properties.
Following the outburst and return to quiescence of MXB 1659-29 (Wijnands et al., 2003, 2004;
Cackett et al., 2008), crust thermal relaxation models of MXB 1659-29 (Brown & Cumming, 2009)
initially indicated that the outer layers in this source had returned to thermal equilibrium with the
core, as suggested by a plateau in the quiescent luminosity after ≈ 1000 days in quiescence. A
recent observation of MXB 1659-29 (Cackett et al., 2013), however, demonstrates that cooling
continued after ≈ 2500 days into quiescence, suggesting that thermal equilibrium with the core
had not yet been established. Cooling at late times near ≈ 2500 days corresponds to cooling layers
deep in the neutron star’s inner crust near a mass density ρ ∼ 1013 g cm−3 , densities predicted for
the formation of nuclear pasta.
In the high density environment of the inner neutron star crust, nuclei are distorted into various
shapes, called nuclear pasta (Ravenhall et al., 1983; Hashimoto et al., 1984). Nuclear pasta has
been studied using quantum molecular dynamics simulations (Maruyama et al., 1998; Watanabe
et al., 2003) and semi-classical molecular dynamics simulations (Horowitz et al., 2004a; Horowitz
& Berry, 2008; Schneider et al., 2013), but the thermal properties of nuclear pasta remain uncertain
(Horowitz et al., 2015). Because the quiescent light curve reveals successively deeper layers with
time, the cooling light curve thousands of days into quiescence is affected by the thermal properties

78

[ ]

140
120
100
80
60
40
101

102

[

103
]

104

Figure 4.11: Crust cooling models of the late time cooling of MXB 1659-29. Cooling models for
MXB 1659-29. The solid gray curve is a model that uses Qimp = 2.5 throughout the entire crust
and T core = 4 × 107 K. The solid blue curve is a model with Qimp = 20 for ρ > 8 × 1013 g cm−3 ,
Qimp = 1 for ρ < 8 × 1013 g cm−3 , T core = 3.25 × 107 K, and using the G08 pairing gap. The
dashed red curve uses the same Qimp as the solid blue curve, but with the S03 pairing gap. The
dotted blue curve is a model with the G08 pairing gap and Qimp = 1 throughout the crust, but
without a low thermal conductivity pasta layer.

79

of the nuclear pasta in the deep inner crust. Horowitz et al. (2015) show that a low thermal conductivity layer may explain the late time cooling of MXB 1659-29, and that spiral defects may lower
the pasta’s thermal conductivity. This prediction is consistent with the low electrical conductivity
of pasta that may explain the cutoff in the spin period distribution of pulsars at P ∼ 10 s (Pons
et al., 2013).
To investigate the impact of a pasta layer on the late time quiescent light curve of MXB 165929, we run a model of the MXB 1659-29 outburst using the best fit parameters from Brown & Cum˙ Edd
ming (2009): an accretion outburst lasting 2.5 years at a local mass accretion rate m
˙ ≈ 0.1 m
˙ Edd ≈ 8.8 × 104 g cm−2 s−1 is the local Eddington mass accretion rate, a neutron star mass
where m
M = 1.6 M , a neutron star radius R = 11.2 km, and a core temperature T core = 2.6 × 107 K.
The impurity parameter is defined for two separate regions: 1) for the crust Qimp = 3.5 following
Brown & Cumming (2009) and only impacts the thermal conductivity in the inner crust where
T < T p , and 2) a pasta layer with a second impurity parameter Qpasta defined at mass densities
ρ > 8 × 1013 g cm−3 and extending to the crust-core transition.
The light curve from our standard model without pasta can be seen in Figure 4.11. In this
model the crust reaches thermal equilibrium with the core by ≈ 1000 days into quiescence, and
so predicts a constant temperature at later times. We also show two models with an impure inner
crust with Qpasta = 25. One of these corresponds to the model in Horowitz et al. (2015) and shows
a long decline in temperature even at times as late as ≈ 2000 days. The other impure crust models,
however do not show a decline at late times, but instead level off. The difference in the light curves
is in the assumed T c (ρ) profiles, which are shown in Figure 4.12. Significant late time cooling only
occurs if there is a normal layer of neutrons at the base of the crust, giving a large heat capacity
there. The pasta layer maintains a temperature difference of ∆T ≈ 3 × 107 K between the inner
crust and core during the outburst. As a consequence, normal neutrons with a long thermal time
80

1010 (a)

pasta

[ ]

109
108

]

107
106

(b)

]

[

105
104
103

(c)

1020

[

1019
1018 11
10

1012

1013
[
]

1014

Figure 4.12: Thermal transport in the inner crust of MXB 1659-29 at the start of quiescence. The
gray vertical lines indicate the neutron drip density and the transition to nuclear pasta. Panel (a):
The temperature profile (solid curve) corresponding to the cooling model in Figure 4.11 with a
Qimp = 20 pasta layer and the G08 pairing gap. The dashed curves show two choices for T c (ρ);
the blue dashed curve corresponds to G08 and the red dotted curve is S03. Panel (b): The heat
capacity profiles for the same models as Figure 4.11. Solid black curve: Qimp = 3.7 throughout
the inner crust. Dashed blue curve: Qimp = 20 for ρ > 8 × 1013 g cm−3 and Qimp = 1 for
ρ < 8 × 1013 g cm−3 using the G08 pairing gap that closes in the crust. Dotted red curve: same as
dashed curve, but with a different choice for T c (ρ) from the S03 pairing gap that closes in the core.
Panel (c): Thermal conductivity profiles for the same models.

81

appear at the base of the crust that cause late time cooling if the neutron singlet pairing gap closes
in the crust. Without normal neutrons at the base of the crust, as is the case if the neutron singlet
pairing gap closes in the core, the crust reaches thermal equilibrium with the core after ≈ 3000 d
and late time cooling is removed.
Some analytic estimates are useful to understand why the late time cooling occurs, and the
crucial role of the normal neutron layer. First, we consider the temperature contrast ∆T between the
inner crust and the core that develops during the accretion outburst (see panel (a) of Figure 4.12).
This is set by the value at which the heat flux through the pasta layer balances the nuclear heating
in the crust (mostly located at shallower densities near the neutron drip region). The heating
rate is nuc = mE
˙ nuc where Enuc ≈ 2 MeV per accreted nucleon. The equivalent heat flux is
Fin ≈ 2 × 1022 erg cm−2 s−1 for an accretion rate of m
˙ = 0.1 m
˙ Edd .
The heat flux through the pasta layer is F ≈ K∆T/H, where K is the thermal conductivity and
H the pressure scale height. Neutrons set the pressure in the inner crust, so that1 H = P/ρg ≈
7 × 104 cm (ρ2/3
Y 5/3 /g14 ) where Yn is the neutron fraction, ρ14 is the mass density in units of
14 n
1014 g cm−3 , and the surface gravity of the neutron star is g = (GM/R2 )(1−2GM/Rc2 )−1/2 in units
of 1014 cm s−2 . As can be seen in Figure 4.13, the thermal conductivity is primarily set by electronimpurity scattering (Equation 2.33) because neutron-impurity scattering only becomes comparable
to electron scattering very close to the crust-core transition. The resulting thermal conductivity
from electron-impurity scattering is given in Equation 2.35. Therefore, the temperature difference
between inner crust and core is
 1/3
 ρ
Qimp ΛeQ
Yn5/3
∆T ≈ 3 × 107 K  14
g14 T 8 (Ye /0.05)1/3
Z

m
˙
0.1m
˙ Edd



 ,


(4.11)

1 In the inner crust, Γ ≡ (∂lnP/∂lnρ) varies with density: at first Γ decreases below 4/3 (the value for degenerate
s
1
1
relativistic electrons) and then it increases for ρ 1013 g cm−3 and approaches Γ1 2 at roughly nuclear density.
For definiteness in computing H, we set Γ1 = 5/3.

82

( )

1021

20
19
18

,

0.7 0.9 1.1 1.3 1.5

[

]

1020
1019
1018
1017

1012

1013

[

pasta

]

1014

Figure 4.13: Scattering frequencies for electrons and neutrons in the inner crust at the beginning of
quiescence for the model with Qimp = 20 at ρ > 8×1013 g cm−3 , Qimp = 1 at ρ < 8×1013 g cm−3 ,
and the pairing gap that closes in the crust (Gandolfi et al., 2008). Subplot: Thermal conductivity K
from electron scattering (dotted red curve), neutron scattering (dashed blue curve), and from both
electrons and neutrons (solid black curve). The mass density ρ is given in units of 1014 g cm−3 .
The region containing nuclear pasta is to the right of the vertical black dotted line.
which is in reasonable agreement with the temperature jumps seen in Figure 4.12, panel (a), between ρ ≈ 8 × 1013 g cm−3 and ρ ≈ 1.5 × 1014 g cm−3 .
Page & Reddy (2012) pointed out that differences in T c (ρ) and the resulting presence or absence
of a layer of normal neutrons at the base of the crust could affect the cooling curves at late times
≈ 1000 days into cooling. We find a much larger effect and on a longer timescale here because the
low thermal conductivity of the nuclear pasta layer keeps the inner crust much hotter during the
outburst. During quiescence, the base of the crust remains at a higher temperature than the core for
≈ 5000 days (see Equation 4.11). The temperature difference between the crust and core results in
a slow decline of the quiescent light curve after

1000 days, as can be seen in Figure 4.11.

Late time cooling in MXB 1659-29 requires that the 1 S0 neutron singlet pairing gap close
83

in the crust. As a result, superfluid neutrons are confined to the inner crust shallower than the
pasta layer at ρ

8 × 1013 g cm−3 where T

T c . By contrast, a recent study of pulsar glitches

suggests that the neutron superfluid extends from the crust into the core continuously (Andersson
et al., 2012). Recent calculations of the neutron effective mass in a non-accreted (Qimp = 0) crust
suggest that mn

mn at the base of the crust (Chamel, 2005, 2012). In this case, a larger fraction

of free neutrons are entrained in the inner crust and the neutron superfluid must then extend into
the core to supply adequate inertia for pulsar glitches (Andersson et al., 2012; Ho et al., 2015). We
note, however, that the above calculation for the neutron effective mass is likely inappropriate for
the impure crust compositions found in the accreting transients studied here. Therefore, we here
assume mn ≈ mn as found in Brown (2013) in the absence of neutron effective mass calculations
in an accreted crust.

84

Chapter 5
Urca Cooling Nuclei Pairs
Neutrino emission through cycles of e− −capture and β− −decay , or “Urca” cycles, was first proposed as a cooling process relevant for stellar collapse (Gamow & Schoenberg, 1940). An Urca
cycle occurs when an e− −capture daughter nucleus undergoes β− −decay on a timescale much
shorter than other possible nuclear reactions, producing the cycle

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

(5.1)

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

(5.2)

where Z is the proton number and A is the mass number of the e− −capture parent nucleus. The
e− −capture parent and the e− −capture daughter thereby form an Urca pair. Neutrinos liberate a
majority ( 2/3) of the e− −capture threshold energy QEC (Gamow & Schoenberg, 1940; Haensel
& Zdunik, 1990) which escape from the star entirely due to the large mean free path of the neutrino.
The strong T 5 dependence of the Urca cycling neutrino emission means that Urca cooling becomes
important in high-temperature systems. For example, Urca cooling has been studied in the context
of white dwarfs interiors at temperatures near ∼ 109 K (Tsuruta & Cameron, 1970). As we will
show, Urca cycles occur in the neutron star interior and become an important cooling mechanism
in accreting neutron stars that reach temperatures T

85

108 K in their oceans and crusts.

5.1

Urca Pair Formation

An accretion outburst in a neutron star transient deposits hydrogen- and helium-rich material into
the neutron star’s envelope. As the envelope temperature rises, accreted hydrogen-rich and heliumrich material in the envelope may eventually ignite unstably, powering a type-I X-ray burst. Burning in type-I X-ray bursts proceeds via the rapid proton-capture process (Fujimoto et al., 1981; Wallace & Woosley, 1981) which produces an array of neutron-rich nuclei and a significant amount of
carbon (Schatz et al., 1999); for example, unstable nuclear burning typically produces nuclei with
masses A ∼ 60–100 in Type I X-ray bursts (Schatz et al., 1998, 2001; Brown et al., 2002; Schatz
et al., 2003; Woosley et al., 2004; Cyburt et al., 2010). The accumulation of more accreted material
in the neutron star’s envelope compresses the carbon-rich layer deeper in the star. Eventually, the
ambient density and temperature are sufficient to ignite the 12 C + 12 C fusion reaction, triggering
a superburst in the envelope (Schatz et al., 2003), and producing ashes in the same mass regime
(Keek & Heger, 2011; Keek et al., 2012). Note that stable nuclear burning in hot neutron star
envelopes

5 × 108 K also produces nuclei of similar masses (Schatz et al., 1999).

The accumulation of material in the envelope by subsequent accretion outbursts compresses the
ashes of nuclear burning deeper into the neutron star. Nuclear reaction network calculations that
follow the compressed ashes of nuclear burning found that ashes contain Urca pairs that become
active in the neutron star crust (Schatz et al., 2014). Urca pairs had previously been considered as
a coolant in neutron stars (Bahcall & Wolf, 1965), but were not found in nuclear reaction network
calculations of the accreted crust. Nuclear reaction network calculations of the accreted crust were
previously done in the zero temperature approximation, which is typically a valid assumption for
e− -captures in the degenerate electron gas where µe

kB T . At zero temperature, e− -captures

occur when µe > |QEC | and β− -decays are blocked because electron phase space is restricted.

86

Urca reaction shells will appear if the network uses a finite temperature instead of the zero
temperature approximation (Schatz et al., 2014). At finite temperature, e− -captures occur when
|QEC | − kB T

µe

|QEC | + kB T and β− -decays are allowed. These reactions occur in an Urca

“shell” with a thermal width µe ± kB T , or (∆R)shell ≈ Ye kB T /mu g (Cooper et al., 2009; Schatz
et al., 2014). A schematic of an Urca reaction layer in a neutron star is shown in Figure 5.1.
Following the approach of Tsuruta & Cameron (1970), in the relativistic limit (µe ≈ EF
me c2 ) the specific neutrino emissivities from e− -captures (+) and β− -decays (−) in the degenerate
electron ocean are
±
ν

≈ me c2

ln 2
ft

F ± n± I ± (EF , T ) ,

(5.3)

where n± is the number density of nuclei, α ≈ 1/137 is the fine structure constant, and the Coulomb
factor can be expressed as F ± ≈ 2παZ/|1 − exp (∓2παZ)|. The electron phase space integrals are

I + (EF , T ) =

I − (EF , T ) =

∞

W W 2 − 1(W − Wm )3 S dW ,

(5.4)

W W 2 − 1(Wm − W)3 (1 − S )dW ,

(5.5)

Wm
Wm

1

where W ≡ Ee /me c2 is the electron energy, Wm ≡ |QEC |/me c2 is the electron capture threshold
energy, S = (1 + exp [(Ee − EF )/kB T ])−1 is the statistical factor, and EF ≈ 3.7 MeV (ρ9 Ye /0.4)1/3
is the electron Fermi energy, and ρ9 ≡ ρ/(109 g cm−3 ). The total neutrino luminosity in the Urca
shell is found by integrating the sum of the specific neutrino emissivities from e− -captures and
β− -decays over the reaction shell,

Lν ≈ 4πR2
shell

( ν+ + ν− ) dz ,

87

(5.6)

Figure 5.1: Schematic of an Urca reaction shell. The red colored region is composed of the parent
nucleus, which electron captures into the daughter nucleus represented by the blue colored region.
The Urca reaction layer has a thermal width QEC ± kB T , where µe ≈ QEC is the center of the Urca
shell.
88

where R is the radius of the neutron star. The Urca cycle occurs over a thin shell defined by the
temperature (∆R)shell ≈ Ye kB T /mu g (Cooper et al., 2009; Schatz et al., 2014) where (∆R)shell
R, Ye ≈ Z/A is the electron fraction, and the surface gravity of the neutron star is g = GM/R2 .
The integral in Equation 5.6 can be solved analytically by transforming the integration variable
dz ≈ (dP/dµe )(dµe /ρg). The approximate neutrino luminosity in the Urca shell is

Lν ≈ L34 × 1034 ergs s−1 XT 95

g14 −1 2
R10 ,
2

(5.7)

where the parameter L34 is a function of nuclear properties only,
 6 
 10 s  56

L34 = 2.15 
ft
A

QEC 5 F ∗
,
4 MeV
0.5

(5.8)

and we define the parameters: T 9 ≡ T/(109 K), g14 ≡ g/(1014 cm s−2 ), R10 ≡ R/(10 km), F ∗ ≡
F + F − /( F + + F − ), and X ≡ Amu n/ρ is the mass fraction of the parent nucleus in the
composition.

5.2

Crust Urca Pairs

A nuclear reaction network calculation that follows compressed X-ray burst and superburst ashes
found Urca pairs form in the neutron star crust (Schatz et al., 2014). Schatz et al. (2014) found
that the presence of Urca pairs is robust for X-ray burst and superburst ashes. Furthermore, at
large enough abundances neutrino cooling from Urca pairs in the crust balances accretion-driven
heating at crust temperatures

2 × 108 K. The Urca pairs in the crust with the strongest neutrino

luminosities are shown in Table 5.1. Because Urca cooling may balance accretion-driven heating,
Urca pairs may leave observational signatures in the cooling light curves of neutron star transients
89

Table 5.1: Urca pairs in the neutron star crust (Schatz et al., 2014).
Urca pair
parent daughter
29 Mg

29 Na

55 Ti

55 Sc, 55 Ca

31 Al

31 Mg

33 Al

33 Mg

56 Ti

56 Sc

57 Cr

57 V

57 V

57 Ti, 57 Sc

63 Cr

63 V

105 Zr

105 Y

59 Mn

59 Cr

103 Sr

103 Rb

96 Kr

96 Br

65 Fe

65 Mn

65 Mn

65 Cr

ρ
[1010 g cm−3 ]
4.79
3.73
3.39
5.19
5.57
1.22
2.56
6.82
3.12
9.45
5.30
6.40
2.34
3.55

µe
[MeV]
13.3
12.1
11.8
13.4
13.8
8.3
10.7
14.7
11.2
7.6
13.3
14.3
10.3
11.7

L34
[X = 1]
4.5 × 104
2.0 × 104
1.6 × 104
1.5 × 104
3.5 × 103
3.0 × 103
3.0 × 103
1.8 × 103
1.7 × 103
1.6 × 103
1.2 × 103
1.2 × 103
1.1 × 103
8.6 × 102

X · L34
[XXRB ]
8.2
95
46
80
—
12
12
17
< 0.1
7.4
< 0.1
0.15
17
13

X · L34
[XSB ]
< 0.1
370
< 0.1
< 0.1
—
4.9
4.9
< 0.1
0
0.51
0
0
< 0.1
< 0.1

and may impact superburst ignition.

5.2.1

Impact on Neutron Star Transients

During an accretion outburst, the crust of a neutron star transient is heated out of thermal equilibrium with the core and can reach temperatures where Urca cooling becomes important. When
accretion halts, a cooling wave propagates inward and the thermal emission powers the quiescent
light curve (Rutledge et al., 2002). Features in the crust temperature profile, such as local temperature dips caused by neutrino cooling, may be revealed in the quiescent light curve. The strength
of Urca cooling depends on the temperature reached during steady state in the neutron star crust.
In particular, the best candidates to possibly confirm the existence of Urca cooling are neutron
star quasi-persistent transients (see Chapter 4) where accretion rates near the Eddington accretion
˙ Edd ≈ 2 × 10−8 M yr−1 ) raise the neutron star’s temperature significantly. Of particurate (M
lar interest is the neutron star transient MAXI J0556-332 (hereafter MAXI; Homan et al. 2014)
90

1038

➤

1037
1036
10

➤

1034

1030

rp

A

=

29

1031

rb
ur

1032

st
-p
as
ro
he
ce
s
ss
as
he
s
br
em
sst
rah
lun
g
pla
sm
on

1033

su
pe

luminosity (erg s−1 )

crust heating
35

1029
107

108

109

temperature (K)
Figure 5.2: Luminosity sources in the neutron star crust as a function of crust temperature. The
red curves indicate the neutrino luminosities from Urca pairs. Blue curves indicate neutrino cooling from other processes in the neutron star crust. The gray shaded region indicates the heating
luminosity from crustal heating.

91

which is the hottest neutron star transient and reaches T ≈ 2 × 109 K in the crust (Deibel et al.,
2015) - making it an ideal testing ground for Urca cooling pairs. Here we detail the work done
in Deibel et al. (2015), where the quiescent observations from MAXI are compared against light
curve models with Urca pairs in the neutron star crust.
We fit the quiescent light curve of MAXI J0556-332 using the thermal evolution code dStar.
The best fit model is discussed in Section 4.2.1, with M = 1.5 M , R = 11 km, and Qshallow =
6.0 MeV per accreted nucleon. The best fit curve can be seen in Figure 5.3. As noted in Sectin 5.1,
Urca cooling pairs have two major impacts on the temperature of the crust during an accretion
episode. First, the Urca pair cools the crust locally through the emission of energy in the form of
neutrinos. Second, the Urca pair thermally decouples layers above and below the reaction layer
because heat is only conducted into the layer, not across the layer. In the case of MAXI, the
thermal decoupling most significantly impacts the shape of the quiescent light curve, as can be
seen in Figure 5.3. The observed MAXI light curve decreases monotonically for hundreds of days
and the cooling front has had time to enter the inner crust. If Urca cooling were present in MAXI,
the crust below the shallow heat source (ρ > ρh ) is thermally decoupled from the shallow heat
source near ρ ∼ 1010 g cm−3 by Urca cooling pairs present between ρ ∼ 1010 − 1011 g cm−3 . As
a result, the light curve “dips” as the cooling front moves deeper than the shallow heat source; that
is, the light curve after ≈ 20 days into quiescence.

5.2.2

Impact on Neutron Star Superbursts

Superbursts are powerful X-ray bursts thought to be triggered by unstable carbon ignition in
the neutron star envelope (Woosley & Taam, 1976; Fujimoto et al., 1981; Cumming & Bildsten, 2001; Strohmayer & Brown, 2002a). Superbursts are observed in neutron stars with ac˙ Edd (Wijnands, 2001) with inferred ignition depths between
cretion rates between ≈ 0.1–0.3 M
92

[ ]

350
300
250
200
150
100
50
100

=
/ = .

101

102
[ ]

103

104

Figure 5.3: Quiescent lightcurve of the neutron star MAXI J0556-332. The crust model is for a
M = 1.5 M , R = 11 km neutron star with Qshallow = 6.0 MeV. The light curve without Urca cooling is shown as a red dashed curve. Model light curves with Urca pairs, with Lν = 1036 ergs s−1 ,
are shown as black curves. From left to right, shell depths are y/yh = 3.4, 6.7, 17, 43, 140, that
correspond to ρ10 = 3, 5, 10, 20, 50, respectively.

93

[ ]

350
300
250
200
150
100
50

=
=
=
=
=

,

,

,
,
,

= , =
= , =
= , =
= , =
= , =

[

]

Figure 5.4: Quiescent lightcurve of the neutron star MAXI J0556-332 with an 33 Al Urca shell.
The crust model is for a M = 1.5 M , R = 11 km neutron star with Qshallow = 6.0 MeV.

94

accretion flow
-9

10

envelope
accumulation and compression

1

ocean (ion liquid)

106
superburst ignition
9

10

100

crust (ion lattice)

1010

Urca cooling pairs

depth [m]

mass density [g cm-3]

103

400
e--capture heating
pycnonuclear fusion heating

1014 core

103

1015

105

Figure 5.5: Schematic of the neutron star’s outer layers and the location of superburst ignition and
Urca cooling pairs.

95

Table 5.2: Properties of the strongest Urca e− -capture parents identified in Meisel & Deibel (2017)
and in Schatz et al. (2014), absent even-A nuclides, excluded by Meisel et al. (2015); Deibel et al.
(2016).
Parent |QEC | (MeV)

log( f t)

L34

XXRB

XSB

XS

5.1

8.2E+3

2.1E-3

1.9E-6

1.6E-4

11.8

4.9

4.2E+3

3.8E-3

4.3E-6

3.4E-4

13.4

5.2

3.7E+4

4.3E-3

4.0E-6

8.8E-5

55 Sc

12.1

4.9

2.4E+3

3.8E-3

1.8E-2

1.3E-3

57 Cr

8.3

11.6

8.6E-5

1.2E-3

1.6E-3

1.7E-3

57 V

29 Mg

13.3

31 Al
33 Al

10.7

4.9

1.2E+3

1.2E-3

1.6E-3

1.7E-3

59 Mn

7.6

11.6

5.2E-5

2.8E-3

3.1E-4

2.3E-3

63 Cr

14.7

14.4

1.1E-6

6.5E-3

6.5E-9

3.6E-3

65 Fe

10.3

11.6

2.1E-4

1.4E-2

4.2E-12

1.6E-2

65 Mn

11.7

11.6

4.1E-4

1.4E-2

4.2E-12

1.6E-2

yign ≈ 0.5–3 × 1012 g cm−2 (Cumming et al., 2006). Here we examine neutrino cooling from Urca
pairs in the crust to determine the impact on carbon ignition conditions in the neutron star ocean.
During an accretion outburst, the accretion-driven crustal heating produces a heat flux into the
overlying neutron star ocean. The ≈ 1–2 MeV per accreted nucleon of heat deposited in the crust,
however, is not sufficient to heat the ocean for superburst ignition at the depths observed. Superburst ignition models require an addition ≈ 1 MeV per accreted nucleon in the crust to provide
≈ 0.1 MeV per accreted nucleon extra heating into the ocean. This occurs because during steady
state ≈ 90% of crustal heating is transported toward the core (Brown, 2000). This extra heating
is needed to have superbursts ignite at the observed ignition depths (Cumming & Bildsten, 2001;
Cumming et al., 2006). Urca cooling pairs in the crust remove much of the heat flux entering the
ocean from accretion-driven crustal heating.
We model the ocean temperature with and without crust Urca pairs to determine their impact
on the ocean temperature and in turn their impact on unstable carbon ignition. We calculate the

96

neutron star ocean’s thermal evolution using the open source code dStar1 which solves the general
relativistic heat diffusion equation using the MESA numerical library (Paxton et al., 2011, 2013,
2015). The model has a neutron star mass M = 1.4 M , neutron star radius R = 10 km, and core
temperature T core = 3 × 107 K. The neutron star accretes until the ocean temperature reaches
steady state. The flux entering the ocean from the crust is F = Qb m,
˙ and we define Qb at y =
˙ Edd )−1 is required for
5 × 1014 g cm−2 . As shown in Cumming et al. (2006), Qb ≈ 0.25(m/0.3
˙
m
superburst ignition at y = 1012 g cm−2 .
Runaway thermonuclear burning of carbon occurs when the local heating rate from carbon
fusion C and the local cooling rate cool satisfy the condition d C /dT > d cool /dT . For the local
heating rate from carbon burning, we use the 12 C+12 C reaction rate given by Caughlan & Fowler
(1988) with screening by Ogata et al. (1993). The local cooling rate is given by cool = ρKT/y2
(Fujimoto et al., 1981) and approximates the heat lost through thermal diffusion. The unstable
carbon ignition curve in an iron ocean with XC = 0.2 is shown in Figure 5.9.
˙ Edd , the minimum acWe run a neutron star model with a local accretion rate m
˙ = 0.3 m
˙ Edd , ≈
cretion rate for ignition in an ocean with XC = 0.2 (Cumming et al., 2006), where m
8.8 × 104 g cm−2 s−1 is the local Eddington mass accretion rate. We place crust Urca pairs with
X · L34 > 1 from Table 5.1 into our superburst ignition model. We test several values of Qshallow
within constraints given by observations. For example, the neutron star transients MXB 1659-29
and KS 1731-260 require Qshallow ≈ 1 MeV during outburst to reconcile thermal evolution models
with quiescent observations (Brown & Cumming, 2009). For the same reason, the hottest neutron star transient to date, MAXI J0556-332, requires Qshallow ≈ 6–16 MeV per accreted nucleon
(Deibel et al., 2015). Specifically, we test Qshallow = 1, 5, 10, 15 MeV per accreted nucleon that
correspond to Qb = 0.1, 0.25, 1.0, 2.0 MeV per accreted nucleon, respectively. The ocean temper1 https://github.com/nworbde/dStar

97

1.2

]

0.8

[

1.0
0.6

(a) X-ray burst ashes

.
.

= .

0.4
0.2

]

0.8

[

1.0
0.6

(b) Superburst ashes

.
.
.
= .

0.4
0.2

]

0.8

[

1.0
0.6

.

(c) No Urca

.
.

= .

0.4
0.2 10
10

.

1011

1012

[

1013

]

1014

1015

Figure 5.6: Ocean temperature profiles during steady state with Urca pairs in the crust for various
values of Qb . Panel (a): Abundances of Urca pairs calculated from X-ray burst ashes. Panel (b):
Abundances of Urca pairs calculated from superburst ashes. Panel (c): Ocean temperature profiles
without Urca cooling in the crust.

98

ature profile for abundances of pairs in X-ray burst ashes (Woosley et al., 2004) and superburst
ashes (Keek & Heger, 2011) is shown in Figure 5.6 for the different values of Qb .
Crust Urca pairs are located deeper than superburst ignition at y
lower than deep crustal heating sources near y

1014 g cm−2 , but are shal-

1016 g cm−2 , as shown schematically in Fig-

ure 5.5. As a consequence, crust Urca pairs impact carbon ignition depths by removing ≈ 80–90%
of the heat flux entering the ocean from crustal and extra heating. As a result, for all values of Qb
tested, carbon ignition occurs deeper than it would otherwise when crust Urca pairs are present. In
a crust composed of X-ray burst ashes, the ocean temperature is limited to T
ignition region and ignition depths are limited to yign

2 × 1011 g cm−2 . In a crust composed of
7 × 108 K and carbon ignition depths

superburst ashes, the ocean temperature is limited to T
are limited to yign

109 K in the carbon

2 × 1012 g cm−2 . Once nuclear uncertainties are taken into account, allowed

carbon ignition depths are yign

1011 g cm−2 (yign

2 × 1011 g cm−2 ) in a crust composed of

X-ray burst (superburst) ashes.

5.3

Ocean Urca Pairs

As first observed by Gamow (1941), some nuclei pairs emit neutrinos more efficiently than others,
even if their abundance is small by comparison. This occurs because the Urca neutrino luminosity
scales as Q5EC and the neutrino luminosity will therefore vary between nuclei.
Here we focus on the complete identification of weaker Urca pairs in the A < 106 mass range
(the highest mass produced in the rp-process) with |QEC |

5 MeV. These Urca cycles occur in

the neutron star ocean, at column depths expected for superburst ignition yign = P/g ≈ 0.5–3 ×
1012 g cm−2 (Cumming et al., 2006). We identify a total of 85 odd-A isotopes with experimentally
determined mass excesses (Audi et al., 2012) which give QEC values accurate within ≈ 1% or

99

0.6

= .

1.2

XRB ashes

1.0

]

]

0.8

[

1.0

(a)

0.8

[

1.2

0.6

0.4

0.6

(b)

= .

0.2

XRB ashes

1.0

]

]

0.8

[

1.0

0.8
0.6

0.4

0.6

(c)

= .

0.2

XRB ashes

1.0

]

]

0.8

[

1.0

0.8
0.6

0.4

0.6

(d)

= .

0.2

XRB ashes

1.0

]

]

0.8

[

1.0

(f)

= .

SB ashes

(g)

= .

SB ashes

(h)

= .

SB ashes

0.4

0.8

[

0.2

SB ashes

0.4

[

0.2

= .

0.4

[

0.2

(e)

0.6

0.4
0.2
1010

0.4
1011

1012

1014

1013

[

0.2
1010

1015

]

1011

1012

1014

1013

[

1015

]

Figure 5.7: Ocean temperature profiles during steady state with Urca pairs in the crust for uncertainties in f t-values. We test the values Qb = 0.1, 0.25, 1, 2 MeV per accreted nucleon. Panels
(a)–(d): Abundances of crust Urca pairs calculated from X-ray burst ashes. Panels (e)–(h): Abundances of crust Urca pairs calculated from superburst ashes. The upper red-dashed curves are with
f t-values enhanced by a factor of 10 and the lower blue-dashed curves are for f t-values reduced
by a factor of 10.

100

better. In addition to ground-state to ground-state transitions, we include e− -captures onto excited
states with excitation energies EX ∼ kB T ∼ 100 keV for the typical temperatures near ∼ 109 K
in the ocean of superbursting neutron stars. Taken together with the study by Schatz et al. (2014)
of Urca cooling pairs in the neutron star’s crust, this study completes a census of Urca cooling
pairs that occur in the neutron star’s outer layers at mass densities below neutron-drip ρ

ρdrip ≈

4 × 1011 g cm−3 .
For each odd-A isotope, we calculate the total neutrino luminosity from e− -capture/β− -decay
cycles using analytic expressions for the reaction rates (Tsuruta & Cameron, 1970), as outlined in
Section 5.1. The f t-values, which measure the transition strength of e− -captures and β− -decays,
are based on experimental data2 when available; otherwise, f t-values are obtained from Table 1 of
Singh et al. (1998), using the experimentally determined ground-state spin-parities J π of the e− capture parent and daughter nuclei (Tuli, 2011). This approach to calculating f t-values provides
more reliable results than purely theoretical calculations using the quasi particle random phase
approximation for the special cases of interest here, where level energies, spins, and parities of all
relevant levels are known experimentally. The 15 strongest Urca pairs near the superburst ignition
depth are shown in Table 5.3.
We find that Urca cooling is present in all odd-A nuclei studied. In Figure 5.8, Urca pairs
are shown relative to the superburst ignition depths constrained from observations between y ≈
0.5–3 × 1012 g cm−2 (Cumming et al., 2006). The strongest cooling pairs occur at depths y
1.1 × 1012 g cm−2 . The two strongest cooling pairs are 23 Na–23 Ne and 25 Mg–25 Na, which occur
near y ≈ 8.6 × 1012 g cm−2 and y ≈ 4.9 × 1012 g cm−2 , respectively.
2 Evaluated Nuclear Structure Data File (ENSDF) — a computer file of experimental nuclear structure data maintained by the National Nuclear Data Center, Brookhaven National Laboratory (www.nndc.bnl.gov) — as of 11
December 2015

101

Table 5.3: Urca pairs active in the neutron star ocean.
Urca pair, ZA X

|QEC | [MeV] y12 1

L34 2

[MeV]

81 Br – 81 Se
1.59
0.1 0.10
34
35
49 Ti – 49 Sc
2.00
0.3 0.12
22
21
65 Cu – 65 Ni
2.14
0.4 0.02
29
28
55 Mn – 55 Cr
2.60
1.0 1.95
24
25
69 Zn – 69 Cu
2.68
1.1 0.78
30
29
57 Fe∗ – 57 Mn
2.70
1.2 3.45
26
25
67 Cu – 67 Ni
3.58
3.7 15.4
29
28
63 Ni∗ – 63 Co
3.66
4.1 14.9
28
27
25 Mg – 25 Na
3.83
4.9 20.3
12
11
81 Se – 81 As
3.86
5.0 4.54
34
33
73 Ga – 73 Zn
4.11
6.5 8.61
31
30
79 As – 79 Ge
4.11
6.5 7.68
33
32
23 Na – 23 Ne
4.38
8.4 42.1
11
10
101 Mo∗ – 101 Nb
4.63
11 8.92
42
41
57 Mn – 57 Cr
4.96
14 17.8
24
25
1 y ≡ y/(1012 g cm−2 ), calculated with g = 1.85.
12
14
2 Calculated with experimental f t-values when possible.

6
5
4 superburst ignition
3
2
1
0
10 20 30

40

50

60

70

80

90 100

Figure 5.8: Depth of Urca cooling pairs of a given mass number A. The size of data points corresponds to the coefficient L34 given in Equation 5.8, and can be found in Table 1. Urca pairs
with L34 ≥ 10 are colored in red. The gray band indicates empirical constraints on the superburst
ignition depth between yign ≈ 0.5–3 × 1012 g cm−2 (Cumming et al., 2006).

102

5.3.1

Observational Signatures

Many Urca pairs appear in the neutron star ocean (Deibel et al., 2016), although they have lower
neutrino luminosities than deeper pairs in the crust. These pairs cool the ocean by neutrino emission near the depths inferred for superburst ignition. In this section, we model the ocean temperature with and without Urca pairs to determine their impact on unstable carbon ignition. The flux
entering the ocean from the crust is F = Qb m.
˙ The flux entering the ocean must be higher than
supplied by crustal heating for models to ignite at observed depths. For example, as shown in Cum˙ Edd )−1 for carbon ignition at y = 1012 g cm−2 .
ming et al. (2006), there must be Qb ≈ 0.25(m/0.3
˙
m
We calculate the neutron star ocean’s thermal evolution using the open source code dStar
(Brown, 2015) which solves the general relativistic heat diffusion equation using the MESA numerical library (Paxton et al., 2011, 2013, 2015). The model has a neutron star mass M = 1.4 M ,
neutron star radius R = 10 km, and core temperature T core = 3 × 107 K. The neutron star accretes
until the ocean temperature reaches steady state.
Runaway thermonuclear burning of carbon occurs when the local heating rate from carbon
fusion C and the local cooling rate cool satisfy the condition d C /dT > d cool /dT . For the local
heating rate from carbon burning, we use the 12 C+12 C reaction rate given by Caughlan & Fowler
(1988) with screening by Ogata et al. (1993). The local cooling rate is given by cool = ρKT/y2
(Fujimoto et al., 1981). The unstable carbon ignition curve in an iron ocean with XC = 0.2 is
shown in Figure 5.9.
˙ Edd , where m
˙ Edd , ≈ 8.8 ×
We run a neutron star model with a local accretion rate m
˙ = 0.3 m
104 g cm−2 s−1 . The heat flux entering the ocean from the crust is F = Qb m
˙ at the location
of the ocean-crust boundary where Γ = 175 (Farouki & Hamaguchi, 1993). We use Qb =
0.25 MeV per accreted nucleon to give ignition at y = 1012 g cm−2 following the relation Qb ≈

103

˙ Edd )−1 (Cumming et al., 2006).
0.25(m/0.3
˙
m
In Figure 5.9 we show the ocean temperature profile with and without Urca pairs for our model.
In an ocean without Urca pairs, superburst ignition occurs at yign ≈ 1.0 × 1012 g cm−2 . When Urca
pairs are located at y > yign at y = 1.2 × 1013 g cm−2 the superburst ignition depth increases (see
panel (a) of Figure 5.9). For an Urca pair with L34 = 10 (L34 = 100) the superburst ignites at a
column y = 1.2 × 1012 g cm−2 (y = 1.2 × 1012 g cm−2 ).
When Urca pairs are located at y < yign the superburst ignition depth increases. In panel (c)
of Figure 5.9, are ocean temperature profiles for an Urca shell located at y = 1.2 × 1010 g cm−2 .
For Urca pairs with X · L34 = 10 (X · L34 = 100) the superburst ignites at a column yign =
1.1 × 1012 g cm−2 (yign = 1.9 × 1012 g cm−2 ).
To determine if cooling from ocean Urca pairs may alter the cooling superburst light curve,
we model superburst cooling following the approach of Cumming & Macbeth (2004). In this approach, unstable carbon burning deposits heat at the carbon ignition depth on a timescale shorter
than a thermal diffusion time. The heat deposition raises the oceans temperature to T

109 K

and sets a steep temperature profile that acts as the initial temperature profile for the subsequent
cooling superburst light curve. The average energy input from superburst burning is E18 ≡
Enuc /(1018 ergs g−1 ) and superburst cooling light curves are typically well fit with E18 ≈ 0.2.
The initial temperature profile is a power law T ∝ yα in column depth and column depth can be
approximated as y = P/g ≈ 3.6 × 1012 g cm−2 ρ94/3 (Ye /0.4)4/3 , where ρ9 ≡ ρ/(109 g cm−3 ), and
superbursts are typically well fit with α ≈ 0.25.
We examine superburst cooling light curves for two ignition depths yign = 3.6 × 1011 and
3.6 × 1012 g cm−2 . We include ocean Urca pairs at y = 0.4 yign with the values X · L34 = 100, 1000
for both ignition depths. The results of these models can be seen in Figure 5.10. Urca cooling
in the ocean only begins to affect the cooling light curve when X · L34
104

100, when the neutrino

]
[
]
[
]
[

1.0
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.9
0.8
0.7
0.6
0.5
0.4
0.3 10
10

(a)

=

(b)

=

(c)

=

1011

1012

[

1013

]

1014

1015

Figure 5.9: Temperature in the ocean as a function of column depth for a neutron star model
˙ Edd , M = 1.4 M , R = 10 km, and core temperature T core = 3 × 107 K. The
with: m
˙ = 0.3 m
dashed-black curves indicate the ignition of unstable carbon burning in a mixed iron-carbon ocean
with XFe = 0.8 and XC = 0.2. The solid-black curves indicate ocean’s without Urca cooling
layers. In each panel, the lower dotted-blue curve is for L34 = 100 and the upper dotted-blue
curve is for L34 = 10, for Urca pairs located at: Panel (a): y = 1.2 × 1013 g cm−2 , Panel (b):
y = 1.2 × 1012 g cm−2 , and Panel (c): y = 1.2 × 1010 g cm−2 .

105

[

]

1038

(

= .

)= .
.

1037
1036
101

102

103

104
[]

105

106

Figure 5.10: Cooling light curves for superbursts with E18 = 0.2, α = 0.25, g14 = 1.6, and
ignition depths yign = 3.2 × 1011 g cm−2 and yign = 3.2 × 1012 g cm−2 . The red dotted curves
are models without Urca cooling. The solid curves are for a cooling source at y = 0.4 yign with
X · L34 = 100, 1000.
luminosity from the Urca pair Lν becomes comparable to the superburst peak luminosity Lpeak . At
realistic abundances of ocean pairs from X-ray burst and superburst ashes, however, all ocean pairs
at these depths have X · L34

1 and Lν

Lpeak . Therefore, we predict that ocean Urca pairs do

not noticeably change superburst cooling light curves.
At large enough abundances ocean Urca pairs have the potential to cool the ocean and increase
carbon ignition depths. At the abundances calculated from X-ray burst and superburst nuclear
burning (Schatz et al., 2014), however, ocean pairs have weak neutrino luminosities (X · L34

1)

and do not impact carbon ignition depths. Furthermore, the neutrino luminosities of ocean Urca
pairs located at y

yign are too weak to cause an early decline of the superburst light curve. That is,

ocean pairs have no effect on the cooling light curve because their neutrino luminosities are small
relative to the superburst peak luminosity (Lν
T

Lpeak ) even in the hot post-superburst ocean near

109 K.

106

Our results suggest that the cooling thermal component of the February 2001 superburst in
4U 1636-536 (Wijnands, 2001; Strohmayer & Markwardt, 2002; Keek et al., 2014a,b) is likely
not caused by Urca cooling in the ocean. The cooling light curve declines much faster after the
superburst peak than model predictions (with uncertainties about the spectral model used to fit the
data, see discussion in Keek et al. 2015). Because the light curve declines after the superburst
peak, Urca cooling would need to be near the carbon ignition depth to explain the cooling trend
(Keek et al., 2015). For the ignition depth near that of the superburst in 4U 1636-536 (yign ≈
2×1011 g cm−2 ), only Urca cooling pairs present at depths y
however, these ocean pairs are typically weak with X ·L34

yign will impact cooling predictions;
1 and do not effect cooling light curve

predictions. Indeed, Urca cooling in the ocean only marginally effects light curve predictions even
for X · L34 = 1000, as can be seen in Figure 5.10.
The shallow ignition depth in 4U 1636-536 near yign ≈ 2 × 1011 g cm−2 (Keek et al., 2015)
provides constraints on the strength of extra heating and the composition of the crust in this source.
A shallow heating source Qshallow ≈ 9 MeV per accreted nucleon is the minimum strength required
to ignite carbon at this depth in a crust composed of X-ray burst ashes; if composed of superburst
ashes the minimum heating strength is Qshallow ≈ 15 MeV per accreted nucleon. A weaker shallow
heat source is consistent with neutron star transients that have similar accretion rates between
˙ Edd , for instance MXB 1659-29 requires ≈ 1 MeV per accreted nucleon (Brown &
m
˙ ∼ 0.1–0.3 m
Cumming, 2009). Moreover, superburst ignition models find that several type-I X-ray bursts occur
during the intervals between the more energetic superbursts (Keek et al., 2012). It is likely that the
crust in 4U 1636-536 is composed primarily of X-ray burst ashes, which accumulate in the crust
over the course of ∼ 100–1000 years in the more frequent type-I X-ray bursts, and cooling from
Urca pairs in these ashes prevails over cooling from pairs in the less abundant superburst ashes.
By contrast, extra heating may be deposited shallower than most crust Urca pairs y
107

2×

1.2

= .

]

0.8

[

1.0
0.6
0.4

.
.

= .

0.2
1010

1011

1012

[

1013

]

1014

1015

Figure 5.11: Steady state ocean temperature as a function of column depth for a neutron star model
˙ Edd , M = 1.4 M , R = 10 km, and core temperature T core = 3 × 107 K. Solid
with: m
˙ = 0.3 m
black curves are steady state temperature profiles for the heating strengths Qshallow = 1, 2, 4. The
dashed black curve indicates the ignition of unstable carbon burning in a mixed iron-carbon ocean
with XFe = 0.8 and XC = 0.2.
1014 g cm−2 as is the case in neutron star transients with extra shallow heating; for example,
KS 1731-260 and MXB 1659-29 (Brown & Cumming, 2009), and MAXI J0556-332 (Deibel et al.,
2015). When deposited at shallower depths, an extra heating strength of Qshallow ≈ 4 MeV per accreted nucleon is required to have carbon ignition at depths consistent with the superburst ignition
depth in 4U 1636-536, as can be seen in Figure 5.11. In this case, no constraints on crust composition are possible because crust Urca pairs do not alter the temperature in the superburst ignition
region. It may be interesting to search for transients that show superbursts during an accretion
outburst because the quiescent cooling may potentially indicate the existence of extra heating and
its depth. A comparison of the strength and depth of extra heating with the inferred superburst
ignition depth will then allow constraints on the crust composition.
The nuclear properties of neutron-rich nuclei are important input for determining the strength

108

of Urca cooling. Although the QEC values for nuclei in the ocean are determined from mass
measurements to within ≈ 1 % (Audi et al., 2012), f t-values and excited state energies are more
uncertain. Shell model calculations of all odd-A nuclei in rp-process burning ashes are needed
to determine excited state energies and transition strengths to higher accuracy. Furthermore, experimentally determined f t-values and mass measurements of neutron-rich nuclei, for instance
at the Facility for Rare-Isotope Beams (FRIB), will allow more accurate determinations of Urca
pair neutrino luminosities and thereby improve constraints on the interior of superbursting neutron
stars.

109

Chapter 6
Conclusion
This dissertation explores the neutron star interior by reconciling numerical models of the neutron
star’s outer layers with neutron star observations. In particular, observations of magnetar giant
flares and cooling neutron star transients can be used to reveal the properties of dense matter in
the neutron star interior that would otherwise be difficult to observe in a terrestrial laboratory. We
build a numerical model of the neutron star’s outer layers to simulate phenomena in these layers.
We then model and analyze a variety of neutron star observations to draw new constraints relevant
for nuclear theory and future astrophysical modeling of neutron stars.
In Chapter 2, we outline the current understanding of the microphysics of the neutron star’s
outermost layers. This numerical model is the foundation for all subsequent investigations of the
interior physics. The model uses theoretical mass models that are fit to experimental nuclear data
to compute the equilibrium composition of the neutron star’s ocean and crust. The composition
is needed to compute oscillation modes of the crust in Chapter 3. Because the composition and
low-density equation of state set the thermal transport properties of the outer layers, our model can
then be used to investigate the thermal evolution of neutron star transients in Chapter 4.
In Chapter 3, we examined crust oscillation modes in highly magnetized neutron stars by
adding a magnetic field to our numerical model of the neutron star crust. The frequencies of
crust torsional oscillations are a compelling match to observed quasi-periodic oscillations in the
tails of magnetar giant flares. Crust oscillation frequencies are uniquely set by the magnetar’s mass
and radius, as well as the crust composition. Only one mass and radius combination can produce
110

the observed oscillation frequencies in a given giant flare. Furthermore, if observed quasi-periodic
oscillations are indeed oscillating crust modes, then models of crust oscillations may reveal the
properties of the densest regions of the neutron star’s crust near the crust-core transition—a density regime where constraints on the nuclear symmetry energy are difficult to obtain.
In Chapter 4, we investigate the quiescent light curves of neutron star X-ray transients. Though
some transients are known to require extra heating during their accretion outbursts to match observed quiescent temperatures, such as KS 1731-260 and MXB 1659-29 (Brown & Cumming,
2009), we find that the hottest transient MAXI J0556-332 requires a surprisingly large amount of
heat deposition during outburst. The quiescent temperature of MAXI J0556-332 requires

6 MeV

per accreted nucleon of heat be deposited during outburst; a factor of 6 larger than previously observed shallow heating in other sources. Interestingly, superburst models require extra heating in
the ocean and shallow crust to provide a heat flux into the envelope sufficient to ignite carbon at
depths consistent with observed light curves (Keek et al., 2008; in’t Zand et al., 2012; Zamfir et al.,
2014). Furthermore, the extra heating is required for superburst models to match observed superburst recurrence times (Cumming & Bildsten, 2001; Strohmayer & Brown, 2002b). The source of
the extra heating, however, remains unknown, and should be the focus of future work. A plausible source of the extra heating may be the rotational energy reservoir in the accretion flow. The
Keplerian energy of the accretion flow is ∼ 80 MeV per accreted nucleon (at the inner-most stable
circular orbit) and the viscous dissipation of rotational energy may occur in an accreted spreading
boundary layer (Inogamov & Sunyaev, 1999, 2010; Philippov et al., 2016).
The late time cooling of neutron star transients should also be investigated further. In Section 4.4, we discuss the late time cooling observations of the neutron star transient MXB 1659-29.
The neutron star appears to resume cooling around ≈ 2500 d which may be indicative of a layer of
normal neutrons with a long thermal time that forms as a consequence of the low thermal conduc111

tivity of nuclear pasta. The normal neutron layer requires that the 1 S0 neutron singlet pairing gap
closes in the crust before the crust-core transition. As a result, constraints on the neutron singlet
pairing gap model are possible through fitting various gap models to the late time light curve of
MXB 1659-29. In particular, the transition density to normal neutrons can be determined, though
the late time cooling will be insensitive to the energy of the gap itself. In a similar manner, theoretical singlet pairing gap models can be examined to determine which of the existing models best
fits the late time cooling observations.
We also outline an observational test to answer an open question in neutron star transient observations (Section 4.3). Observed quasi-periodic oscillations in the transient neutron star Z-sources,
such as Sco X-1, are of unknown origin. Curiously, some observed frequencies are close to the
neutron star ocean’s thermal g-mode frequencies. We show that as the ocean heats up during an
accretion outburst the ocean’s thermal g-mode frequencies increase and when accretion halts, the
predicted g-mode frequencies decrease as the ocean cools. Furthermore, we predict that observed
oscillation frequencies should be larger in the hottest neutron star transients that require anomalous
shallow heating, for example, MAXI J0556-332. These predictions can be tested by further X-ray
observations of neutron star Z-sources along the normal branch.
In Chapter 5 we discuss the recent discovery of Urca cooling in nuclear reaction network calculations of the accreted crust. Schatz et al. (2014) found nuclear pairs that undergo cycles of
e− −capture /β− −decay appear in accreted material compressed into the neutron star’s crust. Urca
cycling pairs cool the neutron star crust and can counteract heating from accretion-driven nuclear
reactions. We find that Urca cooling impacts quiescent light curve predictions of neutron star transients (Deibel et al., 2015; Meisel & Deibel, 2017), and Urca cooling is unlikely to be present in
MAXI J0556-332. We expand upon Schatz et al. (2014) and find that Urca cycling nuclei pairs
appear in the ocean at lower densities as well (Deibel et al., 2016). Although it is unlikely that
112

ocean Urca pairs cool the ocean sufficiently to alter carbon ignition depths in superbursts, we find
that Urca pairs located in the deeper crust potentially alter ignition depths. Urca pairs in the crust
prevent crustal heating from entering the overlying ocean and carbon ignition occurs deeper than
it would otherwise.

113

APPENDICES

114

Appendix A: Local Stability of Non-equilibrium Nuclear Reactions
In the inner crust, an compressed nucleus will undergo non-equilibrium reactions in the degenerate
electron and quasi-free neutron gases therein. For example, a nucleus with proton number Z and
mass number A will undergo an electron capture in an electron gas with chemical potential µe ,
when

Zµe > B(Z, A)− Nmn −Z(m p +me )− B(Z −1, A)+(A−Z +1)mn +(Z −1)(m p +me )+(Z −1)µe , (1)

or
µe > B(Z, A) − B(Z − 1, A) − ∆ ,

(2)

where B(Z, A) is the binding energy of the nucleus and ∆ = mn − m p − me is the binding energy of
the neutron. The same conditions for other possible reactions are:

β decay : µe > −B(Z, N) + B(Z + 1, A) + ∆

(3)

Neutron capture : µn < B(Z, A) − B(Z, A + 1)

(4)

Neutron emission : µn > −B(Z, A) + B(Z, A − 1)

(5)

Combining the capture inequalities and the emission inequalities gives,

µe + µn < 2B(Z, A) − B(Z − 1, A) − B(Z, A + 1) + ∆ ,

(6)

µe + µn > −2B(Z, A) + B(Z + 1, A) + B(Z, A − 1) + ∆ ,

(7)

115

Now combine these inequalities to remove chemical potentials,

1
[B(Z + 1, A) + B(Z, A − 1) + B(Z − 1, A) + B(Z, A + 1)] < B(Z, A) ,
4

(8)

which looks like the nucleus binding energy has to be larger than the average of the binding energies of the surrounding nuclei in order for the nucleus to be stable against non-equilibrium reactions. This is identically expressed as the condition for local stability,
∂2 B ∂2 B
+
< 0 ; ∇2 < 0 .
∂Z 2 ∂A2

116

(9)

Appendix B: Thermal Diffusion
Here we derive the fully general relativistic thermal evolution equation for the neutron star crust
(Equation 2.25) following the approach of Thorne (1977). We begin with Equation (11d) of Thorne
(1977),
1 ∂Π 1 P ∂ρ
1 ∂(LR2 )
=
ε
−
ε
−
+
,
nuc
ν
R ∂t
R ρ2 ∂t
R2 ∂M

(10)

where R = exp(Φ/c2 ) is the gravitional redshift and Φ is the gravitational potential of the neutron
star. the compression term is negligible. Given that C p = ∂Π/∂T for degnerate matter, Equation 10
can be rewritten as
1 ∂(LR2 )
1
∂T
= εnuc − εν −
Cp
2
R
∂t
R ∂M

,

(11)

The LHS of Equation 11 is found by starting with Equation (11h) of Thorne (1977) ,
∂ lnT
∂ lnP
= ∇rad
,
∂Mr
∂Mr

(12)

and Equation (9a) Thorne (1977) ,

∇rad =

κLr P
ε
3
1
+
1
−
,
64π GMr σT 4 HGV
H

(13)

where the correction factors are: H enthalpy, G gravitational-acceleration, V volume, and ε energy.
Plugging Equation 13 into Equation 12 gives

1 ∂T
κLr P
3
1
ε
=
+
1
−
T ∂Mr
64π GMr σT 4 HGV
H
Using Equation (11i) from Thorne (1977) we have the relation

117

∂lnP
.
∂Mr

(14)

∂lnP 1 ∂P
1 GMr
GHV .
=
=
−
∂Mr
P ∂Mr P 4πr4

(15)

Inserting Equation 15 into Equation 14 and noting that H − ε = P/(ρc2 ), gives
1 ∂T
3 κLr 1
GMr GV
=−
−
,
4
4
T ∂Mr
64π acT 4πr
4πr4 ρc2

(16)

GMr GV ∂Φ/c2
=
,
∂Mr
4πr4 ρc2

(17)

1 ∂T
1 ∂Φ
3 κLr 1
+ 2
=−
.
T ∂Mr c ∂Mr
64π acT 4 4πr4

(18)

and using the relation

we can express Equation 16 as

Multiplying by 4T 4 and solving for Lr gives,

Lr = −

∂
4ac
(4πr2 )2 e−4Φ
(e4Φ T 4 ) .
3κ
∂Mr

(19)

∂ R2 ,
Multiplying both sides by 12 ∂M
r
R

1 ∂(Lr R2 )
1 ∂
ac
∂
= 2
−R2 (4πr2 )2 e−4Φ
(e4Φ T 4 )
2
3κ
∂Mr
R ∂Mr
R ∂Mr

(20)

and we have recovered the LHS of Equation 11. Plugging Equation 20 into Equation 11 and

118

simplifying, gives

R−1C p

1 ∂
d
dT
ac
=− 2
(e4Φ T 4 ) + εnuc − εν ,
−R2 (4πr2 )2 e−4Φ
dt
3κ
dMr
R ∂Mr

(21)

The derivative could also be done with respect to the radial coordinate because
∂r
= (4πr2 ρV)−1 ,
∂Mr

(22)

and it can be easily shown that the diffusion equation in terms of the radial coordinate is
1
∂ ¯
∂
dT
=
K(r, T ) (R4 T 4 ) + εnuc − εν ,
dt
C(r, T ) ∂r
∂r

(23)

3C p r2 VR
C(r, T ) =
,
ac

(24)

2
¯ T) = r .
K(r,
κR2

(25)

where

and

119

Appendix C: Axial Perturbations
Here we formulate a version of the torsional perturbation equation appropriate for a highly magnetized neutron star. The ad-hoc form of the perturbation equation will include general relativistic
corrections appropriate for the large neutron star gravity and a finite Alv´en velocity to account for
a large magnetic field in a magnetar. We begin by following the derivation of axial perturbations
from Samuelsson & Andersson (2007) in the Cowling approximation (Cowling, 1941) (i.e., we
ignore changes in the gravitational potential over the thin crust). The perturbation equation is
e2λ
d
ln r4 eν−λ (ε + pt ) v2r ξ + 2
ξ +
dr
vr



 −2ν 2 v2t ( − 1) ( + 2) 
e ω −
 ξ = 0 ,
r2

(26)

where ε is the energy density, ρ is the rest mass energy density, r is the radial coordinate, ω is the
angular frequency, and e2λ ≈ e−2ν ≈ (1 − 2M/(Rc2 ))−1 in the thin crust. Note that the argument of
the logarithm can be made unitless with the appropriate combination of

and c. We will take the

isotropic limit pt = p and vr 2 = vt 2 = v s 2 = µ/ρ. Evaluating the first coefficient,

d
ln r4 eν−λ (ε + pt ) v2r
dr
+

= r4 eν−λ (ε + p)v2s

−1

[

d 4 ν−λ
r (e (ε + p)v2s )
dr

d ν−λ 4
d
d 2 4 ν−λ
e
(r (ε + p)v2s ) +
v (r e (ε + p))] ,
{ε + p} (r4 eν−λ v2s ) +
dr
dr
dr s

(27)

where the result is written to illustrate clearly that the product rule is being used on four functions.
Evaluating each term individually,

d 4
r
dr

= 4r3 ,

d ν−λ
d
2M(r)
dM(r) 1 M(r)
2 M(r) dM(r)
e
=
1−
= −2
− 2 =
−
,
dr
dr
r
dr r
r
r
dr
r
d
dε d p
+
,
{ε + p} =
dr
dr dr
120

d 2
v
dr s

=

dρ
d µ
dµ 1
µ dρ 1 dµ
=
− 2
=
− v2s
,
dr ρ
dr ρ ρ dr ρ dr
dr

where we assume dB/dr = 0, (the magnetic field is not changing with depth in the star). Putting
together all of the pieces in the first coefficient,

A=

+

r4

−1
M(r)
2
(ε + p)v s
[ 4r3 (eν−λ (ε + p)v2s )
1−2
r

2 M(r) dM(r)
−
r
r
dr
+

(r4 (ε + p)v2s ) +

dρ
1 dµ
− v2s
ρ dr
dr

dε d p 4 ν−λ 2
+
(r e v s )
dr dr

(r4 eν−λ (ε + p)) ] .

(28)

(29)

(30)

The differential equation then becomes,





−1
−1
2 (l − 1)(l + 2) 


v2A 
v
M(r)
M(r)
1

s
1 +  ξ + Aξ + 1 − 2
 1 − 2
(1 + v2A )ω2 −
 ξ = 0 ,
2
2
2
r
r
r
vs
vs
(31)






2
2



v 
 e2λ  2λ

e (1 + v2 )ω2 − v s (l − 1)(l + 2) 
1 + A  ξ + Aξ + 
ξ=0.
(32)




A
 v2

r2
v2
s

s

This is the expression used in Deibel et al. (2014). The Newtonian expression in Piro (2005)
can be recovered in the limit e2λ → 1.

121

REFERENCES

122

REFERENCES

Andersson, N. 1998, ApJ, 502, 708
Andersson, N., Glampedakis, K., Ho, W. C. G., & Espinoza, C. M. 2012, Physical Review Letters,
109, 241103
Andersson, N., Jones, D. I., & Ho, W. C. G. 2014, MNRAS, 442, 1786
Audi, G., Wang, M., Wapastra, A. H., Kondev, F. G., MacCormick, M., Xu, X., & Pfeiffer, B.
2012, Chinese Physics C, 36, 2
Baade, W., & Zwicky, F. 1934, Proceedings of the National Academy of Science, 20, 259
Bahcall, J. N., & Wolf, R. A. 1965, Physical Review, 140, 1452
Bahramian, A., Heinke, C. O., Wijnands, R., & Degenaar, N. 2016, The Astronomer’s Telegram,
8699
Baiko, D. A., & Yakovlev, D. G. 1995, Astronomy Letters, 21, 702
—. 1996, Astronomy Letters, 22, 708
Barat, C., et al. 1983, Astron. & Astrophys., 126, 400
Bardeen, J. M., & Petterson, J. A. 1975, ApJ, 195, L65
Baym, G., Bethe, H. A., & Pethick, C. J. 1971a, Nuclear Physics A, 175, 225
Baym, G., & Chin, S. A. 1976, Nuclear Physics A, 262, 527
Baym, G., Pethick, C., & Pines, D. 1969, Nature, 224, 673
Baym, G., Pethick, C., & Sutherland, P. 1971b, ApJ, 170, 299
Belloni, T., Parolin, I., & Casella, P. 2004, A&A, 423, 969
Bildsten, L., & Brown, E. F. 1997, ApJ, 477, 897
Bildsten, L., & Cumming, A. 1998, ApJ, 506, 842
Bildsten, L., & Cutler, C. 1995, ApJ, 449, 800
Bildsten, L., Ushomirsky, G., & Cutler, C. 1996, ApJ, 460, 827
123

Bisnovaty˘i-Kogan, G. S., & Chechetkin, V. M. 1979, Soviet Physics Uspekhi, 22, 89
Brehm, J. J. 1989, Introduction to the Structure of Matter (Wiley)
Broderick, A., Prakash, M., & Lattimer, J. M. 2000, ApJ, 537, 351
Brown, B. A. 2013, Physical Review Letters, 111, 232502
Brown, B. A., Clement, R. R. C., Schatz, H., Volya, A., & Richter, W. A. 2002, Phys. Rev. C, 65,
045802
Brown, E. F. 2000, ApJ, 531, 988
—. 2015, dStar: Neutron star thermal evolution code, Astrophysics Source Code Library
Brown, E. F., & Bildsten, L. 1998, ApJ, 496, 915
Brown, E. F., Bildsten, L., & Chang, P. 2002, ApJ, 574, 920
Brown, E. F., & Cumming, A. 2009, ApJ, 698, 1020
Cackett, E. M., Brown, E. F., Cumming, A., Degenaar, N., Fridriksson, J. K., Homan, J., Miller,
J. M., & Wijnands, R. 2013, ApJ, 774, 131
Cackett, E. M., Wijnands, R., Miller, J. M., Brown, E. F., & Degenaar, N. 2008, ApJ, 687, L87
Casella, P., Belloni, T., & Stella, L. 2006, A&A, 446, 579
Caughlan, G. R., & Fowler, W. A. 1988, Atomic Data and Nuclear Data Tables, 40, 283
Chabanat, E., Bonche, P., Haensel, P., Meyer, J., & Schaeffer, R. 1995, Phys. Scripta, T56, 231
Chadwick, J. 1932, Nature, 129, 312
Chamel, N. 2005, Nuclear Physics A, 747, 109
—. 2012, Phys. Rev. C, 85, 035801
Chamel, N. 2013, Phys. Rev. Lett., 110, 011101
Chamel, N., Page, D., & Reddy, S. 2013a, Phys. Rev. C, 87, 035803
—. 2013b, Phys. Rev. C, 87, 035803
Chen, J. M. C., Clark, J. W., Dav´e, R. D., & Khodel, V. V. 1993, Nuclear Physics A, 555, 59

124

Cooper, R. L., Steiner, A. W., & Brown, E. F. 2009, ApJ, 702, 660
Cowling, T. G. 1941, MNRAS, 101, 367
Cumming, A., & Bildsten, L. 2001, ApJ, 559, L127
Cumming, A., & Macbeth, J. 2004, ApJ, 603, L37
Cumming, A., Macbeth, J., in ’t Zand, J. J. M., & Page, D. 2006, ApJ, 646, 429
Cyburt, R. H., et al. 2010, ApJS, 189, 240
Degenaar, N., Brown, E. F., & Wijnands, R. 2011, MNRAS, 418, L152
Degenaar, N., & Wijnands, R. 2011, MNRAS, 414, L50
Degenaar, N., et al. 2009, MNRAS, 396, L26
—. 2013, ApJ, 775, 48
—. 2014, ApJ, 791, 47
—. 2015, MNRAS, 451, 2071
Deibel, A. 2016, ApJ, 832, 44
Deibel, A., Cumming, A., Brown, E. F., & Page, D. 2015, ApJ, 809, L31
Deibel, A., Cumming, A., Brown, E. F., & Reddy, S. 2017, ApJ, 839, 95
Deibel, A., Meisel, Z., Schatz, H., Brown, E. F., & Cumming, A. 2016, ApJ, 831, 13
Deibel, A. T., Steiner, A. W., & Brown, E. F. 2014, Phys. Rev. C, 90, 025802
Dieperink, A. E. L., & van Isacker, P. 2009, European Physical Journal A, 42, 269
Dubus, G., Kern, B., Esin, A. A., Rutledge, R. E., & Martin, C. 2004, MNRAS, 347, 1217
Duncan, R. C. 1998, Astrophys. J. Lett., 498, L45
Eichler, D., & Cheng, A. F. 1989, ApJ, 336, 360
Farouki, R. T., & Hamaguchi, S. 1993, Phys. Rev. E, 47, 4330
Flowers, E., & Itoh, N. 1976, ApJ, 206, 218

125

Fridriksson, J. K., Homan, J., & Remillard, R. A. 2015, ApJ, 809, 52
Fridriksson, J. K., et al. 2010, ApJ, 714, 270
—. 2011, ApJ, 736, 162
Fujimoto, M. Y., Hanawa, T., & Miyaji, S. 1981, ApJ, 247, 267
Gamow, G. 1941, Physical Review, 59, 617
Gamow, G., & Schoenberg, M. 1940, Physical Review, 58, 1117
—. 1941, Physical Review, 59, 539
Gandolfi, S., Illarionov, A. Y., Fantoni, S., Pederiva, F., & Schmidt, K. E. 2008, Physical Review
Letters, 101, 132501
Gezerlis, A., & Carlson, J. 2010, Phys. Rev. C, 81, 025803
—. 2011, ArXiv e-prints
Giacconi, R., Gursky, H., Paolini, F. R., & Rossi, B. B. 1962, Physical Review Letters, 9, 439
Glampedakis, K., Samuelsson, L., & Andersson, N. 2006, Mon. Not. Royal Astron. Soc., 371, L74
Gudmundsson, E. H., Pethick, C. J., & Epstein, R. I. 1982, ApJ, 259, L19
—. 1983, ApJ, 272, 286
Guillot, S., Servillat, M., Webb, N. A., & Rutledge, R. E. 2013, ApJ, 772, 7
Gupta, S., Brown, E. F., Schatz, H., M¨oller, P., & Kratz, K.-L. 2007, ApJ, 662, 1188
Gupta, S. S., Kawano, T., & M¨oller, P. 2008, Phys. Rev. Lett., 101, 231101
Haensel, P., & Zdunik, J. L. 1990, A&A, 227, 431
—. 2003, A&A, 404, L33
—. 2008, A&A, 480, 459
Hairer, E., Norsett, S. P., & Wanner, G. 1993, Solving Ordinary Differential Equations I. Nonstiff
problems, 2nd edn., Springer Series in Computational Mathematics (Springer-Verlag)
Hashimoto, M., Seki, H., & Yamada, M. 1984, Progress of Theoretical Physics, 71, 320

126

Hasinger, G., & van der Klis, M. 1989, A&A, 225, 79
Hasinger, G., van der Klis, M., Ebisawa, K., Dotani, T., & Mitsuda, K. 1990, A&A, 235, 131
Henyey, L., & L’Ecuyer, J. 1969, ApJ, 156, 549
Hertz, P., Vaughan, B., Wood, K. S., Norris, J. P., Mitsuda, K., Michelson, P. F., & Dotani, T. 1992,
ApJ, 396, 201
Hewish, A., Bell, S. J., Pilkington, J. D. H., Scott, P. F., & Collins, R. A. 1968, Nature, 217, 709
Ho, W. C. G., Espinoza, C. M., Antonopoulou, D., & Andersson, N. 2015, Science Advances, 1,
e1500578
Homan, J., Fridriksson, J. K., Wijnands, R., Cackett, E. M., Degenaar, N., Linares, M., Lin, D., &
Remillard, R. A. 2014, ApJ, 795, 131
Homan, J., Linares, M., van den Berg, M., & Fridriksson, J. 2011, The Astronomer’s Telegram,
3650, 1
Homan, J., et al. 2007, ApJ, 656, 420
—. 2010, ApJ, 719, 201
Horowitz, C. J., & Berry, D. K. 2008, Phys. Rev. C, 78, 035806
—. 2009, Phys. Rev. C, 79, 065803
Horowitz, C. J., Berry, D. K., Briggs, C. M., Caplan, M. E., Cumming, A., & Schneider, A. S.
2015, Physical Review Letters, 114, 031102
Horowitz, C. J., P´erez-Garc´ıa, M. A., Carriere, J., Berry, D. K., & Piekarewicz, J. 2004a,
Phys. Rev. C, 70, 065806
Horowitz, C. J., P´erez-Garc´ıa, M. A., & Piekarewicz, J. 2004b, Phys. Rev. C, 69, 045804
Horowitz, C. J., Schneider, A. S., & Berry, D. K. 2010, Physical Review Letters, 104, 231101
Iida, K., & Oyamatsu, K. 2014, European Physical Journal A, 50, 42
Inogamov, N. A., & Sunyaev, R. A. 1999, Astronomy Letters, 25, 269
—. 2010, Astronomy Letters, 36, 848
in’t Zand, J. J. M., Homan, J., Keek, L., & Palmer, D. M. 2012, A&A, 547, A47

127

Israel, G. L., et al. 2005, Astrophys. J. Lett., 628, L53
Itoh, N., & Kohyama, Y. 1993, ApJ, 404, 268
Jin, R., & Kong, A. K. H. 2016, The Astronomer’s Telegram, 8530
Jonker, P. G., van der Klis, M., Homan, J., M´endez, M., Lewin, W. H. G., Wijnands, R., & Zhang,
W. 2002, MNRAS, 333, 665
Keek, L., Ballantyne, D. R., Kuulkers, E., & Strohmayer, T. E. 2014a, ApJ, 789, 121
—. 2014b, ApJ, 797, L23
Keek, L., Cumming, A., Wolf, Z., Ballantyne, D. R., Suleimanov, V. F., Kuulkers, E., &
Strohmayer, T. E. 2015, MNRAS, 454, 3559
Keek, L., & Heger, A. 2011, ApJ, 743, 189
Keek, L., Heger, A., & in’t Zand, J. J. M. 2012, ApJ, 752, 150
Keek, L., in’t Zand, J. J. M., Kuulkers, E., Cumming, A., Brown, E. F., & Suzuki, M. 2008, A&A,
479, 177
Kuulkers, E., van der Klis, M., Oosterbroek, T., Asai, K., Dotani, T., van Paradijs, J., & Lewin,
W. H. G. 1994, A&A, 289, 795
Lai, D., & Shapiro, S. L. 1991, Astrophys. J., 383, 745
Lau, K. Y. 2012, PhD thesis, Michigan State University
Lee, U. 2014, MNRAS, 442, 3037
Lense, J., & Thirring, H. 1918, Physikalische Zeitschrift, 19, 156
Levin, Y. 2006, Mon. Not. Royal Astron. Soc., 368, L35
Lewin, W. H. G., & van der Klis, M. 2006, Compact Stellar X-ray Sources
Lin, D., Remillard, R. A., & Homan, J. 2009, ApJ, 696, 1257
Link, B. 2012, Mon. Not. Royal Astron. Soc., 422, 1640
Mahmoodifar, S., & Strohmayer, T. 2013, ApJ, 773, 140
Maruyama, T., Niita, K., Oyamatsu, K., Maruyama, T., Chiba, S., & Iwamoto, A. 1998,
Phys. Rev. C, 57, 655
128

Matsumura, T., et al. 2011, The Astronomer’s Telegram, 3102, 1
McDermott, P. N., van Horn, H. M., & Hansen, C. J. 1988, ApJ, 325, 725
McDermott, P. N., van Horn, H. M., & Scholl, J. F. 1983, ApJ, 268, 837
Medin, Z., & Cumming, A. 2011, ApJ, 730, 97
—. 2014, ApJ, 783, L3
—. 2015, ApJ, 802, 29
Meisel, Z., & Deibel, A. 2017, The Astrophysical Journal, 837, 73
Meisel, Z., et al. 2015, Phys. Rev. Lett., 115, 162501
Mereghetti, S. 2008, A&A Rev., 15, 225
Migdal, A. B. 1959, Nucl. Phys. A, 13, 655
Myers, W. D., & Swiatecki, W. J. 1969, Annals of Physics, 55, 395
Nandi, R., Chatterjee, D., & Bandyopadhyay, D. 2012, arXiv:1207.3247
Negoro, H., et al. 2015, The Astronomer’s Telegram, 7943, 1
—. 2016, The Astronomer’s Telegram, 8513
Newton, W., Gearheart, M., & Li, B.-A. 2013, ApJS, 204, 9
Ogata, S., Ichimaru, S., & van Horn, H. M. 1993, ApJ, 417, 265
Oppenheimer, J. R., & Volkoff, G. M. 1939, Physical Review, 55, 374
Oyamatsu, K., & Iida, K. 2007, Phys. Rev. C, 75, 015801
Page, D., & Reddy, S. 2012, ArXiv e-prints
—. 2013, ArXiv e-prints
Paxton, B., Bildsten, L., Dotter, A., Herwig, F., Lesaffre, P., & Timmes, F. 2011, ApJS, 192, 3
Paxton, B., et al. 2013, ApJS, 208, 4
—. 2015, ApJS, 220, 15

129

Pethick, C. J., & Potekhin, A. Y. 1998, Physics Letters B, 427, 7
Pethick, C. J., & Ravenhall, D. G. 1995, Annual Review of Nuclear and Particle Science, 45, 429
Pethick, C. J., Ravenhall, D. G., & Lorenz, C. P. 1995, Nuclear Physics A, 584, 675
Philippov, A. A., Rafikov, R. R., & Stone, J. M. 2016, ApJ, 817, 62
Pines, D., & Alpar, A. 1984, in Birth and Evolution of Neutron Stars: Issues Raised by Millisecond
Pulsars, ed. S. P. Reynolds & D. R. Stinebring, 161
Piro, A. L. 2005, Astrophys. J. Lett., 634, L153
Piro, A. L., & Bildsten, L. 2005, ApJ, 619, 1054
Pons, J. A., Vigan`o, D., & Rea, N. 2013, Nature Physics, 9, 431
Potekhin, A. Y., Baiko, D. A., Haensel, P., & Yakovlev, D. G. 1999, A&A, 346, 345
Potekhin, A. Y., & Chabrier, G. 2000, Phys. Rev. E, 62, 8554
Potekhin, A. Y., Chabrier, G., & Yakovlev, D. G. 1997, A&A, 323, 415
Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 2007, Numerical recipes in
FORTRAN. The art of scientific computing, ed. Press, W. H., Teukolsky, S. A., Vetterling,
W. T., & Flannery, B. P.
Priedhorsky, W., Hasinger, G., Lewin, W. H. G., Middleditch, J., Parmar, A., Stella, L., & White,
N. 1986, ApJ, 306, L91
Rabi, I. I. 1928, Zeitschrift fur Physik, 49, 507
Ravenhall, D. G., Pethick, C. J., & Wilson, J. R. 1983, Physical Review Letters, 50, 2066
Roggero, A., & Reddy, S. 2016, ArXiv e-prints
Ruderman, M. A. 1968, Nature, 218, 1128
Russell, D. M., & Lewis, F. 2016, The Astronomer’s Telegram, 8517
Rutherford, E. 1920, Nature, 105, 500
Rutledge, R. E., Bildsten, L., Brown, E. F., Pavlov, G. G., Zavlin, V. E., & Ushomirsky, G. 2002,
ApJ, 580, 413
Samuelsson, L., & Andersson, N. 2007, Mon. Not. Royal Astron. Soc., 374, 256
130

Sato, K. 1979, Progress of Theoretical Physics, 62, 957
Schatz, H., Bildsten, L., Cumming, A., & Ouellette, M. 2003, Nuclear Physics A, 718, 247
Schatz, H., Bildsten, L., Cumming, A., & Wiescher, M. 1999, ApJ, 524, 1014
Schatz, H., et al. 1998, Phys. Rep., 294
—. 2001, Physical Review Letters, 86, 3471
—. 2014, Nature, 505, 62
Schneider, A. S., Horowitz, C. J., Hughto, J., & Berry, D. K. 2013, Phys. Rev. C, 88, 065807
Schumaker, B. L., & Thorne, K. S. 1983, Mon. Not. Royal Astron. Soc., 203, 457
Schwartz, S. J., et al. 2005, ApJ, 627, L129
Schwenk, A., Friman, B., & Brown, G. E. 2003, Nuclear Physics A, 713, 191
Shapiro, S. L., & Teukolsky, S. A. 1983, Black holes, white dwarfs, and neutron stars: The physics
of compact objects
Shternin, P. S., Yakovlev, D. G., Haensel, P., & Potekhin, A. Y. 2007, MNRAS, 382, L43
Singh, B., Rodriguez, J. L., Wong, S. S. M., & Tuli, J. K. 1998, Nuclear Data Sheets, 84, 487
Skyrme, T. 1959, Nucl. Phys. A, 9, 615
Sotani, H., Iida, K., & Oyamatsu, K. 2015, ArXiv e-prints
Sotani, H., Nakazato, K., Iida, K., & Oyamatsu, K. 2012, Phys. Rev. Lett., 108, 201101
—. 2013a, MNRAS, 428, L21
—. 2013b, MNRAS, 434, 2060
Steiner, A. W. 2008, Phys. Rev. C, 77, 035805
—. 2012, Phys. Rev. C, 85, 055804
Steiner, A. W., Lattimer, J. M., & Brown, E. F. 2010, ApJ, 722, 33
Steiner, A. W., Lattimer, J. M., & Brown, E. F. 2013, Astrophys. J. Lett., 5
Steiner, A. W., Prakash, M., Lattimer, J. M., & Ellis, P. J. 2005, Phys. Rep., 411, 325
131

Steiner, A. W., & Watts, A. L. 2009, Phys. Rev. Lett., 103, 181101
Strohmayer, T., & Mahmoodifar, S. 2014, ApJ, 784, 72
Strohmayer, T., van Horn, H. M., Ogata, S., Iyetomi, H., & Ichimaru, S. 1991, Astrophys. J., 375,
679
Strohmayer, T. E., & Brown, E. F. 2002a, ApJ, 566, 1045
—. 2002b, ApJ, 566, 1045
Strohmayer, T. E., & Markwardt, C. B. 2002, ApJ, 577, 337
Strohmayer, T. E., & Watts, A. L. 2005, Astrophys. J. Lett., 632, L111
—. 2006, Astrophys. J., 653, 593
Sugizaki, M., et al. 2013, PASJ, 65, 58
Sunyaev, R., & Kwant Team. 1989, IAU Circ., 4839
Taam, R. E., & Picklum, R. E. 1978, ApJ, 224, 210
Taam, R. E., Woosley, S. E., & Lamb, D. Q. 1996, ApJ, 459, 271
Thompson, C., & Duncan, R. C. 1995, MNRAS, 275, 255
Thorne, K. S. 1977, ApJ, 212, 825
Titarchuk, L., Seifina, E., & Shrader, C. 2014, ApJ, 789, 98
Tolman, R. C. 1934, Relativity, Thermodynamics, and Cosmology
Tsuruta, S., & Cameron, A. G. W. 1970, Ap&SS, 7, 374
Tuli, J. K. 2011, Nuclear Data Sheets
Turlione, A., Aguilera, D. N., & Pons, J. A. 2015, A&A, 577, A5
Ushomirsky, G., & Rutledge, R. E. 2001, MNRAS, 325, 1157
van der Klis, M. 1989, ARA&A, 27, 517
van der Klis, M., Swank, J. H., Zhang, W., Jahoda, K., Morgan, E. H., Lewin, W. H. G., Vaughan,
B., & van Paradijs, J. 1996, ApJ, 469, L1

132

van Hoven, M., & Levin, Y. 2012, MNRAS, 420, 3035
Ventura, J., & Potekhin, A. 2001, in The Neutron Star - Black Hole Connection, ed. C. Kouveliotou, J. Ventura, & E. van den Heuvel, 393
Wallace, R. K., & Woosley, S. E. 1981, ApJS, 45, 389
Wambach, J., Ainsworth, T. L., & Pines, D. 1993, Nuclear Physics A, 555, 128
Watanabe, G., Sato, K., Yasuoka, K., & Ebisuzaki, T. 2003, Phys. Rev. C, 68, 035806
Watts, A. L., & Strohmayer, T. E. 2006, Astrophys. J. Lett., 637, L117
Wigner, E. 1934, Physical Review, 46, 1002
Wijnands, R. 2001, ApJ, 554, L59
Wijnands, R., Guainazzi, M., van der Klis, M., & M´endez, M. 2002, ApJ, 573, L45
Wijnands, R., Homan, J., Miller, J. M., & Lewin, W. H. G. 2004, ApJ, 606, L61
Wijnands, R., Miller, J. M., Markwardt, C., Lewin, W. H. G., & van der Klis, M. 2001, ApJ, 560,
L159
Wijnands, R., Nowak, M., Miller, J. M., Homan, J., Wachter, S., & Lewin, W. H. G. 2003, ApJ,
594, 952
Wijnands, R., van der Klis, M., & Rijkhorst, E.-J. 1999, ApJ, 512, L39
Wijnands, R. A. D., van der Klis, M., Kuulkers, E., Asai, K., & Hasinger, G. 1997, A&A, 323, 399
Woosley, S. E., & Taam, R. E. 1976, Nature, 263, 101
Woosley, S. E., et al. 2004, ApJS, 151, 75
Yakovlev, D. G., Kaminker, A. D., Gnedin, O. Y., & Haensel, P. 2001, Phys. Rep., 354, 1
Yakovlev, D. G., & Urpin, V. A. 1980, Soviet Ast., 24, 303
Zamfir, M., Cumming, A., & Niquette, C. 2014, MNRAS, 445, 3278

133