Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Relativistic hydrodynamics and other topics in numerical relativity Olabarrieta, Ignacio (Inaki) 2004

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

Item Metadata


831-ubc_2004-99435x.pdf [ 8.71MB ]
JSON: 831-1.0085722.json
JSON-LD: 831-1.0085722-ld.json
RDF/XML (Pretty): 831-1.0085722-rdf.xml
RDF/JSON: 831-1.0085722-rdf.json
Turtle: 831-1.0085722-turtle.txt
N-Triples: 831-1.0085722-rdf-ntriples.txt
Original Record: 831-1.0085722-source.json
Full Text

Full Text

Relativistic Hydrodynamics and other topics in Numerical Relativity by Ignacio (Ifiaki) Olabarrieta Licenciado en Ciencias, Universidad del Pais Vasco, 1998 M . S c , The University of British Columbia, 2000 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 M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F D O C T O R O F P H I L O S O P H Y T H E F A C U L T Y O F G R A D U A T E S T U D I E S • Physics . : • T H E U N I V E R S I T Y O F B R I T I S H C O L U M B I A June 8, 2004 © Ignacio (Ifiaki) Olabarrieta, 2004 Relativistic Hydrodynamics and other topics in Numerical Relativity by Ignacio (Inaki) Olabarrieta Licenciado en Ciencias, Universidad del Pais Vasco, 1998 M . S c , The University of British Columbia, 2000 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 M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F 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 O F G R A D U A T E STUDIES (Department of Physics and Astronomy) We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y O F B R I T I S H C O L U M B I A June 8, 2004 © Ignacio (Ifiaki) Olabarrieta, 2004 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Physics and Astronomy The University Of British Columbia Vancouver, Canada Date . Abstract i i Abstract In this'' thesis I consider three different projects in numerical relativity. The first one is a study of the spherically-symmetric collapse of a scalar field with a potential that mimics the inclusion of angular momentum. This work has been carried out in collaboration with M . W . Choptuik, W . Unruh and J . Ventrella. In this study we found a new family of type II critical solutions which are discretely self similar. The second project involves work I did in another collaboration with M . W. Choptuik, L . Lehner, R. Petryk, F . Pretorius and H . Villegas. Here we study the dynamical evolution of 5-dimensional generalizations of black holes, called black strings, which are known to be unstable to sufficiently long-wavelength perturbations along the string direction. Not only have we been able to dynami-cally trigger the instability, explicitly verifying the results from perturbation theory, we have been able to evolve for sufficiently long times to observe that the system goes through a phase (not necessarily the final end-state) that resembles a series of black holes connected by a thin black string. The third and most extensive part of this thesis is a study of ideal fluids fully coupled to gravity, both in spherical symmetry and in axisymmetry. In this project we have cast both the dynamic and equilibrium equations for general relativistic hydrodynamics in the 2+1+1 formalism and in a way that is tailor-made for the use of high resolution shock capturing methods. In addition, our implementation, for the case of no rotation, is able to evolve discontinuous data and has proven to be convergent. Unfortunately our implementation currently has too much numerical dissipation, and suggests that the use of adaptive methods may be very helpful in achieving long term evolution of star-like configurations. Contents i i i Contents Abstract 1 1 Contents 1 1 1 List of Figures V 1 List of Tables v i i i Acknowledgements 1 X 1 Introduction 1 1.1 3+1 Decomposition 3 1.2 Critical Phenomena 7 1.3 Black Holes and Black Strings 10 1.4 Relativistic Hydrodynamics 13 1.5 Summary of Results 16 2 Scalar Field Collapse with Angular Momentum 18 2.1 Introduction 18 2.2 Equations of Motion 19 2.2.1 Equations 19 2.2.2 Regularity and Boundary Conditions 22 2.3 Results ' 24 2.3.1 Numerics 24 2.3.2 Families of Initial Data 25 2.3.3 Analysis 26 2.3.4 Results 27 2.4 Conclusions 32 Contents 3 Instability of a Black String 38 3.1 Introduction 38 3.2 Equations 39 3.3 Numerical Implementation 41 3.3.1 Evolution 41 3.3.2 Determination of initial data: solving the constraints 43 3.3.3 Finding Event Horizons 45 3.4 Results 47 4 General Relativistic Hydrodynamics in Spherical Symmetry 52 4.1 Introduction 52 4.2 Model/Equations 53 4.3 Numerics 58 4.3.1 Geometric Equations 59 4.3.2 Hydrodynamic Equations 61 4.4 Results 69 4.4.1 New Variables 69 4.4.2 Evolution of T O V solution 72 5 Axisymmetric Hydrodynamics 78 5.1 Introduction 78 5.2 2+1 + 1 Formalism 78 5.2.1 Geometry 81 5.2.2 Fluids 85 5.2.3 Boundary and Regularity conditions 87 5.3 Initial Data for Rotating Stars 88 5.3.1 Integrability Condition 90 5.4 Numerics 92 5.5 Results 97 6 Conclusions and Directions for Future Research 108 6.1 Scalar Field Collapse with Angular Momentum 109 6.2 Instability of a Black String 109 6.3 Relativistic Hydrodynamics 110 Bibliography 112 Contents v A Hydrodynamic equations in conservation form 121 B Scalar fields considered in the 1=1 case 123 C Hydrodynamic equations in the 2+1+1 formalism 127 D Characteristic Structure 136 List of Figures vi List of Figures 1.1 3+1 decomposition of spacetime 6 1.2 Penrose diagram of a spherical symmetric black hole . 11 2.1 Values of log 1 0 (A;) versus 1 28 2.2 Values of log 1 0 (7 ; ) versus 1 29 2.3 Value of max( {maxr {2M(t, r)/r}} in the critical regime as a function of £ 30 2.4 Values of T* versus 1 31 2.5 T*Ai as a function oi I 33 2.6 Fit to periodicity of the times at which maxr {2M/r(t, r)} attains its maximum. . . 34 2.7 Fit to DSS of the times at which max r {2M(t,r)/r} attains its maximum 35 3.1 Values of R m a . x , Rmm and A for the apparent horizon as a function of time 48 3.2 Embedding diagram of the evolution of the apparent horizon 49 3.3 Plots of the apparent horizon and estimates of the event horizon location . 51 4.1 Outgoing bounday conditions for the fluid variables 58 4.2 Control Volume for the spherical symmetric code 62 4.3 Calculation of the Roe flux 66 4.4 ip{°°) calculated using the original variables and a shooting method to integrate the HC 70 4.5 Solution curves for TOV initial conditions 74 4.6 Convergence of the evolution of a TOV solution 76 4.7 Evolution of a TOV solution with varying amounts of numerical flux. . . 77 5.1 Numerical grid for axisymmetric hydrodynamics 96 5.2 Convergence test for a weak gravitating axisymmetric pulse of perfect fluid 98 5.3 Geometry of the setup to study the K H instability 100 5.4 K H instability 101 5.5 Pressure evolution of a weakly self-gravitating pulse of fluid 103 5.6 Evolution of the conformal factor for a weakly self-gravitating pulse of fluid 104 List of Figures 5.7 Evolution of a strongly self-gravitating pulse of fluid. .- 105 List of Tables vi i i List o f Tables 2.1 Families of initial Data 25 2.2 Type of initial data •. . . 25 2.3 Properties Of critical solutions for different values of 1 37 Acknowledgements ix Acknowledgements I want to thank the Government of the Basque Country for providing me with a grant to work on this thesis. I also want to thank my collaborators in the projects discussed and all the current and past members of the numerical relativity group in U B C . Special thanks go to Scott Noble for teaching me so much about relativistic hydrodynamics. Tom Davis, Ph i l Eles, and Bojan Losic for listening to some of my paranoias. M y family and friends in Spain who always have encouraged me to finish this thesis. And finally my supervisor, Matt Choptuik, for all his help, encouragement and patience. Chapter 1. Introduction 1 Chapter 1 Introduction Strong, dynamical gravity plays a crucial role in many interesting astrophysical systems. Especially when the strength of gravity becomes such that the typical length scale, R, of the system is comparable to the gravitational length scale, LQ = GM/c2, denned by the mass, M , of the system, (G and c are, respectively, Newton's constant, and the speed of light in vacuum) or typical velocities are comparable to c, Newtonian gravity fails to completely describe the dynamics of the system. In such cases Einstein's theory of gravitation—general relativity—gives a more accurate description. Some examples of these astrophysical systems are supernovae explosions, neutron stars (including pulsars), and black holes. In Einstein's theory of gravitation the effects of gravity are encapsulated in a Lorentzian metric tensor, gap, denned on a differentiable manifold. This manifold is identified with the region of the universe (i.e. space and time) in which we are interested. This metric tensor tells us how distances and times change from point to point in the manifold and this translates directly to how one measures distances and times in the physical spacetime. The theory provides equations that the metric tensor satisfies, known as the Einstein field equations: Ga0 = —^-Tafj. (1.1) Gap is known as the Einstein tensor and, in general, is a very complicated function of the metric gap as well as its first and second spatial and temporal derivatives. Specifically Ga0 = Ra0 — -ydapR, (1.2) where Rap is the Ricci tensor and R = Tr[Rap] is the Ricci scalar (Tr denotes the trace). On the right hand side of equation (1.1) we have the stress-energy tensor, Tap, which describes the matter content in our spacetime. The tensor form of the Einstein equations hide their great complexity when written out explic-itly in a specific coordinate system. In general they are 10 coupled, non-linear, time-dependent partial differential equations, which for most cases, particularly those of astrophysical interest, are amenable only to numerical solution. Additionally, these equations are invariant under general diffeomorphisms. This invariance is associated with the coordinate invariance or covariance of the Chapter 1. Introduction 2 physical theory. In general we need to fix the coordinate system in some way. This frequently involves solving some extra partial differential equations. The purpose of this thesis is to numerically evolve these equations for three different situations, focusing on the case in which Tap is the stress-energy tensor associated with a perfect fluid so that we must solve the equations of general relativistic hydrodynamics along with Einstein's equations. Spacetimes containing scalar fields as well as vacuum solutions are also considered. In particular in Chap. 2, we discuss a new family of critical solutions that has been discovered in collaboration with M . Choptuik, W . Unruh and J . Ventrella. These results were found for a massless scalar field in spherical symmetry with a potential that tries to mimic the effects of angular momentum. I played a leading role in almost all phases of this project. In Chap. 3, we discuss the main results of a collaboration with M . Choptuik, L . Lehner, R. Petryk, F . Pretorius and H . Villegas in which the instability of a black string, a five dimensional extension to a black hole, was studied dynamically [18]. I focus on the parts of the computations for which I was responsible; this involved solving the initial data constraints and developing some tools for analysis of the physics of the solutions. In the last part of the thesis, Chaps. 4 and 5, we discuss our efforts to study the evolution of a perfect fluid in axisymmetric spacetimes. The equations for<the fluid and the geometry in the so-called 2+1+1 formalism are presented and our numerical implementation of a computer code to solve the system is explained in detail. The remainder of this introductory chapter is organized as follows. First, in Sec. 1.1, we discuss one way to cast equations (1.1) in a way that is appropriate for their solution as an initial-value (dynamical) problem. This formalism is the one used in Chaps. 2, 3, and 4, and it is a crucial part of the formalism used in Chap. 5. Then, in Sec. 1.2, we briefly introduce the topic of critical phenomena in gravitational collapse. A short introduction to black holes and black strings in 5 dimensions follows. In Sec. 1.4, an introduction to, and a summary of, previous work in relativistic hydrodynamics is given. We conclude with a brief summary of the results of the thesis in Sec. 1.5. In the remainder of the thesis we choose units such that G = c = 1, we adopt the Einstein summation convention for tensor-component indices, and we use the signature {—, +, +, +} for the metric tensor gap. Greek/lower-case latin/upper-case latin indices are 4-dimensional (spacetime), 3-dimensional (spatial) and 2-dimensional (spatial) tensor-component indices, respectively. Indices involving letters from the beginning of each alphabet, i.e. a, (3, 7 , . . . , a, b, c,... and A, B, C are generally abstract tensor indices, while those involving letters near the middle of each alphabet, i.e. fi, v, i, j, k,... and / , J, K,... generally refer to tensor components taken with respect to some specific coordinate system. Chapter 1. Introduction 3 1.1 3+1 Decomposition In this section we discuss a way of rewriting equations (1.1) in a form tailor-made for the study of their evolution. There are many people who have worked on the problem of decomposing these equations: [3], [10], [71] and [115] describe some of the main efforts. Here we follow York's development [115], which gives a "3+1" decomposition of the field equations (1.1), and a set of first-order-in-time equations for dynamical variables describing the geometry. Our goals are to choose the values of the free fields denned at an instant of time (constant-* spacelike hypersurface), and to find the solution of the geometry and the matter for the whole subsequent spacetime (future or past development of the initial data). In this decomposition of equations (1.1), it is also possible to identify which geometric fields are to be freely chosen at the initial time—this is a non-trivial issue due to the existence of the constraints (described below). Since we want to set initial data at a constant-* hypersurface, it is natural to foliate the space-time with spacelike hypersurfaces of constant time. This foliation will produce a decomposition of equations (1.1) into equations constraining dynamical quantities intrinsic to each hypersurface, and equations that tell us how to propagate these dynamical variables from hypersurface to hyper-surface. One way of defining the foliation is by introducing a scalar field t with the following property: the level sets defined by this scalar field, i.e. surfaces of constant t, must be spacelike hypersurfaces, meaning that any vector field defined on the hypersurfaces is spacelike. For now, we will assume that these hypersurfaces cover all-of the spacetime (or at least the region of interest) and that they are free of (physical) singularities. Associated with t is the (locally) closed one form, fla, given by fta = V a i . (1.3) Since the slices are spacelike, fla is timelike, and we can therefore write the square of its norm as gaf3nanp = - \ , (i.4) where a is a positive function 1. Notice that the distance (actually the proper time) between hypersurfaces labelled by t and t + dt, as measured by observers moving orthogonally to the slices, is adt. For this reason, a, is usually called the lapse function. It is useful to define the normalized one-form, u>a uja = a f i Q , (1.5) and the unit normal vector, na na = (1.6) 1 Those not familiar w i th tensor analysis should not confuse the a function wi th the a index. Chapter 1. Introduction 4 where the sign has been chosen so that na is future-directed. Using this definition the metric, -yap, induced in the spacelike hypersurface can be written as lap = 9ap + nanp. (1.7) 7 a | a is a purely spatial tensor, i.e. it has no component along na: lap™13 = ga(3np + n a n 0 n p = n a - n a = 0, (1.8) where we have used the fact that n a n a = — 1. It is also important to note that the mixed form l a p = Sap + n a n p , (1.9) is the projection operator onto the spatial hypersurfaces. Associated with the Riemannian met-ric (1.7) there is a 3-metric-compatible covariant derivative, Da, satisfying D a l b c = 0 (1.10) that can also be defined via projection of the spacetime covariant derivative, V a . For example for a spatial vector, Va, we have D5V = rs^pVaV3. (1.11) Using D a , we can calculate the 3-dimensional curvature ^Rabcd, 3-dimensional Ricci tensor ^ R a b , and 3-dimensional Ricci scalar ^R. These quantities describe the internal geometry of each hyper-surface; in addition, we need a description of how the metric 7 t tb changes from slice to slice. This information is encoded in the (symmetric) extrinsic curvature tensor, Kab, which can be defined by: Kab = -^n"lab- (1.12) Here £„« is the Lie derivative along the unit vector field normal to the hypersurface and the factor of —1/2 is a matter of convention. It is worth noticing that Kan is a purely spatial tensor, since Cna commutes with the projection operator (1.9). Now we have all the elements needed to proceed with the 3+1 decomposition of (1.1). Specifi-cally, we consider the following projections of the field equations: n^n^G^ = 8 7 m a r / T Q / 3 , (1.13) l a a n p G a P = 8 n 7 a a n p T a P , (1.14) 7 a a l P b G a f 3 = tot-farfbTaf}. (1.15)' It is relatively straightforward to show that n a v ? G a f i = 1 [WR + {K\)2 - K a b K a b ] , (1.16) Chapter 1. Introduction 5 and 7 ° Q n / 3 G Q ' 3 = Db (Kab - jabKcc) . (1.17) We can now write equations (1.13) and (1.14) as i [WR + (KCC)2 -KabKab] = 8wp, (1.18) Da (Kab - jabKcc) = 8njb, (1.19) where we have made use of the following definitions for the projections of the stress energy tensor: p = nanpTaP, (1.20) f = -y^npT"*. (1.21) Physically, these quantities are the energy and momentum densities, respectively, that time-like observers moving normally to the hypersurfaces would observe. Note that equations (1.18) and (1.19) only involve spatial tensors and spatial derivatives of such tensors. These equations need to be satisfied by the 3-metric, ^ a b , the extrinsic curvature, Kai>, and the matter densities, p and ja, on every hypersurface. In particular, they need to be satisfied at the initial time and thus they constrain the initial data. This implies that not every metric function ^ a b and every component of Kab (12 values in total) are freely specifiable, but that they need to be chosen in such a way that the constraint equations are obeyed. Equation (1.18) is called the Hamiltonian constraint, while (1.19) are usually called the momentum constraints. To develop the other 6 field equations, equations (1.15), it is useful to use a time derivative operator, £jv=, where Na is a vector field dual to the one form fla, meaning that QaNa = 1. It is important to note that this vector field is unique only up to reparametrization of the (spatial) coordinates in each hypersurface. In particular, any vector field ta = (d/dt)a — Na + f3a = ana + (3a, such that (3ana = 0, satisfies £lata = 1 (see Fig. 1.1). We will hereafter assume that we have adopted a so-called "3+1" coordinate system adapted to some such vector field, ta, as illustrated in Fig . 1.1. Since /3Q is a spatial vector, we have f3° = 0 in. such a coordinate system. The vector field f3a encodes the diffeomorphism invariance within the hypersurfaces and, more geometrically, describes the shifting in the spatial coordinates from hypersurface to hypersurface relative to normal propagation (see Fig. 1.1). Hence (3a is known as the shift vector. In terms of £TV ° , the extrinsic curvature tensor can be written as: Kab = --Fr^Nalab = ~ 7. (A°7a6 ~ £/3«7ab) > (1.22) which implies Kab = ~ - A A - Da(3b ) . (1.23) Chapter 1. Introduction 6 t Figure 1.1: 3 + 1 decomposition of the spacetime. The distance (proper time) along the direction of the unit normal vector, n M , between hypersurfaces labelled by t and t + dt is given by adt. Note that 01 and dx% are purely spatial vectors defined on the tangent space of the hypersurface at event (t, x1). In general, in order to move to an event with the same spatial coordinates, x%, on hypersurface t + dt, i.e. in order to move along the vector field (d/dt)^, we have to shift by an amount (3ldt on the future spatial hypersurface relative to normal propagation. The distance between events with coordinates (t, xz) and (t + dt,xl + dx1) is ds and can be calculated by the "Pythagorean theorem" for Lorentzian geometry. Chapter 1. Introduction 7 After some algebra, we can write equations (1.15) as CwKob = C0«Kab - DaDba + a { ( 3 ) iU - 2KacKcb + KabKcc - 8TT Sab - )pab (Scc + p) J , (1.24) where S a b = ^al^bTap, known as the spatial stress tensor, is the total spatial projection of the stress energy tensor onto the hypersurface. Equations (1.18), (1.19) and (1-24) complete the decomposition of equations (1.1). Let us stress that we have: (a) kinematical variables, the lapse and the shift, that encode the coordinate freedom of the theory and that must be specified in some fashion; (b) dynamical vari-ables, jab, Kab, which, loosely speaking, encode the physical information describing the geometry of spacetime. The 3+1 decomposition of the 4-dimensional line-element can be expressed using the 4-dimensional Lorentzian version of the Pythagorean theorem (as described in [71]). Specifi-cally, the square of the spacetime displacement, ds2, between events with coordinates (t,xl) and (t + dt, x% + dx1) can be written as (see Fig. 1.1) ds2 = -(adt)2 + 7y(da: i + f3ldt)(dxj + j3jdt) = -(a2+ PiF)dt2+ 27ijPidtdxj+7ijdxidx:>. (1.25) A n important note is that since only a maximum of 6 of the 10 (second-order) Einstein equations are independent of each other, we have considerable freedom in choosing which specific equations to solve during a dynamical evolution. The initial data for the geometry and matter must always be chosen to satisfy the constraint equations (1.18-1.19), but the evolution equations (1.24), together with the contracted Bianchi identities, Gap',a = 0, then guarantee that the constraints are satisfied at any future or past time. This property of the continuum system of equations is in general lost once the system is discretized; at best the numerical solution only satisfies the unused equations at the truncation order of the approximation. Checking the convergence properties of the residuals of these unused equations thus gives us a very strong check of the validity of our solutions. 1.2 Critical Phenomena The discovery of critical phenomena in gravitational collapse was reported by Choptuik in [17]. He studied the spherically-symmetric collapse to black holes2 of a massless scalar field, </>, minimally 2 B l a c k holes are described in more detail in section 1.3, at this point they can be viewed as regions of spacetime in which even light is trapped and cannot escape to infinity. Chapter 1. Introduction 8 coupled to gravity. The Lagrangian for this problem is given by L = — -4>-a<t>' a, (1-26) 16TT 2v'aV ' v 1 where R is the Ricci scalar (the usual Lagrangian for gravity). At that time, it was already known that, under the same symmetry assumptions, the dynamical evolution of weak configurations of scalar field leads to dispersal of the scalar field to infinity, resulting in Minkowski-like spacetimes [24], and that configurations of sufficiently strongly self-gravitating scalar field give rise to black hole formation [25]. Choptuik's work was the first systematic study of the intermediate regime. In order to study the transition to black hole formation, the initial data were parametrized by a single parameter, p, such that for high parameter values black holes were formed, while for low values, dispersal was the end state. A n example of such a parameter is the initial amplitude of the scalar-field. A t the threshold of black hole formation, p = p*, interesting and unexpected results were found. First of all, as we approach p* from above, the mass of the black hole that is obtained tends to zero following a power law: MBH = C(p-p*)\ (1.27) where M B H is the mass of the black hole, C is a constant of proportionality and 7 « 0.37 is a universal exponent independent of the specific family of initial data used in the calculation. We note that scaling of dimensionful quantities is also seen if the critical solution is approached from below (subcritical) rather that from above. Moreover, for values of p close to p*, the solutions in the strong field regime approach—at least for some finite time and in some finite region of space—a solution called the critical solution, which is universal in the following sense. There exist coordinates (r, p) in which the critical solutions obtained using different families of initial data all take the same form. In [17] so-called polar-areal coordinates, leading to a line-element ds2 = -a(t, rfdt2 + a(t, rfdr2 + r2dQ2, (1.28) r were used to perform the evolutions, with t measuring central proper time. The relationship between these coordinates and (r, p) is given by: r = ln(t*-t) + K, (1.29) p = l n r + K , (1.30) with t* and K depending on the family of initial data used to construct the solution. Not only was the critical solution found to be universal, it also showed an echoing behaviour. This means that in Chapter 1. Introduction 9 the critical regime, any dimensionless function, Z, constructed from the solution has the following symmetry: Z(r-A,p-A) = Z(r,p), • (1.31) where, for the scalar field, the period A « 3.4 is obviously also a universal parameter since all the features of the solution are universal when written in (r, p) coordinates. The transformation given by (1.31) amounts to a discrete rescaling of the original coordinates (t,r), and solutions with this rescaling invariance are thus said to be discretely self similar (DSS). The critical solution is unstable, essentially by construction: small deviations from the precise value p = p* lead to complete dispersal of the field, or to the formation of a black hole. Soon after the publication of [17], similar behaviour was observed in the collapse of axisymmetric gravitational waves [1] and, crucially, in spherically symmetric collapse of perfect fluid with a so-called radiation equation of state [29]. In this last work Evans conjectured that perturbation theory about the continuously self-similar solution that appeared at the black hole threshold for fluid collapse might provide insight into some of the phenomenology that had been observed. Work along precisely these lines was carried out with great success by Koike et al [53] who made the key observation that the "sharpness" of the empirically measured mass scaling law (1.27) for fluid collapse (as well as for the original scalar collapse), strongly suggested that the critical solution, had only a single unstable mode. From this assumption, and using a mathematical devel-opment precisely analogous to that used in the standard treatment of statistical mechanical critical phenomena, it follows that the mass-scaling exponent is simply the reciprocal of the Lyapunov exponent associated with the unstable mode. In follow up work, Gundlach [40] was able to show that this picture also held for the more difficult to treat case of the massless scalar field. Similar behaviour has been found in the collapse to black hole formation for many other types of matter. In general, different matter models behave in two qualitatively different ways at the threshold of black hole formation. In contrast to the massless scalar case, it has been found that the smallest black hole formed for some matter models is non-zero. In analogy with statistical mechanical critical phenomena (first and second order phase transitions) this sort of transition is dubbed "type I", whereas a transition characterized by (1.27) is called "type II". Type I critical solutions exhibit a fundamental length scale; it is this scale that determines the minimum value for the mass of the black hole above threshold. On the other hand, type II critical solutions are generically scale invariant. More specifically, type I critical solutions are usually static or periodic, while type II are self similar, or discretely self similar. Although many matter models have been investigated in spherical symmetry (see [43] for a recent review) very little is known about critical phenomena in less symmetric systems. So far, Chapter 1. Introduction 10 only two systems have been investigated: axisymmetric vacuum, studying the collapse to black hole formation of gravitational waves [1]; and the collapse of a scalar field in axisymmetry [20], [21]. The calculations described in Chap. 2 represent an attempt to investigate the possible effects of angular momentum on critical collapse via an effective potential for a scalar field in spherical symmetry. 1.3 Black Holes and Black Strings The Einstein field equations (1.1) allow for black hole solutions. For example in spherical symmetry the well-known Schwarzschild metric ds2 = g^doT'dx" = - ( l - dt2 + ( l - 2-^pj 1 dr2 + r W , (1.32) is the unique vacuum solution (apart from flat spacetime, which is not a regular limit of Schwarzschild as M —> 0) of the field equations. The region of the spacetime with r < 2M cannot communicate with infinity, i.e. no physical trajectory (timelike or null curve) originating from r < 2M can prop-agate to r —> oo (nor to r > 2M for that matter), and the region r < 2M is thus known as a black hole. In order to describe black holes, it helps to draw diagrams of the corresponding spacetimes that preserve causal structure and which compactify the spacetimes, so that various types of infinity lie at finite coordinate distance. These sorts of pictures are called Penrose diagrams. Figure 1.2 shows one such diagram for the Schwarzschild solution (1.32). Note that each point in the bulk of the diagram represents a 2-sphere manifold of specific radius and that light rays (null geodesies) travel along straight lines inclined at 45°. This diagram has many features worthy of discussion. Starting with the boundary, at the right-most edge of the plot there is a point called i° which represents spacelike infinity and which is the locus where all outgoing spacelike geodesies end. In addition, the point on top, labelled by i+, corresponds to future timelike infinity and is where future-directed timelike geodesies in region I terminate (i~ corresponds to past timelike infinity and is where past-directed timelike geodesies end). The line joining i+ with i° is future null infinity, denoted by I + , and is the place where null, future-directed, geodesies originating in region I terminate (likewise, past null infinity, T~, is the line connecting i° with If a timelike geodesic is within bulk region II it does not end at i+ but instead ends at the dotted line labelled r = 0 (the same happens for null geodesies within region II, they also cannot reach The dotted line corresponds to r — 0 in (1.32) and is, in fact, a physical (crushing) singularity of the spacetime, where the spacetime curvature goes to infinity. The observation that, starting in region II, no timelike or null geodesic can propagate into region I suggests that region II is a black hole. Its boundary is called an event horizon. Chapter 1. Introduction 11 . O r = 0 i Figure 1.2: Penrose diagram of a spherically symmetric black hole. Note that light rays travel at 45 degree angles in this diagram. Region II (shaded) is causally disconnected from I+ and thus lies within a black hole. The boundary of the black hole is called the event horizon. Plotted lines running approximately vertically (horizontally) correspond, in region I, to lines with constant coordinate r (t). We also note that the diagram shown is only half of the full Penrose diagram for the maximally extended Schwarzschild solution, see [71]. Chapter 1. Introduction 12 There are a whole family of spacetimes (Kerr-Newman solutions) with similar characteristics to Schwarzschild, and that can be characterized by just a- few parameters (specifically three —mass, angular momentum and charge). These spacetimes are believed to be the end-states of the evolution of sufficiently massive stars. A t the end of the evolution of a massive star (larger than about 10 M Q ) its structure resembles an onion, with layers of different materials with increasing atomic number towards its centre. The inner most core is composed of iron which is maximally bound with respect to the nuclear force, and hence cannot undergo thermonuclear burning. This core grows from the burning of the outer layers of silicon and sulphur, and when its mass becomes of the order of 1.5 M Q , its pressure cannot counteract the gravitational pull and it collapses. This process releases titanic amounts of energy that are carried out to the outer regions of the star mainly by neutrinos. A highly energetic shock is produced that moves outwards through the envelope. This shock energizes the outer region which undergoes a violent explosion—a supernova—which expels all of the outer layers of the original star. The remnant of this process is a very compact object of a few kilometres radius and, at most, several solar masses. In particular, if the mass of the remnant exceeds about 3 M Q 3 , even degeneracy pressure cannot hold the star up against gravitational collapse [89] and a black hole, whose mass will again be a few M Q , will form. Although there is not yet unequivocal observational evidence for this type of black hole, many candidates have been found in binary systems—for an inventory see [82]. There is also evidence for more massive black holes. Studies of microlensing events and "ultra-luminous" X-ray sources, among others, have provided some hints about the possible existence of black holes with intermediate masses [105], i.e. those with masses between a few solar masses and a few million. Moreover, there are strong indications that super massive black holes, with masses in the range 10 6 — 1 0 9 M Q or greater, can be found in most galactic centres. See [26] and [35] for discussion of the evidence for such a black hole at the centre of the Mi lky Way galaxy. It is a very well known result that black holes in 4 dimensions (3 spacelike, 1 timelike) are stable solutions in general relativity [14]. This means that small perturbations around the solutions tend to die off, and to be radiated away in the form of gravitational radiation. Although this is a very well established result in 4 dimensions, very little is known about the behaviour of solutions with event horizons in a 5-dimensional spacetime (we consider here 4 spatial and 1 temporal dimension). Gregory and Laflamme [37], [38] showed, using linear stability analysis, that at least some classes of black hole solutions in 5 dimensions—those known as black strings, and which have the structure "4-D black hole" x "a line"—are unstable under long wavelength perturbations in the "line" (string) dimension. Given that small perturbations of these solutions tend to grow, a natural question arises: what is 3For some models with differential rotation the upper limit can exceed 4MQ [5]. Chapter 1. Introduction 13 the end state of the evolution of the perturbed black string? There are some arguments supporting the idea that the state will involve some kind of bifurcation of the original solution [37], [38]. Were this true, it would represent a violation of Penrose's cosmic censorship conjecture, generalized to 5 dimensions, which posits that singularities of the type discussed above should generically be hidden from external view by event horizons. On the other hand an argument has been published by Horowitz and Maeda [48] suggesting that bifurcation cannot occur in finite time. In order to study the end stage of the perturbed black string, numerical calculations have been performed [18]. In Chap. 3 we summarize the main results found in that collaborative effort, focusing on the numerical solution of the 5-dimensional analogue of the constraint equations (1.18-1.19). 1.4 Relativistic Hydrodynamics Following [110], we can define a perfect fluid such that in local comoving coordinates the fluid is isotropic. In particular, assuming that the spacetime is Minkowskian (flat), we can write T" = pH, TXX = TYY = TZZ = P, (1.33) where is the proper energy density and P is the hydrostatic pressure. The first generalization of these stress tensor components will be for the case where each fluid element has an arbitrary spatial velocity, vl, with respect to some fixed (lab) frame. This can be achieved via a general boost (spacetime rotation), resulting in = (pa + P) u^u" + Prf". (1.34) Here is the fluid 4-velocity, satisfying vPu^ = -1, and rfv = diag [— 1,1,1,1]. Note that the equations for the conservation of energy and momentum can be written as T^"ILLL = 0, and that' this divergence involves an ordinary derivative since we are in flat spacetime. In order to extend this expression to a generally curved manifold we need only replace the Minkowskian metric by the general Lorentzian metric of the spacetime and partial derivatives with covariant derivatives. Thus, in a general curved spacetime, the stress energy tensor for a perfect fluid is = (pH+P)u'iuu + Pg^ (1.35) with (local) conservation of energy and momentum expressed by ( n ; , = 0- (1.36) Note that this expression is in accordance with the Einstein equations (1.1) since {G^u).^ — 0 by the Bianchi identities. In addition to this conservation law, a fluid can satisfy additional conservation Chapter 1. Introduction 14 laws depending on its composition. In the current case we will assume conservation of particle number (baryon number), which is expressed as (•/").„ =0, (1.37) where J M is the current associated with the particles; i.e. J» = p0u'i, (1.38) with po the rest mass density. As usual, in addition to these equations, we need to specify an equation of state for the fluid that relates the pressure to the energy densities. We also note that the above form for the perfect fluid stress-energy excludes the possibility of thermal or viscous effects—in other words we neglect all the non-adiabatic effects in the fluid (apart from shocks). The proper general relativistic incorporation of such terms is still a matter of open research [72]. In the strong gravity regime, the pressure and stresses are typically so large that we cannot assume the fluid is incompressible. In addition, for highly relativistic configurations, the pressure contributions to the stress tensor can be of the same order as those from the energy density. This makes general relativistic fluids behave very differently from the type of fluids that we encounter in everyday life, where stress energy tensors are dominated by the rest mass density of the fluid, where the assumption of incompressibility is often a very good one, and where characteristic three velocities satisfy General relativistic hydrodynamics (hereafter simply relativistic hydrody-namics, or relativistic hydro) involves the solution of equations (1.36).-(1.37) coupled to equations (1.1). There has been quite a lot of work in relativistic hydrodynamics because it is a good model for a lot of astrophysical systems, and is reasonably tractable computationally. A good review of the work and the methods that have been and are being used to study relativistic hydro can be found in [30]. Here we only briefly touch on some of the main efforts. May & White: The field started with the pioneering work of May and White [68], [69] who studied the collapse of a perfect fluid in spherical symmetry. In order to study the collapse they adopted Lagrangian coordinates (i.e. coordinates comoving with the fluid). This approach was very successful but, due to its Lagrangian nature, it is difficult to generalize to less symmetric spacetimes (for example, axisymmetric solutions) because the matter may follow complicated paths which can "tangle" the coordinates. May and White used a numerical approach based on finite difference approximations, which involves replacement of the derivatives in the equations with finite difference quotients. In order to deal with the shocks that generically arise during the collapse of perfect fluids, the Chapter 1. Introduction 15 continuum equations were modified to include terms that mimic viscosity [108]. The addition of this artificial viscosity tends to smooth out discontinuities over many finite-difference grid zones. This method of solution has also been used by other authors in order to investigate supernovae explosions [98], and neutron star collapse [70], all in spherical symmetry. Wilson formulation: Wilson developed one of the first Eulerian formulations for relativistic hydrodynamics that saw use in numerical work [111]. His approach involved writing the fluid equations (1.36-1.37) as a set of advection equations, i.e. equations of the type: dQi dQ^ dt + dxi - 6 " [ 1 - S 9 } where is again the fluid 3-velocity with respect to the lab (Eulerian) frame. The formulation included the definition of appropriate variables, {Qi(t, x^)}, that allowed the fluid equations to be cast in this specific form. Note that in general Si has terms involving the pressure gradients which are treated as sources. The original implementation of this formalism [111] was used to study the accretion of perfect fluid on to a rotating black hole. This was a "background calculation" in which the fluid's self-gravity was assumed to be negligible, so that the curved geometry acted on the fluid, but not vice versa. Wilson's numerical implementation was again based on finite differences, using so-called upwind derivatives for stability as well as artificial viscosity in regions with shocks. This formulation has been used extensively. Some of the systems that have been simulated using it include stellar core collapse [112], axisymmetric stellar core collapse [85], [28], [95], and coalescence of binary neutron stars using the conformally-flat approximation for the geometry [113]. Nakamura formulation: For the research described in Chap. 4 and especially Chap. 5, it is of crucial importance to review the work due to Nakamura and collaborators [63], [74], [75], [76], [92]. Nakamura and his co-workers used the so-called 2+1+1 formalism for the geometry and the fluid equations. This formalism is a variant of the 3+1 decomposition described in Sec. 1.1 and is designed for use with spacetimes having an axial Ki l l ing vector field (i.e. axisymmetric spacetimes). The formalism was originally developed for the Einstein equations by Geroch [34] and was extended for hydrodynamics by Maeda et al. [63]. In this approach the field equations (and the spacetime) are first decomposed with respect to the axial Ki l l ing vector field. The decomposition is very similar in spirit to a Kaluza-Klein reduction, in which effective matter fields and equations of motion are produced in the quotient 3-dimensional spacetime. After this initial decomposition is complete, the resultant 3-dimensional spacetime is decomposed using an analogue of the 3+1 approach (now, however, it is a 2+1 decomposition, which accounts for the name "2+1+1"). Chapter 1. Introduction 16 Nakamura's numerical implementation of the 2+1+1 equations also used finite difference methods—largely paralleling Wilson's approach—and artificial viscosity to treat shocks. The main application studied was the collapse of axisymmetric stars, including stars with large values of an-gular momentum. Their investigation showed that the collapsing star tends to form a ring structure which, in the case of moderate values of angular momentum, is enclosed by an apparent horizon. On the other hand, for initial stars that were sufficiently rapidly rotating, some evidence for the formation of naked singularities was found [74]. Valencia Formulation: A n important property of the fluid equations that none of the above studies took full advantage of is the fact that the hydrodynamical equations can be written in conservation law form (see App. Here, F^i is the flux associated with Qi along the x J direction and, crucially, Si does not contain any derivatives of the dynamical fluid variables. This was first exploited in the context of relativistic hydrodynamics by the Valencia group in [66]. Casting the equations in the above form allows the use of the Godunov approach (see [57], [58] for a summary of this and other so-called conservative methods) which ensures the correct jumps in values of dynamical variables across discontinuities, as well as the correct shock propagation speeds, (for these reasons such techniques are sometimes called High Resolution Shock Capturing, or H R S C , methods). One advantage with respect to formulations such as Wilson's, is that H R S C methods do not require the use of artificial viscosity for code stability, or to simulate extremely relativistic (v —> c) flows. Some examples of calculations are [91], [78], [31] and more recently [4]. In this thesis we will apply conservative methods to the hydrodynamic equations, and finite difference methods to the geometric equations, both within the 2+1+1 formalism. 1.5 Summary of Results In Chap. 2, we describe a new family of critical solutions resulting from the collapse of a massless scalar field with a particular potential that mimics angular momentum in spherical symmetry. These solutions are parametrized by an angular momentum coefficient, I. We have found that for each value of / we obtain a different discretely self similar critical solution. We also find that A ; decreases approximately exponentially with increasing values of I and, in fact, approaches zero in such a way that the critical solution appears to be periodic in the limit I —» oo. Moreover we have Chapter 1. Introduction 17 found scaling laws similar to (1.27) with ^-dependent mass-scaling exponents, 7J, that also decrease with I. In the collaborative work explained in Chap. 3 we have been able to check, by dynamical evolu-tion of black strings, the results found using perturbative methods by Gregory and Laflamme [37]. Not only have we directly verified that long thin black strings are unstable to perturbations along the string dimension, we have seen some indications that the spacetime evolves to a state (not necessarily the end-state) which can be described as a series of black holes connected by thin black strings. Unfortunately the code crashes at late times due to a coordinate pathology. We are thus unable to make any statements concerning the ultimate end-state of the evolution, since, at the time of the crash, the spacetime is still highly dynamical. Chap. 4 describes a spherically symmetric code for relativistic hydro whose purpose is to test the formalisms and algorithms, both at the continuum and discrete levels, subsequently used in the axisymmetric case. We have seen that, in our coordinate system (which is the natural restric-tion to spherical symmetry of the one used in the axisymmetric code), the standard conservation variables used in the Valencia formalism [91] must be be modified in order to get a well-posed set of geometric constraint equations. In addition, we have found that the Roe approximation used in the computation of numerical fluxes results in a scheme that is too dissipative, at least for grids with constant resolution, to maintain long term evolution of stationary solutions in our chosen coordinate system. In the concluding chapter, we present the equations for relativistic hydrodynamics written in the 2+1+1 formalism and in a way which is amenable for treatment using H R S C methods. The equations for rotational, hydrostatic equilibrium, along with an integrability condition, are also expressed in the same formalism. Finally, we describe our numerical implementation for the case of no rotation. At the current time this code is both too dissipative and too unstable to be able to use it for the study of the long term evolution of stationary solutions. However, we can evolve certain configurations of discontinuous data, and can also demonstrate that the code exhibits second order convergence for situations where the hydrodynamical flow remains smooth. Chapter 2. Scalar Field Collapse with Angular Momentum 18 Chapter 2 Scalar F i e l d Col lapse w i t h A n g u l a r M o m e n t u m 2.1 Introduction In this chapter we explain the main results of a new study of critical phenomena in gravitational collapse. The work that we present here is the result of a collaboration with M . Choptuik, W . Unruh and J . Ventrella. As described in Sec. 1.2, most studies of black hole critical phenomena (or related phenomena in other sets of nonlinear evolution equations) to date have' been performed assuming spherical symmetry as a simplifying assumption (exceptions are [1], [62] and more recently [20], [21]). This simplification has been adopted in most cases because accurate simulation of Type II critical solutions—which exhibit structure at all scales due to their self-similar nature—requires great computational resources. Since spherically symmetric spacetimes do not allow for angular momentum, very little is currently known about the role of angular momentum in critical collapse. For a few cases, most notably the Type I i solutions found in spherically symmetric collapse of a massless scalar field [33], or certain types of perfect fluid [41], [42], perturbative calculations about the spherical critical solutions suggest that all non-spherical modes, including those contributing to net angular momentum, are damped as one approaches criticality. In particular in [33], [41] and [42] using second order perturbation theory it was found that the angular momentum of the black holes produced has the following dependence as a function of the critical parameter p: LBH = Lo{Mp-p*)}(p-p*r, (2.1) where LQ [ln(p — p*)] is some quasi-periodic function that is family-dependent and n. « 0.76 is a universal scaling exponent which is larger than the corresponding 7 (scaling exponent for the black hole mass in (1.27)). These calculations thus suggest that, at least for small deviations from spherical symmetry, the resulting solutions at the verge of black hole formation should remain spherically symmetric in non-symmetric collapse. We also note that an axisymmetric numerical relativity code is currently being developed [21] to study non-perturbatively some effects of angular Chapter 2. Scalar Field Collapse with Angular Momentum 19 momentum in the critical collapse of a scalar field. Here a different approach is taken. Maintaining spherical symmetry, the equations of motion for a massless scalar field are modified by effective terms which mock up some of the effects of angular momentum. As described below, the procedure amounts to performing an angular average over the matter field variables—similar to what is done in [88], [81] and [107]—and results in an entire family of models, parameterized by a principal angular "quantum number", I (we will generally restrict I to take on non-negative integer values, although real-valued l's are also formally possible). We note that since the models remain spherically symmetric, we cannot use them to address the validity of the perturbative calculations mentioned above (e.g. equation (2.1)). Nonetheless, we find interesting results that may shed some light on the effects of angular momentum near the black hole threshold. Some of the main results that have been found are as follows. First, each value of the angular momentum parameter I, apparently defines a distinct critical solution. For I < 10, these solutions are found to be discretely self similar, with values of the echoing exponent A; (see (1.31) and the accompanying discussion) that rapidly decrease (approximately exponentially) as I increases. As a result, for large values of I, and for the time scales for which we are able to dynamically evolve near criticality, the threshold solutions become approximately periodic. In addition, and as expected for Type II solutions, we find that for I < 7 the masses of the black holes formed follow power laws of the type (1.27). As with the echoing exponents, for increasing values of I it is found that the mass-scaling exponent, ji, rapidly decreases, again approximately exponentially in I. The remainder of this chapter is structured as follows. In the following section we describe the recipe used to calculate the effective equations of motion, along with the regularity and boundary conditions imposed in the solution of these equations. In Sec. 2.3 we briefly describe the numerical code, the way the solutions have been analyzed, and then provide a summary of the results obtained for varying values of I. 2.2 Equations of Motion 2.2.1 Equations In order to derive equations of motion, scalar fields of the following form are considered: .$?(t,r,e,<l>)=tl>(t,r)er(e,<t>), m =-1,-1+1,... ,1-1,1 (2.2) where ©["(0, </>) are normalized real eigenfunctions of the angular part of the flatspace Laplacian with eigenvalue 1(1 + 1), and the index m labels the 2^  + 1 distinct orthonormal eigenfunctions for a Chapter 2. Scalar Field Collapse with Angular Momentum 20 given value of I1. By construction, the scalar fields ty™, are not, in general, spherically symmetric and we therefore do not study their collapse directly. Instead, our strategy is find effective equations for the single (t, r)-dependent quantity ip(t,r). To do so, for a specific value of I, we consider the stress-energy tensors, T^Kb, for the 21 + 1 fields 1 a b ~ 87T (2.3) where gab is the metric of the spacetime, V a is the metric-compatible covariant derivative and the non-standard factor of l/(87r) has been introduced to cancel the one in (1.1). Again by construction, the sum of these stress tensors T^ab = J2T(lm)ab, (2.4) m is spherically symmetric and thus depends only on ip(t,r), I, and the metric gab- The net result of this procedure is equivalent to performing an angle average over the individual stress energy tensors T^ab and then summing them, i.e. to computing T^ab = ^ ( T « m U ) (2.5) m where (/ (*> <t>)) = -tf Kd> & s i n (*) Mfy. (2.6) We can now compute the effective equation of motion for the field, ip(t,r), by demanding that the divergence of the total stress energy tensor is zero: V cT<'U = 0. (2-7) The equations for the geometric variables are determined from the 3 + 1 decomposition of the Einstein field equations explained in Sec. 1.1. For the current study we adopt Schwarzschild-like (polar-areal) coordinates, in which the metric (1.25) takes the form: ds2 = -a2(t,r)dt2 +a2(t,r)dr2 +r2d92 + r 2 s in 2 Odxj? . (2.8) Here a(t,r) is the lapse function and a(t,r) is the only non-trivial component of the 3-metric 7ij (both a and a are positive functions). Using this metric, the non zero components of the 1Note that, in general, GJn(6,(/>) will not be eigenfunctions of the azimuthal rotation operator (d/d<j>), since they are real. , Chapter 2. Scalar Field Collapse with Angular Momentum 21 8TT (21 + 1) 2 8ir aa stress-energy tensor for a general value of I are and the stress-energy trace is T(l) = j-(Z) i - i ( n 2 + <i>2) + ^ + i) V>2 ( 2 f + 1) (21 + 1) 8na2 ±(ll2 + $2)-1(1+1) cr n 2 - $ 2 ) , -02 (22 + 1) 8?r -2{U2-*2)-21(1 + 1)^ (2-9) (2.10) (2.11) (2.12) (2.13) In the above expressions, we have made use of the auxiliary variables, $ and II, defined as follows: = ^ , (2.14) Tl(t,r) = a-d4- (2.15) The dynamical equations of motion for these fields, which follow from the definition of $ as well as the wave equation for ip (which in turn can be derived from the vanishing of the divergence of the total stress tensor (2.7)) are then: ~dt dU d fa dr (2.16) (2.17) Note that the dependence of these equations on I is only through the last term in equation (2.17) which is proportional to 1(1 + l)/r2. This term can be thought of as the field-theoretic extension of an analogous term due to the angular momentum potential, l2/r2, in the 1-dimensional reduced problem of a particle moving in a central potential. As mentioned above, equations for the geometric variables result from the 3 + 1 decomposition of the field equations presented in Sec. 1.1, as well as from our choice of coordinates. Specifically, we have the following: 1 da a dr Ida a dr da (21 + 1) 2 (21 + 1) r n 2 + cE>2 + Z(Z + l ) ^ 2 r n 2 + $ 2 - / ( / + i ) ^ v 2 + 2r -2 1 2r — = (21 + l)mn$. dt (2.18) (2.19) (2.20) Chapter 2. Scalar Field Collapse with Angular Momentum 22 Equation (2.18) is the Hamiltonian constraint (1.18), which is used to determine the 3-metric component, a. Similarly, the slicing condition (2.19) fixes the lapse function a at each instant of time, and is often known as the polar slicing condition. It can be derived from the demand that Tr (Kab) = Krr+K6e+I^4> = Krr for all times. The Hamiltonian constraint and slicing condition, with appropriate regularity and boundary conditions, completely fix the geometric variables in this coordinate system. Equation (2.20) is an extra equation derived from the definition of Krr and the momentum constraint (1.19). In our numerical solutions, it is used as a gauge of the accuracy of our simulations, as well as to provide a replacement for the Hamiltonian constraint in certain strong field instances where the numerical constraint solver fails. In addition, we compute the mass aspect function, M(t,r), which serves as a valuable diagnostic quantity in our simulations. The value of this function as r —> oo agrees with the A D M mass (Arnowitt-Deser-Misner mass [3]), and more generally, in a vacuum region of spacetime, measures the amount of (gravitating) mass contained within the 2-sphere of radius r at time t. Moreover, 2M(t, r)/r is useful since its value approaches 1 when a trapped surface is produced and hence (modulo cosmic censorship), a black hole would form in the spacetime being constructed. We note that, as is the case with the usual Schwarzschild coordinates for a spherically symmetric black hole, polar-areal coordinates cannot penetrate apparent horizons, and in fact become singular as they come "close to" black-hole regions of spacetime, where 2M(t,r)/r —> 1. This fact does not present a problem in the study of critical behaviour in our models, since the critical solutions per se have max r [2M(t,r)/r] bounded away from 1. 2.2.2 Regularity and Boundary Conditions In addition to the above equations of motion, appropriate regularity and boundary conditions are needed. At the origin, r = 0, regularity and elementary flatness implies the following expansions: (2.21) l im a(t, r) l + r 2 a 2 ( i )+0 ( r 4 ) (2.22) l im a(t, r) i—>o l im tp(t, r) a0(t)+r2a2(t) + O(r4) (2.23) rlMt)+rl+2^i+2(t)+0(rl+4). (2.24) Chapter 2. Scalar Field Collapse with Angular Momentum 23 This leads to the following boundary conditions: a(t,0) = 1, (2.25) ^(*,0) =. 0, (2.26) ^(t,0) = 0, (2.27) ip(t,0) = Q(rl), (2.28) n(t,0) = O(r'), (2.29) f Ofr '" 1 ) for I > 1, *(*,0) = < (2-30) O(r) for Z = 0. In the continuum, our equations of motion are to be solved as a pure Cauchy problem, on the domain t > 0, r > 0, with boundary conditions at spatial infinity given by asymptotic flatness (i.e. that the matter fields vanish, and that the metric becomes that of Minkowski spacetime, as r —• oo). Computationally, we solve an approximation-to this problem on a finite spatial domain 0 < r < r m a x , where r m a x is some arbitrary outer radius chosen sufficiently large- that we are confident that the numerical results do not depend significantly on its precise value. At the outer boundary, then, the following condition for a is imposed: a(t,rmax) a ( i , r m a x ) = 1. (2.31) This can be viewed as simply providing a convenient normalization for a, since given a solution, a, of the slicing equation (2.19), ka is also a solution, where fc is an arbitrary positive constant. We note that although we have used (2.31) in order to perform the calculations, a different normalization convention—i.e. a different, and time dependent, choice of fc—has been used in order to perform the analysis of the solutions. Specifically, in the analysis we have used central proper time T (t) defined by: T(t)= [ a(i,0)di. (2.32) Jo This definition of time has a natural geometrical interpretation since r = 0 is invariantly defined by the symmetry of the spacetime. For the scalar field variables, II and approximate outgoing-radiation boundary conditions (Sommerfeld conditions) are used: -zr{t,rmax) +—{t,rmax) + = 0, (2.33) Ot OT ''max ^ f t , r m a x ) + ^ f t , r m a x ) + n f t ' r m a x ) = 0. (2.34) ot dr r m a x An important point in the derivation of the equations of motion is the fact that the eigenfuctions in (2.2) are discrete and the allowable values of I are only non-negative integers. Once the equations Chapter 2. Scalar Field Collapse with Angular Momentum 24 are obtained we have relaxed that constraint and have allowed I to take non-negative real values (the physical interpretation of such solutions, if any, is open for discussion). The solutions corresponding to non-integer values of I would have some degree of irregularity at the origin depending on the particular value of I chosen. This implies that only some finite number of derivatives with respect to r will be defined at r = 0. In our particular numerical implementation, which assumes that second derivatives of the variables are defined, we have been able to study the evolution of these systems as long as I > 3. 2.3 R e s u l t s 2.3.1 Numerics We solve equations (2.16), (2.17) for the scalar field gradients, equations (2.18), (2.19) for the geometry, and use (2.14) to reconstruct the field The system is approximated using second order centred finite difference techniques, and coded using R N P L [65]. Numerical dissipation of the Kreiss-Oliger [52] variety was included to damp high frequency modes, and it should be noted that this particular type of dissipation is added at sub-truncation error order, so does not effect the overall accuracy of the scheme as the mesh spacing tends to 0. For the current computations, the damping terms were most useful in regularizing the truncation error estimation procedure that occurs when adaptive mesh refinement ( A M R ) techniques are used. It was also crucial to impose the correct leading-order regularity conditions close to the origin, r — 0 (equations (2.29)-(2.30)) to keep the solution regular during the evolutions. Most of the simulations were done on a fixed uniform spatial grid rj = (j — 1) A r , j = 1,2, • • • , J , J = 1 + r m a x / A r with a typical number of grid points J = 1025, and typical r m a x = 100. In this grid we define discrete values {ipj}, {&j}, {I^} and discrete values of the geometric fields. The regularity conditions are imposed On the dynamical variables by: n i _ f v s a - i / a n , br . - o , [ 0 for I > 1, I I 2 = n 3 /2 2 ( ' " 1 ) for I > 2, (2.36) , 4 / 3 $ 2 - l / 3 $ 3 for I = 1, *i = •{ ' ' (2.37) 0 for Ij: 1, , $ 3 / 2 for I = 2, $ 2 = { (2.38) 1 $ 3 / 2 2 ( ' - 2 ) for l>2. These have been calculated using the regularity conditions (2.29-2.30) and the finite difference Chapter 2. Scalar Field Collapse with Angular Momentum 25 approximation for the first derivative with respect to rq (q taking values 1, I — 1 and / — 2) at r = 0: df_ _ (2-" — 2") h + 29/2 - 2 - « / 3 dr" ~ (2« - 1) A r ? For small values of the angular momentum parameter—specifically for I < 2 based on that described in [17] was used. 2.3.2 Families of Initial Data Our study involved the evolution of 6 different one parameter families of initial data, each defined by an initial profile tp(0,r) as listed in Table 2.1, with specific values of the parameters appearing in the profile definitions as given in Table 2.2. In addition to ip(0\r) we need to provide 17.(0, r) Family Form of initial data, ip(0, r) P (a) A e x p ( - ( r - r 0 ) 2 / a 2 ) A (b) -2A (r - r 0)/CT 2 exp ( - ( r - r 0 ) 2 / a 2 ) A (c) Ar2 (atan(r — ro) — atan(r — ro — cr)) A Table 2.1: Families of initial data and the parameter p that is tuned to generate a critical solution. to complete the specification of the initial data. In all cases we chose 11(0, r) to produce an approximately in-going pulse at the initial time: n(0,r) = $(0,r) = 5W). (2-40) Initial Data Family Parameters 1 (a) ro = 70.0, cr = 5.00 2 (b) ro = 70.0, a = 5.00 3 (c) ro = 70.0, cr = 5.00 4 (a) ro = 40.0, cr = 10.0 5 (a) ro = 40.0, CT = 5.00 6 (a) ro = 70.0, CT = 10.0 Table 2.2: Initial data used in our investigations. The family label is explained in Table 2.1. As previously mentioned, all of the initial data families listed in Table 2.1 have a single free parameter, p, and, as is the usual case in studies of black hole critical phenomena, for any given family we observe two different final states in the evolution, depending on the value of p. For values (2.39) -an A M R algorithm Chapter 2. Scalar Field Collapse with Angular Momentum 26 of p > p* the maximum value of 2M(t, r)/r approaches 1 implying that an apparent horizon is about to form. On the other hand if p < p* the scalar field completely disperses, and leaves (essentially) flat spacetime in its wake. The solution that arises as p —> p* then represents the threshold of black hole formation and, by definition, is the critical solution. We note that these critical solutions are not t —» oo end-states of evolution; rather they persist for only a finite amount of time, and, in fact, are unstable, heuristically representing an infinitely fine-tuned balance between dispersal and gravitational collapse. 2.3.3 Analysis We have calculated p* for the different families of initial data described above, and for different values of I, via bisection (binary search), tuning p in each case to a typical precision of (p — p*) /p ~ 10~ 1 5 (which is close to machine precision using 8-byte real floating point arithmetic). As in the case for I = 0 (where the equations of motion reduce to those for a single, non-interacting massless scalar field, as studied in [17]), the critical solutions for values of / < 9.5 are apparently discretely self similar (DSS). As discussed in the introductory chapter, DSS spacetimes are scale-periodic, meaning that any non-dimensional quantity, Z, obeys the following equation for some specific values of the parameters A and T*: where T is central proper time as defined by (2.32), and T* is the "accumulation time" of the self-similar solution. Note that this is the same scaling invariance as (1.31) but written in the original coordinates given by (2.8). In (2.41) the,integer n denotes the echo number. We also note that due to the discrete ip —> —ip invariance that is exhibited both by the equations of motion as well as the critical solutions themselves, if A is the echoing exponent for which formula (2.41) is satisfied with Z(T,r) = ip(T,r), then the geometric quantities a(T, r), a(T,r), 2M(T,r)/r obey (2.41) with an echoing exponent A / 2 . In order to extract A from our simulations, we use the observation that certain geometric quantities will achieve (locally) extremal values on the spatial domain at discrete central proper times Tn given by where To is the time at which one starts counting the echoes. Specifically, A and T* have been com-puted by a least squares fit for the times Tn at which max r [2M(t, r)/r] achieves a local maximum in time, i.e. by minimizing: (2.41) Tn -T* = (To - T * ) e " A / 2 . (2.42) (2.43) 71=1 Chapter 2. Scalar Field Collapse with Angular Momentum 27 2.3.4 Results Table 2.3 summarizes the values of A ; we have estimated using this procedure; the data are also graphed in Fig . 2.1. Again, note that the reported values for A ; have been calculated using central proper time T instead of proper time at infinity (the parameterization used in the numerical evolutions per se). Also the reported uncertainties have been estimated from the deviations in the A ; values computed across the the six different families of initial data. The first entry in Table 2.3 (I — 0) corresponds to the original case studied in [17]. The second one (I = 1) is apparently the same solution found for the self-gravitating collapse of an SO(3) non-linear cr model, assuming a hedgehog ansatz'[50], [61]. The remainder of the solutions (for the other values of I) are, to the best of our knowledge, new. As was also discussed in the introduction, systems exhibiting type II critical behaviour, where the critical solution is self-similar, generally also exhibit power-law scaling of dimensionful quantities in near-critical evolutions. For example, we can expect the black hole mass, M B H , to scale as M B H = C ( p - P y (2.44) for super-critical evolutions as p —» p*. Here G is a constant that depends on the family of initial data while 7/ is a universal exponent for each value of I, i.e. independent of the specific initial data family used to generate the critical solution. We have observed such scaling in at least some of our computations, but, following Garfinkle and Duncan [32] have found it more convenient to extract 7; by monitoring the maximum value of the trace of the stress tensor, T , which, from the Einstein equations, is proportional to the maximum value of the Ricci curvature. On dimensional grounds T (defined by (2.13)) and R should both scale with an exponent —27. This technique has the advantage of being more precise than a strategy based directly on (2.44) since we can calculate the trace of the stress-energy more accurately than the mass of the black hole formed, and can perform the computation using sub-critical evolutions, where the gradients of field variables generally do not become as large as those in the super-critical cases. The values of 7; as a function of I are listed in Table 2.3 and are plotted in Fig . 2.2. As is characteristic of type-II critical solutions exhibiting discrete self-similarity, 2M(t,r)/r oscillates at higher frequencies and on smaller spatial scales during the course of an evolution in the critical regime. As has already been noted, as I increases, the echoing exponent A ; decreases rapidly. In addition, we observe that the maximum and minimum values between which the spatial maximum of 2M(t, r)/r oscillates increase with I (see Fig . 2.3) indicating that the critical solutions are becoming increasingly relativistic as the angular momentum barrier becomes more pronounced. The amplitude of the oscillations between these extremal values decreases since m i n r [2M(t,r)/r] increases more rapidly than max r [2M(t,r)/r] (see Fig . 2.3). Chapter 2. Scalar Field Collapse with Angular Momentum 28 i i r i i r n r i — i r i i r 0 <1 o 2 A * * I . D . ( i ) I . D . (2) I . D . (3) I . D . (4) I . D . (5) I . D . i (6) 1 1 1 i A i* i A A A i i i 0 2 4 6 8 Figure 2.1: Values of l o g 1 0 (A;) versus I. In this figure we can see that A ; decreases almost exponentially with I. The different lines represent different families of initial data. Assuming universality, the differences between the values calculated for the different families provides one measure of error in our determination of A ; . Chapter 2. Scalar Field Collapse with Angular Momentum 29 o O 1 i 1 1 1 1 1 i i i 1 — • — - • -- -I . D . 1 i i 2 • I . D . ( 2 ) — A I . D . (3) • — • I . D . (4) i x — I . D . (5) 1 A o I . D . (6) 1 A 3 — S • 1 i i 1 i i \ i i 8 I i 0 2 4 6 Figure 2.2: Values.of l o g 1 0 (7;) versus I, where 7; is the scaling exponent defined by (2.44). As for the case of the echoing exponent, A ; , 7; also decreases approximately exponentially with I. We note that due to lack of numerical accuracy we only can reliably compute 7; for I < 6.5 Chapter 2. Scalar Field Collapse with Angular Momentum 30 o 0.1 0.2 0.3 - 0 . 4 0.5 0.6 1 1 1 1 1 1 1 1 1 I 1 • A A A • A A -_ • — — • L o g ( M a x j M a x r | 2 M ( t , r ) / r | | ) -I A L o g ( M i n j M a x r j 2 M ( t , r ) / r | 0 I t i i f 1 i i i 1 i i i 1 i i i 1 i . i i 0 2 4 6 8 Figure 2.3: max t {max r [2M(t,r)/r]} in the critical regime as a function of I (solid line) and the same for mint {max r [2M(t,r)/r]} (dashed line). We see how the maximum and min-imum increase with I. On the other hand the amplitude of oscillation, given by their difference, apparently tends to zero with increasing I. Chapter 2. Scalar Field Collapse with Angular Momentum 31 5 o 4 3 A _L_ A I . D . (1 • I . D . (2 A I . D . (3 • I . D . U - I . D . (5 o I . D . (6 A A i • • • * • * A * A 1 A • A J L 1 * A A J I L O I A I. 6 8 o 5 A A Figure 2.4: Values of l o g 1 0 (T*) versus I. Tf is the "accumulation time" defined by equation (2.41) for a parameter I measured using central proper time. The accumulation time increases almost exponentially showing an increase in the stability of the critical solution as I increases. Chapter 2. Scalar Field Collapse with Angular Momentum 32 Empirically, we have also found that, as we increase I within a family of initial data, although A ; —> 0 and T* —> oo (see Fig. 2.4), the product T ; *A; appears to asymptote to a finite value. Note that ostensibly this product is family dependent (see Fig . 2.5), but that all DSS type-II critical solutions are universal only up to a global scale transformation (r, t) —* (kr, kt), with k an arbitrary positive constant. Choosing k = k(l) for each of the families so that max r [2M(t,r)/r] is attained at some fiducial radius ro, and considering the case I = 10, we find that the normalized asymptotic oscillation frequency, /o, defined by f0 =r0/(T*A) = 4.35 ± 0.01 (2.45) agrees for all families to better than 1%. Again, the quoted uncertainty is estimated from the variation of /o across the different families of initial data. We note that for I = 10 the near-critical solution stays at an near-constant radial position; our spatial resolution is not enough to resolve the small changes associated with the extremely small value of A ; . The radial location of max r [2M(t,r)/r] in this regime is the value of ro that we have used in (2.45). We also note that the observation that /o is apparently well-defined and unique (up to the usual rescalings associated with type-II critical solutions), is consistent with the empirical observation that as I increases, the critical solution becomes ever closer to a periodic solution. In particular, for a periodic solution we have A —> 0, and then Tn - T* = (T 0 - T*) e n A « (T 0 - T*) (1 + n A ) « - (T*A) n-T*, ' (2.46) where To represents the loosely defined time demarking the onset of the critical regime (and whose precise value is clearly irrelevant in the limit T* —> oo) which implies that the maximal value is attained at times Tn: Tn = - (T*A) n. (2.47) As shown in Figs. 2.6 and 2.7, from our simulations for I — 10, we cannot ascertain whether the solution is discretely self similar with A ; very small (< 0.0002), or periodic with period r = T * A . Naively at least, we expect that for I > 10, distinguishing between discrete self similarity and periodicity would become even more difficult. However, it is worth noticing that for I = 20 we have not yet seen evidence for (almost)-periodicity, with period T * A , but have instead seen a more complicated structure near criticality that is not yet understood. 2.4 Conclusions In this chapter we have discussed the results for a model that tries to incorporate some of the effects of angular momentum in the context of critical gravitational collapse. A new family of spherically-Chapter 2. Scalar Field Collapse with Angular Momentum 33 < * r— EH 30 25 20 15 10 1 r 1 r * A J L 8 • 8 I . D . (1) I . D . (2) I . D . (3) I.D. (4) I . D . (5) I . D . (6) • l • . • * * • A A A 1 o i — i — n — | i i r o o o O O o • A A A A A A A A I I 1 . 1 A A A _J L_ l _ 4 6 8 Figure 2.5: T*A/ as a function of Z. The fact that these products remain finite as T* —> co and A( —> 0 is evidence that the critical solutions tend to.a periodic solution in the limit I —> oo. Chapter 2. Scalar Field Collapse with Angular Momentum 34 n Figure 2.6: Fi t of the times Tn at which max r [2M(t, r)/r] reaches its maximum in time (denoted by triangles and scale on the left) assuming a periodic ansatz. Initial data family (1) was used with angular momentum parameter I = 10. We also plot the residuals of each data point with respect to the best fit (denoted by pentagons and a scale on the right) Chapter 2. Scalar Field Collapse with Angular Momentum 35 Figure 2.7: Fi t of the times Tn at which max,. [2M(t, r)/r] reaches its maximum in time (triangles and scale on the left) assuming a self-similar ansatz. Initial data family (1) was used with angular momentum parameter Z = 10. Again, we also plot the residuals of each data point with respect to the best fit (pentagons and scale on the right). Notice that the errors in the fit are of the same order as the errors to a fit assuming periodicity (Fig. 2.6), indicating that from our numerical results we are unable to distinguish between the two types of solutions for I > 10. Chapter 2. Scalar Field Collapse with Angular Momentum 36 symmetric critical solutions, (black hole threshold solutions) labelled by an angular momentum parameter, I, have been found. These solutions have similar properties to those for the I = 0 case originally studied in [17]: specifically, the solutions exhibit discrete self similarity, and have scaling laws for the values of dimensionful quantities in evolutions close to criticality. We have calculated the Z-dependence of the echoing exponents A ; , and the mass-scaling exponents 7;, finding that both decrease rapidly with increasing I, (at least up to I « 10). Moreover, we have argued that as I increases, the critical solution approaches a periodic solution. As we explained in the introduction, we expect that 7; = 1/A; where A; is the Lyapunov exponent associated with the single unstable mode of the critical solution for angular momentum parameter I. Therefore since 7; —> 0 with increasing I, we apparently have A; —> 00. This has the interpretation of increased stability of the critical solution for increasing I, i.e. the period of time that a solution can remain close to criticality (for a fixed amount of fine tuning) increases with I, this can be observed by the increase in Tt* (see Fig . 2.4). We believe that this can be interpreted as an effect of the angular momentum barrier which (partially) stabilizes the collapse against black hole formation. Chapter 2. Scalar Field Collapse with Angular Momentum 37 I Aj li 0 3:43 ± 0.05 0.376 ± 0.003 1 0.460 ± 0.002 0.119 ± 0.001 2 0.119 ± 0.003 0.0453 ± 0.0002 3 0.039 ± 0.001 0.020 ± 0.001 3.5 0.0224 ± 0.0009 0.0127 ± 0.0008 4 0.0132 ± 0.0008 0.0082 ± 0.0008 4.5 0.0077 ± 0.0007 0.0052 ± 0.0006 5 0.0044 ± 0.0007 0.0033 ± 0.0005 5.5 0.0026 ± 0.0006 0.0020 ± 0.0005 6 0.0015 ± 0.0005 0.0013 ± 0.0005 6.5 0.0009 ± 0.0005 0.0008 ± 0.0005 7 0.0006 ± 0.0004 -7.5 0.0004 ± 0.0004 -8 0.0003 ± 0.0004 -8.5 0.0002 ± 0.0003 -9 0.0002 ± 0.0002 -9.5 0.0002 ± 0.0003 -Table 2.3: Summary of the properties of the critical solutions computed for different values of I. Note that both the echoing exponents, A;, and the mass scaling exponents, 7/, rapidly decrease as I increases. Quoted errors have been estimated from the variation in values computed across the different families of initial data. Values of A/ have been calculated using central proper time T normalization of the lapse function, which is the natural normalization for type-II critical behaviour. For I > 6.5 we have not been able to calculate 7/, due to lack of numerical precision. Note that the I = 0 data agree with the original values calculated in [17], and that the I = 1 data agree with values calculated in [61] and [50] using models of completely different origin.-Chapter 3. Instability of a Black String 38 Chapter 3 Instability of a Black String 3.1 Introduction This chapter describes a 5-dimensional problem (4 space plus 1 time dimensions) which involves the dynamical evolution of a perturbed black string. This work was a joint effort with M . Choptuik, L . Lehner, R. Petryk, F . Pretorius and H . Villegas [18]. As mentioned in the introduction, the unperturbed solution (originally proposed by Myers and Perry [73]) has the structure "4-D black hole" x "a line", hence the nomenclature "black string". Gregory and Lafiamme [37] proved that black strings are unstable against long wavelength perturbations in the string dimension. Our work tries to answer the natural question that arises: What is the end state of the evolution of such an instability? In [37] it was conjectured, using entropy considerations,, that unstable black strings fragment producing spatially periodic black hole solutions of the type described in [7]. This process will produce a naked singularity and hence violate (4+1 dimensional) cosmic censorship [47]. However this analysis, which is based on the linearization of the equations around the black string solution, stops being valid once the perturbation grows sufficiently. More recently, Horowitz and Maeda have argued that this proposed bifurcation of the event horizon cannot be achieved in finite time [48]. Moreover they conjectured that the system is likely to evolve to a new stationary solution which is not invariant under translations along the string direction. Following this reasoning, Wiseman solved the equations for equilibrium [114] and found new non-translationally symmetric solutions. The A D M masses, (see [3] for definition of A D M mass) of these solutions are larger than the maximum mass for an unstable black string at the same compactification length, which would appear to rule them out as candidate end-states. In [18] we adopt a different approach, in order to investigate the end state, we dynamically trigger the instability and analyze the subsequent evolution of the spacetime. The organization of this chapter is as follows. In Sec. 3.2 we describe the form of the black string solution, define the coordinates that we use for the evolution, and define the form of the perturbation. In Sec. 3.3, we then explain the specifics of the numerical integration of the equations Chapter 3. Instability of a Black String 39 of motion for the model, paying special attention to the solution of the constraint equations, and the construction of an approximate event horizon finder, which represent my chief contributions to the collaborative effort. Finally, in Sec. 3.4, we summarize the results that were obtained from the study. 3.2 Equations The solution studied by Gregory and Laflamme [37] was first constructed by Myers and Perry [73] and describes a 5-dimensional vacuum configuration with an event horizon. Wi th a specific choice of coordinates—called ingoing Eddington-Finkelstein coordinates—the 5-dimensional line element for the solution takes the form ds2 = - (1 - 2 M / r ) dt2 + AM/rdrdt + (1 + 2M/r) dr2 + dz2 + r2dtt2, (3.1) where M is the mass of the black string, z is the so-called string dimension, and dtt2 = d02 + s in 2 0d(p2 is the metric of a unit 2-sphere. Note that the solution is z-invariant, so that the metric coefficients do not depend on the string dimension and that, paralleling the 4-dimensional Schwarzschild case, we can identify r — 2M with the radius of the black string. In what follows we assume that the string dimension is periodic, or in more mathematical terms, that it has S1 topology, i.e. that z — 0 and z = L are identified, for some constant L. Physically, this choice of topology is not crucial (so long as L is chosen sufficiently large to admit the Gregory-Laflamme instability) but computationally, the choice allows us to sidestep any issues associated with imposing approximate boundary conditions in the string direction. We also note that the solution is spherically symmetric in the sense that (d/dd)a and (d/d(j))a are Ki l l ing vector fields. Gregory and Laflamme found that the above solution is unstable to certain z-dependent perturbations. Specifically, they found that perturbations with wavelengths, Xz, sufficiently large compared to the string radius grow exponentially. We also note that there is some evidence that black strings could have one or more additional unstable modes that preserve the z translational symmetry [104], but if this is the case, then the numerical results discussed below would suggest that such modes grow "more slowly than those that break the symmetry. In order to determine equations of motion for the dynamical evolution of a perturbed black string, we use the decomposition technique summarized in Sec. 1.1, but adapted for a 5-dimensional spacetime (so that we have a 4+1 decomposition). In addition, to minimize computational com-plexity and cost, we retain the spherical symmetry of (3.1), but allow the metric to have both t-and z-dependence in addition to the original r-dependence. Thus, we write the spacetime metric Chapter 3. Instability of a Black String 40 ds2 = ( - a 2 + hAB(3ApB) dt2 + 2hAB(3AdxBdt + hABdxAdxB + hQdtt2, (3.2) where the indices A and B range over the coordinates {r, z}, and the metric functions a, (3A, hAB and hn now generally depend on t, r and z. Our choice of the usual angular coordinates {4>,0}, adapted to the spherical symmetry amounts to fixing 2 of the 5 degrees of coordinate freedom that we have in this problem. The three remaining degrees of freedom are fixed by choosing the lapse function, a, and the two non-trivial components of the shift vector, (3A. For all times, we choose a and f3z to be those associated with the static black string (3.1): a = (1 + 2M/r)~1'2 , 0* = 0. (3.3) The last coordinate choice is used to maintain the condition [hu/r2](t,r) = [ha/r2](Q,r) for all times. This produces the following algebraic constraint on the shift vector component (3T: / r = ^ . (3.4) The main motivation for choosing this particular gauge is to have well behaved coordinates that allow us to penetrate the horizon (as in [93] and [55]). This is crucial for our use of techniques that excise the region of the spacetime containing the singularity. We want to point out that a preliminary approach that fixed /3 r to its black-string value for all times gave rise to coordinate singularities. In particular the radial position of the apparent horizon showed significant variation, approaching zero at late times. In order to decrease this variation, and in the spirit of the so called "horizon-locking" [93] coordinates we chose condition (3.4). The field equations for our model are natural extensions of the 3+1 equations, summarized in Sec. 1.1, for the case of vacuum (vanishing stress-energy tensor). Specifically, the equations we need to solve are H = {4)R + K2 -Ki:jKi:> = 0 , (3.5) Mi ES DjKJ - DiK = 0, (3.6) : -laKij + DjPi + Dipj, (3.7) : a (WRij + KKij^j - 2aKikKkj - DiDja +DipkKkj + DjPkKki + ^ DkKij + aFkiFmjhkmH, (3.8) dKjj dt where Fki = —2SkrSvi. Note that in the last equation, equation (3.8), the Hamiltonian con-straint (1.18) has been added, and that this addition changes the structure of the principal parts of the differential operators involved in the equations. This has been done on the grounds that the Chapter 3. Instability of a Black String 41 equations thus modified have better stability properties than the "bare" equations (see [51], [96], and [56] for reviews). In order to efficiently study the non-linear dynamics resulting from the Gregory-Laflamme insta-bility, we find it convenient to be able to macroscopically "perturb" the black string solution (3.1). We do so by altering hn from its black string form. More specifically, at the initial time we set Here, A is a measure of the strength of the perturbation—in particular, for A = 0 the perturbation vanishes and we recover the black string solution—and q is an integer that controls the spatial frequency in the z-direction. To complete the specification of the initial data for our evolutions, we set all remaining metric components, hab and KAB—except for hrr, KRR and KQQ—to their black-string values. The values of hrr, Krr and K$e are then determined from the constraint equations as described in more detail below. 3.3 Numerical Implementation In this project we chose to perform free evolution, meaning that the constraint equations are only solved at the initial time. One motivation for this choice is that the solution of the Hamiltonian and momentum constraints, which are generally elliptic in nature, is computationally more costly than solving equations of evolution type. This is particularly the case when one wants to imple-ment the numerical solution on massively parallel, distributed memory machines; parallelization of single-grid, finite-differenced evolution equations is straightforward, while parallel treatment of elliptic equations (or other equations with "long-range interactions"), may need to be quite intri-cate. Having solved the constraints at the initial time, we then follow the standard practice in numerical relativity of computing the residuals associated with the discrete constraints as time progresses as one measure of the error in the solution. Here we are exploiting the property that the evolution equations preserve the constraints at the continuum level, and that this property should be preserved by our discrete scheme—to the order of truncation error—provided that the scheme is stable [16]. 3.3.1 Evolution In our numerical implementation we fix the mass of the unperturbed black string to one, i.e. M = 1. A second order finite-difference, Crank-Nicholson treatment for the evolution equations is used, with the resulting algebraic equations for the update values solved using an iterative (3.9) Chapter 3. Instability of a Black String 42 process. This process consists of a constant number of iterations (3) and can be viewed as a specific predictor-corrector scheme for the evolved variables [100]. We discretize on a uniform grid, using coordinates (x, z) where x — r/ (1 + r). Note that since x = 1 corresponds to r = oo, we can actually set boundary conditions at i° (see Figure 1.2). This compactification of the radial coordinate proved crucial to the efficacy of our scheme; our compu-tations done on a grid with finite range of the radial coordinate, r < 0 < r m a x , with the sort of approximate outgoing conditions often used in numerical relativity, lead to spurious results (includ-ing indications of time-independent, z-dependent solutions) which exhibited significant dependence on the specific value of r m a x chosen. The use of numerical dissipation of the Kreiss-Oliger form [52] was also essential to the stability of our numerical scheme, and was particularly important in two regions of the computational domain: close to the horizon of the black hole, as well as for x —> 1 (i.e. close to where the lack of spatial resolution caused outgoing disturbances to ultimately be represented near the Nyquist limit. In this latter case, the Kreiss-Oliger dissipation provides a natural and effective mechanism for "annihilating" the outgoing radiation while minimizing the. amount of artificial reflection back into the interior of the computational domain. Another important aspect of our implementation is the use of black hole excision [101]. As is the case for the 4-dimensional black hole discussed in Sec. 1.3, at r = 0 the black-string spacetime is singular: not only do some curvature terms go to infinity there, some of the metric coefficients blow-up as well. Dealing with such infinities numerically would seem to be fairly hopeless with current techniques, so to circumvent this problem, some region interior to the black string is excised from the computational domain. At the excision surface we do not need to set boundary conditions (this assumes that the time-locus of the interior boundary is either null or spacelike), and in practice we simply use a finite difference approximation of the evolution equations that employs appropriate one-sided difference formulae. The reason that this can work in principle is due to the fact that the interior of a black hole (or black string) by definition cannot influence the exterior, as discussed in Sec. 1.3. Now, for excision to work in practice, the excision surface needs to be chosen to be inside the event horizon. However, the event horizon is a globally defined structure—it cannot be located without knowledge of the entire spacetime. On the other hand the apparent horizon (the outer most marginally trapped surface) can be computed from quantities that are defined at a given instant of time and, assuming cosmic censorship, lies inside the event horizon [47]. In practice then, we locate the apparent horizon periodically and ensure that we are excising within this surface, and thus within the event horizon. The algorithm for locating the apparent horizon is described in detail in [18] and consists of a flow method that corrects the radius of an initial guess for the apparent horizon until the surface has an expansion below some specific tolerance. Specifically, if the radius of the apparent horizon is given by r = R(z), the function R(z) is corrected at every Chapter 3. Instability of a Black String 43 iteration by Rn+i = R n _ Q + A T (3 . 1 0 ) where R n + 1 [Rn] is the value of R(z) at iteration n + 1 [n], 9+ is the outward null expansion at iteration n and AT is the time-step for 'the evolution of the flow. As a final note, following the development of a stable, convergent serial (single processor) code, we constructed a parallel version using the C A C T U S Computational Toolkit [11]. 3.3.2 Determination of initial data: solving the constraints As described previously, once the components of KAB and Kab that we have deemed to be freely specifiable are given at t = 0, we solve (3.5)-(3.6) for the initial values of the remaining geometric variables, hrr, Krr and Kgg. This solution proceeds by iteration—each pass is comprised of three distinct stages, each of which involves the solution of one of the constraints for the appropriate geometric quantity (i.e. one of hrr,Krr,Kg^) treating all other quantities, including the other two constrained functions, as fixed. This process is iterated until the £ 2 _ n o r m (RMS value) of the residuals of all the equations falls below a certain tolerance. To initialize the iteration, we assign (unperturbed) black-string values to hrr, Krr and Kee- A t least for the weak perturbations (small A's in equation (3.9)) considered in our study, this initialization is good enough to yield convergence for the iterative process. We now describe this iterative solution process—in particular, the solution of each individual constraint equation—-in more detail. The constraint equations are discretized on a uniform grid of points {xi,Zj} with xt = (i — 1)A'X, i = 1,...,NX, and ZJ = (j — l)Az, j = 1,...,NZ. The mesh spacings in the x and z directions are A x E E — Xi = 1/(2(NX — 1)) and Az = Zj+i — Zj = L/(NZ — 1) respectively. We typically excise the region £ < 1/2, corresponding to r < M, from the computational domain, i.e. the range of our coordinates is such that x\ = 1/2, X J V ^ . = 1, z\ = 0 and ZNZ — L. We first consider the Hamiltonian constraint (3.5), which in our coordinate system can be viewed as an equation for hrr, and which has the form: F l ^ T + F ^ ^ - + F ^ r 9 - J ^ + (j£-)2 + (hrrf + F 6 h r r = 0 . (3.11) Here, the Fm, m = 1,. . . , 6, are functions that generally depend on all of the metric and extrinsic curvature components and their derivatives except hrr (and its derivatives). We discretize this equation to second order in the mesh spacings using a difference approximation centered at the Chapter 3. Instability of a Black String 44 points (xi+i/2, Zj). The resulting algebraic equations can be written as follows: [hrr]i+l,j — [h* rr\i,j »+l/2,j + t+l/2,j 9 ([^i-rli+l.j + [hrr]i,j) [h rr\i,j-\-l 2[hr rr\i,j — l (F4)i+l/2tj A z 2 ([h rr\i-{-l,j [^rjijj) ~ + [/lrr]i+l,j'+l ~ 2 [ / l r r ] j + l j + [/l rr]i+ [/I 'rrjz+.l,.? —1 + l/2,j 1 / [/irrjj+ij+i — [/lrr], 2 V 2 A z + 2 A z + Az2 [h. rr\i,j — \ 2Az [h rr\i,j — l 2Az -1/2.J rr\i-\-l ,j (3.12) Here {Fm)i+i/2 j a r e second order approximations of the functions Fm at points (2^+1/2,%)- As-suming the values [hrr]ij,j = 1,2, - • • Nz are known, the above system can be viewed as a set of Nz non-linear algebraic equations for Nz unknowns \hrr]i+\^,j = 1, 2, • • • Nz- We can solve this non-linear set of equations using an A^-dimensional Newton-Raphson method. Note that we thus solve the equations "line-by-line" in x^, starting at the inner boundary, i = 1, which is chosen well within the horizon of the string (as mentioned above, typically at r = M which amounts to x\ = 1/2), and where the boundary values, [hrr]ij, j — 1,... Nz, are chosen to be those corresponding to an unperturbed black string. Also note that the algebraic systems obtained in the linearization of (4.38) are: (a) tridiagonal, due to the nearest-neighbor character of our second order finite differ-ence approximations; and (b) cyclic, because of the imposed periodicity in the z-direction. These linear systems can be solved efficiently (in 0(NZ) time) using a cyclic tridiagonal linear solver [86]. We now turn attention to the r component of the momentum constraint (3.6), Mr — 0, which is viewed as an equation for kgg = Kgg/a (the factor of a is introduced to more readily maintain regularity at spatial infinity). In terms of this function, Mr = 0 can be written as dke dx + G2kgg + G 3 = 0, (3.13) where, again, the functions Gm, m = 1,2,3, do not depend on kgg or its derivatives. Note that this equation does not contain any z-derivatives of kgg and thus, for any value of z, is an ordinary differential equation in x, which is moreover first order in x. We discretize (3.13) using + 1/2.J [kg Ax + ( G 2 ) i + 1 / 2 j - \ ([kee]i+i,i + [kee}itj) + (G3)i+1/2J = 0, (3.14) which is a second order approximation centred at the point (xi+1/2, Zj). For any value of Zj, and assuming that the values [kgg]i+i,j are known, the above algebraic equations can be solved for Chapter 3. Instability of a Black String 45 [kge]ij. The boundary conditions [kee]Nx,j, j = at x = 1 (io) are again fixed to their (unperturbed) black string values. The last constraint equation, Mz = 0, is considered to be an equation for krr = r2Krr/a (again, the scaling of the extrinsic curvature component by r2/a is motivated by regularity considerations at x = 1): ff ? £ r + = Q ( 3 1 5 ) dz Once more, the functions Hm, m = 1,2,3, do not involve krr or its derivatives. This equation has the same structure as (3.13), but with the roles of x and z reversed, so that we now have, for any value of x, an O D E in z that we must solve. The second-order discretization used in this case is centred at points (xi, Zj+1/2)'-{ H l ) i . + i / 2 [krr}i,j+l- [krrkj + ^ . . ^ l_ ( [ ^ J . , + 1 + [fc^]. .) + (H3).J + 1 / 2 = 0. (3.16) In order to solve these equations, we again set the boundary conditions [&Vr]i,i, i = 1, • • • Nx to their black string values computed at z = zm\n, then solve for increasing values of j . 3.3.3 Finding Event Horizons As explained before, it is not possible to calculate the intersection of an event horizon with a given spacelike slice of a spacetime without knowledge of the entire spacetime. Specifically, in order to locate an event horizon one must determine the causal past of future null infinity, which in effect means determining the origin of all null geodesies that reach I+. Any region of spacetime not contained in the causal past of X+ (i.e. the "exterior universe") lies within a black hole, by definition, and the surfaces separating black hole interiors from the exterior universe are the event horizons. For the purposes of the black string calculations, it is interesting to attempt to study the actual dynamics of the event horizon—i.e. the time history of the intersection of the event horizon with our spacelike hypersurfaces, so a method that provides a good approximation to the location of the event horizon is needed. Here we describe one technique that we have used, following [49], to do just that. The method involves approximating the location of the boundary of the causal past of some r = constant surface by following radial null rays. We first derive equations for certain null rays in our spacetime. In particular, an appropriate Lagrangian for radial rays (i.e. no motion in the z direction) in our coordinate system is £ = (-a2 + hrr(3r2) (t'f + 2hrr(3rt'r' + hrr (r'f , ' (3.17) where the prime denotes differentiation with respect to some affine parameter A. Since we are interested in null trajectories, we set C = 0. The equation for the radial position of the null rays Chapter 3. Instability of a Black String 46 is then r = ±-jj=-pT=BR(t,r,z)i (3.18) where the plus [minus] sign corresponds to outgoing [ingoing] null rays respectively. Note that this equation is expressed in terms of derivatives with respect to the coordinate time t, which are denoted by an overdot. Once the evolution of the spacetime has been calculated, equation (3.18) for the outgoing case is integrated backwards in time, for all values of Zj, with initial conditions r = ro at the maximum time i m a x achieved in the evolution. We chose ro to be outside the horizon and close to x = 1, as an approximation to J+, or inside the horizon and close to the excision surface. In both cases, at least for our spacetimes, the evolution backwards in time accumulates at the event horizon. The integration is done using a second order Runge-Kutta scheme. Specifically, for a ray with constant coordinate Zj, in order to calculate r™ (the radial position at time tn) from the value r n + 1 (the radial position at time tn+1) with tn = tn+1 — At we use the following approximation: kx = -AtR(tn+1,rn+1,z), (3.19) k2 = -AtR(tn+1 -l/2At,rn+1 +"l/2fei,z), (3.20) rn = r « + i + f c 2 . (3.21) Here, R is defined by the right hand side of equation (3.18). In order to calculate R(tn+1, rn+l, z) and R(tn+1,rn+1,z) we need to determine values of a, (3r and hrr at those coordinate positions. Since the functional form of a is specified a priori, we can use that closed-form specification for the needed lapse values. In order to calculate the required values of /3 r and hrr, we use second-order (bi-linear) interpolation in the (xi,Zj) mesh. Notice that integration of the equations of motion for the rays backwards in time does not give us all of the causal past of r = ro, but only the part of the past that can be reached by purely radial rays. However, we believe that for the spacetimes that we have computed, and particularly since the event horizon is an attractor with respect to backwards integration, tracking these rays allows us to rather accurately locate the horizon. For related discussions of approximate event horizon location in the axisymmetric, four-dimensional case, see [12], [49],[60]. Finally, we should point out that the technique of backwards integration is absolutely crucial to the accuracy and efficiency of this strategy—forward integration of outgoing null rays becomes an increasingly ill-posed problem (especially at the numerical level) for rays emanating from regions closer and closer to the event horizon. Chapter 3. Instability of a Black String 47 3.4 Results Our code was thoroughly tested and showed second order convergence in the mesh spacings as expected. In particular, second order convergence of independent discretizations of the equations of motion was demonstrated, as was second order convergence of the discrete constraint equation residuals. This provides a stringent test of the correctness of our implementation. A l l of the calculations performed in this study have initial data as defined by (3.9) with specific parameter values A = 0.1, q = 1, ro = 2.5 and 5r = 0.5. Using these parameters, we have been able to recover the main perturbative results found by Gregory and Laflamme [37]. Specifically, we found a critical value L c for the string length which is within 2% of the value reported by Gubser [39], L c « 14.3M. Fig . 3.1 shows the rn3.xirnu.rn, Rxn&xi rninirnurn, Rmini values'of the areal radius of the the apparent horizon, as well as the following parameter defined in [39]: In the figure one can see that for a string length marginally larger than the critical value, L c , the apparent horizon gets increasingly distorted as the evolution proceeds, i.e. the maximum value of the areal radius grows while the minimum decreases. On the other hand, for a value of L slightly smaller than L c , the evolution is evidently (physically) stable. In order to most efficiently study the non-linear regime in an attempt to determine the ultimate fate of an unstable black string, we want to have the instability growing as fast as possible. We thus decided to study configurations with a string length L = 1.4LC since perturbation theory predicts that the fastest growing mode has a wavelength close to that value [37]. In Fig . 3.2 we show sequences of embedding diagrams illustrating the evolution of the apparent horizon for a calculation with L = 1.4LC. In this diagram angular dimensions have been suppressed, and new coordinates (f, z) are used so that the coordinate distance along the curve corresponds to proper distance along the apparent horizon. From these plots, we can get some sense of how the perturbation grows, and how there are indications that the late-time solution may be approaching a series of black holes connected by thin black strings. Unfortunately our code crashes soon after the last frame shown in Fig . 3.2 (and such crashes are generic in our late time evolutions of unstable black strings) . Thus we cannot conclude that this chain of black holes connected with a thin black strings is truly indicative of the final state, not least since the spacetime is still highly time-dependent at the time of the crash. Investigation of the behavior of the dynamical variables as £ —> i C rash indicates that the break-down is due to a coordinate pathology. First, we find that i C rash is not significantly dependent on the mesh spacings Aa: and and Az. This suggests that a numerical instability is not responsible for the (3.22) Chapter 3. Instability of a Black String 48 2.03 K 2.01 2 2 g 1.99 d 1.98 1.97 L=0.975LC i L=1.03L c ....... j i . i i I i , i i I i i i . I : 0 100 200 300 t Figure 3.1: Values of the maximum, Rma.x, and minimum, R m i n , areal radius of the apparent horizon as a function of coordinate time for two values of the string length. This calculation has been done for a black string with M = 1. The dotted line corresponds to a string length larger than the critical value, L = 1.03LC: one can clearly see the increase of Rmax and decrease of i ? m i n as a function of time. The continuous line corresponds to a string length smaller than the critical value, L = 0.975LC. In this case the evolution is stable and the small temporal variations observed can be understood in terms of a combination of numerical errors and the "relaxation" of the black string from a slightly excited state induced by the perturbation. 0.015 r< 0.01 0.005 Chapter 3. Instability of a Black String 49 F t=0 H 1—h F t=40 H — I 1 — h F t=80 :—H—I 1 — I — — I 1 1 — I — — I — I — I 1 — — I — I 1 1 —I 1 1—B L t=113 IS-H — I — h - H — I — — I 1 h - M H—i—i——i—i—i—i—\—i—i—i—i——i—i—i—a Figure 3.2: This figure shows snapshots of the apparent horizon, computed in coordinates such that coordinate length of the curve corresponds to proper length along the apparent horizon, and where the two angular dimensions have been suppressed. Note that the figure extends for two periods in z (i.e. the z-span is 2L = 2.8L C ), and that the portion of the curve plotted for negative values of f is included only for better visualization of the horizon dynamics. Chapter 3. Instability of a Black String 50 crash. Second, curvature scalars are calculated, and they appear to remain finite at all events of the evolution indicating that no physical singularity is produced within the computational domain. Finally, to this point in the analysis, only the dynamics of the apparent horizon, and not the event horizon, has been considered. However, Fig. 3.3 shows a comparison of results computed using the approximate event-horizon-locator described in Sec. 3.3.3, and those from apparent horizon location. For the event horizon location, outgoing null rays are traced back in time from 3 different "initial" (actually final) surfaces. Two of these surfaces are defined by 7*1 = 10, r 2 = 4 and lie outside the apparent horizon, while the third is 5 grid points (in our numerical coordinate Xi ) away from the excision surface and inside the position of the apparent horizon. The figure shows how quickly the null surfaces traced backwards from the initial cylinders converge to one another, as well as to the apparent horizon. This indicates that the apparent horizon is indeed a good approximation to the intersection of the event horizon with a given spacelike hypersurface in the spacetimes we have constructed. Chapter 3. Instability of a Black String 51 CO t=150 t=158 t=164 t= 113 t=138 t=144 L t=o : C l ' C2 : AH 1 - - - - - C3 t=40 t = 80 1 1 1 1 1 I I 1 1 1 1 I I . ! 1 1 - 1 0 - 5 0 5 10 Figure 3.3: Plots of the apparent horizon (labeled AH) and estimates of the event horizon location ( C l , C2 and C3) in coordinate space (in contrast to the embedding coordinates used in Fig . 3.2). Here, the C l (C2) curve marks the evolution of the outgoing radial null rays for the final t = 164 surface with ro = 10 (ro = 4 ) . C3 denotes the evolution of outgoing radial null rays, emanating from a surface just inside the apparent horizon at t = 164. Thus, moving backwards in time, these curves should asymptote towards the event horizon of the spacetime. These plots suggest that for most of the evolution (at least), the apparent horizon is an excellent approximation to the event horizon. Chapter 4. General Relativistic Hydrodynamics in Spherical Symmetry 52 Chapter 4 General Relativistic Hydrodynamics in Spherical Symmetry 4.1 Introduction In this chapter we describe the code we have developed for solving the fully coupled equations of general relativistic hydrodynamics in spherical symmetry. From the pioneering work of May and White [68] to more recent codes using H R S C methods, such as that Noble's thesis [79], there have been many different implementations of spherically symmetric general relativistic hydrodynamics that have been used to study a wide variety of problems. These problems include supernova explosions [68], [98], the structure of neutron stars [91], and critical collapse [78], [46], [79]. Overall, these studies have been very successful and many results have been obtained. We view our development of yet another spherically symmetric relativistic code as a logical first step towards our ultimate goal of studying axisymmetric self-gravitating hydrodynamics. The spherically symmetric code serves two main roles: (a) it allows us to experiment with the same formalism, and numerical schemes, as well as the same type of coordinates used in the axisymmetric case, within the context of a much simpler model; and (b) it provides us with the means of computing "benchmark" results that can be used, to test and calibrate an axisymmetric code, provided that spherically symmetric initial data is evolved by the latter. In the spirit of (a), the code has proven to be quite useful since it has allowed us to identify an appropriate set of variables describing the state of the fluid, so that geometric constraint equations actually have solutions in the strong-field regime. We note that existence of solutions is not always guaranteed for the case of non-linear elliptic equations, as emphasized by York [115] for the particular case of the constraints of general relativity. The remainder of this chapter is organized as follows. In Sec. 4.2 and using the 3 + 1 decom-Chapter 4. General Relativistic Hydrodynamics in Spherical Symmetry 53 position explained in Sec. 1.1, we discuss the equations that determine the geometry, as well as the Valencia formalism for the treatment of the fluid variables (see Sec. 1.4 and App. A ) . In the following section, we detail the numerical approach used to solve the equations, focusing on the H R S C method and finite-volume discretization used for the hydrodynamics. In 4.4.1 we explain the problem we encountered using a standard approach to solve the geometric constraint equations to determine initial data, and describe how this problem was resolved through the introduction of new dynamical variables for the fluid. Finally, in Sec. 4.4, we summarize some of the tests that have been performed in order to check the reliability of our numerical implementation. 4.2 Model/Equations As discussed in Sec. 1.4 a perfect fluid has a stress energy tensor of the form = pohu^u" + Pg^, (4.1) where po is the rest mass density, b is the specific enthalpy, P is the pressure and wM is the fluid four velocity. Note that the specific enthalpy h can be written in terms of the specific internal energy, e, as h = 1 + e + P/ po-For this study, we choose so-called maximal/isotropic coordinates in which the spherically sym-metric, time dependent metric takes the 3+1 form ds2 = [-a(t, r) 2 + V(*, r)4f3(t, r)2] dt2 + 2ip4pdtdr + ^ [dr2 + r 2 (dd2 + s in 2 Odcf)} . (4.2) As in our study of scalar collapse in Chap. 2, we adopt angular coordinates 6 and <f> that are adapted to the spherical symmetry. The radial coordinate, r, is fixed by demanding that the 3-metric be conformally flat, and the time slicing is fixed by requiring that the slices be maximal, which means that the trace of the extrinsic curvature vanishes Tr [Kij] = K\ = K\ + K\ + K++ = K\ + 2K% = 0. (4.3) As we will see, this last relation provides an equation for the lapse function that must be solved on each slice as the evolution proceeds. We choose this specific gauge since it is the natural restriction to spherical symmetry of the coordinates used in our axisymmetric implementation, which in turn are based on the coordinates used in [19]. Chapter 4. General Relativistic Hydrodynamics in Spherical Symmetry 54 Following [30], and references within, we now introduce hydrodynamic variables as follows D = poW, (4.4) Sr = p0hW2vr, (4.5) E = pohW2-P, (4.6) r = E-D, (4.7) ^ - ^ o + £ (4.8) auu a W = au* = (1 - vrvr)~1/2 . (4.9) We note that the definitions of D, Sr and r are motivated by our desire to cast the fluid equa-tions into conservation form, as briefly described in Sec. 1.4. Moreover, in our adopted system of coordinates, it is convenient to rescale these variables by an appropriate power of tp ( as we explain in more detail in Sec. 4.4.1). In particular, we will use rescaled variables D = tp6D, Sr = ip6Sr, f = i>&T and P = ip6P. Using the 3+1 formalism outlined in Sec. 1.1, we now derive the equations that will determine the geometric variables. Due to our restriction to spherical symmetry, as well as to our choice of coordinates, we can implement a fully-constrained evolution—wherein all geometric quantities are determined at all times either from the coordinate conditions themselves, or from the constraint equations—and we choose to do so. Given the form of the metric (4.2), and the demand that Kli = 0, the Hamiltonian con-straint (1.18) can be written as 2 3 ' (f + D) i>" + -4>' + —Krr2ijj5 + 27r-^-—-—- = 0, . (4.10) r 16 ip while the r component of the momentum constraint gives ( ^ • ) ' + 3 - ^ - ^ - 8 7 r ^ = 0 . (4.11) In addition, the coordinate conditions give us two more equations. First, the slicing condition, derived from the demand that (4.3) hold for all t is 1 , (z (b + f + ZP) c2 \ ^ K r V ) « ' ] ' " \y (K-f + 4,1 ^ > + * V ( f + \ + j 5 ) j - = 0, (4.12) and fixes the lapse at each time. Second, the requirement that the equations for 7 r r and 700/V2 (which in the 3+1 approach follow immediately from the definition of the extrinsic curvature Chapter 4. General Relativistic Hydrodynamics in Spherical Symmetry 55 components (1-23)) be equal, so that conformal 3-flatness is preserved in time, provides the following O D E (ordinary differential equation) for the single non-trivial shift vector component, (3 = (3r: 2 r fP\' 3 a V T (4.13) These complete the set of geometric equations needed to perform the evolution. In addition, we have an evolution equation for the conformal factor which, again, follows from the definition of the extrinsic curvature component, Krr: iP = - ^ K \ + 1p'/3+^f3'. (4.14) The hydrodynamical equations can be calculated from local conservation of the fluid stress tensor as well as local conservation of the particle number = 0 . (4.15) (4.16) Again, given our restriction to spherical symmetry, two independent equations can be derived from (4.15) i a (Vh~ 1 =gdt i a -g dr \ •06 SR [ Vr - ^ ) + P a 6P + 2 f (f + D + P J S/--U + DY- + 2 P - \ , a \ J a r (4.17) 1 d f Vh\ -T + r=gdt \ r —g dr { V 6 ^ ( V + ^ V } ip6 ) ip4 a 2P? + 3 P + i>6 S? V>4 (f + D + P ) ra (4.18) Here g is the determinant of the metric (4.2), so that y/^g = atp6r2 sin 9. Equations (4.17) and (4.18) represent local conservation of momentum and energy, respectively. From the particle conservation equation (4.16) we get 1 d (Vh~ D v P 0. (4.19) ^g~dt yip6"J ' y/^gdr \ ^ In summary, the complete set of differential hydrodynamical equations that is discretized in Chapter 4. General Relativistic Hydrodynamics in Spherical Symmetry 56 Sec. 2.3.1 is *r 2 dt (r ^) ar2 dr \ a T ^ aj } (4.20) ohdt (r25r) 1 0 arz or a 6 P + 2 ^(f + D + P) i a , o . i s a r J at arz dr .« r. 3 P + + + S r — - ( f + D — + 2 - , i/> a \ / a r (4.21) j a r 2 f ^ r - ^ + A ; r J ip4 a a ra (4.22) tpA (f + D + p) As usual, the above differential equations need to be supplemented with regularity and boundary conditions. At r = 0 the following regularity conditions are imposed: V/ (t,o) = o, (4.23) a' (t, 0) = o, (4.24) K\ (t,0) = o, (4.25) /3(t,0) = • o, (4.26) D'(t,0) = o, (4.27) Sr(t,0) = o, (4.28) r ' ( t ,0) = 0. (4.29) As was the case for our study of scalar collapse, we will approximately solve our equations of motion on a spatially finite computational domain, 0 < r < r m a x . For the geometric variables,.we impose boundary conditions based on the requirement that spacetime be asymptotically flat in the limit r —> oo, and that our time slices be labelled so that coordinate time coincides with proper time at infinity. Specifically, we must then have A(t) l im ip(t,r) 1 + + 0 ( 0 , lim a(t, r) = 1 B(t) + 0(r~% l im P(t,r) C(t) 0(r (4.30) (4.31) (4.32) where A(t), B(t), C(t) are general functions of time (which in practice are not independent of each other). Following standard practice in numerical relativity [116], these fall off conditions (4.30-4.32) Chapter 4. General Relativistic Hydrodynamics in Spherical Symmetry 57 can be expressed as the following mixed (or Robin) boundary conditions, f + fc^=O(r-3)«0, (4.33) or r £ + ^ = O ( r - 3 ) « 0 , (4.34) or r f + ^ = O ( r - 3 ) « 0 , (4.35) or r which are independent of the particular form of the functions A(t), B(t), C(t). This technique allows us to implement (4.30-4.32) using (4.33-4.35) without previous knowledge of these functions. Alternatively, we also have used the following outer boundary condition for a(t,r): -<'•'•'-!:Z1/«v ( 4- 3 6 ) r^oo 1 + M(t, r)/{2r) where M(t, r) = 2r(ip(t, r) — l). This value corresponds to the value of the lapse for the Schwarzschild solution written in isotropic coordinates. For the fluid, we use boundary conditions based on the demand that the flow be purely outgoing at the boundary. We approximate this condition by assuming that the derivative of the conservation variables is zero at and beyond the boundary of the computational domain. At the discrete level we implement this condition using ghost cells, see Fig . 4.1, where we copy the values of the ghost-cell conservation variables from the last physical cell. The hydrodynamical equations derived above do not completely fix the evolution of the fluid; we must close the set of equations by specifying a functional relationship between the pressure on one hand, and the energy and particle densities on the other; i.e. we must fix ah equation of state. In this thesis we restrict attention to the so-called ideal fluid equation of state (EOS) P = (T-l)Poe, (4.37) where T is the adiabatic index that will be taken to be a constant in the range (1,2]. This choice of equation of state, which is an extension of P = (k-Q/m)poT (ICB being the Boltzmann constant) see [13] and [99], admits stationary solutions (in contrast to the ultrarelativistic EOS, P = (r — 1)PH, f ° r example). This is a key feature which makes it a popular choice in relativistic hydrodynamics [30]. In ending this section, we reemphasize that we have introduced so-called conservative variables q = {D,S,T}, in order to cast the fluid equations in conservation law form, i.e. in the form of (1.40). In particular, the conservative variables are in no sense independent of the' primitive variables, p = {po,vR,P} (or equivalently {po,vr,e}), but rather are functionally related to the primitive quantities via equations (4.4)-(4.9). Chapter 4. General Relativistic Hydrodynamics in Spherical Symmetry Nr-2 Nr-1 Figure 4.1: In this diagram we show the cells near the boundary of the computational domain in our numerical scheme. The squares denote the discrete values for the conservation variables and the crosses denote the values of the fluxes. We approximate the outflow boundary condition using ghost cells (shaded on the figure): the values of the conser-vation variables in the ghost cells are identical to the values in the last physical cell. 4.3 Numerics We now describe the algorithm used to solve the coupled system of equations presented in the previous section. We have already noted that due to the symmetry of the problem and to our choice of coordinates, the geometric variables on a given hypersurface Et can be calculated without resort to equations of evolutionary type, assuming that the values of the hydrodynamical variables are known at that time. In order to calculate the solution on a future hypersurface, Ef+At, we adopt an iterative process. The iteration consists of the following main steps: 1. Make an initial estimate for the geometric variables {Gi™ + 1 } at time t + At. 2. Treating the advanced values of the geometric variables as known quantities, evolve the fluid equations to get estimates of the fluid fields {Qjn+1} at t + At. 3. Treating the advanced values of the fluid variables as known quantities, solve the constraint equations to correct the values of { G i ™ + 1 } . Here {Gin+l,i = 1,2,3} = {a,(3,tp} and {Qjn+1,j = 1,2,3} = ib,S„f\, and are evaluated at the advanced discrete time, t = t + At. The two last steps are repeated until the £ 2 norm of the Chapter 4. General Relativistic Hydrodynamics in Spherical Symmetry 59 difference between the values of the fluid variables from two successive iterations is below some tolerance (typically 10~ 1 0 ) . In the remainder of this section, we first discuss the numerical solution of the geometric equa-tions (4.10-4.13), and then the numerical treatment of the hydrodynamical equations (4.17-4.19). The geometric equations (which are all ordinary differential equations in r) are solved using second order finite difference techniques. The fluid equations, on the other hand, are solved using a finite volume approach first developed by Godunov [36], which exploits the fact that the equations are written in conservation law form. 4.3.1 Geometric Equations We first note that, in accord with our 3+1 decomposition of the Einstein equations, neither of the constraint equations (4.10) nor (4.11) involves the lapse function, a or the shift vector component, P (the "kinematical" gravitational variables). Thus,'to solve for the geometric variables, we first view the constraints as a system of 2 coupled ODEs for the quantities ip a n d KrT, then solve that system iteratively, as described below. Once ip and KTr have been determined, (4.12) and (4.13) can be solved for a and p, respectively. In order to finite-difference the geometric equations, we introduce a uniform spatial grid {^i!^2; •••,I'NT.}, where r j + i = rj + A r , A r is the constant mesh spacing, r\ = rmin = 0, and i~Nr = fmax is the outer boundary of the computational domain. Adopting the usual finite difference notation / , = / ( r j ) , we use the following second-order ( 0 ( A r 2 ) ) finite difference approximation of the constraint equations 0, (4.38) 0, i = 2,...,Nr-l, 0, 0, (4.39) 0, i = l,...,Nr-l. These two sets of equations are solved iteratively by first updating the ipi,i — 1,2,-•• ,Nr, as-suming the (Krr)i,i = 1,2,-•• ,Nr are known, then updating the (Krr)i assuming the tpi are known. Fixing the values (Krr)i, equations (4.38) comprise a non-linear system for the unknowns Sipx + 4V>2 - 1p3 2Ar ipi+i - 2ipi + V>i-i. t 2 ipi+1 - ipi-i p 3 , r / 2 5 i n _ _ ( D i + n) Ar2 Ar + T ^ r ) i ^ + 2*- ^ 1pNr - 1 , 3lpNr - 4lpNr-l + 1pNr-2 rNr + 2Ar Ar Ar (ri+1iP2+1+riiP2)'2 -87rS r i + 1 / 2 Chapter 4. General Relativistic Hydrodynamics in Spherical Symmetry 60 tpi. This system is solved using a (global) JV r-dimensional Newton-Raphson method, where at each Newton step all Nr values ipi are simultaneously updated. Each of these updates requires the solution of a tridiagonal linear system, which is accomplished using a standard L A P A C K routine [2]. In the second stage of the constraint-iteration, given the boundary (regularity) condi-tion (4.25), equations (4.39) are solved using a pointwise Newton method for each of the unknowns ( i r r ) i + 1 , i = i , 2 , - - - j v r - i . We note that it would be possible to implement a scheme whereby both systems (4.38) and (4.39) are solved simultaneously using a global Newton iteration. The banded systems (bandwidth ~ 6) that would need to be solved at each Newton step could still be solved in 0(Nr) time. However, we have not explored this option. The slicing condition is treated using a second order discretization similar to that used for the Hamiltonian constraint: —3ai + 4c*2 — 0:3 3 ( a i + i - ai) r2+1/2ip2+1/2 - {ai 2Ar n2 „/,2 (4.40) - i ) r 2 _ 1 / 2 € 1/2 Vi A r ( r i + l / 2 - r ? - l / 2 ) -tf (K\)2 + 4ntf (Di + n + 3P 4) + 4TT (Srfi ONr ( T i + Di+ Pi) 1 - MNJ{2rNr) 0, 0, 2 , . . . , i V r - l , l + MNr/(2rNr) where Mjv r = 2 (i/>jvr — 1) n v r . For the bulk equations we have discretized the expression (l/r2) d/dr as 3d/d(r3). This is a particular instance of a standard technique in numerical relativity (originally due to Evans [28]) whereby terms of the form df(r)/dr with f(r) ~ rp as r —» 0, and for some integer p > 1, are rewritten as p r p _ 1 (d/d(rp))f(r), and then differenced. This approach generally leads to improved behaviour of numerical solutions near r = 0, since the difference scheme is, by construction, consistent with the leading order regularity behaviour of the differentiated function. Equations (4.40) comprise a linear tridiagonal system for the ai which can again be solved using a standard L A P A C K routine. Note that the discrete outer boundary condition used here derives from (4.36). Finally, we need to discretize equation (4.13) for the shift vector component f3. This is done by first introducing a new variable w = /3/r, with discrete representation Wi. In terms of u>, second-order discretization of equation (4.13) results in U>l Wj+i - Wj A r ai+i{K\)i+1 | ai{K\)i n+\ Ti = o, = o, l , . . . , J V r - l . (4.41) Chapter 4. General Relativistic Hydrodynamics in Spherical Symmetry 61 We solve these equations in two stages, first integrating with an arbitrary boundary (initial) con-dition to give a set of provisional values, u>i, then correcting the u>i to produce final values to* that satisfy the true boundary condition. Specifically, we set wi = 0, and integrate outwards. We then set • Wi=Wi + k, (4.42) where A; is a constant value chosen so that WNr - WNr-l + 2 WNr + WNr-i 0. (4.43) A r ' {rNr+rNr-i) This outer boundary condition is derived from (4.35) and the definition of w. Finally, we compute the values of the shift component using = w^i. 4.3.2 Hydrodynamic Equations Treatment of the fluid equations in the discrete domain requires special care. As we have already discussed, the hydrodynamic equations will quite generically develop discontinuities, even if the initial conditions are smooth. In order to handle such discontinuities numerically we have adopted a finite-volume approach, using Roe's approximation for computation of the numerical fluxes. For discussions of these'methods in general see [57], [58], for their application to special relativistic hydro see [67], and for the general relativistic case see [30]. Our approach produces a solver of so-called Godunov type, involving the solution of a Riemann problem at the spatial boundary of each of the discrete volume elements (cells). As we have previously noted (Sec. 1.4), Godunov methods are applicable to any set of hyperbolic evolution equations that has been written in conservation law form. In particular, our fluid equations (4.20)-(4.22) can be written as —2 57 ( r «) + —2 7T a r F = S -arz at arz or (4.44) Here q, F, S are 3-dimensional vectors of dynamical variables, fluxes and sources, respectively: Q = F = D,Sr,r \ , (4.45) (4.46) 0, I 6 P + 2 (f + b + p)]* + sA-(f + p)eL + 2r, a \ / a r .?T*-2P£+{ZP + ip* a a (f + D + Pj I r a (4.47) Chapter 4. General Relativistic Hydrodynamics in Spherical Symmetry 62 In addition we define a vector of primitive variables P = [po,vr,P] (4.48) which can be used toi specify the state of the fluid and that we call primitive variables. These variables, related to q by (4.4-4.9), are useful in order to compute the fluxes F. In our finite-volume approach, we discretize the spacetime region bounded by the hyper-surfaces Et and St+At with a set of uniform-size rectangular cells J O ™ ^ ^ 2 } having vertices { t V n + V i , 7 - i + 1 } , s e e Fig . 4.2. .n+2 .n+1 Ar •< = * >.'-•> S>-» n+1/2 ; • • ? ' / > • /+V/2 ' . ; At i+1/2 i+3/2 i+1 r. i+2 Figure 4.2: Detail of one of the control volumes C^^j2 used in order to discretize the hydrody-namic equations of motion. Note that the vertices of the cells correspond to locations (tn,ri), (tn,ri+i), (tn+1,ri) etc., and are the locations at which the discrete geometric variables are defined. In practice, the cell vertices are the locations at which the geometric variables (which satisfy finite difference equations) are defined. Integrating these equations over any control volume, C™^j2 , we get f — (r2q) droit + f — (ar2F) drdt = [ Sar2drdt, (4.49) lnn + l/2 Ot V ' n n + l/2 dr V ' / ^ n + 1/2 where we have used dCn^}!? = ar2drdt for the infinitesimal 2-volume element associated with the 1+1/2 cell C^yj2. Note that this is not the usual volume element, *J—g = ar2^6, associated with the metric (4.2), because a factor of I(J6 has been absorbed in the source terms S through our definition of the new variables D, Srr, etc. Chapter 4. General Relativistic Hydrodynamics in Spherical Symmetry 63 Using standard integral theorems this last equation can be expressed as (C +V? + 1 / 2Ar - g ? + 1 / 2 r 2 + 1 / 2 A r ) + / - n + l / 2 n+i/2 2 A j . £ , " + 1 / 2 n+1/2 2 AA - hn+1l2 n+1/2 2 A + A ^ (A SfT> (Fi+1 Oj+l % A t _ f i ai r i A t ) = S i + l / 2 « i + l / 2 r . + l / 2 A f A r ' ( 4 ' 5 ° ) where 9™+1)/2, F™+1/'2, and S™+lj2—which are the fundamental discrete hydrodynamical variables— are defined by 9m/2 = Z^-XZ r+1q(C,r)r2dr, (4.51) FT1'" - f + F(t,rt) a dt, (4.52) = 1 T + S(t,r) a r 2 dtdr. (4.53) Q i + l / 2 ri+l/2AtAr J t n J^ Eqns. (4.50), along with the above definitions', comprise the basic discretization adopted for eqns. (4.44). We now schematically write eqns. (4.50) in the following way: Ci% = 9r+i/2 + A£G: + + 1 1 / 2 2 , (4.54) where f t n + l / 2 _ A , (an+X,2r2 F n + 1 / 2 - a n + 1 / 2 r 2 F n + 1 / 2 , \ \ai+1 r i + 1 r t + 1 a- rirl j ~n+\/2 n + 1 / 2 ~ 2 T " + d « + l / 2 ai+l/2 ri+l/2air (4.55) Note that _p" + 1 , / 2 a n c l S\+lj2 depend on the fluid quantities at both the advanced and retarded discrete times, tn+l and tn, respectively. In order to approximate the fluxes, Ft + ^ , we use Roe's approximation [90], which is given by Fi = \{F (pR) + F (pL) - £ \Xa\coar,a) . (4.56) We will explain in detail below the calculation of the different elements that appear in this last expression. For the time being, we note that Xn and rjn are the eigenvalues and right eigenvectors, respectively, of the velocity matrix V — dF/dq. F(pR), F(pL) are computed from equation (4.46) using approximations p R and pL for the primitive variables. In fact, p R and p L are approximations to the primitive variables at the same spatial location—the location of the cell interface, r = rt— but are computed using values defined either to the right or to the left, respectively, of the interface. The process of computing approximations at the cell interfaces is known as reconstruction. In our case we use a "minmod slope limiter" type [106] of reconstruction which results in equations (4.61) Chapter 4. General Relativistic Hydrodynamics in Spherical Symmetry 64 and (4.62) for pL and pR. We note that that the calculation of the primitive variables, p, from the conservative variables, q, is not completely trivial due to the non-linear algebraic relationship between the two sets of quantities. The specific method that we use to compute p(q) is also explained below. Expression (4.56) results in a second order (in space) approximation for the fluxes, that holds in regions of smooth flow, and away from any extrema of the functions being reconstructed. On the other hand, computation of the source term, S^+y^, using the values <7™+1/2 yields only a first order (in time) approximation. To maintain overall second order accuracy in the cell size (again, in regions of smooth flow and away from any maxima in q), we decompose the time step into two sub-steps: = ff+i/a + f ^ + 1 / 2 , (4.57) C++i/2 = 9r + i/2 + A t ^ + + 1 1 / 2 2 . (4.58) G™+y2 corresponds to expression (4.55) evaluated using Ft , Si+y2 and a™, while G™^j2 is com-puted using F™+ 1^ 2, S™+y2 and a™ + 1^ 2 (an interpolation of the lapse function at the half time step). F " + 1 ^ 2 and S™+y2 are calculated from the conservative variables obtained from (4.57) and their corresponding primitive variables. This completes the description of the basic update scheme for the fluid variables q™+y2-We conclude with two additional remarks concerning our numerical scheme. First, we note that the flux for the Sr equation (4.46) is actually split into two distinct pieces, one that contains a term that goes as 2P/r, and the second that absorbs the remaining terms. The first term is manifestly divergent as r —> 0 and directly cancels with the analogous term appearing in the source of (4.47). Second, we observe that difference quotients such as (ai+1r2+1Fi+1 - airfFi) ri+i/2Ar (4.59) (4.60) which appear in G i + 1 / 2 (see (4.55)) are rewritten in the following way 3 (ai+ir2+1Fi+1 - a^Fi) r3 — r3 ' ri+l ri in the same spirit as for the finite difference case explained in the discussion of (4.40). Calculation of the Roe flux We now explain in more detail how to compute the different expressions appearing in formula (4.56). Fig. 4.3 shows the main steps in the calculation of the Roe flux. The reconstructed values pR and pL—from which the quantities F(pR) and F(pL) that appear in the Roe approximation are Chapter 4. General Relativistic Hydrodynamics in Spherical Symmetry 65 calculated—are computed using Pi+l/2 Pi+1/2 = Pi+o-i (ri+i/2 -n), = Pi+i + o-i+x (n+1/2 - n+i) • (4.61) (4.62) in equation (4.46). Note that pt are the primitive values defined at the cell centres, their computa-tion from the conservative variables is explained in detail in the next section. In equations (4-61) and (4.62), cr, is given by o-i = minmod ( s i _ 1 / 2 , s i + i / 2 ) . with i+l/2 — Pi+1 - P i (4.63) (4.64) Finally the minmod function is defined by 0 if minmod(a, b) = j a if b if oh < 0 |a| < |b| and ab > 0 (4.65) \a\ > \b\ and ab > 0. Note that this function is used as a "slope-limiter" in order to decrease spurious oscillations that may appear at discontinuities. If the two slopes <Xj and (Tj + i have the same sign, then the one with smaller absolute value is used to linearly reconstruct the fluid variables. On the other hand if the slopes have differing signs (e.g. at an extremum), a first order reconstruction is performed, i.e. the values of the fluid at the cell centre are assigned to the cell interface. At the extrema, this produces a reduction of the accuracy of the overall scheme from second to first order in the mesh spacing A r . This reconstruction procedure introduces a certain amount of dissipation in the overall scheme [79]. We also note that this is by no means the only viable way of reconstructing; for discussion of other approaches see [58] and [77]. We now explain in detail the characteristic structure of V, which is needed to compute the Roe approximation for the numerical fluxes. The eigenvalues, A Q , and right eigenvectors, r)a, oiV are (see [30]) A 0 A± av 1 — v2c2 vr (l-c2)±cj(l-v2) 1 1 — v2c2) — vrvr (1 — c2) ^ 0 l,hWCT±,hWAr+ - 1 )• (4.68) (4.69) Chapter 4. General Relativistic Hydrodynamics in Spherical Symmetry 66 A P i - 1/2 q ^>p B 3-1 Pi-1 P: P: P P C i i • 1 • 1 • / \ / \ A F i - i F i F i + i Figure 4.3: This figure illustrates the main stages in the computation of the Roe approximation for the numerical flux F. First, in stage A , the primitive variables are computed from the conservative variables at the location of each cell center. In stage B , two approximations for the primitive variables are computed at each cell interface, one from values defined at, or to the left of the interface (pL, equation (4.61)), and the other using values defined at, or to the right of the interface (pR, equation (4.62)). Finally, using these last approximations, the characteristic structures of V, .F(pR) and F(pL) are calculated in stage C, enabling the computation of the Roe fluxes, F, using (4.56). Chapter 4. General Relativistic Hydrodynamics in Spherical Symmetry 67 where we have defined the following quantities: IC = (4.70) K/po - c^ Cl = vr-Vr±, (4.71) v± = 1 7 ^ 4 ' ( 4 - 7 2 ) l/tf - vrvr l/^-vrAr±' Al = \±/a + P/a. (4.74) In addition, we have introduced the sound speed, c s , in the above, defined by c 2 = 1/h (x + P/PQK.), with x — dP/dpo and K = dP/de. The particular values of the primitive variables used in order to compute these characteristic fields are 1/2 (p L + p R ) . The last ingredient we need in order to compute the Roe flux are the u>a; these are the jumps in the characteristic variables associated with the local Riemann problem, and are implicitly defined by q^i/2 - qi+i/2 = X^7?*- (4-75) CL Here (qR, qL) are the values of the conservative variables, which are calculated from the recon-structed primitive variables (p R ,p L ) , using (4.61-4.62). We first reconstruct the primitive vari-ables, then transform to conservative variables, since this approach leads to increased numerical stability relative to direct reconstruction of the conservative variables [79]. Calculation of the primitive variables The final piece of the algorithm for evolution of the discrete fluid quantities involves the calculation D,Sr,f From the of the primitive variables p = [po,v, P] from the conservative variables q = definition of the conservative variables (4.4)-(4.9), as well as the relation h = 1 + e + P/po, we can derive the following equation for the pressure f{P) = D (1 + e) W + P (W2 - 1) - D - r = 0. (4.76) Noting that W(P) = y/Z2/(Z2 — S2), where Z = (r + D + P), and assuming that the equation of state can be cast in the form e = e (p0, P) = (D/W, P), we see that, for given values of D, Sr and r , (4.76) becomes a non-linear equation for P. Given a good initial guess for P, we can find a solution to (4.76) using a Newton-Raphson method, for which we need to be able to compute the derivative df(P)/dP = f'(P). To this end, the following relations are useful: W'{P) = (1 - W2)/VZ2 - S2, (4.77) Chapter 4. General Relativistic Hydrodynamics in Spherical Symmetry 68 Furthermore, in order to calculate (4.78) we use dpo de_ dp D,T,Si de dpo dP p D dW 'W2 ~dP D dW ^w2 ~d~p D dW W2 ~dP + D,r,Si de_ dP de D,T,S{ dPo (-1) + P dp D,r,ST 9P° X , i de_ dP I PO dP_ de + Collecting results, we have f'(P)=D\l + e(^P W + DW D,r,Si W2 K K + W2 - 1 + 2PWW. For the specific equation of state considered here, P = (T — 1) poe, we have f'(P) = W D + 2PW-T W2 T T-1 1. r - I The following relations are also needed at various stages of the update algorithm: PW K = (r - i)po = (r - 1 ) D_ W' x = ( r - i ) e = PW c2 = D ' PWY (r - 1 ) D PWY r - i D (r - 1) + PWT Once the pressure is calculated the rest of the primitive variables can be computed using: 5, (4.79) (4.80) (4.81) (4.82) (4.83) (4.84) (4.85) .T + D + P ' 1 4 « r , PO D\/l — vrvr, P (4.86) (4.87) (4.88) (4.89) ( r - l ) p o The calculation vr using (4.86) can lead to non-physical velocities larger that 1. We note that in order to avoid this problem, an alternative expression to compute vr, which ensures vr < 1, was proposed in [78]. Chapter 4. General Relativistic Hydrodynamics in Spherical Symmetry 69 4.4 Results 4.4.1 New Variables As explained in Sec. 4.2, the variables, D, Sr and f that we have introduced in order to write the hydrodynamic equations in conservative form are not precisely those—D, Sr and T—which are proposed and used by other authors [66], [30]. Although use of either set allows us to write the fluid equations in conservation law form, we will now argue that the set D, Sr and f is a better choice when viewed in the context of providing sources for the geometrical constraint equations. Let us consider the initial data problem for our self-gravitating fluid. In spherical symmetry the gravitational field has no dynamical degrees of freedom, and as a consequence the initial state of the geometry is totally fixed by the state of the matter sources at t = 0. Specifically (and strictly for simplicity of presentation), if we consider time symmetric initial conditions (so that the resulting spacetime has a ( - » - t symmetry about t = 0), then the com-plexity of the differential equations governing our model is considerably reduced. In particular, the radial 3-velocity, vr, of the fluid must identically vanish, and this simplifies or eliminates many of the terms in the constraints that involve the fluid variables. Time symmetry also implies that the entire extrinsic curvature tensor vanishes, so the radial momentum constraint is trivially satisfied. The only non-trivial constraint is the Hamiltonian constraint, which written in terms of D and r (Sr vanishes due to the time symmetry) takes the form with boundary conditions given by (4.33). Here we note that the operator acting on tp is simply the radial piece of the Laplacian written in spherical-polar coordinates, taking into account the fact that ip is spherically symmetric. Given the condition vr = 0, we can view D and r as freely-specifiable quantities that fix the initial state of the fluid. For the sake of concreteness, we consider an initial profile for D given by a gaussian pulse and then compute r (0 , r ) , from the polytropic equation of state P = Kpv. This gives us the condition that P(r) = D{r)r, so that r(r ,0) = P{r)/(T - 1) = D(r)r/(T - 1). We investigate the behaviour of solutions of (4.90) using the following shooting method. For given values of A, ro and A r , and adopting the notation ip(r) = ip(0,r), we choose a trial value for ip(0) (the shooting parameter). Given the second initial condition, ip'(0) = 0 (which follows from regularity, and where a prime now denotes differentiation with respect to r) , (4.90) can be integrated radially outwards to some final radius r = r m a x . The value ip(Q) is then iteratively ^dr (r2dr*P) + 2TT (D + T) V>5 = 0, (4.90) D(0,r) = Aexp -{r-r0)2/A2 (4.91) Chapter 4. General Relativistic Hydrodynamics in Spherical Symmetry 70 refined until the value ^(rmaK) satisfies the discrete boundary condition given in (4.38) to some specified tolerance. For sufficiently small values of the amplitude parameter, A, we encountered no problems in finding solutions of (4.90) using this technique. However, for large values of A the situation was quite different: in such cases, in fact, there seemed to be no values of ip(0) that would produce solutions of (4.90) satisfying the outer boundary condition, lim r^oo ip(r) = 1. In Fig . 4.4 we show the estimates, lim r^oo A) for various values of A, and for T = 1.8. 1 1.1 1.2 1.3 1.4 1.5 Figure 4.4: This figure plots estimated values of tp(oo) as a function of the shooting parameter ip(0), for 10 different values of the amplitude A in the range 0.05 < A < 0.5 and spaced by AA = 0.05. A l l calculations have been done with ro = 0.5, A r = 0.1 (see (4.91)), r m a x = 1 and F = 1.8. ip(oo) is estimated from the assumption that l im^oo ip{r) = tp(oo) + C/r. These results clearly suggest that above a certain threshold amplitude (which in this case is approximately A ~ 0.13), there are no solutions of (4.90) that satisfy the outer boundary condition, lim,.-,;*, ip(r) — 1-We note that the solutions plotted in Fig . 4.4 do not correspond to particularly large values of Chapter 4. General Relativistic Hydrodynamics in Spherical Symmetry 71 the conformal factor tp(r), and that all solutions computed are smooth. We also note that the use of continuation (or homotopic) processes, which use information about the solution for a parameter value A as input (or initialization) for the solution for parameter value A + SA > A, did not help matters. In addition, we verified that the failure to find solutions of (4.90) was not a consequence of the discretization schemes used. In particular, we used several different solution methods and several distinct discretizations and found the same results in all cases. The results of the above experiment agree with an observation made by York [115], who analyzed an equation similar to (4.90). York argues that in order for the full non-linear equation to have a solution, the linearization of the equation with respect to V + Sip, must also have a solution. The claim is that, in general, the linearized equation will not have a solution satisfying Sip —> 0 as r —> oo. The reasoning that York follows involves analysis of the homogeneous equation (i.e the linear equation without sources) which he shows does not satisfy a maximum principle (as one generally wants for elliptic equations), but instead admits solutions have no solution satisfying the boundary conditions at infinity, when the analogue of D + r is freely specified. The argument that the homogeneous equation does not have solutions that asymptotically tend to zero (applied to our Hamiltonian constraint) is based on the fact that (D + r ) ip4 is positive-definite. The sign of this term can be changed by a suitable conformal rescaling of r and D, i.e. by choosing to freely specify conformally related functions D and f defined by D = Dipn and f = ripn for some integer n > 1. In terms of these new variables the linearization of the Hamiltonian constraint is By York's argument, if the second term on the left hand side of the above equation is negative-definite, then the theory of elliptic equations tells us that the solution of the equation exists and is unique. The choice n = 6 in our definitions of D and f thus not only allows us to demonstrate linear stability of the Hamiltonian constraint, it also absorbs the factor ip6 that originates from the determinant of the 4-metric, and which appears in the fluid equations of motion. (Here we note that a factor of tp6 was also introduced in the definitions of Sr and P in order to simplify the equations of motion). —8r (r2dr) + IOTT (D + r) ip4 Sip = -2mp5 (SD + ST) (4.92) that tend to be oscillatory as r —+ oo. York therefore concludes that the full equation will generally (4.93) Chapter 4. General Relativistic Hydrodynamics in Spherical Symmetry 72 4.4.2 Evolution of T O V solution The Tolman-Oppenheimer-Volkoff (TOV) solutions [80], [102], and [103] are the static solutions of the spherically symmetric equations of general relativistic hydrodynamics. These solutions were initially studied by Tolman [102] and [103] and then generalized by Oppenheimer and Volkoff [80] in order to describe neutron stars (stars supported against gravitational collapse by the degeneracy of neutrons). To construct these solutions we adopt a poly tropic equation of state P = Kpl, (4.94) which can be considered as an equation for non-interacting degenerate matter [94]. This equation can be seen as a particular case of the ideal gas equation P = {J? — l)poe in the limit of zero temperature [79]. Let us remind the reader that one of the main uses of the spherically symmetric code explained in this chapter is as a calibrator of our axisymmetric implementation. Ultimately we want to use the axisymmetric code to study the critical collapse of rotating stars. The corresponding spherically symmetric problem is the critical collapse of perturbed T O V solutions, and was studied by Noble [79]. Since the lifetimes of near critical solutions can grow without bound as the critical limit is approached, it is thus crucial that we are able to evolve T O V solutions for long physical times. It is particularly convenient to compute T O V solutions in Schwarzschild-like coordinates (t,r), where the time-independent, spherically-symmetric metric takes the form ds2 = - exp (20(f)) dP ^P-^j 1 df2 + f W (4.95) Using the polytropic equation of state (4.94), the resulting equations for the metric coefficients and pressure are ^? = 4TT?2P, (4.96) df 1 ^ (4.97) df p + P df dP (p + P)(rn + 4irf3P) (4.98) df f (f — 2m) The equations give rise to a continuous family of solutions which can be parametrized by the central pressure P(0). Choosing a particular value for P(0), and with the additional initial con-ditions TO(0) = 0 (regularity) and c4(0) = 0 (normalization of time parameter to central proper time), the set of ODEs (4.96) is integrated outwards from f = 0 using the L S O D A integration package [84]. After the metric coefficients m(f), <f>(f) and the pressure P(f) have been computed Chapter 4. General Relativistic Hydrodynamics in Spherical Symmetry 73 in the Schwarzschild-like coordinates, we must perform a coordinate transformation to the maxi-mal/isotropic coordinates, (t,r), used in our evolution code. Since the T O V solutions are static, we have t(i,f) = i and r(i,f) — r(r). We can determine r(f) by integrating ^ = / 1 _ 2 m ( r ) \ - 1 / 2 r df \ f J f This equation can be obtained by first comparing the angular parts of the two metrics, which yields the relationship tp4(r)r2 = f 2 , then comparing the radial parts, 1 _ 2 m { > ) \ 1 ^ = V>2(r)dr, (4.100) and combining the two results. We now discuss the integration of (4.99). Considering the limit f —* 0 of the equation, and taking into account that m(f) ~ 0 ( r 3 ) by regularity, we obtain the following equation valid for r - » 0 : £ = I (4-101) dr r Therefore the behaviour of r for small values of f is given by r = Ar, with A a constant that will be fixed by demanding that the solution tend to the Schwarzschild form as r —> oo. We thus first integrate the following outwards from r = 0: ur — = 1, if r = 0, (4.102) dr dr = A _ 2 m ( f ) \ - 1 / 2 r ; dr \ r J r which amounts to choosing A = 1. Once this integration is complete, we rescale the solution, r(f), so that it satisfies the correct asymptotic boundary condition by exploiting the fact that if r(f) is a solution to (4.99), then kr(f) (with k constant) is also a solution. In particular, we set r ( r m a x ) = % 1 - + J l - 2 m ( r m a x ) I • ' (4.104) This condition is computed by comparing the metric coefficients of Schwarzschild in the two coor-dinate systems considered, i.e. by comparing d s2 = _ ( 1 _ ^ ) d t 2 + ( i _ ? M i y 1 d r i + r 2 d t f > (4.105) r I \ r 2 / nr \ —2 / 7, r \ 4 and ds2 = -(l-^Ml^j ( l + ^r) d * 2 + ' ( 1 + ^ r ) ( ^ ~ 2 + r W ) , (4.106) where Ms is the A D M mass of the star. Chapter 4. General Relativistic Hydrodynamics in Spherical Symmetry 74 Finally, we note that once r = r(f) is found, tp(r) is calculated using (4.107) As mentioned above, the different T O V solutions can be parametrized by the value of the pressure at the centre of the star or, equivalently, by the central density. Fig . 4.5 shows plots of solution curves for T O V initial conditions. Specifically, the left plot shows the A D M masses of the stars obtained for T = 5/3 (which corresponds to a non-relativistic degenerate fermi gas) and K = 1 as a function of the base-10 logarithm of the stellar radius (we note that the radius of the star, f*, has been defined somewhat arbitrarily by P(f*) < 10~ 1 0 ) . The right plot shows the A D M mass as a function of the base-10 logarithm of the central density. Such plots have been calculated by many authors [44], [79]. Our equilibrium curves qualitatively agree with the ones calculated by [79], and we note that direct comparison is not possible since we are using a different coordinate system. : • i • • • r-p 1 ' i 1 ' 1 i 1 "1 b i I i i i I I i i i I i i i I , i _ d h I i i i i I i i i U I i i i i I i i i i L 0 0.2 0.4 0.6 - 3 -2 -1 0 1 L o g I 0 ( r . ) L o g 1 0 ( p ( 0 ) ) Figure 4.5: These plots show solution curves for T O V initial conditions, computed with T = 5/3 and K = 1. The figure on the left plots the A D M mass of each star against the base-10 logarithm of its radius, f*. The plot on the right shows the A D M mass of the stars as a function of the base-10 logarithm of the central density. In each case, the vertical line indicates a transition between stable and unstable stars. We now proceed to a discussion of typical results computed using our dynamical spherically symmetric hydro code, using T O V configurations as initial data. First, in order to demonstrate Chapter 4. General Relativistic Hydrodynamics in Spherical Symmetry 75 convergence of our code, Fig. 4.6 shows the time evolution of the rescaled central pressure, P(t, 0), from calculations at three different mesh spacings (resolutions). The initial condition for each of the three runs is a T O V solution with a central pressure P(0) = 0.002157, and calculated using a polytropic equation of state with K = 0.1 and F = 5/3. In order to perform the evolution we choose an integration domain with outer radius at r m a x = 5 (approximately 8 times larger than the radius of the star) which we discretize with Nr = 1025, 2049, 4097 points. Note that we use At = XAr, with A = 0.2 held fixed as the spatial resolution is varied, so that each run is characterized by a single discretization scale, Ar. This figure provides evidence that the evolution is second order accurate since the least-squares slope, m ~ dP(t,0)/dt, which should be zero in the continuum limit, is apparently 0(Ar2). Although the code is convergent, it is also quite dissipative (a result of the particular Godunov scheme that we are using), which leads to dispersal of the stars after some time (see Fig. 4.6). This effect seems to be more acute in the current maximal/isotropic coordinates than in the polar/areal coordinate system used, for example, in [79]. In order to decrease the amount of dissipation introduced by the update algorithm, we replace the original equation for the numerical flux (4.56) by Fi = \(F (PR) + F (p L ) ) - I £ \ \ a \ w a t i a , (4.108) where e is a tunable parameter. Fig. 4.7 shows the evolution of P(t, 0) for different values of e. Note that for e = 1 we recover the usual Roe approximation for the numerical flux, whereas for e = 0 the discretization corresponds to a particular finite-difference approximation of the hydrodynamic equations. From the figure we can see that a decrease in e leads to a decrease in the slope of dP(t,0)/dt, but that it also makes the solution increasingly irregular. It is thus evident that at least some amount of the significant dissipativity exhibited by the code can be attributed to our particular computation of the numerical flux. In this chapter we have described a spherical symmetric code to solve the fully coupled equations of general relativistic hydrodynamics. This implementation has not only served as a preliminary step in the construction of the axisymmetric code described in the following chapter, it has allowed us to identify a new set of dynamical variables | . D , Sr, f j with which to describe the fluid. These variables give rise to a well-posed elliptic problem for the constraint equations in the particular coordinate system considered. Moreover, we have also learned that our numerical method is so dissipative that we have difficulty achieving long-term evolution of static, stable, T O V solutions. In order to solve this latter problem, a new numerical implementation that uses adaptive mesh refinement—so that increased grid resolution can be automatically increased where needed—is being developed. This too will be briefly discussed in the next chapter. Chapter 4. General Relativistic Hydrodynamics in Spherical Symmetry 76 Figure 4.6: This plot shows the evolution of the central value of the rescaled pressure, P(t, 0), using fixed initial data, and three separate finite-difference resolutions: Nr = 1025, 2049 and 4097, with A = At/Ar = 0.2. The initial data are computed with a polytropic equation of state (K = 0.1 and T = 5/3) and a central pressure P(0) = 0.002157. The outer boundary of the computational domain for these calculations is r m a x — 5, which is approximately 8 times larger than the radius of the star. We see clear evidence for convergence in the sense that the temporal variation in P(t, 0) decreases as the mesh spacing decreases. More quantitatively, the slopes mo, m i , m 2 (indicated in the figure legend) are tending to zero quadratically in the mesh spacing, as expected for a second order accurate code. Chapter 4. General Relativistic Hydrodynamics in Spherical Symmetry 77 0.0007 A J 1 1 | 1 1 1 | 1 1 1 | l l l \ \ \ s\ \ A i\ e = 0 , 2 5 CD II \ A \ \ A e = 0.5 S-i \ \ \ IP* 0.0006 0.0005 l i i 1 i i i 1 i l i 1 l l \ / 0 200 400 600 t Figure 4.7: This figure shows the evolution of the central value of the conformally rescaled pressure, P(t,0), for a sequence of evolutions starting from the same T O V initial data (P(0) = 0.00046, T = 5/3 and K = 1), but using varying amounts of numerical flux (see (4.108)). It is clear that addition of the numerical flux tends to make the code more dissipative (and hence more stable), and that this in turn can have a significant influence on the long-time evolution of the stars. It is also apparent that the numerical solution becomes increasingly irregular as the numerical flux tends to zero (again, consistent with the stabilizing property of the numerical flux). Chapter 5. Axisymmetric Hydrodynamics 78 Chapter 5 Axisymmetric Hydrodynamics 5.1 Introduction In this chapter we discuss a numerical code that evolves the hydrodynamic equations coupled to the Einstein field equations in axisymmetry. This work makes use of a previously developed code due to Choptuik, Hirschmann, Liebling and Pretorius [19] which solves the Einstein equations in axisymmetry with an optional scalar field matter source. The code described in [19] implements both fully- or partially-constrained evolution (i.e. at least some of the constraint equations are resolved at each time step to fix certain geometric quantities), and has been used to investigate the critical collapse of a massless scalar field in axisymmetry [20]. The work described here involved the following modifications to the existing code: (1) a perfect fluid was incorporated as a source for the geometric evolution and constraint equations, and (2) a routine to integrate the, hydrodynamic equations, using an extension of the numerical method described in the previous chapter was included. This chapter is organized as follows. In Sec. 5.2, a brief introduction to the 2+1+1 formalism is given; this includes a summary of the equations that are to be integrated. Our particular im-plementation is restricted to the case without rotation around the axis of symmetry. Nevertheless, because one of our long-term goals is the study of rotating configurations, we present, in Sec. 5.3, the equilibrium equations for determining initial data describing a rotating star. This system of equations has an integrability condition that we explain in some detail, following the treatment of Bonazzola et al. [8], but recasting their development in the framework of the 2+1+1 approach. Some explanation concerning the numerics is provided in Sec. 5.4. Finally, results, as well as future plans, are summarized in Sec. 5.5. 5.2 2+1+1 Formalism We begin with a brief summary of the 2+1+1 approach, and follow with a discussion of the equations of motion that result from the application of this approach in the context of perfect-fluid-containing spacetimes. The 2+1+1 formalism, originally introduced by R. Geroch [34], exploits Chapter 5. Axisymmetric Hydrodynamics 79 the fact that, by definition, there is a (spatial) axial Ki l l ing vector field in any axisymmetric spacetime. We remind the reader that the standard approach to numerical relativity involves the 3+1 decomposition summarized in Sec. 1.1, and used in previous chapters. The 3+1 approach involves slicing the spacetime into constant-^ spacelike hypersurfaces, and describes the 4-geometry in terms of the 3-geometry intrinsic to each hypersurface, as well as how each slice is embedded in the full spacetime. In the 2+1+1 approach, on the other hand, the spacetime is first decomposed on hypersurfaces orthogonal to the axial Ki l l ing vector field, £ a . In the remaining quotient spacetime, which has 2 spatial dimensions, as well as 1 temporal dimension, a space-plus-time (i.e. 2+1) decomposition is then performed, in complete analogy to the usual 3+1 split. From the 2+1 perspective, the net effects of the init ial decomposition with respect to the "symmetry dimension" are the appearance of additional fields (relative to what one would have for regular 2+1 relativistic gravity), and the appearance of additional source terms in the equations of motion. These fields and source terms are analogous to those that arise in a Kaluza-Klein decomposition (see [83] for a review of this type of decomposition). We now define <f> to be the coordinate associated with the axial symmetry, arid denote the Ki l l ing vector field by £ a = (d/d(j))a. Paralleling the development in Sec. 1.1, where we defined the 3-dimensional spatial metric (1.7) induced on the spacelike hypersurfaces S t , we now define the 3-dimensional space-time metric, jap, la/3 = Sa/3 - J j £ a £ / 3 , (5-1) induced on the 3-dimensional hypersurface orthogonal to £ a , where s = £ a £ Q is the norm of the Ki l l ing field (by symmetry, we need consider only a single 3-dimensional hypersurface in this case). We note that the relative minus sign between the two terms in the definition (5.1) is due to the fact that £ a is spacelike. The mixed form of 7Q /g is the operator that projects onto the hypersurface: 7°/3 = 6ap - 1 ^ , (5.2) (We again remind the reader of our index conventions: 4-dimensional spacetime tensor indices are denoted by greek letters, 3-dimensional ones are denoted by lower case latin characters, and finally 2-dimensional spatial ones are denoted with upper case latin letters.) We reemphasize that in this case the metric on the three dimensional hypersurface is Lorentzian, not Riemannian, as is the case for the standard 3 + 1 decomposition. Therefore the three dimensional indexes ... now take on both spatial and temporal values. The covariant derivative compatible with 7 Q /3 (7a&) is denoted by Da, and can be defined via projection of the 4-dimensional spacetime covariant derivative in the usual manner. Following [20] Chapter 5. Axisymmetric Hydrodynamics 80 the 2+1+1 Einstein equations can be written as: ™ - = - ^ - T C ' W - H - • < 5 3 ) D [ a w 6 ] = 8TTseabcTc<p, (5.4) <3>i*a6 = i l ? a Z ? b S + \s2ZacZhc + 8TT - \labTaa^j , (5.5) where we have introduced the twist vector, u>a coa = jea0XafZx°, (5.6) and the antisymmetric tensor, Zap Zap = da {^Zp} ~ dp (^&) • (5-7) As can be easily verified from the definition (5.6), w a / s 3 is divergence free: Z > ° ( ^ ) = 0 . (5.8) Let us point out at this point that if s = 1, in analogy with Kaluza's original work in reducing 4 + 1 dimensions to 3 + 1, £ a plays the role of a Maxwell field, and then Zab is simply an electromagnetic field strength tensor. Equation (5.3) can be viewed as an evolution equation for the scalar quantity s. Additionally, two of the three equations (5.4) provide evolution equations for uip and u>z while the third is a constraint on those components of the twist vector. Finally, the remaining 3-dimensional field equations (5.5) are to be further decomposed using a space-plus-time split. Performing this split, we can then write the spacetime metric as ds2 = (-a2 + ^ & + H u P 1 ^ dt2 + 2&dfcty + s2d<f>2 + 2 (Hup1 + ^ £ t £ j ) dtdxJ +2£IdxId<f> + (HU + ^iiij^j dx^x3, (5.9) where HAB is the 2-dimensional spatial metric, a is the lapse function, and PA is the 2-dimensional shift vector. In order to describe the perfect fluid, and in analogy with the spherically symmetric Chapter 5. Axisymmetric Hydrodynamics 81 case discussed in Chap. 4, we define the following variables: D = p0W, (5.10) SA = PohW2vA, (5.11) St = p0hW2V4>, • (5.12) E = p0hW2-P, (5.13)' r = E-D, (5.14) = 4 + -- (5-15) a i r a ^ = —V u o' ( 5 - 1 6 ) W = au° = (1 - vAvA -vtv^y1'2 . (5.17) In terms of these quantities the hydrodynamic equations take the form of conservation laws (see section 5.2.2). Note the non-standard, and original to our work, definitions of v$ and , which are made in order to get equations adapted to the 2+1+1 formalism. We now have all of the necessary elements needed to discuss the equations of motion for the geometric and hydrodynamic variables. 5.2.1 Geometry We fix our spatial coordinates by demanding that the 2-metric H A B be conformally flat, so that Hu =tf(t,p,z)fu, (5.18) where fjj = diag [1,1] is the Euclidean metric in cylindrical coordinates (p, z). The time coordinate is fixed by requiring that the trace of the three dimensional extrinsic curvature, ^ K a a , be zero, and where we note that ( 3 ) K = ( 3 ) R i . = _ n i g . ( m s ) + ( 2 ) ^ / 7 ( 5 1 9 ) Here, nl are the components of a unit vector field orthogonal to the constant-t spacelike hypersur-faces. Equations (5.18-5.19), along with appropiate boundary conditions completely determine a and /3 7 and thus exhaust the 2+1+1 coordinate freedom. In order to improve near-axis (p —> 0) regularity of our numerical solutions, we introduce a (roughly) dynamically conjugate pair of variables, a and Cl, defined as follows: pa = l o g ( ^ ) , (5.20) pCi = -2K/-K/. (5.21) Chapter 5. Axisymmetric Hydrodynamics 82 These are evolved in lieu of s and its dynamical conjugate. In addition, we make the following definitions: CT po, (5.22) £* = tfVi/p\ (5.23) UJp = Up/p2, (5.24) Uz = Uz/P3, (5.25) b = 4>6D, (5.26) Si = (5.27) = (5.28) T = (5.29) P = *p6p. (5.30) In the equations presented below, we generally use the variable cr in order to minimize the com-plexity of expressions; the actual variable used in our numerical scheme, however, is o. The variables £5*, wp and wz were originally introduced in the study [21] of a complex scalar field, ^(t, p, z, <f>), endowed with angular momentum via an ansatz *(t,p,z,4>) = <S>(t,p,z)eim*, (5.31) for the specific case m — 1, i.e. with * ( t ) p , z , 0 ) = $(t ) / O ,2 : )e i *. (5.32) The specific powers of p chosen in the definitions (5.23-5.25) facilitate the construction of finite difference schemes whose solutions have good regularity properties. In particular, for the case of ansatz (5.32), it can be shown that the leading order behavior of the twist components is lima;* = p3f(t, z) + 0{p5), (5.33) p-*0 lim up = p3g{t,z) + 0(p4), (5.34) p—>0 \imuz = pAh(t,z) + 0(p5). (5.35) p->0 These expansions then imply that w', w p and u>z are all O(p) as p —> 0—this behavior is more readily maintained in the finite difference domain than 0(p2), 0(p3) etc. A n important point is that these regularity conditions—again, derived for the case of coupling to a scalar field with the ansatz (5.32)—are more stringent than the most general condition in the Chapter 5. Axisymmetric Hydrodynamics 83 fluid case. Forcing the twist vector to have the above leading order behavior forces the angular momentum to go to zero quadratically, i.e.: whereas in the general case we would have S$ = p(t, z) + 0(p2). This restrictive condition is chosen so that we can use equations which are as similar as possible to those described in [19], thereby minimizing the number of required modifications of the code. We also note that in (5.23) a factor of y/j, where 7 is the determinant of the 3-metric, has been factored out so that time derivatives of the lapse do not not appear in the equations of motion. Finally, expressions (5.26-5.30) define conformally rescaled conservative variables, paralleling the definitions made in the spherically symmetric case (see Sec. 4.4.1). We can now summarize the equations for the geometry. We have two equations which are derived from the momentum constraints, and which govern the shift vector components, (3P and lim 5^  = p2q(t,z) + 0(p4) (5.36) p': \ p % + p % + \p%z + \ [I (P% - pp,o) <*,? - {P:„+p%) «,*] + {P%+p%) + i [4 - /3* - paCl) V l P + 6 + /?*) t*] ~ \p* [3fi<r,„ + n,p]_ - | a f l -16TT aSp apuiptf ~ ^ 1 2 e 4 , T = 0 (5.37) P:pp+\p%z - \p%z + \ \ {p% - p%) <*,z - K + p % ) «,P +, -P (P%+P:P)+°,P [P%+P%) + ^ [4 (2p% - 1p% - paCl) ^z + 6 (P% + P*P) tp] ~ \(x* [3ficr, * + fi,,] + 2 (/3* - (3?p) az -16TT aSz ap2wzw 0. (5.38) ip6 •0 1 2e 4 f f The Hamiltonian constraint provides an equation fortp: o 1 2 2 4 — iplPP + ip,zz + -il>,p + 'ip,zV,z+ip,pO-,p +2(otP) +2(<T ) Z ) + -ayP + 2a z z + 2atPP (5.39) Chapter 5. Axisymmetric Hydrodynamics 84 The slicing condition, which fixes the lapse, a, is derived from = 0 and d/dt (^K\) = 0: 1 2 . , , a,pP + a,zz + a z o > + a>cr p + - a p + — {ip,PatP + V,zOt,z) [ P 2 & a + - p d ] +1 [- | (^ ) a -1 ( w 2 - \ (P:P)2 - \ (/?f2)2 -P%P: P+^]} a 5 2 + 5? 52 — + (pWz+u>l)=0. (5.40) + D + P P2e^(f + D + P)\ 2 e ^ The spatial components of the twist, to A ; satisfy the following constraint, calculated from (5.4): 3pCJz + p2Cjz<p - pQPtZ + 167r5^eCT = 0. (5.41) In addition, the twist components satisfy the following evolution equations, also obtained from (5.4): IGTT^S^ + {Cjz<ppz + /3*pu>z) p + 3Cjzpz + + p?pu>„ + 2 ^ + -\Q-Kae' T Sd,Sz p(f + D + P) - (u>* a + w*a i P) p - 3w'a + 4 ijj^paCj* 5 ' (5.42) + -167r^5^ + /3* w, + r ^ z , , + - (CJP,ZPP + P%CJP) 1 ^ 4 16TT-SlpSpaea a ) ' a — (D*a p* ( f + D + P) Expression (5.8) provides an evolution equation for CJ1: i.-.t + 4 V;5 ~.t (5.43) = - 2 £ ' a f t p + (-u>2iZ + 3cZ>z<7i2) a - CJzaz + 2/3* w* + PZCJ]Z + /? p t^ p + - [(3w pa p - w P i P ) a - w p a p + 3 ^ / 3 ^ + - ^ - + + 6 ^ . From the definition of Cl we get an evolution equation for a: a,t = - {Ppo-,PP + PZO-,zP - p2aCl - pPpp + pP) . P The evolution equation for Cl is derived from (5.3): (5.44) (5.45) ~ + ^ +2 ^ a ^p p^ P 2 ^ 4 - 2 (<7 i Pp + 1) a + 2a p p , no-,zaip,z njJ,ppa fia(^,P)2 g2Cj2p2+Cj2p Cjt2pa + 8TT 10 5 2 |_p ( r + D + P ) p 3 ( f + D + P ) pijjS ip5p S2 ,2a (5.46) Chapter 5. Axisymmetric Hydrodynamics 85 Finally, an evolution equation for the conformal factor ip can be derived from ^Ktl = 0: (5.47) This completes the set of equations that we consider for the geometry. Note that not all the equa-tions are independent of each other. In particular, equations'(5.47) and (5.55) are both equations which can be used to update ip. In the case of fully constrained evolution we use (5.55) for that purpose, while (5.47) is used in a partially constrained evolution. 5.2.2 Fluids We now consider the axisymmetric hydrodynamic equations within the 2+1+1 formalism. As usual, the equations governing the evolution of the fluid may be derived from the conservation laws (1.36) and (1.37). Since we want to take advantage of H R S C methods (which have been proven to be very successful in the study of relativistic hydrodynamics) we need to cast the equations in conservation law form. In addition, the equations need to be adapted to the 2+1+1 decomposition of the spacetime. Since the complete derivation of the equations is somewhat lengthy, we simply state them here; details of the derivation are included in App. C. d_ dt (5.48) dt p (5.49) d_ dt a Hi" (3Z a + P5Z, (5.50) (3p + — <^ ae'S* oz (3Z = 0. (5.51) Chapter 5. Axisymmetric Hydrodynamics 86 In the above equations we have introduced the following source terms, noting that they do not contain any explicit derivatives of the hydrodyhamical variables: Si Sf = {Ppe° + tfpe" (f + D + PJ (Sj - 2SJ) ') tf \ P G , J \ P / - - x 3 1 (f + D + P)tf dp" dpz\ pe" fx da - da dp dz ^ [ S D ~ + S.. dp dz pe SpSz tf (f + D + P) Scj, ipApe° (f + D + P^ '-p2 Qae" + ppea 1 fe2 + ~sl) 3 [f + D + Pjtf 4>6 (5.5 2ape" -2a (f + D + P^jtP5 L 6 p a e l p \ d ^ + ptf pe (f + D + p) dpz lbp~dp+bz-dj 2 r.t ip f 9/9 + ape'Tp\^- - pea(f + D) da + Paea + Sz . 1 tfpe'7 I tf -p LO (f + D + P) P{f + D + P) • C^e-" ) S+, (5.53) pe ~ dpp - dp* bp~dz~+bz~dz~ a „2 ~ t tfpe" tf f + D + P (5.54) Exploiting the fact that the coupled Einstein/hydrodynamical field equations presented in this and the previous section are over-determined due to the general covariance of the theory (equivalently the existence of the constraint equations), we can perform a non-trivial check of their overall Chapter 5. Axisymmetric Hydrodynamics 87 consistency. In particular we have demonstrated that the expression resulting from application of {92PP + d2zz) to the right hand side of (5.47) coincides with the expression obtained by taking the time derivative of the right hand side of the Hamiltonian constraint (5.55) written in the form 1p,pp + lf>,zz = 1 Ip ( 2 2 4 — i p t P - i p i Z o - i Z - ip,po-iP - — <2(<rtP) +2(o>) + + 2cr:ZZ + 2o,pp p ° I p +r {\P^ 2+1^  (P% - P%)(tr,) 2+1 w 2 + \ m2+mP+\ w 2 - \ m . (f + b) ti2p + p2ti2 i \ iP2 2e 4 C T V 8 2 V 1 2 e 4 C T f ' l ' Demonstration of the equivalence of the two expressions requires the use of most of the equations of motion. The calculation was carried out using the algebraic manipulation program, Maple [64]. 5.2.3 Boundary and Regularity conditions In order to solve the differential equations presented in the previous section, we need to provide boundary and regularity conditions. Here we restrict attention to the non-rotating case, since we have yet to incorporate angular momentum into our code. The boundary conditions at the outer edges of the (finite) computational domain are a combination of approximate Sommerfeld conditions and relations that follow from asymptotic flatness. More specifically, quantities governed by elliptic equations, i.e. a, ip, PA, satisfy boundary conditions of the form (see [87]) / - /oo + Pip + zftZ = 0, (5.56) where V'oo = «oo = 1 and /3^, = 0. These conditions can be derived from the known large-r fall-off of the metric components for the case of an asymptotically flat spacetime. The variables o and Q, which are radiative in nature, satisfy approximate Sommerfeld conditions at the outer boundary of the computational domain, i.e. asymptotically these variables are assumed to be of the form 9 — 9(t — r)/r, where r = y/p2 + z2. A l l of the dynamical fluid variables obey outflow boundary conditions. Near the z-axis, regularity dictates that a, ip, f3z, D, r and Sz are even functions of p, which implies that their first derivatives vanish at p = 0. On the other hand a, (3P and Sp are odd functions and thus go to zero at least linearly in p as p —> 0. Chapter 5. Axisymmetric Hydrodynamics 88 5.3 Initial Data for Rotating Stars Although the particular implementation to solve the relativistic hydrodynamic equations presented in this thesis does not allow for rotating configurations, one of our long-tetm goals is to study the collapse of rotating stars. For that reason the equilibrium equations for self gravitating perfect fluid configurations with non-trivial angular momentum are presented in this section. A good introduction to this topic is the review paper by Stergioulas [97], where the basic formalism and the main implementations to date are discussed. Here, we will derive the equilibrium equations directly from the dynamical equations introduced in the Sees. 5.2.1 and 5.2.2. In addition to assuming that the time derivatives of the metric and the fluid quantities vanish, we impose the following conditions: (3A = SA = & = 0 = to1 = 0. These choices are made to eliminate any terms that "source" time derivatives in the evolution equations. If upon setting the time derivatives of the dynamical variables to be zero at the initial time, we find that the evolution equations imply that the time derivatives vanish for all times, then the spacetime is stationary and our coordinates are adapted to the time-translational symmetry, i.e. to the timelike Ki l l ing vector field. We thus believe that by imposing the above conditions we have not restricted the type of axisymmetric stationary solutions that can be obtained. From (5.50) we then get the following two equations: V>2 * (5.57) + 6 P (5.58) Equation (5.46) gives the following condition: (5.59) Chapter 5. Axisymmetric Hydrodynamics 89 The equation for the evolution of Co1 (5.44), yields: (—&z,z + 3wzcr z) a - £ z a 2 + - [(3wp0- j P - coPtP) a - LOpatP} 6a yj p tp 0, while the constraint equation for Co A gives: 3pioz + p2Coz,p - pCoPtz + lQirS^e" = 0. (5.60) (5.61) The Hamiltonian constraint, 1p,pp + ->P,zz + ~i>,p + 1p,zO-,z + ll>,p°~,p + 2 (CT pf + 2 (cr, 2) 2 + -aiP + 2CT,22 + 2CTIPP + 1 6 7 T (f + D) Co2 + o2Co2 _j_ p ~ F z _ Q ip2 2e4<*ip8 and the slicing condition, (5.62) 1 a PP + ot,zz + celZo-,z + a<patP + - a i P + — [4>,patP + ^, za, 2) - 4 > ( f + 3 ^ + ^ ) - 4 > p '• v s2 p2e2a (f + JJ + p ) 2 e 4 ^ must also be satisfied. Equations (5.57-5.63) are the equations of general relativistic, hydrostatic equilibrium for a rotating perfect fluid. Following Stergioulas [97], we introduce a new function u(p,z) such that ^ = (s 2 ,w,0,0) . We note that since our coordinates are adapted to the timelike and axial Ki l l ing vector fields, this function has an invariant geometric meaning. Specifically, if na is the timelike Ki l l ing vector field, then we have ^rf=^t=u(p,z), (5.64) where to is a spacetime scalar. Moreover since £t = gt<p—again by our choice of coordinates—co(p, z) is, asymptotically, proportional to the angular momentum of the spacetime. In terms of this new function we can then write the spatial components of the twist vector as 2w 3 / 2 (i>2e°\ tfe pa \ T/W 2w 3 / 2 (ptfe° 1,2 Oo r pa ALO—,—h 2coa z — to p3a to I p p2a 2tfea \LO i>tP LotP P i> 2 (5.65) (5.66) These expressions for the twist vector components automatically satisfy equation (5.60). They also Chapter 5. Axisymmetric Hydrodynamics 90 yield an elliptic equation for LO that can be derived from (5.61): (paip2ea) (aip2ea) ^,pp + W , z z - 2 paip2e 2ecr V,pp + °~,zz + (<J,pf + (<7,z)2 aip2e" ' to -to + 2 la z<7 * + a„On\ h 2 /9 a /9a - 4 [Vvp + + 2iptPo\p + 2ip:Zo:Z] ^  - + 4 [V'.pQ.p + 4>,zCt:Z\ _4JtPf+ (tzf _ = Q (5_67) Defining w = W/O, where u is yet another function that satisfies l i m p ^ o u(p, z) = pu\(z) + 0(p3), we get Wpp +w, 2 Z + ( - ) -PJ o atp2ea p ai/j2ea -u, - 2 V,PP + o,zz + ( c p ) 2 + (a, 2) 2 --a pu + 2 [a 2 C T z + a„o-iP] - + 2 — 4 [^ .pp + 4>,zz + ^ , P ^ , P + \zO\z] — - 8 ^ -p a a p "'' +4 [^,Pa,p + i/>,za,z] — - 4 - ^ (V'.P)2 + (V^f 16 pip"2 = 0. V' "41 p (5.68) We now have a complete set of equilibrium equations. Naively, we might.think that if we provide the form of the angular momentum function S$ we could integrate equations (5.57), (5.58), (5.59), (5.62), (5.63) and (5.68), and obtain an equilibrium configuration. In practice, however, not all functions S<f, would produce a solution. In the next section we investigate the possible forms for the rotation function. 5.3.1 Integrability Condition The equations presented in the previous section have an integrability condition. In particular we have two equations that could be used to compute the pressure—namely (5.57) and (5.58). For our system of equations to be consistent, it is clear that we must obtain the same result irrespective of which of the two is used. This condition restricts the allowable functional form of S^. In theory, we should be able to obtain the appropriate integrability condition, i.e. the condition on S^, by demanding that the derivative with respect to z of the right hand side of (5.57) agrees with the derivative with respect to p of the right hand side of (5.58). This specific procedure is somewhat difficult to carry through, and has not proven to be very illuminating. Instead we will follow Bonazzola et al. [8], but write their results in terms of the variables used in the other sections of this chapter. The four velocity for the case of a stationary star can be written as u» = (u^,W/a,0,-0), (5.69) Chapter 5. Axisymmetric Hydrodynamics 91 where W is defined by equation (5.17) and a is the lapse function. Starting from U * = ^ « " = 6 f + - V = W ^ , (5.70) and assuming that £ t = w (as was done previously), can be expressed in terms of the conservative variables as follows: Again, following [8], this allows us to define a quantity Cl* (not to be confused with Cl): a S s u* s 2 \ p 0 h W 2 J p 2 ^ e 2 ° {f + b + p) (5.72) In terms of Cl*, equations (5.57) and (5.58) take the compact form ^ ^PA . _ £ d + _ " ^ n v (5.73) ( f + L> + P ) ' a W a ( f + D + P ) The above expressions can be further simplified by taking into account the following relationship, ^w 2 = i = j _ = i . ( 5 7 4 ) ( V + £ ) + p ) p 0 ( l + e) + P p 0 h PH+-P' where we have used definitions (5.26-5.30) and the definition of pn introduced in equation (1.33). At least formally, we can define the following functions: M - ITJTTP' (5-75) T _ - ^ (5.7.) a ( f + D + P j This allows us to write (5.73) in the following compact form: (H + \ n a - \ n W ) A = C l * ; A . ' (5.77) It is now easy to see that the demand that the mixed derivative of the left hand side of this expression give the same answer independently of the order of differentiation can be expressed as a condition on the variables appearing in the right hand side, namely F z f l * p - F p C l * z = 0. (5.78) In [8] it is argued that the left hand side of (5.78) can be viewed as the Jacobian of the transfor-mation between (p, z) to [T, Cl*). The fact that the Jacobian is zero then implies that there exists a function $ that relates T and Cl*, i.e. that: $( .F , f i*) = 0. ' (5.79) In this case two possibilities exist: Chapter 5. Axisymmetric Hydrodynamics 92 • If &tjr = 0, then <!?(£}*) = 0 expresses the constancy of fi*; this case is called rigid rotation. • If j 0, then there is a relationship between T and $ that can be written as T = T (fi*); this case is called differential rotation. In the case of differential rotation we can calculate the value of fi* in different parts of the star by fixing the form of ^(fi*) and then solving F ( f i * ) + - 5.80 a2s2 -(s2fi* +OJ)2 for fi*. In both cases (5.77) has a first integral. For the case of rigid rotation we have . H ( p , z ) + l n ( £ ) ( p , z ) = k , (5.81) while for differential rotation the integral is H (P,z) + ln (—J (p, z) + / f(fi*)dfi* = k. (5.82) In both instances A; is a constant. Finally let us point out that the Newtonian limit (pfi* -C 1) of case (5.80) takes the form F(p*) = - p 2 S l * , (5.83) which implies that fl* = fi,* (p2)—hence the terminology "differential rotation". 5.4 Numer i c s In this section we briefly describe several aspects of our numerical implementation. We start with a discussion of the treatment of the geometric equations, follow with a description of the integration' of the hydrodynamic equations, then end with a discussion of issues arising from the coupling of the two systems of PDEs . To date we have only implemented the case of no rotation and therefore will restrict attention to the case = CJ 1 = = 0 for the remainder of the thesis. The code has been implemented using R N P L (Rapid Numerical Prototyping Language) [65] with some specific routines written in Fortran 77. The numerical approximation used for the geometric equations is explained in detail in [19] and [87]. It is based on second order centred finite difference approximations on a uniform grid in the (p, z) plane. More specifically the geometry is computed on a grid of points (see Figure 5.1) denoted .by (pi, Zj) where i = 1 , 2 , N p , j = 1,2,.., such that pj+i = pi + Ap and Zj+i = Zj + Az. Here Ap and Az constants and p\ = 0, p^p = Pmax and z\ — zm-m, ZNZ = z m a x - In practice, we always compute with p m a x = — zm-m = zmax, so that Az = Ap = h (which implies Nz — 1 = 2(NP — 1)). In Chapter 5. Axisymmetric Hydrodynamics 93 addition we choose a discrete time step A t = Xh, where the so-called Courant factor A—which must generally satisfy A < l/y/2 for stability of our numerical scheme—is held constant when we vary the spatial discretization scale, keeping the initial data fixed. Thus, the entire numerical scheme is characterized by the single discretization scale, h, which facilitates convergence testing of the results. The fluid is solved using a finite volume approximation which considers finite difference cells Ci+i/2j+i/2 centred &t (pi+i/2, Zj+i/2)-As we have discussed in Sec. 5.2.1, when we perform a fully constrained evolution the only geometric evolution equations are those for cr and fi, (5.45) and (5.46) respectively. These are discretized using a Crank-Nicholson scheme and second order centred differences for the spatial derivatives as in Chap. 2. We also again apply Kreiss-Oliger style dissipation [52] in order to damp high frequency components which cannot be properly represented at any given resolution, and which tend to result in instabilities in the code, especially near the 2-axis. When we perform partially constrained evolution, we update ip using (5.47) rather than via the Hamiltonian constraint. This evolution equation is also discretized using a Crank-Nicholson scheme. Along with the evolution equations, discrete versions of the constraints (5.37-5.55) and the slicing condition (5.40) must be solved at each time step. These equations, which we assume are always elliptic, are discretized using second order centred finite difference approximations of the spatial derivatives. The resulting discrete systems are solved using an FAS (Full Approximation Storage) multigrid algorithm [9] to determine the advanced values of the discrete lapse and shift components, and, in the case of fully constrained evolution, the discrete conformal factor. The choice of multigrid is motivated by the fact that it is unique among general methods for the solution of finite-differenced non-linear elliptic systems in being able to produce a solution in O(N) time, where N is the number of points in the spatial discretization (N = NZNP in our case). Multigrid algorithms are based on the observation that the decades-old technique of relaxation, while not a very efficient solver of discrete elliptic equations, is often a very efficient smoother of the equations. In particular, for the purposes of illustration, we consider a (scalar) elliptic problem written in the form L u = / , (5.84) where L is some elliptic differential operator (possibly nonlinear), u is the continuum solution, and / is a source function. (Note that u and / are functions of some number of independent variables; the multigrid technique can be applied in any number of spatial dimensions. Also, although the treatment of boundary conditions is an important issue, we will not consider it here since we only wish to illustrate the key ideas underlying the method.) We.discretize (5.84) at some grid Chapter 5, Axisymmetric Hydrodynamics 94 resolution, h, as Lhuh = fh^ (5>85) Here, Lh is the discrete (difference) operator that approximates L, while uh and fh are the discrete solution and source, respectively. We note that at a given resolution, (5.85) will comprise N algebraic equations, where N is the number of grid points, and that we will generally be able to naturally associate one unknown (component of uh) and one equation with any given grid point. We now consider the specific case of Gauss-Seidel relaxation, which, as with all relaxation algorithms, proceeds iteratively. Each iteration, or relaxation sweep, consists of a visit to each grid point (in some prescribed order), where the value of the discrete unknown associated with that point is modified so that, instantaneously, its equation is satisfied. Denote by uh the approximate solution of (5.85) at any stage of the iteration (so that, assuming that the relaxation converges, uh —» uh in the limit of an infinite number of sweeps). Then we define the residual, fh by fh E E Lhuh - fh , (5.86) and the solution error eh by e h = u h - u. (5.87) By the observation noted above, (Gauss-Seidel) relaxation is generally not very efficient at anni-hilating fh and eh, but it is effective at smoothing (i.e. annihilating high frequency components) those quantities. This brings us to another key ingredient of multigrid algorithms—from which their name derives—and that is the use of a sequence of ever-coarser grids that are employed to accelerate the solution process. In particular, once the residual (5.86) has been sufficiently smoothed (this generally requires only a few sweeps, typically 2-4), we can sensibly transfer the discrete problem to a coarser grid, which, from considerations of optimal efficiency, as well as programming convenience, is almost always characterized by a discretization scale 2h. Specifically, on the coarse grid we pose the problem L2hu2h = f2h + fh2h, (5.88) where r2h is computed from r2h = L 2 h I 2 h u h - I 2 h L h u h , . (5.89) and I%h is a so-called restriction operator that transfers a fine grid function to the coarse grid. It can be shown that f2h is an estimate of the truncation error, T%h, of the solution of the coarse grid discrete system L2hu2h = f2h relative to the fine grid problem (5.85). r2h has the property that the solution, u2h, of T2h'2h c2h i J2h (r nn\ L u = } + rh , (5.90) Chapter 5. Axisymmetric Hydrodynamics 95 is precisely the solution of the fine grid problem (5.85), transferred to the coarse grid, i.e. u2h = I2hhuh . (5.91) Thus, r%h (and thus also the approximation f2h) corrects the right hand side of a coarse grid problem, allowing—in principle—solutions with fine grid accuracy to be computed on the coarse grid. The coarse problem (5.88) is significantly less costly to solve than the original fine grid equa-tions (5.85), since in d dimensions it involves N/2d unknowns. More importantly, however, we can apply the above strategy recursively. That is, we perform relaxation sweeps of (5.88) until the corresponding residuals and solution errors are smooth, then transfer to a grid with resolution 4h etc., until we eventually are using a grid that has so few points, that actually solving the discrete equations is very cheap, even if we use a direct method (i.e. simultaneous solution of all of the algebraic equations) rather than relaxation. Once the coarsest-level problem is solved, we begin to work our way back to the fine grid, via a sequence of coarse-to-fine grid transfers. In particular, having (approximately) solved a 2h problem, yielding u2h, we update the finer-grid unknown, uh using uh:=I%h{u2h-I2huh) (5.92) where is a so-called prolongation operator that transfers coarse grid functions to a fine grid. After each of these updates, we again apply a few relaxation sweeps to (5.85) i i i order to ki l l any high frequency components produced by the prolongation operation. Once we are back on the fine grid, we will have completed what is known as a V-cycle, and will generally find that the norm of the residuals and solution errors will have been reduced by some constant factor. Additional V-cycles can then be applied as needed in order to drive the residual below some convergence threshold. This completes the description of the basic operation of a multigrid method. The particular choice of restriction and propagation operators, and some other specifics of the algorithm used in the axisymmetric code are explained in [19],[87]. We now move on to the hydrodynamic equations, (5.48-5.50), which are integrated via a finite volume approximation similar to the one described in Chap. 4. Note that these equations are of the type: | (e'pq) a + ~ (ape'F") + ^ (ae"F*) = S. (5.93) In order to derive a finite volume approximation, the equations are integrated over a control volume defined by C^f^j+1/2 = (tn,tn+1) x (pi,pi+i) x (zj,zj+i), as shown in Fig . 5.1. Chapter 5. Axisymmetric Hydrodynamics 96 In solving the hydrodynamical equations, we also decided to implement a generalization of the H R S C method explained in Sec. 4.3.2. Wi th this approach an approximate Riemann problem is solved at each cell interface to find an expression for the numerical flux. More specifically, the two fluxes Fp and Fz are evaluated via Roe approximations of the type given by expression (4.56). In turn, the Roe fluxes depend on the characteristic structure of equations (5.48-5.50) as described in App. D . The calculation of the fluid quantities at cell boundaries is performed using a one dimensional minmod reconstruction, also, described in 4.3.2. A z • • • • • • • • • • • [i+i • ) (ij •j+1) (i+l/2,j+l) (i+il 1 • j+1/2) (i+l/2,j+l/2) (i+l. • « » • • \ • 1 (i, i i+i< • 2) • • • \ < \ > • A p (i,j) (1+1/2J) o (i+l,j) Figure 5.1: This figure shows a detail of the projection onto the p-z plane of the finite-volume (cell-based) grid used in the numerical solution of the coupled Einstein/hydrodynamical equations in axisymmetry. Note that the fluid variables are computed at points denoted by (i + 1/2, j + 1/2). The p fluxes are computed at a location (i, j + 1/2), and the z fluxes at (i + 1/2, j). Geometric variables, on the,other hand, are computed at locations labelled with open circles. In order to obtain values of the fluid variables at these last locations, or to compute values of the geometric variables at the cell centres, second order (bi-linear) interpolation is used. As was the case in spherical symmetry, because of our choice of (topologically) cylindrical coordinates, one of the fluid equations contains a term which is explicitly divergent as p approaches zero. In particular one contribution to the p flux in (5.50), which governs Sp, contains a term P/p. This term is explicitly canceled, before discretization, by a corresponding term in the source Chapter 5. Axisymmetric Hydrodynamics 97 of (5.50). Finally, the way in which we perform the update of the coupled Einstein/hydrodynamical variables to advance the discrete solution from time t to time t + At is as follows. 1. Initialize the values of the geometry and the fluid at t+At (advanced values) to the previously computed values at time t (retarded values). 2. Improve the estimate of the advanced geometric variables governed by evolution equations by performing a Crank-Nicholson iteration. 3. Improve the estimate of the advanced constrained geometric variables by performing a multi-grid V-cycle on the system of constraint equations. 4. Improve the advanced fluid quantities, qn+1, using a two step method analogous to that described in Sec. 4.3.2. 5. Repeat steps 2-4 until the the norm of the change in qn+1, the norm of the residual of the constraints, and the norm of the residuals of the Crank-Nicholson iteration are below some tolerance. 5.5 Results In order to test the validity of our numerical implementation, we have performed various tests. Setting all the metric coefficients to be those of Minkowski spacetime—namely a = 1, (3A E E 0, a E E 0, tp E E 1—we have tested the part of the code that solves the hydrodynamic equations. Specifically, we have verified that the code is able to properly evolve discontinuous initial data, with surfaces of discontinuity defined by p = const., z = const, or p + z = const. We note this last case (oblique discontinuity) is particularly non-trivial, since our reconstruction to compute the numerical fluxes is one-dimensional at each stage (i.e. at each time step reconstruction is performed along lines of constant p and constant z independently). In addition, convergence tests of smooth initial data have shown second order convergence (see Fig . 5.2), in regions of the computational domain removed from the extrema of the fluid variables. Due to the specifics of the reconstruction procedure explained in Sec. 4.3, the code is, as expected, only first order accurate in the vicinity of local extrema of the solution. Still restricting to the case of flat spacetime, we have been able to evolve more complicated situations. One interesting scenario is the onset of a Kelvin-Helmholtz instability (Figs. 5.3 and 5.4). This type of instability occurs at a tangential discontinuity between two different states of a'fluid (or two different fluids) which are sliding parallel to the discontinuity [54]. A n example Chapter 5. Axisymmetric Hydrodynamics 98 0 1 2 3 4 0 - 2 4 6 8 t t Figure 5.2: Here we show £2 norms of the estimated solution errors computed at two different levels of discretizations (mesh spacings ho = 10/96 and h\ = ho/2 = 10/128). Time development of the estimated errors is shown for the pressure, Pi-ri/2j+i/2, 0 1 1 the left, and for the conformal factor, tpij, o n the right. The solution errors (at each level of discretization) are themselves computed via subtraction of values calculated using two discretization levels, i.e. as fk — fk+i, k = 0,1, where fk is the solution computed on a grid with cell spacing hk, and hk+i = hk/2. For a second order scheme, these differences should be quadratic in the mesh spacing, h, as h —+ 0, so that when we reduce the grid spacing by 2, the differences should decrease by roughly a factor of 4. Thus, we multiply the £2 norm of the finer-grid error estimate by 4 to show more explicitly that our code is second order convergent. We note that the results plotted here came from the evolutions described in more detail in Fig . 5.5 and accompanying text. We also note that at t ~ 3 the pressure becomes discontinuous and we can no longer expect Pj+i/2 j+1/2 to be second order convergent. On the other hand, the conformal factor ostensibly remains second-order convergent throughout the evolution (only the early stages are depicted here), which is a manifestation of the facts that: (a) the constraint equation, being elliptic, tends to smooth discontinuities in its sources; and (b) (related) the gravitational field (particularly the spherical or monopole piece) tends to be responsive to extended, rather than localized, distributions of matter/energy. However, in principle, if we could go to the h —> 0 limit, we would necessarily find that the conformal factor would only be twice differentiable at the locations of shocks, and hence we would observe "loss of convergence" in the geometric variables as well. Chapter 5, Axisymmetric Hydrodynamics 99 of this effect is the generation of waves on the surface of a body of water due to the wind. In our case we have simulated the contact discontinuity (meaning that the pressure remains constant across the discontinuity) of two coaxial cylinders of fluid. The fluid in the inner cylinder is initially at rest with constant density Dj, while the fluid in the inner cylinder is initially moving with a uniform velocity in the z direction—parallel to the discontinuity—and has a constant density Do, with Do < D[. (We have not simulated any situations where any of the fluid is rotating about the symmetry axis). In order to excite the instability, we deform the initial contact discontinuity, so that it has a small amount of curvature along the z-axis. Fig. 5.3 shows a schematic of the initial setup. Fig . 5.4 shows snapshots, from simulations using the same initial data, but two different reso-lutions, during the evolution of the instability. For the specific case shown here, the initial speed of the fluid is v = 1/2. The initial data has been specified using a polytropic equation of state P = Kp0r (like the one used in 4.4.2), with T = 3/2 throughout the fluid, but with distinct values of K chosen for the two cylinders. The initial values of D in the inner and outer cylinders are Di — 1.0 and Do = 0.115 respectively. The limits of the computational domain are pmax = 20, z m i n == — zmax = —10, and the two discretizations used to generate the results shown in Fig. 5.4 are Ap = Az = 20/256 (top figure) and Ap = Az = 20/512 (bottom figure) (for these particular simulation we have not adopted the usual choice pmax = zmax, explained in Sec. 5.4.) Full MPEG animations of this simulation can be viewed at [109]. In the continuum limit, the fluid should become totally turbulent due to the lack of viscosity in our hydrodynamical equations of motion. In practice, there is effective numerical viscosity, due to the fact that the equations are always solved using a finite mesh spacing, but also since, as we have seen in the spherically symmetric calculations, the Roe numerical flux adds dissipation. We find, however, that as the discretization scale h is reduced, our simulations always tend to develop features all of the way down to the particular mesh size being used. We have also performed some evolutions of the hydrodynamic equations fully coupled to gravity. Figs. 5.5 and 5.6 show results from the evolution of a weakly self-gravitating pulse of fluid. The initial conditions in this case are time symmetric with initial distributions of D and r given by 5(0, p,z) = ADexp{- [(p - 3) 2 + z 2 ] } , (5.94) 7(0, p,z) = A r e x p { - [(p - 3) 2 + z 2 ] } , (5.95) with AD = 1.0 x 1 0 - 2 and AT = 1.0 x 10~ 2 . This implies that the initial fluid distribution is toroidal in shape. For the runs described below, the computational domain had pmax = 10.0, - Zmax = — Z m i n = 10.0, NP — 96, NZ = 192, and we used an adiabatic index T = 3/2. The initial maximum value of the pressure is PQ = 4.97 x 1 0 - 3 . In Fig . 5.5 we can see that the Chapter 5. Axisymmetric Hydrodynamics 100 V •Ui J D° V J Figure 5.3: Geometry of the system used to study the Kelvin-Helmholtz instability. Initially, the system state describes two constant-density coaxial fluid cylinders, separated by a con-tact discontinuity. The inner, more dense cylinder (with D — Dj) of fluid starts at rest in the lab frame, while the outer, less dense cylinder (with D = Do) has a uniform velocity in the 2-direction v (along the axis of symmetry). Although not shown in this schematic, in order to trigger the instability the initial discontinuity is perturbed, and is thus not everywhere parallel to the z-axis (or, therefore, to the velocity of the fluid in the outer cylinder). Chapter 5. Axisymmetric Hydrodynamics 101 t=31.25 257 x 257 [0,20], [-10,10] t=31.25 513x513 [0,20], [-10,10] Figure 5.4: This figure shows t = const, snapshots (zooming in close to the surface of discontinuity) from simulations of the Kelvin-Helmholtz instability discussed in the text, starting from the same initial conditions, but using two different resolutions: h = Ar = Az = 20/256 (top) and h = Ar = Az = 20/512 (bottom). Variations in the grey scale highlight regions with different values of D (white to black, low to high density). We can see how the dynamics is highly dependent on the cell spacing, h, which, among other things, sets the amount of numerical viscosity in the discrete equations. Chapter 5. Axisymmetric Hydrodynamics 102 pulse immediately begins to disperse due to hydrostatic pressure. A rough estimation of the force density due to hydrostatic pressure gives P0/S « 5 x 1 0 - 3 where 5 = 1 is the width of the pulse given by (5.94) and (5.95). On the other hand the gravitational force can be estimated by MPo/((T - l)p2) « 8 x 10~ 4 where the mass M is approximated using the spherical symmetric formula M « 2(i/>o — 1), P ~ 3 is the position of the centre of the pulse and tp0 fa 1.1 is the initial maximum of the conformal factor shown in Fig . 5.6. This calculation shows that the force due to hydrostatic pressure is about an order of magnitude more important than the gravitational force at the initial time. The ingoing part of the fluid then propagates towards the axis and produces a shock wave— all the variables describing the fluid become discontinuous—that subsequently propagates outwards. For these initial conditions the fluid eventually completely disperses. F ig . 5.6 shows the evolution of the conformal factor from the same simulation. Note that the self-gravitation of the fluid remains very small throughout the evolution; for example, maxt:P:Z ip(t, p, z) ss 1.1, which is close to the Minkowskian value, tp = 1.0, and which is attained at t — 0. In Fig . 5.7 we show another time-symmetric simulation, with initial profiles for D and f again given by (5.94) and (5.95), but this time with AD — 5.0 x 10~ 2 and J 4 T = 5 x 1 0 - 2 , i.e. with initial amplitudes for the dynamical variables that are 5 times larger than for the weak field simulation just described (max P i 2 P(0, p, z) = 2.5 x 1 0 - 2 ) . Compared to the weak-field case, the evolution is now quite different. The first thing to notice is that, with respect to coordinate time, the fluid evolution unfolds more slowly. This is probably at least partly a coordinate effect, due to the spatio-temporal variation of the lapse function. For example, the minimum value of the lapse at the initial time is m i n P i 2 a(0, p, z) « 0.8 while at the end of the evolution we have m i n P ) 2 a(0, p, z) fa 0.3. The second thing that is apparent from Fig. 5.7 is the fact that the fluid pulse does not spread out due to its pressure in this case. Rather, the fluid is apparently compressed—producing a more compact configuration—and we believe that this is a direct consequence of the self-gravitational interaction. A t late stages in the evolution, this concentrated configuration moves closer to the axis, again, probably due to gravitational self-interaction. Shortly after the time shown in the final frame of the figure, the solver for the primitive variables fails (produces unphysical values) and the evolution cannot proceed. At the current time, the code is prone to numerical instabilities, tending to produce unphysical negative pressures and speeds > 1 in many situations. We believe that at least some of these problems are due to inadequate grid resolution. The maximum resolution that we can obtain is currently constrained by the memory size of the computers on which we run the simulations, as well as on the position of the outer boundaries of the computational domain. In particular, since the boundary conditions explained in Sec. 5.2.3 are only valid for large values of r ~ \J p2 + z2, if we do Chapter 5. Axisymmetric Hydrodynamics 103 t = 0.00 t = 2.50 t = 5.42 t= 17.08 t = 20.00 t = 22.92 Figure 5.5: Time evolution of a weakly self-gravitating pulse of perfect fluid. Here, the initial conditions are time symmetric (i.e. zero initial velocity), with initial profiles for D and T given by (5.94)-(5.95). The computational domain extends from p m i n = 0 to Pmax = 10 and from z m ; n = —10 to zmax = 10. The figure shows the time development of the pressure P(t, p, z) (z axis roughly horizontal in the plots) which has an initial maximum amplitude Po = 4.97 x 10~ 3 . The evolution shows how the pulse initially spreads out due to pressure forces. At t sa 5.0, the ingoing part of the pulse reaches the axis of symmetry and subsequently "self-reflects", producing a shock that marches outwards, and which is most clearly visible in the t = 14.17 and t = 17.08 frames. Chapter 5. Axisymmetric Hydrodynamics 104 t = 0.00 t = 2.50 t = 5.42 t = 8.33 t = 11.25 t = 14.17 t = 17.08 t = 20.00 t = 22.92 Figure 5.6: This figure shows the time evolution of the conformal factor ip(t, p, z) for the weakly self-gravitating pulse described in the text (see also Fig . 5.5). The maximal value, max t ) / 3 ] 2 xp(t, p, z), is about 1.1 and occurs at the initial time. Note that towards the end of the evolution, when the fluid is dispersing to large distances, the conformal factor tends to its Minkowski value, ip = 1.0. Chapter 5. Axisymmetric Hydrodynamics 105 not set the boundaries at large enough distances, we can also get unphysical results. For example the conformal factor, ip(r,t) may actually start growing, rather than falling off, as r —> r m a x . t = 0.0 t = 6.3 t = 8.3 t = 31.3 Figure 5.7: Evolution of a strongly self-gravitating pulse of fluid. Initial conditions are the same as that for the evolution showed in Figs. 5.5 and 5.6 but with amplitude parameters AD and A T two orders of magnitude larger than in the weak-field case. Note that the vertical scale here is very different from that of Fig . 5.5; in particular, in this case, max P ] 2 P(0, p, z) = 2.5 x 1 0 - 2 . We can see how the pressure evolution in this case is very different than that in the weak field simulation. Here the dynamics is apparently dominated by gravity. Instead of dispersing due to its pressure, the fluid compresses as a result of its self-gravity. At later times the torus of fluid collapses towards the axis, and shortly thereafter the calculation of the primitive variables produces unphysical results, and the simulation must be stopped. Despite our concerns about inadequate resolution, and the memory limitations we encounter at the current time when trying to compute at higher resolutions, we wish to stress that, ultimately, additional memory alone is very unlikely to solve the problems we face. For example, one specific long-term goal of this research is to study the critical gravitational collapse of perfect fluids in axisymmetry, including the case of rotating fluids. As explained in Sec. 1.2, In the corresponding spherically symmetric problem it was found that for certain choices of initial data families, the critical solutions are self-similar (see Evans and Coleman [29] and Neilsen and Choptuik for the study of the collapse of ultrarelativistic fluids, Hawke [46] for the collapse of the same type of Chapter 5. Axisymmetric Hydrodynamics 106 fluids in cosmological settings, and Noble [79] for the critical study of the collapse of perfect fluids and, in particular, critical behavior associated with T O V solutions). This strongly suggests that, in the axisymmetric case, evolution of similar types of initial data will also generate structures on all scales, accompanied by the development of large gradients in the dynamical variables. The only plausible way of attacking this problem in a computationally efficient manner is through the development of some sort of mesh refinement strategy to enable resolution to (locally) increase or decrease as necessary. To that end, we have begun the development of a spherical symmetric code that refines the mesh adaptively as required. We are considering different possibilities for the mesh refinement cri-terion. One possibility is to use a Richardson-type estimate (i.e. calculating differences of solutions computed at resolutions hk and kk+i = hk/2 as in the convergence tests discussed previously) for the truncation (solution) error; when this estimate exceeds a specified threshold we refine the mesh. This procedure has been used successfully by Hawke [46], in his study of relativistic fluids, and by LeVeque in [59]. On the other hand, due to the non-analyticity of the solutions (and the general non-smoothness of solutions computed using the H R S C methods employed in this thesis) these estimates are not very accurate and tend to be numerically irregular. This has motivated us to consider the correction to the numerical flux given by equation (4.56) as a measure of the need for refinement. Specifically when the following function: (the different quantities in the above expression are defined in Chap. 4) goes above threshold, we refine. Our current implementation uses a global time step, A i , on all levels, fixed by the finest level of spatial refinement: A i = A A r f i n e s t . The global time-stepping procedure simplifies the treatment of internal boundaries between grids with different mesh spacings, relative to methods such as that proposed by Berger and Collela [6]. We also point out that we have found the use of E N O (Essentially Non Oscillatory) [15] interpolation crucial in our initialization of (new) fine grid values from coarse grid values. We also plan to implement one or more approximations to the numerical flux, other than the Roe approximation described above. We are especially interested in seeing whether the use of an alternate solver helps with some of the strong-field problems that we have encountered in our collapse calculations. In particular, we are considering the use of the Marquina approximation [22] which has been compared with the Roe solver by Hawke [45]; Hawke found that the Marquina method almost always performed better than the Roe one. Another place for improvement is the reconstruction of the conservative variables at the cell interfaces, where we plan to implement interpolations other than minmod. (5.96) Chapter 5. Axisymmetric Hydrodynamics 107 Finally, we have also tried to implement a finite-difference based algorithm to solve the equilib-rium equations for rotating configurations, Sec. 5.3. Our first approach involved the use of multigrid methods of the type described in Sec. 5.4. However, because of our inability to set proper regular-ity conditions on the axis, this approach has not been successful. The fact that most of the codes developed to date to find rotating equilibrium configurations are based on some kind of pseudo-spectral method (see [97] and references within) has inspired us to adopt such a approach. W i t h a spectral technique, regularity conditions can be satisfied by choosing a set of basis (expansion) functions with the appropriate regularity behavior. We plan on developing such an algorithm in the near future. Chapter 6. Conclusions and Directions for Future Research 108 Chapter 6 Conclusions and Directions for Future Research In this thesis we have considered three different problems in numerical relativity. In these projects we have studied self-gravitating dynamics of three distinct types of matter and/or energy: scalar fields; vacuum gravitational energy; and perfect fluids. In two of the chapters, namely Chap. 2 and 4, we studied systems of equations that, due to our restriction to spherical symmetry, were effectively 2-dimensional, i.e. where the dynamic variables depended on 1 spatial dimension as well as time. In Chap. 3 and 5, on the other hand, the problems were effectively 3-dimensional, with variables dependent on 2 spatial dimensions and time. We also made use of 2 different formalisms that cast Einstein field equations in a form suitable for treatment as a Cauchy problem: the so called 3+1 approach and the 2+1+1 formalism. In terms of the numerical analysis, we employed a variety of discretization techniques and solution algorithms in our studies. These included standard centred finite-difference approaches for P D E s of evolution type (scalar and gravitational fields), Godunov/finite-volume/HRSC techniques for the hydrodynamical evolution equations, and multigrid algorithms for the solution of discretized elliptic equations arising from the Einstein constraints. In the remainder of this chapter we summarize our main conclusions. Firstly, in Sec. 6.1, we summarize the results of the collapse of scalar field with a particular potential, as described in Chap. 2. A summary of results from the evolution of a 5-D black string, which we discussed in Chap. 3 follows. Finally we conclude with some of our findings in relativistic hydro. As a continuation of this thesis we plan to continue our numerical studies of the strong field regime in the coupled Einstein equations/hydrodynamical equations. Specifically we plan to incor-porate adaptive mesh refinement into our axisymmetric code. Chapter 6. Conclusions and Directions for Future Research 109 6.1 Scalar Field Collapse with Angular Momentum In Chap. 2 we studied the critical collapse of a scalar field in spherical symmetry. In this collab-orative effort [23], the equations of motion for a massless scalar field included a term that was characterized by an angular momentum parameter, I, and that could be viewed as providing an angular momentum barrier. The additional term is analogous to the angular momentum barrier in the 1-dimensional reduced mechanics problem of a particle moving in a central potential. Following [17], we studied the threshold of black hole formation in our model through parame-terization of the initial data with a single parameter, p. We found evidence for an entire family of critical solutions—one for each value of I—and, as was the case for the original study [17] (equiva-lent to the case 1 = 0), the threshold solutions were found to be discretely self-similar. The critical solutions were thus characterized by two different exponents: (1) A ; , which is the so-called echoing exponent, and which is actually a period when the solution is expressed in coordinates adapted to the scale-symmetry (self-similarity); and (2) a scaling exponent, 7/, that describes how dimension-ful quantities, such as the black hole mass, scale with parameter-space distance from criticality, (p — p*), as p —> p*. We found that both exponents decreased with I in an almost exponential manner. Moreover, we found that as I was increased, the time that a solution remained close to criticality (for fixed (p — p*)/p*) grew, and that the critical solutions appeared to approach periodic configu-rations. We believe that this fact can be understood as the angular momentum barrier having, at least partially, a stabilizing effect against gravitational collapse. 6.2 Instability of a Black String Chap. 3 described a study of the dynamical instability of a particular 5-dimensional relative of the usual 4-dimensional Schwarzschild solution, namely the black string originally, due to Myers and Perry [73]. This was another collaborative effort [18], and in order to study the instability we developed an effective 2+1 dimensional numerical relativity code to evolve the vacuum Einstein equations. The exposition in Chap. 3 focused on our numerical procedure to solve the constraint equations at the initial time, and on our implementation of an algorithm to approximately locate event horizons via "backwards evolution" of null rays, using the full results from the numerical simulations of the spacetimes. The results of our calculations were shown to agree with the predictions from perturbation theory due to Gregory and Laflamme [37]. In addition, we used our code to evolve perturbed Chapter 6. Conclusions and Directions for Future Research 110 black-string spacetimes well beyond the linear regime where perturbation theory is valid. In the non-linear regime, we found evidence that the perturbed spacetime evolves to a configuration that resembles a series of black holes connected by a black string with radius smaller than the progenitor string. Due to the development of a coordinate pathology, our code failed at late times, and we were thus unable to make any statement about the ultimate fate of unstable black strings, particularly since at the time the code crashes, the spacetimes were still highly dynamical. We also observed that, at least for the spacetimes considered, our approximations for the location of the event horizon at any instant of time agreed well with the locations of apparent horizons. 6.3 Relativistic Hydrodynamics In the final two chapters of this thesis, Chaps. 4 and 5, we described the development of codes to simulate, respectively, spherically symmetric, and axially symmetric, general relativistic hydrody-namics. The spherically symmetric code was based on a 3+1 decomposition of the Einstein equations and used coordinates which are the natural restriction to spherical symmetry of the coordinates used in the axisymmetric simulations. This allowed us to experiment with formalisms and numer-ical methods in a simpler setting than that provided by the axisymmetric case. In order to evolve the hydrodynamic quantities, a numerical technique tailored to treat discontinuities stably and cleanly—namely a so-called high resolution shock capturing method—was implemented. Develop-ment of the code also led to the identification of a new set of conservative variables, whose use was crucial in making the solution of the Hamiltonian constraint in the model a well-posed problem. In addition, we observed that—at least at the resolutions used in our simulations—our algorithm was too dissipative to obtain long term evolution of stationary (TOV) solutions. Chap. 5 considered the case of axisymmetric hydrodynamics. The equations of motion for the fluid and the geometry were written using the 2+1+1 formalism [34]. In addition, the fluid equations were expressed in conservation law form, which again permitted the use of high resolution shock capturing methods. We showed that our code converged (at least for the weak evolution shown in Fig. 5.6) as a function of mesh resolution as expected, and that it could stably evolve various configurations of discontinuous initial data. On the other hand, we found that the code was too unstable to permit study of the very strong-field regime, typified, for example, by configurations close to black hole formation. Additionally, as in the spherical case, we found that the algorithm used for the hydrodynamics was too dissipative to permit stable evolution of stars for long times. In order to solve some of these problems we have started implementing a code which adaptively refines the numerical grid of points on the locations where it is required. Some early experiments Chapter 6. Conclusions and Directions for Future Research 111 in spherical symmetry indicate that this is a promising approach. Bibliography 112 Bibliography [I] Abrahams A . M . and C. R. Evans, "Critical Behavior and Scaling in Vacuum Axisymmetric Gravitational Collapse", Phys. Rev. Lett. 70, 2980 (1993). [2] Anderson E . et al, " L A P A C K Users' Guide", Society for Industrial and Applied Mathematics, Philadelphia, (1992) [3] Arnowitt R., S. Deser and C. W . Misner, "The dynamics of general relativity", in Gravitation: An Introduction to Current Research, Ed . L . Witten, J . Wiley & Sons, Inc., New York (1962). [4] Baiotti ai, "Three-dimensional relativistic simulations of rotating neutron star collapse to a Kerr black hole", arXiv:gr-qc/0403029. [5] Baumgarte T. W. , S. L . Shapiro, M . Shibata, "On the maximum mass of differentially rotating neutron stars", Astrophys. J., 528 L29, (2000). [6] Berger M . J . , P. Collela, "Local adaptive mesh refinement for shock hydrodynamics", J. Comp. Phys. 82, 64-84, (1989). [7] Bogojevic A . R. and L . Perivolaropoulos, "Black Holes In A Periodic Universe", Mod. Phys. Lett. A 6, 369, (1991). [8] Bonazzola S., E . Gourgulhon, M . Salgado, and J . A . Marck, "Axisymmetric rotating relativistic bodies: a new numerical approach for "exact" solution", Astron. Astrophys. 278, 421-443, (1993). [9] Brandt A . , "Multi-level Adaptive Solutions to Boundary-value Problems", Math. Comput. 31, 330, (1977). [10] Bruhat Y . , "The Cauchy Problem", in Gravitation: An introduction to current Research, E d . L . Witten, J . Wiley & Sons, Inc., New York (1962). [II] C A C T U S Computational Toolkit, The. [12] Caveny S. A . , M . Anderson and R. A . Matzner, "Tracking Black Holes in Numerical Relativ-ity", Phys. Rev. D 68, 104009, (2003) [arXiv:gr-qc/0303099]. Bibliography 113 [13] Chandrasekhar S., An introduction to the study of stellar structure, University of Chicago Press, Chicago, (1939). [14] Chandrasekhar S., The mathematical theory of black holes, Oxford University Press, Oxford (1983). [15] Chi-Wang Shu, "Essentially non-oscillatory and weighted essentially non-oscillatory schemes . for hyperbolic conservation laws", Technical Report N A S A CR-97-206253 I C A S E Report No. 97-65, Institute for Computer Applications in Science and Engineering, November (1997). [16] Choptuik M . W. , "Consistency of finite-difference solutions of Einstein's equations", Phys. Rev. D 44, 3124-3135, (1991). [17] Choptuik M . W. , "Universality And Scaling In Gravitational Collapse of a Massless Scalar Field", Phys. Rev. Lett. 70, 9, (1993). [18] Choptuik M . W. , L . Lehner, I. Olabarrieta, R. Petryk, F . Pretorius and H . Villegas, "To-wards the final fate of an unstable black string", Phys. Rev. D 68, 044001, (2003) [arXiv:gr-qc/0304085]. [19] Choptuik M . W. , E . W . Hirschmann, S. L . Liebling, and F . Pretorius, " A n Axisymmetric Gravitational Collapse Code", Class. Quant. Grav. 20, 1857, (2003). [20] Choptuik M . W. , E . W . Hirschmann, S. L . Liebling and F . Pretorius, "Critical collapse of the massless scalar field in axisymmetry", Phys. Rev. D 68, 044007, (2003) [arXiv:gr-qc/0305003]. [21] Choptuik M . W. , E . W . Hirschmann, S. L . Liebling and F . Pretorius, "Critical collapse of a complex scalar field with angular momentum", arXiv:gr-qc/0405101. [22] Donat R., and A . Marquina, "Capturing Shock Reflections: A n Improved Flux Formula", J . Comput. Phys. 125, 42-58, (1996). [23] Choptuik M . W. , I. Olabarrieta, W. Unruh, J . Ventrella (in preparation). [24] Christodoulou D., "The problem of a self-gravitating scalar field", Commun. Math. Phys. 105, 337-361, (1986). [25] Christodoulou D., " A mathematical theory of gravitational collapse", Commun. Math. Phys. 109, 613-647, (1987). [26] Eckart A . and R. Genzel, "Stellar Proper Motions In The Central 0.1 Pc Of The Galaxy", Mon. Not. Roy. Astron. Soc. 284, 576, (1997). Bibliography 114 [28] Evans C. R. " A method for numerical relativity: Simulation of axisymmetric gravitational collapse and gravitational radiation generation", Ph.D. Thesis, The University of Texas at Austin, (unpublished), (1984). [28] Evans C. R., " A n Approach For Calculating Axisymmetric Gravitational Collapse", Gravi-tational Radiation, in J . Centrella ed. Dynamical Spacetimes and Numerical Relativity, 3-39, Cambridge, England, Cambridge University Press, (1986). [29] Evans C. R., and J . S. Coleman, "Observation of critical phenomena and selfsimilarity in the gravitational collapse of radiation fluid," Phys. Rev. Lett. 72, 1782 (1994) [arXiv:gr-qc/9402041]. [30] Font J . A . , "Numerical hydrodynamics in general relativity", Living Rev. Rel. 3, 2, (2000) [Living Rev. Rel. 6, 4, (2003)] [arXiv:gr-qc/0003101]. [31] Font J . A . , M . Miller, W . - M . Suen and M . Tobias, "Three Dimensional Numerical General Relativistic Hydrodynamics: Formulations, Methods and Code Tests", Phys. Rev. D 61, 044011, (2000). [32] Garfinkle D. and G. C. Duncan, "Scaling of curvature in sub-critical gravitational collapse", Phys. Rev. D 58, 064024, (1998) [arXiv:gr-qc/9802061]. [33] Garfinkle D. , C. Gundlach and J . M . Martin-Garcia, "Angular momentum near the black hole • threshold in scalar field collapse", Phys. Rev. D 59, 104012, (1999) [arXiv:gr-qc/9811004]. [34] Geroch R., " A Method for Generating Solutions of Einstein's Equations", J. Math. Phys. 12, 918-924, (1971). [35] Ghez A . M . , B . L . Klein, M . Morris and E . E . Becklin, "High Proper Motion Stars in the Vicinity of Sgr A * : Evidence for a Supermassive Black Hole at the Center of Our Galaxy", Astrophys. J. 509, 678, (1998) [arXiv:astro-ph/9807210]. [36] Godunov S. K . , Mat. 56.-47, 271, (1959). [37] Gregory R. and R. Laflamme, "Black Strings And P-Branes Are Unstable", Phys. Rev. Lett. 70, 2837, (1993) [arXiv:hep-th/9301052]. [38] Gregory R. and R. Laflamme, "The Instability of charged black strings and p-branes", Nucl. Phys. B 428, 399, (1994) [arXiv:hep-th/9404071]. [39] GubserS. S., "On non-uniform black branes", Class. Quant. Grav. 19,4825, (2002) [arXiv:hep-th/0110193]. Bibliography 115 [40] Gundlach C , "Understanding Critical Collapse of a scalar field", Phys. Rev. D 55, 695-713, (1997) [arXiv:gr-qc/9604019]. [41] Gundlach C , "Angular momentum at the black hole threshold", Phys. Rev. D 57, 7080, (1998) [arXiv:gr-qc/9711079]. [42] Gundlach C , "Critical gravitational collapse of a perfect fluid with p = kp: Nonspherical perturbations", Phys. Rev. D 65, 084021, (2002) [arXiv:gr-qc/9906124]. [43] Gundlach C , "Critical phenomena in gravitational collapse", Phys. Rept. 376, 339, (2003) [arXiv:gf-qc/0210101]. [44] Harrison B . K . , K . S. Thorne, M . Wakano, and J . A . Wheeler, "Gravitation Theory and Gravitational Collapse", University of Chicago Press, Chicago Illinois (1965). [45] Hawke I., "Computational Ultrarelativistic Hydrodynamics", Ph .D. Thesis, Queen's College, Cambridge, (unpublished), (2001). [46] Hawke I. and J . M . Stewart, "The Dynamics Of Primordial Black Hole Formation", Class. Quant. Grav. 19, 3687, (2002). [47] Hawking S. W. , and G . F. R. Ellis, "The large scale structure of space-time", Cambridge University Press, Cambridge, (1973). [48] Horowitz G. T. and K . Maeda, "Fate of the black string instability", Phys. Rev. Lett. 87, 131301, (2001) [arXiv:hep-th/0105111]. [49] Hughes S. A , C. R. Keeton, P. Walker, K . T. Walsh, S. L . Shapiro and S. A . Teukolsky, "Finding black holes in numerical space-times", Phys. Rev. D 49, 4004, (1994). [50] Husa S., C. Lechner, M . Purrer, J . Thornburg and P. C. Aichelburg, "Type II critical col-lapse of a self-gravitating nonlinear sigma-model", Phys. Rev. D 62, 104007, (2000) [arXiv:gr-qc/0002067]. ' [51] Kel ly B . , P. Laguna, K . Lockitch, J . Pull in, E . Schnetter, D . Shoemaker and M . Tiglio, " A cure for unstable numerical evolutions of single black holes: adjusting the standard A D M equations", Phys. Rev. D 64, 084013, (2001) [arXiv:gr-qc/0103099]. [52] Kreiss H . and J . Oliger, "Methods for the Approximate Solution of Time Dependent Prob-lems", Global Atmospheric Research Programme, Publications Series No. 10. (1973). Bibliography 116 [53] Koike T., T. Hara and S. Adachi, "Critical behavior in gravitational collapse of radiation fluid: A Renormalization group (linear perturbation) analysis", Phys. Rev. Lett. 74, 5170 (1995) [arXiv:gr-qc/9503007]. [54] Landau L . D . and E . M . Lifshifz, Fluid Mechanics, Pergamon Press, London, (1959). [55] Lehner L . , M . Huq and D. Garrison, "Causal Differencing In A D M and Conformal A D M Formulations: A Comparison In Spherical Symmetry", Phys. Rev. D 62, 084016, (2000). [56] Lehner L . , "Numerical relativity: A review", Class. Quant. Grav. 18, R25, (2001) [arXiv:gr-qc/0106072]. [57] LeVeque R. J . , Numerical Methods for Conservation Laws, Basel, Birkhauser-Verlag (1990). [58] LeVeque R. J . , Finite Volume Methods for Hyperbolic Problems, New York, Cambridge Uni-versity Press (2002). [59] LeVeque, R. J . , " C L A W P A C K Version 4.2 User's Guide", (2003). [60] Libson J . , J . Masso, E . Seidel, W . M . Suen and P. Walker, "Event horizons in numerical relativity. 1: Methods and tests", Phys. Rev. D 53, 4335, (1996) [arXiv:gr-qc/9412068]. [61] Liebling S. L . , "Critical phenomena inside global monopoles", Phys. Rev. D 60, 061502, (1999) [arXiv:gr-qc/9904077]. [62] Liebling S. L . , "The singularity threshold of the nonlinear sigma model using 3D adaptive mesh refinement", Phys. Rev. D 66, 041703, (2002) [arXiv:gr-qc/0202093]. [63] Maeda K . , M . Sasaki, T. Nakamura, S. Miyama, " A new formalism of the Einstein Equations for Relativistic Rotating Systems", Prog. Theor. Phys. 63, 719-721, (1979). [64] Maple 8, Waterloo Maple Inc. [65] Marsa R. L . and M . W . Choptuik, "The R N P L User's Guide",^uide/users_guide.html (1995). [66] Mar t i ' J . M . , J . M . Ibanez, and E . Muller, "Numerical Relativistic Hydrodynamics: Local Characteristic Approach", Phys. Rev. D 43, 3794-3801, (1991). [67] Mar t i J . M . , and E . Mueller, "Numerical hydrodynamics in special relativity", Living Rev. Rel. 2, 3, (1999) [arXiv:astro-ph/9906333]. [68] May M . M , and R. H . White, "Hydrodynamic calculations of general relativistic collapse", Phys. Rev. D 141, 1232-1241, (1966). Bibliography 117 [69] May M . M , and R, H . White, "Stellar dynamics and gravitational collapse", Methods Comput. Phys. 7, 219-258, (1967). [70] Miralles J . A . , J . M . Ibafiez, J . M . Mart i , and A . Perez, "Incompressibility of hot nuclear matter, general relativistic stellar collapse and shock propagation", Astron. Astrophys. Suppl. 90, 283-299, (1991). [71] Misner C. W. , K . S. Thorne and J . A . Wheeler, Gravitation, New York, W . H . Freeman and Company (1973). [72] Miiller I. "Speeds of Propagation in Classical and Relativistic Extended Thermodynamics", Living Rev. Rel. 2, 1, (1999). [73] Myers R. C. and M . J . Perry, "Black Holes In Higher Dimensional Space-Times", Annals Phys. 172, 304, (1986). [74] Nakamura T, "General Relativistic Collapse of Axial ly Symmetric Stars Leading to the For-mation of Rotating Black Holes", Prog. Theo. Phys. 65, (1981). [75] Nakamura T., and H . Sato, "General Relativistic Collapse of Non-Rotating, Axisymmetric Stars", Prog. Theor. Phys. 67, 1396-1405, (1982). [76] Nakamura T., K . Oohara and Y . Kojima, "General Relativistic Collapse of Axial ly Symmetric Stars and 3D Evolution of Pure Gravitational Waves", Prog. Theor. Phys. Suppl. 90, 13-101, (1987). [77] Neilsen D . W. , "Extremely Relativistic Fluids in Strong-Field Gravity", Ph .D. Thesis, The University of Texas at Austin, (unpublished), (1999). [78] Neilsen D . W. and M . W . Choptuik, "Ultrarelativistic Fluid Dynamics", Class. Quantum Grav. 17, 733-759, (2000) [arXiv:gr-qc/9904052]. [79] Noble S. C. , " A Numerical Study of Relativistic Fluid Collapse", arXiv:gr-qc/0310116. Ph.D. Thesis, The University of Texas at Austin, (unpublished), (2003). [80] Oppenheimer J . R., and G . M . Volkoff, "On Massive Neutron Cores", Phys. Rev. 55, 374-381, (1939). [81] Olabarrieta I., and M . W . Choptuik, "Critical phenomena at the threshold of black hole formation for collisionless matter in spherical symmetry", Phys. Rev. D 65, 024007, (2002) [arXiv:gr-qc/0107076]. Bibliography 118 [82] Orosz J . A . , "Inventory of Black Hole Binaries", [arXiv:astro-ph/0209041]. [83] Overduin J . M . and P. S. Wesson, "Kaluza-Klein gravity", Phys. Rept. 283, 303, (1997) [arXiv:gr-qc/9805018]. [84] Petzold L . R. and A . C. Hindmarsh, " L S O D A " , Computing and Mathematics Research Divi -sion, 1-316 Lawrence Livermore National Laboratory, Livermore, C A 94550. [85] Piran T. and R. F . Stark, "Numerical Relativity, Rotating Gravitational Collapse and Gravi-tational Radiation", in J . Centrella ed. Dynamical Spacetimes and Numerical Relativity, 40-73, Cambridge, England, Cambridge University Press, (1986). [86] Press W . H . , S. A . Teukolsky, W . T. Vetterling, and B . P. Flannery, "Numerical Recipes in Fortran 77", (2nd edition), Cambridge University Press, New York (1992). [87] Pretorius F . , "Numerical Simulations of Gravitational Collapse", Ph .D. Thesis, The University of British Columbia, (unpublished), (2002). [88] Rein G. , A . D. Rendall and J . Schaeffer, "Critical collapse of collisionless matter: A numerical investigation", Phys. Rev. D 58, 044007, (1998) [arXiv:gr-qc/9804040]. [89] Rhoades C. E . and R. Ruffini, "Maximum Mass of a Neutron Star", Phys. Rev. Lett. 32, 324-327, (1974). [90] Roe P. L . , "Approximate Riemann solvers, parameter vectors and difference schemes", J. Comput. Phys. 42, 357-372, (1981). [91] Romero, J . V . J . M . Ibanez, J . M . Mart i , J . A . Miralles, " A New Spherically Symmetric General Relativistic Hydrodynamical Code", Astrophys. J. 462, 839-854, (1996) [arXiv:astro-ph/9509121]. [92] Sasaki M . , K . Maeda, S. Miyama, T. Nakamura, " A Method of Determining Apparent Horizons in [(2+l)+l]-Formalism of the Einstein Equations", Prog. Theor. Phys. 63, 1051-1053, (1980). [93] Seidel E . , and W . M . Suen, "Towards a singularity-proof scheme in numerical relativity", Phys. Rev. Lett. 69, 1845-1848, (1992). [94] Shapiro S. L . and S. A . Teukolsky, "Black Holes, White Dwarfs, and Neutron Stars", J . Wiley & Sons, Inc., New York (1983). [95] Shibata M . , "Axisymmetric Simulations of Rotating Stellar Collapse in Full General Relativ-ity", Prog. Theor. Phys. 104, 325-358, (2000). Bibliography 119 [96] Shinkai H . A . and G. Yoneda, "Adjusted A D M systems and their expected stability proper-ties", Class. Quant. Grav. 19, 1027, (2002) [arXiv:gr-qc/0110008]. [97] Stergioulas N . , "Rotating stars, in relativity", Living Rev. Rel. 6, 3, (2003) [arXiv:gr-qc/0302034]. [98] Swesty D. , J . M . Lattimer, and E . S. Myra , "The role of the equation of state in the 'prompt' phase of type II supernovae", Astrophys. J. 425, 195-204, (1994). [99] Synge, J . L.,The relativistic gas, Noth-Holland, Amsterdam, (1957). [100] Teukolsky S. A . , "On the Stability of the Iterated Crank-Nicholson Method in Numerical Relativity", Phys. Rev. D 61, 087501, (2000) [arXiv:gr-qc/9909026]. [101] Thornburg J . , "Coordinates and boundary conditions for the general relativistic initial data problem", Class. Quant. Grav. 4, 1119-1131, (1987). [102] Tolman R. C , Relativity, Thermodynamics and Cosmology, Oxford University Press, (1939). [103] Tolman R. C. , . "Static Solutions of Einstein's Field Equations for Spheres of Flu id" , Phys. Rev. 55, 364-373, (1939). [104] Unruh W. , and R. Wald, (in preparation). [105] Van der Marel R. P., "Intermediate-Mass Black Holes in the Universe? - A Review of For-mation Theories and Observational Constraints", arXiv:astro-ph/0302101. [106] Van Leer B . J . , "Towards the conservative difference scheme. V . A second order sequel to Godunov's method", J. Comput. Phys. 32, 101-136, (1979). [107] Ventrella J . F . and M . W . Choptuik, "Critical phenomena in the Einstein-massless-Dirac system", Phys. Rev. D 68, 044020 (2003) [arXiv:gr-qc/0304007]. [108] VonNeumann J . , and R. D. Richtmyer, " A method for the numerical calculation of hydrody-namic shocks", J. Appl. Phys. 21, 232-247, (1950). [109] [110] Weinberg S., Gravitation and cosmology, J . Wiley & Sons, Inc., New York (1972). [ I l l ] Wilson J . R. "Numerical Study of Fluid Flow in Kerr Space", Astrophys. J. 173, 209-219, (1971). Bibliography 120 [112] Wilson J . R. " A numerical Method for Relativistic Hydrodynamics", in L . Smarr ed. Sources of Gravitational Radiation. 423-445, Cambridge, England, Cambridge University Press, (1979). [113] Wilson J . R., G . J . Mathews, and P. Marronetti, "Relativistic Numerical Model for Close Neutron Star Binaries", Phys. Rev. D 54, 1317-1331, (1996). [114] Wiseman T., "From black strings to black holes", Class. Quant. Grav. 20, 1177, (2003) [arXiv:hep-th/0211028]. [115] York J . W . Jr., "Kinematics and Dynamics of General Relativity", in Sources of Gravitational Radiation, E d . L . Smarr, Seattle, Cambridge University Press (1979). [116] York J . W. , T. Piran, "The Initial Value Problem and Beyond", in Spacetime and Geometry, ed. by R. Matzner and L . C. Shepley, U T Press, Austin, (1982). Appendix A. Hydrodynamic equations in conservation form 121 Appendix A Hydrodynamic equations in conservation form The general relativistic hydrodynamic equations can be cast in conservation law form. This fact was first fully exploited in the work of Martiand Muller [66] (also see [30] for a review), who developed what is now known as the Valencia formulation. Here we explicitly show the details of the calculation that puts the fluid equations into such a form. Specifically, our goal is to show that (1.36) and (1.37) can be written as dQj dFh dt dxi Si (A . l ) where the Si do not contain any derivatives of the fluid variables. We begin with (1.36) (T\),^ = 0, and then perform the following manipulations: (A.2) (<K„T"«). 1 = gsu '-g l -g -gT^^ + T5^ u + r { ^ = 0. (A.3) In the above, .we have used equation (4.7.9) of [110] to go from the first to the second line and the Leibniz rule elsewhere. Also note that g is the determinant of the metric gpu. Expression (A.3) can be written as: -j= (^TgT^)^ = g^T^ - T\xT»xgSv. (A.4) This set of equations is thus in the form (A. l ) since the stress tensor for a perfect fluid does not contain any derivatives of the fluid variables. The equation of particle number conservation (1.37) Appendix A. Hydrodynamic equations in conservation form 122 is even simpler to massage into conservation form. Indeed = 0, (A.5) can be written as (A.6) using equation (4.7.7) of [110], and this is also of the desired form. Appendix B. Scalar fields considered in the 1=1 case 123 Appendix B Scalar fields considered in the 1=1 case For illustrative purposes, in this Appendix we explicitly display the scalar fields *}"(£, r, 6, <p) that are considered in the calculations described in Chap. 2, for the specific case I = 1. Note that we want to consider real eigenfunctions of the L2 operator. Moreover, all of the fields J^™ are required to have the same functional dependence on t and r, which we denote ip(t,r). We then have the following 3 fields: = tft,r)Y?{0,<p)=tft,r)SJ^-cos(e), (B.1) * 2 - ^tft,r)[Yl(0^)+Y{-\9,cP)]=tft,r)^sm(e)sm(cp), (B.2) * 3 = ^tft,r)[YlHO,^)-Y{-\e,cP)]=tft,r)]l^-Sm(e)cos(<p), (B.3) where Yfi^O,^) denotes the usual spherical harmonic of degree / and order s. The equations of motion for ip(t, r) can then be derived by taking the divergence df the stress-energy tensor T^~^ which in this case is: T ( i = i ) _ T ( i i ) r ( i 2 ) T ( i 3 ) m 4 x 1ab —±ab + i o f c i " J a b • \ D - ^ ) Appendix B. Scalar fields considered in the 1=1 case 124 where each of the T^" 1 ^ is computed using equation (2.3). Specifically, we have , ( i i ) ( i i ) ( i i ) te T ( i i ) 1t<t> T (n) T, 8?r (n2 + $ 2 ) cos(fl)2 V>2 sin(fl)2 3 aII$cos(<9)2 47T a ' 3 a l t y cos(0) sin(0) 47T = 0, 3 8TT (n2 + $2) cos(0) 2 V 2 a 2 sin(0) 2 , ( i i ) _ - —-$i/> cos(0) sin(0), 47T = 0, ( i i ) _ _3_ 8TT (n2 - $2) r 2 cos(0)2 - V 2 sin(0) 2 ( i i ) _ = 0, - ( H ) _ 8TT sin(0)2 (n2 - $2) r 2 cos(0)2 i/>2 sin(0) 2 (B.5) (B.6) (B.7) (B.8) (B.9) (B.10) ( B . l l ) (B.12) (B.13) (B.14) Appendix B. Scalar fields considered in the 1=1 case 125 (12) T(12) 3 c? ( - r 2 n 2 + r 2 n 2 cos(0)2 - r 2 $ 2 + r2& cos(<ft)2 + V>V - V>2 cos(c6)2a2) cos((?)2 87r . a 2 r 2 3 a 2 (r 2n 2 - r 2 n 2 cos(<£)2 + r 2 $ 2 - r 2 $ 2 cos(</>)2 + </>2 c o s ^ a 2 ) —— • : , (B.15) + 8TT r(i2) = 3 an ( -1 + cos(4>)2) $ cos(g)2 3 an ( -1 + cos(^)2) $ t r 47T a ' 47T a ' 3 an(-l + cos(<?!>)2)sin(6l)i/.cos(6l) , (B.17) 47r a T(i2) = 3 an sin(</>)^  cos(0) cos(6i)2 3 aIIsiii(i/>)V>cos(<ft) . . ^ 4TT a 4TT a 1 ' ' y( i2) = 3 ( - r2 $ 2 + r 2 $ 2 cos^) 2 - r 2 n 2 + r 2 n 2 cos(0)2 - V>2a2 + ip2 cos(</>)2a2) cos(6Q2 8^ r2" 3 r2$)2 _ r2$)2 cos^ )2 + r2U2 _ ^2 cos(^ )2 _ ^2 cos^a2 (B.19) 8TT T r ( e 1 2 ) = ^-$Vsin(</»)2sin(0)cos(0), (B.20) 3 Air Trl2) = *^sin(0)sin(0)Vaw(>), (B.21) T^ 1 2 ) = - A (-V» 2cos(0) 2 + ^ 2cos(0) 2cos(^) 2 + V 2cos(^) 2) 3 2 n 2 cosU)2 sm(e)2 - n2 sin(0)2 - $ 2 cos(^)2 sin(0)2 + sin(0) 2 $ 2 „ , - — r 5 ' (B.22J 87r a J TS2) = ^^2sin(<?!.)cos(6i)cos(0)sin(6i), (B.23) ( 1 2 ) _ 3 ( - r 2 n 2 + r 2 n 2 cos(</>)2 + r 2 $ 2 - r 2 $ 2 cos(0)2 - ip2a2 + j>2 cos(4>)2a2) cos(0)4 T ^ ~ ~ 8 ^ a2 3 (2 r 2 n 2 - 2 r 2 n 2 cos(<£)2 - 2 r 2 $ 2 + 2 r2$2 cos(^)2 + ^>2a2) cos(fl)2 _ _ _ 3 -^ 2 cos (^ ) 2 (a ) 2 - r 2 n 2 +r 2 n 2 cos (< / . ) 2 +r 2 $ 2 - r 2 $ 2 cos (^ ) 2 _ _ . (B.24J 87r az Appendix B. Scalar fields considered in the 1=1 case 126 .(13) _ R ( 1 3 ) + 3 a 2 ( - r 2 n 2 cos(</>)2 - tfo2 - r 2 $ 2 cos(<p)2 + tf cos(cp)2a2) 8TT a2r2 3 a2tfcos(0)2 87 r 2 sin(6>)2 (13) _ te T ( 1 3 ) _ J t 0 — 7^(13) = rr T ( 1 3 ) _ 1re — T ( 1 3 ) _ ±r<t> ~ T ( 1 3 ) _ T ( 1 3 ) _ 3 an$cos(</>)2sin(6>)2 47T a ' 3 an cos(</>)2 sin(g)V> cos(fl) 47r a ' 3 aIIsin((/>)t/>cos(<i) . v ' -^-sm(6>) , 4w a 3 - r 2 $ 2 cos(<?!>)2 + tfa2 - r2U2 cos(cp)2 - tf cos(cp)2a2 "8TT r 2 3 tfa2cos{d)2 " t o r 2 — $ cos(</>)2 s'm(9)ip cos(0), 4TT 3 3 —$sin(<£)V> cos(0) cos(0) - —$sin(0)^cos(0), 47T 47T 3 r 2 n 2 cos(0)2 + tf cos(4>)2a2 - r 2 $ 2 c o s i » 2 - V> 2- 2 8?r a 2 sin(<9)2 • sin(6>)2 + ^ -tfcos(e)2 (2cos (0 ) 2 - l ) , 47T V>2 sin(</>) cos(#) cos(^) sin(0). (B.25) (B.26) (B.27) (B.28) (B.29) (B.30) (B.31) (B.32) (B.33) Appendix C. Hydrodynamic equations in the 2+1+1 formalism 127 Appendix C Hydrodynamic equations in the 2+1+1 formalism In this Appendix we rederive the hydrodynamic equations in the 2+1+1 formalism as originally performed by Maeda et al [63]. In particular the following derivation includes some steps not detailed in [63]. We begin by writing the stress-energy tensor in the following way: T ^ = I r ^ r + + + TijYnvi, ( c i ) where we have used the following projections: r = T^C, (C2) Ti = 7 ^ T M „ , (C.3) nj = lU^- (C :4) (We note that we adopt a different notation here for the various projections of the stress energy tensor than that used in Chap. 5. We make a connection between the two notations at the end of the appendix.) By construction, the projection operator satisfies ev„ = o. (c.5) The fluid equations are derived from equations (1.36) and (1.37). We begin with (1.36): T»vill = 0, (C.6) which can be decomposed in the following way: C(T^)ltl = 0, (C.7) 7i„(rn i A t = 0. (C.8) We then focus on equation (C.7), r (TM,);m = (rr^) ; / i - r"„ (n ; /, = o, (C9) Appendix C. Hydrodynamic equations in the 2+1+1 formalism 128 where the second term on the right hand side is zero since £" satisfies the Ki l l ing equations (£n;v = 0), and T 7 ^ is symmetric. Substituting ( C l ) in (C.9) and performing some manip-ulations, we obtain the following expression: ( ^ " ) + (T^);, = + i ^ + 0"% + r* ( V i ) ; , = 0. (CIO) Here a vertical bar ()|4 denotes the covariant derivative Di in the 3 dimensional quotient spacetime, which is defined by projection of the spacetime covariant derivative, V M , via Di = 7 A i j V M . The first term of (C.10), which is £^(T/S2), vanishes since £ is Ki l l ing . The second term is zero because it is proportional to the trace of the Ki l l ing equations. Thus we have 2. Using the Ki l l ing equation once more, we can express the second term as the derivative of s T* , o x 1 I* 1 2s 2 v >M s v 'I We now turn attention to equation (C.8): 7 i , C ^ " ) ; , , = ( 7 - ^ " ) . ^ - ( 7 „ ) ; M = 0. (C.13) Using ( C l ) and the definition of the projection operator (5.2), we get the following expression: ( ^ 7 ^ 7 ^ ) + ( 7 - 7 ^ 7 V % + T ^ ( | ) = 0 . (C.14) Further manipulation gives ^ ( n ) ; , e + (jh),, + rh ( 7 ^ ) ; / t + ( | ) = 0, (C.15) and combining the second and third terms of the right hand side we find ^ ( ^ r + ; ( ^ ) b . + T ^ (jj) =0. (C16) We now consider the first term of this last equation, which, using the definition of T j , can be written as \ ^ \ ^ = ^^T^)^- ( c - 1 7 ) >i" s Expanding the derivative of the product we have ^ [ ( l \ ) . ^ T x a ^ + l \ C ^ T X a ^ + 1 \ ( T X a ) . ^ } - (C18) We now use the fact that the Lie derivative along the Ki l l ing vector field of the stress energy tensor vanishes, i.e. £ ^ = e (v);A+TV (e).u+TVX (e) M=o, (c.19) Appendix C. Hydrodynamic equations in the 2+1+1 formalism 129 to further manipulate (C.76): 1 (C.20) The second and the third terms in the above sum to 0 (as can be seen by relabelling dummy indexes), and therefore (C.17) can be written as (C.21) We can now introduce this last result into equation (C.16): ( S T ^ ) b . + T ^ (7 / i i ) ; A^ A 0 2 7A* ( ^ ) ; A + ( „2 = 0. (C.22) The expression in brackets can be further simplified using the definition of the projection operator (5.2) Some algebra gives 1 4 ^ ( 6 ) ; A ^ A „ 2 ( ^ ) ; i + ( „ 2 = o, 0. (C.23) (C.24) and then using the. Ki l l ing equations repeatedly, we find J . . S = 0. (C.25) Using the definition of r (C.2) we have -s (--), + r^„ [- (| We now use expression A12 from [34]: V ^ = ^ W / ^ r ^ + 2^2^ V M ] a 2 , to simplify the second term of the equation (C.26): r ^ r = 0 . 2 1 + \ o2 6 W i t h this.result equation (C.26) becomes 1 (C.26) (C.27) (C.28) (C.29) Appendix C. Hydrodynamic equations in the 2+1+1 formalism 130 nlrii = - 1 , = « f cK-)|fe = 9,- (a) /a , = 0,0), - <;•-£>• This concludes the initial decomposition of equations (C.6) in directions parallel and orthogonal to and using quantities adapted,to the projection with respect to the Ki l l ing field. We now proceed to a space-plus-time (2+1) split of equations (C.12) and (C.29). To that end we introduce another projection operator, Hlj, Hij=1ij+71%, (C.30) where nl is the unit-norm, future-directed, orthogonal vector field to the constant time surfaces (which, in an abuse of terminology we will refer to as hypersurfaces). The following relations will be useful in the subsequent derivations: (C.31) (C.32) (C.33) (C.34) (C.35) Note that a is the lapse function, while (31 is the shift vector. In addition, the double vertical bar ()||/ denotes covariant differentiation in the 2-dirriensional spacelike surfaces of our dimensionally reduced spacetime (i.e. differentiation compatible with the 2-metric Hu). We proceed to a decomposition of (C.12) by first splitting T1 into hypersurface-orthogonal and hypersurface-tangential pieces: T* = J ^ n * + S'Wj, (C.36) Jt = -ntf, (C.37) S1 = H'JTK (C.38) Inserting this decomposition into equation (C.12) we obtain -s (sJ+rf),. + -s (sS'lfi),. = 0, (C.39) and expanding the derivatives of the products, we find -s (sJ^n* + J , (n<),. + -s (sS1)^ + S1 ( JP , ) , . = 0. (C.40) It is now convenient to define a new vector field, Nl = an1. Making use of this definition as well as (C.35) we can write (C.40) as _±=£t [yg) -1 ( /? ' )„ , ] + \ (ss% + s1 («v)ti = o. (c.4i) — {sJ^u^ + Jt rv.s l J Appendix C. Hydrodynamic equations in the 2+1+1 formalism 131 Now, using the fact that t = AT' + /?' and rij = (0,0) and regrouping terms, we get the first of our final equations: ?—£t (syfHJ+) + — sJH V V as as \ b — JA,— 0. 11/ We now similarly decompose (C.29): The hypersurface-orthogonal and hypersurface tangential components of the above are n j - (ST\). . - r n ^ + r j e i j k t o k = 0, h'i- (srh)u - TirM + ff'A^'eyfcW* = 0. We begin be manipulating the first term of equation (C.44), Using the decomposition rU = mtfpH + niHjjJJ + nPHuJ1 + H^H^Su, where we have made use of the following definitions we obtain p H = niTijT-1, JJ = - n i H J j T i j , Su = HHHJJT13 , - ^ ( - n V H - ^ j J " 7 ) ] , , . - ^ ^ ) , . . After some algebra and the use of (C.35) we obtain \=\£t (sVHpn) - — as (J1 - pn^-) - r M n l ) i asVH \ / as [ \ a J] ^  ,x '\ Using the following property (C.42) (C.43) (C.44) (C.45) (C.46) (C.47) (C.48) (C.49) (C.50) (C.51) (C.52) (C.53) and the fact that nlKij = 0 and n J ai 7 - = 0 we find that the first term of equation (C.52) can be written l—£t(s^/HPli)-— \as ( V - p H —)] + (P + P H ) \vIvJKIJ-v1^-] + PKrz. (C.54) Appendix C. Hydrodynamic equations in the 2+1+1 formalism 132 Here, we have made use of the fact that the stress energy tensor under consideration is that for a perfect fluid. Therefore we have: Su = (P + pn)v'vJ + HIJP, (C.55) J1 = (P + puW. (C.56) We now shift attention to the remaining terms of (C.44). The second term can be written as (C.57) where we have defined K via K = —suri 1. In order to simplify the third term of (C.44), namely \r3CijkLokn\ we use the following definitions: vJ = -niClH + HiJClJ, Cln = riiLO1, cii = #W. Introducing these into (C.58) we obtain ^ (J^ + S ' f f A ) £ i j k (-nkClH + hkJ) n\ Now, by the antisymmetry of Cijk, the only term that survives in the above is which has only a single factor of ri1. We thus have E1Si, = J^—^Eiv1, where we have made the following definitions: c-ij — TI, ekij, * = 1 2s' E1 = ^-eIJClj. Collecting all of the terms, equations (C.54), (C.57),(C.64), together we find -7=£t (sy^Hpii) - — as (J1 - pu — asVH v J as I \ a (P + PH) a J + PKII + T^ + J^EIVI .= 0 (C.58) (C.59) (C.60) (C.61) (C.62) (C.63) (C.64) (C.65) (C.66) (C.67) Appendix C. Hydrodynamic equations in the 2+1+1 formalism 133 where J1 = (pH + P ) v1, r = J ^ s 2 + Ps2, E$ = KS. Thus our final equation is: 1 £T (SVHPH) + — as{jI-pii as\/H v 7 a s [ 1 -III/ a J " ( J ' + P H ) _p (V, _ + J_ j, (LfyV4> + 2Elv^ = 0. We now turn to (C.45). Manipulation of its first term yields Using the decomposition = ninipH+niHjjJJ + n^HuJ1 + WJHiKSJK, yields J [sH\niHiKJK + HjJH*[SJKHki] b. - ,.. Now, using the definition of the projection operator Hij, property (C.35), and the fact that tl = Nl + Pl, we find, after some algebra, that -^j=£t (Vh) - — {PJ)uj + -—£t (sJi) - — {sJv)^jPJ -a\/H ^ ' a IIj as as 11 Ji + h (asSJi)w+p«nJ + J j (ni)\\j Regrouping terms we have (C.68) (C.69) (C.70) (C.71) (C.72) as ( 5 J 7 - — J/ a —(sy/HJi) + — a s V / T v J as J-±d3J)^ + p ^ { n i ) ^ + JJ ( n 7 ) p . Now using identity (4.7.6) from [110]: (yJ/)NJ = -!= (Vffy.-7/)^  - y J K (HKI)tJ + J ^ r ^ y expression (C.73) becomes ±=£t (SVHJI) + — \ = \asVH (SJ, - ^ j/) (C.73) (C.74) a s \ gJK _ P^jK Jj ([3J) +pHni(nI)lj + JJ (n 7), (C.75) Appendix C. Hydrodynamic equations in the 2+1+1 formalism 134 Expanding the Christoffel symbols, we have —^=£t (sVHJi) + — H gJK _P^jK a as\/hJ asVH~[SJI-—JI a UHiK),j-l(Hu)tK + UHjK)tI ^ { ( / 3 J ) j 7 + \HJLf}K [(HLK)tI + (HLI),K ~ (Hu)tL\ ] + pm3 {ni),- + JJ (n / ) | , j . >\i Regrouping terms, we have rewritten the first term of equation (C.45) as (C.76) as ^=£t (sy/HjA + —^ = las^H (SJ , - —Ji VH V J asVH [ V a - S J " - ( H J K ) T I + a 11 a We now proceed to the last term of (C.45): 1 Using the decomposition and z T b £ijkWk Hl i. Ti = njJs + HjKSK, wk = - f i H n f c + CljHJk, we get -zJ<t>CljQJ j f / K ^ H ^ . (C.77) (C.78) (C.79) (C.80) (C.81) Here we have used the fact that HJKH1 jHkJe^fc = 0 as well as the definition of ejj. Now, using Bj, = 1 / 2 Q H and Ei = l/(2s)ejjQ,J we obtain —2—^J^Ei ^B<pSKejK- (C.82) Again using the fact that we have restricted attention to the case where the stress energy is that of a perfect fluid, we can use S1 = J^v1 which gives -•^J* (2El + -sB<kvKeIK^j . Collecting all the terms together, i.e. collecting expressions (C.76) and (C.83), (C.83) —^=£t (sVHJi) + —\= \asjH (sJj - ^Ji) asy/H V / asy/H [ \ a J a | U + JJ ( n / ) M J - Jj^ - (^2Ej + ^B^ej^j - r 1 a a -SJ"-{HjK)j-c 3 ' (C.84) Appendix C. Hydrodynamic equations in the 2+1+1 formalism 135 and using the fact that r = J^v^ + Ps2 we have ^-==£t (sy/HjA + —^7= \asVH ( 5 J 7 - — J y/H V V a s \ / # L V « - S J * 5 - J ^ + P H ^ + JJ (n7)MJ - V* if 2E, + -B*oV) " V * T - P— , (C-85) Q S a ' ; H J s z " V " s T y T T s J s Now, since €JKCkj = —SJj, BK = cIJdss and JJ (ni)^j = 0, we find our final expression for equation (C.45): l-=£t (SVHJA + —X-= LsVli (sJ! - J 7 -(PR + P) VJVK1-{HJK) J + V J ^ + P H ^ r ~ h3* 2El + eiK {~sB*vK~ \bKv+ .s (C.86) Equations (C.42), (C.68) and (C.86) represent the local conservation of energy, angular momen-tum and linear momentum, respectively. To make connection with the notation used in chapter 5 we make the identifications PH = T + D, (C.87) J l = 5 7 ) (C.88) J<j> = S,/,, (C.89) Si = S^vi, (C.90) Su = SJVJ + HJJP, (C.91) where the quantities on the left are the variables defined and used in this appendix, while those on the right are used in Chap. 5. Finally, we must rewrite the equation of (local) baryon number conservation, (C.92) The above implies sa\/H ^ ' ,M and then using the variables defined in Chap. 5, we have: \ = (saVHpou1*) = 0, —1-F=£T (SVHD) + —!-= \saVHD (v1 - ^-\ say/H V J satfH [ \ a J = 0. (C.93) (C.94) This completes our derivation of the hydrodynamical equations within the 2+1+1 formalism. Appendix D. Characteristic Structure 136 Appendix D Characteristic Structure The characteristic structure for the Jacobian or velocity matrix V?- = dF^/dq1 is given by the following set of eigenvalues: Ao = avp — (3P triply degenerate, (D-1) A± = Y^PcJ {yP ^ ~ ^ ± C s y / ( 1 ~ v 2 ) \ R P P ( 1 " y2°^ ~ V P v P ( 1 " ^} ~ P"' ( D ^ ) and the corresponding right eigenvectors: ro-1 = [ w ^ ' ^ - h w j ' ( D - 3 ) r 0 ,2 = (Wvz,2hW2vpvz,h(rl>i + 2W2vzvz) ,2hW2vzvc/j,Wvz (2hW - 1)) , (D.4) ro,3 = (Wv4>,2hW2vpv4,,2hW2vzv^h(s2 + 2W2v^)iWv4>(2hW-l)), (D.5) r± = (l,hWCi,hWvz,hWv,i„hWAp±-l). (D.6) In the above expressions we have made use of the following definitions K = ( D - 7 ) ii = «/po, (D.8) Cp± = vp-Vp±, (D.9) V£ = (vp-Ap±)/(l/^-vpAp±), (D.10) ^ = ( l / > 4 - «<V) / ( l / > 4 - vpAp±) , ( D . l l ) A£ = X±/a + pp/a, (D.12) where c 2 = l / / i (x + P/PIK), X — 8P/dp0 and K = dP/de. The characteristic structure of the Jacobian matrix in the z direction, V f . = dF*j/dql can be easily calculated from the above results using symmetry arguments, and is explained in [30]. Co-Authorship Statement In this statement I (Ignacio Olabarrieta) describe my role in the dif-ferent collaborations in which I have taken part during my P h . D . studies, and specifically, the collaborative projects which are described in my the-sis, P h . D . thesis "Relativistic Hydrodynamics and other topics in Numerical Relat ivi ty" . In Chapter 2 of my thesis I describe a collaboration wi th M . W . Choptuik, W . Unruh and J . Ventrella. This project involved the study of the crit ical collapse of a scalar field in spherical symmetry wi th a potential that mimics the effects of angular momentum. M y specific contributions included: (1) calculation of the equations of motion for any value of the angular momentum parameter; (2) development of a unigrid numerical code implementing the equations of motion; (3) generation and analysis of the data. Other people in the collaboration had the original idea for the precise model studied, extended my unigrid code to be adaptive, and helped wi th the analysis of the results. In Chapter 3, I describe a collaboration with M . W . Choptuik, L . Lehner, R. Petryk, F . Pretorius and H . Villegas, to dynamically study the instability of a 5 dimensional black string. In this chapter I focus on my main contribu-tions which were: (1) development of a code to solve the constraint equations of general relativity (which was used in order to compute the ini t ia l data) and (2) development of a code to compute the approximate event horizon. Other people were involved in the development of numerical codes to solve the evolution equations (using parallel methods), to locate apparent hori-zons and to compute geometric scalars. Others also performed the runs and collected the data. The work described in Chapter 4 was not done in collaboration wi th anyone. F ina l ly in Chapter 5, I describe my efforts to develop a code to solve the equations for relativistic hydrodynamic in axisymmetry. The work I describe in this chapter is a modification of a code developed by M . W . Choptuik, E . Hirschmann, S. Liebl ing and F . Pretorius to study the collapse of a scalar field. M y work involved the following: (1) obtaining the equations of motion and equilibrium for a perfect fluid in the 2+1+1 formalism and in conser-vation law form, (2) modifying the previously existing code to solve the equations that fix the metric coefficients to include a perfect fluid as a source and (3) wri t ing a code to evolve the hydrodynamic equations using H R S C methods. Ignacio (Inaki) Olabarrieta Research Supervisor: Matthew W . Choptuik Vancouver, 4 August 2004 


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items