Science, Faculty of
Physics and Astronomy, Department of
DSpace
UBCV
Choptuik, Matthew William
Doctor of Philosophy - PhD
University of British Columbia
There is continuing interest in the numerical solution of the field equations of general relativity. Many of the most interesting scenarios involve the generation and propagation of gravitational radiation—this thesis is an investigation of numerical techniques which should be useful for such problems. The study is based on a model system of a single massless scalar field coupled to the gravitational field, with the additional requirement of spherical symmetry. The 3+1 formalism for numerical relativity is used and attention is limited to finite difference techniques. As a preliminary step, the properties of some difference schemes are examined in the context of the wave equation in general coordinates. The accuracy, stability, convergence, dissipation, dispersion and efficiency of the schemes are discussed and three of the methods, which are all centred, second order approximations, are identified as being suitable for radiation problems in general relativity. Following work by Piran, the model problem is solved in three different coordinate systems using variations of the schemes devised for the wave equation. Various initial configurations of the scalar field are evolved, ranging from cases where there is a negligible amount of gravitational self-interaction to those where black holes are formed. Although some numerical difficulties remain with the latter calculations, the basic difference techniques developed for the treatment of the radiative degrees of freedom are found to be satisfactory. In the past it has been argued that there is a fundamental inconsistency in finite difference solutions of Einstein's equations which are generated by algorithms that do not continuously impose the constraint equations which arise from the coordinate freedom of general relativity. The deviations of freely evolved quantities from discrete versions of the constraint equations are supposedly of lower order (in the mesh spacing) than the underlying difference scheme. An argument based on Richardson's hypothesis concerning the expansion of the error of a difference solution in powers of the mesh spacing is made which suggests that there should be no such inconsistency. Numerical results from two fully second order algorithms are presented in support of this argument. The remainder of the thesis deals with the application of local mesh refinement techniques to the model problem. A simplified version of a recent algorithm due to Berger and Oliger is implemented for one particular coordinate system and is used to examine the nature of the scalar field's gravitational self-interaction.
A S t u d y of N u m e r i c a l T e c h n i q u e s for R a d i a t i v e P r o b l e m s G e n e r a l R e l a t i v i t y by M a t t h e w W i l l i a m C h o p t u i k B.Sc, Brandon University, 1980 M.Sc.,The University of British Columbia, 1982 A -Thesis S u b m i t t e d in P a r t i a l F u l f i l l m e n t of the R e q u i r e m e n t s for the Degree of D O C T O R O F P H I L O S O P H Y in T h e F a c u l t y of G r a d u a t e Studies D e p a r t m e n t of P h y s i c s We accept this thesis as conforming W the required standard T h e U n i v e r s i t y of B r i t i s h C o l u m b i a July 1986 © Matthew William Choptuik, 1986 In 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 of the requirements f o r an advanced degree a t the U n i v e r s i t y o f B r i t i s h Columbia, I agree t h a t the L i b r a r y s h a l l make 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 and study. I f u r t h e r agree 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 copying of t h i s t h e s i s f o r s c h o l a r l y purposes may be granted by the head o f my department or by h i s o r her r e p r e s e n t a t i v e s . I t i s understood t h a t copying or 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 not be allowed without my w r i t t e n p e r m i s s i o n . Department of Physics The U n i v e r s i t y of B r i t i s h Columbia 1956 Main Mall Vancouver, Canada V6T 1Y3 August 8, 1986 ii Abstract There is continuing interest in the numerical solution of the field equations of general relativity. Many of the most interesting scenarios involve the generation and propagation of gravitational radiation—this thesis is an investigation of numerical techniques which should be useful for such problems. The study is based on a model system of a single massless scalar field coupled to the gravitational field, with the additional requirement of spherical symmetry. The 3+1 formalism for numerical relativity is used and attention is limited to finite difference techniques. As a preliminary step, the properties of some difference schemes are examined in the context of the wave equation in general coordinates. The accuracy, stabil-ity, convergence, dissipation, dispersion and efficiency of the schemes are discussed and three of the methods, which are all centred, second order approximations, are identified as being suitable for radiation problems in general relativity. Following work by Piran, the model problem is solved in three different coordi-nate systems using variations of the schemes devised for the wave equation. Various initial configurations of the scalar field are evolved, ranging from cases where there is a negligible amount of gravitational self-interaction to those where black holes are formed. Although some numerical difficulties remain with the latter calculations, the basic difference techniques developed for the treatment of the radiative degrees of freedom are found to be satisfactory. In the past it has been argued that there is a fundamental inconsistency in finite difference solutions of Einstein's equations which are generated by algorithms that do not continuously impose the constraint equations which arise from the coordinate freedom of general relativity. The deviations of freely evolved quantities from discrete versions of the constraint equations are supposedly of lower order (in the mesh spacing) than the underlying difference scheme. An argument based o n R i c h a r d s o n ' s h y p o t h e s i s c o n c e r n i n g the e x p a n s i o n of the error of a difference s o l u t i o n i n powers of the m e s h s p a c i n g is m a d e w h i c h suggests t h a t there s h o u l d be no s u c h inconsistency. N u m e r i c a l results f r o m two fully second order a l g o r i t h m s a r e p r e s e n t e d i n s u p p o r t of this a r g u m e n t . T h e r e m a i n d e r of the thesis deals w i t h the a p p l i c a t i o n of l o c a l mesh refinement t e c h n i q u e s to the m o d e l p r o b l e m . A simplif ied version of a recent a l g o r i t h m due to B e r g e r a n d O l i g e r is i m p l e m e n t e d for one p a r t i c u l a r c o o r d i n a t e s y s t e m a n d is used t o e x a m i n e t h e n a t u r e of the scalar field's g r a v i t a t i o n a l self- interaction. Table of Contents A b s t r a c t ii L is t of Tables vii L is t of F igures viii Acknowledgement xv C h a p t e r One. In t roduct ion 1 Chap te r T w o . Theoret ica l F ramework 7 2.1 Review of the 3 + 1 Formalism 7 2.2 Description of the Model and the Equations of Motion 14 2.3 Characteristic Analysis of the Scalar Field 23 2.4 Boundary Conditions at r = 0 24 2.5 Slicing Topology, Apparent Horizons and Black Hole Formation 27 2.6 The Structure of the Radiative Equations 29 Chap te r Three . Difference Schemes for the Scalar F ie ld 32 3.1 Introduction 32 3.2 The 1 + 1 Wave Equation in General Coordinates 33 3.3 Description of the Four Difference Schemes 37 3.4 Von Neumann Stability Analysis 45 3.4.1 Analytic Stability Criterion for the C N Scheme 47 3.4.2 Analytic Stability Criteria for the Angled Scheme . . . . 49 3.5 Numerical Determination of Stability Domains for Al l Four Schemes 51 3.6 Convergence Tests for the Four Schemes 56 3.7 Numerical Accuracy Tests of the Schemes 61 3.8 Numerical Radiation Boundary Conditions 68 Chap te r Four . Some Const ra ined Solutions of the Se l f -Grav i ta t ing Sca lar F i e l d 75 4.1 Introduction 75 4.2 Coordinate Conditions 76 4.2.1 The Maximal Slicing Equation 79 4.2.2 The Polar Slicing Equation 80 V 4.2.3 The Isothermal Equation for 81 4.3 The Structure of the Finite Difference Grid 82 4.4 The Choice of a Time Step and Rescaling of the Lapse 83 4.5 Determining Initial Data 85 4.6 The Maximal/Normal Coordinate System 88 4.7 Weak Field Results in the Maximal/Normal System 91 4.8 The Maximal/Isothermal Coordinate System 100 4.9 Weak Field Results in the Maximal/Isothermal System 105 4.10 The Polar/Radial Coordinate System 109 4.11 Weak Field Results in the Polar/Radial System 114 4.12 Intermediate Field Results for All Three Systems 118 4.13 Strong Field Results and Numerical Problems 128 Chapter Five. Geometric Evolution Equations and the Consistency Problem 139 5.1 Introduction 139 5.2 The Argument for Inconsistency in Numerical Relativity 140 5.3 The Argument for Consistency in Numerical Relativity 143 5.4 The Geometric Evolution Equations in the Polar/Radial System 144 5.4.1 Comparison of the Freely Evolved and Constrained Quantities. ... 150 5.5 Demonstration of Consistency in the Polar/Radial System 159 5.6 The Geometric Evolution Equations in the Maximal/Isothermal System 162 5.6.1 Comparison of the Freely Evolved and Constrained Quantities. ... 166 5.7 Demonstration of Consistency in the Maximal/Isothermal System 177 Chapter Six. A n Application of Local Mesh Refinement 183 6.1 Introduction 183 6.2 Previous Approaches to Grid Adaptation in Numerical Relativity 184 6.3 An Unsuccessful Attempt at Local Mesh Refinement 186 6.4 The Mesh Refinement Algorithm of Berger and Oliger 189 6.5 Adaptation of the Algorithm for the Self-gravitating Scalar Field 191 6.6 Difference Equations for the Mesh Refinement Program 191 6.7 The Mesh Refinement Algorithm in Detail 196 6.8 Results from the Mesh Refinement Program 199 vi 6.9 The Nature of the Gravitational Self-Interaction of the Scalar Field. . . . 208 Chapter Seven. Conclusions 215 Bibliography 219 Appendix A . Basic Difference Techniques for Hyperbolic Equations 227 L i s t of Tables The three coordinate systems used in the study of the self-gravitating scalar field 79 Illustration of regridding overhead in the mesh refinement program 200 v i i i L i s t of F i g u r e s 2.1 Spacetime displacement in the 3+1 formalism 10 3.1 Definitions of difference operators 38 3.2 Von Neumann staggered grid structure 40 3.3 Difference star for the standard scheme 40 3.4 Difference stars for the one-sided (upwind) scheme 41 3.5 Difference star for the Crank-Nicholson scheme 42 3.6 Difference stars for the angled scheme 43 3.7 Difference star for the leap-frog scheme 44 3.8 Von Neumann stability domain of the one-sided (upwind) scheme 52 3.9 Von Neumann stability domain of the Crank-Nicholson scheme 52 3.10 Von Neumann stability domain of the angled scheme 53 3.11 Von Neumann stability domain of the leap-frog scheme 53 3.12 Practical stability domain of the one-sided (upwind) scheme 54 3.13 Practical stability domain of the Crank-Nicholson scheme 54 3.14 Practical stability domain of the angled scheme 55 3.15 Practical stability domain of the leap-frog scheme 55 3.16 Convergence test of the one-sided (upwind) scheme 59 3.17 Convergence test of the Crank-Nicholson scheme 59 3.18 Convergence test of the angled scheme 60 3.19 Convergence test of the leap-frog scheme 60 3.20 Typical /? generated by the coordinate transformation (3.70) 62 3.21 Typical c generated by the coordinate transformation (3.70) 62 3.22 Time series of exact solution ri"+^ for accuracy tests of difference schemes 63 3.23 Time series of Hn+* generated by the sided (upwind) scheme 64 ix 3.24 Time series of generated by the Crank-Nicholson scheme 64 3.25 Time series of generated by the angled scheme 65 3.26 Time series of n"+* generated by the leap-frog scheme. 65 3.27 Comparison of Fourier amplitudes of exact solution and three of the difference solutions 67 3.28 Comparison of accumulated phase error in the four difference schemes. .. 67 3.29 Comparison of various numerical outgoing boundary conditions for spherical wave equation 73 4.1 Staggered mesh structure 83 4.2 Weak-field results in maximal/normal coordinate system: <f> 96 4.3 Weak-field results in maximal/normal coordinate system: a' 96 4.4 Weak-field results in maximal/normal coordinate system: 6' 97 4.5 Weak-field results in maximal/normal coordinate system: a 97 4.6 Weak-field results in maximal/normal coordinate system: Krr 98 4.7 Weak-field results in maximal/normal coordinate system: {rb)' - a 98 4.8 Weak-field results in maximal/normal coordinate system: ab2 99 4.9 Illustration of sensitivity of the solution of the Hamiltonian constraint ...99 4.10 Weak-field results in maximal/isothermal coordinate system: a 106 4.11 Weak-field results in maximal/isothermal coordinate system: Krr 106 4.12 Weak-field results in maximal/isothermal coordinate system: a' 107 4.13 Weak-field results in maximal/isothermal coordinate system: /? 107 4.14 Demonstration of weak-field consistency of the maximally sliced solutions 108 4.15 Weak-field results in polar/radial coordinate system: a 116 4.16 Weak-field results in polar/radial coordinate system: Mass aspect (m) . 116 4.17 Weak-field results in polar/radial coordinate system: a 117 4.18 Weak-field results in polar/radial coordinate system: KTT 117 4.19 Intermediate-field results in maximal/normal coordinate system: <f> 120 X 4.20 Intermediate-field results in maximal/normal coordinate system: a' 120 4.21 Intermediate-field results in maximal/normal coordinate system: a 121 4.22 Intermediate-field results in maximal/normal coordinate system: Krr . . . 121 4.23 Intermediate-field results in maximal/normal coordinate system: {Tb)' - a 122 4.24 Intermediate-field results in maximal/normal coordinate system: ab2 122 4.25 Intermediate-field results in maximal/isothermal coordinate system: <f> .. 123 4.26 Intermediate-field results in maximal/isothermal coordinate system: a .. 123 4.27 Intermediate-field results in maximal/isothermal coordinate system: a .. 124 4.28 Intermediate-field results in maximal/isothermal coordinate system: K\ 124 4.29 Intermediate-field results in maximal/isothermal coordinate system: /? .. 125 4.30 Intermediate-field results in polar/radial coordinate system: <f> 125 4.31 Intermediate-field results in polar/radial coordinate system: a 126 4.32 Intermediate-field results in polar/radial coordinate system: Mass aspect (m) ; 126 4.33 Intermediate-field results in polar/radial coordinate system: a 127 4.34 Intermediate-field results in polar/radial coordinate system: Krr 127 4.35 Strong-field results in maximal/normal coordinate system: <j> 130 4.36 Strong-field results in maximal/normal coordinate system: a 130 4.37 Strong-field results in maximal/normal coordinate system: o 131 4.38 Strong-field results in maximal/normal coordinate system: K9$ 131 4.39 Strong-field results in maximal/normal coordinate system: (rb)' - arbK% 132 4.40 Strong-field results in maximal/normal coordinate system: ab2 132 4.41 Strong-field results in maximal/isothermal coordinate system: a 133 4.42 Strong-field results in maximal/isothermal coordinate system: c 133 4.43 Strong-field results in maximal/isothermal coordinate system: K0g 134 4.44 Strong-field results in maximal/isothermal coordinate system: {rb)' - arbK% 134 xi 4.45 Strong-field results in polar/radial coordinate system: <f> 135 4.46 Strong-field results in polar/radial coordinate system: a 135 4.47 Strong-field results in polar /radial coordinate system: a 136 4.48 Strong-field results in polar/radial coordinate system: Mass aspect (m) 136 4.49 Various functions in polar/radial coordinate system 137 5.1 Instability in polar/radial free evolution 149 5.2 Polar/radial coordinate system: Metric function o updated using Hamiltonian constraint 152 5.3 Polar/radial coordinate system: Metric function o updated using evolution equation 152 5.4 Polar/radial coordinate system: Extrinsic curvature component Krr updated using momentum constraint 153 5.5 Polar/radial coordinate system: Extrinsic curvature component Krr updated using evolution equation 153 5.6 Late-time comparison of freely evolved and constrained a' in polar/radial system 154 5.7 Late-time comparison of freely evolved and constrained Krr in polar/radial system 154 5.8 Deviation of freely evolved a from Hamiltonian constraint in polar/radial system: interior region 155 5.9 Deviation of constrained a from evolution equation in polar/radial system: interior region 155 5.10 Deviation of freely evolved Krr from momentum constraint in polar/radial system: interior region ..156 5.11 Deviation of constrained Krr from evolution equation in polar/radial system: interior region 156 5.12 Deviation of freely evolved a from Hamiltonian constraint in polar/radial system: exterior region 157 5.13 Deviation of freely evolved KrT from momentum constraint in polar/radial system: exterior region 157 5.14 Convergence of deviation of freely evolved a from Hamiltonian constraint in polar/radial system 160 5.15 Convergence of deviation of constrained a from evolution equation in polar/radial system 160 5.16 Convergence of deviation of freely evolved Krr from momentum constraint in polar/radial system 161 5.17 Convergence of deviation of constrained Krr from evolution equation in polar/radial system 161 5.18 Convergence of difference between constrained and freely evolved a in polar/radial system 162 5.19 Maximal/isotherm al coordinate system: (Inc)' = A updated using Hamiltonian constraint 167 5.20 Maximal/isothermal coordinate system: (In a)' = A updated using evolution equation 167 5.21 Maximal/isothermal coordinate system: Extrinsic curvature component Ke$ updated using momentum constraint 168 5.22 Maximal /isothermal coordinate system: Extrinsic curvature component K9g updated using evolution equation 168 5.23 Late-time comparison of freely evolved and constrained (In o)' in maximal/isothermal system 169 5.24 Late-time comparison of freely evolved and constrained Ke$ in maximal/isothermal system 169 5.25 Maximal /isothermal coordinate system: A 170 5.26 Maximal/isothermal coordinate system: Keg 170 5.27 Maximal /isotherm al coordinate system: Freely evolved (In a)' = A near r — 0 prior to scaling of geometric variables by r2 172 5.28 Maximal/isotherm al coordinate system: Freely evolved (In a)' = A near r — 0 after scaling of geometric variables by r2 172 5.29 Maximal/isotherm al coordinate system: Difference between constrained and freely evolved functions of metric components near r — 0 173 5.30 Maximal /isothermal coordinate system: Difference between constrained and freely evolved extrinsic curvature components near r — 0 173 5.31 Deviation of freely evolved a from Hamiltonian constraint in maximal/ isothermal system: exterior region 174 5.32 Deviation of constrained a from evolution equation in maximal/isothermal system: exterior region 174 5.33 Deviation of freely evolved K% from momentum constraint in maximal/isothermal system: exterior region 175 5.34 Deviation of constrained K9 $ from evolution equation in maximal/isothermal system: exterior region 175 5.35 Deviation of freely evolved o from Hamiltonian constraint in maximal/ isothermal system: interior region 176 5.36 Deviation of freely evolved K9g from momentum constraint in maximal/isothermal system: interior region 176 5.37 Convergence of deviation of freely evolved a from Hamiltonian constraint in maximal/isotherm al system 178 5.38 Convergence of deviation of constrained a from evolution equation in maximal/isothermal system 178 5.39 Convergence of deviation of freely evolved K9f from momentum constraint in maximal /isothermal system 179 5.40 Convergence of deviation of constrained K9e from evolution equation in maximal/isothermal system 179 5.41 Convergence of difference between constrained and freely evolved a in maximal/isothermal system 180 5.42 Convergence of difference between constrained and freely evolved K9g in maximal/isothermal system 180 6.1 Grid structure for the BKO refinement scheme 187 6.2 Typical grid structure for the mesh refinement program 192 6.3 "Time reversal" test of difference scheme in polar/radial system: a 195 6.4 "Time reversal" test of difference scheme in polar/radial system: n 195 6.5 Pseudo-code form of mesh refinement algorithm 197 6.6 The three possible positionings of the fine grid and the corresponding segment structures 198 6.7 Mesh refinement test for scalar field alone: $ as determined on coarse grid alone 201 6.8 Mesh refinement test for scalar field alone: Coarse grid representation of $ as determined 3:1 refinement 201 6.9 Mesh refinement test for scalar field alone: Fine grid representation of $ 202 xiv 6.10 Mesh refinement test for self-gravitating scalar field: <j> using coarse grid alone 203 6.11 Mesh refinement test for self-gravitating scalar field: <j> using 5:1 refinement 203 6.12 Mesh refinement test for self-gravitating scalar field: Constrained o using coarse grid alone 204 6.13 Mesh refinement test for self-gravitating scalar field: Constrained a using 5:1 refinement 204 6.14 Mesh refinement test for self-gravitating scalar field: Constrained Krr using coarse grid alone. 205 6.15 Mesh refinement test for self-gravitating scalar field: Constrained Krr using 5:1 refinement 205 6.16 Mesh refinement test for self-gravitating scalar field: Freely evolved o using coarse grid alone 206 6.17 Mesh refinement test for self-gravitating scalar field: Freely evolved a using 5:1 refinement 206 6.18 Mesh refinement test for self-gravitating scalar field: Freely evolved Krr using coarse grid alone 207 6.19 Mesh refinement test for self-gravitating scalar field: Freely evolved Krr using 5:1 refinement 207 6.20 Evolutionary sequence of m' for scalar pulse with negligible self-gravitation 210 6.21a Evolutionary sequence of m' for self-gravitating scalar pulse 211 6.21b Continuation of Fig. 6.21a 212 6.22 Late-time configuration of out-going scalar pulse with negligible self-gravitation 213 6.23 Late-time configuration of previously self-gravitating scalar pulse 213 A . l Finite difference mesh (grid) 228 A.2 Difference star for the leap-frog method 229 A.S Domains of dependence—the CFL condition 232 A c k n o w l e d g e m e n t X V It is my pleasure to thank my research supervisor, Prof. W. G. Unruh, for suggesting the model studied here and for providing direction and encouragement. I also benefitted from discussions with U. Ascher, A. Borde, L. Brewin, N. O Mur-chadha, T. Piran, P. Rastall, J. Thornburg and J. Varah. I am indebted to L. Bowman for helping me prepare the manuscript and providing moral support along with my parents, office mates and members of the department. Finally, I would like to acknowledge the financial support of NSERC and UBC. 1 C H A P T E R 1 . I n t r o d u c t i o n It is customary to preface works in numerical relativity with statements re-garding the great complexity, non-linearity and consequent analytic intractibility of Einstein's equations for any but the simplest physical scenarios. There can be little doubt of the truth of these statements, nor of the conclusion that numeri-cal methods provide the only available means to derive detailed information about the nature of the gravitational interaction for many problems of theoretical and astrophysical interest. At the same time, the application of numerical techniques to general relativity should not be considered an effort of "last resort". As with other branches of computational physics, the essential usefulness of numerical relativity is rooted in the ability afforded by large scale computation to simulate the behavior of a complex system. In effect, the computer is used as a "theoretical laboratory" where the physical consequences of the theory may be studied. Such studies may provide answers to definite questions regarding the behavior of a particular physi-cal system, but there is also the possibility of gaining more global insights into the nature of the theory which may suggest new analytic approaches. Many of the most interesting problems in general relativity involve the gener-ation and propagation of gravitational radiation. The data from the binary pulsar PSR 1913+16 [35,90,91,102] has provided strong indirect evidence for the existence of gravitational waves and there is continuing progress in the difficult experimental program for direct detection of gravitational radiation using resonant bars and laser interferometers [19,79]. The astrophysical situations which are expected to be the best sources of gravitational waves, such as supernovae explosions, stellar collapse, or the interaction of black holes are ideal candidates for numerical study. The main goal of this thesis is to investigate numerical techniques which should be useful for the construction of spacetimes (solutions of Einstein's equations) which contain gravitational radiation. As with most previous work in numerical relativity, the work is undertaken within the framework of the 3+1 formalism for relativity, which is reviewed in Chapter 2, and attention is restricted to finite difference tech-niques. The research is based on a study of a single massless scalar field coupled to the gravitational field, with the additional requirement that the spacetimes be spherically symmetric. This reduces the field equations to a system of partial dif-ferential equations which depend on a single spatial variable and time. This system is sufficiently manageable, in terms of human and computer resources, to allow for extensive numerical experimentation. As discussed in Chapter 2, the scalar field is used as the source of energy for the spacetimes because a spherically symmet-ric solution of the field equations cannot contain gravitational radiation. Thus the scalar field supplies the radiative degrees of freedom for the system and the hope is that the techniques which are developed for treating it will also be applicable to gravitational radiation. Numerical relativity is a young field, and although there have been several successful constructions of radiative spacetimes using finite difference methods [10-12,24,25,67,6 8,71,72,8 1,82,87,8 8], the catalogue of proven techniques is quite limited. In addition, relatively little attention has been devoted to some of the more funda-mental issues involved in the numerical solution of any set of time dependent PDE's. In particular, the issues of stability of the difference approximations used and the convergence properties of the generated solutions are not completely understood. Even for relatively simple cases, such as the one studied here, the complexity of the finite differenced versions of Einstein's equations makes these properties difficult to analyze. To circumvent this problem, an even simpler system, which has the basic structure of general relativistic radiative equations, is introduced in Chapter 3. This model is simply the wave equation (one space, one time dimension) in gen-eral coordinates. The principal advantage of studying this system is that exact solutions can be generated for the purposes of evaluating the numerical methods that are used to approximately solve it. Four difference schemes are devised for the system and their convergence, stability and accuracy are examined. Three of the schemes emerge as suitable candidates for the treatment of radiative degrees of freedom in general relativity, at least for problems with dependence on a single spatial dimension. The choice of coordinates is a central issue in numerical relativity. There is no preferred labelling of events in a generic spacetime and a prescription must be given to fix the coordinate system. One of the most interesting features of Einstein's equations from the point of view of numerical analysis is that the structure of the equations is very much a function of the coordinate conditions which are imposed. In addition, a judicious choice of coordinates can allow some of the basic variables which describe the gravitational field to be eliminated from a numerical scheme, usually at the cost of the appearance of additional differential equations for the variables which determine the coordinate labelling. Moreover, because of the coordinate freedom, some of the field equations are not equations of evolution, but rather equations of constraint which must be satisfied by the geometric and matter variables at all times. Numerically, these constraint equations can be used instead of evolution equations to update some of the geometric variables. The net result of all of this, is that a wide variety of algorithm s can be devised for any given general relativistic problem, quite independently of the particular type of differencing which is used [24,26,68-70]. Following Piran's work on cylindrically symmetric spacetimes [68], Chapter 4 describes the solution of the model in three different coordinate systems, all of which have been used in previous work in numerical relativity. For the most part, the coordinate conditions and constraint equations are used to eliminate the evolution equations for the geometric variables, and the scalar field is treated using variants of some of the difference schemes studied in Chapter 3. Unlike most (possibly all) previous treatments of the 3 + 1 equations, the overall schemes for two of the coordinate systems are designed to be uniformly second order accurate in both space and time. Evidence is presented for the mutual consistency of the results produced by the three programs for the case of weak pulses of scalar radiation. Additional results are presented for cases where there is significant gravitational self-interaction of the scalar field. Finally, attempts are made to produce black hole spacetimes, but these are largely unsuccessful due to various numerical problems which remain to be solved. Analytically, if the constraint equations are satisfied at one instant of time, they are automatically satisfied at all other times. This is not the case for numerical solutions of the field equations. It has been observed in many studies that deviations from discrete versions of the constraints invariably develop unless the constraints are continually enforced. In fact it has been argued that this is indicative of a fun-damental inconsistency of finite difference solutions of Einstein's equations, wherein the order to which freely evolved quantities satisfy the constraints is lower than the approximation order of the evolution scheme itself. This argument is reviewed in Chapter 5 and a counterargument based on Richardson's [73] hypothesis regarding the expansion of the error of a difference solution in powers of the mesh spacing is presented. It is argued that there should be no such inconsistency provided that the solution is sufficiently well behaved. The argument is supported by numerical evidence produced from extensions of the two second order schemes described in Chapter 4, which perform simultaneous free and constrained evolutions of the geometric quantities. These schemes also demonstrate that the methods which were developed for evolving the scalar field can also be used for the geometric variables. In addition several examples of numerical instabilities and irregularities which were encountered and the techniques used to eliminate some of them are discussed. Currently, there is much interest among numerical analysts and computa-tional physicists in the development of adaptive mesh techniques for time depen-dent PDE's [4,16,31,38,47,65,94]. These methods involve tailoring the local level of discretization (grid spacing) to fit the features of the solution, so as to produce uniformly accurate results using a minimum of computational resources. There is considerable potential for the application of such methods to radiative problems in numerical relativity. In order to adequately simulate the dynamics of the best sources of gravitational waves, it is necessary to use a grid spacing in some parts of the solution domain which is orders of magnitude smaller than the distance to which the radiation should be tracked [70,93]. In Chapter 6, following a review of previous approaches to grid adaptation in numerical relativity, and a discussion of an unsuccessful attempt at mesh refinement, an implementation of an algorithm based on recent work by Berger and Oliger [4] is described. The method involves the use of a global coarse grid overlayed with a moving fine grid which tracks a pulse of scalar radiation. The program works quite well, and is used to investigate the nature of the gravitational self-interaction of the scalar field more closely. Chapter 7 contains conclusions and some suggestions for additional research. An appendix describing the basic ideas and techniques of finite differencing and stablity analysis which are used in the body is also included. Some knowledge of general relativity and differential geometry will be useful in reading the following chapter [54,99,101]. In particular the tensor notation and units of [54] are used. Briefly, lower case greek indices (/z, v,...) assume the values (0,1, 2, 3) corresponding to the temporal and three spatial dimensions; latin indices (t, j,...) range over (1,2,3). An Einstein summation convention is used for both types of indices—e.g. gijV* = £3i=i 9%jVx- Finally units are chosen such that G (Newton's gravitational constant) and c (speed of light) are both unity, unless otherwise specified. 7 C H A P T E R 2 T h e o r e t i c a l F r a m e w o r k 2.1 Review of the 3+1 Formalism. General relativity may be regarded as a theory of the geometry of spacetime. The spacetime continuum is represented by a 4-dimensional differentiate manifold, M, endowed with a metric ^g^v which provides, at each event of the spacetime, the means to determine distances to nearby events. Additional fields representing various sources of energy and matter may also be introduced on the manifold, and the total non-gravitational energy content of the spacetime is described by the stress-energy tensor . The fundamental notion that "matter tells spacetime how to curve—the curvature of spacetime tells matter how to move" finds mathematical expression in Einstein's field equations: = KTpv (2.1) Here, the Einstein tensor, G^v, which in general is a highly nonlinear and very complicated function of the metric and its first and second derivatives, contains information about the curvature of the spacetime. K is a constant of proportionality, which in the units used in this work is simply 8~. When written in a coordinate system, (2.1) becomes a set of ten coupled, non-linear, hyperbolic and elliptic partial differential equations for the ten components of the metric tensor. Supplemented with the appropriate general relativistic versions of the equations of motion for the various matter fields, this system provides the means to study the effects of gravitation in any physical scenario where these effects may be considered to be purely classical phenonema. For the purpose of seeking numerical solutions to Einstein's equations, it proves useful to recast the system (2.1) in terms of a space-plus-time decomposition of the spacetime manifold. In this way, the dynamical nature of the gravitational field equations becomes manifest, and the task of spacetime construction may be formulated as an initial-value, or Cauchy, problem. The 3+1 split, which forms the basis for most work in numerical relativity, was first performed by Arnowitt, Deser and Misner [1], and hence is often referred to as the ADM formalism. To some extent, the current work uses the reformulation of the ADM decomposition due to York and various collaborators which is described in detail in [105]. In the 3+1 approach, the manifold, My is to be pictured as an infinite stack, or foliation, of 3-dimensionaI space-like hypersurfaces {£} . Each hypersurface, S, is itself a differentiate manifold with an intrinsic geometric structure induced by the embedding spacetime, and the view of a 4-geometry as the time development of a 3-geometry naturally arises. One must identify appropriate quantities which completely characterize a single slice and the way it is embedded. Once these have been specified on some initial surface, the Einstein field equations then determine the corresponding values on future and past slices. At the same time, a prescription for labelling the events of the spacetime must be given since one needs a coordinate system in which to solve the differential equations and display the results. This may be accomplished by the specification of an additional four quantities per spacetime point. The pair of symmetric 3-tensors, {'3^,y, Kij), form an appropriate set of fundamental dynamica/ variables. The first tensor is simply the intrinsic metric of the hypersurface. The second, which is called the extrinsic curvature tensor, may be thought of as a "velocity" of the 3-metric. (A precise definition will be given below.) The four JcinematicaJ quantities are commonly grouped as a scalar function, a, and a 3-vector, Introducing coordinates i M = {x°, x'} = {t, x1} on the manifold, such that t is constant on any hypersurface, the spacetime displacement between events P and P' having coordinates x^ and x^ + dx*1 respectively, may be written as (Fig. 2.1) ( 4 ) ds 2 = dx" dx" = -a2 dt2 + { 3 ) g i j {dx1 + /3,dt)(<fxJ + pdt) = - (a2 + fiiP*)dt2 + 2/3, dt dx' + ( 3 )j,y dx* dxj (2.2) where /?, = ^g%j&' The role of the kinematical variables can now be elucidated. An observer moving in a direction normal to the hypersurface (from event P to event P" in Fig. 2.1) will measure an elapsed proper time of adt, so a is called the lapse function. The same observer, upon arriving at the nearby hypersurface, will find that he must shift his position by P'dt to arrive at an event which has the same spatial coordinates as the event from which he departed—thus, /?' is termed the shift vector. Given an initial slice with some spatial coordinate system, a prescription for specifying the lapse function and shift vector on that and all other slices will determine the coordinate labelling of the spacetime. The arbitrariness in the choice of Q and /?*, as well as the choice of initial slice and its internal coordinates is the 3+1 realization of the coordinate freedom of general relativity. A given foliation of a spacetime defines a time-like vector field, n**, which is composed of the unit length, future-pointing normals to the spacelike hypersurfaces. Associated with n.*4 is the one-form = ^g^n11. Then the components of the extrinsic curvature tensor are defined by Kij — Kji = -Vyn, (2.3) L(t + dt) /S1 dt 1 *f \ i \ i \ I \ t \ 1 4 a dt \ / xx = constant l / i / \ / \ r — J Figure 2.1. Spacetime displacement in the 3+1 formalism. where V M is the spacetime covariant derivative. Equivalently, the extrinsic curvature may be defined in terms of the Lie derivative of the 3-metric along the vector field, which makes its role as a "velocity" of the 3-metric somewhat more apparent. However, this relationship becomes clear enough once the Einstein equations have been rewritten in 3+1 form, as will be done shortly. As might be expected, it is not convenient to work directly with the compo-nents of the stress energy tensor, T^, in the 3+1 approach. Instead, it is more useful to introduce a scalar p, a 3-vector j*, and a 3-tensor 5,y, defined as follows: p = nllnvT>"/ (2.4) (2.5) 11 These quantities are called the energy density, (3-)momentum density and spatial stress tensor, respectively. Since there is no unique notion of energy, or momentum or stress in general relativity, it should be noted that these terms refer to the fact that observers "momentarily at rest" in the slicing would, for example, determine the local energy density to be p [105]. The derivation of the 3+1 Einstein equations from the generally covariant form (2.1) will not be performed here since the calculation is well documented [78], [105]. The specific form of the equations displayed below is essentially the same as that used by Piran in his study of cylindrically symmetric general relativistic systems [68]. It is necessary, however, to introduce some additional notation and definitions. The 3-metric ^gij will be denoted simply as gij. The associated reciprocal tensor, gl}' satisfies g%3gjk — '^k where S1^ is the usual Kroenecker delta. These two objects are used to raise and lower indices of all 3-tensors. Partial differentiation with respect to time will be abbreviated by an overdot, / ) t is synonymous with df /dx* and a '|' indicates covariant differentiation within a 3-manifold. Four of the Einstein equations do not contain second time derivatives of the metric and therefore represent equations of constraint which must be satisfied by the geometric and matter variables on every hypersurface. They are the Hamiltonian constraint R - K'jK'i + K2 = 16icp (2.7) and the momentum constraints (2.8) where the (3-dimensional) Ricci scalar, R, is the trace of the Ricci tensor R = g'iRij (2.9) and K is the trace of the extrinsic curvature tensor K = g^Ki:j (2.10) An explicit expression for the components Rij which will be needed shortly is p pre pre . pTL p m p n p m 1 1 \ •"•i) — 1 ij,n 1 in,j i 1 n m 1 xj *• jm*- in K^-11-) where the Christoffel symbols Tljk are defined by I","* = \ 9in{9nj,k + 9nk,j ~ 9jk,n) (2.12) In addition to the constraints, there are six equations of evolution for the metric components which are directly derivable from (2.3) 9ij = -2aglhKkj + pkgxj>k + gikfikj + ?fcy/?fc>t- (2.13) and six evolution equations for the extrinsic curvature tensor which come from those of the Einstein equations which do contain second time derivatives. K\ = pKUt + K>kpkti - KkiP\k - a|,-'> + a{R\ + KKJi + 4TT(S - p)6^ - SnSj{) (2.14) The two basic stages in the 3+1 construction of a spacetime may now be summarized. First, initial data must be constructed. This involves deciding what topology the initial hypersurface is to have, the introduction of a spatial coordinate system and a specification of an initial data set {<7,y, if*-; p, ji}, which satisfies (2.7) and (2.8). Because the constraints must be solved, this part of the procedure is, in general, far from trivial. Only 8 of the 12 geometric quantities can be freely chosen, and it is not at all obvious which 8 these should be. An elegant and general treatment of this problem, largely due to York and O Murchada [63,105-107], has provided considerable insight into the theoretical na-ture of the initial data problem, as well as the basis for most of the recent numerical work involving the construction of initial data sets [5,13,45,92,107]. This approach will not be discussed here since the model which is studied in this work is sufficiently simple that reasonable initial data may be generated by an ad hoc scheme. The second stage of a 3+1 generation of a spacetime is the evolution of the initial data. As should be clear from the previous discussion, as well as the struc-ture of Eqs. (2.13) and (2.14), this requires the imposition of coordinate conditions which will determine the values of the lapse, a, and shift j9' everywhere. This specification need not be explicit—the coordinate conditions may be such that a, or /?', or both must satisfy differential equations on each slice. In any case, once the coordinate choices have been made, a and /8' are to be regarded as "given" functions. The right hand sides of (2.13) and (2.14) can then be evaluated on the initial slice, providing the means to determine the geometric quantities on a nearby advanced or retarded slice. Assuming that additional evolution equations for the matter fields are available, the values of the matter variables may also be evolved in this fashion, allowing for a slice by slice construction of the spacetime. In general, such a procedure can be continued indefinitely unless a physical or coordinate sin-gularity develops in the solution. Unfortunately, both of these situations are very real possibilities in general relativity, particularly in strong-field cases. By virtue of the Bianchi identities V G „ „ = 0 (2.15) which arise from the coordinate invariance of general relativity, any initial data set which satisfies the constraint equations, will, when evolved to some future or past hypersurface, automatically satisfy the constraint equations on that slice. Equiv-alently, one might picture the generation of a spacetime in terms of the "free evo-lution" of some minimum number of quantities (depending on the coordinate con-ditions) with the remaining geometric variables automatically adjusting so as to keep the constraints satisfied. The point to note here is that when the procedure of spacetime construction is actually implemented numerically, definite decisions must be made about how each quantity is updated from one time level to the next. For example, one may choose to use one or more of the constraint equations to compute advanced values for some of the geometric variables instead of updating them using the evolution equations. It has been the experience of numerical relativists that the different types of schemes which can be devised sometimes display significant differences in their numerical behavior. This has led to a considerable amount of concern and speculation regarding the precise relationship between the constraint and evolution equations in numerical solutions of the field equations. This prob-lem will be discussed in more detail in Chapter 5 where a probable resolution is suggested. 2.2 Description of the Model and the Equations of Motion. As stated in the introductory chapter, the major goal of this thesis was to investigate and develop techniques which would be useful for radiative problems in numerical relativity. In such a study, it is imperative that the model being exam-ined be sufficiently manageable to allow for experimentation with different solution methods, coordinate systems and levels of discretization. Given the complexity of the field equations and the current availability and cost of computer resources, such demands essentially preclude systems with dependence on more than one spatial di-mension. Today's state-of-the-art numerical relativity codes generally involve two spatial dimensions, but it is clear from the amount of manpower it takes to con-struct such a program, the increasing reliance on the use of supercomputers, and the tendency to use fixed coordinate systems and numerical techniques, that such codes are not the ideal vehicle for performing numerical experiments. The number of spatial dimensions that must be dealt with in a physical system may be reduced by imposing symmetry conditions. Given that the spacetimes to be studied are to depend on time and one space dimension only, there are essentially three choices—the manifolds can have spherical, planar or cylindrical symmetry. Intuitively, a spherically symmetric spacetime, at any instant of time, is invariant under all rigid rotations about some preferred event. Similarly, planar symmetry implies invariance under translations in two independent spatial directions, and cylindrically symmetric systems are invariant under rotations about and transla-tions along some preferred spatial axis. All three types of spacetimes have been studied numerically—a complete catalogue up to the early 1980's is given in [69]. Cylindrically and planar symmetric systems have the advantage that they may con-tain gravitational radiation; spherically symmetric spacetimes may not by virtue of Birkhoff's theorem [54] which states that the only vacuum, spherically symmetric solution of Einstein's equations is the static Schwarzschild solution. (Maxwell's equations for electromagnetism do not admit spherically symmetric radiative so-lutions either. Both of the these results may be understood as a consequence of the transverse nature of electromagnetic and gravitational radiation). At the same time, of the three, spherically symmetric systems are surely the best "first approxi-mation" to many of the most interesting astrophysical scenarios, particularly those involving a single, isolated object. The numerical studies of spherically symmetric spacetimes which have been performed to date [49,52,59,66,76,77,104] have all employed a perfect fluid [54] as the stress-energy source. This involves the solution of the equations of relativistic hydrodynamics along with Einstein's equations. These investigations have long been used to examine various features of stellar collapse and black hole formation, as well as for the development and testing of numerical techniques. The quality of the results has continually improved—good examples of the accuracy of some of the most recent codes can be found in [49] and [66]. There are no radiative degrees of freedom in these hydrodynamical systems, however, and it is not clear (and in fact probably not true), that the numerical methods developed for dealing with a fluid are the best things to use for evolving radiation, of any type, in a general relativistic system. The simplest possible way to introduce radiative behavior into a spherically symmetric spacetime is to use a single, massless scalar field, <j>, as the source of stress-energy [95]. In flat spacetime, such a field satisfies the massless Klein-Gordon equation 4>** = 0 (2.16) which is derivable from the Lagrangian density £# = YLf = -K'YV^ (2.17) where is the Lagrangian scalar, T = y/— with ^g the determinant of the 4-metric and «' is some positive constant. The easiest way to make the scalar field interact gravitationally is simply to add the above density to the gravitational Lagrangian density = r ( < 4 > J Z ( 2 . 1 8 ) where ^R is the 4-dimensional Ricci curvature scalar. Demanding that the action 1= j d*x£, (2.19) be stationary with respect to variations in <f> gives the equation of motion of the scalar field ^ = f ( r ( V ^ ) , , = o (2.20) Similarly, varying the action with respect to the 4-metric components leads to Ein-stein's equations. In particular, the stress energy tensor may be calculated from the following expression [54] r ^ = - 2 J ^ + ( 4 ) ^ (2-2l) which gives r«„ = 2ic'(a\Mf„ - l- ( 4 W ' % ) (2.22) The restriction to spherical symmetry results in an enormous simplification of the Einstein equations. Choosing polar-spherical coordinates (r, 9, <p) adapted to the symmetry, it is well known [54] that the most general 3-metric may be expressed as gij = diag {a2(r,t), r2b2{r,t), r2b2 sin2 9) (2.23) This in turn implies that the extrinsic curvature tensor is also diagonal, with only two independent components Klj = diag(/T r(r,0, K%(r,t), K\) (2.24) Finally, the shift vector has only one non-zero component 0*" = (0r(r, 0,0,0) = (0,0,0) (2.25) 18 as the other two components are annulled by the choice of coordinates. Then the spacetime line element, in its most general form, is [4)ds2 = ( -a 2 + a2/?2) dt2 + 2a2/? dtdr + a 2 dr2 + r2b2 dU2 (2.26) where dU2 = d92 + sin2 0d<f>2 (2.27) Thus, there are only 6 of a possible 16 geometric variables to deal with and in addition, the expressions for quantities such as R are much simpler than in the general case. The connection coefficients Y%jk may be calculated directly from (2.12), al-though the variational approach described in Sec. 14.4 of [54] involves less work. The non-vanishing coefficients are: 1 rr = — 1 60 = ~ o 1 d,i> — Sin 01 $0 a 2al T " = 2 M r V = -sin0cos0 Y+T4 = Y9re r ^ = cot0 (2.28) where ' = d/dr. The calculation of the components of the Ricci tensor jR'y may now be performed using (2.11). For example: RT r — grrRrr = g r r ( - ( r ' r , + r * r # ) ' + ( r ' r , + r*r*)Tr„ - r ' P , 8 - T V ) (2.29) Using the above expressions for the connection coefficients, after some manipulation this becomes aro \ a ) Similarly = = ^ - ( 7 <'*>')') ( 2 S 1 ) The term a.^ which also appears in the evolution equation for the extrinsic curva-ture components is evaluated using a|, b =0yfc(a,-*-rB,-fca>) (2.32) with the result 'R -: ( 7 ) ' ( 2- 3 3» Turning now to the calculation of the stress-energy quantities, is is convenient to introduce the following auxiliary functions: $ EE 4>' (2.35) 11 = -U - M) (2.36) a Also from (2.26), the non-vanishing components of the four metric are given by Wgtt = -o? + a2{l2 W9tr = {%rt = a^ W9rr = a2 <4 W = r2b2 = r2b2 sin2 9 (2.37) and the corresponding contravariant components are found from ^o^ 4 ' ]g p i / — 6^" 1 ' g — — or 'g — 1 'g — pa Wg" = a~2 - /?2«-2 ( V " - r - 2 r 2 ( V ' = r - 2 6 - 2 s i n - 2 0 (2-38) Then 9 = 'g 4>,t + K 'g <P,r n (2.39) and similarly ($2 - n 2) = — 1 (2.40) a' Noting that the components of the one-form are simply (—a, 0,0,0) and taking k' = 1/2, the energy density is computed from (2.4), (2.22) and the last two results 2 rjitt \ Til J P = a T = 2a' ( 2 A l ) Similarly, the only non-zero component of the momentum density, jr, is, from (2.5) jr = <xT\ = (2.42) a and the stress components S'j = </ , f c5 f c ; (2.43) where S^j is defined by (2.6) are ("<) , * (n ! - * 2) Using the above results, the Hamiltonian constraint (2.7) becomes (2.45) $ 2 + n 2 . . = 8TT (2.46) a To compute the single non-trivial momentum constraint, first note that A r\i — r,» + 1 fii-R r ~ 1 r»A n = Krr'+ 2T9r9{KrT - K\) (2.47) K\r = (KTT + 2K99)' (2.48) so that using (2.42), (2.8) yields K99'+[^f{K99-K\) = A« — (2.49) The evolution equations for the metric functions a and b are easily calculated from (2.13) with the result d = -aaKrT + (a/?)' (2.50) b = -abK9, + £{rb)' (2.51) 22 The evolution equation for the extrinsic curvature component Krr is calculated from (2.14) + a {Rrr + KK\ + 4w{S - p) - 87rSrr) (2.52) which, using previously computed results, becomes k% = pK Similarly a \ a / + <* I~ A (—) + KRTr ~ **K ) (2-53) V arb \ a } a1 j k\ = PH',' + - - i - j (r6)')' + «tf JT P (2.54) [rb) a (rb) \ a / Expressions (2.46-54) contain the complete set of 3+1 Einstein equations for the current model. The only remaining equations to be dealt with are the evolution equations for the scalar field. First, solving for <f> from the definition of IT given by (2.36), then differentiating with respect to r gives an evolution equation for $ 6 = + - l l ) ' (2.55) Next, multiplying (2.20) by T and using expressions (2.38) gives ( r ^ " ^ ) , 1 , = ( r ^ 0 , M ) t l + ( r , ' H ' ^ ) i r = (aar 2 6 2 (-a- 2^ + / ? a - y ) ) t + (aar262 ( (a - 2 - /?2oT2) <f>' + /3a - 2^)) = 0 (2.56) Using the definition of II, this becomes (r 2 6 2 n) i t = (r262 ( / ? n + ^ ) ) (2.57) Solving for II, and using (2.51) gives the final result * = ^ (''"' ( ' N + ; * ) ) ' + 2 { A K - " T T ) N <*•"> It should be noted that it is necessary to rewrite the scalar field equations as a pair of coupled first-order equations to avoid the appearance of terms such as d or /?, which are difficult to treat in the 3+1 formalism. 2.3 Characteristic Analysis of the Scalar Field. Equations (2.55) and (2.58) form a first order, quasi-linear system (see, for example, Courant and Hilbert, Vol II, Chap. 5 [18]) for the radiation field. Some insight into the nature of this system may be gained by casting it into the form u t + A u , = B (2.59) where u is the two-component vector ($,II) and A and B are coefficient matrices which depend on u, the geometric variables, r, and t. Specifically A = - U %') <"•» and the characteristic directions t — dr : dt of the system are given by |A - rl | = 0 (2.61) and a simple calculation gives a T = - f i ± -a (2.62) Thus as long as a/a ^ 0, the system has two real, distinct characteristic directions and is therefore hyperbolic. Note that if a is 0 somewhere in the interior of some slice, then the slicing is "bad" from a 3+1 point of view. Such a case indicates either that the foliation is not strictly spacelike, or that the slices are degenerate, which is a type of coordinate pathology. It will always be assumed in the following that a/a > 0—any violation of this inequality will necessarily stop the actual evolution of a spacetime. The quantities given by (2.62) determine the local directions of propagation of (small) disturbances in the scalar field, from which the local domains of influence and dependence may be constructed. This information is of use in understanding the stability conditions of some of the difference schemes which will be devised to solve the model. The above characteristic analysis may seem incomplete in the sense that the evolution equations for the geometric variables have not been treated. Actually, this is not the case, since as noted previously, a spherically symmetric spacetime can possess dynamics only where there is matter present. Thus, "propagating dis-turbances" in the geometric variables will have to be tied to disturbances in the scalar field. From a more fundamental point of view, the 4-metric contains, at each event, all the information necessary to determine the local causal structure of the spacetime. In particular, the past and future pointing null-cones [36] delimit the domains of dependence and influence of the event. Propagation of any information to or from the event must occur within, or on the light cones. Thus another way to determine the characteristic directions given by (2.62) is to look for null directions dr : dt defined by Wds2 = -a2dt2 + a2{dr + 0dt)2 = 0 (2.63) 2.4 Boundary Conditions at r = 0. The Cauchy problem for the self-gravitating scalar field is to be solved on the half plane t > 0, r > 0. Because the scalar field has outgoing characteristics as well as incoming ones, a boundary condition for <f> at r = 0 for t > 0 has to be supplied to make the problem well-posed. Such a condition may be derived simply by demanding that </> be differentiate at r — 0 [96]. In a polar-spherical coordinate basis (f,0,<p), and with spherically symmetric <j>: V<£ = <p'r (2.64) which is not defined at r = 0 unless <j>'{0,t) = 0 (2.65) Similar inner boundary conditions may be derived for the geometric quantities {gij,Klj} using a procedure described in detail by Bardeen and Piran in [2]. Basi-cally, the idea is to demand that all tensor components be regular at r = 0 in the sense that they can be expanded in non-negative powers of "Cartesian" coordinates i , y, and z defined by x — r sin 6 cos <p y — r sin 0 sin ip z — r cos 6 (2.66) Bardeen and Piran note that this imposes a restriction on the coordinate system as well as the geometry of the hypersurface, and that a breakdown in regularity defined in this manner need not imply any physically irregular behavior. The attitude of the above authors, which is also adopted here, is to assume the above regularity initially, and worry about irregular behavior only if it actually develops during the evolution. There is no practical obstacle to constructing regular initial data; the hope is that such initial data sets will generate all or "almost all" of the physically relevant spacetimes. This approach seemed to be adequate for the current work— there was never any evidence of dynamically induced breakdown of regularity, as defined above, at r = 0. The analysis in [2], performed for axisymmetric spacetimes, is considerably more general that what is required for spherically symmetric systems. The only re-sult needed here is that all of the quantities a, b, Krr and K9$ must have expansions in even powers of r near r — 0. This immediately implies the following boundary conditions: a'[0,t) = b'{0,t) = Krr'{0,t) = Kee'{0,t) = 0 (2.67) In general, these relations will be consistent with the evolution equations (2.13), (2.14), only if the kinematical variables satisfy o'(0,t) = 0{O,t) = 0 " (2.68) Some caution must be exercised when imposing coordinate conditions which deter-mine a or /? to ensure that the conditions are compatible with these relations. An additional condition on the metric components at r = 0 arises from the demand that the spacetime be locally flat there. The condition may be obtained by considering the parallel transport of an arbitrary vector around a loop enclosing r — 0 in the equatorial plane of any hypersurface, and demanding that there be no net rotation in the limit that the loop is shrunk to a point. The result is well known [48] and will not be derived here; in the current case, it is simply =0 (2.69) r=0 which can be easily manipulated using expressions from the previous section and the boundary conditions (2.67) and (2.68) to yield o(0,t) = 6(0,*) (2.70) From (2.50-51), it is not hard to see that this and the other boundary conditions then imply Krr{0,t) = K99{0,t) (2.71) In addition to the inner boundary conditions displayed above, some attention will have to be devoted to the problem of devising conditions at an outer boundary at some finite radius, since numerical solution domains cannot be unbounded (unless some compactifying coordinate transformation is employed). Discussion of the these conditions will be deferred to later chapters. 2.5 Slicing Topology, Apparent Horizons and Black Hole Formation. The topology of the spacelike hypersurfaces {£} must be included as part of the specification of the Cauchy problem for the generation of any spacetime. All of the solutions which are constructed in this thesis have the spatial topology R 3 and are asymptotically Hat. Thus the 3-metric on any slice tends to the usual Euclidean metric on R 3 at large distances. Among other things, this allows for the definition of the conserved total mass, (or ADM mass, or ADM energy), M of any of the spacetimes studied. M can be expressed in several forms. For example, if the spherically symmetric line element is written as ds2 = 1>'{dr2 + r2dn2) (2.72) by a suitable choice of radial coordinates, then the mass is given by [6,105] M = - — / V 1' d2Si = lim -2r2V> r (2.73) 2n J r - > o o ' where the integral is taken over the 2-sphere at infinity. Physically this mass could be deduced by observers at large distances from the orbital motion of test particles. The configurations of the scalar field which are evolved in the following are typically comprised of two wave packets—one ingoing and one outgoing—which effectively propagate independently of one another. Because of the spherical sym-metry, it is possible to talk meaningfully about the mass associated with each pulse. For the cases studied, the outgoing wave invariably disperses to infinity, but this need not be the case for the ingoing wave. Using a radial coordinate, r, such that the proper surface area of a r = constant 2-sphere is 4nr2, two length scales can be associated with the ingoing pulse. One is the spatial extent, /, of the pulse—the other is the Schwarzschild radius rs = 2M t r e j 0 , n ? associated with its energy con-tent. For a "weak" pulse, I 3> r s and the pulse will self-reflect off r = 0 and then propagate outwards. However, if / w rg, or smaller, then there is the possibility that the gravitational field will become strong enough to form a black hole—a re-gion in the spacetime from which light rays cannot escape to infinity—as was first demonstrated by Christodolou in 1971 [15]. The boundary of a black hole, known as the event horizon, is a global feature of the spacetime; its location cannot be determined by local measurements. Instead, the trajectories of light rays originating from each event in the spacetime must be traced out to determine which events can communicate with infinity and which can-not. This procedure can be implemented relatively easily in spherically symmetric numerical work [66,76,77], but has not been used here. Instead, the existence of a black hole in a given solution is inferred from the appearance of apparent horizons, also known as marginally trapped surfaces [36]. An apparent horizon is defined as a 2-surface (topologically S2) on which the divergence of the outgoing light rays vanishes, and it can be shown [36] that apparent horizons always lie on or within the event horizon. As Hawking and Ellis [36] point out, the advantage of working with apparent horizons rather that event horizons is that the former depend only on the geometry of the hypersurface in which they are embedded. In component form, the equation for a trapped surface is [25] s'li - K + faJKij - 0 (2.74) where S*1 is the outward pointing, space-like, unit normal to the surface. Special-izing to the case of spherical symmetry, with the line element (2.26) s* = (a_1,0,0) (2.75) so Thus, using K = K\ + 2K% (2.77) (2.74) becomes [rb)' - arbK% = 0 (2.78) The left hand side of this last expression is easily calculated on each hypersurface and its vanishing signals the formation of an apparent horizon, which in turn implies that the spacetime contains a black hole. 2.6 T h e S t r u c t u r e of the R a d i a t i v e E q u a t i o n s . The structure of the 3+1 field equations has been discussed in some detail in the past in an attempt to relate the equations to simpler model systems, and to devise numerical methods [80]. In this section, a few comments in this vein are = 2 (r&T arb (2.76) made regarding the current problem which provide some motivation for the study undertaken in the next chapter. The most crucial issue which must be faced in the finite difference solution of time-dependent problems, and particulary hyperbolic problems, is that of stability. It is a relatively easy matter to construct finite difference operators which approxi-mate the differential operator to an acceptable accuracy, but the discrete equations which result often admit solutions which "blow up" as the calculation proceeds, rendering the scheme computationally useless. It is a general result of the stability theory for finite difference schemes that it is the treatment of the highest order spatial derivatives which governs the stability of the method [74]. For the case of the scalar field equation (2.55) and (2.58), this means that all terms involving or IT' will be important for stability. Rearrangement of the equations to make these terms appear explicitly gives $ = 0 $ ' + - n ' + ••• (2.79) a ll = /m'+ - $ ' + ••• (2.80) a The terms involving /3 in these expressions will be called convective terms in analogy to the simple convection equation u — cu' — 0. They also appear in the evolution equations for the geometric variables and have been responsible for many of the stability problems that have been encountered by numerical relativists, including the author. Part of the problem stems from the fact that the schemes are usually constructed so that it is easy to treat the non-convective terms in a stable and accurate fashion; it is then more difficult to difference the convective pieces to the same order of accuracy while retaining stability. It is always possible to eliminate the convective terms by using the coordinate condition 0 = 0. Historically, this has been a fairly common practice and two of the three coordinate systems used in this thesis employ this choice. However, particularly for multi-dimensional problems, there are strong arguments to be made for reserving the spatial coordinate freedom for purposes other than eliminating these terms. Accordingly, the main purpose of the next chapter is to devise and study some diiference schemes which are suited to the solution of equations of the form (2.79-80). It is expected that variations of these methods will be applicable to more general radiative problems in numerical relativity. 32 C H A P T E R 3 Difference Schemes for the Scalar F i e l d 3.1 Introduction. Given a set of PDE's to be solved numerically, one must first decide how to convert the continuum differential system to a finite algebraic system so that the solution of the latter is an acceptable approximation of the former. Now, a theoretical physicist will probably consider the solution scheme acceptable if he can be reasonably sure that the results it produces represent the physical content of the model as well as the model represents the actual situation he is considering. A numerical analyst, on the other hand, would like to have more detailed information about the behavior of the scheme, including its stability and convergence proper-ties. Finally, the person who actually implements the scheme will not consider it acceptable if it is extremely difficult to program or consumes an excessive amount of computer resources. One might even consider classifying schemes according to these criteria as physically, numerically or computationally acceptable. Obviously, there is considerable overlap among these rather vaguely defined categories. For example, any scheme (and such schemes exist) for which stability and convergence are guaranteed is, in principle, as useful to the theoretical physicist as a method for generating exact solutions. To date, studies in numerical relativity have been almost exclusively geared towards the development of physically and computationally acceptable schemes. This is understandable since most projects have been designed to answer questions of genuine astrophysical interest. The discrete systems of equations involved in such cases are extremely difficult to analyze, and in most instances, researchers are far more concerned with the physical features the solutions display than the detailed behavior of the numerical method. At the same time, a reasonable amount of attention is devoted to the issue of computational efficiency, since these problems are invariably limited by the amount of computer resources which are available. However, there are arguments to be made for devoting more attention to fundamental numerical issues. First, as mentioned above, a numerically accept-able scheme is almost certainly physically acceptable. In addition, a proper un-derstanding of the limitations of the method being used reduces the chances that physical significance is ascribed to solution features which are primarily artifacts of the numerical approximation. This is particularly true for cases when the level of discretization is relatively coarse—as it is for most state-of-the-art relativity calcula-tions. Finally, theoretical knowledge of how a numerical scheme should behave can be extremely useful in the actual construction and verification of computer codes—a feature which even numerical analysts sometimes tend to underemphasize. Keeping this point of view in mind, the remainder of this chapter is devoted to a study of a test problem which captures many of the relevant features of the model described in the previous chapter, but which is much easier to analyze. 3.2 The 1 + 1 Wave Equation in general coordinates. The problem to be examined is the massless Klein-Gordon equation in one space and one time dimension, more commonly known as the wave equation. As-suming that the standard spatial and temporal coordinates are being used, this is an extremely simple problem to solve—both analytically and numerically. This fact will be used to advantage in the following for the purposes of evaluating various difference schemes. However, the methods will be designed to be applicable to the wave equation in general coordinates, which is essentially what is required for the treatment of the self-gravitating scalar field. Consider a flat 2-dimensional manifold with Cartesian coordinates xl = [T, X) (In this chapter, component indices range over the values 0,1 corresponding to time and space, respectively and units are chosen so that c = "speed of light" = 1). Then the metric is simply ~gij = gi]: = diag(-l,l) (3.1) and the scalar wave equation is just -$TT + j>xx = 0 (3.2) Assuming that the solution domain is — oo < X < oo, T > 0, the specification of the problem may be completed by providing initial data in terms of two arbitrary functions f(X) and h{X) representing the ingoing and outgoing parts of the solution at the initial time: 4>(X,0) = f(X) + ~h(X) (3.3) Then, at any later time the solution of (3.2) is simply 4>{X,T) = f{X + T) + h{X-T) (3.4) The wave equation can be rewritten in a symmetric first order form by introducing the auxiliary variables ^ = 4>x n = 4>T (3.5) Then (3.2) becomes $ r = rix r i j — $ x (3.6) 35 with initial data *{X,0) = f'{X) + h'{X) fl{X,0) = f'{X)-h'(X) (3.7) and the complete solution is $(X ,T ) = f'{X + T) + h'{X - T) fl{X,T) — f'(X + T) - h'(X - T) (3.8) where ' denotes ordinary differentiation. Consider now the action of a coordinate transformation on this system. For the studies undertaken below, it is sufficient (and from a practical point of view, relatively necessary) to restrict attention to transformations of the form £' - xi = (x,t) = {X,t{X,T)) (3.9) In these coordinates, the metric components g'3 are U _ t 2 -TT , , 2 -XX * 2 , * 2 g =tT g +tx g - -ty + tx 9tx = gxt = t x g X X = tx gxx = l (3.10) Now, in the 1+1 analogue of the 3+1 formalism, the metric is written as v » - L " - ~ ' . - - % " » . - ) ( 3 n ) corresponding to the line element ds2 = -a2 dt2 + a2{dx + 0 dt)2 (3.12) 36 Identification of components in (3.10) and (3.11) leads to the following expressions for the lapse, shift and "spatial metric component" <* = [tT2-tx2r* (3-13) P = tx{tT2 - ix2)'1 (3.14) •-Ms8)' (si5) Now, a simpler version of the calculation performed in Sec. 2.2 gives the following first order form of the wave equation $ = (/3$ + cn)' (3.16) f[ = (/m + c$)' (3.17) w here 6- 36' * = 6' IT = — — (3.18) C S 7 = t *tXt 1 <3-19) and ' and ' now, and subsequently, denote differentiation with respect to t and x. The relationship between 11} and 11} is determined using elementary rules of coordinate transformations $ = $ - — f l (3.20) n = n - — $ (3.21) Thus, given the inital data (3.7) and the explicit form of the coordinate transfor-mation (3.9), the exact solution in (x,t) coordinates can be determined from these last two expressions and (3.8). This will require the inversion of the coordinate transformation to determine T as a function of t and i. Equivalently, one can de-termine only the initial values $(x,0), II(x,0) and subsequently evolve the system using (3.16-17). This again requires that T(x, t) be known so that tj and tx may be evaluated. In the application of this procedure to be described below, the inversion was performed numerically, as was the evolution of the system (3.16-17). At the same time, the exact solution was generated for comparison purposes from (3.8). It should be noted that the generalization of this technique to the case where x also depends on X and T would be considerably more difficult to implement. Since the transformation of the temporal coordinate alone is sufficient to generate non-trivial lapse and shift functions, which is the whole purpose of this exercise, the additional effort does not seem warranted. At this point, readers who are unfamiliar with the subject of finite difference methods for time dependent problems should read Appendix A for an overview of the basic ideas and techniques which will be used below. In addition, an operator notation for finite difference operations is used extensively in this and following chapters so that the difference equations can be displayed explicitly in a relatively compact form. Fig. 3.1 gives the definitions of all of the operators used in the thesis. 3.3 Description of the Four Difference Schemes. Consider now the problem of constructing difference schemes for the system There are many possible ways of differencing this set of equations, so to narrow the possibilities, some basic "rules-of-thumb" which originate from previous work, as well as from the desire to apply the schemes to general relativistic problems, are (3.22) n = (/3n + e * ) (3.23) Forward and Backward Difference Operators. / . n i l _ »\ A L / ; EE ± V ' A T i ] - / + 0 ( At>) Central Difference Operators. A ^ A i / ; = A*_A*+/;" = f i + l ~ % 2 + f'~l = /"; + 0( As*) Averaging Operators. /^; = ^(/;+/;±,) = /;±'+o(At') Half-step Extrapolation Operator. x«/" = */»_ = /n+Uo( A«») 2 * 2 * > Note: The averaging and extrapolation operators have precedence over other algebraic operations. For example: Figure 3.1. Definitions of difference operators. Note that k and h are sometimes used as synonyms for At and Ax and that, in later chapters, the spatial variable is r, rather than x. employed. The first basic guideline is to keep the schemes centred in time as well as space as much as possible. This generally results in the highest possible order of approximation for a given diameter of the basic difference star. In addition, such schemes very often have a low degree of dissipation [43,75], which is important for hyperbolic systems such as (3.22-23) (as well as for gravitational radiation), where the radiative degrees of freedom do not lose energy via any damping mechanism. Secondly, the diameter of the finite difference stars will be restricted to 3 or fewer mesh points. This constraint greatly simplifies the problem of dealing with spatial boundaries, and for all of the cases studied in this thesis, means that the number of numerical boundary conditions which must be specified is the same as the number of boundary conditions available from physical or analytical considerations. In addition, by demanding that the only schemes with a temporal diameter of 3 be explicit, a maximum of two time levels of any grid function must be stored. This is not crucial for the one dimensional problems studied in this work where storage requirements are not a limiting factor, but it is important for multi-dimensional problems. These considerations essentially limit the schemes to 1st or 2nd order explicit or implicit methods, unless some sort of compact differencing [40] is used. Finally, smoothing techniques or the use of artificial viscosity will be completely avoided in this and subsequent chapters—simply due to a conviction that such methods should not really be necessary for the treatment of scalar, electromagnetic, or gravitational waves propagating in a vacuum. All of the schemes to be studied in this chapter are constructed on the staggered mesh depicted in Fig. 3.2. The use of such a staggered grid structure is common practice in numerical relativity, as well as other areas of computational physics and is apparently originally due to von Neumann [74]. This structure is somewhat useful for constructing centred schemes for second order equations written in first order form, and as will be seen later, can also simplify the treatment of boundaries, particularly in non-Cartesian coordinates. However, it is certainly not necessary to 40 X X X >(C A A A k X X X i A A A X X X X A A A X h X X X Figure 3.2. V o n N e u m a n n staggered g r i d s t r u c t u r e . use a staggered s t r u c t u r e , p a r t i c u l a r l y for o n e - d i m e n s i o n a l p r o b l e m s , a n d its use in this a n d the next two c h a p t e r s is largely a m a t t e r of c o n v e n t i o n . T h e discrete a p p r o x i m a t i o n to $ is defined at the " x " locations a n d p a r t i c u l a r values are denoted i n the u s u a l fashion as $ . t . L i k e w i s e , the values II . 2 are defined only at the "A" l o c a t i o n s . X X Figure 3.3. Difference star for the s t a n d a r d scheme. C o n s i d e r for the m o m e n t , the s y s t e m (3.22-23) for the special case of v a n i s h -i n g 8 a n d c o n s t a n t c. T h e scheme A ' n " - 2 " = c A * * " , (3.24) J 3 ' 2 A' +»; + i=-cA*n; +* (3.25) built on the star of Fig. 3.3 is equivalent to what Richtmeyer and Morton term the "usual finite difference equation" for the wave equation [74]. Here it will be referred to as the standard scheme. (Note that because of the staggered grid, a reference to a star such as the one in Fig. 3.3 is to be understood to include the "conjugate" star x *-> A). The properties of this method are well known. In particular, it has an 0(h2 + k2) truncation error and is stable in the sense of von Neumann (App. A) provided a = ck/h < 1 (which is also the CFL condition (App. A) for the scheme). All of the methods examined here incorporate it in some way or another. A A X X X X A A A A 8 > o 3 < o Figure 3.4. Difference stars for the one-sided (upwind) scheme. The four basic schemes are distinguished by the manner in which the convec-tive terms (/?$)' and {8U.)' are treated and are labelled as such. The first, which will be called the one-sided (OS) or upwind (UP) scheme is built on the stars shown in Fig 3.4. Note that the star which is used at a particular point depends on the sign of 8 there. The difference equations for the case 8 > 0 are A ' + n ; - ^ = A ; (/HI);-* + A * ( c » ) ; + 1 (3.26) A V * ; + i = A + ( / j * ) ; + i + A;(cnj; +» (3.27) The usual technique of Taylor expanding the grid functions about the centers of the difference stars shows that this is an 0[h + k) explicit scheme due to fact that the differencing of the convective terms is only first order accurate. The use of a centred spatial difference applied at the retarded time level is precluded by the fact that such a technique produces a scheme which is unconditionally unstable for the purely convective equations (c = 0, /? = constant). Wilson's scheme [104], which is the most widely used method in 3+1 numerical relativity is partially based on the OS scheme. The use of upwind differencing to stabilize convective terms is common practice in computational physics and rather sophisticated variants of the method which may be useful for numerical relativity are still being devised and studied [34,37,97,98]. A A A X X A A A Figure 3.5. Difference star for the Crank-Nicholson scheme. The second method, called the Crank-Nicholson (CN) scheme, is based on the star in Fig. 3.5. The difference equations are A+nr* =-Ka5pn>r*+as(wr*)+A* {c*)*+h ( 3 . 2 8 ) = \ (A5 + A0* + K (cll)n. + * (3.29) This scheme is implicit and has an 0(h2 + k2) truncation error. Its name de-rives from the resemblance of the star to that used in the original Crank-Nicholson scheme introduced for the heat equation [74]. For a purely convective equation with c o n s t a n t coefficients, it is easy to d e m o n s t r a t e t h a t the m e t h o d is u n c o n d i t i o n a l l y s table. B e c a u s e of the i m p l i c i t n a t u r e of the scheme, al l a d v a n c e d values of the e v o l v e d g r i d f u n c t i o n m u s t be d e t e r m i n e d s i m u l t a n e o u s l y . F o r the o n e - d i m e n s i o n a l s y s t e m s s t u d i e d here, the result ing linear s y s t e m s wil l be t r i - d i a g o n a l . S u c h s y s t e m s are soluble i n 0[N) t i m e (as m e a s u r e d in u n i t s of f u n d a m e n t a l c o m p u t e r o p e r a -t ions) where N is the n u m b e r of u n k n o w n s [62,74], so the scheme is c o m p u t a t i o n a l l y a c c e p t a b l e . S o m e a d d i t i o n a l work is needed to d e t e r m i n e w h e t h e r this scheme c a n be e x t e n d e d t o m u l t i - d i m e n s i o n a l p r o b l e m s in a c o m p u t a t i o n a l l y efficient m a n n e r . A A A A X X X X A A A A F i g u r e 3.6. Difference stars for the angled scheme. T h e t h i r d m e t h o d uses a n idea suggested by R o b e r t s a n d Weiss [75]. It is bui l t on the stars s h o w n in F i g 3.6 and is cal led the angled ( A ) scheme for o b v i o u s reasons. T h e difference equations for the first of the stars i n the figure are A+nr*=\ (a- m)Th+a+ (^ r 1 ) + ( 3 - 3 o ) A + * ; + I = \ (/?*>;:[ + A * (/?*);_,) + A * (en);+* (3.31) L i k e the C N s c h e m e , these equations are i m p l i c i t w i t h 0(h2 + k2) t r u n c a t i o n error. T h e a d v a n t a g e of this scheme is t h a t it is effectively e x p l i c i t — g r i d function values m a y be successively u p d a t e d in the d i r e c t i o n of i n c r e a s i n g (decreasing) j once the first (last) value has been specified. A stabil i ty analysis of the scheme for the case of a simple convection equation indicates stability for all values of a. Roberts and Weiss note, however, that this result essentially assumes that the boundary conditions are set exactly; if they are not, then the errors which are propagated in the sweeping direction will be amplified if the CFL condition is violated and the scheme will be unstable. This point will be discussed in greater detail below. X X X X Figure 3.7. Difference star for the leap-frog scheme. The final method to be studied is the leap-frog (LF) scheme, which is an explicit scheme based on the star in Fig. 3.7. The difference equations are A ^ = AQ [W)]+k + \ ( A * ( C » ) £ [ + Al ( e * ) ^ ) (3.32) A < * ; + , = A g [0*)n.+ i + 1 ( A * (ell) ^ + A ; (en);"*) (3.33) Again, this system has an 0(h2 + k2) truncation error, but unlike the other methods, it is a 2-level scheme. Leap-frog schemes for hyperbolic systems have been the subject of a great deal of study and have been advocated for problems where the conservation of the energy of a solution over long periods of time is important [43,44]. (There are references in the numerical relativity literature to "staggered leap-frog schemes". In the nomenclature of this chapter, these methods are all variants of the standard method. In particular, it is not clear that the convective terms have ever been treated in a true leap-frog fashion.) Note that when /? = 0, the above equations reduce to a system which is equivalent to two steps of the standard scheme. A principal advantage of this type of differencing is that it can easily be extended to multi-dimensional problems with little compromise of computational efficiency 3.4 Von Neumann Stability Analysis. As mentioned in Appendix A, there are a wide variety of techniques available for the stability analysis of finite difference equations. However, for pure initial value problems, and particularly for hyperbolic systems, the method of von Neu-mann, which involves the freezing of coefficients combined with Fourier analysis continues to be the most valuable tool for obtaining practically useful stability criteria [50,74]. In this section, the conditions for von Neumann stability for the four schemes described above are determined. For two of the schemes, this is ac-complished analytically, for the other two methods, a straight-forward numerical technique is used. It should be kept in mind that the conditions derived are only necessary criteria and may not be sufficient to guarantee stability. Some theorems from [53] will be of considerable use in the following. Following Miller, a polynomial n p(z) = a,* 1 o„ z e C 1 (3.34) is called von iVeumann iff its roots m satisfy \m\ < 1 t = l,...» (3.35) and conservative von Neumann iff \m\ = l » = l,...n (3.36) Introducing the notation VN(p) EE p is a von Neumann polynomial. CVN(p) EE p is a conservative von Neumann polynomial. Theorems 6.1 and 6.4 of [53] may be stated as VN(p) if and only if either CVN(p) or (|p*(0)| > |p(0)|) and V N ( P l ) (3.37) CVN(p) if and only if P I E E O and VN(p') (3.38) wh ere P*iz) = ]Ca»-'*'' (3-39) i=0 • M - £ (3.4.) and an overbar denotes complex conjugation, p* (z) is known as the reciprocal polynomial of p(z), and Pi{z), which Miller terms the reduced polynomial, is similar to the Schur transform of p(z) [39]: n-1 T P ( 2 ) = ( a° a« ~ °n an-») 2* i=0 Repeated application of the Schur transform forms the basis of the Schur-Cohn algorithm [39] which can be used to determine whether all zeros of a polynomial lie outside of s o m e closed disk in the c o m p l e x plane. H o w e v e r , this a l g o r i t h m m u s t b e m o d i f i e d w h e n the p o l y n o m i a l has zeros t h a t lie on the u n i t c ircle , a n d this is precisely w h a t M i l l e r has done. 3.4.1 Analytic Stability Criterion for the CN Scheme. T r e a t i n g c a n d /3 as c o n s t a n t s , a n d defining k k oi = c- o2 = f3- (3.42) a fi (3.28-29) m a y be r e w r i t t e n as n" + i = n"-i + °-£ (A;n"+* + Ajn""*) + (J.4S) A p p l y i n g t h e F o u r i e r t r a n s f o r m technique d e s c r i b e d i n A p p e n d i x A , the a m p l i f i c a -t i o n m a t r i x for this s y s t e m is G U ) = ( 2 ° l { e i i " , l ) l , f ) (3.45) y ' \o\z{\ - e~ l f) /w z-oi2s/w2 J K ' w = 1 + 2t'<72 sin £ s = 4 s i n 2 - (3.46) 2 a n d £ = ujh, where u is the F o u r i e r t r a n s f o r m v a r i a b l e . T h e eigenvalues of G ( £ ) a r e t h e n the r o o t s of the c h a r a c t e r i s t i c p o l y n o m i a l p{n) = n2 - (2z - sox2w~2) p + z2 (3.47) a n d for a g i v e n values of the p a r a m e t e r s o\, o~i, the scheme satisfies the v o n N e u -m a n n s t a b i l i t y c r i t e r i a if < 1, for all £, |£| < it. p(fi) is of the general w h e r e w z = — w f o r m 48 p{fi) = ft2 + A\i + B A = S ~ - 2z B = z2 (3.48) (3.49) (3.50) w h e r e A a n d B are, in general , c o m p l e x . T h e n p*{fi) = Bp,2 + Afi+ 1 (3.51) a n d = (l - |B|2) n + A - AB (3.52) N o w clearly, \B\ = \z2\ =• 1 a n d it is easy to check t h a t A = AB. T h u s the r e d u c e d p o l y n o m i a l , px, identically vanishes a n d by (3.38), C V N ( p ) i f a n d only if V N ( p ' ) . N o w whose r o o t fit is p'(li) = 2fi+ - 2z sox 2\w\2 - soxJ 2w2 ~ 2w2 (3.53) (3.54) T h u s \flt\ < 1 o so: 2 w < 1 sol < 4IW|2 =• s i n 2 | o 2 < ( l + Aa\ s i n 2 f) (3.55) w h i c h holds for all £, |£| < w i f a n d only if ax < 1 (3.56) N o t e t h a t this is just the C F L c o n d i t i o n for the "explicit p a r t " of the scheme. T h u s , for a r b i t r a r y 0 2 , the C N amplif ication m a t r i x will have eigenvalues of precisely unit m o d u l u s — s o the scheme is n o n - d i s s i p a t i v e — p r o v i d e d c\ < 1. 3 . 4 . 2 A n a l y t i c S t a b i l i t y C r i t e r i a f o r t h e A n g l e d S c h e m e . U s i n g the same definitions as p r e v i o u s l y , the c o n s t a n t coefficient version of (3.30-31) is n „ + L = n „ _ l ^ / n „ + L n _ L \ „ J 3 2 \ ~ J J / J + li v ' T h e a m p l i f i c a t i o n m a t r i x c a n again be w r i t t e n in the f o r m of (3.45), w i t h the single change (3.59) + f ( , < - , ) W i t h this r e p l a c e m e n t , the calculation proceeds exactly as before. T h e r e d u c e d p o l y n o m i a l vanishes identical ly , a n d d e m a n d i n g that the r o o t s of the derivative of t h e characterist ic p o l y n o m i a l lie w i t h i n or on the unit circle leads to C V N ( p ) s a x 2 < 4\w\2 2 • 2 £ a , sin'' - < 1 + ^ (e*f - l) 02 2 => of s i n 2 £ < 1 + {o2 2 - 2 e r 2 ) . 2e s in — ox2 < (1 - a 2 ) 5 (3.60) Thus, the von Neumann stability condition for this scheme is <J\ < 1 - o-i if o2 < 1 (3.61) °\ < o2 - 1 if o2 > 1 (3.62) One can make some sense of these criteria by noting that the characteristic speeds for the differential system are — 3 — c and — 8 + c with the former corre-sponding to the left moving signal so that (3.61), rewritten as o\ + °2 < 1 (3.63) can be interpreted as a CFL condition on the scheme. (The computational domain of dependence of H*+* does not contain n^*J for any k > 1, for example.) Note that, in contrast to the CN scheme, this method is von Neumann stable for cases where o\ > 1. Presumably, this has something to do with the fact that the stars for the two parts of the scheme share common mesh points here, wheras for the CN method, they are "decoupled". The second part of the condition, (3.62), is harder to explain from this point of view. However, following Roberts and Weiss, note that (3.30) can be rewritten as J } r \ J + l J - l / J - j V 7 P = —,—~—TT T = °LH . (3.65) 2 ( l -a 2 / 2 ) l - a 2 / 2 1 ' with an analagous expression for Now as observed by the above authors in their study of a similar equation, if \p\ > 1, which occurs in this case when o2 > 1, then any errors introduced in setting grid function values (boundary values, for example) will be magnified as the update sweep proceeds. Thus the scheme should be expected to be unstable for a2 > 1, which has been confirmed by the numerical 51 e x p e r i m e n t s to be d e s c r i b e d below. F i n a l l y , note that the angled scheme is also n o n - d i s s i p a t i v e . 3.5 Numerical Determination of Stability Domains for all Four Schemes. T h e u p w i n d a n d leap-frog schemes p r o v e d more difficult to treat a n d a n a l y t i c c o n d i t i o n s for their v o n N e u m a n n stabi l i ty have not been d e r i v e d . H o w e v e r , since b o t h schemes are explicit , a reasonable hypothesis is that they are b o t h C F L l i m i t e d . A s m e n t i o n e d previously , the characterist ic direct ions of (3.22-23) are — /3 — c a n d — 0 + c. S i n c e c > 0, the largest (in a b s o l u t e value) signal speed is given by |/?| + c. T h u s for b o t h schemes, one m i g h t expect the stabil i ty c r i t e r i a to be °\ + |*21 < 1 (3.66) F i g s . 3.8-11 show the results of a n u m e r i c a l d e t e r m i n a t i o n of the stabil i ty d o m a i n s for al l four of the schemes. E a c h dot represents a v o n N e u m a n n stable choice of o\ a n d o i as d e t e r m i n e d by a n explicit c a l c u l a t i o n of the s p e c t r a l radii of the c o r r e s p o n d i n g a m p l i f i c a t i o n matrices for several (« 10 to 20) values of £ r a n g i n g f r o m 0 to 7r [62]. A p o i n t was deemed unstable if any eigenvalue w i t h m o d u l u s greater t h a n 1 + e a p p e a r e d , t was kept very s m a l l a n d was used to d e t e r m i n e w h e t h e r the p o i n t s on the b o u n d a r y of the stabil ity d o m a i n were m a r g i n a l l y stable. T h e s e figures c o n f i r m t h a t the explicit schemes are C F L l i m i t e d a n d also r e p r o d u c e t h e c o n d i t i o n s (3.56) a n d (3.61-62). It s h o u l d also be m e n t i o n e d here t h a t the r e d u c e d p o l y n o m i a l for the leap-frog scheme also vanishes, so that the m e t h o d is n o n - d i s s i p a t i v e whenever it is stable. F i g s . 3.12-15 i l lustrate the results of another set of n u m e r i c a l tests where each d o t now represents an a c t u a l e v o l u t i o n of some initial d a t a u s i n g the difference s y s t e m s d e s c r i b e d a b o v e . F o r each r u n , the spatial g r i d h a d 32 p o i n t s a n d 10 t o 30 t i m e steps were t a k e n . A p o i n t was considered stable if the n o r m of o 0 0 S l -Figure 3 . 8 . V o n N e u m a n n stabil i ty d o m a i n of the one-sided (upwind) scheme. a> E Q> o _C o (/I c o ' ' ' J-'J* vt o - C o • i c o u S I O'O S l -E a> —\ o in c CL ID Figure 3 . 9 . V o n N e u m a n n s t a b i l i t y d o m a i n of the C r a n k - N i c h o l s o n scheme. ,;;!! • • fc> Angled Scheme 1 0 0 SH Figure 3.10. V o n N e u m a n n stabi l i ty d o m a i n of the angled scheme. o E x : o D> O k_ i 1 Q. O <D _ J SH 0"0 S H -Figure 3.11. V o n N e u m a n n stabil i ty d o m a i n of the leap-frog scheme. 54 o a> Q t —: o tn X) c Q-Z> SI 00 SH-F i g u r e 3.12. P r a c t i c a l s t a b i l i t y d o m a i n of the one-sided (upwind) scheme. a> E 0> o -C o to c o b o o z 1 1 c o c_> SI 00 S ' l -ZD F i g u r e 3.13. P r a c t i c a l s t a b i l i t y d o m a i n of the C r a n k - N i c h o l s o n scheme. 55 o P tT Angled Scher SI 00 SH-Figure 3.14. Practical stability domain of the angled scheme. a> o E a> _c o 60 u 1 > 1 CL 0 a> S'l 00 S"l-ZD Figure 3.15. Practical stability domain of the leap-frog scheme. t h e final solution was less by some factor, 6, t h a n the n o r m of the init ial d a t a . (Il/jlloo = maxj \fj\). 6 was 2.0 for the three n o n - d i s s i p a t i v e scemes a n d 1.01 for t h e u p w i n d scheme w h i c h , as wil l be seen, dissipates rather badly. N o t e that the o n l y s ignif icant change f r o m the p r e v i o u s d i a g r a m s is the absence of a n y region of s t a b i l i t y for the angled scheme w h e n ai > 1, as was e x p e c t e d . A g a i n , this is a n o t h e r r e m i n d e r t h a t the von N e u m a n n c o n d i t i o n is, i n general , only necessary for s t a b i l i t y . 3.6 Convergence Tests for the Four Schemes. R e c a l l f r o m Sec. 3.3 t h a t the one sided scheme has an 0(h + k) t r u n c a t i o n er-r o r w h i l e the other schemes are 0[h2 + k2) accurate. In general, w i t h o u t a stronger n o t i o n of s t a b i l i t y t h a n is b e i n g used here, it is not a t r i v i a l m a t t e r to p r o d u c e r i g -o r o u s results c o n c e r n i n g the precise rate of convergence of the solution of difference e q u a t i o n s to the exact solution f r o m t r u n c a t i o n error estimates. S u c h results wil l n o t be p u r s u e d here, a l t h o u g h they c o u l d u n d o u b t e d l y be o b t a i n e d — a t h o r o u g h d i s c u s s i o n of the techniques i n v o l v e d is given by K r e i s s a n d O l i g e r [44]. In practice, w h e n d e a l i n g w i t h p r o p e r l y c o n s t r u c t e d schemes a n d sufficiently s m o o t h funct ions, one comes to expect t h a t the solution error will vanish as r a p i d l y as the t r u n c a t i o n e r r o r . T h i s h y p o t h e s i s , w h i c h w i l l be discussed in greater detail in C h a p t e r 5, c a n be verified n u m e r i c a l l y for the schemes d e s c r i b e d above by p e r f o r m i n g convergence tests. S u c h a test c a n be d e s c r i b e d schematically in the fol lowing way. A s s u m e t h a t a differential s y s t e m Lf = 0 (3.67) w h e r e L is some differential o p e r a t o r a n d / is the solution vector is discretized on s o m e m e s h s u c h t h a t the grid spacings in each of the c o o r d i n a t e directions are O(l) multiples of a basic increment h. Then the discrete solution / satisfies Lhfh = 0 (3.68) where L h is the difference approximation to L. Now, provided that the exact solution / is known, the convergence rate of the difference scheme can be tested by solving the problem (3.68) on grids with varying mesh sizes. In particular, assume that the discrete solution is calulated using grid spacings of h and 2h. Then, one can define a convergence factor, c ,^ as for some discrete norm || • ||. Then, if the error in the computed solution is 0{h), as would be expected from a scheme with O(h) truncation error, this quantity will approach 2 as h —+ 0. Similarly, if the solution error converges quadratically in the mesh spacing, c& will tend to 4. As described in Sec. 3.2, an exact solution for the system (3.22-23) can be generated by performing a coordinate transformation on a solution to the ordinary wave equation. The transformation chosen for the convergence tests, as well as for the experiments in subsequent sections, was x — X t = cx tanh(c2T) cos(c3A:) (3.70) where c\, c 2, and C3 are parameters which are chosen so that the quantities 3 and c which result assume reasonable values during an evolution. In addition, c3 is chosen so that periodic boundary conditions can be applied on the spatial boundaries of the grid. Graphs of the typical behavior of 3 and c are shown in Figs. 3.20-21. Each of the convergence tests used four grids which covered the spatial interval 0 < x < L. The mesh spacings used were hn = 2kn = 2~[*+n)L n = 0,1,2,3 (3.71) with L — 60.0, so the coarsest spatial increment was roughly 1 and the finest ~ - 1 * o about 1/8. The exact values rL 2 , and, in the case of the leap-frog scheme, 11^ . 2 , $ + L , were supplied as initial values. The tests ran for 2 + time steps except for the case of the one sided scheme which was advanced for 2 6 + ' steps. Every 2 steps the £2 norm of the error in IL 2 was computed and output. Once this process had been completed on all four grids for a given scheme, the error norms were combined to calculate the convergence factors c ^ , c 2 f t , and c 4 f t as a function of time. The results of these tests are shown in Figs. 3.16-19. In all of the graphs, the chain-dashed line indicates the expected asymptotic convergence factor. In the case of the upwind scheme, the convergence rate drops as a function of time, then levels off near the expected value of 2. This is due to the fact that 3 grows from an initial value of 0 everywhere, so that the second order part of the scheme dominates at early times. The results for the other schemes provide a good indication that the difference solutions are converging quadratically to the exact solution. The initial transient behavior which appears (most noticeably in the CN and LF schemes) is presumably due to the use of exact initial data for all runs. This is inconsistent with the O(h) or 0{h2) variation in the grid function which would be expected at t = 0 if the solution had been begun at an earlier time. In addition to demonstrating that the solutions of the various schemes can be expected to converge as quickly as the truncation error, this process of convergence testing was a very useful check of the routines 59 F i g u r e 3.16. C o n v e r g e n c e test of the one-sided (upwind) scheme. q o o c c o 8h : 4h 4h : 2h 2h : h Crank-Nicholson Scheme F i g u r e 3.17. C o n v e r g e n c e test of the C r a n k - N i c h o l s o n scheme. Angled Scheme X 0 S 10 15 * Figure 3.18. Convergence test of the angled scheme. 8h 4b 4h 2h ______ 2h h Leap-frog Scheme 0 5 10 15 « Figure 3.19. Convergence test of the leap-frog scheme. which were used to solve the equations. Incorrect code would usually produce a wrong (unexpected) rate of convergence. 3.7 Numerical Accuracy Tests of the Schemes Figs. 3.22-26 show the exact and calculated values of the grid function (at selected times) which were generated using the parameters ci — 2.0, c-i — 7r/l5 and C3 = 0.1 in the transformation (3.70). 64 grid points were used with L = 60.0, so the mesh was fairly coarse. The tests ran for 256 time steps with k — 1/8 h. The initial data was / W = / W = „ P ( - ( £ ^ ! > ) ' ) h{X) = h{x)=0 (3.72) which produces a single left moving wave packet that starts near the right edge of the grid. Figs. 3.20-21 show the functions 0 and c at selected times during the evolution. Several features of these plots are noteworthy. First, there is considerable variation in the (exact) function II . 2 as the pulse propagates in the transformed coordinates which indicates that the difference schemes are being subjected to a non-trivial test. The dissipative nature of the one sided scheme is very evident from Fig. 3.23. Although the basic shape of the waveform is reproduced reasonably well, much of the energy in the pulse has been lost by the end of the run. The last time slice plotted also shows a definite elongation of the wave which is mainly a dispersive effect—the higher frequency components of the solution travel too slowly. The results for the second order schemes are clearly niuch better, and virtually identical. Deviations between the different numerical solutions are difficult to detect on the scale of Figs. 3.24-26. There is very little evidence of dissipation as the flattening 62 Figure 3.21. T y p i c a l c generated by c o o r d i n a t e the t r a n s f o r m a t i o n (3.70) Figure 3.22. Time series of the exact solution II. 2 for accuracy tests of difference schemes. Figure 3.24. Time series of IIn+2 generated by the Crank-Nicholson scheme. 65 Figure 3.25. Time series of IT . 2 generated by the angled scheme. Figure 3.26. Time series of IT. 2 generated by the leap-frog scheme. of the peak of the latest curve is accompanied by a deepening of the trough in all three graphs. Again, this is indicative of a dispersive effect. Another method of comparing the difference solutions with each other and with the exact solution is to compute (fast) Fourier transforms (FFT's) [56] of the various grid functions to see what is happening with the discrete Fourier amplitudes and phases. Fig. 3.27 shows a plot of amplitude versus wavenumber for the exact and three sets of computed values of 11^. 2 using the data from the runs described above. Again, although the one sided scheme reproduces the basic shape of the spectrum, the attenuation of even the low frequency components is obvious. The second order schemes are much better and produce very similar results. The spectrum generated by the angled scheme is not included in the figure since it cannot be distinguished from the C N scheme. To qualitatively illustrate the dispersive effects in the schemes, another run was performed with the same initial data on twice as fine a grid (128 x 512). The phase of each Fourier component was accumulated and at the end of the run, for each wave number q, the quantity — ; — r ^ — (3-73) \gexact v ' 1 was calculated, the logarithm of this quantity for the four different methods is shown in Fig. 3.28. Only the lower half of the wavenumbers are included in the plot, since the results for the upper half are unreliable, and as can be seen from the previous figure, there is very little power in the upper half of the spectrum. The one sided scheme is clearly more dispersive than the second order schemes which, again, exhibit very similar behavior. None of the methods does a particularly good job at high frequencies; as noted by Roberts and Weiss, this is the most serious problem which remains with any centred, second order, non-dissipative scheme o Wave number Figure 3.27. Comparison of Fourier amplitudes of exact solution and three of the difference solutions. O 0 4 8 12 16 20 24 28 32 Wave number Figure 3.28. Comparison of accumulated phase error in the four difference schemes. In s u m m a r y , all of the three second o r d e r schemes a p p e a r to have the p r o p e r -ties w h i c h are required to h a n d l e a r a d i a t i v e p r o b l e m in general c o o r d i n a t e s , at least for t h e case of a single spatial d i m e n s i o n . T h e m e t h o d s have reasonable stabi l i ty c r i t e r i a , a low degree of d i s s i p a t i o n ; fairly low d i s p e r s i o n except at high frequen-cies, a n d are d e m o n s t r a b l y second order a c c u r a t e in b o t h space a n d t i m e . T h e u p w i n d s c h e m e is clearly i n a d e q u a t e in c o m p a r i s o n . A l l m e t h o d s require 0[rig) c o m p u t a t i o n a l t i m e , where ng is the n u m b e r of p o i n t s i n the s p a c e t i m e g r i d , a n d are therefore c o m p u t a t i o n a l l y acceptable. F i n a l l y as m e n t i o n e d previously , because the s c h e m e s are designed for the wave e q u a t i o n in general c o o r d i n a t e s , they c a n be e x p e c t e d to work for the self -gravitat ing field. 3.8 Numerical Radiation Boundary Conditions. T h i s sect ion, w h i c h is essentially i n d e p e n d e n t of t h e rest of the c h a p t e r , deals w i t h the issue of i m p o s i n g a p p r o x i m a t e outer b o u n d a r y condit ions in the n u m e r i c a l t r e a t m e n t of a r a d i a t i v e p r o b l e m . A t t e n t i o n will be r e s t r i c t e d to the case of spherical s y m m e t r y , a n d the wave e q u a t i o n in the u s u a l c o o r d i n a t e s (r, t) c2 4>tt= -j[r2<Pr)r (3.74) w i l l be a d o p t e d as a m o d e l . ( N o t e t h a t c is a c o n s t a n t here.) M a t h e m a t i c a l l y , the c o n s t r u c t i o n of a spacetime using, for e x a m p l e , the 3+1 f o r m a l i s m is a pure init ial -value ( C a u c h y ) p r o b l e m . In other words, the spacelike h y p e r s u r f a c e s i n v a r i a b l y have no outer edges, so there is no need to i m p o s e outer c o n d i t i o n s o n a n y r a d i a t i v e q u a n t i t y , be it scalar, e l e c t r o m a g n e t i c or g r a v i t a t i o n a l r a d i a t i o n . N u m e r i c a l l y , of course, there wil l always be an outer edge of the n u m e r i c a l d o m a i n unless some c o m p a c t i f i c a t i o n scheme is u s e d . A s s u m i n g , as is the case here, t h a t c e n t r e d differences are b e i n g used, it is necessary to provide b o u n d a r y c o n d i t i o n s or modify the difference equations at, or near, the outer edge. C o n s i d e r the requirement that the amount of radiation which enters the region covered by the mesh is negligibly small. This is a commonly used assumption in numerical relativity which is probably physically reasonable in most cases. Then, given that the desired solution has no incoming radiation, the goal is to ensure that the treatment of the outer numerical boundary does not generate such radiation and that the amount of reflection of an outgoing pulse off the boundary is minimized. For general multi-dimensional radiative systems this is quite a difficult prob-lem and has received considerable attention (see for example, [41] and the references contained therein). Fortunately, for the case of spherically symmetric waves it is relatively easy to come up with acceptable boundary conditions. It should be noted that (numerical) boundary conditions have to be taken into account when assess-ing the overall stability of a difference scheme. The stability theory for the mixed initial-boundary value problem is fairly well-developed (see [32,42] for some early work and [43,44] for more complete discussions), but, as mentioned in App. A, is generally more difficult to apply than von Neumann analysis. Roughly speaking, the stability analysis involves looking for "normal modes" of the complete set of difference equations (boundary conditions included) and establishing that none of these modes grow "too rapidly" with time. This type of analysis has not been per-formed in this thesis, but empirical evidence suggested that the boundary conditions described below did not cause instabilities. The general incoming and outgoing solutions to (3.74) are 4>m — F{r + ct) (3.75) r 4>out — G{r - ct) (3.76) r where F and G are arbitrary functions. If the outer edge of the grid is at r = R, then the demand that H)t + cM r =0 (3.77) is equivalent to the statement that no radiation enters the grid from r > R. This will be called an outgoing radiation boundary condition. This equation can then be finite differenced to provide conditions that can be used numerically. Since (3.77) is an exact expression for the case of spherically symmetric waves, its effectiveness in a numerical scheme is limited only by the order to which it is differenced. This assumes, of course, that the overall scheme remains stable. Follow-ing Piran [68], one possible expression may be derived directly from the expression (3.76) for an outgoing wave. Let 6 — ck — c At. Then Tj4>nj+l = R4>{R, t + At) = {R - 6)<f>{R - 6, t) (3.78) <f>nj+1 = (1 - 6/R) {9{R,t) - 69r{R,t)+0{62)) Hil-S/rjUfi-Sbl+j) (3-79) Because a first order backwards difference is used to approximate (f>r(R,t), this is a first order condition. It is the spherical wave equivalent of the expression used by Piran in his study of cylindrical systems, and is also very similar to the conditions that have been used in some recent axisymmetric numerical relativity codes [28]. A second order condition which can be used in conjunction with some of the interior schemes described previously can be derived by discretizing (3.77) using second order backward differences. At r = R (3.80) A L ( 3 ^ + 1 - + O(At') = -c (^f- + A r _ ( 3 ^ + 1 -+ 0(Ar 2) (3.81) so (3.82) Note that the value Sn is required. This presented no problem in the programs constructed in this study since two full time levels were always maintained. Another method which can be used to minimize the effects of wave reflection off an artificial boundary involves the introduction of terms to the differential system which damp, or filter the waves near the boundary. Although such a technique will not be used in subsequent chapters, a brief description of one method is included here because it may be useful in multidimensional numerical relativity calculations. The idea, due to Israeli and Orszag [41], is called a sponge filter . Equation (3.74) is modified by the introduction of a term which allows outgoing waves to propagate while filtering incoming radiation: <t>tt = £ ( r 8 * r ) r - u(r) ((r6)r + (r<f>)t) (3.83) The Newtonian cooling coefficient, v[r) > 0, is chosen to vanish in the region 0 < r < r 0 , so the width of the filter is R - r0. In the interval r0 < r < R, the authors suggest the choice u{r) = A{r - rf)n{R - r){R - rf)-n-2(n + 1)(» + 2) (3.84) The constants, A and n can be used to "tune" the filter and [41] contains a plot of the relative performance on a model problem for various values of the parameters. In order to evaluate the effectiveness of these various boundary conditions, another program which solves (3.83) using finite differences was written. Again, introducing the auxiliary variables $ = Sr, II = <pt, (3.83) is rewritten as $t = FIr (3.85) nt = l ( r 2 $ ) r - i / ( ^ + r($ + n)) (3.86) and discretized as - + ( + M+n;-*)) (3.88) A^; = n;n-" (3.89) The intitial data was of the form ^(r,O) = <£ 0exp^- ( !L^£) j (3.90) n(r,0) = 0 (3.91) which produces one ingoing and one outgoing pulse. Some of the outgoing pulse gets reflected off the outer boundary and propagates inward. By monitoring 6{0,t), for example, at later times, a measure of the amount of reflection that has occurred is obtained. Such results are shown in Figs. 3.29. For the first two plots on each graph, the function u was identically 0 and the outer values of <f> were updated using (3.79) and (3.82) respectively. For the sponge filter runs, a cooling coefficient of the form (3.84) was used with A — 1.0 and n — 2. The other parameters for all runs were R - 25.0, ry - 12.5, 60 = 0.05, r0 - 10.0, A - 2.0 and h = Ar = 2A<. Note that the vertical scale of the first plot is twice that of the second. In halving the mesh spacing, the maximum peak to trough variations in the reflected waves Figure 3.29. Comparison of various numerical outgoing boundary conditions for the spherical wave equation. are reduced by factors of very close to 2 and 4 using the first and second order conditions, respectively. It is more difficult to quantitatively assess the convergence rate for the case of the sponge filter, but the effectiveness of the technique is clear. Finally, it is estimated that an amplitude of ±0.05 in these plots represents about 1% reflectance of the outgoing wave. 75 C H A P T E R 4 S o m e C o n s t r a i n e d Solutions of the S e l f - G r a v i t a t i n g Scalar F i e l d 4.1 Introduction. As discussed previously, a major motivating factor for the numerical study of a general relativistic system with one spatial degree of freedom is the opportunity it presents for experimentation with the method used to find discrete solutions of the equations. This fact has been recognized and used to advantage by several researchers since the advent of numerical relativity in the early 1970's. The most notable of the one dimensional studies is Piran's work on cylindrical general rel-ativistic systems [68]. From the point of view of a researcher who is new to the field, this work is very useful in that it clearly describes the basic considerations that must go into any numerical scheme based on the 3 + 1 formalism. In terms of "experimentation", Piran solved the cylindrically symmetric vacuum Einstein equa-tions using several coordinate systems and various combinations of the evolution and constraint equations for the geometric quantities, and, to a certain extent, was able to verify the results with analytic solutions. Among other things, the results demonstrated the efficacy of the basic numerical evolution scheme which had origi-nally been developed for problems in relativistic hydrodynamics by Wilson [104]. The work described in this chapter is very much in the spirit of Piran's study.-In particular, numerical solutions of the model described in Chapter 2 are investi-gated in three different coordinate systems and evidence for the mutual consistency of the results is presented. In addition, weaknesses and limitations in the numerical techniques used are identified and discussed. For the most part, evolution equations for the geometric variables are not used in this chapter. The applications of versions of the schemes devised in the previous chapter to the geometrical quantities will be discussed in the next two chapters. 4.2 Coordinate Conditions. The question of the choice of coordinates has always been a central issue in numerical relativity and the literature contains several thorough discussions of the topic, particularly for the case of axisymmetric spacetimes [2,70,81]. Here, some of the main conclusions of these studies which are pertinent to the current work are briefly reviewed and the various coordinate conditions employed in this thesis are then described. It is generally accepted that the specification of the lapse function is the most crucial part of choosing a coordinate system for the purposes of constructing a spacetime, particularly when the gravitational field is expected to be strong. Work by Penrose, Hawking and Geroch [36] in the 1960's clearly established the fact that singularities are a generic feature of solutions of Einstein's equations. In particular, any spacetime which contains one or more black holes will contain singularities which can develop from perfectly regular initial data—such as the data representing a sufficiently massive star prior to its collapse. The problem is that once a time slice of a numerically generated spacetime intersects a singularity, the evolution must necessarily cease (with current techniques). This may leave some interesting and physically relevant portions of the spacetime (part of the exterior of a black hole, for example) "unexplored". In order to avoid this scenario, a great deal of effort has gone into the study of slicing conditions (choices of the lapse function) which "avoid singularities". A singularity is avoided by having the lapse function tend to zero sufficiently rapidly 77 in the vicinity of the singularity, which has the effect of "freezing" the evolution there. The best known, and certainly the most widely used, such condition is called maxima/ slicing [23,27,57-61,66,68,76,77,105]. Here one demands that the trace of the extrinsic curvature vanish on all slices of the foliation As will be seen below, this requirement leads to an "elliptic" equation which deter-mines the lapse on each hypersurface. The name "maximal slicing" derives from the fact that such a slice has maxi-mal volume with respect to variations of the surface within the enveloping spacetime. This condition also has the advantage that it simplifies the 3+1 equations and can be used to eliminate one of the diagonal extrinsic curvature components from the equations. A second slicing condition, which also appears to have the singularity avoid-ance property, has been recently advocated by Bardeen and Piran [2] and imple-mented by Stark and Piran [70,71,87,88] in the context of axisymmetric stellar collapse. Again, the slicing is defined through a condition on K: This relation leads to a parabolic equation which the lapse must satisfy on each slice. Bardeen and Piran refer to the slices generated by this condition as polar hypersurfaces. As with maximal slicing, (4.2) can be used to eliminate one of the K%j, and when combined with an appropriate choice of spatial coordinates, polar slicing can lead to a considerable simplification of the 3+1 equations. Although some effort has been expended in the investigation of spatial coor-dinate choices which have some special geometric (and possibly physical) signifi-cance [9,85,86,10 5], the spatial degrees of coordinate freedom are most often used K = gijK11 = 0 (4.1) K = K\ (4.2) to produce simplification in the system of equations to be solved. The most obvious possibility is to use a zero shift vector which results in spatial coordinate trajectories which are normal to the hypersurfaces. Such coordinates eliminate the convective derivative terms which appear in the evolution equations (Sec. 2.6), and thus sim-plify the task of producing stable difference schemes which are properly centred in time. This condition has been successfully employed in many numerical relativity codes to date [10-12,24,26,5 7-61,68,82]. Another general approach to the choice of spatial coordinates, due to Wil-son [104], involves imposing algebraic relationships among some of the components of the 3-metric. Such conditions will generally lead to parabolic or elliptic equa-tions which must be satisfied by the components of the shift vector on each slice. Care must be taken to ensure that the conditions are compatible with the slicing choice and with other considerations such as regularity at r = 0 in spherical polar coordinates. The two such conditions which are used in this thesis are o = b (4.3) and 6 = 1 (4.4) (refer to (2.26) for the definition of c and b). The first of these has become known as the isothermal condition, and has proved quite useful in spherically symmetric collapse calculations [104,76,7 7]. The second choice is the specialization to spherical symmetry of the radial condition, introduced by Bardeen and Piran [2]. In this case 2-spheres at a radius r have a proper surface area 4nr2, so that the radial coordinate has a direct geometrical significance. Both of these conditions reduce the complexity of the equations derived in Chapter 2 and eliminate at least one variable from the set [a,b,KTr, K9g, /?}. Table I summarizes the three coordinate systems that will be used in the following study. These systems will be referred to as the maximal/normal, maximal/isothermal a n d polar/radial systems, respectively, T h e v a r i o u s differential e q u a t i o n s w h i c h the lapse a n d shift satisfy in these c o o r d i n a t e s y s t e m s are d e r i v e d i n the fol lowing subsect ions. S ection Slicing Spatial Coordinate Condition Condition 4.6 maximal normal [8 — 0) 4.8 maximal isothermal (a = 6) 4.10 polar radial [b = 1) T a b l e I. T h e three c o o r d i n a t e s y s t e m s used i n the s t u d y of the self-grav-i t a t i n g scalar field. 4.2.1 T h e M a x i m a l S l i c i n g E q u a t i o n . A s s u m i n g t h a t the init ia l surface has K = 0, a m a x i m a l foliation is generated b y d e m a n d i n g t h a t K = 0 (4.5) T a k i n g the trace of (2.14), this e q u a t i o n b e c o m e s K = K\ = 8K'- alili + a[R +K2+ 4n{S-3p)) =0 (4.6) N o w a n d f r o m (2.44-45) *2 Using K — K' — 0 along with the last two expressions in (4.6) yields a second order homogenous equation for the lapse function It has been frequently observed that the scalar curvature, R, generally requires many operations to evaluate numerically so that it is usually advantageous to solve for R from the Hamiltonian constraint (2.7), and substitute the result into the maximal slicing equation [70,76,77]. Using the fact that K99 = ~^Krr in a maximally sliced spacetime, (4.9) then becomes For a spherically symmetric spacetime, the difference in numerical efficiency which results from using (4.10) instead of (4.9) is not very significant, and (4.9) was most frequently used in this work since the scalar curvature was calculated for diagnostic purposes anyways. 4.2.2 T h e P o l a r S l i c i n g E q u a t i o n . Since, for the case of spherical symmetry, K = KTT + 2K9g, the polar slicing condition K = KrT implies that K9g = 0 everywhere. Again, assuming that this holds on the initial hypersurface, the condition (4.9) a = 0 (4.10) K% = 0 t > 0 (4.11) must be imposed to generate the polar foliation. Using (2.54) and K9g = K9/ = 0, this may be written as (4.12) 81 which is a first order homogenous equation for the lapse function. As shown in Table I, this slicing is only used in conjunction with the spatial coordinate choice b — 1, in which case (4.12) becomes a ' + ( 1 : r ! - 7 ) a = 0 ( 4 I 3 ) 4.2.3 The Isothermal Equation for /?. The isothermal condition is imposed by choosing initial data such that a — b and then demanding d = b t>0 (4.14) Identifying the right hand sides of equations (2.50) and (2.51), and using a — b, a' — b' gives the following first order equation for the shift vector component /? fi'-£ = r(P)' = a(K'r-K't) (4.15) Since this coordinate choice is used in combination with maximal slicing, the equa-tion may be further simplified to r(j)'= \«K\ (4.16) It should be clear that all three of the equations derived above for a and /? must be supplemented by boundary conditions. These will be discussed in the following descriptions of the various numerical schemes used to solve the equations. 4.3 The Structure of the Finite Difference Grid. As in the previous chapter, the numerical schemes described below are based on the staggered mesh structure depicted in Fig. 4.1. Assume there are J + l grid points in the radial direction. Then the numerical domain G may be defined as G 0 U G A U G x where G 0 = { (Jj - i) Ar , n At ) | j = 0,..., J; n = 0,. .. } G A = { ( ( ' " A r , (n- ^ A t ) | / = 0 , . . . , J ; n = 0,... } G x = {( jAr,nAt) | j = 0,..., J - 1; n = 0,... } The various grid functions which will be used are defined on the different subgrids as follows. G i n Ln n j.n \ 0:\a.,b,a,<p. r • *U} on on One further feature of the grid structure needs explanation. G 0 and G ^ contain points which have a radial coordinate less than 0 (the left-most points in Fig. 4.1 ). As discussed in section 2.4, all quantities / defined on these subgrids are assumed to be expandable in power series in r2 about r = 0. In this case, it can easily be demonstrated that by maintaining the identity = / " or f^ + ^ — in the course of a solution, central difference expressions such as AQ/" or A ^ A l / " will retain second order accuracy. In addition, by displacing G 0 and away from r = 0, the immediate need to deal with terms which are singular there (e.g. l/r) is eliminated, although as will be seen, the procedure is definitely not a panacea for numerical difficulties near r = 0. o * o X O X o X O O A j A A A A At A O X O X O X O X O X O A J A A A A A O X O X O X O X O X O A ! A A A A A O * O X O X O X O X O A | A — Ar -— A A A A O x O X O X O X O X O r = 0 F i g u r e 4 . 1 . Staggered m e s h s t r u c t u r e for s o l u t i o n of the sel f -gravitat ing scalar field. 4 . 4 T h e C h o i c e of a T i m e S t e p a n d R e s c a l i n g of t h e L a p s e . In section 3.3 it was shown t h a t the m a x i m u m l o c a l speed of signal p r o p a g a -t i o n for the scalar field is C m . , = - + | 0 | (4.17) o w h i c h w i l l general ly lead to the fol lowing C F L l i m i t a t i o n on the t i m e step of an e x p l i c i t s c h e m e . A* < (4.18) c m o i T h i s result is well k n o w n [28,68,70] a n d t h e analogy w i t h the case of h y d r o d y n a m i c s is f r e q u e n t l y m a d e . F o r a n expl ic i t h y d r o d y n a m i c a l scheme, one m u s t have (4.19) where Vs is the local speed of sound and V is the velocity of the fluid with respect to the coordinate system. The way in which inequality (4.18) is usually maintained in numerical relativity codes is to modify the time step At , so that A* = !-ljL\L (4.20) Here fg is some fixed positive constant less than unity. The problem with this procedure is that, because the time step is changing, simple schemes such as those described in the previous chapter are no longer guaranteed to be second order in time. This problem was avoided here by making use of the fact that both of the slicing conditions, (4.9) and (4.13), are homogeneous in a. Thus, any lapse function which satisfies the maximal slicing, or polar slicing equation on some slice can always be multiplied by some constant factor to give a new lapse function which still maintains the required coordinate condition. In particular, given the grid functions a", a", and 0n = \ir_Qn. i, consider the quantity 3 3 3 J+2 a= min (l - \0n\) (4.21) o<j<J an V 1 J '/ ^ ; 3 If the lapse is rescaled by this factor then the maximum local signal speed will always be unity and a uniform time step can be used. This procedure was used in all of the schemes discussed in this and the next two chapters. (Note that a will become non-positive if | Bj \> 1 for some j, but this situation never occurred in the cases studied). In practice, cmax was invariably attained on the outer edge of the grid, so this process is equivalent to a specification of the outer value of the lapse which is necessary anyway for the slicing equations to be well-posed. The usual procedure is to set a = 1 and argue that the outer edge of the grid is in a "nearly flat" region of the spacetime. If this is true, then a = 1 corresponds to the natural demand that coordinate time be a direct measure of the proper time at infinity. Since o —• 1, 0 —• 0 as r —• co in all of the coordinate systems used here, this feature is recovered using the rescaling process, in the limit r —* oo, provided that the rescaling occurs at the outer edge of the mesh. 4.5 Determining Initial Data. Chapter 2 contained a very brief discussion of the initial value problem for general relativity. There it was noted that, in general, because of the constraint equations, the determination of appropriate initial data is not a trivial problem. In this study, full advantage is taken of the spherical symmetry of the problem, along with additional simplifying assumptions, in order to make the task of choosing initial data as easy as possible [68]. Initial values are required for a, 6, Krr, K9$, $ and II at t = 0. These six functions must satisfy the Hamiltonian constraint (2.46) and the momentum constraint (2.47). In addition, the data must be of a form which is compatible with the coordinate conditions which are to be imposed during the evolution, or some prescription for transforming it to such a form must be given [2,76]. To avoid the extra work involved with the latter possibility, is most convenient to consider determining initial data for the maximal/normal and maximal/isothermal systems separately from the polar/radial system. In fact, the process is sufficiently straightforward that it is only necessary to describe it in detail for the first case. A particularly easy way to generate initial data which satisfies the momentum constraint is to require that the initial data surface represent a moment of time symmetry in the spacetime [54]. This means that the spacetime possesses a t —+ — t symmetry, so that the evolution "backward in time" from t — 0 will be identical to the evolution "forward in time" from the initial surface. Initial data of this type has already been used in the numerical tests of the outgoing radiation boundary conditions for the scalar field alone described in the last section of the previous chapter. A time symmetric slice has an identically vanishing extrinsic curvature tensor; in a d d i t i o n the t i m e derivat ives of the f u n d a m e n t a l m a t t e r variables m u s t also v a n i s h . It is clear f r o m (2.49) t h a t if KTT{r,0) = K99{T,0) = n(r,0) = 0, then the m o m e n t u m c o n s t r a i n t is t r i v i a l l y satisfied. It is also clear t h a t this choice is c o m p a t i b l e w i t h b o t h the m a x i m a l a n d p o l a r slicing c o n d i t i o n s . T h e r e m a i n i n g init ial d a t a is d e t e r m i n e d by freely specifying <j>(r,0) a n d de-m a n d i n g t h a t a ( r , 0 ) = 6 ( r , 0 ) . T h e H a m i l t o n i a n c o n s t r a i n t (2.46) then d e t e r m i n e s the single r e m a i n i n g f u n c t i o n , a. A s in Sec. 3.8, the i n i t i a l configuration of the field is of the f o r m [68] *(r,0) = <f>0exp ^ - j (4-22) w h e r e <po, ro a n d A are parameters. W i t h the c o n d i t i o n s a = b, Krr = i f % — JI — 0, the H a m i l t o n i a n c o n s t r a i n t m a y be w r i t t e n as \ ( r V ) ' + w * V = 0 ( 4- 2 3) r 2tf(r) = o(r,0)* (4.24) T w o b o u n d a r y condit ions are required for the solution of this equation. O n e follows i m m e d i a t e l y f r o m the regulari ty c o n d i t i o n a'(0, r) = 0 V>'(0) = 0 (4.25) T h e o t h e r , w h i c h is i m p o s e d at the outer b o u n d a r y of the n u m e r i c a l g r i d , is an a p p r o x i m a t e r e l a t i o n due to Y o r k a n d P i r a n [107] w h i c h is based on the k n o w n a s y m p t o t i c b e h a v i o u r of ip. A s y m p t o t i c flatness of the spacetime requires that ^ „ 1 + f l + O (4.26) for sufficiently large r. Differentiating the above expression and solving for c\ gives ci = - r V + o Q ) (4-27) Substituting in (4.26) and rearranging yields *'+^! = o(-L) (4.28) which is then replaced by the approximate condition 0 - 1 0 ' + (4.29) r=R where R is the radius of the outer edge of the numerical mesh. An approximate solution to the problem defined by (4.23), (4.25) and (4.29) is readily obtained using finite differences. Working on one row of the grid G 0 described above, the discrete equations are -^A r _ ( r » + , A ^ , - ) + TT *° -. = 0 j — 1, .. . ,J - 1 (4.30) 00 = 01 (4.31) 0 7 1 - 1 *o*j-i + — = <> (4.32) r J - l Equations (4.31) and (4.32) are used to eliminate 0 Q and 0^ from (4.30), resulting in a tridiagonal linear system of J — 1 equations which is solved using a packaged routine [62]. The solution for 0y is then squared to give a°. Due to the use of the staggered grid structure, initial values for the "momen-tum" variables, Krr, K$g, and f l are needed at t — — A t , rather than at t — 0. In o r d e r to preserve second order a c c u r a c y in t i m e , these q u a n t i t i e s s h o u l d be c a l c u -l a t e d to 0( A t 2 ) . T a y l o r e x p a n s i o n a b o u t t = 0, c o m b i n e d w i t h the t i m e s y m m e t r y c o n d i t i o n gives IC.'* = --AtK° + 0{ A t 2 ) (4.33) 3 2 ^ w h e r e K represents a n y of the m o m e n t u m variables. K°_. m a y be c a l c u l a t e d f r o m differenced f o r m s of the right h a n d sides of (2.53), (2.54), a n d (2.68), once all of the d a t a at f = 0 has been d e t e r m i n e d . A c t u a l l y (4.33) has 0 ( A t 3 ) t r u n c a t i o n error since all even t i m e d e r i v a t i v e s of the m o m e n t u m variables v a n i s h as a consequence of t i m e s y m m e t r y . 4.6 The Maximal/Normal Coordinate System. T h e f u n d a m e n t a l n u m e r i c a l variables used in this c o o r d i n a t e s y s t e m are since 0 v a n i s h e s a n d the other extrinsic c u r v a t u r e c o m p o n e n t c a n be evaluated as KTr n. - — —2Kl 2 f r o m the m a x i m a l slicing c o n d i t i o n . B e c a u s e this was the first coordinate system studied, the n u m e r i c a l treatment of the equations is somewhat less u n i f o r m t h a n for the other two systems. In p a r t i c u l a r , following P i r a n [68], the m a x i m a l sl icing a n d c o n s t r a i n t e q u a t i o n s are solved u s i n g a f o u r t h order R u n g e -K u t t e ( R K ) r o u t i n e [8]. A p a r t f r o m the fact t h a t second order accuracy is the best t h a t c o u l d p o s s i b l y be achieved in the overall scheme (so t h a t there is really no p o i n t in u s i n g a f o u r t h order t e c h n i q u e ) , the R K t e c h n i q u e does not generalize to t h e case of a d d i t i o n a l spatial d i m e n s i o n s . T h e success of its a p p l i c a t i o n here relies o n the fact t h a t all three of the e q u a t i o n s c a n be solved as initial value p r o b l e m s . The lapse is determined from the maximal slicing equation by introducing the auxiliary variable 8 = a' and rewriting (4.9) as the first order system ct — 8 6' = fit + ha (4.34) a fr2b2\' . , , 2 , , f l = -TW\-a-) h = a 2 R - M 2 (4.35) The initial conditions are n ^ = ^ A r c ^ c ^ (4.36) where car& is an arbitrary constant, which is typically 1. The expression for 8* follows from a Taylor expansion of a' about r = 0 to 0( Ar 2 ) and the use of the regularity conditions a'(0,t) — $(0,t) — 0. The coefficient functions j\ and h a r e evaluated using centred differences, and the system (4.34-35) is integrated outwards using the R K routine. The lapse is then rescaled as described in section 4.4. In a similar fashion, introducing the variables z = rb and £ = z', the Hamil-tonian constraint (2.46) becomes *' = e i' =-\- + hZ + U*+- (4.37) Z Z Z h = — h = -a a I U = «2 (K\K\ + V / ) - 2TT (n 2 + <_) (4.38) 9 0 w i t h i n i t i a l c o n d i t i o n s *i = fl = ^Ar&J (4.39) H e r e , use is m a d e of the fact t h a t since 8 — 0 a n d K9 ${0,i) —• 0, (K9g(0,t) — Krr(0,t) f r o m regularity, K9g{t,0) = -\Kr r{t,0) f r o m m a x i m a l s l ic ing) , (2.51) gives 6(0, t) — 0 so that to 0{ Ar 2 ) , 6* is constant . O n c e (4.37) has been i n t e g r a t e d , 6^ is c a l c u l a t e d by d i v i d i n g 2" by rj. F i n a l l y , the m o m e n t u m c o n s t r a i n t (2.49) is a single first order equation for K ; K9e = hK99 + f7 (4.40) h = -Z—T h = 4TT 4.41 rb a w h i c h is s o l v e d using the R K r o u t i n e w i t h the i n i t i a l c o n d i t i o n / ^ + 2 " = 0 (4.42) It is k n o w n [68] that the c o m b i n a t i o n of m a x i m a l s l ic ing a n d zero shift v e c t o r i m p l y t h a t the d e t e r m i n a n t of the three-metric is a conserved quantity . T h i s is easily verified i n the c u r r e n t case, using the evolut ion e q u a t i o n s (2.50-51). = [ab2) = 6 2 d + 2a66 \r£ sin 9 J t tl = - a a 6 2 [KrT + 2K99) = -aab2K = 0 (4.43) T h i s fact c o u l d be used to e l i m i n a t e the metric function a f r o m the n u m e r i c a l scheme or at least to p r o v i d e a simple way of u p d a t i n g a j" = a j ° ( 6 ; / 6 ; ) 2 (4.44) F o l l o w i n g P i r a n , r a t h e r than a d o p t i n g this p r o c e d u r e , the e v o l u t i o n e q u a t i o n (2.50) is used to u p d a t e a n d the c o n s t a n c y of the q u a n t i t y a™b2™ serves as an i n d i c a t i o n of the q u a l i t y of the numerical s o l u t i o n . T h e discret ized f o r m of (2.50) is w h i c h is not quite second order in t i m e since the lapse is not evaluated at the center of the c o m p u t a t i o n a l star. F i n a l l y , the equations for the scalar field variables are, f r o m (2.55) a n d (2.58). i n \ la. i \ (4 .52) + 2 a ; * ; ; - * M + n ; - * (4 .53) T h i s s y s t e m is a s t r a i g h t f o r w a r d generalization of the s t a n d a r d scheme for the wave e q u a t i o n discussed in C h a p t e r 3. A g a i n , note t h a t these a p p r o x i m a t i o n s are not n o t c o m p l e t e l y s e c o n d order in t i m e since the a/a t e r m in (4.52) is not properly c e n t r e d in t i m e , nor is the K9$ t e r m in (4.53). T h e second order o u t g o i n g r a d i a t i o n condit ions derived in 3.8 are used to c o m p u t e <f>n,+l, n1 + 2; the other values of.<pn+1 are calculated f r o m the defining J J J r e l a t i o n $ n + ! = Ar,d>n+1. 4.7 Weak field results in the Maximal/Normal System T h e results of a solution of the above equations for a relat ively weak pulse of scalar r a d i a t i o n are displayed in F i g s . 4 .2-8. T h e init ia l c o n f i g u r a t i o n of the scalar field is given by (4.22) w i t h = 0.01, ro = 20.0, a n d A = 5.0, c o r r e s p o n d i n g to an i n i t i a l t o t a l m a s s of 0.066. T h e o u t e r b o u n d a r y is at r = 50.0, 300 grid points are used in the radial direction and At = | Ar . Every sixth and twelfth point in the radial and temporal directions, respectively, is plotted in the figures, and because the plotting routine connects points with straight line segments, some smoothness is "lost" in regions with high gradients. (Some remarks concerning the surface plots displayed in this and following chapters are in order. First, the viewpoint changes from plot to plot to avoid in order to make the most interesting features of each graph visible. The reader should always check the horizontal axes to determine the orientation of the ( r , t) coordinate system. In addition, the self-scaling feature of the graphics routine [20] which was used to produce the surface plots sometimes leads to a mislabelling of the ordinate values on the t or r axes. This occurs because the routine wants the end of the axes to correspond to "round values" of the coordinates. Places where this may lead to confusion will be pointed out.) The low energy of the pulse of radiation is reflected in the fact that the evo-lution of <j> (Fig. 4.2) is virtually indistinguishable from a solution in flat spacetime. The field produces small changes in the geometry of the spacetime, but these changes do not significantly affect the propagation of the field. The small disturbances in the geometry which are "tied" to those of the scalar field are most evident in the extrinsic curvature component Krr (Fig. 4.6), although they can also be seen in the spatial derivatives of the metric functions, a' and b' (Figs. 4.3-4). Because it is the derivatives of the scalar field which couple to the gravitational field, there are two peaks in the gravitational disturbances versus the single peaks of the scalar waves. The behaviour of the lapse function a (Fig. 4.5) during the time interval when the ingoing scalar wave is passing through r = 0 is already indicative of the collapsing phenomenon which slows the evolution in regions of a stronger gravitational field. An interesting feature of the plots of a' and b' is the appearance of the sta-tionary ripples in the region where the scalar waves originate. This is a result of the particular radial coordinate being used, which is completely determined by the 93 initial data condition a(r,0) = b{r, 0) and the choice (3 — 0. As Unruh pointed out, the 3-metric ds2 = a2 dr2 + r2b2 dU2 (4.54) is flat if there exists a radial coordinate R{r), which puts the metric into the form ds2 = dR2 + R2 dU2 (4.55) which requires that {rb)' - a = 0 (4.56) This quantity is displayed in Figure 4.7. It is clear from the plots of this function and the lapse that the part of the spacetime within the sphere encompassed by the reflected outgoing wave is flat—as it should be. The "coordinate waves" in a' and b' could be eliminated by an appropriate radial coordinate transformation on the initial hypersurface. In particular, if a new coordinate, f, defined by f 3 = 3 f a3{f,0)f2df (4.57) Jo were used, the 3-metric on the initial slice would assume the form f 2 ds2 = a2 df2 + — dn2 (4.58) a with \/g = i2 sin 9 (4.59) as in the case of flat spacetime in the usual spherical polar coordinates. In fact, the form (4.59) would be preserved for all time (at least in the analytic case) because of the conservation of yjg (Eq. (4.43)). It is then easy to argue, using (4.56) and the assumption of regularity, that the unique flat-spacetime solution of (4.58) has o = 1. T h i s p r o c e d u r e has not yet been t r i e d n u m e r i c a l l y . In a n y case, the " r i p p l e s " are a g o o d e x a m p l e of features in a n u m e r i c a l l y c o n s t r u c t e d s p a c e t i m e w h i c h develop because of the c o o r d i n a t e choices. S u c h features m a y be difficult to d i s t i n g u i s h f r o m " p h y s i c a l l y - r e a l " effects, or p r o g r a m m i n g errors. T h e b u m p s w h i c h a p p e a r i n c' a n d b' near r = 0 d u r i n g the p e r i o d of self-reflection of the i n g o i n g scalar pulse seem to be due to inadequacies i n the n u m e r i c a l p r o c e d u r e . F i g . 4.8 shows t h a t ab2 r e m a i n s quite c o n s t a n t before a n d after this p e r i o d , b u t t h a t a noticeable p e r t u r b a t i o n is i n t r o d u c e d . T h e m a g n i t u d e of this defect, as well as the size of the b u m p s in a' a n d b' decrease as the m e s h is refined. T h e w a y i n w h i c h the H a m i l t o n i a n c o n s t r a i n t is solved m a y be the source of this effect. In p a r t i c u l a r , because g r a d i e n t s in all functions are greatest i n the v i c i n i t y of r = 0, for a u n i f o r m m e s h s p a c i n g , the a c c u r a c y of any q u a n t i t y c a n be e x p e c t e d to be worst there. S i n c e the c o n s t r a i n t e q u a t i o n is i n t e g r a t e d o u t w a r d s , e r r o r s i n c u r r e d i n the first few values are p r o p a g a t e d t h r o u g h o u t the s o l u t i o n . T h i s results in the oscil lations w h i c h c a n be seen in ab2 i n the p e r i o d 15 < t < 30, w h i c h also a p p e a r i n 6, b u t not in a. It m a y also be t h a t c o n s t r a i n i n g 6" t o have a c o n s t a n t value is not the best t h i n g to do, since the "best" 0 ( A r 2 ) a p p r o x i m a t i o n to the exact s o l u t i o n w o u l d al most certainly have a fluctuating In order to test these hypotheses, several a t t e m p t s were m a d e to solve the H a m i l t o n i a n c o n s t r a i n t as a b o u n d a r y value p r o b l e m using the regulari ty c o n d i t i o n 6' = 0 at the i n n e r b o u n d a r y a n d a value calculated f r o m the e v o l u t i o n equation for b at the o u t e r edge of the g r i d . T h i s is the p r o c e d u r e used i n the s o l u t i o n of the c o n s t r a i n t e q u a t i o n in the i s o t h e r m a l c o o r d i n a t e s y s t e m (Sec. 4.8) a n d it works well there. T h e s i t u a t i o n here is far less satisfactory. A s an i l lustrat ion of the difficulties e n c o u n t e r e d , i n one of the a t t e m p t s , t h e c o n s t r a i n t equation was w r i t t e n as -j{r3q')' + jxq' + f2q+ hq 3 = 0 (4.60) with 2 a' a g = 6 3 /, = /, = _ a 2ri h = 2, (II2 + * 2 ) - a 2 ( K \ K \ + \K\2) - - + (4.61) \ Z / ra LT which was differenced as Ar_ ( r / + l A r + g y ) + / ^ A j g , . + / 2 y 9 y + hflfk = 0 (4.62) y 2 with boundary conditions 9o = ? i «j = *J 1 (4-63) where ft" is to be calculated from a differenced form of the evolution equation for 6. A routine which solved the above system using a global Newton iteration was implemented and a test solution b with corresponding coefficient functions f\, fl and fz were generated. (Analytic expressions for 6, a, II, $ and K% were specified, then (4.60) was solved for Krr.) Plots of the exact solution and two different solutions of (4.62) are shown in Fig. 4.9. The systems were solved on the interval [0,1] with Ar = 1/64. The dotted line is the solution obtained with 3 qj — bj^ — cosh(l) = 1.543... , while the dashed line has qj — {bj + 2 The value e = —3.11 x 10 - 3 was chosen experimentally to minimize the global error. The Newton iteration converges in both cases and the 0(Ar 2) truncation error of (4.62) when applied to b has been verified using a sequence of geometrically decreasing mesh sizes. Furthermore, additional experimentation indicated that, to a good approximation, c £ « c Ar 2 (4.64) 96 Figure 4.3. Weak-f ield results in m a x i m a l / n o r m a l c o o r d i n a t e s y s t e m : a' Figure 4 . 4 . Weak-field results in maximal/normal coordinate system: 6' Figure 4 . 5 . Weak-field results in maximal/normal coordinate system: a 08 Figure 4.7. Weak-field results in maximal/normal coordinate system: (rb)' - a •0 5o. 0 Figure 4.8. Weak-field results in maximal/normal coordinate system: ab ID 0.1 0.2 0.3 0.4 r Figure 4.9. Illustration of sensitivity of the solution of the Hamiltonian constraint to supplied outer boundary value. w h e r e € o p t r - m i n i m i z e s the global e r r o r a n d c is a c o n s t a n t . T h e error in the p e r t u r b e d s o l u t i o n a p p e a r s to converge q u a d r a t i c a l l y , while the u n p e r t u r b e d case exhibits slower convergence at the mesh sizes used (up to A r = 1/1024), a l t h o u g h it m a y b e c o m e a s y m p t o t i c a l l y q u a d r a t i c . A g o o d e x p l a n a t i o n for this b e h a v i o r has not b e e n f o u n d yet, b u t it is a p p a r e n t f r o m the sensit ivity of the s y s t e m (4.62) to the s u p p l i e d outer value, t h a t this is n o t a g o o d way to n u m e r i c a l l y solve the c o n s t r a i n t e q u a t i o n . 4 . 8 The M a x i m a l / I s o t h e r m a l C o o r d i n a t e System. T h e p r i m a r y n u m e r i c a l variables used here are n n+ -T h e c o o r d i n a t e c o n d i t i o n s are used to u p d a t e b. a n d Krr . 2 . It is convenient v 3 3 to c o m p u t e a n d store the grid f u n c t i o n as well as QN L. T h e scheme w h i c h 3 + 2 3+2 is d e s c r i b e d below is, w i t h the e x c e p t i o n of the o u t e r b o u n d a r y c o n d i t i o n s on the scalar field v a r i a b l e s , completely second order in space a n d t i m e . In contrast to the p r e v i o u s scheme, the c o n s t r a i n t a n d c o o r d i n a t e c o n d i t i o n equations are solved using c e n t r e d difference techniques w h i c h w o u l d generalize to the case of a d d i t i o n a l s p a t i a l v a r i a b l e s . T h e lapse is again d e t e r m i n e d f r o m the m a x i m a l slicing e q u a t i o n (4.9) w h i c h is differenced as A ; A r _ a J + fx^a. + f2jan. = 0 (4.66) A ^ a n C . 3 >>, = ~ + T3 " 2 pn (4.67) where a = b has been used. Eq. (4.66) can be solved explicitly for ot"+i, s ° the values of can be determined in order of increasing j (3 < j < J) given initial values % = «! = Cart a!>c a r 6 ( l+^ArV;X) (4-68) where, as previously, cor& is an arbitrary constant and the expression for « 2 is derived using a Taylor series technique. The lapse is then rescaled to maintain the constancy of the maximum local signal speed. The isothermal condition for the shift component is written in the form , 3aKrr 7 = ~ — (4.69) P 1 = -r and differenced as Given the boundary conditions 7 7 = 7 7 » = 0 (4.72) 2 2 these systems are easily integrated outwards. The shift components are then calcu-lated from Ct - (v! - Ci) '4-74' which corresponds to choosing the constant of integration in (4.69) such that the value of 8 at the outer edge of the mesh is 0. This condition is only correct in the limit r, —* oo, whenever there is matter in the numerical domain since, in this for large r. However, the use of an outer condition based on this asymptotic behavior had no discernible effect during the periods when waves were not passing through the outer boundary. When there is a pulse leaving the grid, (4.75) is no longer valid. The Hamiltonian constraint is solved in much the same fashion as described in the section on the determination of initial data, except that now the equation is nonlinear. Again, introducing rp = y/a, the differenced form of (2.46) is coordinate system [76] (4.75) - L A : (( M + r , - ) 8 A ; * ; ) + / , , . + fij4>n; = o 3 (4.76) (4.77) with the boundary conditions (4.78) where c" is determined from a discrete version of (2.50) n - 1 (4.79) which assumes that the (/?a)' term is negligible at the outer boundary. The system (4.76-7 8) is solved using a global Newton method [64] using as an initial estimate values generated from another discretized version of (2.50). Each Newton step involves the solution of a tridiagonal system, and convergence to a given tolerance is (usually) achieved with a number of iterations which is independent of the mesh size. Thus the solution is obtained using O(J) operations. The momentum constraint is the same as in the previous coordinate system. However, the differenced version used here is, = - 3 ,t ihn ..r Ttn+2 x'+*" • « ; n ; The boundary conditions are K};+± = K';+l = 0 (4.82) and (4.80) is easily integrated outwards. The equations for the scalar field variables are differenced using a general-ization of the Crank-Nicholson scheme discussed in the previous chapter. They are where z = ra. As can be seen, rather extensive use of averaging and extrapolation operators is required to center the equations in the interest of maintaining second order accuracy. The inner boundary conditions are $ L = 0 2 and, to facilitate the solution of (4.83) and (4.84) using the tridiagonal solver, outer boundary conditions are also supplied. The outer value for the spatial derivative of the field is computed from AV#;_ 1 = A ; ( x « + ( H ) ; _ i n r _ t ) <4.») n+-which assumes the convective term is negligible at the outer boundary. 11^ 2 is calculated using the first order outgoing radiation condition (3.79). The second order condition (3.82) cannot be used without destroying the tridiagon ality of (4.84). The angled scheme, on the other hand, could be used in conjunction with the second order condition with no loss in efficiency, provided that the unknowns are consistently determined in the direction of increasing j. The value 4r.+i is also computed from a first order outgoing approximation, then the values <f>^+ , j < J are determined from <f>n,+1 and <&n • 4.9 Weak Field Results in the Maximal/Isothermal System. Figs. 4.10-13 show the results of a solution of the above equations using the same initial data and mesh as described in Sec. 4.7. Note the similarity between the plots of a and Krr here with the ones in the maximal/normal coordinate system; the plot of <p, which is not displayed, is indistinguishable from Fig. 4.2. Because the slicing conditions and initial data are identical for the two sets of results presented so far, it is fairly easy to provide evidence that the same physical solution is being generated by the two different schemes. In fact, since the outer value of 3 is always 0 in the isothermal coordinate system, the outer edge of the two numerical grids trace out the same path through the spacetime. This, combined with the fact that a and Krr are both scalars with respect to spatial coordinate transformations, implies that the outer value of o should be the same, for all time, in both systems. Finally, since the rescaling of the lapse for the two runs is always such that a/a is 1 at the outer boundary, the outer values of the lapse functions should also be the same. This is observed numerically to a high degree of precision—the quantity an — 1 is the same in both solutions to 0.2%; a"i itself, is the same to J J within 2 parts per million. It is therefore clear that the sets of t = constant slices generated during the two evolutions are virtually identical. The difference between the two systems is, of course, that different radial coordinates are being used. However, on any given slice in the normal solution, which has ds2 = a2df2 + f 2"b2dU2 (4.87) o.0 F i g u r e 4.10. Weak-field results in m a x i m a l / i s o t h e r m a l c o o r d i n a t e system: a 107 F i g u r e 4.13. Weak-field results in m a x i m a l / i s o t h e r m al coordinate system: 8 108 F i g u r e 4.14. D e m o n s t r a t i o n of weak-field consistency of the two m a x i m a l l y sliced solutions. T h e t o p plot was generated by the m a x i m a l / i s o t h e r m a l p r o g r a m ; the b o t t o m plot f r o m the m a x i m a l / n o r m a l p r o g r a m v i a a n u m e r i c a l c o o r d i n a t e trans-f o r m a t i o n . it is s t r a i g h t f o r w a r d to p e r f o r m the c o o r d i n a t e t r a n s f o r m a t i o n f —> r ( f ) w h i c h puts the 3-metric i n t o i s o t h e r m a l f o r m ds 2 = a2 [dr2 + r2dU2) (4.88) T h e two basic equations are a df = a dr (4.89) f b = rc (4.90) w h i c h give a n equation for r df dr (4.91) w h i c h is n u m e r i c a l l y integrated in the m a x i m a l / n o r m a l p r o g r a m , w i t h the b o u n d -a r y c o n d i t i o n rj = fj. O n c e r has been d e t e r m i n e d , a is calculated f r o m (4.90). A n i n t e r p o l a t i o n o p e r a t i o n is then p e r f o r m e d on a to p r o d u c e a set of values cor-r e s p o n d i n g to a u n i f o r m g r i d in the ( c o m p u t e d ) r c o o r d i n a t e . F i n a l l y this set o f values is differenced to p r o d u c e an a p p r o x i m a t i o n to a'. T h e results of this process are i l lustrated in F i g . 4.14. T h e top plot shows a' as c o m p u t e d by the i s o t h e r m a l p r o g r a m . ( T h i s is a p o r t i o n of the plot of F i g . 4.12, seen f r o m a different v i e w p o i n t . ) T h e b o t t o m plot is a' as c a l c u l a t e d f r o m the m a x i m a l / n o r m a l p r o g r a m . T h e q u a l i -t a t i v e a g r e e m e n t is clearly very g o o d ; q u a n t i t a t i v e l y , the average relative difference between the two functions is a b o u t 3 % . T h i s result provides evidence t h a t , at least for weak pulses, b o t h n u m e r i c a l schemes are p r o d u c i n g the same solutions. 4.10 The Polar/Radial Coordinate System. T h e s i m p l i f i c a t i o n of the 3+1 equations w h i c h results f r o m the use of this c o o r d i n a t e s y s t e m i n a spherical ly s y m m e t r i c spacetime is r e m a r k a b l e . Since the s p a t i a l c o o r d i n a t e condition is b = 1, a n d the polar slicing c o n d i t i o n implies K% = 0, the evolution equation for 6, (2.51), implies that 3 — 0, so these coordinates are also spatially normal. This leaves the set of variables to be dealt with numerically. Not surprisingly, this coordinate system has proven useful in analytic studies of the model—see for example Berger et al [3], and Un-ruh [95]. In a vacuum region of the spacetime, the coordinates are simply those of the familiar Schwarzschild line element « ) d s 2 = _ f\ _ ^ df2 + ^ _ *M_yl d r 2 + r2 d n 2 ( 4 Q 2 ) where M is the mass contained in the interior of the vacuum region. The lapse is determined from the polar slicing equation (4.13), which is dif-ferenced as K°"+fij+i =» (4-93) (4.94) The initial condition is a" = a? = cari (4.95) with c a rj arbitrary. Eq. (4.93) gives a * + l i n terms of o^, the integration proceeds outwards, then the lapse is rescaled as usual. Note that because of this rescal-ing, which typically sets a™ — instead of = l/a" the time coordinate used in the numerical procedure becomes the "true" Schwarzschild coordinate only in the limit r j —> oo. However this simply means that the slices of foliation are rela-belled slightly which does not lead to any undue difficulty in interpreting the results physically. One other point about the solution of (4.93) is worth noting. In their discussion of the polar slicing condition for the case of axisymmetric spacetimes [2], Bardeen and Piran state that the equation for the lapse is a parabolic equation (with the radial coordinate playing the role of time, and x = cos0, the other in-dependent variable) which must be integrated inwards—presumably because the problem would not be well posed for outwards integration. However, in the current case, it can easily be argued that it doesn't matter which way (4.93) is integrated. For example, in the case of outwards integration, (4.93-94) give J - 1 i , n n T T 1 + °~j a , = a, II -J 1 1 1 1 - CT, y=i 3 where 1 A S 3 2 o+h while, for inwards integration, l . n re T T 1 ~ ~ °i a i = a r II j=J-l 3 Now, independent of the integration direction, the lapse is (almost) always rescaled so as to fix the value a^. Clearly, the order in which the factors (1 -I- o~j)/(l — Qj), or their reciprocals, are multiplied together can have no bearing on the "stability" of the difference system (which might be monitored, for example, using the ratio The metric component, o J , is updated at each time step using the Hamiltonian constraint. In the polar/radial coordinate system, (2.46) becomes a first order equation for o a' + ha + ha3 = 0 (4.96) / 2 = _ ( l + 2„r(n> + *')) / ,= ! (4.97) T h e b o u n d a r y c o n d i t i o n for this e q u a t i o n is a(0,t) — 1 w h i c h follows i m m e d i a t e l y f r o m the e l e m e n t a r y flatness c o n d i t i o n a(0, t) = 6(0, t). N u m e r i c a l l y , it is convenient t o i n t r o d u c e the new v a r i a b l e A = In a, a n d r e w r i t e (4.96) as A' + f3e2A + f2=0 (4.98) w i t h the b o u n d a r y c o n d i t i o n A(0, t) — 0. A s e c o n d order difference a p p r o x i m a t i o n is given by A i ^ " + / 3, + i e x p ( 2 l . r + A ; ) + / 2 ) - + i = 0 (4.99) (4.101) w i t h the init ia l c o n d i t i o n A n 0 = A H x = 0 (4.102) (4.99) is integrated o u t w a r d s using a point-wise N e w t o n i t e r a t i o n to treat the n o n l i n -earity. Initial est imates for An. are generated u s i n g a discretized version of the evo-l u t i o n e q u a t i o n for a. It s h o u l d be noted t h a t (4.99) m u s t be integrated o u t w a r d s . A s i m p l e p e r t u r b a t i o n analysis a b o u t the flat solution An = 0, ( /2.J .1 = — / 3 , j _ l ) s h o w s t h a t a s m a l l error i n a value suppl ied at t h e outer edge of the n u m e r i c a l g r i d is a m p l i f i e d d u r i n g a n inwards integration, a n d t h a t the overall a m p l i f i c a t i o n factor diverges as t h e mesh is refined. It is useful to define, in analogy w i t h the line element (4.92), the m a s s aspect function m(r,t) [2,70] m = - r ( l - a ~ 2 ) (4.103) which may be unambiguously interpreted in any vacuum region of the spacetime as the amount of mass/energy contained within the 2-sphere of proper surface area 47rr2 at time t. As a brief aside, using the Hamiltonian constraint (4.96), it is easy to calculate the spatial derivative of the mass aspect rn' = 2TT—* = 47rr2p (4.104) a1 which is clearly non-negative. This, with the condition m(0,t) = 0, provides an easy proof of positivity of energy for this simple model which is reminiscent of Geroch's early positive energy theorem for the case of a maximal slice [29]. The extrinsic curvature component, KRT . 2 , is computed from the momentum constraint (2.49), which is algebraic in this coordinate system K\ = -ATCT (4.105) a This is simply discretized as M r _ ( x + * " + . ) nn+l> Finally, the difference equations for the scalar field variables are A ' + * ; + J = A ; ( , 4 ( f ) ; n; + ' ) A ' n ; - = ^ A i ( , ; ( ^ ) n * ; + , ) (4,08) with boundary conditions = 0 (4.109) n*+2" = n"+2" (4.110) o n the inner b o u n d a r y a n d second order outgoing r a d i a t i o n condit ions on the outer b o u n d a r y . T h e entire scheme has second order t r u n c a t i o n error in b o t h space a n d t i m e . B e c a u s e the c o o r d i n a t e c o n d i t i o n b — 1 is i n c o m p a t i b l e w i t h a = b for any case e x c e p t flat s p a c e t i m e , t h e init ial d a t a generated by the p r o c e d u r e d e s c r i b e d in Sec. 4 . 5 c a n n o t be used d i r e c t l y here. It w o u l d be possible to n u m e r i c a l l y t r a n s f o r m the a — b d a t a to the r e q u i r e d f o r m , but this has not been done. Instead, t h e scalar field is freely specified in the f o r m (4 .22) a n d the i n i t i a l hypersurface is again a s s u m e d to be an i n s t a n t of t i m e - s y m m e t r y . E q u a t i o n s (4 .93 ) a n d (4 .96 ) c a n t h e n o o — — be solved for the init ia l values a . a n d a.. I I . 2 is c a l c u l a t e d f r o m a T a y l o r series _ i_ e x p a n s i o n a b o u t t — 0 a n d the difference e q u a t i o n ( 4 . 1 0 8 ) , a n d K* ^ 2 is c o m p u t e d f r o m ( 4 . 1 0 6 ) . T h e i n i t a l d a t a has, at worst, second order t r u n c a t i o n error. G i v e n t h e s a m e i n i t i a l specification of the scalar field, (i.e. for cf> a given function of r a d i a l c o o r d i n a t e ) the t o t a l mass of the solution tends to be less for the case b — 1 t h a n for a — b, p a r t i c u l a r l y w h e n the mass gets large. 4.11 W e a k f i e l d r e s u l t s i n t h e P o l a r / R a d i a l S y s t e m . T h e results s h o w n in F i g s . 4.15-18 were generated using the same grid struc-t u r e a n d i n i t i a l p a r a m e t e r s for the scalar field as the previous two solutions. A g a i n , a plot of <p is not s h o w n since it looks i d e n t i c a l to F i g . 4 . 2 . A l t h o u g h the m a s s here is s l ightly less t h a n p r e v i o u s l y , to two digits it is still 0 . 0 6 6 . T h e results in t h i s c o o r d i n a t e s y s t e m look p a r t i c u l a r l y c lean n u m e r i c a l l y — t h i s c a n be p a r t l y at-t r i b u t e d to the relat ive s i m p l i c i t y of the equations a n d the fact t h a t b o t h of the " n o n - t r i v i a l " geometric variables are u p d a t e d using the c o n s t r a i n t e q u a t i o n s . N o t e t h a t , i n c o n t r a s t to the case of m a x i m a l slicing, p r o p a g a t i n g d i s t u r b a n c e s are clearly v i s i b l e i n t h e lapse f u n c t i o n , a n d the overall response of a to the presence of m a t t e r is m u c h m o r e l o c a l t h a n for m a x i m a l slicing. A s an i n d i c a t i o n of the precision of the s o l u t i o n , the outer value of the mass aspect function r e m a i n s c o n s t a n t to w i t h i n 0 . 0 2 % in the p e r i o d prior to the a r r i v a l of the initial outgoing pulse. A f t e r the wave has left the g r i d , varies by less t h a n 0.01%. F u r t h e r m o r e , it has been observed t h a t refining the mesh (both spatial ly a n d t e m p o r a l l y ) by a factor of 2 reduces the f l u c t u a t i o n in by very nearly a factor of 4, w h i c h is a g o o d i n d i c a t i o n that the overall s c h e m e has been i m p l e m e n t e d p r o p e r l y a n d is second order accurate. T h e resemblance between the m e t r i c function a s h o w n in F i g . 4.15, a n d the p l o t of the q u a n t i t y used to d e t e r m i n e w h e n the m a x i m a l / n o r m a l solution is s p a -tially flat ( F i g . 4.7), is quite noticeable. A g a i n , this c a n be largely u n d e r s t o o d in t e r m s of the different radial c o o r d i n a t e s w h i c h are used in the two s o l u t i o n s . T h e t r a n s f o r m a t i o n r e l a t i n g the line element in the m a x i m a l / n o r m a l s y s t e m ds2 = - a 2 di2 + a2 df2 + f2b2 dU2 (4.111) to t h a t of the p o l a r / r a d i a l s ys tem ds2 = - a 2 dt2 + a2 dr2 + r2 dll2 is given by - a 2 = f t 2 a 2 - i , f 2 d 2 a 2 = f r 2 a 2 - i / a 2 (4.113) r 2 = f 2 b2 T h e c l a i m here is t h a t the derivatives ft, t r are s m a l l , a n d c a n be neglected in c o m p a r i s o n to f r , ttt w h i c h are of order unity. T h e n a, dr o dx dr w (rb)' df (4.114) (4.112) 116 F i g u r e 4.16. Weak-field results in p o l a r / r a d i a l c o o r d i n a t e s y s t e m : M a s s aspect ( m ) Figure 4.18. Weak-field results in polar/radial coordinate system: KTr where ' denotes differentiation w i t h respect to f. T h u s , a (4.115) B u t since, for the weak field case a — I + a l> a , <C 1 b - 1 + & p bt « : 1 (4.116) it follows t h a t ( l + o i ) 7 « 1 + d i -1 + (f^) = 1 + a - (4.117) so that to a first a p p r o x i m a t i o n a — 1 is j u s t the q u a n t i t y w h i c h is p l o t t e d in F i g . 4.7. T h i s p r o v i d e s m o r e evidence for the m u t u a l c o n s i s t e n c y of the n u m e r i c a l results o b t a i n e d in the three c o o r d i n a t e systems. 4.12 Intermediate Field Results for all Three Systems. In this section, some results f r o m the three different schemes for cases where there is significant g r a v i t a t i o n a l self- interaction of the scalar field are presented. T h e first set of plots, s h o w n in F i g . 4.19-24, were generated in the m a x i m a l / n o r m a l c o o r d i n a t e s y s t e m . T h e init ia l d a t a parameters are <J>Q — 0.053, ro = 20.0 a n d A = 5.0, g iving an initial t o t a l mass of 1.93. T h e outer b o u n d a r y of the grid is a g a i n at r = 50.0, b u t 600 mesh p o i n t s are used in the r a d i a l direct ion a n d the e v o l u t i o n consists of 2000 t i m e steps w i t h At = | A r . In this case, the p r o p a g a t i o n of the scalar field is clearly affected by the c h a n g i n g geometry; this b e h a v i o r will be e x a m i n e d m o r e closely in C h a p t e r 6. T h e variations in the geometric quantities are much greater than previously, and al-though the spacetime appears to be flat to the interior of the final outgoing pulse (Fig 4.23), the problem with the non-conservation of ab2 has become significantly worse. The second set of plots—Figs. 4.25-29—were produced, using the same ini-tial data and run parameters, with the isothermal/maximal program. Note that, qualitatively, the plots of (f>, a, and Krr are very similar to the preceeding graphs. One rather unfortunate feature of isothermal coordinates is that gradients tend to become even more enhanced near r — 0 than they do in normal coordinates, for example, and this single fact explains most of the differences which can be observed between the two sets of plots for a and KTT. There is still fairly good quantitative agreement (a few percent) between the results from the two programs, as deter-mined by comparing the maximum or minimum value of a quantity like KrT on corresponding slices. Notice that the shift (Fig 4.29) attains fairly high values (0 is comparable to a/a), so that the convective terms in the scalar field evolution equa-tions are important. This illustrates the success of the Crank-Nicholson scheme in maintaining a stable and non-dissipative evolution of self-gravitating radiation. It should be noted that the angled scheme has also been used in the program with virtually identical results. The most noticeable numerical difficulty which has ap-peared is the "kink" in <f>, and to a lesser extent in Krr, near r = 0. The precise cause of this behavior is unknown, but the defect does diminish as the mesh is refined. Finding the source of the problem is complicated by the fact that the ir-regularity develops essentially simultaneously in most of the numerical variables (or their derivatives) in the regime where the difference equations are most difficult to analyze. The last set of intermediate field results, from the polar/radial program, is shown in Figs. 4.30-34. The same set of run parameters was used as for the other two sets, but in this case, the initial mass is 1.72. Thus, this is a weaker field than 120 F i g u r e 4.19. Intermediate-f ield results in m a x i m a l / n o r m a l c o o r d i n a t e s y s t e m : <j> 0-0 Figure 4.21. Intermediate-field results in maximal/normal coordinate system: a 122 F i g u r e 4.27. Intermediate-f ield results i n m a x i m a l / i s o t h e r m a l c o o r d i n a t e system: a 125 Figure 4.30. Intermediate-field results in polar/radial coordinate system: <j> 126 Figure 4.31. Intermediate-field results in polar/radial coordinate system: a Figure 4.32. Intermediate-field results in polar/radial coordinate system: Mass aspect (m) 1 2 7 Figure 4.34. Intermediate-field results in polar/radial coordinate system: K rr previously, which accounts for the reduced "overshoot" in <f> at r = 0, as well as the fact that the final outgoing pulse is further out at the time that the evolution is stopped. There are no significant numerical difficulties here—in particular the representation of KrT is actually much smoother that it appears in Fig. 4.34. 4.13 Strong Field Results and Numerical Problems. This last section reports on the behavior of the three codes when the space-times being constructed contain black holes. None of the programs is successful in producing an evolution of a black hole spacetime which can be continued indefi-nitely; in all cases, numerical problems develop which cause the codes to crash, and additional research is needed to remedy the difficulties. Again, the first set of results, shown in Figs 4.35-40, is from the maxi-mal/normal program. The initial data parameters were <f>0 — 0.10, TQ = 15.0 and A — 4.0, giving an initial total mass of 5.36. 600 grid points were used in the radial direction, the outer edge of the grid was at r = 50.0 and 4000 time steps with At = | Ar were taken. The behavior of the lapse (Fig. 4.36) is typical of what has been observed in previous studies of strong gravitational fields [57-61,70,71,76,77,82,87,8 8]. Once the ingoing pulse has reached the central region of the spacetime, the lapse quickly "collapses" and effectively freezes the evolution of the scalar field, as can be seen in Fig. 4.35. The value of a at r = 0 at the end of the run ( t = 167, not 200 as labelled on the plot) is « 7.0 x 10~10. Coincident with the collapse of a, apparent horizons appear in the solution. Recall from Sec. 2.5 that the quantity {rb)'-arbK0e (4.118) vanishes at the location of such a horizon, and that apparent horizons cannot He outside of an event horizon. It is clear then, from Fig. 4.39 that a black hole has formed. However, the problem with the solution of the Hamiltonian constraint for 6, which was described previously, has become much more severe. It appears that the outwards integration of the constraint becomes increasingly less accurate in the vicinity of the apparent horizon where steep gradients form. This leads to the behavior which can be seen in Fig. 4.40 where the value of ab2 is apparently diverging with time outside of the horizon, whereas it should be constant. Again, it would seem that solving the Hamiltonian constraint as a boundary value problem would help but, as discussed in Sec. 4.7, attempts to do this have failed so far. Figs. 4.41-44 show some results in the maximal/isothermal coordinate system using the same initial data and grid structure. In this case, however, the evolution can be advanced for only about 400 time steps before the program crashes due to a divergence of the Newton iteration in the Hamiltonian constraint solver. As can be seen from Fig. 4.44, an apparent horizon has not yet formed at this time. Although quite steep gradients have developed in most quantities, this does not seem to be the sole cause of the problem with the constraint solver since the time at which the code halts is relatively independent of the spatial or temporal mesh sizes used. It is possible that the Jacobian matrix of the nonlinear system for is becoming singular, but this has not been investigated. This problem does not seem to have been encountered by Shapiro and Teukolsky [76,77] in their study of spherically symmetric perfect fluid collapse in the same coordinate system. They also employed a global Newton method to solve the constraint equation, using the type of initial estimate that is used here. They were able to continue the evolution of a black hole essentially indefinitely—the limiting factor in their code is the fact that the entire spatial coordinate system will eventually "fall into the black hole". The final set of plots, Figs. 4.45-49, were generated in the polar/radial co-ordinate system, again with the same initial data parameters and mesh spacings. The total initial mass in this case is 3.62, and the program ran for about 1200 time steps before crashing. It is known [2], that apparent horizons cannot form in this coordinate system in spherically symmetric spacetimes; this is particularly obvious 130 Figure 4 . 3 6 . Strong-field results in maximal/normal coordinate system: a 131 132 Figure 4.40. Strong-field results in maximal/normal coordinate system: ab 133 Figure 4.41. Strong-field results in maximal/isothermal coordinate system: a Figure 4.42. Strong-field results in maximal/isothermal coordinate system: a 134 Figure 4.44. Strong-field results in maximal/isothermal coordinate system: (rb)'—arbK Figure 4.45. Strong - f i e l d dinate system: i> 136 Figure 4.47. Strong-field results in polar/radial coordinate system: a from expression (4.118) which is clearly identically unity for b — 1, Ke$ — 0. Nu-merical experience [70,71,87,88) has shown that the lapse function collapses "even more rapidly" than in the case of maximal slicing, but one can argue that the polar foliation covers all of the spacetime which is outside of the event horizon (as in the case of standard Schwarzschild coordinates). Furthermore, it would seem that a collapsed lapse is a fairly clear indication that the spacetime being generated does contain a black hole. Since the radial coordinate is areal (measures proper surface area), the outer boundary of the grid does not get "dragged inwards" as in the case of the other two coordinate systems, so by waiting for a sufficiently long time, any radiation which will get out to infinity will eventually be seen escaping from the grid. Conversely, if little or no radiation is seen crossing the outer boundary after a sufficient integration time has elapsed, the spacetime almost certainly must contain a black hole. in o If) Figure 4.49. Various functions in polar/radial coordinate system shortly before the strong-field evolution crashes. The collapsing behavior of the lapse is evident in Fig 4.47. The innermost value at the end of the evolution is w 2.0 x 10~3. The evolution of the ingoing pulse has essentially halted (Fig 4.45) by this time. The mass of this pulse is 2.08, with a corresponding Schwarzschild radius of 4.16, and very steep gradients and sharp peaks form in all quantities in the vicinity of this radius (Fig. 4.49). As an example of the type of numerical difficulty which was encountered, consider the differenced polar slicing equation (4.93), solved explicitly for (the temporal superscript is suppressed here): 1 + a , \ ^ h j + L (4.120) For any fixed Ar , it was found that the maximum value of / , i_ (for a point near r = ' + 2 -Mingoing) would eventually become large enough to make a • > 1, at which point further computed values of the lapse would be negative! Needless to say, this quickly stopped the program. Special schemes can be devised to circumvent this problem, but it is clear that when such behavior arises, the basic underlying assumptions of smoothness on the scale of the grid, which are made prior to differencing the equations, are no longer valid. The solution from any difference technique is liable to be suspect in such a case. At the same time, one cannot argue on physical grounds for the introduction of artificial viscosity (as has been used in previous strong field calculations in numerical relativity [57-61,76,77,82]), since the gradients which are developing are due mainly to the coordinate system being used, not the dynamics of the scalar or gravitational fields. Considerations such as these provided much of the impetus for the study of mesh refinement techniques which is described in Chapter 6. Finally, it should be noted that Piran and Stark [71] have also reported numerical problems with this coordinate system, which prevented them from carrying out evolutions of collapsed axisymmetric gravitational radiation indefinitely. C H A P T E R 5 G e o m e t r i c E v o l u t i o n E q u a t i o n s a n d the C o n s i s t e n c y P r o b l e m 5.1 Introduction. The existence of the constraint equations has proven to be a mixed blessing for numerical relativists. On the one hand, experience has shown that it is generally easier to develop a stable scheme if some, or all, of the constraints are used in lieu of evolution equations [24,68]. In addition, comparison of constrained and freely evolved solutions generated from the same initial data has usually shown the former to be numerically "cleaner", and therefore, presumably more accurate. Still, it is possible to construct viable schemes which do not use the constraint equations, other than for the determination of initial data. It is well established that the solutions produced by these schemes do not, in general, satisfy the (discretized) constraint equations on any but the initial slices. This has caused a great deal of concern, both to practitioners of numerical relativity and to outside observers of the field. This concern seems to derive from the general consensus that the only good spacetime is a properly constrained spacetime. In fact, it has been argued [70,89] that there is a fundamental inconsistency in the solutions produced by free evolution schemes which employ standard finite differencing techniques. It is claimed that in such schemes, the discrete constraint equations will be satisfied to a lower order of accuracy than the order of the differenced evolution equations. One of the main purposes of this chapter is to present an argument, supported with numerical evidence, which suggests that the issue of "evolution vs. constraint equations" is less problematic than it has been made out to be. In particular, the results presented strongly suggest that there should be no basic inconsistency of the type outlined above in properly constructed difference schemes for the solution of Einstein's equations. Successful finite differencing of the geometric evolution equations was, of course, a necessary prerequisite for generating these results, as well as one of the principal goals of the thesis, so the schemes used are discussed in some detail. 5.2 The Argument for Inconsistency in Numerical Relativity. The basic argument which has been presented for the inconsistency of free evolution finite difference schemes is quite simple [70,89]. Let Q and T denote the sets of geometric and matter variables for some solution of the field equations. As-sume, for simplicity, that the problem has a dependence on only a single spatial coordinate x and time. Then, for any particular metric or extrinsic curvature com-ponent g E Q, let and g^ denote the exact and approximate values, respectively, on some grid with fixed mesh spacings In addition, assume that all differential expressions are replaced with second order finite differences. Then a typical evolution equation for a geometric quantity At = cQh, Ax = cxh (5.1) c 0, C l = O(l) 9 = HQ* T) (5.2) becomes, upon discretization Bf»;-L*($;,f;) (».>) where D*f and L f c are 0(h?) approximations to d/dt and L. Now, from the definition •of the order of a difference scheme, = L f c ( 5 ; , f ; ) + o ( f c » ) ( 5 . 4 ) where the 0(fe2) quantity, f*, is just the local truncation error. Provided that the scheme is stable, and all functions involved are sufficiently smooth, it is reasonable to expect * ; = y ; + e ; = , ; + o ( f c 8 ) ( 5 . 5 ) where t \ is the error in the computed solution. (In principle one should be able to establish the validity of ( 5 . 5 ) from the difference and differential equations and some additional assumptions (such as smoothness), but for a typical scheme for Einstein's equations this is a very messy task, and it is never done in practice.) Now, the argument continues, consider one of the constraint equations ff(S,T)=0 ( 5 . 6 ) and assume, as will often be the case, that this equation includes terms involving the first and second spatial derivatives of g. Thus, for the sake of argument, write ( 5 . 6 ) as Hg = P (9, 7) g" + Q (9, T) g' + R (Q, T) = 0 ( 5 . 7 ) and difference it to second order as H k * ; = P * ( ^ f ; ) A * A ' y » + Q* (9n3,fn.) A*gn. + /2 f e (^ n , f ; ) =0 (5.8) where is another approximation of g™, which in general will not coincide with g". Now, if the substitutions gn. —• gn, y. —+ Gn. and T" —* T™ are made in the last expression, the right hand side becomes another 0(h2) quantity which is just the local truncation error, f ?, of the discretized constraint equation. Then replacing g^ by g* in (5.8) and using (5.5), one gets H*f. = f; + p h (s*, f;) A * A i e ; + Q * (<?;, t ; ) A 0 * e ; + o ( * » ) (5.9) where there are additional 0(ri2) terms which may arise because the coefficient functions Ph, Qh, and Rh in (5.8) can depend on g .. Now, the crucial claim in the inconsistency argument is that since e? is an 0(h2) quantity, and the operators A+Al and A p involve the factors h~2 and h~l respectively, it follows that A * A l e ; = 0(l) (5.10) A*0en.=O(h) (5.11) Thus, it is concluded, the discrete constraint equation, when applied to the freely evolved second order approximation is satisfied to zeroth (or, at best, first) order and consistency is lost. 5.3 The Argument for Consistency in Numerical Relativity. The proposed counterargument to the above analysis is equally simple. The key point goes back to Richardson's [73] 1909 hypothesis that "the errors of the in-tegral [solution of the difference equations] and of any differential equations derived from it, due to using ... central differences instead of differential coefficients, are of the form" h2f2{x,t) + h*f4{x,t) + ---It is important to stress that the functions f2, ft, . •. are independent of the mesh spacing. Richardson illustrated the validity of this hypothesis by explicit numerical calculation of a problem with an analytic solution and then suggested the technique, now known as Richardson extrapolation, whereby solutions of the same difference equations using different mesh sizes are combined to eliminate terms in the error expansion to produce a more accurate solution. The validity of the Richardson hypothesis for a given differential system and difference scheme may be established by substituting expansions of the form S} = <*3 ' e 2 . + h* e 4 . +•••) (5.12) for all relevant functions into the difference scheme. The difference operators are then replaced by their corresponding differential operators and truncation expan-sions. The zeroth order terms in the resulting expression will be eliminated by virtue of the original differential equations (consistency of the difference scheme) and additional systems of differential equations for the error functions e2, t\, ... will result. It must then be argued, order by order, that these systems are well posed so that the error functions are bounded. (See Chapter 10 of Kreiss and Oliger [44] for a simple example and Marchuk [50] for a complete discussion of the procedure). This process is not trivial for non-linear systems such as Einstein equa-tions, but accepting the Richardson proposal, is probably not much more an act of faith than accepting the hypothesis that the field equations themselves have a unique, well-behaved solution. Once the Richardson hypothesis has been accepted, it is easy to see why there should be no consistency problem. For example, the term A+Ale^ becomes Az+Ax_en. = A* A l [h2e2n.+h4 etn. + • • • ) = h2 ( A ; A l e 2 ; + /i 2 A*A x _e 4 n . +•••) = h2(e2"" + 0{h2)) (5.13) Similarly Ax0en.^h2{e2'n.+O(h2)) (5.14) So provided e2 and e2" are well behaved, these results show that all terms on the right hand side of (5.9) are 0(h2), in contradiction to previous claims. This type of argument also holds for the case of discrete evolution equations being applied to constrained quantities which, has also been purported to lead to the same kind of inconsistency. Both of these claims (and others which appear in the literature) are based on the incorrect notion that a numerical differentiation must automatically reduce the accuracy of any quantity by one order. 5.4 The Geometric Evolution Equations in the Polar/Radial System. The equations to be solved here are the same as in section 4.10 with the addition of the evolution equations for a and Krr. Using the coordinate conditions 6 = 1, /? = 0 and K = K'\ in (2.50) and (2.52) these are a=-aaKTr (5.15) K\ = 6+ aK\2 (5.16) where . „ _ i ( = : j ' . 0 ( i ( i y + ^ ) a \ a / \ar\a/ c r / (5.17) In analogy with the treatment of the scalar field, instead of evolving a directly, the quantity It should be pointed out that there is no compelling reason to evolve A or a' rather than a. In fact, for the case of second order centred schemes, to which the current study is limited, virtually the same difference equations can be developed using either form. This will not necessarily be the case for higher order schemes, however, where it may be advantageous to keep the pair of equations for the metric and extrinsic curvature components as symmetric as possible (with respect to the spatial derivatives) so that the same.methods can be used to difference both equations. In addition there is some case to be made for evolving a quantity such as A whose value fluctuates about 0, versus c which has a nominal value of 1, from the point of view of numerical precision. On the other hand, the increased complexity which results from differentiating a general evolution equation may limit the usefulness of the technique. As it is written, the evolution equation for Krr does not explicitly contain the second spatial derivative of o (first derivative of A). Such terms are present, how-ever, as can be seen by using the slicing condition (4.13) to replace the occurrences A = (In a)' (5.18) is introduced. The evolution equation for A is then A — - (aKrr)' (5.19) of a' and a" in (5.17). After some manipulation, (5.17) becomes (5.20) where z = a 2 - 1 (5.21) With s rewritten in this form, it is tempting to adopt the simple pair of equations A = -aKrr' (5.22) K\ = -aA' (5.23) as a model system to facilitate stability analyses of possible difference schemes for (5.15-16). This approach is of questionable merit, however. For example, if a is treated as a constant, then a stability analysis of an explicit scheme for (5.22-23) would surely result in a CFL condition aA« < Ar (5.24) instead of - A t < Ar (5.25) a which is the expected condition, since as argued in Chapter 2, (and as seen in the results of the previous chapter) the disturbances in the gravitational field must propagate with the same speed as those in the scalar field. A second, related problem is that such an analysis ignores the influence of the scalar field. Yet without the scalar field, there can be no dynamics in the geometric variables, so it would seem that any proper treatment of stability would have to be done in conjunction with the equations for the matter variables. The main difficulty is due to the fact that, unlike the case in any other dynamical system, the geometric variables determine the causal structure of the "background" on which the disturbances propagate. As a result, the form of the evolution equations strongly depends on the basic variables which are used and one must be very cautious about trying to relate a given system to some model system. The approach which is taken here is to assume that the same basic schemes which work for the scalar field should also work for the gravitational field since the structure of the two sets of equations is similar. The evolution equation for A is not difficult to treat numerically. The differ-ence scheme used is the obvious one Au ; r - A;(xW^; + l) (5-26) The inner boundary condition is AnL = 0 (5.27) which, because of the grid structure, does not have to be explicitly imposed after the first time step. Likewise, no outer boundary condition is required. Following the computation of of ^ 4 " + J , the metric function an+1 is updated from n+l I A An+1 I n+1 . „ . °, = e x P I 2 ^ A r A k + L ) ' % = 1 (5 2 8 ) .fc=o 2 The evolution equation for Krr, which proved more difficult to handle, is differenced as This is a set of decoupled nonlinear equations for the quantities K^+2 which may be solved using a (single variable) Newton method with initial values supplied from the 0( At) difference equation < i f ; ; ^ 8 ; + tt;(^;-')2 (5.so) The inner boundary conditions are JT;; + » = jr; ; +» =0 (5.31) At the outer edge of the grid a first order outgoing radiation condition is used which is based on the large r behaviour K r f i t K r ^ + o ( l ) (5.32) This is most easily deduced from the constraint equation (4.105), together with the known behaviour of II and a, and the assumption that the spacetime is nearly flat where the condition is imposed. In the first attempt at evolving the extrinsic curvature with this scheme, the source term s was differenced from (5.17) as 1 / A r _ a n \ - °* (^A+ (•/ "-<)+8* ('+*W'lY) (5-33) This procedure led to an instability which can be seen in Fig. 5.1 (In this, and following plots, the superscripts E and C are used to identify quantities which are updated using evolution and constraint equations, respectively). The stationary disturbance grows fairly slowly at first, then rapidly explodes, halting the solution. 149 Figure 5.1. Instability in polar/radial free evolution when the differ-ence expression (5.33) is used in the scheme. The same behaviour appears in the extrinsic curvature component. This does not seem to be a "classic" instability in the sense that some Fourier components are being magnified since a reduction in the time step does not appreciably affect the absolute (as a function of t) rate of growth. At the same time, decreasing the spatial mesh spacing does slow the growth rate. This suggests that expression (5.23), combined with the slicing equation for the lapse, (4.93), may in effect be a "better" approximation to something like s + Arp, where p involves higher order derivatives of o, rather than s itself. It is plausible that the differential system, modified in this way, would admit the growing solution observed in Fig. 5.1 [100]. Although the cause of the instability is almost certainly due to inappropriate differencing of s, the problem can be fixed by using the constraint equation (4.105) to eliminate Krr from the evolution equation for A. The solution then remains stable until the ingoing pulse reaches r = 0, at which point irregular behaviour in KTr quickly develops. Apart from this last problem, this technique does not produce a "true" free evolution scheme, since the direct coupling of Krr to A, through the evolution equation, has been removed. It was at this point that the form (5.20) for s was derived, the idea being that it would be best to start from a differential expression in which the dependence of Krr on a and its derivatives is as explicit as possible. The differenced form of (5.20) used is -«;(»;(«;-i)/'*+»"('i#;-t/•;)')' (5-3<) The difference equations for A and KTT were incorporated into the scheme described in Sec. 4.10. In effect, the program was constructed so that completely free and completely constrained evolutions of the geometric variables, using the same initial data, were performed simultaneously. In general, the constrained metric functions were used in the equations for a, $ and II. However, the freely evolved quantities could also be used without affecting the overall stability or convergence of the program. 5.4.1 Comparison of the Freely Evolved and Constrained Quantities. Typical results of this scheme, for the same initial data and grid structure used in Sec. 4.10, are shown in Figs. 5.2-5.7. The agreement between constrained and evolved quantities is very good—the differences which can be observed in the last two plots are most likely due to dispersion of the evolved functions. The scheme remains stable up to the CFL limit, the solution is well behaved near r = 0, and the radiation conditions allow the waves in the geometric variables to propagate cleanly off the mesh. These comments continue to hold as the mass of the solution is increased—the cause of the ultimate demise of the code remains the inability to handle the very large gradients which develop near r = 2Mingoing for those spacetimes which contain black holes. During the course of the evolution, the deviations of the freely evolved geomet-ric quantities from the discretized constraints and the deviations of the constrained functions from the discrete evolution equations are also computed. Suppressing grid indices, these defects are denoted as 0 a , 0 a , 0 K\" and 0 KTr , where C E to avoid overcomplicating the notation, the definition of the operators 0 and 0 depends on what they are operating on according to 0Ca° = 0 o E a E = o O C K ; C = O 0 E K R R E = 0 (5.35) Because the momentum constraint is algebraic in this coordinate system the quan-ta E tity 0 Krr can be evaluated simply as O C K ; E = K ; C - K [ E (5.36) C E C E Quantities analogous to 0 a and 0 KTr have been computed and displayed previously by Piran in his study of cylindrical gravitational waves [68]. In addition, Nakamura and his coworkers, [57-61], in studying various aspects of axisymmetric stellar collapse, have adopted the general strategy of freely evolving all of the geo-metric variables and using the measured deviations from the constraint equations as an indication of the reliability of the solution. Centrella [10-12] used a similar E C technique in her study of plane-symmetric cosmologies. Quantities such as 0 a , which measure how well the constrained geometric variables satisfy discrete versions of the evolution equations, do not appear to have been computed before. 152 Figure 5.3. Polar/radial coordinate system: Metric function a updated using evolution equation. Figure 5.5. Polar/radial coordinate system: Extrinsic curvature component K r updated using evolution equation. 154 m Figure 5 . 6 . Late-time comparison of freely evolved and constrained a! in po-lar/radial system. Figure 5 . 7 . Late-time comparison of freely evolved and constrained KTr in po-lar/radial system. Figure 5.9. Deviation of constrained a from evolution equation in polar/radial system: interior region. 156 Figure 5.11. Deviation of constrained KTr from evolution equation in polar/radial system: interior region. 157 Figure 5.13. Deviation of freely evolved K rr from momentum constraint in po-lar/radial system: exterior region. Plots of the four deviations defined above are displayed in Figs. 5.8-13. These were calculated in the course of the same run which produced the data of the previous figures. Because these functions are more oscillatory than a or Krr, the loss in resolution in plotting the entire domain is unacceptable, so the first four graphs cover the inner third of the grid for the first half of the evolution. All of the functions display the same qualitative behaviour. The deviations track the scalar pulse inward, are magnified near r = 0, and then follow the reflected pulse outward. There is a noticeable overall amplification in all quantities as a result of the interaction near the origin, but other than that, the deviations do not grow a great deal during the evolution. The initial growth rate is certainly no more rapid than that of the momentum variables. The large ripples which appear in 0 a during the period of self-reflection, but not in the other graphs may be understood c in terms of the "non-locality" of the difference system for a . A perturbation introduced near r = 0 is propagated outwards by the solution process. Since this perturbation is largest and most rapidly varying when the pulse is near the origin the observed behaviour is not very surprising. The other two plots show the outer third of the grid at later times. The largest disturbances are due to the originally ingoing pulse. The other features are due to the outgoing part of the initial pulse—the radiation condition imposed on KTr is not perfectly absorbing and an inwardly propagating disturbance forms. These results are quite different from those displayed by Piran in [68]. His plots of the deviations show much more erratic behaviour and it is difficult to detect the propagation of the wave packets in them. In addition his defect quantities did not return to a quiescent state after the waves had left the grid, as the functions plotted above do (except for pulses reflected from the outer boundary). It also seems that these differences cannot be attributed solely to the use of different mesh spacings, or the fact that the spacetime generated here is nearly flat, since the qualitative features of Figs. 5.8-13 are visible in both strong field and coarse mesh cases. 5.5 Demonstration of Consistency in the Polar/Radial System. In order to demonstrate the evolution/constraint consistency of the difference scheme, a series of runs which evolved the same initial data on grids with varying mesh spacings were performed. The initial data parameters were <f>o = 1.0 x 10~3, r0 = 50.0 and A = 5.0 giving a mass of 3.9 x 10 - 3 and the outer edge of all grids was at r = 100. The initial data and number of time steps were chosen so that there was no interaction with either boundary before the runs terminated. This was done so that, as much as possible, the numerical experiments would be unaffected by boundary effects, since the arguments for and against inconsistency are based solely on the interior difference equations. Grids with Ar = 2 At = ^ 0 < / < 4 (5.37) 2 were used, and all runs stopped at t = 25. Thus the coarsest grid was 100 x 50 / C E C E E C and the finest 1600 x 800. Every 2 steps the l2 norms of 0 a , 0 Krr , 0 a , E C C E 0 KTr > and a — a were computed and output. The results of the various runs were then combined to calculate convergence factors for the deviations. The level I convergence factor c L of any of the above quantities, Q, is the ratio (5.38) which can be computed for 1 < n < 4. The results of this process are shown in Figs. 5.14-18 (Note that in these plots, h is the finest radial spacing used and K = K r r ) . The dotted line in each graph indicates the ideal convergence factor of 4. Although the convergence rates 160 O iri ui O * o Q> O q <u * o c P S CD > o q 0.0 Vri*! •*! ft ."ui, _ iii. _ _ ^ ^ r A , ^ , -a 16h : 8h -A 8h : 4h H H 4h : 2h -*—x 2h : h 5.0 10.0 15.0 20.0 25.0 t Figure 5.14. Convergence of deviation of freely evolved a from Hamiltonian con-straint in polar/radial system. Figure 5.15. Convergence of deviation of constrained a from evolution equation in polar/radial system. o i n - ^ _ _ X T O ^ q o c > C O 9 O ^ m oi 0.0 c V YY YT YYT YT Y T Y Y ¥ ¥¥' * *'Y * Y'* * Y'Y w ¥ ¥'¥ WW* ¥ ¥ * « A A A A A A A A A A A A A A ,», » > » » » I'I I'I A u ^ ji, ^ ^ ^ 5.0 10.0 15.0 20.0 25.0 t Figure 5.16. Convergence of deviation of freely evolved KTT from momentum constraint in polar/radial system. o iri in j U & A * * » x V.'iT ¥¥'¥¥'¥ YY ¥ Y ¥ ¥ Y ¥ ¥¥ W f f f l C3—B—• 16h : 8h & — A — A 8h : 4h H 1 H 4h : 2h x — x — x 2h : h 0.0 5.0 25.0 Figure 5.17. Convergence of deviation of constrain ed KTT from evolution equation in polar/radial system. 162 O i r i o m 1 "t o O o o CD O C in CD O) i_ > C o o o CN 0.0 • — B — • 16h : 8h - — A — A 8h : 4h H 1 H 4h : 2h x—x—x 2h : h 5.0 10.0 15.0 20.0 25.0 t Figure 5.18. Convergence of difference between constrained and freely evolved a in polar/radial system. deteriorate with time on the coarser grids, there can be little doubt that all of the deviations are 0(h2) quantities, so there is no inconsistency. 5.6 The Geometric Evolution Equations in the Isothermal/Maximal System. The scheme described in this section involves the equations of Sec. 4.8 and the additional evolution equations (2.50-52) and (2.54). Although two of these equations are unnecessary because of the coordinate conditions a — b and KTr = —2K9$, it is instructive to examine all of the equations before deciding on which pair to evolve. As in the previous section, auxiliary variables A = (In a)' B = (In 6)' (5.39) (5.40) are introduced. Then (2.50-51) lead to A= {PA)'- (aKrr)' + 0" (5.41) B = (/?_)'- (aK',)'+ (^j (5.42) and (2.52) and (2.54) become, with K = 0 Now, as can be seen from the isothermal condition for the shift vector, (4.69), the term 0" on the right hand side of (5.41) involves spatial derivatives of the extrinsic curvature components. Likewise, from (4.9), the term involving a" in the evolution equation for Krr involves both first and second derivatives of the metric functions. If the right hand side of (4.69) is differentiated and used to eliminate 0" from (5.41), then (5.42) is recovered, since one is just undoing the procedure that was used to derive the condition for the shift in the first place. For the same reason, if the maximal slicing equation is used to eliminate a" from (5.43) and the substitution KTT —• —2K\ is made, the result is simply (5.44). In view of the results of Sec. 5.4 this suggests that it would be best to difference the pair (5.42) and (5.44) and numerical experiments confirm this hypothesis. Schemes constructed by applying the usual centred differencing techniques to (5.41) and/or (5.43) demonstrated non-propagating, growing solutions very similar to the one shown in Fig. 5.1. This type of behaviour also seems to have been encountered by Nakamura et al [59], in a study of spherically symmetric stellar collapse. There, radial coordinates were used in conjunction with maximal slicing and they encountered problems with an evolution equation which contained a term proportional to 0' where 0, in turn, was proportional to an extrinsic curvature component. They went further than using the coordinate condition to eliminate 3', however. The momentum constraint was also used to eliminate the resulting spatial derivative of the extrinsic curvature component. This stabilized their scheme, but it was not made clear that the last step was strictly necessary. In any case, it would appear that it is advisable to use the coordinate conditions to ensure, as much as possible, that the spatial derivatives of the metric and extrinsic curvature components appear explicitly in the evolution equations prior to differencing. The other major problem which was encountered in the construction of a free evolution scheme in this coordinate system was irregular behaviour in the geometric variables in the first few grid points near r = 0. An example of this type of behaviour has already been seen in the previous chapter in Fig. 4.25. The situation with the evolution equations for A (substitute B — A in (5.42)) and K9 g is much worse however, since irregularities are generated even for weak pulses. Part of the problem seems to be due to the presence of terms in (5.44) which individually diverge but collectively cancel as r —• 0. (The right hand side of (5.44) must vanish at r = 0 to be compatible with K9f(0,t) — 0.) Originally, (5.44) was differenced essentially "as is"; following discussions with Ascher, O Murchadha and Unruh, the equation was rewritten so that the terms proportional to r _ 1 were grouped together before differencing. This was not sufficient to regularize the behaviour. Numerical experiments suggested that there was no problem in differencing the expression to 0( Ar 2) everywhere, but that near the origin, the spatial derivatives of K9g would be poorly approximated. This type of problem has occurred before in numerical relativity codes [26,104]—the most successful approach to the problem so far seems to be the procedure suggested by Nakamura [57] of using "regularized variables" such as r~2(K rr — K9 g) as the basic quantities for a numerical solution. However, this method also involves the use of x = r2 as the spatial coordinate, which does not seem appropriate for radiation problems, since, for example, an outgoing wave packet will be represented on fewer and fewer mesh points. The best results obtained in this research involved the introduction of new variables A = r2A (5.45) K\ = r2K\ (5.46) as the basic quantities to be evolved. The decision to scale A and K% by r2 was motivated by a somewhat vague notion that this would inhibit the irregularity. The evolution equations for these variables are derived from (5.41-42) A = 3A - {aK09)' + (^B' -^JA + yK99 + rfi'- 3 (5.47) -i ' 23 ~ „ rlaa)' + (r2aa')' K e = BK\ - -?-K90 - 1 7 ; f- (5.48) where the coordinate condition a — b has been used to obtain the last expression. Angled differences (Sec. 3.3) are used to properly centre the convective terms. (A variation of the Crank-Nicholson scheme can also be used.) The difference equations are -3 f A A. i + A.A. i - A l y , a . Kl 2 3+2 — V+X-e j J+ 2 ..t „ 0 n + 2 - (r, A; (ao); + A!, ( ML (ra) J A l a ' ) ) / a 3 * (5.50) With the inner boundary conditions AL = 0 (5.51) 2 k°nQ+1* = K 9 ^ = 0 (5.52) Both of these systems can be solved explicitly using an outwards sweeping direction. Outer boundary conditions must also be supplied since the expressions , J + 2 and Ar+Kg™ 2 are not defined at the last grid point. Again, first order outgoing radiation conditions based on the leading order asymptotic behaviour A~f{t-r) (5.53) \{^A ~9{t-r) (5.54) are used. After has been calculated, the metric function an+1 is updated using J+ 2 3 n+l — exp / j-i n + l no, - <?{';U*'*Z\j) <»•«> where the outer boundary value is determined from equation (4.79). K$*+2 is computed from *;; +* = r,»*;;+* (5.56) 5.6.1 Comparison of the Freely Evolved and Constrained Quantities. As with the scheme described in the Sec. 5.4, the above difference systems were incorporated into the program described in Sec 4.8 so that a direct comparison between the freely evolved and constrained quantities could be made. Figs. 5.19-26 Figure 5.20. M a x i m a l / i s o t h e r m a l c o o r d i n a t e system: (In a ) ' e v o l u t i o n e q u a t i o n . = A u p d a t e d using 168 leg o ^ I I I I L_ 10 20 30 40 r Figure 5.23. L a t e - t i m e c o m p a r i s o n of freely evolved a n d constrained (In a ) ' in m a x i m a l / i s o t h e r m a l system. Figure 5.26. Maximal/isothermal coordinate system: K $. This is the other function which is evolved. show some results from a run made with the same initial data parameters and grid structure used earlier. Again, the agreement between the constrained and evolved function appears very good. Figs. 5.25-26 show the quantities which were actually evolved. Note that A looks very much like the negative of the mass aspect function shown in Fig. 4.16. This is to be expected, since in isothermal coordinates , M , * ~ - - Y (5.57) for small M , in a vacuum region. This behaviour may also be verified by performing the coordinate transformation from radial to isothermal coordinates under the same assumptions used in Sec. 4.11. Figures 5.27-28 show the effect of the introduction of the variables A and K9g. The top plot shows the quantity A which was evolved directly in conjunction with K9 , using a centred scheme very similar to (5.49-50). The bottom graph is the same function calculated from r~2A. The kink in the second radial value of A in Fig. 5.27 is typical of the irregular behaviour which was observed. Although the use of the variables A and K6\ appears to have improved the situation somewhat, examination of the deviations between constrained and evolved quantities, shown in Figs. 5.29-30 suggest that there is still a problem. The functions plotted in these figures are defined by A (In a)' = (in a°)' - r2A (5.58) AK$ =K99C - r2k99 (5.59) Note that the deviations in K9g are largest near r = 0 and a small cusp appears after the pulse has left. In addition, A (In a)' is not smooth near the origin when the quantity is less than 0. As the mass of the spacetime increases for a fixed mesh spacing, the cusping gets worse although it always remains localized near the axis and does not severely 172 173 Figure 5.30. Maximal/isothermal coordinate system: Difference between con-strained and freely evolved extrinsic curvature components near r = 0. 174 Figure 5 . 3 2 . Deviation of constrained a from evolution equation in maximal/isothermal system: exterior region. 175 imal/isothermal system: exterior region. mal/isothermal system: exterior region. 176 Figure 5.36. Deviation of freely evolved K% from momentum constraint in max-imal/isothermal system: interior region. affect the exterior s i tuation. T h e precise cause of the effect has yet to be d e t e r m i n e d , a n d this r e m a i n s a n o u t s t a n d i n g p r o b l e m w i t h this scheme. D e v i a t i o n s of the evolved q u a n t i t i e s f r o m the c o n s t r a i n t e q u a t i o n s a n d vice v e r s a are s h o w n in F i g s . 5.31-36. T h e first four plots exclude the i n n e r m o s t p o r t i o n o f the solution d o m a i n since the d e v i a t i o n s near r = 0 are m u c h larger. T h e same q u a l i t a t i v e features a p p e a r here as i n the previous s c h e m e except t h a t the defects w h i c h track the reflected outgoing pulse fall off w i t h i n c r e a s i n g distance. Ingoing pulses in the d e v i a t i o n are again generated at the outer b o u n d a r y p r i m a r i l y because the r a d i a t i o n c o n d i t i o n s on the geometric variables are not perfectly a b s o r b i n g . T h e (J E C ft E b e h a v i o u r of 0 a a n d 0 Kg near r = 0 is s h o w n in F i g s . 5.35-36. T h e r a p i d g r o w t h a n d lack of s m o o t h n e s s in t h e first few grid p o i n t s are clear indicators t h a t i m p r o v e m e n t is needed in this area. 5.7 Demonstration of Consistency in the Maximal/Isothermal System. T h e consistency of this scheme was tested i n the same fashion as d e s c r i b e d p r e v i o u s l y . In this case, the initial d a t a p a r a m e t e r s were the same as those w h i c h g e n e r a t e d the results d i s p l a y e d in the previous sect ion, a n d the o u t e r edge of all grids was at r — 50. F o u r grids were used with Ar = 4 A* = ^ 0 < I < 3 (5.60) a n d the r u n s e n d e d at t w 6. T h e results of the test are shown in F i g s . 5.37-42. A g a i n it seems clear t h a t there is no basic inconsistency between the c o n s t r a i n e d a n d E ff C a C ff E e v o l v e d quanti t ies . T h e only p r o b l e m s seem to be w i t h 0 Kg a n d Kg -Kg . T h e convergence factor of the former a p p e a r s to t e n d to a value less t h a n 4, while t h a t of the latter begins to deteriorate at the highest level of refinement. T h i s suggests t h a t there m a y be some p r o b l e m w i t h the way the m o m e n t u m c o n s t r a i n t is s o l v e d . It is also possible t h a t the rather extensive use of e x t r a p o l a t i o n in this s c h e m e c o n t r i b u t e s to this effect a n d it m i g h t be i n s t r u c t i v e to m o d i f y the scheme O N o Q> 0 9 <D * o c CD > C O O O in in 0.0 • — B — • 8h : 4h £ s — & — A 4h 2h - I 1 1- 2h h °" -^ IHH'A' * A'* ^ AAAA'AAA AA&£** 1.0 2.0 3.0 t 4.0 5.0 6.0 F i g u r e 5.37. C o n v e r g e n c e of deviation of freely evolved a f r o m H a m i l t o n i a n c o n -s t r a i n t in m a x i m a l / i s o t h e r m a l s y s t e m . o in m CN 0.0 -——— 4h : 2h — I h 2h : h 1.0 2.0 f $ a A l A A A A A A A A A A A ji AA A A S a n a • • • • • • • • • • • • • 3.0 t 4.0 5.0 6.0 F i g u r e 5.38. C o n v e r g e n c e of deviat ion of c o n s t r a i n e d a from evolution equation i n m a x i m a l / i s o t h e r m a l s y s t e m . 0.0 ^ • • • • • • • • • • • • • • • • • • • D O • — B — E 3 8 h : 4h A A A 4h 2h H 1 1- 2h h 1.0 2.0 3.0 t 4.0 5.0 6.0 F i g u r e 5 . 3 9 . C o n v e r g e n c e of deviation of freely evolved K9$ f r o m m o m e n t u m c o n s t r a i n t i n m a x i m a l / i s o t h e r m al system. q i n U l *5 9 o c CD > c o q 0.0 • — B — B 8 h : 4h A — A — A 4h : 2h H 1 1- 2h : h 1.0 2.0 S B B D • • • • • • • • • • • • • 3.0 t 4.0 5.0 6.0 F i g u r e 5 . 4 0 . C o n v e r g e n c e of d e v i a t i o n of c o n s t r a i n e d K\ f r o m evolution equation i n m a x i m a l / i s o t h e r m a l s y s t e m . q in U l O m o O o m CD O c CD O) V_ CD £ ° in CN Q — B — B 8 h : 4h A — A — A 4h 2h H 1 ^ 2h h 0.0 1.0 2.0 3.0 t 4.0 5.0 6.0 Figure 5.41. Convergence of difference between constrained and freely evolved a in maximal/isothermal system. q in m 1 o ^. o M— O CD O c in CD K) O) i _ CD AU o O O m CN 0.0 B—B—B 8 h : 4h A — A — A 4h : 2h 4 2 h : h 1.0 2.0 3.0 t 4.0 5.0 6.0 Figure 5.42. Convergence of difference between constrained and freely evolved K9g in maximal/isothermal system. by a d d i n g an inner i teration at each t i m e step to e l iminate the need for e x t r a p o l a -t i o n . In fact, several c o d i n g errors a n d i n c o r r e c t difference expressions have been d i s c o v e r e d in the p r o g r a m w h i c h generated these results by e x a m i n i n g q u a n t i t i e s s u c h as the deviat ions s h o w n in F i g s . 5.37-42, or the changes in a grid f u n c t i o n as the m e s h is refined. C o n v e r g e n c e tests have been used i n n u m e r i c a l r e l a t i v i t y , b u t t h e rate of convergence is s e l d o m , if ever, m e a s u r e d . M o r e often, t h e tests are q u a l i t a t i v e in nature, a n d a solution is c o n s i d e r e d convergent if it " d o e s n ' t change m u c h " fol lowing a m e s h refinement. H o w e v e r , if one wishes to c o n s t r u c t a fully s e c o n d or der scheme, as was the goal here, it seems i m p o r t a n t to p e r f o r m a q u a n -t i t a t i v e m e a s u r e m e n t of the convergence rate. A m i s t a k e m a d e in a second o r d e r difference expressions can easily lead to an expression w h i c h is first order accurate. S u c h a n error wil l a lways show u p in a q u a n t i t a t i v e test b u t m a y easily go u n n o t i c e d in a q u a l i t a t i v e one. In c o n c l u d i n g this c h a p t e r , a few a d d i t i o n a l r e m a r k s c o n c e r n i n g the i n t e r p l a y of the c o n s t r a i n t a n d evolution e q u a t i o n s i n n u m e r i c a l r e l a t i v i t y seem in order. T h e results f o u n d here suggest t h a t the technique of m o n i t o r i n g the q u a l i t y o f a freely e v o l v e d solution using the d e v i a t i o n s f r o m the c o n s t r a i n t s is quite reasonable. T h i s p r o c e d u r e has been c r i t i c i z e d by Stewart [89], but his objections are based on the a s s u m p t i o n of inconsistency. It is not o b v i o u s , however, t h a t these d e v i a t i o n s c a n be used d i r e c t l y to p r o v i d e a q u a n t i t a t i v e estimate of the error in the s o l u t i o n as is s o m e t i m e s done [10,59]. T h e fact t h a t the c o n s t r a i n t equations exist a n d can be used i n this m a n n e r as a consistency check for a free e v o l u t i o n is often t o u t e d as a s p e c i a l feature of E i n s t e i n ' s equations. T h i s is true to a certain extent, b u t there is p r o b a b l y no m o r e i n f o r m a t i o n to be g a i n e d f r o m e x a m i n i n g c o n s t r a i n t d e v i a t i o n s t h a n f r o m e s t i m a t i n g the l o c a l t r u n c a t i o n error of the difference scheme u s i n g a m e s h r e f i n e m e n t t e c h n i q u e . T h e latter m e t h o d c a n be a p p l i e d to any system of difference e q u a t i o n s . In a d d i t i o n , researchers w h o c o n s t r u c t c o n s t r a i n e d schemes s h o u l d note 182 that deviations from the evolution equations will also provide a measure of the quality of the solution. 183 C H A P T E R 6 A n A p p l i c a t i o n of L o c a l M e s h Refinement 6.1 Introduction. A generic feature of al l of the s o l u t i o n s of the field e q u a t i o n s w h i c h have b e e n d i s p l a y e d so far is the relat ively s m a l l p r o p o r t i o n of the s p a c e t i m e grid in w h i c h there is significant spatial or t e m p o r a l v a r i a t i o n in the d y n a m i c variables. O f c o u r s e , this is largely due to the fact t h a t the spacet imes c o n t a i n only one or two pulses of scalar r a d i a t i o n , whose spatial extent is fairly s m a l l c o m p a r e d to that of the c o m p u t a t i o n a l d o m a i n . H o w e v e r , f r o m the c u r r e n t t h e o r e t i c a l u n d e r s t a n d i n g of s t r o n g g r a v i t a t i o n a l fields a n d the p r o d u c t i o n of g r a v i t a t i o n a l r a d i a t i o n [93], it seems clear that this will also be a feature of m a n y of the most interesting astrophysical a p p l i c a t i o n s of n u m e r i c a l relat ivi ty . P a r t i c u l a r l y for those cases w h i c h involve one or m o r e black holes (or other col lapsed objects) the strong field region, where the d y n a m i c s are m o s t difficult to c o m p u t e , w i l l have a length scale m a n y t i m e s smaller t h a n the distance to w h i c h the p r o d u c e d r a d i a t i o n should ideally be followed. T h i s s i t u a t i o n is not u n i q u e to relat ivi ty . M a n y other fields of c o m p u t a t i o n a l p h y s i c s i n v o l v e t i m e d e p e n d e n t p r o b l e m s w h i c h are significantly m o r e difficult to a p p r o x i m a t e in s o m e portions of the s o l u t i o n d o m a i n t h a n in others. T h i s has led t o r e s e a r c h , w h i c h has been a c c e l e r a t i n g recently, d i r e c t e d towards the d e v e l o p m e n t of techniques where the mesh s t r u c t u r e of the grid is a d a p t e d to the features of the s o l u t i o n . T h e m a j o r m o t i v a t i n g factor for this research is the gain in c o m p u t a t i o n a l efficiency which stands to be made. This is particularly true for problems in several spatial dimensions—Hedstrom and Rodrique [38] have expressed the view that some multi-dimensional calculations may not even be possible without the development of adaptive-grid methods. Moreover, if the degree of reliability of numerical solutions of complicated systems of PDE's is someday to approach that currently attained in the solution of ODE's, for example, then the development of such methods would seem to be essential. To date, there have been basically two approaches taken to the construction of adaptive-grid algorithms for time dependent PDE's [38]. The first method involves tracking solution features by means of a coordinate transformation, which is usually performed numerically. The discrete equations are then solved on a uniform mesh in the transformed coordinate system. The major problem which must be solved is how to choose the transformation so as to put the mesh points where they are needed. An example of the use of such a technique for the study of self-gravitating (Newtonian) spherically symmetric and axially symmetric gas flows is discussed by Tscharnuter and Winkler in [94]. The other approach is to keep the coordinate system fixed and introduce finer meshes in regions where the solution is relatively difficult to approximate. Recently, Berger and Oliger [4] have presented a general algorithm based on this idea, suitable for the treatment of hyperbolic systems in one or more spatial dimensions. This chapter describes an implementation of a variant of their algorithm for the solution of the self-gravitating scalar field in the polar/radial coordinate system. The results of the study are promising and suggest that mesh refinement techniques could play an important part in the development of numerical relativity. 6.2 Previous Approaches to Grid Adaptation in Numerical Relativity. Much of the work in numerical relativity has dealt with collapsing or collapsed objects in spherical or cylindrical coordinates, and there has been a general need for increased resolution in the central areas of the solution domain. One method which has been used quite frequently is to employ a non-uniform mesh in the radial direction. Thus the mesh spacing becomes a function of mesh index; typically a "geometrically stretched" grid which satisfies Ar y = q A r y _ , (6.1) for some constant q > 1 is used [28,57-61,70]. There are a few problems with this approach, however. Unless q satisfies ? = l + ?oAr 0 (6.2) for some other constant go as Aro —• 0, then the usual centred difference formu-lae will no longer have C?(Ar2) truncation error. Second order accuracy may be retained by using larger computational stars, but this is rarely, if ever done in prac-tice. Piran has previously suggested [70] that a condition like (6.2) should be used to determine the "stretching factor" q. Secondly, if an explicit difference scheme is used, then the maximum allowable time step is determined by the finest spacing Aro- Thus, small time steps are taken in regions where larger steps would presum-ably suffice and the gain in computational efficiency is not as complete as possible. Finally, this type of grid structure is limited in its ability to adapt to the solution. In particular, the highest density of grid points may not always be needed at the smallest radii, as is evident from the strong field results displayed in Chapter 4. In addition, in traversing across the coarsening mesh, outgoing pulses of radiation will suffer increasingly from dispersive effects and features which developed on the scale of the well-resolved interior region are likely to be lost. Another technique which is often used in hydrodynamical numerical relativity codes is to introduce a grid velocity to allow the mesh points to adapt to the matter flow. However, unless the velocity field is constant, the mesh spacing will n o t be u n i f o r m a n d a c c u r a c y wil l be lost, unless condit ions a n a l o g o u s to (6.2) are m a i n t a i n e d f r o m slice to slice. It s h o u l d be n o t e d t h a t the i d e a of u s i n g c o o r d i n a t e t r a n s f o r m a t i o n s for grid a d a p t a t i o n , w h i c h was m e n t i o n e d previously , wil l p r o b a b l y b e of l i m i t e d use i n n u m e r i c a l r e l a t i v i t y in the near future. S u c h t r a n s f o r m a t i o n s w i l l general ly be i n c o m p a t i b l e w i t h the c o o r d i n a t e c o n d i t i o n s i m p o s e d to s impli fy t h e e q u a t i o n s . F i n a l l y , M a n n has done some work on grid a d a p t a t i o n in c o n j u n c t i o n w i t h his a p p l i c a t i o n of finite element techniques to spherically s y m m e t r i c , perfect fluid col lapse [49]. O n e significant advantage of finite element m e t h o d s is t h a t a u n i f o r m p l a c e m e n t of nodes (grid points) is not as i m p o r t a n t for a c c u r a c y or ease of i m p l e -m e n t a t i o n as it is for finite difference m e t h o d s . M a n n used a fixed n u m b e r of nodes i n his r u n s , a l lowing t h e m to m o v e f r o m one t i m e step to a n o t h e r , a n d the m o t i o n of the i n f a l l i n g m a t t e r was followed quite well. 6.3 An Unsuccessful Attempt at Local Mesh Refinement. T h e first a t t e m p t s to use a local m e s h refinement technique for the self-g r a v i t a t i n g scalar field were based on a m e t h o d suggested by B r o w n i n g , K r e i s s a n d O l i g e r ( B K O ) [7]. T h e basic s t u c t u r e of the grid in the spatial d i r e c t i o n is shown i n F i g u r e 6.1 for the case of a 3:1 ref inement (H — 3h). T h e coarse g r i d is defined by G H = {Xj+(J -l-)H, .7 = 0,1,...} (6.3) w h e r e X1 is the interface c o o r d i n a t e a n d the fine grid is g f t = { X j - [j- l)fc, y = 0 , l , . . . } (6.4) 187 x x x x x x x ^ x l x ^ X X F i g u r e 6.1. G r i d s t r u c t u r e for the B K O refinement scheme. N n N o t e t h a t the grid p o i n t s of the fine grid are n u m b e r e d r i g h t to left. L e t U a n d u . J J d e n o t e the coarse a n d fine g r i d a p p r o x i m a t i o n s of a n y f u n c t i o n w h i c h is governed by a n e v o l u t i o n e q u a t i o n . P r o v i d e d t h a t the difference scheme is based o n a s y m m e t r i c s t a r w i t h d i a m e t e r < 3, t h e n a d v a n c e d values U*f+1 a n d tt™ + 1, for J > 0, j > 0 c a n be d e t e r m i n e d in the usual fashion. T h e p r o b l e m r e m a i n s to d e t e r m i n e the b o u n d a r y values U^+1 a n d «Q+ 1- C o n s i d e r the case w h e n the two grids are al igned i n t i m e as they are in F i g u r e 6 . 1 . T h e n B K O suggest the m a t c h i n g c o n d i t i o n s TTN+1 jTN+l n+l n+l . . U 0 + U 1 = uo + u i ( 6 - 5 ) TTN+1 TTN+l n+l n+l U — U u — u 1 ff ° - ° h 1 (6 -6) w h i c h are based on the c o n t i n u i t y of u a n d its first spatial derivat ive at Xj. Clearly these c o n d i t i o n s d e t e r m i n e the b o u n d a r y values once the interior values have been c a l c u l a t e d . If the g r i d is refined in the same m a n n e r i n t i m e as i n space, t h e n this p r o c e d u r e c a n only be used every H/h steps on the fine m e s h . B K O suggest that at t h e other steps, some f o r m of i n t e r p o l a t i o n f r o m p r e v i o u s l y c o m p u t e d coarse grid values c a n be used to d e t e r m i n e « o + 1 - T h e a u t h o r s established the s t a b i l i t y of the c o m b i n e d difference e q u a t i o n s for the case of leap-frog differencing of h y p e r b o l i c s y s t e m s w h e n the refinement was m a d e only in space. T h e usefulness of the tech-n i q u e for a p r o b l e m w i t h a b o u n d a r y layer was d e m o n s t r a t e d , a n d finally it was s h o w n t h a t there w o u l d be p r o b l e m s if a wave w h i c h was p o o r l y represented on the c o a r s e m e s h p r o p a g a t e d across the interface. T h e basic m e t h o d o u t l i n e d a b o v e was used to c o n s t r u c t a t r i a l p r o g r a m to solve the wave e q u a t i o n in spherical s y m m e t r y w i t h an a r b i t r a r y n u m b e r of levels of r e f i n e m e n t , w i t h the p r o v i s o t h a t the mesh size was a n o n d e c r e a s i n g function of r . T h e refinement was p e r f o r m e d in t i m e as well as space a n d each level of m e s h was a c t u a l l y one of the usual t e m p o r a l l y a n d spatial ly staggered grids. T h e c o m b i n a t i o n of t i m e refinement a n d staggering significantly c o m p l i c a t e d the task of d e t e r m i n i n g c o n d i t i o n s for b o u n d a r y values at the interface, a n d the flow of the p r o g r a m was quite i n v o l v e d . H o w e v e r , the m e t h o d seemed to work q u i t e well for t h e wave e q u a t i o n so it was i n c o r p o r a t e d into the schemes d e s c r i b e d i n Sees. 4.6 a n d 4.10 w h i c h solve the full m o d e l p r o b l e m in the c o o r d i n a t e systems where 3 = 0. A d d i t i o n a l c o m p l i c a t i o n s were i n d u c e d by the increase in n u m b e r of f u n c t i o n s to be d e a l t w i t h a n d the fact t h a t some of the quantit ies were d e t e r m i n e d by i n i t i a l value or b o u n d a r y value e q u a t i o n s rather t h a n evolut ion equations. T h e results, even for weak pulses, were less t h a n satisfactory. T h e represen-t a t i o n of the scalar field c o n t i n u e d to be g o o d everywhere, b u t the metric functions a n d extrinsic c u r v a t u r e c o m p o n e n t s d e v e l o p e d irregularit ies at the m e s h interfaces w h i c h s o m e t i m e s led to d i s c o n t i n u o u s b e h a v i o u r . T h e exact cause of the trouble w a s never l o c a t e d , as the c o m b i n a t i o n of the v a r i o u s u p d a t i n g p r o c e d u r e s used to d e t e r m i n e the interface values was quite c o m p l i c a t e d . S o m e a t t e m p t s were m a d e to r e m e d y this p r o b l e m by h a v i n g the coarse a n d fine m e s h overlap rather t h a n s i m p l y a b u t , b u t these were largely unsuccessful. H o w e v e r , a p a r t f r o m the p r o b l e m s at the interface it was clear t h a t a m e s h refinement p r o c e d u r e w o u l d be useful in better r e s o l v i n g the b e h a v i o u r in the v i c i n i t y of r = 0. In p a r t i c u l a r , i n d i c a t i o n s such as t h e a m o u n t of v a r i a t i o n in a q u a n t i t y , w h i c h should have been c o n s t a n t at large r d u r i n g the p e r i o d of self-reflection of the i n g o i n g pulse, were e x a m i n e d . A mesh r e f i n e m e n t r e d u c e d the fluctuations f r o m the case of a single coarse grid to a level comparable to that achieved by using a fine grid everywhere. At the same time, the local refinement solution required considerably less computer resources than a fine grid treatment. 6.4 The Mesh Refinement Algorithm of Berger and Oliger. Apart from the difficulties encountered in trying to implement the scheme of BKO for the general relativistic problem, the method suffers from several of the shortcomings outlined in Sec. 6.2 in the discussion of the use of a geometrically stretched grid. In particular, the grid does not fully adapt to the solution and the coarse grid must be sufficiently fine to adequately represent the initial ingoing pulse, unless a refinement extends to include the support of the initial data. In addition, it was clear that the irregular behaviour which developed in the geometric variables at the interfaces was due to deficiencies in the algorithm for updating boundary values when disturbances were passing through the interface. A rather obvious potential solution for all of these difficulties was to attempt to move the regions of refinement along with the scalar waves so as to minimize the propagation of signals across the mesh boundaries. Fortunately, the feasibility of such an approach had been established by Berger and Oliger [4] in an article which describes a very general algorithm for adaptive mesh solutions of hyperbolic systems. The basic idea of their method is to employ a single uniform coarse mesh which covers the entire computational domain, over-layed with uniform fine meshes in those areas requiring refinement. The refinement is performed in time as well as in space. Any grid can admit a properly contained subgrid and there can be an arbitrary number of levels of refinement, each charac-terized by some unique mesh spacing. The interior equations are evolved separately on each grid. After a time step has been taken on all grids at some level, evolution at that level is suspended until the equations on every grid at finer levels have been integrated to the same time. Boundary values may be determined from previously computed values in the case of overlapping grids at the same level, or from inter-polation from coarser grid values as in the previously described scheme. Whenever the level / systems have been updated to a time which aligns with the coarser, level I — 1 grid, those function values on the coarse grid which lie in refined regions are replaced with values transferred or interpolated from the fine grid. This process, known as injection, ensures that the representation of any function on any grid is as accurate as possible. This is useful for the purposes of displaying the results, but more importantly, it is clear that without doing this, any numerical difficulties which arise when a single coarse mesh is used to compute the solution would also arise on the coarse mesh during the refinement solution. If these problems do not remain localized, then there is the danger that the entire solution could become corrupted. In effect, it is not really necessary to evolve coarse grid equations at points which are in the interior of a refinement, but the programming is simpler if all equations are evolved, and the waste of computational effort is essentially negligible, particularly for multi-dimensional problems. The placement of refinements in Berger and Oliger is done in a very general and automatic fashion. The idea is to periodically generate local error estimates using a technique based on the Richardson hypothesis discussed in the previous chapter. Regions where a new refinement is needed, or where a previously in-troduced fine mesh is no longer needed are identified and the spatial domain is regridded by deleting and adding meshes from the grid structure. Initial values on a new fine mesh may be interpolated from a coarse grid or transferred from another fine grid prior to its deletion. This process is quite easy for problems in one space dimension. For the case of two or more dimensions, the authors have developed a very interesting algorithm for clustering points into arbitrarily oriented rectangular grids, in order to minimize the total area of the refinement. 6.5 Adaptation of the Algorithm for The Self-gravitating Scalar Field. Even for the case of a single spatial dimension, a general implementation of Berger and Oliger's algorithm demands a substantial programming effort beyond that which is required to implement the basic difference scheme. In addition, there were additional issues which had to be addressed in the solution of the general relativistic problem which were not treated in [4]. As a result, the program which is described in the remainder of this chapter can currently handle only a single refined mesh overlaying the global coarse mesh. In addition, the positioning of the fine mesh is not determined using a local error estimation technique. Instead, the initial data is constructed so that there is only a single ingoing pulse, and the center of the refinement tracks the extremal value of the wave packet. These restrictions simplify the implementation of the algorithm significantly but it should be evident from the following that they are relatively unimportant to the success of the application. 6.6 Difference Equations for the Mesh Refinement Program. The three computational examples of the mesh-refinement algorithm given in [4] all use one-level difference schemes which, among other things, keeps the amount of interpolation which must be done at regridding time to a minimum. At the same time, it was clear that the use of a fully staggered mesh structure in conjunction with the algorithm would unduly complicate the program as it had in the previous attempt at mesh refinement described in Sec. 6.3. Therefore, it was decided to retain the spatial staggering, but to define all functions at the same time level. An example of the mesh structure for the case of 3:1 refinement is shown in Fig. 6.2. It is necessary to use an odd (3,5,7,...) refinement factor to get a staggered fine mesh which is properly aligned with the coarse mesh. Now, it is not an easy matter to design an efficient, fully second-order scheme for Einstein's equations based on two time levels of this type of structure. Such a scheme must obviously be implicit and the large, sparse (but not tri-diagonal) non-linear system which would result would surely require the equivalent computational work of many time steps of an explicit method. This suggested that the possibility of using a two-step method for the scalar field/Einstein equations should be investigated. o X o °x° O X O » 0 » X 0 X o X ° O X O "©» X 0 x ° x ° O X O o X o o X o o x o °x° O X O X 0 X x © x X 0 x O X O °x° O X O x 0 x xQx x 0 x O X O °x° O X O o X o 0 X 0 O X O °x° O X O X 0 X x 0 x x 0 x O X O °x° O X O x 0 x x 0 x x O x O X O °x° O X O 0 X 0 o X o O X O °x° O X O X O X x 0 x X O X O X O °x° O X O x 0 x x 0 x x 0 x O X O °x° O X O o X o 0 X 0 o x o •X' X O X xQx O X O °x° X o x x 0 x O X O °x° 0 X 0 Figure 6 .2 . Typical grid structure for the mesh refinement program for the case of a 3:1 refinement. Such a method is easily constructed in the polar/radial coordinate system used in sections 4.10 and 5.4. Because the spatial staggering is retained, the same difference equations can be used to solve the slicing equation for a and the Hamil-tonian constraint for a. The evolution equations for the scalar field variables are differenced as K n ; ^ A i [ ^ ( ^ ) \ U ) ,6.8, The advanced values of the field, <f>™+1 are then computed using the second order outgoing radiation condition and $ n + [ = Ar+<f>n+1. Similarly, the discrete evolution equations for the geometric variables are (6.10) (6.9) where, again, A = (ln a)' and sn. is given by ( 5 . 3 4 ) . a n + is updated using ( 5 . 3 8 ) , and all of the boundary conditions are the same as is previous calculations. The momentum constraint is again trivially soluble for K rr. Note that the evolution equation for K rr is easier to solve here than on the fully staggered mesh since the A linearized stability analysis of ( 6 . 7 ) and ( 6 . 8 ) suggests that the system will be Von Neumann-stable provided which is also the CFL condition. This condition on the time step is twice as restric-tive as it was previously, but is still acceptable, and is the price that must be paid for retaining the spatial staggering. The scheme is also non-dissipative. Time sym-metric initial data is readily determined for this system using essentially the same technique described in Sec. 4.10. However, because the difference method used here is a two-level scheme, initial values at t = — At , as well as at t = 0 are needed. For the scalar field variables, these values are computed using a Taylor series expansion about t — 0, and the difference equations (6.7) and (6.8). The geometric variables at t — — At can then be calculated from the slicing and constraint equations as they are on the t = 0 slice. The time symmetric initial data is acceptable for the purposes of testing the overall difference scheme, but for the mesh refinement problem, initial values which produce an ingoing-only pulse are required. In the case of flat spacetime, such a advanced value, KTT n + appears linearly in (6.10). « . 1 * - At < - Ar a ~ 2 (6.11) pulse in the spherically symmetric scalar field may be generated by freely specifying the initial value of <f> <f>{r,0) = S ( r ) (6.12) and demanding that rj> - (r<p)' = 0 (6.13) This then determines the initial value of II n(r,0) = ^(r,0)= *(r/) (6.14) This same procedure can be used in the self-gravitating system. Expressions (6.12) and (6.14) are used to specify $° L and 11° and the rest of the initial data is calculated in the same manner as for the time symmetric initial data. The above difference equations, along with the slicing and constraint equations from Sec. 4.10, were implemented in a program designed solely to test the overall difference scheme. The results were very satisfactory and quite similar to those generated by the program described in Sec. 4.10. As an illustration of the stability and non-dissipativity of the code, Figs. 6.3 and 6.4 show plots of the scalar field momentum and lapse function generated in a run where the sign of the time step was flipped after half the requested time steps had been taken. In the ideal case, the initial data should be recovered at the end of the run. Even though the mesh is fairly coarse (Ar = 4 At = 1.0) and the ingoing pulse is significantly self-gravitating, the final data agrees with the initial values to within ~ 1%. (Only the inner half of the actual computational domain is displayed in the figures. The time reversal test will fail if the outgoing pulse reaches the outer boundary since the finite differenced radiation boundary conditions are not invariant under At —• — At.) The program was also used to verify the stability condition (6.11) for a range of field strengths. Figure 6.3. "Time reversal" test of difference scheme in polar/radial system: a. Note that the grid used is not temporally staggered. 6.7 The Mesh Refinement Algorithm in Detail. The flow of the refinement algorithm is reasonably straightforward—a pseudo-code form is given in Fig. 6.5. The major difference between this program and those that are discussed in [4] is that boundary value equations, as well as evolution equa-tions, have to be solved here. Since the solution to a boundary value problem cannot be determined locally, it is slightly trickier to deal with the Hamiltonian constraint and slicing equation in the refinement scheme than the evolution equations. In order to treat these equations, it is convenient to view the spatial domain at any time step as being composed of two or three segments. Fig. 6.6 depicts the three possible types of segment structure for the case of a single refinement; in each in-stance, the shaded area represents the region of refinement. The program maintains a representation of this structure which includes, for each segment, the mesh size, a pointer to the beginning of the segment, and the number of grid points contained. The slicing equation and Hamiltonian constraint are always integrated outwards segment by segment, with the final values of one segment being used to provide the initial value for the next segment. This procedure produces the most globally accurate solutions to the boundary value problems. It is also necessary to do this so that the numerical difficulties described in Sec. 4.13 do not crash the code on the coarse mesh. A global solution of the slicing and constraint equations is only possible when there is temporal alignment of the fine and coarse meshes (step 14 of Fig. 6.5). Otherwise, the equations are partially integrated on the coarse mesh only if the first segment is coarse (step 8), and on the fine grid (step 12). In both cases the lapse function is fixed by specifying the last value in the segment. This value is determined by interpolation in step 12, and by extrapolation of previous coarse grid values in step 8. Following a global solution, the lapse is rescaled in the usual fashion. It should be mentioned here that linear or bilinear interpolation generally proved adequate for determining fine grid boundary values. 1: Input run parameters. 2: Compute extent of mesh refinement. S: Initialize grid structure. 4: Compute initial data on both levels. 5: LOOP for requested number of steps. 6: Evolve coarse grid quantities. 7: IF first segment is coarse T H E N 8: Solve constraint and slicing equations on first segment. 9: E N D IF 10: LOOP refinement factor times. 11: Evolve fine grid quantities. 12: Solve constraint and slicing equations on fine segment. 13: END LOOP 14: Solve constraint and slicing equations on all segments. 15: Inject fine grid values into coarse grid. 16: IF refinement should be moved T H E N 17: Move fine grid and reinitialize fine grid functions. 18: END IF 19: END LOOP Figure 6.5. Pseudo-code form of mesh refinement algorithm. The evolution equations (steps 6 and 11) are essentially treated as suggested by Berger and Oliger. The auxiliary quantities A™ and , , never require J + 2 3+2 boundary conditions, but the corresponding primary variables, a" and <f>n , , as 3 3+2 well as II™ and K l n, do. On the coarse grid, these conditions are the ones used for 3 3 the case of a single mesh. On the fine grid, the boundary values are always (because of the spatial staggering—see Fig. 6.2) determined from interpolation using coarse grid values unless the positioni ng of the mesh is such that a regularity or outgoing radiation condition should be used. 198 First segment Second segment Third segment fine coarse Figure 6.6. The three possible positionings of the fine grid and the corresponding segment structures. The run parameters (step 1) include the usual initial data parameters for the field, the mesh refinement factor p, the number of time steps to be taken on the coarse mesh and a quantity which determines the extent of the fine grid. The fine mesh is centred at the peak of the initial pulse, (r = rQ) and extends a distance w/2 to the left and right such that <5(r0±w;/2,O) < c,*(r 0,0) (6.15) The size of the refinement stays fixed during the computation. Periodically (step 16) the grid index, m, of the maximum value of \<f>™\ is located and, if necessary, the grid is translated so that the center of the fine grid is at r « rm (Step 17). The reinitialization of the fine grid quantities is accomplished through a combination of shifting storage array values, and linearly interpolating values from the coarse grid. 6 . 8 Results From the Mesh Refinement Program. The overall flow of the mesh refinement program was first tested using only the scalar field variables. An example of the results generated is shown in Figs. 6.7-9. The coarse mesh has 20 points, with Ar = 4 At = 1.0 and a 3:1 refinement which covers about 50% of the coarse grid is used. The first two plots show $ as calculated without and with the refinement respectively. The improved quality of the solution in the refined case is clear. The ringing behavour and dispersion of the packet evident in the first graph is essentially eliminated, as is the small outgoing pulse due to large truncation errors in the differenced Taylor series expansion used to calculate the initial data. The third plot shows the fine grid representation of output every few steps on the fine mesh. The four basic stages of the evolution are discernible in the qualitative changes in the character of the plot. In the first stage the mesh is moving inward, and the periodic regridding leads to the corrugated appearance in Once the left edge of the fine grid reaches r = 0, it stops and the 200 Refinement CPU time per Ratio grid point (ms) 1 : 1 0.097 3:1 0.106 5 : 1 0.100 7 : 1 0.099 9 : 1 0.097 Table II. Illustration of regridding overhead in the mesh refinement program. function looks smooth for a period. The refinement then moves outward until it reaches the outer boundary; the pulse then propagates off of it. There were no fundamental difficulties encountered in incorporating the geo-metric variables into the scheme, although it did take some care to get all the various cases of boundary conditions working properly. Results obtained using a mesh re-finement are compared with those generated on a coarse grid alone in Figs. 6.10-19. In each case the top and bottom plots are from the unrefined and refined runs, respectively. There are significant self-gravitational effects in this test, the mass of the spacetime is « 0.92 and at times this energy is contained within a sphere of about 4 Schwarzschild radii. The coarse grid has mesh sizes Ar = 4 At = 1.0 and a 5:1 refinement which covers about 1/2 of the spatial domain is used. The improvement which is achieved by using the refinement is particularly vivid for the freely evolved geometric variables. The refinement results are indistinguishable from a solution generated on a single fine mesh, but the computational cost of the former is about 25% of the latter. In fact, there is not much overhead incurred in the regridding procedure. Table II shows a list of timing estimates, per space-time grid point, for various levels of refinement using an Amdahl 5860, the IBM FORTRAN-H Compiler and 8-byte real arithmetic. (These timing estimates are very nearly independent of the particular initial data used since, apart from the 201 Figure 6.8. Mesh refinement test for scalar field alone: of $ determined using 3:1 refinement. Coarse grid representation 202 Figure 6.9. Mesh refinement test for scalar field alone: Fine grid representation of 203 Figure 6.10. Mesh refinement test for self-gravitating scalar field: <f> using coarse grid alone. Figure 6.11. Mesh refinement test for self-gravitating scalar field: <j> using 5:1 refinement. 204 Figure 6.12. Mesh refinement test for self-gravitating scalar field: constrained a using coarse grid alone. 205 Figure 6.15. Mesh refinement test for self-gravitating scalar field: constrained Krr using 5:1 refinement. 206 *0.0 Figure 6.17. Mesh refinement test for self-gravitating scalar field: freely evolved a using 5:1 refinement. Figure 6.19. Mesh refinement test for self-gravitating scalar field: freely evolved K rr using 5:1 refinement. iterative solution of the Hamiltonian constraint, the scheme is completely explicit.) The maximum overhead is about 10% and by the time the grid has been refined by a factor of 9, the extra operations needed to translate the refinement are essentially negligible in comparison to the work expended in solving the difference equations. Finally, the fluctuations in the mass aspect function at the outer edge of the grid, prior to the arrival of the pulse, is reduced by a factor of 24.5 by using the 5:1 refinement, which is a good indication of the second order convergence rate that is to be expected from the scheme. 6.9 The Nature of the Gravitational Self-Interaction of the Scalar Field. The mesh refinement program was used to examine more closely the features of the gravitational self-interaction of the scalar field which was described briefly in Sec. 4.12. Of all of the functions which were evolved or computed by the program, the spatial derivative of the mass aspect proved to be the most useful for displaying the qualitative features of the interaction. Recall from (4.104) that this quantity is also the radial energy density of the scalar field. Fig. 6.20 shows a series of plots of this quantity at successive instants of time during the period of self-reflection of a pulse with negligible gravitational self-energy. The evolution proceeds from left to right, top to bottom. The initial configuration for 6 was of the form The use of the fourth power in the exponential has the effect of flattening the top of the pulse and producing a better separation between the two bumps in the energy profile. At the beginning of the sequence, the leading edge of the pulse is undergoing m 2 = 47rr p — 27rr 2n 2 + $2 a2 (6.16) (6.17) self-reflection at r = 0, while the trailing edge is still incoming. As the evolution continues, the energy peaks pass through each other and, roughly speaking, the net effect is the transferral of some energy from the trailing edge to the leading edge of the pulse. It is readily verified from analytic considerations that this is the expected behaviour for a weak pulse. The situation for a self-gravitating packet is markedly different. Fig. 6.21 shows the time development of m' for a spacetime with a mass of ss 0.92. The coarse mesh extends to r = 100 and, again Ar = 4 At = 1.0. A 9:1 refinement which covers 40% of the spatial domain is used; thus the radial axes of the plots extend from r = 0 to r ss 43M. The initial stages of the sequence are quite similar to the weak field case. However, as the reflected leading edge of the pulse passes through the incoming trailing edge, a third peak appears in the energy profile (last three frames of Fig. 6.21a). This peak then appears to travel inwards, combining with the peak from the trailing edge to produce the final outgoing profile seen in the last few frames. Thus, in this case, one of the overall effects would seem to be that energy is transferred from the leading edge to the trailing edge of the pulse. The small bump which shows up in the final frames displayed is also quite curious. As an indication of the trustworthiness of these results, the mass as measured at the outer edge of the grid is conserved to within 0.1%. In addition there is very good agreement between the freely evolved and constrained geometric variables. Figs. 6.22 and 6.23 show the late time waveforms of the scalar field for the two cases discussed above. The initial spatial extent and positioni ng of both pulses were identical and the plots represent the same instant of proper time as measured at the outer edge of the grid to within 1%. The clearly visible lag in the position of the strong pulse as compared to the weak pulse is partly due to the fact that the former has to "travel farther" in the curved geometry. The second peak in the pulse (which looks more irregular than it actually is, due to the injection of the fine grid representation onto the coarse mesh) is probably due to the "backscattering" 210 Figure 6.20. Evolutionary sequence of m' (radial energy density) for scalar pulse with negligible self-gravitation during period of self-reflection off r = 0. Evolution proceeds left to right, top to bottom. Figure 6.21a. Evolutionary sequence of m' (radial energy density) for self-gravitating scalar pulse during period of self-reflection off r = 0. 2 1 2 A - x r ^ x A ^ x A^-x ^ ^ /v ^ A -A ^ A ^ A ^ A ^ A , A . A . A , A . A . A A Figure 6.21b. Continuation of Fig. 6.21a. 213 Figure 6.22. Late-time configuration of outgoing scalar pulse with negligible self-gravitation. The initial configuration was of the form given by Eq (6.17). Figure 6.23. Late-time configuration of previously self-gravitating scalar pulse. The initial configuration was of the same form as that of the weak pulse. effect. Although the long tail which has developed seems to be genuine, the small bump at r — 0 is due to the fact that the refinement starts moving outwards a bit too soon and a small amount of noise is generated at the refinement interface. Similarly, there is a small amount of "ringing" at r — 0 in the weak pulse, which was evolved without the use of a refinement. If the initial amplitude of the scalar field is increased by a few percent, the resulting spacetime almost certainly contains a black hole. In this case, the plots of m' look very similar to those of Fig. 6.21, except that the combination of the "backscattered" bump and the pulse from the trailing edge of the original packet produce a gravitational field sufficiently strong that they do not escape. This is a rather novel way to form a black hole since it would appear that neither of the original pulses in the energy profile would form one on its own. Although this interaction effect is not well understood yet, the fact that the third pulse appears at fairly large distances—r « 10M in the case above—suggests that one can rule out explanations based on scattering off of the peak in the effective potential for massless particles [54] which occurs at r = 3M in a Schwarzschild spacetime. The creation of pulses of radiation due to the interaction of other packets has been observed in the study of cylindrical gravitational waves [70,72], but it is not known whether these two phenomena are at all related. Finally, being able to employ a single level of refinement is insufficient to cure the numerical problems encountered in the polar/radial coordinate system, for the case of black hole spacetimes. Although the program runs longer with a refinement than without one, it still crashes due to the development of very steep gradients. However, the initial success of the refinement technique suggests that it may not be too difficult to keep refining the mesh in the vicinity of r = 2M in order to maintain smooth representations of the functions. Unfortunately, this might lead to a very rapid growth in the total number of (spacetime) grid points as a function of integration time, which would render the method computationally unacceptable. C H A P T E R 7 C o n c l u s i o n s There are several clear conclusions that can be drawn from this work regard-ing the evolution of Einstein's equations in general, and radiative spacetimes in particular. First it seems established that it is possible to evolve self-gravitating radiation in a uniformly second order, non-dissipative fashion through consistent use of central differencing. To a certain extent, this has been known since Eppley and Smarr's work on the two black hole collision [24,82]. However, the convergence of their scheme was never established, and for the most part, they worked in normal coordinates so the convective terms were eliminated. The results of Chapters 4 and 5 in the maximal/isothermal coordinate system show that both the Crank-Nicholson and angled schemes provide viable alternatives to Wilson's method [68,104]. Al-though a direct comparison with the latter has not been made yet, it is probable that the second order centred schemes are preferable for the treatment of radiation, since Wilson's method is always dissipative [28], and becomes essentially one-sided where steep gradients form. At the same time, the suitability of these two centred methods has only been established for a problem in one dimension. Additional work will be needed to see if they can be extended to multi-dimensional problems, while still maintaining 0(ng) run time, where ng is the number of grid points on a slice. In particular, there is the possibility that they could be used in conjunction with some sort of splitting technique [46,50], whereby a time step of a multi-dimensional calculation is performed as a sequence of updates, each of which involves only a single dimension. However, the leap-frog scheme, which was not applied to the self-gravitating scalar field, should be directly extendable to multi-dimensional problems, so it would be of interest to attempt to use it in conjunction with the maximal/isothermal coordinate system. The results of Chapter 3, where the three second order methods were seen to behave very similarly for the test problem, suggest that there should not be major difficulties in doing this. Also, in view of the success of the basic scheme used in Chapter 6, it may be worthwhile to consider constructing schemes on an "unstaggered" mesh if full second order accuracy is desired. This would eliminate the need for most of the averaging and extrapolation operations which were necessary in the methods discussed in Chapters 4 and 5. It should also simplify the treatment of some of the non-linearities in the evolution equations. The usefulness of the preliminary study of finite difference techniques for the scalar wave equation in general coordinates which was described in Chapter 3, cannot be overemphasized. When a complicated system of difference equations is implemented as a computer program and "instabilities" appear, it is very easy to suspect the basic difference scheme. Very often, however, the source of the problem will be one or more coding errors in the program. Once the stability and convergence properties of the basic difference schemes discussed in Chapter 3 had been established, it was possible to approach the full model problem with confidence that it would be possible to get a stable evolution using the centered methods. Another principal observation which has been made is that the errors in differ-ence solutions to Einstein's equations should be expected to have Richardson-type expansions in powers of the mesh spacing. This immediately resolves the "incon-sistency question", and also provides a very useful way of checking the implemen-tation of a scheme. Although this would hardly come as a surprise to most nu-merical analysts—there are many well known and widely used algorithms based on Richardson extrapolation—it may remove some of the misconceptions and doubts that various relativists have expressed about 3+1 numerical relativity. In particu-lar, it is probably reasonable to assume that any stable, consistent free evolution scheme will converge to a properly constrained solution of the field equations. If this is the case, then fundamentally, it should not matter which of the many conceivable algorithms, based on varying usage of the constraint and evolution equations, is employed. Of course, "fundamentally" has to be interpreted as "in the limit as the mesh spacings tend to 0". The problem is that, at any given time, the amount of computa-tion involved in answering the most interesting questions usually precludes working anywhere near such a limit. This has certainly been the case in most previous work in numerical relativity, and even in this research. Even if it does not matter, as a point of principle, whether the constraints are imposed or not, it still seems that constrained schemes are easier to construct (from the point of view of getting a stable evolution) and probably more accurate (Compare Figs. 6.15 and 6.19 for ex-ample) at a given level of discretization. Such observations have been made before [24,26,68] but good explanations of the behavior have not been provided. This is an interesting problem, because if one could figure out exactly why a constrained solution has better properties, the information might be useful in constructing bet-ter evolution schemes. The major difficulty with constrained schemes is that they usually require the solution of elliptic equations at every time step, and even with 0(ng) methods such as multi-grid [14,33], are apt to run several times slower than a corresponding free evolution method. In addition, state-of-the-art gravitational radiation calculations tend to use coordinate systems where certain combinations of geometric variables are directly identifiable (at least asymptotically) with the "real" radiative degrees of freedom. These quantities are invariably treated via evo-lution equations, and it is imperative that the corresponding difference schemes be as accurate as possible. From the point of view of developing more efficient and accurate methods for numerical relativity, the results of the previous chapter are very encouraging. At least for problems in one spatial dimension, it should not be overly difficult to construct complete implementations of algorithms such as Berger and Oliger's for the Einstein equations. There will be a continuing need for reliable one-dimensional programs in the future, if for no other reason than to provide a testing facility for more general codes. Although the biggest gains in efficiency would result from the application of mesh refinement to multi-dimensional problems, it would appear from [4] that some additional work needs to be done to allow for the solution of elliptic equations, as well as evolution equations. It is interesting to note that the primary data structure used in [4] could also be used with the multi-grid method, and it may be feasible to combine the algorithms. The two most serious unresolved problems in this work are the irregularities near r = 0 encountered in isothermal coordinates and the difficulties in the strong field regime, particularly with the Hamiltonian constraint. The first problem is particularly troublesome since it arises even for weak fields and thus seems to be indicative of a basic flaw in the difference scheme in the vicinity of r = 0. As a general rule, the application of special techniques near the origin was avoided, but perhaps more available information about the behavior of the various functions (regularity conditions) should be incorporated into the schemes. The problems with the Hamiltonian constraints in the maximally sliced coordinate systems deserve ad-ditional attention. As the results of Sec. 6.9 indicate, there will be some interesting physics going on with the scalar field as a black hole forms and maximally sliced spacetimes have been shown to be suitable for long time evolution of black hole spacetimes [66,76,77]. Despite these problems, I feel that this work represents a step towards the development of useful and reliable "computational laboratories" for the study of gravitational phenomena. B i b l i o g r a p h y [l] Arnowitt, R., S. Deser and C. W. Misner. "The Dynamics of General Relativ-ity.", in Gravitation-An Introduction t o Current Research, L. Witten, ed. New York: Wiley, 1962. [2] Bardeen, J. M. and T. Piran. "General Relativistic Axisymmetric Rotating Systems: Coordinates and Equations." Phys. Rep. 96 (1983), 205-250. [3] Berger, B. K., D. M. Chitre, V. E. Moncrief, and Y. Nutku. "Hamiltonian Formulation of Spherically Symmetric Gravitational Fields." Phys. Rev. D. 5 (1972), 2467-2470. [4] Berger, M. J. and J. Oliger. "Adaptive Mesh Refinement for Hyperbolic Par-tial Differential Equations." J. Comp. Phys. 53 (1984), 484-512. [5] Bowen, J. M., J. Rauber, and J. W. York. "Two Black Holes with Axisymmet-ric Parallel Spins: Initial Data." Class, and Quant. Grav. 1 (1984), 591-610. [6] Brill, D. R. "General Relativity: Selected Topics of Current Interest." Supp. Nuovo Cimento 2 (1964), 1-56. [7] Browning, G., H.-O. Kreiss and J. Oliger. "Mesh Refinement." Math. Comp. 27 (1973), 29-39. [8] Burden, R. L., J. D. Faires and A. C. Reynolds. Numerical Analysis. Boston: Prindle, Weber and Schmidt, 1978. [9] Carfora, M. and N. Tzoitis. "Some Remarks on Radiation Gauges in General Relativity." Phys. Lett. 88A (1982) 443-446. [10] Centrella, J. "Interacting Gravitational Shocks in Vacuum Plane-Symmetric Cosmologies." Astrophys. J. 241 (1980), 875-885. [11] Centrella, J. "Numerical Cosmologies.", in The Origin and Evolution o f Galaxies. B. J. T. Jones and J. E. Jones, eds. Boston: D. Reidel, 1983. 41-64. [12] Centrella, J. and R. A. Matzner. "Plane-Symmetric Cosmologies." Astrophys. J. 230 (1979), 311-324. [13] Choptuik, M. "A Study of Numerical Techniques for the Initial Value Problem of General Relativity.", Univ. of British Columbia M.Sc. thesis (unpublished), 1982. [14] Choptuik, M. and W. G. Unruh. "An Introduction to the Multi-grid Method for Numerical Relativists." To appear in Genera/ Relativity and Gravitation. [15] Christodoulou, D. Princeton Univ. thesis, 1971 (unpublished), cited in Unruh [95]. [16] Ciment, M. "Stable Difference Schemes with Uneven Mesh Spacings." Math. Comp. 25 (1971), 219-227. [17] Courant, R., K. Friedrichs and H. Lewy. "Uber die Partiellen Differenzen-gleichurgen der Mathematischen Physic." Math. Anna/en 100 (1928), 34-74. English translation in IBM J. (1967), 215-234. [18] Courant, R. and D. Hilbert. Methods of Mathematical Physics. Vbiume II: Partial Differential Equations. New York: Wiley and Sons, 1962. [19] Deruelle, N. and T. Piran, eds. Gravitational Radiation. Amsterdam: North-Holland, 1983. [20] D ISSPLA User's Manual. Integrated Software Systems Corporation, 1981. [21] Dupont, Todd. 'Mesh Modification for Evolution Equations." Math. Comp. 39 (1982), 85-107. [22] Eardley, D. M. "Theoretical Models for Sources of Gravitational Waves.", in Gravitational Radiation. N. Deruelle and T. Piran, eds. Amsterdam: North-Holland, 1983. 257-296. [23] Eardley, D. M. and L. Smarr. "Time Functions in Numerical Relativity: Marginally Bound Dust Collapse." Phys. Rev. D. 19 (1979), 2239-2259. [24] Eppley, K. R. " The Numerical Evolution of the Collision of Two Black Holes", Princeton Univ. Ph.D. thesis, 1975. [25] Eppley, K. R. "Evolution of Time-Symmetric Gravitational Waves: Initial Data and Apparent Horizons." Phys. Rev. D. 16 (1977), 1609-1614. [26] Eppley, K. R. "Pure Gravitational Waves.", in Sources of Gravitational Radiation. Cambridge: Cambridge University Press, 1979. L. Smarr, ed. 275-291. [27] Estabrook, F., H. Wahlquist, S. Christensen, B. DeWitt, L. Smarr, and E. Tsiang. "Maximally Slicing a Black Hole." Phys. Rev. D. 7, (1973), 2814-2817. [28] Evans, C. R., L. Smarr, and J. R. Wilson. "Numerical Relativistic Collapse and Collisions with Spatial Time Slices." (preprint, 1983). [29] Geroch, R. "Energy Extraction." Ann. NY Acad. Sci. 224 (1973), 108-117. [30] Griffiths, D. F. and A. R. Mitchell. The Finite Difference Method in Partial Differential Equations. New York: John Wiley and Sons, 1980. 221 [31] Gropp, W. D. "A Test of Moving Mesh Refinement for 2-D Scalar Hyperbolic Problems.", SIAM J. Sci. Stat. Comput. 1 (1980), 191-197. [32] Gustafsson, B., H.-O. Kreiss and A. Sundstrom. "Stability Theory of differ-ence Approximations for Mixed Initial Boundary Value problems. II.", Math. Comp. 26 (1972), 649-686. [33] Hackbusch, W. and U. Trottenberg, eds. Lecture Notes in Mathematics 960 Multigrid Methods. New York: Springer-Verlag, 1981. [34] Harten, A., P. D. Lax, and B. van Leer. "On Upstream Differencing and Godunov-Type Schemes for Hyperbolic Conservation Laws.", SIAM Rev. 25 (1983), 35-61. [35] Haugan, M. P. "Post-Newtonian Arrival Time Analysis for a Pulsar in a Bi-nary System." Astrophys. J. 296 (1985), 1-12. [36] Hawking, S. W. and G. F. R. Ellis. The Large Scale Structure of Space-Time. Cambridge: Cambridge University Press, 1973. [37] Hawley, J. F., L. Smarr and J. R. Wilson. "A Numerical Study of Nonspherical Black Hole Accretion II. Finite Differencing and Code Calibration." (preprint, 1983) [38] Hedstrom G. W. and G. H. Rodrique. "Adaptive-Grid Methods for Time-Dependent Partial Differential Equations.", in Lecture Notes in Mathe-matics 960 Multigrid Methods Ed. W. Hackbusch and U. Trottenberg. New York: Springer-Verlag, 1981. 474-484. [39] Henrici, P. Applied and Computational Complex Analysis. Vol I. Toronto: John Wiley and Sons, 1974. [40] Hirsh, R. S. "Higher Order Accurate Difference Solutions of Fluid Mechanics Problems by a Compact Differencing Technique.", J. Comp. Phys. 19 (1975), 90-109. [41] Israeli, M., and S. A. Orszag. "Approximation of Radiation Boundary Condi-tions." J. Comp. Phys. 41 (1981), 115-135. [42] Kreiss, H.-O. "Stability Theory for Difference Approximations of Mixed Initial Boundary Value Problems. I."Math. Comp. 22 (1968), 703-714. [43] Kreiss, H.-O. Numerical methods for Solving Time-dependent Prob-lems for Partial Differential Equations. Montreal: Led Presses de L'Universite de Montreal, 1978. [44] Kreiss, H.-O. and J. Oliger. Methods for the Approximate Solution of Time Dependent Problems. Garp Publications Series No. 10, 1973. 222 [45] Kulkarni, A. D. "Time-Asymmetric Initial Data for the N Black Hole Problem in General Relativity." J. Math. Phys. 25 (1984), 1028-1034. [46] LeVeque, R. J. and J. Oliger. "Numerical Methods Based on Additive Split-tings for Hyperbolic partial Differential Equations." Math. Comp. 40 (1983), 469-497. [47] Lucier, B. J. "A Stable Adaptive Numerical Scheme for Hyperbolic Conser-vation Laws." SI AM J. Numer. Anal. 22 (1985), 180-203. [48] MacCullum, M. et al; Schmutzer, E., ed. Exact Solutions of Einstein's Field Equations. Cambridge: Cambridge University Press, 1980. [49] Mann, Patrick J. "Applications of the Finite Element Method to General Relativity." Oxford Ph.D. thesis, 1982. [50] Marchuk, G. I. Methods of Numerical Mathematics . New York: Springer-Verlag, 1982. [51] Matzner, R. A., and L. C. Shepley, eds. Spacetime and Geometry: The Alfred Schild Lectures. Austin: University of Texas Press, 1982. [52] May, M. M. and R. H. White. "Hydrodynamic Calculations of General-Relativistic Collapse." Phys. Rev. 141 (1966), 1232-1241. [53] Miller, J. J. H. "On the Location of Zeros of Certain Classes of Polynomials with Applications to Numerical Analysis." J. fnst. Maths. Applies. 8 (1971), 397-406. [54] Misner, C. W., Thorne, K. S. and Wheeler, J. A. Gravitation. San Francisco: W. H. Freeman and Co., 1973. [55] Mitchell, A. R. and D. F. Griffiths. The Finite Difference Method in Partial Differential Equations. Chichester: Wiley, 1980. [56] Moore, C. ed. Discrete Fourier Transforms with Applications to Fourier Series and Convolution Integrals. University of British Columbia Com-puting Center Documentation, 1984. [57] Nakamura, T. "General Relativistic Collapse of Axially Symmetric Stars Lead-ing to the Formation of Rotating Black Holes." Prog. Theor. Phys. 65 (1981), 1876-1890. [58] Nakamura, T. "Numerical Relativity." Supp. Prog. Theor. Phys. 70 (1981), 202-214. [59] Nakamura, T., K.-I. Maeda, S. Miyama and M. Sasaki. "General Relativistic Collapse of an Axially Symmetric Star. I." Prog. Theor. Phys. 63 (1980), 1229-1244. [60] Nakamura, T., K.-I. Maeda, S. Miyama and M. Sasaki. "Numerical Integration of Axially Symmetric Systems.", in Proceedings of the Second Marcel Grossmann Meeting on General Relativity. R. Ruffini, ed. New York: North Holland, 1982. 675-677. [61] Nakamura, T. and H. Sato. "General Relativistic Collapse of Non-Rotating, Axisymmetric Stars." Prog. Theor. Phys. 67 (1982), 1396-1405. [62] Nicol, T. ed. U B C M A T R I X A Guide to Solving Matrix Problems. University of British Columbia Computing Center Documentation, 1982. [63] O Murchadha, N. and J. W. York Jr. "Initial Value Problem of General Rela-tivity. I. General Formulation and Physical Interpretation." Phys. Rev. D. 10 (1974), 428-436. [64] Ortega, J.M. and W. C. Rheinboldt. Iterative Solutions of Non-Linear Equations in Several Variables. New York: Academic Press, 1970. [65] Osher, S. and R. Sanders. "Numerical Approximations to Nonlinear Conser-vation Laws with Locally Varying Time and Space Grids." Math. Comp. 41 (1983) 321-336. [66] Petrich, L. I., S. L. Shapiro and S. A. Teukolsky. "Oppenheimer-Snyder Col-lapse with Maximal Time Slicing and Isotropic Coordinates." Phys. Rev. D. 31 (1985), 2459-2469. [67] Piran, T. "Gravitational Waves and Gravitational Collapse in Cylindrical Sys-tems.", in Sources of Gravitational Radiation. L. Smarr, ed. Cambridge: Cambridge University Press, 1979. 409-422. [68] Piran, T. "Numerical Codes for Cylindrical General Relativistic Systems." J. Comp. Phys. 35 (1980), 254-283. [69] Piran, T. "Numerical Relativity.", in Proceedings of the Second Marcel Grossmann Meeting on General Relativity. R. Ruffini, ed. New York: North-Holland, 1982. [70] Piran, T. "Methods of Numerical Relativity.", in Gravitational Radiation. N. Deruelle and T. Piran, eds. Amsterdam: North-Holland, 1983. 203-256. [71] Piran, T. and R. F. Stark. "Gravitational Radiation, Gravitational Collapse and Numerical Relativity." (Talk given by T. Piran at the XII Texas Sympo-sium, December 1984, Jerusalem, Israel.) (Preprint, 1985.) [72] Piran, T., P. N. Safier and R. F. Stark. "A General Numerical Solution of Cylindrical Gravitational Waves." (Preprint, 1985). [73] R i c h a r d s o n , L . F . " T h e A p p r o x i m a t e A r i t h m e t i c a l Solution by F i n i t e Differ-ences of P h y s i c a l P r o b l e m s i n v o l v i n g Differential E q u a t i o n s , w i t h a n A p p l i c a -t i o n to the Stresses in a M a s o n r y D a m . " Phil. Trans. Roy. Soc. 210 (1910), 307-357 [74] R i c h t m y e r , R . D . a n d K . W . M o r t o n . D i f f e r e n c e M e t h o d s f o r I n i t i a l -V a l u e P r o b l e m s . Second E d i t i o n . N e w Y o r k : W i l e y , 1957. [75] R o b e r t s , K . V. a n d N . O. Weiss. " C o n v e c t i v e Difference S c h e m e s . " M a t h . Comp. 20 (1966) 272-299. [76] S h a p i r o , S. L . a n d S. A . T e u k o l s k y . " G r a v i t a t i o n a l C o l l a p s e to N e u t r o n S t a r s a n d B l a c k Holes: C o m p u t e r G e n e r a t i o n of S p h e r i c a l S p a c e t i m e s . " Astrophys. J. 235 (1980), 199-215. [77] S h a p i r o , S. L . a n d S. A . T e u k o l s k y . "Relativist ic Stellar D y n a m i c s on the C o m p u t e r . I. M o t i v a t i o n a n d N u m e r i c a l M e t h o d . " Astrophys. J . 298 (1986), 35-57. [78] S m a r r , L . " T h e S t r u c t u r e of G e n e r a l R e l a t i v i t y w i t h a N u m e r i c a l I l l u s t r a t i o n : T h e C o l l i s i o n of T w o B l a c k H o l e s . " U n i v . of T e x a s at A u s t i n P h . D . thesis, 1975. [79] S m a r r , L . ed. S o u r c e s of G r a v i t a t i o n a l R a d i a t i o n . C a m b r i d g e : C a m b r i d g e U n i v e r s i t y P r e s s , 1979. [80] S m a r r , L . " B a s i c C o n c e p t s in F i n i t e Differencing of P a r t i a l Differential E q u a -t i o n s . " , i n S o u r c e s of G r a v i t a t i o n a l R a d i a t i o n . L . S m a r r , ed. C a m b r i d g e : C a m b r i d g e U n i v e r s i t y P r e s s , 1979. 139-159. [81] S m a r r , L . " G a u g e C o n d i t i o n s , R a d i a t i o n F o r m u l a e a n d the B l a c k Hole C o l l i -s i o n . " , in S o u r c e s o f G r a v i t a t i o n a l R a d i a t i o n . L . S m a r r , ed. C a m b r i d g e : C a m b r i d g e U n i v e r s i t y P r e s s , 1979. 245-274. [82] S m a r r , L . " S p a c e - t i m e s G e n e r a t e d by C o m p u t e r s : B l a c k Holes w i t h G r a v i t a -t i o n a l R a d i a t i o n . " A n n . N.Y. Acad. Sci. 302 (1977), 569-604. [83] S m a r r , L . " C o m p u t a t i o n a l R e l a t i v i t y : N u m e r i c a l a n d A l g e b r a i c A p p r o a c h e s . R e p o r t of W o r k s h o p A 4 " , in G e n e r a l R e l a t i v i t y a n d G r a v i t a t i o n , B . B e r t o t t i et. al . , eds. B o s t o n : D . R e i d e l , 1984. 163-183. [84] S m a r r , L . , A . C a d e z , B . D e W i t t , a n d K . E p p l e y . " C o l l i s i o n of T w o B l a c k H o l e s : T h e o r e t i c a l F r a m e w o r k . " Phys. Rev. D. 14 (1976), 2443-2452. [85] S m a r r , L . a n d J . W . Y o r k , J r . " R a d i a t i o n G a u g e in G e n e r a l R e l a t i v i t y . " Phys. Rev. D. 17 (1978), 1945-1956. [86] S m a r r , L . a n d J . W . Y o r k , J r . " K i n e m a t i c a l C o n d i t i o n s in the C o n s t r u c t i o n of S p a c e t i m e . " Phys. Rev. D. 17 (1978), 2529-2551. [87] S t a r k , R . F . a n d T . P i r a n . " G r a v i t a t i o n a l W a v e E m i s i o n f r o m R o t a t i n g G r a v -i t a t i o n a l C o l l a p s e . " ( p r e p r i n t , 1985.) [88] S t a r k , R . F . a n d T . P i r a n . " A N u m e r i c a l C o m p u t a t i o n of the G r a v i t a t i o n a l R a d i a t i o n f r o m R o t a t i n g G r a v i t a t i o n a l C o l l a p s e . " ( P a p e r d e l i v e r e d at the F o u r t h M a r c e l G r o s s m a n n M e e t i n g on G e n e r a l R e l a t i v i t y , J u n e 1985, i n R o m e , Italy; to be p u b l i s h e d in the P r o c e e d i n g s of the M e e t i n g . ) [89] S t e w a r t , J . M . " N u m e r i c a l R e l a t i v i t y . " , in C l a s s i c a l G e n e r a l R e l a t i v i t y . W . B . B o n n o r , J . N . I s h a m a n d M . A . H . M a c C a l l u m , eds. C a m b r i d g e : C a m -b r i d g e U n i v e r s i t y P r e s s , 1984. [90] T a y l o r , J . H . , L . A . F o w l e r a n d P. M . M c C u l l o c h . " M e a s u r e m e n t s of G e n e r a l R e l a t i v i s t i c Effects i n the B i n a r y P u l s a r P S R 1913 + 16." N a t u r e 277 (1979), 437-440. [91] T a y l o r , J . H . a n d J . M . W e i s b e r g " A N e w T e s t of G e n e r a l R e l a t i v i t y : G r a v i -t a t i o n a l R a d i a t i o n a n d t h e B i n a r y P u l s a r P S R 1913+16." Astrophys. J. 253 (1982), 908-920. [92] T h o r n b u r g , J . C o o r d i n a t e s a n d B o u n d a r y C o n d i t i o n s f o r t h e G e n e r a l R e l a t i v i s t i c I n i t i a l D a t a P r o b l e m U n i v . of B r i t i s h C o l u m b i a M . S c . thesis ( u n p u b l i s h e d ) , 1985. [93] T h o r n e , K . S. " T h e T h e o r y of G r a v i t a t i o n a l R a d i a t i o n : A n I n t r o d u c t o r y R e -v i e w . " , i n G r a v i t a t i o n a l R a d i a t i o n . N . D e r u e l l e a n d T . P i r a n , eds. A m s -t e r d a m : N o r t h - H o l l a n d , 1983. 1-57. [94] T s c h a r n u t e r a n d K . - H . W i n k l e r . " A M e t h o d for C o m p u t i n g S e l f g r a v i t a t i n g G a s F l o w s w i t h R a d i a t i o n . " Comp. Phys. Comm. 18 (1979), 171-199. [95] U n r u h , W . G . " N o t e s on B l a c k - H o l e E v a p o r a t i o n . " Phys. Rev. D. 14 (1976), 870-892. [96] U n r u h , W . G . P e r s o n a l c o m m u n i c a t i o n . [97] v a n L e e r , B . " T o w a r d s the U l t i m a t e C o n s e r v a t i v e Difference Scheme. II. M o n o -t o n i c i t y a n d C o n s e r v a t i o n C o m b i n e d in a S e c o n d - O r d e r S c h e m e . " , J. Comp. Phys. 14 (1974), 361-370. [98] v a n L e e r , B r a m . " T o w a r d s the U l t i m a t e C o n s e r v a t i v e Difference S c h e m e III. U p s t r e a m - C e n t e r e d F i n i t e - D i f f e r e n c e S c h e m e s for Ideal C o m p r e s s i b l e F l o w . " J. Comp. Phys. 23 (1977), 263-275. [99] W a l d , R . M . G e n e r a l R e l a t i v i t y . C h i c a g o : U n i v e r s i t y of C h i c a g o P r e s s , 1984. [100] W a r m i n g , R . F . a n d B . J . H y e t t . " T h e M o d i f i e d E q u a t i o n A p p r o a c h to the Sta-b i l i t y a n d A c c u r a c y A n a l y s i s of Finite-Dif ference M e t h o d s . " J. Comp. Phys. 14 (1974), 159-179. [101] W e i n b e r g , S t e v e n . G r a v i t a t i o n a n d C o s m o l o g y : P r i n c i p l e s a n d A p p l i -c a t i o n s of t h e G e n e r a l T h e o r y of R e l a t i v i t y . N e w Y o r k : J o h n W i l e y a n d S o n s , 1972. [102] W e i s b e r g , J . M . a n d J . H . T a y l o r . " O b s e r v a t i o n s of P o s t - N e w t o n i a n T i m i n g Effects in the B i n a r y P u l s a r P S R 1913+16." Phys. Rev. Lett. 52 (1984), 1348-1350. [103] W i l s o n , J . R . " N u m e r i c a l S t u d y of F l u i d F l o w s in a K e r r Space." Astrophys. J . 173 (1972), 431-438. [104] W i l s o n , J . R . " A N u m e r i c a l M e t h o d for R e l a t i v i s t i c H y d r o d y n a m i c s . " , in S o u r c e s of G r a v i t a t i o n a l R a d i a t i o n . L . S m a r r , ed. C a m b r i d g e : C a m b r i d g e U n i v e r s i t y Press, 1979. 423-44 5. [105] Y o r k , J . W . " K i n e m a t i c s a n d D y n a m i c s of G e n e r a l R e l a t i v i t y . " , in S o u r c e s of G r a v i t a t i o n a l R a d i a t i o n . L . S m a r r , ed. C a m b r i d g e : C a m b r i d g e U n i v e r s i t y Press, 1979. 83-126. [106] Y o r k , J . W . " T h e Initial V a l u e P r o b l e m a n d D y n a m i c s . " , in G r a v i t a t i o n a l R a d i a t i o n . N . Deruel le a n d T . P i r a n , eds. A m s t e r d a m : N o r t h - H o l l a n d , 1983. 175-201. [107] Y o r k , J . W . a n d T . P i r a n . " T h e Initial V a l u e P r o b l e m a n d B e y o n d . " , in S p a c e t i m e a n d G e o m e t r y : T h e A l f r e d S c h i l d L e c t u r e s . R. A . M a t z n e r a n d L . C. Shepley, eds. A u s t i n : U n i v e r s i t y of T e x a s P r e s s , 1982. 147-176. A P P E N D I X A B a s i c Difference Techniques for H y p e r b o l i c E q u a t i o n s This appendix reviews the basic ideas of finite differencing and stability anal-ysis for time dependent PDE's which are used in the body of the thesis. The discussion is limited to problems in one space and one time dimension and the simplest hyperbolic equation ut — aux = 0 is used for illustrative purposes. Most of this material can be found in one form or another in any of a number of text books and monographs. In particular Richtmeyer and Morton [72], Kreiss and Oliger [41], Kreiss [42] and Mitchell and Griffiths [53] are all very useful. The review by Smarr [78], which is directed specifically at relativists may also be of interest. Finite difference methods for the solution of PDE's are based on the familiar notion of replacing a continuum of function values by a countable set of values de-fined at discrete points which constitute a mesh or grid. The differential operators are then replaced by finite difference quotients, resulting in a set of algebraic equa-tions whose solution is an approximation of the solution of the continuum system. Whenever possible, it is convenient to work with a uniform grid, such as the one depicted in Fig. A . l , which is characterized by unique mesh spacings in each of the coordinate directions. In the case shown in the figure, the grid covers the half-plane — oo < x < oo, t > 0 and the mesh points are the set {(xj,tn)\xj — 228 jh,tn = nk) for integral j and n, with n > 0. For any function / , the following notation is understood: /*? = /(xy,tn). ri • • • • • k a • • • i i • • [_ • • • • • • • i i i [p • — h — • • • i i • E> •• E3- - a - - -x Figure A . l . Finite difference mesh (grid). For a given differential system, there are many possible finite difference ap-proximations and many ways of deriving them; no attempt is made here to sum-marize the available techniques. Instead, one particular scheme for the solution of (A.l) is examined in some detail to illustrate basic concepts and tools of analysis. Introducing the difference operators (as in Chap. 3) AQ, A 0 defined by xn+l _ fn—1 A D / ? = ' u j (A3) the leap-frog scheme for (A.l) is A 0 v ; - a A ^ = 0 [AA) where vj denotes the approximation to the exact value u™. For any given differ-ence scheme, the set of locations of the values referenced at a mesh point (ij, tn) constitute the difference star of the method at that point. The star of the leap-frog scheme is shown in Fig A.2. The diameter of a difference scheme will be defined as the maximum number of mesh points which is covered by the star in either of the coordinate directions, so the diameter of the leap-frog method is 3. It is also common to refer to a scheme with a temporal diameter of q + 1 as a q-level or q-step method. A difference scheme such as the current one, whose star contains only a single point at the most advanced time level can be solved explicitly for the unknown at that point, and is therefore called an explicit scheme. Otherwise, the unknowns at the advanced level will be coupled and the scheme is implicit. Finally, a g-level method requires the specification of initial values • • • v ' - 1 , . after which the solution may be repeatedly advanced in time using the difference equations. • • • • Figure A.2. Difference star for the leap-frog method. Consider a general differential system Lu(x,t) - 0 (A.5) where L is some differential operator, and a corresponding difference system Lk'kv? = 0 (A.6) where Lh'k is a finite difference operator. Then the local truncation error, r", of the difference scheme will be defined by = Lh'ku] (A.l) and provides a measure of how well the difference operator locally approximates the differential operator. The difference scheme is consistent if ||r?|| -> 0 as h,k^0 (A.S) where CO l l /yl l= £ Uh (^9) j=-oo Secondly, the difference scheme is 0(h°l + k*2) accurate if < c(nk)(hSl + k52) as h, k 0 (A.10) for some function c(t) which is bounded on any finite time interval. For the case of the leap-frog scheme for (A.l), it easy to calculate k2 h2 A 0 < - aAgu? - — uut - a—uxxx + 0(h4 + k4) (A.ll) J 1 6 6 by performing Taylor series expansions of u" + 1 , u" _ 1 , and u"_j about the center of the star. Thus provided the derivatives uttt and uxxx are bounded, the scheme is 0(h2 + k2) accurate. Finally, the notion of strong stability is often introduced to capture the notion that the norm of the solution of a difference scheme 231 should not grow "too rapidly" in time. The argument goes that nature permits, at most, exponential growth in any physical quantity, so a difference scheme is called strongly stable if there are constants k, a such that | | „ ; | | < fce"*||t,}|| (A.12) It can be quite difficult to establish this type of stability for complicated systems of difference equations. It is however, a useful notion since, at least for linear constant coefficient problems, it allows bounds to be placed on the solution error e" = u"-v" in terms of the truncation error. A theorem due to Lax [72] basically states that, for these problems, given consistency, strong stability implies convergence (e" —» 0 as h, k —• 0) and vice versa. In their landmark 1928 paper, Courant, Friedrichs and Levy (CFL) made the observation that the solution of a difference system for a hyperbolic equation cannot converge to the exact solution if the physical domain of dependence of a grid point is not contained in the numerical domain of dependence [17]. For explicit schemes, this places a restriction on how large k can be, for given o and h, before the method becomes unstable. The situation is illustrated for the model system (with a — l) in Fig. A.3. In both cases, the physical domain of dependence of the advanced grid point is indicated by the dashed line; the numerical domain is simply the rest of the computational star. The second case is unstable because the CFL condition k a = o - < 1 (A.13) h is violated. A CFL condition is not usually considered much of a restriction since, for hyperbolic problems, accuracy considerations suggest taking k w h. The general theory of stability for difference approximations of hyperbolic problems (including mixed initial/boundary value problems) is well-developed [32,40], but the task of applying this theory to determine whether a particular scheme is 232 \ \ • \ • > • • \ a < 1 stable cr > 1 unstable Figure A . 3 . Domains of dependence—the CFL condition. stable can become very difficult. In fact, many different notions of stability, and different techniques of analysis have been introduced for the purposes of examining the properties of difference methods. In practice, the technique due to von Neu-mann of looking at the properties of the difference scheme in the frequency domain continues to be the most useful way of determining practical stability conditions for hyperbolic problems. The basic idea is to write the method being analyzed in the form \{x,t + k) = Qv(x,*) where v is, in general, a vector of unknowns (multi-step schemes can always be written in this form by the introduction of auxiliary variables) and Q is the finite difference transition matrix. Assuming Q does not depend on t or x, in Fourier space (A. 14) becomes vB + 1(w) = G (0 vn{u) (A.15) for each t, where +oo e-'wx v(i) dx (A.16) Here, £ = u>h and G ( £ ) is called the amplification matrix. The von Neumann condition demands that all of the eigenvalues, fij, of G ( f ) satisfy < l + 0{k), |e| < Tr (A.17) It is clear that this is a necessary condition for stability; if it is violated for some values of £, then the corresponding Fourier components will grow more than expo-nentially as the solution proceeds. For problems where it is known that there is no exponential growth in the exact solution, the condition may be strengthened to i < i , lei < i r (A.IS) and it is this form which is used in Chap. 3. Consider the application of this technique to the model scheme the auxiliary variable = r™, (A.4) may be written as where, again, o — ak/h. Under Fourier transformation OAQ -* 2tosin $ (A.20) so the amplification matrix is G(0=(2''Tf J) (A-21) with eigenvalues given by the roots of the characteristic polynomial . Introducing (A.19) p{fi) — M2 — 2tcT sin £u — 1 (A.21) T h e s e r o o t s are t o s i n 1 - a2 s i n 2 t: (A.22) T h u s , c o n d i t i o n (A.18) is satisfied p r o v i d e d t h a t a < 1, w h i c h is just the C F L c o n -d i t i o n . A l s o , i n this case, the eigenvalues of the amplif icat ion m a t r i x have precisely u n i t m o d u l u s for a l l £; a scheme w i t h this p r o p e r t y is said to be non-dissipative since there will be no a t t e n u a t i o n of the F o u r i e r c o m p o n e n t s as the solution proceeds. In p r a c t i c e , the p r o b l e m of d e t e r m i n i n g w h e t h e r the roots of the characteristic poly-n o m i a l lie w i t h i n the unit circle in the c o m p l e x p l a n e c a n b e c o m e q u i t e difficult, p a r t i c u l a r l y for m u l t i - s t e p m e t h o d s or w h e n there are several u n k n o w n s i n v o l v e d . In this r e g a r d , the t h e o r e m s d u e to M i l l e r [52] (see also H e n r i c i [39]) c a n be quite useful since t h e y p r o v i d e a way of r e d u c i n g the task of l o c a t i n g the r o o t s o f an n - t h degree p o l y n o m i a l to the p r o b l e m of f i n d i n g the r o o t s of a p o l y n o m i a l of degree n — 1. U n f o r t u n a t e l y , a l t h o u g h this gives a definite a l g o r i t h m for d e t e r m i n i n g the neces-s a r y c o n d i t i o n s for h a v i n g \fij\ < 1, the p r o c e d u r e c a n still become algebraical ly c o m p l i c a t e d . S t r i c t l y s p e a k i n g , von N e u m a n n s t a b i l i t y analysis is only a p p l i c a b l e to l inear, p u r e init ia l value p r o b l e m s w i t h c o n s t a n t coefficients. H o w e v e r , n u m e r i c a l experi -ence has s h o w n t h a t it c a n also be usefully a p p l i e d to m o r e general p r o b l e m s by l i n e a r i z i n g the s y s t e m (if necessary) a n d freezing coefficients prior to F o u r i e r trans-f o r m i n g . T h e v o n N e u m a n n c o n d i t i o n is t h e n establ ished for all conceivable values o f t h e coefficients to d e t e r m i n e the p r a c t i c a l stabi l i ty criteria. C o n s i d e r the solution of ( A . l ) w h e n the init ial d a t a is a single F o u r i e r c o m -p o n e n t w i t h w a v e n u m b e r q u(i, 0) = e tqx (A.23) T h e n u(x,t) = e i q [ x + a t ) = e * " ( « x + w ' ) (A.24) which is a wave moving to the left with phase velocity a. Now for this simple equation, (but also for scalar, electromagnetic or (weak) gravitational radiation propagating in vacuum), the phase velocity is independent of q so there will be no dispersion of a general wave packet. However, this property is usually not retained in difference solutions to (A.l), or other non-dispersive hyperbolic systems. For example, consider the leap-frog semi-discretization of (A.l) vt — OAQV (A.25) and look for solutions of the form v = e*rtx+B,t> (A.26) Then one finds , sin qh , a = —f- (A.27) qh so the difference scheme is dispersive and in fact a' —+ 0 at the highest frequencies qh —• 7T. Similarly, it can be shown that the group velocity of the difference solution is c — a cosqh (A.28) which tends to 0 even sooner. This means that there are severe problems with this scheme in the numerical transport of disturbances having wavelengths comparable to the mesh spacing. This problem is inherent to most finite difference approxi-mations, although there has been much effort directed towards the development of methods with particularly low levels of dispersion [34].
Thesis/Dissertation
eng
Physics
Vancouver : University of British Columbia Library
For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use.
Graduate
A study of numerical techniques for radiative problems in general relativity
Text