UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Neutron star metallurgy Hoffman, Kelsey Llyn 2011

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.

Item Metadata

Download

Media
24-ubc_2012_spring_hoffman_kelsey.pdf [ 937.75kB ]
Metadata
JSON: 24-1.0072379.json
JSON-LD: 24-1.0072379-ld.json
RDF/XML (Pretty): 24-1.0072379-rdf.xml
RDF/JSON: 24-1.0072379-rdf.json
Turtle: 24-1.0072379-turtle.txt
N-Triples: 24-1.0072379-rdf-ntriples.txt
Original Record: 24-1.0072379-source.json
Full Text
24-1.0072379-fulltext.txt
Citation
24-1.0072379.ris

Full Text

Neutron Star Metallurgy by Kelsey Llyn Hoffman B.Sc., The University of Alberta, 2004 M.Sc., McGill University, 2006 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in The Faculty of Graduate Studies (Astronomy) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) October 2011 c© Kelsey Llyn Hoffman 2011 Abstract The crust of a neutron star plays an important role in the emission observed from it. The thermal emission generated in the core of the neutron star passes through the crust, thus it is important to know what is in the crust in order to understand how the emission is shaped and altered. The crust itself may be responsible for the observations of glitches from neutron stars and also as a source of gravitational waves. This thesis is two-fold. The first goal is to calculate the composition of the neutron star crust of a non-accreting neutron star. The second is to use the calculated crustal compositions in molecular dynamics simulations in order to determine the shear modulus and breaking strain of the crustal material. The composition of the crust is found to be dependent on how the neu- tron star cooled. Nuclear reactions within the crust are quenched as the star cools. The composition of the crust, envelope, and atmosphere are calcu- lated after the nuclear reactions are quenched. With the settling timescales of the various isotopes in the crust, some of these isotopes are able to float up to the neutron star surface and form the atmosphere. Three different cooling methods were used in these calculations – modified Urca cooling, a thick crust and a thin crust – each produces different atmospheric and crustal compositions. The calculated crustal abundances are then used as initial conditions in molecular dynamics simulations. A shear force is introduced by de- forming the simulation box. The shear modulus and breaking strain are calculated for the three different crustal compositions as well as for per- fect pure face-centered cubic (FCC) and body-centered cubic (BCC) sys- tems. The upper limit, from the perfect crystal lattice structure, on the breaking strain is found to ∼ 0.11 − 0.12 and the shear modulus is found ii Abstract to be 6.5 × 1030 dyne/cm2. These properties predict glitch amplitudes of ∆Ω/Ω∼10−3. The gravitational wave strain amplitudes for PSR J2124- 3358 are also predicted to be greater than the observed upper limits. This indicates that the neutron star crust is not a perfect BCC lattice which deformed to 10% of the maximum. iii Preface A version of Chapter 2 has been published as Hoffman, K. and Heyl, J. (2009) Compositional freeze-out of neutron star crusts Monthly Notices of the Royal Astronomical Society 400, 1986. I conducted all the simulations and analysis and wrote the majority of the manuscript, with comments from Jeremy Heyl. iv Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Preface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . v List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi Glossary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Neutron Star Structure . . . . . . . . . . . . . . . . . . . . . 1 1.2 Neutron Star Flavours . . . . . . . . . . . . . . . . . . . . . . 4 1.3 Neutron Star Atmospheres . . . . . . . . . . . . . . . . . . . 6 1.4 Crustal Cracking . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.4.1 Glitches . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.4.2 Bursts from Soft Gamma-ray Repeaters . . . . . . . . 11 1.4.3 Gravitational Waves . . . . . . . . . . . . . . . . . . . 12 1.5 Current Crustal Understandings . . . . . . . . . . . . . . . . 13 1.6 Methods Employed . . . . . . . . . . . . . . . . . . . . . . . 14 1.6.1 The Yukawa Potential . . . . . . . . . . . . . . . . . . 15 1.6.2 Nuclear Reaction Network . . . . . . . . . . . . . . . 15 v Table of Contents 1.6.3 Molecular Dynamics . . . . . . . . . . . . . . . . . . . 16 1.7 Thesis Goals . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2 Compositional Freeze-Out of Neutron Star Crusts . . . . . 18 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.2 Programs Utilized . . . . . . . . . . . . . . . . . . . . . . . . 20 2.2.1 Nuclear Statistical Equilibrium . . . . . . . . . . . . . 20 2.2.2 Nuclear Reactions . . . . . . . . . . . . . . . . . . . . 21 2.3 Cooling Curves . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.4 Crustal Mass Fraction Calculations . . . . . . . . . . . . . . 24 2.5 Timescales: Freeze-Out . . . . . . . . . . . . . . . . . . . . . 27 2.5.1 Settling Timescale . . . . . . . . . . . . . . . . . . . . 27 2.5.2 Reaction Timescale . . . . . . . . . . . . . . . . . . . 28 2.6 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.6.1 Case 1: Modified Urca . . . . . . . . . . . . . . . . . 29 2.6.2 Case 2: Thick Crust . . . . . . . . . . . . . . . . . . . 31 2.6.3 Case 3: Thin Crust . . . . . . . . . . . . . . . . . . . 34 2.7 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 3 Molecular Dynamics . . . . . . . . . . . . . . . . . . . . . . . . 39 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.2 Molecular Dynamics Simulations . . . . . . . . . . . . . . . . 39 3.3 LAMMPS . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 3.4 Characterization of the Output . . . . . . . . . . . . . . . . . 44 3.4.1 Radial Distribution Function: RDF . . . . . . . . . . 45 3.4.2 Parameters of Example Simulations . . . . . . . . . . 46 3.4.3 RDF: One Atomic Type . . . . . . . . . . . . . . . . 48 3.4.4 RDF: Two Atomic Types . . . . . . . . . . . . . . . . 51 3.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 4 Cut-Off Radius and Total Energy . . . . . . . . . . . . . . . 57 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 4.2 Simulation Set-Up . . . . . . . . . . . . . . . . . . . . . . . . 58 4.3 Results: FCC Crystal . . . . . . . . . . . . . . . . . . . . . . 58 vi Table of Contents 4.4 Results: BCC Crystal . . . . . . . . . . . . . . . . . . . . . . 58 4.5 Discussion and Conclusions . . . . . . . . . . . . . . . . . . . 60 5 The Effect of Electron Screening on Melting Temperature 64 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 5.2 Electron Screening Length . . . . . . . . . . . . . . . . . . . 66 5.3 Expected Melting: Robbins et al. (1988) . . . . . . . . . . . 67 5.4 Simulation Parameters and Calculations . . . . . . . . . . . . 68 5.5 Results: Melting Temperature . . . . . . . . . . . . . . . . . 69 5.6 Results: Simulation Length and Temperature . . . . . . . . . 71 5.7 Results: Fitting for Melting Temperature . . . . . . . . . . . 74 5.7.1 Linear Fitting . . . . . . . . . . . . . . . . . . . . . . 75 5.7.2 Dealing with Hockey Sticks . . . . . . . . . . . . . . . 78 5.8 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 6 Mechanical Properties of Pure FCC and BCC Crystals . 82 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 6.2 Applying Shear . . . . . . . . . . . . . . . . . . . . . . . . . 83 6.2.1 Deformation . . . . . . . . . . . . . . . . . . . . . . . 83 6.2.2 Simulation Parameters . . . . . . . . . . . . . . . . . 84 6.2.3 Calculation of Material Properties . . . . . . . . . . . 85 6.3 Breaking Strain of Pure Crystals . . . . . . . . . . . . . . . . 86 6.3.1 Box Size Effects . . . . . . . . . . . . . . . . . . . . . 87 6.3.2 Pure FCC Crystal . . . . . . . . . . . . . . . . . . . . 88 6.3.3 Pure BCC Crystal . . . . . . . . . . . . . . . . . . . . 88 6.4 Second Break . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 6.5 Forward and Reverse . . . . . . . . . . . . . . . . . . . . . . 95 6.6 Temperature Dependence . . . . . . . . . . . . . . . . . . . . 98 6.7 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 7 Mechanical Properties of Neutron Star Crusts . . . . . . . 105 7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 7.2 Simulation Set-up . . . . . . . . . . . . . . . . . . . . . . . . 106 7.3 Melting Simulations . . . . . . . . . . . . . . . . . . . . . . . 109 vii Table of Contents 7.4 Stress-Strain Relationships . . . . . . . . . . . . . . . . . . . 112 7.5 Second Break . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 7.6 Forward and Reverse . . . . . . . . . . . . . . . . . . . . . . 124 7.7 Temperature Dependence . . . . . . . . . . . . . . . . . . . . 126 7.8 Discussion and Conclusions . . . . . . . . . . . . . . . . . . . 131 8 Application of Breaking Strain . . . . . . . . . . . . . . . . . 134 8.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 8.2 Glitches . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135 8.3 Soft Gamma-ray Repeaters . . . . . . . . . . . . . . . . . . . 137 8.4 Gravitational Waves . . . . . . . . . . . . . . . . . . . . . . . 138 8.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141 9 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143 9.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143 9.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . 145 9.2.1 Atmosphere Thermal Evolution . . . . . . . . . . . . 145 9.2.2 Bursts from Soft Gamma-ray Repeaters . . . . . . . . 146 9.2.3 Planetary Interiors . . . . . . . . . . . . . . . . . . . 146 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148 viii List of Tables 2.1 Summary of the atmosphere and crust compositions . . . . . 30 3.1 Unit conversions for the reduced units . . . . . . . . . . . . . 44 3.2 Size of the simulation boxes . . . . . . . . . . . . . . . . . . . 45 3.3 Pair coefficients for the RDF example simulations . . . . . . . 48 4.1 Total energy with cut-off radius for a FCC crystal . . . . . . 59 4.2 Total energy with cut-off radius for a BCC crystal . . . . . . 60 5.1 Melting temperatures for the various cut-off radii . . . . . . . 74 5.2 Melting temperatures and simulation length . . . . . . . . . . 75 5.3 Linear fits with no errors . . . . . . . . . . . . . . . . . . . . 77 5.4 Linear fits using two data points . . . . . . . . . . . . . . . . 78 6.1 Breaking strain for the FCC crystal at various strain rates . . 93 6.2 Breaking strain for the BCC crystal at various strain rates . . 93 6.3 Temperature dependence for FCC and BCC crystals . . . . . 103 7.1 Isotope number fractions . . . . . . . . . . . . . . . . . . . . . 109 7.2 Parameters for the impure simulations . . . . . . . . . . . . . 110 7.3 Pair coefficients for the modified Urca composition . . . . . . 110 7.4 Pair coefficients for the thick crust composition . . . . . . . . 111 7.5 Pair coefficients for the thin crust composition . . . . . . . . 111 7.6 Breaking strain: modified Urca composition . . . . . . . . . . 119 7.7 Breaking strain: thick crust composition . . . . . . . . . . . . 121 7.8 Breaking strain: thin crust composition . . . . . . . . . . . . 121 7.9 Temperature dependence of FCC and BCC crystals . . . . . . 131 ix List of Tables 8.1 Gravitational wave strain amplitudes . . . . . . . . . . . . . . 140 x List of Figures 1.1 Diagram of the regions of a neutron star . . . . . . . . . . . . 3 1.2 Diagram of a pulsar glitch . . . . . . . . . . . . . . . . . . . . 9 2.1 Three cooling curves . . . . . . . . . . . . . . . . . . . . . . . 25 2.2 Lightest isotopes for modified Urca cooling . . . . . . . . . . 32 2.3 Freeze-out for modified Urca cooling . . . . . . . . . . . . . . 33 2.4 Lightest isotopes for a thick crust . . . . . . . . . . . . . . . . 35 2.5 Lightest isotopes for a thin crust . . . . . . . . . . . . . . . . 36 3.1 Diagram of neighbour list construction . . . . . . . . . . . . . 42 3.2 RDF for FCC and BCC crystals . . . . . . . . . . . . . . . . 49 3.3 RDF for the melted FCC and BCC crystals . . . . . . . . . . 50 3.4 RDF for a FCC crystal with two types of atoms . . . . . . . . 53 3.5 RDF for a BCC crystal with two types of atoms . . . . . . . 54 3.6 RDF for a melted FCC crystal with two types of atoms . . . 55 3.7 RDF for a melting BCC crystal with two types of atoms . . . 56 4.1 Determining the cut-off radius required for a FCC crystal . . 62 4.2 Determining the cut-off radius required for a BCC crystal . . 63 5.1 Example melting temperature determination . . . . . . . . . 70 5.2 Energy ratio as a function of the screening parameter . . . . 72 5.3 Melting temperature as a function of the screening parameter 73 5.4 Melting temperature and cut-off radius regimes . . . . . . . . 76 5.5 Comparing the calculated to expected melting temperatures . 79 6.1 Illustration of box deformation and strain calculation . . . . . 84 xi List of Figures 6.2 A schematic of a stress-strain relationship . . . . . . . . . . . 87 6.3 Breaking strain and the size effects of a FCC crystal . . . . . 89 6.4 Breaking strain and the size effects of a BCC crystal . . . . . 90 6.5 Breaking strain of a FCC crystal . . . . . . . . . . . . . . . . 91 6.6 Total energy with FCC crystal deformation . . . . . . . . . . 92 6.7 Breaking strain of a BCC crystal . . . . . . . . . . . . . . . . 94 6.8 Investigating a second break of a FCC crystal . . . . . . . . . 96 6.9 Investigating a second break for a BCC crystal . . . . . . . . 97 6.10 Reversing the strain applied to a FCC crystal . . . . . . . . . 99 6.11 Reversing the strain applied to a BCC crystal . . . . . . . . . 100 6.12 Temperature dependence on the breaking strain for a FCC crystal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 6.13 Temperature dependence on the breaking strain for a BCC crystal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 7.1 RDF of the modified Urca composition for all particles . . . . 113 7.2 RDF of the modified Urca composition for each isotope . . . 114 7.3 RDF of the thick crust composition for all atoms . . . . . . . 115 7.4 RDF of the thick crust for each isotope . . . . . . . . . . . . 116 7.5 RDF of the thin crust for all atoms . . . . . . . . . . . . . . . 117 7.6 RDF of the thin crust for each isotope . . . . . . . . . . . . . 118 7.7 Breaking strain for the impure lattices at ṡ = 20× 10−6 . . . 120 7.8 Comparing the perfect and imperfect thick crust crystals . . . 122 7.9 Breaking strain of a thick crust imperfect crystal at ṡ = 20× 10−6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 7.10 Investigating a second break for the impure crusts . . . . . . 125 7.11 Reversing the strain applied to the impure crusts . . . . . . . 127 7.12 Temperature dependence of the breaking strain of the modi- fied Urca composition . . . . . . . . . . . . . . . . . . . . . . 128 7.13 Temperature dependence of the breaking strain of the thick crust composition . . . . . . . . . . . . . . . . . . . . . . . . . 129 7.14 Temperature dependence of the breaking strain of the thin crust composition . . . . . . . . . . . . . . . . . . . . . . . . . 130 xii Glossary Name Description AXP Anomalous X-ray Pulsar BCC Body-Centerd Cubic FCC Face-Centered Cubic LAMMPS Large Atomic/Molecular Massively Parallel Simulator RDF Radial Distribution Function SGR Soft Gamma-ray Repeater M Mass of the Sun Ω Angular Frequency ρcore Core Density φm Maximum Breaking Strain NVE Constant Number Volume and Energy µ Shear Modulus Qmax22 Maximum Quadrupole κ Inverse Screening Length α Fine Structure Constant λ Screening Length Parameter a Characteristic Length Uo Characteristic Energy n Number Density rc Cut-Off Radius (Ze)2 Pair Coefficient ṡ Strain Rate h Gravitational Wave Strain Amplitude xiii Acknowledgements I would first like to thank my supervisor Jeremy Heyl for suggesting this project and for helping make my time at UBC an enjoyable experience. I would also like to acknowledge my supervisory committee, Harvey Richer, Jorge Rottler, and Ingrid Stairs, for all the suggestions along the way. I would like to thank my examiners Chad Sinclair, Jaymie Matthews, and Charles Horowitz, for their comments and suggestions during the examina- tion process. I would like to thank my office mates, past and present, Jaime, Nicole, Laura, and Sam for joining in on all the tea runs. An extra shout out goes to Jaime for all her help near the end of this process. Thanks to Aaron, Chris, Cindy, Marjorie, Thomas, and Tyron for always being ready to take a break and file something in the mullet drawer. A huge thank you to Jason Rowe for supporting me along the way and always being ready with words of encouragement and a mai tai. I would finally like to thank my family, Nina, Marvin, Ken, Donna, and Lori, for all their support. xiv Dedication For my mother. xv Chapter 1 Introduction 1.1 Neutron Star Structure Metallurgy is the study of metals and their properties. The periodic table has many metals, including iron. A handful of iron, a volume of 60 cm3, weighs less than half a kilogram. Taking the same volume of neutron star material in one hand would be a bone-breaking mass of over 1010 tonnes! Neutron stars are objects which are all about extremes and in this thesis the material properties of neutron stars are examined. Neutron stars are born in the aftermath of a Type ii supernova, resulting from the core collapse of a star with a mass greater than 8 M1 (e.g., Lat- timer & Prakash, 2004). The upper limit on the mass of the progenitor star that would form a neutron star after its nuclear fuel is consumed depends on the metallicity of the star; for a metal-poor star the upper limit is 25 M. This upper limit increases with metallicity (Heger et al., 2003). The newly formed neutron star is hot, with temperatures in the core reaching in excess of 1010 K (Page et al., 2006). After the neutron star is formed, it quickly cools by neutrino emission and after about one day the temperature drops to 109 − 1010 K, with neutrino cooling dominating for the next ∼103 years (e.g., Shapiro & Teukolsky, 1986). Neutron stars are extremely dense objects, with masses on the order of ∼1.5 M and radii of ∼12 km (Lattimer & Prakash, 2004). The average neutron star density is ∼4 × 1014 g cm−3, greater than the density of an atomic nucleus, 2.3 × 1014 g cm−3. The masses of some neutron stars have been determined from binary pulsar systems. The most precisely measured neutron star masses range from 1.25 to 1.44M (Thorsett & Chakrabarty, 1One M is the mass of the sun, 1.9889×1033 g. 1 1.1. Neutron Star Structure 1999). Recently, the highest neutron star mass has been measured to be 1.97 ± 0.04 M in the binary millisecond pulsar system PSR J1614-2230 (Demorest et al., 2010). The radii of neutron stars are harder to measure. Observing the thermal spectra of neutron stars results in a measurement of the gravitational red- shift, which depends on both the radius and the mass of the star. With one observation the mass and radius can not become untangled. In recent work, different spectroscopic observations of the same source have been used to independently determine its mass and radius (Özel et al., 2010a). Neutron stars are not only extreme in terms of their density, they also have extreme magnetic fields, reaching more than 1015 Gauss (Kaspi, 2010). For this discussion the structure of a neutron star is divided into six regions: inner core, outer core, inner crust, outer crust, envelope, and at- mosphere. These regions are displayed schematically in Figure 1.1 and are briefly described below. Inner Core: (ρcore) The majority of the neutron star mass is found within the core (Lattimer & Prakash, 2004). The behaviour of the inner core is not well understood as it is at super-nuclear densities (Lattimer & Prakash, 2007). The inner core is thought to be composed of exotic particles such as pions or kaons, or even strange quark matter (Lat- timer & Prakash, 2004). The mass measurement of PSR J1614-2230 indicates that a transition to quark matter can occur within the core (Özel et al., 2010b). Outer Core: (∼2×1014 g cm−3 ≤ ρ ≤ ρcore) The outer core is a superfluid and superconductor containing a mixture of nucleons, electrons and muons (Lattimer & Prakash, 2004). Inner Crust: (4.3 × 1011 g cm−3 ≤ ρ ≤ 2 × 1014 g cm−3) The inner and outer crust together extend 1-2 km below the neutron star surface (Lattimer & Prakash, 2004). The inner crust is bounded on the lower density side by neutron drip, the point at which neutrons leaked out of the atomic nuclei, at ρ = 4.3× 1011 g cm−3. Electrons, protons, nuclei 2 1.1. Neutron Star Structure Inner Core Outer Core Surface: Atmosphere and Envelope Outer Crust Inner Crust Neutron Drip ρ=106 ~10km 2×10 ~100m ~1km 11 14 4.3×10 Figure 1.1: A schematic of the different regions of a neutron star. The density of the regions, in units of g/cm3, are listed on the left side of the diagram and the boundaries of the regions are labelled on the right. The atmosphere and envelope contain a small fraction of the neutron star mass. The outer crust is a lattice of nuclei. The inner crust begins at neutron drip with free neutrons increasing in number as the density increases. The outer core is a superfluid and superconductor. Details of the inner core are unknown. 3 1.2. Neutron Star Flavours and free neutrons exist together within the inner crust. As the density increases, so does the number of free neutrons (Shapiro & Teukolsky, 1986). Outer Crust: (106 g cm−3 ≤ ρ ≤ 4.3 × 1011 g cm−3) The outer crust is a solid lattice of nuclei, without any free neutrons (Shapiro & Teukolsky, 1986). In the cold catalyzed matter hypothesis the outer crust is a body-centered cubic (BCC) of iron, 56Fe (Chamel & Haensel, 2008). Envelope: The envelope is involved with thermal transport of energy (Lat- timer & Prakash, 2004) and is also referred to as the neutron star ocean. This layer of the neutron star is in the liquid phase and is composed of electrons and nuclei (Potekhin, 2010). Atmosphere: The atmosphere shapes the emerging thermal spectrum (Lat- timer & Prakash, 2004). The bottom of the atmosphere is where the optical thickness is close to unity (Potekhin, 2010). The composition of the atmosphere depends on whether the neutron star has gone through any accretion events. More on the expected composition of neutron star atmosphere is discussed in Section 1.3. It is the upper regions of the neutron star – the atmosphere, the envelope, and especially the outer crust – which are the focus of this thesis, with regards to the composition of each layer, as well as the mechanical properties of the crust. 1.2 Neutron Star Flavours Neutron stars were first theorized in 1934 as the transition from an ordinary star via a supernova explosion (Baade & Zwicky, 1934). The high temper- ature of their formation led to the expectation that neutron stars would be observed in the form of X-ray sources (Chiu & Salpeter, 1964). However, it was in the form of radio pulsars that neutron stars were first observed in 1967 (Hewish et al., 1968). Since these first observations many different 4 1.2. Neutron Star Flavours flavours of neutron stars have subsequently been observed. These differ- ent types of neutron stars have been reviewed by Kaspi (2010), Mereghetti (2001), and Popov (2008). The focus of this thesis is on young neutron stars. The types of stars that will be discussed include radio pulsars, magnetars – which includes anomalous X-ray pulsars (AXPs) and soft gamma-ray re- peaters (SGRs) – X-ray dim isolated neutron stars, and strange or quark matter stars. Each type is summarized below. Radio Pulsars: These neutron stars are observed at radio frequencies. They have regular periodic pulses, with the prototypical radio pulsar being the Crab pulsar (Kaspi, 2010). Pulsars are observed to slow their rotation period over time. This ‘spin-down’ is a result of the magnetic field slowing down the rotation, releasing energy (Kaspi, 2010). Magnetars: Magnetars are slowly rotating neutron stars with extremely strong magnetic fields, B > 1014G (Woods & Thompson, 2006). The high magnetic field is the source of the observed emission (Mereghetti, 2008). Magnetars can be divided into two subclasses: Anomalous X-ray Pulsars (AXPs) and Soft Gamma-ray Repeaters (SGRs). The AXPs are named for their pulsed X-ray emission. The SGRs are known for their short bursts, but in three sources giant flares have also been observed. These giant flares emit more energy than the short bursts, and after this massive release of energy the light curves display a decay into a softer energy tail with pulsations (Woods & Thompson, 2006). X-ray Dim Isolated Neutron Stars: There are seven confirmed X-ray dim isolated neutron stars. These neutron stars are young and, be- cause of their low luminosity, they are nearby with distances of ∼ 500 pc, or closer (Kaplan, 2008). X-ray dim isolated neutron stars have thermal X-ray spectra, with no indications of a non-thermal compo- nent, and they do not have radio counterparts (Kaspi, 2010). Strange Stars: Strange stars may form if at high pressures there is a phase transition from normal matter to quark matter(Alcock & Olinto, 1988). Unlike normal neutron stars, where the equation of state has 5 1.3. Neutron Star Atmospheres pressure vanish as the density goes to zero, strange stars, like most solid matter, are self-bound and the pressure vanishes at a finite den- sity (Lattimer & Prakash, 2001). The surface of a strange star could either be bare or have a thin crust of material which is the same as the outer crust of a neutron star, where there are no free neutrons (Alcock & Olinto, 1988). Currently, there is no strong evidence for the existence of strange stars. The composition of the upper regions of the neutron star is calculated in this thesis for the case where the neutron star has not undergone any accretion events, from either fallback or the interstellar medium. This is an unphysical situation whose closest approximation is a X-ray dim isolated neutron star. The compositions are calculated using various cooling curves, including a cooling curve appropriate for strange matter stars. Though the case of a completely non-accreting neutron star is strictly speaking unphys- ical, the mechanical properties of such a crust can be used as a limiting case when considering radio pulsar glitches, SGR giant flares and gravitational wave emission. 1.3 Neutron Star Atmospheres The thermal emission of a neutron star is embedded with the signature of the composition of the atmosphere. In order to understand the thermal emission, these observations have been fit with black body models. Analysis of the fits have indicated that the observations deviate from a pure black body spectrum (Zavlin & Pavlov, 2002). Subsequently, atmosphere compositions such as hydrogen, carbon, as well as silicon have been used to model the observed thermal emission. The composition of the neutron star atmosphere depends on the details of accretion, as well as nuclear reactions occurring on the surface. In a situation where there has not been many accretion events, an iron atmosphere would be expected for a young hot neutron star (Rajagopal et al., 1997). After a neutron star is formed, light elements such as hydrogen 6 1.4. Crustal Cracking and helium may be found on the surface. These lighter elements could then undergo nuclear reactions, depleting their abundance. Thus, without additional accretion or convective mixing, an iron atmosphere would be produced (Chiu & Salpeter, 1964). Due to gravitation settling, the atmosphere of a neutron star which has undergone accretion would be composed of the lightest isotope accreted (Za- vlin & Pavlov, 2002). For example, as hydrogen is the lightest element, any accreted mixture containing hydrogen, to an abundance for it to be optically thick, would result in a hydrogen atmosphere. The neutron star atmosphere can change with time as nuclear reactions occur. The hydrogen in the at- mosphere can diffuse to a lower depth and become depleted due to proton capture reactions (Chang et al., 2010). This diffusive nuclear burning could lead to a helium atmosphere, but recently a neutron star with a carbon at- mosphere has been observed (Ho & Heinke, 2009), indicating that helium can be consumed in a similar process as hydrogen. Heavier atmospheres have been examined, such as silicon models. Silicon atmosphere models have been used to describe the thermal spectra of the isolated neutron star RX J185635-3754 (Pons et al., 2002; Walter et al., 2004). This silicon at- mosphere could result from nuclear reactions within the crust, as will be explored in Chapter 2. 1.4 Crustal Cracking The outer crust of a neutron star is at subnuclear densities and as a result is composed of a lattice of nuclei. For this study of the properties of the neutron star crust, magnetic fields are not applied and accretion is ignored; material with these properties is known as cold catalyzed matter. In the case of cold catalyzed matter the outer crust would be a body-centered cubic (BCC) lattice of 56Fe (Chamel & Haensel, 2008). Stresses, either due to rotation effects or the magnetic field, applied to the neutron star crust could induce crustal cracking events, which are known as ‘starquakes’, or ‘crustquakes’. Observational phenomena which have been attributed to starquakes as a 7 1.4. Crustal Cracking possible explanation include pulsar glitches and SGR bursts. The properties of the crust are also important for the emission of gravitational waves from the neutron star. Observations of events associated with starquakes have been found to share similar statistical properties to terrestrial earthquakes (Cheng et al., 1996). In order to determine the role that crustal cracking may play on these observed phenomena the shear modulus and breaking strain of the crustal material must be determined. The breaking strain indicates the degree of deformation at which the material yields. The shear modulus is the ratio of the applied shear stress to the degree of shear strain in the linear regime. The shear modulus is also dependent on the direction in which the shear stress is applied. The mechanical properties of the crust not only have implications for starquakes but also for the possibility of mountains forming on the neutron star. These mountains on the surface could be a source of gravitational waves. The neutron star phenomena of glitches, SGR bursts, as well as gravitational waves are discussed in Sections 1.4.1, 1.4.2, and 1.4.3, respectively. 1.4.1 Glitches Glitches are phenomena that have been observed in timing observations of radios pulsars and have also been observed in anomalous X-ray pulsars (AXPs). The rotation of neutron stars slows down due to magnet dipole emission. A glitch is an event where there is a sudden increase in the rotation frequency and then the period decays back toward the pre-glitch rotation period, but the period derivative changes (Lyne et al., 2000). A schematic of an idealized glitch is shown in Figure 1.2. In the case of AXP glitches, the glitch has on occasion been observed to coincide with a burst event (Woods & Thompson, 2006). Two models which have been used to explain glitches are the starquake model and the other is the superfluid vortex model (Chamel & Haensel, 2008). These two models are discussed below. In the starquake model the glitches are due to a change in the moment of inertia of the neutron star resulting from the crust cracking and readjusting its shape (Chamel & Haensel, 2008). The discussion of starquakes as applied 8 1.4. Crustal Cracking Glitch Time R ot at io n Fr eq ue nc y Figure 1.2: A schematic of a pulsar glitch. The rotation frequency of a neutron star typically decreases with time. The glitch is the point where there is a sudden increase the rotation frequency, but eventually the neutron star returns to the pre-glitch rotation frequency. 9 1.4. Crustal Cracking to glitches follows Ruderman (1969). The strain, at the point of cracking, of the neutron star surface is given by: φ = 7R5 8GM (ω2I − ω2) sin2 2θ, (1.1) where ωI is the initial angular velocity, ω is the current angular velocity, θ is the polar angle, R is the radius and M is the mass of the neutron star. As the neutron star rotation slows, the shape of the star changes from an oblate to a more spherical shape. This induces a stress on the crust until the breaking strain is reached. By approximating the interior of the neutron star as having constant density the change in the surface of the neutron star as the maximum breaking strain, φm, is reached due to the rigidity, or shear modulus of the crust, is given by: δR(θ) R ∼ 95µφmR 7GMρ [ 1− ( Ri R )7] (1− 3 cos2 θ), (1.2) where Ri is the inner radius of the crust, µ is the shear modulus, and θ is the polar angle. After the crust yields to the applied stress there is a decrease in the moment of inertia, and thus an increase in the angular frequency. This change in the angular frequency, ∆Ω, is related to the change of the moment of inertia through: ∆Ω Ω = −∆I I = 2 δR(θ = pi/2) R , (1.3) where ∆I is the change in the moment of inertia. With the starquake model of neutron star glitches, it is possible to predict the time between subsequent glitches. The glitch amplitude is related to the time between glitches as: ∆Ω Ω <  δt τ , (1.4) where  is proportional to Ω2, τ = Ω̇/2Ω, and δt is the time between glitches (Chamel & Haensel, 2008). The starquake model has been found to be con- sistent with some observations, for example the Crab pulsar, but it fails to 10 1.4. Crustal Cracking predict the time interval of glitches for the Vela pulsar (Shapiro & Teukolsky, 1986). Another interpretation of the pulsar glitches considers not only the solid outer crust region, but also the superfluid interior of the neutron star, the superfluid vortex model. In this interpretation of pulsar glitches angular momentum is transferred between these two regions (Chamel & Haensel, 2008). The following discussion follows that in Lyne et al. (2000). The superfluid within the neutron star is rotating in quantized vortices which are pinned to the crust. As the rotation of the neutron star changes with time, the superfluid region rotates at a different rate than the crust. With these two regions rotating at different rates, the pinned vortices become unpinned to the crust. This unpinning releases angular momentum to the crust, resulting in a glitch. It is also possible that the unpinning of the vortices could lead to the crust cracking. 1.4.2 Bursts from Soft Gamma-ray Repeaters The SGRs are classified by the observations of short bursts of low-energy gamma-rays which can reach a peak luminosity of 1041 erg s−1 (Woods & Thompson, 2006). There have been eleven observed SGRs, 7 confirmed and 4 candidates2. Three of these 11 sources have been observed to produce giant flares. The giant flares are events which emit more energy than the short bursts and are observed to decay into a softer-energy tail which has pulsations (Woods & Thompson, 2006). The first giant flare was observed on March 5, 1979 from SGR 0526-66 (Mazets et al., 1979). This giant flare had a peak luminosity of around 5× 1044 erg s−1(Helfand & Long, 1979). In the decay of the burst, oscillations with a periodicity of ∼8 s were discovered in the gamma-ray emission (Terrell et al., 1980). The second giant flare to be observed was from SGR 1900+14 on August 27, 1998; this flare had a periodicity of 5.126 s in the decay tail (Hurley et al., 1999). The giant flare on December 27, 2004 from SGR 1806-20 was the third giant flare observed and at its peak was ∼100 times more energetic than the August 1998 flare 2Data from the McGill SGR/AXP online catalog, http://www.physics.mcgill.ca/ pul- sar/magnetar/main.html 11 1.4. Crustal Cracking and the tail had a periodicity of 7.56 s (Hurley et al., 2005). Both the giant flares with the subsequent oscillations and the bursts are believed to be associated with crustal events. These bursts from SGRs can be interpreted as being due to stresses placed on the crust from the magnetic field inducing crustal cracking (Chamel & Haensel, 2008). The magnetic field required to produce a burst is related to the breaking strain as, B = 1015 ( ESGR 1041erg )−1/2 ( l 1km )( φm 10−3 ) G, (1.5) where ESGR is the energy of the burst and l is the length of the fracture (Thompson & Duncan, 1995). As on Earth, such a large event on the surface would likely cause vibrations to propagate throughout the star (Chamel & Haensel, 2008). 1.4.3 Gravitational Waves A neutron star symmetric about its rotation axis does not radiate gravita- tional waves, but a deformation on the surface of the neutron star would result in the emission of gravitational waves (Chamel & Haensel, 2008). An- other signature of an asymmetrically shaped neutron star would be in the free precession of the star (Stairs et al., 2000). In the context of the me- chanical properties of the crust, the emission of gravitational waves due to a deformation on the crust would be the focus of this study. The amplitude of the quadrupole moment, Qmax22 , depends on the breaking strain of the crust, Qmax22 ≈ 1038g cm2 ( φm 10−2 ) , (1.6) where φm is the breaking strain of the crust (Haskell et al., 2006). The gravitational waves emitted from the neutron star would be observed at Earth to have a strain amplitude of h = 16 5 ( pi 3 )1/2 GQ22 Ω2 d c4 , (1.7) 12 1.5. Current Crustal Understandings where Ω is the angular frequency and d is the distance to the neutron star (Ushomirsky et al., 2000). The strain amplitude of a neutron star can be predicted based on the mechanical properties of the crust. Observation up- per limits have been placed on the gravitational wave strain amplitude of 78 radio pulsars using the Laser Interferometer Gravitational Wave Observa- tory (LIGO) (Abbott et al., 2007). A detection of gravitational waves from a neutron star, with a comparison to predicted strain amplitudes, would put constraints on the strength and structure of the crust. 1.5 Current Crustal Understandings The strength of the crust has been examined using various methods: ex- amining the conditions required for shear waves, using molecular dynamics simulations and other methods including perturbation theory. The break- ing strain of a neutron star has been predicted to range between φm = 10−5 and φm = 10−2, and the shear modulus has been predicted to be µ = 1030dyne cm−2 (Smoluchowski, 1970). The breaking strain prediction in Smoluchowski (1970) is based on arguments that the neutron star crust would be expected to be impure in composition and also contain various defects. As a result, the crust would be weaker than theoretical calculations of perfect crystals of less dense, terrestrial matter (see Kittel, 1976). The shear modulus prediction results from treating the crust as a Coulomb lattice which has the approximate relationship of µ ∼ (Ze)2n4/3, where the charge Z ranges between 30 and 50, and n is the number density (Smoluchowski, 1970). In order to understand what happens to the crust at the point in which it yields to the applied stresses, Jones (2003) compared the nature of the neutron star fault planes to that of tectonic earthquakes. In Jones (2003) the shear modulus was assumed to be µ = 1030 dyne cm−2, near neutron drip for a BCC lattice. The pressure exerted by the neutron star is found to be orders of magnitude larger then the shear modulus, and as a result the crust can not support a void for a long enough time-span. Jones (2003) finds that brittle fractures cause by Maxwell stresses are not expected to 13 1.6. Methods Employed occur in the neutron star crust, though this work was calculated for a case of a pure crust. Recently, molecular dynamics simulations have been performed in order to understand the crust of an accreting neutron star, which would contain impurities. The composition of an accreting crust was calculated by Gupta et al. (2007). As a result of these simulations, it was found that chemical separation occurs in the crustal regions, with a liquid ocean being enhanced with low atomic elements compared to the solid crust. The chemical sep- aration could then have an impact on the shear modulus (Horowitz et al., 2007). This chemical phase separation analysis was followed up with semi- analytic Monte Carlo calculations that were able to reproduce the same chemical compositions in the liquid and solid regions of the crust (Medin & Cumming, 2010). Molecular dynamics simulations have also been used to calculate the breaking strain of an accreted neutron star crust in Horowitz & Kadau (2009). The breaking strain of this type of crust was found to be close to φm = 0.1 and could support a mountain on the neutron star which would produce gravitational waves detectable by LIGO. It has also been found by using molecular dynamics simulations that the breaking stress of the crust is dependent on the temperature of the simulation (Chugunov & Horowitz, 2010) The shear modulus has also been calculated using thermodynamic per- turbation theory, as in Baiko (2011). With this type of calculation, both the classical and quantum regimes of ion motion are considered, whereas molec- ular dynamics simulations only consider the classical regime. The quantum effects are found to only become important for light elements which are at high densities. 1.6 Methods Employed In order to investigate the properties of the neutron star crust, various meth- ods have been employed in this thesis. The isotopes in the crust interact with each other via a Yukawa, or screened Coulomb, potential. A nuclear reaction 14 1.6. Methods Employed network was used in the calculation of the isotopes of the crustal compo- sition. In order to investigate the mechanical crustal properties, molecular dynamics simulations were employed. The Yukawa potential, nuclear reac- tion networks and molecular dynamics are briefly discussed below. 1.6.1 The Yukawa Potential Due to the high densities and pressures of the neutrons in the crust, the ions are pressure ionized and interact via a screened Coulomb or Yukawa poten- tial. The electrons are relativistic and degenerate. The Yukawa potential is given as: Vij = ZiZke 2 r e−κr, (1.8) where Zi and Zj are the ion charges, r is the pair separation and κ is the inverse electron screening length. The inverse screening length can be given by the inverse Fermi screening length, κ = 2α1/2kF /pi 1/2, where kF is the electron Fermi momentum and α is the fine structure constant. The inverse screening length is parameterized by introducing the parameters λ = κa, where a = n−1/3 is the characteristic particle spacing. Yukawa systems have been studied extensively and the work of Robbins et al. (1988) will be referred to frequently in this work. In Robbins et al. (1988) the phase diagram and dynamics of a Yukawa system were studied using molecular dynamics. The study explored the melting transition of the Yukawa system, as well as the transition from face-centered cubic (FCC) to body-centered cubic (BCC) crystal phases. The melting transition in Robbins et al. (1988) will be compared to current calculations in Chapter 5. 1.6.2 Nuclear Reaction Network In Chapter 2, a nuclear reaction network is used to calculate crustal compos- tions. Such nuclear reaction networks are discussed in Hix & Meyer (2006). A nuclear reaction network is composed of a system of first-order differen- tial equations. The differential equations describe the abundance of each isotope considered in the network. Sparse matrices are often used in these 15 1.7. Thesis Goals types of calculations in order to decrease memory usage and computational time. The sparse matrices are a result of neglecting reactions between cer- tain species and instead focusing on those reactions which occur at faster rates, such as the capture or release of neutrons, protons, α-particles, or photons. These differential equations are integrated with the isotope abun- dances updated for the next time-step of the calculation. 1.6.3 Molecular Dynamics This thesis uses molecular dynamics simulation to investigate the mechani- cal properties of the crust. The simulations are performed with the software LAMMPS, the Large Atomic/Molecular Massively Parallel Simulator, which is discussed in Chapter 3. The technique of molecular dynamics is discussed in depth in Allen (2004). The term, molecular dynamics, refers to solving Newton’s equations of motion for a system of particles. The particles in the systems of interest for this thesis interact via the Yukawa potential. With an initial configuration of particle positions the particle motions are calculated with a velocity-Verlet algorithm (Swope et al., 1982). For each particle, the position, velocity, and acceleration is calculated, and with this data the thermodynamic properties of the system can be determined. The molecular dynamics simulations completed in this work make use of neighbour lists in order to decrease the computational time for each calculation. Instead of calculating the potential between all the particles in the system, the neigh- bour list identifies the particles within a specified distance from the particle of interest. It is only the particles within the neighbour list for which the potential is calculated. The parameters for the cut-off length used for the neighbour list are examined in Chapter 4. 1.7 Thesis Goals The main focus of this thesis is to investigate the properties of a non- accreting neutron star crust. The first part of this investigation of neutron star properties is the calculation of the composition of the non-accreting 16 1.7. Thesis Goals crust. This calculation of the crust composition is discussed in Chapter 2, which also discusses the composition of the atmosphere. These crustal abundances are used with molecular dynamics for an investigation of the me- chanical properties of the crust. The molecular dynamics simulation package LAMMPS and the characterization of the output are discussed in Chapter 3. The cut-off radius used in the simulations is covered in Chapter 4. The phase transition of a Yukawa system is calculated in Chapter 5 and com- pared to that of Robbins et al. (1988). Before performing the molecular dynamics simulations on the crustal compositions, the shear modulus and breaking strain of a pure FCC and BCC crystal are calculated in Chapter 6. The same types of simulations are then applied to the crustal compositions in Chapter 7. The results of the impure crust calculations are compared to observations linked to crustal cracking in Chapter 8 in order to determine if fracturing is a workable model for some of the observed phenomena in neutron stars. Finally Chapter 9 summarizes and concludes what has been found in this thesis, and what the future may bring. 17 Chapter 2 Compositional Freeze-Out of Neutron Star Crusts3 2.1 Introduction In order to determine the detailed composition of the envelope and atmo- sphere, the chemical evolution as nuclear reactions abate are calculated for three different cases of neutron star crusts. For all of these calculations, the crustal properties are investigated for neutron stars without fallback accre- tion. The three cases are distinguished by the method of cooling: neutron star cooling via a modified Urca process without the thermal influence of the crust, and direct Urca process with both a thick crust and a thin crust. The Urca process is a method of neutron star cooling via neutrino emission due to beta decay (Shapiro & Teukolsky, 1986). The neutron star crust extends 1-2 km below the surface inward to a density of around 1014 g/cm3. The neutron star envelope is defined to be the upper layer of the crust that throttles the heat flow running from a density of 107 g/cm3 outward (Hernquist & Applegate, 1984; Heyl & Hern- quist, 2001). The atmosphere lies at relatively a low density and comprises a column density of about 1 g/cm2. Though the envelope and atmosphere contain a negligible amount of the neutron star mass (or the mass of the crust for that matter), they both play a crucial role in shaping the observed properties of neutron stars. In order to fully interpret the observed emission, it is important to understand the crustal composition, especially the lightest trace elements that can float to the neutron star surface and form the en- 3A version of this chapter has been published. Hoffman, K. & Heyl, J. (2009) Compo- sitional freeze-out of neutron star crusts MNRAS 400, 1986. 18 2.1. Introduction velope and atmosphere. As the atmosphere shapes the emergent spectrum, understanding the composition could yield predictions of spectral features in the neutron star thermal emission. Since the envelope influences the trans- port and release of thermal energy, a change in its composition also changes the thermal conductivity and inferred surface temperature of the neutron star (Lattimer & Prakash, 2004). Trace light elements among the calculated isotopes may have enough time to float to the surface of the neutron star before the layer crystallizes and form the atmosphere or envelope. Without fallback from the supernova or other accretion, a neutron star is expected to have an atmosphere of iron-group elements (Chiu & Salpeter, 1964). An example of a family of neutron stars which are closest to this idealized situation are the X-ray dim isolated neutron stars. The isolated neutron star RX J185635-3754 has been studied extensively. Its surface composition was estimated by fitting the X-ray spectra to various model at- mospheres, including models attributed to fallback accretion. These models have included: blackbody, hydrogen, helium, iron, silicon-ash atmospheres (Pons et al., 2002), and later extended to two blackbodies, pure silicon, low- iron silicon ash and magnetic hydrogen atmospheres (Walter et al., 2004). An atmosphere model deviates from a blackbody resulting a different deter- minations of the effective temperature of the star and the different atmo- sphere compositions also have different opacities (Zavlin & Pavlov, 2002). The model atmospheres that best fit the spectrum of RX J185635-3754 are those with heavy elements, though the predicted absorption lines associated with the atmosphere composition are not seen (Pons et al., 2002; Walter et al., 2004). After the star has cooled to the point where the nuclear reactions are quenched, the mass fractions in the crust are calculated. These crustal abun- dances have an effect on the atmospheric composition. The nuclear reactions quench when the cooling timescale is shorter than the inverse of the reaction rate. The mass fractions are calculated by using the 489 isotope reaction network torch (Timmes et al., 2000), for each of the cooling cases. The results from the modified Urca case are compared to analytic calculations of the freeze-out of the neutron star crust using the nuclear statistical equi- 19 2.2. Programs Utilized librium program nse (Seitenzahl et al., 2008). The compositions calculated in this chapter are used later in the study of the mechanical properties of the neutron star crust in Chapter 7. 2.2 Programs Utilized In order to calculate the composition of the neutron star crust the programs nse (Seitenzahl et al., 2008) and torch (Timmes et al., 2000) were used. Both of these codes were written by Frank X. Timmes and are available on his website4. The nse code calculates the abundances for a system in nuclear statistical equilibrium. The torch code is a nuclear reaction network which calculates the abundances of isotopes of a system that is not necessarily in nuclear statistical equilibrium. Each of the codes uses the cooling curves and the density of the neutron star crust as input. These two codes are discussed below. 2.2.1 Nuclear Statistical Equilibrium At high temperatures and densities, the nuclei may be in nuclear statistical equilibrium (Shapiro & Teukolsky, 1986). The code nse (Seitenzahl et al., 2008) calculates the abundances of isotopes using the nuclear statistical Saha equations. Nuclear statistical equilibrium relations have the form: (Z,A) + p ⇀↽ (Z + 1, A+ 1) + γ (2.1) (Z,A) + n ⇀↽ (Z,A+ 1) + γ, (2.2) where Z is the atomic charge and A is the atomic mass. These reactions maintain a balance. Here, only reactions of the simplified form: A ⇀↽ B + C, (2.3) 4http://cococubed.asu.edu/code pages/codes.shtml 20 2.2. Programs Utilized will be discussed. For a reaction of the type A → B + C the nuclear Saha equation is nBnC nA = 1 h3 ( 2pi mB mC kT mA )3/2 GBGC GA e−Q/kT , (2.4) where ni is the number density of the atom, mi is the mass, Q is the dif- ference of the binding energies: Q = (mB + mC − mA)c2, and Gi is the statistical weight: Gi = ∑ (2Ir + 1)e −Er/kBT , (2.5) where Ir is the spin of the r th excited state and Er is the energy above the ground state. The nse code contains the reaction rates for 47 isotopes up to nickel. The code takes the temperature, density and electron fraction, Ye, as initial arguments and calculates the abundances of the isotopes in these conditions. 2.2.2 Nuclear Reactions The torch (Timmes et al., 2000) code is a nuclear reaction network which contains the reaction rates for 489 isotopes, up to technetium. The abun- dances of the various isotopes are calculated at a specified density, tempera- ture and burning time. The nuclear reaction network construction is covered in detailed in Timmes et al. (2000) and is summarized below. The reaction network begins by determining the mass fraction, X, of an isotope, i: Xi = ρi/ρ and the corresponding molar abundance of the isotope: Yi = Xi/Ai, where ρi is the mass density of the specific isotope, ρ is the total mass density, and Ai is the mass number of the isotope. From the continuity equation of the isotope, a set of partial differential equations is constructed: dYi dt +∇ · (YiVi) = Ṙi (2.6) where Ṙi is the total reaction rate and Vi is the mass diffusion velocity, which is set to zero. As a result the reaction network is composed of a set 21 2.3. Cooling Curves of ordinary differential equations: dYi dt = Ṙi. (2.7) This system of ordinary differential equations is integrated using a Bader- Deuflhard method coupled with MA28 sparse matrices. Three methods for integrating and eight different matrix packages are compared in Timmes (1999). This code was adapted to read cooling curve files, which contain the temperature and time step, for a specified density. The calculation starts in nuclear statistical equilibrium, but subsequent steps use the abundances calculated in the previous step as initial conditions. 2.3 Cooling Curves Three different cooling curves were used in order to calculated the crustal abundances. These cooling curves included a modified Urca process, cooling of a thick crust, and the cooling curve for a thin crust. Each of the cooling curves start near a temperature of 1010 K and cool until the nuclear reactions are quenched. In the case of the modified Urca process the nuclear reactions are quenched before a year has past. The core temperature and the cooling time are de- termined by the modified Urca equation (Shapiro & Teukolsky, 1986): ∆t = 1yrT−69 (f) { 1− [ T9(f) T9(i) ]6} , (2.8) where T9(f) is the temperature of the outer core (Tc) in units of 10 9 K. For simplicity the temperature of the crust is assumed to be isothermal at densities above 107 g/cm3 (e.g. Heyl & Hernquist, 2001). As the crust of the neutron star cools via the diffusion of heat to the interior in the form of a cooling wave, the surface of the star is expected to remain hotter than the core until the heat reservoir of the crust has been exhausted (Lattimer et al., 1994). The thick and the thin crust cooling models take into account this thermal disconnect between the crust and 22 2.3. Cooling Curves core, whereas the modified Urca model does not. The cooling curve for a thick crust is based on the model in Lattimer et al. (1994). In this model the core is cooling rapidly and the cooling wave takes 15 years to diffuse through the crust. This model is appropriate for a normal neutron star. The nuclear reactions in the crust are quenched before the cooling wave passes through the crust, thus the core cooling does not affect the reactions in the crust. The cooling curve used in this work is taken from Figure 3 in Lattimer et al. (1994), where the hottest parts of each curve set the age and temperature of the neutron star. The thin crust cooling curve is appropriate for a strange star and the expression for the cooling curve is described next. In the thin crust case the core is cooling via the direct Urca process, where the temperature re- lation is t ∼ 20 T−49 s (Lattimer et al., 1991). The crust of the thin crust neutron star is cooled by crustal bremsstrahlung. The cooling of the crust is determined from the crustal luminosity and the total thermal energy of the crust. The luminosity of crustal bremsstrahlung goes as Lbrems ∼ (5× 1039 erg/s)(Mcr/M), where Mcr is the mass of the neutron star crust and the total thermal energy of the crust is Ucr = 3 2kT Mcr Amu (Shapiro & Teukolsky, 1986), where A is the average atomic mass and mu is the mass unit. The crust of a strange star reaches a density of a few times 1010 g/cm3, instead of neutron drip in the case of a normal neutron star (Stejner & Mad- sen, 2005). For a crustal density of 4.8×1010 g/cm3 the equilibrium nucleus is 80Zn (Shapiro & Teukolsky, 1986). Using Figure 1 of Brown et al. (1998), the time scale for the cooling wave to pass through the crust is on the order of a month. As a result, for the thin crust case the crustal cooling is domi- nated by the crustal bremsstrahlung for the first month. After a month has passed the direct Urca process dominates. The temperature of the crust can thusly be described as T = Tbrems(e −(t/month)2) + TDU (1− e−(t/month)2), (2.9) where Tbrems is the temperature of crustal bremsstrahlung and TDU is the temperature at a specific time step using the equation for direct Urca cooling. 23 2.4. Crustal Mass Fraction Calculations The Tbrems starts at a temperature of 10 10 K and decreases by 0.1% each step. The time, or age, of the star is given by the crustal bremsstrahlung equation for a thin crust t = 3 10 kb M Amu 2× 10−31 s K erg T−59 [ 1− T 510 ] (2.10) for a neutron star starting at a temperature of 1010 K, where T10 is the temperature in units of 1010 K. The cooling curves for the three cases are shown in Figure 2.1. For densities above 107 g/cm3 the crust is assumed to be isothermal and Figure 2.1 reflects the cooling curves appropriate for densities at 107 g/cm3 and greater. For densities below 107 g/cm3 the input temperature is interpolated between the core (Tc) and the surface (Ts): log(T ) = [ log(Tc)− log(Ts) 7 ] log(ρ) + log(Ts), (2.11) where the surface temperature is given as: Ts = (10 K 1/2Tc) 2/3 (Shapiro & Teukolsky, 1986). 2.4 Crustal Mass Fraction Calculations For each of the densities investigated, the system within torch is started in nuclear statistical equilibrium. The temperature and the length of time at which the system burns is determined by the cooling curve. After the reactions have quenched the abundances at a specific density are calculated and the isotopes which are sufficiently abundant to compose the atmosphere and envelope are determined. In order for an isotope to form the atmosphere it needs to have a surface density greater than 1 g/cm2 ∼ σT /mp, the ratio of the Thompson cross- section to the mass of the proton. The first step in determining which isotopes rise to the surface is to calculate the total pressure at the density of the calculation. For densities greater than 106 g/cm3 the pressure is given 24 2.4. Crustal Mass Fraction Calculations Cooling Curves Log(Time [sec]) L og (T em p [K ]) 0 1 2 3 4 5 6 7 8 9 7. 5 8. 0 8. 5 9. 0 9. 5 10 .0 Modified Urca Thick Crust Thin Crust Figure 2.1: The three different cooling curves used for the abundance calcu- lations. The methods for cooling include the modified Urca, thick crust and a thin crust. These curves are appropriate for densities at 107 g/cm3 and above. For densities less than 107 g/cm3, the curves would be interpolated between the surface temperature and the temperature at 107 g/cm3. The modified Urca equation is given in Shapiro & Teukolsky (1986). The thick crust is taken from Lattimer et al. (1994) and is appropriate for a normal neutron star. The cooling curve for the thin crust is presented in this work and is appropriate for a strange star. 25 2.4. Crustal Mass Fraction Calculations by: P = 1.2435× 1015 µ 4/3 e ( ρ 1 g/cm3 )4/3 dyne/cm2, (2.12) where it is assumed that relativistic electrons dominate the pressure and that µe = 2. For the case where a density of 10 6 g/cm3 is used in the mass fraction calculation the pressure is given by: P = 1.42180× 1025φ(x) dyne/cm2, (2.13) where φ(x) = 1 8pi2 [x √ 1 + x2(2x2/3− 1) + ln(x+ √ 1 + x2 )] (2.14) and x = pf/mec, the ratio of Fermi momentum (pf ) to the product of the mass of the electron and the speed of light (Shapiro & Teukolsky, 1986). With the pressure calculated the next step is to calculate the column density between the particular density and the neutron star surface: P/gNS, where gNS is the surface gravity given by: gNS = GM R2 ( 1− 2GM c2R )−1/2 . (2.15) For a 10 km and 1.4 M neutron star the surface gravity is 2.43×1014 cm/s2. Finally, the minimum mass fraction required for an isotope to rise to form the atmosphere is given by the ratio of 1 g/cm2 to the calculated column density. For a mass fraction to be sufficient to form the envelope the isotopes which have abundances on the order of parts per million over the entire crust are of interest. On the other hand, in the atmosphere the isotopes need to have abundances on the order of parts per billion or less. Using a typical column density of about 4× 109 g/cm2 and 1 g/cm2 for the envelope and atmosphere, respectively, the minimum mass fractions required for an isotope to float to form the envelope is calculated by scaling according to column densities. 26 2.5. Timescales: Freeze-Out 2.5 Timescales: Freeze-Out To calculate the neutron star freeze-out three different timescales need to be considered: the cooling time (τc, for example Eq. (2.8) for the modified Urca case), the settling time (τs), and the nuclear reaction timescale (τrxn). The composition of the neutron star depends on the timescales; for the case τrxn < τc < τs, the particular species are expected to be in nuclear statistical equilibrium (NSE); for τc < τrxn < τs the particular nuclear reactions are quenched – the star is cooling faster than the reactions can occur; for the case τs < τrxn < τc the isotopes can float up faster than the reactions bring them to NSE. For the case τs < τfreeze the light isotopes can float all the way up to the top before the layer freezes. The time when the particular layer of the neutron star freezes is deter- mined by comparing the potential energy between the ions that compose the crust to their kinetic energy (Shapiro & Teukolsky, 1986), Γ = Potential Energy Kinetic Energy = (Z e)2 ra T kb , (2.16) where e is the electron charge and ra is the ion radius such that the product of (4/3)pir3a and the ion number density is unity. When Γ > 180 the layer freezes out, or crystallizes, and the light species can no longer float upward. 2.5.1 Settling Timescale As well as determining the lightest isotopes which are optically thick, it is important to calculate if the isotopes can even reach the surface before the layer freezes. The gravitational settling could also be so efficient that the light isotopes could float up before the layer cools and the nuclear reactions quench. Thus, it is important to estimate the settling time of the isotopes of interest. 27 2.6. Results Brown et al. (2002) calculated the sedimentation or settling timescale for a neutron star atmosphere, τs ≈ 105s [ Z3.91 Z 0.3 2 ρ 1.3 5 A1.81 g14 2T 0.37 (A2Z1 −A1Z2) ] . (2.17) This is an estimate of the time for a nuclide to settle down over a pres- sure scale height — a negative value means that the nuclide ascends. This timescale is typically around a few years for the envelope and about 106 yr for the outer crust (ρ < 1012 g cm−3) for Silicon-28 in a background of Iron-56. In particular, until the bulk of the Nickel-56 has decayed the settling time for Silicon-28 is much larger because A2Z1 − A1Z2  1; consequently, the nucleons differentiate gravitationally after the nuclear reactions effectively cease, i.e. τrxn < τs. 2.5.2 Reaction Timescale In the modified Urca cooling case the reaction rate timescales were calcu- lated using the subroutines in torch. Each of the reaction rate calculations depends only on the input temperature and the densities of the various species. As there are many different ways to make a specific isotope, e.g. 28Si, the rates which lead to the creation of the isotope are added together to get the timescale of the reaction rate, τrxn. In order to calculate the abundance of the alpha particles and the other species the nuclear statistical Saha equations as implemented in nse are used. A particular species is assumed to freeze out of equilibrium when the reaction timescale exceeds the cooling timescale. The abundance in nuclear statistical equilibrium at the freeze out temperature gives an alternative estimate of the final abundance of the nuclides. 2.6 Results In order to determine the expected composition of the neutron star crust and atmosphere, in the cases of the modified Urca and the thick crust, a 28 2.6. Results density of 107g/cm3 was examined. For the thin crust a density of 107g/cm3 would crystallize before any of the isotopes had time to reach the surface, so a density of 106g/cm3 was examined. The higher density regions of the crust freeze before the lower density, thus the isotopes would have more time to float to the surface from the lower density regions of the crust. The results from the nuclear reaction network are compared with those of a semi-analytic freeze-out calculation in the modified Urca case. Due to gravitational stratification the lightest isotope which can reach the surface would form the atmosphere. Each of the three cases are discussed below, with Table 2.1 summarizing the lightest isotopes which could form the at- mosphere as well as the top three isotopes with the highest mass fractions which thus compose the crust. 2.6.1 Case 1: Modified Urca Mass Fractions At a density of 107 g/cm3 the corresponding pressure is: 1.1×1024 dyne/cm2. At this pressure the column density to the surface is: 4.4 ×109 g/cm2. The resulting minimum mass fraction required for an isotope to be optically thick on the surface is: 2.3×10−10 at this density. Isotopes with a mass fraction greater than 2.3×10−10 will have a surface density of 1 g/cm2. The lightest elements to be optically thick on the surface and have time to reach the surface before crystallization of the layer occurs are shown in Figure 2.2, where the horizontal line indicates the minimum mass fraction required to be optically thick on the surface. The lightest elements to rise to the surface which are optically thick are: 28Si, 32S, 34S, and 36Ar. In particular the abundance of 28Si is about 3 × 10−9 so a layer ∼10 g/cm3 of silicon lies on the surface of the star. Due to the overabundance of silicon and the effects of gravitational stratification other isotopes would not be visible. The top three isotopes of the crust at this layer are 56Fe, 54Fe, and 60Ni. 29 2.6. R esu lts Structure Modified Urca Thick Thin Isotope X Isotope X Isotope X Atmosphere 28Si 2.695×10−9 50Cr 6.108×10−10 40Ca 2.859×10−8 32S 5.956×10−9 53Mn 3.228×10−6 50Cr 5.967×10−5 34S 1.164×10−9 54Fe 7.474×10−6 53Mn 1.337×10−4 36Ar 2.639×10−9 55Fe 1.751×10−5 54Fe 6.477×10−1 — — 57Co 1.492×10−6 55Fe 5.470×10−3 Crust 56Fe 0.559 56Fe 0.9286 54Fe 0.6477 54Fe 0.3187 60Ni 0.03138 58Ni 0.2209 60Ni 0.05175 52Cr 0.02485 56Fe 0.1018 Table 2.1: A summary of the atmosphere and crust compositions along with the mass fraction, X, of each isotope for the three different types of cooling. 30 2.6. Results Freeze-Out Using the steps outlined in Section 2.5 the cooling (τc), settling (τs), and the nuclear reaction (τrxn) timescales for the case of 28Si are calculated. The results of these calculations are displayed in Figure 2.3. These calcu- lations compare the age of the neutron star to the settling, crystallization temperature, creation and destruction timescales of 28Si for two densities: 107 g/cm3 and 1012 g/cm3. The temperatures at which the layers crystal- lize are 4.8 × 107 K and 2.2 × 109 K, for the densities of 107 g/cm3 and 1012 g/cm3, respectively. The creation and destruction rates are given as d lnXi/dt, where Xi is the abundance of the 28Si isotope. These rates are calculated by using the reaction rate routines in torch to determine the energy release per unit mass, these are then multiplied by the abundances calculated from the nse code. The abundances were also calculated using the nuclear Saha equation and are the output of the nse routine. The two different methods for calculating the relative abundances, the output from torch and using nse, are displayed in Figure 2.2. It is clear that the results from torch do not follow nuclear statistical equilibrium precisely. With Figure 2.3 the quenching temperature for the silicon reactions are determined. The reactions become quenched when the reaction rate time becomes longer than the age of the neutron star. From the plot the quench- ing temperature for the reactions which destroy silicon are 2.12× 109 K and 2.15×109 K for the densities of 1012 g/cm3 and 107 g/cm3, respectively. The reactions which create silicon are quenched at 3.05 × 109 K for both of the densities. These approximately give the temperature range over which the abundance of silicon levels out in Figure 2.2. 2.6.2 Case 2: Thick Crust In the case of the thick crust, the density of 107 g/cm3 was examined for the isotopes which can rise to the surface. The pressure, column density to the surface and the minimum mass fraction required to be optically thick are the same for the modified Urca case. For the thick crust the lightest isotopes which are optically thick and can rise to the surface are 50Cr, 53Mn, 54Fe, 31 2.6. Results ρ = 107g/cm3 Temperature[K] L og (R el at iv e A b u n d an ce ) 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0×109 −9 .5 −9 .0 −8 .5 −8 .0 −7 .5 −7 .0 −6 .5 −6 .0 −5 .5 −5 .0 28Si 32S 34S 36Ar Figure 2.2: Lightest isotopes for which the mass fraction abundance would be great enough that the isotope will be optically thick on the neutron star surface for modified Urca cooling. All of these isotopes have time to reach the surface before the layer crystallizes. These are the mass fractions for the density of 107 g/cm3. The corresponding pressure and column density for this neutron star density are 1.1 × 1024 dyne/cm2 and 4.4×109 g/cm2, respectively. The horizontal line indicates the minimum abundance required for an isotope to have a surface density of 1 g/cm2 to be optically thick. The abundances for 28Si, 32S, and 36Ar in nuclear statistical equilibrium using nse are depicted for comparison by the nearly dashed vertical lines, the colours correspond to the isotope type. 32 2.6. Results Freeze-Out log(Temperature [K] lo g( τ [s ]) 8.4 8.6 8.8 9.0 9.2 9.4 9.6 9.8 10.0 0 2 4 6 8 10 12 14 Age Settling Create 28Si Destruction 28Si Figure 2.3: Comparing the cooling, settling, and nuclear reaction timescales of 28Si for the neutron star densities of 107 g/cm3 (solid lines) and 1012 g/cm3 (dotted lines) for modified Urca cooling. The black vertical line denotes the temperature at which the layer at 1012 g/cm3 crystallizes, the temperature for the crystallization of the layer at 107 g/cm3 is cooler than 108.5 K and so lies to the left of the plot. The layer of 1012 g/cm3 crystallizes at 2.2×109 K, and the layer at 107 g/cm3 crystallizes at 4.8× 107 K. Note that the cooling and the creation reaction timescales overlap for the two densities. 33 2.7. Conclusions 55Fe, and 57Co, as shown in Figure 2.4. The three most abundant isotopes of this crust were found to be 56Fe, 60Ni, and 52Cr. 2.6.3 Case 3: Thin Crust For the thin crust the isotopes which had the possibility to be optically thick could not reach the surface before the layer crystallized for densities at 107 g/cm3 and higher; consequently, the layer at a density of 106 g/cm3 was examined. At this density the corresponding pressure is: 2.3×1022 dyne/cm2. The column density to the surface at this pressure is 9.6 × 107 g/cm2, re- quiring a minimum mass fraction of 1.0×10−8 for an isotope to be optically thick on the surface. The isotopes which are not only optically thick, but can also rise to the surface are shown in Figure 2.5. These isotopes include 40Ca, 50Cr, 53Mn, 54Fe, and 55Fe. The isotopes with the top abundances for the crust at this layer were found to be 54Fe, 58Ni, and 56Fe. 2.7 Conclusions In the case of cooling via the modified Urca process, without the thermal influence of a crust, the results (e.g. Fig. 2.3) show that silicon has sufficient time to float to the top from a density of 107 g/cm3 before the layer freezes or the neutron star is observed. However, the settling time from 1012 g/cm3 is found to be too long for light species to float up before the crust freezes; therefore, the atmosphere in this case is likely to be composed of silicon but the envelope is likely to be composed of iron-group elements that have been chemically separated by gravitational settling. Deeper layers are unlikely to be chemically separated, at least by gravity. The torch code has been used in order to calculate the lightest isotopes which would rise to the neutron star surface and these results were compared to semi-analytic calculations. Using the torch code the mass fractions at a density of 107 g/cm3 were calculated. The lightest isotopes found to rise to the surface and be optically thick in the atmosphere are: 28Si, 32S, 34S, and 36Ar. On the other hand, there is not sufficient time for light elements 34 2.7. Conclusions Thick ρ = 107g/cm3 Temperature [K] L og (R el at iv e A b u n d an ce 2.0 2.5 3.0 3.5 4.0 4.5 5.0×109 −1 0 −9 −8 −7 −6 −5 −4 −3 −2 50Cr 53Mn 54Fe 55Fe 57Co Figure 2.4: The lightest isotopes for which the mass fraction abundances are large enough that the isotopes which can reach the surface are optically thick for a thick crust. The mass fractions are calculated for a density of 107 g/cm3. The horizontal line indicates the minimum abundance required for an isotope to have a surface density of 1 g/cm2 and be optically thick. 35 2.7. Conclusions Thin ρ = 106g/cm3 Temperature [K] L og (R el at iv e A b u n d an ce ) 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5×109 −8 −7 −6 −5 −4 −3 −2 −1 0 40Ca 50Cr 53Mn 54Fe 55Fe Figure 2.5: The lightest isotopes which can reach the surface and are also optically thick in the case of a thin crust. These isotopes are calculated for a density of 106 g/cm3. The corresponding pressure and column density at this density are 2.3× 1022 dyne/cm2 and 9.6× 107 g/cm2, respectively. The horizontal line indicates the minimum mass fraction required in order for an isotope to be optically thick. 36 2.7. Conclusions to percolate from high density, such as a density of 1012 g/cm3, to form the envelope. However, calculations of freeze-out at higher densities are required to determine the precise composition of the envelope. Semi-analytic calculations of the freeze-out of 28Si were performed, for the modified Urca case, assuming local thermodynamic equilibrium until the cooling rate exceeds the reaction rate. The rates from torch and the abundances from the code nse were used in order to calculate the rates of nuclear reactions involving 28Si; the reactions become quenched when the reaction time is longer than the age of the neutron star. The creation reac- tions are quenched at a temperature of 3.05× 109K for both the densities of 107 g/cm3 and 1012 g/cm3. The reactions which destroy silicon are quenched at 2.12× 109 K and 2.15× 109 K at a density of 1012 g/cm3 and 107 g/cm3, respectively. The calculated quenching temperatures of the silicon reactions agree with the results calculated using the torch code directly. In the case of the thick crust, the atmosphere could be formed by 50Cr rising to the surface from a density of 107 g/cm3. For a neutron star with a thin crust and direct Urca cooling the layers at 107 g/cm3 and higher crystallize before the isotopes could have time to reach the surface. The atmosphere in the thin crust case could be formed from 40Ca raising to the surface from a layer at a density of 106 g/cm3. The neutron stars with crusts stay hotter longer than a star cooling via modified Urca, thus there is more time with these hotter stars to use up the lighter elements in nuclear reactions. Unless there has been significant accretion from either the supernova debris, the interstellar medium or a companion, neutron star atmospheres are unlikely to be composed of iron, helium or hydrogen. In the case where some accretion occurs, due to the reactive surface on these neutron stars, diffusive nuclear burning could destroy the accreted material (Chang et al., 2010). The type of isotope composing the atmosphere depends on the cool- ing mechanism at work in the star: 28Si for modified Urca, 50Cr for a thick crust, and 40Ca for a thin crust and direct Urca cooling. The formation of these atmospheres can provide additional justification, with fallback ac- cretion, for isolated neutron stars fit with intermediate atmosphere models. 37 2.7. Conclusions Understanding how this novel composition of the atmospheres affects the neutron star emission may provide new insights on the observed spectra of neutron stars. The model thermal spectrum as well as the predicted spec- tral lines for the calculated compositions could be calculated in the future and compared to the emission observed from isolated neutron stars. The composition of the neutron star crust for each of the three cases were also calculated and will be used in subsequent chapters in order to investigate the mechanical properties of each of these crusts. 38 Chapter 3 Molecular Dynamics 3.1 Introduction In the previous Chapter it was described how various crustal compositions are computed. This Chapter explores how mechanical properties are de- termined using the Large Atomic/Molecular Massively Parallel Simulator (LAMMPS). LAMMPS is a software that computes molecular dynamics simula- tions. The output of these simulations can not only be used to determine the breaking strain of the material but also the crystalline structure of the sample material. The Radial Distribution Function (RDF) is used to char- acterize the structure of the simulation as either liquid or solid. In the case where a solid is formed, the type of crystal configuration can also be deter- mined, such as FCC or BCC. Chemical separation in the crust can also be investigated with the calculated RDF of the different types of atoms in the system. Molecular dynamics simulations in general and the LAMMPS software are discussed below. This is followed by a discussion of the calculation of the RDF, which is applied to four example simulations. 3.2 Molecular Dynamics Simulations Two main computational methods which could be used to calculate the mechanical properties are either Monte Carlo or molecular dynamics simu- lations of a system. The Monte Carlo method calculates the probability of a system being in one state versus another state at any specific time. With molecular dynamics simulations the thermodynamic properties are derived by integrating Newton’s equations of motion. Both computational meth- ods calculate the properties of a system of many particles. With molecular 39 3.2. Molecular Dynamics Simulations dynamics a specific force can be applied collectively to the system, instead of atoms being randomly selected, as in Monte Carlo, in order to calculate the effect of the applied force. Molecular dynamics calculations are able to provide a direct determination of the dynamical properties of the system of intrest and will be used in order to understand the neutron star crust. This discussion of molecular dynamics follows Allen (2004). As men- tioned above, molecular dynamics solves Newton’s equations of motion for a system of particles: dri dt = vi dvi dt = − 1 mi N∑ j=1,i 6=j ∂ ∂ri V (|ri − rj |), (3.1) where mi is the mass of atom i, ri and vi are the position and velocity vectors, and V is the potential. The potential with which the particles interact in this work is the screened Coulomb or Yukawa potential: Vij = ZiZje 2 rij e−rijκ, (3.2) where Z is the charge of the particle, rij is the distance between the two interacting particles is rij , and κ is the inverse screening length. With the positions of each of the particles interacting via a defined poten- tial, the equations of motion define the motion of the particles at a moment of time. A typical algorithm used to follow the evolution of particle motion is the velocity-Verlet (Swope et al., 1982), which calculates the positions, velocities and accelerations of all the particles with time. The velocity- Verlet algorithm is similar to the Leapfrog algorithm (Hockney & Eastwood, 1988), both being variations of Verlet integration (Verlet, 1967), a method to integrate Newton’s equations of motion. With Leapfrog integration the position and velocity calculation are offset by half a time-step, whereas for velocity-Verlet integration the position and velocity are calculated at the same time-step. The velocity-Verlet algorithm results from an expansion of 40 3.2. Molecular Dynamics Simulations the particle position with time, r(t+ ∆t), r(t+ ∆t) = r(t) + v(t)∆t+ a(t) 2 ∆t2 + ... (3.3) and is time reversible. As ṙi = ṗi/mi and ṗi = fi the velocity-Verlet algo- rithm can be expressed as: pi(t+ 1 2 ∆t) = p + 1 2 ∆t fi(t) ri(t+ ∆t) = ri(t) + ∆t pi(t+ 1 2 ∆t)/mi pi(t+ ∆t) = pi(t+ 1 2 ∆t) + 1 2 ∆t fi(t+ ∆t). (3.4) The velocity-Verlet algorithm works by evaluating the above equations, as well as calculating the force for each of the interactions of the particles in the system. The thermodynamic properties of the system are determined by averag- ing the properties of all the particles in the system. For example, consider the calculation of the kinetic energy. The kinetic energy for each particle is given by: 1/2mv2. An average of the kinetic energy of all the particles would result in the measured kinetic energy thermodynamic property of the system. In order to decrease the number of calculations performed and thus short- ening the time required for the simulations to complete, neighbour lists can be constructed. With neighbour lists the particles which are within a cut- off radius from a particular particle are listed and it is only these particles which have their interactions calculated. The neighbour list construction is schematically displayed in Figure 3.1. The determination of the cut-off radius for a simulation is discussed in Chapter 4. An extra ‘buffer’ region is also defined in order to keep track of particles which could cross over into the cut-off radius boundary during each time step. The frequency at which the list is constructed is defined for each of the simulations. Molecular dynamics simulations typically utilize periodic boundary con- ditions. The use of periodic boundary conditions decreases the surface effects 41 3.3. LAMMPS Figure 3.1: An illustration, based on Figure 6 from Allen (2004), of the neighbour list construction. The potential cut-off radius is indicated by the solid circle, and the buffer region is indicated by the dashed circle. The white particles would have the potentials calculated, the gray particles are kept track of in the neighbour list, but do not have the potential calculated. The black particles are not included in the neighbour list at all. which result from small simulation sizes. With periodic boundary conditions particle images are created such that a particular atom interacts with the neighbouring particles, as well as the particle images. In the situation where an atom leaves the simulation box, the potential between the particle of in- terest and the incoming particle are considered in the calculation. Periodic boundary conditions are used in all the simulations discussed within this thesis. 3.3 LAMMPS The molecular dynamics simulations presented in this thesis were computed with the software LAMMPS (Large Atomic/Molecular Massively Parallel Sim- ulator). LAMMPS is an open-source software made available by Sandia Lab- oratories5 and is fully documented on the LAMMPS website6. LAMMPS is con- 5http://lammps.sandia.gov 6http://lammps.sandia.gov/doc/Manual.html 42 3.3. LAMMPS structed to run on and take advantage of parallel computing systems (Plimp- ton, 1995), which can significantly decrease the time it takes for a simulation to complete. Neighbour lists are also utilized within LAMMPS in order to in- crease computational efficiency. All the simulations presented in this thesis used the velocity-Verlet algorithm to integrate the equations of motion. All the simulations use reduced units. The reduced units are based on the Yukawa potential, Vij = ZiZke 2 r e−κr, (3.5) where r is the distance between particles, Z is the charge and κ is the inverse screening length. For the reduced units, characteristic quantities are now introduced: a characteristic length scale, energy, mass and charge. The characteristic length scale, a, is parameterized by a = n−1/3, where n is the number density. The characteristic energy scale, Uo, is given by Uo = (Z e)2 a , (3.6) which is expressed with respect to a charge of Z = 26. The characteristic mass, ma, is expressed in units of the mass of iron, ma = 56/A, thus ma = 1 for iron. The characteristic charge is expressed in units of the iron charge, Z = 26, thus Za = Z/26. The reduced units are related to the physical quantities through these characteristic quantities; these relationships are listed in Table 3.1. For example, the characteristic quantities for 56Fe at a density of ρ = 1011 g/cm3 are a characteristic length of a = 10−11 cm and a characteristic energy of Uo = 10 −5 erg. Using the reduced units the Coulomb parameter can be expressed as Γ = (Ze)2 rkbT = (〈Z〉 Z26 )2 1 T ∗ a r , (3.7) where 〈Z〉 is the average charge of the system and r is the ion spacing. In the simulations the number density is set to unity, thus both a and Uo are unity for the simulations. When Γ = 180 a transition between liquid and solid occurs. 43 3.4. Characterization of the Output Unit Reduced Unit Temperature T* = T kbUo Energy E* = EUo Distance r*= ra Pressure P* = P a 3 Uo Time t* = t/(ma 2 2Ua )1/2 Table 3.1: The reduced unit system used in the LAMMPS simulations. Pa- rameter values quoted in this thesis from the simulations are in reduced units. The physical quantities are determined by using appropriate values the characteristic distance, a, and energy, Uo, parameters for the material. Different versions of LAMMPS were used for the work in this thesis. The bulk of the calculations of the mechanical properties the LAMMPS version 21 Jul 2010 was used. When other versions are were used in calculations, the version of LAMMPS is stated. All simulations were performed with pe- riodic boundary conditions. As the number density is set to unity for all simulations, the number of unit cells dictates the box size of each of the simulations. The number of particles and the box lengths are listed in Table 3.2 for the five sizes of simulation box sizes used in this thesis for FCC and BCC crystals. The time for a calculation to be performed increases with the number of particles in the simulation box. In Section 6.3.1 the simple shear results from the different box sizes are examined, though a thorough test of the dependence on particle number and the calculated properties would require a larger range of simulation box sizes. There are always twice as many particles in the FCC boxes as the BCC boxes. 3.4 Characterization of the Output The output of the molecular dynamics simulations are then characterized by calculating the radial distribution function (RDF) of the sample. With the RDF it is possible to examine correlations between all the particles of the sample, particles of the same type and cross correlations between particles of different types. The state of the system – liquid or crystal – can also be 44 3.4. Characterization of the Output Unit Cell BCC FCC Volume Particles Box Length Particles Box Length 10× 10× 10 2000 12.599 4000 15.874 12× 12× 12 3456 15.119 6912 19.049 20× 20× 20 16000 25.1984 32000 31.748 25× 25× 25 31250 31.498 62500 39.685 35× 35× 35 85750 44.0972 171500 55. 559 Table 3.2: The number of particles and box length in reduced units for BCC and FCC crystals for the different box sizes used in this thesis. These values are for a number density set to unity. examined. In this section the method for calculating the RDF is discussed. The procedure for the RDF is applied to example simulations: FCC and BCC crystals, of one atom type and two, which are subsequently melted. The RDF is calculated both in the initial solid configuration and in the final, liquid, configuration. 3.4.1 Radial Distribution Function: RDF The RDF is the same as the pair correlation function. The calculation indicates at which distances scales there is structure, or a lack thereof. Peaks in the RDF are indicative of structure occurring on those length scales, while valleys indicate voids of structure. The distance scales considered for the RDF calculations are from r = 0 to half the box length size, L/2, where L is the total box length. The parameters required for the RDF calculation are a list of all the particle positions and the box length, L. The box is divided into 1000 shells of thickness dr, where dr = (L/2)/1000. The RDF is calculated by essen- tially creating a histogram of the number of pairs, g[k], within a distance interval or shell r + dr. In this histogram, k = r/dr, which determines the bin number, ensuring the g[k] array begins at g[0]. The distance between each pair of particles is calculated as: distance = √ (∆x)2 + (∆y)2 + (∆z)2. (3.8) 45 3.4. Characterization of the Output Due to the periodic boundary conditions there are two possible distances between particle pairs. In order to account for the pair distance degeneracy, the smallest distance between particles is used in the RDF calculation. For example, if ∆x is greater than L/2, the distance ∆x is subtracted from the box length L: ∆x2 = L − ∆x1, and it is this ∆x2 which is used in the calculation of the pair distance. As long as this pair distance is within L/2 the pair is used in the calculation g[k]. The actual RDF, g(r), is then calculated by normalizing, with respect to density, the histogram g[k]: g(r) = 2 g[k] (4pi/3) n N (r3upper − r3lower) , (3.9) where n is the number density and N is the total number of atoms in the sample investigated, rlower and rupper are the bounds of the shell, with rlower = k dr and rupper = dr+rlower, and r is the midpoint of the shell. This normalization ensures that g(r) = 1 for an isotropic distribution, though the cross-correlation RDF calculations do not normalize to unity. The RDF was calculated for four example simulations, in crystal and liquid state. The ex- amples include simulations with two atom types in which cross correlations are also investigated. 3.4.2 Parameters of Example Simulations In order to demonstrate the RDF calculation, four example systems were examined. These example simulations consist of FCC and BCC crystals with one atom type, and FCC and BCC crystals with two different types of atoms. These systems began in a crystalline state and then were heated to form a liquid. There were 62500 atoms used in the FCC crystal case and 31250 atoms in the BCC crystal case. The mass of each atom is the same value of 1.0, corresponding to an atomic mass of A = 56, and the initial temperature was 0.001. Each simulation had a time step of 0.09 and ran for 1 million steps. The atoms interacted via the Yukawa potential. The temperature increase was controlled by a Langevin thermostat. With a Langevin thermostat the system of particles in the simulation are within a 46 3.4. Characterization of the Output heat bath. From collisions of the simulation particles with the heat bath a frictional force as well as a random force, which imparts a random velocity kick, controls the temperature of the system (Schneider & Stoll, 1978). The Langevin thermostat just controls the simulation temperature, the time in- tegration was performed with the nve fix: constant number, volume, and energy, which accounts for the modified forces due to the Langevin thermo- stat. The 7 July 2009 version of LAMMPS was used in these simulations. The inverse screening length values appropriate for FCC and BCC crys- tals were taken from Robbins et al. (1988). A discussion of the cut-off radius determination is covered in Chapter 4. For the FCC crystal simulations the atoms interacted via the Yukawa potential with an inverse screening length of κ = 3.74 and a cut-off radius of interaction of rc = 3.07 reduced units. The temperature was increased from 0.001 to 0.004. In the FCC crystal case of one atom type the pair coefficient, Z2e2, was set to 1.0 for all atoms. For the two atomic types simulations, there were 15625 atoms of type 1 and 46875 atoms of type 2, 1/4 and 3/4 of the total number, respectively. The type 1 atoms interacted with a pair coefficient of (Ze)2 = 1.0, the type 2 atoms interacted with a pair coefficient of (Ze)2 = 1.1598 and the type 1 and type 2 particles interacted with each other with a pair coefficient of (Ze)2 = 1.0769. The two atom types correspond to 56Fe and 56Ni for type 1 and 2, respectively. These pair coefficients for the two types of atoms simulations are listed in Table 3.3. The mass of all atoms, even of different types, was the same value at 1.0 mass units. For the BCC crystal simulations the atoms interacted via the Yukawa potential with an inverse screening length of κ = 0.88 and a cut-off radius of interaction of rc = 9.0 reduced units. In order to melt the crystal the tem- perature was increased from 0.001 to 0.04. For the BCC crystal simulation of one atomic type, all the particles interacted with the same pair coefficient of (Ze)2 = 1.0. In the case of the BCC crystal simulation with two types of atoms, there were 15625 atoms of type 1 and 15625 atoms of type 2, half of the total number each. The pair coefficients used and the mass of the particles in the simulations are the same in the BCC case as the FCC case. 47 3.4. Characterization of the Output Type 1 2 1 1.0 1.0769 2 1.0769 1.1598 Table 3.3: Pair coefficients, (Ze)2, used in the simulations containing two atom types, for both the FCC and BCC lattice example simulations. The type 1 atoms correspond to 56Fe and the type 2 atoms correspond to 56Ni. 3.4.3 RDF: One Atomic Type The RDF was calculated of the initial structure of the simulation and the structure after melting the FCC and BCC crystals. In this section only the results of the FCC and BCC crystals containing one atomic type, where the pair coefficient is set to 1.0 for all particles, are investigated. The distance scale in the RDF plots are to a distance of r = 6 in order to display the short- length scale structure. The initial crystalline structures of the two different crystals are compared in Figure 3.2. The FCC and BCC crystals present different spacings in the peaks of the RDF. Both crystals were heated until a liquid formed and the RDFs of the liquid FCC and BCC structures are compared in Figure 3.3. Both the FCC and BCC curves display structure on the small distances, but at large distances both FCC and BCC curves go to unity, indicating a lack of large distance scale structure in the liquid form. 48 3.4. Characterization of the Output g( r) FCC and BCC, Crystal 0 10 20 30 radius 0 1 2 3 4 5 6 10 20 30 FCC BCC Figure 3.2: The RDFs of the initial FCC and BCC crystal structures. The atoms in both crystal types interacted via a Yukawa potential. The FCC atoms interacted with an inverse screening length of κ = 3.74 and a cut-off radius of rc = 3.07. In the case of the BCC crystal, an inverse screening length of κ = 0.88 was used with a cut-off radius of rc = 9.0. The two different crystals have the same initial peak, but different peak spacings after the initial peak. 49 3.4. Characterization of the Output g( r) FCC and BCC, liquid 0. 0 0. 5 1. 0 1. 5 2. 0 2. 5 radius 0 1 2 3 4 5 6 0. 5 1. 0 1. 5 2. 0 2. 5 FCC BCC Figure 3.3: The RDFs of the melted FCC and BCC crystals. The atoms interacted via a Yukawa potential with an inverse screening length of κ = 3.74 and a cut-off radius of rc = 3.07 for the FCC crystal and κ = 0.88 with rc = 9.0 for the BCC crystal. The temperature was increased from a value of T ∗ = 0.001 to T ∗ = 0.004 in the FCC case and T ∗ = 0.001 to T ∗ = 0.04 in the BCC case. Both curves demonstrate that there is structure on small distance scales, but not at the larger distances. The FCC simulation contained twice the number of atoms as the BCC simulation. 50 3.4. Characterization of the Output 3.4.4 RDF: Two Atomic Types Simulations of melting FCC and BCC crystals composed of two types of atoms were also performed in order to investigate the structure and cross- correlations in these simulations. The type 1 atoms correspond to 56Fe and the type 2 atoms correspond to 56Ni, with the pair coefficients listed in Table 3.3. The two atom type FCC crystal consisted of 15625 type 1 atoms, with a number density of 0.25, and 46875 type 2 atoms, with a number density of 0.75. The FCC unit cell has four atoms, the type 1 atoms are placed at one point of the unit cell and the type 2 atoms at the other three. In the BCC crystal case, of a total 31250 atoms in the simulation, half were type 1 and the other half of the atoms type 2, for a number density of 0.5 for each type. The BCC unit cell has 2 atoms, as a result each type is at one of the unit cell points. The RDF was calculated for the type 1 particles, the type 2 particles, and for the cross-correlations of the two types of particles in both the FCC and BCC crystal cases. These different RDFs are compared in Figures 3.4 and 3.5 for the initial FCC and BCC crystalline structure, respectively. For the FCC crystal, Figure 3.4, the two atom types have different peak spacings, due to the location of each atom type in the unit cell. For the BCC crystal, Figure 3.5, the number densities of each atom type are the same and the RDFs of the two atom types have the same spacings. As in the one atomic type FCC and BCC crystals, the two atom type crystals were heated to melting and RDFs were taken of the liquid state of the two atom types as well as for the cross-correlations between the two types of atoms. The liquid state RDFs of the FCC crystal are shown in Figure 3.6. The three curves – type 1, type 2, and cross-correlation – display similar shapes and at large distances a lack of structure. The pair correlations of just the type 1 and just the type 2 atoms go to unity at large distances. The cross-correlations between type 1 and type 2 atoms do not normalize to unity at large distances, but the flat curve does indicate a lack of structure. The melted BCC crystal curves are compared in Figure 3.7, in the liquid state all – type 1, type 2, and cross-correlation – indicate that at large distances 51 3.5. Conclusions there is no structure evident. These different RDFs are not typically used in astronomy studies and these examples demonstrate the different types of situations which will be examined with regards to the neutron star crust. The possibility of chemical phase separation can be investigated with the use of these example RDFs. A signature of chemical phase separation in the RDFs would be, for example, one type of atom to be in the liquid state and another atom type in the solid state. Chemical phase separation in the neutron star crust is investigated in Section 7.3. 3.5 Conclusions Molecular dynamics is a useful computational tool in order to determine not only the mechanical properties of a material, but also the structure of the sample. Molecular dynamics simulations solve Newton’s equations of motion for a system of particles interacting with a defined potential. The molecular dynamics simulations in this work are performed using the software LAMMPS. The structure of the results to these calculation can be characterized using the RDF. FCC and BCC example simulations have been calculated and the output characterized with the RDF. These example simulations are useful for comparison to simulations in subsequent chapters. The RDFs of the FCC and BCC crystals have been shown to have different spacings of the peaks. Calculating RDFs of a solid is a method which can be used to discriminate between a FCC or BCC global structure. By computing RDFs it is also possible to investigate the possibility of chemical phase separation in a situation where there are different types of particles. With chemical phase separation certain atom types may be found in the liquid state or in the solid state, which would be indicated in the RDFs of individual atom types. The radial distribution function calculation will be used in following chapters in order to characterize the structure of the neutron star crust. 52 3.5. Conclusions g( r) FCC: two types, Crystal 0 10 20 30 0 10 20 30 radius 0 1 2 3 4 5 6 5 10 15 type 1 type 2 cross Figure 3.4: The RDFs of the initial two atomic type FCC crystal. The simulation contains 15625 type 1 atoms and 46875 type 2 atoms. The atoms interacted with a screening length of κ = 3.74 with a cut-off radius of rc = 3.07. The pair coefficients between the two atom types are listed in Table 3.3. The top figure is the RDF of just type 1 atoms, the middle figure is the RDF of type 2 atoms, and the bottom figure is the cross-correlations of the two atoms types. Each of the RDFs show different spaces of the peaks. 53 3.5. Conclusions g( r) BCC: two types, Crystal 0 10 20 30 0 10 20 30 radius 0 1 2 3 4 5 6 10 20 30 type 1 type 2 cross Figure 3.5: The RDFs of the initial two atomic type BCC crystal. The simulation contains 31250 type 1 atoms and 31250 type 2 atoms. The atoms interacted with a screening length of κ = 9.0 with a cut-off radius of rc = 0.88. The pair coefficients between the two atom types are listed in Table 3.3. The top figure is the RDF of just type 1 atoms, the middle figure is the RDF of type 2 atoms, and the bottom figure is the RDF of the cross- correlations of the two atom types. The peak spacings of the type 1 and type 2 RDFs are the same in this case, unlike the FCC crystal. This is to be expected as the BCC unit cell only has two atoms. 54 3.5. Conclusions g( r) FCC: two types, Liquid 0 1 2 0 1 2 radius 0 1 2 3 4 5 6 0. 2 0. 4 type 1 type 2 cross Figure 3.6: The RDFs of the melted two atomic types FCC crystal. The temperature was increased from T ∗ = 0.001 to T ∗ = 0.004. The simulation parameters are the same as for Figure 3.4. The top figure is the RDF of just type 1 atoms, the middle figure is the RDF of type 2 atoms, and the bottom figure is the RDF of the cross-correlations of the two atom types. The curves, which behave in similar fashions, have a lack of structure for large distances indicated by the flat curve. 55 3.5. Conclusions g( r) BCC: two types, Liquid 0. 0 0. 5 1. 0 1. 5 0. 0 0. 5 1. 0 1. 5 radius 0 1 2 3 4 5 6 0. 2 0. 4 type 1 type 2 cross Figure 3.7: The RDFs of the melted two atomic types BCC crystal. In order to melt the crystal the temperatur was increased from T ∗ = 0.001 to T ∗ = 0.04. The parameters for the simulation are the same as for Figure 3.5. The top figure is the RDF of just type 1 atoms, the middle figure is the RDF of type 2 atoms, and the bottom figure is the RDF of the cross-correlations of the two atom types. The curves are all very similar with g(r) becoming flat at large distances, indicating a lack of structure. 56 Chapter 4 Cut-Off Radius and Total Energy 4.1 Introduction In molecular dynamics simulations the interactions between all the vari- ous particle pairs are calculated. However, since the potential decreases with increasing separation, the impact of particle-particle interactions be- comes unimportant for large relative distances. This property can be used to speed up the time it takes for each calculation, with a cut-off radius uti- lized in conjunction with a neighbour list. By employing a cut-off radius only the interactions between particle pairs within the specified radius are calculated, those particles outside the cut-off radius are not included in any pair-interaction calculations. In order to determine the optimal cut-off ra- dius to use, the total energy at zero temperature is calculated for varying cut-off radii. This calculation of the total energy at zero temperature results in just the energy due to pair interactions via the Yukawa potential. The total energy at zero temperature as a function of cut-off radius is an indi- cation of how much energy due to Yukawa interactions is being neglected and optimal cut-off radius which balances speed and energy. An arbitrarily selected difference of 1% between the expected total energy and the total energy is used as an indication of the optimal cut-off radius to be used in the molecular dynamics simulations. This choice for the optimal cut-off ra- dius was calculated for both FCC and BCC crystal configurations and the results of these calculations were compared to the cut-off radius employed in Robbins et al. (1988). 57 4.2. Simulation Set-Up 4.2 Simulation Set-Up The calculations of the total energy of the FCC and BCC crystals were performed on an Opteron-Beowulf cluster using the July 21, 2010 version of LAMMPS. For the FCC energy calculations 62500 atoms were used. The BCC energy calculations contained 31250 atoms. See Table 3.2 for the relationship to box size. Only one calculation was required for each of these simulations: the total energy at the start of the run where the temperature is zero. At zero temperature the non-thermal energy component is determined. For each of these runs the cut-off radius was varied in order to determine the effect of the cut-off radius on the total energy at zero temperature. For all the simulations periodic boundary conditions were used. 4.3 Results: FCC Crystal For the FCC crystal simulations an inverse screening length of κ = 3.74 was used. The cut-off radius for calculating particle interactions ranged from 1.0 to 10.0 in reduced distance units. The results of the total energy calculations at zero temperature as a function of cut-off radius are listed in Table 4.1. At a cut-off radius of 7.0 the total energy converges to the 10−9 level, as displayed in Figure 4.1. The point where the total energy is within 1% of the maximal total energy occurs between a cut-off radius of 2.0 and 3.0; the horizontal line in Figure 4.1 indicates this radius value. A cut-off radius of 3.0 is found to be within 0.04% of the maximal total energy. 4.4 Results: BCC Crystal For the BCC crystal calculations an inverse screening length of κ = 0.8835, appropriate for a charge of Z = 26, where κ = (24α3/2 pi1/2 nZ)1/3 with a number density of 1.0, was used. The cut-off radius was varied from 1.0 to 15.0 for these calculations. The results of the total energy calculations for the 31250 particle BCC crystal are listed in Table 4.2, and graphically in Figure 4.2. In the FCC case the total energy was found to converge to 58 4.4. Results: BCC Crystal rc Energy % Discrepancy 1.0 0 100 1.5 0.080317068 11.6 2.0 0.089598473 1.34 3.0 0.090773608 0.042 4.0 0.090808517 0.0036 5.0 0.090811738 7.2×10−5 6.0 0.090811801 2.2×10−6 7.0 0.090811803 0 8.0 0.090811803 0 9.0 0.090811803 0 10.0 0.090811803 0 Table 4.1: The calculated total energy of a κ = 3.74, 62500 particle, FCC crystal at zero temperature for various cut-off radii. The percent discrepancy listed in the table is the discrepancy of the simulation energy from the converged energy value of 0.090811803 energy units. the 0.01% level, whereas the energy calculations for the BCC crystal do not converge to the same degree for the cut-off radii investigated. Another method is employed in order to determine the total energy expected, and also the radius at which there is a 1% difference between the total and expected energy. In order to determine the total energy expected, the number of pairs, NP (r), within a shell of ∆r thickness, and the energy at the shell midpoint r is calculated as: E(r) = NP (r) NT e−κr r , (4.1) where NT is the total number of particles for the crystal. The energy at the specific radius is a summation of all energies calculated at smaller radii. This E(r) is plotted with the total energy from the simulations in Figure 4.2. The two energies appear to be in agreement with each other. The 1% difference from the maximum pair energy is indicated by the horizontal line. At a cut-off radius of rc = 8 the calculated energy is within the 1% level. 59 4.5. Discussion and Conclusions rc Energy % Discrepancy 1.0 0 0 1.5 2.1803033 68.8 2.0 2.877923 58.9 3.0 4.9195019 29.7 4.0 5.9864294 14.4 5.0 6.4834227 7.3 6.0 6.7404808 3.6 7.0 6.8817533 1.6 8.0 6.9400936 0.79 9.0 6.9730391 0.32 10.0 6.9862956 0.13 11.0 6.9920647 0.051 12.0 6.9948678 0.011 13.0 6.9961824 0.0076 14.0 6.9967485 0.016 15.0 6.9969933 0.019 Table 4.2: The total energy calculations for a κ = 0.8835, 31250 particle, BCC crystal at zero temperature for various cut-off radii. The percent discrepancy listed in the table is the discrepancy of the energy calculated from the simulations to the maximum energy calculated using the RDF of the particle positions. The maximum energy from the RDF calculations is 6.995648, note that the simulations at larger cut-off radii have a larger energy than this value. 4.5 Discussion and Conclusions In order to determine the ideal cut-off radius for calculating the pair in- teraction, simulations at zero temperature for a varying cut-off radius were performed. Both FCC and BCC crystals were investigated. For the FCC simulation, using a cut-off radius of rc = 3.0, or κrc = 11.48, the total en- ergy is within 1% of the maximum total energy. For the BCC crystal the total energy is within 1% of the calculated maximum total energy with a cut-off radius of rc = 8.0, or κrc = 7.068. In Robbins et al. (1988) for a λ = κa ≥ 3, where a = n−1/3, a cut-off ra- dius of rc = 3.07a was used, or that κrc ≥ 8 was required to achieve accurate 60 4.5. Discussion and Conclusions molecular dynamics simulations. At zero temperature the crystal structure is BCC below λ = 1.72 and FCC above this value (Robbins et al., 1988). This κrc ≥ 8 limit from Robbins et al. (1988) is based on molecular dynam- ics simulations which demonstrate the force decreasing at large distances, as well as the pair-correlation function flattening for distances greater than rc = 8/κ. The results in these current simulations are in agreement with the cut-off radius as defined in Robbins et al. (1988). As a result, for λ ≥ 3 a cut-off radius of 3.07a is used in future simulations, and a cut-off radius of rc = 8/κ is used for smaller λ values in the simulations. 61 4.5. Discussion and Conclusions l l l l l l l l l l l FCC Crystal Cut off Radius T ot al E n er gy 1 2 3 4 5 6 7 8 9 10 0. 01 0. 02 0. 03 0. 04 0. 05 0. 06 0. 07 0. 08 0. 09 0. 10 FCC-1% Figure 4.1: The total energy at zero temperature for a 62500 particle FCC crystal as a function of cut-off radius. At a cut-off radius of 7 the total energy calculated converges, see Table 4.1 for energy values. The horizontal line labelled ‘FCC-1%’ indicates the total energy which is within 1% of the converged total energy. 62 4.5. Discussion and Conclusions l l l l l l l l l l l l l l l l l BCC Crystal Comparison Cut off Radius T ot al E n er gy 0 2 4 6 8 10 12 14 1 2 3 4 5 6 7 8 BCC Data Pair Energy Pair Energy-1% Figure 4.2: Comparing the total energy from the simulations, labelled ‘BCC Data’, to the total energy calculated from the pair distribution, labelled ‘Pair Energy’. The BCC simulations did not converge to the 0.01% level as the FCC simulations did. The horizontal line ‘Pair Energy-1%’ indicates where the simulation energy agrees to the maximal pair energy to the 1% level. 63 Chapter 5 The Effect of Electron Screening on Melting Temperature 5.1 Introduction Two different methods were used in order to determine the appropriate cut- off radius to use in the molecular dynamics simulations. In Chapter 4 the cut-off radius was determined using the total energy at zero temperature. In this chapter another method is used to determine appropriate simulation parameters – calculating the melting temperature of the system of particles and then comparing the results to those in Robbins et al. (1988). In Rob- bins et al. (1988) the phase diagrams and dynamical properties of a Yukawa system of particles were investigated. This studied resulted in the melting transition line of as a function of screening length, as well as the transition line from FCC to BCC crystals. Comparing the calculated melting temper- ature at various cut-off radii to the phase diagram in Robbins et al. (1988) provides not only a test of the cut-off radii as determined in Chapter 4 but also sets a benchmark for the molecular dynamics simulations as performed with LAMMPS. In the limit where the inverse screening length, κ, goes to zero the Coulomb limit of the Yukawa system is reached. This Coulomb limit is also known as a one component plasma (Robbins et al., 1988). In the Coulomb limit the melting temperature can be characterized by the Coulomb coupling parameter Γ, which is the ratio of the potential and thermal energy (see Eq. 64 5.1. Introduction 2.16). When Γ becomes large enough, crystallization occurs (Shapiro & Teukolsky, 1986). This transition, between liquid and solid, has been in- vestigated via Monte Carlo simulations; in Stringfellow et al. (1990), for example, the phase transition was found to depend on the structure of the material, with ΓBCC = 178 and ΓFCC = 192. The properties of dense ion- ized matter have also been covered in a series of papers. In particular, in Pollock & Hansen (1973) Monte Carlo simulations were performed in which a melting parameter of Γm = 155±10 was determined. Potekhin & Chabrier (2000) found the phase transition to be Γm = 175.0 ± 0.4 in the Coulomb limit when an analytic equation of state was employed. The melting temperature in a Yukawa system with a non-zero inverse screening length has been investigated many times, for example in Meijer & Frenkel (1991), Robbins et al. (1988), and Vaulina et al. (2002). In such systems the melting temperature can be characterized by the Lindemann criterion, where an ion lattice melts when the thermally-induced root-mean squared fluctuations of the ion position compared to the original ion position ( √〈δr〉2/r), reach a high enough value (Shapiro & Teukolsky, 1986). Melting was found to occur at values of 0.19 and 0.16 for BCC and FCC structures, respectively (Meijer & Frenkel, 1991). In Robbins et al. (1988) the melting of Yukawa systems with varying screening lengths was found to occur near 0.19 for both BCC and FCC structures. The inclusion of electron screening effects increases the temperature at which a crystal melts (Robbins et al., 1988). In Vaulina et al. (2002) a single parameter was found to determine the melting line of a multicomponent plasma within the Yukawa regime. This parameter was a modified Coulomb coupling parameter. In this chapter LAMMPS is used in order to determine how the melting temperature depends on the value of the inverse screening length, κ. These melting temperature results are then compared to the phase transition re- sults of Robbins et al. (1988). The goal of this comparison is the first check the determination of the cut-off radius to be used in future simulations and second, using the Robbins et al. (1988) as a benchmark as indication if there are problems with the LAMMPS simulations. The phase diagram of a Yukawa system has been extensively studied, any major deviations found from this 65 5.2. Electron Screening Length would be an indication of possible errors within the simulation. The version of LAMMPS used for these simulations is 7 Jul 2009. Initial structure of the crystal in the simulation was determined by the BCC to FCC transition as discussed in Robbins et al. (1988). 5.2 Electron Screening Length In the case of degenerate, relativistic electrons, the electron screening length can be determined by the Thomas-Fermi approximation. In this approxi- mation the inverse screening length is given as κ = 2α1/2kF /pi 1/2, where kF is the electron Fermi momentum. The Fermi momentum is given by: kF = ( 3pi2 ρ mu Z A )1/3 , (5.1) where ρ is the mass density, mu is the atomic mass unit, A is the mass number, and Z is the charge number. Substituting the Fermi momentum into the expression for inverse screening length results in an expression for κ, κ = 2α1/2 pi1/2 ( 3pi2 ρ mu Z A )1/3 , (5.2) which depends on mass density and the ratio of Z/A. Following Robbins et al. (1988), the unitless parameter λ is now intro- duced, where λ = κa. As in Robbins et al. (1988), a is the characteristic inter-particle distance, a = n−1/3, where n is the number density. Us- ing the λ parameter a relationship between the inverse screening and the charge Z can be determined. The Fermi momentum can be expressed as kF = (3pi 2ne) 1/3, where ne is the electron number density. The electron num- ber density is related to the number density as ne = 〈Z〉n. Thus, the Fermi momentum can be expressed as: kF = (3pi 2〈Z〉n)1/3. As a result the λ pa- rameter is related to the Fermi momentum with λ = κa = 2α1/2kF n −1/3/pi1/2 or λ = κa = (24α3/2pi1/2〈Z〉)1/3. (5.3) 66 5.3. Expected Melting: Robbins et al. (1988) With this equation it is thus possible to relate the inverse screening length to the charge of the system of particles, and this equation is used to specify the inverse screening length for the pure and impure neutron star crusts,which is presented in Chapters 6 and 7, respectively. 5.3 Expected Melting: Robbins et al. (1988) Simulations are performed in which a crystal is heated past its melting point. These melting temperatures calculated from the LAMMPS simulations are used in comparison to the results of Robbins et al. (1988). In Robbins et al. (1988) the melting temperature as a function of λ was found to have the following relationship: T̃ = 0.00246 + 0.000274λ, (5.4) where T̃ is defined to be kbT/(ma 2ω2E) and ωE is the Einstein frequency and ω2E = 〈ω〉2, where 〈ω〉2 is the mean squared phonon frequency. In order to compare the melting temperature (equation 5.4) to the results of the simulations, both T̃ and the simulation results are expressed in terms of the unitless energy ratio kBT/Ua, where Ua = Uo e −λ and the energy Uo = (Ze) 2/a. The parameter T̃ is related to the energy ratio using the values of mω2Ea 2/Uo as a function of λ found in Table 1 of Robbins et al. (1988). Setting the value of mω2Ea 2/Uo found in the table as a parameter β(λ), T̃ can be rewritten as: T̃ = kbT ma2ω2E = kbT βλUo = kbT β(λ)Ua eλ , (5.5) with Uo = Uae λ. As a result, the expression for kbT/U in terms of the T̃ parameter is: kBT Ua = T̃ βeλ. (5.6) In this manner the expression for the melting temperature can be compared to the melting temperature from the LAMMPS simulations. The β(λ) value 67 5.4. Simulation Parameters and Calculations used to convert T̃ to the energy ratio is dependent on the structure, either FCC or BCC. For values of λ = 3 and smaller the structure is BCC and FCC for the larger λ values. 5.4 Simulation Parameters and Calculations The melting temperature simulations were run with the 7 July 2009 version of LAMMPS. The simulations were performed with constant Number, Volume, and Energy (NVE) and the temperature was controlled with a Langevin thermostat (Schneider & Stoll, 1978). The simulations contained one type of atom with all particles interacting with a Yukawa potential with the pair-coefficient, (Ze)2, set to unity. The atom mass was was also set to unity. The FCC simulations used 4000 atoms and 2000 atoms were used in the BCC simulations. The melting temperature values are calculated at various inverse screening lengths, with the inverse screening length was varying from κ = 0.5 to κ = 6.45. For simulations with κ > 3 the initial lattice structure is FCC, and BCC below this value. Four different cut-off radii were investigated; rc = 16/κ, 8/κ, 4/κ, and 3.07a. For three of the cases the cut-off radius is dependent on the value of the inverse screening length, while the fourth case sets the cut-off radius to a value based on the particle spacing. Each simulation was heated over 2× 106 steps with a 0.09 time step. In every case the crystal had melted by the end of the simulation. The mean squared displacement (MSD) of the particles was calculated as the temperature increased and is used in determining the melting temper- ature. The melting temperature in this work is defined as the temperature where the largest jump in the mean squared displacement occurs. Another method which could of been used to determine the melting temperature would be to equate the free energy of the liquid and solid phase of the sys- tem. Examples of how the melting temperature is determined are shown in Figure 5.1, which plots the logarithm of the MSD as a function of the reduced temperature for simulations with κ = 0.5 and κ = 1.29. The crystal was found to melt for the κ = 0.5 case at a temperature of T ∗ = 0.0104433. The data for the κ = 1.29 case shows multiple jumps in the MSD, but the 68 5.5. Results: Melting Temperature melting occurred at the largest jump, at a temperature of T ∗ = 0.0106. The reduced temperature from the simulations can be related to both the energy ratio, kbT/Ua and the T̃ parameter from Robbins et al. (1988). The reduced temperature is defined as T ∗ = kbT/Uo, and since Ua = Uo e−λ, T ∗ = kbT Uo = kbT Ua eλ . (5.7) As a result, the energy ratio is simply kbT/Ua = T ∗ eλ. In order to relate T ∗ to T̃ , ma2ω2E is equated to 2/3λ 2Ut, where Ut is the total potential energy (Robbins et al., 1988), thus T̃ = kBT ma2ω2E = 3kBT 2λ2Ut . (5.8) As kbT = T ∗Uo and the total potential energy is the kinetic energy sub- tracted from the total energy (ET ), then in reduced units Ut = E ∗ T − 3/2T ∗. Substituting the expression for the total energy into the above equation for T̃ gives the relationship: T̃ = 3T ∗ 2λ2(ET − 3/2T ∗) , (5.9) since the characteristic energy, Uo = 1 for these simulations. 5.5 Results: Melting Temperature The crystal melting simulations were performed for crystals of inverse screen- ing length κ = 0.5 to κ = 6.45 with cut-off radii rc = 16/κ, 8/κ, 4/κ, and 3.07a, where a is the characteristic spacing of the atoms. The calculated melting temperatures for each of the cut-off radii are presented in Table 5.1, along with the initial structure of the crystal. The entries with asterisks indicate simulations which had mean squared displacements with multiple jumps. The empty entries occur in the cases where the cut-off radius is close to or less than the characteristic particle spacing; in these cases the particles would not interact with any nearby particles. 69 5.5. Results: Melting Temperature l ll lll l lll llll llll l ll lll llll lll lllll llll llllll llllll lll llllll l ll l ll l l l l l l l l l l l l l l ll ll ll ll ll lll ll l lll ll lll l lll l ll ll lll l ll l l ll l lll l ll l l llll lll lll ll llll lllll ll lll lll l ll l lll lll ll l ll lll lll llll lll lll lll lll llllll lll l ll ll lll l lll l l lll ll ll l ll ll l lllll ll ll l ll ll llll lll l ll l lll ll lll lllll ll l ss ss s s s s ssssssssssssss s sss ssssssssssssssss sss s s s s s s s s s s ss ss ss ss ss s ss sss sss ssss sss s ss ss ssss sss ss s sss s s ssss s s sss ss sss s ss ss s sss ss s ss sss s s ssss ss ss ss ss sss s ssss s sss ss ss ss s ssss sss ss ss ssss ss ss s s sss sss ss ss s sss sss sss ss s ss sss ss s sss sss s ssssssssss s ss ss s ssss ss ss sss ss ssss l s Melting Temperature Examples T∗ L og 10 (M S D ) 0.005 0.010 0.015 0.020 0.025 0.030 0.035 0.040 −2 −1 0 1 2 3 κ = 0.5 κ = 1.29 Figure 5.1: Two examples of how the melting temperatures from the simula- tions are determined. The melting temperature is defined as the temperature where the jump in the logarithm of the mean squared displacement (MSD) is the largest. In the example of κ = 0.5 with a rc = 8/κ the melting oc- curs at T ∗ = 0.0104433. In the κ = 1.29 example, which has a rc = 4/κ, there are multiple jumps in the MSD, but the largest occurs at T ∗ = 0.0106, indicating the melting temperature. 70 5.6. Results: Simulation Length and Temperature Each of the melting temperatures calculated at the different inverse screening lengths and cut-off radii are compared to the Robbins et al. (1988) melting temperatures. The calculations at each cut-off radius are compared to the expected melting temperatures in terms of kBT/Ua and T̃ as shown in Figures 5.2 and 5.3, respectively. At small values of λ all the energy ratios in Figure 5.2 appear to agree, except for the cut-off radius of rc = 3.07a. At larger values of λ the calculated melting temperatures are larger than the Robbins et al. (1988) result. In comparing the T̃ melting values in Figure 5.3 differences between the expected result and the simulated results are more apparent. At a cut-off radius of rc = 3.07a there is a large degree of discrepancy from the expected melting temperature for values of λ < 3. This indicates that for κ < 3, a cut-off radius of rc = 3.07a is inappropriate. The melting temperatures calculated using a cut-off radius of rc = 4/κ, like for the cut-off radius of rc = 3.07a, are found to differ from the expected melting temperature for certain values of λ. The melting temperatures cal- culated with the larger cut-off radii are found to be closer to the the expected melting temperature. At larger values of λ the melting temperatures show better agreement with the expected melting temperature at the various cut- off radii. This result is not unexpected considering the results from Chapter 4, where for κ = 4.74, which corresponds to λ = 3.74, a cut-off radius of rc = 3.0 was found to be appropriate. At a κ = λ = 0.8835 a larger cut-off radius of rc = 8.0, or rc = 7.068/κ was required. 5.6 Results: Simulation Length and Temperature An additional investigation was undertaken in order to determine if the calculated melting temperature is sensitive to the rate at which the crystal is heated. This is a test to see if the number of steps used in the simulations are appropriate. For this, one case was considered; a FCC λ = 6 crystal with a cut-off radius of 3.07a, which corresponds to rc = 18.42/κ. With the simulations discussed in Section 5.4, 2 × 106 steps were taken in order melt the crystal. Additional simulations of 1 × 106 and 3 × 106 steps were performed on a FCC crystal with κ = 6.0 and rc = 3.07a in order to compare 71 5.6. Results: Simulation Length and Temperature l ll l l l l l l l l l l l s ss s s s s s s s s s s s n nnn n n n n n : ::: : : : : : : : : : : 6 6 6 6 6 6 6 l s n : 6 λ = κa, a=1 (k b T m el t) /( U a ) 1 2 3 4 5 6 7 0. 05 0. 10 0. 15 0. 20 0. 25 0. 30 0. 35 0. 40 0. 45 16/κ 8/κ 4/κ 3.07a Robbins 1988 Figure 5.2: Comparing the energy ratio at the melting temperature for various cut-off radii to the results from Robbins et al. (1988) (labelled Rob- bins 1988). The label indicates the cut-off radius used in the calculation, rc = 16/κ, 8/κ, 4/κ, or 3.07a. Note that there are only values for rc = 4/κ for λ ≥ 3.23, as above this value of λ the cut-off radius was too short for there to be interactions between particles. 72 5.6. Results: Simulation Length and Temperature l lll l l l l l l l l l l s sss s s s s s s s s s s n n nn n n n n n : : : : : : : : : : : : : : l s n : λ = κa, a=1 T̃ m 1 2 3 4 5 6 7 0. 00 2 0. 00 3 0. 00 4 0. 00 5 0. 00 6 0. 00 7 0. 00 8 0. 00 9 0. 01 0 0. 01 1 0. 01 2 16/κ 8/κ 4/κ 3.07a Robbins 1988 Figure 5.3: Comparing the T̃ melting value of the different cut-off radii melting temperature calculations to the Robbins et al. (1988) result of T̃m = 0.00246 + 0.000274λ. Simulations with smaller values of λ require a longer cut-off radius in order to agree with the expected melting temperature. Note that for a cut-off radius of rc = 4/κ at λ = 3.23 the cut-off radius is rc = 1.23. 73 5.7. Results: Fitting for Melting Temperature Structure κ T∗m T∗m T∗m T∗m rc = 16/κ rc = 8/κ rc = 4/κ rc = 3.07a BCC 0.5 0.0110142 0.0104433 0.0101333* 0.02* BCC 0.71 0.0105465* 0.0107639 0.0125488* 0.017* BCC 0.81 0.0107369 0.0105282 0.00950195* 0.014* BCC 0.92 0.0103464 0.00990401 0.00893095 0.013* BCC 1.08 0.0100772 0.0100009 0.0125845 0.012* BCC 1.29 0.0094542 0.00945241 0.0106* 0.010* BCC 1.61 0.00870435 0.0088535 0.0074933* 0.00889075* BCC 2.15 0.00712177* 0.00719909 0.0074511 0.00734446 FCC 3.23 0.00444924 0.00428558 0.00613765* 0.00439987 FCC 3.74 0.00323955 0.00326031 – 0.00329566 FCC 4.5 0.00206131 0.00205862 – 0.00211135 FCC 5.5 0.0010634 0.00100772 – 0.00106841* FCC 6.0 0.000742646 0.000736148 – 0.000740376 FCC 6.45 0.000538555 0.00063504 – 0.000541194 Table 5.1: Melting temperatures for the various cut-off radii. Those entries with an asterisk (*) beside the temperature are simulations in which there are multiple jumps in the mean squared displacement. These results are also displayed in Figures 5.2 and 5.3. the melting temperature of the three different heating rates. The melting temperature, T ∗m, and the energy ratio, kBT/Ua, for the three simulations are compared in Table 5.2. These calculations indicate that the melting temperature is not strongly dependent on how the fast the crystals were heated. This is a limited range of heating rates and larger range of heating rates may indicate at which rate there is a strong dependence on the number of steps used in the simulation. 5.7 Results: Fitting for Melting Temperature The calculation of the interactions between particles is cut-off at a prescribed radius in order to decrease the time it takes for a simulation to run. In an ideal situation each of the particles would interact with all the other parti- cles; thus, the cut-off radius would be infinite. A method to approximate 74 5.7. Results: Fitting for Melting Temperature Steps T ∗m kBT/Ua 1× 106 0.000753143 0.30383957 2× 106 0.000740376 0.298688996 3× 106 0.000744994 0.33055203 Table 5.2: Dependence of the melting temperature and the energy ratio on how quickly the crystal is heated. These melting temperatures are calculated for a FCC crystal with an inverse screening length of κ = 6.0 and a cut-off radius of 3.01a. the value of the melting temperature at an infinite cut-off radius would be to extrapolate to the value. The melting temperature has been calculated for various cut-off radius values listed in Table 5.1. The intercept of the melting temperature as a function of 1/rc, as well as 1/(rc κ), provides an estimate of the melting temperature at 1/rc = 0 or rc = ∞. The relationship of the melting temperature as a function of 1/(rc κ) is displayed in Figure 5.4. The data in this figure also shows the data sets at cut-off radius regimes used for fitting for the melting temperature at an infinite cut-off radius. The linear fitting results of the data are discussed below in Section 5.7.1 and the fitting of the more difficult data sets are discussed in Section 5.7.2. 5.7.1 Linear Fitting Linear fitting was performed on all of the sets of data, each set is indicated in Figure 5.4. Since the melting temperatures calculated for λ = 0.5 to 1.61 for a cut-off radius of rc = 3.07 disagree from the temperature calculations for the other cut-off radii, see Figure 5.3, these temperatures are not included in the fits. The cut-off radius at 4/κ for λ = 3.23 is also discrepant and is similarly not included in the fitting. The data was fit to a straight line using χ2 minimization, assuming equal errors for all data points (Press et al., 1992). The fits to the melting temperature data are listed in Table 5.3. The intercept of the linear fit of the melting temperature as function of 1/rc is the melting temperature that the crystal would have if all the particles in the simulations were interacting with each other. 75 5.7. Results: Fitting for Melting Temperature l l l s s s n n n : : : 6 6 6 l l l s s s n n n n : : : 6 6 6 l l l ss s n n n : : : l s n : 6 l s n : 6 l s n : 1/κrc T * m 0.05 0.10 0.15 0.20 0.25 0.30 0. 00 2 0. 00 4 0. 00 6 0. 00 8 0. 01 0 0. 01 2 κ=0.5 κ=0.71 κ=0.81 κ=0.92 κ=1.08 κ=1.29 κ=1.61 κ=2.15 κ=3.23 κ=3.74 κ=4.5 κ=5.5 κ=6.0 κ=6.45 Figure 5.4: The melting temperature of the various screening lengths for the different cut-off radius regimes. These are the data sets used for the fitting of the melting temperature at infinite cut-off radius. Where the data intercepts the y-axis indicates the expected melting temperature if all particles interacted with each other. 76 5 .7 . R esu lts: F ittin g for M eltin g T em p eratu re κ intercept inter err slope slope err χ2 goodness 0.5 0.0111692 0.00027227 -0.0087625 0.00329307 4.94208×10−8 1 0.71 0.00965406 0.000441867 0.0157652 0.00376363 1.30166×10−7 1 0.81 0.01125 0.0001993 -0.00841777 0.00148797 2.648×10−8 1 0.92 0.0108329 2.88825×10−5 -0.00824202 0.000189852 5.56131×10−10 1 1.08 0.00878542 0.000895635 0.0133468 0.00501508 5.3478×10−7 1 1.29 0.0088804 0.000376806 0.00507717 0.00176646 9.46544×10−8 1 1.61 0.00864728 7.68895×10−5 0.000807452 0.000336397 2.87792×10−9 1 2.15 0.00701236 6.06399×10−5 0.000842561 0.000174097 5.11456×10−9 1 3.23 0.00461823 9.22479×10−5 -0.000773071 0.000287032 1.70762×10−9 1 3.74 0.00324459 8.2951×10−5 6.01302×10−5 0.000233286 1.5094×10−9 1 4.5 0.0021084 7.33018×10−5 -8.03213×10−5 0.000179266 1.46898×10−9 1 5.5 0.00112122 2.40744×10−6 -0.000165155 4.99476×10−6 2.07394×10−12 1 6.0 0.000745889 3.26158×10−6 -1.2749×10−5 6.27981×10−6 4.24711×10−12 1 6.45 0.000464251 2.0473×10−5 0.000209779 3.69968×10−5 1.82231×10−10 1 Table 5.3: The linear fitting parameters of each of the sets of inverse screening length. The intercept gives the extrapolated melting temperature of the specific inverse screening length as if all the atoms in the simulation interacted with each other. The linear fitting of the data was performed without errors attached to the entries. 77 5.8. Conclusions 5.7.2 Dealing with Hockey Sticks For some of the simulations there is a ‘hook’ in the data, or rather the data set takes on the shape of a hockey stick. In these cases the data was refit using only two of the data points which are at the largest rc values from the data sets. The data sets which were refit in this manner did not have errors calculated with the linear fit, as the fitting routine requires three data points for the error calculation. The refit data along with the resulting intercept and slope are listed in Table 5.4. The final melting temperatures including this refit data, corresponding to an infinite cut-off radius, are compared to the results from Robbins et al. (1988) in Figure 5.5. The melting temperatures are taken from Table 5.4 for the κ values that had to be refit; all other melting temperatures are taken from Table 5.3. The comparison of the extrapolated melting temperature to that of Robbins et al. (1988) – see in Figure 5.5 – shows good agreement. κ intercept slope 1.08 0.0101535 -0.00113034 1.29 0.00945599 -2.2199 3.74 0.00309697 0.00609982 4.5 0.00174493 0.0011249 6.0 0.000725367 4.60792×10−5 6.45 0.000552303 -3.41039×10−5 Table 5.4: Linear fitting parameters from using only two data points from the set. This table only includes data files which were refit due to the hook in the data, other fits are the same as in Table 5.3. This fit was performed without attaching errors to the data. 5.8 Conclusions In order to test that the LAMMPS simulations were performing as expected and also to confirm the required cut-off radius values found in Chapter 4, the melting temperature was calculated for different values of inverse screening length. The cut-off radius value was also varied for each value of inverse 78 5.8. Conclusions l l ll l l l l l l l l l l s s s s s s sl s λ = κa, a=1 (k b T m el t) /( U a ) 0 1 2 3 4 5 6 7 0. 05 0. 10 0. 15 0. 20 0. 25 0. 30 0. 35 Fitted Data Robbins 1988 Figure 5.5: Comparing the the energy ratio of the fitted melting tempera- tures found from the LAMMPS to the results of Robbins et al. (1988). The two are in good agreement, especially at the smaller values of λ. 79 5.8. Conclusions screening length. The crystal structure used in the simulations depended on the value of inverse screening length: for simulations with κ > 3 the initial structure used was FCC and for inverse screening lengths below, the structure used was BCC. Each initial crystal was heated and the melting temperature was determined. The melting temperature in these simulations was defined as the temperature at which there is the largest jump in the MSD. In order to approximate the conditions where all the particles interact with each other, the cut-off radius was varied, the value of the melting temperature at infinie cut-off radius can be found by extrapolating the data set for a specific screening length. The melting temperature found in the extrapolation was then compared to the results from Robbins et al. (1988). The simulation and Robbins et al. (1988) results were found to be in good agreement, especially at smaller values of inverse screening lengths. The determination of the melting temperature can be affected by the simulation parameters and can be measured with different methods. The effect of the rate at which the crystal was heated on the melting temperature was investigated by changing the number of steps for the heating to occur over. The range of rates used in this chapter did not indicate a dependence of the melting temperature with the heating rate. The heating rates, with a larger range of values than investigated, may start to indicate the number of simulation steps at which the heating rate becomes important for the melt- ing temperature. The inverse screening length and the corresponding cut-off radius used may also change the dynamics of the system, and would need to be further investigated in the future. The melting temperature could also be calculated in a different manner, such as equating the free energy of the liq- uid and solid phases of the simulation. The method used in this chapter for determine the melting temperature via the MSD could be double checked by recrystallizing the melted material and calculating the crystallization temperature. If there is no hysteresis than the melting and crystallization temperatures would be in agreement. With these considerations, from the comparisons fo the results of the melting temperatures calculated in this chapter the Robbins et al. (1988), indicates that the melting criterion be- tween the two is not very different. Comparing the melting temperatures 80 5.8. Conclusions at the different cut-off radii indicated that at small values of λ, or small inverse screening lengths, a larger cut-off radius is required than for those simulations with larger inverse screening lengths. With these simulations the results from Chapter 4 are confirmed in terms of the cut-off radius required at a value of κ and the method of implementing the molecular dynamics using LAMMPS is also confirmed. 81 Chapter 6 Mechanical Properties of Pure FCC and BCC Crystals 6.1 Introduction The mechanical properties of a crystal, such as the breaking strain and shear modulus, can be determined using molecular dynamics simulations. In the context of understanding a neutron star crust, molecular dynam- ics simulations have previously been performed to understand the breaking strain (Horowitz & Kadau, 2009) and also the breaking stress (Chugunov & Horowitz, 2010) of both an accreted neutron star crust composition as well as pure BCC lattices. A pure crystal was simulated in Horowitz & Kadau (2009) as well as simulations which contained impurities, defects and grain boundaries. The various types of defects in the simulated samples were found to only slightly reduce the measured breaking strain to a value around 0.1. In Chugunov & Horowitz (2010) the maximum value of the stress-strain relationship, the breaking stress, was explored. For these sim- ulations a sample containing only one type of ion was calculated. With this sample the breaking stress at various inverse temperatures and strain rates was examined. Within the timescales simulated in Chugunov & Horowitz (2010) the breaking stress was found to have a strong dependence on the temperatures. In this chapter the tools which are used to examine the neutron star crust are developed and applied to pure FCC and BCC crystals. In Chapter 7 the techniques discussed here are used to characterize the nature of a non- accreting neutron star crust. From these simulations the shear modulus and 82 6.2. Applying Shear breaking strain are calculated by applying a shear to the sample crystal. This chapter is divided in the following manner: the process for applying shear in a simulation, the calculation of the thermodynamic and mechanical properties, and the simulation results of applying a simple shear to pure FCC and BCC crystals. Studies were also conducted to investigate additional properties of the crystals through simulations looking for a second break, running the simulation forward and then reversing the direction of shear, and the temperature dependence on the shear modulus and breaking strain. The simulations discussed in this chapter are primarily used to create a molecular dynamics tool kit which will be applied to the neutron star crust material in Chapter 7. 6.2 Applying Shear In Horowitz & Kadau (2009) two different methods were used to introduce the shear to the simulations. The first method was to deform the periodic boundaries of the simulation box. The second method was to create three layers of atoms, top and bottom boundary layers where the atoms do not interact and an active region, the two boundary move opposite to each other to introduce the shear (see Horstemeyer et al., 2001). For the simulations in this chapter a shear was introduced using LAMMPS version 21 Jul 2010 by deforming the simulation box containing the crystal lattice, either FCC or BCC. The particles in the box continued to interact via the Yukawa poten- tial while the deformation was applied. For each simulation calculations of the thermodynamic and mechanical properties were performed. The defor- mation, simulation parameters, and calculations are discussed fully below. 6.2.1 Deformation In order to deform a crystal using LAMMPS the command fix deform was invoked. With this command the box was deformed in the X-direction. A velocity was also defined which specifies the rate at wich the tilt of the box, or strain, was changed. This is the same as specifying a strain rate with 83 6.2. Applying Shear ṡ = v/l, where ṡ is the strain rate, v is the velocity and l is the original box length in the direction perpendicular to the direction of shear. A schematic of the box deformation is shown in Figure 6.1, where the box is deformed in the x-direction and the ratio ∆x/l is the strain. l ∆x Figure 6.1: The box is deformed in the x-direction. The strain is calculated by taking the ratio of ∆x to l, where ∆x is the amount the top of the box has moved, and l is the length of the box in the y-direction. Strain measures the degree of deformation. 6.2.2 Simulation Parameters In both FCC and BCC cases pure crystals were simulated and the mass of all atoms in the system was set to unity. The pair coefficient, (Ze)2, was also set to unity for both cases. An inverse screening length of κ = 3.74 was used in the FCC simulations. With this κ value a cut-off radius of rc = 3.07a is appropriate. In the BCC simulations a κ = 0.8835 was used with a cut-off radius of rc = 9.055, or rc = 8/κ. Before the box was deformed, the simulation first underwent a stage of thermal equilibration. This was done by giving the particles an initial velocity, which was below the melting point of the crystal. The simulations were performed with NVE time integration – constant number, volume and energy. Thermal equilibration of the system was performed for 200 steps, with a time step size of 0.01. Constant temperature during the deformation 84 6.2. Applying Shear was maintained by rescaling the temperature to a specified value when it exceeded set bounds in a manner similar to Horowitz et al. (2007) rather than a Langevin thermostat because the temperature rescaling allows for direct and rigid temperature control. 6.2.3 Calculation of Material Properties Temperature and Energy The system temperature and energy were both calculated within the LAMMPS simulations. The temperature of the system of particles was calculated as: 3/2N k T = KE, where KE is the sum of the kinetic energy of the particles, KE = ∑ 1/2mv2. Constant temperature was maintained and controlled through temperature re-scaling. The deformation of the simulation box contributes an additional velocity to the simulation, where particles near the top of the box have a larger additional velocity than the particles near the bottom of the box. This is a non-thermal contribution and was subtracted from the temperature measurement. The total energy of the simulations is the summation of the kinetic and potential energies of the system. Pressure and Stress The pressure and stress tensors were calculated for all particles in the sim- ulation within LAMMPS. The pressure tensor includes a kinetic energy term and the calculation of these values did not include the non-thermal contri- bution due to deformation of the box. Pressure and stress tensors differ in the direction: pressure has the opposite direction of the stress, or the pressure is the negative of the stress. Along with the pressure tensor calcu- lation, the total pressure of the system of particles was also calculated. In order to calculate the stress tensor, the stress tensor for each atom was first calculated. The global stress tensor of the system was then calculated by taking a summation of the per-atom values which was then divided by the volume of the simulation. 85 6.3. Breaking Strain of Pure Crystals Stress-Strain Relationship The amount of strain is a measure of the degree of deformation of the simu- lation box. Strain is measured by taking the ratio between the displacement in the direction of applied shear, in this case the displacement is in the X- direction, and the original length of the box. Thus, the strain is measured as ∆x/l, as displayed schematically in Figure 6.1. In order to ensure the sim- ulation was not being overstrained, the rate of strain, ṡ, is chosen such that ṡ/ωp < 1, where ωp is the plasma frequency. The plasma frequency of a ma- terial is: ωp = (4piZ 2e2ni/mi) 1/2. In the simulations described here, mass, mi, number density, ni, and charge, Z 2e2 are all set to unity for pure FCC and BCC simulations. Thus, in the simulation units, the plasma frequency is simply ωp = (4pi) 1/2. Note that in Chugunov & Horowitz (2010) ratios of ṡ/ωp between 10 −4 and 10−7 were used in their work. The strain rate is the velocity of the shear applied divided by the length of the box: v/l, thus velocity of the shear in the simulations is kept much below v = l (4pi)1/2. The stress-strain relationship, the amount of stress or pressure compared to the degree of strain, can give an indication of how strong the material is and when the material yields. A schematic of the expected stress-strain relationship is shown in Figure 6.2. The shear modulus is the slope of the linear region. At the point where the stress-strain relationship is no longer linear in nature, the material has yielded to the applied stress. The breaking strain in this work is the strain at which the material has yielded to the applied shear forces. The breaking stress is the amount of stress at the point of yielding. 6.3 Breaking Strain of Pure Crystals The breaking strain and shear modulus of both pure FCC and BCC crystals were calculated. The effects of the simulation size was first investigated by examining five different FCC and BCC box sizes which all had a 20× 10−6 strain rate applied. A shear was applied to the FCC and BCC simulations until a strain of ∆x/l = 0.2, 20% deformation, was reached. To determine 86 6.3. Breaking Strain of Pure Crystals Yield St re ss Strain Linear Figure 6.2: A schematic of a stress-strain relationship. The shear modulus of the material is the slope of the linear region. The breaking strain is the point at which there is a change from the linear increase of stress with strain, this is labeled ‘Yield’. the breaking strain and shear modulus a range of strain rates were applied to the crystals in order to determine if the simulations had converged to the same values of breaking strain and shear modulus. 6.3.1 Box Size Effects The effect which the size of the simulation box has on the stress-strain results was tested by comparing simulations with unit cell volumes of 10× 10× 10, 12 × 12 × 12, 20 × 20 × 20, 25 × 25 × 25, and 35 × 35 × 35, for both the FCC and BCC structures. The corresponding number of particles used in the simulations are listed in Table 3.2. The simulations with smaller unit cell volumes correspond to fewer particles in the simulation. A strain rate of 20 × 10−6 was applied in all cases. The stress-strain relationships of the 87 6.3. Breaking Strain of Pure Crystals five sizes of the FCC and BCC crystals are compared in Figures 6.3 and 6.4, respectively. For the five sizes of the FCC and BCC crystal simulations the shear modulus is found to be the same. The figures indicate differences in the breaking strain for the different sizes. The FCC crystals with the largest box sizes yield before the smaller boxes. The BCC crystal with the largest box sizes also yield first, but the two smaller box sizes yield at nearly the same point. The difference in the breaking strain is likely due to the number of particles in the simulation. Simulations of pure crystals with a simulation box size of 25× 25× 25 unit cells are presented next. 6.3.2 Pure FCC Crystal The FCC simulations with various applied strain rates contained 62500 atoms. Strain rates of 20, 5, 2.5, and 1.25×10−6 were applied to the crys- tal. The breaking strain and shear modulus determined for each strain rate are listed in Table 6.1, as well as the stress-strain relationships for the different strain rates are shown in Figure 6.5. This figure indicates that the shear modulus for the four strain rates are simular, but the breaking strain has a dependence on the strain rate applied. The largest strain rate of ṡ = 20 × 10−6 is found to yield after the smaller strain rates, whereas the smaller strain rates yield at approximately the same strain. The total energy as function of strain is displayed in Figure 6.6 for an FCC crystal. This plot displays the same yielding events, with the smaller scale plastic events occurring after the crystal fails the first time. The breaking strain is approximately 0.1 and the shear modulus is approximately 0.086 in reduced pressure units in the FCC case. 6.3.3 Pure BCC Crystal The BCC simulations contained 31250 particles. This crystal case represents a perfect 56Fe BCC crystal. Strain rates of 20, 5, and 2.5×10−6 were applied to the perfect crystal. The breaking strain and shear modulus are listed in Table 6.2 for the three different strain rates, with the stress-strain relation- ships compared in Figure 6.7. The shear modulus is the same for the three 88 6.3. Breaking Strain of Pure Crystals FCC Sizes Strain σ x y [P re ss u re ] 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.200 .0 00 0. 00 1 0. 00 2 0. 00 3 0. 00 4 0. 00 5 0. 00 6 0. 00 7 0. 00 8 0. 00 9 10x10x10 12x12x12 20x20x20 25x25x25 35x35x35 Figure 6.3: The stress-strain relationship of a pure FCC crystal comparing five different sizes of simulations. The number of particles range from 4000 in the smallest simulation box to 171500 in the largest, the particle numbers are listed in Table 3.2. For the same applied strain rate, ṡ = 20× 10−6, the five sizes share the same shear modulus, but different breaking strains. The larger simulations yield before the smaller. 89 6.3. Breaking Strain of Pure Crystals BCC Sizes Strain σ x y [P re ss u re ] 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20 0. 00 5 0. 01 0 0. 01 5 0. 02 0 0. 02 5 0. 03 0 10x10x10 12x12x12 20x20x20 25x25x25 35x35x35 Figure 6.4: The stress-strain relationship of a pure BCC crystal comparing five different sizes of simulations. In the smallest simulation box there are 2000 particles and 85750 particles in the largest simulation box. The range of the number of particles in the simulation boxes are listed in Table 3.2. For the same applied strain rate of ṡ = 20 × 10−6 the five sizes of simulations share the same shear modulus, but the breaking strain is different. The larger simulation boxes are found to yield before the smaller. 90 6.3. Breaking Strain of Pure Crystals FCC Strain σ x y [P re ss u re ] 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.200 .0 00 0. 00 1 0. 00 2 0. 00 3 0. 00 4 0. 00 5 0. 00 6 0. 00 7 0. 00 8 20 5 2.5 1.25 Figure 6.5: The stress-strain relationship for a 62500 particle, pure FCC crystal with an inverse screening length of κ = 3.74 with applied strain rates of 20, 5, 2.5, and 1.25 ×10−6. The shear modulus is the same for all the strain rates and the breaking strain is similar at small strain rates. 91 6.3. Breaking Strain of Pure Crystals FCC Strain T ot al E n er gy 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20 0. 09 38 0. 09 39 0. 09 40 0. 09 41 0. 09 42 0. 09 43 20 5 2.5 1.25 Figure 6.6: The energy corresponding the the deformation of a pure FCC crystal with an inverse screening length of κ = 3.74 for the strain rates of 20, 5, 2.5, and 1.25 ×10−6. After the initial breaking point, at ∼0.1 smaller yield events occur. 92 6.4. Second Break Strain Rate ṡ/ωp Breaking Shear (×10−6) (×10−7) Strain Modulus 20 56.42 0.104669 0.08570 5 3.979 0.0998332 0.08569 2.5 1.989 0.100611 0.08576 1.25 0.995 0.0998332 0.08573 Table 6.1: The breaking strain and shear modulus for a 62500 particle, pure FCC crystal, with κ = 3.74 for strain rates of 20, 5, 2.5, and 1.25×10−6. The ṡ/ωp ratio is below unity for all strain rates. values of ṡ and the breaking strain converges for smaller values of the strain rate. The breaking strain of a perfect BCC crystal was found to be approx- imately 0.11 and the shear modulus was determined to be approximately 0.27 in reduced pressure units. Strain Rate ṡ/ωp Breaking Shear (×10−6) (×10−7) Strain Modulus 20 56.42 0.114508 0.27217 5 3.979 0.112896 0.27219 2.5 1.989 0.111395 0.27220 Table 6.2: The breaking strain and shear modulus for a pure BCC crystal, with κ = 0.8835 for the three different strain rates. The value of ṡ/ωp was kept below unity for all strain rates investigated. 6.4 Second Break After the first break in some of the breaking strain simulations the pressure began to once again increase linearly with strain. This is likely to be an indication that the crystal is undergoing a second major yielding event. In order to investigate the possibility that the system undergoes a second yielding event, or second break, the simulations are run with a strain rate of 20 × 10−6 for a total strain of ∆x/l = 0.4, twice as much strain as the previous simulations. The FCC and BCC crystal cases were both examined 93 6.4. Second Break BCCFe Strain σ x y [P re ss u re ] 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20 0. 00 5 0. 01 0 0. 01 5 0. 02 0 0. 02 5 20 5 2.5 Figure 6.7: The stress-strain relationship of a 31250 particle pure BCC crystal with an inverse screening length of κ = 0.8835 with applied strain rates of 20, 5, and 2.5×10−6. The shear modulus for the three strain rates is the same at ∼0.27 reduced pressure units. The breaking strain is ∼0.11. 94 6.5. Forward and Reverse for the occurrence of a second break. The stress-strain relationship for the FCC crystal is shown in Figure 6.8. The FCC simulation did not show any evidence for a second break occurring for the strain rate of ṡ = 20 × 10−6. The strain-strain relationship for the BCC crystal is presented in Figure 6.9. The BCC simulations did show an indication of a second break occurring. The shear modulus for the first and second break were different with a shear modulus of µ∼ 0.27 before the crystal first yields and a shear modulus of µ∼ 0.11 before the second break. The shear modulus of the second break was measured to occur between a strain of 0.133 and 0.186. 6.5 Forward and Reverse Simulations were also performed where the direction of applied shear was changed after the material yielded. These simulations were conducted to investigate the effect changing the direction of shear has on the shear mod- ulus of the crystal and if the yield event is reversible. If the yielding event was reversible, the same shear modulus would be found in both directions. These simulations were performed on both the FCC and the BCC crystals. The crystals were deformed at a rate of 20 × 10−6 until the fracturing oc- curred and then the same rate of strain was applied in the reverse direction. For plotting purposes the strain measurements in the figures continues to increase after reversing the shear direction as (step size×∆t × 20 × 10−6). For the FCC crystal, a strain of 0.106 was reached and then the deforma- tion of the crystal was reversed. The stress-strain relationship of the FCC crystal is shown in Figure 6.10 with both the forward and reverse directions. The shear modulus changes from µ∼ 0.087 in the forward deformation, to a shear modulus of µ∼ 0.079 when the deformation is reversed. The shear modulus in the reverse direction was measured between a strain of 0.165 to 0.188 for the FCC crystal. In the BCC case, the initial shear modulus was measured to be µ∼ 0.27. The box deformation was reversed at a strain of 0.12, and a shear modulus of µ∼ 0.17 was measured between a strain of 0.175 and 0.2115 in the reverse direction. The stress-strain relationship of the BCC crystal is shown in Figure 6.11. As the shear modulus is not the 95 6.5. Forward and Reverse FCC: Second Break Strain σ x y [P re ss u re ] 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.400 .0 00 0. 00 1 0. 00 2 0. 00 3 0. 00 4 0. 00 5 0. 00 6 0. 00 7 0. 00 8 Figure 6.8: The stress-strain relationship of a 62500 particle, pure FCC crystal with an inverse screening length of κ = 3.74 for a strain rate of 20 × 10−6. This simulation was run to a larger deformation in order to investigate a second yielding event. There is no indication that a second break occurred at this strain rate for the FCC crystal. 96 6.5. Forward and Reverse BCCFe: Second Break Strain σ x y [P re ss u re ] 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0. 00 5 0. 01 0 0. 01 5 0. 02 0 0. 02 5 Figure 6.9: The stress-strain relationship of a 31250 particle, pure BCC crystal with an inverse screening length of κ = 0.8835 for a strain rate of 20 × 10−6. The simulation box was deformed to a strain of ∆x/l = 0.4 in order to investigate the occurrence of a second material failure. A second break occurred for the BCC crystal between a strain of 0.133 to 0.186 with a shear modulus of µ ∼ 0.11 in reduced pressure units. 97 6.6. Temperature Dependence same in both directions, the yielding event is not reversible for either case of crystal structure. 6.6 Temperature Dependence The temperature dependence of the shear modulus and breaking strain was also investigated with the molecular dynamics simulations. Pure FCC and BCC crystals were run at hotter and colder temperatures relative to the original simulations. The temperature for the hotter simulation was held at 0.002, and the temperature of the colder simulations were held at 0.0005. The temperature of the original simulation was held at 0.001 in reduced temperature units. All these temperatures are below the melting tempera- tures of the two crystal types. For the FCC crystal the melting temperature is approximately 0.003 and the melting temperature is approximately 0.01 for the BCC crystal. All temperature simulations were run at a strain rate of 20 × 10−6. The stress-strain relationship for all three temperatures are compared in Figure 6.12 for a FCC crystal and in Figure 6.13 for a BCC crystal. The simulations for both FCC and BCC crystals have different shear moduli and breaking strains at each temperature, as listed in Table 6.3. The BCC simulations indicate small differences for the temperature dependence. In both the FCC and BCC simulations the shear modulus was largest for colder simulations. The temperature dependence on the breaking strain for FCC and BCC crystals are different; for FCC crystal the hottest simulation yielded first, whereas the coldest simulations in the BCC case yielded to the shear first. 6.7 Conclusions A molecular dynamics toolkit was developed in this chapter for understand- ing neutron star crust material and applied to two sample cases: a pure FCC and a pure BCC crystal. The crystal structures within the simula- tion boxes were deformed with the molecular dynamics simulations in order to determine the shear modulus and breaking strain of the crystals. For 98 6.7. Conclusions FCC: Forward and Back Strain σ x y [P re ss u re ] 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20 −0 .0 02 0. 00 0 0. 00 2 0. 00 4 0. 00 6 0. 00 8 Flipped Figure 6.10: The stress-strain relationship of a pure FCC crystal with an inverse screening length of κ = 3.74 and a strain rate of 20×10−6. After the material yielded, the box deformation is reversed. The line labelled ‘Flipped’ indicates the stress applied in the reverse direction. The shear modulus in the forward direction was ∼ 0.087 and in the reverse direction the shear modulus was ∼ 0.079 measured between a strain of 0.165 and 0.188. 99 6.7. Conclusions BCCFe: Forward and Back Strain σ x y [P re ss u re ] 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20 0.22 0.24 −0 .0 05 0. 00 0 0. 00 5 0. 01 0 0. 01 5 0. 02 0 0. 02 5 0. 03 0 Flipped Figure 6.11: The stress-strain relationship of a pure BCC crystal with an inverse screening length of κ = 0.8835 and a strain rate of 20 × 10−6. The deformation of the box was reversed after the material yields. In the forward direction the shear modulus was µ ∼ 0.27, after the box deformation was reversed the shear modulus is µ ∼ 0.17. The shear modulus was measured between a strain of 0.175 and 0.2115. 100 6.7. Conclusions FCC: Temperature Strain σ x y [P re ss u re ] 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.200 .0 00 0. 00 1 0. 00 2 0. 00 3 0. 00 4 0. 00 5 0. 00 6 0. 00 7 0. 00 8 0. 00 9 0. 01 0 T = 0.002 T = 0.001 T = 0.0005 Figure 6.12: The stress-strain relationships for FCC crystals at different tem- peratures. The melting temperature for the FCC crystal is approximately 0.003. All simulations had a strain rate of ṡ = 20× 10−6 applied. The shear moduli differ for all three temperatures with the colder simulations having the largest shear modulus. The hottest simulation broke first. 101 6.7. Conclusions BCCFe: Temperature Strain σ x y [P re ss u re ] 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20 0. 00 0 0. 00 5 0. 01 0 0. 01 5 0. 02 0 0. 02 5 0. 03 0 T = 0.002 T = 0.001 T = 0.0005 Figure 6.13: The stress-strain relationship for BCC crystals at different temperatures. At a temperature of approximately 0.01 the BCC crystal melts. A strain rate of ṡ = 20 × 10−6 was applied in all three cases. For each temperature the shear modulus and breaking strain differ, though T = 0.001 and T= 0.0005 were found to be very similar. The shear modulus of the colder simulation was slightly larger then the warmer simulations. The colder simulation breaks before the hotter simulations, unlike the FCC simulations. 102 6.7. Conclusions Structure Temperature Temperature Breaking Shear T ∗ K Strain Modulus FCC 0.002 6.9×108 0.0738188 0.05185 0.001 3.5×108 0.104669 0.08570 0.0005 1.7×108 0.11673 0.08973 BCC 0.002 6.5×108 0.11712 0.26069 0.001 3.5×108 0.114508 0.27217 0.0005 1.7×108 0.112285 0.27821 Table 6.3: The temperature, breaking strain and shear modulus for a pure FCC and BCC crystal. The temperature and shear modulus are in reduced units. Note that the FCC crystal melts at T ∗∼0.003 and the BCC crystal at T ∗∼0.01. The physical temperature is at the base of the crust at ρ = 1014 g/cm2. For a FCC crystal the difference with temperature for both the breaking strain and the shear modulus was greater than those for the BCC simulations. simulations run at a temperature of T ∗ = 0.001 the FCC crystal is found to have a shear modulus of µ∼ 0.086 and a breaking strain of ∼ 0.1. The BCC crystal at this same temperature had a shear modulus of µ∼ 0.27 and yielded to the applied shear at ∼ 0.11. Only one crystal orientation was investigated in all the various simulations, the same shear force applied on a different orientation of the crystal could result in a different stress-strain relationship. Additional simulations were performed in order to characterize the ma- terials. Both simulations were deformed to a strain of ∆x/l = 0.4, twice as much as the first set of simulations, at a strain rate of ṡ = 20 × 10−6 in order to investigate the possibility of a second yielding event. The FCC crystals did not show any indications of a second break, but the BCC crys- tal did. The shear modulus for the second break of the BCC lattice had a different shear modulus compared to the first break: µ∼ 0.11 for the sec- ond break versus µ∼ 0.27 for the first break. Simulations were performed where the crystal was deformed to the point of yielding, at which point the deformation was reversed. These simulations showed a change in the shear modulus when the simulation was reversed, indicating that the system was 103 6.7. Conclusions not reversible after the yielding event occurred. With a rate of strain of ṡ = 20 × 10−6 the FCC forward shear modulus was found to be µ∼ 0.087 but µ∼ 0.079 on the reverse direction. For the BCC case, with a strain rate of ṡ = 20 × 10−6, the shear modulus was µ∼ 0.27 in the forward direction and µ∼ 0.17 in the reverse direction. The temperature dependence on the shear modulus and the breaking strain was tested by running simulations at the strain rate of ṡ = 20× 10−6 but at three different temperatures, T ∗ = 0.002, 0.001, and 0.0005, for both FCC and BCC crystals. Note that at T ∗∼ 0.003 the FCC crystal melts and at T ∗∼ 0.01 the BCC crystal melts. The shear modulus was larger for the cooler simulations for both the FCC and BCC cases, though in the BCC case it is not significant as found in the FCC simulations. The shear modulus was expected to decrease linearly with increasing temperature (Nadal & Le Poac, 2003). The breaking strain in the FCC case was larger at colder temperatures. The BCC crystal was found to fail earlier with cooler temperatures. The differences in the breaking strain are small as the material cools and this could have implications for SGR burst activity. Due to the small changes in the breaking strain with temperature for these perfect crystals the bursting activity would not be expected to change for a SGR as the crust cools. These different types of simulations will be used on the impure crystals which represent the different neutron star crust compostions in order to probe the mechanical properties of the crust. The pure BCC case from this chapter will be compared to the impure cases in the next chapter. The breaking strain for the BCC crystal simulation of ∼ 0.11 is in agreement with the results found by Horowitz & Kadau (2009). Continuing to deform the simulation box past the first yielding event and observing that a second yielding event can occur is important to do further studies on, as it may have implications on the energetics of subsequent SGR bursts. 104 Chapter 7 Mechanical Properties of Neutron Star Crusts 7.1 Introduction With the tool kit developed using the pure crystal cases from Chapter 6 the properties of impure crystals are now investigated. The impure crystals used are representative of different neutron star crust compositions. The crystals have a perfect BCC structure, which sets an upper limit on the shear modulus and breaking strain of the neutron star crust. The addition of defects to the perfect crystal would be expected to decrease the breaking strain and shear modulus, though Horowitz & Kadau (2009) found this to only be a modest effect in the case of the composition of an accreting neutron star crust. Also in the case of the accreted neutron star crust material chemical phase separation was found to occur, where the liquid ocean was enriched with elements with low atomic numbers (Horowitz et al., 2007). In the case of a non-accreting neutron star crust, which is the focus of this chapter, the crustal composition is not pure iron, but includes some impurities. There are three cases for which the crustal abundances were calculated for a non-accreting neutron star. The abundances used for the crusts were calculated using the reaction network torch, as discussed in Chapter 2. The abundance calculations were based on the cooling curves of three cases of non-accreting neutron stars: a neutron star which cooled via the modified Urca process, a neutron star with a thick crust, and a neutron star with a thin crust. Each of the differently cooled neutron star crusts resulted in 105 7.2. Simulation Set-up different compositions, these crustal compositions are listed in Table 2.1. For each of the three cases the melting temperature, occurrence of chemical separation and stress-strain relationship were investigated. The simulations include neutron star crust compositions which have previously not yet been explored. As well, simulations were run where the box is strained beyond the first yielding event and also where the strain is reversed, these types of simulations have also not been previously studied. The stress-strain rela- tionships of these impure crystals were also compared to the pure crystal from Chapter 6. 7.2 Simulation Set-up In order to add impurities to the simulation, a data file was created for a perfect BCC lattice structure and the different atoms were randomly as- signed to lattice locations, weighted by the number abundance of each of the isotopes. A FCC structure was not considered since the inverse screening length for the crustal compositions were appropriate for a BCC structure. Each of the simulations contained 31250 atoms. The number of each type of atom was determined from the mass fractions of the isotopes. These mass fractions were calculated by torch, as discussed in Chapter 2. The nuclear reaction network software torch produced abundances for 489 isotopes, but only the top three of these were used in the simulations. As a result the mass fractions from torch were rescaled so that the mass fractions of the three isotopes add up to unity. For a neutron star cooling via the modified Urca process the top three most abundant isotopes were found to be 56Fe, 54Fe, and 60Ni. The three most abundant isotopes for a neutron star with a thick crust were 56Fe, 60Ni, and 52Cr. A thin crust had 54Fe, 58Ni, 56Fe as the three most abundant isotopes. The fraction of the number of types was found by starting with mass fraction: Xi = ρi/ρ, where ρi = nimuAi and ni = Ni/V , which is the number density. Substituting the last two expressions into the equation for 106 7.2. Simulation Set-up Xi gives: Xi = ρi ρ = NimuAi ρ V , (7.1) and ρ = Mtotal/V = N 〈A〉mu/V , then Xi is: Xi = NimuAiV V N 〈A〉mu . (7.2) Re-arranging and solving for the number fraction gives: Ni N = Xi 〈A〉 Ai . (7.3) The average mass is given by 〈A〉, the definition of which is described below. The total number of isotopes in the simulation is the sum of the different types of atoms: Ni +NJ +Nk = N ; using the equation for number fraction, the total number is given as N〈A〉(Xi/Ai + Xj/Aj + Xk/Ak) = N . Thus, the average mass is: 1 〈A〉 = Xi Ai + Xj Aj + Xk Ak . (7.4) The number density of electrons is given as ne = 〈Z〉nion, thus the average charge is the ratio of the two number densities: ne/nion. The total number of electrons is Ne,total = Ne,i +Ne,j +Ne,k. As the number density of an ion type is ni = Xi ρ/(muAi), the number density of electrons for a type of atom is: ne,i = ZiXi ρ muAi . (7.5) The number density is n = N/V and the density is ρ = mtotal/V , now: ne,i = Ne,i V = ZiXi muAi mtotal V = ZiXi muAi N mu 〈A〉 V , (7.6) re-arranging to solve for the number of electrons results in: Ne,i = ZiXiN 〈A〉 Ai . (7.7) 107 7.2. Simulation Set-up As a result the total number of electrons is then: Ne = N〈A〉 ( ZiXi Ai + ZjXj Aj + ZkXk Ak ) . (7.8) The average charge is found from the ratio of electron density to isotope density, 〈Z〉 = ne/nion = Ne/Nion, thus: 〈Z〉 = 〈A〉 ( ZiXi Ai + ZjXj Aj + ZkXk Ak ) , (7.9) which was used as the charge in the calculation of the inverse screening length, κ = (24α3/2 pi1/2 〈Z〉)1/3. The average charge of the isotopes dictate the inverse screening length of the system. The degree of impurity of a sample can be quantified by the impurity factor, Qimp. The impurity factor is Qimp = 〈Z2〉 − 〈Z〉2, or (Itoh & Kohyama, 1993): Qimp = n −1 ion ∑ i ni (Zi − 〈Z〉)2, (7.10) where Qimp = 0 for a sample without impurities. The larger the value of Qimp the larger the amount of impurities in the sample. The top three most abundant isotopes of 489 calculated from torch were used for the three different cooling curve cases: the modified Urca, a thick crust, and a thin crust. Table 7.1 summarizes the top three isotopes for each method of cooling along with the original mass fraction of each isotope from the torch calculations. The rescaled mass fraction, the number fraction of each isotope, and the mass value used in the molecular dynamics simulations is also listed in the table. All the particles in the simulation interact via a Yukawa potential. The parameters used for determining the interactions between particles for each composition are listed in Table 7.2, this includes 〈Z〉, which is used to cal- culate the inverse screening length κ, the cut-off radius, and the Qimp pa- rameter are also included in the table. The pair coefficients, (Ze)2 for the pair interactions of each of the crustal compositions are listed in Table 7.3 for the modified Urca crust, Table 7.4 for the thick crust, and Table 7.5 for 108 7.3. Melting Simulations Isotope A Z Mass Fraction Rescaled X Number Fraction Mass Modified Urca Composition 56Fe 56 26 0.559 0.6014 0.5961 1.0 54Fe 54 26 0.3187 0.3429 0.3524 0.964 60Ni 60 28 0.05175 0.0557 0.0515 1.071 Thick Crust Composition 56Fe 56 26 0.9286 0.9429 0.9431 1.0 60Ni 60 28 0.03138 0.0319 0.0298 1.071 52Cr 52 24 0.02485 0.0252 0.0271 0.929 Thin Crust Composition 54Fe 54 26 0.6477 0.6675 0.6808 0.964 58Ni 58 28 0.2209 0.2276 0.2161 1.036 56Fe 56 26 0.1018 0.1049 0.1031 1.0 Table 7.1: The top three isotopes from the torch calculations for three different cooling curves: modified Urca, thick crust, and thin crust cool- ing. The original mass fraction from the torch calculations, as well as the rescaled mass fraction (X) are included. The number fractions determine the number of each type of isotope used in the molecular dynamics simulations. The mass value is a scaled mass for the molecular dynamics simulations, where the mass is unity for A = 56. the thin crust. The value of the pair coefficient was scaled depending on the value of Z. In the case where Z = 26 the pair coefficient (Ze)2 = 1. 7.3 Melting Simulations In each of the crustal composition cases the impure BCC crystals were sub- jected to an increase of temperature and then cooled down to the starting temperature. These simulations allowed for a calculation of the melting tem- perature and to confirm the structure of the system after re-crystallization. The crystals were heated from T ∗ = 0.001 to a value of T ∗ = 0.02 over 106 steps. Note that a pure BCC crystal was found to melt at T ∗∼0.01 in Chapter 6. Once melted, the liquid was allowed equilibrate by fixing the temperature at T ∗ = 0.02 for an additional 106 steps. During the equilib- rium portion of the simulation the atoms were able to cross the length of the 109 7.3. Melting Simulations 〈A〉 〈Z〉 κ rc Qimp mod Urca 55.501 26.103 0.8847 9.043 0.195 thick 56.011 26.005 0.8836 9.054 0.228 thin 55.071 26.432 0.8884 9.005 0.678 Table 7.2: The inverse screening length and cut-off radius used in the simu- lations. The parameters 〈A〉 and 〈Z〉 were used in the calculation simulation parameters. The Qimp value indicates the degree of impurity of the sample. type 1 2 3 1 1.0 1.0 1.0769 2 – 1.0 1.0769 3 – – 1.1598 Table 7.3: Pair coefficients, (Ze)2, of the Yukawa potential for the modified Urca composition. Type 1 atoms are 56-Fe, type 2 atoms are 54-Fe and the type 3 atoms represent 60-Ni. simulation box over ∼ 500 times, allowing the system to equilibrate. The simulation was next cooled down to the original temperature over 2 × 106 steps, resulting in a re-crystallization of the solid. It was important to know the structure of the final crystal for future simulations. For example if the sample was BCC or FCC. Also, the pair correlation functions of each type of atom would indicate if there was chemical separation occurring in the simulations. The melting temperatures of the crusts were calculated in the same man- ner as described in Chapter 5, by defining the melting temperature as the temperature corresponding to the largest jump in the mean squared dis- placement. Following this method, the melting temperature of the modified Urca crust was found to be T ∗ = 0.0103312. The thick crust had a melting temperature of T ∗ = 0.010131087. In the case of the thin crust the tran- sition from solid to liquid occurred at a temperature of T ∗ = 0.0104601. The three different crusts were observed to melt near the same temperature with the thin crust melting at the highest temperature. The melting tem- peratures for the impure crusts are also similar to that of the pure BCC 110 7.3. Melting Simulations type 1 2 3 1 1.0 1.0769 0.9231 2 – 1.1598 0.9941 3 – – 0.8521 Table 7.4: Pair coefficients, (Ze)2, of the Yukawa potential used in the simulations for the thick crust compositions. Type 1 atoms are 56-Fe, type 2 atoms are 60-Ni and type 3 atoms are 52-Cr. type 1 2 3 1 1.0 1.0769 1.0 2 – 1.1598 1.0769 3 – – 1.0 Table 7.5: Pair coefficients, (Ze)2, of the Yukawa potential used in the simulations for the thin crust composition. Type 1 atoms are 54-Fe, type 2 atoms are 58-Ni and type 3 atoms are 56-Fe. crystal of Chapter 6. At a density of ρ = 1011 g/cm3 these melting tem- peratures correspond to ∼ 7 × 108 K. In the simulations where the crystal was deformed the temperature of the simulations was always less than the melting temperature. The final structure and the possibility of chemical separation were quan- tified via pair correlations between the particles of the sample. The radial distribution function (RDF) of all the particles was calculated for the re- crystallized crust at the end of the simulation. The RDF of all the particles at the end of the simulation was compared to the RDF of a BCC crystal in all three cases. The RDF was also calculated for each isotope type, in order to look for signatures of chemical separation. The RDF of all isotope types in the modified Urca crustal composition was compared to that of a BCC lattice as shown in Figure 7.1. The am- plitude of the peaks, g(r), was scaled down by a factor of 5 in order to visually compare the spacings of the two structures. The spacings of g(r) for the modified Urca crustal composition correspond to those of the BCC lattice structure. The individual pair correlations for each type of isotope 111 7.4. Stress-Strain Relationships are compared in Figure 7.2. Evidence of chemical phase separation would be apparent if different types of isotopes had different RDFs, such as one type of isotope being in the solid state and another in the liquid state. The RDFs of each of the isotopes did not exhibit a signature of chemical separation. In the case of the thick crust composition, the final structure from the cooling of the crustal melt was compared to a BCC lattice as presented in Figure 7.3. As in the modified Urca case the amplitude of the peaks are scaled down. Comparing the spacings of the peaks of the thick crust to that of a BCC lattice indicates the thick crust lattice corresponds to a BCC lattice. Figure 7.4 compares the RDFs of individual isotope types in order to investigate chemical separation, but there are no signatures evident. For the thin crust composition the RDF of all the atom types was com- pared to a BCC structure and is shown in Figure 7.5, where again the am- plitude of the BCC RDF is scaled down for visual comparison. The RDFs of each of the atom types in the thin crust are compared in Figure 7.6. The structures of each of the different types of isotopes did not indicate chem- ical separation. The lack of evidence for chemical phase separation in all three cases could be due to the final temperature of the simulations, perhaps the final temperature was low enough that all atom types crystallized and thus different phases were not evident. At hotter temperatures chemical separation may become evident. 7.4 Stress-Strain Relationships Following the procedure outlined in Chapter 6 the mechanical properties, such as breaking strain and shear modulus, were calculated for the three dif- ferent crustal compositions. For the deformation simulations, perfect BCC lattices were used with impurities randomly assigned to lattice locations. The results of the impure BCC crystals were also compared to the pure BCC crystal case. Strain rates of 20× 10−6, 5× 10−6, and 2.5× 10−6 were used in these simulations with the lattice deformed to a 0.2 strain. As per- fect BCC lattices were used in these simulations and it is unlikely that a neutron star would have a perfect BCC lattice crust, the results of the shear 112 7.4. Stress-Strain Relationships Mu3 End: All types radius g( r) 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 1 2 3 4 5 6 Mu3 data BCC Figure 7.1: The RDF of all the atomic types in the recrystallized modified Urca composition compared to that of a BCC structure. The BCC structure is scaled down in amplitude in order to compare the spacings of the peaks. The spacings of the modified Urca agree with those of the BCC lattice. 113 7.4. Stress-Strain Relationships g( r) Mu3 End: Different Types 0 2 4 6 0 2 4 6 radius 0 1 2 3 4 5 6 2 4 6 56Fe 54Fe 60Ni Figure 7.2: The RDFs of individual atom types within the recrystallized modified Urca composition. Differences in the RDFs, such as one indicating a liquid structure, would be an indication of chemical phase separation. There is no evidence of chemical phase separation in this crustal composition. 114 7.4. Stress-Strain Relationships Thick End: All types radius g( r) 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 1 2 3 4 5 6 7 Thick data BCC Figure 7.3: The RDF of all atomic types in the recrystallized thick crust composition compared to that of a BCC structure. The BCC structure is scaled down in amplitude in order to compare the spacings. The spacings of the thick crust RDF correspond to those of the BCC lattice. 115 7.4. Stress-Strain Relationships g( r) Thick End: Different Types 0 2 4 6 8 0 2 4 6 8 radius 0 1 2 3 4 5 6 2 4 6 8 56Fe 60Ni 52Cr Figure 7.4: The RDFs of individual atoms types within the recrystallized thick crust composition. Differences in the peak spacings of the RDFs would be an indication of chemical separation. There is no evidence of chemical separation in this composition. 116 7.4. Stress-Strain Relationships Thin End: All types radius g( r) 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 1 2 3 4 5 6 Thin data BCC Figure 7.5: The RDF of all atomic types in the recrystallized thin crust composition compared to that of a BCC structure. The BCC structure is scaled down in amplitude in order to compare the spacings. The spacings of the thin crust composition correspond to those of a BCC lattice. 117 7.4. Stress-Strain Relationships g( r) Thin End: Different Types 0 2 4 6 0 2 4 6 radius 0 1 2 3 4 5 6 2 4 6 54Fe 58Ni 56Fe Figure 7.6: The RDFs of individual atom types in the recrystallized thin crust composition. Differences in the RDFs between the types of atoms would be an indication of chemical phase separation. There is no evidence of chemical phase separation for the thin crust composition. 118 7.4. Stress-Strain Relationships modulus and the breaking strain set upper limits, as adding in defects would likely work to weaken the structure, decreasing the breaking strain. The calculations of the breaking strain and the shear modulus for the three strain rates are compiled in Table 7.6 for the modified Urca compo- sition, Table 7.7 for the thick crust composition and Table 7.8 for the thin crust composition. These results were also compared to those of the pure BCC lattice from Chapter 6. An example of a stress-strain relationship of the three impure BCC lattices compared to the pure BCC lattice at a strain rate of ṡ = 20 × 10−6 is displayed in Figure 7.7. The modified Urca and thick crust display similar properties to the pure BCC lattice. The thin crust is found to yield to the stress at a larger degree of deformation than the other compostions. For all three strain rates, ṡ = 20 × 10−6, 5 × 10−6, and 2.5× 10−6, the modified Urca and thick crust compositions shared sim- ilar shear moduli and breaking strains to the pure BCC lattice; the thin crust composition consistently has a higher breaking strain. All three com- positions have shear moduli of ∼ 0.27. The breaking strain for the modified Urca and thick crust occurred at φm∼ 0.11, and at φm∼ 0.12 the thin crust yielded to the applied stress. Strain Rate ṡ/ωp Breaking Shear (×10−6) (×10−7) Strain Modulus 20 56.42 0.115620 0.2735 5 3.979 0.113396 0.2736 2.5 1.989 0.113174 0.2735 Table 7.6: The breaking strain and shear modulus for a perfect BCC lattice with a composition of the modified Urca crust. As mentioned earlier the calculated shear modulus and breaking strain for the above perfect crystals set an upper limit for the neutron star crust as the neutron star crust is unlikely to be a perfect crystal. Applying a strain to the recrystallized structure found in Section 7.3 places a limit on the stress-strain relationship for an imperfect crystal. The imperfect crystal was created by using the results from the melting simulation runs. For the 119 7.4. Stress-Strain Relationships BCC Rate 20 Strain σ x y [P re ss u re ] 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.200 .0 00 0. 00 5 0. 01 0 0. 01 5 0. 02 0 0. 02 5 0. 03 0 mu3 thick thin bccfe Figure 7.7: The stress-strain relationship for the three crustal compositions: modified Urca (mu3), thick crust (thick), and thin crust (thin), compared to a pure BCC lattice (bccfe). A strain rate of 20 × 10−6 was used in all four simulations. The shear modulus was found to be similar for all the simulations. The thin crust yielded to the deformation at a higher strain than the other three simulations. As these results use a perfect BCC lattice structure, these simulations set an upper limit for the shear modulus and breaking strain of a neutron star crust. 120 7.4. Stress-Strain Relationships Strain Rate ṡ/ωp Breaking Shear (×10−6) (×10−7) Strain Modulus 20 56.42 0.113285 0.27179 5 3.979 0.111895 0.27178 2.5 1.989 0.111673 0.27179 Table 7.7: The breaking strain and shear modulus for a perfect BCC lattice structure with a thick crust composition. Strain Rate ṡ/ωp Breaking Shear (×10−6) (×10−7) Strain Modulus 20 56.42 0.122846 0.26644 5 3.979 0.120122 0.27620 2.5 1.989 0.119455 0.27711 Table 7.8: The breaking strain and shear modulus for a perfect BCC lattice containing a thin crust composition. melting simulations a perfect BCC crystal with impurities is melted and recrystallized. Though the recrystallized material has a global structure of BCC, in the cooling process defects, such as grains are introduced, creating an imperfect crystal. The stress-strain relationship for the thick crust in the perfect and imperfect crystal cases are compared in Figure 7.8. The two cases show large differences with regards to the shapes of the curves. In Figure 7.9 the imperfect crystal is shown in more detail. Unlike perfect crystals, imperfect crystals display plastic activity occurring at small strains. The properties such as shear modulus, breaking strain, and breaking stress are different between the perfect crystal and the imperfect crystal containing defects. The figure is for a thick crust, but all three impure crystals displayed this same type of behaviour for the imperfect crystals. This result is different than that observed in Horowitz & Kadau (2009). 121 7.4. Stress-Strain Relationships Thick Strain σ x y [P re ss u re ] 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20 0. 00 0 0. 00 5 0. 01 0 0. 01 5 0. 02 0 0. 02 5 Perfect Imperfect Figure 7.8: Comparing the stress-strain relationships of a perfect and an imperfect crystal containing defects for the thick crust composition. A strain rate of ṡ = 20×10−6 was applied in both simulations. The two crystals show different behaviours in their stress-strain relationship. 122 7.4. Stress-Strain Relationships Thick Strain σ x y [P re ss u re ] 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20 0. 00 0 0. 00 1 0. 00 2 0. 00 3 0. 00 4 0. 00 5 0. 00 6 0. 00 7 0. 00 8 0. 00 9 0. 01 0 Figure 7.9: The imperfect crystal, which contains defects, for the thick crust composition with a strain rate of ṡ = 20 × 10−6 applied. At small strain values small yielding events can be observed. 123 7.5. Second Break 7.5 Second Break Simulations were performed in order to investigate the possibility of a second yielding event occurring after the first in the case of perfect BCC crystals. These simulations were run to a strain of ∆l/l = 0.4, which is twice as large as the original simulations. A strain rate of ṡ = 20 × 10−6 was applied in all of the cases. The stress-strain relationship of the three different crustal compositions is compared in Figure 7.10. In all three cases the shear modulus of the second fracturing event decreases as compared to the first event. For the modified Urca composition the shear modulus decreased from µ∼ 0.27 to µ∼ 0.21 for the second break, measured between 0.13 to 0.17 strain. In the case of the thick crust composition the shear modulus decreased from µ∼ 0.27 before the first break to µ∼ 0.12, measured between a strain of 0.14 to 0.2, for the second break. Finally, for the thin crust composition the shear modulus was µ∼ 0.28 prior to the first break and µ∼ 0.05 prior to the second break, determined between 0.13 to 0.16 strain. 7.6 Forward and Reverse Simulations where the crystal is deformed to the point where the lattice yields, and then the deformation is reversed were performed in order to investigate the reversibility of the yield event. These simulations were per- formed on all three impure compositions with perfect BCC lattices at a strain rate of ṡ = 20 × 10−6. When the breaking point was reached the same strain rate was applied but in the reverse direction. The stress-strain relationships for the three different crustal compositions are compared in Figure 7.11. For plotting purposes, the strain measurements on the x-axis of the figure continue to increase after reversing the shear direction as (step size×∆t× 20× 10−6). For both the modified Urca and thick crust composi- tions a strain of 0.12 was reached before reversing the direction. For the thin crust composition a strain of 0.1265 was reached before reversing direction. In all three cases the shear modulus decreased after reversing the direction of shear. For the modified Urca composition the shear modulus changed 124 7.6. Forward and Reverse BCC Second Break Strain σ x y [P re ss u re ] 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0. 00 0 0. 00 5 0. 01 0 0. 01 5 0. 02 0 0. 02 5 0. 03 0 mu3 thick thin Figure 7.10: The stress-strain relationship of the three impure compositions. The simulations were deformed to a total strain of 0.4 at a strain rate of 20× 10−6 in order to search for a second yielding event. The shear modulus was found to decrease from the first break to the second. In the case of the thin crust the shear modulus dramatically decreased. 125 7.7. Temperature Dependence from µ∼ 0.27 in the forward direction to µ∼ 0.24, measured between a total strain of 0.19 and 0.21, in the reverse direction. The shear modulus was found to be µ∼ 0.27 in the forward direction and µ∼ 0.23, between a total strain of 0.2 and 0.24, in the reverse direction for the thick crust composi- tion. In the thin crust composition the shear modulus was µ∼ 0.28 in the forward direction and µ∼ 0.2 in the reverse direction, measured between 0.13 to 0.20 total strain. 7.7 Temperature Dependence The temperature dependence of the breaking strain and shear modulus was investigated in the perfect BCC crystals with crustal compositions, as was done for the pure crystal cases discussed in Chapter 6. The original sim- ulations were performed at a temperature of T ∗ = 0.001 and additional simulations were performed at a hotter temperature of T ∗ = 0.002 and a colder temperature of T ∗ = 0.0005. The melting temperature of the three different compositions was found to be T ∗∼ 0.01, thus all the additional temperatures are below the melting temperature of the crystal. A strain rate of ṡ = 20 × 10−6 was applied in these simulations. Each of the three crustal compositions were simulated in this investigation of the temperature dependence of the stress-strain relationship. The stress-strain relationships of the three different temperatures are compared for the modified Urca com- position in Figure 7.12, for the thick crust composition in Figure 7.13 and finally for the thin crust in Figure 7.14. For all three different compositions the colder simulation, T ∗ = 0.0005, had a larger shear modulus and break- ing stress. For the modified Urca and the thick crust simulations the colder simulations have breaking strains smaller compared to the hotter simula- tions, which is the same behaviour observed for the pure BCC case. For the thin crust simulations, the shear modulus was observed to differ for the three simulation temperatures, however the breaking strain was found to be very similar. Table 7.9 compares the breaking strain, shear modulus and the peak energy for the different temperatures and compositions. 126 7.7. Temperature Dependence BCC Forward and Reverse Strain σ x y [P re ss u re ] 0.00 0.05 0.10 0.15 0.20 0.25−0 .0 15 −0 .0 10 −0 .0 05 0. 00 0 0. 00 5 0. 01 0 0. 01 5 0. 02 0 0. 02 5 0. 03 0 mu3 thick thin Figure 7.11: The stress-strain relationship at a strain rate of 20 × 10−6 for three different crustal compositions. In these simulations the box is deformed to the breaking point and then the direction of deformation is reversed, dotted lines indicate the reverse direction. The shear modulus in the reverse direction is less than in the forward direction. 127 7.7. Temperature Dependence MU3: Temperature Strain σ x y [P re ss u re ] 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20 0. 00 5 0. 01 0 0. 01 5 0. 02 0 0. 02 5 0. 03 0 T = 0.002 T = 0.001 T = 0.0005 Figure 7.12: The stress-strain relationship at a strain rate of 20 × 10−6 for different temperatures of a perfect BCC lattice with a modified Urca composition. The colder temperature run has the largest shear modulus and breaking stress. 128 7.7. Temperature Dependence Thick: Temperature Strain σ x y [P re ss u re ] 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.200 .0 00 0. 00 5 0. 01 0 0. 01 5 0. 02 0 0. 02 5 T = 0.002 T = 0.001 T = 0.0005 Figure 7.13: The stress-strain relationships for a BCC lattice with a thick crust composition at three different temperatures. The simulations all had a strain rate of 20× 10−6 applied. The melting temperature of a thick crust is approximately 0.01, above the temperatures of the runs. The coldest simulation is found to have a steeper slope, or larger shear modulus than the other two temperatures. 129 7.7. Temperature Dependence Thin: Temperature Strain σ x y [P re ss u re ] 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20 0. 00 5 0. 01 0 0. 01 5 0. 02 0 0. 02 5 0. 03 0 T = 0.002 T = 0.001 T = 0.0005 Figure 7.14: The stress-strain at three temperatures for a BCC lattice con- taining the composition of a thin crust. A strain rate of 20×10−6 was applied in all three cases. All three simulations have similar breaking strains, though the coldest simulation has the largest shear modulus and breaking stress. 130 7.8. Discussion and Conclusions Structure Temperature Breaking Shear T ∗ Strain Modulus Modified 0.002 0.116064 0.26253 Urca 0.001 0.115620 0.27351 0.0005 0.112952 0.29558 Thick 0.002 0.115397 0.26052 Crust 0.001 0.113285 0.27179 0.0005 0.112396 0.27789 Thin 0.002 0.122846 0.26644 Crust 0.001 0.121901 0.27637 0.0005 0.122401 0.28060 Table 7.9: The temperature, breaking strain and shear modulus for the three different crustal compositions. The three different crust compositions melt at approximately 0.01, all temperature simulations are below the melting temperature. At a density of ρ = 1014 g/cm3 the corresponding physical temperaturs are 6.9×108 K for T ∗ = 0.002, 3.5×108 K for T ∗ = 0.001 and 1.7× 108 K for T ∗ = 0.0005. 7.8 Discussion and Conclusions Using the compositions calculated for a non-accreting neutron star, the me- chanical properties of these crusts have been investigated. The compositions, as discussed in Chapter 2, are dependent on how the neutron star cools. A neutron star cooled via the modified Urca process, a neutron star with a thick crust and one with a thin crust were considered in this study. The melting temperature and structure of recrystallization have been investigated in all three cooling cases. The structure and chemical phase separation of the three cases were investigated through the RDFs of all the particles, as well as each type of particle. With using the RDFs, chemical phase separation would be evident if RDFs suggesting a solid state structure and a liquid state were present for different isotopes. Horowitz et al. (2007) found for an ac- creted crust composition lighter Z isotopes to be preferentially in the liquid phase. None of the simulations of the non-accreted crusts showed any evi- dence for chemical separation. This lack of evidence of chemical separation could be due to the temperature at which this calculation was computed; at 131 7.8. Discussion and Conclusions a higher temperature chemical separation may occur. The isotopes in the non-accreted crust also have charges which are very similar and this could have also had an impact on the occurrence of chemical phase separation. The calculation of the static structure factor (see Horowitz & Berry, 2009) is left for future work. The three compositions were used in simulations to determine the shear modulus and breaking strain of the crustal type and compared that of a pure iron BCC crystal. The modified Urca and the thick crust were found to have shear moduli very similar to the pure iron BCC crystal. A change of compositions in Horowitz & Kadau (2009) was also found to have little effect on the stress-strain relationship. The simulations used perfect BCC lattices with the impurities randomly placed at lattice locations. With the perfect crystal structure upper limits can be placed on the neutron star crust, as defects, if anything, would be expected to decrease the rigidity of the crust. The thin crust has a breaking strain which is larger then the two other impure simulations and the pure BCC simulation. Imperfect crystals with the same compostions compared to perfect crystals indicate the two systems behave differently with the same applied stress. This result is different from the results observed in Horowitz & Kadau (2009), where there were 6 grains introduced to the simulation box and only moderately effected the stress-strain relationship. The effect of the number of grains and the resulting stress-strain relationship is something to be investigated in futre work. Note that the results presented here for the non-accreted crust compositions have a different applied strain rate than those in Horowitz & Kadau (2009). Simulations were also performed to have the perfect crystal deform past the first point of breaking in order to search for a second breaking event. The shear modulus was found to be less for the second breaking event as compared to the first time the crystal yielded to the deformation. Simula- tions which deformed the crystal to the point of breaking and then had the deformation reversed indicate that the shear modulus in the reverse direc- tion is less compared to the shear modulus in the forward direction, thus the yielding event is not reversible. The temperature dependence on the 132 7.8. Discussion and Conclusions shear modulus and the breaking strain was also investigated for all three crustal compositions of the non-accreting crusts. Hotter temperatures were found to have the smallest shear modulus for all three crustal compositions. The failure of the crystals were found to occur at smaller strains at colder temperatures. The breaking stress was found to be greatest for the colder temperatures. As in the pure crystal simulations the change in breaking strains with temperature is small, suggesting that SGR burst activity would also not change much as the neutron star cools. Simulations of the imperfect crystals at various temperatures would help to constrain the expectations for SGR burst activity. The results of the upper limits placed on the shear modulus and breaking strain for the different compositions are applied to neutron star observations in the next chapter. 133 Chapter 8 Application of Breaking Strain 8.1 Introduction The observational consequences of crustal cracking are investigated in this chapter using the calculated shear modulus and breaking strain. The shear modulus and breaking strain from the perfect crystal simulations for the three different crustal compositions and pure iron crust are considered in context of crustal cracking. The amplitude of glitches, elastic energy of Soft Gamma-ray Repeaters, and gravitational wave detection are considered in terms of the results of the simulations. The modified Urca and thick crust compositions have shear moduli and breaking strains close to that of a pure iron BCC lattice, so it is expected that these three cases will have similar observational consequences. The shear modulus calculated from the simulations is measured in re- duced pressure units and is scalable with density. From Table 3.1 the reduced pressure is expressed as: P ∗ = P a3 Uo , (8.1) where P is the physical pressure, Uo is the characteristic energy and a = n −1/3 is the characteristic distance. By substituting the values of Uo and a into the above equation for the system of interest, the physical pressure can be determined. The value of the shear modulus was found to be the same for the three impure and the pure iron crusts, 0.27 in reduced units. The shear modulus in physical units is given by the expression µ = 0.27Uo/a 3 = 0.27 (Ze)2/a4, using the unit conversions discussed in Section 3.3. 134 8.2. Glitches The observational consequences of crustal cracking are typically investi- gated at the bottom of the neutron star crust, at a density of ρ = 1014 g/cm3. The results for the simulations of the crusts from Chapter 7 are appropriate for densities ranging from ∼106 g/cm3 to ∼1011 g/cm3, but in order to com- pare to previous work, the simulation results are extrapolated to the greater density. At the bottom of the crust the neutron fraction is higher than at lower densities which are closer to the surface. At the bottom of the crust the neutron fraction is Xn = 0.8, also the charge can be set to Z = 20 and the atomic mass can be set to A = 88 (Ushomirsky et al., 2000). With the above considerations the number density is related to the neutron fraction as: n = a−3 = ρ m (1−Xn) = 10 14 g/cm3 88mu × 0.2, (8.2) where mu is the atomic mass unit. At the density at the bottom of the crust, ρ = 1014 g/cm3, the characteristic length is a = n−1/3 = 1.94 × 10−12 cm or 19.4 fm. The characteristic energy at the bottom of the neutron star crust is Uo = (Ze) 2/a = 4.76× 10−5 erg, which is ∼30 MeV, for a charge of Z = 20. As a result the shear modulus at the bottom of the neutron star crust is µ = 6.5× 1030 dyne/cm2, which is close to the previously used value of 1030 dyne/cm2 (for example Ruderman, 1969). Note that at a density of ρ = 1011 g/cm3 the shear modulus would be µ ∼ 1027 dyne/cm2. The maximal breaking strain is a unitless quantity, thus the simulation values are the physical values with a breaking strain of φm ∼ 0.11− 0.12 used in order to compare the simulation results to observations. The breaking strain found for the non-accreted neutron star crust is in agreement of the breaking strain calculated in Horowitz & Kadau (2009) of close to φm = 0.1, for an accreted crustal composition. 8.2 Glitches Using the shear modulus and breaking strain as calculated from the crustal simulations it is possible to now calculate the corresponding glitch ampli- tude, ∆Ω/Ω, where ∆Ω is the change in angular frequency. As discussed in 135 8.2. Glitches Chapter 1, following Ruderman (1969) the glitch amplitude is related to the shear modulus and breaking strain as ∆Ω Ω = 2 ( 95µφmR 7GMρ [ 1− ( ∆R R )7]) , (8.3) where µ is the shear modulus, φm is the breaking strain, and ∆R is the crust thickness. Following Ruderman (1969) a density of 1014 g/cm3 is also used to calculate the glitch amplitude. The shear modulus calculated from the simulations at a density of ρ = 1014 g/cm3 is µ = 6.5 × 1030 dyne/cm2. With the calculated breaking strain of φm ∼ 0.1 and for a neutron star with a crustal thickness of 1 km, a radius of 10 km and a mass of 1.4 M, the corresponding calculated glitch amplitude is ∆Ω/Ω ∼ 10−3 for glitches at the bottom of the crust. Glitch amplitudes have been observed to range primarily from 10−9 to 10−6 (Lyne et al., 2000), with PSR B2334+61 having the largest observed glitch amplitude of ∆Ω/Ω = 20.5× 10−6 (Yuan et al., 2010). The calcu- lated glitch amplitude of ∆Ω/Ω∼10−3, using the value of breaking strain and shear modulus for a crust with perfect crystal structure, is much larger than has been observed from neutron stars. Two possible explanations for this difference between the calculated and observed glitch amplitude are, first a question of the structure of the crust and second, the origin of glitches. The shear modulus and breaking strain used in the glitch amplitude calcu- lations are from a perfect crystal. These perfect crystal parameters place an upper limit on the glitch amplitude, a smaller breaking strain would result in smaller glitch amplitudes. The neutron star crust is unlikely to be a perfect crystal. In Section 7.4 an imperfect crystal neutron star crust demonstrated different behaviour in the stress strain curves than the perfect crystals. The breaking strain of the neutron star crust is likely dependent on the structure of the crust and how close the structure is to a perfect crystal. The differences found in the breaking strain and shear modulus in Section 7.4 for the perfect versus the imperfect would likely lead to a prediction of smaller glitch amplitudes for an imperfect crystal neutron star crust. In terms of glitch origins, glitches are also unlikely to involve starquakes alone, 136 8.3. Soft Gamma-ray Repeaters they might also require vortex unpinning in the neutron star superfluid and subsequent transfer of angular momentum (Lyne et al., 2000). 8.3 Soft Gamma-ray Repeaters A total of eleven Soft Gamma-ray Repeaters (SGRs), seven confirmed and four candidates, have been observed7. The characteristic bursts from these SGRs have reached peak luminosities of 1041 erg/s (Woods & Thompson, 2006). Three of the total eleven SGRs have been observed to also exhibit gi- ant flare events, which have been observed to reach a luminosity of 1044 erg/s (Helfand & Long, 1979). The bursts can be interpreted as being triggered by crustal cracking due to stresses placed on the crust from the magnetic field (Chamel & Haensel, 2008). The strength of the magnetic field required to crack the neutron star crust and release the observed burst energy, is dependent on the breaking strain of the material: B = 1015 ( ESGR 1041erg )−1/2 ( l km )( φm 10−3 ) G, (8.4) where ESGR is the energy of the burst, l is the length of the crustal fracture and φm is the breaking strain (Thompson & Duncan, 1995). A limit on the size of the crust fracture can be estimated using the observed energy emitted during the burst, the neutron star’s magnetic field strength, and the breaking strain calculated from the simulations. The estimated size of the crustal fracture can be used as an indication of if the bursting event can be attributed to crustal events alone or if the burst is related to another region of the neutron star. If the fracture length is too large, the burst can not be attributed to the crust alone. Using Equation 8.4, the length of the fracture can be expressed as: l = ( B 1015 G )( ESGR 1041 erg )( 10−3 φm ) km. (8.5) 7Data from the McGill SGR/AXP online catalog, http://www.physics.mcgill.ca/ pul- sar/magnetar/main.html 137 8.4. Gravitational Waves In the case of typical SGR bursts, with burst energies of ESGR = 10 41 erg, assuming a magnetic field of B = 1015 G and using the perfect crystal break- ing strain value, φm ∼ 0.1, the size of the fracture is l = 10 m. But, in the case of giant flares, the energy associated is 1000 times larger than the typ- ical bursts with ESGR = 10 44 erg. Using the giant flare energetics and the same values for the magnetic field and breaking strain as for the typical bursts, the fracture length is l = 300 m. The crust of the neutron star extends to a depth of approximately 1 km, with the outer crust itself extending ∼100 m. Whether the SGR bursts can be attributed to crustal cracking depends on if the fracture length associated with observed energy is within the bounds of the crust. With a fracture length of l = 10 m, for the lower energy SGR burst events, the typical SGR bursts could be due to crustal cracking, for crusts which have a breaking strain of φm = 0.1. The giant flare events require a much larger fracture length, l = 300 m. A fracture of this size would encompass ∼15 − 30% of the crust depth, so the fracture could propagate along the surface or vertically. The fracture length required to reach the giant flare energetics is within the size constraints of the neutron star crust, and thus the giant flares could be attributed to crustal cracking. These crustal fracture lengths were calculated using the breaking strain found for perfect crystal configurations. The perfect and imperfect crystals were found to display different mechanical properties. For a neutron star crust with a smaller breaking strain, the fracture length required for a burst would increase. Depending on how close the crust is to a perfect or imperfect structure would dictate either how the fracture propagates, along the surface or vertically with depth, or if crustal cracking could be ruled out as a mechanism for SGR bursts. 8.4 Gravitational Waves A neutron star perturbed to be asymmetric about its rotation axis may have a varying quadrupole moment and may be a source of gravitational waves. A neutron star with a varying quadrupole moment could result from either the star undergoing precession or from a mountain on the star surface. The 138 8.4. Gravitational Waves free precession of a neutron star is a possible result from an asymmetrically shaped star (Stairs et al., 2000). From these gravitational wave emitting sources the observed strain amplitude of the wave at Earth is expressed in terms of a sinusoidally varying quadrupole: h = 16 5 ( pi 3 )1/2 GQ22 Ω2 d c4 , (8.6) where Q22 is amplitude of the varying quadrupole moment, Ω is the angular frequency and d is the distance to the neutron star (Ushomirsky et al., 2000). With respect to the study of the mechanical properties of the crust, it is the gravitational waves emitted due to mountains on the neutron star surface which are of interest. The strength, or breaking strain, of the crust dictates the mountain size and thus, the maximum quadrupole which a neutron star can support. The maximum quadruple is found to depend directly on the breaking strain of the crust. The calculation of the maximum quadrupole is covered in Ushomirsky et al. (2000). Following this calculation, the quadrupole can be related to both the breaking strain and shear modulus. In Ushomirsky et al. (2000) the maximum quadrupole is expressed as Q22 = γ φm I, where γ is a factor which depends on the mass and radius of the neutron star, φm is the breaking strain, and I is an integral which is directly proportional to the shear modulus, as a result Q22 is also directly proportional to the shear modulus, µ. For a shear modulus of µ = 1030 dyne/cm2, the quadrupole is expressed as: Q22 = 10 38 g cm2(φm/10 −2). Introducing the directly propor- tional dependence of Q22 on the shear modulus results in the relationship: Q22 = 10 38 g cm2 ( φm 10−2 )( µ 1030dyne/cm2 ) . (8.7) For the simulation results of the perfect crystal structure with a breaking strain of φm = 0.1 and at a density of ρ = 10 14 g/cm3, where µ = 6.5 × 1030 dyne/cm2, the quadrupole is Q22 = 6.5× 1039g cm2, assuming the crust is maximally deformed. With the calculated quadrupole, it is possible to compare calculated 139 8.4. Gravitational Waves limits of gravitational wave strain amplitudes to the observed limits. Ob- servational upper limits using LIGO have been placed on the gravitational wave strain amplitudes of 78 pulsars in Abbott et al. (2007). The observa- tional limits of three of the 78 pulsars are compared to the results of the simulations. The three pulsars examined are: PSR J1603-7202, PSR J2124- 3358, and PSR J0534+2200, the Crab pulsar. Two of the three pulsars are recycled and as such have undergone periods of accretion, as well, the Crab pulsar is in a supernova remanent and may have accreted some mat- ter. Even though the three pulsars investigated are not necessarily pristine, the breaking strain results found for the non-accreting neutron star crust composition do agree with that of the accreted crust reported in Horowitz & Kadau (2009). The strain amplitudes are calculated using Equation 8.6 and the quadrupole calculated for a maximally deformed, perfect crystal crustal structure. The LIGO observational limits of the gravitational wave strain amplitude, h, are compared to the strain amplitudes predicted for each pulsar based on the simulation results in Table 8.1. The distance, d, and the spin frequency, ν, are also listed in the table. Pulsar d ν LIGO Predicted kpc Hz Maximum log(h) Maximum log(h) PSR J1603-7202 1.64 67.38 -24.58 -25.2 PSR J2124-3358 0.32 202.71 -24.31 -23.5 PSR J0534+2200 2 29.80 -23.51 -25.0 Table 8.1: The distance, d, spin frequency, ν, the observational and predicted gravitational wave strain amplitude, h, for the three pulsars of interest. The predicted gravitational wave strain amplitudes are below the observational LIGO upper limits for two of the cases. For pulsar PSR J2124-3358 the predicted strain amplitudes indicate that if the crust was perfectly BCC in structure and maximally deformed, than gravitational waves should have been detected from this source. The distances of each pulsar were taken from the ATNF Pulsar Catalogue (www.atnf.csiro.au/people/pulsar/psrcat (Manchester et al., 2005). The spin frequency and LIGO upper limits are from Abbott et al. (2007). The predicted maximum strain amplitudes are calculated using the distance and spin frequency of each of the pulsars, as well as the quadrupole as calculated using the simulation results. 140 8.5. Conclusions The comparison of the LIGO observational limits and the predicted strain amplitudes in Table 8.1 indicate which neutron stars would be ex- pected to have gravitational waves detected, if the crust is perfectly BCC and maximally deformed. Of the three pulsars examined, PSR J2124-3358 has a predicted strain amplitude greater than the observational limit. This means if PSR J2124-3358 had a perfect BCC lattice and the crust was maximally deformed, or deformed to even ten percent of the maximum, gravitational waves would have been detected from this source. With the reasonable assumption that the crustal strain is uniformly distributed up to the breaking point, the non-detection of the gravitational waves suggests that the neutron star crust is not a perfect BCC crystal. The non-detection of gravitational waves could also result from the crust not being maximally deformed. The search for gravitational waves from pulsars could be used as a method to determine the structure, strength, or degree of deformation of the neutron star crust. 8.5 Conclusions The observational consequences of crustal cracking were investigated us- ing the shear modulus and breaking strain calculated using the molecular dynamics simulations. The observation consequences of crustal cracking considered included glitches, SGR bursts, and gravitational waves. The shear modulus from the simulations scales with density and at a density of ρ = 1014 g/cm2, the bottom of the crust, the shear modulus is µ = 6.5× 1030 dyne/cm2, for a perfect BCC crystal configuration. The breaking strain used for comparison to observation was φm ∼ 0.1, for a perfect BCC crystal for the three different crustal compostions as well as the pure iron crust. These two parameters were used in all the comparisons to observa- tions. The glitch amplitude for a perfect crystal at a density of ρ = 1014 g/cm2 was found to be ∆Ω/Ω ∼ 10−3, which is much larger than the largest ob- served glitch. The differences between the calculated and observed glitches could result from the actual structure of the crust, as well as the glitch ori- 141 8.5. Conclusions gins. The crust may not be a perfect BCC crystal as thus have a different breaking strain. Also, the glitches may not be due to crustal cracking alone, but may also require vortex unpinning. In terms of SGR bursts, the frac- ture length required in order to release the observed energy was calculated. Crustal cracking could be a possible explanation for the SGR bursts if the fracture length is within the constraints of the crust. For the typical SGR burst of ESGR = 10 41 erg the fracture length was found to be 10 m, small enough to occur in the crust depth. The energy associated with giant flares require fracture lengths of ∼300 m, which would encompass a large fraction of the crust. For the gravitational waves, the strain amplitude for a max- imally deformed, perfect crystal crust, was compared to the observational limits placed by LIGO. The maximal gravitational wave strain amplitude predicted for PSR J2124-3358 was found to be ten times larger than the upper limits determined from LIGO. If PSR J2124-3358 had perfect BCC crystal lattice crust and was deformed to ten percent of its maximum, than gravitational waves would have been detected from this pulsar. The non- detection of gravitational waves suggests that a neutron star crust is not a perfect BCC crystal lattice or that the crust is not deformed to even ten percent of the maximum. 142 Chapter 9 Conclusions 9.1 Summary The outer layers of the neutron star: atmosphere, envelope, and crust all play an important role in the observations of neutron stars. The atmosphere shapes the emerging thermal spectrum, while the envelope regulates the thermal transport of energy (Lattimer & Prakash, 2004). The strength of the crust and its ability to fracture has been tied to observations of glitches, SGR bursts and the possibility of eventually observing gravitational wave emission from a neutron star. In order to understand these upper regions of the neutron star, the compositions of the regions were first calculated for three different methods of cooling. With the calculated crustal compositions upper limits on mechanical properties such as the shear modulus and breaking strain of the neutron star crust were also calculated. Several observational consequences of these upper limits were also calculated: predicted glitch amplitudes, energy from SGR bursts, and predicted strain amplitudes of gravitational waves. The composition of the atmosphere and the crust was calculated using the nuclear reaction network torch. For the composition calculations a non- accreting neutron was considered, as a result the atmosphere of such a star would be dependent on isotopes produced due to nuclear reactions in the crust floating to the surface of the neutron star before the crust freezes. Three methods of cooling were considered in these calculations: a neutron star cooling via the modified Urca process, a neutron star with a thick crust, and one with a thin crust. Each of the different methods of cooling resulted in a different atmosphere composition. For a neutron star which cooled by the modified Urca process the atmosphere would be 28Si. A neutron star 143 9.1. Summary with a thick crust would have a 50Cr atmosphere. Finally a thin crust would be expected to have a 40Ca atmosphere. A non-accreting neutron star would be expected to have an iron crust. The composition calculations indicated that this was not the case and the three crusts contained various impurities. The isotopes which were calcu- lated to compose the crust were used in molecular dynamics simulations in order to determine the shear modulus and breaking strain of the systems. In each case a perfect BCC lattice structure was initially used. With the per- fect crystal structure an upper limit could be placed on the shear modulus and breaking strain. The three impure cases were also compared to a pure perfect BCC lattice. All three impure crusts share a simular shear modulus to the pure lattice, µ = 6.5× 1030 dyne/cm2 at a density of ρ = 1014 g/cm3. The breaking strains were also similar to the pure system, with the modi- fied Urca and thick crust having a breaking strain of φm∼ 0.11, which is the same as the pure system. The thin crust composition was found to have a breaking strain of φm∼ 0.12. These perfect crystals set upper limits to the properties of the neutron star crust. The behaviour of the imperfect versus the perfect crystals is quite different. The imperfect crust simulation boxes were created from re-crystallized melted perfect BCC crystals. The imperfect crystals display multiple plastic events at small amounts of strain before a larger yielding event. The perfect crystals do not show the small events before the crystal fails. More simulations will need to be performed in order to fully understand the properties of an imperfect crust. The observational consequences of the calculated upper limits were con- sidered in terms of glitches, SGR bursts, and gravitational waves. The glitch amplitudes are predicted to be on the order of ∆Ω/Ω∼ 10−3, which is much larger than the largest observed glitch (see Yuan et al., 2010). For the SGRs, the size of the fracture required for the energy release in the burst to be as- sociated with crustal cracking was calculated. For the typical SGR bursts of E = 1041 erg the fracture would be on the order of ∼ 10 m. In the case of giant flares with E = 1044 erg the required fracture size is ∼ 300 m. In terms of gravitational wave emission the gravitational wave strain amplitude was 144 9.2. Future Work calculated for three different pulsars. For PSR J2124-3358 the predicted maximal strain amplitude is ten times greater than the upper limit from LIGO (see Abbott et al., 2007). This indicates that PSR J2124-3358 does not have a crust which is a perfect BCC lattice that has been maximally deformed, or deformed to 10% of the maximum. 9.2 Future Work The work in this thesis can be used as a starting point for several different extensions. The atmosphere has been calculated for three different neutron star cooling methods, but the thermal evolution of the atmosphere can also be examined. With the molecular dynamics tool kit developed in this the- sis, the role the neutron star crust plays in giant flares from SGRs can be explored further. Other extreme matter besides neutron star crusts can be investigated using molecular dynamics, such as the matter in the interiors of planets. 9.2.1 Atmosphere Thermal Evolution With the calculation of the atmospheric composition, the three different cooling curves resulted in different atmospheres. These calculations were based on the material being fully ionized. As the neutron star cools there is a possibility that the atmosphere composition could change. As the neutron star cools, the material may no longer be fully ionized, which results in a change of the settling timescales of the isotopes. As the settling timescales change, the atmospheric compositions could change over time without ac- cretion occurring. This introduces the possibility that self-consistent cooling models could be created, where the type of atmosphere which fits the ob- served thermal emission would indicate the cooling processes within the neutron star. The surface composition for the different neutron stars were found not to be iron, but to be isotopes with nuclear free energy. This re- active surface would also change the expected surface composition after an accretion event of matter onto the star (Chang et al., 2010). 145 9.2. Future Work 9.2.2 Bursts from Soft Gamma-ray Repeaters The study of the shear modus and breaking strain of neutron star crusts can be extended to understanding what happens to the neutron star before and after a burst occurs. The crust may become more brittle or stronger after a bursting event. Subsequent bursts would be affected by changes in the structure of the crust, which could be characterized via the static structure factor. The mechanical properties of the crust are thought to be associated with the observed quasi-periodic oscillations observed in the tails of giant flares from SGRs. In order to predict the frequencies at which the oscillations occur astroseismology models can be developed. These types of models would give further understanding of neutron star structure. 9.2.3 Planetary Interiors Neutron stars are not the only astrophysical objects which can be further understood with molecular dynamics, these types of simulations can also be performed in order to understand planetary interiors. Calculating the melting temperature of carbon and iron at extreme pressures is important for giant planets and Earth-massed planets. It has been suggested that Uranus and Neptune contain a layer of carbon (Ross, 1981). Molecular dynamics simulations have been performed to calculate the melting curve of carbon at extreme pressures and temperatures (Correa et al., 2006). These simulations contained 128 and 256 atoms and size effects were noted in the results of the simulations. Larger simulations could be performed in order to avoid these size effects. For exoplanets with Earth-like masses the question of whether they have a liquid or solid core has implications for the strength of a magnetic field. Calculating the melting curve for iron at extreme pressures would help to answer this question. Chemical separation may also occur in giant planets (Helled et al., 2011). Simulations of the giant planets would be used for understanding the distribution of material in the planet interior. Performing molecular dynamics simulations on planets would help to create a self-consistent model for the evolution of planets. In this thesis the outer layers of the neutron star have been studied. Due 146 9.2. Future Work to nuclear reactions in the crust, a non-accreting neutron star has been found to not have an iron atmosphere, but instead a reactive surface. The strength of the crust had also been explored with molecular dynamics simulations. Using perfect crystal lattices, upper limits have been placed on the shear modulus and breaking strain of the crust. With the tools developed in this thesis there are many different areas for future study, including the evolution of the neutron star atmosphere, as well as the extreme matter of neutron star crusts and also the interiors of both giant and rocky planets. 147 Bibliography Abbott, B., Abbott, R., Adhikari, R., Agresti, J., Ajith, P., Allen, B., Amin, R., Anderson, S. B., Anderson, W. G., Arain, M., & et al., 2007, Phys. Rev. D 76(4), 042001 Alcock, C. & Olinto, A., 1988, Annual Review of Nuclear and Particle Science 38, 161 Allen, M. P., 2004, in Computational Soft Matter: From Synthetic Polymers to Proteins, Lecture Notes, Vol. 23, pp 1–28, NIC Series 2004 Baade, W. & Zwicky, F., 1934, Proceedings of the National Academy of Science 20, 259 Baiko, D. A., 2011, ArXiv e-prints Brown, E. F., Bildsten, L., & Chang, P., 2002, Ap.J. 574, 920 Brown, E. F., Bildsten, L., & Rutledge, R. E., 1998, Ap.J. 504, L95 Chamel, N. & Haensel, P., 2008, Living Reviews in Relativity 11, 10 Chang, P., Bildsten, L., & Arras, P., 2010, Ap.J. 723, 719 Cheng, B., Epstein, R. I., Guyer, R. A., & Young, A. C., 1996, Nature 382, 518 Chiu, H.-Y. & Salpeter, E. E., 1964, Physical Review Letters 12, 413 Chugunov, A. I. & Horowitz, C. J., 2010, MNRAS 407, L54 Correa, A. A., Bonev, S. A., & Galli, G., 2006, Proceedings of the National Academy of Science 103, 1204 148 Bibliography Demorest, P. B., Pennucci, T., Ransom, S. M., Roberts, M. S. E., & Hessels, J. W. T., 2010, Nature 467, 1081 Gupta, S., Brown, E. F., Schatz, H., Möller, P., & Kratz, K., 2007, Ap.J. 662, 1188 Haskell, B., Jones, D. I., & Andersson, N., 2006, M.N.R.A.S. 373, 1423 Heger, A., Fryer, C. L., Woosley, S. E., Langer, N., & Hartmann, D. H., 2003, Ap.J. 591, 288 Helfand, D. J. & Long, K. S., 1979, Nature 282, 589 Helled, R., Anderson, J. D., Podolak, M., & Schubert, G., 2011, Ap.J. 726, 15 Hernquist, L. & Applegate, J. H., 1984, Ap.J. 287, 244 Hewish, A., Bell, S. J., Pilkington, J. D. H., Scott, P. F., & Collins, R. A., 1968, Nature 217, 709 Heyl, J. S. & Hernquist, L., 2001, MNRAS 324, 292 Hix, W. R. & Meyer, B. S., 2006, Nuclear Physics A 777, 188 Ho, W. C. G. & Heinke, C. O., 2009, Nature 462, 71 Hockney, R. W. & Eastwood, J. W., 1988, Computer simulation using particles Horowitz, C. J. & Berry, D. K., 2009, Phys. Rev. C 79(6), 065803 Horowitz, C. J., Berry, D. K., & Brown, E. F., 2007, Phys. Rev. E 75(6), 066101 Horowitz, C. J. & Kadau, K., 2009, Physical Review Letters 102(19), 191102 Horstemeyer, M. F., Baskes, M. I., & Plimpton, J., 2001, Acta Materialia 49, 4363 149 Bibliography Hurley, K., Boggs, S. E., Smith, D. M., Duncan, R. C., Lin, R., Zoglauer, A., Krucker, S., Hurford, G., Hudson, H., Wigger, C., Hajdas, W., Thompson, C., Mitrofanov, I., Sanin, A., Boynton, W., Fellows, C., von Kienlin, A., Lichti, G., Rau, A., & Cline, T., 2005, Nature 434, 1098 Hurley, K., Cline, T., Mazets, E., Barthelmy, S., Butterworth, P., Marshall, F., Palmer, D., Aptekar, R., Golenetskii, S., Il’Inskii, V., Frederiks, D., McTiernan, J., Gold, R., & Trombka, J., 1999, Nature 397, 41 Itoh, N. & Kohyama, Y., 1993, Ap.J. 404, 268 Jones, P. B., 2003, Ap.J. 595, 342 Kaplan, D. L., 2008, in Y.-F. Yuan, X.-D. Li, & D. Lai (ed.), Astrophysics of Compact Objects, Vol. 968 of American Institute of Physics Conference Series, pp 129–136 Kaspi, V. M., 2010, Proceedings of the National Academy of Science 107, 7147 Kittel, C., 1976, Introduction to solid state physics Lattimer, J. M. & Prakash, M., 2001, Ap.J. 550, 426 Lattimer, J. M. & Prakash, M., 2004, Science 304, 536 Lattimer, J. M. & Prakash, M., 2007, Physics Reports 442, 109 Lattimer, J. M., Prakash, M., Pethick, C. J., & Haensel, P., 1991, Physical Review Letters 66, 2701 Lattimer, J. M., van Riper, K. A., Prakash, M., & Prakash, M., 1994, Ap.J. 425, 802 Lyne, A. G., Shemar, S. L., & Smith, F. G., 2000, MNRAS 315, 534 Manchester, R. N., Hobbs, G. B., Teoh, A., & Hobbs, M., 2005, A.J. 129, 1993 150 Bibliography Mazets, E. P., Golentskii, S. V., Ilinskii, V. N., Aptekar, R. L., & Guryan, I. A., 1979, Nature 282, 587 Medin, Z. & Cumming, A., 2010, Phys. Rev. E 81(3), 036107 Meijer, E. J. & Frenkel, D., 1991, J. Chem. Phys. 94, 2269 Mereghetti, S., 2001, in F. Giovannelli & G. Mannocchi (ed.), Frontier Objects in Astrophysics and Particle Physics, pp 239–+ Mereghetti, S., 2008, Astronomy & Astrophysics Reviews 15, 225 Nadal, M.-H. & Le Poac, P., 2003, Journal of Applied Physics 93, 2472 Özel, F., Baym, G., & Güver, T., 2010a, Phys. Rev. D 82(10), 101301 Özel, F., Psaltis, D., Ransom, S., Demorest, P., & Alford, M., 2010b, Ap.J. (Letters) 724, L199 Page, D., Geppert, U., & Weber, F., 2006, Nuclear Physics A 777, 497 Plimpton, S., 1995, Journal of Computational Physics 117, 1 Pollock, E. L. & Hansen, J. P., 1973, Phys. Rev. A 8, 3110 Pons, J. A., Walter, F. M., Lattimer, J. M., Prakash, M., Neuhäuser, R., & An, P., 2002, Ap.J. 564, 981 Popov, S. B., 2008, Physics of Particles and Nuclei 39, 1136 Potekhin, A. Y., 2010, Physics Uspekhi 53, 1235 Potekhin, A. Y. & Chabrier, G., 2000, Phys. Rev. E 62, 8554 Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P., 1992, Numerical recipes in C. The art of scientific computing Rajagopal, M., Romani, R. W., & Miller, M. C., 1997, Ap.J. 479, 347 Robbins, M. O., Kremer, K., & Grest, G. S., 1988, J. Chem. Phys. 88, 3286 151 Bibliography Ross, M., 1981, Nature 292, 435 Ruderman, M., 1969, Nature 223, 597 Schneider, T. & Stoll, E., 1978, Phys. Rev. B 17, 1302 Seitenzahl, I. R., Timmes, F. X., Marin-Laflèche, A., Brown, E., Magkotsios, G., & Truran, J., 2008, Ap.J. (Letters) 685, L129 Shapiro, S. L. & Teukolsky, S. A., 1986, Black Holes, White Dwarfs and Neutron Stars: The Physics of Compact Objects, Black Holes, White Dwarfs and Neutron Stars: The Physics of Compact Objects, by Stuart L. Shapiro, Saul A. Teukolsky, pp. 672. ISBN 0-471-87316-0. Wiley-VCH , June 1986. Smoluchowski, R., 1970, Physical Review Letters 24, 923 Stairs, I. H., Lyne, A. G., & Shemar, S. L., 2000, Nature 406, 484 Stejner, M. & Madsen, J., 2005, Phys. Rev. D 72(12), 123005 Stringfellow, G. S., Dewitt, H. E., & Slattery, W. L., 1990, Phys. Rev. A 41, 1105 Swope, W. C., Andersen, H. C., Berens, P. H., & Wilson, K. R., 1982, J. Chem. Phys. 76, 637 Terrell, J., Evans, W. D., Klebesadel, R. W., & Laros, J. G., 1980, Nature 285, 383 Thompson, C. & Duncan, R. C., 1995, MNRAS 275, 255 Thorsett, S. E. & Chakrabarty, D., 1999, Ap.J. 512, 288 Timmes, F. X., 1999, Ap.J. (Suppl) 124, 241 Timmes, F. X., Hoffman, R. D., & Woosley, S. E., 2000, Ap.J. (Suppl) 129, 377 Ushomirsky, G., Cutler, C., & Bildsten, L., 2000, MNRAS 319, 902 152 Bibliography Vaulina, O., Khrapak, S., & Morfill, G., 2002, Phys. Rev. E 66(1), 016404 Verlet, L., 1967, Physical Review 159, 98 Walter, F. M., Pons, J. A., Burwitz, V., Lattimer, J. M., Lloyd, D., Wolk, S. J., Prakash, M., & Neuhäuser, R., 2004, Advances in Space Research 33, 513 Woods, P. M. & Thompson, C., 2006, in Lewin, W. H. G. & van der Klis, M. (ed.), Compact stellar X-ray sources, pp 547–586, Cambridge: Cambridge Univ. Press. Yuan, J. P., Manchester, R. N., Wang, N., Zhou, X., Liu, Z. Y., & Gao, Z. F., 2010, Ap.J. (Letters) 719, L111 Zavlin, V. E. & Pavlov, G. G., 2002, in W. Becker, H. Lesch, & J. Trümper (eds.), Neutron Stars, Pulsars, and Supernova Remnants, pp 263–+ 153

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            data-media="{[{embed.selectedMedia}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0072379/manifest

Comment

Related Items