@prefix vivo: . @prefix edm: . @prefix ns0: . @prefix dcterms: . @prefix dc: . @prefix skos: . vivo:departmentOrSchool "Science, Faculty of"@en, "Physics and Astronomy, Department of"@en ; edm:dataProvider "DSpace"@en ; ns0:degreeCampus "UBCV"@en ; dcterms:creator "Tonita, Aaryn"@en ; dcterms:issued "2009-01-05T16:06:11Z"@en, "2009"@en ; vivo:relatedDegree "Master of Science - MSc"@en ; ns0:degreeGrantor "University of British Columbia"@en ; dcterms:description """The production of axisymmetric initial data for distorted black holes at a moment of time symmetry is considered within the (3+1) context of general relativity. The initial data is made to contain a distorted marginally trapped surface ensuring that, modulo cosmic censorship, the spacetime will contain a black hole. The resulting equations on the complicated domain are solved using the piecewise linear finite element method which adapts to the curved surface of the marginally trapped surface. The initial data is then analyzed to calculate the mass of the space time as well as an upper bound on the fraction of the total energy available for radiation. The families of initial data considered contain no more than few percent of the total energy available for radiation even in cases of extreme distortion. It is shown that the mass of certain initial data slices depend to first order on the area of the marginally trapped surface and the gaussian curvature of prominent features."""@en ; edm:aggregatedCHO "https://circle.library.ubc.ca/rest/handle/2429/3325?expand=metadata"@en ; dcterms:extent "727289 bytes"@en ; dc:format "application/pdf"@en ; skos:note "Initial Data for Axially Symmetric Black Holes With Distorted Apparent Horizons by Aaryn Tonita B.Sc, Acadia University 2005 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF Master of Science in THE FACULTY OF GRADUATE STUDIES (Physics) The University of British Columbia May 2009 c© Aaryn Tonita, 2009 Abstract The production of axisymmetric initial data for distorted black holes at a moment of time symmetry is considered within the (3+1) context of general relativity. The initial data is made to contain a distorted marginally trapped surface ensuring that, modulo cosmic censorship, the spacetime will contain a black hole. The resulting equations on the complicated domain are solved using the piecewise linear finite ele- ment method which adapts to the curved surface of the marginally trapped surface. The initial data is then analysed to calculate the mass of the space time as well as an upper bound on the fraction of the total energy available for radiation. The families of initial data considered contain no more than few percent of the total energy available for radiation even in cases of extreme distortion. It is shown that the mass of certain initial data slices depend to first order on the area of the marginally trapped surface and the gaussian curvature of prominent features. ii Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Motivation for studying the problem . . . . . . . . . . . . . . . . . . 2 2 The Physical Problem . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.1 Conventions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2 ADM formalism and York’s conformal transverse traceless procedure 11 2.3 Black hole horizons and asymptotic flatness . . . . . . . . . . . . . . 13 2.4 Final formulation of the problem . . . . . . . . . . . . . . . . . . . . 16 3 The Tools . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3.1 The finite element method . . . . . . . . . . . . . . . . . . . . . . . . 20 3.2 Adaptive polygonisation of parametric curves . . . . . . . . . . . . . 31 3.3 The boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . 34 4 The Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 4.1 Convergence tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 iii 4.1.1 Convergence of independent residuals . . . . . . . . . . . . . 40 4.2 Mass and radiation analysis . . . . . . . . . . . . . . . . . . . . . . . 46 4.2.1 Ellipsoidal marginally trapped surfaces . . . . . . . . . . . . . 54 4.2.2 Head-on-collision-like data . . . . . . . . . . . . . . . . . . . . 63 4.2.3 Hemispherical shells . . . . . . . . . . . . . . . . . . . . . . . 67 4.2.4 Cylindrical marginally trapped surfaces . . . . . . . . . . . . 75 4.2.5 The sphere-del-torus . . . . . . . . . . . . . . . . . . . . . . . 83 5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 5.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 5.2 Future directions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 iv List of Tables 4.1 Linear regression of independent residuals . . . . . . . . . . . . . . . 44 4.2 Coefficients of r = 0.5 ellipsoid mass model as a function of l/r. . . . 55 v List of Figures 2.1 Typical problem domain. . . . . . . . . . . . . . . . . . . . . . . . . 17 3.1 Uniform discretization of domain. . . . . . . . . . . . . . . . . . . . . 24 3.2 Non-structured discretization of domain. . . . . . . . . . . . . . . . . 25 3.3 Illustration of adaptive polygonisation. . . . . . . . . . . . . . . . . . 33 3.4 Vertices needed for derivatives and step sizes. . . . . . . . . . . . . . 36 4.1 Independent residual analysis of sphere-del-torus initial data with r1 = 5.3, r2 = 1.325, r3 = 2.45, t2 = 0.707pi, t1 = sin −1 ( r1−r3 r1 sin t2 ) , and r0 = 600 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 4.2 Convergence of masses with domain size for collision-like data for r1 = 0.2, r2 = 0.4, z1 = 0.25, z2 = −0.5, and w = 2.2. . . . . . . . . . 51 4.3 Convergence of ADM mass with decreasing grid size for χ = 2−3, r = 0.13 andρ = 1.0 rounded cylinder initial data. . . . . . . . . . . . 52 4.4 Convergence of horizon mass with gridsize for χ = 2−3, r = 0.13, and ρ = 1.0 rounded cylinder initial data. . . . . . . . . . . . . . . . . . . 53 4.5 Mass measurements of r = 0.5 ellipsoids with varying l/r. . . . . . . 57 4.6 Contour plot of l/r = 0.56, r = 0.5 ellipsoidal initial data. . . . . . . 58 4.7 Contour plot of Ψ = |∇ψ| for l/r = 0.56, r = 0.5 ellipsoidal initial data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 vi 4.8 Contour plot of l/r = 0.21, r = 0.5 ellipsoidal initial data. . . . . . . 60 4.9 Contour plot of l/r = 1.96, r = 0.5 ellipsoidal initial data. . . . . . . 61 4.10 Contour plot of Ψ = |∇ψ| for l/r = 0.21, r = 0.5 ellipsoidal initial data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 4.11 Contour map of collision like initial data with r1 = 0.2, r2 = 0.4, z1 = 2.7× 10−2, and z2 = r1z1/r2 . . . . . . . . . . . . . . . . . . . . 65 4.12 Mass estimates of pseudo-collisional surfaces with r1 = 0.2, r2 = 0.4, z1 = 0.25, z2 = z1 − d and w = 2.2. . . . . . . . . . . . . . . . . . . . 66 4.13 The generating function of the hemispherical shell. . . . . . . . . . . 68 4.14 Mass estimates for r = 0.53 hemispherical shells with varying χ. . . 69 4.15 The maximum radiation loss of r = 0.53 hemispherical shells with varying χ. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 4.16 A closeup of ∆E of r = 0.53 hemispherical shells. . . . . . . . . . . . 72 4.17 Contour plot of χ = 0.05, r = 0.53 hemispherical shell initial data . . 73 4.18 Contour plot of χ = 0.8, r = 0.53 hemispherical shell initial data. . . 74 4.19 Mass estimates for r = 0.53 superellipsoids . . . . . . . . . . . . . . . 76 4.20 Mass estimates for r = 0.53 cylinders with rounded corners. . . . . . 77 4.21 The maximum radiation loss of r = 0.53 superellipsoids. . . . . . . . 78 4.22 The maximum radiation loss for r = 0.53 cylinders with rounded corners. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 4.23 Contour plot of χ = 0.15, r = 0.53 rounded cylinder initial data. . . 80 4.24 Contour plot of χ = 0.15, r = 0.53 rounded cylinder initial data gradient . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 4.25 Contour plot of p = 0.15, r = 0.53 superelliptical initial data . . . . . 82 4.26 The circles creating the initial sphere-del-torus. . . . . . . . . . . . . 83 vii 4.27 The arcs removed from the sphere-del-torus. . . . . . . . . . . . . . . 84 4.28 The bezier curves completing the sphere-del-torus. . . . . . . . . . . 84 viii Acknowledgements I would like to acknowledge the academic support of my supervisor, Matthew Choptuik, in helping me understand the problem in clear and simple terms and for his ability to demonstrate the ease and utility of numerical methods, Bill Unruh for his meticulous and expeditious help with the final editing, the Argelander Institut für Astronomie and particularly Peter Schneider for providing me with the opportunity to work without sacrificing my family, and finally the financial support of NSERC. Aaryn Tonita The University of British Columbia May 2009 ix Dedicated to my mother, for sticking by me during some tough times. x Chapter 1 Introduction General relativity is the modern theory of the gravitational field. Black holes are a natural consequence of this theory and provide a model for studying the dynamics of the gravitational field alone. In order to study the dynamics of the gravitational field, the initial state must be specified. This data must in general conform to a set of constraint equations which complicate the specification of initial states. This thesis generates initial data for axisymmetric spacetimes which contain non-spinning black holes by specifying that the initial data contain a marginally trapped surface with a specified “shape” in coordinate space. This is done by solving the conformally decomposed constraint equations subject to initial time symmetry. Numerical methods are used to solve these equations, namely the finite element method. The following subsection provides the background and motivation for this problem and the methods used to solve it. 1 1.1 Motivation for studying the problem Numerical analysis seeks to solve those problems which are intractable by other forms of analysis. Specifically, numerical relativity seeks to solve problems in general relativity in the dynamic and strong field regimes where other forms of analysis generally fail. The most natural way to accomplish this is to first specify the initial state of the gravitational field and then use the laws of physics to describe its subsequent evolution. This defines the initial value problem of general relativity. There are various ways to specify initial data and arrive at physically rea- sonable spacetimes that satisfy the Einstein equations, for instance characteristic approaches [26] which specifies initial data on a null hypersurface and Cauchy for- mulations (such as the 3+1 decomposition) where the initial data is specified on a spacelike hypersurface [10] [9]. Early and influential Cauchy formulations were introduced by Arnowitt, Deser and Misner [20] and Choquet-Bruhat [8] and is sim- ilar to other Cauchy initial value problems in that it is assumed that a field can be known for all space at some initial time. The ADM formulation yields a system of second order hyperbolic evolution equations and four elliptic constraint equations. In this approach, the physical variables chosen to represent the gravitational field in the ADM picture are, roughly, a spatial metric for each spatial hypersurface and the time derivative of that metric. The extrinsic curvature is used as a natural expres- sion of the instantaneous time derivative of the gravitational field on the Cauchy slice [25]. Counting the components of the spatial metric and the extrinsic curva- ture, there are twelve unknown quantities of which four are fixed by the constraint equations and another four are arbitrary buy fixed by any particular choice of coor- dinates leaving four free quantities but no unique way of knowing which quantities should be specified to arrive at a particular and distinct physical spacetime. 2 This problem was studied by Lichnerowicz [19] and then later in detail by York [28] who developed what is now known as the conformal transverse traceless decomposition. Later, the Thin-Sandwich decomposition was also developed to give an approach which connected the extrinsic curvature with certain physical quantities [30]. The conformal decomposition treats the initial spatial metric as a conformal transformation of a chosen background metric, and uses conformal formulas to easily calculate derivative terms. This corresponds to a choice of five dynamical variables and simultaneously a coordinate choice on the initial spatial hypersurface. The conformal factor is then used to satisfy the Hamiltonian constraint. Due to the simplicity of conformal transverse traceless decomposition constraint equations for time symmetric initial data, this thesis uses this approach for decomposing the constraint equations and choosing the free data. A point of time symmetry is a time where for t′ = −t one has for any quantity u, ∂u/∂t = ∂u/∂t′. Asymmetric points are all other points. The reason for producing time symmetric initial data is essentially due to the automatic satisfaction of three of the four constraint equations and a simpler remaining constraint. For these reasons, time symmetric initial data is produced in this thesis. Black holes play a prominent role in the field of general relativity as the physical manifestation of gravitational fields so strong that the pressures supplied by the other forces of nature were not enough to prevent the collapse of a star into a singularity. Because of this fact, the fact that they can exist without the presence of further mass distributions, their significance as sources of gravitational radiation, and simple curiosity, their dynamic behaviour has been studied for decades. These reasons motivate this study of distorted black holes; in order to study directly their behavior in general situations. By definition, a black hole is a region in spacetime 3 from which no particle can escape to infinity, even photons; technically speaking it is the boundary of the causal past of future null infinity [25][16]. This definition demands knowledge of the entire spacetime in order to determine if a black hole is present or not. Knowing only the metric and its time derivative at a single instant of time, it is impossible to know exactly where the boundary of the interior of a black hole occurs. However, due a theorem by Penrose [22] it is known that initial data possessing a marginally trapped surface or an apparent horizon will develop singularities and under certain assumptions a black hole. Under general circumstances it is necessary to assume cosmic censorship, that all singularities in General Relativity are hidden behind event horizons, however, for the conformally flat, time symmetric initial data studied in this thesis it was proven by Jang and Wald [18] that the singularity would be hidden by an event horizon, and so the spacetime will contain a black hole. Since a black hole is a surface out of which no physical influence, even light, can propagate, one had the intuition that the physical data on the interior of the horizon could be ignored. This has led to the development of excision techniques, where the spacetime whithin the horizon (in fact the trapped surface) is ignored. While it is true that no physical information can escape the horizon, coordinate information can escape with no limit to the speed of propagation. Thus an alternative approach to excision was developed: the puncture technique. Coupled to gauge conditions which rapidly ensure that spacetime region outside the horizon is populated by coordinates. Further, initial data is described for any variable u in the system as u = ubounded + usingular (1.1) and a particular expression for usingular is used so that ubounded is genuinely bounded and is subject to modified constraint equations [5]. While the puncture method has 4 advantages, we are primarily interested in looking at the effects of changing the shape of the apparent horizon on the initial data for the black hole. This is much easier to do in the “excision” picture, and therefore the excision approach is used in this thesis. The marginally trapped surface has now been extensively used in numerical relativity to excise the singular regions of black holes from the numerical evolution experiments. A quasi-equilibrium (a Cauchy slice is said to be in quasi-equilibrium if some of its instantaneous time derivatives vanish, which ones and how many differ by author and method) method for generating initial data with arbitrary spin black holes has been developed and tested by Cook and Pfeiffer [12][11]. In the applications of these initial data they used apparent horizons which were coordinate spheres even though the boundary condition defined in these papers enables one to use arbitrary shapes of the excised region: so long as a conformal transformation is found which takes the actual shape of the surface to a sphere, a Killing vector representing the spin of the black hole in the spherical coordinates can then be transformed back onto the deformed surface as a conformal Killing vector which is all that is required to specify the spin for this type of initial data. Because of the difficulty in computing conformal transformations for arbitrary distortions, the distorted black hole initial data produced in this thesis is non-spinning. Some previous studies have been carried out on the dynamics of single, dis- torted black hole spacetimes. The initial data in Bernstein et al. [3] was specified in a way so that the spacetime contains two asymptotically flat universes connected by a coordinate sphere throat superimposed on a Brill type gravitational wave back- ground. The distortion in this spacetimes is generated by the shape and amplitude of the gravitational wave and is aimed at studying the interaction between black holes 5 and gravitational waves. Some static solutions have been investigated by Geroch and Hartle [15] where the black hole was distorted by external matter sources. Husa [17] studied black hole spacetimes generated by specifying a toroidal marginally trapped surface. Husa [17] looked at some particular distorted marginally trapped surfaces. To create initial data representing distorted black holes, Teukolsky [24] use perturbation theory to generate small deviations away from known, highly symmet- ric solutions. For strong distortions, with the hope of creating strong gravitational wave signals, perturbation theory cannot be used. The number of multipole mo- ments needed to be computed for some initial data is simply too large. For this reason numerical methods are used in this thesis to study distorted black hole ini- tial data. The preceding has outlined the mathematical problem to be solved in this thesis, namely the problem of producing initial data for black hole spacetimes using the conformal decomposition of the constraints, choice of time symmetry, and the use of a marginally trapped surface to ensure the space time will be black hole initial data. The remaining problem to solve is still quite difficult, namely because of the arbitrary shape that will be used for the marginally trapped surface. To solve this problem, this thesis uses numerical methods. The most common numerical method implemented to solve a partial differ- ential equation in general is the finite difference method, a numerical method of solving partial differential equations by approximating derivatives with differences between nearby points. This method becomes more complicated when the domain is irregular because the grid points of uniform square meshes will not generally intersect curved boundaries and finite difference formulae become complicated for non-uniform meshes, meshes which can be easily adapted to a curved domain. The 6 boundary conditions for uniform rectangular meshes become quite complicated [13] for curved domains. With a curved domain, the finite element method is the natural computational method of choice [4] for solving partial differential equations. The domain of distorted black hole initial data with an excised region will in general be curved and so this thesis uses the finite element method to solve the constraint equations. The finite element method casts the elliptic constraint equations into an in- tegral formulation known as the weak form by multiplying the equations by a test function and integrating by parts to transfer one derivative operator onto the test function. In this way the solution need only be once differentiable and the solu- tion is found by solving a variational problem: finding the function such that the integral equation is satisfied for all test functions. The domain is then discretised and piecewise polynomial basis functions are defined with compact support over this discretised domain. The function space spanned by these basis functions is a subspace of the space of once differentiable functions of which the actual solution is an element. By solving the variational problem in this limited subspace an approx- imation to the true solution is found, the approximation being the function which minimizes the distance to the actual solution in a norm defined by the weak form [4]. The first choice in specifying a finite element method is in choosing how the domain is discretised. The lowest order discretisation of a domain is a triangulation, and is used in this thesis due to the straightforwardness of triangulating a curved domain. The next choice of specifying a finite element method is in choosing the basis functions. For piecewise constant basis functions, each triangle becomes asso- ciated with a single degree of freedom (the coefficient of the basis function in the 7 linear combination of basis functions representing the solution) and the resulting approximation is only first order accurate in the average mesh size. For a piece- wise linear, continuous polynomial basis functions each node of the triangulation becomes associated with a coefficient of a basis function (usually simply the value of the approximation at the vertices) and the resulting method is second order. This thesis uses the second order accurate piecewise linear finite element method. Much research has been done on the finite element method. Of particular interest is the use of superconvergent gradient recovery techniques to estimate the a posteriori error associated with a given triangular subdomain [21]. In general a numerical method will have a given order of convergence, but in finite element solutions one can show that at certain points inside each element, the order of con- vergence is up to one order higer. This phenomenon is known as superconvergence and when it is used to estimate the gradient, it is called superconvergent gradi- ent recovery. The a posteriori errors this method yields are ideal for determining where a mesh needs refinement to locally achieve a desired accuracy. This type of refinement is used frequently, for instance in the publicly available finite element solver PLTMG [2]. This tool is ideal for the solution of the problem at hand: it uses adaptive mesh refinement on a triangularization of a general two dimensional domain to find a piecewise linear approximation to the smooth solution. Due to the public availability of PLTMG and its level of sophistication in mesh refinement, it is used in this thesis to study the problem at hand. An overview of this thesis is as follows: using the piecewise linear finite element method package PLTMG 9.0, initial data representing axisymmetric black hole spacetimes at a point of time symmetry is generated in the 3+1 decomposition of Einstein’s equations with a marginally trapped surface inner boundary. The theory 8 outlining the constraint equations and their formulation for this particular problem (Chapt. 2) as well as the theory of the finite element method and other numerical methods (Chapt. 3) is reviewed. In Chapt. 4 the initial data is tested for convergence and the masses and maximum radiation limits of various spacetimes is presented. None of the initial data produced permits more than 6% of the mass contained in the spacetime to be radiated, as computed through the difference between the horizon mass and the ADM mass. It is seen that some surfaces chosen as marginally trapped surfaces are not the apparent horizons (the outermost marginally trapped surface), another surface in these spacetimes surrounding the chosen inner boundary is the actual apparent horizon. In this thesis, axial symmetry, maximal slicing and time symmetry are used to decompose the constraint equations. A curve is used to generate a surface a revolution on which a boundary condition is imposed by forcing it to be a marginally trapped surface. A large radius is chosen at which the domain is truncated and a boundary condition is imposed on conformal metric here to enforce asymptotic flatness. Finally, the finite element method is used to solve the resulting boundary value problem for the conformal factor. Conclusions and the results of thesis are presented in Chapt. 5. 9 Chapter 2 The Physical Problem In this chapter the mathematical framework required to study the specific problem of generating axisymmetric initial data with specified marginally trapped surfaces is developed. The final section contains the exact mathematical statement of the initial value problem to be solved. 2.1 Conventions Units will be adopted such that G = c = 1 and metric signature (−,+,+,+) will be used. Greek indices {µ, ν} will denote four dimensional quantities, latin indices near the beginning of the alphabet {a,b,...} will be used to denote two di- mensional and indices {i,j,k} will represent three dimensional spatial quantities; spatial quantities are distinguished by orthogonality to the timelike vector, nµ de- scribed in the following section. Equations will always be written in terms of a coordinate basis. In some cases it will be necessary to distinguish between the three dimensional quantitiy and a two dimensional quantity and this will be accomplished by prepending a (2) superscript. 10 2.2 ADM formalism and York’s conformal transverse traceless procedure In order to foliate spacetime into a sequence of spacelike hypersurface a scalar function t is defined over the entire spacetime. The timelike vector associated with this coordinate is the vector field tµ which is tangent to the xi =constant curves. The infinite number of level sets of the time function define the hypersurfaces Σt. There is a unit, future directed, timelike vector orthogonal to these hypersurfaces, denoted nµ, related to the vector field of the time function as tµ = βµ+αnµ, where βµ and α are called the lapse vector and shift function [20]. In terms of these quantities the line element of the complete spacetime is ds2 = −(α2 − βiβi)dt2 + 2βidxidt+ hijdxidxj , (2.1) where hij denotes the spatial metric of the hypersurfaces. It describes the gravita- tional field on a hypersurface, and can also be written as hµν = gµν + nµnν , (2.2) which can be treated as a completely spatial tensor since it is of non maximal rank and so will now be written as hij . Initial data consists of two symmetric, spatial tensors hij and Kij , the ex- trinsic curvature, chosen to represent the first derivative of the gravitational field because of its behaviour under spatial coordinate transformations, in the 3+1 de- composition is is given by Kij = −1 2 Lnhij , (2.3) where L is the Lie derivative. This contains a total of twelve independent compo- nents. In general, a maximum of eight of these components can be specified and the 11 other four are fixed by the constraints. In general, it is not known what is the “physically relevant” way of choosing which of the twelve components are freely specifiable, it is mathematically arbitrary. Although choosing coordinates and the constraints themself in some way constrain some of them. Various useful ways of choosing these free quantities have been developed [29][30] which can be used to decouple the equations, both attributable largely to York. The first aspect of these methods is to assume that the metric of the spatial slice (hij) is conformally related to a freely chosen metric (h̄ij), hij = ψ 4h̄ij . (2.4) With this, one is able to easily write down the covariant derivative in terms of the co- variant derivative of the chosen metric and extra terms from the conformal factor, as well as to write down the Ricci scalar in terms of the same. The final choice involved in these methods involves treating the trace of the extrinsic curvature K = Kii as a freely specifiable quantity and the methods differ by how they accomplish this. The major difference between the conformal transverse traceless decomposition and the conformal thin sandwich approximation is how they decompose the conformal extrinsic curvature. Because of the assumptions made in this thesis, outlined below, the methods become equivalent. Here we choose the maximal slicing condition K = 0 and time symmetry so that Kij = 0, and a conformally flat spatial metric h̄ij = δij . The choice of maximal slicing allows the momentum constraint to decouple from the Hamiltonian constraint, while the choice of conformal flatness simplifies the derivative operators. The choice of flat metric fixes the coordinate degrees of freedom other than the conformal factor which is fixed by the Hamiltonian constraint and an overall unit of distance measure. The choice of maximal slicing and time symmetry fixes the 12 dynamical degrees of freedom to contain zero momentum and zero spin and the only remaining degree of freedom is a quadrupole and higher multipole moment which can be specified by the choice of marginally trapped surface. The only remaining constraint equation after these choices is the Hamiltonian constraint which is given by ∇̄2ψ + 1 8 K̄ijK̄ ijψ−7 = 0, (2.5) where ∇̄ represents the flat space derivative operator. These are the equations that conformally flat, maximally embedded, instan- taneously time symmetric initial data in the ADM-York picture must satisfy. Any initial space which satisfies these equations will generate a spacetime consistent with general relativity, but not necessarily a black hole space time. In order to ensure that there is a black hole within the spacetime, we need proper boundary conditions, which we describe next. 2.3 Black hole horizons and asymptotic flatness In order to have a black hole within our spatial slice we will excise a region and specify a boundary condition which guarantees the existence of a black hole. As noted in Sec. 1.1 this is merely one method for constructing black hole spacetimes chosen over the puncture technique because of its adaptability to a distorted horizon. The problem with this approach is that a black hole has a global definition: the region from which no light can escape to infinity. Since we are generating initial data, we have no way of tracing light rays to infinity to ensure that a surface is an event horizon. To ensure that the initial spatial slice describes a black hole spacetime, we impose that a freely specified surface be a marginally trapped surface. If the 13 spatial unit normal to the surface is denoted sµ then the condition that must be imposed on the marginally trapped surface is that the divergence of the outgoing null geodesics must vanish on the surface. Physically this means that neighboring null geodesics whose tangent is normal to the marginally trapped surface will converge to a point shortly after being continued off the marginally trapped surface [22]. Assuming cosmic censorship, this then ensures that the marginally trapped surface is contained in the event horizon. The outgoing null vector, kµ, orthogonal to the trapped surface can be written in terms of the spatial normal vector sµ and the timelike vector normal to the Cauchy surface, nµ, as kµ = nµ + sµ; (2.6) using this equation and the fact that gµν = hµν + nµnν , the condition that the surface be a marginally trapped surface, ∇µkµ = 0, becomes hµν∇µkν + (kµ − sµ)(kν − sν)∇µkν = 0. (2.7) Then kµ∇µkν = 0 because kµ is tangent to a null geodesic and kν∇µkν = ∇µ(kνkν)− kν∇µkν = −kν∇µkν = 0, (2.8) where the last equality follows since raising and lowering an index commutes with covariant differentiation. So we have, hµν∇µkν + sµsν∇µkν = hµν∇µ(nν + sν) + sµsν∇(µkν) = 0, (2.9) because of the symmetry of sµsν . Now note that Kµν = ∇(µnν) and sµsν∇µsν = 0, so that we have hµν∇(µnν) + hµν∇µsν + sµsν∇µnν = Kµµ +Dµsµ + sµsνKµν = 0, (2.10) 14 where Dµ is the spatial covariant derivative. Now, keeping in mind that all the above quatities have vanishing temporal components and that Kij = 0, s i = ψ−2s̄i and ∇isi = ∇̄isi+6si∇̄i lnψ where ∇̄i is the flat space derivative operator [25], the boundary condition imposed on the marginally trapped surface becomes ∇̄is̄i + 4s̄i∇̄i(lnψ) = 0, (2.11) in the flat coordinate system. This is the boundary condition prescribed on the inner boundary used to construct black hole initial data. We want to study a physical system in isolation, so at a sufficient distance from the black hole we want the space to appear flat. This means that the conformal factor of the previous section must approach unity (in “inertial coordinates”) as the distance from the centre approaches infinity. Perhaps more realistically for a physical system which will never be totally isolated, we can demand that the spacetime have the correct asymptotic behaviour, that is asymptotic flatness. Asymptotic flatness ensures that far from the origin, the conformal factor will have the asymptotic expansion in non-rotating non-boosted coordinates, [27] ψ = 1 + E 2r +O(r−2). (2.12) for some as yet, undetermined E. Differentiating this condition, solving for E and resubstituting the expression to eliminate E yields the boundary condition. ∂ψ ∂r + ψ − 1 r = O(r−3). (2.13) By treating the right hand side of the above equation as if it were zero and imposing this condition at a finite radius, it describes an approximate boundary condition to enforce asymptotic flatness. This ensures that sufficiently far from the black hole, the conformal factor drops off at the correct rate and approaches unity. 15 Satisfying the first equation on an inner boundary and the second equation on an outer boundary ensures that the initial data will represent an asymptotically flat spacetime containing a black hole. There is a certain freedom in choosing both of these surfaces, although the proper physical outer boundary condition in the inertial coordinates used here is a Dirichlet boundary condition specified at infinity; the Robin type (mixed) boundary condition of Eq. (2.13) is merely an approximation chosen for the computational simplicity of having an outer boundary at finite radius. A sphere or cube is perhaps the most natural choice for the outer boundary to ease computations, but for the marginally trapped surface we have freedom to use our intuition and imagination to construct black hole spacetimes which contain a marginally trapped surface with a specified shape in the conformally flat coordinates. Next we choose a coordinate system and simplify the problem, mathemati- cally and computationally, by imposing axisymmetry. 2.4 Final formulation of the problem Although fully three dimensional evolution problems are cutting edge, and initial data for these problems are routinely generated, they still present great chal- lenges. By imposing axisymmetry, one can reduce the complexity of the problem and still study problems of interest. For the purpose of computational simplicity, this thesis studies problems with this symmetry imposed. With these assumptions, we choose axisymmetric coordinates (r, z, φ), where ∂φ is a Killing vector and consequently, ∂φψ = 0. The computational domain then is Ω ⊂ R2 and there are three distinguishable parts of the boundary ∂Ω: the curve which generates the marginally trapped surface, ΓA, the axis of rotation, ΓR, and the boundary, Γ∞, where space is truncated and boundary conditions derived from 16 Figure 2.1: Typical problem domain. Shown is a typical domain for the initial value problem. The outer, semi-circular boundary, Γ∞, represents the truncation of the domain. The inner, irregular bound- ary, ΓA, represents the marginally trapped surface. The z-axis, ΓR, is used as the axis of symmetry. asymptotic flatness are prescribed, (see Fig. 2.1). If sa is the outward pointing unit normal to the domain then the equations which govern this physical problem are simply: ( ∂2r + 1 r ∂r + ∂ 2 z ) ψ = 0 on Ω, (2.14) sa∇aψ + ψ∇as a 4 = 0 on ΓA, (2.15) sa∇aψ + ψ − 1√ r2 + z2 = 0 on Γ∞, (2.16) sa∇aψ = 0 on ΓR. (2.17) As pointed out in personal correspondence by Choptuik, in analogy to the 17 (3+1) decomposition, one has on the 2-surface that ∇asb = ∇(asb) +∇[asb], (2.18) and the first term of the right hand side is merely the extrinsic curvature of the two surface, Hab ≡ ∇(asb). (2.19) Then taking the trace of the above derivative results in the antisymmetric part dropping out, yielding ∇asa = Haa = H. (2.20) Where the conformal background of the problem is merely euclidean space, the trace of the extrinsic curvature for a surface of revolution is known: the sum of the gaussian curvature of the generating curve, κγ and the rotational curvature, κφ. Useful expressions for these two quantities are κγ = dsa ds ta, (2.21) κφ = s(r) r , (2.22) where ta is the tangent vector associated with the curve ta = dγds , s represents a parameterisation by arc length so that tatbδab = 1, and s (r) is the r component of the normal vector. The marginally trapped surface boundary condition becomes sa∇aψ + ψ 4 ( dsa ds ta + s(r) r ) = 0. (2.23) Eqns. (2.14-2.17) together with Eq. (2.23) provide the full description of the problem of finding axisymmetric, time symmetric initial data in the York-ADM method for black holes with an initially specified non-spherical marginally trapped surface. The goal is to find solutions of interest and to study physical properties 18 like the mass and possible radiation content (both to be described in Chapter 4) of these snapshots of spacetime. 19 Chapter 3 The Tools Now that the physical problem has been specified in a precise mathematical language it remains only to solve the problem. Due to the potential complexity of the domain, it is only practical to seek numerical solutions to the equations as outlined in Sec. 2.4. A natural choice for a numerical solution of a differential equation on a non rectangular domain is the finite element method. It is this tool that will be used to solve the Hamiltonian constraint equations and seek solutions of the distorted black hole problem. 3.1 The finite element method As a simple example of how the finite element method works, consider the boundary value problem ∂2u(x) ∂x2 − g(x) = 0 (3.1) subject to u(0) = u(1) = 0. (3.2) 20 One multiplies by a test function v(x) satisfying the boundary conditions and inte- grates by parts to achieve ∫ 1 0 ( ∂u ∂x (x) ∂v ∂x (x) + g(x)v(x) ) dx = 0, (3.3) and the problem is then to find u(x) such that the above integral is satisfied for all possible v(x). Because this integral only requires the function be once differentiable, and because it is a variational problem, it is distinguished from the actual differential equation and called the weak form of the problem. One cannot evaluate this integral for every function which satisfies the boundary conditions because there are infinitely many and is so not amenable to computer automation, but if we look only in a finite dimensional subspace of the total function space we can find the solution to the corresponding discrete problem computationally. We construct a piecewise linear subspace in the following way: discretise the domain into n equal intervals of length h = 1/n and at the endpoint of interval j (for j = 1..n− 1) we define a “hat” basis function φj(x) =   0 for |jh− x| > h (x− (j − 1)h)/h for x < jh ((j + 1)h− x)/h for x > jh (3.4) and its first derivative dφj dx (x) =   0 for |jh− x| > h 1/h for x < jh −1/h for x > jh (3.5) which is of course undefined at the interfaces. The span of these basis functions define the subspace that the approximation will belong to. That is, we solve a modified variational problem where we search for the function uh in this subspace 21 such that the exact same intregral, Eq. 3.8, is satisfied for all test functions in this subspace. The solution and test functions are expanded in this basis as uh(x) = n−1∑ j=1 ujφj(x) (3.6) and v(x) = n−1∑ j=1 vjφj(x), (3.7) where the test functions now belong to the finite dimensional subspace. The integral is then evaluated and yields v · (Ahuh +G) = 0. (3.8) Above v and uh are R n−1 vectors whose j’th element is given by vj and uj respec- tively; Ah is a matrix whose elements are given by (Ah)ij = ∫ 1 0 dφi dx (x) dφj dx (x)dx, (3.9) it is a tridiagonal matrix with diagonal elements given by 2/h2 and off diagonal elements given by −1/h2; and G is a Rn−1 vector whose elements are given by Gj = ∫ 1 0 g(x)φj(x)dx, (3.10) and is in practice approximated using the midpoint rule as Gj = ∫ jh+h jh−h g(x)φj(x)dx ≈ g(jh). (3.11) Because the vector v is arbitrary, the term in brackets in Eq. 3.8 must vanish, by bringing the matrix product to the right hand side one acquires a simple tridiagonal matrix equation which one may recognize as being identical to the second order centred finite difference approximation to the differential equation. 22 This procedure outlines the finite element method in brief. One casts the problem in weak form, discretises the domain, defines a basis with a finite number of elements to approximate the function space, and then solves the system of coupled coefficients which results from evaluating the weak form. The rest of this section elucidates how to do this for the distorted black hole initial data problem. As shown in Fig. 2.1, the domain of the problem is made up of a half-disk with an arbitrarily curved boundary replacing some portion of the symmetry axis. Attempting to cover such a domain with a uniform finite difference grid results in the boundary generally intersecting the grid in between grid points, as shown in Fig. 3.1. The uniform grid points almost never coincide with the curved boundary and the use of accurate finite difference formulae becomes complicated. A natural way to solve this problem is to use the finite element method on an unstructured mesh adapted to the marginally trapped surface, as shown in Fig. 3.2. The following development is based upon the normal Ritz-Galerkin treatment of the finite element method for elliptic problems as outlined in [6] or [4]. To use the finite element method one first needs to cast the differential equation in its weak form. The weak form of a differential equation is an integral form that eliminates second derivatives. For notational simplicity we introduce a vector, a(ψ), defined by a(ψ) ≡   ∂ψ/∂r ∂ψ/∂z   , (3.12) and a scalar function, f , defined by f(ψ) ≡ 1 r ∂ψ ∂r . (3.13) Note that f(ψ) = a1(ψ)/r, this duplicate notation is used because it mimics stan- dard notation in weak forms and makes the input of these functions into PLTMG 23 Figure 3.1: Uniform discretization of domain. Shown is the uniform discretization of the domain shown in Fig. 2.1. The original boundaries are shown in red. It is seen that the grid points may lay close, but generally not on the boundary with the exception of the z-axis. This phenomenon causes various problems. 24 Figure 3.2: Non-structured discretization of domain. Shown is a non-structured discretization of the domain shown in Fig. 2.1 produced by PLTMG. The grid vertices now intersect the boundary exactly which is now treated as a polygonisation of the original smooth curves. simpler. We now use notation such that ∇ ≡   ∂r ∂z   . (3.14) Then, leaving s as the outward pointing unit normal to the domain, the differential Eq. (2.14) and boundary condition Eqs. ( 2.23, 2.16, 2.17) can be rewritten as −∇ · a(ψ)− f(ψ) = 0 in Ω, (3.15) s · a(ψ) = −ψ 4 ( dsa ds ta + s(r) r ) ≡ gA on ΓA(r, z, ψ), (3.16) s · a(ψ) = 1− ψ√ r2 + z2 ≡ g∞ on Γ∞(r, z, ψ), (3.17) s · a(ψ) = 0 ≡ gR(r, z, ψ) on ΓR. (3.18) The following function spaces are then introduced H1c(Ω) = {φ ∈ H1(Ω)|φ is continuous on ∂Ω} (3.19) 25 with H1(Ω) = { φ ∈ L2(Ω) ∣∣∣∣ ∂φ∂xi ∈ L2(Ω), i = 1, 2 } (3.20) and L2(Ω) is the space of L2 integrable functions on Ω. Then multiply Eq. (3.15) by v ∈ H1c integrate over Ω and use integration by parts to arrive at the following variational form of the differential equations called the weak form of Eqs. (3.15– 3.18) A(ψ, v) ≡ ∫ Ω (a(ψ) · ∇v + f(ψ) v ) dr dz − ∫ ∂Ω g(r, z, ψ) v ds = 0 (3.21) with g(r, z, ψ) defined as the piecewise function g(r, z, ψ) =   gA(r, z, ψ) (r, z) on ΓA g∞(r, z, ψ) (r, z) on Γ∞ gR(r, z, ψ) (r, z) on ΓR . (3.22) By constuction any function satisfying the original differential equations will also satisfy the variational form, but it is possible that some functions satisfying the weak form will not satisfy the differential equation. Clearly those solutions which are not twice differentiable will not satisfy the differential equation, thus the weak form is seen as a relaxation of the requirements on differentiability. The weak formulation of the axisymmetric, time symmetric distorted black hole initial data problem is to find ψ ∈ H1c such that for all v ∈ H1c we have A(ψ, v) = 0. The Sobolev space H1c is still infinite dimensional, and the problem is not yet amenable to a computational solution. The finite element method remedies this by looking at a finite dimensional subspaces of H1c, designated H1c,h (here h is an abstract quantity specifying the general scale of discretization of a particular domain), and defining the finite element problem as: find ψh ∈ H1c,h such that A(ψh, v) = 0 for all v ∈ H1c,h. It can be shown that the ψh found in this way is 26 the closest possible approximation to ψ in the chosen subspace, that is it minimizes |ψ − ψh| in an appropriate norm [6]. By converting an infinite dimensional vector space to a finite dimensional one, a numerical solution is easily found through linear algebra techniques. Now we need to define that subspace. On a uniform grid of a rectangular domain, the piecewise linear finite ele- ment method can be converted to a linear system of equations which is identical to the linear system of equations formed for finite difference method using stan- dard centred, second order difference formulae. Conceptually however, the finite element method approximates the solution on the entire domain while the finite difference method approximates the solution only at the grid points. The finite el- ement method approximates the solution ψ of the differential equation, while the finite difference method approximates the differential equation itself. The basis functions of the linear subspace H1c,h are chosen so they are non vanishing only on a few neighbouring elements; that is, the basis functions have compact support. In this way the integral Eq. (3.21) is converted into a sparse matrix system and typical sparse matrix methods can be used to quickly find the approximation. In order to achieve second order accuracy of the approximate solution (where order of accuracy denotes the power with which the error converges in relation to the average grid spacing), ψh, piecewise linear basis functions are used. A triangulation of the do- main is the most flexible and easily adaptable method for discretizing the domain in conjunction with piecewise linear elements. The basis functions will ensure that the approximation will be linear over the given triangular subdomains. The following sets are defined, the triangulation: Th, the set of edges: Eh and the set of vertices (nodes): Nh. Here, the label h will typically represent an appropriate measure of the grid spacing, such as maxEi∈Eh(Ei) or maxTi∈Th √∫ Ti dA. In what follows we 27 will assume that h = maxTi∈Th √∫ Ti dA. Because the domains are unstructured the individual triangles of the domain can only be implicitly defined. The edges are of course the individual line segments of the boundaries of the triangle and the nodes are the endpoints of the edges. The triangulation of the domain can be implicitly defined by the following properties 1. if Ti ∈ Th then Ti ⊂ Ω and ⋃ Ti = Ω, 2. if Ti 6= Tj ∈ Th and Ek ∈ Eh and nl ∈ Nh then Ti ∩ Tj = ∅ or Ti ∩ Tj = nl or Ti ∩ Tj = Ek. This essentially means that all the triangles cover the domain and that individual triangles overlap only by an edge or a vertex at most and is embodied by Fig 3.2. The inner curved boundary will be polygonated, so the domain is to be interpreted as an outer semicircle, rotational axis and polygonated inner boundary. The triangles on the outer boundary are not strictly speaking triangles since they will contain arcs of triangles in fact, but will be referred to as such in any case. Then define the particular subsets Th ⊃ T i = {T ∈ Th|ni ∈ Nh and ni ∈ T} (3.23) over which the basis functions will be continuous and nonvanishing and will vanish outside. These domains are merely the set of triangles which have ni as a vertex. The set of points contained in this domain is defined as Nh ⊃ N i = {n ∈ Nh|n ∈ T i}. (3.24) Then associated with this neighbourhood and central node ni the basis function φi : Ω 7→ R is defined by four characteristics: 28 1. if x 6∈ T i then φi(x) = 0, 2. if T ∈ T i then φi|T ∈ P1, the space of linear polynomials over T , 3. φi(ni) = 1 for ni ∈ N i, and 4. if n ∈ N i and n 6= ni then φi(n) = 0. This says in plain language 1. the basis function is vanishing outside of its associated domain (element), 2. the basis function is linear over each triangle in the element and merely con- tinuous across edges, 3. the basis function has a value of one on the central vertex of the element, 4. the bassis function takes a value of zero at all other vertices in the element. Where item 2 and 4 together with 3 imply that over each triangle in the domain, the basis function decreases linearly to zero. This also implies that the basis function will be vanishing on the outer boundary of its associated element. Everything is identical on the boundaries except for the curved outer boundary where the only difference is that the basis function will not vanish on this edge. With these definitions in place, any f ∈ H1c,h can be written as a linear combination of the basis functions. Defining N as the number of nodes in Nh, the specific linear combination of the approximation is ψh = N∑ i=1 ψiφi, (3.25) and more generally v = N∑ i=1 viφi, (3.26) 29 where the ψi and vi are coefficients. With this description in hand, the integral in Eq. (3.21) can be evaluated in terms of the basis functions and transformed into a system of linear equations in the unknown coefficients. We define vectors ψ and v ψ ∈ RN with 〈ei,ψ〉 = ψi, (3.27) v ∈ RN with 〈ei,v〉 = vi, (3.28) where the set {ei} are the Euclidean basis vectors of RN and 〈·, ·〉 is the Euclidean inner product on RN . Then the integral Eq. (3.21) is approximated by vTLψ = 0, (3.29) where the elements of the N ×N matrix L are determined by Lij = A(φi, φj), (3.30) where A(u, v) is the weak form functional given in Eq. 3.21. This then defines the piecewise linear finite element method on the bulk of Ω. For each ni ∈ ∂Ω one row of L is formed from an analogous approximation of a boundary condition equation. Then, due to the arbitrariness of v, the linear system to solve is Lψ = 0, (3.31) which has a unique solution so long as the rank of the matrix is maximal. By con- struction Lii = A(φi, φi) 6= 0 so the only way it could not be of maximal rank is if columns were not linearly independent. The linear independence of the columns of L follows directly from the linear independence of the basis functions φi. Therefore, the linear system has a unique solution. Note that, because nowhere have the lo- cal sizes of the elements been specified (only, implicitly, their maximum size), it is, 30 in principle, possible to have a mesh with any given triangulation. This will only result in a different form for L whose specific values are not discussed here. Fur- thermore, the triangulation of irregularly shaped domains can easily be performed in an adaptive way if the tesselation is allowed to be unstructured. This section has surveyed the necessary background to understand the piece- wise linear finite element method on triangular tesselations of the domain. There are many freely available programs of varying sophistication that are able to solve such problems. The program PLTMG 9.0 by Randolph Bank is one such program [2]. It solves general two dimensional elliptic problems using piecewise linear finite element methods with adaptive discretization. PLTMG solves the resulting linear system using the multigrid method. This thesis uses PLTMG to solve the initial value problem for axisymmetric black holes described above. 3.2 Adaptive polygonisation of parametric curves As described previously, our goal is to solve the hamiltonian constraint for the conformal factor on a domain that has an inner boundary which is an arbitrarily specified curve in the (r, z) plane. A key obstacle that arises in numerically dealing with such a curved boundary is in polygonising it. PLTMG needs a boundary specified either by line segments or arcs of circles, so the curved inner boundary must be polygonised in some way. To do this, a discrete number of points are chosen with which to represent this curve, assuming a parametrization. The easiest way to accomplish this is to simply choose a number of points and have them equally spaced in the curve’s parameter. However, this may cause problems because a general parameterisation of the curve will not have points equally spaced along the curve ∂Ω. Also, to correctly capture the shape of the curve the (open) polygon 31 must sample the curve more where the curvature is higher. For this reason the polygonisation should adapt and increase the number of points sampled at areas of high curvature. This section outlines one known method of adaptively sampling a parametric curve[14]. The curve is denoted γ : [0, 1] 7→ R2 and the parameter is denoted λ ∈ [0, 1]. Because γ(λ) ≡ (x1(λ), x2(λ)) is a point of R2, γ(λ1)− γ(λ2) is a Euclidean vector. The polygon will be defined by a sequence of parameters Λ = {λi ∈ [0, 1]|i, j ∈ N, ; i, j < Nc; if i < j then λi < λj}, (3.32) Nc will thus be the number of vertices in the polygon. The opening angle at a given vertex, γ(λi), is defined for i = 1, 2 . . . , Nc − 1 as θi = cos −1 ( 〈γ(λi−1)− γ(λi), γ(λi−1)− γ(λi)〉 |γ(λi−1)− γ(λi)| |γ(λi−1)− γ(λi)| ) , (3.33) with 〈·, ·〉 the Euclidean inner product in R2 and | · | its associated norm. This is evaluated using the natural inner product of the flat metric space because we are specifying the shape of the apparent horizon in the flat space and not the confor- mally related actual hypersurface. The definition is merely the dot product of two consecutive segments in the polygon, centred on a vertex under consideration: the angle measures the angle between two consecutive line segments in terms of positions on the curve. The polygonisation algorithm also needs a temporary set containing the current end points of the curve Σ = {σi ∈ [0, 1]|i, j ∈ N; i, j < Ne; if i < j then σi > σj}. (3.34) One controls the polygonisation by specifying a critical angle. For a given vertex there are vectors pointing to the previous vertex and the subsequent vertex, see Fig. 3.3. If the angle between these two vectors is greater than the critical angle 32 Figure 3.3: Illustration of adaptive polygonisation. Here, the algorithm proceeds from the left to the right as we adaptively poly- gonise a semi-circle (shown in blue) parameterised by distance along the z-axis: γ = ( √ 1− z2, z). One sees that the first midpoint was already added to the curve, this is therefore the second iteration of the algorithm loop. Fig. a: adaptive poly- gonisation always occurs between endpoints λNc and λNe . Fig. b: the midpoint in parameter space, λNc+1 has been added to the curve and the opening angle, θNc+1, between the two new line segments has been computed. Fig. c: the opening angle between the two line segments was too small so the endpoint was added to the set Σ and the other vertices renumbered. The next possible adaptive polygonisation will again be between λNc and λNe . (closer to pi), it is considered a straight line and information (the vertex) is not needed there. If the angle is smaller than the critical angle, it is considered to be curved and the vertex is kept. After choosing a critical opening angle, θc, the algorithm for producing the set proceeds as: 1. Begin with λ1 = 0 and σ1 = 1 and Nc, Ne = 1. 2. Consider the opening angle of the midpoint θNc+1 with λNc+1 = λNc+σNe 2 and let λNc+2 = σNe , if θNc+1 < θc then let σNe+1 = λNc+1 and Ne → Ne + 1, else let λNc+1 = σNe , Nc → Nc + 1 and Ne → Ne − 1. 3. If Ne = 0 stop, otherwise repeat step 2. The procedure illustrated in Fig. 3.3. This process will terminate once all the midpoints between the current vertices would have opening angles greater than the critical angle. That is to say, the process terminates when adding the midpoint 33 would be adding a vertex onto a curve which is essentially a “straight line” in terms of the critical angle. For θc sufficiently close to pi the straight line segments will provide a suitably accurate representation of the underlying curve. Note that it is necessary to choose θc < pi in order for the algorithm to terminate, since in the limit that θc → pi, Λ→ [0, 1]. With the polygonal representation of the curve fully adapted to the curva- ture, all that remains is to compute the actual curvature of the curve to evaluate the boundary condition. This all important step amounts to numerical differentiation on the polygonal curve. This differentiation does not necessarily have access to Λ al- though it must have access to Γ = {γi ∈ ∂Ω|γi = γ(λi), λi ∈ Λ}. The arclength can be approximated as the Euclidean distance between points or as integrated distance along the polygon. This process is described next. 3.3 The boundary conditions The regularity boundary condition, Eq. (2.17), and the asymptotic flatness condition, Eq. (2.16), are straightforward to implement exactly in the finite ele- ment scheme; one simply chooses boundary basis funcitons which automatically satisfy these conditions. However, approximating at the marginally trapped surface Eq. (2.23) is problematic since the equation depends heavily on the geometry of the surface. Specifically, it depends heavily on derivatives of the position of the curve γ(s). The function that must be computed is K(γ) ≡ 4gA(ψ) ψ = dsa ds ta + s(r) r . (3.35) where gA was the function defined in Eq. (3.16). There must be at least be C 1 smoothness along the entire polygonised curve to evaluate this equation. The normal 34 vector will be computed as a pi/2 rotation of the tangent vector and will have the same order of error. Due to interpolation theory the expected error between the curve and its representative polygon is O(|λi − λi+1|2), so a finite difference formulae to the same accuracy is desired in order to be consistent with the interior approximation. In order to find the tangent to O(h2) at a position x ∈ ΓA along the polygon, three vertices are necessary. The distances along the polygonised curve to the nearest two points are labeled in increasing order h1 and h2. These measure distance to the endpoints of the line segment containing γ(x), the point the derivative is desired at. On the first line segment of the polygon or on the second half of a general line segment the third vertex chosen wil come from the succeeding line segment. Similarly, on the final interval or on the first half of a general line segment the vertex chosen will come from the preceeding line segment. In both cases the distance along the polygon to this curve will be denoted h3. The first endpoint of the line segment will be denoted γ1 ≈ γ(x ± h1), the second endpoint γ2 ≈ γ(x ± h2) and the spare vertex γ3 ≈ γ(x ± h3); the plus or minus signs depend on the condition described above and illustrated in Fig. 3.4. In the case that γ3 comes from the succeeding endpoint, let h1 denote the distance to the point γ2, h2 to the point γ1 and h3 denote the distance to the third vertex, as shown in the first part of Fig. 3.4. In practice, these values are measured using the distance along the polygon which roughly corresponds to true arc length. 35 Figure 3.4: Vertices needed for derivatives and step sizes. In both figures λ increases to the right, and the derivative is to be computed at point γ(x). Fig. a shows the case where the nearest endpoint, γ1, has a greater value of λ than γ2; the third vertex, γ3 comes from the endpoint of the succeeding (higher λ) line segment. In Fig. b, γ(x) comes from the first half of a line segment; γ1 has a lower value of λ than γ2. Here γ3 must come from an even lower value of λ. In both cases h3 = h+ h1. For the case of the γ3 coming from the succeeding line segment the coefficients a = h3 − h2 (h3 − h1)(h1 + h2) , (3.36) b = − h1 + h3 (h3 + h2)(h1 + h2) , (3.37) c = h2 − h1 (h3 − h1)(h3 + h2) (3.38) are used to compute the tangent vector δ+s γ = dγ ds + e+ = aγ2 + bγ1 + cγ3. (3.39) The error term to leading order is again given as e+ = 1 6 d3g ds3 (h1h2 + h2h3 − h1h3), (3.40) where g : [0, 1] − R is any component of γ. In the case that γ3 comes from the preceeding line segment hi denotes the distance from x to γi (see Fig. 3.4 and the coefficients take the form a = − h2 − h3 (h1 − h3)(h1 + h2) , b = h1 + h3 (h3 + h2)(h1 + h2) , c = h2 − h1 (h1 − h3)(h3 + h2) , (3.41) 36 and the corresponding finite difference formula is δ−s γ = dγ ds + e− = aγ1 + bγ2 + cγ3, (3.42) with the leading order error term given by e− = 1 6 d3g ds3 (h1h2 + h2h3 − h1h3), (3.43) The normal and tangent vectors created in this way will be normalized to O(h2) because of their definition with respect to arc length, that is 〈δ±γ, δ±γ〉 = 1+O(h2). However, they can be normalized fully in order to ensure that dds〈δ±γ, δ±γ〉 = 0. The tangent vector is therefore defined as ta ≡ δ ± s |δ±s | , (3.44) using the standard Euclidean norm, to ensure precise normalization. For the derivative of the normal vector, it is neceessary to consider the actual description of the resulting vector in terms of the original curve. Denoting the curve as γ(s) = (γr(s), γz(s)) T and because the normal vector is a pi/2 rotation of the tangent vector, the derivative of the normal vector is therefore dsa ds =   d 2γz ds2 −d2γr ds2   . (3.45) To achieve an O(h2) finite difference formula at any point of the polygon, the formu- lae need four vertices, the endpoints of the containing line segment and the endpoint of both the preceeding and succeeding line segment (unless we are in the first or last line segment). If associated with the point γi is the signed distance si from x then the general finite difference formula for a function at four separated points is δ2sg(s)|s−1(x) = d2f ds2 + e2 = aγ1 + bγ2 + cγ3 + dγ4, (3.46) 37 with coefficients given by a = −2 s2 + s3 + s4 (s1 − s2)(s1 − s3)(s1 − s4) , (3.47) b = −2 s1 + s3 + s4 (s2 − s1)(s2 − s3)(s2 − s4) , (3.48) c = −2 s1 + s2 + s4 (s3 − s1)(s3 − s2)(s3 − s4) , (3.49) d = −2 s1 + s2 + s3 (s4 − s1)(s4 − s2)(s4 − s3) . (3.50) This formula can be used to calculate the derivative of the normal vector to second order by correctly inserting γz and γr for g in the above formula. The error of this formula to leading order is e2 = − 1 12 d4g ds4 (s1s2 + s1s3 + s1s4 + s2s3 + s2s4 + s3s4). (3.51) These are the formulae that can successfully be applied to a polygonal repre- sentation of a parametric curve to correctly account for the curvature and calculate the marginally trapped surface condition properly. The above formulae are used to create a program which automatically applies the correct boundary conditions for a given curve without knowing the closed form expression of its Gaussian curvature. As an added benefit, this feature can be used to study marginally trapped surfaces which are only known numerically; this was done for the collision-like initial data family(Sec. 4.2.2). 38 Chapter 4 The Results After finding initial data for axisymmetric distorted black holes, a few phys- ical properties of these systems can be examined. It is possible to place an upper bound on the energy of the gravitational radiation that could be emitted during evolution of the data. The mass of the black holes can also be analysed in rela- tion to their distortion. First, however it is necessary to analyse the initial data for convergence and accuracy of the numerical solutions. In the following section, the solutions are shown to converge, and evidence is given that strongly indicates that there is no error in the implementation. Physical properties of the systems are discussed in Sec. 4.2. 4.1 Convergence tests Convergence testing is a crucial aspect of computational methods and nu- merical analysis. The idea is to estimate the error of the approximate solution in relation to the continuous solution. In the context of the calculations in this thesis, there are basically two things that can go wrong 39 1. the computational method is not converging, or 2. the computational method is not converging towards the continuum solution. In the first case the computational method is ill posed or unstable. In the second case, there has been a mistake in the implementation of the solution to the discrete equations. With no knowledge of the continuous solution, testing for these errors is dif- ficult. However, there are ways around this difficulty, essentially by showing that a sequence of solutions on finer meshes represents a Cauchy sequence in an appro- priate norm. That is enough to show that the computational method is converging to a particular solution. By using an independent method to compute residuals derived from the differential equation, and showing convergence of these residuals, the existence of the limit to the continuous solution is justified. In effect, the com- puted approximate solution is directly shown to be a convergent representation of the continuum solution. 4.1.1 Convergence of independent residuals Using the original differential equation, Eq. (2.14), and some finite difference approximation (FDA) of it, the finite element solution is interpolated to a uniform rectangular grid and a residual calculated with respect to the FDA operator. The residual computed in this way is independent of the method used to find the approx- imate solution and provides proof via an alternate method that the correct solution is being obtained; it eliminates the possibility of a flaw in the implementation of the solution algorithm. Using operater notation, L = ∂2r + 1r∂r + ∂2z . Eq. (2.14) is rewritten as Lψ = 0. (4.1) 40 Then a corresponding second order finite difference operator is denoted by Ldiff and the approximate solution ψdiff to the approximate differential equation Ldiffψdiff = 0, (4.2) will have a Richardson expansion[7] ψdiff = ψ − h2ediff +O(h4) (4.3) where ediff is an h independent function of position and h is the spacing between two adjacent grid points. If the finite element solution is converging towards the continuous solution then it can be written that ψelem = ψ + eelem, (4.4) where e ∈ H1c and |eelem| tends towards zero as the discretization of the finite element method is refined. So using Eqs. (4.2,4.3), the application of Ldiff to the finite element approximation ψelem yields Ldiffψelem = Ldiff(ψ + eelem) = Ldiff(h2ediff +O(h4) + eelem) (4.5) The crucial assumption for the efficacy of the use of independent residuals to demon- strate convergence in this case is that Ldiffeelem ≪ h2Ldiffediff (4.6) i.e. that the error of the approximation is negligible even after being operated on by Ldiff and the measured residual is dominated by the error of the finite difference operator. If this were not the case the residual would not converge in the way expected for the finite difference operator since it would have a contribution based on the unknown error coming from the approximation. In order to ensure that this 41 assumption holds, the approximate solution used for this analysis, the approximate solution is taken to have as fine a mesh as possible and the finite difference mesh is taken to be a few orders of magnitude larger. When this condition applies, the residual becomes r ≡ Ldiffψelem ≈ h2Ldiffediff , (4.7) where r is the independent residual vector. If there are η nodes in the uniform mesh discretization of Ω then r ∈ Rη and using the norm (|v| = √ 1/n ∑ v2i ) to compute the average residual gives r = |r| = √ 1 η ∑ (Ldiffediffh2)2 < h2 supLdiffediff = kh2 (4.8) for constant k = sup (Ldiffediff) which is also O(1) so long as ediff is bounded and smooth. In order for the foregoing derivation to be accurate, it is crucial that Eqn. (4.6) hold. Because Ldiff is an approximate differential operator and eelem is only piece- wise smooth, this condition is non trivial. Standard a priori error estimates show that ∇ψelem converges to first order but ∇2ψelem is no longer square integrable and local measurements vanish or diverge. However, one expects that the real solution for the initial data problem studied in this thesis be smooth (at least more than just once differentiable). So the lack of twice differentiability can be considered an artifact of the finite element method, a similar phenomenon occurs in the evolution of the vacuum electromagnetic field using the finite element method. For these rea- sons one must use a reconstruction method in postprocessing ψelem if one wishes to analyse its second derivatives and here a bicubic interpolant is used: x = (r, z) ∈ Ω P (x) = 3∑ i=0 3−i∑ j=0 pijr izj . (4.9) 42 The ten coefficients {pij} are determined by requiring that P (x) be an exact in- terpolant for the vertices on the nearest ten nodes (not the neighbouring ten, but the ten points out of all points which are closest to this node). This interpolation smooths the piecewise linear approximation so that the second derivative of the approximate solution is no longer locally vanishing. The residual we compute is then r = LdiffP (x), (4.10) where P is defined separately for each x. With the numerical differentiation and the interpolation to the required points now defined, the independent residuals can be computed. On initial data whose outer boundary lies at radius (truncation radius) r0, a geometric sequence of step sizes, h(l), is used, where the parameter l ∈ R represents the discretization level: h(l) = r0 2 2−l. (4.11) Then two levels are compared simultaneously with d = l2−l1 > 0. The ratio between the components of the residual vectors at locations which coincide on both grids is calculated and averaged. The relationship between the average ratio of residuals at the two levels, if the above approximation was valid, will be e(d) = ∣∣∣∣r(l1)|xr(l2)|x ∣∣∣∣ = C (2−l1)2(2−l2)2 = C22d. (4.12) By taking logarithms of the above equation, the relation is linearised as ln(e(d)) ln(2) = ln(C) ln(2) + 2d, (4.13) so the normalized logarithm of the error ratio will scale linearly with the difference in the levels of discretization and have a slope of 2. With initial data representing 43 regression of slope All data 2.07± 0.01 Upper bound 2.41± 0.01 Table 4.1: Linear regression of independent residuals The results of linear regression applied to the independent residual data plotted in figure 4.1. The regressions performed are for the entire data set and the outlying upper points, defined as, for a fixed d, the point with the largest value (upper bound). The slope of 2 represents the expected scaling behaviour of the residuals under the assumption of small source error and second order finite differencing, as described by equation 4.13. The upper bound may be an outlying superconvergent branch. a sphere with a torus removed from the equator (Sec. 4.2.5), the convergence of the independent residuals is shown in Fig. 4.1 and the results of the regression are summarized in Table 4.1. The scatter in the plot is due to the ill posed nature of differentiation on scattered data points. For higher order finite element methods, each derivative of the approximate solution decreases the order of convergence by 1. Naively carrying this behavior over to the piecewise linear method, one would presume the error to be O(1), here it is remedied by taking as accurate a solution as possible and smoothing the result with an interpolation. A more sophisticated approach would be to use the superconvergent points mentioned in the introduction to compute an interpolation of higher order. However, the straightforward interpolation used here was seen to be sufficient for proof of convergence. Because the upper lying branch of points (denoted as the upper bound in the figure) in the plot of residuals seems to be outlying, its convergence rate was measured to determine if it was a set of badly converging points. Its rate of convergence was instead measured to be of higher order than expected and may be due to unexpected correlations with the 44 0 1 2 0 2 4 6 Figure 4.1: Independent residual analysis of sphere-del-torus initial data with r1 = 5.3, r2 = 1.325, r3 = 2.45, t2 = 0.707pi, t1 = sin −1 ( r1−r3 r1 sin t2 ) , and r0 = 600 Lines represent least squares fit to uppermost points (short dash) and all points (solid line). There is still some oscillation in the derivatives. The error shows second order dependence on the difference of discretization level, justifying the assumption of equation 4.6. The initial data represents a sphere with a torus removed, (Sec. 4.2.5), recall that r0 represents the truncation radius. 45 superconvergent points. Its exact origins are not necessary for showing the evidence of convergence and it is merely fit and displayed here due to its outlier quality. To summarize, using a higher order interpolation than that defined by the finite element basis functions, independent residuals were computed on a nontrivial family of initial data. It is seen that these residuals converge to second order in the discretization. This implies that the finite element solution as computed is in fact converging to the correct solution to the initial value problem. In the follow- ing sections, some of the physical properties of some initial data families will be presented. 4.2 Mass and radiation analysis For vacuum initial data describing black hole spacetimes, measures of the mass of the gravitational system are of specific interest. In the current case there are two such measures that are of particular interest. The first is the mass defined by Arnowitt, Deser and Misner when they first considered the problem of dynamics in general relativity [1]: this quantity is computed as a surface integral for asymp- totically flat spacetimes, such as those currently under consideration. Specifically, the ADM mass is given by MADM = lim ρ→∞ 1 16pi ∑ i,j ∮ [ ∂γij ∂xi − ∂γii ∂xj ] sjdA, (4.14) where sj is the outward pointing spatial normal to the surface on which the integral is computed. Computing this in an axisymmetric, conformally flat spacetime and 46 on a coordinate semi-circle one finds MADM = lim ρ→∞ MADM(ρ), (4.15) MADM(ρ) = −ρ2 ∫ pi 0 (ψr sin θ + ψz cos θ)ψ 5 cos θdθ, (4.16) where ρ is the radius of the coordinate sphere on which the mass is computed, θ is an azimuthal angle measured from the positive z axis in the axisymmetric coordinate system, and MADM(ρ) will be referred to as the truncated ADM mass. This mass is found below to monotonically decrease in a way which is inversely proportional to the radius at which it is measured: MADM(ρ) ≈M∗ + k ρ (4.17) for some unknown M∗ and k. The second mass measure comes from the the area of the marginally trapped surface. Marginally trapped surfaces are always internal to the apparent horizon (or coincident with it) and, assuming cosmic censorship, the apparent horizon is internal to the event horizon. Thus, at any instant of time, the area of a marginally trapped surface is smaller than the area of the intersection of the event horizon and the t = const. surface. Because the area of the event horizon is related to the mass of the black hole in a specific way, one can both bound the mass of the black hole and, in many cases, determine a reasonable estimate of the mass by using the area of the marginally trapped surface. The area of a surface of revolution in a conformally flat spacetime is given by A = 2pi ∫ r(s)ψ4(r(s), z(s))ds (4.18) where s represents an arclength parameterization of γ(s) ≡ (r(s), z(s)). This area 47 is then used to compute the “apparent horizon mass” MAH = √ A 24pi , (4.19) which is just the Hawking energy of the apparent horizon in a vacuum spacetime: the amount of mass-energy surrounded by the apparent horizon [23]. Jang and Wald [18] previously showed that for time symmetric initial data where the apparent horizon can be continuously deformed to a 2-sphere at infinity that the mass determined by the area of the apparent horizon is always less than the ADM mass. The two mass measures defined above clearly measure two different quanti- ties, with the horizon mass only being an approximation of the event horizon mass and we note that there are additional mass like quantities which may be computed but we do not[23]. After an evolution of the system, the ADM mass (at infinity, not truncated) of the system will be the same as the original ADM mass but the mass associated with the apparent horizon can increase as the black hole absorbs gravitational radiation. The difference between the two quantities upon reaching equilibrium in a vacuum spacetime can be interpreted as the total energy emitted by the system in the form of gravitational radiation. On an initial data slice, the dif- ference in the two quantities can be interpreted as the maximum amount of energy available for radiation[3]. Because the initial data is computed on a subdomain of the entire spacetime only the truncated ADM mass can be calculated and it is necessary to quantify how the measurement of the mass depends on the radius at which the computational domain is truncated. This relationship is shown in Fig. 4.2. The data is calculated with a fixed ratio of the average finite element edge length to the truncation radius. The two mass measurements, MADM andMAH, are computed approximately 48 using the trapezoid integration rule. The estimate of MADM is computed as MADM = −12ρ2 η∞−1∑ i=1 ψ(xi) 5 ρ−1 (ψr(xi)ri + ψz(xi)zi + ψr(xi+1)ri+1 + ψz(xi+1)zi+1)∆θ/2, (4.20) where xi ≡ (ri, zi) are the ordered vertices comprising the endpoints of the arcs defining the circular outer boundary Γ∞, η∞ is the number of such points and ρ is the radius of the outer boundary; ∆θ is computed as ∆θ = tan−1 ( ri+1 zi+1 ) − tan−1 ( ri zi ) . (4.21) The area of the outermost marginally trapped surface is computed as MAH = 2pi ηA−1∑ i=1 ( riψ(xi) 4 + ri+1ψ(xi+1) 4 ) ∆s 2 , (4.22) where xi ≡ (ri, zi) are the points actually contained in ΓA and chosen to approximate it, ηA is the number of such points, and ∆s is computed using the Euclidean norm (recall that xi ∈ R2): ∆s = |xi − xi−1| . (4.23) These mass measurements explictly depend on the size of the mesh as well as the radius of the outer boundary. Since the error in the solution due to the finite sized mesh is converging to second order and the gradient of the data converges to first order, these masses will converge to first order. This behaviour is explicitly shown in Fig. 4.3 and Fig. 4.4, which, respectively, show the linear dependence of the ADM mass and the horizon mass on the basic scale of discretization; the data used in these figures is derived from initial data where the apparent horizon is cylindrical with rounded edges as described in Sec. 4.2.4. The radius of curvature of the rounded edges is one eighth of the radius of the cylinder. 49 With a solid understanding of the limitations of the masses measured from the computational data, we can now explore the dependency of the mass on the “level of distortion” of the black hole. Each distortion is controlled through a di- mensionless quantity which is specific to the parameterisation of the curve; the quantity is generally unphysical as it can be changed by reparameterisation. How- ever, the general behaviour of certain families of initial data can be described for a given parameterisation. One must keep in mind all the sources of error that occur in the computation. When derivatives are involved in post processing, the error in the post processed quantities will generally not converge to zero at second order in the number of ver- tices. For the mass estimate derived from the area of the marginally trapped surface, uncertainty arises because the ∆s used in the approximation to the integral is only an approximation of the actual change in arc length between the two points and the error is roughly proportional to the stepsize. The accuracy of this approximation depends on the critical angle used in the polygonization procedure. Intrinsic to PLTMG is the possibility that the initial unstructured mesh will form triangles tiny in area but possessing large edge sizes; this seems to be correlated with a highly sampled marginally trapped surface. When these triangles occur, an unphysical ridge occurs along the long edges, linear behavior where is should be curved. Fi- nally, changing any parameter of the simulation will alter the entire mesh in an unpredictable fashion, producing these errors in a similarly unpredictable fashion, demonstrating a statistical dispersion about the true values of order h for any small change in a given parameter. 50 Figure 4.2: Convergence of masses with domain size for collision-like data for r1 = 0.2, r2 = 0.4, z1 = 0.25, z2 = −0.5, and w = 2.2. The ADM mass and the horizon mass (MADM and MAH respectively) calculated from head on collision-like initial data (see Sec. 4.2.2) and plotted against the radius at which the computational domain was truncated. The initial data used is meant to represent the post-collision state of two black holes, with a mass ratio of 2:1, after a common apparent horizon has formed. 51 0.002 0.004 0.006 0.008 1.01 1.02 1.03 Least squares fit. Figure 4.3: Convergence of ADM mass with decreasing grid size for χ = 2−3, r = 0.13 andρ = 1.0 rounded cylinder initial data. The ADM mass plotted versus 1/ √ Nh, the inverse of the square root of the total number of vertices, which is a measure of the average grid spacing. The initial data was the rounded cylinder initial data (see Sec. 4.2.4) with χ = 2−3, where χ is a parameter controlling how sharp the edges of the cylinder are. The data is fit using a least squares regression showing the first order dependence of the average edge length of the mass. The equation of the line isMADM = (−3.69±0.05)/ √ Nh+(1.03) with negligible error in the intercept at these significant figures. 52 0.002 0.004 0.006 0.008 0.402 0.404 0.406 0.408 Least squares fit. Figure 4.4: Convergence of horizon mass with gridsize for χ = 2−3, r = 0.13, and ρ = 1.0 rounded cylinder initial data. Same as Fig. 4.3 excep here the horizon mass is plotted versus 1/ √ Nt. The equation of this line is Mhorz = (−1.04± 0.01)/ √ Nh + (0.41079± 0.00005). 53 4.2.1 Ellipsoidal marginally trapped surfaces A simple, yet interesting class of initial data is created by choosing the ap- parent horizon to have the shape of an axisymmetric ellipsoid. This is represented in the computational region by half an ellipse and is generated from the parametric curve, γ(t) = (r sin t, l cos t) = r ( sin t, l r cos t ) . (4.24) Increasing r for a fixed value of l/r will increase the overall mass in a linear fashion and is equivalent to a coordinate rescaling. A convenient way to investigate the effects of distortion is to fix r and vary the dimensionless parameter l/r. The mass estimates for this family of initial data are plotted in Fig. 4.5. Plotted as well are least square fits for the ADM mass and the marginally trapped surface mass following the models MADM = aADM + l r bADM + r l cADM, (4.25) MAH = aAH + l r bAH + r l cAH, (4.26) where aADM, bADM, ..., cAH are coefficients whose estimated values are given in Table 4.2. Should these formulae remain valid as l/r → 0 and l/r → ∞ then the maximum radiation loss can be extrapolated to extreme values: lim l/r→0 ∆M ≈ 1− cAH cADM = (5.6± 0.3),% (4.27) lim l/r→∞ ∆M ≈ 1− bAH bADM = (4.4± 0.1).% (4.28) From the plot of the mass of the ellipsoids versus the l/r ratio, there is a minimum to the mass when the ratio of the semimajor axis to semiminor axis is about l/r = 0.56. As the ratio is decreased from unity the total volume of the surface decreases and this effect initially dominates the energy content, but after the 54 Mi ai bi ci MADM 0.32± 0.02 0.52± 0.01 0.14± 0.01 Mmts 0.34± 0.01 0.49± 0.01 0.14± 0.01 Table 4.2: Coefficients of r = 0.5 ellipsoid mass model as a function of l/r. The results of a two dimensional linear regression to fit MADM and Mmts for ini- tial data with ellipsoidal marginally trapped surface. The error terms bound the calculated coefficients in a statistical 95% confidence interval. critical value the energy required to form such a distorted marginally trapped surface increases the energy over the decrease in volume. From the contour map of the conformal factor of this initial data, Fig. 4.6, it is probable that the conformal factor of the initial data has its maximum forming off the z-axis at some (rmax, 0), and since a marginally trapped surface always contain now or in the future a singularity it may be that this maximum represents, in the continuum limit, a singularity with non-pointlike topology or a forming singularity which may maintain a non-pointlike topology. This is further corroborated by the contour map of the magnitude of the gradient of the conformal factor, Ψ = |∇ψ| = |∇iψ∇iψ|, (4.29) see Fig. 4.7. The contours of the gradient show a local minimum along the marginally trapped surface at the axis, indicating the formation of a local saddle point at this intersection. One may speculate that the central singularity associated with near- spherically-symmetric initial data has vanished. This phenomenon is even more evident from the contour plots for initial data with l/r = 0.21, where the gradient is actually decreasing in magnitude along the axis as it meets the marginally trapped surface (see Fig. 4.8 and Fig. 4.10). Fig. 4.9 illustrates typical initial data from the other region of parameter 55 space where the marginally trapped surface is an oblate ellipsoid characterized by l/r > 1. In this case, because of the behaviour of Ψ one may hypothesise the existence of two coaxial point-like singularities. To summarize the ellipsoidal initial data, it has been shown that a highly distorted marginally trapped surface can increase the energy in the spacetime rela- tive to an apparent horizon with constant curvature (a sphere) in the conformally flat spacetime. The behaviour of the two masses shows that the greatest increase in MAH relative to MADM occurs for increasing l/r; however, a study holding MADM constant and varying l/r would be a more stringent demonstration of this. The total gravitational energy that could possibly be emitted has been estimated to be bounded to a small fraction of the total energy. Evidence is seen that indicates the possibility of richer singularity structure in spacetimes with distorted marginally trapped surfaces; however, because the singular region is not part of the computa- tional domain, this work is unable to prove this point. 56 Figure 4.5: Mass measurements of r = 0.5 ellipsoids with varying l/r. This plot shows the ADMmass and the mass measurement of the marginally trapped surface plotted versus the axial to radial diameter ratio. The lines represent the least square fit to equations 4.25 and 4.26, there is no guarantee that these models are accurate representations of the true functions. Some radiation is estimated by these calculations even for l/r = 1 due to the use of a finite computational domain as well as intrinsic variability due to the truncation error. 57 Figure 4.6: Contour plot of l/r = 0.56, r = 0.5 ellipsoidal initial data. This initial data represents the minimum of the mass curve shown in figure 4.5. The peak in the conformal factor at the edge of the disc becomes increasingly prominent for smaller values of l/r, the extra energy which is stored in this region of the gravitational field begins to dominate the mass as the l/r ratio is decreased. 58 Figure 4.7: Contour plot of Ψ = |∇ψ| for l/r = 0.56, r = 0.5 ellipsoidal initial data. This initial data represents the minimum of the mass curve shown in figure 4.5. The peak in the conformal factor at the edge of the ellipsoid becomes increasingly prominent for smaller values of l/r. The value of the conformal factor in this region of the domain begins to increasingly dominate the mass measurement as the l/r ratio is decreased. 59 Figure 4.8: Contour plot of l/r = 0.21, r = 0.5 ellipsoidal initial data. The circular contours near the edge of the disc suggest that there may be a singular- ity near there if the initial data was continued into the marginally trapped surface’s interior. From this data alone it is impossible to tell. 60 Figure 4.9: Contour plot of l/r = 1.96, r = 0.5 ellipsoidal initial data. This initial data set represent a contour plot of the conformal factor for the greatest value of l/r plotted on the mass plot, figure 4.5. 61 Figure 4.10: Contour plot of Ψ = |∇ψ| for l/r = 0.21, r = 0.5 ellipsoidal initial data. The circular contours near the edge of the disc suggest that there may be a singular- ity near there if the initial data was continued into the marginally trapped surface’s interior. 62 4.2.2 Head-on-collision-like data We now consider a family of initial data meant to superficially look like a collision of two black holes at some time after a single apparent horizon associated with with final black hole has formed. As before, we assume that the initial data correspond to a moment of time symmetry. We first discuss how the marginally trapped surface is constructed. Adopting conformally flat, cylindrical coordinates (r, z), one finds that the conformal factor for the Schwarzschild metric is ψS = 1 + m 2 √ r2 + z2 . (4.30) Since the radius of the event horizon is related to the mass as RBH = m/2, this can also be expressed as ψCS = 1 + RBH√ r2 + z2 . (4.31) One then defines a function by summing two of these conformal factors with distinct coaxial centres z = z1 and z = z2: ω(r, z) = ( 1 2 + R1√ r2 + (z − z1)2 ) + ( 1 2 + R2√ r2 + (z − z2)2 ) . (4.32) The function ω has various level sets, that is collections of points (r, z) where ω(r, z) = w for some constant w called the weight. Because of the smoothness of ω, the level sets form a one parameter family of planar curves, that can be spec- ified to be maginally trapped surfaces. The factors of 1/2 are included in ω so that when z1 = z2 then ω becomes the conformal factor for a single Schwarzschild black hole with RBH = R1 + R2. The level sets of this function are used to produce marginally trapped surfaces. The motivation for this being that one might expect that the apparent horizons of two colliding black holes to slowly merge and using the more naive approach of superimposing two spheres would generate creases. 63 In order to keep the data somwhat centred, instead of choosing both z1 and z2, one may choose z1 and then solve for z2 by computing the Newtonian centre of mass associated with two point masses and having this centre of mass lie on the origin. That is, we choose z1 and z2 is given by the relation z1R1 + z2R2 = 0. (4.33) If the separation of the centers, d = |z1 − z2| for a given weight is too large, then the level set will be two disjoint curves. There is also a level set before separa- tion where a point lying on the axis of symmetry has a discontinuous tangent. This point has a singular Gaussian curvature, and consequently the boundary condition is singular there. This level set, where the gaussian curvature becomes singular, is avoided to prevent computational difficulties. In Fig. 4.11 a contour plot of a member of this family of initial data is shown, its parameters are described in its caption. The mass measurements are plotted in Fig. 4.12. The study of this family of initial data is hampered by the fact that the curve is only known implicitly which introduces further numerical errors into the solution process. This is seen in Fig. 4.12. Much like the ellipsoidal marginally trapped sur- face, this family of initial data points to the possibility of richer singularity structure. The large statistical variance in mass measurements could possibly be solved with a smoother interpolation algorithms since these problems do not occur for the curves without numerical error. 64 Figure 4.11: Contour map of collision like initial data with r1 = 0.2, r2 = 0.4, z1 = 2.7× 10−2, and z2 = r1z1/r2 . This figure shows a zoom in onto the marginally trapped surface formed from the level set of two confomal factors. The position of the second conformal factor is chosen to place the Newtonian centre of mass at the origin. and the masses have a 2 : 1 ratio, with the largest mass being 0.4 (m). The conformal factor achieves a global maximum on the surface near the smaller conformal factor’s origin. 65 0 0.1 0.2 0.46 0.48 0.5 0.52 0.54 0.56 d Figure 4.12: Mass estimates of pseudo-collisional surfaces with r1 = 0.2, r2 = 0.4, z1 = 0.25, z2 = z1 − d and w = 2.2. The parameter d represents the separation between the centers of the two conformal factors d = |z1 − z2|. The initial data used is meant to represent the post-collision state of two black holes, with a mass ratio of 2:1, after a common apparent horizon has formed. As the separation parameter approaches the value at which the level set will split, the mass sharply increases due to the increasing Gaussian curvature of the marginally trapped surface. One sees the O(h) dispersion with small changes of the parameter in this plot. 66 4.2.3 Hemispherical shells This section uses hemispherical shells as marginally trapped surfaces: hemi- spheres whose interiors have been removed. They are studied merely as an object of curiosity. The surfaces are characterized by an external radius, r0 and a thickness ratio χ; the edge of the shell is rounded off by a torus of diameter equal to the thickness of the shell. We first define two different curves parameterized by t ∈ [0, 1] γ0(r, t) = r (sin (pi/2t) ,− cos (pi/2t)) , (4.34) γ1(r, χ, t) = r (1− (1− χ)/2, 0) + r(1− χ)/2 (cos (pit) , sin (pit)) . (4.35) Note that γ0(r, 1) = γ1(r, χ, 0) and γ0(χr, 1) = γ1(r, χ, 1). The curve γ0 will be used to form both the inside (radius r) and outside (radius χr) of the shell, the curve γ1 is used to round the edge of the shell (radius (1 − χ)/2r). The arclength of each curve is S[γ0](r) = r pi 2 , (4.36) S[γ1](r, χ) = r 1− χ 2 pi, (4.37) and the total arclength of each curve is ST = S[γ0](r) + S[γ1](r, χ) + S[γ0](χr) = rpi. (4.38) Thus, the following curve is parameterized by arclength, γ(r, χ, tST ) =   γ0 ( r, tSTS[γ0](r) ) for 0 ≤ tST < S[γ0](r) γ1 ( r, χ, tST−S[γ0](r)S[γ1](r,χ) ) forS[γ0](r) ≤ t < S[γ0](r) + S[γ1](r, χ) γ1 ( χr, ST (1−t)S[γ0](χr) ) forST − S[γ0](χr) ≤ tST ≤ ST (4.39) 67 and is a continuous curve that begins and ends on the r = 0 axis. This curve generates the hemispherical shell marginally trapped surfaces studied in this section and displayed in Fig 4.13. Figure 4.13: The generating function of the hemispherical shell. The hemispherical shell generating curve: two concentric hemicircles capped by a semicircle. The outer radius, r can be factored out of the equations with a coordinate rescaling and has no impact on the physics and non-trivial variations in physics must be induced via changes in dimensionless parameters. The thickness parameter χ is the parameter of interest here. There are two singular values of this parameter when the radius of the various circles vanishes and the Gaussian curvature (1/r) blows up: when the thickness (1 − χ)r approaches zero the shell becomes infinitesimally thin and the radius of the circle at the end of the curve goes to zero, and as (1−χ) 68 approaches unity the radius of the interior circle goes to zero. Mass estimates for hemispherical shells with 0 < χ < 1 are plotted in Fig. 4.14. 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 Figure 4.14: Mass estimates for r = 0.53 hemispherical shells with varying χ. The domain represents the thickness of the shell in terms of the fraction of the total external radius. The MTS mass quickly overtakes the ADM surface, indicating the presence of an external apparent horizon and signaling the unfeasibility of black holes with this particular shape. Plots of the conformal factor, Fig. 4.17 and Fig. 4.18 show that a local and sometimes global minimum of the conformal factor forms on the interior of the shell. However, the plots of the maximum radiation loss show that at χ = 0.8 the mass estimated from the marginally trapped surface is greater than the ADM mass, 69 indicating that there is another marginally trapped surface enveloping the striking features of this initial data, which has not been searched for in this thesis. 70 0.2 0.4 0.6 0.8 -0.2 -0.1 0 Figure 4.15: The maximum radiation loss of r = 0.53 hemispherical shells with varying χ. The rapid drop of the maximum radiation loss (∆E = (MADM −MAH)/MADM) to negative values is indicative of the formation of an external apparent horizon which forms as the curvature on the inside of the hemispherical shell becomes singular. 71 0.2 0.4 0.6 0.8 0 0.01 0.02 0.03 Figure 4.16: A closeup of ∆E of r = 0.53 hemispherical shells. The domain here is truncated to show the range in which ∆E can be accurately computed before the formation of the seperate apparent horizon. It can be seen that hemispherical shells can produce a few percent of the total mass of the system in gravitational radiation. 72 Figure 4.17: Contour plot of χ = 0.05, r = 0.53 hemispherical shell initial data This initial data set represent a contour plot of the conformal factor for the smallest value of χ plotted on the mass plot, figure 4.14. 73 Figure 4.18: Contour plot of χ = 0.8, r = 0.53 hemispherical shell initial data. This initial data set represents a contour plot of the conformal factor for the largest value of χ plotted on the mass plot, Fig. 4.14. The negative value of the maximum radiation loss shows that another marginally trapped surface exists in this initial data set, so the interior region of this data will not necessarily be visible to observers at infinity. 74 4.2.4 Cylindrical marginally trapped surfaces There is a special class of marginally trapped surfaces which poses a special challenge: namely, the class of surfaces with kinks. In particular, we consider here two families of marginally trapped surfaces that both limit to a cylinder, i.e. to a curve that contains right angles in the computational domain. The first family is generated from the rotation of the superellipse γSE(t) = r(sin (pit)| sin (pit)|p−1, cos (pit)| cos (pit)|p−1), (4.40) for t ∈ [0, 1]. The parameter p controls the shape of the marginally trapped surface and creates a sphere when p = 1 and a cylinder when p = 0: for intermediate values superellipses are generated. The second class of marginally trapped surfaces is generated by a piecewise function which represents a cylinder of radius, r, whose edges have been rounded with a radius of curvature χr, with χ ∈ (0, 1]. When χ = 1 the marginally trapped surface is a sphere and when χ = 0 it is a cylinder. The masses for the superellipses and rounded cylinder family are plotted in Fig. 4.19 and Fig. 4.20 respectively. As the edges become sharper, a second marginally trapped surface forms outside the one used as the boundary. This marginally trapped surface is the true apparent horizon which would also have a lower area than the one specified as the inner boundary, subsequently the horizon mass computed from the boundary marginally trapped surface overtakes the ADM mass and the maximum radiation loss estimate becomes negative. This is evidenced in Figs. 4.21-4.22. 75 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 p Figure 4.19: Mass estimates for r = 0.53 superellipsoids The domain represents the power p parameterising the superellipsoids as in Eq. (4.40) The mass increases at a fairly slow rate until the marginally trapped surface is nearly cylindrical, the mass then increases more rapidly. 76 0.2 0.4 0.6 0.8 1.2 1.4 1.6 1.8 2 p Figure 4.20: Mass estimates for r = 0.53 cylinders with rounded corners. 77 0.2 0.4 0.6 0.8 1 0.025 0.03 0.035 0.04 0.045 p Figure 4.21: The maximum radiation loss of r = 0.53 superellipsoids. The maximum radiation loss reaches a peak at a value of about 4.9% of the total mass of the system, and then the mass estimated from the marginally trapped surface drops as presumably an apparent horizon forms elsewhere. 78 0.2 0.4 0.6 0.8 0.025 0.03 0.035 0.04 0.045 p Figure 4.22: The maximum radiation loss for r = 0.53 cylinders with rounded corners. 79 Figure 4.23: Contour plot of χ = 0.15, r = 0.53 rounded cylinder initial data. This initial data set represents a contour plot of the conformal factor of cylinder whose edges are rounded by a torus of radius χ = 0.15 the radius of the cylinder. 80 Figure 4.24: Contour plot of χ = 0.15, r = 0.53 rounded cylinder initial data gradient This initial data set represents a contour plot of the magnitude of the gradient of the conformal factor of a cylinder whose edges are rounded by a torus of radius χ = 0.15 the radius of the cylinder. 81 Figure 4.25: Contour plot of p = 0.15, r = 0.53 superelliptical initial data This initial data set represents a contour plot of the conformal factor of a p = 0.15 superellipse. 82 r r r 2 3 1 Figure 4.26: The circles creating the initial sphere-del-torus. The sphere-del-torus generating curve begins as these circles with the shown dimen- sions: an outer circle of radius r1 an inner circle of radius r2 located a distance r3 away from the axis of rotation. 4.2.5 The sphere-del-torus A final example of a non-spherical marginally trapped surface is the sphere with a torus removed from it. One begins with a semi-circle of radius r1 centered at the orgin and a circle of radius r2 laying at (r3, 0) in the r > 0 halfplane. This is diagrammed in Fig. 4.26. One chooses angles t1 and t2; measuring from the r- axis towards the z-axis with origin laying at the center of the respective circles, one removes the arc from t1 to −t1 from the circle with radius r2 and arc from t2 to −t2 from the circle with radius r1; this is diagrammed in Fig. 4.27. One then connects these two broken circles with a Bezier polynomial whose extra control points are located a distance r = r1− r2− r3 along the tangent to the circles at the endpoints, this is diagrammed in Fig. 4.28. The sphere-del-torus is the volume of revolution of this curve. Initial data generated from this curve has the unfortunate property of con- 83 tt 2 1 Figure 4.27: The arcs removed from the sphere-del-torus. The arcs have been removed from the two circles for the second step in constructing the sphere-del-torus generating curve. An arc of 2t1 is removed from the outer circle and an arc of 2t2 is removed from the inner circle. Figure 4.28: The bezier curves completing the sphere-del-torus. The two circles are complete by Bezier curves so that the entire sphere-del-torus curve is continuous and differentiable. 84 taining negative values of the conformal factor. This occurs because the inner circle of the sphere-del-torus curve has a negative divergence of its normal vector. This negative divergence forces the slope of the conformal factor to be positive in the or- thogonal direction to the marginally trapped surface and the conformal factor must decrease as one approaches the marginally trapped surface. Calculations show that the marginally trapped surface mass estimate is greater than the ADM mass; as discussed previously, this indicates that this marginallly trapped surface is not an apparent horizon. For these reasons, physical analysis of this family of initial data is not reported here and this section is included only to describe the non-trivial family of initial data used to compute the independent residuals discussed in Sec. 4.1.1; 85 Chapter 5 Discussion 5.1 Conclusions Use of the marginally trappped surface equation as an inner boundary con- dition makes it possible to study a large class of initial data. The finite element method can be used to cope with the curved geometries of arbitrary marginally trapped surfaces with arbitrary coordinate shapes as has been shown. Axisymmetric initial data representing single black hole spacetimes has been created such that a few percent of the total energy of the system is conceivably available for radiation. Because the physical energy available for radiation is always related to the scale of the system, the total energy can be increased by increasing the overall size of the system. The relative energy available for radiation for a given mass must be related to dimensionless parameters and this behavior has been studied for various families of initial data. For one case where the behaviour of the mass function was modelled - the ellipsoidal initial data - the extrapolations to the extreme values suggested that still only a few percent of the total energy would be 86 available for radiation, comparable with previous initial data families [3]. In terms of the analysis presented in this thesis, the main drawback in using the finite element method with a non- uniform mesh arises when the solution data is numerically differentiated during postprocessing, especially when the curve is not smooth. In this case the second derivatives will in general have an O(1) statistical dispersion for small parameter changes. A smoothing method needs to be developed in order to reduce this source of error in post-processing: although high order inter- polation somewhat decreases the variance of the error in numerical differentiation, it does not completely eliminate it. This was observed while studying the independent residuals, because the lower order interpolations yielded residuals consistent with noise, they were not displayed. As first mentioned in Sec. 4.2, due to the theorem by Jang [18] the initial data created in this thesis must have Mmts < MADM for Mmts measured on the apparent horizon. As the Gaussian curvature of a maginally trapped suface approaches a singular value, the measured value of Mmts surpasses MADM indicating another marginally trapped surface, the actual apparent horizon, forms outside the first. These apparent horizons were not searched for since the goal of this thesis was to study the range of spacetimes where the apparent horizon could be specified. This behaviour is independent of the family of curves which limits to the singular surface. This was seen in the two families of surfaces with cylindrical limits as well as the hemispherical shell initial data. It was not seen for the ellipsoidal data but smaller values of l/r could show this phenomenon. This thesis has used the apparent horizon boundary condition to create black hole spacetimes which contain a marginally trapped surface with a chosen geometry at a given snapshot of time. It has been shown that this method gives only limited 87 control over specifying the spacetime as sometimes an outer marginally trapped sur- face will exist; one cannot use this method and guarantee that the apparent horizon will have a specific geometry. Initial data representing various shapes have been cre- ated and analysed. This thesis has demonstrated the first usage of the marginally trapped surface condition to generate initial data for black hole spacetimes from generic axisymmetric surfaces. 5.2 Future directions Using the quasiequilibrium framework for blackhole spacetimes [11] or exci- sion boundary conditions[12] this work could be extended to include rotating space- times. New software would be necessary to solve the momentum constraint in addi- tion to the Hamiltonian constraint solved here because PLTMG is unable to solve systems of equations. In addition, the basic idea used here could be extended to more general deformations of the marginally trapped surface if a three dimensional finite element solver was available. If the generating curve of the marginally trapped surface was expressed in terms of the axisymmetric spherical harmonics, the exact form of those basis func- tions could be used to yield exact boundary conditions on the marginally trapped surface as opposed to the approximate finite difference curvatures used here. One could also use a spectral decomposition of the boundary in terms of the spherical harmonics if the exact expression for the curve was not known. This would elim- inate any error coming from the implementation of the boundary condition which currently is controlled by a parameter independent of the discretization level. It is always possible to transform a distorted object into a sphere, so long as the topology is the same. Therefore, for conformally flat initial data with a given 88 shaped apparent horizon, there would always be a coordinate system where the same initial data on the same hypersurface would have the apparent horizon of a sphere. The difference would arise in the derivative operators due to the conformal metric no longer being flat. The necessity of using the finite element method would be eliminated and high accuracy spectral type methods could be used in that case. The resulting equations would be much more complicated in this approach. The resulting spacetimes would be the same however; indeed, the only real difference that can be made in the choice of the conformal metric is by making a suitable choice to eliminate “junk” radiation. Finally, if the spacetimes were evolved, the amount of gravitational radiation and its form could be studied. Furthermore, by studying the lightlike geodesics of the system originating from behind the black hole and propagating forward to an observer, one could study the effect of distortion on visible matter. Evolving the spacetime would need a dynamic finite element program or a better analysis of the error associated with applying finite difference formulae to the finite element solution. It is possible to using PLTMG, to introduce an error function based on the finite difference residual discussed in Sec 4.1.1. Using this error function PLTMG could adapt its mesh to minimize this residual[21], in a more sophisticated fashion. 89 Bibliography [1] R. Arnowitt, S. Deser, and C.W. Misner. Gravitation: An Introduction to Current Research, chapter ”The dynamics of general relativity”, page 227265. Wiley, 1962. [2] R. E. Bank. PLTMG: A Software Package for Solving Elliptic Partial Differ- ential Equations Users’ Guide 8.0. Philadelphia: Society for Industrial and Applied Mathematics, 1998. [3] David Bernstein, David Hobill, Edward Seidel, and Larry Smarr. Initial data for the black hole plus brill wave spacetime. Phys. Rev. D, 50(6):3760–3782, Sep 1994. doi: 10.1103/PhysRevD.50.3760. [4] Dietrich Braess. Finite Elements. Theory, Fast Solvers, and Applications in Solid Mechanics. Cambridge University Press, Cambridge, 1997. [5] Steven Brandt and Bernd Brügmann. A simple construction of initial data for multiple black holes. Phys. Rev. Lett., 78(19):3606–3609, May 1997. doi: 10.1103/PhysRevLett.78.3606. [6] Zhangxin Chen. Finite Element Methods and their Applications. Springer, Berlin, Germany, 2005. 90 [7] M. W. Choptuik. Consistency of finite-difference solutions of einstein’s equa- tions. Phys. Rev. D, 44(10):3124–3135, Nov 1991. doi: 10.1103/Phys- RevD.44.3124. [8] Y. Choquet-Bruhat and R. Geroch. Global aspects of the cauchy problem in general relativity. Commun. Math. Phys., 14:329–335, 1969. [9] Y. Choquet-Bruhat and J.W. York. The cauchy problem. In A. Held, editor, General Relativity and Gravitation: One Hundred Years After the Birth of Albert Einstein, volume 1, pages 99–172. Plenum, New York, U.S.A., 1980. [10] Greg Cook. Initial data for numerical relativity. Living Reviews in Relativity, 3(5), 2000. URL http://www.livingreviews.org/lrr-2000-5. [11] Gregory B. Cook. Corotating and irrotational binary black holes in quasi- circular orbits. Phys. Rev. D, 65(8):084003, Mar 2002. doi: 10.1103/Phys- RevD.65.084003. [12] Gregory B. Cook and Harald P. Pfeiffer. Excision boundary conditions for black-hole initial data. Physical Review D (Particles, Fields, Gravitation, and Cosmology), 70(10):104016, 2004. doi: 10.1103/PhysRevD.70.104016. URL http://link.aps.org/abstract/PRD/v70/e104016. [13] Gregory B. Cook, Matthew W. Choptuik, Mark R. Dubal, Scott Klasky, Richard A. Matzner, and Samuel R. Oliveira. Three-dimensional initial data for the collision of two black holes. Phys. Rev. D, 47(4):1471–1490, Feb 1993. doi: 10.1103/PhysRevD.47.1471. [14] L. de Figueiredo. Adaptive sampling of parametric curves, 1995. URL citeseer.ist.psu.edu/article/defigueiredo93adaptive.html. 91 [15] Robert Geroch and James B. Hartle. Distorted black holes. J. Math. Phys., 23 (4):680–692, 1982. URL http://link.aip.org/link/?JMP/23/680/1. [16] S.W. Hawking and G.F.R. Ellis. The Large Scale Structure of Space-Time. Cambridge Monographs on Mathematical Physics. Cambridge University Press, Cambridge, U.K., 1973. [17] S. Husa. Initial data for general relativity containing a marginally outer trapped torus. Phys. Rev. D, 54(12):7311–7321, Dec 1996. doi: 10.1103/Phys- RevD.54.7311. [18] Pong Soo Jang and Robert M. Wald. The positive energy conjecture and the cosmic censor hypothesis. Journal of Mathematical Physics, 18(1):41–44, 1977. doi: 10.1063/1.523134. URL http://link.aip.org/link/?JMP/18/41/1. [19] A. Lichnerowicz. L’integration des équations de la gravitation relativiste et le problème des n corps. J. Math. Pures Appl., 23:37–63, 1944. [20] Charles Misner, Kip Thorne, and John Wheeler. Gravitation. W. H. Freeman, 1973. [21] Jeffrey S. Ovall. Asymptotically exact functional error estimators based on superconvergent gradient recovery. Numerische Mathematik, 102(3):543–558. [22] Roger Penrose. Gravitational collapse and space-time singularities. Phys. Rev. Lett., 14(3):57–59, Jan 1965. doi: 10.1103/PhysRevLett.14.57. [23] Lszl B. Szabados. Quasi-local energy-momentum and angular momentum in gr: A review article. Living Reviews in Relativity, 7(4), 2004. URL http://www.livingreviews.org/lrr-2004-4. 92 [24] Saul A. Teukolsky. Rotating black holes: Separable wave equations for gravita- tional and electromagnetic perturbations. Phys. Rev. Lett., 29(16):1114–1118, Oct 1972. doi: 10.1103/PhysRevLett.29.1114. [25] Robert Wald. General Relativity. The University of Chicago Press, 1984. [26] Jeffrey Winicour. Characteristic evolution and matching. Living Reviews in Relativity, 8(10), 2005. URL http://www.livingreviews.org/lrr-2005-10. [27] J. W. York, Jr. and T. Piran. The Initial Value Problem and Beyond. In Spacetime and Geometry, pages 147–+, 1982. [28] James W. York. Gravitational degrees of freedom and the initial-value prob- lem. Phys. Rev. Lett., 26(26):1656–1658, Jun 1971. doi: 10.1103/Phys- RevLett.26.1656. [29] J.W. York Jr. Covariant decompositions of symmetric tensors in the theory of gravitation. Ann. Inst. Henri Poincare A, 21(4):319–332, 1974. [30] J.W. York Jr. Conformal “thin-sandwich” data for the initial-value problem of general relativity. Phys. Rev. Lett., 82:1350–1353, 1999. 93"@en ; edm:hasType "Thesis/Dissertation"@en ; vivo:dateIssued "2009-05"@en ; edm:isShownAt "10.14288/1.0066863"@en ; dcterms:language "eng"@en ; ns0:degreeDiscipline "Physics"@en ; edm:provider "Vancouver : University of British Columbia Library"@en ; dcterms:publisher "University of British Columbia"@en ; dcterms:rights "Attribution-NonCommercial-NoDerivatives 4.0 International"@en ; ns0:rightsURI "http://creativecommons.org/licenses/by-nc-nd/4.0/"@en ; ns0:scholarLevel "Graduate"@en ; dcterms:title "Initial data for axially symmetric black holes with distorted apparent horizons"@en ; dcterms:type "Text"@en ; ns0:identifierURI "http://hdl.handle.net/2429/3325"@en .