UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Magnetohydrodynamics of neutron star interiors Easson, Ian Whiteman 1975

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
831-UBC_1975_A1 E28_6.pdf [ 8.62MB ]
Metadata
JSON: 831-1.0084835.json
JSON-LD: 831-1.0084835-ld.json
RDF/XML (Pretty): 831-1.0084835-rdf.xml
RDF/JSON: 831-1.0084835-rdf.json
Turtle: 831-1.0084835-turtle.txt
N-Triples: 831-1.0084835-rdf-ntriples.txt
Original Record: 831-1.0084835-source.json
Full Text
831-1.0084835-fulltext.txt
Citation
831-1.0084835.ris

Full Text

MAGNETOHYDRODYNAMICS OF NEUTRON STAR INTERIORS by IAN WHITEMAN EASSON B.Sc. U n i v e r s i t y of B r i t i s h Columbia, 1969 M.Sc. U n i v e r s i t y of B r i t i s h Columbia, 1971 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY i n the Department of PHYSICS We accept t h i s thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA A p r i l , 1975 I n p r e s e n t i n g t h i s t h e s i s i n p a r t i a l f u l f i l m e n t o f t h e r e q u i r e m e n t s f< a n a d v a n c e d d e g r e e a t t h e U n i v e r s i t y o f B r i t i s h C o l u m b i a , I a g r e e t h a t h e L i b r a r y s h a l l m a k e i t f r e e l y a v a i l a b l e f o r r e f e r e n c e a n d s t u d y . I f u r t h e r a g r e e t h a t p e r m i s s i o n f o r e x t e n s i v e c o p y i n g o f t h i s t h e s i s f o r s c h o l a r l y p u r p o s e s may b e g r a n t e d b y t h e H e a d o f my D e p a r t m e n t o r b y h i s r e p r e s e n t a t i v e s . I t i s u n d e r s t o o d t h a t c o p y i n g o r p u b l i c a t i o n o f t h i s t h e s i s f o r f i n a n c i a l g a i n s h a l l n o t b e a l l o w e d w i t h o u t my w r i t t e n p e r m i s s i o n . D e p a r t m e n t o f T h e U n i v e r s i t y o f B r i t i s h C o l u m b i a V a n c o u v e r 8, C a n a d a D a t e 6^1 , /c?7\ i ABSTRACT The dynamics of the charged p a r t i c l e s i n the f l u i d i n t e r i o r of rotat i n g magnetized neutron stars (pulsars) i s investigated. I t i s shown that a magnetohydrodynamic approach i s v a l i d under a wide v a r i e t y of conditions. The small amplitude waves that can propagate i n the "charged" < - 4 f l u i d are sound waves (period ^> 10 sec), i n e r t i a l waves (- 1 sec) and hydromagnetic-inertial waves (£ months). Generally, the most e f f e c t i v e damping mechanism i s v i s c o s i t y . Viscous damping times f o r hydromagnetic-i n e r t i a l waves can be as long as hundreds of years. Several normal sound mode and normal i n e r t i a l mode c a l c u l a t i o n s are performed. The motion of the charged f l u i d over the time scales of i n e r t i a l waves and hydromagnetic-inertial waves i s simulated on a computer. Mode-mode coupling, boundary layers, and the ac c e l e r a t i o n of the crust due to charged f l u i d motion are also studied. The most important conclusion i n th i s thesis i s that the time the charged f l u i d as a whole takes to respond to a sudden disturbance such as a change i n the angular v e l o c i t y or acc e l e r a t i o n of the crust i s the time f o r a long wavelength hydromagnetic-inertial wave to cross the st a r . Contrary to an assumption which has been made i n models of pulsar post-g l i t c h behaviour, the charged f l u i d response time i s not nec e s s a r i l y small compared to p o s t - g l i t c h r e l a x a t i o n times. In view of t h i s , models of p o s t - g l i t c h behaviour which assume that the charged f l u i d response time i s small compared to p o s t - g l i t c h r e l a x a t i o n times should be re-examined. TABLE OF CONTENTS A b s t r a c t .' T a b l e Of C o n t e n t s L i s t Of T a b l e s . . . L i s t Of F i g u r e s . . . . . . . Acknowledgments CHAPTER I: INTRODUCTION . Rev iew Of N e u t r o n S t a r P r o p e r t i e s The S o l i d C o r e The F l u i d I n t e r i o r The C r u s t • • • — The Magnetosphere Rev iew Of P u l s a r T i m i n g F e a t u r e s G e n e r a l T i m i n g F e a t u r e s T i m i n g I r r e g u l a r i t i e s Charged F l u i d Dynamic s : What Has Been Done And What Needs To Be Done What Has Been Done • What Needs To Be Done CHAPTER I I : THE MODEL OF THE CHARGED FLUID M a g n e t o h y d r o s t a t i c E q u i l i b r i u m Of A R o t a t i n g S t a r The E q u a t i o n s Of M a g n e t o h y d r o s t a t i c E q u i l i b r i u m The E q u a t i o n Of S t a t e . . . . . A N u m e r i c a l Example The M a g n e t i c F i e l d S t r u c t u r e •••• Charged F l u i d E q u a t i o n s Of M o t i o n The Range Of V a l i d i t y Of A F l u i d D e s c r i p t i o n The Mass E q u a t i o n The Momentum Equations The Gravitational Source Equation The Magnetic Induction Equations : The Boundary Conditions The Velocity Boundary Conditions The Magnetic Field Boundary Conditions ........ The Crust And Core Equations ;. ... CHAPTER III: TYPES OF CHARGED FLUID WAVES ........ Dispersion Relations Sound Waves * Inertial Waves And Hydromagnetic-inertial Waves Inertial Waves Hydromagnetic-inertial Waves Damping Due To Viscosity And Superfluid Drag .... Wave Damping By Heat Conduction Resonant Mode-Mode Coupling Boundary Layers Simplified Charged Fluid Equations Sound Time Scales Inertial Time Scales Hydromagnetic-inertial Time Scales CHAPTER IV: SOUND MODES CHAPTER V: INERTIAL MODES AND INERTIAL FLOW ....... Free Inertial Modes Of Oscillation Cylindrical Geometry ...<.-Inertial Flow First Flow Simulation Second Flow Simulation i v CHAPTER VI: MAGNETOS TROPHIC FLOW .139 Equations For The Velocity Field...... 142 The Charged Fluid Response Time Problem .. 145 The Magnetic Field Structure 147 Flow Simulation 152 Basic Time Scales Of The Flow 153 Details Of The Flow *. . 155 CONCLUDING REMARKS . 181 BIBLIOGRAPHY .... , 187 APPENDIX A: NUMERICAL METHODS FOR SOUND WAVES 190 APPENDIX B: NUMERICAL METHODS FOR INERTIAL FLOW 192 Spatial Mesh Systems 192 The Coupled Pressure And Potential Equations 195 Time Differencing And Time Centering 196 APPENDIX C: NUMERICAL METHODS FOR MAGNETOSTROPHIC FLOW 198 V LIST OF TABLES I. Overall Structure Of Neutron Stars 2 II. Timing Parameters For NP0532 10 III. Non-Dimensional Parameters 20 IV. Some Microscopic Parameters Of The Fluid Region 22 V. The Parameters Of A Rotating Neutron Star Model ... 29 VI. Forces Acting On The Charged Fluid ...... 37 VII. Eigenfrequencies Of Some Standing Sound Waves 87 VIII. Eigenvalues - A 2 For Some Axisymmetric Inertial Modes 99 IX. Eigenfrequencies Of Some Buoyancy-Inertial Modes 102 X. Charged Fluid Waves And Their Characteristics ............... 183 XI. Convergence Of Sound Wave Eigenf requencies 191 v i LIST OF FIGURES 1. A Schematic Picture Of A Rotating Magnetized Neutron Star . . . . . . 3 2. Charged Fluid Mass Fraction As A Function Of Total Density x In The Fluid Region 57 3. Logs Of Squares Of Sound Speeds As A Function Of Total Density In The Fluid Region 59 4. And As Functions Of Radius In A Rotating Neutron o o Star Model 61 5. Toroidal Magnetic Field Perturbation Required To Balance The Quasi-Steady Poincare Force 63 6. Radial Parts Of Some £ •<* 0 Charged Fluid Sound Modes . . . . . . . . . . 88 7. Radial Parts Of Some JL = 1 Charged Fluid Sound Modes 90 8. Radial Parts Of Some J! = 2 Charged Fluid Sound Modes 92 9. Radial Parts Of Some Axisymmetric Inert ia l Modes . . . . . . . . . . . . . . 114 10. Fluid Flow In The Highest Frequency m = 1 Buoyancy-Inertial Mode ,. 116 11. Fluid Flow In The Highest Frequency m =2 Buoyancy-Inertial Mode 118 12. Fluid Flow In.The Highest Frequency m = 3 Buoyancy-Inertial Mode 120 13. Numerical Grids Used In The Inert ia l Flow Simulation 122 14. Inert ial Flow After 1.98 Rotations In The F i r s t Simulation 124 15. Inert ial Flow At Various Times In The Second Simulation 126 16. Core And Crust Angular Acceleration Due To Shear Stress In Second Inert ial Flow Simulation 135 17. Power Spectrum Analysis Of Crust Angular Acceleration In Second Inert ial Flow Simulation 137 18. Perturbation In Crust Angular Acceleration Due To Viscous Torque 159 19. Power Spectrum Analysis Of Crust Angular Acceleration Data In Figure 18. 161 20. Magnetostrophic Flow Of Charged Fluid Induced By Crustal Deceleration 163 v i i ACKNOWLEDGMENTS I would like to thank Prof. M.H.L. Pryce for his constant encouragement and guidance. He was particularly helpful in keeping me from going on wild goose chases, and was always available with good advice whenever I needed i t . I thank my external examiner and a l l my committee members, and in particular Greg Fahlman, who gave me useful advice on power spectrum analysis. I thank the staff of the U.B.C. Computing Centre and the faculty of the Computer Science Department, including Carol Bird and Jim Varah, for useful discussions on numerical methods. For the same reason I thank Jeanie Eilik in astronomy, who had numerical problems similar to mine. I have had useful conversations with many people at the University of Illinois, including Gordon Baym, Ron Eisner, Don Lamb, Chris Pethick, Dave Roberts, Roger Smith and particularly with Fred Lamb. This research was supported by a National Research Council of Canada post-graduate scholar-ship. Finally, and most of a l l , I would like to thank my wife Elaine for typing the thesis and for keeping my spirits up a l l the way along. CHAPTER I: INTRODUCTION Very soon a f t e r the discovery of pulsars i t was r e a l i z e d that there i s probably only one kind of object, which can display the observed pulsar c h a r a c t e r i s t i c s . That object i s a r o t a t i n g neutron star which possesses a large magnetic f i e l d . This thesis i s concerned with one p a r t i c u l a r feature of r o t a t i n g , magnetized neutron s t a r s : the magnetohydrodynamics of the charged f l u i d i n t h e i r i n t e r i o r s . Very l i t t l e work has hitherto been done i n this area. Accordingly, the f i r s t part of the thesis i s devoted to a c a r e f u l i n v e s t i -gation of various aspects of the charged f l u i d dynamics and to estimates of t h e i r importance. The second part contains d e t a i l e d numerical calcu-l a t i o n s of some possible charged f l u i d flows and modes of o s c i l l a t i o n . This introductory chapter begins with a short review of some theoret-i c a l l y predicted properties of neutron s t a r s , with emphasis on those properties which are most relevant to the charged f l u i d dynamics. 1 Next, some features of pulsar observations are reviewed, with emphasis on those aspects which one hopes w i l l reveal something about the i n t e r i o r . Previous work on the charged f l u i d dynamics i s c r i t i c a l l y reviewed. F i n a l l y , a number of important questions about the dynamics are posed. In t h i s thesis I w i l l attempt to answer these questions. 1 For more de t a i l e d information on neutron stars and pulsars, see the review a r t i c l e s by ter Haar(1972), Ruderman(1972) and Ginzberg(1971). Review of Neutron Star Properties The neutron star structures and properties described i n t h i s section are based on the t h e o r e t i c a l c a l c u l a t i o n s of many people and not for the most part on d i r e c t observational evidence, because while there i s a wealth of information contained i n pulsar observations, i t i s not easy to extract from these observations information about neutron star i n t e r i o r s . Never-theless, a consensus has emerged over the l a s t several years on many aspects of neutron star structure. This view i s presented below i n general terms. Many points are discussed more thoroughly and q u a n t i t a t i v e l y i n Chapter I I . A general overview of neutron star structure i s given i n Table I. Table I Over a l l Structure Of Neutron Stars Radius Composition Mass Average density Temperature Surface magnetic f i e l d Rotational period Inner structure Surrounded by about 10 km mainly neutrons about one solar mass roughly nuclear ( I O 1 4 gm cm - 3) 10 6 to 10 8 K 1 0 1 0 to I O 1 4 gauss 30 msec to a few seconds s o l i d core (?) f l u i d i n t e r i o r s o l i d crust magnetpsphere Figure 1 shows a schematic picture of a neutron st a r . It may help the reader i f he r e f e r s to Table I and to Figure 1 while reading the present s e c t i o n . I t i s convenient f o r the purposes of t h i s thesis to d i s t i n g u i s h three regions i n the i n t e r i o r of a neutron s t a r : a s o l i d core, a f l u i d i n t e r i o r , and a s o l i d crust. These three regions and the region immediately surrounding the s t a r , the magnetosphere, w i l l now be discussed i n turn. 3 Figure 1 . A schematic picture of a rotating magnetized neutron star. NEUTRON SUPERFLUID AND CHARGED FLUID 5 The S o l i d Core > A s o l i d core may or may not e x i s t at the centre of some neutron s t a r s . . The uncertainty a r i s e s mainly from the unknown nature of the i n t e r p a r t i c l e 15 -3 i n t e r a c t i o n at the high d e n s i t i e s there (about 10 gm cm ). If present, the s o l i d core i s expected to consist of a neutron c r y s t a l . Detailed c a l c u l a t i o n s of the c r y s t a l properties have been made by Canuto and Chitre 36 (1974). They f i n d that the material moduli of the c r y s t a l are about 10 -2 dyne cm. , which implies a very r i g i d core. The l a t e s t most refined c a l c u -l a t i o n s i n t h i s area are those of Pandharipande and.Smith (1975), who f i n d that s o l i d cores are possible i n stars with masses within .1 M Q of the maximum mass f o r a stable star (Smith, p r i v a t e communication). For t h e i r two d i f f e r e n t p o t e n t i a l s , t h i s maximum mass i s 1.8 and 2.3 M . I t i s uncertain i f neutron stars i n t h i s mass range can be formed. The F l u i d I n t e r i o r More than 90% of the p a r t i c l e s i n the f l u i d i n t e r i o r are degenerate n o n - r e l a t i v i s t i c neutrons which form a s u p e r f l u i d of the Bardeen-Cooper-S c h r i e f f e r type. The s u p e r f l u i d "rotates" i n the sense that there are v o r t i c e s i n whose cores the neutrons are normal. The most widely accepted view i s that the s u p e r f l u i d i s approximately uniformly r o t a t i n g , with aligned v o r t i c e s (Ruderman, 1970; Ruderman and Sutherland, 1974). The r e s t of the f l u i d i n t e r i o r consists of charged p a r t i c l e s , p r i m a r i l y of equal numbers of protons and electrons. The protons are degenerate, non-r e l a t i v i s t i c and probably form a type II superconductor. The electrons are degenerate and r e l a t i v i s t i c but not superconducting (Baym et_ al_. , 1969b) . Because of t h e i r strong tendency towards l o c a l charge n e u t r a l i t y , i t i s natural to treat these charged p a r t i c l e s as a sin g l e f l u i d - a "charged" f l u i d . The dynamics of t h i s charged f l u i d i s the subject of t h i s t h e s i s . The only d i r e c t i n t e r a c t i o n between the sup e r f l u i d and the charged f l u i d comes about by the s c a t t e r i n g of electrons off those few neutrons -18 which are not i n a s u p e r f l u i d state (about 10 of a l l neutrons). The small f r a c t i o n of protons which i s not i n the superconducting state can also scatter off the normal neutrons. The f l u i d i n t e r i o r can thus be thought of as being composed of two interpenetrating, very Weakly i n t e r -acting f l u i d s . Some of the transport c o e f f i c i e n t s of these f l u i d s depend s e n s i t i v e l y on the temperature of the i n t e r i o r . For example, the v i s c o s i t y calculated by Flowers and Itoh (1975) i s <* T . Their paper contains the most c a r e f u l and up-to-date c a l c u l a t i o n s of transport c o e f f i c i e n t s i n neutron stars, so I w i l l use t h e i r r e s u l t s i n t h i s t h e s i s . 9 The temperature i n the i n t e r i o r may vary from about 10 K for neutron stars a few hours old to perhaps as low as 10~* K f o r extremely old ones (Tsuruta and Cameron, 1965). Most observed pulsars are so o l d that t h e i r g i n t e r i o r temperatures are £ 10 K. The thermal conductivity of the e l e c -24 8 trons i s so high (about 10 erg/cm-deg-sec at T = 10 K; Flowers and Itoh, 1975) that the f l u i d i n t e r i o r i s p r a c t i c a l l y isothermal. Estimates of the I n t e r i o r temperature of a neutron star of a given mass and age are i n a very confused state, because the temperature i s s e n s i t i v e to a large number of factors such as magnetic f i e l d strength, s u p e r f l u i d i t y , pion condensation, and so on. Greenstein and McClintock (1974) have recently reviewed the problem of neutron star temperatures, both from a t h e o r e t i c a l and an observational point of view. There i s a strong magnetic f i e l d . i n the i n t e r i o r . Just how strong i t i s and what i t s structure i s are not known. There i s f a i r l y general agree-ment, however, that t h i s f i e l d i s a remnant of a seed f i e l d which was greatly amplified i n the collapse which formed the neutron s t a r . A somewhat s i m p l i s t i c argument i l l u s t r a t e s this idea. If a star of 11 2 10 cm radius which possesses a p o l o i d a l f i e l d of 10 gauss collapses to form a neutron star of 10^ cm radius, then according to the argument, con-servation of magnetic f l u x Implies the r e s u l t i n g f i e l d should have a strength 2 11 6 2 12 of 10 x ( 1 0 / 1 0 ) =10 gauss. Two reasons why t h i s post-collapse value i s u n r e l i a b l e are that the choice of i n i t i a l f i e l d strength i s somewhat a r b i t r a r y and that the a m p l i f i c a t i o n of magnetic f i e l d by d i f f e r e n t i a l r o t a -t i o n during collapse i s ignored. Such d i f f e r e n t i a l r o t a t i o n can produce f i n a l f i e l d strengths which are orders of magnitude larger than those suggested by the f l u x conservation argument (LeBlanc'and Wilson, 1970). Adding to the uncertainty i s the p o s s i b i l i t y that j u s t a f t e r collapse, while the star i s completely f l u i d , large magnetic f i e l d s i n the i n t e r i o r can be transported to the surface by f l u i d i n s t a b i l i t i e s , where the f i e l d s can decay much more r a p i d l y than they could i n the i n t e r i o r . After the star cools s u f f i c i e n t l y f or a s o l i d crust to form, which may take about an hour (Tsuruta and Cameron, 1965), t h i s process i s no longer possible. A f t e r crust formation, the magnetic f i e l d i n the i n t e r i o r can only decay by r e s i s t i v e d i f f u s i o n . A magnetic f i e l d with length scale X w i l l 2 2 decay i n a time - 4 i r a X /c , where a i s the e l e c t r i c a l conductivity. How-29 -1 ever, the conductivity i s so high i n the i n t e r i o r (about 10 sec for a temperature T ='10** K; Baym e t ' a l . , 1969b) that the decay time for \ - 10^ cm i s much longer than the age of the universe. An estimate of the i n t e r i o r magnetic f i e l d strength i n the neutron star forming the Crab pulsar was given by Ostriker and Gunn (1969). Their 8 magnetic dipole r a d i a t i o n model of pulsar slowing down predicted that the Crab pulsar i s - 25% older than i t i s known to be. They a t t r i b u t e d the age discrepancy to the neglect of g r a v i t a t i o n a l quadrupole r a d i a t i o n which i s continually emitted because of a small equatorial e l l i p t i c i t y . They estimated that an i n t e r i o r p o l o i d a l f i e l d of - 10''''' gauss would produce an e l l i p t i c i t y s u f f i c i e n t to eliminate the discrepancy. However, Gunn and Ostriker assumed that the i n i t i a l angular speed fi(0) of the pulsar was very much greater than i t s present value £3(t) . They 2 require t h i s assumption i n order to drop terms of order [fi(t)/fi(0)] from the i r age equation, which then allows them to predict an age for the pulsar which i s independent of. n(0). If they had used an £2(0) which i s 2.3fi(t), rather than a number >> fi(t), there would have been no discrepancy at a l l . There are several good reasons why fi(0) = 2.3fi(t) i s a reasonable assump-ti o n (Ruderman, 1972). I t i s thus quite unnecessary to invoke an i n t e r n a l f i e l d of 10'"^  gauss i n order to s a t i s f a c t o r i l y explain the age discrepancy. Equally important, t r i a x i a l i t y due to the s o l i d crust or to a possible s o l i d core could also produce the required equatorial e l l i p t i c i t y . Whatever the structure of the i n t e r i o r f i e l d and whatever i t s exact strength, i t i s almost c e r t a i n l y strong enough to have an important influence on the charged f l u i d dynamics. The Crust The crust i s a coulomb l a t t i c e composed of neutron r i c h n u c l e i bathed i n an electron gas. I t i s between a few hundred meters and several km thick. 29 31 -2 The bulk and shear moduli of t h i s s o l i d range from 10 to 10 dynes cm 14 -3 The density v a r i e s from about 2 x 10 gm cm at the f l u i d i n t e r i o r - c r u s t interface to near zero at the surface. At densities higher than about 11 ' -3 10 gm cm , the solid is expected to be permeated with a neutron superfluid. The electrical resistivity n = c /4ira due to electron-phonon scattering 2 2 -1 '. in the crust falls steadily from about 10 cm sec near the surface to -9 2 -1 7 about 10 cm sec near the crust-fluid boundary, for T = 10. K (Flowers and Itoh, 1975). The resistivity due to electron-phonon scattering is <* T. The resistivity due to impurity scattering is temperature independent, and 2 is proportional to <(AZ) > X ^ » where x ^ is the impurity concentration and 2 <(AZ) > is the mean square charge difference between the impurities and the -3 2 surrounding nuclei. For x. , ~ 10 and <(AZ) > - 1, the impurity scattering resistivity is roughly equal to the electron-phonon resistivity at T = 10^  K throughout the crust. Taking a simple model in which ri varies inversely with distance from 3 - 1 the surface at a rate a cm sec , dimensional considerations suggest that the decay time x for the magnetic field in a crust of thickness X is -3 "~A 3 ™"1 X2 T - X /a. For a = 10 cm sec and X = 1km, T is about 10 years. All pulsars whose ages are known or have been estimated are younger than several million years. A somewhat more realistic model of the dependence of n on distance from the surface will not change the conclusion that the magnetic field in the crust will not decay significantly over the lifetime of the pulsar. Thus, the field can be considered to be effectively "frozen" into the crust. 10 The Magnetosphere The magnetosphere is the plasma outside the surface of the star. It owes its existence to the large magnetic field frozen into the highly elec-trically conducting rotating star. In most models, a portion of the mag-netosphere co-rotates rigidly with the crust, and pulses of electromagnetic radiation are emitted in time with the rotation of the crust. The signif-icance of this for studies of the interior is that pulsar timing allows one to determine the angular speed of the crust, from which i t may be possible to infer something about the interior. Review of Pulsar Timing Features General Timing Features Pulsars are very steady clocks. However, if the pulse arrival times for a pulsar are recorded for a period of several months, i t is observed that on the average the pulsar slows down. As an illustration of timing data, the timing parameters for pulsar NP0532 (the Crab) are given in Table II. The Crab is the fastest pulsar known, and primarily for this reason i t is also the only pulsar for which ft is measurable. Table II Timing Parameters For NP0532 Q = h = 1.9 x 102 sec - 1 -2.4 x 10"9 6.9 x IO - 2 0 sec - 3 —? sec ^  ec" The slowing down of a pulsar can be understood i n terms of the braking action of the electromagnetic torque acting on the s o l i d surface of the st a r . The s i z e of the surface magnetic f i e l d can be estimated as follows: assume a magnetic f i e l d configuration (e.g., a "rot a t i n g d i p o l e " ) . Calculate the Maxwell st r e s s tensor at the surface and, by i n t e g r a t i o n , the e l e c t r o -magnetic torque. Assume a p e r f e c t l y r i g i d star and c a l c u l a t e the angular deceleration = torque/(moment of i n e r t i a of s t a r ) . Comparing this number with the observed deceleration allows the surface magnetic f i e l d s i z e to be estimated. Two comments about t h i s procedure are i n order. The f i r s t i s that due to large uncertainties and inconsistencies i n the present models of the magnetosphere, the deduced values of the magnetic f i e l d are uncertain by a few orders of magnitude. For example, the model of Ostriker and Gunn (1969) implies that surface magnetic f i e l d s f o r most pulsars are a few times 12 10 gauss, whereas the dense magnetosphere model of Roberts and Sturrock (1972) implies surface magnetic f i e l d s about an order of magnitude smaller than t h i s . The second comment i s that neutron stars are not p e r f e c t l y r i g i d . The preceding view of a uniformly slowing down s t a r i s an o v e r s i m p l i f i c a t i o n . It might be expected that there are i r r e g u l a r i t i e s i n pulsar timings a r i s i n g from the lack of r i g i d i t y of neutron s t a r s . Timing I r r e g u l a r i t i e s If the pulse a r r i v a l times f o r a pulsar are recorded f o r about a year or longer and the r e s u l t s f i t t e d by any smooth function such as a cubic or quartic polynomial, there are sometimes s i g n i f i c a n t d i f f e r e n c e s between the observed and f i t t e d a r r i v a l times. These d i f f e r e n c e s , c a l l e d timing r e s i d u a l s , 12 are measurable for the Crab pulsar (Drakej 1972), the Vela pulsar PSR0833 (Nelson et al.,1970; Reichley, 1972), and several other pulsars (Manchester et al., 1974; Manchester and Taylor, 1974). There appear to be two distinct kinds of irregularities. In the Crab, Vela and PSR1508+55 pulsars, the residuals show large "glitches", that i s , sudden changes in the rotational period of the star. Glitches are followed by a relaxation phase, lasting a few days for the Crab and about a year for the Vela, in which the sudden change in period relaxes exponentially. These glitches have been interpreted as due to cracking in the solid crust or in a solid core (Ruderman, 1969; Baym and Pines, 1971; Shaham et al., 1973)'. Very simple models of the superfluid motion following such glitches are able to account well for the relaxation behaviour immediately following the glitch (Baym et al., 1969c), and this agreement has been interpreted as evidence for an I n t e r i o r which is at least partly superfluid. This post-glitch theory assumes that both superfluid and charged fluid can be treated as rigidly rotating, with the charged fluid co-rotating with the solid parts of the star, and that the observed relaxation time i s the time i t takes the superfluid to adjust to the glitch. In support of the assumption that over the relaxation time scale i t is a good approximation to treat the charged fluid as if i t rigidly co-rotates with the solid parts of the star, Baym et al. (1969c) argue that within roughly the time i t takes 12 an Alfv6n wave to cross the star, about 100 sec for a 10 gauss magnetic field, the charged fluid can readjust itself to changes, in the angular velocities of the solid parts of the star. The other kind of timing irregularity is present in some pulsars in which large glitches have not yet been seen, as well as in those in which large glitches have been seen. The most complete residual data are for 13 the Crab i n the o p t i c a l (Boynton et a l . , 1972; Groth, 1975a,b; Groth, private communication). The r e s i d u a l s are quasi-sinusoidal, with a period between about one and two thirds the length of the data, for any given data set. The quasi-sinusoidal period thus grows longer as more data i s c o l l e c t e d . A power spectrum analysis of the residuals shows that the power at frequency -4 OJ i s roughly proportional to OJ . The amplitude of the residuals also grows with the length of the data, from several 100 ysec for about half a year of data to a few msec f o r a few years of data. Two kinds of explanations f o r the quasi-sinusoidal behaviour seem possib l e . The f i r s t (Boynton e t ' a l . , 1972) assumes that some random process i s occuring. This assumption i s consistent with the growth df the period of the quasi-sinusoidal residuals with the length of the data. Furthermore, i f the process i s a random walk i n the frequency of the crust a power -4 spectrum « GO would r e s u l t , i n agreement with observations. Such a f r e -quency walk could be produced by micro-glitches. The second type of explanation i s that the neutron star i s o s c i l l a t i n g simultaneously i n a large number of long l i v e d normal modes. For example, Ruderman (1972) assumed that the s u p e r f l u i d was o s c i l l a t i n g mainly i n i t s longest.period Tkachenko mode. His theory was able to account for the amplitude and period of the quasi-sinusoidal r e s i d u a l s known at that time. Since then, the amplitude and period have increased as more data has been gathered so the s u p e r f l u i d o s c i l l a t i o n model has become untenable. This example i l l u s t r a t e s one basic problem of the second kind of explanation: the longest period normal mode must have a period ^ the length of the data i n order f o r the residuals to appear to be the r e s u l t of noise. Data c o l l e c t i o n ceased i n A p r i l , 1974 (Groth, p r i v a t e communication). A -4 second basic problem Is to account f o r the w trend of the power spectrum i n a nautral way. 14 Charged F l u i d Dynamics: What Has Been Done And What Needs To Be Done What Has Been Done Not very much research has been done on.the charged f l u i d dynamics. One l i n e of research has been the i n v e s t i g a t i o n of magnetohydrodynamic i n s t a b i l i t i e s i n neutron s t a r s . This work was i n i t i a t e d by Vandakurov (1972) who noted that for the highly degenerate material i n neutron stars the pressure depends only very weakly on the temperature,' and who speculated that as a r e s u l t an a r b i t r a r i l y weak magnetic f i e l d could cause convection. This was investigated further by Tayler (1973b) , who considered the possible e f f e c t s of the very small thermal corrections to the equation of sta t e . He concluded that i t i s necessary f o r the magnetic f i e l d B to be ^ 5 x 10^'T gauss i n order for the star to be convectively unstable. Tayler's (1973b) work was c r i t i c i z e d by Chanmugam (1974), who pointed out that Tayler had t a c i t l y assumed that the growth time x of the i n s t a b i l i t y was >> the time x^ for convected matter to reach nuclear equilibrium. Since nuclear equilibrium i s established by the weak i n t e r a c t i o n , Chanmugam argued that x << x^ i s more l i k e l y . Assuming t h i s ordering of time scales, he found that a necessary condition f o r convective i n s t a b i l i t y i s B ^ 10"*"^  gauss, independently of the temperature. There are two major f a c t o r s that have not been taken into account i n a l l t h i s work on convective i n s t a b i l i t y . One i s r o t a t i o n , which generally tends to s t a b i l i z e f l u i d s (Acheson and Hide, 1973) and, which as Hide (1971) points out, i s very important i n the dynamics of the f l u i d core of neutron s t a r s . The other i s the highly conducting s o l i d crust through which magnetic l i n e s of force thread. The l i n e s of force tend to be pinned at the s o l i d - f l u i d i n t e r f a c e . This " l i n e - t y i n g " tends to greatly s t a b i l i z e 15 a magnetic f i e l d configuration ( K r a l l and T r i v e l p i e c e , 1973). Vandakurov (1970) attempted to f i n d some of the normal modes and eigenfrequencies of a completely f l u i d r o t a t i n g neutron star with a t o r o i d a l magnetic f i e l d that vanished at the surface. His ca l c u l a t i o n s were performed for values of a dimensionless parameter which, i f the star i s rot a t i n g at the rate of the Crab pulsar, imply he i s considering f i e l d s of about 1 0 ^ gauss. He i n c o r r e c t l y neglected bouyancy forces which are the same order of magnitude as other e f f e c t s which he included, despite claims to the contrary i n h i s paper. This could be a serious omission f o r the non-axisymmetric modes he considers. He also i n c o r r e c t l y assumed, i n e f f e c t , that the neutrons as w e l l as the charged f l u i d can respond to the J x B magnetic body f o r c e . As a r e s u l t , the unperturbed d e n s i t i e s which appear i n h i s equations are too large by about a factor = 100. The d e t a i l e d c a l c u l a t i o n s of Vandakurov (1970) thus seem to be i n c o r r e c t , although his conclusion that the t o r o i d a l f i e l d he considered i s unstable may be correct. Heinzman je_t aJ-. (1973) attempted to determine the reaction of ro t a t i n g s t a r s , both magnetized and unmagnetized, to externally applied torques. An e s s e n t i a l s i m p l i f i c a t i o n i n t h e i r paper i s based on t h e i r claim that for isochoric r o t a t i o n a l motion, i n which V • v = 0, the non-linear term v • V v i n the f l u i d momentum equation written i n an i n e r t i a l frame i s cancelled by gradient terms. This i s equivalent to saying that i n a frame of reference which i s co-rotating with the mean angular v e l o c i t y of the s t a r , not only i s the c e n t r i f u g a l force conservative (which i t i s ) , but so i s the C o r i o l i s force (which i t i s not). This fundamental error r e s t r i c t s the v a l i d i t y of t h e i r c a l c u l a t i o n s to the case of neutron st a r s which to a f i r s t approxi-mation do not rotate. 16 Hide (1971), i n probably the most important work i n t h i s area to date, pointed out that r o t a t i o n has an extremely important influence on the dynamics of the charged f l u i d . As a r e s u l t , two possible types of waves that can e x i s t i n the charged f l u i d are i n e r t i a l waves and hydromagnetic-i n e r t i a l waves. Hide did not inv e s t i g a t e these waves i n any d e t a i l . What Needs To Be Done Li s t e d below are some of the important questions about the charged f l u i d dynamics which need to be answered and the reasons why I consider them important. I t i s these questions that I attempt to answer i n t h i s t h e s i s . 1. What i s the range of v a l i d i t y of the f l u i d d e s c r i p t i o n of the charged p a r t i c l e s i n the i n t e r i o r ? Within t h i s regime, to what extent are the usual equations of magnetohydrodynamics v a l i d ? The motivation f o r these questions i s p r a c t i c a l because a f l u i d approach to a plasma problem i s usually considerably easier than a k i n e t i c approach. Since the highly degenerate plasma i n the i n t e r i o r i s i n many respects d i f f e r e n t from non-degenerate laboratory plasmas, i t i s necessary to determine under what conditions, i f any, the usual magnetohydrodynamic equations are applicable to the charged f l u i d . 2. What types of small amplitude waves can e x i s t i n the charged f l u i d ? What can damp them? How much are they damped? How do the answers to these questions depend on the gross parameters of the problem, such as fi, B, and T? I m p l i c i t i n these questions i s the approximation that the charged f l u i d can be treated as uniform and i n f i n i t e i n extent. Hide (1971) provided a p a r t i a l answer to the f i r s t of these questions. 17 3. What kinds of mode-mode coupling can occur? Can this weakly non-linear effect lead to heavy damping of a wave which is only lightly damped in the linear approximation? The answers to these questions are important in determining to what extent the calculated linear damping rates of waves are accurate. 4. What are some normal modes of oscillation of the charged fluid? This question is a more realistic version of question 2 in that i t takes into account the Inhomogeneity of the fluid and the effects of boundary conditions. 5. How do charged fluid oscillations affect the crust? What is the relationship, if any, between the charged fluid and the quasi-sinusoidal residuals? The answer to these questions are important because i t is only via the observed crust angular velocity that one can learn anything observationally about the dynamics of the interior. 6. How does the charged fluid respond under a variety of in i t i a l conditions and external forces? In particular, how does i t respond and how quickly does i t respond to changes in the crust angular velocity? The question of the extent to which the charged fluid can rigidly slow down with the crust is important. Furthermore, there is a published treatment of the problem (Heinzman et al., 1973) which is incorrect. How rigidly the charged fluid acts following a glitch is extremely important to the question of whether or not post-glitch relaxation constitutes strong evidence in favour of a partly superfluid Interior. 18 CHAPTER I I : THE MODEL OF THE CHARGED FLUID In t h i s chapter, the basic assumptions and approximations that w i l l be made about the charged f l u i d and i t s i n t e r a c t i o n with the other parts of the star are outlined. The equations of the charged f l u i d dynamics are derived, and order of magnitude estimates of the s i z e of various e f f e c t s are given. Magnetohydrostatic Equilibrium Of A Rotating Star A neutron star i s formed by collapse. Immediately a f t e r formation i t i s i n a very agitated state (LeBlanc and Wilson, 1970). Within roughly the time i t takes f o r a sound wave to cross the star i t should be o s c i l -l a t i n g about hydrostatic equilibrium. These o s c i l l a t i o n s are damped within seconds by the emission of g r a v i t a t i o n a l r a d i a t i o n and neutrinos. The star cannot be d i f f e r e n t i a l l y r o t a t i n g to any great extent, p a r t l y because of the i n h i b i t o r y e f f e c t of v i s c o s i t y , but more importantly because of the e f f e c t of the magnetic f i e l d . Any shear i n the charged f l u i d v e l o c i t y w i l l r a p i d l y "wind up" the magnetic f i e l d l i n e s and amplify the magnetic f i e l d within an e-folding time - l / ( v e l o c i t y shear). For a large shear, t h i s i s approximately the r o t a t i o n a l period of the s t a r . A f t e r several rotations, the magnetic f i e l d w i l l be magnified to such an extent that the magnetic forces are s u f f i c i e n t l y large to e f f e c t i v e l y r e s i s t d i f f e r e n t i a l r o t a t i o n . Shortly a f t e r d i f f e r e n t i a l r o t a t i o n has e f f e c t i v e l y stopped, the star should be near magnetohydrostatic equilibrium. The youngest pulsar known i s the Crab, which i s 921 years old. Even t h i s young pulsar should be close to r i g i d r o t a t i o n i n magnetohydrostatic equilibrium, assuming the absence of continual v i o l e n t e x c i t a t i o n . The 19 proper point of departure f o r an i n v e s t i g a t i o n of the i n t e r n a l motions of neutron stars which are not. extremely young i s therefore the problem of the magnetohydrostatic equilibrium of a r i g i d l y r o t a t i n g neutron sta r . This problem should best be treated within the framework of a r e l a t i v -i s t i c theory of g r a v i t a t i o n . There are at l e a s t two problems with any r e l a t i v i s t i c approach, however. The f i r s t i s one of p r i n c i p l e : which of the several v i a b l e r e l a t i v i s t i c theories of g r a v i t a t i o n should be used? The second i s a p r a c t i c a l problem. For consistency, the same g r a v i t a t i o n a l theory should be used i n the magnetohydrostatic c a l c u l a t i o n as w i l l be used i n the c a l c u l a t i o n of the f l u i d dynamics. The l a t t e r c a l c u l a t i o n requires a large number of g r a v i t a t i o n a l v a r i a b l e s (e.g., s i x for Einstein's theory) as compared to the one i n Newton's theory. I therefore decided to use Newton's theory of g r a v i t a t i o n throughout my c a l c u l a t i o n s . As a r e s u l t , such equilibrium quantities as the pressure and density w i l l be i n error by a maximum of about pressure/(density x c ), which i s < 1/3. The upper l i m i t of 1/3 occurs only at the centre of the most massive neutron s t a r s . For most neutron s t a r s , the error i s probably more l i k e 10%. The most important dynamical consequence of the non-r e l a t i v i s t i c approximation occurs i n the C o r i o l i s term i n the f l u i d dynamic equations. The neglect of the r e l a t i v i s t i c "dragging of i n e r t i a l frames" (Cohen, 1971) produces an error of about the same s i z e as i n the equilibrium equations. There i s one case i n which the n o n - r e l a t i v i s t i c approximation would produce q u a l i t a t i v e l y wrong r e s u l t s . Newtonian g r a v i t a t i o n a l theory i s not capable of p r e d i c t i n g some of the i n s t a b i l i t i e s of massive neutron s t a r s . Such c r i t i c a l l y massive neutron stars are not considered here. 20 The Equations of Magnetohydrostatic Equilibrium ?--*Po-Po*<'o - ? l V ? l ? > V ^ V x ( V - x B J The n o n - r e l a t i v i s t i c magnetohydrostatic equations are [II.1] [II.2] [II.3] o K o p„ = P„(p J . o o o In these equations, p Q i s the t o t a l equilibrium density, p Q i s the t o t a l equilibrium pressure, * Q - l s the equilibrium g r a v i t a t i o n a l p o t e n t i a l , • B i s the equilibrium magnetic f i e l d , and ^  i s the angular v e l o c i t y . . The temperature w i l l be denoted T Q . Three of the most v a r i a b l e or uncertain parameters i n both the e q u i l i b -rium problem and i n the problem of the charged f l u i d dynamics are ft , B and o o T . The angular speeds of known pulsars span roughly two decades, the magnetic f i e l d strength i s uncertain by at l e a s t two orders of magnitude, and the temperature of neutron stars of varying ages probably spans two or three decades. To f a c i l i t a t e making order of magnitude estimates, and to make manifest how ph y s i c a l quantities depend on ftQ, B q and T q , non-dimensional versions of these parameters are introduced i n Table I I I . A non-dimensional-ized length scale i s also introduced i n that table. Table, III Non-Dimensional Parameters T (°K)/10 7 11 '11 B q(gauss)/10 tiQ (rad/sec)/10" A, = A(cm)/10 - . 21 The value chosen to non-dimensionallze ftQ i s such that pulsar angular v e l o c i t i e s l i e i n the range .01 <v Q.^ ^  2. The surface magnetic f i e l d s for most pulsars as deduced from the dense magnetosphere model (Roberts and Sturrock, 1972) correspond to B.^ - 1. Values of T 7 in.the range .1 £ £ 10 are probably reasonable for neutron stars neither younger than the Crab nor extremely o l d . A length X, = 1 corresponds to 10 km, which i s a t y p i c a l radius f o r a neutron s t a r . For l a t e r reference, the values of a number of microscopic parameters of the charged f l u i d are c o l l e c t e d i n Table IV. . ' • Several approximations besides the n o n - r e l a t i v i s t i c one are already present i n [ I I . 1 - 3 ] . In [II.3] the temperature dependence of the pressure has been neglected because of the degeneracy of the s t e l l a r material - the Fermi energy of the neutrons i s about 200 Mev, so k T lz -rn B o Fn 5 x 10" 6 T ? « 1. Because of the presence of a magnetic f i e l d , the material pressure tensor should be anisotropic rather than i s o t r o p i c as i n [ I I . 1 ] , However, Ostriker and Gunn (1969) estimated that the e f f e c t of the pressure aniso-tropy on the s t e l l a r structure should be very much smaller than the e f f e c t of the J x B magnetic force, which i s i t s e l f extremely small as. i s shown below. Because of the small e f f e c t of the magnetic f i e l d , f l u c t u a t i o n s i n the permeability of the s t e l l a r material due to the de Haas -Van Alphen e f f e c t have been neglected. These f l u c t u a t i o n s tend to be thermally smeared out i n any case. In [ I I . l j , the s i z e of the c e n t r i f u g a l term r e l a t i v e to the p V$ term o o -4 2 i s = 1 0 which i s << 1. In making t h i s estimate, I have used $ q = GM/R 20 2 -2 33 - 10 cm sec , where M i s the mass of the star - 10 g and R i s the s t e l l a r radius = 10^ cm. The r e l a t i v e s i z e of the magnetic term i s (B O/4TTP O)/* O - 10 B^. I t i s therefore an excellent approximation 22 Table IV Some Microscopic Parameters Of The F l u i d Region Name Symbol Approximate Value 36 —3 Electron number density n 5 x 10 cm J e • f i i-\36 —3 Proton number density n 5 x 10 cm 3 P 38 -3 Neutron number density n 5 x 10 cm n Electron Fermi energy e F e 110 Mev Proton Fermi energy e-c 5 Mev *P Neutron Fermi energy e F n •200 Mev Electron Lorentz f a c t o r y 220 Elect r o n - e l e c t r o n -m~10 ^ . T 10 T-. sec sca t t e r i n g time e 7 -2 Electron mean free path X 3 T_ cm e 7 21 -1 Electron plasma frequency 0) 8 x 10 sec 21 -1 Ion plasma frequency t o ^ 3 x 10 sec 16 „ -1 pe p i sec Electron cyclotron i n - L b T> frequency ce 11 Ion cyclotron frequency B^^ sec ^  -13 Electron Debye length 3 x 10 cm 10 -2 2 -1 Charged f l u i d kinematic v 3 x 10 cm sec shear v i s c o s i t y c 25 -1 Thermal conductivity K 10 T ? erg/deg-cm-sec " ' 23 to completely neglect the effect of the magnetic field on the stellar structure and to work to only first order in the rotation. In this approximation, all variables can be expanded as the sum of a spherically symmetric term and a term with an angular dependence of P^Ccose). A spherical coordinate system (r,6,<f>). is used, with 9 = 0,TT as the rotation axis. The second term arises from the rotation, and is smaller than the 2 2 first by about flQR For example, the density may be expanded as [II.4] P (r,9) = P^0)(r) + p^ 2 )(r)P 9(cos6) . o o o 2. Under this expansion scheme, equations [II.1-3] reduce to a set of coupled ordinary differential equations. The Equation of State An equation of state must be chosen before these equations can be solved. There are by now a large number of equations of state for neutron st a r s . , 14 -3 All are in substantial agreement below densities of 2 x 10 gm cm , that is, in the crust whose properties are fairly well understood. Above that densityj and particularly near the centre of the star, the strong inter-particle interactions are not well understood. As physicists have begun to understand these interactions better, more and more refined equations of state have been devised. Unfortunately, this refining process has not yet converged - theoretically calculated neutron star structures continue to fluctuate by about the same amount as each new equation of state is proposed. In such a situation i t i s best to use the simplest equation of state which takes into account the features in which one is interested. The simplest equation of state which 14 -3 takes into account charged particles at densities greater than 2 x 10 gm cm 24 i s the Harrison-Wheeler equation of state (Harrison et a l . , 1965), so I chose i t for the work i n t h i s t h e s i s . The d e t a i l s of th i s equation of state i n the c r u s t a l region are not of i n t e r e s t f o r the purposes of t h i s t h e s i s ; howevery i t s behaviour i n the f l u i d region i s . The Harrison-Wheeler equation of state assumes that the f l u i d region i s composed of non-interacting degenerate neutrons, protons, and electrons. The r e l a t i v e composition at a given density i s fi x e d by the requirements of charge n e u t r a l i t y and 8-equilibrium. I t follows from charge n e u t r a l i t y that the average number de n s i t i e s of electrons and protons are equal, and thus t h e i r Fermi momenta are as w e l l : [II.5] m c sin h ( t /4) = m c sinh(t /4) , e e p p where m and t are the electron mass and dimensionless energy parameter, e e resp e c t i v e l y . Similar notation i s used f o r the protons and the neutrons. In 3-equilibrium the reactions [II.6] n —> p + e + v e [II.7] p + e —> n + v e proceed at equal rates. Since the neutrinos escape from the star , 8-. equilibrium implies that the chemical p o t e n t i a l of the neutrons i s equal to the sum of the chemical p o t e n t i a l s of the electrons and the protons: 2 2 2 [II.8] m c cosh(t /4) = m c cosh(t /4) + m c cosh(t /4) n n e e p p In [II.8], the r e s t mass energies are included i n the chemical p o t e n t i a l s . The t o t a l density i s the sum of the p a r t i a l d e n s i t i e s [II.9] p = D[(m /m ) 4 ( s i n h t - t ) + (m /m ) 4 ( s i n h t - t ) o e n e e p n p p + (sinh t - t )] , n n which can be written as the sum of the charged f l u i d density P c q and the uncharged f l u i d density PUQ' [11.10] p = p + p o co uo In [II.9], D i s the characteristic density: [11.11] D = 7rc 3m 4/4h 3 = 5.7 x 10 1 4 gm cm"3 . The pressure is the sum of the partial pressures: [11.12] p = P[(m /m ) 4(sinh t - 8 sinh(t /2) + "3t ) o e n e e e + (m /m ) 4(sinh t - 8 sinh(t /2) + 3t ) P n p e p + (sinh t - 8 sinh(t /2) + 3t ] , n n n which may be written as the sum of the charged f l u i d pressure .p . and the uncharged f l u i d pressure P U Q-[11.13] p = p + p . O C O uo In [11.12], P is the characteristic pressure: [11.14] P - | Dc 2 = 1.7 x 10 3 5 dynes cm"2 . Given any one of the three energy parameters, e.g. t^, the other two can be calculated using [II.5] and [II.8]. Then p (t ) , p (t ), p (t ), o p o p c o p p (t ), etc. can be calculated from [II.9] and [11.12]. This procedure co p was programmed on a computer to produce the equation of state P 0(P Q) f ° r use in equilibrium calculations and also to produce two other, functions " 2 useful in charged f l u i d dynamical calculations: p (p ) and s (p ). The " J C O O CO o latter function is defined by [11.15] s2nn = dp /dp . C O C O C O 26 S c q will be referred to as the charged fluid sound speed; i t appears as the propagation speed in the wave equation for sound-like waves in the charged fluid. An analogous quantity for the fluid as a whole is [11.16] s 2 = dp /dp o o o Strictly speaking, i t is not proper to refer to S q as the sound speed of the fluid interior because the neutron fluid and the charged fluid are weakly coupled and therefore do not oscillate together. The energy parameters in the density range 2 x 10"^ £ P q £ l O ^ vary over 20 "v t " V 25, .2 £ t ^ .5 and 1 £ t £ 2. Consequently in this density e p n range the electrons are very relativistic, the protons are non-relativistic, and the neutrons vary from non-relativistic to somewhat relativistic/ The main contribution to p arises from the neutrons, so [II.9] implies o that for t" « 1, n * [11.17] p ' = \ J D . o b n Similarly the main contribution to the charged fluid mass density arises from the protons, so [11.18] p * \ t H , c o o p and thus the fraction of charged fluid by mass over this density range is [11.19] p^/p. = (t It ) 3 - 10"2 . co o p n The exact numerical results for p /p as a function of p given in Figure • c o o o 2 bear out the rough estimate [11.19]. The main contribution to the total pressure arises from the neutrons, so for t a 1, n 27 For the charged fluid, the relativistic electrons provide most of the pressure, so after using [II.5] to eliminate t in favour of t , i t follows e p that • y [11.21] p '* ~ t 4 P . co 32 p It follows from [TI.17] and [11.20] that the equation of state is:very 5/3 nearly that of an ideal non-relativistic neutron gas, p « p , when t - 1. o o n When t is closer to two the gas is becoming relativistic, and the exponent should become smaller than 5/3. From [11.18] and [11.21], the charged fluid 4/3 equation of state p (p ) behaves like p - p . It also follows that co co co co [11.22] s 2 Is1 = 4t It1 - .5 . co o P n These approximate pressure-density relations are born out by the exact 2 2 numerical calculations. In Figure 3, s and s are plotted as functions o co 2 of p Q . The change in slope of s^ at high densities is due to the changeover from a non-relativistic to a slightly relativistic regime for the neutrons. 2 2 The s curve is straighter than the s curve because the protons remain co o 2 2 non-relativistic over almost the entire density range. Finally, s c 0 / s 0 * s about .55 and .60 throughout, in agreement with [11.22]. A Numerical Example The coupled differential equations which result from the expansion scheme [II.4] were solved on a computer.1 Two parameters define a particular model. These parameters are the central density p^°^(0) and the non-^ dimensional angular speed i ^ * Spherically symmetric functions such as p^°^ (r) do not depend on 9,^ (2) 2 in the approximation considered, and functions such as p Q (r) scale as Thus, from a particular solution obtained for given values of p^ 0^(0) and fl^* the family of solutions corresponding to this central density and any angular velocity may be obtained by scaling. The expected range of values for p^°^(0) i s not great. Most neutron star models of moderate mass have central densities which l i e within a 15 -3 factor of two either way of 10 gm cm In addition, there is a well-known scaling law (§109; Landau and Lifschitz, 1969) which allows one to generate a class of approximate solutions parametrized by p^°^(0) from a solution with a given p^°^(0). The scaling law states that p^°^(r)R^ i s a function only of r/R. It is exact in the absence of rotation and when the stellar material is entirely composed of a single non-relativistic degenerate gas. The latter condition was previously shown to be satisfied to a good approximation in the f l u i d interior. The value of p ^ ( 0 ) chosen for the numerical calculations was o 15 . ' - 3 10 gm cm ; ft- was taken to be 1.9, the angular speed of the Crab pulsar (see Table IT). The results for p^0"* (r) and p^ 2^ (r) are displayed in Figure 4. The scaling methods just mentioned may be used to generate a large class of approximate solutions from these results. 1 The routine STIFF from the U.B.C. Computing Centre was used. 29 The spherically symmetric part of the density drops smoothly from 10 gm cm at the centre to 2 x 10 gm cm at r "= R £ - 8.9 km, the radius of the crust. There are several kinks in the graphs of and for p^°^ < 2 x 10r** gm cm ^. These arise from changes in the form of the equation of state in the crust. The density drops to near zero at the stellar radius r = R - 15.8 km. The main features of this particular model are summarized in Table V. Table V The Parameters Of A Rotating Neutron Star Model Central density Angular speed Mass Radius Inner crustal radius Crust moment of inertia 10 1 5 gm cm 190 sec - 1 1.71 x IO 3 3 gm (about .9 M ) 15.8 km 8.9 km 4.06 x lO 1^ gm cm The Magnetic Field Structure The equations of magnetohydrostatic equilibrium place constraints on the form of the magnetic f i e l d . Taking the curl of 1/Pq times [II.1] and using the dependence of P q only on pQ produces [11.23] V x {[B x (V x B ) ] /4TP } = 0 , o o o which must be satisfied in addition to the Maxwell equation [11.24] i = 0 o It is an excellent approximation to neglect the dependence of p on S in o o [11.23], and to use the density as determined in a hydrostatic equilibrium calculation, since as previously shown, the magnetic f i e l d has only a minor effect on the stellar structure. ' 30 • There are a very large number of magnetic f i e l d structures consistent with [11.23].." One large class i s formed by magnetostatic f i e l d s which s a t i s f y $ x B Q = (5 ; a more general c l a s s i s formed by force-free f i e l d s , which s a t i s f y V x B Q .« B q. I t i s not presently possible to decide on the basis of theory alone which of these many solutions are a c t u a l l y r e a l i z e d i n neutron star i n t e r i o r s . It i s quite probable that the magnetic f i e l d i n the i n t e r i o r i s magnetohydrodynamically stable, considering the very long time a v a i l a b l e fo r i t to a t t a i n such a stable state. As was indicated i n Chapter I, the subject of magrietohydrodynamic s t a b i l i t y i n s i d e neutron stars i n not yet well developed.. Further consideration of t h i s question i s deferred u n t i l Chapter VI. However, some information about the i n t e r i o r magnetic f i e l d structure may be gleaned from the observed slowing down rate of pulsars. The angular v e l o c i t y of the crust i s a l t e r e d by torques exerted by the magnetic stresses at both i t s outer and inner surfaces. There i s no p a r t i c u l a r reason to suppose that the torque exerted on the inner surface i s very much smaller that that on the outer. I t i s thus reasonable to assume that the torque on the inner surface i s roughly equal to that e x t e r i o r torque required to produce the observed pulsar slowdown. From t h i s condition one may obtain as follows an estimate of B B at P t the inner surface, where B and B are the sizes of the p o l o i d a l and t o r o i d a l P t parts of B q, res p e c t i v e l y . The components of magnetic stress at the surface 2 r = R which can exert torques are = B B /4TT, S O the force i s ~ B B R , c M p t p t c' 3 . and the torque i s - B^B^R^ Equating t h i s torque to TO, where I i s the moment of i n e r t i a of the crust, implies that 31 [11.25] B B '* Ift/R 3 p t c at the surface of the crust. It i s not nec e s s a r i l y true that B p B t has roughly the value [ 1 1 . 2 5 ] throughout the f l u i d i n t e r i o r . However, there i s a large class of magnetic f i e l d s f o r which t h i s i s the case. For example, i f the magnetic f i e l d i s symmetric about any axis and the angular v a r i a t i o n s i n P q are n e g l i g i b l e , i t follows from [11.23] that even though B^ and B^ may separately vary by many orders of magnitude throughout the star, the product B p B t * s roughly constant. The values of I and R f o r the p a r t i c u l a r model given i n Table V can c be used to obtain a very crude estimate of B B at the inner c r u s t a l surface P t for any pulsar with known Q. For the Crab (see Table I I ) , the r e s u l t i s 1 ft 9 B B = 10 gauss . P t If one believes that B £ 1 0 ^ gauss, then i t follows from [ 1 1 . 2 5 ] that o for a l l known pulsars the magnetic f i e l d at the inner c r u s t a l surface i s either almost purely p o l o i d a l or purely t o r o i d a l . Charged F l u i d Equations Of Motion In t h i s section, the equations govering the motion of the charged f l u i d are derived. There are four sets of equations: the mass equation, the f l u i d momentum equations, the g r a v i t a t i o n a l source equation, and the magnetic induction equations. F i r s t l y , however, the range of length scales and time scales over which the charged p a r t i c l e s may be treated as a f l u i d are examined. The reason f o r wanting to treat the charged p a r t i c l e s as a f l u i d i s that the charged f l u i d equations are much simpler than the k i n e t i c equations. 32 A great deal of information about the system is lost in going to a fluid picture, however. For example, within a purely fluid picture, i t is not possible to treat collisionless damping of waves or any velocity space microinstabilities, both of which are often very important in- plasmas. The Range Of Validity Of A Fluid Description I assume that both the protons and neutrons are superfluid (Baym et al., 1969a), so that the most effective scattering mechanism in establishing local thermodynamic equilibrium is electron-electron scattering. The scat-tering time T -has-been estimated by Heinzman and Nitsch (1972): -~\I7 7 - 1 0 - ? [11.26] x = a / ( L I ) = 10 T, sec , e e Bo 7 where a is the fine structure constant and U £ is the chemical potential of the electrons, which is about 110 Mev (Ruderman, 1972). The value given by Flowers and Itoh (1975) for the viscosity n g of the electron fluid whose relativistic mass density is p g may be used in the approximate kinetic relation [11.27] x p > 3n Jo v* , e e e le to obtain an independent estimate of T G . Assuming that the electron Fermi velocity v„ - c, which is an excellent approximation, the result given by re [11.27] agrees with that of [11.26] to within a factor of two. Two necessary conditions for the validity of a fluid description of the charged particles are that the time scale T and the length scale A of any disturbances satisfy - i n -? [11.28] x » T « 10 T, sec e 7 and [11.29] X » A = V T , T - 3 T~2 cm . e Fe e 7 Here X Is the electron mean free path. The most extreme case occurs in e -2 very old pulsars where T^  may get as low as 10 , in which case X >> .3 km is necessary for a f l u i d description. A necessary condition for the validity of a single fl u i d description of the electrons and the protons, as opposed to a two-fluid description, is that [11.30] T » uf"!" ~ 3 x 10~ 2 2 sec , p i where to^ i s the ion (proton) plasma frequency, defined by (/ 2X1/2 4irn e \ .. E—1 = 3 x 10 sec" , mp / 36 —3 where n^ is the number density of protons, which is = 5 x 10 cm . Note that because the protons are non-relativistic, i t is a good approximation to use the proton rest mass in [11.31]. For the r e l a t i v i s t i c electrons, however, the plasma frequency i s (, 2\l/2 — * 8 x 10 sec" , 2 where Y i s the electron Lorentz factor, which is - e„ /m c - 220. 1 Fe e Another necessary condition is that length scales be much longer than electron or proton Debye lengths X^e and p^^ > e.g«> [11.33] X » X D e = vFea,-J = 3 x 1 0 _ 1 3 cm . 34 A comparison:of;[11.28] with [11.30] and [11.29] with [11.33] shown that if a fluid description is valid, then quasi-neutrality is an excellent approx-imation. The Mass Equation If the charged fluid is perturbed away from the magnetohydrostatic equilibrium, its density will in general change: [11.34] p. = pnn + p' . C CO c In what follows, only small departures from equilibrium are considered, so the density perturbation is relatively small: [11.35] |pVP c ol « 1 • The equation describing changes in is [11.36] 8o'/8t + V • (p cv c) = - P V T N , where v c is the charged fluid velocity. The term ~pV T N represents the creation or destruction of charged fluid by the weak interaction, which tends to create charged fluid from neutrons if p^_ < 0, and produces neutrons from charged fluid i f p^ > 0. As », the equation expresses the conser-vation of charged fluid mass. It is important to estimate because only in the limit of large x N can the charged fluid be treated as a distinct fluid with its own equation of state. If the charged fluid density is perturbed away from equilibrium, the electron and proton Fermi energies change by <Se-pe and ^ e^, respectively. For relatively small Fermi energy perturbations, the amount of energy E per unit volume which is available for radiation away by neutrinos is 35 Fe Fp . Fe The smallness of .|6e_ I compared to |6e„ I , which r e s u l t s from e_ << e„ , 1 Fp 1 r 1 Fe' Fp Fe was used to drop the second term on the right-hand side of [11.37]. The rate of the si n g l e neutron 3-decay re a c t i o n [11.38] n p + e + v . e insid e the f l u i d region of neutron stars i s very low, because momentum conservation requires the emission of electrons and protons with momenta small compared to e^ e and e F p » r e s p e c t i v e l y . The p r o b a b i l i t y of holes being a v a i l a b l e i n the electron and proton Fermi seas f o r these emitted - e F n / k B T o p a r t i c l e s i s = e , which i s << 1. Conservation of momentum i s e a s i l y s a t i s f i e d i n URCA reactions such as [11.39] n + n —> n + p + e + v g . The luminosity L a r i s i n g from such reactions was calculated by Bahcall and Wolf (1965) f o r the case i n which the a v a i l a b l e energy i s s o l e l y thermal. 8 The luminosity f o r that, case i s « (k BT) . However, when the a v a i l a b l e B energy can arise.both from thermal e f f e c t s and from density induced changes 6e F i n Fermi energies, the expression for L contains a serie s of terms i n (kgT) 8, ( k B T ) 6 ( 6 e F ) 2 , (kgT)^\6e^)^, etc (Pethick, private communication). Since the 6e F's considered here are << kgT, the most important contribution to L which includes the e f f e c t s of small density perturbations i s 6 2 (kgT) (6 E f) . One may then take Bahcall and Wolf's (1965) expression 8 6 2 and replace (kgT) by (kgT) ( 5 e F e ) to get a f a i r l y good approximation to the luminosity a r i s i n g from small charged f l u i d density perturbations. The r e s u l t i s [11.40] L = I O 2 3 T^ ( 6 e F e ) 2 Mev cm"3 y e a r " 1 , , 36 if 6e F e is measured in Mev. In [11.40], a term « (^e^) is neglected, because f 6e_ I .« |6e„ I. 1 Fp1 ' Fe1 The decay time is thus [11.41] x„ = f * 10 1 2 T i 6 years. N L 7 In the calculation of Bahcall and Wolf (1965) i t was assumed that the neutrons are not superfluid. If superfluidity is taken.into account, the decay time 12 " —6 of 10 T^  years should not be shortened, so [11.40] can be looked upon as a lower limit for T ^ . The effects of charged fluid creation or destruction are thus completely negligible over the time scales of interest, and mass conservation holds: [11.42] + V • ( p v ) = . d t C C The Momentum Equations The equations describing the acceleration of the charged fluid have the form UP v ) • [11.43] — — .+ V • ( p v c v c ) = (sum of a l l forces/unit volume 0 on the charged fluid}. The various forces are given in Table VI. In that table, the terms which support the charged fluid when i t is in magnetohydrostatic equilibrium with the uncharged fluid have been subtracted off, although no linear-izations have yet been made. The frame of reference which is used co-rotates with the solid crust at angular velocity Q: [11.44] ft(t) = ^+^'(t) . 37 Table VI Forces Acting On The Charged F l u i d C entrifugal P c v ( i | S x ? | 2 , -,Jik ft x o r| ) C o r i o l i s - p 2 ^ x v c c Poincare - p fi x r Gr a v i t a t i o n a l - (p V$ - p v V ) C CO o Pressure - V • p' Magnetic ( 1 / 4 * ) $ • [B B - i l -o o t < * 2 2 i - > - , Viscous - V • a' Superfluid drag P„(v - v )/x c s c D The small v a r i a t i o n s ^' i n are due mainly to electromagnetic torques on the crust which cause i t s general slowdown. In a l l c a l c u l a t i o n s i n th i s t h e s i s , the z axis i s taken to l i e p a r a l l e l to ft . i o The various terms l i s t e d i n Table VI are now discussed. Each w i l l be l i n e a r i z e d i n the sense that terms which are second order i n perturbed quantities w i l l be dropped. This l i n e a r i z a t i o n i s an excellent approximation for the small amplitude perturbations that are considered i n t h i s t h e s i s . A discussion of some of the expected e f f e c t s of the dropped non-linear terms i s given i n Chapter I I I . Centrifugal force To f i r s t order i n primed v a r i a b l e s , the c e n t r i f u g a l force can be written as p'V(-^-|ft x r l 2 ) + p V(-|-|^ x r | 2 - -rift x r | 2 ) . The second term i s conser-c .2 o 1 co 2 2' o ' vative and can be balanced by pressure gradients within roughly the time i t - 4 takes a sound wave to cross the s t a r , about 1 0 sec. Over time scales •••• " 38 shorter than t h i s , fi cannot change appreciably. Thus no matter what the time scale, the second term may be neglected. The f i r s t term i s discussed l a t e r along with the g r a v i t a t i o n a l term. 2. C o r i o l i s force The C o r i o l i s force i s extremely important i n the charged f l u i d dynamics, A f u l l e r discussion of i t s importance i s delayed u n t i l Chapter I I I , which deals with the types of waves which can propagate i n the charged f l u i d . The l i n e a r i z e d f orm of the C o r i o l i s term i s —2fi x p v . The non-o co c l i n e a r convective term V • (P V cv" r) In the f l u i d momentum equation [11.43] i s smaller that the C o r i o l i s term by about a f a c t o r V /fi'A, where the order c o of magnitude of v^ i s denoted V , and A represents the length scale over -*• which v c v a r i e s . The condition that the f l u i d i s close to r i g i d r o t a t i o n , [11.45] V /fi R « 1 , c o i s equivalent to the smallness of the convective term compared to the C o r i o l i s term f o r long wavelength (A = R) disturbances. 3 . Poincare" force The P o i n c a r i f o r c e 1 i s a f i c t i t i o u s force which arises from the a c c e l -eration of the coordinate frame attached to the crust. If the crust i s precessing the o s c i l l a t i n g Poincare force which r e s u l t s can excite the charged f l u i d at the precession frequency. If the crust i s not precessing but only gradually slowing down, the r e s u l t i n g Poincare for c e i s almost constant i n time. In such a s i t u a t i o n , the f l u i d has ample time i n which to balance t h i s f o r c e . Because the force 1 The force was so named by Malkus (1970), who considered the p o s s i b i l i t y that i t plays an important r o l e i n the dynamics of the f l u i d core of the earth. i s non-conservative, i t cannot be balanced by a pressure gradient. However i t can be balanced by the magnetic body forces. This w i l l be shown' l a t e r i n t h i s chapter. 4. G r a v i t a t i o n a l force To f i r s t order i n primed q u a n t i t i e s , the g r a v i t a t i o n a l f o r c e consists of two terms: - p % Q and -PJ$<, w h e r e the g r a v i t a t i o n a l p o t e n t i a l pertur-bation $' i s defined by [ 1 1 . 4 6 ] $ = $ + $ o The f i r s t term represents the buoyancy force on an element of charged f l u i d , and may be conveniently combined with the unbalanced p o r t i o n of the c e n t r i f u g a l force to produce -p'V($ - •^ •1^  x r | 2 ) . c o 2 o The second term represents the g r a v i t a t i o n a l force on the charged f l u i d due to perturbations i n the g r a v i t a t i o n a l p o t e n t i a l . These perturbations a r i s e from changes i n pc, whereas $ q a r i s e s from both the charged f l u i d and su p e r f l u i d . The second term i s thus about p /p - 1% the s i z e of the f i r s t co o It i s therefore a f a i r l y good approximation to neglect -p vV. co 5. The pressure term In the absence of a magnetic f i e l d , the perturbation*!?' i n the charged c f l u i d pressure tensor would be i s o t r o p i c , and -V could be written as c -Vp c > In the presence of a magnetic f i e l d , P^ can be a n i s o t r o p i c . I t i s important to determine whether or not the anisotropy can be neglected. On the c l a s s i c a l l e v e l the conditions f or a nearly i s o t r o p i c pressure tensor are ( K r a l l and T r i v e l p i e c e , 1973) [11.47] ( a . A ) 2 « 1/(TU> c 1) « 1 , • 4 0 where a i s t h e i o n ( p r o t o n ) c y c l o t r o n r a d i u s , to ' . - " i s t h e i o n c y c l o t r o n i c i • J f r e q u e n c y , and X and T a r e m a c r o s c o p i c l e n g t h and t i m e s c a l e s , r e s p e c t i v e l y , to ^ i s g i v e n by [ 11 .48 ] tu , = eB /m c = 1 0 1 5 B., s e c " 1 . c i o p 11 N o t e t h a t a n o n - r e l a t i v i s t i c f o r m s u c h as [ 11 .48 ] f o r i s a good a p p r o x i -m a t i o n . F o r l a t e r r e f e r e n c e , t h e e l e c t r o n c y c l o t r o n f r e q u e n c y i s [ 11 .49 ] oi = eB /ym c - 1 0 1 6 B,, s e c " 1 . ce o e 11 A c o m p a r i s o n o f to 1 w i t h t h e e l e c t r o n - e l e c t r o n s c a t t e r i n g t i m e x i n [11.261 c i ° e r e v e a l s t h a t i n t h e f l u i d r e g i m e where x >> x , t h e c o n d i t i o n l / ( x u .) << 1 e c i i s obeyed . By q u a s i - n e u t r a l i t y , t h e F e r m i momenta o f t h e e l e c t r o n s and protons are 9 - 1 e q u a l , so v_,. .- e_ /cm = 6 x 10 cm sec , and t hu s F i F e p [ 11 .50 ] a , = v_.w~] " = 6 x 1 0 ~ 6 B7* cm . l ( l c i l l A c o m p a r i s o n o f [ 11 .29 ] and [11 .50 ] shows t h a t when X >> X e , ( a i / X ) << 1. 2 The answer t o t h e q u e s t i o n o f w h e t h e r o r n o t (aj\) i s << l / ( x t o c i ) depends on t h e r e l a t i o n s h i p o f x and X, I t i s shown i n C h a p t e r I I I t h a t t h e r e a r e t h r e e c l a s s e s o f waves w h i c h c an p r o p a g a t e i n t h e c h a r g e d f l u i d : sound wave s , i n e r t i a l waves and hydromagnet i c -^ i n e r t i a l wave s . F o r e ach c l a s s o f wave , t h e r e i s a d i f f e r e n t r e l a t i o n be tween x and X. F o r sound 2 wave s , t h e r e l a t i o n be tween x and X i s x - X/s , so ( a ,/X ) << l/ (xto ..) co i c i 2 - 5 - 1 i f X >> a.to ,/s = 1 0 B,, cm, w h i c h i s s a t i s f i e d i n t h e f l u i d r e g i m e , i c i c o 11 ° -1 1/2 F o r i n e r t i a l wave s , x - Q i n d e p e n d e n t l y o f X, so i f X >> (to ./ft ) a . = o r c i o x 1/2 2 200(B i ; L/ft 2) cm t h e n (aJX) « l / ( x t o c i ) . F o r h y d r o m a g n e t i c - i n e r t i a l — 6 **2 2 2 waves , x = 10 ^ 2 B 1 1 X ' s o ( a ± / X ) K < 1 ^ T ( J J C i ^ f o r a 1 1 w a v e l e n g t h s i f - 1 3 - 2 n2*ii > > 4 x 1 0 ' "'41. On the quantum level there is an additional restriction, which is that the electron Landau levels are thermally smeared out: [11.51] k„T » -no) . 1 5 o ce -1 -2 which can be written as T^B^ » 10 In what follows, i t is assumed that these conditions are met so that the pressure tensor can be treated as isotropic. The pressure perturbation p^ can arise from both density and temperature perturbations: U 1 - 5 2 1 • However, the charged fluid is highly degenerate, and in the limit of extreme degeneracy (kgT /e p e —> 0), both factors (3p /8T) and T' in the second 0 6 0 pc term of [11.52] separately vanish. In this limit, p^ becomes independent of T and the heat transfer equation becomes decoupled from the other equations of motion. This point is discussed more quantitatively later. For the moment we shall ignore the small thermal corrections to the charged fluid equation of state p = p (p ) calculated in the first section of this c c c chapter, so to first order in primed quantities [11.53] p' = s 2 p' . C CO c Equation [II.53] can be used to eliminate p^  from the momentum equations. Using [11.53] and the magnetohydrostatic equation [11*1] with the negligible magnetic term dropped, the buoyancy term can be written as [11.54] -P5($ - x r| 2) = - p ' ^ , C O z'o ' C O where the dimensionless quantity ^ q which depends only on the hydrostatic solution is defined by 42 i s n n [11.55] The l i m i t s o f i n t e g r a t i o n on ^ a r e a r b i t r a r y , because ^ appea r s o n l y as a g r a d i e n t , c o l l e c t i n g the buoyancy and p r e s s u r e g r a d i e n t t e r m s , we: have [ 1 1 . 5 6 , ::.^^"-p;»(. 0-i |S o x i \ h - - ( * + W o c 6. The m a g n e t i c body f o r c e s To f i r s t o r d e r i n p r i m e d q u a n t i t i e s , t h e m a g n e t i c body f o r c e s a r e - ( l / 4 i r ) B o x (7 x B ' ) - ( 1 / 4 T T ) $ ' X (V" x B )', where the magne t i c p e r t u r b a t i o n i s d e f i n e d by [11.57] B - B + P . 0 As i n t he m a g n e t o h y d r o s t a t i c e q u a t i o n s [11.23] f o r B , any de Haas -Van A l p h e n o o s c i l l a t i o n s i n t he p e r m e a b i l i t y a r e n e g l e c t e d because o f the t h e r m a l smear ing [11.51] o f e l e c t r o n Landau l e v e l s . As p r e v i o u s l y m e n t i o n e d , t h e m a g n e t i c body f o r c e s c a n b a l a n c e the q u a s i - s t e a d y Po incareV f o r c e due t o c r u s t a l s l o w i n g down. The c o n d i t i o n f o r such a f o r c e b a l a n c e i s [11.58] - ~[S" x(v" x P) + B' x (V x f ) ] - p^^n x r - Vp/ - p ^ = fj C C O From [ 1 1 . 5 8 ] , t he r e l a t i v e s i z e o f m a g n e t i c p e r t u r b a t i o n r e q u i r e d to do t h i s i s [11.59] B ' /B . p m 2 / ( B 2 / H ) , O CO o - 5 - 2 w h i c h i s about 10 B n f o r t he C r a b . The m a g n e t i c f i e l d has a v e r y l o n g t ime i n wh i ch to a d j u s t t o t he Po incare" f o r c e due t o c r u s t a l s l o w i n g down. ' 43 Thus i t seems reasonable that given any B q satisfying [11.23], a small B ' will be present such that [11.58] is satisfied. However, this is not possible for every magnetohydrostatic For example, in the axially symmetric case, the (J) component of [11.58] is [11.60] ( B • V ) ( r sine B ' ) - p- tiv2 sin 29 . o <p CO The right-hand side of [11.60] is </0, so i f B q has any closed field lines within the fluid region, [11.60] has no single-valued solution for B ^ . Magnetic fields B of this closed field line type and perhaps other types o as well cannot balance the quasi-steady Poincare force by making only small changes B * ' in their structures, and are thus in this sense unstable. ->• '-->•'•*> A simple example of a stable B is a uniform field B = B z, where r o o p V B Q - O " . The B ^ required to balance the Poincare force is displayed in Figure 5. For this calculation, B . ^ = 1, ti is that of the Crab (see Table 15 —3 ' II) and the total central density is 10 gm cm (the model of Figure' 4)• • ' This result for B ' can be easily scaled for other values of B . . , ti and <p II P<°><0>. For the remainder of this thesis, only those B * 's which are stable in o the sense just explained are considered. It is then convenient to absorb the B ' required to balance the quasi-steady Poincare" force into the definition of B q and to drop the quasi-steady Poincare force from the equations. The resulting equations are valid only over time scales short compared to the time over which ti varies. For the Crab this time is roughly ti/ti - a thousand years (see Table II). Any Poincarg force due to precession with a period « ti/ti should be retained in the equations as a term of the form -p ti x r if i t is not small compared to other retained terms, co p However, precession will not be considered in this thesis. 44 . The Viscous term The charged fluid experiences viscous effects due to electron-electron scattering. The charged fluid shear viscosity n in the absence of a magnetic field has been calculated by Heinzman and Nitsch (1972). Flowers and Itoh (1975) calculated the electron fluid shear viscosity n e which is related to n by n lp = n lp . The two results for n are in good agreement. The charged fluid kinematic shear viscosity V £ is defined as follows: 1 9 1D -9 9 _1 [11.61] v - n lp «.±r c = 3 x 10 Tn cm sec . c c c 3 e 7 -2 Note that the shear viscosity is a T . On the other hand, the Landau "theory 2 of Fermi liquids predicts that the bulk viscosity is « T (Abrikosov and Khalatnikov, 1958), which vanishes in the extreme degeneracy limit. Accord-ingly, the bulk viscosity is neglected. Since V £ depends on T, which fluctuates slightly due to the fluid flow, v c experiences slight fluctuations. The effect of these viscous fluctu-ations on the flow is negligible i f the temperature fluctuations are small. The size of the temperature fluctuations is discussed in Chapter III. The estimate [11.61] for is based on the assumption that the magnetic field vanishes. This is a good assumption if the electron cyclotron fre-quency io is much smaller than the electron-electron collision frequency T" 1. From [11.26] and [11.49] i t follows that [11.62] oo T = 106 B-.T"2 » 1 , ce e 11 7 so i t is a poor approximation to neglect the effect of the magnetic field on the viscosity. This was noted by Heinzman et al.(1973). 45 In the limit <«> T °°> which applies to the charged fluid in view of.[11.62], the isotropic viscosity tensor with non-zero components - 'n for the B = 0* case is replaced by an anisotropic tensor. This tensor depends on the direction of B but not on its magnitude. Its non-zero components are s t i l l = n c (Chapman and Cowling, 1939). Thus, the order of magnitude 2 of the viscosity term is equal to n V /X in both the u' x ~ > 0 limit and the OJ x —> °° limit. The ratio of the size of the viscosity term to the ce e J size of the Coriolis term is to within a factor of two the Ekman number E: [11.63] E = v /Q X 2 3 x io-4 n~ 1T" 2x~ 2 c o 2 / 0 where the dimensionless length scale X, is used (see Table III). Since o E is very small for X, - 1 and T., ^  .1, the viscosity should have only a o 7 small effect on large scale charged fluid flow within pulsars which are not extremely did. In such cases where the viscosity has l i t t l e effect the anisotropy in the viscosity will not have an Important effect on the motion. The superfluid drag Electrons in the charged fluid can scatter off the normal neutrons in the superfluid vortex cores, as can the small fraction of protons which are not superconducting. As a result, there is a drag force which tends to reduce the relative velocity v - v between the neutron superfluid and the charged fluid. For small relative, velocities, this drag force may be represented by a term P c ( v g ~ V c ) / T n t^ i e charged fluid momentum equation. Feibelman (1971) estimated the relaxation time required to reduce the relative velocity between the electrons and the neutrons, and showed that 46 the electron drag is larger than the proton drag. Baym et al.(1975) discuss the relationship of this relaxation time with the pqst-glitch relaxation time of pulsars. For reasonable values of the adjustable parameters, two of which are the temperature and the neutron energy gap size A, the theore-tically predicted post-glitch relaxation times can be the order of days (as in the Crab) or a year (as in the Vela). The fluid momentum equation has been discussed term by term, and a number of approximations have been made and justified. Given below are the resulting linearized equations: [11.64] |-(p v ) = -2ft xp v - Vp' - p'vV 3t co c o co c "c c o; .-p v*' - (1/4TT)B x(V, x B ' ) - (1/4TT)B' x(V x B ) CO o o -v" •"'"a' + p (v - v ) / - r n . co s c D The Gravitational Source Equation The gravitational potential perturbation satisfies [11.65] V2*' = 4TTGP' . c This equation is coupled to the momentum equation through the term -pcQV$', which as previously noted is only about 1% of the size of the buoyancy term for long wavelength motions. In that case", i t is a fairly good approximation to neglect the term -p^V*' and not to solve [11.65]. This is done in a l l the free oscillation calculations in this thesis. 47 • - ^  v The Magnetic Induction Equations The magnetic field evolves according to [11.66] I? = -cv" x E ° t subject to the constraint [11.67] V • B = 0 . For processes happening so slowly that electromagnetic radiation can be neglected, the current J is given by the pre-Maxwell equation -h -> A-jr n T» ' 1 [11.68] V x B = - J c To close these equations, a relationship must be found between the electric field E and the other variables in the problem, such as and p' , i.e., a generalized Ohm's law appropriate to the charged fluid must be found. A standard procedure (Krall and Trivelpiece, 1973) used in forming a generalized Ohm's law begins with the momentum equations for the electron and proton fluids. Not a l l terms in these equations will be written out; those that are omitted are analogous to terms in the one-fluid equations which are small. A term that appears in the two-fluid equations which has no analogue in the one-fluid equations is the electron-proton collision term. The resulting linearized two-fluid equations are 3v [11.69] n m - ~ = -2ft x n m v - Vp' - n'm V$ oa a 9t o o a a a a a a o v ^ 6v + n q ( E ' + - ^ x B ) + n m - ^ coll. pa a c o oa a fi.t In [11.69], a = {e,p} labels electrons and protons. Subscripts 'o' refer to unperturbed quantities and primes refer to perturbations, as in the 4 8 o n e - f l u i d e q u a t i o n s . The symbols n , v , q and m r e f e r t o t he number d e n s i t y , v e l o c i t y , cha r ge and e f f e c t i v e mass, r e s p e c t i v e l y , o f p a r t i c l e s o f t ype a . Because t h e e l e c t r o n s a r e r e l a t i v i s t i c and the p r o t o n s n o n -* * * r e l a t i v i s t i c , m = ym and m - m . I n terms o f t h e s e v a r i a b l e s t h e p e r t u r -e e p p * b a t i o n s i n J , i n t h e cha r ge d e n s i t y and i n the c e n t r e o f mass v e l o c i t y a r e d e f i n e d by [11.70] J * =\q n v / ^ot oct a a [11.71] p ' = ^ n ' q q /_ a a [11.72] v = ^ \ n m c p / ~"oct~q V a c a * ft M u l t i p l y i n g [11 .69] by <l a /m a and summing o v e r a p r o d u c e s an e q u a t i o n f o r t h e t i m e e v o l u t i o n o f 3 ' . T h i s e q u a t i o n becomes c o n s i d e r a b l y s i m p l e r i f terms o f o r d e r m g / m p - .1 a r e n e g l e c t e d compared t o one . I t w i l l be shown t h a t t h i s a p p r o x i m a t i o n does n o t p r o d u c e comparab le e r r o r s i n the g e n e r a l i z e d Ohm's l aw. The e q u a t i o n f o r dJ'/3t i s [I.I.73] H = x J + \ V p . ' C M - P 'V$ at o * e q o m n e 2 + + ^ (P + f x l ) - J' x.f. m* c ° cm* 0 6 t s e c o l l . y CM The term Vp^ i s t he g r a d i e n t o f the p e r t u r b a t i o n i n the e l e c t r o n p r e s s u r e as v i e w e d i n t h e c e n t r e - o f - m a s s f r a m e . S i n c e a s m a l l r e l a t i v e v e l o c i t y between the e l e c t r o n and p r o t o n f l u i d s p r o d u c e s an e x t r e m e l y l a r g e c u r r e n t wh i ch tends t o r e d u c e the r e l a t i v e v e l o c i t y , v - v = v and thus V p ' e p c e e q u a l s V p e t o a f a i r a p p r o x i m a t i o n . W i t h the s t a n d a r d a p p r o x i m a t i o n t h a t the c o l l i s i o n term ( 6 ^ ' / 6 t ) c o l l i s « J ' , [11.73] may be r e w r i t t e n as t he g e n e r a l i z e d Ohm's law 4 9 V [11.74] E' = - x B + {- + ^ [^ + (2ft - co B ) x ]}P c o a 2 9t o ce o pe • i . 4irp1 en re 2 o e to pe 21 - i where w is the electron plasma frequency, which is = 10 sec , B Q is a unit vector in the direction of B , and a is the conductivity in the o absence of a magnetic field. The -v x B* term is the induced emf. If i t is the only significant c o term in the Ohm's law, the magnetic field evolves according to the infinite conductivity induction equation [11.85]. There are four terms linear in J'. One is the usual Ohm term J'/o. Another is the Hall term « B x J ' which o represents the magnetically induced anisotropy in the electrical conduction. Another anisotropic term « ft x J ' arises from the rotation of the star. o There is a term a 9?/9t which represents the electric field set up by very rapidly fluctuating currents. Electric fields are also set up by gradients in the electron pressure and by the action of gravity on fluid regions where the perturbation in the charge density does not vanish. Because of quasi-neutrality one might expect the term a P^®^ to be -> -v small. This is indeed the case. Using V • E' = 4TTP^ to estimate Pq ~ E'/4irA, i t follows that 2 U i r ' p ' v * / o i I cm I'I [H.75] V P 6' - IO" 2 8 A"1 « 1 if A » IO' 2 8 Such short wavelengths are, of course, unphysical. It will be shown in Chapter III that sound waves are the highest frequency waves that the charged fluid can support. Because they have the highest frequency, the pressure perturbation in them is higher for a given Vc than for other kinds of waves. Thus |Vp'| £ p V /x where the sound wave period x = A/s . c co c co Since |pj % |p;|,V 50 Vp'/en [II.76] '• — — 2_ < co co B -15 -1 x - l i - • en XB D H • e o 1 * — v x B 1 c c o which is « 1 if X » 10~ 1 5 B~* cm. The four terms linear in 5'stand in the ratio - : •—- — 4ir 4TT a V 2 x " 2 2 V OJ OJ pe pe 0) Using the conductivity a given by Baym et al.(1969b), and the oj P e fact that the shortest period oscillations are sound waves, the ratios are roughly 10~ 3 1 T 2 : 10 - 3 7 X ^ 1 : 10 3 9 &2 : 10~ 2 5 Bu.!-'The Hall term is dominant except for some extreme conditions such as"a very low magnetic field or a very high temperature. To estimate the size of the Hall term one needs an estimate of J * . From [11.68] [11.77] J' ~ cB'/4irX . • • •' . • ' ' i ' -| ' An assumption about the origin of the perturbations B must be made in order to estimate B'. The assumption I make is that if" v c = fi, then B' = fi, i.e., any perturbations in B originate in the action of charged fluid motion on B^ . That being the case, i t is reasonable to suppose that the induced emf is not very much smaller than the Hall term. Equation [11.66] can then be used to estimate B1 in terms of v : c [11.78] B' =VB T / X . c o The relative size of the Hall term is thus | (4lT0J /OJ ) B X J ' l OJ C T c o H I . , , , ^ Pe . _ c ^ _ . 1 0-5 2 , _ . — v x B OJ X I c e o pe For sound waves, T - X/s . For such waves the Hall term is negligible if co [11.80] X » 10~ 1 5 B n cm (sound). 51 F o r i n e r t i a l waves , x - ft i n d e p e n d e n t l y o f A. F o r them the c o n d i t i o n i s [11.81] A » 1 0 " 3 ( B 1 1 / f t 2 ) l / 2 cm ( i n e r t i a l ) . 2 — 6 - 2 F o r h y d r o m a g n e t i c - i n e r t i a l waves , x/A - ' 1 0 n 2 B l I ' ' 8 0 o v e r t h e t i m e s c a l e s o f t h e s e waves, the H a l l te rm i s n e g l i g i b l e i f [11.82] Q 2 B 1 1 < K ( h y d r o m a g n e t i c - i n e r t i a l ) . The c o n d i t i o n s [ 11 .80 -82 ] a r e e a s i l y s a t i s f i e d when t h e cha r ged p a r t i c l e s can be t r e a t e d as a f l u i d ( see e q u a t i o n s [ 1 1 . 2 8 - 2 9 ] ) . Under such c o n d i t i o n s , t he l i n e a r i z e d f o rm of t he g e n e r a l i z e d Ohm's law r e d u c e s to v [11.83] , E ' + — x S = 0" , C O w h i c h g i v e s the l i n e a r i z e d i n d u c t i o n e q u a t i o n [11.84] = v" x (v x $ ) , dt C O of w h i c h the n o n - l i n e a r f o rm i s [11.85] ff: = v" x ( v c x B) . I t f o l l o w s f r om [11.84] t h a t a n e c e s s a r y c o n d i t i o n f o r the v a l i d i t y of a l i n e a r i z e d t r e a t m e n t i s [11.86] V « A/x . c The c o n c l u s i o n f r om t h i s a n a l y s i s i s t h a t f o r a w ide v a r i e t y o f A and x wh ich encompass t he e n t i r e f l u i d r e g i m e , the m a g n e t i c f i e l d e v o l v e s a c c o r d i n g to the i n f i n i t e c o n d u c t i v i t y i n d u c t i o n e q u a t i o n . 52 The Boundary Conditions The charged fluid is enclosed within two surfaces which are roughly concentric spheres (or within one if there is no solid core). The tran-sition between crust and fluid is a sudden one (Baym et al. , 1971). The transition from the fluid to any solid central core probably also involves a phase change (Pandharipande and Smith, 1975). It is necessary to specify the behaviour of the charged fluid velocity and the magnetic field at these boundaries in order to be able to calculate the fluid motion. The Velocity Boundary Conditions Since the charged fluid has a non-zero viscosity, the appropriate velocity boundary conditions at a solid surface Z are "no-slip" conditions, that is, [11.87] v| E = fi . It is assumed in [11.87] that the frame of reference is chosen so that the surface is not moving. That is not the case with the core. For i t , the boundary condition is [11.88] vj = Afl x r , where Afi is the difference in angular velocity between crust and core. Although [11.87-88] are the correct velocity boundary conditions, i t is not clear that they are an adequate and consistent set of boundary conditions when one formally lets the viscosity tend to zero, which may be a desirable approximation in certain circumstances. • 53 Evan i n t h e i n v i s c i d l i m i t , t he s u r f a c e i s impermeable to the cha r ged f l u i d , so t he component o f v e l o c i t y no rma l t o i t i s z e r o : [11.89] v • nlj, = % , where n i s a u n i t v e c t o r no rma l t o L I t f o l l o w s f r om M a x w e l l ' s e q u a t i o n s that the t a n g e n t i a l component o f I i s c o n t i n u o u s a t the i n t e r f a c e . J u s t i n s i d e the c r u s t E v a n i s h e s i n the i n f i n i t e c o n d u c t i v i t y a p p r o x i m a t i o n . Thus , [11.90] n x I|E' - - J .' . Ohm's law [11.83] and [11 .89 -90] i m p l y t h a t [11.91] (B • n ) ( n x v ) | E 4 , so t h a t [11.92] v| E = ^ i f B • n | E 1 u" . Thu s , i f t h e r e i s a n o n - z e r o component o f m a g n e t i c f i e l d no rma l to Z a t a p o i n t , t hen the f l u i d s a t i s f i e s n o - s l i p boundary c o n d i t i o n s t h e r e i even i n t h e i n v i s c i d l i m i t . In t he d e r i v a t i o n o f the boundary c o n d i t i o n [ 1 1 . 9 1 ] , i t was i m p l i c i t l y assumed t h a t the cha r ged f l u i d i s s t r o n g l y c o u p l e d to the magne t i c f i e l d . I t w i l l be shown i n C h a p t e r I I I t h a t ove r t ime s c a l e s ^ ft ^ , t he two a r e o n l y v e r y weak ly c o u p l e d . Under such c o n d i t i o n s , i t i s a v e r y good a p p r o x i m a t i o n t o c o m p l e t e l y n e g l e c t t h e e f f e c t o f t he m a g n e t i c f i e l d on the f l u i d f l o w , and to u se the i n v i s c i d f l u i d e q u a t i o n s w i t h no m a g n e t i c body f o r c e . t e r m . I f such a m a t h e m a t i c a l a p p r o x i m a t i o n i s made, o n l y one boundary c o n d i t i o n can be imposed on a t Z. T ha t c o n d i t i o n i s the " s l i p " boundary c o n d i t i o n [ 1 1 . 8 9 ] . 5 4 The Magnetic Field Boundary Conditions The penetration depth 6 of variable magnetic fields with frequency to into the crust is (§45; Landau and Llfshitz, 1966) [11.93] 6 = c/^ooj , where a is the crust conductivity. For the lowest frequency waves that can propagate in the charged fluid (see Chapter III), 6 is tens of cm, which is considerably smaller than the thickness of the crust. It is thus an excel-lent approximation to assume that the solid parts of the star are infinitely conducting. Therefore the magnetic field in these regions cannot decay, and since we are assuming these regions are also perfectly rigid, the magnetic field is static there, as seen in the crustal frame. From Maxwell's equations, B* • n is continuous across a boundary Z. However, on the solid side of E, B • n is static, and thus B* • n on the fluid side is also. In this sense, the magnetic field lines are "pinned" at the boundaries. However, in the infinite conductivity limit the tangen-tia l component of B, n x B, is not necessarily continuous across I. What actually happens Is that if the magnetic field lines in the fluid are oscillating with frequency u, n x B* can vary a large amount in the small distance 6(to) given in [11.93]. 55 The C r u e t And C o r e E q u a t i o n s The s o l i d c r u s t and the s o l i d c o r e ( i f one e x i s t s ) a r e v e r y r i g i d . In t h i s t h e s i s , t h e y w i l l be t r e a t e d as p e r f e c t l y r i g i d . A n e c e s s a r y c o n d i t i o n f o r the a p p r o x i m a t i o n to be a good one i s t h a t t he s t r e s s on a s o l i d p a r t o f the s t a r i s s m a l l e r than i t s m a t e r i a l m o d u l i t imes the maximum s t r a i n i t can t ake w i t h o u t f r a c t u r i n g . F o r examp le , i f t h e i n n e r c r u s t can t a k e a - 2 s t r a i n o f 1 0 . , t h e maximum s t r e s s t h a t i t c o u l d t ake w i t h o u t f r a c t u r i n g 2 7 2 9 - 2 would be about 1 0 t o 1 0 dynes cm I f the cha rged f l u i d i s i n s h e a r i n g m o t i o n , then t h e r e a r e v i s c o u s s t r e s s e s s e t up i n i t . These shea r s t r e s s e s p r o d u c e t o r q u e s on t h e s o l i d p a r t s o f t he s t a r , and thus p e r t u r b t h e i r a n g u l a r v e l o c i t i e s . The re a r e a l s o v a r i a b l e e l e c t r o m a g n e t i c s t r e s s e s on the s o l i d b o u n d a r i e s wh i ch can e x e r t t o r q u e s . Fo rmu lae f o r t h e p e r t u r b a t i o n iir i n t h e a n g u l a r v e l o c i t y o f the c r u s t caused by the p e r t u r b a t i o n s i n the m a g n e t i c s t r e s s a r e g i v e n be low. S i m i l a r f o r m u l a e can be d e r i v e d f o r a p o s s i b l e s o l i d c o r e . F o r s i m p l i c i t y , the r e s u l t s a r e g i v e n f o r t h e s p e c i a l c a s e o f an a x i s y m m e t r i c m a g n e t i c f i e l d whose a x i s c o i n c i d e s w i t h the r o t a t i o n a l a x i s . F o r t he same r e a s o n , the p e r t u r b a t i o n s In the magne t i c s t r e s s a r e l i n e a r i z e d . I t i s a s t r a i g h t -f o r w a r d m a t t e r t o g e n e r a l i z e t h e s e f o r m u l a e to more c o m p l i c a t e d f i e l d g e o -m e t r i e s and t o t h e n o n - l i n e a r c a s e . S i n c e the e f f e c t s o f r o t a t i o n on the - 4 2 / s t r u c t u r e o f the s t a r a r e o f o r d e r 1 0 Q^, the s l i g h t l y n o n - s p h e r i c a l shape o f b o u n d a r i e s i s n e g l e c t e d . The t o r q u e s e x e r t e d by f l u i d m o t i o n s l i e a l o n g the z a x i s . The v i s c o u s t o r q u e on the c r u s t i s g i v e n by [ 1 1 . 9 4 ] N ' = 2 , T R 3 \ 0\ s i n 2 8 de , Jy L where ia the (r,<f>) component of the fluid shear stress. The magnetic torque on the crust due to perturbations in the magnetic stresses at the inner crustal surface is [ 1 1 . 9 5 ] N' = 2 T T R 3 J sin 2 e de f B K f B [Bl / or 1 f = ZTTR 3 / sin26 de c Jz In [ 1 1 . 9 5 ] , the symbol [B'] stands for the jump in B' at E going from the <P <P fluid side £ to the solid side E + . In the infinite conductivity limit B ; ' E I ( B ; IE- + B ; I E + ) = \ B ; ' Z - -I[b;]' s i n c e > ; i z + • w n i s h e 8 i n t h i s limit. The total crust angular acceleration is the sum of a term ftQ due to quasi-steady external torques and any precessional motion, and a term ft' due to the shear stresses and the perturbations in the magnetic stresses communicated by the fluid interior: [ 1 1 . 9 6 ] ft = $ + ft' . o If the crust is treated as a spherical rigid body with a scalar moment of inertia I, then [ 1 1 . 9 7 ] ft' = (N^ + N^/I . The power spectrum of the perturbations ft'(t) in the angular acceleration of the crust due to charged fluid flow will be used in Chapters V and VI to determine the periodicities in some computer simulated flows. 57 Figure 2. Charged f l u i d mass f r a c t i o n as a function of t o t a l density i n the f l u i d region. The d ensity i s i n grams. 58 59 Figure 3. Logs of squares of sound speeds as a fu n c t i o n of t o t a l density i n the f l u i d region. The change i n slope of the s 2 graph above I O 1 5 gm r e l a t i v i s t i c . 15 -3 ° 10 cm i s due to the neutrons becoming more 60 CO o L d LU Q_ CO =) O CO Li_ O CO LU CC < ZD o CO LL. O 20.4, 20.0 19.6 19.2 18.8 CO O 18.4 - j 13.0 1 1 1 13.5 14.0 14.5 LOG DENSITY J 15.0 15.5 6 1 Figure 4. p*'"'' and pK£1 as functions of radius i n a r o t a t i n g neutron o o star model. P Q ° ^ ^ s t n e s p h e r i c a l l y symmetric part of the t o t a l mass density, p^ i s the r a d i a l part of the density perturbation due to rota t ion. The radius i s measured i n km. 62 63 Figure 5. Toroi d a l magnetic f i e l d perturbation required to balance the quasi-steady Poincare" f o r c e . The unperturbed magnetic f i e l d i s assumed to be a uniform one p a r a l l e l to the r o t a t i o n axis 11 and to be 10 gauss. The r o t a t i o n axis i s v e r t i c a l ; the equator i s ho r i z o n t a l . The t o r o i d a l f i e l d i s given i n a contour p l o t . Contour plots with s i m i l a r l a b e l i n g conventions are used throughout the thesis f o r the azimuthal parts of B' and v . The maximum of the absolute value of the quantity ( i n t h i s case, max|B^|, which i s 4.93 x 10 5 gauss) i s given on each p i c t u r e . A value of +10 i s assigned to the point where t h i s maximum occurs i f the quantity i s p o s i t i v e there. Otherwise, -10 i s assigned to i t . MAX B A = 4 .93 x IO5 GAUSS 1 65 CHAPTER I I I : TYPES OF CHARGED FLUID WAVES In the l a s t chapter, the equations of the charged f l u i d dynamics were derived. The present chapter uses these equations to f i n d the general c h a r a c t e r i s t i c s of charged f l u i d waves. Dispersion Relations When the length scale A of a free o s c i l l a t i o n i n the charged f l u i d i s small compared to the distance R over which the equilibrium s o l u t i o n v a r i e s , d i s p e r s i o n r e l a t i o n s may be derived r e l a t i n g the frequency w and the l o c a l ->-wavenumber k of the o s c i l l a t i o n . Even when A - R, the dispersion r e l a t i o n s derived upon neglecting Inhomogeneities i n the equilibrium s o l u t i o n can serve as us e f u l guides i n determining the r e l a t i o n between the length and time scales of the o s c i l l a t i o n s . In the present section, such inhomoge-n e i t i e s are neglected and the charged f l u i d i s treated as i n f i n i t e i n extent. The motivation for t h i s somewhat crude treatment of the waves i n the f l u i d i s that i t enables one to determine e a s i l y the possible types of waves and th e i r prime c h a r a c t e r i s t i c s . I t i s assumed throughout t h i s thesis that the length scales A and time scales T of charged f l u i d disturbances are such that the conditions derived i n the previous chapter for the charged p a r t i c l e s to be treated as a f l u i d are obeyed. These conditions are s a t i s f i e d f o r an extremely wide range of A and T. In p a r t i c u l a r , they are e a s i l y s a t i s f i e d i n the case A = R. To a large extent, the work i n th i s section i s an adaptation to neutron stars of the general theory of waves i n r o t a t i n g e l e c t r i c a l l y conducting f l u i d s i n the presence of a magnetic f i e l d , which has recently been reviewed by Acheson and Hide (1973). When i t i s l i n e a r i z e d and s p a t i a l g r a d i e n t s i n p c Q a r e n e g l e c t e d , the mass c o n s e r v a t i o n e q u a t i o n becomes [ I I I . l ] s " 2 | | — + V . v = 0 . co 9 t c The n o t a t i o n P ' = p ' / p has been u sed In [ I I I . l ] f o r t he r e d u c e d p r e s s u r e C CO p e r t u r b a t i o n . I s h a l l n e g l e c t s e v e r a l terms i n t he l i n e a r i z e d momentum e q u a t i o n . One i s t h e Po incare " t e r m , w h i c h i s o m i t t e d b e c a u s e o n l y f r e e o s c i l l a t i o n s a r e o f i n t e r e s t h e r e ; . I s h a l l a l s o n e g l e c t -p 7*' and -p'v"i|i . The g r a v i t a t i o n a l ° CO C O term - p c o $ * V was shown i n t he p r e v i o u s c h a p t e r t o be f a i r l y s m a l l . A l s o , s i n c e * ' s a t i s f i e s an e l l i p t i c e q u a t i o n , t h i s g r a v i t a t i o n a l te rm canno t be r e s p o n s i b l e f o r any new t y p e s o f waves . In t h e buoyancy te rm " P ^ v ^ * ^ Q v a r i e s r o u g h l y l i n e a r l y w i t h d i s t a n c e to t he c e n t r e o f t he s t a r ( see C h a p t e r I V ) . T h u s , buoyancy can be p r o p e r l y t r e a t e d i n t h e d i s p e r s i o n r e l a t i o n a p p r o a c h o n l y i f A << R. However, i n t h i s s h o r t w a v e l e n g t h l i m i t , ~ P ^ ^ 0 i s r o u g h l y A/R t imes t h e s i z e o f t he p r e s s u r e g r a d i e n t term - V p ' and i s thus n e g l i g i b l e . T h e q u e s t i o n o f the c e f f e c t s o f buoyancy on modes w i t h l e n g t h s c a l e s A = R i s a d d r e s s e d i n C h a p t e r s IV t h r o u g h V I f o r each c l a s s o f o s c i l l a t i o n i n d i v i d u a l l y . The damping terms ( v i s c o s i t y and s u p e r f l u i d d rag ) a r e t e m p o r a r i l y n e -g l e c t e d ; they a r e r e - i n t r o d u c e d i n a l a t e r s e c t i o n where t h e y a r e shown t o damp most o s c i l l a t i o n s o n l y l i g h t l y . The r e s u l t i n g h i g h l y s i m p l i f i e d momentum e q u a t i o n s a r e 3v [III.2] -r-S- = -2&n x £ - VP ' - T - T — B x (V x f ••)•• . 3t o c Atrp o co The l i n e a r i z e d i n d u c t i o n e q u a t i o n i s [ I I I .3] |f- = V x ( v , x t ) . at c o 67 - > • - > ' • If a l l v a r i a b l e s vary i n time and space as r , then [III.4] -iw s ~ 2 P' + i l k • v = 0 C O c [III.5] -iojv = -2ft x v - ikP' - -.—-—- B x ( i k x B * ) c o c 4irp o co [ I I I . 6 ] -ioJB' = i k • B v - i t • v B O C C O Sound Waves An equation f o r P' may be derived by taking the sc a l a r product of k with [III.5] and using [ I I I . 4 ] : [III.7] ( k 2 - u) 2s" 2)P' = -it • [-2ft x v - . 1 B x (ik x B')] . V co o c 4irp o • co ••; • Note than i n the l i m i t ft , B —> 0, [III.7] reduces to the Fourier analyzed o o form of a wave equation f o r P'; i n the absence of r o t a t i o n and a magnetic f i e l d , the only kind of wave which can propagate i s a sound wave.with d i s -persion r e l a t i o n [III.8] a) - s k (ft = B = 0) . co o o If the magnetic body force and the C o r i o l i s force are s u f f i c i e n t l y small, sound waves with the dispersion r e l a t i o n [III.8] can propagate. The t o t a l a cceleration, C o r i o l i s a c c e l e r a t i o n and magnetic acceleration have 2ft 2 2 — 2 sizes i n the r a t i o 1 : : V • k co , res p e c t i v e l y , where V. i s the Alfv£n OJ A r A speed: 68 If to = s kj then the three accelerations are in the ratio 1 : 2ft Is k : co o co 2 2 -2 -12 2 V ^ / S C Q - 1 : 10 ^2*6 : ^11' S° t* i e c o n ^ i t i o n s necessary for the propagation of sound waves which are not much affected by rotation or the magnetic field are' [III.10] ft >, « 102 ;, B.r « 106 . These conditions are well satisfied unless the magnetic field or. the angular velocity are so large that they have a substantial effect on the equilibrium structure of the star. Such cases are not considered in this thesis. Chapter IV is.devoted to a fuller examination of sound waves in the charged fluid. Inertial Waves. And Hydromagnetic-inertial Waves 2 If ID « s k, then the left-hand side of [III.7] becomes k P' to a co 2 2 2 good approximation. Dropping to Is In comparison to k is equivalent to co replacing the wave equation for P' with an elliptic equation in which pressure effects are instantaneously propagated. In such an approximation, equation [III.4] becomes the incompressibility condition [III.11] ik • v - 0 c and [III.6] simplifies to [III.12] -iwf' = ik • B v . . o c Multiplying [III.5] by ik and using [III.11-12] produces [III.13] [(-ioo)2 - (V. • i t ) 2 ] it x v = (2ft • it)(-ito)v . A c o c Since the l i n e a r operation ( k x ) has the eigenvalues 0, l i k , of which the f i r s t i s excluded by the i n c o m p r e s s i b i l i t y condition, the dispersion r e l a t i o n which follows from [III.13] i s : (2ft • k) [III.14] O J Z + ~ OJ - (V • kr = 0 , which i s equation (4.11a) of Acheson and Hide (1973). A c h a r a c t e r i s t i c dimensionless parameter associated with [III.14] i s (2ft • k) , , [III.15] 5 = — ° — - 10 ft9Aft B..7 . 1 VT t 2 o i i kV. • k A Note that the estimate of the magnitude of t% i n [III.15] i s not accurate i f k and ft*Q are almost perpendicular. However, the estimate Is accurate f o r almost a l l d i r e c t i o n s of lc. In order to save space, statements such as these about the dependence of the s i z e of various quantities and e f f e c t s on the d i r e c t i o n of k w i l l be omitted f o r the remainder of t h i s chapter. If £ >> 1, then the roots of [III.14] s p l i t into two branches, widely separated i n frequency. The high frequency branch i s formed by i n e r t i a l waves: 2ft" • £ [III.16] co = ± — ° k ^ 2 0 0 Q 2 s e c " The lower frequency branch i s formed by hydromagnetic-inertial waves: (V. * t ) 2 k , 0 _ 9 _•• _-, [III.17] OJ = t — = 10. B,, \, ft0 sec ^ . 2ft • t 1 1 6 2 o Hide (1971) made an important contribution to the understanding of the charged f l u i d dynamics when he noted that £ i s generally >> 1 for o s c i l l a t i o n s with length scales A - R i f the magnetic f i e l d i s not very 12 much larger than 10 gauss. As a r e s u l t , two modes of the charged f l u i d are the i n e r t i a l modes and the hydromagnetic-inertial modes. 70 The parameter E, can be « 1 i f the magnetic f i e l d i s extremely strbng 4 -4 ^ B l l > : > 1 0 } o r i f t h e w a v e l e n 8 t n i s extremely small ' ( A << 10 ). Under such conditions r o t a t i o n has only a minor e f f e c t on the waves, which are e s s e n t i a l l y pure hydromagnetic (Alfveh) waves. In this thesis, I am pr i m a r i l y concerned with the less extreme conditions i n which £ •>>'.!. I n e r t i a l Waves The periods of i n e r t i a l waves are t y p i c a l l y h a l f the r o t a t i o n period P, or s l i g h t l y longer: [III.18] " T j . X \ p • _4 Since the periods of sound waves are £ R / S C Q ~ 10 s e C > there i s no overlap between the frequency range of sound waves and that of i n e r t i a l waves i n neutron s t a r s . It follows from [III.12] that over the time scales T T the r a t i o of the magnetic energy of these waves to t h e i r k i n e t i c energy i s ; extremely small: [IH.19] - f ^ f - T 2 ~- IO" 8 ft!2 A " 2 B 2 - P v 2 6 11 2 co c The r a t i o of the magnetic body force to the C o r i o l i s force over time scales _2 i s also equal to $ . Under these circumstances the e f f e c t of the mag-netic f i e l d can be completely neglected. Hydromagnetic-Inertial Waves The periods of these highly d i s p e r s i v e waves can be extremely long: [III.20] = B" 2 A 2 ft2 months. o If £ >> 1, there i s no overlap between the periods of these waves and the periods of i n e r t i a l waves. Unlike pure Alfv£n waves, the group v e l o c i t y 3to/3k of hydromagnetic-inertial waves i s not n e c e s s a r i l y p a r a l l e l to The magnetic energy i n hydromagnetic-inertial waves i s very much larger than the k i n e t i c energy: I m . 2 1 ] I B ^ / | , , ? 2 „ 1 0 8 fi-2 n 2 x 2 _ 2 V c i n contrast to the i n e r t i a l waves [III.19] where t h i s r a t i o i s reversed. Over time scales T , the magnetic body forces are roughly the same size as the C o r i o l i s f o r c e : | [B x (V x B ' ) ] / 4 T T | [III.22] ° — ; « 1 , 12ft x p v I . ' O C O c 1 whereas the t o t a l force on a f l u i d element i s smaller than either by the 2 factor £ . Damping Due To V i s c o s i t y And Superfluid Drag For order of magnitude estimates of viscous damping times, the use of an i s o t r o p i c v i s c o s i t y tensor i s s u f f i c i e n t . For incompressible flow, iso-t r o p i c v i s c o s i t y contributes an extra term v £V to [ I I I i 2 ] . The e f f e c t of t h i s term on the dispersion r e l a t i o n i s to change equation [III.14] to 2ft" • t + ±\ c [III.23] to2 + + . i v k 2 to - (V A • k ) 2 = 0 . For i n e r t i a l waves, i t follows from [III.23] that „->- ->-2ft • k . [III.24] to = + — \ i v k 2 . .. k c The viscous damping time f o r i n e r t i a l waves i s thus [III.25] T,(TI) = 1/v k 2 = 10 A 2 T 2 sec. V c 6 7 It i s straightforward to show that the viscous damping time f or sound waves 2 i s also about l / v c k . For hydromagnetic-inertial waves, the dispersion r e l a t i o n which follows from [III.23] i s ( v . - k ) 2 k(V. • k ) 2 _ . [III.26] co = — ~ 2 «' + „ A V - U " v k Z • +(2fi • k)/k + i v k 2ft • k c o c o (HI) 2 (I) The viscous damping time f o r these waves i s thus about £ times T y : [III.27] T T ( 7 H I ) = 5 2/v k 2 - 100 ft2 \t B 7 2 T 2 years. The very long viscous damping time f or hydromagnetic-inertial waves when 2 -2 ^ i s >> 1 r e s u l t s from the f a c t that only a small f r a c t i o n £ of the energy i n these waves resides i n f l u i d motion. Equation [III.23] can also be used to determine the viscous damping times f o r waves with wavelengths so small that £ < K 1, i . e . when k >> ft^/V^. -4-Under such conditions the fa c t o r l(2ft • k)/k i n [III.23] may be dropped. o 2 The viscous damping time f o r such Alfven waves i s - 2/v ck , which for 2 2 — 6 2 2 — 2 k >> ft /VA i s << 2V./v ft = 10 B , . T., ft. sec. Such waves are thus damped o A A c o 11 7 2 i n a very short time. To analyze the e f f e c t s of su p e r f l u i d drag on the three basic types of waves, I assume that v = 0 i n the ro t a t i n g frame. These e f f e c t s may then 2 be determined by formally replacing v^k i n [III.23-27] by 1 / T Q - For example, the i n e r t i a l wave damping time due to su p e r f l u i d drag i s j u s t [III.28] T D I } = T D • For hydromagnetic-inertial waves the superfluid drag damping time is [III.29] T < H I ) = ? 2x D . Since x D can be days to years, T ^ H ' 1 " ^ can be up to 108 years if £ = 104. 2 When T q >> l/v ck , which is true for a wide range of conditions, superfluid drag is a less important damping mechanism than is viscosity. Sound waves can be damped by collisionless processes such as Landau damping, but for wavelengths longer than the electron mean free path, collisionless damping is insignificant (Pippard, 1955). Wave Damping By Heat Conduction Since the charged fluid is not completely degeneratej the mass density p- depends on the temperature. If the temperature is perturbed by T', there will be a density perturbation (9p /9T) T', which is acted upon by gravity to produce a thermal buoyancy force -(9p /9T) T'V$ .Because temperature C Pc ° perturbations which are created by charged fluid motion can be damped by heat conduction, the thermal buoyancy force can damp waves in the charged fluid. The appropriate form of the heat flux equation (Acheson and Hide, 1973) when linearized is [III.30] ct | r ' + v • VT I = <V2T' . l o t c oi a is the heat capacity per unit volume, and K is the thermal conductivity. The viscous heating term a' • V v £ does not appear in [III.30] because i t is second order in small quantities. An Ohmic dissipation term of the form ) does not appear because of the assumed infinite c c o electrical conductivity. Because the protons are superfluid, the heat capacity arises almost entirely from the relativistic electrons (§58; Landau and Lifschitz, 1969) [III.31] a = _ ( 3 T T2 ) 2 / 3 ,2 2/3 16 3cii "B e k n T = 10 T_ erg/deg-cm The thermal conductivity as calculated by Flowers and Itoh (1975) is not strongly density dependent in the fluid region, and is roughly equal to [III.32] 25-1 K - 10 T^  erg/deg-cm-sec. The corresponding thermal diffusivity x is [111.33] 9 —2 2 —1 X = </ct * 10 T7 cm sec . If the gradient in T q is treated as roughly constant, a Fourier analysis in space and time of [III.30] yields [III.34] •T' = 2 vc VT ico - x k The temperature gradient in the fluid region can be roughly estimated by equating the temperature gradient heat f lux JKVTJ in the interior with the 4 flux CTSBTS of blackbody radiation from the stellar surface at temperature Tg. a g B is the Stephan-Boltzman constant. The result is [III.35] a T 4 -fP^ 10~ 1 3 cm"1 for T s .- 106 K, which is very much larger than R_1. The mass density of the charged fluid is due primarily to the protons, so the temperature dependence of p at c temperatures low compared to the proton Fermi temperature is [III.36] P (T ) = p (T = 0) C O C O Thus the thermal buoyancy force i s [III.37] - l -^ l T'V* = 0 I-3-2] S i ° - v » 3T o 'co U 2 V ° \ £Fp / (iu - x kV o 2 - 3 - 2 - 2 -1 For low frequency waves such that to << xk - 10 Tn A , sec , ito can be 2 2 neglected ,in comparison to yk i n [III.37]. For o £ xk , which i s probably the case f o r both sound and i n e r t i a l waves, io> may not be neglected i n 2 comparison with xk l n estimating the magnitude of the thermal buoyancy force. However, our i n t e r e s t i n [III.37] l i e s i n the damping e f f e c t that 2 2 i t has. This damping ar i s e s only from the xk term i n ito - X k . Thus i f thermal buoyancy damping time x i s defined by the condition that the thermal buoyancy force with to set equal to zero be about p V / T i n co c magnitude, then i t follows, from [III.37] that [III.38] T ( B ) * -^P-\ k..T B o^  2 T 2 .2 o R z _n4 - 4 - 2 I 0 I The associated damping time f o r sound waves and i n e r t i a l waves i s T , and f o r hydromagnetic-inertial waves i s E, T . Wave damping due to heat conduction seems to be less important than viscous damping i f Ty £ 1, but f o r » 1, i t could be the dominant damping mechanism. Resonant Mode-Mode Coupling I t was previously concluded that hydromagnetic-inertial waves i n the charged f l u i d tend to have extremely long damping times. One of the e s s e n t i a l assumptions that led to t h i s conclusion was that a l l waves are of such small amplitude that non-linear e f f e c t s could be neglected. In t h i s section, some of these non-linear e f f e c t s are b r i e f l y examined. 76 The question i s : can the energy i n a hydromagnetic-inertial wave be transferred to other waves at such a rate that i t can be dis s i p a t e d i n a time much less than the energy damping time predicted by a l i n e a r theory? For the moment, neglect the l i n e a r damping e f f e c t s and consider the simpler question of whether energy can be resonantly transferred from a hydromagnetic-i n e r t i a l wave with r e a l frequency co^  and wavenumber to two waves with r e a l frequencies and wavenumbers io 2, k"2; to^j-tc^. The resonance conditions are (Sagdeev and Galeev, 1969): [ I I I .39] cuj = u> + a> [III.40] k x = £ 2 + t3 . Several p o s s i b i l i t i e s are immediately excluded by [ I I I .39-40] and by the dispe r s i o n r e l a t i o n s [ I I I .8, 16 and 17], A hydromagnetic-inertial wave cannot beat with i t s e l f , because i t s di s p e r s i o n r e l a t i o n i s non-linear i n k. Two hydromagnetic-inertial waves cannot undergo fusion into either a sound wave or an i n e r t i a ! wave because of the frequency mismatch. For the same reason, a hydromagnetic-inertial wave cannot decay into an i n e r t i a l wave and a sound wave. The. coupling between a hydromagnetic-inertial wave and either two i n e r t i a l waves or two sound waves should be extremely weak: as previously shown, a l l but a small f r a c t i o n £ of the energy i n a hydromagnetic-inertial wave resides i n the magnetic f i e l d , whereas f o r i n e r t i a l waves and sound waves almost a l l the energy resides i n the k i n e t i c energy of the wave. However, even i f the coupling were not extremely weak, energy cannot be resonantly transferred from a hydromagnetic-inertial wave to either a p a i r of i n e r t i a l waves or sound waves. This may be demonstrated as follows. Because the damping times f o r these high frequency waves are much shorter than the periods of hydromagnetic-inertial waves, i t s u f f i c e s to consider the case i n which i n i t i a l l y almost a l l the energy i s i n the low frequency wave. Resonant transfer of energy i s possible only i f the wave which i n i t i a l l y possesses most of the energy has a higher frequency than the other two waves (Sagdeev and Galeev, 1969). However, the frequencies of hydromagnetic-inertial waves are very much lower than the frequencies of either i n e r t i a l waves or sound waves. In view of these arguments, i t seems u n l i k e l y that the energy i n a hydromagnetic-inertial wave of period T can be resonantly transferred Hi to sound waves or. i n e r t i a l waves and then diss i p a t e d i n a time << T . . n i Boundary Layers Acheson and Hide (1973) review the theory of hydromagnetic boundary la y e r s . In t h i s s e c t i o n , t h e i r general r e s u l t s are applied to the problem of boundary layers at the f l u i d - s o l i d i n t e r f a c e s within neutron s t a r s . In the charged f l u i d near the boundary surface E with unit normal n, a l l p hysical v a r i a b les <J> are s p l i t into two parts, a part <j> which i s roughly constant throughout the boundary layer and a part <J> which v a r i e s s i g n i f i c a n t l y throughout the boundary layer: [ I I I . 41] <f> • = 4 + <(> . If the usual boundary layer approximations are made and the curvature of E i s neglected (Acheson and Hide, 1973), t h e , l i n e a r i z e d magnetohydrodynamic equations may be Fourier analyzed, i n time and space. The time and space v a r i a t i o n s of a l l boundary layer quantities are assumed to be of the form [III.42] <J) oc exp(icot + yz) , 78 where z measures the distance from E along the unit normal n. The r e s u l t i n g algebraic equations have n o n - t r i v i a l solutions only i f a disp e r s i o n r e l a t i o n i s s a t i s f i e d : [III.43] [ ( i u ) ( i u - v Y 2) " (v\ • n)2y2]2 - (2fl • n ) V = 0 . C A O The boundary layer thickness A (to) i s defined by [III.44] A = -1/Re(y) . There are two c h a r a c t e r i s t i c frequencies which appear i n the disp e r s i o n r e l a t i o n [ 1 1 1 , 4 3 ] . One Is 2ft Q • n and represents the influence of r o t a t i o n on the boundary layer structure. The other c h a r a c t e r i s t i c frequency i s ->• -> 2 (V • n) /v , which i s the rate at which viscous d i f f u s i o n can proceed across A C the boundary layer along the f i e l d l i n e s . I t was shown i n Chapter II that the magnetic f i e l d at the c r u s t a l E i s probably almost e n t i r e l y p o l o i d a l or almost e n t i r e l y t o r o i d a l . In the p o l o i d a l case, [III.45] (V. •• n) 2/v = V 2/v = 1 0 " 2 B?. T 2 s e c " 1 . A ., C A C XX I In the t o r o i d a l case, |V. • n| i s << V , and the c h a r a c t e r i s t i c frequency (V. • n ) 2 / v i s « 1 0 " 2 B 2 T 2 s e c " 1 . A c 11 7 For waves with frequencies to i n the range [III.46] • 1-2^ • n| > OJ » (V A • n ) 2 / v , o 1 A c the magnetic term i n the dispersion r e l a t i o n [III.4 3 ] i s n e g l i g i b l e . The boundary layer thickness f o r such waves i s ,1/2 [III.47] A = ' 2, c 1 0 4 Q~1/2 T' 1 cm, + \2& • n| o 1 which i s about a hundredth of R £ i f &2 - T ? = 1. These boundary layers are 79 the Ekman boundary layers often found at the boundaries of rotating fluids. In the low frequency range [III.48] OJ « (V. • n)2/v , A c the viscous term in the dispersion relation [III.43] is negligible. The boundary layer thickness for such frequencies is V • n [III.49] A - A . /2QQ • nco It follows from [III.49] that i f the charged fluid is oscillating in a hydromagnetic-inertial mode whose length scale is X and if the magnetic field is predominantly poloidal at E, then [III.50] A - X , i.e. there is effectively no boundary layer. If the magnetic field is predominantly toroidal near E, then [III.51] A « A . The presence of a boundary layer in the toroidal case can be understood in terms of a simple physical picture in which the field lines near E can easily "slip" by one another parallel to E. In the poloidal case, tension in the field lines acts to oppose any such motion and thus irihibits the presence of a boundary layer. 80 Simplified Charged Fluid Equations 2 It has been shown in this chapter that for £ >> 1, there are three basic time scales in the charged fluid dynamics corresponding to the periods of sound Waves, inertial waves and hydromagnetic-inertial waves. For each time scale, one or more of the terms in the charged fluid equations are negligible. In this section, a simplified set of charged fluid equations is derived for each of the three time scales by dropping the negligible terms. . There are great advantages to the use of such simplified equations. Not only do they make the mathematics easier, but they also enable one to clearly separate the effects of different time scales. All mode-mode coupling is neglected in such an approach. The remaining three chapters in this thesis are devoted to a more thorough examination of these approxi-mate equations and some of their solutions. Sound Time Scales Over the time scales of sound waves, i t has been shown that the magnetic body forces are negligible and the Coriolis force is small. A simplified set of equations for these time scales is thus [III.52] 3 — + V • (p v ) = 0 t CO c [III.53] | r (p v ) « -Vp' - p'V> . d t C O C r C C O The viscous term could be added to [III.53] i f viscous damping is of interest In the absence of i t , the boundary condition on v"c is the slip boundary condition [11.89]. 81 I n e r t i a l Time Scales For time scales longer than those of sound waves, i t was shown that 9P'/ot was n e g l i g i b l e . However, i t was assumed that gradients i n P c q could be neglected, which i s not necess a r i l y the case when A - R. This suggests that the i n c o m p r e s s i b i l i t y condition v"*v = 0 i s not a good approximation over these time scales, but that the a n e l a s t i c condition (Batchelor, 1953) V • (p„rtv ) = 0 i s . C O c The a n e l a s t i c condition i s a good approximation i f [III.54] .|v" • (p v )I « p V /A . 1 C O C 1 C O c From the mass conservation equation i t follows that [ i i i . 5 5 ] |v • (p v )| = | a P 1 / a t| * p ' / s 2 T ' co c 1 c 1 c co From the momentum equation f o r time scales T ft \ i t follows that o [III.56] p' = 2Aft p V . c o co c Together., [III.54] through [III.56] imply that i f A 2ft [111.57] - y ^ « 1 , S T . C O then the a n e l a s t i c condition [III.58] v" • (p v ) = 0 C O c i s a good approximation. Condition [III.57] i s e a s i l y s a t i s f i e d f o r a l l T £ fl The use of [III.58] mathematically " f i l t e r s " out a l l sound waves from the solutions to the f l u i d equations, j u s t as the i n c o m p r e s s i b i l i t y condition V*v c= 0 does f o r f l u i d s of nearly constant density. 82 The magnetic body forces are n e g l i g i b l e over the time scales of i n e r t i a l waves, so the s i m p l i f i e d momentum equation i s [ I I I . 5 9 ] | r (P v ) = -2$ x p v - vp' - . p ' v f ' • 9t co c o co c *c c o The viscous term may be added i f desired, and the s l i p boundary condition on V £ replaced by a n o - s l i p condition. The charged f l u i d dynamics over the time scale of i n e r t i a l waves i s examined more thoroughly i n Chapter V. Hydromagnetic-inertial Time Scales Since [III.57] i s s a t i s f i e d for hydromagnetic-inertial time scales as w e l l as i n e r t i a l time scales, the an e l a s t i c condition [III.60] V • (p v ) = 0 C O c i s a good approximation on t h i s time scale also. Over these time scales, i t was shown that the 9(p v )/9t term i n the C O c f l u i d momentum equations i s smaller than the C o r i o l i s and magnetic body force 2 terms by a fac t o r £ . Dropping t h i s r e l a t i v e l y small term from the momentum equations r e s u l t s i n [III.61] ft = -2^ x p v - v"p' - p'V> - r - B x(V x B')- 7 - B 1 x(V x B o co c c c r o 4TT O 4TT C Dropping the 3 ( p C O V c ) / 9 t term mathematically f i l t e r s out i n e r t i a l waves i n the same manner as dropping the 9p^/9t term i n the mass equation f i l t e r s out sound waves. I t was previously shown that resonant coupling of hydro-mag n e t i c - i n e r t i a l modes and the higher frequency modes i s very small so the use of equations [III.60-61], which s t r i c t l y eliminates such mode coupling, should not lead to q u a l i t a t i v e l y i n c o r r e c t r e s u l t s i n th i s regard. 83 A l l e x p l i c i t time evolution over these time scales resides i n the magnetic induction equation [III.62] 11^- = V x (v x l ) . dt C O The appropriate v e l o c i t y boundary conditions are the n o - s l i p conditions i f 13 • n|j, 4 0, and s l i p boundary conditions i f B • n|^, = 0. The viscous term may be included i n the momentum equation [III.61], but i t a f f e c t s the f l u i d flow only i f B • n|^ =0. The charged f l u i d dynamics over the time scale of hydromagnetic-inertial waves i s examined more thoroughly i n Chapter VI. 84 CHAPTER IV: SOUND MODES In the preceding chapter, i t was shown that the highest frequency oscillations of the charged fluid are sound waves. The present chapter explores the properties of these waves in greater detail. One reason for studying these waves in more detail is to determine the effect of the buoyancy force -p^ViJJo on them. This force cannot be properly treated in the dispersion relation approach of Chapter III. The mass conservation equation is 3p* [IV. 1] s 1 T - ^ + v" • (p v ) = 0 . C O d t C O C The simplified momentum equation is From these two equations follows 2 32p' ? . 9 [IV.3] s —TT- = vp' + • Vp' + V ty p' . co 2 . c ro rc ro*c d t Equation [IV.3] differs from the usual wave equation in that buoyancy terms are present and the sound speed is position dependent. Standing wave solutions of [IV.3] satisfy (cf. equation [III.7]) [IV.4] s 2 {V2p» + v\|i • Vp1 + V2ij) p'} = -o>2p' C O c o rc o c c In [IV.4], distances are measured in units of the inner crustal radius Rc> and sound speeds in units of S C Q(0)> t n e sound speed at the centre of the star. The frequency co is measured in units of s (0)/R . co c Because the magnetohydrostatic solution departs from spherical symmetry 2 2, - 4 2 only by a factor ftQR /$o - 10 ti^, i t Is a very good approximation to take ; 85 2 i|»o and S c q to be functions of r only. In this approximation, [IV.4] is separable in spherical polar coordinates. Solutions of the form [IV.5] p^ (r,6,<|>) = f^(r)Rm(cose)eim* were sought, fg satisfies the ordinary differential equation [iv.6] s 2 o{f- +. (f + *;>f. + o;« + -m±21 ) f } = r where ' means d/dr. In the absence of gravitational effects and variation in the sound speed, [IV.6] becomes the spherical Bessel equation, with solutions j (ti^-.r) and n^Cu^r). However, the gravitational terms may be neglected only for r << 1, since i|> is of order 1. The sound speed vari-o ation may be neglected only if the fluid region is very thin. Because of the neglect of viscous effects, slip boundary conditions should be used at the outer boundary r = 1 and at any inner boundary r = a, representing a,solid central core. The r-component of the momentum equation then implies 3p' dip [IV.7] 0> - r — ^ - p' -r-°-3r *c dr at the boundaries. In terms of f^ , [IV.7] becomes [IV.8] f j (a) + f^ (a)^ »(a) = 0 [IV.9] f ^ U ) + f ^ C D ^ ' d ) « 0 . If there is no inner core, the appropriate boundary conditions at the centre r = 0 are that the pressure perturbation and its derivatives are continuous: [IV.10] f'(0) = 0 o [IV. 11] f (0) - 0 if £> 0 . 86 The eigenvalue d i f f e r e n t i a l equation for f was solved numerically, subject to the appropriate boundary conditions, for various values of Z and a. 1 The functions s (r) and :ty (r) used i n the c a l c u l a t i o n were based CO o on the magnetohydrostatic equilibrium model of Chapter I I . They were determined by l e a s t squares f i t t i n g 2 s i x t h order polynomials to the spher-2 i c a l l y symmetric parts of s (r,9) and ty ( r , 6 ) . This f i t t i n g produced CO o [IV.12] s 2 (r) -.1 - .8697 r 2 + .3277 r 4 - .0650 r 6 co and [IV.13] ^ o ( r ) = 1.2168 r 2 - .0129 r 4 + .0593 r The f i v e lowest p o s i t i v e eigenf requencies to^ , for each of £ = 0, 1 and 2 are given i n Table VII, f o r a v a r i e t y of core s i z e s . For every eigenfrequency tog > 0 there i s another one -tog. Note that a l l the eigenf requencies OJ, are r e a l ; the buoyancy force — p Q neither damps the waves nor renders them unstable. The corresponding eigenmodes f^ f o r some core sizes are plotted i n Figures 6, 7 and 8, r e s p e c t i v e l y . The zero frequency eigenmodes i n the s p h e r i c a l l y symmetric case represent the pressure diff e r e n c e between two equilibrium configurations with d i f f e r e n t c e n t r a l pressures. I t can be seen from the figures that the eigenmodes are q u a l i t a t i v e l y l i k e s p h e r i c a l Bessel functions. Neither the buoyancy force nor the gradients i n the sound speed has any q u a l i t a t i v e e f f e c t on the eigenmodes. See Appendix A for d e t a i l s of the numerical method used. 2 Using the program OLQF from the U.B.C. Computing Centre. Table VII Elgenfrequencies Of Some Standing Sound Waves JL a = 0.0 a = .25 a - .50 a = .75 . 0.0 0.0 0.0 0.0. 3.74 3.96 5.07 8.98 0 6.56 7.24 9.81 17.8 9.30 10.6 14.6 26.7 12.0 14.0 19.4 35.6 1.74 1.68 1.44 1.14 4.94 4.76 5.35 9.05 1 7.78 7.73 9.95 17.9 10.5 10.9 14.7 26.7 13.3 14.2 19.5 35.6 2.64 2.63 2.43 1.97 5.99 5.88 5.88 9.21 2 8.91 8.66 10.2 1.7.9. 11.7 11.6 14.9 26.8 14.5 14.7 19.6 35.7 88 Figure 6. Radial parts of some 1=0 charged f l u i d sound modes. In th upper diagram, i t i s assumed that there i s no s o l i d c e n t r a l core. In the lower diagram, a s o l i d c e n t r a l core of .25 t the inner c r u s t a l radius i s assumed to be present. Similar comments apply to the next two f i g u r e s , f i s p l o t t e d a g a i n s t r f o r the f i r s t f our modes. e imes 90 Figure 7 . Radial parts of some £ = 1 charged f l u i d sound modes. r 92 Figure 8. Radial parts of some ^ = 2 charged f l u i d sound modes. CHAPTER V: INERTIAL MODES AND INERTIAL FLOW In Chapter I I I , i t was shown that the charged f l u i d i s capable of o s c i l l a t i n g at frequencies equal to or somewhat less than twice, that of the r o t a t i o n a l frequency of the st a r . The present chapter explores some features of these i n e r t i a l normal modes, as well as charged f l u i d flow over these time scales. Free I n e r t i a l Modes Of O s c i l l a t i o n The equations to be solved are the an e l a s t i c condition [V.I] V" • (p v ) - 0 C O c and the s i m p l i f i e d momentum equation '3(p v ) [V.2] - f 0 - 5 - - -2$ x p v - Vp' - p'Vij, . 3t o co c c e o When Fourier analyzed i n time, equations [V.l-2] can be reduced to an eigenvalue equation f o r p \ To make t h i s reduction, i t i s convenient to define new va r i a b l e s s and n": [V.3] s = p v CO c [V . 4 ] h = -Vp' - p'ViJ, , c c ro i n terms of which the equations become [V.5] V • s = 0. [V.6] | f + 2^ x l = S 3t o Measuring distances i n units of R^ * s i n units of 2fi QR c and frequencies i n units of 2fto> [V.6] becomes, a f t e r Fourier analyzing i n time, -+• " ->- -> [ V . 7 ] -icos + z x s = h . Equation [ V . 7 ] can be solved f o r s: 2 2 2-> ~ ->• [V .8] co (to - l ) s = -co h + itoz x h + h z z Taking the divergence of [V . 8 ] and using the a n e l a s t i c condition produces 2 . 3h [V.9] 0.- -to V" • £ - icoz • v" x h + -T-2- . d z If an azimuthal dependence f o r p' of the form e l m <' > i s assumed, where c y l i n -c d r i c a l coordinates (z,x,<f>) are used, [V.9] becomes 9 -1 3IJJ ' [V.10] v V + h ' v P ' •+ [ v > + S i i L — T - 0 - ] p r c o c o x 3x c -2 = to 3 2p' 3ib 3p' 3 2 4) c T o c o , . 2 3z 3z ' 2 P c 3 z 3 z In the l i m i t ViJ^—> fi,-[V.10] reduces to an equation used to study the i n e r t i a l o s c i l l a t i o n s of a uniform density confined f l u i d (Greenspan, 1968) In the present case the buoyancy terms, which are not n e g l i g i b l e i f A = R^, render [V.10] inseparable i n both s p h e r i c a l polar and c y l i n d r i c a l co-ordinates. C y l i n d r i c a l Geometry Equation [V.10] i s separable i n c y l i n d r i c a l coordinates i f , instead of taking to be a function of r only, i t i s taken to be a function of x only. The r e s u l t i n g s i t u a t i o n i s that of an external g r a v i t a t i o n a l f i e l d everywhere perpendicular to the r o t a t i o n axis. This i d e a l i z a t i o n w i l l be .. • 96 made here. In addition, the boundaries w i l l be changed from two concentric spheres to two coaxial cylinders of f i n i t e length. The replacement of a .spherical geometry by a cylindrical geometry for reasons of mathematical simplicity is frequently used in the study of normal modes in rotating bodies, both magnetized and unmagnetized (Acheson and Hide, 1973). The resulting cylindrical eigenmodes and eigenfrequencies, while not quantitatively accurate, are usually found to possess many of the qualitative features of the spherical, eigenmodes and eigenfrequencies. A cylindrical idealization is in any case more r e a l i s t i c than a dispersion relation approach such as was used in Chapter III because the former takes into account, even i f only approximately, the effects of boundary conditions and inhomogeneities in the equilibrium solution, while the latter does not. In the cylindrical idealization, equation [V.10] becomes 3z 9x d 2 i p „ . -1. dty 2 dx x This equation has solutions of the form [V.12] P ; = Z(z) f(x) e±m* , where [V.13] Z(z) = cos(kz) or sin(kz) . The function f(x) satisfies the ordinary differential equation -1 2 [V.14] f " 4- (— + »')f' + {fr" + ( 1 + m" ^ y•.+ k2(co~2 -1) - \ } f = 0 x o o x o 2 X where ' means d/dx. S i n c e ty^(r) i s t o a v e r y good a p p r o x i m a t i o n j u s t p r o p o r t i o n a l t o r ( s ee [ I V . 1 3 ] ) ; t h e ^ q ( X ) c ho sen i n t h e p r e s e n t i d e a l i z a t i o n i s [V.15] ^ o ( x ) = | x 2 , where A i s a p o s i t i v e c o n s t a n t o f o r d e r 1. F o r t h i s c a s e , f s a t i s f i e s [V .16] f " + ( i + A x ) f 1 - f = - X Z f , x 97 2 where [V.17] A 2 = k 2(u> 2 - 1) + A (2 + moj " 1 ) . 2 To e v e r y e i g e n v a l u e -A t h e r e c o r r e s p o n d s a c o u n t a b l e i n f i n i t y o f n o r m a l 2 modes. The e i g e n f r e q u e n c i e s u ( k ) c o r r e s p o n d i n g t o an e i g e n v a l u e - A a r e g i v e n by [V.18] co(k)-•• ^ * ^ ~ 2 A + k ^ J 1 / 2 2 ( A 2 - 2A + k ) From [ V . 1 8 ] , i t i s p o s s i b l e t o l e a r n s e v e r a l t h i n g s about t h e e i g e n f r e q u e n c y s p e c t r u m . F o r e x a m p l e , i n t h e m a t h e m a t i c a l l i m i t A = 0 o f no g r a v i t a t i o n , [V .18 ] i m p l i e s t h a t |tu| < 1, o r , i n u n n o r m a l i z e d u n i t s , j co | < 2 f i Q , a r e s u l t w h i c h was g i v e n i n C h a p t e r I I I . A n o t h e r i n f e r e n c e i s t h a t i f A i s r e a l and [V .19 ] A 2 > 2A - k 2 - A 2 m 2 / ( 4 k 2 ) , t h e n t h e r e i s no damping o r g r o w t h o f t h e i n e r t i a l wave s . S i n c e v i s c o s i t y has been n e g l e c t e d , t h e p r o p e r bounda ry c o n d i t i o n s f o r t h i s p r o b l e m a r e s l i p bounda ry c o n d i t i o n s . I t f o l l o w s f r o m [ V . 4 , 6, & 12] ' t h a t t h e v a n i s h i n g o f s a t t h e uppe r and l o w e r s u r f a c e s o f t h e c y l i n d e r s , z w h i c h a r e t a k e n a t z = ~t 1, i m p l i e s t h a t dZ/dz v a n i s h e s t h e r e . 98 The permissible k-values are thus [V.20] k = n for cos, and [V.21] k - (n •+ -|)IT for sin, where n = 0,-1,-2 The vanishing of s at the inner cylinder surface x = a < 1 and at the x outer cylinder surface x = 1 implies [V.22] f*(a) + (Aa - m L o _ 1/a) f(a) =0 and [V.23] f'(l) + (A - mo)"1) f(1) = 0 . If there is no inner cylinder, the proper boundary conditions are that the pressure and its derivatives are continuous on the rotation axis: [V.24] f'(0) =Q if m = 0 [V.25] :•• f (0) =0 if m 4 0 . The boundary conditions for axisymmetric modes are similar to those of the sound modes in Chapter IV, so the differential equation for f(x) was solved numerically in this case in the same way that the sound wave equations were solved.1 The five eigenvalues -X which are smallest in magnitude for this m = 0 case are given in Table VIII, for four values of a. A value of 2 x 1.2168 = 2.4336 for A was used (see equation [IV.13]). See Appendix A. 99 2 Table VIII Eigenvalues -X For Some Axisymmetric I n e r t i a l Modes a = 0.0 a = .25 a = .50 a = .75 -4.86 -4.86 -4.85 -4.84 , -20.4 -25.8 -47.7 -169. -56.0 -80.6 -169. -654. -112. -171. -37 2. -1460 -188. -298. -656. -2590 From Table VIII i t can be seen that there i s no damping or growth of these modes i n the i n v i s c i d approximation. The eigenfunctions f f o r the four eigenvalues of smallest magnitude are plotted i n Figure 9. The eigenfunctions f are q u a l i t a t i v e l y l i k e l i n e a r combinations of Bessel functions J (Ax) and N (Ax), which are the m = 0 solutions of [V.16] o 0 for A = 0. I t follows that, as with the sound modes, buoyancy has no qual-i t a t i v e e f f e c t on axisymmetric i n e r t i a l modes. While t h i s conclusion i s based on the use of a c y l i n d r i c a l geometry, I can see no reason why i t should not also hold i n a s p h e r i c a l geometry. The v e l o c i t y f i e l d s corresponding to a given f can be formed by multiplying f by cos(kz) or sin(kz) and using [V.8] to f i n d s. One general feature of the v e l o c i t y f i e l d s i n a l l these axisymmetric modes i s that, since s, and s d i f f e r by a purely imaginary m u l t i p l i c a t i v e factor when cp X m = 0, the meridional and azimuthal flow are out of phase i n time by T T / 2 . The boundary conditions f o r the axisymmetric modes are independent of 2 the eigenvalue -A . However, the boundary conditions f o r the non-axisymmetr 2 2 modes do depend on -A . Furthermore, the dependence of -A i s i n general 100 non-linear. As a r e s u l t , the non-axisymmetric mode problem i s quite d i f f i c u l t to solve. There i s , however, one class of non-axisymmetric modes for which the problem i s t r a c t a b l e , and moreover, which i s of some i n t e r e s t . When k = 0 (no z-dependence of p^), the boundary conditions can be written i n a form 2 which i s l i n e a r i n -X : [V.26] -aAf'(a) - [ ( a A ) 2 + 2A]f(a) = - X 2 f ( a ) [V.27] . - A f ' ( l ) - [ A 2 + 2A]f(1) = -X 2f(1) . Because these boundary conditions take the same form as the eigenvalue equation [V.16], I t i s possible to use the methods outlined i n Appendix A f o r the s o l u t i o n of the problem. I t follows from [V.18] that f or these k = 0 modes, [V.28] OJ = Am/(X2 - 2A) or 0 . One branch corresponds to steady ( t o = 0) flow, while the other branch corresponds to o s c i l l a t o r y (OJ ^  0) motion. In the mathematical l i m i t A 0, the motion i n the l a t t e r branch also becomes steady. In t h i s l i m i t (which is,not relevant to neutron s t a r s ) , a l l buoyancy forces vanish, and the Taylor-Proudman theorem (Acheson and Hide, 1973) holds. That theorem states that steady flow i n r o t a t i n g , non-magnetized constant density f l u i d s i s •'two-dimensional" i n the sense that 3s/3z = 0. This theorem can be derived by assuming that OJ = 0 and V^ o = fi, and taking the c u r l of [V.7]. When V^ li f fi, the theorem i s not necess a r i l y true because the buoyancy o force -PcViJ>o gives r i s e to a torque <* Vp^ x which vanishes i f and only i f p' i s constant on surfaces of constant . Since 3s/9z i s * Vp' x v\jj c o c o when OJ = 0, 3s/3z need not vanish. In p a r t i c u l a r , i t does not vanish f or non-axisymmetric modes. 101 Conversely, non-axisymmetric two-dimensional (k =0) inertial modes need not necessarily be steady, as shown by [V.28]. The physical reason for this is the torque exerted on the fluid by the buoyancy force- Because of the importance of the buoyancy force for the non-axisymmetric oscillatory two-dimensional inertial modes, these modes, will be called buoyancy-inertial modes. To my knowledge, these modes have not been, previously discussed in the literature. The eigenfrequencies for some m = 1, 2, and 3 buoyancy-inertial modes are given in Table IX. The eigenfrequencies are a l l real, so these modes are a l l stable and undamped. Note that a l l the eigenfrequencies satisfy |to| < 1, or in ordinary units, |co| < 2fi Q, just as for ordinary inertial waves. However, the highest frequency buoyancy-inertial mode has a period which is more than five times longer than that of the highest frequency axisymmetric inertial mode. The motion, in the buoyancy-inertial modes is such that vanishes, i.e., the motion is everywhere parallel to the equator, since 9s/3z =0 and s vanishes at z = 1 1. The s field associated with these waves can z be calculated from [V.8].. At a given time, i t is of the form [V.29] s a [ ( f + Axf) - uf 1 - f] cos(m<j>) [V.30] S l « [af 1(f' + Axf) - - f] sin(m(j)) . cp x These fields for the highest frequency m - 1, 2, and 3 modes are plotted in Figures 10, 11, and 12, for a = 0 and a = .25. The cellular flow pattern arises from the e form for p^; the number of cells is 2m. The buoyancy force can have a qualitatively important effect on m 4 0 modes, in the sense that motion which would be steady in the absence of i t is oscillatory in its presence. While this conclusion was arrived at on the basis of a simplified cylindrical geometry, i t is true even in the Table IX Eigenfrequencies Of Some Buoyancy-Inertial Modes m a = 0.0 a = .25 a = .50 a = .75 .139 .107 .055 . .015 .047 .032 .015 .004 1 .023 .015 . 007 .002 .014 .009 .004 .001 .009 .006 .003 .001 .166 .155 .098 .029 .066 .056 .029 .008 2 .036 .028 .013 .003 .022 .017 .008 .002 .015 .011 .005 .001 .167 .164 .124 .044 .075 .070 .041 .011 3 .043 .038 .020 .005 .028 .023 .011 .003 .020 .016 .007 .002 103 s p h e r i c a l case. This can be seen most e a s i l y i n equation [V.10], which applies to any geometry. If the flow i s axisymmetric (m = 0) or i f the buoyancy forces vanish (3iJ>o/3x = 0), then [V.10] i s l i n e a r i n co , as a . r e s u l t of which the two branches of o s c i l l a t i o n are the same up to a change i n the sign of to. If the flow i s non-axisymmetric (ra ^  0) and the buoyancy force perpendicular to the r o t a t i o n axis does not vanish everywhere (9*|W8x f 0) , then there i s an a d d i t i o n a l non-zero term a co 1 i n [V.10]. As a r e s u l t of th i s term, the two branches of o s c i l l a t i o n can have quite d i f f e r e n t eigenmodes and eigenfrequencies. I n e r t i a l Flow In the f i r s t section of t h i s chapter some of the free i n e r t i a l modes of o s c i l l a t i o n of the charged f l u i d were investigated. In the present section, the charged f l u i d flow which develops from chosen i n i t i a l conditions i s examined. Since the f l u i d motion w i l l be followed only over the time scales of i n e r t i a l waves, I c a l l t h i s flow i n e r t i a l flow. In. order to f i l t e r put sound waves from the solutions to the equations, the a n e l a s t i c condition [V.31] V • (p v ) = 0 c c i s used. The s i m p l i f i e d momentum equation i s [V.32] 3(p v )/3t = -Vp' - p'V> - p vV - V - 2Q x p v + Q . c c r c c o c c c In equation [V.32], Q represents any unbalanced body forces which are roughly constant over these time scales. Such forces might occur, f o r example, i f the magnetic body force and the quasi-steady Poincare' force are not balanced at the i n i t i a l time. I t i s convenient to include such a term i f one i s interested i n determining the response of the charged f l u i d to an applied force. , Taking the divergence of [V.32] and using [V.31] produces an equation for p^ which i s e l l i p t i c , rather than hyperbolic, because of the neglect of sound waves: [V.33] V 2p' + V> . • Vp' + V2if/ p' = -v" • [2^ x p v - Q + p" vV + V • *a] c o c o c c c c (cf. equation [ I I I . 7 ] ) . Equation [V.33], when combined with the momentum equation, ensures that the f l u i d flow s a t i s f i e s the a n e l a s t i c condition for a l l times t > 0 i f i t i s s a t i s f i e d at t = 0. As an a l t e r n a t i v e to [V.31-32], i t i s possible to solve [V.32-33] i f [V.31] i s used as an i n i t i a l constraint on v . • c The g r a v i t a t i o n a l p o t e n t i a l perturbation s a t i s f i e s [V.34] v 2 * ' = ( 4 T T G / S 2 )p' . C O c Equation [V.34] and the pressure equation [V.33] form a pair of coupled e l l i p t i c equations. If the v e l o c i t y f i e l d and Q are known at a p a r t i c u l a r time, [V.33-34] allow the pressure and g r a v i t a t i o n a l p o t e n t i a l perturbations at that time to be calculated. The proper boundary conditions on the v e l o c i t y are the no-sli p conditions [11.87-88], since v i s c o s i t y i s taken into account i n [V.32]. For reasons of numerical s i m p l i c i t y , an i s o t r o p i c viscous tensor was used i n place of an anisotropic tensor which depends on the d i r e c t i o n of the magnetic f i e l d . The use of an i s o t r o p i c tensor i s ju s t i f i a b l e , b e c a u s e the Ekman number i s so small compared to 1 that v i s c o s i t y , whether i s o t r o p i c or anisotropic, does not have a very important e f f e c t on the f l u i d flow. 105 The magnetic induction equation i s not included i n the equations [V.31-34] which are to be solved because the changes i n B over the time scale fi. 1 are • o so extremely small that time v a r i a t i o n s i n the magnetic body force are neg-l i g i b l e , as was: shown i n Chapter I I I . However, changes i n B can cause changes i n the magnetic torques on the s o l i d parts of the star (see [II.95]). For the induction equation to be decoupled from the other equations, i t i s also necessary that these v a r i a t i o n s i n the magnetic torque be small compared to the viscous torque. From equations [11.78, 94, and 95], the r a t i o of the v a r i a t i o n i n the magnetic torque N^ to the viscous torque N' i s o V [V.35] N B / N V = V A T / V c ' for f l u i d motion with time scale T. For i n e r t i a l flow, T - fi \ so o [V.36] N'/N' « V 2 / ( f i o V c ) ~- 10" 5 B ^ fi'1 T 2 . Thus, the v a r i a t i o n s i n the magnetic torque are r e l a t i v e l y n e g l i g i b l e and the induction equation may be completely dropped unless B i s quite large compared to I Q 1 1 gauss or i f the neutron star i s so young that i t s temper-ature i s much higher than 10^ K. I t i s assumed i n the d e t a i l e d c a l c u l a t i o n s i n t h i s chapter than N' << NA, so the induction equation i s not solved. I t i s necessary to make d e f i n i t e choices for the values of the inde-pendent parameters of the problem, such as p (0), fi and T , i n order to o o o perform d e t a i l e d c a l c u l a t i o n s such as those reported below. Since no one set of parameters i s t y p i c a l of neutron stars as a c l a s s , any p a r t i c u l a r set may be c r i t i c i z e d as being somewhat a r b i t r a r y . However, the use of s c a l i n g laws such as those mentioned i n Chapter II allows one to generate a large class of approximate solutions from any p a r t i c u l a r s o l u t i o n . 106 The equilibrium model used i n the ca l c u l a t i o n s below was the one calculated i n Chapter II - see Table V and Figure 4. The value chosen for 11 2 -1 the kinematic v i s c o s i t y was v = 2 x 10 cm sec , whose associated c Ekman number f o r X - R i s so small compared to 1 (see [11.63]) that the re s u l t i n g f l u i d flow should be almost completely independent of the precise value of v c > The temperature corresponding to t h i s value of i s T^ - .4 . F i r s t Flow Simulation The charged f l u i d motion was simulated on a computer by solving very accurate f i n i t e d ifference approximations to the momentum equations and the coupled p^, $' equations. This was done as f o l l o w s . 1 At the i n i t i a l time t = 0, a v e l o c i t y f i e l d v (t = 0) was chosen which s a t i s f i e s the c anel a s t i c condition (or, more p r e c i s e l y , the f i n i t e d i f f e r e n c e approximation to i t ) . The coupled p^, $' equations were solved for p ^ ( t = 0) and * ' ( t = 0), which were then used i n the momentum equations to get 3 ( p c V c)/3t(t = 0 ) . This time.derivative and v c ( t = 0) were used to generate an approximation to V £ ( t = At > 0), where At i s a time i n t e r v a l small compared to a l l charac-t e r i s t i c time scales i n the problem. The procedure was then repeated at times t = At, 2At, etc. Each phy s i c a l v a r i a b l e was defined on one of four staggered grids which are uniform i n sp h e r i c a l coordinates . These grids are displayed i n Figure 13. The s l i g h t non-sphericities i n the s o l i d boundaries were neglected. It can be seen from the fi g u r e that the grids are coarsest near the crus t . Because a sp h e r i c a l coordinate system becomes singular at the o r i g i n , i t was numerically necessary to assume that there i s a s o l i d c e n t r a l core i n order See Appendix B f o r further d e t a i l s of the numerical methods used. to avoid the s i n g u l a r i t y . A very small core of radius Rg = 1 km was chosen. The moment of i n e r t i a of t h i s s p h e r i c a l core was determined to 40 2 be 2.48 x 10 gm cm by numerical i n t e g r a t i o n of the model i n Table V. Because of computer storage r e s t r i c t i o n s , only axisymmetric f l u i d flows could be simulated, i n which a l l quantities such as v , v., v., etc are r - o <p independent of azimuthal angle <p. Two of the computer simulations which were performed are reported here. The f i r s t simulation was performed i n order to see how the f l u i d responds to a suddenly applied force. This force was taken to be the Poincare" force due to c r u s t a l slowing down: [V.37] Q(t) = - P c & 0 x r f o r t - 0. The c r u s t a l deceleration ft. chosen was that observed for the Crab (see " o Table I I ) . For s i m p l i c i t y , the i n i t i a l charged f l u i d v e l o c i t y was assumed to vanish i n the co-rotating frame: [V.38] v (t = 0) = fj . c Due to the constant Poincare" force i n the -Hp d i r e c t i o n , v, i n i t i a l l y <P increased l i n e a r l y with time, as did.the perturbation ft' i n the c r u s t a l a c c e l e r a t i o n , which depends only on the <p component of.v . One might expect that the C o r i o l i s force would cause the meridional flow to grow i n i t i a l l y q u a d r a t i c a l l y with time, and become comparable i n magnitude with v^ a f t e r about one r o t a t i o n . This does not happen. A f t e r one r o t a t i o n , |v^ g l / l ^ l -3 i s only about 10 The reason for the r e l a t i v e l y small meridional flow i s that pressure gradients are f a i r l y successful i n balancing the C o r i o l i s force i n meridional planes. This can be seen by considering a model i n v i s c i d , constant density f l u i d with no g r a v i t a t i o n a l f i e l d present. The equations describing the 108 acceleration of such a f l u i d by the Poincare force are [V.39] 3v /3t = -3P'/3z z [V.40] 3v /3t = -3P'/3x + 2ft v, X o <p [V.41] 3v A/3t = -ft x - 2ft v (j) o o x where P' i s the perturbation i n the kinematic pressure. The s o l u t i o n to these equations which s a t i s f i e s the i n i t i a l conditions [V.38] i s [V.42] v = v = 0 Z X [V.43] = -ft xt cp o [V.44] P' = -ft ft x 2 t . o o Thus i n t h i s model, pressure gradients are completely e f f e c t i v e i n suppressing meridional flow. The computer generated s o l u t i o n of [V.32-34] for a f t e r about two rotations i s displayed i n Figure 14. The v e l o c i t y f i e l d i n that f i g u r e i s s p l i t i n t o a meridional or p o l o i d a l part v^ and an azimuthal or t o r o i d a l part v^. The meridional flow i s represented by a vector p l o t ; the azimuthal flow i s represented by a contour p l o t . For ease of v i s u a l i z a t i o n , the vectors v^ are attached to a uniform 15 x 15 g r i d i n (z,x) coordinates. I t would have been extremely hard to draw and v i s u a l i z e the flow using the grids i n s p h e r i c a l coordinates, as can be seen from the high density of these g r i d points hear the core (see Figure 13). It i s necessary t o . i n t e r -polate the vectors v^ before they be displayed on a uniform (z,x) g r i d . However, no matter what g r i d v^ i s displayed on, i n t e r p o l a t i o n i s necessary because v r and v Q are computed on two d i f f e r e n t grids (see Appendix B). 109 The contour p l o t shows that the computed s o l u t i o n i s e s s e n t i a l l y r i g i d r o t a t i o n i n the +<j> d i r e c t i o n with an extremely large v e l o c i t y shear near the crust, p a r t i c u l a r l y near the equator. Af t e r only two rotation s , v i s c o s i t y has not yet had time to reduce the shear. Apart from the shear which a r i s e s from the presence of no - s l i p boundaries, the computed s o l u t i o n i s p r a c t i c a l l y i d e n t i c a l to the s o l u t i o n [V.42-44] of the simple model equations [V.39-41], which do not take boundary conditions into account. This agreement i s consistent with the soundness of the numerical method. Because the numerical c a l c u l a t i o n s included buoyancy whereas the simple model [V.39-41] does not, the agreement between the two i s also consistent with the conclusion reached e a r l i e r i n thi s chapter that the buoyancy force has l i t t l e q u a l i t a t i v e e f f e c t on axisymmetric i n e r t i a l flow. The undulations which are apparent i n the contour i i n e s near the crust i n Figure 14 a r i s e because the numerical g r i d i s coarsest at these large r a d i i (see Figure 13) and because c o i n c i d e n t a l l y the shear i n v^ i s extremely large i n the same region. These are the most extreme conditions possible. Thus the s i z e of these a r t i f i c i a l wiggles gives some i n d i c a t i o n of the maximum numerical errors which can possibly a r i s e i n the numerical simulations. The c a l c u l a t i o n was not ca r r i e d further than two rotations because i t was apparent that nothing of i n t e r e s t was going to happen for a very long time. The question of what happens to the charged f l u i d i f the Poincare force i s "turned on" suddenly i s returned to i n Chapter VI, where the flow i 8 i s followed f o r a time - 10 times longer than i t was here. 110 Second Flow Simulation A second simulation was performed to see what type of charged f l u i d flow develops from a simple i n i t i a l v e l o c i t y condition. In th i s simulation, i t was assumed that the quasi-steady Poincare force was balanced by magnetic body forces, as explained i n Chapter I I : [V.45] Q = fi . The i n i t i a l v e l o c i t y f i e l d was chosen to be a simple azimuthal flow which vanishes on the crust and core surfaces: [V.46] V x(r,e,t=0) « (r - R ) ( r - R ) sine , <p s c The r e s u l t i n g flow at eight d i f f e r e n t times i s given i n Figure 15. I t can be seen that because of the C o r i o l i s force, meridional flow develops out of the i n i t i a l l y purely azimuthal flow. The maximum bf |v | never exceeds about a t h i r d of m a x l v J , and the o s c i l l a t i o n s i n time of maxlv I ' <P ' . P have about the same amplitude as the o s c i l l a t i o n s of max|v^|. Thus, the flow which develops consists of two parts. The f i r s t part i s a roughly steady flow i n the +<p d i r e c t i o n . If the c a l c u l a t i o n had been c a r r i e d forward 2 for a time >> R /v^, v i s c o s i t y would have slowed t h i s steady flow. The second part of the flow i s o s c i l l a t o r y i n character, with the amplitude of the meridional and azimuthal o s c i l l a t i o n s equal to about a half of the magnitude of the steady azimuthal flow. At .15 of a r o t a t i o n period a f t e r t = 0, a large scale meridional flow has developed, consisting of a large s w i r l centered near r = .75Rc> 6 = 3ir/8, plus a small counterswirl near the core. The azimuthal flow i s developing a large shear near the crust. At .31 ro t a t i o n s , the counterswirl has grown. By .47 rotatio n s , the o r i g i n a l s w i r l i n the meridional flow has disappeared. I l l What remains i s the counterswirl, which i s now very large. The shear i n the azimuthal flow near the crust has begun to decrease, and by .78 r o t a t i o n s , the flow i s much l i k e i t was sh o r t l y a f t e r t = 0. The o s c i l l a t o r y part of the flow thus has a period of roughly .7 of a r o t a t i o n period. From the large s p a t i a l scale of the O s c i l l a t i o n , i t can be reasonably concluded that the o s c i l l a t i o n i s e i t h e r the fundamental mode or another one very close to i t i n frequency. The flow patterns at about 11 and 16 r o t a t i o n periods give no i n d i c a t i o n of the f l u i d switching to another mode. To show more c l e a r l y the gross features of the flow, the angular accelerations of the crust and the core due to f l u i d shear stress are plotted i n Figure 16. The accelerations of the s o l i d parts of the star have two components, one a steady p o s i t i v e a c c e l e r a t i o n due to the steady part of the flow i n the +<j> d i r e c t i o n , and the other an o s c i l l a t o r y a c c e l e r -ation due to the o s c i l l a t o r y part of the flow. Figure 17 shows a power spectrum analysis of the c r u s t a l a c c e l e r a t i o n graph. 1 There i s a sharp peak i n the power spectrum at a l i t t l e l e s s than .75 r o t a t i o n s , with a much smaller peak at about one r o t a t i o n . The o s c i l l a t o r y motion thus consists of the f l u i d o s c i l l a t i n g mainly i n a normal mode with a period of - .75 P, and to a smaller degree i n a normal mode with period - P. Both these periods > are - .5 P, which they must be because the modes are i n e r t i a l . One feature of the flow pattern which can be seen most predominantly a f t e r 16 rotations i s some sinuous r i p p l e s i n both the meridional stream l i n e s and the azimuthal contour l i n e s near the r o t a t i o n axis above the s o l i d core. This feature can be seen developing at e a r l i e r times and i s c l e a r l y associated with the presence of the core. A f t e r only .15 r o t a t i o n s , In the power spectrum a n a l y s i s , the routine FOURT from the U.B.C. Computing Centre was used to Fourier analyze the data. 112 one r i p p l e has formed j u s t above the core, and i s connected with the small counterswirl i n the meridional flow i n the same region. As time goes on, more r i p p l e s grow further and further up from the core. These r i p p l e s are not a numerical a r t i f a c t . There are a number of strong reasons to support t h i s statement. F i r s t l y , these r i p p l e s are larger and more pronounced than the r i p p l e s near the crust i n the f i r s t simulation (Figure 14) which d e f i n i t e l y are numerical a r t i f a c t s . Since the l a t t e r r i p p l e s a r i s e from a combination of the most severe conditions the computer program can handle, any larger r i p p l e s such as the ones under consideration should not be a r t i f a c t s . Secondly, these r i p p l e s do not a r i s e from a g r i d point to grid point v a r i a t i o n i n the numerical s o l u t i o n . Even on the scale that the s o l u t i o n has been pl o t t e d , the flow i s i n no sense turbulent; rather, the flow near the r o t a t i o n axis above the core i s a coherent but r i p p l y one. The scale of the r i p p l e s i s several times the distance between the g r i d points (see Figure 13). Thi r d l y , i f the r i p p l e s were i n any way the r e s u l t of too coarse a numerical g r i d , one might reasonably expect that t h i s would r e s u l t i n numerical i n s t a b i l i t y . There i s no i n d i c a t i o n whatsoever of any such i n s t a b i l i t y . F i n a l l y , and most importantly, the wavelength of the r i p p l e s i s almost p r e c i s e l y equal to the radius of the core, which i s i n no way simply related to the gr i d spacing. Since the r i p p l e s are not a numerical a r t i f a c t , what are they? Consider the graph of the core angular a c c e l e r a t i o n i n Figure 16. Neglecting the steady acceleration, we see that the core i s o s c i l l a t i n g about the r o t a t i o n axis. Because the moment of i n e r t i a of the core i s so much smaller than the c r u s t a l moment of i n e r t i a , the core angular a c c e l e r a t i o n i s about 113 a hundred times larger than the c r u s t a l angular a c c e l e r a t i o n . Thus, the core i s o s c i l l a t i n g at such a rate that i t could have a s i g n i f i c a n t e f f e c t on the f l u i d flow. What i s the r e s u l t of moving a s o l i d body such as the core i n an enclosed r o t a t i n g f l u i d such as the charged f l u i d ? If the body i s moving slowly or i s o s c i l l a t i n g at a rate much slower than 2ft, the f l u i d which l i e s above the body p a r a l l e l to the r o t a t i o n axis moves r i g i d l y with the body and forms a Taylor-Proudman column (Greenspan, 1968). There can be large shears at the Interface between the f l u i d column and the rest of the f l u i d . If the frequency of o s c i l l a t i o n of the s o l i d body i s comparable to but not greater than 2ft, a columnar structure i s also formed. However, the column i s not t r a n s l a t i o n a l l y i n v a r i a n t along the r o t a t i o n axis as i t i s when OJ « 2ft ; structure appears i n the column which p e r i o d i c a l l y repeats i t s e l f along the r o t a t i o n axis. The r i p p l y f l u i d flow above the o s c i l l a t i n g s o l i d core i s very much l i k e the preceding d e s c r i p t i o n of Taylor-Proudman columnar flow f o r the case i n which ft i s comparable to 2ft. I conclude that the o s c i l l a t i o n of the s o l i d core has set up a Taylor-Proudman column, which has r i p p l y structures i n i t because the core i s o s c i l l a t i n g at frequencies comparable to 2ft. 114 Figure 9. Radial parts of some axisymmetric i n e r t i a l modes. The upper pi c t u r e shows the eigenfunctions f which s a t i s f y [V.16] f o r 2 the four a = 0 eigenvalues -X which are smallest i n s i z e (see Table V I I I ) . The lower p i c t u r e i s s i m i l a r , except a - .25. 116 Figure 10. F l u i d flow i n the highest frequency m = 1 buoyancy-inertial mode. The flow f i e l d s", which i s rel a t e d to the flow v e l o c i t y by s = P c o v c > i s displayed at an instant of time i n a vector p l o t . The r o t a t i o n axis points out of the paper. The z -> component of s vanishes. What i s shown i s (s , s , ) , which x cp i s independent of z. The length of each l i n e segment i s proportional to the magnitude of s, and i s oriented along s" I n the upper p i c t u r e , a = 0 (no inner c y l i n d e r ) ; i n the lower, a = .25. Similar comments apply to the next two f i g u r e s . 1 1 8 Figure 1 1 . F l u i d flow i n the highest frequency m = 2 buoyancy-inertial mode. 119 120 Figure 12. F l u i d flow i n the highest frequency m = 3 buoyancy-inertial mode. 121 122 Figure 13. Numerical grids used i n the i n e r t i a l flow simulation. The r o t a t i o n axis i s v e r t i c a l ; the equator i s h o r i z o n t a l . Both are indicated by heavy l i n e s , as are the s o l i d core and s o l i d crust surfaces. Each p h y s i c a l v a r i a b l e i s defined on one of the four overlapping g r i d s . See Appendix B for a d e t a i l e d d i s c u s s i o n of the g r i d s . 123 GRID LEGEND SOLID CORE RADIUS o + ( i , j + l / 2 ) • ( i + 1/2 , j ) • ( i + 1/2 , j + 1/2) INNER CRUSTAL RADIUS 124 Figure 14. I n e r t i a l flow a f t e r 1.98 rotations i n the f i r s t simulation. The v e l o c i t y f i e l d v^ i s s p l i t into a meridional or p o l o i d a l part Vp, which i s represented i n a vector p l o t , and an azimuthal part v^, which i s represented i n a contour p l o t . The contour l a b e l i n g conventions are those of Figure 5. 126 Figure 15. I n e r t i a l flow at various times i n the second simulation. This f i g u r e extends over several pages. Each p i c t u r e i s labeled by the time i n r o t a t i o n periods at which the flow i s represented. N 3 VO 130 132 MAX V P = 0.90 x 10 cm/sec MAX = 0.95 x I0" 6cm/sec U> LJ MAX V P = 0.20 x I0" 6cm/sec MAX V, = 0.77 x I0" 6cm/sec LO -IN 135 Figure 16. Core and crust angular a c c e l e r a t i o n due to shear stress i n second i n e r t i a l flow simulation. Note that the core ac c e l e r a t i o n i s much larger than the crust a c c e l e r a t i o n , and that each graph consists of a steady component plus a smaller o s c i l l a t o r y component. 136 CORE ANGULAR ACCELERATION (IO 1 0 rod/sec 2) CRUST ANGULAR ACCELERATION (IO"12 rod/sec2) 137 Figure 17. Power spectrum analysis of crust angular a c c e l e r a t i o n i n second i n e r t i a l flow simulation. 0.12 CVJ b 0.08 i UJ ^0.04 CL 0 1 1 0 0.25 0.50 139 CHAPTER VI: MAGNETOSTROPHIC FLOW It was shown in Chapter III that the longest period oscillations in the charged fluid are hydromagnetic-inertial waves. This chapter is concerned with the dynamics of the charged fluid over time scales comparable to the periods of these waves. Acheson and Hide (1973) call fluid motion over these time scales magnetostrophic flow. The simplified equations over these long time scales were derived in Chapter III. They are the induction equation ''-»••.' [VI. 1] f f = V x (v c x B) , the Maxwell equation [VI.2] V • B = 0 , the anelastic condition [VI.3] V • (p cv c) = 0 , and the momentum equation [VI.4] 5 = -2« x p v - ~ [ t x (V x B) - B x (V x B ) ] - Vp' - p'V> c c 4ir o o c c o The gravitational term -p^JQ 1 is not included in [VI.4]. This term was included in the inertial flow simulations of the previous chapter. However, i t had l i t t l e effect on the flow since as I showed in Chapter II, it is smaller than the buoyancy term -p^Ji by a factor of about P Q / P c o ~ 100. There is no viscous term in [VI.4]. I showed in Chapter II that i t is very much smaller than the Coriolis term, unless T_ << 1 or A , << 1. 7 o Furthermore, over time scales T , viscosity has a negligible effect on the HI flow in the vicinity of solid boundaries, unless the normal component of B at the boundary is very much smaller than its.tangential component. Since i t would be extremely difficult to take into account viscous effects in the type of computer simulations reported in this chapter, the computations are restricted to conditions in which such effects are negligible. A superfluid drag term is not present in.[VI.4]. In Chapter III, I showed that the damping effect of superfluid drag on hydromagnetic-inertial waves is generally very much less important than the damping effect of viscosity. Since the-latter is very small for a wide range of physical conditions (see [III.27]) and has been neglected here, i t is consistent to neglect superfluid drag. The quasi-^steady Poincare" term does not appear in equation [VI.4] becaus I assume that the B field in [VI.4] is the sum of the magnetohydrostatic field and the small perturbation required to balance ~ P C ^ 0 x r • The fluc-tuations in the crust angular velocity induced by fluid motion are not taken into account in [VI.4]. For the particular fluid flow simulation reported in this chapter, i t will be shown that the fluctuating Poincare force is relatively small. In general i t need not be small. The appropriate boundary conditions on are no-slip conditions (see [11.89 and 92]), since I have restricted attention in this chapter to magnetic fields which generally have non-zero normal components near solid boundaries. Thus [VI.5] v J E =3 . As explained in Chapter II, the boundary condition on B in the infinite electrical conductivity limit is that its normal component is static as viewed in a frame of reference co-rotating with the boundary. In terms of the perturbation B' in the magnetic field, this boundary condition is 141 [VI.6] B' • n|z = 0 . In the inertia! flow simulations i t <was necessary to assume that there was a solid central core because of numerical difficulties associated with the choice of a spherical coordinate system. Since the calculations in the present chapter were performed in cylindrical coordinates, there is no such necessity here. Accordingly, i t was assumed that there is no solid central core. Thus, [VI.5-6] apply at the inner crustal surface. As explained in Chapter II, the crust is subject to a fluctuating magnetic torque N' and to a torque N' due to viscous shear stress, both of which arise from the charged fluid motion. The general formula for the relative size of N' and N' is given by [V.35]. Combining that result with the estimate [III.20] o V for the time scale T „ t of hydromagnetic-inertial waves produces H i [VI.7] N^/Nl » 103 n2 A* T 2 . Thus, unless T-. << 1 or A, << 1, the fluctuations in the magnetic torque on the crust are larger than the viscous torque on the crust. The torques on the crust arising from the fluid flow change the angular velocity of the crust. Can the quasi-sinusoidal timing residuals observed in the Crab (Boynton et al., 1972) be due to this effect? It seems unlikely. The period of the lowest frequency hydromagnetic-inertial mode would have to be longer than the length of the data train for the Crab, about five years. It follows from [III.20] that for the period to be five years or longer the mean interior magnetic field strength would have to be v 2 x 10^ gauss, which is somewhat low compared to the value expected on the basis of the naive flux conservation argument (see Chapter I). Furthermore, one must find a process for exciting hydromagnetic-inertial waves in such a way that the predicted timing power spectrum tP(m) has the -4 observed form P^ co with the correct power coefficient (Boynton et al., 1972). 142 I have not been able to think of any excitation process which will do this. One can probably f i t the observed ^ ( t u ) by choosing the amplitudes of the hydromagnetic-inertial modes properly. However, the number of arbitrary parameters in such a procedure is equal to the number of normal modes, whereas a frequency walk model has only one parameter, namely P q . Thus, a frequency walk model is far less arbitrary than a hydromagnetic-inertial mode model. Equations For The Velocity Field In magnetostrophic flow, the fluid instantaneously adjusts so that the total force on any fluid element vanishes (see [VI.4]). This adjustment is analogous to the instantaneous adjustment of the pressure in incompressible flow in such a way that the incompressibility condition is satisfied. The adjustment, times in both cases for realistic flow are of course not zero, but are much smaller than the time scales of interest. The instantaneous adjustment of the velocity and pressure in magneto-strophic flow is reflected in the mathematical structure of the magnetostrophic equations [VI.1-4]. At a given i n i t i a l time t = 0, p^ and V £ cannot be chosen arbitrarily, but must satisfy the constraint equations [VI.3-4], just as B must satisfy [VI.2]. The induction equation [VI.1] can be used to calculate B at the next "instant" t = At. Equations [VI.3-4] must then be used to determine v £(t = At), and so on. In this dynamical problem, i t is necessary to be able to solve [VI.3-4] for the velocity field, given an arbitrary Lorentz force. For that purpose, it is useful to introduce the variables [VI.8] u = 2Slp(yc , i n terms o f w h i c h [ V I . 3 - 4 ] become [V I .9] V" • u = 0 [V I .10] '• • ft « x u - (v" + v\|> )p» + ft t where ft i s t he body f o r c e : [V I . 11] ft > - 7 H B x (V" x B) - B x (V x B ) ] . HV O O E q u a t i o n s [V I . 9 -10 ] must be s o l v e d f o r u , g i v e n ft. Because o f computer s t o r a g e r e s t r i c t i o n s , o n l y a x i s y m m e t r i c f l u i d f l o w s c o u l d be s i m u l a t e d . F o r t he same r e a s o n , a t t e n t i o n was r e s t r i c t e d to f l o w s w h i c h p o s s e s s e q u a t o r i a l symmetry o f t he f o l l o w i n g t y p e i n c y l i n d r i c a l c o o r d i n a t e s (z,x,<p), where z = 0 i s t he e q u a t o r : [V I .12] v z ( - z ) = - v z ( z ) [V I .13] v x ( - z ) = v x ( z ) [V I .14] Because o f t he assumed ax i s ymmet ry , t he <p component o f [V I . 10] i s [V I .15] u = Q. . x <p E q u a t i o n [V I .15] e x p r e s s e s t he b a l a n c e i n the fp d i r e c t i o n between the C o r i o l i s f o r c e and the m a g n e t i c body f o r c e . The z component o f u can be d e t e r m i n e d f r o m the a n e l a s t i c c o n d i t i o n [ V I . 9 ] : [ V I - 1 6 ) • and f r o m t h e v a n i s h i n g of u a t z = 0 ( see [ V I . 1 2 ] ) . , 1 4 4 The pressure perturbation terms can be eliminated from [VI.1 0 ] by -> -> applying the operator "(V + ViJ> ) x" to them. The result is o [VI.17] (v" + V> ) x (-z x u) = -(v" + Vifi ) x Q . o o Note that although p^ does not appear in [VI.17], the effects of buoyancy are s t i l l taken into account through the terms involving v\ji . . The <p compo-nent of [VI.17] yields an equation for u,: [vi . 1 8 ] I — + — j U(j).- ^  + — j Qz - ^  + a7"j Qx • As an aid in selecting solutions to [VI.18], I require that u, vanishes when .-> Q vanishes. This eliminates the arbitrariness associated with the possible presence of geostrophic flow (Greenspan, 1968). This arbitrariness arises from dropping the 9(p c V c)/3t term in the momentum equation, and is analogous to the arbitrariness in the pressure solution in incompressible hydrodynamic flow. In incompressible flow, an arbitrary constant may be added to the pressure solution. There is a simple yet instructive example of a solution to [VI.15-18] for the case in which 5 q is a uniform magnetic field in the z direction, and in which the magnetic body force is linearized. For this case, the charged fluid velocity is given by B \ [VI.19] p v CO C \87Tft , O/ V x B which clearly satisfies the anelastic condition. The flow field [VI.19] has a number of unusual features. The fluid velocity is everywhere parallel to the current J ' . The linearized Coriolis force 2ft x p v and the o co c linearized magnetic body force - ^ Q x (V x exactly cancel so that in the linear approximation the pressure perturbation p^ vanishes everywhere. 145 These features are possible because B q and fiQ are parallel, and will probably not be present in the non-aligned case. -v ->-Note that for an arbitrary B', the charged fluid velocity v c given in [VI.19] does riot necessarily satisfy the no-slip boundary conditions [VI.5]. Thus there can be a boundary layer, whose existence seems to be in contra-diction with statements made in Chapter III that there are no oscillatory hydromagnetic-inertial boundary layers where B • n[^ ^ 0. However, there /is no contradiction because any such boundary layers are transient. An - . -y extremely large shear in v c near E will cause B' to grow in such a way that the shear in v^ is rapidly reduced. The Charged Fluid Response Time Problem One of the essential simplifying approximations which have been made in two-component models of post-glitch timing behaviour of pulsars (Shaham et al,, 1973; Baym et al., 1975) is that the charged fluid co-rotates with the crust following a glitch. In support of this approximation, i t has been argued (Baym et al., 1969c) that within the time i t takes ah Alfve*n wave to cross the star,, the charged fluid can adjust to the new angular velocity of the crust. Since this crossing time is small compared to the post-glitch relaxation time, i t is argued, one may safely neglect the internal motions in the charged fluid and simply assume i t co-rotates with the crust. However, I showed in Chapter III that unless the magnetic field is -1 -A so large that &2 is « 10 (i.e. that E, « 1 for A - R), the waves which couple strongly to the magnetic field are hydromagnetic-inertial waves, which behave like Alfven waves only i f A « R. Such very short wavelength Alfven waves are generally damped by viscosity in very much less than a second, and so cannot propagate across the star. Thus, the assumption which has been made in two-component models that the charged fluid as a whole responds to a disturbance in the Alfven wave crossing time is not well founded, unless B q is ^  10 gauss. What then is the charged fluid response time? By definition, i t is the time for the fluid as a whole to adjust after the magnetic field lines have been sheared. Thus, i t is probably the time scale of the longest wave-length oscillations that are strongly coupled to the magnetic field, i.e., roughly the period T„ of hydromagnetic-inertial modes with X R. This Hi period can be the order of a month, for example, if B^ is roughly - 1 within the star. The observed post-glitch relaxation time for the Crab is several days (Boynton et al., 1972). Thus, the charged fluid response time need not be small compared to the post-glitch relaxation time. The magnetostrophic equations can be used to determine the behaviour of the charged fluid in response to changes in the angular velocity of the crust. One possible computer simulation involves a sudden change in the angular velocity of the crust (a glitch). A second possible simulation involves a sudden change in the angular acceleration of the crust. In this chapter, the results of a computer simulation of the second kind are reported. The choice of this simulation over a glitch simulation was based on numerical convenience, and not because of a lack of interest in the details of charged fluid post-glitch flow. However, the charged fluid response time, which is of great interest, should be almost the same in the chosen simulation as in a glitch simulation. 147 The Magnetic Field Structure The charged fluid flow which develops following a disturbance depends on the general structure and strength of the magnetic field in the Interior. The choice of magnetic field structure and strength in any magnetostrophic calculations must be guided by a number of considerations including pulsar observations, theories of neutron star formation, magnetohydrodynamic stability and simplicity. One can learn very little about the interior magnetic field from pulsar observations. Its mean value < B > is probably << 10^ gauss: if < B > were o ° o ^ 10 gauss, the star would be so distorted that the main energy loss mechanism would be the emission of gravitational radiation, in disagreement with timing observations (Ruderman, 1972). As regards the field structure, > 10 I deduced from pulsar slowing down rates that if B q . ^ 10 gauss at the inner crustal surface, the magnetic field there is almost purely poloidal or almost purely toroidal (see [11.25]). Theories of neutron star formation are of little help in determining the interior field structure and strength. The value of <Bq> predicted by the simplistic flux conservation argument mentioned in Chapter I is uncertain by several orders of magnitude. In addition the flux conservation argument makes no prediction of the post-collapse magnetic field structure. Ruderman and Sutherland (1973) have attempted to remove some of the uncertainty regarding the strength of the pre-collapse magnetic field. They argue that the stellar matter which later collapses to form a neutron star forms a convective core of radius - 4 x 10^  cm and that the magnetic field energy there is equal to the energy of convective motions. This argument 9 leads to a pre-collapse B of about 3 x 10 gauss. If the naive flux conservation argument is then applied, final magnetic field strengths of - -• • ' 148' 12 about 4 x 10 gauss are predicted to occur. As noted in Chapter I, the flux conservation argument completely ignores the amplification of a magnetic f i e l d by thie shearing action of differential rotation during the collapse forming a rotating neutron star. However, the most serious problem in applying theories of neutron star formation to the question of the structure and strength of the magnetic f i e l d inside neutron stars arises from present ignorance of the magneto-hydrodynamics of neutron stars prior to crust formation. Such newly formed stars are probably subject to magnetohydrodynamic i n s t a b i l i t i e s with growth times much less than the crust formation time. As a result, the magnetic f i e l d structure and strength in a neutron star a few days old could very well bear l i t t l e resemblance to what they were at formation. Thus, on the basis of formation arguments alone, one cannot say anything with any degree of certainty about the structure of the magnetic f i e l d inside neutron stars. If the magnetic f i e l d structure inside a neutron star i s unstable, the magnetic flux might be quickly redistributed in roughly the time T : through-Hi out the interior so that the structure becomes stable. If this is not possible, which is the case i f the in s t a b i l i t y arises from the topology of the f i e l d lines (Markey and Tayler, 1973), f i n i t e amplitude oscillations w i l l occur. The topology of the f i e l d lines can be changed in the resistive 2 diffusion time - A. /n, which for very short wavelength ( A << R) i n s t a b i l i t i e s can be hundreds of years or less depending on the temperature. The ultimate result of these resistive changes w i l l probably be the alteration of the unstable part of the f i e l d structure to a stable form. Thus i t is reasonable to r e s t r i c t ones attention at f i r s t to magnetic fields which are magneto-hydr©dynamically stable. .^ -; 149 Unfortunately, l i t t l e work has been done on magnetic f i e l d s t a b i l i t y inside neutron stars (Vandakurov, 1970, 1972; Tayler, 1973b; Chanmugan, 1974). This work has primarily been concerned with convective inst a b i l i t y . Fujimoto and Murai (1972) have considered the buoyant instability of an isolated toroidal magnetic flux tube in a gas layer just outside the surface of a neutron star; Although the concept of magnetic buoyancy may be applied to the question of magnetic s t a b i l i t y inside neutron stars, the details of the work of Fujimoto and Murai (1972) on the exterior f i e l d are irrelevant to the present discussion. Convection occurs i f an element of f l u i d which rises adiabatically and which is in (material plus magnetic) pressure balance with i t s surroundings is lighter than them. The most recent and complete work on adiabatic convection in neutron stars i s by Chanmugan (1974), who assumed that the nuclear relaxation time x ^ is much longer thaii the growth time of magnetic i n s t a b i l i t i e s . Both he and Tayler (1973b) argued that this growth time < 15 is an A l f v i n wave crossing.time. If B ^10 gauss, the growth time is more likely to be the crossing time t D of a hydromagnetic-inertial wave. HI Comparing the estimate [11.41] for x with the estimate [III.20] for x , N HI we see that indeed x , >> x U T , in agreement with Chanmugan's assumption that N Hi x ^ is larger than instability growth times. Chanmugan (1974) concluded that B must be ^ lO 1^ gauss for magnetic buoyancy to be able to drive convection inside neutron stars. Thus, convective i n s t a b i l i t i e s are probably unimportant inside neutron stars. Non-convective magnetic i n s t a b i l i t i e s in non-rotating stars have been considered by Tayler (1973a) and Markey and Tayler (1973). The basic cause of such i n s t a b i l i t i e s i s the curvature of the f i e l d lines, from which energy can be extracted, causing small perturbations to grow. Tayler (1973a) considers t o r o i d a l magnetic f i e l d s and shows that i n the absence Of r o t a t i o n they may be unstable to the interchange of f l u x tubes i f there i s a non-zero current density on the axis. However, when ro t a t i o n i s present i t tends to suppress these i n s t a b i l i t i e s . Somewhat, s i m i l a r conclusions about the s t a b i l i t y of t o r o i d a l f i e l d s have been reached by Acheson and Hide (1973). They point Out that i n a r a p i d l y r o t a t i n g f l u i d the interchange of two material rings which conserves t h e i r angular momenta requires an input of energy greatly i n excess of the energy a v a i l a b l e from an unstable t o r o i d a l f i e l d structure. More p r e c i s e l y , a t o r o i d a l f i e l d configuration B^(x) i s stable to a x i -symmetric perturbations i f and only i f 4ft 2 , ' . ^ t V I - 2 0 ] IT ~ & ( V A / X > " ° ' which i s equation (4.33) of Acheson and Hide (1973). Thus, If the magneto-hydrostatic t o r o i d a l f i e l d v a ries on a length scale X, where [VI.21] X i 50 B n Q'1 cm, the structure i s stable to small axisymmetric perturbations. However, i n a non-axisymmetric disturbance, the t o r o i d a l f i e l d l i n e s are twisted and the r e s u l t i n g Lorentz force can o f f s e t the s t a b i l i z i n g influence of the C o r i o l i s force through a magnetostrophic balance. Markey and Tayler (1973) discuss the i n s t a b i l i t i e s of p o l o i d a l f i e l d s i n non-rotating s t a r s . Their main conclusion i s that i f there are any closed f i e l d l i n e s w i t h i n the s t a r , the f i e l d i s subject to i n s t a b i l i t i e s s i m i l a r to those i n a pinched gas discharge. Closed p o l o i d a l f i e l d l i n e s occur i f the f i e l d i s much stronger near the centre of the star than e l s e -where. In Chapter I I , I concluded that an axisymmetric f i e l d B whose p o l o i d a l part B - B^ cp contains closed f i e l d l i n e s cannot balance the 151 quasi-steady Poincare force. Thus i t seems u n l i k e l y that the p o l o i d a l part of B i s very much stronger i n the deep i n t e r i o r than near the crust. At present, the uncertainty about magnetohydrbdynamic s t a b i l i t y i n s i d e r o t a t i n g neutron stars i s great; the only magnetohydrostatic f i e l d which i s d e f i n i t e l y known to be stable i s a uniform f i e l d p a r a l l e l to the r o t a t i o n axis. I t i s stable because there are no magnetic pressure gradients and because there i s no curvature of the f i e l d l i n e s . F i n a l l y , s i m p l i c i t y should be a guide i n choosing a f i e l d structure to use i n the magnetostrophic simulation. The simplest magnetohydrostatic f i e l d i s the uniform one j u s t mentioned. If the f i e l d strength i s << 10"^ gauss, such a f i e l d i s compatible with pulsar observations. I t t r i v a l l y s a t i s f i e s the condition of being e i t h e r predominantly p o l o i d a l or predom-inan t l y t o r o i d a l near the inner crust. I t contains no closed f i e l d lines', and as j u s t mentioned, i s the only magnetohydrostatic f i e l d configuration i n s i d e r a p i d l y r o t a t i n g neutron st a r s which i s known to be s t a b l e . For these reasons, I chose a uniform magnetic f i e l d aligned with the r o t a t i o n axis as the basic f i e l d for the magnetostrophic simulation. To this basic f i e l d I added the small t o r o i d a l f i e l d displayed i n Figure 5, which i s required to balance the quasi-steady Pbincar^ force. The a c t u a l value of B q chosen was 1 0 1 1 gauss. This choice i s not r e s t r i c t i v e , because a whole cla s s of solutions f o r other B values can o be formed by using the s c a l i n g of v a r i a b l e s on B ^ indicated throughout the t h e s i s . The angular speed fi and uniform deceleration fi were chosen o % o to be those of the Crab (see Table I I ) . The r o t a t i n g hydrostatic equilibrium model was chosen to be the one calculated i n Chapter IT. 152 Flow Simulation The magnetostrophic flow which develops following a sudden change in the angular acceleration of the crust was simulated on a computer in order to see how quickly the charged fluid responds to the disturbance and to determine some of the flow characteristics. The magnetostrophic equations [VI.1, 15, 16, and 18] were replaced by finite difference equations. A l l physical variables were defined on a uniform 15 x 15 grid ih the section 2 2 2 z - 0, x - 0 Of the <p = 0 plane, with the region z + x > R£ excluded.1 The simulation ran for 125 days of simulated time. There are two possible descriptions of the i n i t i a l disturbance which are equivalent in the sense that the flow which follows is the same in both cases. In the first description, the magnetic field is uniform and the star is not decelerating prior to the in i t i a l time t =0. At t = 0 a steady external torque is suddenly applied to the crust, i.e., the Poincare force is "turned on". In the second description, the star is decelerating and the magnetic field B q is nearly uniform prior to t = 0. The non-uniformity in B q arises from the relatively small toroidal field required to balance the quasi-steady PoincarS force. This toroidal field was calculated in Chapter II and is displayed in Figure 5. At t = 0, the small balancing toroidal field is "turned off", that is, the i n i t i a l perturbation B ' ( t = 0) in the magnetic field is the negative of the field in Figure 5. This i n i t i a l perturbation B ' ( t = 0) clearly satisfies the Maxwell equation [VI.2] and the magnetic field boundary condition [VI.6]. non-linearity in the subsequent flow is of order B ' / B , which is o ee [11.59]). If terms of this size had been neglected compared to 1 Appendix C for further details of the grid and the numerical methods terms of order 1, which they were not i n the actual c a l c u l a t i o n , the r e l a t i o n -ship between vV_ and B' would be that given i n [VI. 19]. Thus to a f i r s t approximation, the properties of the flow f i e l d v c i n [VI.19] hold for the simulated flow. In p a r t i c u l a r , the pressure perturbation p^ ar i s e s only from the small n o n - l i n e a r i t i e s , and thus buoyancy has very l i t t l e e f f e c t on the f l u i d flow. Basic Time Scales Of The Flow In order to show more c l e a r l y some of the basic time scales of the flow, the perturbation i n the c r u s t a l angular a c c e l e r a t i o n due to viscous stress i s displayed i n Figure 18. An i s o t r o p i c viscous tensor was chosen, as i n 11 2 -1 Figure 16, and V £ was taken to be 2 x 10 cm sec , which corresponds to a temperature T_ - .4. Any other value of v could have been chosen, or the 7 c anisotropy i n the v i s c o s i t y taken into account, or the magnetic torque f l u c t u a t i o n s included, but none of these would have changed the time scales. I t w i l l be shown l a t e r that the most prominent mode i n the flow simula-t i o n has a wavelength X^ = .3. Using ti^ = 1.9 and T^ = .4, equation [VI.7] can be used to estimate that N'/N' i s - 20 for t h i s flow. The t o t a l pertur-B V bation ft' i n the.crust angular a c c e l e r a t i o n i s therefore - 20 times the -13 -2 acceleration shown i n Figure 18, i . e . , ft' i s - 4 x 10 sec , which i s -4 . about 10 times the s i z e of ft .. Thus, the f l u c t u a t i n g Poincard force -*•. -y -P cft x r i s r e l a t i v e l y n e g l i g i b l e , and i t was a good approximation not to include i t i n the f l u i d momentum equation [VI.4]. The most s t r i k i n g feature i n Figure 18 i s the o s c i l l a t i o n s i n ft'(t). There appears to be one predominent period of about a day. On closer inspection, i t can be seen that there are two quite d i f f e r e n t phases to ft'(t). The f i r s t phase, which I c a l l the charged f l u i d response phase, 154 lasts a few days after t = 0. It is distinguished by a positive time average of ft'. During this phase the charged fluid responds to the crustal slowing down in such a way as to oppose i t . The charged fluid response time in this simulation is thus a few days, which is comparable to the observed post-glitch relaxation time for the Crab pulsar (Boynton et al., 1972). The important point is not that the two times are almost equal, but that the charged fluid response time is not very much smaller than the post-glitch relaxation time. ' • • • * . The charged fluid response phase is followed by a phase in which fi' fluctuates about a mean of zero. In the latter phase the fluid oscillates about a state of co-rotation with the decelerating crust. The oscillations are undamped because viscosity and superfluid drag have been neglected in the fluid equations. A power spectrum analysis of the crustal acceleration graph in Figure 18 was performed in order to examine more quantitatively this periodicities in the flow.1 The resulting power spectrum is given in Figure 19, which shows a large number of well defined peaks. These peaks are not the result of fluctuations at a single point in the discrete Fourier analysis, because each is composed of at least several points (there; are 2,000 points in the discrete frequency domain). There can be no doubt that the peaks are not the result of numerical "leakage" from the main peak at a period of 1.2 days. If there were a significant amount of such leakage, one would expect i t to produce peaks in the frequency region > 25 rad/day, where there are no peaks. In fact, there should be no peaks in this region because waves of that high a frequency cannot propagate on the computational grid, whereas waves whose frequencies are lower than this can. 1 The program EOURT from the U.B.C. Computing Centre was used in the Fourier analysis. The original 50,000 data points were averaged in blocks of 25 prior to Fourier analysis. The power spectrum was smoothed with a running average of 3 points. I thank G. Fahlman for recommending this procedure to me. It is thus quite certain that these peaks are not art i facts of the Fourier analysis. The peaks are due to the osc i l la t ion of the charged f lu id in hydromagnetic-inertial modes with frequencies equal to the frequencies at the peak maxima. Some of the prominent modes have periods of about .30, .60, .96, 1.2, 2.4, 3.8, 4.9, 6.3, 10 and 13 days. These normal mode periods were computed for, values of = ^» ^2 = ^"* 9 , a n < * ^c = ^ a n" The periods for other magnetic f i e ld strengths, crustal angular velocit ies and inner crustal rad i i may be obtained by scaling, because i t follows from [III.20] -2 2 that a l l normal mode periods are « £1^  t R c C'^111)/8.9 ] . Details Of The Flow Pictures of the charged f lu id velocity V £ and the perturbation B' in the magnetic f i e l d at a number of different times are given in Figure 20, which extends over many pages. Each picture is labeled with the time in days at which the situation is depicted. The meridional and azimuthal parts of v £ are separately displayed, as are the poloidal and toroidal parts of B' It was previously noted that B'(t = 0) is the negative of the toroidal f i e ld depicted in Figure 5, so i t is unnecessary to provide a separate picture of i t in Figure 20. The i n i t i a l azimuthal flow v^Ct = 0) vanishes, so only the meridional component of v (t = 0) i s displayed. Since the magnetic body force Q at t = 0 is equal to ~ P C ^ 0 x r by the the choice of i n i t i a l conditions, i t follows from [VI.15] that [VI.22] v x ( t = 0) . o The z component of v ' ( t = 0) can be determined from [VI.16]. Because is 156 smaller at larger r a d i i , v (t = 0) i s larger near the crust than near the z centre. These features of v £ ( t = 0) are displayed i n the f i r s t picture of Figure 20. Because v c ( t = 0) does not satisfy the no-slip boundary condition, a thin transient boundary layer is formed at the crust * After .06 days, an azimuthal flow v^ has begun to develop through the action of the Coriolis force on the meridional flow v . This azimuthal flow is in such a direction as to oppose the crustal slowing down. At this time v, i s s t i l l small in comparison to the meridional flow. Because of the large shear near the crust in the i n i t i a l meridional flow, the developing azimuthal flow has a similar shear. The magnetic f i e l d perturbation has a non-zero poloidal part B*^  at .06 days. It can be seen that the component of B^  normal to the crust vanishes as required by the magnetic f i e l d boundary condition [ V I .6]. On the other hand, i t can be seen from the pictures of B' in Figure 20 that the tangential component of B at any point on the crust i s not fixed, because the boundary conditions in the limit a -*-'«> do not place any restrictions on n x B' ] ^. BV l i e s nearly parallel to the unperturbed magnetic f i e l d because 9B' 9v 9v . [VI.23] r — - B -r—^- - B -—— z , 9t o 9z o 9z ' where the dependence of v^(t = 0) only on x (see [VI.22]) has been used in the last step of [VI.23]. As a result, |9B'/°t| i s << |9B'/9t| for some X z time after t =0. Since B z = B x = 0 i n i t i a l l y , |Bx| is << |B^| during this time. After a quarter of a day, and B^  have grown larger. The largest growth i n v^ i s near the crust at the equator since the unperturbed magnetic f i e l d is tangent to the crust there and as a result large shears in the x direction cannot be effectively resisted there. 157 The r e g i o n i n w h i c h v , and B l have l a r g e s h e a r s i s l a r g e r than i t was a t .06 d a y s . T h e r e i s a s m a l l r e g i o n n e a r the r o t a t i o n a x i s i n wh i ch v , i s • <P n e g a t i v e . These f e a t u r e s a r e e a s i l y i n t e r p r e t e d i n terms o f the d i s p e r s i v e n a t u r e o f t he h y d r o m a g n e t i c - i n e r t i a l waves . The sudden on se t o f t he c r u s t a l d e c e l e r a t i o n g e n e r a t e s h y d r o m a g n e t i c - i n e r t i a l waves o f v a r i o u s w a v e l e n g t h s near t he c r u s t . Because the phase v e l o c i t y o f h y d r o m a g n e t i c - i n e r t i a l waves i s i n v e r s e l y p r o p o r t i o n a l t o t h e i r w a v e l e n g t h ( see [ I I I . 1 7 ] ) , the s h o r t e r w a v e l e n g t h d i s t u r b a n c e s t r a v e l more q u i c k l y than the l o n g e r w a v e l e n g t h ones . A s h o r t w a v e l e n g t h d i s t u r b a n c e i s r e s p o n s i b l e f o r t he s m a l l r e g i o n i n wh i ch v , < 0. By .5 d a y s , somewhat l o n g e r w a v e l e n g t h waves have r e a c h e d t h i s a r e a , and t h e r e g i o n has grown. In t he meant ime, the l o n g e s t w a v e l e n g t h waves a r e s t i l l n e a r t he c r u s t , g r a d u a l l y w i d e n i n g the shea r r e g i o n . By t h r e e q u a r t e r s o f a day , t he L o r e n t z f o r c e s have grown t o s u c h an e x t e n t t h a t i n t h e ' l a r g e shea r r e g i o n n e a r the c r u s t a t t he e q u a t o r , v^ has become n e g a t i v e , and ft' i s < 0. A f t e r one d a y , ft' i s s t i l l < 0 and the r e g i o n o f n e g a t i v e v^ near the r o t a t i o n a x i s has grown. The l i n e s o f c o n s t a n t v , and B ' a r e t e n d i n g to l i n e up a l o n g the d i r e c t i o n o f t he u n p e r t u r b e d m a g n e t i c f i e l d , so t h a t |av^/3z| and thus j 1 1 a r e as s m a l l as p o s s i b l e . By 1.25 d a y s , t he a z i m u t h a l f l o w nea r the c r u s t i s a g a i n i n t he + <f> d i r e c t i o n , and ft' i s > 0. The p o l o i d a l and t o r o i d a l components o f each o f v £ and B ' a r e by now comparab le i n m a g n i t u d e . A f t e r 1.5 d a y s , t h e cha r ged f l u i d i s n e a r i n g the end o f i t s r e s p o n s e phase bec au se the l o n g e r w a v e l e n g t h waves have begun t o r e a c h the c e n t r e . The s i t u a t i o n a t t he l a s t day o f the s i m u l a t i o n , day 125, d i s p l a y s c h a r a c t e r i s t i c s w h i c h were p r e s e n t i n t h e f i r s t two days o f t h e f l o w , i n c l u d i n g the s m a l l n e s s o f |B ' | compared to |B ' | and the a l i g n m e n t o f the c o n t o u r s o f e q u a l v ^ and B^ a l o n g t h e d i r e c t i o n o f B * q . The l e n g t h s c a l e 158 for changes in v c and B' is about Rc/3, so that the period of the oscillations 2 '" should be about (1/3) times the period for longest wavelength waves. This is the order of days, which is compatible with the calculated predominant period of 1.2 days. At 125 days, there are indications of waves with, wavelengths which are very small compared to R^ . Such small scale structures would not be present after 125 days if viscosity had been taken into account and if T^  were ^  1, because the damping times of such waves can be the order of a few days or (HI) less (see [111.27]). Since is an extremely sensitive function of A, the large scale structures present at 125 days would probably persist for several years. The results of the computer simulation can be summarized and interpreted as follows. The sudden onset of the crust deceleration shears the magnetic field lines near the crust, particularly near the equator; This generates a variety of hydromagnetic-inertial waves. Because the shorter wavelength waves travel more quickly, they reach the stellar interior after only a fraction of a day, at which time the longer wavelength waves are s t i l l near the crust, widening the shear region. After several days, the longer wave-length waves reach the interior. At this time, the charged fluid is oscil-lating about a state of rigid rotation with the decelerating crust. The oscillations persist for some time thereafter. 159 Figure 18. Perturbation i n crust angular a c c e l e r a t i o n due to viscous torque. The perturbations i n the crust angular a c c e l e r a t i o n due to v a r i a t i o n s i n the magnetic torque induced by the f l u i d flow are larger than the v i s c o s i t y induced perturbations. However, i t i s not the amplitude of the perturbations but only t h e i r character as a function of time that i s of i n t e r e s t here. There are two d i s t i n c t phases apparent i n s i ' ( t ) . One i s a charged f l u i d response phase l a s t i n g a few days i n which the time average of fi'(t) i s p o s i t i v e . This i s followed by a second phase, i n which the time average of fi'(t) i s zero. 161 Power spectrum analysis of crust angular a c c e l e r a t i o n data i n Figure 18. The zero of the v e r t i c a l scale i s a r b i t r a r y . Peaks i n the spectrum correspond to frequencies to of hydro-ma g n e t i c - i n e r t i a l normal modes of the charged f l u i d . The 2 -1 -2 frequencies to are °= tt^ [R (km)/8.9] , and these were computed f o r B.. - 1, ft = 1.9, and R (km)/8.9 = 1. The 11 z c lack of structure at frequencies above 25 rad/day i s a numerical a r t i f a c t due to the f i n i t e g r i d s i z e . 163 \ Figure 20. Magnetostrophic flow of charged f l u i d induced by c r u s t a l d eceleration. In these p i c t u r e s , the charged f l u i d v e l o c i t y v^ and the perturbation B' i n the magnetic f i e l d are displayed at various times. Both v"c and B*' are s p l i t i nto p o l o i d a l parts v , B*' and t o r o i d a l parts v., B', which are separately p p 9 <P displayed. The conventions f o r the vector p l o t s and contour p l o t s are as i n Chapter V. Time i s measured i n days. Note that the notation MAX B stands f o r the maximum value that P |B^| a t t a i n s at a given time. Similar notation i s used f o r v , v, and B'. The i n i t i a l v e l o c i t y f i e l d i s purely meridional, P 9 9 so v , ( t = 0) i s not displayed. The i n i t i a l B' vanishes, and 9. P B^Ct = 0) i s equal to the negative of the f i e l d i n Figure 5. M A X V p = 0 . 2 8 x I O * c m / s e c MAX V P = 0.28 x 10 cm/sec MAX V^ = 0.34 x 10 cm/sec ON MAX Bp = 4.38 x IO4 GAUSS MAX = 4.82 x IO5 GAUSS MAX V p = 0 .28 x 10 4 em/sec MAX V^ = 0.10 x 1 0 4 cm/sec <3\ 00 MAX Vp =0.29 x IO*"4 cm/sec | j \ TIME I I I I * I 1 4 1 * * 0.500 MAX V^ =0.13 x Id 4 cm/sec TIME = 0.500 c4 MAX V p = 0.29 x 10 cm/sec TIME MAX V^ = 0.15 x I0 4 cm/sec TIME = 0.751 MAX Bp = 1.93 x l O 5 GAUSS MAX Bx = 4.13 x I O 5 GAUSS ro 173 MAX Bp = 2.61 x IO 5 GAUSS MAX B, = 3.94 x IO5 GAUSS 175 MAX Bp = 3.30 x IO 5 GAUSS MAX B^ = 3.80 x IO 5 GAUSS MAX Bp = 4.00 x IO 5 GAUSS MAX B^, = 3.70 x IO 5 GAUSS 00 o 181 CONCLUDING REMARKS The conclusions I have reached are presented i n the form of answers to the sets of questions which I posed at the end of Chapter I. For ease of reference, the questions are repeated here. 1. What i s the range of v a l i d i t y of the f l u i d d e s c r i p t i o n of the charged p a r t i c l e s . i n the i n t e r i o r ? Within t h i s regime, to what extent are the usual equations of magnetohydrodynamics v a l i d ? The charged p a r t i c l e s may be treated as a f l u i d i f the length scale X -2 of a disturbance i s >> 3 T^ cm (the electron mean free path) and i f the -10 -2 time scale T i s >> 10 T^ sec (the electron mean-free time). Under such conditions, q u a s i - n e u t r a l i t y i s an excellent approximation. The magnetic f i e l d evolves i n time according to the i n f i n i t e conductivity induction equation. The creation or destruction of charged f l u i d by the acti o n of the weak i n t e r a c t i o n on density perturbations i s completely n e g l i g i b l e , and the usual mass conservation equation holds. Off-diagonal terms i n the material pressure tensor due to f i n i t e Larmor radius e f f e c t s are r e l a t i v e l y small. Other r e l a t i v e l y small terms i n the f l u i d momentum equation are the viscous term and the s e l f - g r a v i t a t i o n a l term. To a very good approximation, the heat f l u x equation i s decoupled from the other equations. To my knowledge, there has been no previous attempt to answer t h i s f i r s t set of questions. What types of small amplitude waves can e x i s t i n the charged f l u i d ? What can damp them? How much are they damped? How do the answers to these questions depend on the gross parameters of the problem, such as Q, B and T? ' ' 182 The answers to these questions are given i n Table X. The most s i g n i f i c a n t r e s u l t s are the long periods of hydrbmagnetic-iriertial waves and t h e i r extremely long damping times. Two a d d i t i o n a l damping mechanisms not l i s t e d i n the table are heat conduction and e l e c t r i c a l r e s i s t i v i t y , which are s i g n i f i c a n t only i f the star i s so young that T^ i s >> 1. The possible existence of i n e r t i a l waves and hydromagnetic-inertial waves ins i d e neutron stars was noted by Hide (1971). However, he did not consider the question of damping or investigate further the magnetohydro-dynamics of the charged f l u i d . 3. What kinds of mode-mode coupling can occur? Can t h i s weakly non-linear e f f e c t lead to heavy damping of a wave which i s only l i g h t l y damped i n the l i n e a r approximation? The second of these questions concerns hydromagnetic-inertial waves. If such a wave can resonantly couple to two other waves, one of which i s a sound wave or an i n e r t i a l wave, the p o s s i b i l i t y should be considered that energy can be transferred to the sound or i n e r t i a l wave and dissipated i n a time short compared to the period of the hydromagnetic-inertial wave. I found that the only such resonant coupling involves a hydromagnetic-i n e r t i a l wave and e i t h e r a p a i r of i n e r t i a l waves or a p a i r of sound waves. However, I showed that resonant energy transfer was not possible i n t h i s s i t u a t i o n . 183 Table X Charged F l u i d Waves And Their C h a r a c t e r i s t i c s TYPE PERIOD VISCOUS DAMPING TIME SUPERFLUID DRAG DAMPING TIME SOUND 10~ A X£ sec o 2 2 10 X, T^ o 7 sec days INERTIAL1" 1 r o t a t i o n 2 period 2 2 1 0 X 6 T 7 sec days HYDROMAGNETIC-INERTIAL1" B U X 6 " n 2 months 100 B~ 2 ft2 years 4 2 X6 T7 106 ft2 x2 v years NON-DIMENSIONALIZED VARIABLES NOTATION WAVELENGTH X, = X(cm)/10 6 CRUST ANGULAR SPEED Q2 = ft(rad/sec)/102 MAGNETIC FIELD B l l "-^(sauss)/^ 1 INTERIOR TEMPERATURE T ? = T(°K)/10 7 SUPERFLUID DRAG TIME* x 1 = T^day)/!. It i s assumed above that.B • ft. X, >> 10 11 2 6 TD d e f i n e d by;drag f o r c e / u n i t volume on charged f l u i d c s c D 184 4. What are some normal modes of oscillation of the charged fluid? I calculated some of the sound modes and their eigenfrequencies in the inviscid approximation. Neither the variations in the charged fluid sound speed nor the buoyancy force had any qualitative effect on the modes. I calculated some inertial eigenmodes and eigenfrequencies in a crude representation of the interior consisting of two co-axial rotating cylinders of finite length with a gravitational field perpendicular to the rotation axis. Viscosity was neglected. A l l axisymmetric modes were stable and relatively unaffected by the gravitational buoyancy force. I also calcu-lated a number of non-axisymmetric modes which are independent of the distance along the rotation axis. The frequencies of these buoyancy-inertial modes, as I call them, become zero if the gravitational field is turned off 5. How do charged fluid oscillations affect the crust? What is the relationship, i f any, between the charged fluid and the quasi-sinusoidal residuals? The motion of the charged fluid can change the angular velocity of the crust either through the viscous torque or through variations in the magnetic torque. If the fluid is oscillating at frequencies of inertial waves or higher, the viscous torque is the more effective of the two mechanisms. At hydromagnetic-inertial frequencies, the fluctuations in the magnetic torque are generally more effective. It appears unlikely that charged fluid oscillations can account for the residuals. A mean interior magnetic field of ^  2 x lO 1^ gauss would be required. Furthermore, there seems to be no natural way that hydro--4 magnetic-inertial modes can produce an to timing residual power spectrum. 185 6. How does the charged f l u i d respond under a v a r i e t y of i n i t i a l conditions and external forces? In p a r t i c u l a r , how does i t respond and how quickly does i t respond to changes i n the crust angular v e l o c i t y ? ; I simulated the charged f l u i d flow over two quite d i f f e r e n t time scales on a computer to see how the f l u i d responds. Two simulations were performed i n which the f l u i d motion was followed over the time scales of i n e r t i a l waves. In one of these simulations, the charged f l u i d i n i t i a l l y has a small azimuthal flow as seen i n a reference frame co-rotating with.the cru s t . A meridional flow develops through the action of the C o r i o l i s f orce, and the f l u i d begins to o s c i l l a t e p r i m a r i l y i n an i n e r t i a l mode with a period - .75 times the r o t a t i o n period. One i n t e r e s t i n g feature of the flow i s a type of Taylor-Proudman column p a r a l l e l to the r o t a t i o n axis above a s o l i d core which o s c i l l a t e s r e l a t i v e to the crust. I also simulated the f l u i d motion over the long time scales of hydro-magn e t i c - i n e r t i a l waves. I t i s only over these time scales that the magnetic f i e l d and the charged f l u i d are strongly coupled. P r i o r to the i n i t i a l time t = 0, the magnetic f i e l d is. a uniform 11 10 gauss f i e l d p a r a l l e l to the r o t a t i o n axis. The star i s not decelerating, At t = 0, a steady external torque i s applied to the crust so i t decelerates s t e a d i l y . The angular v e l o c i t y and deceleration chosen were those of the Crab pulsar-Although the simulation performed i s not a g l i t c h simulation, the response time of the charged f l u i d should be roughly the same i n both cases. It takes several days for the longer wavelength (A = R /3) hydromagnetic-i n e r t i a l waves formed at the crust near the equator to reach the centre of the s t a r . At the end of t h i s time, the charged f l u i d i s o s c i l l a t i n g about a state of co-rotation with the decelerating crust. The observed p o s t - g l i t c h 186 r e l a x a t i o n time f o r the Crab pulsar i s sdveral days. The major conclusion of my thesis i s that the charged f l u i d response time i s not necessa r i l y very much shorter than p o s t - g l i t c h r e l a x a t i o n times. In view of t h i s conclusion, simple two-component models of p o s t - g l i t c h behaviour which assume that the charged f l u i d response time i s much shorter than post-g l i t c h r e l a x a t i o n times need to be re-examined. 187 BIBLIOGRAPHY Abrikosov, A.A. and Khalatnikov, I.M. Sov. Phys. Usp. 1, 68 (1958) Acheson, D.J. and Hide, R. Rep. Prog. Phys. _36, 159 (1973) Bahcall, J.N. and Wolf, R.A. Phys. Rev. 140, B1452 (1965) Batchelor, G.K. Quart. J . Roy. Meteorol. Soc. 29, 224 (1953) Baym, G., Lamb, D.Q., and Lamb, F.K. "Influence of External Torques on the the Wobble Motion of Two-Component Neutron Stars", preprint (1975) Baym, G., Pethick, C , and Pines, D. Nature 224, 673 (1969a) . Baym, G., Pethick, C , and Pines, D. Nature 224, 674 (1969b) Baym, G., Pethick, C., and Pines, D. Nature 224, 872 (1969c) * Baym, G., Pethick, C , and Sutherland, P. Ap. J . 170, 299 (1971) Baym, G. and Pines, D. Ann. Phys. 66_, 816 (1971) Boynton, P.E., Groth, E.J., Hutchison, D.P., Nanos, J r . , G.P., Partridge, R.B., and Wilkinson, D.T. Ap. J. 175, 217 (1972) Canuto, V. and Chitre, S.M. Phys. Rev. D 9, 1587 (1974) Chanmugan, G. M.N.R.A.S. 169, 353 (1974) . ' Chapman, S. and Cowling, T.G. The Mathematical Theory Of Non-Uniform Gases, Cambridge U n i v e r s i t y Press (1939) Cohen, J.M. pp. 87-95, R e l a t i v i t y And G r a v i t a t i o n , Gordon and Breach (1971) Drake, F. Talk at 61"*1 Texas Symposium, New York (1972) Feibelmah, P.J. Phys. Rev. D 4_, 1589 (1971) Flowers, E. and Itoh, N. "Transport Properties of Dense Matter", preprint (1975) Fujimoto, M. and Murai, T. Publ. Astron. Soc. Japan 24, 269 (1972) Ginzberg, V.L. Sov. Phys. Uspekhi 14, 83 (1971) Greenspan, H.P. The Theory Of Rotating F l u i d s , Cambridge U n i v e r s i t y Press (1968) Greenstein, G. and McClintock, J.E. Science 185, 487 (1974) Groth, E.J. "Timing of the Crab Pulsar I I . Method of Analysis", preprint (1975a) 188 Groth, E.J. "Timing of the Crab Pulsar I I I . The Slowing Down and the Nature of the Random Process", pr e p r i n t (1975b) Harlow, F.H. and Welch, J.E. Phys. FI. 8_, 2182 (1965) Harrison, B.K., Thorne, K.S., Wakano, M., and Wheeler, J.A. G r a v i t a t i o n a l  Theory And G r a v i t a t i o n a l Collapse, Univ. of Chicago Press (1965) Heinzman, H., Kundt, W., and Schrufer, E. Astron. and Astrophys. 27, 45 (1973) Heinzman, H. and Nitsch, J . Astron. and Astrophys. 21, 291 (1972) Hide, R. Nature Phys. S c i . 229, 114 (1971) K r a l l , N.A. and T r i v e l p i e c e , A.W. P r i n c i p l e s Of Plasma Physics, McGraw-Hill (1973) Landau, L.D. and L i f s h i t z , E.M. Electrodynamics Of Continuous Media, Pergamon Press (1966) Landau, L.D. and L i f s h i t z , E.M. S t a t i s t i c a l Physics, Pergamon Press (1969) Le Blanc, J.M. and Wilson, J.P. Ap. J . 161, 541 (1970) Malkus, W.V.R. Mathematical Problems In The Geophysical Sciences, p. 207, Amer. Math. Soc. (1970) Manchester, R.N. and Taylor, J.H. Ap. J . 191, L63 (1974) Manchester, R.N., Taylor, J.H., and Van, Y.Y. Ap. J . 189, L119 (1974) Markey, P. and Tayler, R.J. M.N.R.A.S. 163, 77 (1973) Nelson, J . , H i l l s , R., Cudaback, D., and Wampler, J . Ap. J . 161, L235 (1970) Ostriker, J.P. and Gunn, J.E. Ap. J . 157, 1395 (1969) Pandharipande, V.R. and Smith, R.A. Nuc. Phys. A237, 507 (1975) Pippard, A.B. P h i l . Mag. 46, 1104 (1955) Reichley, P.E. Talk at 6 t h Texas Symposium, New York (1972) Roache, J.P. Computational F l u i d Dynamics, Hermosa Pu b l i c a t i o n s , Albuquerque New Mexico (1972) 189 Roberts, D.H. and Sturrock, P.A. Ap. J . 173, L33 (1972) Ruderman, M. Nature 223, 597 (19o9) Ruderman, M. Nature 225, 619 (1970) Ruderman, M. Ann. Rev. Astron. Astrophys. 10, 427 (1972) Ruderman, M. and Sutherland, P.G. Nature Phys. S c i . 246, 93 (1973) Ruderman, M. and Sutherland, P.G. Ap. J . 190, 137 (1974) Sagdeev, R.Z. and Galeev, A.A. Nonlinear Plasma Theory, W.A. Benjamin (1969) Shaham, J . , Pines, D., and Ruderman, M. Ann. of the N.Y. Acad. S c i . 224, 190 (1973) Tayler, R.J. M.N.R.A.S. 161, 365 (1973a) Tayler, R.J. M.N.R.A.S. 162, 17 (1973b) ter Haar, D. Physics Reports 3_, 57 (1972) Tsuruta, S. and Cameron, A.G.W. Nature 207, 364 (1965) Vandakurov, Yu.V. Ap. L e t t . 5., 267 (1970) Vandakurov, Yu.V. Sov. Astr. 16, 265 (1972). 1 9 0 APPENDIX A;, NUMERICAL METHODS FOR SOUND WAVES In b r i e f , the d i f f e r e n t i a l equation f o r was solved by approximating i t by a f i n i t e matrix eigenvalue problem. The i n t e r v a l [a,l] was divided into n + 1 equal i n t e r v a l s of length Ar: [A.l] ( i ) r = a + iAr so that r ( 1>,..., r ( n> are the i n t e r i o r points, and r ( o ) and r ( n + 1 > a r e the boundary points. The function f, and i t s f i r s t and second d e r i v a t i v e s were approximated as follows: [A.2] [A.3] [A.4] , f f ( r ( 1 ) ) - > f, ( 1> f ; ( r - o -> (f (i+l) _ f ( i - l ) ')/2Ar f ' ' ( r ( i > ) - > ( f J i + 1 > . 2 f i 1 > + f C i - D ) / ( A r ) 2 The d i f f e r e n t i a l equation f o r f^ then becomes a set of n l i n e a r equations i n n + 2 unknowns, the i t h equation r e l a t i n g " f ^ 1 " 1 ^ , f ^ and f . The boundary conditions were handled by solving t h e i r d i f f e r e n c e analogues f o r fj°^ and f j - n + 1 ) j a n ( j then s u b s t i t u t i n g these expressions t i l into the f i r s t and n equations. The r e s u l t i s a matrix eigenvalue problem: [A.5 ] A • F = -O/F , where A i s a non-symmetric, t r i - d i a g o n a l , n x n r e a l matrix, and _F i s the column vector (f . . •, f ^ ) • Equation [A.5] was solved with the use of routines FIGI2 and IMTQL2 of the subroutine package EISPACK.1 The c a l c u l a t i o n s were performed i n 1 Obtained from the U.B.C. Computing Centre. 191 double p r e c i s i o n on an IBM 370 machine. To test convergence of the answers, the equations were solved with n = 25, 50 and 100. The dependence of the lowest f i v e eigenf requencies of the 1, a = 0 case on n i s displayed i n Table X. I t can be seen from Table XI that the numerical method i s converging w e l l . The f i n a l r e s u l t s were a l l calculated with n = 100. Table XI Convergence Of Sound Wave Eigenfrequencies n = 25 n = 50 n = 100 1.750 1.742 1.738 5.018 4.966 4.937 7.888 7.823 7.781 10.65 10.59 10.55 13.33 13.32 13.28 APPENDIX B: NUMERICAL METHODS FOR INERTIAL FLOW 192 The equations were solved by replacing d e r i v a t i v e s by di f f e r e n c e s , and solving the r e s u l t i n g algebraic equations. S p a t i a l Mesh Systems A set of four staggered grids i n the upper right-hand quadrant of the z = 0 plane was used, s i m i l a r i n many respects to the MAC mesh system of Harlow and Welch (1965). These four grids can be generated by the product of a p a i r of r grids and a pa i r of 6 g r i d s : [B.l] r . = R g + ( i - |)Ar [B.2] r i + l / 2 = R s + ( i " [B.3] 6 j ' - (j - |)AG [B .4] 3 j + l / 2 = ( j " 1 ) A 9 where i - {1,...,^} , j = {1,...^} , R g i s the radius of the s o l i d c e n t r a l core, and [B.5] Ar = (R - R ) / ( N - 2) c s r [B.6] A6 = ( T T/2)/(N - 2) . 6 The four meshes may be labeled by ( i , j ) , ( i , j + -|), ( i + - | , j ) , and ( i + -T, j +4). In the ca l c u l a t i o n s reported i n t h i s thesis i both N and N 2 2 r x i were taken to be 15. A representation of the four grids i s given i n Figure 13. 193 The three components ^ v r* V0» v (j )) °f t n e v e l o c i t y were placed on the ( i ( i , j +-|"), and ( i , j ) meshes, re s p e c t i v e l y . The magnetic f i e l d components were s i m i l a r l y placed. The va r i a b l e s p' and $' were placed on c the ( i , j ) mesh. Of the stress tensor components, a , a . . , , and a,, were r r ; • 00 <pcp placed on the ( i , i ) mesh; a . on the ( i + i k i ) mesh; a.,'on the (i,1 +4) rep 2. . 0<p L mesh; and a on the ( i + \>i + ^ ) mesh. This choice of meshes and v a r i a b l e placements was designed so that a l l d e r i v a t i v e s can be replaced by second order accurate c e n t r a l d i f f e r e n c e s . For example, consider the reasons for the placement of p' r e l a t i v e to v c r and v •.: The r a d i a l component of the f l u i d momentum equation contains the terms [B-7] ( p c v r ) - - - ^ . With the chosen placement of v a r i a b l e s , [B.7] becomes C B - 8 ] ' . . I F : ^ c V i + i / 2 l j = "'- [ K h + i t i - to'chj'*'''- • I t can be seen from [B.8] that the placement of ( p ' ) . . between (v ). , *c i , j r x 1-1/2,.] 2 and ( v r ) ^ + ] l / 2 j produces a truncation error of 0((Ar) ) i n evaluating 9pV3r, 2 which i s four times smaller than the truncation error 0((2Ar) ) which would have ar i s e n i f only one g r i d had been used. S i m i l a r l y , the 6 component of the f l u i d momentum equation contains the terms ^ J£ (>cV = • * ' - 7 ^ " v ' which becomes [B.10] 8 a l ( p c v 6 > i , j + i / 2 - f ^ i . j + l ' < P c > i , ^ / ( r i A 0 ) ' I t i s apparent from [B.10] that the formal truncation error i n evaluating 3p'/39 i s four times smaller than i t would have been i f p' and v. had been c - c 0 placed on the same g r i d . Not only p', v and v„, but a l l v a r i a b l e s were placed so as to achieve c r 9 a f o u r - f o l d decrease i n s p a t i a l d e r i v a t i v e truncation errors compared to that achievable with a s i n g l e g r i d with the same N^, . In terms of truncation error, the staggered 15 x 15 grids used were equivalent to the use of a s i n g l e 30 x 30 g r i d , but with no more computer storage space than with a si n g l e .15 x 15 g r i d . This choice of staggered grids also allows boundary and symmetry conditions to be e a s i l y s a t i s f i e d . For example, consider the symmetry conditions on the equator 9 = TT/2. These are 3v [ B . l l ] jf- (0 = TT/2) .= 0 [B.12] . , v Q ( 9 = TT/2) = 0 [B.13] ? T ( 8 = * , 2 ) = 0 ' whose f i n i t e d i f f e r e n c e versions are simply [ B- 1 4 ] ^rVl / a . N - ( v r>i+l/2,N -1 0 o [B'15] ( V i , N - 1/2 = ° o Occasionally, variables are required at points where they have not been placed. In that case, averaging of nearby values was used. For example, we define 195 The Coupled Pressure And P o t e n t i a l Equations Given the v e l o c i t y f i e l d and the magnetic f i e l d at some time, the two coupled e l l i p t i c equations f o r p^ and $' were solved i n a " s e l f - c o n s i s t e n t " manner; that i s , a t r i a l s o l u t i o n f o r $' was assumed, and the equation for p' solved. These new values f o r p' were then used to c a l c u l a t e new c *c values, and so on, u n t i l convergence was reached. At each stage i n t h i s i t e r a t i v e process, i t i s necessary to solve an e l l i p t i c equation with given right-hand side. For example, the Poisson equation f o r i s [B.18] r ( r 3r ) + 2 . n 36 ( s i n 6 36 ) = ~T K ' 1_ 3 ', 2' 3$' 2 .3 r sin6 co which must be solved f o r a s s u m i n g that p' i s known. The f i n i t e d i f f e r e n c e c version of [B.18] i s [B.19] •r ± / i-1/2 (*') i - l , j + /Ar VZ 8 i n Bi + l / 2 , / A r \ ^ T A B J sine ^ +\TKe) 1+1/21 + ' r-2 S i n 9 i - l / 2 s i n 6". r. I r. ^ A 6 l I sine sin6 + - J - l/2 smt (*')• = 4uG(Ar) Z ^ 1»J which, along with f i n i t e d i f f e r e n c e forms f o r the $' boundary conditions, constitute a set of l i n e a r equations i n N N = 225 unknowns. The f i n i t e r e d i f f e r e n c e version of the e l l i p t i c equation f o r p' i s not given because i t C i s much more complicated than [B.19], 196 The use of successive overrelaxatibn (SOR) methods to solve e l l i p t i c equations l i k e [B.19] i s very time consuming. Since [B.18] has to be solved several times on each time step, the computing time and cost which would a r i s e i f SOR methods were used would be completely p r o h i b i t i v e . . I theref ore adapted the EVP method of Roache (1972) , which had been designed f o r the Poisson equation i n Cartesian coordinates, to the present problem. The EVP method has two great advantages over SOR methods. F i r s t l y , i t i s not I t e r a t i v e as SOR methods are. Thus, one does not need to worry about.whether or not the s o l u t i o n has converged. Secondly, i t i s about two orders of magnitude f a s t e r than SOR methods. The actual c a l c u l a t i o n s of the p^, $' solutions were performed i n double p r e c i s i o n on an IBM 370 computer. Time Di f f e r e n c i n g And Time Centering There are e s s e n t i a l l y two kinds of time d i f f e r e n c i n g : i m p l i c i t and e x p l i c i t . In an i m p l i c i t method, a set of coupled algebraic equations must be solved i n order to c a l c u l a t e v a r i a b l e s l i k e v at the advanced time. In c an e x p l i c i t method, the advanced values of v a r i a b l e s are given e x p l i c i t l y i n terms of t h e i r values at the present and past time steps. Because of i t s r e l a t i v e s i m p l i c i t y , an e x p l i c i t method was used i n the present c a l c u -l a t i o n s . The p a r t i c u l a r e x p l i c i t timing d i f f e r e n c i n g method used was the Leapfrog method (Roache, 1972). This method can be i l l u s t r a t e d by the f i n i t e d i f f e r e n c e analogue of [B.7:]. If a time g r i d [B.20] t k = kAt i s used, then [B.7] becomes 197 [B..21] / \ k + l , , k-1 ( p c V i + i / 2 . i - frcVi+w,.-! 2At Ar and so, neglecting p' compared to p , c C O [B.22] ( V r } i + l / 2 , j - ( V r ) i - f l / 2 , j 2At [ (P:)L , - ( p ' ) k The Leapfrog method i s thus a centred second order accurate method. It was found that the time s p l i t t i n g phenomenon (Roache, 1972) appeared i f a l l steps were leapfrog. This was cured by taking a forward time step about every 40 steps. The viscous terms were retarded by one time step i n order to avoid an i n s t a b i l i t y of the Richardson type. The s i z e of the time step At was constrained by several e f f e c t s : convection, d i f f u s i o n , and i n e r t i a l e f f e c t s . These constraints are roughly 12 [B.23] At < Ar/V - 10 sec (Courant condition) c [B.24] [B.25] At < (Ar) / v = 1 sec c At < 2TT/fi = 10 2 s ec The most r e s t r i c t i v e of these constraints on the step s i z e arises from the C o r i o l i s term i n the f l u i d equations. The At chosen for these c a l c u l a t i o n s ,„-3 was 10 sec. 198 APPENDIX C: NUMERICAL METHODS FOR MAGNETOSTROPHIC FLOW The magnetostrophic equations were replaced by an approximate set of f i n i t e d i f f e r e n c e equations, which were solved on an IBM 370 computer. A l l v a r i a b l e s were defined on a s i n g l e uniform g r i d (z^ ,x^) i n c y l i n d r i c a l coordinates: [C.l] • x = ( i - l)Ax [C.2] z. = (j - l ) A z , for i , j = 1, 2, .... , N, where [C.3] Ax = Az = R /(N - 1) c For the c a l c u l a t i o n s reported i n t h i s t h e s i s , N was 15. Although the choice of c y l i n d r i c a l coordinates rather than sp h e r i c a l coordinates makes the equations considerably simpler, i t complicates the boundary conditions. The representation of the inner surface of the crust i n the (z^ ,x_^ ) mesh i s one with many sharp edges. A sin g l e c y l i n d r i c a l g r i d (z^,x^) was chosen i n preference to a set of staggered c y l i n d r i c a l grids analogous to the staggered s p h e r i c a l g r i d s used i n the i n e r t i a l flow c a l c u l a t i o n s (see Appendix B) because the computer "bookkeeping" involved i n representing the e s s e n t i a l l y s p h e r i c a l c r u s t boundary i n a number of staggered c y l i n d r i c a l grids i s extremely complicated. In f a c t , a preliminary version of the computer program which used such staggered grids could not be made to work f o r t h i s reason, and was abandoned i n favour of the present program. One consequence of the use df a sin g l e g r i d i s that the s p a t i a l truncation errors are about four times larger than those i n the I n e r t i a l flow c a l c u l a t i o n s . In an e f f o r t to determine what influence these truncation errors had on the computer generated s o l u t i o n , a simulation was performed with N = 20 but with the same time step s i z e as i n the N = 15 c a l c u l a t i o n s . Because the computer cost per time step i n the former simulation i s about twice that of the l a t t e r , the N = 20 simulation was terminated a f t e r only a few days of simulated time; Within that simulated time i n t e r v a l , the N = 20 simulation gave p r a c t i c a l l y the same r e s u l t s as the N =15 one. Thus, neither truncation error nor the i r r e g u l a r representation of the crust boundary i s a serious problem i n the N = 15 c a l c u l a t i o n s . In t h i s g r i d , s p a t i a l d e r i v a t i v e s were replaced by c e n t r a l differences. For example, K B ' ) , ^ - (B'). , [c 4] '(v" x B') — > ""x^'i*1 v " x ' i , j - i ^ r v i + i . j ' v l V i - r , j -2Az Also, the s p a t i a l d i f f e r e n c i n g of the right-hand side of the induction equation [0.5] | f - ^ = V x (v x B) d t C was designed so that V • B was numerically l o c a l l y conserved i n time. I t i s a well-known numerical phenomenon (Roache, 1972) that although the - » - - * • conservation of a quantity such as V • B i s guaranteed by a set of p a r t i a l d i f f e r e n t i a l equations such as [C.5], not a l l d i f f e r e n c i n g techniques lead to the conservation of the quantity i n a numerical simulation. For example, i f c e n t r a l d i f f e r e n c i n g had been applied to the right-hand side of the expanded form 9B' ± ->• ->.->- ->• [C-6] -TT~'= B • V v - (v • V)B - B(V • v ) dt C C C of the induction equation, the di f f e r e n c e version of V • B would d e f i n i t e l y not have been conserved as the simulation proceeded. Spurious magnetic monopoles would have appeared i n the s o l u t i o n . This problem was avoided by c a l c u l a t i n g v c x B at a l l g r i d p o i n t s , and then applying c e n t r a l d i f f e r e n c i n g to the form [C.5] of the i n d u c t i o n equation. This method i s thus numerically f l u x conserving. Time d e r i v a t i v e s were approximated by Leapfrog time d i f f e r e n c i n g (see Appendix B). Thus, i f we define [C.7] g = v c x B , the d i f f e r e n c e form of [C.5] i s z i , j z i , j x/lAx cp i + l , j *-tp i - l , j J [C.9] ( B ' ) k + 1 = ( B ' ) J - 1 - | ^ [ ( £ . ) * , + 1 - , , 1 x 1,3 x i , j 2Az w < p i , j + l 9 i , j - l t c i o ) A Von Neumann s t a b i l i t y a n a l y s i s (Roache, 1972) of t h i s d i f f e r e n c i n g method was c a r r i e d out. I t i n d i c a t e s that the method i s numerically s t a b l e i f At i s chosen to be l e s s than 1//2 times the period of a hydromagnetic-i n e r t i a l wave of wavelength Az. In the present, c a l c u l a t i o n s the numerical f a c t o r was taken to be 1/20, rather than l//2~ , to be conservative. The s t a b i l i t y a n a l y s i s also i n d i c a t e s that the method has zero numerical damping, a d e s i r a b l e feature that i s of t e n associated with Leapfrog time d i f f e r e n c i n g . A t o t a l of 50,000 time steps were taken i n the c a l c u l a t i o n . There was no i n d i c a t i o n whatsoever of numerical i n s t a b i l i t y even a f t e r t h i s very l a r g e number of steps, i n agreement with Von Neuman s t a b i l i t y a n a l y s i s . I conclude that the numerical method i s a sound one. 

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.831.1-0084835/manifest

Comment

Related Items