Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Critical collapse of collisionless matter in spherical symmetry Olabarrieta, Ignacio (Iñaki) 2000

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

Item Metadata


831-ubc_2001-0090.pdf [ 2.31MB ]
JSON: 831-1.0089848.json
JSON-LD: 831-1.0089848-ld.json
RDF/XML (Pretty): 831-1.0089848-rdf.xml
RDF/JSON: 831-1.0089848-rdf.json
Turtle: 831-1.0089848-turtle.txt
N-Triples: 831-1.0089848-rdf-ntriples.txt
Original Record: 831-1.0089848-source.json
Full Text

Full Text

CRITICAL COLLAPSE OF COLLISIONLESS M A T T E R IN SPHERICAL SYMMETRY By Ignacio (Inaki) Olabarrieta Licenciado en Ciencias, Universidad del Pais Vasco, 1998 A THESIS S U B M I T T E D IN P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F M A S T E R O F S C I E N C E in T H E F A C U L T Y O F G R A D U A T E S T U D I E S D E P A R T M E N T O F PHYSICS A N D A S T R O N O M Y We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y O F BRITISH C O L U M B I A September 2000 _) Ignacio (Ihaki) Olabarrieta, 2000 In presenting this thesis in partial fulfillment 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 6224 Agricultural Road Vancouver, B.C. , Canada V6T 1Z1 Date: lt> ~ 12- 'XooG Abstract We perform a numerical study of the critical regime for the general relativistic collapse of collisionless matter in spherical symmetry. The evolution of the matter is given by the Vlasov equation (or Boltzmann equation) and the geometry by Einstein's equations. This system of coupled differential equations is solved using a particle-mesh (PM) method. This method approximates the distribution function which describes the matter in phase space with a set of particles moving along the characteristics of the Vlasov equation. The individual particles are allowed to have angular momentum different from zero but the total angular momentum has to be zero to retain spherical symmetry. In accord wih previous work by Rein, Rendall and Schaeffer, our results give some indications that the critical behaivour in this model is of Type I (the smallest black hole in each family has a finite mass). For the families of initial data that we have studied it seems that in the critical regime the solution is a static spacetime with non-zero radial momentum for the individual particles. We have also found evidence for scaling laws for the time that the critical solutions spend in the critical regime. ii Table of Contents Abstract ii Table of Contents iii List of Tables v List of Figures vi Acknowledgements vii 1 Introduction 1 1.1 Gravitation and Critical Phenomena 1 1.2 Collisionless Matter (A little bit of history) 3 2 The Equations 5 2.1 The 3+1 Formalism of General Relativity 7 2.1.1 3+1 Formalism in Spherical Symmetry 10 2.1.2 Maximal Areal Coordinate System 13 2.2 Stress-Energy Tensor 15 2.3 Evolution Equations 17 2.4 Conserved Quantities 19 3 The Code 22 3.1 Particle-Mesh Methods 22 3.2 The Field Equations 24 iii 3.3 The Evolution Equations 26 3.4 Interpolation and Restriction 27 3.5 Initial Set of particles 29 4 Checks 30 4.1 Mass Conservation 30 4.2 Einstein Clusters 32 4.3 Numerical Errors 34 5 Results 39 5.1 Is the critical solution Static ? 39 5.2 Time Scaling 44 5.3 Different Families of Initial Data 46 5.4 Conclusions and Future Work 49 Bibliography 51 A Summary of Equations in Maximal-Areal Coordinates 53 A.l Field equations 53 A. 2 Evolution Equations 54 B Direct Approach To Solutions of Einstein-Vlasov System 55 B. l Field Equations 55 B.2 Vlasov Equation 56 C Summary of Equations in Polar-Areal Gauge Coordinates 58 D Pseudo-code 60 iv List of Tables Values of time scaling exponent for the different families (t.s. stands for "almost time symmetric"). The errors above are assumed to be the same as that of the t.s. family List of Figures 3.1 Kernel for the first order interpolation (NGP) 27 3.2 Kernel for the second order interpolation (CIC) 28 4.1 Value of 2M(r, t)/r at r = r m a x as a function of time 31 4.2 Radial coordinate of 7 particles as a function of time for a cluster with gaussian energy density profile 33 4.3 Statistical Error 35 4.4 Truncation (dotted lines) and statistical error estimates (dashed lines) at t=0 36 4.5 Truncation (dotted lines) and statistical error estimates (dashed lines) at t=20. . 37 5.1 Snapshots of da/dt 40 5.2 Three ensembles of da/dt at t = 156 41 5.3 maxr (2M(t, r)/r) as a function of time 42 5.4 Shift function as a function of time 43 5.5 Evolution of Sr(t,r) 44 5.6 Scaling law for the time 45 5.7 Critical solution for different families of initial data 47 5.8 r 2 5 a a ( r ) for the different families during the critical regime 48 5.9 Scaling behaviour for different families of initial data 49 vi Acknowledgements I want to thank the rest of the graduate students under the supervision of Matt Choptuik, especially Jason Ventrella, with whom I have had many useful discussions. I also want to thank Bojan Losic, Noureddine Hambli, Luis Lehner and Bill Unruh who have been there for any question I ever had. I am grateful to all the people who have encouraged me from Spain, my family and friends, and the people of the department of Theoretical Physics and History of Science in the University of the Basque Country who allowed me to work there for some time during the development of this thesis. Finally, my supervisor, Matt Choptuik, for keeping on asking me where are the particles? and for being a friend in addition to a supervisor. vii Chapter 1 Introduction 1.1 Gravitation and Critical Phenomena Einstein's theory of gravitation connects the geometry of spacetime to its matter content via the field equations1: Gab = 8irTab. (1.1) Here Gat, is the Einstein tensor, whose coefficients are complicated functions of the metric, gab, and its first and second derivatives. Ta<, is the stress energy tensor, which depends on the matter, and also, in general, on the metric. In a given coordinate system, these equations are non-linear partial differential equa-tions for the metric coefficients and can give rise to very interesting solutions for the spacetime. Some of the most interesting phenomena occur at the threshold of black hole formation as was first discovered numerically by Choptuik [1], who studied the collapse of a massless scalar field in spherical symmetry. For initial data characterized by one parameter p (the maximum amplitude of the scalar field for example), Choptuik found a critical value p* such that for values of p, p > p*, the evolution gave rise to the formation of a black hole, while for p < p* the scalar field dispersed to infinity. Near p = p* the space-time approached a universal solution, independent of the particular choice of the initial shape of the pulse (e.g. gaussian, tanh, etc). It was found that the near-critical dynamics was characterized by self-similar will use units where c—1 and G = 1 (c being the speed of light, and G Newton's gravitational constant) 1 Chapter 1. Introduction 2 oscillations with scaling factor eA. Choptuik also found that the masses of the black hole formed obeyed a scaling law: MBHcx\p-p^ (1.2) where 7 was a universal exponent, again, universal in the sense of being independent of the choice of the family of initial data. This equation shows that, in principle, one can form black holes with arbitrarily small masses in the model. The transition to black hole formation in this case is thus continuous in the mass of the black hole. In analogy with the critical behavior in statistical mechanics this was called a Type II critical solution. A lot of work has been done in the last few years trying to find the critical solutions for different types of matter. Type I transitions (i.e. transitions where the mass of the smaller black hole is finite) have been also found [2]. In this case, the solutions that have been found are either static or periodic [3]. Here it was found that the time that the solution stays in the critical regime scales like: r ~ — cr In \p — p*\ (1.3) where, again, a is a universal constant2. The process of tuning p to p* can be understood as choosing initial data such that at the threshold of black hole formation, we tune-out any growing component. Since we can tune away the growing components using only one parameter, the solution apparently has exactly one unstable mode. Linear perturbations of the critical solutions [3] agree with that observation. In this thesis we study the critical regime for collisionless matter in spherical symme-try. If we think of the matter as being composed of individual particles, allowing these particles to have angular momentum different from zero (keeping the net angular momen-tum equal to zero to remain in spherical symmetry) can give us some insight concerning 2 Note that this constant depends on the units used. Chapter 1. Introduction 3 the importance of angular momentum in critical phenomena, without needing to go to spacetimes with less symmetry. 1.2 Collisionless Matter (A little bit of history) Collisionless matter (or dust) has been widely studied in general relativity. Einstein [4] used this type of matter to investigate the physical significance of the singularity at r = 2 M for the Schwarzschild solution. He studied clusters of particles rotating in circular orbits (and spherical symmetry) and proved that for clusters in equilibrium the Schwarzschild mass function cannot reach M(r) = r/2 (r being the areal coordinate), therefore these configurations cannot have an event horizon nor a singularity. He postu-lated that this result also would hold for systems not in equilibrium and concluded that in "physical reality" there is no configuration of matter such that M(r) = r/2 (i.e. black holes cannot form). Although he came to the wrong conclusion, he studied one of the first models for relativistic stars. In the case when the angular momentum of each particle is zero, and the spacetime is spherically symmetric, Oppenheimer and Snyder [5] showed that a homogeneous distri-bution of this kind of matter generally forms a singularity. Using comoving coordinates they found the solution for the spacetime and deduced that singularity formation would happen in finite proper time (at least from the point of view of the comoving observers). Tolman [6] and Bondi [7] generalized this result for non-homogeneous configurations, restricting to the case where no particle overtakes any other during the evolution. Eardley and Smarr [8] considered the same model as Tolman and Bondi and studied the existence of maximal (K = 0) and constant mean curvature (K = KQ) slicing condi-tions. Shapiro and Teukolsky [9], [10], [11], [12], [13] have studied this system in detail in the context of models of stellar clusters. They developed a code very similar to ours Chapter 1. Introduction 4 to study properties of relativistic clusters, such as stability and collapse to a black hole. They also extended the code to treat more general spacetimes , i.e. with axial symmetry, [14] to study naked singularities. Kendall [16], [17], and Rein [15] have a series of papers addressing general properties for the Einstein-Vlasov system (the Cauchy problem, existence of static solutions in spherical symmetry, existence of general asymptotically flat solutions, etc.). More recently, Rendall, Rein and Schaeffer [18] have published the first study of critical behavior (black-hole threshold behavior) for this type of matter. Indeed, in this last paper Rendall et al argue that the critical solution in the case of collisionless matter is generically Type I. Our results agree with theirs in that we also find evidence for Type I transitions for certain initial data families which we have studied in this thesis. In addition to that result, we argue here that the critical solutions in this model are static, and we present evidence for scaling laws of the form (1.3). Finally we have preliminary investigations aimed at determining whether the critical solution is unique (universality). Chapter 2 The Equations The dynamical state of dust can be described by a distribution function, / : f(xa,Pa) = dN/dVp (2.1) where N is the number of particles and Vp is the volume in phase space. In our case, and since the matter is collisionless, the volume in phase space is a conserved quantity during the evolution of the system (Liouville theorem). This implies that the distribution function is also a conserved quantity: This is the collisionless Boltzmann or Vlasov equation. This equation plus Einstein's equations (1.1), all restricted to spherical symmetry, i.e.: form the system of equations that we want to solve numerically. Numerical solutions of PDEs generally involves discretization of the continuum prob-lem, and here we can consider at least two approaches: 1) Consider a set of particles which approximates the distribution function at the initial time and evolve the particles along the characteristics of equation (2.2). The geometry is computed at each time step using the energy densities derived from the positions and velocities of the particles. 2) Discretize the Vlasov equation for the distribution function in phase space and solve for (2-2) f(t,x\Pi) = f(t,Rx\RPi) with Re SO{3) i = 1,2,3 (2.3) 5 Chapter 2. The Equations 6 it directly in phase space. In this case the energy densities, required in the update of the geometry, are calculated by integrating the distribution function over momentum space. In the particle approach, we have the problem that a given set of particles is just one of the infinite number of possible realizations for given density and velocity profiles, therefore we introduce statistical errors. In the second case we have to solve the Vlasov equation, a partial differential equation in phase space. Although we have taken the particle approach due to the simplicity of the evolution equations, we include the equations for the direct integration in Appendix B. As with any problem in "3+1" (or ADM) numerical relativity, we want to be able to specify initial data on a spatial hypersurface and then evolve these data in time. To do this, we need to split the equations (1.1) into a set of constraint equations (equations that must be satisfied at each instant of time) and dynamical or evolution equations (equations that tell us how to evolve the geometric quantities in time). In the most general case, we will have four constraint and six second-order-in-time evolution equations. In order to do this splitting, we will make use of the 3 + 1 formalism due to Arnowitt, Deser and Misner (ADM), for a review of this formalism see [19]. We have followed a summary by Choptuik [20] which is itself based on a review by York [24]. In the remainder of this chapter, we derive the system of equations that we solve. We start with a summary of the 3+1 formalism and its restriction to the spherical symmetric case, and a definition and description of the maximal-areal coordinate system. In section 2.2 we explain how the Einstein's equations are coupled to the matter, and in 2.3 we explain the particle evolution equations in our coordinates. The chapter ends with some calculations of conserved quantities. Chapter 2. The Equations 7 2.1 The 3+1 Formalism of General Relativity. We consider a spacetime manifold with metric gab defined on it, and assume that we can slice the manifold into spacelike hypersurfaces, defined by t = constant, where t is our time coordinate. At least locally then, the space-like surfaces can be described in terms of a closed one-form, f^1, which is just the gradient of t na = dt = Vat. (2.4) The norm of this one-form is: gabnaQb = -a~2, (2.5) where a is a scalar function commonly known as the lapse function. As long as our slices are spacelike, a will be strictly positive. We introduce a normalized one-form, na, defined by na = -aQ,a - -aVat. (2.6) The contravariant vector, na = gabnb, which is the future-directed normal to the surfaces, can be viewed as the four velocity field of observers moving orthogonally to the slices. In this section, again following York, we will describe the 3+1 split in terms of 4-dimensional tensors. At the same time, we will need to decompose various tensors (in-cluding the Einstein and stress energy tensors appearing in (1.1)) into pieces parallel to na ("timelike" pieces) and pieces orthogonal to na ("spacelike" pieces). Thus we define Wn = -Wana, (2.7) where the sign convention is again York's, while for covariant vectors we define Wn = Wana. (2.8) 1 Indexes in this section are abstract indexes as in [26] Chapter 2. The Equations 8 For the projection onto the hypersurfaces we use the projection operator: 7a6 = 6% + nanb (2.9) The projection of the space-time metric to the hypersurfaces is: 7a6 = l°aldb 9cd = 9ab + "a « 6 (2-10) and is, in fact, the metric induced on the hypersurfaces by the 3+1 splitting. (Note that 7 a 6 is not necessarily the inverse of jab, as follows immediately from the definition of jab, and in this section we raise and lower all indexes with the four dimensional metric gab and its inverse). We can now introduce the natural derivative operator on the space-like surfaces, by projection of the four dimensional covariant derivative: £»a = 7"aV 6 . (2.11) It is easily shown that this derivative is compatible with the spatial metric: D a l b c = 0. (2.12) Having defined a derivative operator on the hypersurfaces, we can compute the asso-ciated Riemann tensor 3Habcd which measures the intrinsic curvature of the hypersurface. For example, for a spatial one form uia (uana — 0) we have (DaDb - DbDa) uc = 3Kabcdud. (2.13) We can then also define the spatial Ricci tensor, 3TZac = 3TZabCb and the spatial Ricci scalar, 3Tl = 3Tlaa. Finally, we define the extrinsic curvature tensor, Kab: Kab ' -Danb = V Q n 6 - naab, (2.14) Chapter 2. The Equations 9 where ab = naDanb is the dual to the four acceleration field of the observers moving with the slicing. The geometry of the space-time is completely defined by the spatial metric, jab, which describes the geometry of each slice, and the extrinsic curvature tensor which tells us how each slice is embedded in the four dimensional space-time. We can write Einstein's equations in terms of the tensors defined above, as well as the following projections of the stress energy tensor: p = Tabnanb (2.15) ja = 7 V T V ) (2-16) Sab = 7 a c 7 6 d r c d (2.17) These are interpreted as the local energy density, momentum density and spatial stress tensor, respectively, for an observer moving orthogonally to the space-like hypersurfaces. The Hamiltonian constraint is found by contracting Einstein's equations with na twice, Gabnanb = Tabnanb, yielding 31l + K 2 - K \ K b a = 16?rp (2.18) where K — gab Kab = 7 a 6 Kab is the trace of the extrinsic curvature and 3H is the spatial Ricci scalar defined previously. If we contract Einstein's equations once with na, and then project the resulting vector onto the hypersurface, jacGcbnb = ^acTcbnb, we get the momentum constraint: DbKab - DaK = 87rf. (2.19) These two constraint equations involve only spatial tensors and spatial derivatives of spatial tensors, and must be satisfied on each spacelike slice, including the initial slice. Chapter 2. The Equations 10 In order to derive the evolution equations we use the Lie derivative along the vector field, Na: Na = ana + Pa = d/dt (2.20) where /3a is an arbitrary spatial vector field, commonly known as the shift vector. Note that Na satisfies: NaQa = 1 (2.21) It can be shown that Lie differentiation along Na commutes with the projection operator. We can now write the rest of equations (1.1) as two different sets: (1) the evolution of the spatial metric: Caab = -2aKab + C0jab (2.22) which can also be viewed as a definition of Kab, and (2) the evolution equations for the extrinsic curvature: CtK\ = CPK\ - DaDba + a3Kab + KK\ + 8TT {^P\ (S - o) - 5 a 6 ) (2.23) where S is the trace of tensor given by equation (2.17). 2.1.1 3+1 Formalism in Spherical Symmetry We now restrict our attention to spherical symmetry. In contrast to the previous section, where our tensor expressions involved abstract 4-dimensional indices, we will be mostly concerned with the components of specific 3-tensors in this section. Thus, Latin indices i,j,k,- • • range over the spatial values 1,2 and 3, and all indices of the spatial tensors are raised and lowered with 7 U and 7;J respectively, where 7tJ' jjk — Slk-The most general 3-metric in spherical symmetry can be written as: V = a2 (t, r) dr2 + r 2 b2 (t, r) dtf (2.24) Chapter 2. The Equations 11 where dQ,2 is the metric on the unit sphere, dQ2 = d92 + sin2 9 d(f>2. From the 3-metric, we can compute the associated connection coefficients or Christoffel symbols: ^ jk — g') ' ( 7 n j , f c + 7 n f c , j ljk,n) In our coordinates, the non-zero connection coefficients are: (2.25) T r r r = a!/a rree = - ((r262)') / (2a2) T*H = - sin 2 9 (r2b2)' / (2a2) Terg = ((r262)') / (2r262) r % = - sin 9 cos 6 = (r 26 2)'/(2r 26 2) r % = cot0 We also need to compute the components of the Ricci tensor: 3 n pn pn I pn nm pn pm " i j — i ij,n ^ in,j > L nm ij 1 jm 1 in from which we find the non-zero components: 3 R r r = - [(r2b2)' I (r2b2)}' + a' (r2b2)' / (ar2b2) - \ [(r2b2)' / (r2b2)\' * R e e = - [(r2b2)' / (2a2)]' + 1 - [a' (r2b2)' / (2a3)] (2.26) (2.27) (2.28) 3R^ = sin2 9 Ree The mixed components are given by 3R% = -2 [(rb)'/a]'/ (arb) lRe e — R^s — a - (rb{rb)' /a)'] / (a(rb)2) (2.29) Chapter 2. The Equations 12 Another quantity which we need is , where "|" denotes covariant differentiation with respect to the 3-metric. <x\}j = l j k { a i k - T n i k a j n ) (2.30) which has the following non-zero components: = {a'I a)' /a a , / = a' (r2b2)' / (2r W ) a , / = (rb)' a'/ (rba2) (2.31) Finally, the Ricci scalar, 3R — ^ PHJ, is given by: (2.32) , 2 3R = .2fMY + 4 ,_M \ a I rb \ az At this point we have calculated all the geometric objects intrinsic to the spatial slices that we need. Now, we embed these slices into the most general spherically symmetric spacetime whose metric can be written: ds2 = (-a2 + a2 ft2) dt2 + 2a2ftdtdr + a2dr2 + r2b2d£l2 (2.33) where a, ft, a and b are all functions of r and t. Note that due to the spherical symmetry, the shift vector has only a radial component; ft1 = (ft,0,0), so ft will be called the shift function. In these coordinates n M = (—a, 0,0,0) and n M = (1/ct, —ft/a,0,0). The extrinsic curvature can be written in terms of the 4-metric coefficients: From this we compute the non-vanishing components for the current spherically sym-metric case: KrT = ^ ( ( « / ? ) ' - d ) (2.35) Kee = ^(P(rb)'-rb) (2.36) KH = s i n 2 ^ 0 (2.37) Chapter 2. The Equations 13 The corresponding components with mixed indexes are: K \ = - L ( ( a / ? ) ' - d ) (2.38) aa v ' K \ = ^ (0 (rb)' - rb) (2.39) K% = K \ (2.40) and the trace of the extrinsic curvature is Finally, the only non-zero component of DiK = Ku is KV = K,r = - — (- + (a/3' - d) + — ((a/3)" - d') - \ c l (ft (r6)' _ ^ rb 2.1.2 Maximal Areal Coordinate System. Einstein's equations allow for coordinate freedom that we need to fix. We have chosen maximal-areal coordinates (defined below), but have also included the equations for the polar-areal system in Appendix C. The main advantage of maximal-areal coordinates is that, in contrast to the polar-areal case, the slices used can penetrate apparent horizons. As the name suggests, in the maximal-areal coordinate system, the radial coordinate is areal, so that the proper area of 2-spheres with radius r is 47rr2. In terms of the general spherically-symmetric 4-metric (2.33), this means that b(r,t) = 1. The time coordinate is fixed by demanding that the 3-slices be maximal, i.e. that K (r, t) — 0. This leads to a condition on the lapse function a (r, t) (the so-called slicing condition), which must be satisfied at each instant of time. Thus, in maximal-areal coordinates the 4-metric (2.33) takes the form: ds2 = ( - a 2 + a2p2) dt2 + 2a2pdtdr + a2dr2 + r2dtf (2.42) Chapter 2. The Equations 14 Using the expressions we derived in the previous section, we compute the relevant geo-metrical quantities we need in order to write Einstein's equations in terms of these metric coefficients: R = 2 [-2 (1/a)' + a/r (a - 1/a2)] /(ar) R\ = -2 ( l /o ) ' / (ar) R% = B** = [a - (r/a)'] / (ar 2 ) K% = ((a B)' -a) /(a a) K \ = K% = 8/ (r a) (2-43) K\\i = Krr,r + 2 (K% - Kee) /r K]r = 0 a\/ = (a1/a)' /a a , / = ct\J* = a'/ (r a2) We can now specialize equations (1.1): Hamiltonian Constraint: a' 3 \a2rKee2 + Airra2p + ^ ( l - a2) (2.44) a 2 2r v ' Momentum Constraint: Kee = - - K d e - 4 n j r (2.45) r Evolution of the 3-metric: a = 2aaKee + (aB)' (2.46) B = arK°e (2.47) Slicing Condition: a» = a' l°L _ 2\ + 2 « + fl2 _ ^ + W a ( 5 _ 3 p ) ( 2 4 g ) \a rJ rl \ a ) Chapter 2. The Equations 15 where in deriving the last formula, we have made use of the slicing condition (K = 0), and the evolution equations for the extrinsic curvature components, equations (2.23). As we will discuss in more detail below, we have chosen to do a constrained evolution, which, in this case means that we use the constraint equations, rather than evolution equations, to update a and Keg. 2.2 Stress-Energy Tensor In this section we explain how we calculate the energy densities p, jr and S that appear in equations (2.44-2.48). We approximate the distribution function (2.1) by a set of N "spherical particles". Since the particles only interact with each other gravitationally, we have N T^ = J2 Tr (2.49) i=l where T/ i l / is the stress energy tensor for a single particle. For a point particle we can write: 7 T = —<Kr - r i ( i ) ) . (2.50) rrii Here, pf is the p component of the 4 momentum of the z-th particle, rm is its rest mass, and fi(t) is the spatial position of the particle at time t. Using equations (2.15-2.17) for p, S and jT specialized to our coordinates, we get: \jr]t = a [T*r]. , (2.53) Now we relax the point particle approximation and assume that each particle is a spher-ically symmetric shell of mass, uniformly distributed over a region A r of space2. These 2 Later on this A r will be the mesh spacing used in the finite difference solution of the geometrical equations. Chapter 2. The Equations 16 shells of matter are averages of the shells at that point with magnitude of the angular momentum I pointing in all posible directions. Therefore the angular momentum of the average is zero, I = 0, but I2 ^ 0. The proper volume that each particle occupies is then: t 2 t Vi = ^ - A r [ ^=g~ d(f)d6 = 4irAraa7-tPi (2.54) rrii J rrii and we can approximate the delta function that appears in (2.50) by 1 /V^. This yields: ^ AirAra r2 ^ ^ Si = + {S\)t = 4 f f A r 0 3 [ ^ j 1 ' r ? + 4)rArar/^'] j ( 2 ' 5 6 ) w . = s L W 1 ( 2 ' 5 7 ) Here a, a and B are evaluated at r = r-j, [p^ is defined as [p4]- = ^[p^, and [I]? = \pe\i* + [p</.]j2/ sin2 6i is the square of the magnitude of the angular momentum of the i-th particle. We then introduce quantities which do not explicitly depend on the geometry (since after updating the particle positions we want to solve for the geometry): [p]i = a [p]t, [S\\. = a 3 [SMi, [S\]. = a [S a a], and = a [jr]v We interpolate the one-particle quantities to the continuum and sum over all the particles to find the total values: Np f = T,fW(r-n) (2.58) i=i where fi is any of the barred quantities defined above for each particle, / is the corre-sponding quantity in the continuum case, and W(r — rt) is an interpolation function (see section 3.4 for further explanation of how we define this function). Having defined (2.58) we can now write equations (2.44-2.48) as a' 1 — a2 3 , TS6 2 , , / o r n \ — = — h -ralK\ + 4?rarp (2.59) a 2r 2 Chapter 2. The Equations 17 Kee' = A K o g _ 4 7 r t (2.60) r a a» = a > ( * _ l ) + ^ ( a 2 _ l + 2r-)+4naa(^ + S\-3p) (2.61) a r J r2 \ a J \ a2 0 = arK\ (2.62) 2.3 Evolution Equations The equations of motion for the particles are just the geodesic equations (the character-istics of the Vlasov equation). These can be expressed in terms of the four momentum of the particle: paVapb = 0, (2.63) so in our coordinate system we have: % = -r\t(pt)2-2T\Tprpt-rKT{ff-rt96^ (2.64) ^ = -rtt (pf - 2 P „ P y - rrr {prf - r M £ (2.65) ^ = - 2 r ^ P y - r V / (2 .66 ) where the T's are the Christoffel symbols defined by (2.25), and I2 = pg2 +p^,2/ sin 2 6 is the square of the magnitude of the angular momentum as before. We want to consider the 4-momentum as a function of our time coordinate and not the proper time of each particle. In order to do the transformation we can use the chain rule (d/dr = (dt/dr) (d/dt) and pl = dt/dr): % = - r \ t p t - 2 T \ r f - r t r r ^ - r t 0 e 4 i (2-67) dt pl plr* drf (r>r)2 I2 ^ = - r a P ' - 2 r l r p ' - r „ ^ J - - r m — t ( 2 . 68 ) « £ = _ 2 r » r , ^ - r » « < (2.69) dt pl pl Chapter 2. The Equations 18 Furthermore it is convenient to express these equations in terms of pr rather than pr: Pr = 9r^ = a2(t, r)p(t, r)pl + a2(t, r)pr (2.70) Solving for pr: p r { t ) = a J M - ^ r ) p t { t ) ( 2 ' 7 1 ) We can now compute the total derivatives with respect to coordinate time as follows: A = 1. + _ _ _ _ _ _ = j _ , P___ / 2 72) dt dt drdtdr dt pl dr { ' ' Applying this operator to equation (2.71) we get: d£_ = \_dpr_ _2f_(da fda\ _ M _ (dp f_dp\ dt a2 dt a*\dt ptdr) P dt P \dt p< dr J { ' ' Substituting equations (2.71) and (2.73) into equation (2.68) we obtain: dpr da t dp 1 dap2 I2 . dt or dr aA dr pl plr6 which is the evolution equation for pr. To derive the evolution equation for r we use the definition of pT (pT = dr/dr) which after some manipulation yields: dr pr P (2.75) dt a2pt To find the time component of the 4 momentum we make use of the normalization condition p^p^ = —m2: <** = f ^ f ^ V 2 ( 2 - 7 6 ) It is also convenient, as previously mentioned, to use p* = apl rather than pl itself. Then the resulting geodesic equations in these coordinates are (we have included equations for p ,^ p* and dp^/dt because use we use them sometimes for visualization purposes): = _ ^ + _^ + ______; + ___ ( 2 7 7 ) dt dr dr r a3 dr p* p*r3 Chapter 2. The Equations 19 dr ~di ~ Pl = dp0 2 / apr dt ~ r yaPp1 dd ape dt ~ pt <t>2 1 (I2 sin2 6\r4 d(f> ap* ~dt pt (2.78) + ^ (2-7 9) (2.81) P92) (2-82) (2.83) 2.4 Conserved Quantities Here we include the calculation of two quantities that we have used to check the code. The first quantity that we discuss is a conserved quantity in flat spacetime related to the energy of the particles. This gives us a check of the evolution when we fix the geometry to be Minkowskian (flat), which was useful in the development of the code. The second quantity that we calculate is a mass aspect function that coincides with the A D M mass at infinity. The conservation of the value of this function at the outer limit of our range of integration gives another check in the case when we solve the fully coupled problem. In flat space-time, the geometry is static, and therefore we have a time-like Killing vector field. Choosing coordinates such that: ds2 = -dt2 + dr2 + r2 d92 + r2 sin 20d0 2, (2.84) we can write this Killing vector field as £" = d/dt. If we contract £„ with the stress energy tensor T^, we get a vector field whose divergence is zero: V M (T"%) - ( V M T H + T"" ( V ^ „ ) = 0. (2.85) Chapter 2. The Equations 20 The first term on the right hand side is zero by the conservation of the stress-energy tensor, and the second is zero by the assumption that f satisfies the Killing equation3, and by the symmetry of T 7"'. Integrating this divergence over the four volume we get: / V^T^^^dd'x = [ T^f^JWg d3x = 0 (2.86) Jv JdV where the second integral is an integral over the 3-hypersurface dV (boundary of V). We choose this 3-hypersurface to be composed of two different hypersurfaces of constant t, t = ti and t = t2, joined by a timelike-hypersurface at spatial infinity. We assume that the stress-energy tensor vanishes at spatial infinity, so that it does not contribute to the integral. Let us also assume that n M = (—1,0,0,0) for the surface with ti, and nu = (1,0,0,0) for t2. In this case: f TuJWgd3x = f TttJ(Vgd3x (2.87) JdVH V JdVt2 V and since t\ and t2 are arbitrary this means that the integral is constant for any value of the time and the integrand calculated using (2.50) C = T^- (2.88) is a conserved quantity 4 . The other quantity that we calculate in this section is the mass aspect function: (2.89) To derive (2.89), we consider the transformation relating the metric (2.42) to the Schwarzschild metric (in a region of vacuum): 2M\ , , / 2M\ - l dsl = _ ( i _ — ) dtl + ( i _ — j dr] + r2sdtf (2.90) 3 V M 6 + V ^ M = 0 4 N o t only (2.88) is a conserved quantity but, in flat spacetime, pl for each individual particle is conserved. Chapter 2. The Equations 21 Here the V subscript stands for Schwarzschild and we wish to stress the point that the Schwarzschild time and our time coordinate are different, while the radial coordinates are the same. Following [9] we get: dt, = C{dt + Ddr) (2.91) drs = dr (2.92) where the last equation holds because both radial coordinates are areal. Using these formulas in (2.90) we get: d$l = - ( l - Cfdf - 2 ( l - C M * _ (i _ V±\ C2DV + ( l - ^ . y 1 d r 2 + r W (2.93) It is easy to see that in order for expressions (2.42) and (2.93) to be equal we need: 2 = _ - a 2 + a2p2 . ^ l-2M/r V D = fl2/3 0 0 (2.95) If we now compare the dr2 terms of (2.42) and (2.93), we get an equation for the mass: fl4/?2 • 1 - a 2 , (2.96) -a2 + a2B2 ' \-2M/r which after some manipulation can be rewritten as: 2M 1 B2 , n . = 1 - ~2 + ~2 • 2 - 9 7 This last equation is interesting in its own, because when this quantity, 2M(t,r)/r, becomes equal to 1, we know that a marginally trapped surface has been formed. We can see this by computing the expansion of the null geodesies (see for instance [25]), which in this coordinates can be written as: l-.(«,r)r^.(«,r) = !-*<«• yt«'r>. (2.98) cn(t, r) Therefore, if the expansion is zero, 1/a2 = B2/a2, and 2M(t,r)/r = 1. Chapter 3 The Code In this chapter we explain the details of the method we have used to solve the equations which were presented in the last chapter. 3.1 Particle-Mesh Methods We have used a method of calculation known as "particle-mesh" (PM). For a very com-plete review of this and other particle methods see [21]. This technique solves the coupled Einstein-Vlasov equations by splitting each time step into two stages: 1) solution of the field equations on a mesh, and 2) updating of the positions of the particles via their equations of motion. To provide some motivation for the use of P M methods let us first consider the N body problem in Newtonian gravity (where we assume N is very large). We could calculate the forces experienced by each particle Fj directly using Once the acceleration at each particle is known, we could integrate the equations of motion over some small time interval At giving the new positions and velocities of the particles. This method is denoted particle-particle, and computationally is very expensive since we have to perform 0(N2) calculations per time step. The "particle-mesh" (PM) approach to the problem involves calculation of the mass density due to the ./V particles on a mesh (assuming each particle occupies a finite region 22 Chapter 3. The Code 23 in space). On this mesh, we can use finite difference approximation (or some other discretization technique, like a spectral method) to solve the Poisson equation for the gravitational potential: V 2 $( t , r ) = 4nGp{t,r). (3.2) Once the potential is found, we can interpolate the value of its derivative to the actual particle positions to find the forces. A time step is then concluded by solving the evolution equations for each particle to provide updated positions and velocities. P M methods are most effective when the close interactions between the particles are not important for the evolution of the system. The main reason for this is that smoothing of the particles in a finite region of the space produce a large error in the density close to the particle positions. However, in cases when the dynamics is dominated by a mean field the particle-mesh method is much faster than direct integrations like particle-particle. In our case we want to solve the Vlasov equation (equation (2.2)) by approximating the distribution function with a set of point particles. Further approximation is necessary to couple to the geometric equations, since a point particle gives rise to infinite densities. To solve the geometric equations we thus have to find the stress-energy tensor for the particles, where we consider each particle to be smoothed out over a finite region of space. That is precisely the P M method defined above. In the next sections of this chapter we will explain the specifics of the integration of the field equations, the evolution equations for each particle and the interpolation scheme that we have used. We also include pseudo-code for the main routine of our program in Appendix D. Chapter 3. The Code 24 3.2 The Field Equations We first explain how the equations (2.59-2.62) for the geometry are solved numerically, assuming we know the quantities p, J r , S£, Sg 1 . The first two equations, equations (2.59-2.60), are integrated from the origin, r = 0, using the LSODA [22] integrator. The boundary conditions are given by the spherical symmetry of the space-time, and by the demand that the spacetime be locally flat at r = 0. They are a(t, 0) = 1, and Kee(t,0) = 0. We compute the values of the functions <2j and K6Qi on a uniform grid of Nr points at positions r*: r 0 = 0, n = h, r2 = 2h,..., rNr = r m a x . In order to compute the values at r = ri+x, we supply to LSODA the values of the functions at r = and the derivatives computed using equations (2.59-2.60) at r = ri+i/2, using p and J r averaged between the i-th point and the (i-t-l)-th point: [ P W = \{[p]i + \PUi) (3-3) IJrli+l/a = \ (&r]i + (3-4) Once we have calculated a we can solve the slicing equation: n1 9 2rv n' a" = a'(- - -) + ^ ( a 2 - 1 + 2 r - ) + 47raa(5rr + Saa - 3p), (3.5) a r r2 a with the boundary conditions: a'(t,0) = 0 (3.6) a(t,oo) = 1. (3.7) Here the first condition follows from (3.5) demanding that a be regular at r = 0 and the second one follows from asymptotic flatness, and the demand that t measure proper 1we will explain how to compute these quantities in section 3.4 Chapter 3. The Code 25 time at infinity. We use a second-order finite difference approximation on the previously defined mesh of points {r;}: ai+i - ai_i fai+i - aj-i 2 ai+i — 2a>i + ttj-i h~2 2h 2ha{ f # l - r - i + + Rearranging this equation gives us: where: fi = ^ + 2 | ) a ' - 1 - ( ^ + f t ) a i + U - 4 ) a < - 1 = 0 2hdi ft = ^ f ^ - l + 2 n ^ t i ) + 4 , 0 , [ I ^ + [ 5 ^ - 3 ^ In addition to these equations we have the boundary equation at r = 0: l / ^ + / 2 / ( 2 / Q \ / - 2 A 2 - P 2 \ " D + I/ /* 2 - / 2/(2/i) j 0 1 + ^ + l/h? - h/(2h)) C * 2 = 0 (3.8) (3.9) (3.10) (3.H) (3.12) (3.13) (3.14) which can be derived from the forward finite difference approximation to a' = 0 at r = 0: —3ai + 4o!2 — CH3 2h 0, (3.15) and equation (3.11) with i = 2. We have also the boundary condition at r = rm a x: OiNr = (2M^ K (3.16) Nr This is an approximation to (3.7) at a finite value of r = rm a x, where M is the mass aspect function defined by equation (2.89). These equations form a linear system of algebraic equations that can be solved using a tridiagonal solver (we have used the L A P A C K [23] routine dgtsv). Chapter 3. The Code 26 3.3 The Evolution Equations To evolve the particles' positions and momenta we integrate the geodesic equations (2.77-2.78). The values of the coefficients in these equations (basically products and quotients of a, a, 8, a', a' and 8') must be calculated at the particles' positions using the values obtained at the mesh points {rj}. The mesh values are interpolated to the particle positions using the same operator kerner used to produce mesh values from particle quantities (this procedure is explained in section 3.4). These equations are integrated using the LSODA integrator as in the case of the Hamiltonian and momentum constraints (2.59-2.60). For a given particle at position rn and momentum p™ at time step n, we calculate the new position rn+l and momentum p™+ 1 at t = tn+l by suplying to LSODA the values of r", pn and their derivatives at time step t = tn. This procedure has an accuracy O(At) (since the derivatives are calculated using the metric coeficients at t = tn) where At = tn+l — tn. In our program we chose a value of At proportional to h, i.e. At = A h, where usually A = 1.0. We need to take special care if a particle leaves the range of integration (r > rm a x) or if it reaches the origin. In the first case we stop considering the particle and we change the total number of particles to (N — 1). When a particle p reaches the origin, in other words when rp < 0, we reflect the particle, i.e. we set: [r]p = • - M , (3.17) H = • - fa], (3-18) WP => Mp (3-19) As mentioned previously, equations (2.82), for p^, and (2.83), for <fi, are also sometimes solved for visualization purposes and is also integrated in the same way. Chapter 3. The Code 27 3.4 Interpolation and Restriction In this section we explain how we calculate the energy distributions on the mesh from a given set of particles. In other words we will discuss the different possibilities for the interpolation operator W(r—ri) appearing in equation (2.58). At the end of the section we explain how to use the same operators to restrict the coefficients of the geodesic equations calculated on the mesh (i.e. the geometric quantities), to the particles' positions. We can define a hierarchy of interpolation operators which can be classified by how many derivatives of the resulting density (interpolant) are continuous. The lowest-order interpolation (the resulting density is discontinuous, piecewise constant), called nearest grid point interpolation, assigns the density of each particle to its nearest grid point (NGP). Thus the interpolation kernel can be written as: 1 : \x-xp\< h/2 1 P l (3.20) 0 : otherwise where xp is the position of the particle and h is the mesh spacing of the grid on which W(x — xp) W(x-p) h / 2 i-X i i+1 Figure 3.1: Kernel for the first order interpolation (NGP). we solve the field equations. Figure 3.1 displays a typical NGP kernel operator for a particle at position p, and a grid with grid points at r = n-i, n, r m , etc. This kernel will produce a density Chapter 3. The Code 28 N distribution given by: p(x) = YJPPW(x-xp) (3.21) i=l A commonly used higher order interpolation is called "cloud in cell" (CIC) interpo-lation which distributes the density of each particle into two grid cells. In this case the resulting density is piecewise linear and the kernel takes the form: 1 — \x — xv\lh : \x — xv\ < h i 1 p l - (3.22) 0 : otherwise Figure 3.2 is a graph of the kernel of the CIC interpolation. In our calculations we have W{x -xp) = { Figure 3.2: Kernel for the second order interpolation (CIC). used the CIC interpolation operator. To restrict the Christoffel symbols (and other geometric quantities) calculated on the mesh to the particles' positions we use the same kernel in the following way: F{xp) = £ F{xi)W(xi - xp) (3.23) i=l where CCp IS the position of the particle, the grid points, and F is any of the coefficients that appear in the evolution equations. The coefficients F are products and quotients of the metric coefficients and their derivatives. In order to calculate derivatives we use the following finite difference approximation on the mesh: [q\ = (qi+l - /h + 0(h2) (3.24) Chapter 3. The Code 29 and then use equation (3.23) to find an approximate value of q' at the particle's position. 3.5 Initial Set of particles. To initialize the set of particles which we evolve, we specify the particle distribution (number of particles per unit of areal coordinate) and the velocity distribution (specifi-cally the number of particles per pr and I). This correspond to a distribution function in the sense of (2.1) which is separable in the r, pr and / variables. In other words: Moreover, instead of specifying P as a function of p r we specify P = P{pr) where Pr = Pr/a. This allow us to calculate the value of pl = yjm2 + p* + l2/r2, and therefore p, jr and S, for each particle without knowing the geometry in advance. This allows us to decouple the tasks of specify initial data for the particles, and ensuring that the initial data satisfies the constraint equations. We use 1-dimensional Monte Carlo techniques applied to each of the functions R(r), P(pr), L(l) to get an specific ensemble of the particles. In theory the statistical error that we produce goes like 0 ( i V _ 1 / 2 ) , where N is the number of particles that we use to sample the distribution function (see section 4.3). f(r,pr,l)=R(r)P(pr)L(l) (3.25) Chapter 4 Checks In this chapter we explain different checks we have performed to test that we have properly implemented the algorithm. The first of these test was checking if the quantity given by equation (2.88) was conserved in the flat spacetime case. In our code this value was conserved up to the tolerance of the LSODA integrator for initial data which was gaussian in both the radial and velocity distributions. We have also checked that a test particle moving in a fixed background follows the geodesic of the background, and that its energy is conserved. The rest of the checks that we describe here are for the fully coupled problem. In section 4.1 we discuss the conservation of the mass aspect function at the outer boundary (r = r m a x ) . We also have simulated certain clusters of particles in circular motion known as Einstein clusters; results from these simulations are discussed in 4.2 in section 4.2. At the end of the chapter we study how the numerical errors scale with the number of particles we use to sample the distribution function / , and as well as with the mesh spacing h. 4.1 Mass Conservation The value of the mass function given by equation (2.89) at infinity should be a conserved quantity if the spacetime is asymptotically flat (see for instance [19]). As an approxi-mation to the value at infinity we can monitor the value at the maximum value of the radial coordinate, M ( t , r m a x ) . This check will make sense as long as no matter leaves the 30 Chapter 4. Checks 31 range of integration (note that energy cannot escape in the form of radiation since that is forbidden in spherical symmetry). In Figure 4.1 we plot M(t,rmax) for three different mesh scales (NT = 64, Nr = 128 and 7Yr = 256) and with a fixed number of particles, N = 10000. 0.8 0.6 x (0 E t-, I  I-- 0.4 <o a 0.2 0 10 20 30 40 50 t Figure 4.1: Value of 2M(t, r)/r at r = r m a x as a function of time. We observe that the value of M is conserved until t « 40 which s the time when the particles begin to leave the range of integration. In the inset we see a detail of the function between t = 0 and t = 40 showing that the conservation is better than 1% for Nr = 64 and improves somewhat as we increase the number of grid points. Theoretically (and in the limit N —> oo) we should observe that the variation of the value of mass should decrease linearly with the number of points, i.e. first order. Although the convergence test here is not as definitive as it tends to be for pure finite-difference methods, the inset of Fig. 4.1 provides fairly good evidence that our code's mass conservation would converge linearly with h in the continuum limit. Chapter 4. Checks 32 4.2 Einstein Clusters Probably the strongest check of the reliability of our code is its ability to simulate clusters of particles that only have rotational motion (zero radial component of the 4-momentum). Let us calculate what initial conditions we must set in order to obtain such a cluster. Assuming we know what the energy density of the cluster of particles is, let us cal-culate the velocity profiles required in order to produce a cluster in equilibrium. In the process of obtaining the velocity profiles, we will also find the geometry of the space-time created by the cluster. The metric coefficient a(r) can be calculated using equation (2.89). If the spacetime is static then /3(r) = 0 and the resulting equation is: a t r ) = ( 1 _ 2 ^ ) ) - , i (4.1) where we calculate M(r) by integrating the density p(r) over the radial coordinate, M(r) = / A-xr2 p{r)dr. To calculate the lapse function we can make use of the evolution equation for the extrinsic curvature (2.23), which in this case reduces to: ^—^ = aa — Anaar2p. (4.2) After a little bit of manipulation this equation can be written as: ol I - a 2 M(r) 1 (4.3) a 2r r 2 ( l - 2 M ( r ) / r ) ' We want to find the angular momentum l(r) needed to get the particles rotating with zero radial component of the 4 momentum. Imposing that condition (pr = 0 and dpT/dt = 0) on the geodesic equation for pr, equation (2.77), we obtain: da _t t l2a Taking into account that or r^p1 f = . L a + i l , (4.5) Chapter 4. Checks 33 and using equation (4.3), we find the equation for l(r): j2/ \ 2 2 (M(r)y l*(r) = mlrl 1 1 (4.6) r J ( l - 3 M ( r ) / r ) This is the same expression as the angular momentum for a particle in circular orbit around a Schwarzschild black hole with mass M(r) (see [26]). In Figure 4.2 we have plotted the areal coordinate for some particles in a cluster with gaussian energy density profile and total mass M « 0.57 as a function of time. Further we can see that the particles retain roughly the same radial coordinate for a time longer than 500M (2 com-plete orbits). Also, the motion appears to be stable since the particles oscillate radially with roughly constant amplitude (see inset). We found that for clusters with gaussian <- 4 5.8 5.6 E=-"5.4 5.2 5 50 100 150 200 100 200 300 Figure 4.2: Radial coordinates of 7 particles as a function of time for a static cluster with gaussian energy density profile. energy densities and initial maximum 6 M / r = 1.01, the cluster becomes unstable and formed a black hole, whereas for clusters with initial maximum 6 M / r = 0.97 the cluster Chapter 4. Checks 34 was stable1. This result agrees with the result that Einstein clusters become unstable for 6M/r = 1 (see [9] and references therein). We have also simulated distorted static clusters for which we specify (1 — e) times the angular momentum fixed by (4.6); this allows us to control the deviation from the ideal static case. We have observed results similar to the case just described where the perturbation is effectively introduced by the numerical approximation itself; i.e. our distorted clusters tend to be stable when max (6M(i, r)/r) < 1. (4.7) t,r 4.3 Numerical Errors. There are two principal sources of numerical errors in our code: 1) sampling of the initial distribution function with a finite number of particles (statistical error), and 2) approxi-mation of the derivatives in the field equations by finite-difference operators (truncation error). In the continuum limit — where the number of particles, iV, approaches infinity, and the mesh spacing, h, goes to 0 — we need to check that these numerical errors go to zero to establish that, in this limit, we are solving the continuum system of equations (Vlasov equation coupled with Einstein's equations). First let us study how the statistical error varies with the number of particles. We can estimate the statistical error in one function by subtracting the functions computed using two different initial sets of particles (two different sets, each sampled using a Monte Carlo method and the same distribution function). Here, the specific function which we monitor is 2M(t,r)/r, and we take the L2 norm of the difference to be an estimate of the statistical error at each time. In other words, if ( 1 ) (2M(t,r)/r) and ( 2 ) (2M(t,r)/r) ^ h i s cluster actually achieves gets 6 M / r > 1 during certain parts of the evolution. Chapter 4. Checks 35 are calculated from two different sets of particles (generated from the same distribution function), then the statistical error estimate at a given time tQ is computed by: &s(t0,r) = < x > (2M(t0,r)/r) - ™ (2M(t0,r)/r) (2) (4.8) Furthermore, we take the L2 norm of this function: As{t0) = Nr (4.9) In Figure 4.3 we have plotted these estimates as a function of the number of particles 2.5 3 Log(N) Figure 4.3: Statistical Error. at the initial time tQ = 0 and for fixed mesh spacing, h. We have displayed 9 different estimates for the statistical error, and the shaded area of the graph shows bounds for the ensemble. The average slope is m = —0.44 ± 0.08 where the uncertainty is the standard deviation calculated from the individual slopes. We see that these results roughly agree with the expectation that the statistical error should scale as N'1/2. Let us turn our attention now to the truncation error. We assume a Richardson Chapter 4. Checks 36 _ 4 5 I i i I i i i I i i i I i i i I i i i! - 2 -1 .8 -1 .6 -1 .4 Log(h) Figure 4.4: Truncation (dotted lines) and statistical error estimates (dashed lines) at t=0. expansion for a mesh function qh(t,r), calculated with a given mesh spacing h: qh(t,r) = q(t,r) + Aqs(t,r) + hpep(t,r) + 0(hP+1) (4.10) where Aqs is the statistical error. We estimate the truncation error in qh(t,r) via qh(t,r)-qh'(t,r) (4.11) where h! < h, specifically we have used log(h') « —2.4). In this way we use the function calculated with h' as a better approximation (a reference) to the analytic solution than the one with larger h. It is important that the statistical error A o s <C hpep otherwise this recipe to calculate the estimate for the truncation error is meaningless since qh{t,r) — qh' (t, r) would be dominated by the statistical error. In Figure 4.4 we show the base ten logarithm of the L2 norm of the estimates of both the truncation and the statistical errors at t = 0. We have plotted the truncation error estimates with dotted lines, and the statistical ones with dashed lines. Chapter 4. Checks 37 In these calculations we have not fixed the total number of particles, but rather the number of particles per grid cell (the value we use in this particular case is 500 particles per grid cell on average). The reason for fixing the number of particles per cell, rather than the total number of particles, is that the statistical error scales with h (for the total number of particles fixed) roughly as h~l. Increasing the number of particles as we decrease h, we compensate for this bias by maintaining roughly constant and small statistical error. In the part of the Figure where the statistical error is much smaller than the truncation error we have calculated the average slope for the truncation error (the solid line in the Figure). The resulting slope is m = 2.02 ± 0.08 where again the uncertainty is the standard deviation. At t = 0 we thus observe second order convergence which is in agrement with the scheme that we have used to compute the geometric quantities. Log(h) Figure 4.5: Truncation (dotted lines) and statistical error estimates (dashed lines) at t=20. Figure 4.5 shows the results of our error analysis at t = 20. We see that in this case the slope is not 2 but smaller, m = 1.36 ± 0.08. This is consistent with the observation Chapter 4. Checks 38 made in section 3.3 that the numerical scheme that we have used for the evolution is only first order accurate in h. Chapter 5 Results In this chapter we summarize the main results we have obtained. In order to find the critical solutions (solutions sitting at the threshold of black hole formation) we have performed bisection searches using the total rest mass (sum of the rest masses of the individual particles) of the cluster, Ma, as the search parameter p. In section 5.1 we discuss some characteristics of the critical solutions. In 5.2 we give some indications for the existence of a scaling laws, r ~ — oTn \p — p*\, and we discuss the possibility of universality in the model in section 5.3. The chapter ends with some conclusions and some ideas for future work. 5.1 Is the critical solution Static ? In this section we present evidence that, for initial data with non-zero angular momentum, as we approach criticality the solution for the space time approaches a static solution. We will first focus on one family of initial conditions and afterwards we will try to explain what we observed for different families. The initial distributions on equation (3.25) are: R(r) oc r 2 e ( - ( r - r o ) 2 / A ' ) 0(r) (5.1) P(pr) oc e ( - ( f t - * » ) a / A * a ) (5.2) L{1) oc e H - ' o ) W ) 0 ( / ) (5.3) where pr = pT/a, and 0 is the step or Heaviside function. For r0 — 5, A r = 1, pTo = 0, 39 Chapter 5. Results 40 Ap> = 2, l0 = 12, A; = 2 (we will call this family "almost" time symmetric because the gaussian for the radial momentum is centered around pr = 0) we found the critical solution for a cluster rest mass about M0 « 1.3. In Figure 5.1 we show a few snapshots from the evolution of da/dt resulting from initial data which is close to criticality but which eventually disperses. At early times, da/dt oscillates, but between t — 117 to t = 195 it seems to get close to zero. In the last snapshot (t = 234) we observe how da/dt -A,—sfi-— t=9.5 t=29 t=78 t=117 • v J t=195 t=234 0 5 10 15 20 r Figure 5.1: Snapshots of da/dt. has become negative, corresponding to dispersal of the particles. We observe similar behavior for the time derivatives of all the metric coefficients. In order to see how close to zero da/dt becomes, we show in Figure 5.2, a detail of a(t, r) at t — 156 for three different sets of particles. The statistical error is roughly of the same order as the function itself (this is not true from r « 6 to r « 8 where the three ensembles agree on a non-zero value for the derivative, however we suspect that this feature will decrease if we tune closer to the critical solution and if we use greater <0 0.02 0 -0.02 -0.04 _t=156 . . . I . . . I . . . I Chapter 5. Results 41 resolution). This shows that the derivative of the function could be zero as we anticipate but it doesn't give us a definite answer. The value of da/dt could be smaller that the statistical error. In order to definitively answer this question we will need to have better resolution and decreased statistical errors. As we have commented in section 2.4 an 0.002 F~T 10 r 156 » Ensemble 1 - Ensemble 2 • Ensemble 3 Figure 5.2: Three ensembles of da/dt at t = 156. interesting function we monitor is 2M(t, r)/r. This function is a scale invariant quantity which tells us a how close we are to trapped surface formation. We plot 2M{t,r) max r (5.4) in Figure 5.1. In the inset we show a detail for the period of time when the maximum seems to have a roughly constant value. Again, between t = 140 and t = 190 the statistical errors are of the same order as the fluctuations in 2M(t,r)/r If we accept that the metric coefficients are independent of our time coordinate then we have a stationary space-time. If, in addition, the vector Na = d/dt introduced in Chapter (2) is orthogonal to the spatial hypersurfaces, then the space-time is static. Chapter 5. Results 42 i—'— 1—•—•—=f I* • of — I • I 0 100 200 300 t Figure 5.3: maxr (2M(t,r)/r) as a function of time. This then implies that Qa = Na, or that the shift vector, and in particular the shift function /3(t,r), is zero. In Figure 5.4 we show the evolution of the shift function. During the period when the time derivatives of the metric coefficients are close to zero, the shift function /3(£, r) also is close to zero, or at least of the same order as the statistical fluctuations as before. Thus the solution for the space-time sitting at the threshold of black hole formation could be a static solution (if it is not static, the amplitude of the derivatives and the shift function have small amplitudes relative to the statistical error). We also note that the total current density jr also tends to be of the order of the statistical fluctuations, but that the S£ component of the stress energy tensor is not zero in the critical regime, as we can see in Figure 5.5. Therefore, on average there are the same number of particles with positive (outwards) and negative (inwards) r components of the 4-momentum but the absolute value of pr is not zero on average. u.o A Ensemble 1 * Ensemble 2 * Ensemble 3 0.2 0.76 0.758 0.756 0.754 I 1 1 1 i 100 120 140 160 180 20' Chapter 5. Results 43 t=29 1=0 t=9.5 t=49 t=78 1=117 1 ' 1 1 _t=156 t=195 t=234 0 5 10 15 20 r Figure 5.4: Shift function as a function of time. This result and the fact that the maximum value of 2M(r)/r for this critical solution is ~ 0.76 indicate that this cannot be one of the Einstein clusters we have considered pre-viously since there are no equilibrium clusters (either stable or unstable) with maximum 2M(r)/r larger than 2/3. Chapter 5. Results 44 0.0027 0.0018 0.0009 0 -0.0009 t=0 t=9.5 k Ii t=29 t=49 ft JI t=78 1 *» * i Ii t=117 k a ii l i • I • 1 I i • 1 • • • 1 I • I - t=156 • k -- ft L / t t=195 1 A t=234 4 A , . 1 , . . 1 , , . 1 , , r-0 5 10 15 20 Figure 5.5: Evolution of Sr(t,r). 5.2 Time Scaling A Type I critical solution usually approaches a static or a periodic solution. As was argued in the previous section, we believe the critical solution in the current case is static. As we approach p* the dynamical solution spends more and more time close to this putative static solution. We have calculated the time t that each of the dynamical solutions spends inside some radius r = r0 (in particular r0 = 6). Figure 5.6 shows the plot of T = t — tc versus ln(p) (natural logarithm) where t is the time that the solution with search parameter p spends inside r = r0, and tc is the same quantity for the solution closest to criticallity. We show three different ensembles of particles for the gaussian family given by equa-tions (5.1-5.3) and again using "almost time symmetric" initial data. We also have nor-Chapter 5. Results 45 T i i i [ i r l n | p - p ' | Figure 5.6: Scaling law for the time. The error for each value of o is the standard deviation of the slope computed using a least squares fit. malized the critical solution in such a way that it has A D M mass equal to 1 at infinity1, i.e.: r -> r / M c ( f , r m a x ) (5.5) t t / M c ( f , r m a x ) , (5.6) where M°(t*,rmax) is the value of the mass aspect function at r = r m a x for the solution closest to criticality during the critical regime. We can see that there is a rough linear relation between the time r and ln|p—p*\ with slope — ( 5 . 2 ± 0 . 2 ) (the error is an estimate of the statistical error between different sets of particles): T ~ - ( 5 . 2 ± 0 . 2 ) l n | p - p * | (5.7) Qualitatively this agrees with other type I critical solutions. 1This will be useful when we compare with other families of initial data. Chapter 5. Results 46 5.3 Different Families of Initial Data The results we showed in the previous section are for Gaussian initial data given by equations [5.1- 5.3] with pro = 0 ("almost time symmetric" initial data) and l0 = 7. In this section we show that similar results are found for other families of initial data. We have chosen initial data with pT0 = —4 and l0 = 3, l0 = 5, l0 — 7, l0 = 12, (with same r0 = 5, A r = 1, Ap r = 2, A ( = 2 as before) and find that the critical solution for each of these cases appears to approach a static solution. For smaller values of the angular momentum the mass in the critical solution gets concentrated closer to the origin, which makes it more difficult to resolve the solution with a uniform finite-difference grid. We also have evolved an initial family with the following one particle distributions: R(r) oc ( l - tanh ((r - r 0 ) / A r ) 2 ) 0(r) (5.8) P(pr) ex (1 - tanh((p r -pro)/&Pr)) (5.9) L{1) a (l - tanh ((/ - / 0 ) /A,) 2 ) 9(Z) (5.10) with r0 = 5,A r = 1, pr0 = —4, Apr = 2, l0 = 7, A; = 2. For this distribution we also find that the solution with small p — p* (p being the total rest mass of the cluster as before) spends some time close to an apparently static solution. In Figure 5.7 we show profiles for 2M(r)/r for the different families, each separate profile being selected from the corresponding period of near-critical evolution. Since different initial conditions set different scales for the problem we have normalized the results such that the total A D M mass of the solutions to which they approach is 1 at infinity (rescaling given by equations (5.5-5.6)). We can see that after normalization, the peak of 2M(r)/r is roughly at the same r* = 2.3. We can also see that the better a solution is resolved, the closer it conforms to the best resolved solution (the Gaussian with lQ — 7 and "almost time symmetric" initial data). This is an indication that the critical solution might be universal although, Chapter 5. Results 47 0 2 4 6 8 10 r Figure 5.7: Critical solution for different families of initial data. Chapter 5. Results 48 we again need better resolution and more particles to be able to give a definitive answer. We are also interested in seeing the importance of the distribution of angular momen-tum on the critical solutions. In Figure 5.8 we show r2 Saa(r) where Saa(r) is the function defined by equations (2.58) and (2.56) for the different families of initial data during the critical regime. This function, r2Saa(r), is a dimensionless quantity which measures the square of the angular momentum of the distribution of particles. We have rescaled the radial coordinate (and time) in the same way as in Figure 5.7. We can see that there is i 1 1 1 1 1 1 1 1 r Figure 5.8: r2Saa(r) for the different families during the critical regime. no apparent agreement between the functions calculated from different families of initial data. More work needs to be done in order to see what is the dependence of the critical solution with respect the angular momentum distribution. Moreover, estimates for the truncation and statistical errors for each calculation of the critical solutions must be Chapter 5. Results 49 calculated to properly compare the solutions. (7 ,=5 .0±0.1 a„=5.0±0.1 -20 h h -40 •60 h a=5.7±0.1 Gauss., 1 =3 CT,=5.1±0.1 <72=5.3±0.1 <7,=5.2±0.1 Gauss., 1„=13 (t.s 15 10 <7=5.5±0.1 Gauss., 1 =5 CT=5.9±0.2 Tanh. 1 =7 ln|p-p*| Figure 5.9: Scaling behaviour for different families of initial data . We have also estimated o defined by r = — crln|p — p*\ ( 5 . H ) for the different families. In Figure 5.9 we show the scaling behaviour for the different families of initial data that we have studied in this thesis. The error for each value of a in the Figure is, as beforre, the standard deviation of the slope computed using a least squares fit. We have collected the values that we have obtained for a in Table 5.1. 5.4 Conclusions and Future Work We have studied critical behavior at the threshold of black hole formation for collisionless matter. Our results indicate that, for families with non-zero angular momentua, the critical solution could be static with non-zero radial particle momenta. We have found Chapter 5. Results 50 Family a Gaussian, l0 = 7 (set 1) 5.0 ± 0 . 2 Gaussian, l0 = 7 (set 2) 5.0 ± 0 . 2 Gaussian, l0 = 3 5.7 ± 0 . 2 Gaussian, l0 = 5 5.5 ± 0 . 2 Gaussian, Z0 = 12 4.9 ± 0 . 2 Gaussian, l0 — 12 t.s. (set 1) 5.1 ± 0 . 2 Gaussian, l0 = 12 t.s. (set 2) 5.3 ± 0 . 2 Gaussian, l0 = 12 t.s. (set 3) 5.2 ± 0 . 2 Tanh, l0 = 7 5.9 ± 0 . 2 Table 5.1: Values of time scaling exponent for the different families (t.s. stands for " almost time symmetric"). The errors above are assumed to be the same as the statistical error for the t.s. family (see equation (5.7)). evidence for a life-time scaling law which is to be expected for Type I critical solutions. In order to answer all these questions more rigorously we would need some way to increase our resolution where it is needed. Some kind of adaptive code (such as the ones used in finite difference studies [1]) would be useful. It may be that the development of a finite difference code to solve the Vlasov equation directly in phase space would be the best way to attack this problem. This would allow us to incorporate the adaptivity that we need, as well as providing us with better-understood convergence properties. Bibliography [1] M.W. Choptuik, "Universality and scaling in the gravitational collapse of a scalar field". Phys. Rev. Lett, 70, 9-12, 1993. [2] M.W. Choptuik, P.Bizon, and T. Chmaj, "Critical Behaviour in Gravitational Col-lapse of a Yang-Mills Field". Phys. Rev. Lett, 77, 424-427, 1996 [3] C. Gundlach, "Critical phenomena in gravitational collapse". Adv. Theor. Math. Phys. 2 (1998) 1 [gr-qc/9712084]. [4] A. Einstein, "On Stationary System with Spherical Symmetry Consisting of Many Gravitating Masses". Ann. Math., 40, 922-936, 1939 [5] J.R. Oppemheimer and H. Snyder, "On Continued Gravitational Contraction". Phys. Rev., 56, 455-459, 1939 [6] R.C. Tolman, "Effect of Inhomogeneity on Cosmological Models". Proc. Nat. Acad. Sci, 20, 169-176, 1934 [7] H. Bondi, "Spherically Symmetrical Models in General Relativity". Mon. Not. Roy. Astr. Soc, 107, 410-425, 1947 [8] D .M. Eardley, "Time functions in numerical relativity: Marginally bound dust collapse". Phys. Rev. D, 19, 2239-2259, 1979 [9] S.L. Shapiro and S. Teukolsky, "Relativistic Stellar Dynamics on the Computer. I. Motivation and Numerical Method". Ap. J., 298, 34-57, 1985 [10] S.L. Shapiro and S. Teukolsky, "Relativistic Stellar Dynamics on the Computer. II. Applications". Ap. J., 298, 58-79, 1985 [11] S.L. Shapiro and S. Teukolsky, "The Collapse of Dense Star Clusters to Supermas-sive Black Holes: The Origin of Quasars and AGNs". Ap. J., 292L, 41-44, 1985 [12] S.L. Shapiro and S. Teukolsky, "Relativistic Stellar Dynamics on the Computer. IV. Collapse of a Star Cluster to a Black Hole". Ap. J., 307, 575-592, 1986 [13] F .A. Rasio, S.L. Shapiro, S.A. Teukolsky, "Solving The Vlasov Equation in General Relativity". Ap. J., 344, 146-157, 1989 51 Bibliography 52 [14] S.L. Shapiro, S.A. Teukolsky, "Formation of Naked singularities: The violation of Cosmic Censorship". Phys. Rev. Lett, 66, 994-997, 1991 [15] G. Rein, "Static solutions of the spherically symmetric Vlasov-Einstein system", (unpublished) [gr-qc/9304028]. [16] G. Rein, A.D. Rendall, "Global existence of classical solutions to the Vlasov-Poisson system in a three-dimensional, cosmological setting", (unpublished) [gr-qc/9306021]. [17] A.D. Rendall, "An introduction to the Einstein-Vlasov system", (unpublished) [gr-qc/9604001]. [18] G. Rein, A.D. Rendall and J. Schaeffer, "Critical collapse of collisionless matter: A numerical investigation". Phys. Rev. D58 (1998) 044007 [gr-qc/9804040]. [19] C.W. Misner, K.S. Thorne, J.A. Wheeler, "Gravitation". W.H. Freeman, San Francisco, 1973 [20] M.W. Choptuik, "The 3+1 Einstein Equations", (unpublished notes), 1998 http:// Phy387N/Notes.html [21] R.W. Ffockney, J.W. Eastwood, "Computer Simulation Using Particles". McGraw-Hill Inc.,1981 [22] L.R. Petzold and A . C . Hindmarsh, " L S O D A " . Computing and Mathematics Re-search Division, 1-316 Lawrence Livermore National Laboratory, Livermore, C A 94550. [23] E. Anderson et al, "LAPACK Users' Guide". Society for Industrial and Applied Mathematics, Philadelphia, 1992 [] [24] J. W. York, "Kinematics and Dynamics of General Relativity", in "Sources of Gravitational Radiation", L. Smarr, ed. Cambridge University Press, 1979 [25] M.W. Choptuik, "Numerical Techniques for Radiative Problems in General Rela-tivity". Univ. of British Columbia Ph.D. thesis (unpublished), 1986 [26] R.M. Wald, "General Relativity". The University of Chicago Press, 1984 Appendix A Summary of Equations in Maximal-Areal Coordinates A . l Field equations a a Ki' = I - a 2 3 + ^-ra2Kg2 + Anarp 2r r a a P = 2, a r' arKl = a'{- - -) + ^ ( a 2 - 1 + 2r-) + 4naa& + 3f - 3/5) a r r2 a a1 or 2 (A.12) (A.13) (A.14) (A.15) where and i APS 47rAr r, 51 = 1 A P I 2 1 * I2 4?rAr fr[ r\p\ 47rAr fr{ r\ v W ( x - Xi) = { 1 — | a; — Xi\/h 0 |a; — Xi\ < h otherwise (A.16) (A.17) (A.18) (A.19) (A.20) 53 Appendix A. Summary of Equations in Maximal-Areal Coordinates 54 A.2 Evolution Equations dpr da_t 88 adapr2 l2a . . = --s-p + + ~i^-^r + -=7-3 ( A - 2 1 ^ dt or or a6 or j r prr4 ± = B%L-I, (A.22) dt alpr $ = , / m 2 + ^ + i ! (A.23) V az rz % - -'te-'Wab^-'") (A'24) d0 _ ap* dt ~ d(j> _ ap* ~dt ~ ~]f (A.25) P*2 = -^T-A-i-P9) (A.26) surd \ r 4 / (A.27) Appendix B Direct Approach To Solutions of Einstein-Vlasov System B . l Fie ld Equations The distribution function is defined by: dN f(t,r,pr,l) = —, (B.28) where N is the number of particles and V is the volume in phase space. Then the stress energy tensor can be calculated via: T^ = J p»p»f(t,r,pr,l)±dVp (B.29) where the volume element in momentum space, Vp, is restricted by - m 2 = &*pV (B.30) to the following expression (see [13]): = Jpdpedp, We can introduce momentum coordinates adapted to the symmetry: I2 = Pl + A - n (B-32) sin 9 ib = tan P4> i (B.33) and pl given, as before, by: p ' = P + f , 4 ( B - 3 4 ) 55 Appendix B. Direct Approach To Solutions of Einstein-Vlasov System 56 The volume element dVp in these variables is: dprdl2dip . . The stress-energy tensor components are: Tt = -TilP'fdPrdl2 + [prfdprdl2 (B.36) alrl J aazrz J Tr = -TT- I PrfdPrdl2 (B.37) 2 /prfdprdl2 + / \fdPrdl2 (B.38) T ao?r2 J aHr* J p Te_T<i> _ 1 n rl2 1 0 ~ 2 - f i I %dpTdl2 (B.39) e r r 4 J p-We are able then to construct the source terms: / 2 72 9 = a^If{t>r'PrJ)h2 + ^ 2+72dPrdl2 (R40) jr = —7 / f(t,r,pr,l)prdprdl2 (B.41) a r z 7 5; = J L / P 2 / ( t , r , P r , Q 2 a 3 r 2 7 ^ / m 2 + p 2 / a 2 + j 2 / r 2 Se = JL[ l2H^,pr,l) M 2 ar*J Jm2 + p2/a2 + / 2 / r 2 B.2 Vlasov Equation If the system is collisionless, we have Liouville's theorem concerning the conservation of the volume in phase space, V, and we also know that the particle number N is a conserved quantity. Therefore we can write the conservation of the distribution function, / = N/V, as: ^ = 0 (B.44) dr Appendix B. Direct Approach To Solutions of Einstein-Vlasov System 57 This equation is known as the "collisionless Boltzmann's equation" or the "kinetic equa-tion" . Expanding the total derivative we have: +____v=0 (B45) dr dxa dr dpa where xa are the space-time variables, and pa are the momentum variables. In spherical symmetry this equation reduces to: dl + drdl + dpLdl = Q dt dt dr dt dpr where dr/dt and dpr/dt are given by equations (2.77) and (2.78). This gives: df apr df da_t dp a dap2 al2 df , , dt a2p* dr dr dr a? dr p* pr* dpr Appendix C Summary of Equations in Polar-Areal Gauge Coordinates In this section we present the equations for a second choice of coordinates. Specifically we again adopt areal radial coordinates (implying that the proper surface area of the spheres centered at r — 0 and with radius r will be Airr2), but fix the time coordinate using polar slicing. This slicing condition is implemented by demanding: TiK = K\ = K \ (C.48) which implies K6e — K^^ = Keg = K^^ = 0. Using expression (2.39) we can express K% as: K99=±(0(rb)'-rb) = O, (C.49) we have: P(rb)' = rb (C.50) Choosing areal coordinates, 6=1 ,6 = 0, the above equation implies: p = 0. (C.51) The resulting metric is: ds2 = -a2(t,r)dt2 + a2(t,r)dr2 + r2d£l . (C.52) 58 Appendix C. Summary of Equations in Polar-Areal Gauge Coordinates 59 we then have the following geometrical quantities: R = 2 [-2 (l/o)' + o/r (1 - 1/a2)] /(ar) R\ = 2a'/ (ra3) R% = R^^ = 1/ (ar2) (a - (r/a)') KrT = K = -a/ (aa) K% = K++ = 0 (C.53) = — (<V ( a Q0)' — 2d/ (raa) A 'lr = i ^ r r ( r = — (d/ (aa)) a| r l r = 1/a (a1/a)' a , / = a , ^ = a'/ (a2r) Briefly we summarize the resulting equations of motion: Momentum Constraint: a = —4irraajr (C.54) Hamiltonian Constraint: a' = 4irra3p-^{a2-l) (C.55) Evolution of the 3-metric 9 i ] = -2oc9ikK) (C.56) Evolution of the extrinsic curvature k\ = - - ( - ) + ' ^ - + 4TTa(S-p)-SiraSrr (C.57) a \ a J ra4 K% = 0 = — < ^ + —Aa--+r^\+ATia{S-P)-^aS66 (C.58) alr arz \ a a1) k% = 0 = - ^ - + ^ r f a - - + r 4 l+47ra(5 -p ) -87ra5^ (C.59) alr arz \ a az J where the two last equations are zero by the choice of slicing. Appendix D Pseudo-code routine vlasov # n := position of the i-th particle # \Pn\i : = rnomentum of i-th the particle # '•— T£ component on the j grid point # a,j, aj and Bj are the metric coefficients # Mj := mass aspect function # Generate the initial positions r; and velocities [pM]j for the particles read parameters of the distribution function for i = 1 . . . N p a r i i c i e s Generate rj and using Monte Carlo for a given distribution function end for t = 0 # Time loop while t < tmax for j = 1... Ngridpoints Calculate [T£]j from rj and [pM]i end for Solve for aj, Bj and aj on the mesh # Check apparent horizon formation MIMjJTj = 1 stop Ouput aj, Bj and aj for i = l.-Nparticl es Calculate [T^]j at t = t + 1/2 di from a,-, /3j and aj end for t = t + dt for i — ^••••Npariicies Integrate geodesic equations to find new r-j and [p^ ]* if n > r m a x Kill i-th particle, Npartides = Npartides — 1 if r-j < 0 Reflect the particle if Npartides = 0 Stop end for end while end routine 60 


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



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


Related Items