Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Relative stability in critical collapse of scalar fields with angular momentum Akbarian Kaljahi, Arman 2010

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

Item Metadata

Download

Media
24-ubc_2011_spring_akbariankaljahi_arman.pdf [ 534.83kB ]
Metadata
JSON: 24-1.0071512.json
JSON-LD: 24-1.0071512-ld.json
RDF/XML (Pretty): 24-1.0071512-rdf.xml
RDF/JSON: 24-1.0071512-rdf.json
Turtle: 24-1.0071512-turtle.txt
N-Triples: 24-1.0071512-rdf-ntriples.txt
Original Record: 24-1.0071512-source.json
Full Text
24-1.0071512-fulltext.txt
Citation
24-1.0071512.ris

Full Text

Relative Stability in Critical Collapse of Massless Scalar Fields with Angular Momentum by Arman Akbarian Kaljahi  B.Sc., Sharif University of Technology, 2008 M.Sc., The University of British Columbia, 2010  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in The Faculty of Graduate Studies (Physics)  THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) December 2010 c Arman Akbarian Kaljahi 2010  Abstract We study numerically the gravitational collapse of two massless scalar fields in spherical symmetry using an approach which incorporates some of the effects of angular momentum. Each field is characterized by a distinct angular momentum parameter, l, with l = 0, 1, 2, . . ., and has distinct equations of motion. We focus on black critical behaviour, in which we look for solutions, and the properties thereof, that represent the onset of black hole formation in parametrized families of initial data. Although many different critical solutions have been discovered since the early 1990’s, almost all work has involved a single matter source, and by coupling more than one type of matter to the gravitational field, new and interesting questions arise. In particular, we consider the issue of the relative stability of the critical solutions in our model, in which we look for signs of instability of one field in the presence of the critical solution associated with the other. Our main result is that fields corresponding to larger values of l appear to be unstable in the context of lower-l critical solutions. This implies that in a two-field universe, the field with larger l would always dominate at precise criticality.  ii  Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  ii  Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  iii  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  v  List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  vi  Abstract  List of Tables  Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  1  1.1  Gravitational Collapse and Critical Phenomena  . . . . . . . . . . . . . . .  1  1.2  Relative Stability of Critical Solutions . . . . . . . . . . . . . . . . . . . . .  3  2 Mathematical Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . .  7  2.1  3 + 1 Formulation of General Relativity . . . . . . . . . . . . . . . . . . . .  2.2  Einstein’s Equations in Spherical Symmetry  2.3  Massless Scalar Fields with Angular Momentum in Spherical Symmetry . . 16  2.4  Diagnostic Quantities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18  2.5  Summary of Equations  3 Numerical Methods 3.1  3.2  8  . . . . . . . . . . . . . . . . . 13  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22  Solving the Scalar Field Equations . . . . . . . . . . . . . . . . . . . . . . . 22 3.1.1  Finite Differencing . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22  3.1.2  Boundary Conditions  . . . . . . . . . . . . . . . . . . . . . . . . . . 25  Solving the Geometric Equations . . . . . . . . . . . . . . . . . . . . . . . . 26  iii  Table of Contents 3.2.1  Boundary Conditions  . . . . . . . . . . . . . . . . . . . . . . . . . . 27  3.2.2  Black Hole Creation . . . . . . . . . . . . . . . . . . . . . . . . . . . 28  3.2.3  Logarithmic Coordinates  . . . . . . . . . . . . . . . . . . . . . . . . 28  4 Code Tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 4.1  Redundant Equation Evaluation . . . . . . . . . . . . . . . . . . . . . . . . 32  4.2  Conservation of Energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33  4.3  Convergence Test  5 Results  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39  5.1  Critical Collapse: Single Field Case  5.2  Critical Collapse: l1 = 1, l2 = 2 . . . . . . . . . . . . . . . . . . . . . . . . . 41  5.3  Critical Collapse, l1 = 1, l2 = 3 . . . . . . . . . . . . . . . . . . . . . . . . . 43  6 Conclusions  . . . . . . . . . . . . . . . . . . . . . . 39  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51  Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53  iv  List of Tables 5.1  Excitation of the l2 = 2 Critical Solution . . . . . . . . . . . . . . . . . . . . 44  5.2  Excitation of the l2 = 3 Critical Solution . . . . . . . . . . . . . . . . . . . . 44  v  List of Figures 1.1  Critical Curve . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  5  2.1  Foliation of Spacetime . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  8  2.2  Coordinate System in the 3 + 1 Decomposition . . . . . . . . . . . . . . . . 12  3.1  Finite Differencing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23  3.2  Black Hole Formation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30  4.1  Residual r − t . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35  4.2  Residual θ − θ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36  4.3  ADM Mass . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37  4.4  Convergence Test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38  5.1  Mass Gradient at the Critical Regime . . . . . . . . . . . . . . . . . . . . . 45  5.2  Geometry at Critical Regime . . . . . . . . . . . . . . . . . . . . . . . . . . 46  5.3  Excitation in l2 = 2 Critical Solution . . . . . . . . . . . . . . . . . . . . . . 47  5.4  Snapshots of Two Field Critical Evolution . . . . . . . . . . . . . . . . . . . 48  5.5  Critical Curve for l1 = 1 and l2 = 2 . . . . . . . . . . . . . . . . . . . . . . . 49  5.6  Excitation of l2 = 3 Critical Solution . . . . . . . . . . . . . . . . . . . . . . 50  vi  Acknowledgements I am mostly grateful to my supervisor, Dr. Matthew Choptuik, for all the experience and knowledge I have learned during my studies. I also want to thank Dr. Jeremy Heyl and Dr. Bill Unruh for their help and insightful comments. Special thanks to my friend Sanaz Vafaei who has been a great support during writing this thesis, and I am also thankful to the rest of the graduate students in the numerical relativity group with whom I had many useful discussions.  vii  Chapter 1  Introduction 1.1  Gravitational Collapse and Critical Phenomena  Einstein’s equations describe how the presence of matter alters the geometry of spacetime through ten nonlinear second order partial differential equations. Written in tensor form these are Gµν = 8πTµν ,  (1.1)  where Gµν is the Einstein tensor, that includes first and second order derivatives of the metric components, and Tµν is the total energy-momentum tensor for any matter sources [1]. Over the years, many solutions have been found for these equations in various cases. The most relevant for the current work started with numerical investigations by Choptuik [2], who studied the collapse of a massless scalar field in spherical symmetry. These studies led to the discovery of very interesting solutions for the spacetime at the threshold of black hole formation. Specifically, for any initial configuration of collapsing matter that was characterized by one parameter, p, (such as an overall amplitude of the radial profile of the scalar field), Choptuik found a critical value, p∗ , such that for p > p∗ the evolution of matter resulted in collapse and black hole creation, while for p < p∗ the matter dispersed to infinity. Near p = p∗ the spacetime exhibited a very interesting evolution—approaching a universal solution that was independent of the initial shape of the matter distribution. Choptuik also found that the near-critical evolution of spacetime was characterized by discrete self-similarity of the scalar field profile, in the sense that after a certain time δt the field profile had the same shape, but on a length scale smaller by a factor of e∆ , and again after δt/e∆ this so-called echoing behaviour repeated by contraction of the length scale of  1  1.1. Gravitational Collapse and Critical Phenomena the dynamic solution by another factor of e∆ . The discrete self-similarity of the critical solutions is reflected in the periodic behaviour of the scalar field profile (as well as all the quantities derived from it, such as components of the energy momentum tensor) when plotted in the logarithm of spatial coordinate, while the echoing exponent, ∆, is directly related to the period of the oscillations in the logarithmic scale. The echoing exponent ∆ was also found to be universal in the sense that it did not depend on the specifics of the initial data. Furthermore, Choptuik discovered that in the super-critical regime, the mass of the black holes that formed were well described by a scaling law: MBH ∝ |p − p∗ |γ ,  (1.2)  where γ was again a universal exponent. Because black holes of infinitesimal mass can form in this case, the critical solution is called Type II, since the behaviour is analogous to a second order phase transition in statistical mechanics, where the black hole mass is viewed as an order parameter. After the initial discovery of critical phenomena, there has been much additional research done on the behaviour, and the reader is referred to [3] for a comprehensive review of the subject. We note that the subsequent studies have led to the identification of another type of critical solution, known as Type I [4], in which the mass of the black hole starts from a finite value as a function of the tuning parameter, p. In this case, the critical solutions are either static or periodic. Furthermore as the tuning parameter, p, approaches the critical value, p∗ , the length of time the solution remains in the critical regime scales as: τ ∼ σ ln |p − p∗ | .  (1.3)  Here σ is also a universal constant for a particular matter model. Critical solutions of Einstein’s equations are unstable, and for any small perturbation, the system tends to evolve to either a collapsed or dispersed state. However, as described in full detail in [3], within perturbation theory such solutions tend to have only a single growing mode. Furthermore, the growth factor associated with this unstable mode is  2  1.2. Relative Stability of Critical Solutions directly related to the empirically measured quantities, γ and σ, for the type II and type I cases, respectively. Since the original work involving a massless scalar field, critical phenomena have been discovered using a wide variety of different matter models in spherical symmetry: these include a perfect fluid [5–7], a Yang-Mills field [4], and non-linear sigma models [8]. In addition there have been a few studies carried out in axisymmetry, including the collapse of pure gravitational waves [9], a massless scalar field [10], and a complex scalar field with angular momentum [11]. In this thesis we study the critical collapse of two massless scalar fields with angular momentum in spherical symmetry. The scalar fields interact via the gravitational field, and to maintain spherical symmetry while having an effective angular momentum we follow the method introduced in [12]. Using scalar fields with distinct angular momentum as our model, we explore the critical regime where there are two different matter sources coupled to spacetime. The study of such a system gives us some insight regarding the stability of one of the matter sources in the context of the critical solution associated with the other source.  1.2  Relative Stability of Critical Solutions  In scenarios where one type of matter collapses, there is only one tuning parameter, p, and one critical value, p∗ , that controls and distinguishes collapse and dispersal solutions. However, in a case with two different matter sources, with corresponding control parameters p and q, one would expect a “critical curve” to separate dispersal from collapse: this is schematically illustrated in Fig. 1.1. For the case we consider in this thesis, the two different matter sources are scalar fields with different angular momenta: we denote these fields by ψ1 and ψ2 . Again, we assume that we have parametrized the initial data for the two fields such that, individually, they have critical parameter values p∗ and q ∗ , respectively. To be specific, p∗ is the critical parameter value for the first scalar field, ψ1 , in the absence of ψ2 , while q ∗ is the critical parameter value for the second field ψ2 when the first field identically vanishes throughout 3  1.2. Relative Stability of Critical Solutions spacetime. Then as shown in Fig. 1.1, the critical curve connects the points (0, p∗ ) to (q ∗ , 0). Here we introduce the notation (¯ q , p¯)  1  to represent the set of critical points along  the curve. We again emphasize that even though the critical solutions for many matter types have been studied previously, there has been relatively little work done investigating the nature of the solutions along such critical curves. This provides much of the motivation for the current work. The study of a scenario with two types of matter provides a model with an intrinsically broader solution space than encountered in investigations of the critical collapse of a single matter source. First, as just outlined, we will generically have two-dimensional parameter spaces to consider, rather than single-dimensional ones. Since each of the scalar fields have different characteristics, such as the echoing exponent, ∆, at criticality, it is natural to inquire about the spacetime structure along the critical curve in Fig. 1.1. In particular, according to our standard understanding of critical collapse [3] each of the critical solutions, (0, p∗ ) and (q ∗ , 0) are one-mode unstable but with different growth factors. One question that then arises is the stability of the field ψ2 in the presence of the critical solution for pure ψ1 evolution and vice versa. This type of study was initiated by Choptuik in [13], where it was called the “relative stability of critical solutions”. In Choptuik’s case, ψ1 was a Yang Mills field (within the context of the so-called purely magnetic ansatz), ψ2 was a massless scalar field, and the calculations were performed in spherical symmetry. His results suggest that for anything but pure Yang-Mills evolution, the critical solution is the one associated with the massless scalar field. This also indicates that the massless scalar field is intrinsically unstable in the context of the critical YangMills solution. We follow Choptuik’s methodology closely in what follows: this includes studying the limits (¯ q , p¯) → (0, p∗ ) and (¯ q , p¯) → (q ∗ , 0) in an attempt to determine the behavior of small perturbations of one type of matter in the presence of the critical solution associated with the other. 1  We note the ordering (¯ q , p¯) may seem unnatural. However, we prefer it to (¯ p, q¯) since the control  parameter in our numerical experiments will always be q, and there is thus a good argument for making it the first member of the ordered pair.  4  1.2. Relative Stability of Critical Solutions  p✻ p∗..r.......................  .......... .......... .......... ......... ......... ......... ........ ........ ........ ....... ....... ....... ....... ....... ....... ....... ...... ...... ...... ...... ...... ...... ...... ...... .... .... .... .... .... .... .... .... .... .... .... .... .... ... .... .... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...  Black Hole r  (¯ q , p¯)  No Black Hole  r  q∗  ✲  q  Figure 1.1: Evolution with two distinct matter sources generically results in one of two final fates. For cases where both sources are sufficiently weak, all matter will disperse (no black hole); in instances where the gravitational interaction induced by the combined effect of the two sources is sufficiently strong, the evolution results in a black hole. The two regions in the (q, p) parameter space are separated by a curve that we denote the critical curve, and on which the system exhibits critical behaviour. Clearly this curve connects the points p∗ ≡ (0, p∗ ) and q ∗ ≡ (q ∗ , 0).  5  1.2. Relative Stability of Critical Solutions More precisely, we first start with pure evolution of ψ1 , and using a binary search we tune the amplitude p to the critical value p∗ (up to machine precision). Then we turn on the second field with a fixed amplitude q¯ and re-tune the amplitude of the first field to determine the critical point (¯ q , p¯). Finally we study the phenomenology along the critical curve in Fig. 1.1, as well as the limits q¯ → 0 and q¯ → q ∗ . Our results suggest that the type II critical solution associated with the scalar field with larger angular momentum dominates at criticality. Specifically, along the critical curve we consistently find evidence for growth of the scalar field with higher angular momentum, even for very small initial amplitudes (relative to the second field). We thus conjecture that the relative stability seen in our model is analogous to that seen in Choptuik’s work, where the fields with high and low angular momenta correspond to the scalar and Yang-Mills fields, respectively. We now briefly outline the remainder of this thesis. We start with a discussion of the mathematical formulation underlying our model in Chap. 2. In Chap. 3 we describe the numerical methods that were used to approximately solve the equations of motion governing that model. Chap. 4 is dedicated to the various tests performed to validate and calibrate our numerical code. Results are presented in Chap. 5, and we finish with some brief conclusions and suggestions for further research in Chap. 6. We also note that we use units in which c = G = 1, and that we choose the metric signature (−, +, +, +). We also adopt Penrose’s abstract index notation as used, for example, in Wald [15]. In particular, letters from the beginning of the Latin alphabet, {a, b, c, ...}, will generally denote abstract indices. For spacetime components of tensors we use Greek indices {µ, ν, ...} that range over {0, 1, 2, 3} (0 being the time index), while for spatial components ranging over {1, 2, 3}, we use the set {i, j, k, l, m, n}.  6  Chapter 2  Mathematical Formulation In this chapter we derive and describe the equations that govern our system of two massless scalar fields with angular momentum that interact through the gravitational field. This system consists of two sets of coupled partial differential equations: Einstein’s equations which govern the metric functions, and modified wave equations (on a curved spacetime) for the scalar fields. In the first section, we start with a review of the “3 + 1”, or Arnowitt-Deser-Misner (ADM), formulation of Einstein’s equations [20]. This formalism provides an evolutionary picture of the intrinsic and extrinsic geometries of three-dimensional spatial slices in time, thereby providing a description of the whole spacetime. The ADM decomposition yields four constraint equations—namely the Hamiltonian and momentum constraints—and six evolution equations which contain second time derivatives of the metric coefficients. The constraint equations need to be satisfied on each slice, including the initial time slice. The task of determining initial conditions for the geometry and matter fields that are consistent with the constraints is known as the initial data problem in relativity. Once initial data have been set, and coordinates have been chosen, the evolution equations provide the system of partial differential equations (PDEs) to evolve the metric functions in time. In Sec. 2 we simplify the ADM equations assuming spherical symmetry and a particular choice of coordinates. The matter content of the system is discussed in Sec. 3. Following the work done by Olabarrieta et al [12], we describe how to retain spherical symmetry while allowing the scalar fields to “mock up” some of the effects of angular momentum. The result of this procedure is modified wave equations (Klein-Gordon equations) for the massless scalar fields. Finally, we provide expressions for the needed components of the energy-momentum tensor. This completes the equations for the coupled system of  7  2.1. 3 + 1 Formulation of General Relativity  ✄✄✗n  µ ...................................................................................... ..................... ................... .............. ..................... ..................... ......... ................................. .... .... ................................................................................. .... .... .... .... .... .... .... .... ... ........................................................................................................ ... . . . . . . . . . . . . . . . . . . . . . . . . . . . . ... . . . . . . . . . ................... ... .......... . . . . . . . . . . . ... . . . . . . . . . . . . . . . . . . . . . . . t2 ........ ..................... ... .......... ........ . . . . . . ... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .................................................................................. .... .... ... ... .... . . . .... .... .... .... .... . . .... .... .... .................................................................. . . . . . . . . .... . . . . . . . . . . . . . . .... ............... .... .. ......... ..... ............ . .. ............................................................................................................................ ... ..... .......... . . . . . . . . . . . . . . . . . . . . . . . . . . . . ...... . . . . . . . . . . . . . . . . ................... .......... ................ ... ............ ... . . . . . . . . . . . . . . . . . . . . . . . ... . . . . . . . . . . . . . . .......... .................... ... ... t1 ......... .......... ......... . . . . . .......... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .............................................................................................. ... .... .... ............... .... .............. . . .... . . . .... . . . ...................... .... .................... .... .... . . . . . . . . . . . . . . . . . . ............................................. .... .... .... .... .............................................................................. .... .... .... .... .............. ................... .... ... .............. ............ .... ... ..... ... .......... ...... ..................... ... ...... .......... ... ............ ... ... . .......... ... ... .......... ........... . . . ................... ............ ... . . . . . . . . ................. .... ................................... ................................................... ...... .... ...... .... ................................................. .... ............................. ................ .... .... ................ ............ .... ..... ............ ........... ..... .......... ...... .................... ...... .......... ....... .... .......... .......... . ............ . ........ . . . . ............. . . . . . . . ...... ................. ......................................................................................  Σt3 ........................................................................  ✄  ✄  Σ  Σ  t3 t2 t1  Figure 2.1: Foliation of spacetime into 3 dimensional space-like surfaces. geometric and matter fields.  2.1  3 + 1 Formulation of General Relativity  We assume that the spacetime manifold can be foliated into a family of space-like three dimensional surfaces, Σt , corresponding to constant slices of a scalar function t, where t can be interpreted as a global time function. The geometry of the 4-dimensional spacetime is then translated into the intrinsic geometry of three-surfaces (the “3-metric”) and a description of how these slices are embedded in the four-space (the extrinsic curvature tensor). In terms of these quantities, Einstein’s equations then become a set of constraint and evolution equations for the 3-metric and the extrinsic curvature. For a foliation (t,Σt ) we can define a one-form, Ωa , by Ωa = ∇a t .  (2.1)  The 4-metric gab then allows us to calculate the norm of Ω as ||Ω||2 = g ab ∇a t∇b t ≡ −  1 , α2  (2.2)  where α is known as the lapse function, since it measures how fast proper time elapses along the normal vector Ωa between two slices Σt and Σt+dt . We next define a unit-normalized vector, na , by na ≡ −αg ab Ωb .  (2.3) 8  2.1. 3 + 1 Formulation of General Relativity Then na is a time-like, future-pointing vector normal to the spatial hypersurfaces: na na = α2 g ab Ωa Ωb = −1 .  (2.4)  Here we assume that α > 0, and the negative sign in (2.2) is chosen so that na is timelike and the 3-surface Σt is spacelike everywhere. Given the normal vector we can now construct the projection tensor, γba , that can be used to project arbitrary tensors into the hypersurfaces: γba = δba + na nb .  (2.5)  Calculating the operation of γba on a vector W b yields: γba W b = γba (W⊥b + W b ) = (δba + na nb )((W⊥b + W b ) = W⊥a + W a + na (nb W⊥b + nb κnb ) = W⊥a + κna + na (nb κnb ) = W⊥a .  (2.6)  Here W⊥a and W a denote the components of W a that are parallel and perpendicular, respectively, to nb . Then, by definition, W⊥a na is zero and W a ≡ κna . Equation (2.6)  shows that γba is indeed the desired projection tensor. Furthermore we can apply the projection operator to other tensors. For instance, we define: ⊥Tab ≡ γac γbd Tcd .  (2.7)  Then using the fact that γba nb = 0, it is clear that the contraction of ⊥Tab with any vector parallel to na is zero, and that ⊥Tab is therefore a spatial tensor. The spatial 3-metric is introduced by projecting the spacetime metric onto the hypersurface: γab ≡ γac γbd gcd = (δac + na nc )(δbd + ndb )gcd = gab + na nb .  (2.8) 9  2.1. 3 + 1 Formulation of General Relativity γab defined in this fashion encodes the intrinsic geometry of the hypersurfaces, Σt . Given the projection operator, the spatial derivative operators are defined by projecting the four dimensional covariant derivative onto the hypersurfaces: Da = γab ∇b .  (2.9)  It is then easy to show that Da is compatible with the 3-metric: Da γbc = 0.  (2.10)  Once the derivative operator on the hypersurfaces is in hand, we can define the spatial Riemann tensor associated with Da : (Da Db − Db Da )ωc ≡ 3 Rabc d ωd .  (2.11)  In precise analogy to the 4-dimensional case, we can also define the spatial Ricci tensor: 3  Rac = 3 Rabc b ,  (2.12)  and the spatial Ricci scalar: 3  R = 3 Ra a .  (2.13)  Finally, we define the extrinsic curvature as (up to a sign) the projection of the gradient of the normal vector onto the spatial hypersurface: Kab ≡ −γad γbd ∇c nd = −Da nb .  (2.14)  Equivalently, Kab can be expressed in terms of the Lie derivative of the 3-metric: 1 Kab = − Ln γab ⇒ ∂t γab = −2αKab + Da βb + Db βa . 2  (2.15)  The 3-metric and extrinsic curvature (γab , Kab ) completely determine the spacetime geometry and are to be considered as dynamical variables evolving in time. On any hypersurface they define the “instantaneous” state of the gravitational field, with the 3-metric and extrinsic curvature corresponding to position and velocity (momentum) variables, respectively, if we consider the geometry of spacetime as a dynamical system. 10  2.1. 3 + 1 Formulation of General Relativity We can now write the Einstein equations in terms of these dynamical variables, (γab , Kab ). Contracting the field equations with na twice, Gab na nb = 8πTab na nb , yields the Hamiltonian constraint: 3  R + K 2 − Kab K ab = 16πρ .  (2.16)  Here K is the trace of the extrinsic curvature and ρ is defined by ρ ≡ na nb T ab .  (2.17)  Similarly, contracting Einstein’s equations once with na , and projecting onto the hypersurface, γda Gdc nc = 8πγda T dc nc , gives the momentum constraints: Db K ab − Da K = 8πj a .  (2.18)  j a ≡ γda T dc nc .  (2.19)  In the above, j a is defined by  The quantities ρ and j a can be interpreted as the local energy and momentum densities, respectively, that are measured by an observer moving orthogonal to the hypersurfaces, Σt . Equations (2.16) and (2.18) are called constraint equations since they only involve spatial tensors and spatial derivatives of such tensors and, in particular, do not contain second time derivatives of dynamical variables. Nonetheless, they must be satisfied on each time slice, including the initial hypersurface. Appropriate initial data for the 3-metric and extrinsic curvature is thus restricted to a set (γab , Kab )(t=0,x) that satisfies the Hamiltonian and momentum constraints. In order to derive the evolution equation, we first define a coordinate system. Naturally, the function t used for slicing the spacetime is adopted as the time coordinate. The vector field, ta , tangent to t = const. trajectories thus satisfies: t a ∇a t = 1 ,  (2.20)  but is not necessarily orthogonal to the hypersurfaces. In fact, we can express ta as a linear combination of its component, λna , normal to the 3-surface and a purely spatial 11  2.1. 3 + 1 Formulation of General Relativity  ................................................................................................................... Σt+dt ................... ........................ ........... .....................  Σt  .............. ..................... ...... ..................... ......... ................................. ..... .... ................................................................................. .... .... i .... .... x a .... .... β dt .... .... ... i ... ... dx ... ... ........... ... ... ... ... µ ... ... αn dt .... .... i i .... .... x + dx .................................................................. . . . . . . . . .... . . . . . . .... . . . . . . ............... ..... a ............................................................................. . . . . . . ............ ............................... t dt ..... ..................... .... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ...... . . . . . . . . . . . . . . ................... .......... ................ .. ................... .......... . ....................... . . . . . . . . . . . ...... ........................ds .......... .......... ......... . ...... . . . . . . . . . . . . .......... . . . . . . . . . . . . . . . . .............................................................................................. .... .... ....... . . . . . . .............. . . . . .... . . . . . ...................... .................... .... . . . . . . . . . . . . . . . . . ............................................. .... .......................................................... .... .... . .............................. .... ... ... xi ... ... ... ... ... ... ... ... ... .... .... .... .... ................................................. . . . . . . . . . . . . . . . . . . .... . . . . . . . . . . .... ................ ........... . . .... . . . . . .... . . . . . . . . . . . ............ ..... ...... ..... .......... ...... .................... ...... .......... ....... .......... ..... .......... .......... . ............ ......... . . . . . . . . . ............... . . ......... . . . . ..................... . . . . . . . . . . . . . . ......................................................  ✏❍ ✏✶ ✏ ✁✕ ❍ ❳❳ ❳ ✻ ✁ ❳❳❍ ❳❍ ❳❍ ✁ ❥ ❍ ❳ ③ ❳ ✶ ✏ ✏ ✁ ✏✏ ✏ ✁ ✏✏ ✁✏✏  t + dt  t  Figure 2.2: The vector field ta is not necessarily perpendicular to the spacelike hypersurfaces. The shift vector, β a , which is a spatial (hypersurface-intrinsic) vector, encodes the difference between the normal direction and ta , while the lapse measures the ratio of proper time to coordinate time, relative to propagation in the orthogonal direction. vector, β a . Then by virtue of (2.3) and (2.20), we can deduce that the coefficient λ is in fact the lapse function: ta = αna + β a .  (2.21)  Equation (2.21) can be considered as the definition of the “shift vector”, β a , which describes how the spatial coordinates, xi , shift from slice to slice (relative to normal propagation) as illustrated in Fig. 2.2. By fully projecting Einstein’s equations onto the 3-surface, γca γdb Gcd = 8πγca γdb T cd , we get the evolution equation for extrinsic curvature: Lt Kab = − Da Db α + α( 3 Rab − 2Kac Kbc + KKab ) 1 + 8πα(Sab − γab (S − ρ)) + Lβ Kab . 2  (2.22)  Here Sab , known as the spatial stress tensor, is defined by Sab ≡ γac γbd Tcd ,  (2.23)  S is the trace of Sab , and L again denotes the Lie derivative. For our purposes, the following  12  2.2. Einstein’s Equations in Spherical Symmetry version of the above evolution equation that uses mixed indices is more convenient Lt K a b = − Da Db α + α( 3 Ra b − 2K a c Kbc + KK a b ) 1 + 8πα(S a b − γ a b (S − ρ)) + Lβ K a b . 2  (2.24)  Together, equations (2.15) and (2.22) completely determine the evolution of the gravitational field (γab , Kab ), or, equivalently, and from now on, (γab , K a b ). Here we emphasize once more that the initial data for the system must satisfy the Hamiltonian and momentum constraints. We note that, in general, there are no evolution equations for the lapse and shift: this can be understood in terms of the fact that these functions encapsulate the coordinate degrees of freedom in the 3 + 1 decomposition (i.e. they represent the general covariance of the theory). Thus, in constructing a spacetime using the 3 + 1 approach, the lapse and shift must be completely specified in some fashion in order to fix the coordinatization of spacetime. We finish this section by noting that we can relate the components of the 4-metric to those of the 3-metric, as well as the lapse and shift components, using a four-dimensional version of the Pythagorean theorem. Specifically, from Fig. 2.2 we see that the spacetime line element can be written as ds2 = −α2 dt2 + γij (dxi + β i dt)(dxj + β j dt) .  2.2  (2.25)  Einstein’s Equations in Spherical Symmetry  In this section we specialize the 3 + 1 equations to the case of spherical symmetry. We start with a brief discussion of the most general form of the spherically-symmetric line element, then further specialize to a particular coordinate system known as “polar-areal” which we have used for all of the calculations described in this thesis. We choose spherical polar coordinates (t, r, θ, φ) adapted to the spherical symmetry. Then a suitably general form of the line element is ds2 = (−α2 + a2 β 2 )dt2 + 2a2 βdtdr + a2 dr2 + r2 b2 dΩ2 ,  (2.26)  where α, β, a and b are all functions of t and r alone, dΩ2 is the line element on the unit 2-sphere, and the shift vector β i has components (β r , 0, 0) ≡ (β, 0, 0). 13  2.2. Einstein’s Equations in Spherical Symmetry For the components of the extrinsic curvature tensor, it can be shown that we must, in general, have K ij  = diag K r r (t, r), K θ θ (t, r), K φ φ (t, r) = diag K r r (t, r), K θ θ (t, r), K θ θ (t, r) .  (2.27)  That is, there are only two independent components of the extrinsic curvature tensor in spherical symmetry. We now completely fix the spatial coordinate, r, by demanding that it measure proper surface area (i.e. that it be “areal”). With respect to the general line element (2.26), this means that we must have b(t, r) ≡ 1. Using b(t, r) = 1, we can then compute the non-vanishing components of extrinsic curvature: Krr =  1 aα  (aβ) −  K θ θ = K φφ =  ∂a ∂t  β . rα  (2.28) (2.29)  We can now fix the time coordinate by requiring that the extrinsic curvature have only an r-r component. We thus have K θ θ = K φφ = 0 .  (2.30)  From (2.29), and assuming a non-vanishing lapse, we then find β(t, r) ≡ 0 .  (2.31)  This choice of time is known as polar slicing. Our general spherically symmetric line element (2.26) has thus become ds2 = −α2 dt2 + a2 dr2 + r2 dΩ2 , while the 3-metric has the simple form:  a2 0 0   γij =  0 r2 0  0 0 r2 sin2 θ      .   (2.32)  (2.33)  14  2.2. Einstein’s Equations in Spherical Symmetry Finally, the extrinsic curvature is     Krr 0 0 −a/(αa) ˙ 0 0         K ij =  0 0 0 = 0 0 0 ,     0 0 0 0 0 0  (2.34)  where we have adopted the notation A′ ≡ ∂A/∂r and A˙ ≡ ∂A/∂t.  Having the 3-metric, we can compute the non-vanishing components of the spatial Ricci tensor: a′ a3 r  3  Rr r = 2  3  Rθ θ = Rφ φ =  (2.35) a′ 1 + 2 3 ra r  1−  1 a2  (2.36)  Using these last results, as well as (2.34) in (2.16), we find, after some manipulation, that the Hamiltonian constraint can be written as the following ordinary differential equation for the 3-metric component, a:  1 − a2 a′ = − 4πra2 T t t . a 2r  (2.37)  Turning to the momentum constraints (2.18), only the radial component is non-trivial, and, using (2.15) it can be written as  a˙ = 4πra2 T r t . a  (2.38)  We note that in these last two equations T denotes the total energy momentum tensor for the two scalar fields. To ensure that the choice (2.30) of a polar time coordinate is maintained in time, it is sufficient to demand that K˙ θ θ = 0 .  (2.39)  Using this requirement in the θ-θ component of (2.22), along with other results from above, we get:  α′ a2 − 1 = + 4πra2 T r r . α 2r  (2.40)  This is to be viewed as an ordinary differential equation for the lapse function (the so-called slicing condition), which must be satisfied on each hypersurface. 15  2.3. Massless Scalar Fields with Angular Momentum in Spherical Symmetry  2.3  Massless Scalar Fields with Angular Momentum in Spherical Symmetry  We now turn to the derivation of the equations of motion for the massless scalar fields, and the computation of the associated components of the energy-momentum tensor. Following Olabarrieta et al [12] we separate the angular and radial dependencies of a generic scalar field, Ψm l , having angular momentum by writing: l Ψm l (t, r, θ, φ) = ψ (t, r)Qlm (θ, φ) .  (2.41)  Here Qlm are the real eigenfunctions of the angular part of the flat-space Laplacian with eigenvalue l(l + 1), l = 0, 1, 2, . . . while m encodes the 2l + 1 distinct eigenfunctions corresponding to m = −l, −(l − 1), . . . , l − 1, l. The Qlm can be expressed in terms of the usual spherical harmonics, Ylm :   Yl0    √1 (Ylm + (−1)m Yl−m ) Qlm = 2     √1 Y |m| Y l|m| − (−1) l−|m| i 2  for  m = 0,  for  m > 0,  for  m < 0,  (2.42)  The energy-momentum tensor corresponding to Ψm l is: (lm)  1 m c m m Tab = ∇a Ψm l ∇b Ψl − gab ∇ Ψl ∇c Ψl , 2  (2.43)  and since Ψm l has angular dependence, this tensor would also depend on θ and φ in general. Therefore, without further restriction, the coupled equations of motion for the matter and geometry are not spherically symmetric in this form. In order to keep our calculations spherically symmetric while retaining at least some of the effects of angular momentum, we consider an energy-momentum tensor corresponding to a sum over all admissible m for any given l l (l)  (lm)  Tab =  Tab .  (2.44)  m=−l  Then, by construction,  (l) T  ab  is spherically symmetric and only depends on a single radial  profile, ψ l (t, r). It is these ψ l that will constitute our dynamical matter fields in the spacetimes we will be constructing. 16  2.3. Massless Scalar Fields with Angular Momentum in Spherical Symmetry To derive the equation of motion for the field, ψ l (t, r), which we denote simply by ψ in the remainder of this section, we can use the fact that the divergence of the energy momentum-tensor,  (l) T  ab ,  must be zero: g ac ∇c  (l)  Tab = 0.  (2.45)  ∂ψ , ∂r  (2.46)  a ∂ψ , α ∂t  (2.47)  Introducing the following auxiliary variables Φ≡  Π≡  and using the form (2.32) for the metric, we find the following equations of motion for a scalar field with angular momentum l: ∂ ∂Φ = ∂t ∂r  α Π , a  α 1 ∂ ψ ∂Π r2 Φ − l(l + 1)aα 2 . = 2 ∂t r ∂r a r  (2.48)  (2.49)  Furthermore, the corresponding nonzero components of the energy momentum-tensor are (l)  T tt = −  (2l + 1) 1 2 ψ2 2 (Π + Φ ) + l(l + 1) , 8π a2 r2 (2l + 1) 2α ΠΦ, 8π a3  (2.51)  ψ2 (2l + 1) 1 2 2 (Π + Φ ) − l(l + 1) , 8π a2 r2  (2.52)  (l)  (l)  T rr =  (l)  (2.50)  T θθ =  T rt =  (l)  T φφ =  (2l + 1) 2 (Π − Φ2 ) . 8πa2  (2.53)  The total energy-momentum tensor is simply the sum of the two energy-momentum tensors associated with each of the fields:  T µν =  (l1 )  T µ ν [Π1 , Φ1 , ψ1 ] +  (l2 )  T µ ν [Π2 , Φ2 , ψ2 ] .  (2.54) 17  2.4. Diagnostic Quantities We note that the term l(l + 1)aαψ/r2 that appears in (2.49) plays the role of an effective angular momentum barrier, in analogy with the 1-dimensional reduced problem of a particle moving in a central potential [12].  2.4  Diagnostic Quantities  In this section we describe various quantities that have been used to check the code that we constructed to approximately solve the equations derived above. First, in order to monitor the energy (mass) content of the system, we introduce the mass (or mass aspect) function as follows: M (t, r) ≡  r 2  1−  1 a2  .  (2.55)  The motivation for this definition can be seen by comparing the form of the metric in polar-areal coordinates (2.32) with the usual Schwarzschild metric: 2M ds = − 1 − r 2  2  −1  2M dt + 1 − r 2  a (t, r) →  2M 1− r  −1  .  dr2 + r2 dΩ2 ,  (2.56)  (2.57)  In the limit r → ∞ the mass function approaches the Arnowitt-Deser-Misner (ADM) mass  lim M (t, r) = MADM  r→∞  (2.58)  which is a strictly conserved quantity. Physically, M (t, r) measures the amount of mass inside a sphere with radius r at time t: in vacuum regions it is thus (locally) conserved and this fact can be used to advantage in the testing of a spherically symmetric code. The other quantities that are used as a gauge of the accuracy of our numerical calculations come from the fact that, due to the general covariance of general relativity, Einstein’s equations are an overdetermined set. In the current case, for example, only one of the two equations (2.37) and (2.38) is needed to find the metric function, a. Now, one can think of the Hamiltonian constraint (2.37), as the t-t component of Einstein’s equations, and the momentum constraint, (2.38) as the t-r component. Due to the assumption of 18  2.4. Diagnostic Quantities spherical symmetry, the geometry of each spatial slice is completely determined by the energy-momentum tensor on that slice. The field equations are called “fully constrained” in such an instance, since only the matter content—and not the history of the gravitational field—is needed to compute the geometric functions on each slice. The absence of gravitational waves in spherical symmetry is also intimately related to this observation. Since we have two equations for a that are equivalent to each other, one should be trivially satisfied if the other one is solved correctly. As we will describe in Chap. 3, we choose to compute a by solving the Hamiltonian constraint at each time. Therefore, as an accuracy test, we directly monitor the degree to which the t-r component of the Einstein equation,  Gr t ≡  2∂a/∂t = 8πT r t , 3 ra  (2.59)  is satisfied. Here T r t is given by (2.51) and (2.54). Similarly, the θ-θ component of Einstein’s equations. Gθ θ ≡ − −  1 ∂α ∂a ∂α ∂a + α3 + α2 r −α2 a 3 3 rα a ∂r ∂r ∂r ∂r 2 2 ∂ α ∂ a ∂α ∂a = 8πT θθ , α2 ar 2 + a2 αr 2 − a2 r ∂r ∂t ∂t ∂t  (2.60)  where T θ θ is given by equations (2.54) and (2.53), can be viewed as an evolution equation for the extrinsic curvature component K r r . We will use it, however, as a check of the validity and accuracy of our numerical code. Since (2.60) involves all of the geometric and scalar field variables, it forms the basis for a strong test of our numerical solutions.  19  2.5. Summary of Equations  2.5  Summary of Equations  Our equations of motion for a system of two massless scalar fields with angular momentum in spherical symmetry, and interacting via the gravitational field are given by: ∂ ∂Φ1 = ∂t ∂r  α Π1 , a  (2.61)  ∂Π1 ψ1 1 ∂ α = 2 r2 Φ1 − l1 (l1 + 1)aα 2 , ∂t r ∂r a r ∂Φ2 ∂ = ∂t ∂r  (2.62)  α Π2 , a  (2.63)  α ψ2 1 ∂ ∂Π2 r2 Φ2 − l2 (l2 + 1)aα 2 . = 2 ∂t r ∂r a r  (2.64)  Here a and α are components of the spherically symmetric metric in polar-areal coordinates: ds2 = −α2 dt2 + a2 dr2 + r2 dΩ2 ,  (2.65)  and satisfy the (overdetermined) set of equations  − −  a′ 1 − a2 = − 4πra2 T t t , a 2r  (2.66)  a˙ = 4πra2 T r t , a  (2.67)  a2 − 1 α′ = + 4πra2 T r r , α 2r  (2.68)  1 ∂α ∂a ∂α ∂a −α2 a + α3 + α2 r rα3 a3 ∂r ∂r ∂r ∂r 2 2 ∂ α ∂ a ∂α ∂a α2 ar 2 + a2 αr 2 − a2 r = 8πT θθ . ∂r ∂t ∂t ∂t  (2.69)  In this last set of equations the energy-momentum components, T µ ν , are given by: T tt = − −  ψ12 (2l1 + 1) 1 2 2 (Π + Φ ) + l (l + 1) 1 1 1 8π a2 1 r2 (2l2 + 1) 1 2 ψ22 2 (Π + Φ ) + l (l + 1) 2 2 2 8π a2 2 r2  ,  (2.70) 20  2.5. Summary of Equations  T rt =  (2l2 + 1) 2α (2l1 + 1) 2α Π1 Φ1 + Π2 Φ2 , 3 8π a 8π a3  T rr = +  T θ θ = T φφ =  (2l1 + 1) 1 2 ψ12 2 (Π + Φ ) − l (l + 1) 1 1 1 1 8π a2 r2 ψ2 (2l2 + 1) 1 2 (Π2 + Φ22 ) − l2 (l2 + 1) 22 2 8π a r  (2.71)  ,  (2l1 + 1) 2 (2l2 + 1) 2 2 (Π − Φ ) + (Π2 − Φ22 ) . 1 1 8πa2 8πa2  (2.72)  (2.73)  21  Chapter 3  Numerical Methods In this chapter we describe the various numerical methods we have used to solve the system of coupled partial differential equations (2.61-2.73). First, the finite differencing scheme used to solve the equations of motion for the scalar fields is discussed, along with the treatment of boundary conditions for those fields at r = 0 and large r (r → ∞). In the second section, the numerical solvers for the geometric quantities a and α are described, again with the boundary conditions that are imposed to determine unique solutions. We then discuss the criterion we use to determine when a black hole is forming in a calculation, and end with a brief discussion of a change of radial coordinate that we have implemented in order to provide more spatial resolution in the central region of the computational domain.  3.1  Solving the Scalar Field Equations  To approximately solve the evolution equations for either of the scalar fields, we first use a second order finite-difference discretization to convert the continuum problem into a set of algebraic equations. We then solve the algebraic system to some tolerance using an iterative technique.  3.1.1  Finite Differencing  As shown in Fig. 3.1, we use the notation fin for the discretized value of a function f (t, x) at discrete time, tn , and discrete radial location ri . We note that, in principle, our differential equations should be solved on the domain 0 ≤ r < ∞: however, since computer resources are finite, our finite difference mesh cannot extend to infinity unless we employ spatial compactification or an equivalent technique. We thus truncate the spatial grid at some 22  3.1. Solving the Scalar Field Equations  fin+1 ✉  ✉  n fi−1  fin  ✉  tn+1  ✉  n fi+1  ✉  ✉  tn  ✉  tn−1  fin−1 ✉  ✉  ri−1  ri  ri+1  Figure 3.1: Finite difference discretization of the function f (t, r) on the r-t plane finite radius, rmax , with the understanding that we are to ensure that our numerical results do not depend significantly on the location of the outer boundary. Once a finite difference grid has been introduced we can replace the derivative operators with the following second order (in the temporal and spatial mesh spacings) algebraic operators:  n f n − fi−1 ∂f n (t , ri ) = i+1 + O(∆r2 ) , ∂r 2∆r  (3.1)  f n+1 − fin ∂f n+ 1 (t 2 , ri ) = i + O(∆t2 ) . ∂t ∆t  (3.2)  We also define the second order time-averaging operation: 1 1 f (tn+ 2 , ri ) = (fin + fin+1 ) + O(∆t2 ) . 2  (3.3)  Using this last operator, we can approximate the equations of motion for the scalar 1  fields, (2.61 -2.64), with a difference scheme centred at (tn+ 2 , ri ). Such a discretization is known as a Crank-Nicholson scheme. For instance, (2.61) is discretized by φ1 n+1 − φ1 ni i ∆t  =  1 2  n+1 n+1 n+1 n+1 n+1 αi+1 Π1 n+1 i+1 /ai+1 − αi−1 Π1 i−1 /ai−1 2∆r  +  1 2  n Π n /an − αn Π n /an αi+1 1 i+1 i+1 i−1 1 i−1 i−1 2∆r  (3.4)  where the left and right hand sides in the above are second order approximations of ∂φ1 ∂t  1  (3.5)  (tn+ 2 ,ri )  23  3.1. Solving the Scalar Field Equations and ∂ ∂r  α Π1 a  ,  1  (3.6)  (tn+ 2 ,ri )  respectively. Furthermore, we relate the mesh spacing in the time direction, ∆t, to the spatial mesh size, ∆r, through the so-called Courant condition: ∆t = λ∆r .  (3.7)  Here, the Courant number λ will generally be held fixed when we change the resolution of the finite difference mesh, so that the overall difference scheme is characterized by a single discrete scale. We will sometimes denote this scale by h and identify it with the spatial spacing ∆r. For difference schemes and/or solution methods of the discrete equations such as those we adopt, there will often be restrictions on the size of λ in order that the scheme be numerically stable. For low values of l we typically choose λ ≈ 0.3, but smaller values are needed as l increases. Once we have replaced all the continuum variables with discrete equivalents and derivative operators with corresponding second order accurate finite differences, we can view the scheme as constituting algebraic equations for each of the variables f at the grid point (tn+1 , ri ). We write that scheme as Dfin+1 = 0  (3.8)  where the details of the operator D will depend on the particular discrete variable and, to a certain extent, on the spatial location being considered. We solve the system of equations (3.8) using a single-step, pointwise Newton-GaussSeidel relaxation scheme [21] as follows. We compute the value of the variable at the advanced time step, fin+1 , iteratively, and denote the approximation to this value at any stage of the iteration by f˜in+1 . As an initial guess we set f˜in+1 = fin . The single-step iteration then proceeds by updating the unknown using Rf˜n+1 i . f˜in+1 = f˜in+1 − Jf˜n+1  (3.9)  i  Here, Rf˜n+1 is the residual of the difference equation evaluated at the current iteration, i  i.e. Rf˜n+1 = Df˜in+1 , i  (3.10) 24  3.1. Solving the Scalar Field Equations and Jf˜n+1 is the diagonal element of the Jacobian matrix associated with (3.8): i  Jf n+1 = i  ∂D , ∂fin+1  (3.11)  again evaluated using the latest iterate, f˜in+1 . As the name suggests, and assuming that it converges, the single-step, pointwise Newton-Gauss-Seidel technique solves the discrete equations (3.8) (for all values of i) by performing one step of Newton’s method for the unknown fin+1 while keeping the values of all other unknowns fixed. A complete iteration—known as a relaxation sweep—then consists of taking such steps for all variables and at all spatial locations. Relaxation sweeps continue until some convergence criterion, typically based on the size of the residuals, is attained.  3.1.2  Boundary Conditions  In addition to discretization of the equations of motion for the scalar fields, we need to deal with boundary conditions for these fields at r = 0 and r = rmax . The boundary conditions at the origin come from the demand that the solution be regular there, while the outer boundary conditions are derived from a requirement that there be no incoming radiation (waves) at r = rmax . Equivalently, we require that the solution at r = rmax describe purely outgoing waves. Regularity of the scalar fields at the origin is enforced via ψ(t, 0) = O(rl ) ,  (3.12)  Π(t, 0) = O(rl ) ,  (3.13)    O(rl−1 ) Φ(t, 0) =  O(r)  for l ≥ 1,  (3.14)  for l = 0.  Such boundary conditions are implemented numerically by imposing continuity of the value of the scalar field, or its derivatives, depending on the value of the angular momentum, l.  25  3.2. Solving the Geometric Equations At r = rmax we wish the scalar field solutions to describe outgoing-only waves: this approximates so-called Sommerfeld conditions which would be achieved in the limit r → ∞. As previously discussed, we must choose rmax sufficiently large that our results do not depend significantly on its magnitude, and the artificial reflection of outwards propagating radiation off the outer boundary is a chief concern in this respect. However, in our calculations, choosing rmax sufficiently large has not been a major problem, and discretized versions of the following conditions, which are based on 1/r fall off of the scalar fields and an assumption of flatness of the spacetime at the outer boundary, work well:  3.2  ∂Φ Φ(t, rmax ) ∂Φ (t, rmax ) + (t, rmax ) + = 0, ∂t ∂r rmax  (3.15)  ∂Π Π(t, rmax ) ∂Π (t, rmax ) + (t, rmax ) + = 0. ∂t ∂r rmax  (3.16)  Solving the Geometric Equations  We recall that the spacetime metric in polar-areal coordinates is given by ds2 = −α2 dt2 + a2 dr2 + r2 dΩ2 ,  (3.17)  and that the spacetime geometry is thus completely determined by knowledge of the two functions a(t, r) and α(t, r). We choose to use the Hamiltonian constraint (2.66) to compute a, while the lapse function is fixed by the polar slicing condition (2.68). Now, equation (2.66) is a first order ordinary differential equation (ODE) of the form 1 ∂a = f (r) + g(r) a2 (r) , a ∂r  (3.18)  where f (r) and g(r) are given functions of the scalar fields, the radial coordinate and the values of the angular momenta for each of the fields. Defining a new variable, A, by A = ln(a) ,  (3.19)  the equivalent differential equation for A becomes ∂A(r) = f (r) + g(r) exp(2A(r)). ∂r  (3.20) 26  3.2. Solving the Geometric Equations To solve (3.20), we use the following second order implicit finite-difference scheme 1 1 Ai+1 − Ai = (fi+1 + fi ) + (gi+1 + gi ) exp(Ai+1 + Ai ) . ∆r 2 2  (3.21)  Here we have again adopted the notation fi ≡ f (ri ) ≡ f (i∆r), and we note that the scheme is centred at the midpoint, r = ri+ 1 . Assuming that Ai is known, (3.21) can be 2  solved for Ai+1 using Newton’s method. Once a has been determined, we use a similar second order finite difference technique to integrate the ODE (2.68) for the lapse function explicitly (i.e. Newton method’s is not required). The finite differencing of (2.68) for B ≡ ln(α) yields Bi+1 − Bi Hi+1 + Hi = , ∆r 2  (3.22)  where H(r) is given by the right hand side of (2.68): H(r) =  a2 − 1 + 4πa2 T r r . 2r  (3.23)  Furthermore we perform the integration for the Bi from r = rmax to the origin, since, as described in the next sub-section, we use a boundary condition at r = rmax to fix the solution.  3.2.1  Boundary Conditions  To solve (2.66) we need one boundary condition which is provided by the requirement that the spacetime be locally flat (regular) at the origin: a(t, 0) = 1 .  (3.24)  Since the equation for α is also first order, we need one boundary condition for it as well. Even though the slicing condition is coupled to the equation for a, it is linear in α: therefore for any given solution, α(t, r), the function kα(t, r), where k is an arbitrary positive constant is also a solution. This reflects the fact that we are free to relabel the polar slices via t → F (t)  (3.25)  27  3.2. Solving the Geometric Equations without changing the form (2.32) of the spacetime metric. Our specific choice uses α(t, rmax )a(t, rmax ) = 1 ,  (3.26)  which can be motivated by a requirement that our metric approach the Schwarzschild form (2.65) as r → ∞. This choice also ensures that, as r → ∞, t measures proper time.  3.2.2  Black Hole Creation  For initial data with sufficiently large amplitude(s) of the scalar fields, or which otherwise develops strong gravitational interactions, black hole formation is possible. In such instances, and as the matter collapses toward the origin, the value of 2M (t, r)/r asymptotes to 1 at some r = rBH , and as t → ∞. Here we recall that M (t, r) is the mass function given by 2M (t, r) = r  1−  1 a2 (t, r)  .  (3.27)  We use this behaviour as our criterion for the creation of a black hole with radius rBH , or equivalently mass MBH = rBH /2. Here we must point out that due to our coordinate choice, we cannot detect the formation of apparent horizons—which are often used in numerical relativity calculations as proxies for black holes—in our calculations. Nonetheless, past experience with computations in polar-areal coordinates [2] have shown that black hole masses can be estimated quite well using the technique just described. We also point out that when black hole formation is imminent the values of a(t, r) become quite large in the vicinity of r = rBH , and the spatial derivatives of a even more so. This can lead to a failure of the discretized Hamiltonian constraint, in which case we use a discretized version of the momentum constraint (2.67) to update a. A typical profile of a(t, r) from a black-hole-forming computation is illustrated in Fig. 3.2.  3.2.3  Logarithmic Coordinates  To conclude this chapter we note that once we had implemented a code based on the finite difference equations described above, we found that for near-critical evolutions— and particularly for small l—we could not use sufficient resolution (small enough ∆r) to  28  3.2. Solving the Geometric Equations compute accurate solutions at a reasonable computational cost. We therefore introduced a new radial coordinate, ζ, defined by r → ζ = log(r + ǫ) ,  (3.28)  in which ǫ was chosen to be of order ǫ ≈ 10−2 rmax . After transforming our system of equations to the (t, ζ) coordinate system, we applied finite difference techniques on a mesh uniform in ζ (i.e. with constant ∆ζ). As will be seen from the results in Chap. 5, this approach was quite successful in providing increased resolution where it was needed, but without being prohibitively costly.  29  3.2. Solving the Geometric Equations  Figure 3.2: Value of the geometric function, a(t, r), as a function of radius, at an instant when the scalar fields are compressed close to the origin and have created a very strong gravitational field. The resulting spacetime contains a black hole with an (approximate) radius r = rBH as indicated by the dashed line. As discussed in text, in such scenarios the Hamiltonian solver often fails, due to the large gradients in a(t, r), and we use the momentum constraint to continue the numerical solution.  30  Chapter 4  Code Tests In this chapter we discuss various methods we have used to check the design and implementation of our numerical solver for the system of equations (2.61–2.73). We also explain some diagnostic tools that were employed to verify convergence of the solutions to the continuum limit. As we discussed in Chap. 2, and provided that a suitable choice of coordinates is made (and the polar-areal system is suitable), only two of the four independent components of Einstein’s equations are needed to determine the geometry of spacetime in spherical symmetry. Thus, provided that our numerical solution is correct, the two equations that are not used to compute the geometric variables should be automatically satisfied up to the accuracy of our numerical method. This yields a natural way to at least partially verify our solutions, which we call redundant equation evaluation, and our specific methods based on this technique are described in Sec. 3.1. We can also test out code by monitoring quantities that should be conserved in the continuum limit. In particular, we describe a test of the conservation of total mass (energy) from a typical code run in Sec. 3.2. Finally in Sec. 3.3, we present the results of a convergence test, in which solutions generated from fixed initial data but using different mesh scales are compared. The evidence from this test, combined with the other checks, suggests that our numerical solutions converge to their continuum counterparts at the expected second-order rate. We note that all of the tests performed in this chapter used the same initial data: namely two ingoing gaussian pulses, as described in the next chapter, for two scalar fields with l1 = 1 and l2 = 2, respectively, and where the pulses were chosen to overlap.  31  4.1. Redundant Equation Evaluation  4.1  Redundant Equation Evaluation  As we have discussed, our numerical approach is based on the solution of the t-t and r-r components of Einstein’s equations. This means that (2.67) and (2.69), which derive from the t-r and θ-θ components, respectively, must be satisfied if we are solving the t-t and r-r equations correctly. We thus define the following quantities Rt r = Gt r (a, α) − 8πT t r (ψ1 , ψ2 , a) ,  (4.1)  Rθ θ = Gθ θ (a, α) − 8πT θ θ (ψ1 , ψ2 , a) ,  (4.2)  which, for a numerical solution, should be residual in nature. That is, if we replace the right hand sides of the above with finite difference approximations which are, for example, second order in the mesh spacing, h, then as h → 0 we should find Rt r = O(h2 ) and Rθ θ = O(h2 ). Here · is any suitable discrete norm, such as the ℓ2 norm, ·  2,  defined  for any grid function, fi , by fi  2  =  N 2 i=1 fi  N  ,  (4.3)  where N is the total number of spatial grid points. In practice we coded redundant equation evaluators (REEs), as described above, independently of the main code. We note that compared to the solver itself, the REEs are relatively straightforward to develop, so are much less likely to suffer from implementation errors. Moreover, since the chance of having two independent codes which are in agreement (in a convergence sense), but which both have bugs is extremely unlikely, the REEs—combined with the convergence tests described at the end of this chapter—can be considered as constituting a robust test of our code. Results from a typical test using the REEs are shown in Fig. 4.1. Displayed is the time evolution of the l2 norm of the residual Rt r from a sequence of calculations that used the same initial data and mesh spacings h, h/2 and h/4. Convergence of the residuals is evident, and closer examination reveals that the convergence is second order, as expected. Fig. 4.2 shows the analogous plot for the residual Rθ θ , which again displays second order convergence as the resolution is increased. We therefore conclude that the original 32  4.2. Conservation of Energy solutions—based on the t-t and r-r components of Einstein’s equations—are such that the other two components are automatically satisfied in the continuum limit. As can be seen from both Fig. 4.1 and Fig. 4.2, the REE norms show more variation around time t ≈ 10, which is when the scalar fields are most centrally concentrated, and most influenced by the effective angular momentum barrier. During this period of time the evolution of the system is more dynamic, the scalar fields have larger gradients, and this results in the relatively large fluctuations. However, it should be stressed that there is still strong evidence for convergence of the REEs in this epoch, and the REE norms are, in fact, smaller than at earlier and later times. We note that, as will be seen in Sec. 4.3, analogous behaviour is also observed in the basic convergence tests we performed on the dynamical variables. However, it is not seen in our convergence test of the conserved mass, which is computed using the value of M (t, rmax ), and which is thus relatively insensitive to the interior dynamics.  4.2  Conservation of Energy  Since our spacetimes are asymptotically flat, the value of the ADM mass, MADM , defined by MADM = lim M (t, r) r→∞  (4.4)  must be a conserved quantity (see [1]). We can therefore monitor the value of the mass function at the outer boundary, r = rmax , of the computational domain as a test of our code’s accuracy. We note that M (t, rmax ) is, of course, only an approximation to the ADM mass, which is only defined at infinity. However, as long as the none of the scalar field matter has left the solution domain region, M (t, rmax ) should be conserved. Fig. 4.3 shows a plot of M (t, rmax ) as a function of time during a period when there is no scalar field leaving the computational domain. We observe that the value of M (t, rmax ) is roughly conserved, and that, more importantly, as the resolution increases we have convergence to conservation. This last fact is illustrated by the inset in the figure, which shows the deviation of M (t, rmax ) from its time-averaged value, also as a function of t: closer examination of the deviation shows that it is exhibiting second order convergence 33  4.3. Convergence Test to zero.  4.3  Convergence Test  In this section, we describe the convergence tests we have performed on the metric functions and the scalar fields. These tests suggest that all of the functions that are evolving in time are second order convergent as h → 0. Together with the fact that the redundant equation residuals are also converging as expected, this provides very strong evidence that our numerical solutions are actually converging to continuum solutions of (2.61-2.73) as desired. To perform a basic convergence test, one needs to run the code with the same initial data, but using three different mesh spacings, say h, h/2 and h/4. Let us denote our numerical solution for any of the evolving fields by χh (t, r), with a corresponding exact (continuum) solution, χ⋆ (t, x). We then have the following relation χh (t, r) = χ⋆ (t, r) + τ (t, x) ,  (4.5)  where τ is the solution error. Assuming that our finite difference discretization is second order, then following Richardson [22], we expect that τ can be written as: τ (t, r) = h2 e2 (t, r) + h4 e4 (t, r) + . . .  (4.6)  Here it is important to stress that the various functions e2 , e4 , etc. that appear in the above are h independent. We now define the convergence factor, Q[χ](t), associated with χ by Q[χ](t) =  ||χh − χh/2 ||2 . ||χh/2 − χh/4 ||2  (4.7)  Using (4.6), we can then easily see that the convergence factor should asymptote to 4 in the continuum limit: lim Q[χ](t) =  h→0  ||e1 (t)||2 |h2 − (h/2)2 | = 4. ||e1 (t)||2 | (h/2)2 − (h/4)2 |  (4.8)  Fig. 4.4 plots the value of the convergence factor as a function of time for the metric functions a and α, as well as one of the scalar fields ψ1 . The displayed results suggest that our numerical method is second order convergent for all of the variables. 34  4.3. Convergence Test  Figure 4.1: l2 norm of the residuals Rt r defined in (4.1) as a function of time, from a sequence of calculations using the same initial data, and mesh spacings h, h/2 and h/4. Convergence of the residuals is evident.  35  4.3. Convergence Test  Figure 4.2:  l2 norm of the residuals Rθ θ , defined by (4.2), versus time, and at three  distinct resolutions. The displayed data provides strong evidence that the θ-θ component of Einstein’s equations is satisfied as h → 0.  36  4.3. Convergence Test  Figure 4.3: Value of M (t, r) at r = rmax as a function of time for three different resolutions. The inset shows the deviation of M (t, rmax ) from its time-averaged value.  37  4.3. Convergence Test  Figure 4.4: Value of the convergence factor, Q(t), for three different field variables, a, α, and ψ1 , as a function of time. Convergence of all three quantities is apparent.  38  Chapter 5  Results In this chapter we present the results we have obtained. Sec. 5.1 summarizes calculations of single field critical collapse for the cases l = 1 and l = 2. These calculations are not original, but our results are in agreement with computations previously reported in [12]. In Sec. 5.2 we discuss the critical collapse of two scalar fields with angular momentum l1 = 1 and l2 = 2, focusing on the relative stability problem defined in Chap. 1. Sec. 5.3 is similarly devoted to the l1 = 1 and l2 = 3 case. Our main result is the evidence that fields with higher angular momenta are unstable in the context of a lower-l critical solution.  5.1  Critical Collapse: Single Field Case  In this section we summarize the results of our studies of critical collapse when the matter content is a single scalar field with angular momentum, l. We choose the following gaussian form for the initial scalar field profile: ψ(0, r; p) = p e  (r−r0 )2 δ2  ,  (5.1)  where p, r0 and δ are parameters. We note that this form thus describes a spherical shell of scalar field, centred at r = r0 , with a thickness that is controlled by δ. Here and in the following we will always use the overall amplitude p as the tuning parameter, keeping r0 and δ fixed as we tune to criticality. We choose the following values for the fixed parameters: r0 = 20.0 ,  (5.2)  δ = 3.0 ,  (5.3)  while the outer boundary is rmax = 50. To complete the specification of the initial data 39  5.1. Critical Collapse: Single Field Case we choose the time derivative of the scalar field via ∂ψ ∂t  = t=0  ∂ψ ∂r  ,  (5.4)  t=0  so that, to a good first approximation, the scalar pulse is purely ingoing at the initial time. Our search for a critical solution begins by determining values p< and p> such that for p = p< and p = p> , a black hole does not form or does form, respectively. We then tune p using a bisection search until we have determined a critical value, p⋆ , to about machine precision (relative error of ≈ 10−14 ). In studying our critical solutions—both here and in subsequent sections—useful diagnostic quantities are the mass gradient functions, dmi /dr, associated with the scalar fields. These are given by dmi (t, r) := r2 ρi = r2 nµ nν Tiµν dr  (i = 1, 2) ,  (5.5)  where i is the scalar field label, and ρi is the mass density associated with the i-th field. With respect to the mass (mass aspect) function, M (t, r), defined by (2.55), one can take the spatial derivative of (2.55) and then substitute the derivative a′ , computed from (2.66), in the result to find dM dm1 dm2 = + . dr dr dr  (5.6)  We can thus view dmi /dr, i = 1, 2 as encoding the mass content due to each of the corresponding fields: in particular we can use the mass gradients as a measure of the interaction and possible energy transfer between the two. Finally, we note that from (5.6) we have  ∞ 0  dm1 dm2 + dr dr  ∞  dr = 0  dM dr = MADM , dr  (5.7)  where, again, MADM is the conserved ADM mass. Fig. 5.1 shows the mass gradient for a near-critical, single scalar field calculation with angular momentum l = 1. As with the majority of the graphs in this chapter, the plot uses a (natural) logarithmic radial coordinate, and displays the mass gradient profile at a late time, T , when the system is still strongly gravitating on the smallest scale. The discrete self-similarity (echoing) associated with the critical solution is apparent, and from 40  5.2. Critical Collapse: l1 = 1, l2 = 2 the plot we have estimated the echoing exponent to be ∆l=1 = 0.50 ± 0.05. Within the level of the estimated error in our calculations, this result is in agreement with that of Olabarrieta et al [12] who cite ∆l=1 = 0.460 ± 0.002. The echoing behaviour can also be seen in the metric function a(T, r), which is displayed in Fig. 5.2. We performed a similar critical search for the single field case with l = 2. We again observed a discretely self-similar solution at criticality, with an estimated echoing exponent ∆l=2 = 0.15 ± 0.05. Once more, this result is in agreement with that of [12], which quotes ∆l=2 = 0.119 ± 0.003.  5.2  Critical Collapse: l1 = 1, l2 = 2  We now study a two field case with l1 = 1 and l2 = 2. The initial profiles for the scalar fields are given by ψ1 (0, r; p) = pe  (r−r01 )2 2 δ1  ψ2 (0, r; q) = qe  (r−r02 )2 2 δ2  ,  (5.8)  ,  (5.9)  where we choose r01 = 20.0 ,  (5.10)  δ1 = 3.0 ,  (5.11)  r02 = 22.0 ,  (5.12)  δ2 = 3.0 ,  (5.13)  so that there is overlap between the two fields at the initial time, as well as during the subsequent evolution. As before, we set the time derivatives of the scalar fields so that the initial pulses are ingoing. In addition, the parameters r01 , r02 , δ1 and δ2 are held fixed throughout our numerical experimentation, while the amplitude parameters p and q are varied. We now set the amplitude of ψ2 —which has angular momentum l2 = 2—to q = q¯ = 1.0 × 10−5 , and then tune the overall amplitude, p, of the first field (l1 = 1). We find a critical value p¯ ≈ 5.345 × 10−2 . For this critical parameter pair, namely (¯ q , p¯) = 41  5.2. Critical Collapse: l1 = 1, l2 = 2 (10−5 , 5.345 × 10−2 ), the ratio of the mass content of ψ2 to ψ1 is of the order 10−8 . Thus ψ2 is very weakly gravitating in this case. However, as shown in Fig. 5.3, examination of the late-time behaviour of the mass gradient of the field with higher angular momentum reveals a very interesting effect. As can be seen from the plot, the amplitude of the oscillations associated with the discretely self-similar solution for the l = 2 field—although small—are growing in time (right to left). This suggests that the l = 2 field is inherently unstable in the context of the l = 1 critical solution, analogously to the behaviour previously observed by Choptuik, where the massless scalar field was found to unstable in the context of the Type II Yang-Mills solution [13]. However, we note that in the current case the induced oscillations apparently have the same period, ∆, as that of the l = 2 solution. In Choptuik’s computations, if the scalar field is initially perturbative (very weak), the period of the excitation seems to be related to that of the Yang Mills solution. If we start with different initial amplitudes q = q¯, the observed behaviour persists and, in fact, as shown in Fig. 5.4, even in the case that the initial mass contributions of the two fields are comparable, the oscillations in ψ2 increase in time and eventually dominate the mass content of the spacetime (on the relatively small spatial scales on which the critical solution is still strong-field). In order to characterize the excitation of the oscillations of the l2 = 2 critical solution, we define Am i to be the maximum value of dmi /dr achieved over any given echo, where the index m enumerates the maxima (from right to left) in Fig. 5.3, and i = {1, 2} labels the two scalar fields. We then define the relative growth, ∆A2 /A2 , in the amplitude as Am+1 − Am ∆A2 := 2 m 2 . A2 A2  (5.14)  As can be seen from Fig. 5.3, the amplitude of oscillations grows in an exponential fashion. βm , where C and β are positive Therefore, by fitting a function of the form Am 2 = Ce  constants, we can determine the growth factor via ∆A2 Ceβ(m+1) − Ceβm = = eβ − 1 . A2 Ceβm  (5.15)  We then calculate ∆A2 /A2 by computing the average of β over a range of about 10 oscillations. The resulting values, which are compiled in Table 5.1, show that the relative 42  5.3. Critical Collapse, l1 = 1, l2 = 3 increase in the size of the oscillations is almost independent of the size of the initial data. This indicates that the excitation of the l2 = 2 critical solution would happen no matter how small the ψ2 field is initially. Combined with the fact that the log-period of the oscillation appears to be the same as that for the pure l = 2 critical solution, it is also plausible that the conjectured instability may be characterized by the same, single growing mode that characterizes that solution. Finally, in Fig. 5.5 we plot the critical curve as defined in Chap. 1 for the case l1 = 1 and l2 = 2. Here we have renormalized the tuning parameters so that they are the total mass contributions, mi , defined by ∞  mi (t, r) = 0  dm dr , dr  (5.16)  for each of the fields. We emphasize that for all cases except m2 = 0 (where the l = 2 field vanishes identically) our results suggest that the critical solution is the one associated with the l = 2 field.  5.3  Critical Collapse, l1 = 1, l2 = 3  We have also performed an analysis identical to that described above, but for the case l1 = 1 and l2 = 3 (except that we now use rmax = 100). Again we find evidence that the scalar field with larger angular momentum gets amplified in the context of the critical solution with lower l. Additionally, the relative growth of the mass gradient amplitude in one oscillation, ∆A2 /A2 , appears largely independent of the size of the initial data. Fig. 5.6 illustrates the growth in dm2 /dr for (¯ q , p¯) = (1.0 × 10−5 , 5.344−2 ), and the growth factor, ∆A2 /A2 is listed for a wide range of q¯ in Table 5.2. Finally, we ran a second set of numerical experiments, both for l1 = 1, l2 = 2 and l1 = 1, l2 = 3, but where the nature of the overlap of the two fields was changed. Specifically, we exchanged the centre-points of the two gaussian profiles by setting r01 = 22.0 and r02 = 20.0. Although this naturally results in different critical values (¯ q , p¯), we still observe that the field with higher angular momentum is excited in the presence of the lower-l critical solution. This provides some evidence that the effect we are seeing is not due to the particulars of the initial data we are using. 43  5.3. Critical Collapse, l1 = 1, l2 = 3  Table 5.1:  q¯  p¯  M2 /M1  ∆A2 /A2  10−7  5.345 × 10−2  7 × 10−13  0.10  10−6  5.345 × 10−2  7 × 10−11  0.11  10−5  5.345 × 10−2  7 × 10−9  0.12  10−4  5.345 × 10−2  7 × 10−7  0.12  10−3  5.344 × 10−2  7 × 10−5  0.11  10−2  5.261 × 10−2  7 × 10−3  0.08  2.0 × 10−2  4.977 × 10−2  3 × 10−1  0.05  Excitation of the l2 = 2 critical solution. This table lists various critical  parameter pairs (¯ q , p¯) (first and second columns), the ratio of total initial masses associated with the l1 = 1 and l2 = 2 fields (third column), and the estimated growth factor of the second field’s mass content (fourth column). Note that the control parameter is q¯, the overall initial amplitude of the second scalar field, ψ2 . We observe that especially for small values of the control parameter, the growth factor of ψ2 ’s mass contribution is largely amplitude-independent.  Table 5.2:  q¯  p¯  M2 /M1  ∆A2 /A2  10−6  5.344 × 10−2  1.1 × 10−9  0.15  10−5  5.344 × 10−2  1.1 × 10−7  0.14  10−4  5.344 × 10−2  1.1 × 10−5  0.15  3.0 × 10−4  5.344 × 10−2  9.5 × 10−5  0.15  10−3  5.344 × 10−2  1.1 × 10−3  0.14  10−2  5.308 × 10−2  1.1 × 10−1  0.04  Excitation of the l2 = 3 critical solution. Data columns are as in Table 1.  Again, the observed growth factor ∆A2 /A2 is almost independent of q¯ for a wide range of that parameter.  44  5.3. Critical Collapse, l1 = 1, l2 = 3  Figure 5.1: Mass gradient function, dm/dr(T, log(r +ǫ)), for a single field calculation with angular momentum, l = 1. Here T is a late time when the evolution is near-critical and the field is strongly self-gravitating on the smallest scale (left-most part of the plot), and ǫ is chosen to be e−1 for this computation. Here, and in rest of the plots in this chapter, we use the natural logarithm, (denoted by log) in the abscissa. The (discrete) self-similarity of the solution is reflected in the log-periodic oscillations of the mass gradient function.  45  5.3. Critical Collapse, l1 = 1, l2 = 3  Figure 5.2:  Metric function, a(T, log(r + ǫ)). in near-critical regime. Plot details are  analogous to Fig. 5.2. Again, the discrete self-similarity of the near-critical solution is evident, especially in the inset.  46  5.3. Critical Collapse, l1 = 1, l2 = 3  Figure 5.3:  Excitation of l2 = 2 field in presence of l1 = 1 critical solution. Plotted  is dm2 /dr(T, log(r + ǫ)), where T is again some late time in the critical regime, and ǫ is a small parameter: the growth in the function from right to left indicates a growth (instability) of ψ2 as a function of time in the context of the l1 = 1 critical solution. Moreover, the log-period (∆) of the oscillations is consistent with that of the l = 2 critical solution.  47  5.3. Critical Collapse, l1 = 1, l2 = 3  Figure 5.4:  This figure shows snapshots of the evolution of dmi /dr(t, log(r + ǫ)) from  a near-critical evolution with l1 = 1, l2 = 2, and for a case where the initial total mass content, mi of each of the two fields is comparable. Evolution proceeds left-to-right, topto-bottom. The plot provides further evidence for the dominance of the field with higher angular momentum at criticality.  48  5.3. Critical Collapse, l1 = 1, l2 = 3  Figure 5.5: Plot of the critical curve along which the system exhibits critical behaviour for two massless scalar fields with angular momentum l1 = 1 and l2 = 2. We conjecture that, apart from the point (m ¯ 2, m ¯ 1 ) = (0, m⋆1 ), the critical solution is the one associated with the l = 2 field.  49  5.3. Critical Collapse, l1 = 1, l2 = 3  Figure 5.6: Excitation of l2 = 3 field in presence of l1 = 1 critical solution. Plot details are as in Fig. 5.3. As before, this figure provides evidence for the amplification of the scalar field with higher angular momentum in the context of the critical solution associated with the lower-l field.  50  Chapter 6  Conclusions In this thesis we have studied the critical phenomena associated with the spherically symmetric collapse of two massless scalar fields whose equations of motion mock up some of the effects of angular momentum. Each field is characterized by a distinct angular momentum parameter, l, l = 0, 1, 2, . . ., distinct equations of motion, and distinct critical solution characteristics. We wrote and tested a numerical code to solve the system of PDEs describing the coupled Einstein/scalar equations, and then used the code to investigate critical collapse in several cases. Our main result is strong evidence that the type II critical solution of the scalar field with larger angular momentum gets amplified in the presence of the critical solution associated with the lower-l field. Thus, for our model, whenever the evolution involves two matter fields which physically overlap during critical evolution, we can expect the latetime critical solution (i.e. the solution on arbitrarily small spatial and temporal scales) to be identical to that corresponding to the higher-l field (and with the lower-l field being driven to 0 on those scales). Our results also lend support to a conjecture previously made by Choptuik [13] that the critical solutions associated with various types of matter may constitute some sort of a hierarchy. Clearly then, directions for future work could include the study of the relative stability phenomena for additional matter models (pairs, triples, etc.) interacting via the gravitational field. The possibility of nontrivial interaction between two or more types of matter in the presence of inherently unstable critical solutions could provide a rich phenomenology to investigate. In addition, the studies to date could be extended by relaxing the restriction to spherical symmetry, thus considering more generic spacetimes such as those with an axial symmetry, or even with no symmetry at all.  51  Chapter 6. Conclusions Another way in which the current work could be extended would involve incorporating adaptive mesh refinement (AMR) techniques into our finite difference code. This would allow us to better resolve the fine structure that develops in the type II critical solutions encountered in our model. Particularly noteworthy in this regard is the algorithm due to Berger and Oliger [23] in which the local finite difference resolution is modified dynamically through the monitoring of truncation error estimates. This method has previously been shown to work very well for the l = 0 case of our current model [2], so we expect that we could use it to significant advantage in our studies as well.  52  Bibliography [1] C. W. Misner, K. S. Thorne, and J. A. Wheeler. Gravitation. W.H. Freeman and Company, New York, (1973). [2] Matthew W. Choptuik. Universality and scaling in gravitational collapse of a massless scalar field. Phys. Rev. Lett., 70:9–12, 1993. [3] Carsten Gundlach.  Critical phenomena in gravitational collapse.  Phys. Rept.,  376:339–405, 2003. [4] Matthew W. Choptuik, Tadeusz Chmaj, and Piotr Bizon. Critical Behaviour in Gravitational Collapse of a Yang- Mills Field. Phys. Rev. Lett., 77:424–427, 1996. [5] Charles R. Evans and Jason S. Coleman. Critical phenomena and self-similarity in the gravitational collapse of radiation fluid. Phys. Rev. Lett., 72(12):1782–1785, 1994. [6] David W Neilsen and Matthew W Choptuik. Critical phenomena in perfect fluids. Classical and Quantum Gravity, 17(4):761, 2000. [7] Patrick R. Brady, Matthew W. Choptuik, Carsten Gundlach, and David W. Neilsen. Black hole threshold solutions in stiff fluid collapse. Class. Quant. Grav., 19:6359– 6376, 2002. [8] Steven L. Liebling, Eric W. Hirschmann, and James Isenberg. Critical phenomena in nonlinear sigma models. J. Math. Phys., 41:5691–5700, 2000. [9] A. M. Abrahams and C. R. Evans. Critical behavior and scaling in vacuum axisymmetric gravitational collapse. Phys. Rev. Lett., 70:2980–2983, 1993.  53  Bibliography [10] Matthew W. Choptuik, Eric W. Hirschmann, Steven L. Liebling, and Frans Pretorius. Critical collapse of the massless scalar field in axisymmetry. Phys. Rev., D68:044007, 2003. [11] Matthew W. Choptuik, Eric W. Hirschmann, Steven L. Liebling, and Frans Pretorius. Critical collapse of a complex scalar field with angular momentum. Phys. Rev. Lett., 93:131101, 2004. [12] Ignacio Olabarrieta, Jason F. Ventrella, Matthew W. Choptuik, and William G. Unruh. Critical Behavior in the Gravitational Collapse of a Scalar Field with Angular Momentum in Spherical Symmetry. Phys. Rev., D76:124014, 2007. [13] M. W. Choptuik.  Relative stability of black hole threshold solutions.  http://laplace.phas.ubc.ca/People/matt/Doc/Talks. [14] J. W. York, Jr. Kinematics and dynamics of general relativity. In L. Smarr, editor, Sources of Gravitational Radiation. Seattle, Cambridge University Press, (1979). [15] R. M. Wald. General Relativity. University of Chicago Press, Chicago IL, (1984). [16] M.  Choptuik  and  R.  L.  Marsa.  Rnpl  reference  manual.  http://laplace.physics.ubc.ca/People/matt/Rnpl/index.html. [17] Ignacio Olabarrieta and Matthew W. Choptuik. Critical phenomena at the threshold of black hole formation for collisionless matter in spherical symmetry. Phys. Rev., D65:024007, 2002. [18] Jason F. Ventrella and Matthew W. Choptuik. Critical phenomena in the Einsteinmassless-Dirac system. Phys. Rev., D68:044020, 2003. [19] M.  W.  Choptuik.  The  3+1  einstein  equations.  http://laplace/People/matt/Teaching/98Spring/Phys387N/Notes.html. [20] R. Arnowitt, S. Deser, and C. W. Misner. The dynamics of general relativity. In L. Witten, editor, Gravitation: An Introduction to Current Research. New York, Wiley, (1962). 54  Bibliography [21] Richard S. Varga. Matrix Iterative Analysis, 2nd Ed. Springer, Berlin, (2000). [22] L. F. Richardson. The approximate arithmetical solution by finite differences of physical problems involving differential equations, with an application to the stresses in a masonry dam. Philosophical Transactions of the Royal Society of London. Series A, 210:307–357, (1911). [23] Marsha J Berger and Joseph Oliger. Adaptive mesh refinement for hyperbolic partial differential equations. Journal of Computational Physics, 53(3):484 – 512, 1984.  55  

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.24.1-0071512/manifest

Comment

Related Items