Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Axisymmetric boson stars in the conformally flat approximation Rousseau, Bruno 2003

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

Item Metadata

Download

Media
831-ubc_2004-0105.pdf [ 4.16MB ]
Metadata
JSON: 831-1.0091568.json
JSON-LD: 831-1.0091568-ld.json
RDF/XML (Pretty): 831-1.0091568-rdf.xml
RDF/JSON: 831-1.0091568-rdf.json
Turtle: 831-1.0091568-turtle.txt
N-Triples: 831-1.0091568-rdf-ntriples.txt
Original Record: 831-1.0091568-source.json
Full Text
831-1.0091568-fulltext.txt
Citation
831-1.0091568.ris

Full Text

Axisymmetric Boson Stars in the Conformally Flat Approximation by Bruno Rousseau B . S c , Universite de Montreal, 2001 A T H E S I S S U B M I T T E D IN P A R T I A L F U L F I L M E N T O F T H E REQUIREMENTS F O R T H E D E G R E E OF M A S T E R OF SCIENCE in The Faculty of Graduate Studies (Department of Physics and Astronomy)  We accept this thesis as conforming to the required standard  T H E UNIVERSITY OF BRITISH C O L U M B I A November 12, 2003 © Bruno Rousseau, 2003  In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission.  Department of Physics and Astronomy The University of British Columbia Vancouver, Canada Date ^  D<?c€^ber  ~1<^~B>  Abstract  ii  A b s t r a c t  In this thesis, boson stars are studied under the conformally flat approximation for the metric. Numerical simulations are used to extract solutions. The geometry is solved for using the constraint equations of general relativity and a generalization of the Klein-Gordon equation is used to evolve the matter fields. Type I critical phenomena is studied i n the context of this approximation in axisymmetry by perturbing a boson star with a real massless scalar field.  Contents  iii  C o n t e n t s  Abstract  ii  Contents  iii  List of Tables  vi  List of Figures  vii  Acknowledgements  xii  1 Introduction 2  Analytical Work  2.1 Conventions and Formalism 2.2 Constraint Equations 2.2.1 The Extrinsic Curvature 2.2.2 Gauss Weingarten equations 2.2.3 Gauss Codazzi equations 2.2.4 The Equations of Constraint 2.3 Dynamical Equations for the Geometry 2.4 Matter Content 2.4.1 Lagrangian and Hamiltonian Formulation 2.4.2 Stress Energy Tensor and Matter Distribution 2.4.3 Noether Charge 2.5 Approximation and Assumptions 2.6 Equations 2.6.1 Conformal Flatness 2.6.2 The Extrinsic Curvature 2.6.3 The Lapse function 2.6.4 The Ricci Scalar 2.6.5 The Explicit Equations 2.6.6 Boundary Conditions 2.6.7 Spherical limit 2.7 Initial Data 2.7.1 Single Static Spherical Boson Star 2.7.2 Transformation to Isotropic Coordinates 2.8 Apparent horizons  1 3  3 4 5 -. 7 7 8 9 12 12 15 15 17 18 18 19 21 21 22 24 26 • 27 27 29 30  Contents 2.9  iv  A D M Mass  34  3 Numerical Implementation 3.1 3.2  3.3  3.4  3.5  4  44 44 44 46 47 49  Multigrid Algorithm 4.1.1 Necessary Condition for Relaxation Time Evolution System 4.2.1 Amplification Matrix and Von Neumann Condition . . . . 4.2.2 Limiting cases  Numerical Testing  5-2  5.3 5.4 5.5  7  36 36 36 37 38 39 39 40 40 40 41 42  Scheme Analysis  5.1  6  - • •  4.1 4.2  5  36  Initial Data Definitions 3.2.1 Two Dimensional Domain Discretization 3.2.2 Discretization and Finite Difference Operators Solving the Elliptic System: Multigrid Methods 3.3.1 Discretization on the Grid 3.3.2 Relaxation 3.3.3 Transfers and the F A S Algorithm Solving the Hyperbolic System: Time Evolution 3.4.1 Interior Difference Scheme 3.4.2 Boundary Conditions Obtaining Solutions  Solution Error and Tolerance 5.1.1 Definitions 5.1.2 Solution Error Estimates 5.1.3 Values of the Tolerance Convergence and Independent Residuals 5.2.1 Convergence 5.2.2 Independent Residuals Stability of Boson Stars 5.3.1 Time Evolution of Boson Stars Spherically Perturbed Boson Stars Newtonian boson stars collision  52 '.  52 52 53 53 53 54 57 61 63 66 69  Critical Phenomena  71  6.1 6.2 6.3  71 72 74  Introduction Numerical Setup Results .  Conclusion  .  82  Bibliography  83  A Perturbed Boson Stars Data  85  Contents  B Grid Function Slicing for Critical Phenomena  v 102  List of Tables  L i s t  5.1  o f  vi  T a b l e s  Approximation to the truncation error for all grid functions at level 6,7. The value is taken at T=5.0, when the effects due to the use of the same initial data have disappeared. The tolerance for level 6 will thus be set to 2.0 x 10~ , which is one order of magnitude smaller than the smallest error above 64 Comparison between the frequencies obtained from the code and the values obtained by Hawley and Choptuik by perturbation theory. o is obtained by taking the square root of the values quoted in the article and multiplying by ao as obtained from the initial data code. The uncertainties are approximate. Values quoted come from level 7 data. The outer radii used for the domain are: R x — 64 for <f> = 0.06,0.1 and R = 48 for (p = 0.14 and <p = 0.18. The lack of resolution, due to numerical cost constraints, render the code unable to produce valid solutions at cf>o = 0.22, which is too close to the unstable branch 67 6  5.2  p e r t  ma  0  0  0  m a x  List of Figures  L i s t  2.1  o f  vii  F i g u r e s  Schematic of a trapped surface. S is the trapped surface, s° is the spacelike normal to S, n is the timelike normal to £ and Z° is the outgoing null vector from S a  3.1 3.2  5.1 5.2 5.3 5.4  30  Domain in the pz plane Coarse gridpoint (circle) surrounded by its fine grid neighbors (squares and triangles). The full-weighted injection operator attributes a value to a grid function at the coarse point by summing its values over the fine points with weights \ for the fine point at the position of the coarse point, | for the triangle points and ^ for the square points. The half-weighted injection operator attributes half the fine grid value of the grid function at the circle point as the coarse grid value of the grid function at the circle point. Note that this is equivalent to the full-weighted injection operator in the limit that the values are zero at the triangle points and equal at the circle and square points. The bilinear interpolation operator, which allows to transfer information from a coarser grid to a finer grid, distributes the coarse value at the circle point on the finer grid with following weights: 1 for the circle point, \ for the triangle points and \ for the square points. The actual value of a grid function thus transmitted at a given fine point is  37  the sum of all the coarse point values distributed on that point. .  43  The central value of the lapse as a function of time for the convergence testing trial data. The geometry is significantly curved. Convergence factor for trial initial data at level 5,6,7 for all the grid functions Convergence factor for trial initial data at level 6,7,8 for all the grid functions The convergence factor of da/dp is shown here. As can be seen, the convergence is much better than the convergence of a itself. This is possibly due to the fact that the boundary conditions on the lapse are only correct asymptotically, in the limit that the outer boundary of the computational domain tends to infinity. .  55 56 57  58  viii  List of Figures  5.5  Convergence factor calculated at level 5,6,7 for all the independent residuals for the trial initial data. Some of the 1/p terms are badly resolved on the axis at level 5; this is probably why such bad convergence is observed in certain grid functions 5.6 Convergence factor calculated at level 6,7,8 for all the independent residuals for the trial initial data 5.7 £2 norm of the independent residuals for a and ip. Clearly, these residuals are converging to zero. The other grid functions behave similarly 5.8 Universal quantities as functions of (p , as obtained from the initial data code. The mass parameter attains a maximum of 0.633 at the value cp n ~ 0.27. The radius of the star is defined as the radius at which the amplitude of the scalar field is 1% of its central value. In spherical symmetry, stars with cp smaller than cp n are stable and stars with <p bigger than the critical parameter are unstable 5.9 Norm of the scalar field, Noether charge and A D M mass as functions of coordinate time for <po = 0.1, Rmax — 32. A s can be seen in the Figure, results appear to deteriorate around t = 300, where the scalar field attains about 1% of the central value on the outer boundaries at the coarsest level 5.10 Norm of the scalar field, Noether charge and A D M mass as functions of coordinate time for (po = 0.1, R = 64. The figures show that the star is indeed stable and that the outer boundaries need to be placed further away than R =32 5.11 Amplitude of the scalar field for a truncation error perturbed boson star with parameter r/> = 0.18 using the spherically symmetric code briefly described in the text. Computations at levels 7 and 10 and shown. Here, R x = 32 which is the same value that leads to an unstable evolution at level 7 with the 2-D code. As can be seen on the plot, the level 7 spherical calculation also appears unstable whereas the level 10 data shows that the star is in fact stable. . 5.12 Newtonian collision between two light boson stars. The value of (po is 0.02 and the stars are initially 120 Planck units apart. The mass of each star is 0.28 Planck units. The central value of | $ | reaches a maximum at t ~ 2310 Planck units  59 60  61  0  CI  0  CT  0  62  65  max  max  66  0  ma  6.1  6.2  Typical form of the initial data. The norm of the complex scalar field and the real massless scalar field are superimposed in the figure Central value of a for diverse values of the family parameter for supercritical level 5 data. As can be seen in the figure, the closer the parameter is to its critical value, the longer the solution behaves like the critical one  68  70  73  76  List of Figures 6.3 6.4  6.5  6.6  6.7  ix  Central value of a for diverse values of the family parameter for supercritical level 6 data Lifetime spent in the critical regime as function of the family parameter, A, for level 5 data. Proper time at the center of the grid is used; the time variable must be rescaled using r = J adt (the shift is negligible). The value of the critical parameter, denned in equation (6.4), is given by 7 = 9.7, where the quoted uncertainty comes from the uncertainty on A*. The value of 7 is obtained by making a linear fit to the data Lifetime spent in the critical regime as function of the family parameter, A, for level 6 data. Proper time at the center of the grid is used; the time variable must be changed using r = / adt (the shift is negligible). The value of the critical parameter, defined in equation (6.4), is given by 7 — 3.7, where the quoted uncertainty comes from the uncertainty on A*. The value of 7 is obtained by making a linear fit to the data Lifetime spent in the critical regime over comparable ranges of the family parameter, A, for level 5 data and level 6 data using coordinate time. The values of 7 are obtained by making linear fits to the data. The disagreement between level 5 and level 6 is much smaller here than for the computations of 7 that were made using central proper time. A D M mass for supercritical and subcritical solutions for level 5 data. Once an apparent horizon forms, the value of the A D M mass rapidly increases, probably indicating that the scheme quickly becomes less accurate due to the steep gradients arising from the impending coordinate pathology. When the solution oscillates in the critical regime, though, the bulk of the matter is concentrated around the origin and the computed mass yields a sensible result. Clearly, the solution does not "shed mass" while in the critical regime (type II behavior) and the black holes formed must have a finite mass, which is typical of type I behavior  Amplitude of the scalar field as a function of time. After t ~ 150, boundary effects generate high frequency noise. Level 8 data was only generated up to t=250 due to high computational costs. Parameters: (f> = 0.06, A = 0.02, R =64 A.2 Convergence factor for some of the grid functions at level 6,7,8 for 50 < t < 150. Parameters: <f> = 0.06, A = 0.02, R =64 A.3 Convergence factor for some of the grid functions at level 6,7,8 for 50 < t < 150. Parameters: <f> = 0.06, A = 0.02, R =64 A.4 Convergence factor for some of the grid functions at level 6,7,8 for 50 < t < 150. Parameters: 4> = 0.06, A = 0.02, R =64  77  78  79  80  81  A.l  0  0  max  86  0  max  87  0  max  88  max  89  List of Figures  x  A.5 A D M mass at level 6,7,8. Boundary effects become important around t ~ 800 and the solution is probably untrustworthy after that time. The mass is only computed up to t = 250 at level 8 because of large computational cost. Parameters: = 0.06, A = 0.02, R =64 A.6 Noether charge at level 6,7,8. Boundary effects become important around t ~ 800 and the solution is probably untrustworthy after that time. The charge is only computed up to t = 250 at level 8 because of large computational cost. Parameters: 4> = 0.06, A = 0.02, i ? =64 A.7 Amplitude of the scalar field as a function of time. Level 8 data was only generated up to t=150 due to high computational costs. Parameters: 0 = 0.10, A = 0.02, R =64 A.8 Convergence factor for some of the grid functions at level 6,7,8 for 50 < t < 150. Parameters: <f> = 0.10, A = 0.02, R =64 A.9 Convergence factor for some of the grid functions at level 6,7,8 for 50 < t < 150. Parameters: = 0.10, A = 0.02, R =64 A . 10 Convergence factor for some of the grid functions at level 6,7,8 for 50 < t < 150. Parameters: (f> = 0.10, A = 0.02, R =64 A . l l A D M mass at level 6,7,8. Parameters: 4> = 0.10, A = 0.02, R =64 A.12 Noether charge at level 6,7,8. Parameters: 0 = 0.10, A = 0.02, i ? =64 A.13 Amplitude of the scalar field as a function of time. Only level 7 data is presented. A t lower resolution, the star is unstable and no oscillating solution could be produced. Parameters: </> = 0.14, A = 0.02, i ? =32 A . 14 A D M mass and Noether charge at level 7. Boundary effects render the solution invalid around t=800. Parameters: = 0.14, A = 0.02, R =32 A . 15 Amplitude of the scalar field as a function of time at level 7. The star is perturbed by truncation error. Parameters: <j> = 0.18, A = 0.0, = 32 A.16 (f> = 0.18: A D M mass and Noether charge at level 7 for truncation error perturbation. Parameters: <p =. 0.18, A = 0.0, i?m = 3 2 max  0  m a x  O  max  0  max  max  0  max  0  max  O  0  90  91  92  93  94  95 96  m a x  97  m a x  98  max  0  99  100  0  0  ax  101  List of Figures B.l  xi  p = 0 (dotted line) and z = 0 (full line) spacelike slices of the lapse for \A - A*\ ~ 1 0 ~ at level 5. A t initial times, the lapse is highly aspherical because of the perturbation, but once the solution is in the critical regime (t > 400) the radial and axial slices become close to identical. The critical solution's axial slice past t ~ 500 seems off centered by one grid point to the left with respect to the radial slice. This effect converges away, as will be seen at level 6 15  B.2  p = 0 (dotted line) and z = 0 (full line) spacelike slices of the conformal factor for \A - A*\ ~ 1 0 at level 5. The same features present i n the lapse data can be observed. 104 p = 0 (dotted line) and z = 0 (full line) spacelike slices of the amplitude of the complex scalar field for \A — A*\ ~ 1 0 ~ at level 5. Even though the boson star parameter is 4>o = 0.27, the presence of the perturbation changes the amplitude of the boson star significantly. The star appears to be i n a aspherical excited mode, as the nodes on the axis slice indicate. The star seems to be dissipated away before the end of the lifetime of the critical solution, but a look at the boson star data at a relative scale by reducing the range of the vertical axes indicates that massive scalar field is still present in the vicinity of the origin for the duration of the critical regime. Clearly, the star is not well resolved; it occupies roughly a subgrid of 7 by 15 points 105 p = 0 (dotted line) and z =.0 (full line) spacelike slices of the amplitude of the real scalar field for \A — A*\ ~ I O at level 5. The data presented is super-critial, i. e. it leads to the formation of a black hole 106 p = 0 (dotted line) and z = 0 (full line) spacelike slices of the lapse for \A — A*\ ~ 10~ at level 6. The asymmetry observed at level 5 is now gone 107 p = 0 (dotted line) and z = 0 (full line) spacelike slices of the conformal factor for \A - A* | ~ 10~ at level 6 108 p = 0 (dotted line) and z = 0 (full line) spacelike slices of the amplitude of the complex scalar field at a late time for | A — A* | ~ 1 0 at level 6. The star is better resolved than at level 5, but is still clearly very poorly resolved 109 - 1 5  B.3  15  B.4  - 1 5  B.5  9  B.6  9  B.7  103  - 9  Acknowledgements  xii  A c k n o w l e d g e m e n t s I would like to thank all the members of the Numerical Relativity group at U B C , especially my supervisor Matthew Choptuik and Inaki Olabarrieta, for their help and their wise counsel. I would also like to thank all my fellow graduate students at St. John's College-UBC for their support and friendship.  Chapter  C h a p t e r  1.  Introduction  1  1  I n t r o d u c t i o n  The gravitational interactions between heavenly bodies is described by Einstein's masterpiece theory, general relativity. This theory, although conceptually simple and very elegant, is difficult to apply to generic situations. Indeed, one cannot escape its inherent non-linearities i n the strong gravitational regime. Thus, to study situations of real astrophysical interest, one is naturally led to numerical relativity, the study of Einstein theory through numerical simulations. Especially interesting to the general relativity community is the question of how infalling neutron star binaries behave before they coalesce. These systems are believed to be good gravitational radiation source candidates and a quantitative understanding of their interactions would go a long way towards interpreting future measured data from the big gravitational wave detectors. In an effort to reduce the computational work required to solve the equations of general relativity, Wilson, Matthews and Marronetti [7] were led to propose the Conformally Flat Condition ( C F C ) approximation, which states that when gravitational radiation is small , the spatial metric can be approximated to a good precision by a flat space metric multiplied by a conformal factor. If correct, this approximation can potentially yield approximate solutions at a fraction of the numerical cost of solving the full Einstein equations. Another problem of interest i n general relativity is the question of critical phenomena. These phenomena, as first observed by Choptuik and described in Gundlach's review article [13], describe the behavior of spacetime at the threshold of black hole formation. Basically, by considering a one-parameter family of spacetimes (family parameter p), where large values of p lead to black hole formation and small values lead to matter dispersal, one can "tune" p to p*, the critical value above which black holes are formed. The very high cost of tuning the parameter down to machine precision has confined numerical studies of this phenomena to spherical symmetry: very few investigations have been performed in more than one spatial dimension. This lack of higher dimensional investigation combined with the proposed C F C approximation has inspired this work. In this thesis, spherical boson stars perturbed by an axisymmetric massless scalar field pulse are studied in the C F C approximation and the main questions asked are whether the approximation will allow for critical phenomena and if the obtained behavior agrees with the spherical solutions obtained by Hawley and Choptuik [6]. This problem is technically important because if critical phenomena can be accurately studied within the C F C approximation, it will be possible to explore critical phenomena away from spherical symmetry at a fraction of the numerical cost of the full relativistic equations.  Chapter 1. Introduction  2  In the second chapter of this thesis, most of the analytical work leading to the expressions of the equations within the chosen approximation is presented, following standard references. Furthermore, the equations relevant to the generation of initial data and expressions for conserved quantities are obtained. In the third chapter, the numerical implementation and solution process is described; in Chapter 4 some numerical stability analysis is performed and a Von Neumann stability condition is obtained. In Chapter 5, the numerical solutions are tested for convergence and conservation of mass and Noether charge. Furthermore, the perturbation frequencies of stable boson stars are obtained from the axisymmetric numerical implementation as a non-trivial test of the correctness of the code. Finally, critical phenomena are studied in Chapter 6. It was found that the C F C indeed allows for critical phenomena, but that the obtained scaling exponents are ambiguous; the result obtained at a given resolution is found to agree very well with the spherical result of Hawley and Choptuik, but the same exponent computed at a finer resolution is found to be quite different. Further work would be necessary to clarify this discrepancy.  Chapter  C h a p t e r  2. Analytical  Work  3  2  Analytical  W o r k  Unless otherwise stated, most of the derivations in this section will follow the standard derivations of Misner, Thorne & Wheeler [1] and York [3]. Throughout this thesis, the "god given" units will be used: G = c = h=l  (2.1)  This, in effect, means that all quantities will be measured in units of the Planck quantities, the Planck mass, Planck time and Planck length:  2.1  Conventions and Formalism  Spacetime is considered as a 4 dimensional manifold endowed with a Riemannian geometry. The A D M 3+1 formalism will be used to describe the split of spacetime into a foliation of spacelike hypersurfaces, { E } , identified by a timelike coordinate. Abstract tensorial notation will be used when dealing with intermediate equations, i n which Latin letters ranging from a to f used as indices will indicate ranks of tensors but will not represent actual components. When dealing with equations from which one expects to compute numerical quantities, tensorial components will be used, for which Greek indices will range from 0 to 3 (spacetime indices) and Latin indices from i to n will range from 1 to 3 (space indices). Two vector bases will be used extensively; (n ,{ef}) and ( £ , { e " } ) , where {ef} is the set of spacelike vectors consistent with the given form of the metric. The covariant derivative in the 4 dimensional manifold will be denoted V whereas it will be denoted D on a given spacelike slice E . V i = V ? will denote the covariant derivative in the direction of the vector e\ and = V « represents the covariant derivative in the direction of the timelike unit normal to E . The space covariant derivative is related to the spacetime covariant derivative by projection: for example, for a vector A lying on £ , D A = (S +n n )(S + n n )V A . 4 Christoffel symbols will be denoted rjj;„ when they relate to the 4 dimensional covariant derivative and Tj when the relate to the 3 dimensional one. Further, R b and K b will be used to denote the Ricci tensor and the extrinsic curvature on a spacelike slice respectively. a  a  a  a  e  n  b  b  a  b  d  d  c  k  a  a  c  c  a  a  b  d  Chapter  2. Analytical  Work  4  The study of general relativity is dominated by the Einstein equation: Gab = 8nT b a  The Einstein tensor is given by:  1  (2.3)  Gab — %-ab ~ ^9ab^-  where % b is the Ricci tensor and "R is the associated Ricci scalar in the 4 dimensional manifold. In the A D M formalism, the metric has the form: a  -{a -p p ) 2  • Pi  k  k  or 2  oc . . ~2 2  I  a  The geometry of the spacetime is thus expressed through the following variables: a is called the lapse and measures the time interval between two spacelike hypersurfaces, for observers moving normally to the slices; P is the shift vector and measures how the space coordinates move about from one slice to another; finally, 7„t, is the 3-metric of the spacelike hypersurface. The unit normal to £ is denoted n . In terms of this normal, the tangent vector to the spacetime curve of constant space coordinate, i.e. the vector pointing in the coordinate time direction, is given by t" = an + p . In terms of components in the A D M coordinates, n = (l/a;—P /a) and n = ( - a ; 0 ) . Further, ^ b = g b + n rib as can be shown by projection along all the basis vectors. The stress-energy tensor can also be split into components lying in E and normal to it. The form of the tensor, in all generality, is given by: a  a  a  M  a  l  M  a  T b = S b + J( nb) + P a  a  a  a  a  (2.4)  Here S b is the spacelike stress tensor, J is the matter current density and p is the matter energy density, all as measured by observers moving normally to the slices. a  2.2  a  m  Constraint Equations  Not all components of the Einstein equation are related to the dynamics of the geometry. In fact, 4 components of the equations only involve data on a given spacelike slice E ; these equations are called constraint equations and they are related to the initial value problem. These equations must be initially satisfied before the dynamics of the Einstein equation can be studied. Differential geometry shows how these components of the Einstein equation are amenable to a formulation involving only the geometry of a given slice. First, the extrinsic curvature must be denned.  Chapter 2. Analytical Work  2.2.1  5  The Extrinsic Curvature  n being the unit normal to the slice E , its covariant derivative along a spacelike direction must also lie in E in order to preserve its length. Hence, the components of the extrinsic curvature tensor are defined by: a  V>  =-K/e°  a  •  (2.5)  The extrinsic curvature tensor is understood to be lying entirely i n E and can be shown to be symmetric:  =  =  -{V )e)g  =  -Vi{n e))  =  n ( r £ > ) = n ( r£e>)  =  Kn  a  in  ab  -{V i )e  b  iT b  +n Vie)  b  b  o 4  4  6  e  b  (2.6)  Here the symmetry of the Christoffel symbols and Vj<? (, = 0 were used. The extrinsic curvature also has the the following useful properties. Consider the Lie derivative of the space metric along the unit normal vector field: 0  £filab  where a = V n . a  n  =  £n(9ab + n Tl )  =  £h(g b)  =  Vn  a  a  a  b  + £h{n )n a  +  b  + V n +n a  b  b  a  a  n £ (n ) a  n  b  + na  b  b  (2.7)  a  Clearly:  a  nVn  =  a  b  a  V (n n )  -  -n V n  =  a  b  .=  a  = na  a  a  b  a  a  a  b  -n V n a  b  a  0  (2.8)  =  a  nVn  n n V n =n {n \/ n ) a  b  b  b  a  a  b  a  = 0  "  (2.9)  Projecting the Lie derivative along various directions:  n n £ "f a  b  n  ab  n e £nlab a  b  i  =  n n (V n  =  na  =  0  a  + V6na + na  b  a  b  b  a  + na a  b  — na b  a  + na)  b  b  a  — na a  b  a  (2.10)  =  n e (V n  =  e\a - e a  =  0  a  + Vn  b  a  b  b  a  + na a  b  + na) b  a  b  b  b  (2.11)  Chapter e"eji'ft7a!,  2. Analytical  Work  6  =  e"eJ(V ni, + V ; , n + n a\, + n a )  =  efrim  =  -2K  a  a  a  b  a  + elVjUa (2.12)  tj  One thus obtains the geometrical, coordinate free result: Kab = -^£hlab  (2.13)  The extrinsic curvature can also be expressed v i a projection. Consider -7a 7(, V n C  d  (c  d)  =  -^(S +  n n )(S +n n )(V n  =  ~]^{^anb + ^bn +n n V n  +  +n n V n  2n nbn n V n )  c  c  a  d  a  + Wn)  d  b  b  c  d  d  d  a  c  a  d  + n ra Vf,n +  c  a  b  c  b  =  -^(V n  =  -\£hlab  a  0  + V n  &  6  a  d  b  a  +n a  a  nnVn c  c  c  d  a  d  c  d  +n a )  b  b  a  (2.14)  where equation (2.7) was used. Thus, using equation (2.13): K  ab  = -7aVV  n  ( c  (2.15)  d )  This result can be taken further. Consider:  wb a  =  -7a 7ft V[ n  =  -\~ialb {Vcn -  C  d  c  d]  V n )  d  d  d  (2.16)  c  and project this equation along all the basis vectors: nnw a  =  b  ab  = n e\  a  Wab  e e\w a  ab  b c la  d  b  [c  (2.17) -n eh l V n a  c  a  d  b  [c  =  0  =  -e?e$7 VV 7i  d]  (2.18) [ c  Q  = =  d]  0  =  a  {  -n n j V n  d ]  -e ie W n C  d  [c  - -ele {\/ n l  d  c  d  1 {epin 2~ K  d  j  d]  - V n ) c  d  - e?V,-n ) c  {K -K ) 2' l  ji  0  ii  (2.19)  Chapter Hence w  ab  2. Analytical  7  Work  vanishes identically. This implies K  ab  =  -7o 76 V n  =  -V n(, - na  C  d  c  a  a  d  (2.20)  b  and that therefore V n a  2.2.2  = -K  f t  - na  ab  a  (2-21)  b  Gauss Weingarten equations  As it turns out, the spacetime covariant derivative of a space vector along a space direction can be expressed using the extrinsic curvature and the space covariant derivative.  V »= r?.e£  (2.22)  4  ie  Projecting this equation along all the basis vectors yields the desired result: .n„VieJ  =  V i K e " ) -e?Vi(n„) = -e?Vi(n„).  = Kij OabelVie" = =  (2.23)  5 „ e £ e ? V e ? = e e{ 7 S / V e ^ d  6  c  lbd  e e{D e)  =  d  lbd  k  f  k  0  c  e%Die)  lah  = ^ele^ = r k 7w  .(2.24)  Comparing the above results with the general form of V i e " V i e ? = C$e% + C n  a  2  (2.25)  the Gauss Weingarten equations are obtained: V i e " = T%e% - K n  a  i3  2.2.3  (2.26)  Gauss Codazzi equations  Some components of the 4 dimensional Riemann tensor can be expressed i n terms of the extrinsic curvature and the spatial Riemann tensor. The Gauss Codazzi equations describe how to extract those components. B y definition: X(e«,e )et b  j  = (V V -V V )et i  j  j  i  (2.27)  In the coordinate basis used, this implies = (V*V,- - V , - V i K  (2-28)  Chapter  2. Analytical  Work  8  Using the Gauss Weingarten equations and the definition of extrinsic curvature one finds ViV^  =  ,  Vi{T e1-K n ) l  a  ik  jk  =  d (T )ef  +  =  d (T )et  + T] (rfle -K n )-d (K )n  =  {$(1%) + Tf T'  l  i  jk  l  i  T V e^-d (K )n -K W n l  a  jk  a  jk  k  a  m  k  i  a  jk  jk  i  +  a  il  i  jk  + K K>}e?  im  i  - {di(K )  jk  jk  +  T K }n l  a  jk  u  (2.29) Permuting indices i and j in the last expression and subtracting, one gets the Gauss Codazzi equation:  =  {0i(i%)-0;(rU) + -{di(K )  - dj(K )  jk  =  {R  l kij  r ^ + T K«  + K K\  -  l  ik  jk  - KitK'^ef  jk  - {D K t  - / J ^ K  jk  (2.30)  From this, one can deduce 14 of the 20 independent components of the 4 dimensional Riemann curvature by projection: %iki  =  9a e^ el  —  Rijki + KjiK  =  DjKik-DkKij  b  jkl  — Kj K  ik  k  (2.31)  it  (2.32)  The coordinate independent Gauss Codazzi equations are thus: -L^afccd  =  Rabcd + K -K  -L^na6c  =  DbKac ~ D K  ba  c  ac  ~ K Kd bc  a  (2.33) (2-34)  ab  where " ± " means that all free indices are projected onto S.  2.2.4  The Equations of Constraint  The constraint equations are obtained by projecting the Einstein equation. The Hamiltonian constraint is obtained by: nnG a  b  ab  n n (Jl -^g X) a  8nn n T a  b  ab  ab  =  Sirp  Xhn + \%  =  87rp  b  ab  =  m  ro  (2.35)  Here, equations (2.3) and (2.4) have been used. The left-hand side is easier to compute with indices in the ( n ; {e?}) basis. Accounting for the symmetries of a  Chapter 2. Analytical Work  9  the Riemann tensor and using the Gauss Codazzi equations: 2 1 2  =  liR + K^K^-KijK^) 2  (2.36)  The momentum constraint is given by: n e\G  =  a  ab  1  ^n e)T a  ab  n e*(:R , - ^g X)  =  -8TT J*  Xni  =  -8nJi  a  ai  ab  (2.37)  Using the second set of Gauss-Codazzi equations, one obtains: •^m  T'' '•^•jhki  =  =  -^"(DkKij-DiKjk)  (2.38) i  Thus, the constraint equations are given by: R DK  ah  b  2.3  =  -K\K\  + K K  =  j DK  + &TrJ  ac  b  c  + 167rp  ah  ab  (2.39)  m  (2.40)  a  b  Dynamical Equations for the Geometry  The remaining components of the Riemann tensor can be related to the dynamical evolution of the geometry. In particular, it will be useful to relate some of these components to the Lie derivative of the extrinsic curvature along the coordinate time direction. First, the Lie derivative of the extrinsic curvature along the unit normal field is given by: £K n  = nVK  + K Vn  c  ab  c  + K Vn  c  ab  cb  (2.41)  c  a  ac  b  Next, using equation (2.21), one obtains:  V V n a  ft  c  = V (-K a  bc  -n a ) b  c  =  -V K  - V (n )a  =  -V Kbc  + K a  a  bc  a  Hence, in geometrical language:  a  ab  b  c  c  - n W {a ) b  + naa a  b  c  a  c  - n V (a ) b  a  c  (2.42)  Chapter 2. Analytical Work  Xabc n  a  d  =  (V V - V V )n  =  VK  d  c  c  c  cd  -V K d  d  - Ka  bd  6  - naa  b  c  +K a  bc  cd  10  b  + nWa  d  d  +naa  b  d  b  c  b  - nVa  c  c  d  (2.43)  b  Projecting one index of the above along the normal direction, one finds Rabcd.n n a  d  =  nVK  - Va - nVK  =  -n V K  d  c  c  d  bc  bd  b  b  bc  + K S/ n  cd  d  c  d  +K Vn  b  cd  b  - Wa  d  =  -£ K  - KK  - V a - nnVa  b  c  b  bc  c  6  c  - naK b  d  b  b  - V a  d  bd  -n n V a  - aa  d  b  d  d  c  c  6  - aa  d  d  nnVa  c  - K Vn  c  d  cd  c  b  d  c  £K n  bc  - aa  =  n  d  - K Vn  d  d  b  d  nnVa c  - aa -  d  bd  d  b  b  (2.44)  c  Projecting the remaining indices onto E then yields —  j^bacdn rX  ^•ifijh  =  £K  + KikK) + a  =  £K  +K Kj  h  ij  n  %Va  d ilb  + Djdi + a aj  k  ij  + e\e  c  iaj  ik  t  a  a  d  (2.45)  This result can be further simplified i n an A D M coordinate system. Indeed, in such coordinates, = (—a;0). Thus V (^)  =-  v  4  r ^ ) = V (^)  • (2.46)  M  so one has  V„n = V n„ + ^(n V„a - n^V^a) M  M  M  (2.47)  which leads to  =>  fli  =  n' (V n + ^(n V„« - n V^a))  =  —V„a + - n „ n V„a  /  p  1/  1„ a  M  1 a  v  „,_,  = A(ln(a))  (2.48)  This further implies: Djdi + aiCtj  =  DjDi(\n(a)) + D ( l n ( a ) ) £ » j ( l n ( a ) )  =  —DiDjCt a  i  (2.49)  Chapter  2. Analytical  Work  11  Thus, the dynamical equation in A D M coordinates becomes: £ Kij  = y. j  h  ih  - KK  - ^-DiDja  k  h  ik  j  (2.50)  Next, we turn to the Einstein equation which relates the Ricci tensor to the matter stress-energy tensor. The Einstein equation can be rewritten: Xab = 8n{T - ^g g T )  (2.51)  cd  ab  ab  cd  Projecting onto S , using equation (2.4) and the Gauss Codazzi equations, we find g^ninj+^ikj -R  + Rij + K K  - KK  k  ifljh  i3  k  k  ik  j  =  MSio +  \liAPm-S\))  =  8Tv(Sij +  (p - S )) k  m  k  From this is obtained: £K fl  =  ij  - ^ D i D j a + Rij + K K  k  {j  - M ^  k  -  2K K  k  ik  j  +p;(Pm-.S\))  (2-52)  In order to turn this into a Lie derivative in the coordinate time direction, consider: £K t  ab  =  (an +n^cK  =  a£ K  + £pK  =  a£ K  + £pK  + K V (an +p )+Kcb^a(an  c  c  ab  h  h  ab  ac  + n°(K 'V a  ab  ab  c  c  b  ac  b  +  +  n  K V a) cb  a  (2.53)  ab  Similarly: £tlab  =  Ct£nJ  + £piab  ab  (2-54)  Both of the above results rely on the spacelike nature of the space metric and the extrinsic curvature, i. e. their vanishing product with the unit normal. Furthermore, the Lie derivative of the inverse space metric is given by: £n(l lcb) aC  = =  ^£nl  ab  'Ycb£fi1  = =  where equation (2.13) was used. derived:  0 ~~ t 7 £n~Ycb  -l l £nlcd aC M  2K  ab  (2.55)  Finally, the following useful result can be  Chapter 2. Analytical  £ (K\) t  =  a£ (K )  +  =  a£ (Y )K  =  2aK K  =  - DDa  i  fl  i  i  0  - Y'DtDja  ij  i  £ (K\) 0  + a{R + K^K^  - 2K K  i  j  m  i  S\))  m  +£ {K\)  i  j  - 47r(3p -  ii  ii  + a{R + K\K -4ir{?>p -S ))  ii  1  i  + aryVXfiKij +  ij  ii  12  £ (K )  j  n  Work  (2.56)  0  Using the Hamiltonian constraint, equation (2.39), the above result reduces to: £ (K\)  =  t  2.4  - DDa  + a(K K ^Aix{  ii  1  i  + S )) + £ (K )  i  i  id  i  Pm  i  i  0  i  (2.57)  M a t t e r Content  A massive complex scalar field is chosen as a matter source primarily because it is the simplest type of matter that allows for the existence of a static star-like configuration, i.e. a self gravitating lump of matter.  2.4.1  Lagrangian and Hamiltonian Formulation  The matter content is described by the scalar field: $ = 0 i + i(i>2  (2.58)  where <j>\ and 4>2 are real-valued. The Lagrangian density associated with this field is L  *  =  -^-(ff V $V $*+m $$*)  =  - ^ ( 3  =  i  o 6  0  1  a  + Z,  (2.59)  2  Q  t  V  6  0  ^ V ^ + m V ) - ^(  a 5  6  V 0 V^ a  2  2  +mW) (2.60)  02  where m is the mass parameter of the scalar field. We thus obtain the following action for the matter i n the curved spacetime: I  = jV=9(i0.  + Lfr)dV  (2.61)  Here, dV is the coordinate volume element and y/—g is the square root of the determinant of the metric. Hence, ^/^gdV is the proper volume element in the curved spacetime. Extremizing this action with respect to each component of the scalar field, we get the Klein-Gordon equation. However, from the 3+1 point of view, it is more useful to consider the Hamiltonian formulation of the  Chapter  2. Analytical  Work  13  dynamics of the scalar field. Explicitly, i n an A D M coordinate system, the Lagrangian is given by:  +  ~ ( - % + 2%/rc^ 0 07T \ cr or  - ^ d r f a d r f a + m  2  Y di<f> d i> j  a  ^  2  j(  a  (2.62)  )  where 4> is the coordinate time derivative of the a-th component of the scalar field (i.e., a takes on the values 1 or 2). We define the conjugate momentum field as the derivative of the effective Lagrangian (i.e. the expression that appears in the action) with respect to the coordinate time derivative of the scalar field: a  5(^/ gL ) ::  4>a  Ha  (2.63) This then implies: 4TTO:  2  <$>a = —7=O  a  +  (2.64)  P di(f> l  a  W i t h this i n hand, a Hamiltonian density can now be constructed:  =  ^(L  &i4>i + <72</>2  0  1  +L  0  2  )  (2.65)  We thus have  2-na 2  =  V~ 9  2  + audita  + V-J-^idifadrta °K  +^ - r n ^ OT!" 2  2  (2.66)  The Hamiltonian is then given by:  Di = J H$dv  (2.67)  where the integration is over a 3-volume. The associated dynamical equations are given by: <Pa =  J7ocr a  5-K Ha  The functional derivative obeys the following axioms:  (2.68)  Chapter  2. Analytical  Work  14  ^(y)=S(x-y)  (2.69)  J' F{y)cj>{y)dy = F(x) Scf>(x)  (2.70)  Further, the functional derivative is distributive and when it acts on the derivative of the field, partial integration must be used, noting that boundary terms vanish:  64>(x) These properties generalize to an arbitrary number of dimensions. Using these properties, we get, for example: S_ 54>  + =  I  § - Y ^ l ^ ) d  -2di(-  V  y  (2.72)  8TT  Using similar derivations, the dynamical equations for the matter can be deduced. In order to simplify the expressions, we define an alternate conjugate momentum to the matter field: a = 4iro a  (2.73)  a  In terms of these new fields, the dynamical equations are given by:  4>a O-a =  (2.74)  '-9 d (p O )+d (y/^ d 4 )-^m <f> i  i  i3  a  i  1  2  j  >a  a  (2.75)  These expressions reduce to the flat spacetime Klein-Gordon equation for the case of a flat geometry, as they should. The Lagrangian density can also be expressed i n terms of the new conjugate momentum: (2.76)  Chapter  2.4.2  2. Analytical  Work  15  Stress Energy Tensor and Matter Distribution  From the action principle for the geometry of spacetime, the stress energy tensor is given by: T  = - 2 ^ | + g L^  ab  (2.77)  ab  Noting that: SL^  _  n^d^a  =  1  (d^Ma)  -(t»  (2.78)  - p^d^  =  (2.79)  and using equation (2.4): the following components of stress energy tensor for the complex scalar field can be deduced,  i =  Sij  2  2  ^ T ^ i ^ f a + ^ ^ M a + r n  — Tij  2  ^ )  (2.80)  o=l  (2.82)  1  =  2  ^ E (  o=l  2.4.3  2  3  7 ^  v  C  3  T  a - 7  i  J  ^  a  d ^  a  - 3 m y ^  (2.83)  '  Noether Charge  The celebrated Noether's theorem allows one to compute conserved quantities from infinitesimal symmetries of the Lagrangian density. Following derivations found in Peskin & Schroeder [17] and in Wald [2], in curved space time one obtains Noether's theorem as follows. Consider a Lagrangian density L(4>, VM(/>)  Chapter 2. Analytical Work  16  and an infinitesimal transformation cf> -» <j> + 9A<f>—where 9 is an infinitesimal parameter and A</> is a given function—leaving this Lagrangian density unchanged. Then:  L(^,V 0) =  L((j> + 0A<f>,V  M  =  +6 A$))  L(<f>, V„0) + 9 A ^ + 0 V ( A 0 ) - | ^ — + 0(9 ) 2  M  (2.85)  This implies: o  =  A4£  + V „M V( A 0r  d  L  at • ' "* 'a(v 0) M  =  V  ^ (  ^ ^ )  A  +  A  K | -  V  s  ^ ^ ) ) ) v  =0!  '  A conserved current is thus obtained:  V j"  =  M  0  (2.88)  If more than one function A<f> is considered, all individual currents of the above form must be summed. In order to obtain a conserved charge, the concept of a tensor density is needed. First, define the current density: r  = V=93»  (2-89)  and the covariant derivative of the determinant:  V  ^ =5 v^- r^y=5 4  M  M  (2-90)  The rationale for this definition is that, even though the determinant of the metric is not a tensor and its covariant derivative makes no formal sense, the covariant derivative of the Levi-Civita tensor entering the definition of the element of proper volume reduces to the above expression. Using the above definitions, we obtain:  =  d^r  (2.9i)  From this, a conserved charge can be deduced: Q = J 3°d V 3  (2.92)  Chapter  2.  Analytical  Work  17  Indeed, using the vanishing of the current at spatial infinity, we have Q dt  d  _ = ~  '  f *0J3 jd 3°d "  =  - j  =  0  a  0  3  aa'dPv  (2.93)  It is clear from equation (2.59) that the transformation $ -> e' $ leaves the Lagrangian density, and thus the equations of motion, invariant. For 9 infinitesimal, the above transformation reduces to: 9  <j>! -> A0i  =  (pi- 9(p  (p2 -> 4>2 + 9 (pi  -4>2  Ah  2  = (pi  (2.94) (2.95)  Using equation (2.87), the form of the conserved current is deduced: 3  .— =  Ji +32  ^((h^drfa-Wdufa)  (2.96)  Thus, using the form of the A D M metric, the conserved charge is given by: 47T az 1  47T  a  z  =(0lCT -02CTl)  (2-97)  2  Q =  K  j{(PKJ2  - 4>2Tl)d V 3  (2.98)  Here K is an arbitrary constant and d V is the element of coordinate volume. 3  2.5  Approximation and Assumptions  The goal of this thesis is to investigate the dynamics of a massive complex scalar field in axisymmetry under the conformally flat condition ( C F C ) . This approximation is discussed by Wilson, Matthews & Marronetti [8]. A cylindrical coordinate system (t,p,z,(p) is used on the spacelike hypersurface and azimuthal symmetry is assumed. The C F C prescribes the form of the metric on the 3surface: / ip 0 0 4  Hi =  0  \ o  V o 4  o  Heuristically, this approximation is thought to minimize free radiation, i.e. the emission of gravitational disturbances that are not directly correlated with disturbances with in the matter fields. The validity of this approximation remains  Chapter 2. Analytical Work  18  to be tested in order to confirm whether or not it is of any relevance to situations of astronomical interest. The C F C greatly reduces the number of functions that need to be computed. Moreover, the full Einstein equations are not necessary to determine those functions that remain in the evolution scheme: indeed the geometric quantities can be computed from the constraints, as well as from the specific choice of time coordinate: only the matter fields need to be updated using equations of evolution type. To fix the time coordinate, the maximum slicing condition will be used. This gauge condition demands that the trace of the extrinsic curvature vanishes at all times, i.e. = £K = 0; its use generates an elliptic equation for the lapse. Hence, using equations (2.39), (2.40) and (2.57), the system describing the geometry is given by: i  t  R  =  DiK  =  ij  KylC*+  i  16irp  m  8irJ  j  Because of the symmetry, = 0, and we are left with the following set of functions that completely characterize the geometry and the matter field: a  =  ct{t, p, z)  1> =  i>(t,p,z)  P>  n  2.6  =  P (t,p,z)  =  P (t,P,z)  p  z  =  (j>i(t,p,z) +i<t>- (t,p,z  =  oi(t,p,z)  2  +  ia (t,p,z 2  Equations  The above approximation and assumptions, in conjunction with the Einstein equation and the equations of motion for the matter fields, yield two systems of equations: the time evolution equation for the matter field and the constraints and slicing condition for the geometry variables. The exact form of the equations must be derived in the chosen coordinate system, and we note that the C F C leads to expressions that are quite simple compared to those that would be needed i n a fully general-relativistic treatment.  2.6.1  Conformal Flatness  Again, following Wilson, Matthews and Marronetti [7], the spatial metric is given by: la = Tp ia A  i  ij  = 4>- Y 4  j  (2-99)  where 7^ is the flat 3-metric in cylindrical coordinates. A l l geometric quantities in the curved space can be expressed i n terms of their flat space equivalents:  19  Chapter 2. Analytical Work  the latter will be denoted with a For instance, the relationship between the curved and flat Christoffel symbols is given by:  =  \l  {dan + djlu -  =  f * + | {8 D^  kl  ) + | (Stdrf + 618*1, -  + 5 Dii, -  k  vtf'M)  %^D,fl>)  k  = ft+4  (2-100)  where the last expression defines the A j. This, i n turn, makes it possible to express the curved space covariant derivative i n terms of the flat covariant derivative: k  DiV = DiV + A k  k  (2.101)  k j jV  Further,using standard flat space calculus notation, we obtain: fibiDjf  =  V /  =  Vf-Vg  Y DifD j  j9  (2.102) (2.103)  2  The following properties of A j will be useful: k  4  =  m  (- )  4  = ^Dtf  (2.105)  7 4  =  (2-106)  7«4  = - ^ 7 * W  (2-107)  H  A  2  73,4 = 8 i 4 + f f 2.6.2  r o  ^-f^.-f^  104  (2.108)  The Extrinsic Curvature  The coordinate time derivative of the space metric, given that i n the A D M coordinates f* = (1,0), is given by: £tlij  =  t V„7y+7<*V,-t* + 7jyv"i**  =  Vo7i +7i r* +7fc/rf  M  4  j  fc  0  0  = iij On the other hand, recalling equation  =  (2-109) (2.13),  -2aK  iS  we have  + Difij + Djfii  (2.110)  Chapter 2. Analytical Work  20  This leads to: 7y = -2oK  + D$} + Djpi  ij  (2.111)  Since the space metric is assumed to be conformally flat for all time, the above equation can be turned into an equation for K^. Indeed, since the space metric is flat, its time derivative can be written:  7«=4^V7«  (2.112)  The time derivative of the conformal factor is solved for by taking the trace of the two above expressions. Then, using Y^Kij = 0, the form of the extrinsic curvature can be deduced: K  ij  = - L ( « * Z J ^ +y'*Z? ^ - | I>*0*)  (2.113)  y  7  f c  fc  7  This equation can be recast i n terms of the flat space covariant derivative. Using equations (2.101), (2.105) and (2.106) we find K  ij  =  ^(? Dk0 +¥ Dk0 -l'r Dkl3 ) k  1  k  i  ii  k  1  cop ^ 4  k  3  m  k m  >  "v  V  =  i  p  '  ^f*«  ( -114) 2  The covariant derivative of the extrinsic curvature can also be expressed in terms of flat space operators: DiK*  =  DiK  ij  + A K^ k  ki  + A{ K  (2.115)  ik  k  Using the properties of A j, this can be further simplified: k  A K k  i3  ki  =  ^(DpP)K  =  ^(4 f k  aip = Thus:  (2.116)  i3  m  D  p -\Ai f D k  m  k  k  m  r)  ^(DnP)K  3 i3  (2.117)  Chapter 2. Analytical Work  =  21  -L.bi^kV)  (2.118)  Finally, the trace of the square of the extrinsic curvature can also be expressed in terms of flat space operators: KijKV  =  =  l i m l j n  K  ^ ^ ( ¥  m n  k  K^  b  k  ^ +  ^ b p -^rb p )x k  i  k  ( A/3" + 7 'A/? n  ra(  k  k  ro  h D,p ) mn  7  l  o  (2.119)  2.6.3  The Lapse function  It will be useful to know how to express the double covariant derivative of the lapse i n terms of the flat space operator. Using equations (2.102), (2.103) and (2.107) we find  = =  2.6.4  ^(Dibja-A^b.a) ^  2  a  +  ^ -  V  a  )  ( - °) 2  1 2  The R i c c i Scalar  The Riemann tensor must also be expressed i n terms of flat space operators. Using equations (2.104) to (2.108):  =  d ^ + A ^  +  =  d ^ + A ^ e l  +  =  (diT^+t^e^  ^ + A ^ D ^ it^+Af^f^+A^ + idiAfM  + A^A^e^  .  +{t? A% + TZA? )e* k  =  k  Dibjet + ibiAf^ +{%)Kk  m  +  + 2t A^)e n  a  k(i  m  A^A^ (2.121)  Chapter  2. Analytical  Work  22  Permuting the indices i and j and subtracting, the Riemann tensor is obtained: Kij  =  R +2D AY +2A» A%  (2.122)  n  k ij  [i  ]k  k[j  n  Contracting the indices appropriately, the Ricci tensor is obtained. through the algebra yields:  Carrying  + ^D il>D fl> - ~^%T b ^b ^  (2.123)  n  i  j  m  n  Further, the Ricci scalar is given by the contraction of the Ricci tensor:  B u t R is the curvature scalar of flat space, which vanishes! Hence, the curvature scalar of the conformally flat 3-geometry can be written: R = ~V 1>  (2-125)  2  2.6.5  The Explicit Equations  Combining all of the above results, the actual equations to be solved can finally be deduced. Only the momentum constraint still needs to be manipulated i n order to yield equations for the non-vanishing components of the shift vector. Only the p and z components of the momentum constraint equation are of interest since the axial symmetry implies that the <j> component vanishes identically. As already noted i n equation (2.118), the covariant derivative of the extrinsic curvature is given by: DiK'1  =  ^D^KV)  =  (f D F  + ^Dtf*  k  k  -  l^DkP")]  (2-126)  For the flat cylindrical 3-metric the only non-vanishing Christoffel symbols are given by: fJ  3  = -P,  = i  (2.127)  It can be shown that i n this case the Riemann tensor vanishes and thus the flatspace covariant derivatives commute. Using this fact, the momentum constraint equations thus become: 167n/> J> 10  =  ^ ( ^ ) [ v a IL  7  ^  f  c  ^  +  ^  f  c  ^ _ | f ^ * ] 3  Chapter 2. Analytical Work  OL  =  L  23  3  Di ( ^ ) [ f D p  + l^DkP -  +^\f DiD 0i  + lfJDiD p ]  k  j  k  k  \i D p ] ij  k  k  (2.128)  k  k  k  ft L  O  Using the above values for the Christoffel symbols and noting that P^ vanishes by symmetry, one gets for i, j 6 {1,2}:  -  =  dB dP P ^~ + ^ - + — dp dz p  (2.129)  Dip  =  dtp '  (2.130)  Dp  =  tL.  (2.131)  p  t  Dp  k  k  j  3  3  z  p  3  P  This, i n turn, leads to:  f'DiD.pJ  = f (diD pi  -f? D pJ  k  k  rD p  +t{ D p ) n  n  n  =  f did p  -f'tl^p  =  f did p  +drP -5 ^  =  k  iDk  k  k  3  k  +s f fl b p  3  k  j  l3  3  3  (2.132)  lj  P  di[^  3  3  j  k  k  2  + ?£ + ^}  (2-133)  Using the above results, and standard flat space operator notation, the elliptic system for the geometry is given by: ,5  V V  =  -^{KijK^  Va  =  -jVlp-Va  =  -\dj(y • P) + ~P" 6 p  2  2  Vp 2  j  2,  + 167rp )  (2.134)  m  + a^iKijK +4Tr(p  + S ))  13  i  m  1>  i  (2.135)  + m/> (167r^') 4  -a*(in (^)) ( V * ^ + 6' d p k  i  k  | < F ( V • p))  (2.136)  with KijK given by equation (2.119). Noting that in the chosen coordinate system = ctpip , the explicit equations can now be written down. The matter evolution equations are (i = 1,2): lJ  6  £  ,  » ^  (  +  ^  -api> m <l>i e  2  )  +  »  (  /  , .  I  +  <  ^ » ) (2.138)  Chapter  2. Analytical  Work  24  The geometric equations are: dhp_  dhp_  dz  dp  +  2  +  CT? +a  lfhp  2  p dp  2  4p ip 2  ~i[{-dp~) <Pa dz 2  +  cPa dp  +  2  Ida p dp  dz  4rd pp_  1 dpP  3 L dp  p dp  2  +  2  \-dz-)  dz dz J  a  s  dp  (o\+e\)-a^m\<p\  i a /? 2  p  + <p )  (2.140)  2  4a 50i pt/> V " dp  2  3 5p<9z  2  V dz  +  ipldp 2  2  y~dp~) da dxpi  2  - - r ^ i t i + 4>l) (2-139) 4  +  2 rda dtp  + p ip 9^  7  6  dp J  1  (2.141)  dp 4 dp 3 dz 2  z  dP dp  1 a /3" 3 dpdz  2  2  1 dp" 3p dz  2  p dp  2  2 o  +  4a / d<pi pip' CTi dz —  ,i , a/j'  a/?  6  3az^  2  r  r  a ^ [ dp  l n (  2  p  dp  a  2  yS^-  p i CJZ r- + —  dPaz \  dP L  d<p "  hdz CT-2-  z  dz  dp J  (2.142) (2.143)  Kyi?'  1  4/dP 2o? 3 V dp  dp  2.6.6  dP <9z  p  z  +  dz  3Vp  V dz  )( <9p  +3Vp  dp J  (2.144)  az  Boundary Conditions  The above system of equations is valid for p > 0. It will also be useful to know what form this system takes on the axis. B y the regularity condition, i.e. demanding that the equations be smooth and finite on the axis, we get the following conditions on the axis: A(p,z,t)  =  A (z,t)  B(p,z,t)  =  pB (z,t)  0  l  + p A (z,t) 2  2  + p B (z,t) 3  3  + O{p )  (2.145)  4  +  0(p ) 5  (2.146)  Chapter  2. Analytical  Work  25  with: A  G  a,i>,p ,4> 4>  (2.147)  B  6  /3 ,a a  (2.148)  z  u  2  p  u  2  Hence, the system on the axis becomes:  ^  w  +  +  - To.  =  P  8  ~ r  c^a dp  c^a dz  2  2 +  19a pdp  _ ~  _2dadip i>dzdz 2a  /3 3 dz  2  dp  2  K--K  13  =  p  a  ^ s > i  0  )  ( 2  -  1 4 9 )  /<90 \ 1  2  2  2  t j  ^  i *™  Ki K  + "*) - ^ (<t>\ a  + 02)  m2  (2-150)  0  p dp  (2.151) 3  3p dz  1  pip V 6  /?"\  2  2Q2  i  V>IY<90l\  . +  +  +  2>,  / 2 ,  1  w  3 V dp  c>2  1  +  a 2  ^  /  ^dpPdp dp dz  2  dp dz z  3\ p )  V dp +  (2.153)  with the understanding that all terms involving 1/p remain finite on the axis, due to the regularity conditions. Finally, the outer boundary conditions can be obtained using the following argument: assuming the matter distribution has compact support, the solution inside must be connected to a vacuum solution outside the matter distribution. The system is easily solvable for a spherically-symmetric, vacuum spacetime. Far enough from the central region the matter distribution will seem approximately spherical and the inside solution should be joined to a spherical vacuum solution. In spherical symmetry, the system of equations in vacuum becomes: da dr <Pip dr 2  2  2  2 da r dr 2 dip r dr  2 da dip ip dr dr =  0  (2.154) (2.155)  Chapter  2. Analytical  Work  26  The general vacuum spherical solution, taking into account that the geometry should be flat at infinity, is thus: ^(r) = 1 + r a  ( )  =  r  1  (2.156)  r+ A  +  ( 2  1 5 ?  )  Here A and B are integration constants. Eliminating these integration constants from the general vacuum solution, we get the following conditions:  %  +  -  0  P"  =  0  (2.160)  P  =  0  (2.161)  z  '  <  2  -  1  5  8  >  B y imposing these conditions on the interior solutions at the outer boundary, the solutions will match smoothly with the vacuum solutions.  2.6.7  Spherical limit  It will be useful to know what the above system reduces to in the case of spherical symmetry. B y making the following change of coordinates, where (j> is the spherical polar angle and Oi is the conjugate momentum to the scalar field obtained from the Lagrangian in spherical symmetry: p — r sin 4> z = r cos 4> P =P sin(f> p  P  r  = P cos 4>  z  r  r sin <p A  and requiring that all grid functions depend only on r and t, the spherical equations are deduced:  ^  +  ?^ = - | £ ^  +  Ir(f -7 ) :  :  2  +^ & + & 5V  25V  V  5  -  »  (dP  r  2  m 2  ^ ^ +  ' Py  1  r  + ( f )  2  ]  (2-162) V'  5  2  2  o  r  ,  2  +  ,  2  ,  (-3)  Chapter d^T dr 2  2dp_ r dr  =  jr_ _ 3 a r V / ^ _ ^ v dr  6  2  +  f  2. Analytical / v \ / r  1  Work  27  dfa cdfo\ dr dr ) l _ 9 a _ 6 5V>\ J Va dr ip dr ) ta  -  <»»>  m ctip 4>i 2  (2.166)  6  with boundary conditions: r = 0:  |=0.  £ = 0 , ^ = 0.  f  =0,  r ^ + V - l = 0 ,  n/>^+a-1=0,  6=0  (2.167)  /T = 0  (2.168)  The extrinsic curvature is given by: 2 V > \dp 4  j3  r  4"t> =  r ~Y  Krr  (2.170)  K  sm 4>K^  (2.171)  2  K  ee  =  2  Of course, this system is exactly what the constraint equations and evolution equations reduce to in spherical symmetry when the conformally flat condition is imposed; in this case, it is not an approximation but a perfectly valid coordinate choice.  2.7  Initial D a t a  The simplest self gravitating system that can be studied is the static spherical boson star. It is thus an appropriate start for the numerical study of the C F C approximation. Further, it is critical that the initial data used in the collisions be that of stable boson stars. Indeed, if the initial data was simply specified as random lumps of scalar fields, these lumps could simply fly apart due to poor initial data, and not due to actual gravitational interaction between the two bodies.  2.7.1  Single Static Spherical Boson Star  Boson stars are discussed in the article by Jetzer [18]. A complex massive scalar field governed by the Klein-Gordon equation can be molded into a static body  Chapter 2. Analytical Work  28  with compact support called a boson star; for such a configuration of the field, the dispersive nature of the field is exactly counterbalanced by its gravitational self interaction. Boson star solutions are most readily computed i n Schwarzschild-like coordinates. The metric i n this coordinate system is given by: ds = -a(r) dt 2  2  + a(r) dr + r d f i  2  2  2  (2.172)  2  and the matter field is described by the following ansatz: $ = 4>(r)e  (2.173)  iwt  This ansatz, due to its harmonic time dependence, results in a time-independent stress tensor and thus a time-independent geometry. The Einstein equation plus the standard expression for the Lagrangian and the stress energy tensor yield the following system: da dr da. dr d^ dr 2  2  ,2^3^,2 a a 2r~~2r aa a. ~2r ~ 2r fa 1 -(—+ - -m ra 4> )^+m a 4>(l--^) Vr r J dr V 3  +  =F('+^)+?(g)  2  <"">  2  2  2  2  2  2  (2.176)  2  ma / z  l  where a = wa was used. As previously stated (see equation (2.145)), regularity at the origin demands that all first order derivatives vanish at r = 0. Further, it is clear that for the equations to remain finite at the origin the condition a | o = 1 must be respected in order to avoid a conical singularity. Hence, close to the origin, the functions must have expansions: r=  a(t,r)  =  1 + a (t)r + 0(r )  a(t,r)  =  a {t) + a (t)r  0(t,r)  =  Mt)+<t>2(t)r + 0(r )  2  (2.177)  4  2  + 0(r )  2  0  (2.178)  4  2  2  (2.179)  4  In order for the scalar field configuration to represent a boson star, it must effectively have compact support, i.e. the amplitude of the scalar field must die off exponentially. This results in an eigenvalue problem for the value of a\ =oIndeed, only a discrete set of values for this parameter will yield a boson star configuration. Demanding asymptotic flatness in the geometry turns this into an eigenvalue problem for w. Indeed, the system of equations outside the matter distribution becomes: r  ± dr da d  =  f(l-a ) 2r ' a „. 2  v  y  (2.180) '  Chapter  2. Analytical  Work  29  Here C is an integration constant. In the limit r —> co, asymptotic flatness, as well as imposing the coordinate condition equating coordinate time and proper time at infinity demands that both a and a be equal to 1. This implies C — 1/w and w = lim (2.183) ]—yoo act Once this eigenvalue is obtained, a single spherical static boson star is obtained. The parameters describing the stars are the amplitude of the scalar field at the origin <p and the time frequency eigeinvalue w. 0  2.7.2  Transformation to Isotropic Coordinates  The above boson star was expressed in the Schwarzschild coordinate system; the result must be transformed to the isotropic coordinate system i n order to agree with the coordinate system chosen for the dynamical evolution. Preserving spherical symmetry, the two metrics compare: ds  2  =  -a(r) dt  =  -a{R) dt +ip(R) {dR  2  2  2  + a(r) dr 2  2  A  2  2  + r dfl  ,  (2.184)  + R dtt)  (2.185)  2  2  Here R is the new radial coordinate, and a and ip and the lapse and conformal factor respectively. Assuming the coordinate change is of the form R = R(r) and equating the angular dependence in both coordinate systems, the following equations are obtained: r  =  Rip  (2.186)  adr  =  ip dR  (2.187)  2  2  This leads to: adr dR — =  (2-188)  In order to turn this into the sought relation, R = R(r), boundary conditions are needed. We demand that the two coordinate system agree at the asymptotically flat infinity, i.e. Wmr^^R/r = 1. Outside the matter distribution, the form of a is known: ^outside = ( ! - — ) Here M is the Schwarzschild mass of the star. relation between R and r is given by: r  R  am  (2-189) Using equation (2.188), the  r a(r)dr  'Ro =>R(r)  =  E o e x p f - ^ ^ ^ ]  (2.190)  Chapter  2. Analytical  Work  30  Figure 2.1: Schematic of a trapped surface. S is the trapped surface, s is the spacelike normal to S, n ° is the timelike normal to E and l is the outgoing null vector from S. a  a  Considering coordinates outside the star and using a f" r  R(r t) ou  —  Ro exp  R(r ut) Tout 0  Here r 00:  o u t  1  'r0  VV  - ^  Ro  2  0UtS  jde  :  df  i  -  2Mfi  + 11  2  (2.192)  ' out - ^ro + 1  To  (2.191)  is a value of the coordinate outside the star. Taking the limit r t —> o u  Ro = \ (ro ~ M + ^r  2  - 2Mr )  (2.193)  0  Thus, choosing r outside the star and knowing the value of M, the relationship between R and r can be established: 0  R(r) = \ ( r - M + ^ r 0  2  - 2 M r „ ) exp [ - £ "  (2.194)  Knowing the form of this relationship, $(R),a(R) and ip(R) can immediately be deduced. Finally, this spherical configuration can easily be expressed in a cylindrical coordinate system by noting that R = \/P + z . 2  2.8  2  Apparent horizons  Assuming cosmic censorship, singularities of spacetime are hidden behind event horizons. Unfortunately, it is impossible to determine the position of the event horizon unless the full solution of the Einstein equations is known; the event horizon is a global entity. What one can do in a dynamical spacetime containing a singularity is to pinpoint the position of the apparent horizon, which is the outermost trapped surface on a given time slice. This essentially means looking for an instantaneous approximation (something that can be calculated on a given time slice) to the event horizon. If an apparent horizon is present, and again assuming cosmic censorship, it is then known that there must exist an event horizon, and that the spacetime thus contains a black hole.  Chapter  2. Analytical  Work  31  Following the presentation due to Pretorius [19], an apparent horizon is defined as follows. Consider a spacelike closed surface S with spacelike unit normal s on a given time slice E. The outgoing null vector from this surface is given by: a  l  =  ll  =  a  a  a  +s  a  (2.195)  a  n  0  (2.196)  Define the outgoing expansion parameter K+ as the projection on the surface S of the covariant derivative of the outgoing null vector: K+ = (  - s s )V l  Q b  a  7  (2.197)  b  a b  This, i n effect, is a measure of how lightlike geodesies leaving the surface towards the outside diverge from one another. A trapped surface is a surface for which K+ < 0, which indicates that outgoing geodesies converge! The outermost marginally trapped surface, for which K = 0, is called the apparent horizon. The presence of an apparent horizon is a local property ( i n time) and (again, assuming cosmic censorship) signals the presence of an event horizon. It should be noted that an event horizon may intersect a t = const, surface without the existence of an apparent horizon on the slice. Setting K+ to zero, a differential equation denning the apparent horizon (or more properly, a marginally trapped surface) can be derived: +  =  h " -s s )V (n  =  {-f -s s )(-K -n a  + Vs)  =  ( ° -s s )(-K -n a  + Ws)  =  -K + s s K ,+  =  -K + s s K  a  a  a  ab  a  +s)  b  b  b  b  ab  &  a  a  b  a  b  7  ab  a  a  a  b  7  a  b  a  b  (2.198)  a  b  b  s s )W s  +Ds  ab  a  b  ,( * -  ab  a  b  a  Here equation (2.21) and the fact that s is a space vector were used. Hence, since the trace of the extrinsic curvature vanishes due to our choice of time slicing, we have: a  Ds  a  a  +ssK  ab  a  b  =0  . (2.199)  In order to solve the above equation, i t is useful to view the apparent horizon surface S as the 0-level surface of a scalar function, F i.e.:  S: l y e E I F f y ) = 0}  ,  (2.200)  Then, the spacelike unit normal to the surface is given by the normalized gradient: 1  s*  =  ^  =  N  =  y/6 d Fd F mn  m  n  (2 201) (2.202)  Chapter  2. Analytical  Work  32  Since the spacetimes we are concerned with are axially symmetric, the projection of the apparent horizon surface, for anyfixedvalue of the azimuthal angle, is described by a curve starting and finishing on the axis in the p z-plane. The apparent horizon is then the surface of revolution of this curve about the axis. Regularity on the axis demands that the tangent to the curve be in the p direction on the axis. The one dimensional curve in the pi;-plane is most easily expressed through the change of coordinate:  r = v V + (z- z )  9 = tan  2  0  -1  (—?—) \ Z —  (2.203)  ZQ I  The surface is thus described by the function:  F =>r  = r - R{9) = 0 = R(9)  (2.204) (2.205)  Using the above change of coordinates, the partial derivatives of F are given by: F  p  F  tZ  R' =• sin(0)-cos(0)— R R' R = cos(9)+sm(9) — 2am(«)coe ^ W  2  cos(0) R „ - ,^  (2.206) (2.207) - A i ( l - - )  +  cos (0)^ , R" (l+^) R ^ ' R , ™ ' sin(6>) / 2  r  R  (2-208) R"\  2  R  R  V  R )  Here R! and R" are the first and second derivatives of the function R(9) with respect to 9. Using these results, equation (2.199) becomes: 0  =  D^  +  =  ds + T s i  SiSjK^ i  i  j  ij  +  s K  ij  iSj  = HV^ w '> ^ ' ' ' s  Nip  2  +  T  dkF+  >didjF - ^ W d j F N ip 2  2  1  3  -  ^ ^ A F + ^ d i F d j F  Nip  K> d Fa F  3  ^d^djF (2.211)  Chapter  2. Analytical  Work  33  The Christoffel symbols for the conformally flat space are easily calculated we then find  -J^^Nip+ ^S^di^diF  -H^TijdkF = Nip 2  Nip p  %3  2  3  3  (2.212) '  v  Premultiplying by Nip , the differential equation becomes: 2  0  =  \S' didjF - —diNdjF +^\  + [-5 ditl>djF + jjrIC>diFdjF\  3  t3  =/  =T  (2.213) To proceed further, the derivative of N must be known: N  =  yjF% + F%  (2.214)  1  ftij  ^—diNdjF  =  ^[F F  + 2F F F +F F ]  2  (2.215)  2  pp  z  p  pz  zz  Use of equations (2.206) through (2.210) reduces the above expression to:  Thus, the term independent of the geometry in the differential equation becomes: I  =  5 d d F+^-—d Nd F l3  i  1  r,  #"i  i  i  1 r  ,  i  1  /m^'i  r (T/)  12  2  t + (f I'] [2* - ™t(e)R'} + ^ 1  R +R 2  j  i  r,  - ii"]  (2-217)  To compute the source term, T, i.e. the term involving the geometry of the spacetime under consideration, the extrinsic curvature must be written down explicitly. Using equation (2.114) we find: "  K  R  Z  Z  =  = =  l ^ l W - K - j ] d  ^  2^K  -  *  + / ?  '  <J  -  ^  (2-218) ( 2  '  2 1 9 )  ( 2  "  2 2 0 )  Chapter  2. Analytical  34  Work  Thus, the source term becomes: T  [Itf'diiltdjF+^KVdiFdjF]  =  •4 (20* - 0P - ^)F* + F, F {p% + 0> )1 p  z  (2.221)  tP  Hence, the final expression for the differential equation describing the apparent horizon is given by: R"  =  [l + ( ^ ) ' l [2R 2  R  cot(9)R'}  +  R  ^-  R  +  [R  2  (2.222)  + R»]?  The boundary conditions are given by the regularity condition on the axis: R'(9 = 0) =• 0 R'(9 = TV) = 0  (2.223) (2.224)  In spherical symmetry, the trapped surface has to be spherical, and thus the normal is given by: t  _  I 9j Fgp/t 13  \jY^ d F ptid F pf n  n  =  S  m  S  l  p  (2-225)  This leads to the apparent horizon equation: 0  =  DiS  =  dis'  +  {  SiSjK  ij  rij+SiSjK*'  +  dr\ip ) 2  2_chp ip dr 3  +  +  ip \ip 2  dr  +  r)  +  ip  A  (  J _ / 6 dj)_ 2\  _2_/d/^ _ F\  ip \ipdr  3a V dr  2  rJ  '  '  ( 227) 2  r /  Thus, in spherical symmetry, when the function defined by the right hand side of the last expression has a zero crossing, an apparent horizon is present.  2.9  A D M Mass  There are many ways of defining the "mass" of spacetime. The ADM mass, as defined by R. Wald [2] is given by:  Chapter  2. Analytical  35  Work  = 22. is 1 1 [Sf - 0 ] " ' * *  <- > 2 228  t,i=i  where the integration is to be performed on a closed surface in the asymptotically flat region that contains all of the matter in its interior. In this case, the ADM mass expression reduces to:  MADM  -  -  -  j  ip ^,zpdp+ 5  I J Z  il> ip, pdz 6  P  m  i „  ,p=Pmax  j  ip i{),zpdp 5  (2.229)  Chapter 3. Numerical Implementation  36  Chapter 3  Numerical Implementation The numerical solution of the problem at hand is divided into three parts: generating the initial data, which i n this case will involve the spherical boson stars, solving the elliptical system for the geometry by using finite difference techniques and multigrid methods at each time step, i n effect generating a time evolution for the geometry, and solving the time evolution equations for the matter fields, also using finite difference techniques. ' •  3.1  Initial Data  The initial data problem is straightforward. For spherical static boson star initial data, the problem is one dimensional. The Schwarzschild system in spherical coordinates is put into first order form and the value of a | o is used as a shooting parameter (see equations (2.174), (2.175) and (2.176)). Tuning this parameter produces a star centered at the origin with the scalar field asymptotically reaching zero amplitude. A decaying exponential is then attached to the trailing scalar field to ensure it smoothly goes to zero. The spherical data thus obtained can then be transfered to the cylindrical grid. A single star can be considered for stability testing or two identical stars can be put on the grid a distance apart to study head-on collisions. r =  3.2  Definitions  3.2.1  Two Dimensional Domain Discretization  Once initial data is given, the two dimensional problem can be approached. The domain is rectangular, given by 0 < p < R and — - R m a x < z < Rmax- The numerical domain is a discretization of the domain with NI = N = 2 +1 points i n the p direction and by N2 = N = 2 + 1 i n the z direction. Here, "level" is a parameter that measures how many discrete points are used to approximate the domain . For example, at level = 3, the grid will be approximated by 9 x 17 points and at level 7, it will be approximated by 129 x 257 points. The grid points are identified by two integers, typically i and j. Hence, the grid point designated by (i, j) is located at p(i) = ih, z(j) = — i ? + jh, where h = R ax/{N — 1) is the mesh spacing. A s can be seen by the above definitions, the mesh spacing is the same i n both directions. m  a  x  l e v e l  p  l e v e l + 1  z  m a x  m  P  Chapter 3. Numerical Implementation  37  Rmax  P  N2  Rmax >  -Rmax NI Figure 3.1: Domain in the pz plane  3.2.2  Discretization and Finite Difference Operators  Following Ascher et al. [14] [15] [16] , the differential equations must be transformed into finite difference equations i n order to be solved numerically. To do so, all the differential operators are transformed into finite difference operators which agree with the differential operator up to second order i n a Taylor expansion, i.e. all schemes used will be of second order. A l l the finite difference operators will either be centered or one-sided, the one-sided operators being used exclusively on the boundaries. Let u be any grid function of interest and p, z the cylindrical coordinates. The centered operators used are given by: du -^^D (u) p  =  1 — [ (i + l,j)-2u{i,j)+u(i-l,jj]  (3.1)  0 - ^ ( 1 * )  =  ^[u(i,j  (3.2)  2  2  u  + l)-2u(i,j)  + u(i,j-l)]  [ ^+u +i  ) - ^ - y  1  gjit^ %M = D  w  -u(i + l,j-l)  + u(i-l,j-1)]  + i) (3.3)  ^^D (u)  =  ^[u(i + l,j)-u(i-l,j)]  (3.4)  ^^D (u)  =  l.[ (i,j  (3.5)  p  z  u  + l)-u(i,j-lj]  The forward one-sided operators used are given by: DIBDYM = -^[u(4,j)-4u(3,j) + 5u(2,j)-2u(l,j)]  (3.6)  Chapter 3. Numerical Implementation  38  KBDYM  =  - ^ [ « ( i 4 ) - 4 u ( i , 3 ) + 5u(i,2)-2u(i l)]  D ,BDYI(U)  =  -^[u(3,j)-4u(2,j)  + 3u(l,j)]  (3.8)  D ,BDYI(U)  =  -^-[u(i,3)-4u(i,2) 2h  + 3u(i,l)}  (3.9)  P  Z  >  >  (3.7).  (3.10) and the backward one-sided operators used are given by: DIBDY2(U)  1  =  -^[u(Nl-3,j)~4u(Nl-2,j)  h  2  +5u(JVl-l,j)-2u(JVl,j)] D ,BDY2(U)  =  2  Z  (3.11)  -^[u(i,N2-3)-4u(i,N2-2)  h  2  +5u(i, N2-1)-  2u{i, N2)}  (3.12)  D ,BDY2(U)  =  ^[u{Nl-2,j)-4u(Nl-lJ)+3u(Nl,j)]  (3.13)  D ,BDY2(U)  =  ^[u(i,N2-2)-4u{i,N2-l)  (3.14)  P  Z  + 3u(i,N2)]  The terms involving ^ but known to remain finite on the axis are put i n finite difference form via: -> -^[u(3,j)-8u(2,j)  ^Uxis -\ s p  axi  + 7u(l,j)]  -»• - ^ K 3 , i ) - 8 « ( 2 , j ) ]  (3.15) (3.16)  These expressions are derived by assuming the following/behavior near the axis: u  =  A + Bp + Cp + ...  (3.17)  v  =  p(A + Bp + Cp ...)  (3.18)  In particular, u € {a,tp,P ,&}, z  I0 =° ^=° lp  2  4  2  v € {/3 ,Oi}. p  4  Finally:  =  -  - i ^ ^ ( 3 . i + l ) - 8 W , i + l) -(3 (3,j-l) p  3.3  + 8(3 (2,j-l)] p  (3.19)  Solving the Elliptic System: Multigrid Methods  For a given matter distribution, the geometry of a spacelike slice of the spacetime can be determined by solving the elliptic system presented i n the previous chapter. The algorithm of choice for a problem of this nature is the multigrid method, first introduced by A . Brandt [4]. The algorithm relies on relaxation of the discretized system and on transfers of grid functions between grids with different mesh spacing.  Chapter 3. Numerical Implementation  3.3.1  39  Discretization on the G r i d  The equations are discretized to second order using mostly centered difference operators. Ghost points are introduced all around the boundary and they are eliminated from the difference equations by using the discretized boundary conditions obtained from equations (2.158) and (2.159). Ghost points are a virtual extension of the existing numerical grid; when the lattice defined by the discretization of space is extended beyond the original domain, the ghost points are the first virtual points encountered directly outside the domain boundary. The only troublesome points are the corners: the inner corners (on the axis) have enough conditions to eliminate all ghost points (the regularity condition and the outer boundary condition). This idea cannot be applied to the outer corners, however, because there is only one equation, the outer boundary condition, available but two ghost points next to the corner, which is why backward/forward differencing is used i n the z direction. Finally, terms i n 1/p known to remain finite on the axis are discretized using the operators defined in the last section.  3.3.2  Relaxation  The numerical grid is split into a "red" grid and a "black" grid, in analogy with a chess board. If i + j is even, a point is said to belong to the red grid and if i + j is odd, the point belongs to the black grid. A full red-black relaxation step consists in a relaxation sweep on the red subgrid followed by a relaxation sweep on the black subgrid. In the context of multigrid, the role of relaxation is to smooth the error and Gauss-Seidel relaxation with red-black ordering is a good way to reach that goal. A t each grid point, the discretized equations can be viewed as a set of four (three on the axis) algebraic relations i n the grid functions at that point. Hence, at each point we want to solve: P(v)  =  0  (3.20)  v  =  (a(i,j),i,(i,j),p»(i,j),p (i,j)) z  (3.21)  Following a strategy suggested by Choptuik [22], the Newton method is employed to approximate a solution to this equation. This method consists in applying the following iteration until a certain tolerance is attained (in practice, a single iteration step is enough). v"  +1  J  =  tT-J-ip"  = «w.fl.*.fl>  (3.22)  (3.  23)  d{vi,v ,v ,v ) 2  3  i  Here J is a "jacobian matrix" and J is its inverse. As is well known, given a good initial guess this method reduce the error in the approximation quadratically. Gauss-Seidel red-black relaxation reduces high frequency error modes down very efficiently, but leaves lower frequency error mode virtually untouched. If relaxation was the only solution process used, it would take an inordinate amount -  1  Chapter  3. Numerical  Implementation  40  of operations(i.e. 0(N ) operations where N is the total number of grid points) to reach a decent approximation to the solution. Multigrid, on the other hand, can solve these problems in O(N) operations. 2  3.3.3  Transfers and the F A S Algorithm  The described relaxation is a good smoother, but a bad solver, since low frequencies are slow to converge. Fortunately, the meaning of "smooth" depends on the mesh size. B y transferring the smooth data from a fine mesh to a coarse mesh, what was low frequency error becomes high frequency on the coarse grid and can efficiently be relaxed away. This is exactly what the multigrid algorithm relies on to solve the difference equations rapidly. A lengthy discussion of multigrid methods would take us too far off course, but we note that since the problem solved is essentially non-linear, the F A S method is used, meaning that a full approximation to the solution is being manipulated at each multigrid level. The coarsening is done by simply doubling the mesh size i n both dimensions. Following Brandt [4], the grid functions themselves are transmitted from fine grid to coarse grid using a full-weighted injection operator. The defects, which are zero on the black grid and vary smoothly on the red grid, are transmitted using a half-weighted injection operator. Finally, the approximations to the errors are transmitted from a coarse grid to the finer grid using the bilinear interpolation operator. A l l these operators basically perform weighted averages over neighboring points and are defined in the caption of figure (3.2). Using the multigrid method applied to the discretized elliptic system, the geometry of the spacetime can be determined efficiently from the matter distribution.  3.4  Solving the Hyperbolic System: Time Evolution  Ideas relating to the numerical solutions of hyperbolic equations are presented in Ascher et al. [14][15][16] . The time evolution of the matter field is given by equations (2.137) and (2.138). the index i = 1,2 for the matter fields will be suppressed, with the understanding that the following applies to both the real and the imaginary part of the scalar field. Further, the function names will now designate their gridfunction approximations, and not their actual continuum counterparts.  3.4.1  Interior Difference Scheme  As suggested by Choptuik [22], the Crank-Nicholson scheme with Kreiss-Oliger dissipation is used to obtain a finite difference version of the hyperbolic system. This scheme is second order in both time and space. Kreiss-Oliger dissipation is equivalent to introducing a term of the form h -^{u + u ) to the equations,, where e is a small parameter and u is (j) or a as the case may be. 4  pppp  zzzz  Chapter 3. Numerical Implementation  41  The Crank-Nicholson scheme can be summarized as follows. For a continuum equation of the form du T t  =  £  u  the discrete system is given by „ " + i _ -Tn  - n 4- -Tn  n+1  n  n  where L is the discrete equivalent of the operator £, k is the mesh spacing in time, h the mesh spacing in space and n denotes the time level. Define the following operators: KO{u)  = -^[u{i + 2,j)-4u{i +  l,j)+6u{i,j)-4u{i-l,j)  +u(i - 2,j) + u(i,j + 2) - 4u(i,j + 1) +6u(i,j) - 4u(i, j - 1) + u(i,j + 2)] AVG{u)  = |[u (i,j)+u n  n+1  (3.24)  (i,J')]  (3-25)  Then the following implicit finite difference equations are obtained: 4> {i,j) = 4> {i,i) n+1  n  +kAVG[^ o (i,j) n+1  =  + P"D (cf>) + p D {cj>) - KO(<j>)]  (3.26)  z  p  z  a (i,j) n  +kAVG[/3»D (a) p  +  +o{D {p") + D (0*))  P D (CJ) Z  p  z  z  +ai> {pDfa) + pD (cf>) + D (cf>)) - ap^m 4> 2  2  2  z  p  +2 nl,(D (il>)D (<l>)+D {1>)D (<t>)) af  p  p  z  x  +p^ {D {a)D {4>) + D (a)D (ct>)) - KO(a)] 2  p  3.4.2  p  z  z  (3.27)  Boundary Conditions  On the axis, regularity conditions prescribe what form the boundary conditions should take: dp a  = 0  By introducing ghost points next to the axis, and using the regularity condition to solve for the value of the scalarfieldat that point, it is possible to use the inside equation on the axis by simply setting the derivative in the p direction  Chapter 3. Numerical Implementation  42  to zero. Hence, the boundary equations on the axis are given by:  4> \l,3) n+  <p (l,j) n  +kAVG[^ cr  0  +0 D (4>) -KO{4>)} Z  Z  (3.28) (3.29)  The equations for the outer boundaries are the same as the equations inside the domain but the one-sided boundary operators replace the centered operators where it is appropriate. It is understood that the KO operator is only used in the direction parallel to the boundary, i.e. only in the direction in which it is well defined.  3.5  Obtaining Solutions  The time evolution code was generated using R N P L , the higher level language developed by Choptuik and Marsa [23]. R N P L was designed to solve hyperbolic equations and it is thus used to solve the time evolution system for the matter. The geometry is updated using an add-on multigrid routine designed and programmed by the author. The principal parameters used by the program are the tolerance (the maximum value of the defects tolerated by the solver when iteratively solving the equations), the dimensions of the domain, the value of e, the Courant factor (the ratio of the time meshsize and the space meshsize, k/h) and the maximum number of time steps to be taken.  Chapter 3. Numerical Implementation  /  /\  r /  43  \  J / \  /\  \  Figure 3.2: Coarse gridpoint (circle) surrounded by its fine grid neighbors (squares and triangles). The full-weighted injection operator attributes a value to a grid function at the coarse point by summing its values over the fine points with weights \ for the fine point at the position of the coarse p o i n t , | for the triangle points and ^ for the square points. The half-weighted injection operator attributes half the fine grid value of the grid function at the circle point as the coarse grid value of the grid function at the circle point. Note that this is equivalent to the full-weighted injection operator i n the limit that the values are zero at the triangle points and equal at the circle and square points. The bilinear interpolation operator, which allows to transfer information from a coarser grid to a finer grid, distributes the coarse value at the circle point on the finer grid with following weights: 1 for the circle point, | for the triangle points and \ for the square points. The actual value of a grid function thus transmitted at a given fine point is the sum of all the coarse point values distributed on that point.  Chapter 4. Scheme Analysis  44  Chapter 4  Scheme Analysis Finite difference systems, in a sense, are richer than differential systems. Even if a differential system is well behaved, a bad finite differencing of it can lead to spurious results. The following analysis, following Ascher [14] [15] [16], provides some necessary conditions for the chosen schemes to be stable and convergent.  4.1  Multigrid Algorithm  The elliptic system governing the geometry is quite complicated. Second order differencing has been used throughout, so we expect the truncation error to be of second order in the mesh spacing. A sharp and comprehensive analysis of the finite difference system would be very arduous, if at all possible, and outside the scope of this thesis. There is, however, a useful necessary condition that can be derived with respect to relaxation.  4.1.1  Necessary Condition for Relaxation  Let us consider the flat space case, i.e. a matter-free initial condition. Then, the lapse equation only involves the lapse function and the problem becomes fairly simple. The matter fields all vanish, and the geometry is given by a •= 1, ip = 1 and /3 = 0. It is crucial that relaxation be numerically stable with respect to this simplified problem, or else there is no hope for it to be a good solver for the general case. W i t h this simplifying assumption the elliptic system becomes: (4.1) (4.2) and the outer boundary conditions are the same as before, i.e.:  0 Clearly, the ip equation is uncoupled and can be solved on its own, without requiring reference to a. Once ip is known, a can be deduced.  Chapter  4. Scheme Analysis  45  Using centered differencing, relaxation inside the domain is given by: ^n+x  {i  j  }  V(> + 1,7) + ^ " ( » ~ 1,J) + ^ " ( » , J + 1) + V>"(M - 1)  =  + ^(r{i  + lj)-r(i-l,3))  (4-3)  q " ( i + 1, j ) + q " ( i - 1, j ) + ct (i,j n  n + 1  (V*^ ^ +  -  1  - i,i)))  4 ^ ^ ( ^ M • (^  n + 1  + 1) + a"(t, j - 1)  (t,j + l ) - ^  + l)-«"(M-l))x  n + 1  (i,i-l)))  (4-4)  Note that the superscript n in the above expressions denotes the relaxation sweep iteration, not the time level. We define the errors in the grid functions, and e™ v i a  V  =  a"  =  1+  4  (4-5) (4.6)  l + el  and rewrite the above iteration procedure in terms of these defects. Working to first order i n the errors, we get:  e  „+i (  e (i + 1, j ) + e"(i - l,j) + e (i,j + 1) + e"(i, j - 1) n  j  }  =  n  + ~ ( e ( i + l,j)-e (i-l,j)) n  n  (4-7)  which is valid for both e^, and e . If we consider a given Fourier mode of the form: e" = e e ^ ' \ we get ( do not confuse the grid index i and i = the index always appear i n grid functions): Q  n  l  , p + ? i 2  Clearly, for /o(i) > h, i.e. inside the domain, the norm of the amplification factor is smaller than 1 and relaxation reduces the error. The same analysis can be carried through at the outer boundary, but there's a twist: the ghost points introduced i n the centered differencing of the inside  Chapter 4. Scheme Analysis  46  equations must be solved for in the boundary conditions. Let's consider the P-— Pmax boundary for definiteness. Using the boundary conditions, we get for the ghost points: </>(iVl + l , j )  =  ^(Nl-l,j)  +  -^ (l-il (Nl,j)) )  )  z(i) p(Nl) {i>(Nl,j + l)-TP(Nl,j-l)) U  (4.9)  1  NI (a(Nl,j  p(Nl)  + l)-a(Nl,j  -1))  (4.10)  Now, we substitute the above result into the interior equation. Focusing on the ip equation, we get:  =  4p + 2hp + h ( 2  2  ~  2r{Nl  h  j  ) +  >  r{N1  j  +  1  }  +r(N,j-i)) 2 ( 4 p  +  2  y  2  l p  +  ^ ) (  2  /  t  - ^ " ^ ' ^  1  )  -i; (N,j-l)))  (4.11)  n  B y introducing the error and looking at a specific Fourier mode, we get:  = M  MNl?  Ap(Nl)  2  (cos(M+e-^)  + 2hp(Nl) + h  2  z(2p(Nl) + h) h + 2p(Nl)(2p(Nl) + h) 2  sin(£,))  (4.12)  We see that the second term above contains a term of the form z/(2p). If our scheme is to be stable, this term must be smaller than 1, or else the amplification factor will be greater than unity and the error will blow up! We thus obtain the following necessary condition: < 2p„ This condition prohibits the use of a very narrow grid where z is why the domain used is 0 < p < R , -R x < z < R the above condition.  max  max  4.2  ma  max  > p ax • This which respects m  T i m e Evolution System  The discretized system commanding the time evolution can be analysed using standard techniques. We will consider the geometry as given and we will perform "frozen coefficient" analysis.  Chapter  4. Scheme  Analysis  47  Scheme Order and Truncation Error The truncation error is calculated by substituting the continuum solution in the difference scheme and then performing Taylor expansion. As discussed i n the previous chapter, for the continuum equation f  =' < «  the Crank-Nicholson scheme is of the form /  „+i _k  Lfn+1  =  f  n  k  +  Lfn  where / is the function whose approximation we seek and L is the discretization of the differential operator £ . We thus get: r =  -W  + *>P>*) +  W>P>*)  l (f(t + k,p, z) + f(t, L  +  Knowing that L is second order i n the mesh size, we deduce: T  =  -(/ + kf + - f f a + iftu) + f l_ k  t  r = ~(ft +  \fu  +  (  £  +  0  (  f  t  2  )  )  (  (  /  +  f  c  /  t  +  ^  /  t  t  )  +  /  )  +  0  (  f  c  3  )  + £(f + \ft + *fu) + 0(h ) + 0(k ) 2  +  r = (-ft + £(f)) + \(~ftt + £(ft))  3  + 0(h ) + 0(k ) 2  2  Now, using the original differential equation, we have  0(K ) + 0(k ) 2  r =  2  and we thus see that the Crank-Nicholson scheme is of order (2,2).  4.2.1  Amplification M a t r i x and V o n Neumann Condition  We now proceed to a Von Neumann analysis of the Crank-Nicholson scheme— such analyses generally give rise to necessary conditions for stability of a difference method applied to a time dependent problem. We first define the following coefficients: Cl  =W  C  2  =  C  3  =  C  -  (4 13)  p"  (4.14) (4.15)  4  =  apV  C  b  =  a<p +pD (a)ip  C  6  =  pD (a)ip  C  7  =  -m apij  8  =  D (f3») + D (p )  C  (4.16)  2  2  2  p  2  z  2  p  +2poal>D (ilf) p  +2patpD (iP) z  (4.17) (4.18) (4.19)  6  z  z  (4.20)  Pf  z))  Chapter 4. Scheme Analysis  48  In terms of these coefficients, our discrete system is given by  ,  n+l  . / , \n (4.21)  where  k( 2  C D + C D - KO 2  V CV  + C D 5  P  3  Z  + C £> + C4CD + 7J ) 2  2  P  6  2  Ci C + C D + C D - KO 2  8  P  3  Z  (4.22) B y going to Fourier space, we can diagonalize the difference operators. applying the operators to single Fourier modes, we get  (4.23)  ~* h  D p  D  h  z  By  shife)  (4.24)  1  (4.25) (4.26)  KO  |({l-cosfc)) + {l-costo)f) 2  (4.27)  We thus obtain ( M =  -KO  + {(C  2  sin(£ ) + C sin(&)) p  3  C -KO  C + ^ ( s i n ( ^ ) + sin (^)) 2  2  7  V  +|(C sin(e ) + 5  p  8  C sin(e )) 6  + | ( C s i n ( £ ) + C sin(£ ))  2  2  p  3  2  / (4.28)  where KO is the Fourier transform of the KO operator described above. Thus, in Fourier space, we have the equation: 4>  n+l  -(,-*)-(,,*)(#)"  and therefore an amplification matrix given by (4.29) Let A be an eigenvalue of G and A an eigenvalue of M. Then we have: Gv  =S> ( J - M ) ~ (l + M)v  =  Av  = Av  Chapter 4. Scheme Analysis • » (i + M^v  =  ( l + A^Mv  =  > Mv  49  (i-M^Av (A-l)u  =  => A  =  1 1+ A 1+A 1- A  (4.30)  Physically, we don't expect the "norm" of the matter field,i.e. the amount of matter i n the spacetime, to increase with time; heuristically, conservation of energy prohibits it. Applying the Von Neumann necessary stability condition, which in this case demands that the spectral radius of G be smaller or equal to 1, we obtain:  1+A  2  1— Aj  < 1  => Re(A) < 0 Hence, the V o n Neumann stability condition reduces to demanding that the real part of the eigenvalues of M be smaller than 0. Although the form taken by the eigenvalues is quite complicated, it can be calculated easily by the quadratic root formula. We define the following quantities: A  =  a  =  b =  C - (sin (|)+sin (|)) 4  8  / (C 2  l  4  e  2 8  + 4C C )+4C C (sin (|)+sin (|)) 2  1  7  1  2  4  4Ci/j(C sin(£ ) + C sin(£ )) 5  p  6  z  W i t h these quantities now defined, the Von Neumann condition reduces to :  la + Va + b 2  A<  2  2h  2  We notice right away that k, the time step, is unrestricted! In practice, however, due to the way R N P L solves the Crank-Nicholson equations (point-wise NewtonGauss-Seidel iteration) there actually is a restriction on k. 4.2.2  Limiting cases  The above Von Neumann condition does not tell us much in its current form, which is a direct consequence of the fact that our equations are so complicated. However, we can consider some limiting cases:  Chapter  4. Scheme Analysis  50  1) low frequency modes (4> = £ = 0) z  A  ->  C  a  ->  /i (C|+4CiC )  b  ->  0  s 2  7  Recalling the definitions of the coefficients, we have that C1C7 = — m a and C = r>(/3") +Z> (y3 ). Typically, /? <C 1 and m a = 0(1), which implies that a < 0. Hence, the condition turns into something of the form Cs < 0; this is not strictly true everywhere, but we expect Cs to be many orders of magnitude smaller than 1, so it is hoped that this slight violation will not have dire consequences. 2  2  s  p  2  2  2  z  2) high frequency modes (£ = £ = ir — 29,9 < 1) p  z  Let us set Cs = 0 for simplicity. We then have, expanding to first order in 9: A  ->  -2e  a  -¥  4h C C  b  ->  8hC {C  +  2  l  7  1  f)  8C C (l-9 ) 2  l  4  + C )9 6  The result will depend on the sign of a. We define Ti  =  4/i CiC +8C C  T  2  =  8C C =a ^--  T  3  =  8/iCi(C +C )  2  7  1  = a ( - 4 / i m + -^) 2  4  2  2  2  1  4  5  =8h-^(D (aprp )+D (apip )) 2  6  p  2  z  W i t h these definitions, the Von Neumann condition becomes:  ' ^ ( ^  ^  '  +  | •|>-( ^  1 +  K I  )  ^  ' '  +  w •  ,  )  ,  Thus  Hence, we see that if T i < 0, the scheme obeys the Von Neumann condition (i.e. is dissipative) for high modes for e > 0. On the other hand, if T\ > 0, the scheme  Chapter 4. Scheme Analysis  51  obeys the condition only if e > F /2h. For typical values a = rjj = h = m = l, this yields e > 2. 1  Of course, what has been done cannot be taken at face value because we are not dealing with a constant coefficient problem. Hence, Fourier analysis is not "rigorously" justified, but experience of others has shown that it is still a good "ball park" tool. As it turns out, the results obtained are somewhat pessimistic: the scheme seems stable even for small values of e.  Chapter 5. Numerical Testing  52  Chapter 5  Numerical Testing 5.1 5.1.1  Solution Error and Tolerance Definitions  Consider the differential system d <p=£4>  (5.1)  t  and its discretized version D 4> = L4>  (5.2)  t  where <j) is the exact discrete solution. The solution error is given by e = 4>-4>  (5.3)  and, of course, is a measure of how "good" the discrete solution is compared to the exact solution. Since the problem at hand is essentially non-linear, the numerical solutions are obtained by an iterative procedure. Thus, these solutions are never exactly equal to the discrete solutions, but would converge to them in the limit of infinite number- of iterative steps. The residual is defined as the amount by which a given estimate of the discrete solution fails to satisfy the discrete equations, ie by Residual — —D <t> t  num  + L<j>  num  (5.4)  where 4> is the numerical approximation to the discrete solution. The difference between the discrete solution and the numerical approximation is not readily available. Here, we assume that the residuals and errors are of the same order; this assumption seems to lead to good results in this case. The numerical solution is obtained by iteratively driving the residual below a certain tolerance. Thus, the difference between the discrete solution and the numerical approximation is : num  4> — 4>num — tol  (5.5)  It is of great importance to know what the tolerance should be set to since it determines how much effort it takes to obtain the approximation. Manipulating the above expressions, one obtains: <t> ~ 4>num — e + tol  (5.6)  Chapter 5., Numerical Testing  53  It is clear that setting the tolerance much smaller than the solution error leads to wasted effort, since the numerical approximation, relative to the continuum solution, has intrinsic error e.  5.1.2  Solution Error Estimates  Since in general exact solutions are not known, the value of the solution error cannot be calculated exactly. However, it can be approximated by using different levels of discretization. Let h be the mesh size of level of ultimate interest. Then e = <f> - j> h  (5.7)  h  setting the tolerance to a very small number, the following estimates hold:  0~4«m,  4> -4>num h  '  (5-8)  leading to:  e ~ Kum - 4>Lm h  (5-9)  To be on the safe side, the tolerance should then be set to one or two orders of magnitude smaller than the truncation error estimate.  5.1.3  Values of the Tolerance  The above arguments are easily generalized to the problem at hand. To find the value of the tolerance at each level for a specific problem, one solves the equations using a very small tolerance (say 1 0 ) , at the desired level and at one level higher. Subtracting the two solutions, one obtains an approximation to the truncation error for every grid function and the tolerance should be set to one or two orders of magnitude smaller than this approximation. In effect, since the scheme used is of second order, the tolerance at higher levels can simply be deduced from the tolerance at lower levels by dividing the low level tolerance by four for every level between the low level and the high level. - 1 0  5.2  Convergence and Independent Residuals  The scheme used is of second order in all gridfunctions. In order to be confident the features observed in the numerical results are indeed features of the exact solution, it must be established that the code indeed present this second order convergence behavior. Furthermore, it must be established that the right equations are being solved by using an independent residual evaluation.  Chapter 5. Numerical Testing  5.2.1  54  Convergence  The convergence factor can be used to measure how well a code is converging. For a second order scheme, we define the convergence factor, Cf via  ~ \\P -f \\  C f  h  h  ( 5 - 1 0 ) 2  where ||-|| is the £2 ( R M S ) norm. In the limit ft-tOwe should find Cf ->• 4 for a second order accurate scheme. To test the convergence properties of the code, the following initial data was evolved 2  0t  =  A e x p [ - ^ t ^ ]  (5.11)  4,2 =  0  (5.12)  01  =  0  (5.13)  02  =  2pAexp[- ^  f 2  -]  (5-14)  with parameters A = 0.1, 6 = 8, level € {5,6,7,8}, T = 5.0 and domain 0 < p < 32, - 3 2 < z < 32, leading to h e {1.0,0.5,0.25,0.125}. The Courant factor was set to 0.1 implying that 50 time steps are taken at the coarsest resolution. As can be seen in Figure (5.1), the initial data specified leads to large general relativistic effects, as evidenced by values of the lapse that are substantially smaller than 1. Figures (5.2), (5.3) and (5.4) display the convergence factor as a function of time for level 5,6,7 and then level 6,7,8. Considering the resolution limitations of the code, the results seem to indicate second order convergence. The chosen data is generic; one can be reasonably sure that the fact that the code converges for this initial data implies that it converges for essentially arbitrary initial data. max  Chapter 5. Numerical Testing  55  Figure 5.1: The central value of the lapse as a function of time for the convergence testing trial data. The geometry is significantly curved.  Chapter 5. Numerical Testing  56  Figure 5.2: Convergence factor for trial initial data at level 5,6,7 for all the grid functions  Chapter 5. Numerical Testing  0  1  2  3 time  4  5  0  57  1  2 time  3  4  Figure 5.3: Convergence factor for trial initial data at level 6,7,8 for all the grid functions  5.2.2  independent Residuals  In order to ensure that the right equations are being solved, an independent residual evaluation must be performed, as suggested by Choptuik [22]. This consists i n taking a solution and calculating its residuals with respect to a distinct discretization of the original continuum equations. These residuals, if indeed the solution solves the right set of equations, should behave like the truncation error, i.e. in this case should converge to zero to second order (assuming that the independent discretization is also second order). The specific independent discretization used involves second order finite difference operators that are onesided. The domain is split in four quadrants and the finite difference operators are defined i n every quadrant such that they always point inward, towards the center of the domain. Results relevant to the trial initial data are presented in Figures (5.5), (5.6) and (5.7).  5  Chapter 5. Numerical Testing  58  3.8 3.  6  I  0  1  1  1  2  L _  3  :  1  4  5  time  Figure 5.4: The convergence factor of da/dp is shown here. A s can be seen, the convergence is much better than the convergence of a itself. This is possibly due to the fact that the boundary conditions on the lapse are only correct asymptotically, in the limit that the outer boundary of the computational domain tends to infinity.  Chapter 5. Numerical Testing  59  (  time  time  Figure 5.5: Convergence factor calculated at level 5,6,7 for all the independent residuals for the trial initial data. Some of the 1/p terms are badly resolved on the axis at level 5; this is probably why such bad convergence is observed in certain grid functions.  Chapter 5. Numerical Testing  0  60  1  2  3  4  5  3  4  5  time  0  1  3  2 time  4  0  1  2 time  Figure 5.6: Convergence factor calculated at level 6,7,8 for all the independent residuals for the trial initial data.  Chapter 5. Numerical Testing  61  L2 norm of the independent residuals of a  L2  L2 norm of the independent residuals of ip  L2  Figure 5.7: li norm of the independent residuals for a and ip. Clearly, these residuals are converging to zero. The other grid functions behave similarly.  5.3  Stability of Boson Stars  Static boson stars have been studied extensively i n the past and their stability properties are well understood i n spherical symmetry, as discussed by Jetzer [18]. The boson star solutions form a one parameter family where (po, the central value of the magnitude of the scalar field, is a convenient parameter. Although it appears as though rh, the Klein-Gordon mass parameter, is also a freely specifiable quantity, it only in fact defines the mass scale. Indeed, it can be easily checked that the system described by equations (2.174) to (2.175) is invariant under the change: m  Am,  w -> Xiv,  r r -¥ —  (5.15)  Chapter 5. Numerical Testing  62  Figure 5.8: Universal quantities as functions of cf> , as obtained from the initial data code. The mass parameter attains a maximum of 0.633 at the value 0 jt ~ 0.27. The radius of the star is defined as the radius at which the amplitude of the scalar field is 1% of its central value. In spherical symmetry, stars with cf>o smaller than </> j are stable and stars with 0 bigger than the critical parameter are unstable. 0  cr  cr t  O  Chapter 5. Numerical Testing  63  for A a dimensionless parameter. Hence two solutions differing only by the value of the m parameter are in fact the same solution, but expressed in different units. Thus, universal (independent of the value of m) dimensionless quantities describing the boson stars can be constructed: m^M_ MM' pl  m R M L  M m  star  pl  pl  pl  '  pl  l  p  l  W  i &  -  i b j  Here M is the A D M mass of the star, R tar is its radius and w is its frequency eigenvalue. The rather counter intuitive fact that Mm is independent of m is derived by using the geometric definition of M. Using linear perturbation theory, it is possible to predict the behavior of the spherical stars under small radial perturbations. Since the perturbation equations are linear in the perturbed fields, it is possible to Fourier decompose the perturbations into oscillatory "modes" with time dependence exp(—ryi), with 7 the frequency of the mode. Demanding that such perturbations conserve the Noether charge and vanish at infinity leads to an eigenvalue problem for 7. For 7 > 0, the star is stable with an oscillatory perturbation about the boson star solution, while for 7 < 0, the star is unstable. In the latter case, perturbations increase exponentially and the system leaves the linear perturbation regime in a finite time. Boson stars are known to be stable under radial perturbation as long as the parameter is lower than a certain critical value, ^> rit- This critical value corresponds to the maximum of the A D M mass, M x s  2  2  C  m a  5.3.1  Time Evolution of Boson Stars  When studied numerically, the boson star solutions are always in a slightly perturbed state. Indeed, the solution error appearing in the generation of the initial data, in conjunction with the fact that this data is obtained by solving a system of ordinary differential equations to a much higher degree of accuracy than could be obtained via a corresponding solution in two dimensions, as well as further perturbations generated by the time evolution of this data are bound to take the boson star configuration away from the exact analytical solution. Thus the numerical time evolution of a boson star never yields a precisely static solution. A first example is considered in order to establish how the numerics should be set up. The first boson star studied has parameter 4>o = 0 . 1 , well within the stability region with respect to spherical perturbations. From Figure (5.8),. it is seen that the radius of such a star is about 20 Planck lengths for m = M ; . A preliminary run is launched to determine what value the tolerance should have. The parameters chosen are: p  e = 0.001,  A = 0.1,  i?  m a x  =  tolerance =  32, 10"  /euef = 6,7,8 1 0  T  m a x  = 25.0 (5.17)  Chapter 5. Numerical Testing  W/ -/*)(xio- ) h  5  a  TP  P"  P  28  13  2  2  0i 37  z  64 02 10  (71  0-2  59  77  Table 5.1: Approximation to the truncation error for all grid functions at level 6,7. The value is taken at T=5.0,'when the effects due to the use of the same initial data have disappeared. The tolerance for level 6 will thus be set to 2.0 x 10~ , which is one order of magnitude smaller than the smallest error above. 6  where A is the Courant factor. The approximation to the truncation error can be found in table (5.1). The tolerance is thus set to: .  level tolerance  6 2.0 x I O "  6  7 5.0 x 1 0 -  7  8 1.0 x I O "  7  (5.18)  Next, knowing what the tolerance should be set to, the initial data is evolved for a longer time, with A = 0.5. As can be seen in Figure (5.9), the solution appears to change behaviour around t = 300. Charge and mass conservation begin to become seriously violated and the scalar field starts behaving in a way that is inconsistent with a stable perturbed boson star. The scalar field reaches an amplitude of 1% of its central value on the outer boundaries at t ~ 400. The boundary conditions were derived under the assumption that the scalar field would vanish at the outer boundaries and thus the solutions become untrustworthy when the scalar field amplitude becomes large there. Thus, to see whether the star is really unstable or if the features observed are only boundary effects, the same initial data is evolved, but using R = 64. As can be seen on Figure (5.10), the bad behavior observed at i Z = 32 is a boundary effect. As a general rule of thumb, the outer boundary should be set three boson star radii away from the origin. T h i s seems to be the best compromise between distant outer boundaries and good resolution, since the computational resources available are limited. max  m a x  Chapter 5. Numerical Testing  65  0.566  500  50  100 150 200 250 300 350 400 450 500 time  Figure 5.9: Norm of the scalar field, Noether charge and A D M mass as functions of coordinate time for </> = 0.1, R = 32. A s can be seen in the Figure, results appear to deteriorate around t = 300, where the scalar field attains about 1% of the central value on the outer boundaries at the coarsest level. 0  max  Chapter 5. Numerical Testing  66  0.545 MADM 0.544  |$(0)|  2  0.012 0.011 0.01 0  100  200  300  400  500  600  700  time  Figure 5.10: Norm of the scalar field, Noether charge and A D M mass as functions of coordinate time for </> = 0.1, R ax = 64. The figures show that the star is indeed stable and that the outer boundaries need to be placed further away than Rm = 32. 0  m  ax  5.4  Spherically Perturbed Boson Stars  As a non trivial test of the numerical implementation, spherically perturbed boson stars are studied. In particular, i n a recent article Hawley and Choptuik [6] quote the frequencies of the fundamental perturbation mode as calculated in linear perturbation theory. Their values must be modified slightly i n order to match the system of units used in this thesis. Their perturbative results are compared to those obtained from the current numerical implementation in table 5.2 . The perturbation used for <p = {0.06,0.10,0.14} is given by: 0  1*1 =  4> -+<!> +8 </>  8<f>  A exp [  p +z 2  2  A  2  (5.19)  with A = 0.02 and A = 2, leading to a spherical perturbation that is substantially larger than the truncation error, thus allowing convergence testing. The  Chapter 5. Numerical Testing  4>o 0.06 0.10 0.14 0.18 0.22  0J 0.922 0.873 0.826 0.781 0.738 O  O-pert  (Xl0~ ) 4  154 ± 3 225 ± 2 275 ± 1 293 ± 1 267 ± 1  O-num ( X l 0 ~ 161 ± 17 242 ± 25 273 ± 26 279 ± 20 ???  67  4  )  Table 5.2: Comparison between the frequencies obtained from the code and the values obtained by Hawley and Choptuik by perturbation theory. o~p t is obtained by taking the square root of the values quoted in the article and multiplying by ao as obtained from the initial data code. The uncertainties are approximate. Values quoted come from level 7 data. The outer radii used for the domain are: Rmax = 64 for 4> = 0.06,0.1 and R = 48 for 0 =0.14 and 0 = 0.18. The lack of resolution, due to numerical cost constraints, render the code unable to produce valid solutions at (p = 0.22, which is too close to the unstable branch. er  0  m a x  O  O  0  central values of different grid functions as a function of time and the convergence properties of the solutions can be viewed in Appendix A . The frequencies are calculated from the graphs of |3>(0)| vs. time; by measuring the position of successive mimimae (within uncertainty), one can deduce a value for the period of oscillation and thus the frequency. This is possible because the amplitude of the oscillation is dominated by the amplitude of the lowest mode, the fundamental mode we want to consider. Unfortunately, the resolution used is inadequate to properly study perturbed solutions close to the boson star unstable branch, i.e. solutions with (po close to 0.27. As it turns out, the perturbed 4>o = 0.14 solution is unstable at level 6 and the (po = 0.18 solution is unstable at level 7. Stars with a (po value higher than 0.14 are thus considered without additional perturbation. The (po = 0.18 star behaves well at level 7 under truncation error perturbation. The program crashes for a (po = 0.22 star, on the other hand; the only way to keep evolution going for a long time is to add so much dissipation that the star "evaporates" into a lighter configuration. To ensure that this indeed is caused by lack of resolution, the equivalent situation i n spherical symmetry is considered. Using a one dimensional spherical code, the solutions can be studied with much higher resolution. The spherical code is constructed i n the exact same way that the axisymmetric code was built, with the exceptions that 1) the code is 1 dimensional, and 2) the boundary conditions are used directly to solve for the boundary points, instead of indirectly through the use of ghost points. A s can be seen in Figure (5.11), spherical calculations lead to what appears to be an unstable star at level 7, but turns out to be stable at a finer resolution. It is thus not surprising that the two dimensional code should exhibit the same properties. 2  Chapter 5. Numerical Testing  i  1  1  68  :—r  1  0  50  100  150  200  250  300  350  400  450  500  time  Figure 5.11: Amplitude of the scalar field for a truncation error perturbed boson star with parameter (p = 0.18 using the spherically symmetric code briefly described in the text. Computations at levels 7 and 10 and shown. Here, R x = 32 which is the same value that leads to an unstable evolution at level 7 with the 2-D code. A s can be seen on the plot, the level 7 spherical calculation also appears unstable whereas the level 10 data shows that the star is i n fact stable. 0  ma  Chapter 5. Numerical Testing  5.5  69  Newtonian boson stars collision  In the limit of light stars moving very slowly, the conformally flat approximation, like the full Einstein equations, reduce to Newtonian gravity. A non trivial test for the code is to see if it will generate collision between stars in the Newtonian limit in a time comparable to the timescale predicted by Newtonian gravity (assuming that the stars are point masses). The Newtonian equation for the gravitational force:  can easily be integrated to yield the time two point masses should take to collide. Choosing mi = m = M and do as the initial distance between the two stars, the time to collision i n Planck units is given by: 2  star  (5.21) Results from a sample calculation are presented in Figure (5.12), where two light boson stars (<fo = 0.02, M DM = 0.28) are initially placed 120 Planck units apart. The Newtonian prediction for the collision time is given by t i ~ 1950 and the level 6 result show a collision time of 2310. Considering that the situation is not purely Newtonian and that the stars are not point masses, this seems to indicate that the solution reduces to the Newtonian case, as it should. A  co  Chapter 5. Numerical Testing  70  t=0.00  am  0.00e+00  2.00e-02  1=1800.00  6.76e-07  2.03e-02  t=2310.00  9.35e-07  8.07e-02  Figure 5.12: Newtonian collision between two light boson stars. The value of 0 is 0.02 and the stars are initially 120 Planck units apart. The mass of each star is 0.28 Planck units. The central value of |#j reaches a maximum at t ~ 2310 Planck units. O  Chapter 6. Critical Phenomena  71  Chapter 6  Critical Phenomena 6.1  Introduction  The subject of critical phenomena has been extensively reviewed by Gundlach [13]; here we present only a brief overview of some key results. Consider an initial data set for time evolution of a general relativistic system characterized by a parameter p such that for large enough values of the parameter, the spacetime contains a black hole and that for small enough values of p all the matter dissipates to infinity. Then there must be a critical value of the parameter, p*, at the threshold of black hole formation. The question of what the critical solution, i.e. the spacetime lying on the edge between dissipation and collapse, looks like was first studied by Choptuik [20]. The critical solution is defined as the limit solution for the spacetime when p —> p*, and subsequent work suggested that it has only a single growing mode in perturbation theory. According to the sign of the coefficient in front of that growing mode, the near critical spacetime will either collapse or dissipate all its matter to infinity. Choptuik found unexpected and interesting behavior for diverse families of initial data leading to spherically symmetric collapse. Depending on the matter model used, two types of symmetry have been observed for solutions close to criticality: the critical solution can be static or periodic (type I solution) or the critical solution can have scale-translation symmetry, where the solution is either continuously or discretely self-similar (type II solution). Working with a massless scalar field, Choptuik [20] found type II solutions. B y tuning the parameter p closer and closer to p*, it was found that the solution exhibited "discrete self-similarity" (DSS), i.e.: Z,(r,t) =  Z (e r,e t) nA  nA  t  n = {1,2,3,...},  (6.1)  A ~ 3.44  where Z* is any dimensionless dynamical variable computed at criticality, and where the time coordinate t must be defined appropriately. It was found that the masses of the black holes formed for p close to p* scale as  M~c(p- *y P  (6.2)  The value found by Choptuik for 7 was 0.37. Other models were also investigated and it was found that these initial data sets led to continuous self-similar (CSS) critical solutions: (6.3)  Chapter 6. Critical Phenomena  72  These types of critical solutions were called "Type II DSS" and "Type II CSS". Other matter models exhibit a different behavior, dubbed "Type I". For these solutions, for p > p* black hole formation is switched on with a finite mass. The scaling appears in the lifetime of the solution, i. e. solutions close to criticality behave like the critical solution for a time: AT  ~  —7 In \p — p* I + const.  (6.4)  The geometric functions for the critical solution are either periodic or static rather than self similar. In particular, this is the type of behavior observed when a spherical boson star is appropriately perturbed by a spherical pulse of real scalar field [6]. The non-spherical case has received little attention in the past; Abrahams and Evans [24] were the first to study the collapse of gravitational waves in vacuum in axisymmetry. In this thesis, a single spherical boson star perturbed by an incoming axisymmetric wave of massless real scalar field is used to investigate critical phenomena in a non-spherical case within the context of the conformally flat approximation.  6.2  Numerical Setup  Many technical factors influence which one-parameter family of initial data will be studied. The competing demands of highest possible resolution for lowest possible numerical costs greatly constrain the kind of situations practically available for study. The situation chosen is given by a single spherical boson star that is perturbed by an incoming pulse of massless real scalar field. The massless scalar field, described by </> andCT ,obeys the exact same equations as the components of the complex scalar field with m = 0 (see equations (2.137) and (2.138)). The code is thus naturally extended to deal with these fields; the matter source terms for the geometry are modified by introducing the new real field and the evolution code is modified as well. Since more dynamics at the boundaries is expected, all the matter fields are endowed with outgoing boundary conditions at the outer boundary. This amounts to sssuming that the field will behave approximately like a flatspace outgoing spherical wave when it reaches the outer boundary. The form of the perturbation at the initial time is given by: 3  ^ *a  =  Aexp[-^-]  Wex [-4^----«] P  =  3  a  op  \z\<z  0  \z\>z  "  (6 0  5)  (6.6)  The initial data for the scalar field basically represents a cylindrical gaussian pulse of finite extent in the z direction. This finite extent is chosen such that the bulk of the pulse is initially contained in the domain considered (see figure  Chapter 6. Critical Phenomena  73  t=0.00  4.79e-10  2.70e-01  Figure 6.1: Typical form of the initial data. The norm of the complex scalar field and the real massless scalar field are superimposed in the figure.  74  Chapter 6. Critical Phenomena  (6.1)). The initial data for the conjugate variable a% is chosen to make the pulse initially ingoing. The initial data describing the boson star is taken to be the single boson star solution. Care must be taken not to specify too large a perturbation, or else p + S as described by equation (2.84), can become negative, leading to a non-physical a > 1. The parameters chosen are: l  m  <j> = 0.27,  i:  R =40, <5 = 8, p = 32, z = 35 6 = 10, level = 5,6, A = 0.1, e = 0.25  0  max  0  0  (6.7)  and A is the tunable parameter.  6.3  Results  A , the family parameter, was tuned to high precision at level 5 and 6- In both cases, near critical solutions appear exhibit type I behavior. The critical values of the parameters obtained are: A*  =  0.3185799460534935 ± 5 x 1 0 ~  A%  =  0.2117047263 ± 1 0 "  b  16  (6.8) (6.9)  1 0  Solutions close to criticality all exhibit the same behavior: initially, the geometry is non-spherical due to the incoming pulse but it "shakes off" the non-spherical perturbation and the solution becomes periodic for a time before either dispersing most of the matter to infinity, or collapsing to a black hole, depending on the sign of A — A*. The nature of the distribution of the matter left at the origin by a subcritical evolution is unclear, but a very large fraction of the initial matter is dispersed to infinity. Evidence of the spherical nature of the critical solution can be viewed in Appendix B . The black holes formed by supercritical solutions seem to have a finite mass and the value of the scaling exponent, as defined by equation (6.4) is given by 9.7 at level 5 and 3.7 at level 6. The value cited by Hawley and Choptuik [6] for spherical simulations is 9.2. The large discrepancy at level 6 motivates further investigation. The relatively naive fitting method used for level 5 data probably is not appropriate at level 6 because the critical solution was not tuned to the same precision. To obtain a measure of just how different the level 6 result is from the level 5 result, two techniques were used: 1) a 2-parameter fit was used on level 6 data, fitting 7 and A*, to the experimental points dropping the last one, and 2) a similar analysis was performed on a subset of the level 5 data that mimics the gathered data at level 6. Using the gnuplot data fitting package, the 2-parameter fit on the subset of the level 6 data excluding the final point yielded: A*  6  =  0.211705 ± 4 x 1 0  7e  =  3.7 ± 0.3  - 5  Chapter 6. Critical Phenomena  75  This is similar to the first naive one-parameter fit with the value of A$ stated earlier. Using the same algorithm on the full data set at level 5 yielded: A*  5  7  / u 5  "  =  0.31858  =  9.4  ±  1 x  1 0 ~  5  ± 0 . 3  This is slightly different from the value obtained previously for 7 . Using only a truncated set approximating the level 6 set, the values obtained were: Al -ytrunc  =  0.31858  =  8.4  ±  5 x  1 0 "  5  ± 0 . 8  As can be seen from these results, the level 5 truncated set does yield a result that is substantially different than the full set value, but the value is still very different from what is observed at level 6 . We can also compute another scaling exponent, 7, which for the case of a spherical static critical solution should be related to 7 via 7 = 7/ao where ao is the central value of the lapse during critical evolution. Specifically 7 is calculated via the same basic technique used above in the computation of 7 , but using coordinate time rather than central proper time. As can be seen in figure 6 . 3 , the difference between the estimates of 7 for level 5 and level 6 calculations appears to be of the order of 5 0 % instead of the roughly 2 0 0 % the previous discussion seems to imply. This would indicate that the method employed to compute 7 is somewhat imprecise and that the discrepancy between the results at the two resolutions is not that large. Furthermore, the critical solutions at the two discretization scales appear to be somewhat distinct (due, presumably, to the substantial discretization errors on the relatively coarse meshes used), which may well account for the remaining discrepancy. Indeed, the minimum values of the lapse function in the critical regime are significantly different: at level 5 , the mimimum value of the lapse is a 5 | i = 0 . 1 7 9 and at level 6 it is a 6 | i n = 0 . 0 9 5 . This indicates that the solutions considered are indeed different, which no doubt contributes in part to the different values of 7 that were computed. m  n  m  Although the obtained values are of the same order of magnitude as the value cited by Hawley and Choptuik, it remains unclear whether the observed data in axisymmetry behaves i n the same way as the spherically symmetric case. If higher resolutions could be achieved it would be interesting to see if the axisymmetric value of 7 converges to the spherical one.  76  Chapter 6. Critical Phenomena  =6.0 x 10~  6  k u l J J J yy j Jyy y y j yy j y jyyJ jjyJoillilLi = 1.3 x I O "  1 1  \A - A*\ = 5.0 x I O "  1 6  w  0  500  1000  1500  Figure 6.2: Central value of a for diverse values of the family parameter for supercritical level 5 data. A s can be seen in the figure, the closer the parameter is to its critical value, the longer the solution behaves like the critical one.  2000  Chapter 6.  VdUu J u VdUJJU j u JJ V J U U U JJu  0  i 200  400  77  Critical Phenomena  \A -A*\  = 1.5 x 10~  6  \A-A*\  =8.1  9  \A-A*\  600  x 10"  = 1.1 x  10-  10  800  Figure 6.3: Central value of a for diverse values of the family parameter for supercritical level 6 data.  1000  Chapter 6. Critical Phenomena  -30  -25 ln\A-A*\  -20  78  15  -10  Figure 6.4: Lifetime spent in the critical regime as function of the family parameter, A, for level 5 data. Proper time at the center of the grid is used; the time variable must be rescaled using r = J adt (the shift is negligible). The value of the critical parameter, defined in equation (6.4), is given by 7 = 9.7, where the quoted uncertainty comes from the uncertainty on A*. The value of 7 is obtained by making a linear fit to the data.  79  Chapter 6. Critical Phenomena  -20  -19  -18 ln\A-A*\  -17  16  -15  -14  -13  Figure 6.5: Lifetime spent in the critical regime as function of the family parameter, A, for level 6 data. Proper time at the center of the grid is used; the time variable must be changed using r = j adt (the shift is negligible). The value of the critical parameter, defined i n equation (6.4), is given by 7 = 3.7, where the quoted uncertainty comes from the uncertainty on A*. The value of 7 is obtained by making a linear fit to the data.  80  Chapter 6. Critical Phenomena  At  1150 1100 1050 1000 950 900 850  —O  i  I  i  I  i  i  75 = 35 ± 6  i  O i  -23  i -22  1  650  i -21  I  i -20  -19 ln|yl -A*\ I  -18  i -17  o  I  l  i -16  I  76 = 23 ± 3  O  600 At  i  550 O  500 450  i  -23  -22  i  i  i  -21  -20  -19  i -18  i  "1>  -17  -16  -  ln\A-A*\  Figure 6.6: Lifetime spent in the critical regime over comparable ranges of the family parameter, A, for level 5 data and level 6 data using coordinate time. The values of 7 are obtained by making linear fits to the data. The disagreement between level 5 and level 6 is much smaller here than for the computations of 7 that were made using central proper time.  Chapter 6. Critical Phenomena  81  Figure 6.7: A D M mass for supercritical and subcritical solutions for level 5 data. Once an apparent horizon forms, the value of the A D M mass rapidly increases, probably indicating that the scheme quickly becomes less accurate due to the steep gradients arising from the impending coordinate pathology. When the solution oscillates in the critical regime, though, the bulk of the matter is concentrated around the origin and the computed mass yields a sensible result. Clearly, the solution does not "shed mass" while in the critical regime (type II behavior) and the black holes formed must have a finite mass, which is typical of type I behavior.  Chapter 7. Conclusion  82  Chapter 7  Conclusion In this thesis, boson stars studied in axisymmetry under the conformally flat approximation have been shown to behave similarly to the spherical solutions of the Einstein-Klein-Gordon equations under small perturbation. The program devised to solve the set of equations has been shown to converge well, to conserve the A D M mass and the Noether charge and independent residual evaluations have confirmed that indeed the right equations were being solved. Additional results, though preliminary, indicate that the model also admits Type I critical behavior at the black hole threshold when a boson star is per- turbed by a massless scalar field pulse. Results indicate that within this approximation the critical solution seem to be spherical, as evidenced in appendix B , even though the perturbations are highly non spherical. Furthermore, as p —> p*, the solutions in the critical regime exhibited periodicity that probably can be interpreted as the "ringing" of an unstable boson star. On the other hand, a definitive value of the scaling exponent 7 has not yet been obtained, since the current results display significant dependence on the mesh spacing used. A major drawback suffered by numerical simulations in axisymmetry is the high computational cost for adequately resolved calculations, particularly relative to spherically symmetric simulations. Although there is no "magic resolution" and solution features should be observable even at coarse discretizations, convergence and conservation testing are less crisp than they would be in a highly resolved spherical (one dimensional) numerical investigation. Furthermore, the desire for the highest resolution possible interferes with the demand for large numerical domains in order to accommodate the boundary conditions, which are only valid in the asymptotically flat region. Bigger domains mean coarser resolution or longer solving time, which can become problematic in a tuning problem. There still remains the difficult problem of quantifying the validity of the conformally flat approximation. This question is of high interest to numerical relativists, for solving the conformally flat system of equations is much simpler than tackling the full Einstein equations. In this thesis, this approximation seems to yield solutions consistent with the expected behavior of a full Einstein solution at a fraction of the numerical price, so a definitive exposition of its validity would be profitable to the numerical relativity community.  Bibliography  83  Bibliography [1] C. Misner, K . Thorne, J . Wheeler: "Gravitation", W . H . Freeman, New York, N Y , 1973 [2] R. Wald: " General Relativity", The University of Chicago Press, Chicago, IL, 1984 [3] J . York: "Kinematics and Dynamics of General Relativity", in Source of Gravitational Radiation, L . Smarr, ed. Cambridge University Press, Cambridge, 1979 [4] A . Brandt: "Guide to Multigrid Development", in Multigrid methods: proceedings of the conference held at Koln-Porz, Nov. 23-27 1981, Springer Verlag, Berlin, 1982 [5] W . Briggs: " A Multigrid Tutorial", S I A M , Philadelphia, P A , 2001 [6] S. Hawley and M . Choptuik: "Boson Stars Driven to the Brink of Black Hole Formation", Phys. Rev. D, 62, 104024, 2000 [7] J . Wilson, G . Mathews and P. Marronetti: "Relativistic Numerical Model for Close Neutron Star Binaries", Phys. Rev. D, 54, 1317-1331, 1996 [8] J . Wilson and G . Mathews: "Binary Induced Neutron Star Compression, Heating, and Collapse", Astrophys. J., 482, 929-941,1997 [9] G . Mathews, P. Marronetti and J . Wilson : "Relativistic Hydrodynamics in Close Binary Systems: Analysis of Neutron Star Collapse", Phys. Rev. D, 58, 043003, 1998 [10] G . Cook, S. Shapiro, S. Teukolsky: "Testing a Simplified Version of Einstein's Equations for Numerical Relativity", Phys. Rev. D, 53, 5533-5540, 1996 [11] M . Miller and Wai-Mo Suen: "Towards a Realistic Neutron Star Binary Inspiral", 2003, (unpublished) [gr-qc/0301112] [12] M . Choptuik, E . Hirschmann, L . Liebling and F . Pretorius: " A n Axisymmetric Gravitational Collapse Code", Class. Quant. Grav., 20, 1857-1878, 2003 [13] Carsten Gundlach: "Critical Phenomena in Gravitational Collapse", Phys. Rept, 376, 339-405, 2003  Bibliography  84  [14] U . Ascher and L . Petzold: "Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations", S I A M , Philadelphia, 1998 [15] U . Ascher, S. Ruuth, and R. Spiteri: "Implicit-Explicit RUnge-Kutta Methods for Time-Dependent Partial Differential Equations", Applied Numerical Mathematics, 25, 151-167, 1997 [16] U . Ascher, S. Ruuth, and B . Wetton: "Implicit-Explicit Methods for TimeDependent Partial Differential Equations", SIAM J. Numer. Anal, 32, 797-823, 1995 [17] Michael Peskin and Daniel Schroeder, " A n Introduction to Quantum Field Theory",Perseus Books, Cambridge, M A , 1995 [18] P . Jetzer: "Boson Stars",Phys. Rept, 220, 163-227, 1992 [19] F . Pretorius, "Numerical Simulations of Gravitational Collapse", P h . D . thesis, The University of British Columbia (unpublished), 2002. [20] M . W . Choptuik, "Scaling and Universality in Gravitational Collapse of a Massless Scalar Field" Phys. Rev. Lett, 70, 9-12, 1993. [21] M . W . Choptuik," Numerical Techniques for Radiative Problems i n General Relativity",Ph.D. thesis,The University of British Columbia (unpublished), 1986 [22] M . W . Choptuik, personal communication, 2002 [23] R . L . Marsa,"Radiative Problems i n Black Hole Spacetimes", P h . D . thesis, The University of Texas at Austin (unpublished), chapters 8 & 9, 1995 [24] A . M . Abrahams and C . R . Evans,Phys. Rev. Lett, 70, 2980 (1993)  Appendix A. Perturbed Boson  Stars Data  85  Appendix A  Perturbed Boson Stars Data Here is displayed the data relevant to the calculation of the frequency of the fundamental mode for a single spherical boson star under spherical perturbation as well as evidence that the code behaves as expected and converges properly. The static boson star solutions are labeled by two parameters: 0o, the maximum value of the scalar field's amplitude at the origin, and m , the mass parameter in the Klein-Gordon equation, m is always set to 1 and 0o is chosen between 0.06 and 0.18, which always puts the static boson star solution in the stable region of parameter space (see figure 5.8). Two types of perturbation are considered. Truncation error perturbation, which is the perturbation imposed on the boson star solution by the fact that the continuum solution doesn't satisfy the discretized equations exactly, drives the solution away from staticity. This perturbation is always present and is the only form of perturbation acting when no other perturbation is introduced "by hand". A second type of perturbation is also considered. The initial data is modified by adding a spherical gaussian pulse to the initial value of the scalar field. This perturbation is parametrized as follows:  4>i ->  0i +  dcp  (A.l)  ->  02  (A.2)  0"1  -»  Ol  (A.3)  0  ->•  <7  02  2  +  2  pH  (A.4) (A.5)  where $(* = 0) = 0 i + i<f> is the initial value of the scalar field and U(t = 0) = o i + io is the momentum conjugate field to The particular form chosen respects the demand that the conjugate momentum vanish on the axis. The perturbation is given by: 2  2  2  2  50 = .4 exp [ - ^ ^ - ]  (A.6)  Hence the perturbation is described by two parameters: A and A, with the understanding that A = 0 represents truncation error perturbation. Finally, the numerical solving process also involves freely adjustable parameters: . R a x , which defines the position of the outer boundary, e, the KreissOliger dissipation coefficient, h and k, the meshsizes in space and time, and m  Appendix A. Perturbed Boson Stars Data  86  0.0065 0.006 0.0055 h 0.005 |$(0)|  2  0.0045 0.004 0.0035 0.003 100  200  300  400  500  600  700  800  900 1000  time  Figure A . l : Amplitude of the scalar field as a function of time. After t ~ 150, • boundary effects generate high frequency noise. Level 8 data was only generated up to t=250 due to high computational costs. Parameters: (f> = 0.06, A = 0.02, i ? = 64 0  m a x  level, which measures the coarseness of the mesh. A s discussed in Chapter 3, all these parameters are not independent:  ^ A  +1  =  =  (  \  ^  (A.8)  where A, the Courant factor, has been introduced. Thus, a sufficient set of parameters to describe a particular implementation is composed of level, R x, A and e. Thus, i n summary, a numerical solution will be labeled by the following parameters: ma  0o,  A,  A = 2,  level = 6,7,8,  R , max  A = 0.25,  e = 0.01  (A.9)  where the numerical values indicate that these parameters are the same for all the graphs presented in this appendix.  Appendix A. Perturbed Boson Stars Data  a  87  tp  time  Figure A . 2 : Convergence factor for some of the grid functions at level 6,7,8 for 50 < t < 150. Parameters: 4> = 0.06, A = 0.02, R = 64 0  max  Appendix A. Perturbed Boson Stars Data  01  88  02  C  60 '  80  100 time  120  140  Figure A . 3 : Convergence factor for some of the grid functions at level 6,7,8 for 50 < t < 150. Parameters: 0 = 0.06, A = 0.02, # = 64 O  m a x  Appendix A. Perturbed Boson  Stars Data  89  Appendix A. Perturbed Boson Stars Data  0.467 0.4668 0.4666  90  level 6 level 7 level 8  0.4664 0.4662 0.466 MADM  0.4658 h 0.4656 0.4654 h 0.4652 0.465 b 0.4648 0  100  200  300  400  500  600  700  800  900 1000  time  Figure A . 5 : A D M mass at level 6,7,8. Boundary effects become important around t ~ 800 and the solution is probably untrustworthy after that time. The mass is only computed up to t = 250 at level 8 because of large computational cost. Parameters: <j> = 0.06, A = 0.02, i ? = 64 0  m a x  Appendix A. Perturbed Boson Stars Data  -0.921  n  1  1  1  r  91  level 6 level 7 -H level 8  -0.922 -0.923 Q -0.924 |-0.925 -0.926 -0.927  j  0  i  i  i  i  100 200 300 400 500 600 700 800 900 1000 time  Figure A . 6 : Noether charge at level 6,7,8. Boundary effects become important around t ~ 800 and the solution is probably untrustworthy after that time. The charge is only computed up to t = 250 at level 8 because of large computational cost. Parameters: cf> = 0.06, A = 0.02, R = 64 0  max  Appendix A. Perturbed Boson Stars Data  0.026 level 6 0.024 - level 7 level 8 0.022  i  I  1  0.02  92  1  /  \  l\~  $(0)|fJ.018 0.016 0.014  0.01  \.,\  /  0.012 "  0  /  1  1  i  i  200  400  600  800  \ 1000  time  Figure A.7: Amplitude of the scalar field as a function of time. Level 8 data was only generated up to t=150 due to high computational costs. Parameters: (f> = 0.10, A = 0.02, i ? = 64 0  m a x  Appendix A. Perturbed Boson Stars Data  93  Appendix A. Perturbed Boson Stars Data  94  Appendix A. Perturbed Boson Stars Data  95  0"!  60  80  100 time  120  140  Figure A.10: Convergence factor for some of the grid functions at level 6,7,8 for 50 < t < 150. Parameters: 0 = 0.10, A = 0.02, i i = 64 O  m a x  Appendix A. Perturbed Boson Stars Data  Figure A . l l : ADM mass at level 6,7,8. Parameters: <f> = 0.10, Q  A = 0.02,  R  max  = 64  96  Appendix A. Perturbed Boson Stars Data  Figure A.12: Noether charge at level 6,7,8. Parameters: (f> = 0.10, A = 0.02, 0  jR  m a x  = 64  97  Appendix A. Perturbed Boson Stars Data  0  200  400  600  98  800  1000  time  Figure A.13: Amplitude of the scalar field as a function of time. Only level 7 data is presented. A t lower resolution, the star is unstable and no oscillating solution could be produced. Parameters: </> = 0.14, A = 0.02, P ^ a x = 32 0  Appendix A. Perturbed Boson Stars Data  M  0  6  4  5  99  Q  1000  1000  Figure A.14: ADM mass and Noether charge at level 7. Boundary effects render the solution invalid around t=800. Parameters: <f> = 0.14, A = 0.02, i ? = 32 0  m a x  Appendix A. Perturbed Boson Stars Data,  100  Appendix A. Perturbed Boson Stars Data  101  Appendix B. Grid Function Slicing for Critical Phenomena  102  Appendix B  Grid Function Slicing for Critical Phenomena Here is presented some time evolution data for single static boson star perturbed by an axisymmetric massless scalar field pulse. The data shows evidence that, although the perturbation is highly non-spherical, the resulting close-to-critical solutions appear to be close to spherical. A l l the data presented is super-critical, i.e. A > A* and thus the solution eventually develops a black hole.  Appendix B. Grid Function Slicing for Critical Phenomena  ^  / / /  t=0.00  t=50.00  t=i^o^  t=200.00  103  V=ioo.oo  t=500.00 -40-30-20-10 0 10 20 30 40 Figure B . l : p — 0 (dotted line) and z = 0 (full line) spacelike slices of the lapse for \A — A* \ ~ 1 0 ' at level 5. A t initial times, the lapse is highly aspherical because of the perturbation, but once the solution is in the critical regime (t > 400) the radial and axial slices become close to identical. The critical solution's axial slice past t ~ 500 seems off centered by one grid point to the left with respect to the radial slice. This effect converges away, as will be seen at level 6. 1 5  Appendix B. Grid Function Slicing for Critical Phenomena  t=0.00  t=50.00  .•'^t=150^  ^t^20C3^  t=400.00  t=500.00  104  t=100.00  2 1.8  t=600.00  1.6 1.4 1.2 1 ' i i i i^i^i r -40-30-20-10 0 10 20 30 40 Figure B.2: p — 0 (dotted line) and z = 0 (full line) spacelike slices of the conformal factor for \ A — A* \ ~ 1 0 at level 5. The same features present in the lapse data can be observed. - 1 5  Appendix B. Grid Function Slicing for Critical Phenomena  ."^  t=0.00  t=50.00  t=150.00  t=200.00  A  0.16 0.14 0.12  t=400.00  i  i  i  \  i  i  i t=300.00  . t=500.00  i  105  t=600.00  • x  -40-30-20-10 0 10 20 30 40  |*(0)|  3.5e-07 3e-07 2.5e-07 2e-07  .-.A  :  1.5e-07 le-07 5e-08 0 -40  t=1282.50 i -30  i -20  i • -10  i V i 0  10  i 20  i 30  40  Figure B.3: p = 0 (dotted line) and z = 0 (full line) spacelike slices of the amplitude of the complex scalar field for \A — A*\ ~ 1 0 ~ at level 5. Even though the boson star parameter is <j>o = 0.27, the presence of the perturbation changes the amplitude of the boson star significantly. The star appears to be in a aspherical excited mode, as the nodes on the axis slice indicate. The star seems to be dissipated away before the end of the lifetime of the critical solution, but a look at the boson star data at a relative scale by reducing the range of the vertical axes indicates that massive scalar field is still present in the vicinity of the origin for the duration of the critical regime. Clearly, the star is not well resolved; it occupies roughly a subgrid of 7 by 15 points. 15  Appendix B. Grid Function Slicing for Critical Phenomena  t=50.00  1.5 1 0.5 0 (O) 0  ^  •  106  \ ^ t = 100.00  t=150.00  t=200.00  ^^=300^0(3  J^40(l.00  . V^soo.oo  t=550.00  3  -  I  I  I  I  I  I  I  -40-30-20-10 0 10 20 30 40 Figure B.4: p = 0 (dotted line) and z = 0 (full line) spacelike slices of the amplitude of the real scalar field for \A - A*\ ~ 1 0 at level 5. The data presented is super-critial, i.e. it leads to the formation of a black hole. - 1 5  Appendix B. Grid Function Slicing for Critical Phenomena  t=0.00  107  t=50.00  1  " , , , ' /  ,t=300 00  '  r  t=350.00  -40-30-20-10 0 10 20 30 40 Figure B.5: p = 0 (dotted line) and z = 0 (full line) spacelike slices of the lapse for \A - A*\ ~ 1 0 at level 6. The asymmetry observed at level 5 is now gone. - 9  Appendix B. Grid Function Slicing for Critical Phenomena  t=0.00  t=50.00  \  108  t=100.00  — — -  • • i • ' i '  i  -40-30-20-10  \  t=150.00  -^t=250^  \  t=300.00  t=400.00  i  i  i  i  0 10 20 30 40  Figure B.6: p = 0 (dotted line) and z = 0 (full line) spacelike slices of the conformal factor for \A — A*\ ~ I O at level 6. - 9  Appendix B. Grid Function Slicing for Critical Phenomena  109  |*(0)|  Figure B.7: p = 0 (dotted line) and z = 0 (full line) spacelike slices of the amplitude of the complex scalar field at a late time for — A* | ~ 10~ at level 6. The star is better resolved than at level 5, but is still clearly very poorly resolved. 9  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

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

Comment

Related Items