{'fraction_non_alphanumeric': 0.03936831315266833, 'fraction_numerical': 0.01652892561983471, 'mean_word_length': 4.263572753329961, 'pattern_counts': {'":': 0, '<': 38, '<?xml version=': 0, '>': 20, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Numerical Simulation of General Relativistic Shock Waves by a Locally Inertial Godunov Method Featuring Dynamical Time Dilation

We introduce what we call a locally inertial Godunov method with dynamical time dilation, and use it to simulate a new one parameter family of general relativistic shock wave solutions of the Einstein equations for a perfect fluid. The forward time solutions resolve the secondary reflected wave (an incoming shock wave) in the Smoller-Temple shock wave model for an explosion into a static singular isothermal sphere. The backward time solutions indicate black hole formation from a smooth underlying solution via collapse associated with an incoming rarefaction wave. As far as we know, this is the first numerical simulation of a fluid dynamical shock wave in general relativity.


Introduction
This thesis presents a numerical study of a canonical family of spherically symmetric solutions to the Einstein equations of General Relativity, such that shock wave formation is simulated in forward time and black hole formation is indicated in backward time. To start the simulations, we introduce a new one parameter family of initial data satisfying the constraints of the Einstein equations, and to perform the simulations, we develop a locally inertial Godunov method incorporating a new dynamical time dilation feature to account for the relativistic effects of curvature.
Our point of departure is the work of Groah and Temple [GT04], who gave the first existence theory for shock wave solutions of the Einstein equations, starting from general initial data of bounded variation in the density and velocity. The setting is restricted to spherically symmetric spacetimes in standard Schwarzschild coordinates. For the analysis, Groah and Temple introduced the locally inertial Glimm scheme, a numerical method based on approximating spacetime by flat Minkowski space in each grid cell, such that the method was amenable to detailed mathematical estimates sufficient to prove convergence to weak (shock wave) solutions of the Einstein equations. The curvature of spacetime in the locally inertial Glimm scheme is accounted for by coordinate transformations between the locally inertial approximations in each grid cell. This technique is pedagogically interesting because it parallels the so called correspondence principle, the physical principle that general relativity should reduce to special relativity in sufficiently small neighborhoods of spacetime. Groah and Temple concludes in [GT04] that the curvature of spacetime necessarily becomes discontinuous at shock waves, resulting in solutions solving the Einstein equations in the weak sense of the theory of distributions. In this thesis, we develop these locally inertial formulations of the Einstein equations into a viable numerical method, what we refer to as a locally inertial Godunov method with dynamical time dilation, and we introduce a new one parameter family of interactive solutions on which we demonstrate the convergence of the method.
Our numerical demonstrations are backed up by a general theorem which we prove.
Namely, we prove if a sequence of approximate solutions, generated by our locally inertial Godunov method converges pointwise almost everywhere to a function along with the total variation of the fluid variables remain uniformly bounded under the limit, (exactly what we can demonstrate numerically), then the limit solution is an exact weak solution to the Einstein equations. By this general theorem, we need only demonstrate numerical convergence to a limit, with bounded oscillations, in order to conclude our simulated solutions accurately represent exact (weak) solutions of the Einstein equations.
The one parameter family of spacetimes on which we test the convergence of the method are interesting in their own right. This one parameter family is constructed by matching initial data from a critical (k = 0) Friedmann-Robertson-Walker (FRW) spacetime to initial data for a Tolmann-Oppengeimer-Volkoff (TOV) spacetime at a point where the spacetime is only Lipschitz continuous, the smoothness required for our convergence theorem. Since our numerical method is based on standard Schwarzschild coordinates, this requires a coordinate transformation from FRW coordinates to standard Schwarzschild coordinates. Such a transformation requires an integrating factor, and we use the integrating factor derived by Smoller and Temple in [ST95], (we call the resulting standard Schwarzschild coordinates version of FRW the FRW-1 spacetime). We also derive a new integrating factor that leads to a different standard Schwarzschild coordinates version of FRW, (which we call FRW-2).
The initial discontinuity that separates FRW from TOV in the standard Schwarzschild coordinates initial data then generates a region of interaction between the two spacetimes, a region for which there is no closed form solution, and can only be constructed by numerical simulation. The ultimate goal of this thesis is then to simulate this region of interaction, and demonstrate convergence of the locally inertial Godunov method.
The numerical simulation of the region of interaction between FRW and TOV is a problem motivated by the original work of Smoller and Temple in [ST95], further developed in [ST03,ST04]. In [ST95], the first exact shock wave solution of the Einstein equations is constructed by matching an FRW spacetime to a TOV spacetime across a shock interface, such that the resulting matched solution, given by exact formulas, is a true weak solution of the Einstein equations. The resulting matched spacetime (we denote as the FRW/TOV metric) with a shock wave was a simple model for an explosion into a static singular isothermal sphere, relevant to models of star formation [ST95]. Based on Smoller and Temple's work, we built a computer visualization [JSV] to demonstrate the qualitative behavior of their model. The Smoller and Temple (Sm/Te) construction requires taking a different sound speed on the FRW and TOV sides of the shock wave finely tuned to reduce the region of interaction between FRW and TOV to a single pure shock wave. As a result, they only required the existence of, not an explicit formula for, the integrating factor. In our simulation, we take the same equation of state on both sides, a more realistic assumption, and as a consequence, a region of interaction, for which there are no exact formulas, between the lightlike curves emanating from the initial point of connection between FRW and TOV metrics is generated. Based on this background, we interpret our numerical FRW/TOV model simulation in forward time as resolving the secondary reflected wave in the Sm/Te model of an explosion into a static singular isothermal sphere. In particular, our simulation confirms for the first time the reflected wave is a shock wave, not a rarefaction wave. We also simulated the time reversal of the FRW/TOV solution. In this case our simulation indicates the formation of a black hole in the resulting solution, which consists entirely of rarefaction waves. Based on our simulation, one can only expect the observed strong rarefaction wave evolving toward the coordinate center will form a black hole. Note that a black hole cannot be computed exactly in standard Schwarzschild coordinates because black holes generate coordinate singularities in this coordinate system. A proposal for future work would be to simulate the black hole all the way by transforming the solution to Eddington-Finkelstien or Kruskal coordinates [Wei72].
We start our discussion by briefly introducing the Einstein equations in Chapter 2. In this introduction, we formulate these equations in the simplest setting for shock waves in General Relativity; this setting is the spherically symmetric Einstein equations with the stress energy tensor for a perfect fluid and an equation of state p = σρ, ρ = const. In this setting, the Einstein equations reduce to a conservation law with source terms coupled with a system of ODEs.
Chapter 3 sets out the background information for the FRW, the TOV, and the FRW/TOV matched models. In particular, we construct a one parameter family of initial data that agrees with the FRW and TOV metrics in standard Schwarzschild coordinates on either side of a single point where they match continuously. The continuity of the metric at the matching point implies the initial data meets the constraints of the Einstein equations introduced by Groah and Temple in [JGB06]. More specifically, the initial data is built by connecting a FRW metric with a TOV metric across an initial discontinuity in the fluid variables, such that the metric components are matched Lipschitz continuously. This one parameter family of initial data enables us to construct the one parameter family of interacting spacetimes that solve the Einstein equations. To be simulated by our locally inertial Godunov scheme, the metric must be in standard Schwarzschild coordinates, and while the TOV is already in this coordinate system, the FRW metric needs to be mapped over. We demonstrate two mappings of the FRW metric over to standard Schwarzschild coordinates, denoted as FRW-1 and FRW-2. With these two coordinate transformations, the one parameter family of shock wave solutions, the FRW/TOV metric, has two different forms based on the coordinate system used. These two forms are used as a further test of the correctness of our numerical simulations, by verifying the simulations produce the same spacetime in two different coordinate systems.
In Chapter 4, we introduce the locally inertial Godunov method. This chapter begins with a discussion of the solution to the Riemann problem [Smo83] for relativistic compressible Euler equations in Minkowski spacetime. The Riemann problem is the fundamental building block of the Godunov method, [LeV92]. We outline the solution to this Riemann problem as demonstrated by Smoller and Temple in [ST93]. In order to incorporate the Riemann problem in the Godunov method for simulating non-flat spacetimes, we augment the solution to the Riemann problem with the time dilation feature. To this end, we determine the effects of time dilation on the Godunov step and incorporate it into our locally inertial Godunov method. The chapter ends with a layout of the locally inertial Godunov method and all the steps necessary to implement it. This method is a fractional Godunov scheme that defines the solution inductively and contains four major steps: a Riemann problem step, a Godunov step (with time dilation), an ODE step, and an update step.
In Chapter 5, we prove the main theorem that establishes the consistency of our locally inertial Godunov method with time dilation. It states, assuming the convergence of our approximate solution constructed using our method with a total variation bound at each time step, the convergent solution is indeed a weak solution to the Einstein equations. The two main difficulties of this proof are handling the jump in the fluid variables across each time step and showing the jump in the metric across each space step is accounted for by the term added to the ODE step. The assumptions to this theorem are obtained numerically by the simulation results. Thus, our theorem is perfectly tailored to our numerical simulation: we numerically establish convergence and a total variation bound, and this main theorem proves that assuming these conditions, we can conclude convergence to a weak solution of the Einstein equations.
In Chapter 6, we begin the numerical study of the locally inertial Godunov method with dynamic time dilation by simulating the pure non-interacting FRW and TOV spacetimes individually. Each of these continuous models provide us with a different scenario to test and fine tune our method. These models are also important for correctly handling the non-interaction regions and the ghost cells for the subsequent shock wave model. Within this chapter are graphs of the solutions and tables of the convergence results, along with the explanations behind them.
In Chapter 7, the first glimpse of shock wave solutions to the Einstein equations is shown by simulating the FRW-1/TOV model. We demonstrate this one parameter family converges and produces quantitatively different solutions as we vary the free parameter.
Regardless of the parameter chosen, the solution always contains two shock waves, an incoming and outgoing wave, enclosing a region of higher density, denoted as the interaction region. Over and above the convergence and graphical results, we show the interaction region converges to the cone of sound (the extent of the sound like information emanating from the initial discontinuity), and we also show outside this interaction region the FRW-1 and TOV metrics are preserved under the simulation. Changing coordinate systems, the FRW-2/TOV model simulation is similar to the FRW-1/TOV model. We discover and numerically confirm the FRW-2/TOV model is the same solution as the FRW-1/TOV model that differs by a non-linear coordinate transformation, and this second model provides a pedagogically interesting numerical confirmation of the covariance of the Einstein equations in standard Schwarzschild coordinates.
In Chapter 8, we examine the time reversed FRW/TOV model. The numerical simulation shows this solution consists of two rarefaction waves, one outgoing and one incoming, encompassing a region of lower density. In particular, the simulation indicates there are no shock waves present. Therefore, this simulation shows the underlying solution is smooth and thus a strong solution to the Einstein equations, unlike the forward model. We record the convergence of this solution and the convergence of the interaction region to the cone of sound along with the preservation of the non-interaction regions. The simulation convinces us the incoming rarefaction wave will form a black hole in finite proper time. That is, we extend the time of the simulation up until extreme relativistic effects of time dilation, characteristic of a black hole, appear. We cannot simulate all the way up to black formation in standard Schwarzschild coordinates because black holes are coordinate singularities in standard Schwarzschild coordinates [ST97b], but the metric component A is monotonically decreasing to a small value before time dilation makes our numerical calculations impractical. Indeed, we numerically demonstrate the mass function near the black hole formation gets within 1.084 ≤ 1.125 = 9/8 of the Schwarzschild radius of that mass; thus, we have numerically confirmed the solution gets within the Buchdal statiblity limit of 9/8ths the Schwarzschild radius. Buchdal's theorem states that whenever the mass of the star gets within 9/8ths of the Schwarzschild radius for that mass, there is no static configuration with a finite pressure at the center capable of holding the mass up [ST97b]. The implication then is a star that gets within 9/8ths the Schwarzschild radius cannot be prevented from collapsing into a black hole; thus, we view our results here as a numerical physical proof that black hole formation must occur in the time reversed FRW/TOV model. Thus, the simulation indicates a black hole forming out of a smooth solution to the Einstein equations. With the success of running the forward and reverse FRW/TOV models, we have simulated general relativistic analogs of both the 1 and 2 family of shock and rarefaction waves, that is, the analogs of all the elementary waves of conservation laws. This fact gives us confidence in the ability and accuracy of our method to simulate any solution to the Einstein equations for a perfect fluid in standard Schwarzschild coordinates.
In Chapter 9, we put dimensions back into the problem, giving our simulations of the earlier chapters a physical context. For convenance and concreteness, we set the Einstein's In Appendix A, we outline the code used to perform the simulations in this paper. A CD is attached that contains all this code, which is approximately 8,000 lines. The appendix is designed to provide the reader an overview of the organization of the code and to help the reader navigate to specific areas of interest.

Preliminaries
This chapter introduces the background to the topic of this thesis, the locally inertial Godunov method with dynamical time dilation. This is a fractional step numerical method that approximates the convective part of the fluid dynamics in each grid cell by the solution of a Riemann problem in a (flat) Minkowski background spacetime, such that each cell is endowed with a time dilation factor. Averages of the fluid variables are taken, and a second fractional step then evolves these averages according to evolution by a dynamical system that accounts for the undifferentiated source terms. The main issue is at shock waves the gravitational metric has a jump discontinuity in the derivative, so the curvature of the gravitational metric becomes discontinuous at shocks, and the purpose of this thesis is to confirm our method is effective at numerically simulating general relativistic shock waves.
The locally inertial Godunov method with dynamical time dilation is based on the locally inertial formulation of the Einstein equations for spherically symmetric spacetime metrics in standard Schwarzschild coordinates and was first derived by Groah and Temple in [GT00]. Groah and Temple provided an exposition of their work along with the work of Smoller and Temple in [JGB06]. All the equations in this chapter are taken from the joint work of Groah, Smoller, and Temple in [JGB06]. In this study, we consider the stress energy tensor a perfect fluid with an equation of state p = σρ 1 , where σ ≡ Const, the simplest possible setting for shock wave propagation in general relativity. In this case, the particle number density and entropy decouple from the energy and momentum equations, and the fluid equations close at the level of the Einstein equations alone. This equation of state models flow at constant temperature, and in the case σ = c 2 /3, it models the important cases of free particles in the extreme relativistic limit, as well as the case of pure radiation, c.f. [Wei72]. In [JGB06], Groah and Temple show that the four standard Schwarzschild coordinates Einstein equations are weakly equivalent to the system of PDE's obtained by 1 C.f. [JGB06] where σ 2 is used in place of σ

Introduction to Einstein Equations
The fundamental equations of general relativity, the Einstein field equations, is given by where G is the Einstein curvature tensor, T is stress energy tensor (the source of gravity), and κ is Einstein's coupling constant where G is Newton's gravitational constant, and c the speed of light. For concreteness and convenience, we take the speed of light as c = 1 and Newton's gravitational constant as G = 1 throughout this paper. The consequences of setting these constants to one is discussed in detail in Chapter 9. These equations describe how the mass/energy of an object curves space and how the curvature, in turn, stretches or squeezes matter in the three spatial directions.
The component form of the Einstein equations is denotes the individual components of the Einstein curvature tensor. We let i, j ∈ {0, 1, 2, 3} refer to the individual components in a given coordinate system. Throughout, the Einstein summation convention is used where repeated up-down indices are summed from 0 to 3, unless otherwise noted.
Modeling the sources of the stress energy tensor by a perfect fluid, T is given by where w denotes the unit 4-velocity vector of the fluid (tangent vector to the world line of the fluid particle), ρ denotes the energy density (measured in the inertial frame moving along with the fluid), and p denotes the fluid pressure.
For simplification, the only metrics under consideration are spacetime metrics g which are spherically symmetric. A general spherically symmetric metric takes the form where the metric components A, B, C, D are functions of the time and space variables (t, r) only, dΩ 2 = dθ 2 + sin 2 θdφ 2 represents the standard line element of the unit 2-sphere, and x ≡ (x 0 , . . . , x 3 ) ≡ (t, r, θ, φ) is the spacetime coordinate system. A spherically symmetric metric is in standard Schwarzschild coordinates 2 if it is written in the more simple form of (2.7) ds 2 = −B(t, r)dt 2 + 1 A(t, r) dr 2 + r 2 dΩ 2 .
A classic argument transforming a spherically symmetric metric (2.6) into standard Schwarzschild coordinates (2.7), which is thoroughly explained in Chapter 3.
In this paper, the simplest setting for shock wave propagation in General Relativity is assumed: a spherically symmetric metric in standard Schwarzschild coordinates (2.7) where the sources are modeled by a perfect fluid (2.5) with an equation of state where σ is assumed to be constant and √ σ represents the sound speed. Using MAPLE to put the standard metric (2.7) into the Einstein equations (2.1) gives us the following system of four coupled partial differential equations, where the quantity Φ is where "prime" represents ∂ ∂r and "dot" represents ∂ ∂t . The mass function is M ≡ M (t, r) is defined through the identity (2.14) which is interpreted as the total mass inside radius r at time t. Using this total mass, Assuming a perfect fluid (2.5), the components of the stress energy tensor T ij satisfy the following relations x (T 00 (2.26) Here, {} ,0 and {} ,1 represent the derivative with respect to the time and space variable, respectively. The other two equations (2.9) and (2.11) are rearranged as The variable x is written instead of r when the equations are expressed as a system of conservation laws consistent with the literature on the subject. This conservation form is an advantageous formulation because equations (2.25) and (2.26) form a system of conservation laws with source terms, where the conserved quantities u = (T 00 M , T 01 M ) are independent of the metric (c.f. (2.21) and (2.22)). Moreover, all the time derivatives of the metric (2.10) and (2.12) cancels out, and the source terms to the conservation laws (2.25) and (2.26) are independent of any derivatives, resulting in a solvable system. The resulting system (2.25)-(2.28) can be written compactly as (2.29) Here, the conserved quantities are the Minkowski energy and momentum densities and the metric components are written as The flux function becomes The source terms of the conservation laws, (2.25) and (2.26), are (2.35) The LHS of the metric ODEs, (2.27) and (2.28), are Thus, for weak solutions to the Einstein equations, the coupled PDEs (2.9)-(2.12) are equivalent to a conservation law with source terms paired with a system of ODEs (2.33). In this thesis the convergence of the locally inertial Godunov method (defined in Chapter 4) is shown to be a consistent and effective first order method for computing solutions of system (2.33), equivalent to the Einstein equations when shock waves are present and the gravitational metric is C 1,1 . Such solutions only solve the Einstein equations in the weak sense of the theory of distributions. Our general system (2.33) applies to the case of spherically symmetric metrics in standard Schwarzschild coordinates, and in the next chapter, a list of candidates meeting this criteria is assembled.

One Parameter Family of Shock Wave Solutions
This chapter presents the theory behind the one parameter family of shock wave solutions of the Einstein equations whose numerical simulation is the subject of this thesis.
Our point of departure of this chapter is the paper [ST95] where the metric components A, B, C, D are functions of the time and space variables (t, r) only. The expression dΩ 2 = dθ 2 + sin 2 θdφ 2 represents the standard line element of the unit 2-sphere. A metric in standard Schwarzschild coordinates takes the form where again the metric components are functions of the time and space alone. Throughout this paper, the barred coordinates (t,r) are reserved for metrics in a standard Schwarzschild coordinates system (3.2) while the unbarred coordinates (t, r) are for metrics that are only spherically symmetric (3.1). This FRW metric transformation in not only necessary for the matching, but it is also a prerequisite for the locally inertial Godunov scheme since it is only capable of simulating metrics in standard Schwarzschild coordinates. When constructing these shock wave solutions, Sm/Te proved the existence of the necessary coordinate transformation, but they did not require any detailed information about it. Here, we do not have this luxury since this information is essential to run the simulation. Our main goal of this chapter is to uncover these details on transforming the spherically symmetric FRW metric into standard Schwarzschild coordinates. We discover there exists at least two of these transformations of the FRW metric, denoted as FRW-1 and FRW-2. In the subsequent chapters, these FRW metrics expressed in different standard Schwarzschild coordinates are used to test our locally inertial Godnuov method, by simulating these continuous metrics.
In Section 3.1, we briefly discuss the FRW metric. The FRW-1 metric is obtained by providing a candidate for the transformation of the FRW metric into standard Schwarzschild coordinates and proving the validity of this transformation. The fluid variables also get mapped over as functions of the new coordinates. In Section 3.2, we review the standard procedure for mapping a spherically symmetric metric to standard Schwarzschild coordinates, and the FRW-2 metric is constructed using this procedure. The other metric in the matching, the TOV metric, is discussed in Section 3.3. This metric is already in standard Schwarzschild coordinates, but the main aspects of it needed in this paper are explained in this section. We conclude with the details of matching the FRW and TOV metrics continuously across a shock wave interface in Section 3.4. It is in this section, we make precise the one parameter family of shock wave solutions by defining the shock wave surface.

The FRW-1 Metric
The conformally flat (k = 0) FRW metric takes the form The fluid density, as a function of time t, associated with this metric takes the form is Einstein's coupling constant (2.2). To determine the fluid velocity, the fluid in this metric is comoving relative to the background. More precisely, its velocity vector, w = (w 0 , w 1 , w 2 , w 3 ), has the form w i = 0, for i = 1, 2, 3, or, in words, the movement of the fluid is limited to the time coordinate and none of the space coordinates. The FRW metric g being diagonal and the velocity vector having unit length implies This velocity vector is expressed in terms of the classical coordinate radial velocity v ≡ dr/dt of the particle paths of the fluid as v = 0.
Recently, Smoller and Temple derived a new set of equations with solutions describing a two parameter family of expanding wave solutions of the Einstein equations which includes the critical Friedmann universe in the standard model of cosmology [ST]. In deriving this set of equations, they discover a coordinate map for the FRW metric (3.3) over to standard Schwarzschild coordinates (3.2) for σ = 1 3 , and where R(t) = √ t is the scale factor. We state their result and derivation in the following theorem goes over to the following metric in standard Schwarzschild coordinates where the fluid velocity v is related tor/t by Proof. To map the FRW metric tensor (3.3) from the unbarred coordinates (t, r) over to the barred coordinates (t,r), we use the tensor transformation law More specifically, at each point, g transforms by the matrix transformation law (3.11)ḡ = J T gJ for a bilinear form since the Jacobian matrix J = ∂x i ∂x α transforms the vector components of thex-basis ∂ ∂x α into their corresponding coordinates of the x-basis ∂ ∂x i . The inverse of the Jacobian matrix is easily computed by taking partial derivatives of the coordinate transformation equations (3.7), resulting in So the Jacobian matrix is where (3.14) The metric in barred coordinatesḡ is computed as (3.15) Using (3.14) to solve for gives us Note, the new metric is a function ofr/t, and it is useful to define a new variable η as As a function of η, the metric becomes or equivalently it is written in standard Schwarzschild coordinates form (3.2), (3.20) The metric is also a function of the fluid velocity v, shown by finding the relationship between v and our variable η. To this end, the old velocity vector w is related to the new velocity vectorw using the tensor transformation law so for our comoving velocity vector w = (1, 0, 0, 0) this law simplifies to Using (3.23) with our coordinate transformation (3.7), the new velocity vector coordinates are For a spherically symmetric metric in standard Schwarzschild coordinates (3.2) the radial velocity v ≡ dr/dt is computed as For our metric (3.20), A = B −1 , implying √ AB = 1. By (3.18), the relationship between v and η is We substitute this relation (3.26) into (3.20) to prove (3.8). By defining the variable ξ =r/t, the coordinate transformation (3.7) implies proving (3.9).
The FRW metric in standard Schwarzschild coordinates (3.2) (3.28) is a function of ξ alone. We consider now the fluid variables (ρ, v). The fluid velocity v is implicitly a function of ξ by equation (3.9), but through manipulation the fluid velocity v is a function of ξ explicitly where the minus sign is taken to ensure v < 1. The density ρ is a function of the FRW time t (3.4) and is transformed into a function of ξ by multiplying the density function (3.4) bȳ r 2 to obtain (3.31)r 2 ρ = 3 4r 2 κt 2 = 3 4 where ξ =r/t.

The FRW-2 Metric
In this section, we show the existence of another mapping of the FRW metric (3.3) into standard Schwarzschild coordinates. As opposed to the a priori knowledge of the FRW-1 transformation in Theorem 3.1.1, the FRW-2 transformation is developed through construction. There is a classic argument in [Wei72] to build a coordinate transformation (t, r) → (t,r) taking a general spherically symmetric metric (3.1) over to standard Schwarzschild coordinates (3.2), and this argument is repeated to build the FRW-2 transformation. In particular, we take the spherically symmetric FRW metric (3.34) ds 2 = −dt 2 + R 2 (t){dr 2 + r 2 dΩ 2 }, over to the standard Schwarzschild coordinate form We first match the spheres of symmetry (3.36)r 2 dΩ 2 = R 2 r 2 dΩ 2 , so the first component of the coordinate mapping (t, r) → (t,r) is defined by the following (3.37)r ≡r(t, r) = R(t)r.
Rewriting (3.37) as we differentiate and square to obtain the following differential form relationships, Substituting (3.40) into (3.34), results in the FRW metric in (t,r) coordinates In order to complete the coordinate transformation (t, r) → (t,r),t =t(t, r) must be defined to eliminate the cross term dtdr in (3.41). Since it is already matched, ther 2 dΩ term is ignored in our subsequent calculations to simplify notation. We show the procedure for eliminating the cross term for a generalized metric of the form (3.42) ds 2 = −C(t,r)dt 2 + D(t,r)dr 2 + 2E(t,r)dtdr, and handle our specific case afterwards. To eliminate the cross term in (3.42), an integrating factor Ψ = Ψ(t,r) is chosen to satisfy the equation For our specific case, we have and substituting these into the general form (3.47) gives us For the pure radiation phase, where σ = 1/3, the cosmological scale function becomes R(t) = √ t, and the resulting metric is The definition of η (3.18) allows us to rewrite the metric in the barred coordinates as (3.51) Notice this metric (3.51) closely resembles the FRW-1 case (3.20) except for inclusion of Ψ, indicating our ability to recover the FRW-1 metric by choosing the correct integrating factor. With this in mind, we look for solutions of the integrating factor equation (3.43).
In light of the functions C(t,r) and E(t,r) defined in (3.48), equation (3.43) becomes

THE FRW-2 METRIC 25
The following lemma expresses the existence of two solutions for this integrating factor equation Lemma 3.2.1. The following PDE for the function Ψ(t,r) has at least two solutions, a constant solution Ψ 1 (t,r) = Ψ 0 and a dynamical solution where Ψ 0 is some constant.
Proof. The constant function, Ψ 1 (t,r) = Ψ 0 , solves the PDE (3.53) since To find the other solution for the integrating factor Ψ 2 (t,r), a simplification of (3.53) is obtained by noticing so Ψ 2 satisfying (3.53) is equivalent being a solution to the following PDE We proceed by constructing a function Ψ 2 satisfying (3.57). Suppose Ψ 2 has the form where f (η) is an unknown function of the predefined variable η =r/t. Taking the partial derivatives of 3.58 leads to where the prime represents differentiation with respect to the variable η. Substituting η into equation (3.57) gives us Plugging in the partials (3.59) into the PDE (3.60) and simplifying results in the following .
The second integrating factor Ψ 2 , the dynamical solution, has the corresponding 1-form dr .
Again, since this 1-form is exact, there exists a functiont(t,r) such that To solve fort we integrate the ∂t ∂r equation to obtain for some function g(t). Taking ∂ ∂t of (3.71) leads to g(t) = C, the constant function. By setting g(t) = 0, the coordinate transform corresponding to the second integrating factor (3.54) is The radial velocity v ≡ dr/dt is computed as where we used (3.51) that √ AB = 1/Ψ. Remarkably, this relationship matches the FRW-1 case (3.26). Using this relationship, the FRW-2 metric (3.51) is rewritten in terms of the v as (3.76) Unlike the FRW-1 metric, the fluid variables are not a function of the ratio ξ =r/t, so the FRW-2 metric relies on the unbarred coordinate time t. The coordinate transformation (3.72) is used to find t as a function of the new coordinates (t,r), All these results for the FRW-2 metric are recorded in the following thereom The fluid variables (ρ, v) corresponding to this metric are where the unbarred time coordinate is the following function of (t,r) This section ends with the lemma stating any Ψ satisfying (3.43) makes dt (3.44) into an exact 1-form.
Let C be a curve in (t,r) space parameterized by α, denoted by C = (t(α),r(α)), going from the point P 1 = (t(α 1 ),r(α 1 )) to the point P 2 = (t(α 2 ),r(α 2 )), shown in Figure 3.1. Now take the line integral of the 1-form dt along the curve C to obtain Using the chain rule along with (3.85), we have df dα dα = f (t(α 1 ),r(α 1 ))) − f (t(α 2 ),r(α 2 ))). Hence, C dt is only dependent on the endpoints of the curve C which means C dt is independent of path; therefore, the 1-form dt is closed. Since in subsequent chapters we only consider a convex domain D = {r min ≤ x ≤ r max , t ≥ 0}, Poincare's Lemma [Rud76] states a closed 1-form on a convex domain is exact, proving the claim.

The TOV Metric
For the outer TOV solution in our shock wave simulation, we use the general relativistic static isothermal sphere with the equation of state p = 1 3 ρ for the pure radiation phase as derived by Smoller and Temple in [ST95]. For completeness we summerize the results required here. This TOV metric has the form (3.88) ds 2 = −B(r)dt 2 + 1 1 − 2GM (r) r dr 2 +r 2 dΩ 2 .
The time metric component B takes the form where √ σ is the speed of sound, and the mass M function is given as which agrees with equation (3.4) in [ST95]. The fluid variables in the TOV metric are where the velocity is zero since the TOV metric is static.
Recall for both FRW metrics, FRW-1 and FRW-2, the relationship v =r/2t holds. Using this relationship along with the density equation (3.4), the square of the velocity becomes Substituting (3.97) into (3.96) and rearranging produces By taking c = 1 and using equation (3.5) to replace κ, this equation simplifies to The equation (3.99) defines the shock surfacer =r(t) implicitly in the (t,r) coordinates.
To obtain the corresponding curve in the original FRW coordinates (t, r), the substitution ofr = R(t)r into (3.99) is made to define the shock surface r = r(t). Interestingly, equation (3.99) is independent of the integrating factor Ψ, meaning the shock surface is the same curve in (t,r) space regardless of the integrating factor chosen.
To obtain Ψ, the other metric component must be matched at the shock surfacer(t), (3.100) This matching allows us to finish the definition of Ψ by defining the integrating factor constant Ψ 0 once the TOV time scale factor B 0 is chosen. In the simulation setup, we reverse the roles of these constants; we choose Ψ 0 such that the coordinate speed of light

Locally Inertial Godunov Method
This chapter is dedicated to the algorithm for the locally inertial Godunov method featuring dynamical time dilation. Although in this thesis, the method is used to simulate the family of shock waves identified in Chapter 3, the method can be applied to simulate general spherically symmetric flows. This method is a modification of the locally inertial Glimm method by Groah and Temple [JGB06,GT04], and many of the equations within this chapter are taken from their work. In their paper, Groah and Temple devise the fractional Glimm method to prove the existence of shock wave solutions to the spherically symmetric Einstein equations for a perfect fluid. This method is a technique to evolve solutions from initial profiles for the conserved quantities u(t, x) and the metric A(t, x) in standard Schwarzschild coordinates satisfying the Einstein equations. Recall, standard Schwarzschild coordinates take the form (4.1) ds 2 = −B(t,r)dt 2 + 1 A(t,r) dr 2 +r 2 dΩ 2 .
Here, we modify the locally inertial Glimm method by replacing the random choice Glimm step in Groah and Temple's work with an averaging Godunov step to obtain more consistent and less jagged solutions, with the ultimate goal of simulating these shock wave solutions in General Relativity. Also, the purpose of the Groah and Temple construction is to prove an existence theorem; therefore, the development lacks some details needed to numerically construct the solution. A goal of this chapter is thus to fill in these details, like an algorithm to find the middle state to the Riemann problem. Our method can be interpreted as a locally inertial scheme, in the sense it exploits the locally flat character of spacetime. Each grid cell is considered a locally flat frame, and we handle the time dilation between frames, by choosing a reference frame relative to which time can be synchronized. The frame chosen is the one in which the factor √ AB is one. We denote this reference frame as a unitary frame. We also denote the factor √ AB as the coordinate speed of light to represent the rate at which light travels relative to the unitary frame. Throughout this paper, we will refer to it as the speed of light, not to be confused with the speed of light constant c = 1, and it will be clear from context which one is being used. Note that for any single frame, the metric component B (4.1) can be rescaled by a change in the time coordinate to make the (coordinate) speed of light one in that frame, so the reference frame is quite arbitrary and not so important in tracking time. However, it is important what this reference frame's speed of light is relative to the other frames around it, causing time dilation between the frames. These ideas will be explored and expanded throughout the chapter.
In Section 4.1 we discuss the initial value problem in Special Relativity. In this section, we study solving the Riemann problem in Minkowski spacetime, which is a unitary frame.
Solving the Riemann problem is at the heart of the Godunov step of our method, and its importance cannot be overstated. Section 4.2 explains time dilation and how it affects the averaging of our Godunov step, extending our ability to perform this step on nonunitary frames. With the background material set, we conclude in Section 4.3 by stating the algorithm of the locally inertial Godunov method. This method formulates the solution inductively and has four major steps: a Riemann problem step, a Godunov step (with time dilation), an ODE step, and an update step.

The Initial Value Problem in Special Relativity
In this section we develop solutions to the relativistic compressible Euler equations in flat Minkowski spacetime for the case p = σρ with a constant σ together with the initial conditions Note, this equation (4.2) is our conservation law ((2.25) and (2.26)) without the source term where √ AB = 1. This problem is a specific case of the initial value problem for a general system of nonlinear hyperbolic conservation laws in the sense of Lax [Smo83], (4.4) In our case, we have To distinguish between the two sets of variables, we refer to the pair (ρ, v) as the fluid variables and to the pair (u 0 , u 1 ) in (4.5) as the conserved quantities. Both of these variables play an important role in the locally inertial Godunov method. The conserved quantities are needed in implementing the Godunov step, and the fluid variables are needed to solve the Riemann problem along with giving us more physical meaning behind our simulation.
Fortunately, Groah and Temple showed [JGB06] there is a 1 -1 correspondence between the fluid variables and the conserved quantities, giving us the following result Proposition 4.1.1. The mapping (ρ, v) → (u 0 , u 1 ) is 1 -1, and the Jacobian determinant of this mapping is both continuous and non-zero in the region ρ > 0, |v| < c.
In the locally locally inertial Godunov method, we constantly transfer back and forth between the conserved quantities and the fluid variables, so the inversion of the mapping (ρ, v) → (u 0 , u 1 ) defined in (4.5) is necessary to find the fluid variables as a function of the conserved quantities, stated in the following corollary to Proposition 4.1.1 Proof. Using the definition of the conserved quantities (4.5), both are solved as equa- Setting these equations, (4.9) and (4.10), equal to each other and simplifying gives us the In order to determine which solution to use, we normalize the speed of light (i.e. c = 1) and rewrite (4.12) as Since u 0 > u 1 and σ + 1 > σ, the following inequalities hold u 0 (σ+1) 2u 1 σ > 1 and 4(u 1 ) 2 σ 2 (u 0 ) 2 (σ+1) 2 < 1. In order to keep v less than the speed of light (i.e. v < 1), the minus sign is taken in (4.12) to obtain (4.7). Substituting this v into (4.10) gives us ρ as (4.8).
The main focus of this section is to solve the Riemann problem for system (4.2). The Riemann problem is the initial value problem with initial data u 0 (x) of a pair of constant states separated by a jump discontinuity at x = 0, (4.14) Note, in light of Proposition 4.1.1, the conserved quantities u L and u R are uniquely determined by the fluid variables (ρ L , v L ) and (ρ R , v R ), and the symbol u can be interpreted as either set of variables throughout this paper.
A general theorem of Lax [Smo83] states for any system of conservation laws (4.4) which is strictly hyperbolic and genuinely nonlinear in each characteristic field, the Riemann problem for this system has a unique solution in the class of elementary waves if u R and u L are sufficiently close. The proof for this result is a constructive one, relying on the structure of the state space. More specifically, given a point u L in u-space, there exists a family of i-wave curves connecting the point u L to any other sufficiently close point u R by traversing a 1-wave curve followed by a 2-wave curve, where the middle state u M is the intersection of these two curves. Depending on the direction along the 1-wave curve, a shock or rarefaction wave is between the states u L and u M , and similarly for the 2-wave curve, the direction taken determines whether a shock or rarefaction wave is between the states u M and u R . As a side remark, each i-wave curve has second order contact at the point u L (i.e. their first two derivatives are equal at this point) [Smo83]. Figure 4.1 gives us a sample of the web-like structure for these family of curves for a 2 system of conservation laws referred to as the (non-relativistic) p-system [Smo83], for p = 1. In this figure, the red dot represents the left state u L . The 1-shock curve is represented by the blue/red graph, and the 1-rarefaction is represented by the red/cyan graph. The 2-shock and 2-rarefaction curves, originating from the 1-wave curve, are represented by green and yellow graphs, respectively, with the red stripped curves representing the ones emanating from the left state. One can interpret the coordinate system of wave curves defined for each u L as giving us a road map providing directions from our starting point u L to our destination u R : the 1-wave curve is the first street traveled and 2-wave curve next street traveled, and the middle state u M is where the streets traveled intersect. The solution becomes the left state connected by a 1-wave to the middle state which is connected by a 2-wave to the right state. This result is only a local one because for certain conservation laws (including γ-law gases), not all the points in u-space can be connected by this web structure, due to the formation of the so-called "vacuum" states [Smo83].
For other systems, this local connectivity of points in u-space can be extended to a global connectivity of points, where the "vacuum" states do not appear. Fortunately, for and Moreover, the solution is given by a 1-wave followed by a 2-wave, satisfies ρ > 0, and all speeds are bounded by c. This solution is unique in the class of rarefaction and admissible shock waves. The There is a typo in the statement of the Riemann invariants in the Groah and Temple paper [JGB06] in equations (2.5.73), (2.5.74), (4.2.12), and (4.2.13). As well as moving freely between the conserved quantities and the fluid variables, the ability to go back an forth between the fluid variables and the Riemann invariants is needed within the locally inertial Godunov method.
The fluid variables (ρ, v) are recovered from the Riemann invariants (r, s) by algebraic manipulation of (4.17) and (4.18). More specifically, subtracting (4.17) from (4.18) to solve for the density and adding (4.17) to (4.18) to solve for the velocity To solve the Riemann problem, it is necessary to know the speed of the shock waves in the solution to determine which side of the discontinuity the point of interest is located.
The speed for the 1-shock is the following function of beta β while the speed for the 2-shock is These shock speeds, (4.22) and (4.23), are calculated in a frame where the particle velocity v is zero. To obtain these quantities in an arbitrary frame, the Lorentz transformation law for velocities must be applied. The transformation law is given as follows: if in a Lorentz transformation, the barred frame moves with velocity v measured in the unbarred frame with s i as the speed of the i-shock wave measured in the barred frame, and if s denotes the speed of the shock in the unbarred frame, then The eigenvalues λ i (ρ, v) of the system (4.2) are used to determine the speed of the rarefaction waves. These eigenvalues are λ 1 = − √ σ and λ 2 = √ σ when the particle velocity is zero. Using the Lorentz transformation law (4.24) for velocity v, the eigenvalues become (4.25) To find the solution of the Riemann problem within a rarefaction wave, we need to solve for the fluid variables based on the eigenvalues and the Riemann invariants. More precisely, formulas for v as a function of the eigenvalues, ρ as a function of (s, v), and ρ as a function of (r, v) are necessary to solve for a state in a rarefaction wave. A quick calculation on (4.25) and (4.26) shows v is a function of the first eigenvalue or is a function of the second eigenvalue Another quick calculation on (4.17) and (4.18) provides us with ρ as a function of (r, v) and ρ as a function of (s, v) 4.1.2. Finding the Middle State. This subsection describes a numerical algorithm to find the middle state and the associated waves to the Riemann problem for our system (4.2). To accomplish this task, not only does one need to determine in which region the right state is located, governing the set of curve equations to use, but one also needs to solve these curve equations which cannot be solved explicitly. With this in mind, finding this middle state poses quite a challenge with the fluid variables (ρ, v). Fortunately, this problem is easier in the Riemann invariant coordinate system or the rs-plane. This coordinate system simplifies the i-wave curves in the state space because the rarefaction curves become straight lines, resulting in a clearer segregation of the state space. Another advantage is the i-wave curves are geometrically invariant across the rs-plane, so the shape of the wave curves are independent of the base point (r L , s L ), implying an i-wave curve at a point in this plane can be mapped by rigid translation onto any other point. Within this subsection, we deal with the Riemann invariants as our variables, keeping in mind the goal is to obtain the middle state in the fluid variables.
To distinguish between the different coordinate systems, we denote the fluid variables (ρ, v) by lower case u and the Riemann invariants (r, s) by upper case U . More specifically, where the transformations (4.17) and (4.18) are used. In the rs-plane, the 1-shock curve S 1 for the system is given by the following parametrization with respect to the β, 0 ≤ β < ∞: (4.32) and the 2-shock curve S 2 is given by: . In this coordinate system, both rarefaction curves are straight lines parallel to the coordinate axises along the positive directions in the rs-plane. More precisely, the 1-rarefaction curve R 1 is: (4.36) and the 2-rarefaction curve R 2 is given by: for the parameter β, 0 ≤ β < ∞. Notice within the equations for the shock and rarefaction curves (4.32)-(4.37) the differences ∆r and ∆s along the curves only depend on the parameter β, proving the geometric invariance mentioned earlier. For notation purposes, we define the following sets based at a point U L for the shock curves and the rarefaction curves (4.39) Using this knowledge, a simpler family of curves is built to solve for the middle state U M , as seen in Figure 4.2.
Using these curves, the rs-plane is segregated into four disjoint open regions I, II, III, and IV as depicted in Figure 4.3. In each region, a possible right state is placed and denoted by U i for i = 1, 2, 3, 4. This figure also shows the 1-wave and 2-wave curves emanating from U L as the thick lines with the 2-wave curves associated with the points U i for i = 1, 2, 3, 4 emanating from the 1-wave curve as dashed lines. Region I corresponds to connecting the point U L to the point U 1 by 1-rarefaction wave followed by a 2-shock wave. Region II corresponds to the 1-shock/2-shock wave case, while region III and region IV corresponds to 1-shock/2-rarefaction and 1-rarefaction/2-rarefaction wave cases, respectively. Since the wave curves are geometrically invariant, we only consider the change from the right to left states, or the quantities ∆r ≡ r R − r L and ∆s ≡ s R − s L where the right state is U i for i = 1, 2, 3, 4, depending on the region under consideration. The middle state U M for all Family of curves for our system in the rs-plane these cases is found by considering each quadrant in the rs-plane separately. After the quadrant breakdown, we discuss the technique used to find the correct β(s) needed to solve the shock wave equations (4.32) and (4.33). Note that the shock curves cannot be solved exactly, and a threshold of ǫ must be predefined, with the goal of solving the shock curve equations with a maximum error of ǫ.
We start with the point U 4 in +/+ quadrant, which is determined by ∆r > 0 and ∆s > 0. This quadrant is the same set as region IV, the 1-rarefaction/2-rarefaction case.
Since the rarefaction curves are straight lines, the middle state becomes r M = r R and s M = s L .
Next we turn our attention to the point U 3 in the -/+ quadrant, where ∆r < 0 and ∆s > 0. This quadrant is a subset of the region III, the 1-shock/2-rarefaction case. Our goal is to find β such that (r R , s * ) ∈ S 1 (U L ). After finding β, the middle state will be The point U 1 in the +/-quadrant, where ∆r > 0 and ∆s < 0, is similar to U 3 in the -/+ quadrant. This quadrant is a subset of the region I, the 1-rarefaction/2-shock case. Due to the geometric invariance of the 2-wave curve, we find a β such that the point (r * , s R ) ∈ S 2 (U L ), and the middle state will be (r M , s M ) = (r L + (r R − r * ), s L ).
Finally, consider the point U 2 in the -/-quadrant, where ∆r < 0 and ∆s < 0. This quadrant is a superset of region II, the 1-shock/2-shock case, so a mechanism is needed to determine when U 2 is not in region II but in region I or III as discussed below in the two shock algorithm. If U 2 is in region III, we find β 1 such that U M = (s M , r M ) ∈ S 1 (U L ) and β 2 such that (s R , r R ) ∈ S 2 (U M ). The middle state U M is, the intersection of the 1-shock curve S 1 (U L ) and the 2-shock curve S 2 (U M ) crossing the right state U R .
When solving one shock curve equation, as in the +/-or -/+ quadrant, a standard bisection method [BF97] is used to find β. We show how we implement the bisection method for the point U 3 , where a 1-shock is followed by a 2-rarefaction. In this case, the bisection method is used to find the β that satisfies the equation ∆r(β) ≡ r R − (r L + S r 1 (β)) = 0, and our goal is to find a β such that ∆r(β) < ǫ. The bisection method requires a starting interval [β min , β max ] where there exists a β within this interval such that ∆r(β) = 0. The shock curves, S r i (β) and S s i (β) for i = 1, 2, are monotone decreasing functions of the variable β. To find our interval, an initial guess, like β = 10 5 , is chosen and ∆r(β) is computed. If ∆r(β) > 0, our guess is too small (i.e. r R > r L + S r 1 (β)), and β is decreased. The power of our initial guess β = 10 k is decreased (i.e. k = 4, 3, 2, 1, 0, −1, . . .) until a k is found such that ∆r(10 k ) < 0, which means the guess is too big now. With this k, we set β max = 10 k+1 and β min = 0. On the other hand, if for the initial guess ∆r(β) < 0, it is too big (i.e. r R < r L + S r 1 (β)), and β needs to be increased. The power of our initial guess β = 10 k is increased (i.e. k = 6, 7, 8, 9, . . .) until k is found such that ∆r(10 k ) > 0, which means our guess is too small now. With this k, we set β max = 10 k and β min = 10 k−1 . Either way, with the interval [β max , β min ] established, the bisection method is implemented until a β is found where ∆r(β) < ǫ. For the -/+ quadrant, the same algorithm is used, but for the is sought out, and our goal is to find (β 1 , β 2 ) such that ∆r(β 1 , β 2 ) < ǫ and ∆s(β 1 , β 2 ) < ǫ. To perform the bisection algorithm for both β 1 and β 2 , intervals [β 1 min , β 1 max ] and [β 2 min , β 2 max ] are needed, just like solving for one shock equation. Since the equation (4.40) is dominated by the parameter β 1 , in order to get our interval [β 1 min , β 1 max ], we assume β 2 = 0 simplifying (4.40) into ∆r(β 1 ) ≡ r R − (r L + S r 1 (β 1 )) = 0 and repeat the same procedure as above, increasing or decreasing the power k in our guess 10 k , for the parameter β 1 . Similarly, we repeat this process to find [β 2 min , β 2 max ], by setting β 1 = 0 in (4.41) and working with the equation ∆s(β 2 ) ≡ s R −(s L +S s 2 (β 2 )) = 0. If either initial guess 10 k gets too small, like k < −20, then we assume the point U 2 is in the wrong region. Depending on which parameter it is, β 1 or β 2 , the U 2 must be in region I or region III, respectively. In either case, our problem boils down to only solving one shock curve equation, and it is handled according to the above procedure. Equipped with our initial guess (β 1 , β 2 ), we compare the errors ∆r(β 1 , β 2 ) and ∆s(β 1 , β 2 ) against each other. If ∆r(β 1 , β 2 ) < ∆s(β 1 , β 2 ), we perform a bisection method step on β 1 ; otherwise, we perform the step on β 2 . We repeat this process until ∆r(β 1 , β 2 ) < ǫ and ∆s(β 1 , β 2 ) < ǫ.
4.1.3. Solving the Riemann Problem. Being able to find the middle state with the connecting waves enables us to build the solution to our Riemann problem in its entirety.
In particular, for an arbitrary point (t, x), the solution u(t, x) can be determined for any Riemann problem. This subsection covers the details of finding this solution at an arbitrary point. We also use this process to build a Riemann problem simulator.
Suppose we have a Riemann problem (4.14) to our conservation law (4.2) and a point (t, x), and we want to solve for the state u * ≡ u(t, x) = (ρ * , v * ) as a solution to this problem. Converting u L and u R over to Riemann invariants, using (4.17) and (4.18), and implementing the above procedure, the middle state U M is determined along with the 1wave connecting to the left state and the 2-wave connecting to the right state. Converting Without loss of generality, we consider the 1-shock/2-rarefaction case, where the other cases are handled similarly. Since the solution to the Riemann problem is self-similar, the ratio of our point of interest s ≡ x/t is compared to the wave speeds, usually done left to right. We first compare this ratio against the 1-shock speed s 1 , computed by the transformation law (4.24) with v = v L and the speed (s i ) computed by (4.22) using the β found in the middle state algorithm. If s < s 1 , then the state proceeds the 1-shock, and u * = u L is the solution. Otherwise, we compute the speeds for the left and right sides of the 2-rarefaction fan. The left speed s L 2 is computed by (4.26) with u M as the input, and the right speed s R 2 is computed similarly using u R instead. If s 1 < s < s L 2 , then the state is between the shock and rarefaction waves, and u * = u M . If s R 2 < s, then the state is beyond the rarefaction wave, and u * = u R . If s L 2 < s < s R 2 , then the state is within the rarefaction wave, and more work needs to be done to find the solution. We use the fact that within the rarefaction wave the solution varies smoothly, and every state u between u M and u R moves at the speed λ 2 (u) [Smo83]. Since the speed for our state is λ 2 (u * ) = x/t, the inversion of the second eigenvalue (4.28) gives us the fluid velocity v * . With the fluid  as a yellow dot, by the 1-shock curve. Also, the middle state is connected to the right state u R , shown as a red dot, by the 2-rarefaction curve. The other two panels show the density and velocity profiles after one unit of time. Here one can see, going left to right, the left state, the 1-shock wave, the middle state, the 2-rarefaction wave, and the right state in these graphs.

Time Dilation Between Space Time Cells
With the capability of solving the Riemann problem for the relativistic compressible Euler equations in a unitary frame, like Minkowski spacetime, this section is dedicated to only changes the speeds of the waves, and when considering solutions to the Riemann problem, it has no effect on the states themselves [JGB06]. Since the locally inertial Godunov method has a fixed spatial distance for all the frames, this factor determines how time is sped up or slowed down relative to the unitary frame. One way to view this phenomenon is this factor stretches or shrinks the height of our spacetime cell relative to other frames, as illustrated in To proceed, we need to distinguish between two different spacetime cells. In the last section, we dealt exclusively with a cell containing the whole Riemann problem posed, which we will refer to this cell as a Riemann cell. Now, we consider a cell containing an entire    where λ = ∆t ∆t < 1 is the ratio between the two times.
Proof. Suppose we have a solution to both Riemann problems u(t, x) in the Godunov cell, along with a maximum time ∆t (4.42) and a shorter time ∆t < ∆t. Since each half of the Godunov cell contains a disjoint self-similar Riemann problem, we prove the claim on one half of the grid cell. Without loss of generality, we study the left half (i.e. 0 < x < ∆x/2) with the origin located at the bottom left of this cell. We start by constructing the relationship between the function u(t, x) at the two times. The fastest speed in the Godunov cell (i.e. the speed of light) is referred to as α = √ AB. By the self-similarity of the Riemann Problem, the following relationship holds The average at the time ∆t is directly calculated as (4.45) Since ∆x 2 = α∆t, this implies and (4.45) becomes proving the claim.
Since the Godunov step only considers the zero speed states u L * and u R * and the factor √ AB has no affect on scaling a speed of zero, we expect the affect on the Godunov step of this time dilation to be minor. A quick calculation for the left half of the Godunov cell (4.48) Hence, the only effect of a reduced time is to shorten the input time in the Godunov step.

Locally Inertial Godunov Method
Equipped with all the necessary tools, we state the algorithm of the locally inertial Godunov method. This method is a fractional step scheme started by choosing some parameters. In particular, we choose a minimum radius r min , a maximum radius r max , the number of spatial gridpoints n, and a start time t 0 . In our simulations, the number of spatial grid points n is chosen to be a power of two (i.e. n = 2 k for some k). From these parameters, the mesh width ∆x is determined to be and is fixed throughout the scheme. Let (x i , t j ) represent a mesh point in an unstaggered grid defined on the domain The spatial points are defined as (4.51) Unlike the mesh width, the time step or the mesh height, ∆t, changes from one time step to the next because there is no way to determine beforehand the smallest ∆t satisfying the where the minimum is taken over all the spatial gridpoints at time t j of the metric A ij = (A ij , B ij ), where these entries are defined shortly. Starting at t 0 , our temporal mesh points are defined as We assume at our current time t j for j ≥ 0 there exists a solution u(t j , x) and A(t j , x) This solution is either provided as the starting solution at t 0 or from the last iteration of the locally inertial Godunov scheme constructed inductively. To implement the method, this solution is discretized into piecewise constant states. Discretizing the conserved quantities u(t j , x), let u ∆x be given by piecewise constant states u ij at time t = t + j as follows: For notational convenience, we denote x i+ ≡ x i+ 1 2 and x i− ≡ x i− 1 2 throughout this paper.
We define the grid rectangle R ij so the mesh point (x i− , t j ) is in the bottom center of it, which is diagrammed in Figure 4.8. Each grid rectangle is a Riemann cell, containing a solution to a distinct Riemann problem. We are limited to solving Riemann problems within Riemann cells having a constant speed of light, as discussed in the last section.
To this end, the metric source A = (A, B) must be approximated by a constant value, denoted A ij , in each Riemann cell R ij throughout the simulation. These constant values are established by setting This approximation makes A ∆x discontinuous along each line x = x i , i = 1, . . . , n, at each To implement the Godunov step, we need boundary profiles at the left and right boundaries along with the initial profiles at time t 0 because the Godunov step is a three point method, and the points x 1 and x n need left and right partners, respectively, to pose the boundary Riemann problems. These boundary profiles are used to implement the boundary Riemann problems and are referred to as ghost cells. The left and right ghost cells, located at the points x 0 and x n+1 , respectively, must be consistent with our numerical solution to the Einstein equations around these boundaries. Any inconsistency in these boundary conditions would result in errors propagating into our solution, corrupting the data; therefore, data is needed for the left ghost cell u 0,j and A 0,j along with the right ghost cell u n+1,j and A n+1,j that are solutions to the Einstein equations synchronized with the data close to the boundary. Figure 4.9 displays the location of these ghost cells.
The discontinuities of the metric A ∆x are staggered relative to the approximate solution u ∆x as illustrated in (4.57) Let u RP ij (t, x) denote the solution of (4.57) within the Riemann cell R ij , and define as the Riemann problem step of the fractional step scheme.
Equipped with the solutions to the Riemann problems in each Riemann cell R ij , we implement the Godnunov step to obtain the average of fluid variables u RP ∆x across the intervals [x i− , x i+ ] at the next time step t j+1 . Since the metric A is different on both sides of x i , separate averages must be taken over the left and right half cells and combined to obtain the true average. In particular, letū L ij andū R ij be the average on the left and right half cells, respectively. Also, let u L x i+ ) represent the zero speed states left and right Riemann problem, respectively, as shown in Figure 4.6.
To perform the Godunov step on the left half cell, we compute and do the same for the right half cell where the 2 accounts for the half cell calculations. Taking the average of these results leads to (4.61)ū ij = 1 2 {ū L ij +ū R ij }, defining our Godunov step of the method.
We proceed to define the ODE step. Letû(t, u 0 ) denote the solution to the following initial value problem We define the approximate solution u ∆x (t, x) and A ∆x (t, x) analytically to derive the piecewise formulas used to update the numerical scheme and to be used in the convergence proof of the next chapter. The conserved quantities are defined by the formula Therefore, u ∆x (t, x) is equal to u RP ∆x (t, x), the solution to the Riemann problems, plus a correction term from the ODE step of the method. The metric is derived from the definition of the mass In terms of these equations, define the metric as T 11 M (u ∆x (r, t)) dr.
Finally, in order to update the metric and conserved quantities, we use the Riemann problem averagesū ij to replace the Riemann problem solution u RP ∆x (t, x) and perform numerical integration on the analytical equations (4.65)-(4.68). This process leads us to define The mass is and the metric becomes Note that since the metric is staggered relative to the conserved quantities, we use the in between values, like x k− and u ∆x (x k− , t j+1 ) in the update step. Let A i,j+1 = (A i,j+1 , B i,j+1 ) denote the constant value for A ∆x on R i,j+1 . This concludes the update step and completes the definition of the approximate solution u ∆x and A ∆x by induction.
To summarize the method, after the setup of the gridpoints, Riemann cells, initial profiles, and ghost cells, the locally inertial Godunov method constructs the solution inductively with four major steps: a Riemann problem step, a Godunov step (with time dilation), an ODE step, and an update step. The Riemann problem step is described in equations (4.57) and (4.58). Formulas (4.59)-(4.61) denote the Godunov step. The ODE step is detailed in

Convergence of the Method
The focus of this chapter is proving the main theorem of this thesis. This theorem states that if a solution (u ∆x , A ∆x ) → (u, A) using the locally inertial Godunov scheme converges and has a total variation bound at each time step, then (u, A) are weak solutions to which is weakly equivalent to the Einstein equations (2.9)-(2.12). This proof is a modification of the Groah and Temple argument using the locally inertial Glimm scheme [JGB06], with a few differences. The main difference is the solution update at each new time step.
In this paper, an average of the fluid variables is taken verses a random sampling, but the steps leading up to this update step are the same. Another difference is the assumed total variation bound is used to bound the Riemann problem solutions, as opposed to Groah and There are two main points to the proof. The first point is to show the discontinuities in the metric A along the boundary of Riemann cells are accounted for by the inclusion of the term in the ODE step (4.62). The second point is to show the jump in the approximate solution u ∆x along the time steps are of order ∆x. In their work [JGB06], Groah and Temple did not need the convergence and total variation assumptions because with the Glimm scheme, these assumptions are proven as long as there exists a total variation bound on the initial data, a truly remarkable feature of the scheme. In this thesis, these assumptions are shown numerically by the simulation results in the subsequent chapters. Thus, our theorem is perfectly suited to our numerical simulation: we numerically establish convergence and a total variation bound, and the theorem of this section proves that assuming these conditions, we can conclude convergence to a weak solution of the Einstein equations.

Convergence to a Weak Solution
The main theorem of this thesis is the following Theorem 5.1.1. Let u ∆x (t, x) and A ∆x (t, x) be the approximate solution generated by the locally inertial Godunov method starting from the initial data u ∆x (t 0 , x) and A ∆x (t 0 , x) for t 0 > 0. Assume these approximate solutions exist up to some time t end > t 0 and converge to a solution (u ∆x , A ∆x ) → (u, A) as ∆x → 0 along with a total variation bound at each any time t such that t j ≤ t < t j+1 . By the self similarity of the solution to the Riemann problem, this result also implies a total variation bound for any space coordinate within the Riemann cell. More specifically, we have the following bounds: for any x and t within the Riemann cell R ij .
All the functions f , G, and g derived in [JGB06] are smooth, and it is the metric that is only Lipschitz continuous. The smoothness of these functions is used throughout this proof.
Let This variable is created to unify all the time steps, and more importantly, used to calculate the maximum number of time steps needed to go from t 0 to t end .
We now follow the development of Groah and Temple in [JGB06]. Recall, u RP ∆x (t, x) denotes the collection of the exact solutions in all the Riemann cells R ij for the Riemann problem of the homogenous system So u RP ∆x (t, x) satisfies the weak form of this conservation law in each Riemann cell where ϕ is a smooth test function with Supp(ϕ) ⊂ [t 0 , t end ) × [a, b] for a < r min < r max < b.
Remember,û(t, u 0 ) denotes the solution to the ODÊ (5.9) Therefore, Also, recall u ∆x denotes the approximate solution obtained using the fractional step method.
Since our fractional method takes the Riemann problem solution and feeds it into the ODE step, u ∆x is defined on every Riemann cell R ij as This expression implies the error between the approximate solution and the Riemann problem solution is on the order of ∆x; a fact that is repeatedly used throughout the proof.
Define the residual ε = ε(u ∆x , A ∆x , ϕ) of u ∆x and A ∆x as the error of the solution in satisfying the weak form of the conservation law (5.1) by where (5.13) (5.14) The expression i=n+1 i=1,j denotes a double sum where the index i runs across all the spatial gridpoints, and the index j runs across all the temporal gridpoints. Remember, n is the number of spatial gridpoints, and there are n + 1 Riemann cells as depicted in Figure   4.9. Our goal is to show ε(u ∆x , A ∆x , ϕ) = O(∆x) because if the approximation converges (u ∆x , A ∆x ) → (u, A) as ∆x → 0, then the limit function satisfies the condition of being a weak solution to the Einstein equations ε(u, A, ϕ) = 0.
Substituting (5.11) into (5.12) gives us Plugging the weak form of the conservation law (5.8) of each grid rectangle into (5.15) gives where the number of time steps is proportional to T /∆t and the number of space steps is O(1/∆x) by (4.49).
Since u RP ∆x (t + j , x) = u ∆x (t + j , x), the following sum is rearranged to become where the term u ∆x (t j , x) is added and subtracted to isolate the jump in the solution u ∆x across the time step t j . We define this jump ε 1 = ε 1 (u ∆x , A ∆x , ϕ) as and this definition allows us to rewrite (5.17) as But the last sum is rearranged to cancel the boundary conditions as follows: and similarly (5.25) Note that the resulting double sum in (5.23) lost a term, resulting in only n terms.
To simplify the I 1 ij term, we add and subtract a term deviating from it by an order of ∆x, use integration by parts on the new term, and with the result add and subtract another term to reduce the expression further. To this end, let (5.26) From the total variation bound on the Riemann problems and the smoothness of f , this term is bounded by (5.27) and the above procedure reduces the term to where (5.29) and (5.30) Again by smoothness and the total variation bound, we have (5.31) Substituting (5.23) and (5.28) into (5.22) along with using (5.11) as an identity leaves us The second sum represents the jump in the flux function f , resulting from the discontinuities in the metric A, and the first sum is the addition to the ODE step (5.9) specifically designed to cancel these jumps in the flux.
To see how the cancelation works, we perform a Taylor expansion on the test function, and we add and subtract terms deviating by order ∆x. The first sum in (5.32) is expanded (5.34) Plugging this result (5.34) into (5.32) gives us where one term remains due to the mismatch in the number of terms in the spatial sum.
Clearly, this last term is O(∆x).
So the residual boils down to with all that remains to show is To estimate ε 1 , we break up the sum by each time step t j and define Recall, the approximate solution for the new time step t + j is computed by the Godunov step, using averages at the top of each Riemann cell R ij . In particular, the solution at each new time step is To finish the proof, a lemma is needed, which is proven at the end of this chapter. This lemma states the difference of the ODE step taken on an average verses the solution to the Riemann problem across the top of the Riemann cell is bounded by the total variation of the Riemann problem.
Lemma 5.1.1. Let u RP ∆x represent the solution of the Riemann problem in the Riemann cell R i,j−1 andū ∆x (t) denote the average of the Riemann problem solution across Riemann cell. Letû be the solution obtained by the ODE step (5.9) and ϕ be a smooth test function.
Then the following bound holds for some constant C.
Using Lemma 5.1.1, (5.38) is rewritten as solutions to the ODE step (5.9) and bounded by By the total variation bound on u ∆x (t j , ·), the residual is bounded by Therefore, ε = O(∆x) and the proof is complete.
To prove Lemma 5.1.1, a preliminary result is needed: given a function on a set of points the difference of the function between any point and the average is bounded by the total variation of that function on the set. This result is provided by the following x i+ be the average of u on this set. Then we have Proof. The second inequality is true by the definition of the total variation (5.46) sup To prove the first inequality, we assume it is false to obtain a contradiction, so suppose Relabel the u-coordinates by an isometry ϕ : u → v that maps the point u(x * ) to the origin in the v-coordinates (i.e. ϕ(u(x * )) = 0), and the vectorū − u(x * ) in the direction of the 1st coordinate v 1 , as show in Figure 5.1.
Since the average of a collection of points is independent of the coordinate system in which they are labeled in, we have The following inequality holds by transforming equation (5.47) over to v-coordinates This inequality |v| < |v| is an obvious contradiction, proving the first inequality in (5.45). Now we prove the lemma used in the proof of Theorem 5.1.1 Proof of Lemma 5.1.1 Recall the solution to the ODE step has the form: This solution implies the LHS of (5.41) is written out as x i+ where the test function in the first term is approximated by a Taylor expansion. By the definition of the average functionū, the first term is zero. By the smoothness of G, the bound (5.41) is proven by where Lemma 5.1.2 is used to bound the difference between the average and the solution to the Riemann problem.

Continuous Models
Before simulating the shock wave models, we briefly explore simulating more simple ones, which we denote as the continuous models. The continuous models are the pure FRW-1, FRW-2, and TOV metrics developed in Chapter 3, and these models are not as complex as the discontinuous shock wave models. Our purpose in simulating these models is to use exact formulas to test the numerical convergence and the accuracy of our locally inertial The models we consider in this paper have infinite extent, with the radii going from zero to infinity, and we refer to this infinite region of space as the universe of our simulation.
Even though theoretically space can extend out to infinity, the computer can only simulate a finite region of spacetime in the model, so we need to demarcate the minimum radiusr min and the maximum radiusr max of the simulated region, along with the number of gridpoints n to simulate it. Also, we need to decide when to begin the simulation of our models, so a start time t 0 is chosen. We choose to set these parameters as This chapter is organized into four sections. The first three sections, Sections 6.1-6.3, are dedicated to simulating the FRW-1, TOV, and FRW-2 metrics, respectively. These sections start out by reiterating the equations developed in Chapter 3 that are needed to run the simulation. More precisely, these equations are used to build the initial profiles and ghost cells. With this data set, we use our locally inertial Godunov scheme to simulate these models and provide a glimpse of them with various pictures. In Section 6.4, we discuss how the errors and convergence rates are computed, and we record the numerical convergence in the simulation of all three continuous models.

Simulating the FRW-1 Model
Recall, using the coordinate transformation where the variable ξ is defined as ξ =r/t and is a function of v by the following relation (6.6) ξ = 2v 1 + v 2 .
The integrating factor constant is suppressed, which is equivalent to setting Ψ 0 = 1.
With this parameter set, the initial profiles and ghost cells for the simulation are constructed by equations (6.3)-(6.6) along with the standard parameters (6.1). The resulting profiles are illustrated in Figure 6 The simulation is run for one unit of time, t end = t 0 + 1, with the results depicted in Since a positive velocity indicates matter is moving outward, we expect this decrease in density corresponding to a decrease in the mass function. The metric has little change, and moreover, the solution as a whole has a very similar shape from the initial profiles seen

Simulating the TOV Model
Remember, the TOV metric takes the form (6.7) ds 2 = −B(r)dt 2 + 1 1 − 2GM (r) r dr 2 +r 2 dΩ 2 , Setting the extra parameter B 0 = 1 and using the standard parameters (6.1), the initial profiles and ghost cells are built from these equations (6.7)-(6.10) (with σ = 1/3), and the resulting profiles are shown in Figure 6 Running the simulation for one unit of time, the evolution of the TOV metric along this time frame is shown in Figure 6.4. This metric is static, so we expect it to remain unchanged throughout this simulation. The density, metric B, and the mass profiles are unchanged. The constant profiles of the velocity and metric A appear to have changed. This appearance is due to the small numerical errors and the auto-zoom feature of the graphing tool which magnifies these errors.

Simulating the FRW-2 Model
Remember, using the coordinate transformation r = √ tr, t = Ψ 0 2 4t 2 +r 2 t , (6.11) the FRW metric in standard Schwarzschild coordinates (the FRW-2 form) is (6.12) with the metric components (6.13) the integrating factor (6.14) Ψ(t,r) = Ψ 0 t 4t 2 +r 2 , This relationship is explained in more detail in the next chapter and is a rearrangement of equation (7.49). The initial profiles and ghost cells for the simulation are created by equations (6.12)-(6.16), using the standard parameters (6.1) and Ψ 0 set in (6.17). These profiles are pictured in Figure 6.5. In this figure there is more curvature in the graphs of FRW-2 than the FRW-1 model (Figure 6.1). This curvature produces higher values in all the increasing graphs (density, velocity, metric B, and mass) and lower values in the decreasing graphs (metric A).
Allowing the simulation to run for one unit of time, the evolution of the FRW-2 metric is shown in Figure 6.6. Again, the FRW-2 metric follows the same overall pattern as the FRW-1 metric (Figure 6.2), and the fluid variables are decreasing as time progresses. We expect these similarities since both metrics, FRW-1 and FRW-2, are derived from the same FRW metric. In the FRW-2 universe, the (coordinate) speed of light changes uniformly across the simulated region from 1 at t 0 to 1.0667 at t end . Thus, the FRW-2 model has a dynamic (coordinate) speed of light. In contrast to the FRW-1 and TOV models, our simulation of the FRW-2 model tests the locally inertial Godunov method in a setting in which the (coordinate) speed of light is dynamic.

Convergence Results
This section discusses the numerical convergence results for all three of the continuous models. These results are recorded for the FRW-1, TOV, and FRW-2 metrics in Table   6.1, 6.2, and 6.3, respectively. The tables are organized to show, from left to right, the number of gridpoints used, the density results, the velocity results, the metric A results, and the metric B results. The variable results are partitioned into two values: the onenorm error and the convergence rate. Since we possess exact formulas for the solution of all three models, the 1-norm error is numerically computed between the numerical solution and this exact solution. The 1-norm is a natural norm to use for conservation laws because it requires integrating the function, and the weak form of the conservation law gives us information about these integrals. The 1-norm is our chosen method for computing the error and showing numerical convergence in this paper. The convergence rate is computed by taking the log 2 of the ratio in successive errors, enabling us to measure the decrease in error relative to the increase in the number of gridpoints. For example, a rate of 1 means that using twice the number of gridpoints reduces the error by half. A rate less than 1 means the error is reduced by less than a half, while a rate greater than 1 means the error is reduced by more than a half. Looking at the tables, all the errors are decreasing as the number of gridpoints increase with convergence rates around 1. The convergence rate of 1 is expected because we are implementing a first order method on continuous solutions.
These results indicate the simulation is producing an accurate numerical representation of all three models, giving us the green light to simulate the discontinuous models in the following chapters.

Shock Wave Model
In this chapter we perform the numerical simulation of a spherically symmetric general relativistic outgoing shock wave, together with the secondary incoming reflected wave. This simulation demonstrates the secondary reflected wave is also a shock wave. To reiterate, we To simulate the secondary reflected wave and develop a picture of these solutions, we assume both TOV and FRW spacetimes satisfy the same equation of state, which is p = c 2 3 ρ, for the pure radiation stage of the early universe. To run the simulation, we use an exact formula for the FRW spacetime in standard Schwarzschild coordinates first constructed in [ST], detailed in Chapter 3. This transformation puts the FRW metric into the "same coordinate system" as the TOV metric, enabling us to match these metrics together, with the FRW metric on the inside, to form the FRW/TOV matched metric. Since we have two transformations of the FRW metric into standard Schwarzschild coordinates, denoted as FRW-1 and FRW-2, we have two matched models, FRW-1/TOV and FRW-2/TOV, at our disposal to simulate. This matching gives us exact expressions for the initial data (a point of continuity between FRW and TOV metric components, such that the density, pressure, and velocity suffer a jump discontinutiy) and boundary data. This data is the information necessary to run the locally inertial Godunov scheme, developed in Chapter 4, to perform the numerical simulation of the interaction region between the two spacetime metrics. In particular, since light speed is a speed limit for propagation of signals, outside this lightcone the solution should remain FRW (on the inside) and TOV (on the outside). We call this region contained within the light cone that emanates from the initial discontinuity the region of interaction. Besides simulating shock wave formation, we also compute the cone of light and cone of sound, the light and sound information that propagates outward from the initial discontinuity. Both the cone of light and the cone of sound emanate and expand away in both directions from the initial discontinuity. Because light travels faster than sound, the cone of sound must be contained in the cone of light. Since we expect from [JGB06] that there is no lightlike propagation in spherically symmetric spacetimes, the cone of sound should be the true limit of propagation of information; therefore, outside the cone of sound we expect the solution should remain FRW (on the inside) and TOV (on the outside). Figure 7.1 illustrates all these expectations. This prediction is an important result we demonstrate in our numerical simulation. In fact, since the numerical simulation involves integration of the metric components at each time step, we know of no mathematical proof that the weak solutions constructed by the Groah and Temple theorem [JGB06] actually have the sound speed as a propagation speed for the information in the solution. Just like the continuous models, we consider the entire infinite region of the FRW/TOV model as the universe of our simulation, and we need to define initial parameters to run the simulation. Recall, our chosen parameter to the one parameter family of shock wave solutions is the radius of the initial discontinuityr 0 separating the two metrics. We choose to set these parameters as (7.1)r min = 3,r max = 7,r 0 = 5, n = 2 14 = 16, 384, in all the simulations in this chapter, unless otherwise stated. The units to these are explained in Chapter 9. Also, throughout this chapter, when we discuss mesh points, we use the space coordinate x to match the notation of the locally inertial Godunov scheme as denoted in Chapter 4, and when we discuss points within our universe, we use the radial coordinate r to explicitly express these points as a radial coordinate of the metric. When we discuss the derivative of a function in our solution, we compute these derivatives numerically using a standard three point method [BF97].

Region of Interaction
In Section 7.1 we introduce material needed for testing the accuracy and certainty of computing the correct solution in our numerical simulation. We start by computing the (coordinate) speed of light and sound. We move to tracking the borders for both the FRW and TOV side of the simulation, and we conclude this section by finding the boundary condition for the TOV side. We setup the FRW-1/TOV matched model simulation in Section 7.2, determining the initial profiles and boundary conditions we can discretize and feed into the locally inertial Godunov scheme. We use this setup in Section 7.3 to run the simulation and obtain results. We record the convergence of the entire simulation, using the successive mesh refinement technique introduced by Colella and Woodward in [WC84].
We determine the interaction region is not only contained in the cone of sound, but it is exactly synchronized with the cone of sound region. We also show convergence of the non-

Preliminaries
In our simulation, we track the outer boundary of the lightlike and sound like information emanating from the initial discontinuity and compare them to the trajectory of the shock waves. In order to make this comparison, we need to determine how fast light and sound propagates through our universe. Since space and time are bent differently across the universe, we expect these speeds to be dependent on the point of interest. We first compute the (coordinate) speed of light. To simplify the calculation, we suppress ther 2 dΩ 2 term and take the 2-dimensional metric in standard Schwarzschild coordinates Consider a lightlike curve, γ(ξ) = (t(ξ),r(ξ)), which has the velocity vector Since γ(ξ) is a lightlike curve, the length of this curve is zero or Solving this equation for the coordinate speed, we obtain the speed light travels as taking the plus or minus sign for light traveling in the positive or negative direction, respectively.
To find the speed at which sound travels in (t,r)-coordinates, we first note the sound speed is tied to the equation of state p = σρ as being √ σ. This speed is relative to the background fluid at rest (i.e. v = 0). To find the speed of sound relative to fluid moving at a speed v, we need to apply the Lorentz transformation law for velocities (4.24). More specifically, the speed √ σ transforms into Again the plus/minus sign determines the direction in space being traveled. In order to determine the (t,r)-coordinated speed, we magnify this speed by the factor √ AB, and the speed of sound becomes Notice how the speed of sound is not only dependent on the metric (A, B), but also on the movement of the fluid v, indicating the dependence on the medium the sound travels through. where we set the speed as s = l ± or s = s ± for the light or sound information, respectively, choosing the plus/minus sign depending on the direction traveled.
For the cone of sound verses the shock wave position comparison and error estimates, we want to know where the FRW metric stops and the interaction region begins; therefore, we need a mechanism to determine the border between the two. Also, we need to find a similar border between the TOV metric and the interaction region. Theoretically, the position of the shock wave can be used as the borders, but since the solution is numerically simulated with a first order method, these shock waves are smeared by numerical diffusion. This diffusion bleeds into the metrics, causing numerical error in testing this region against the known model. In order to handle these errors, we develop a criteria, based on studying the numerical solution, to determine where the numerical diffusion ends and the metric begins.
The border criteria we choose for the FRW side is where the spatial derivative of the fluid velocity first changes sign. In terms of notation, the border is the first positionr such that (7.9) ∂v ∂r (r) ∂v ∂r (r + ∆r) < 0 This criteria is sufficient because the velocity is an increasing function ofr for the FRW metric in standard Schwarzschild coordinates and only decreases when it hits the shock.
To implement this criteria, we start the gridpoint under consideration at the minimum radius x i = x 1 =r min , which is in the FRW metric because of the boundary condition, and increase the index i until the numerical derivative changes sign, giving us the point where the numerical diffusion first takes place.
Next, we explore the border criteria for the interaction region and the TOV metric.
Again, the fluid velocity is used as our indicator. Since the fluid velocity is theoretically zero, we would ideally find the point where it first becomes non-zero, but in the simulation it is only approximately zero, due to numerical error. Instead, the border criteria we choose for the TOV side is where the absolute value of the spatial derivative of the fluid velocity first becomes greater than 0.01. More specifically, the border is the first positionr such that (7.10) ∂v ∂r (r) > 0.01.
This criteria detects the first significant change in the velocity and ignores small changes occurring from numerical error. To implement this strategy, we start the gridpoint at the maximum radius x i = x n =r max , which is in the TOV metric, and decrease the index i until the numerical derivative satisfies (7.10).
To run our simulation, we need boundary data from both the FRW and TOV metrics to maintain a consistent solution at the edges of the simulated region of the universe.
Much to our surprise, the TOV metric boundary condition changes although the TOV metric itself is time independent. Even though the theory predicts there is a TOV metric matched continuously to the right side of the interaction region, the TOV metric allows for arbitrary changes of the time coordinatet, and we do not know which of these is matched in the solution except at the initial matching by solving for the correct B 0 . More precisely, consider an arbitrary time coordinate map ϕ : τ →t applied to the TOV metric (7.11) ds 2 = −B(r)dt 2 + 1 1 − 2GM (r) r dr 2 +r 2 dΩ 2 .
Since this metric is time independent, the effect oft = ϕ(τ ) is just a scaling of the B metric component, resulting in the transformed metric of (7.12) In (τ,r)-coordinates, the metric is no longer time independent because the time metric which TOV metric is used as the boundary condition. This rematching also enables us to pick out the correct TOV metric model to test against the simulated solution in order to prove convergence of the TOV side of our universe, which will be discussed in more detail later.
To explain this in different words, the theory in [JGB06] tells us the reason we must

FRW-1 and TOV Matched Model Setup
In this section, we cover the details of building the initial profiles and the boundary data required to setup the FRW-1/TOV model simulation. We restate both metrics for ease of reference. Using the coordinate transformation (7.14) the FRW metric in standard Schwarzschild coordinates (the FRW-1 form) is with the metric components (7.16) and the fluid variables where the variable ξ is defined as ξ =r/t and is a function of v by the following relation Whereas, the TOV metric is written where the parameter γ is a constant dependent on σ (7.22) γ = 1 2πG Recall, for this simulation we let σ = 1/3.
Assuming an initial start timet 0 > 0 to be determined later, we choose an initial radius for the discontinuity ofr 0 , and we want to match the metric components (A, B) continuously at the starting point (t 0 ,r 0 ) for the initial discontinuity in the fluid variables. We start by matching the metric A component at this point Let v 0 = v(r 0 /t 0 ) represent the fluid velocity at the discontinuity so (7.23) implies where we substituted (7.22) for γ. Take note that v 0 is independent of the our free parameter r 0 , it is quite astounding the velocity of the fluid at the discontinuity remains the same regardless of the placement of the discontinuity. Equipped with the value of v 0 , we use (7.18), rewritten as (7.25)r to find the unknown starting timet 0 as The independence of v 0 fromr 0 along with (7.26) implies the initial start time is proportional to the initial radius of the discontinuity. Findingt 0 enables us to build the initial profile of the FRW-1 metric for any radial coordinater <r 0 by computing ξ =r/t 0 and using equations (7.15)-(7.17).
To compute the TOV metric, the A metric component is already determined beforehand by a constant (7.20). To find the other metric component, match it at the discontinuity forcing the constant B 0 to take the form With the TOV time scale B 0 , we can build the TOV metric for any radial coordinater >r 0 by using the equations (7.19)-(7.21).
Combining all this data together, we build the following functions v init (r), ρ init (r), A init (r), and B init (r) to use as the initial data at timet 0 (7.26), depending on the free parameterr 0 . Because the other functions are based on v in the FRW-1 space, we start by stating the fluid velocity where ξ =r/t 0 . Based on this function v init , the initial density profile is One last piece of information needed to run the simulation is the boundary conditions or the ghost cells. Recall, we need both boundary conditions to evolve according to the Einstein equations to maintain consistency at the edges of our simulated universe. For the FRW-1 side, at the gridpoint x 0 , the fluid velocity is a function of timet j where the variable ξ is defined as ξ = x 0 /t j . With the fluid velocity, the fluid density becomes (7.34) Since the metric components are staggered relative to the fluid variables, we need to compute the half gridpoint and use it to find the corresponding velocity We use this velocity in the following computation for the metric components, The boundary condition for the TOV side is easier to implement. Since the TOV metric is independent of time, we set the data values during the initial setup of the function profiles, and they remain the same. These static values work for the fluid variables and the metric component A, but the function B, specifically the time scale B t (7.13), changes and must be rematched during each time step, as discussed in the last section. Using the above criteria for the TOV border (7.10), let x i be the gridpoint position of this border. We rematch the time scale at timet j by the following formula where B(t j , x i ) is the simulated solution at the coordinate (t j , x i ). In this section, we look at results of the simulation for the FRW-1/TOV matched model.
We use the initial profiles (7.29)-(7.32) along with the boundary conditions (7.33)-(7.38) developed in the last section to run the simulation. Figure 7.2 shows the initial profiles with these parameters for the fluid variables (ρ init , v init ) and the metric A init , along with the mass. By selecting the initial discontinuity atr 0 = 5, equation (7.26) gives the initial With these initial profiles, we run the simulation for one unit of time (i.e.t end =t 0 + 1). Figure 7.3 depicts the evolution of the fluid variables (ρ, v), giving us a frame by frame view for the evolution of the fluid variables across this time frame, evenly distributed from the left frame att 0 to the right frame att end . After the initial timet 0 , two shock waves are formed, the stronger shock moving out toward the TOV side and the weaker shock moving in toward the FRW-1 side, creating a pocket of higher density expanding and interacting with both the FRW-1 and TOV metrics; therefore, the secondary wave to the strong shock for this solution of the Einstein equations is another shock wave, reflecting back in.
Next, we focus our attention to the resulting solution at the end time,t end . Figure  The convergence of this solution, tested by successive mesh refinements, is shown in Table 7.1. Since we are implementing a first order method, we expect a convergence rate of 0.5 for the discontinuous fluid variables, while we expect a convergence rate of 1 for the continuous metric. Looking at the results in the table, the convergence rate of the fluid variables start around 0.5 for the big mesh sizes, as expected, and improve to 1 as we mesh refine, which is a little surprising. In contrast, the convergence rate for the metric components start around 1 and on average continue to stay around 1 under mesh refinement, regardless of the high/low swing in the convergence rate for the metric B component. We believe this change in the rates is due to numerical error from integrating the metrics up across the universe.
Next we study the preservation of the FRW-1 and TOV metrics outside of the interaction region. Figure 7.5 shows the numerical solution after one unit of time for a mesh with less On the other side, the TOV border criteria is sufficient because the fluid velocity is almost constant until the smeared shock wave causes a significant change. Even though it looks   Table 7.2. Table 7.2 shows, from left to right, the number of gridpoints, the position of the left edge to the cone of sound, the position of the FRW-1 border, the error between the two, and the rate of convergence associated with the error. The corresponding results are displayed for the TOV side in Table 7.3, where we record the convergence between TOV border and the right edge of the cone of sound. Notice how the error between the cone of sound and both borders are decreasing, approaching zero. Moreover, the convergence rate is increasing, approaching 1, a linear rate which we would expect from a first order method.
These results gives us confidence in our method for determining the edges to the cone of sound and the FRW and TOV borders.
Equipped with a method for determining the position of the FRW-1 and TOV borders, these borders are used to show convergence of the metrics outside the interaction region.
These borders are indicators on where to stop computing the error between the simulated and model solutions. For example, we consider the FRW-1 side error computation for n = 1024. Using Table 7.2, the numerical diffusion on the FRW-1 side ends at r * = 4.8495, Figure 7.5. Smearing of the shock waves for solution with less gridpoints (n = 2 10 = 1024) so when computing the error, we only consider mesh points x i such that x 1 ≤ x i ≤ r * . As we mesh refine, the point r * increases, getting closer to the cone of sound as shown in Table   7.2, and the region where we are performing the error calculation expands, converging to the edge of the interaction region. Using this procedure, the results for the convergence for the FRW-1 and TOV side are listed in Tables 7.4 and 7.5, respectively. These tables reveal both sides are converging to their respective model solution at the expected rate of 1. We find it remarkable that the TOV metric is preserved, regardless of all the mix-up in the interaction region. In particular, the TOV metric components remain intact even though we are integrating through the interaction region to obtain the metric.
Since our free parameter is the position of the initial discontinuity, we have a one parameter family of solutions to the Einstein equations based on the initial radius, and we briefly explore changing this parameter. Figure 7.7 shows the initial profile for shock position atr 0 = 95, withr min = 93 andr max = 97. We make a couple of observations  All these profiles are just stretched out over a longer region of space, causing them to look more like straight lines. There is less density and greater mass in this region of the universe since we are farther from the center of it. Figure 7.8 shows the solution after one unit of Notice for this solution the weak shock has a positive speed as opposed to the former case, as seen in Figure 7.4. The difference in speed from the earlier simulationr 0 = 5 indicates we are dealing with a quantitatively different solution from before. We explore changing this parameter many times and determine two conclusions in each case. One is the resulting solution always has a region of higher density surrounded by two shock waves, a strong shock on the TOV side and a weak shock on the FRW-1 side, and the other is the shock waves have different speeds, resulting in quantitatively different solutions. Hence, we truly have a one parameter family of quantitatively different shock wave solutions to the Einstein equations.

FRW-2 and TOV Matched Model Setup
After exploring the FRW-1/TOV model simulation, we cover the details in the setup of the FRW-2/TOV model. We recall the FRW-2 metric, so using the coordinate transformationr = √ tr, t = Ψ 0 2 4t 2 +r 2 t , (7.39) the FRW metric in standard Schwarzschild coordinates (the FRW-2 form) is Remember, unlike FRW-1, the FRW-2 metric relies on the FRW time coordinate t, which is the following function of (t,r) We match this metric to the TOV metric (equations (7.19)-(7.22)), and we follow a procedure similar to the FRW-1/TOV matching. Again, we choose the initial shock positionr 0 , assuming an initial start timet 0 > 0, and we want to match the FRW-2 and TOV metrics continuously at the point (t 0 ,r 0 ). Like before, we match the metric A component as in equation (7.23) to obtain where v 0 represents the initial velocity at the interface. Unlike before, we do not have a relationship between v 0 andt 0 (7.18) to solve for the initial timet 0 . Instead, we possess a relationship between v 0 and the FRW time coordinate t 0 in equation (7.43) to find the FRW time coordinate at the discontinuous interface In order to find the initial timet 0 , we solve for the integrating factor constant Ψ 0 , requiring us to determine the integrating factor. In the FRW-1 matched model, the integrating factor constant was suppressed, which is equivalent to setting Ψ 0 = 1, causing the (coordinate) speed of light √ AB on the FRW-1 side to be one. We follow the same paradigm with the FRW-2 matched model by choosing Ψ 0 such that the speed of light is one on the FRW-2 side, which is equivalent to setting Ψ = 1. More specifically, we choose Ψ 0 such that at the discontinuity (7.47) Ψ = Ψ 0 t 0 4t 2 0 +r 2 0 = 1, implying the integrating factor constant must be Now we use (7.39) to solve for the initial time withr 0 , t 0 (7.46), and Ψ 0 (7.48) as at the coordinate (t 0 ,r 0 ) provides us with Since Ψ = 1, the time scale constant B 0 becomes exactly the same as the FRW-1 case (7.28). With B 0 , we can compute the TOV metric for any radial coordinater >r 0 by using the equations (7.19)-(7.21).
With the matching complete, we build the functions v init (r), ρ init (r), A init (r), and B init (r) to use as the initial data at timet 0 . These initial profiles are similar to their FRW-1 matching counterparts (7.29)-(7.32), except the velocity is a function of η instead of ξ, and we rely on the FRW time coordinate t(t,r) defined in (7.44). Again, due to the dependence on the function v, we state the fluid velocity first where η =r/t. The initial density profile is Based the function v init , the metric components are giving by Since the boundary conditions are derived from the model equations, they contain the same discrepancies between the FRW-1 and FRW-2 matchings as the initial profiles. For the FRW-2 side, at the gridpoint x 0 at timet j , we use (7.44) to obtain the FRW time t j (t j , x 0 ) for computing the fluid velocity and density. We find the fluid velocity and the fluid density Since the metric components are staggered relative to the fluid variables, we need to compute the half gridpoint as before (7.35) and use it to find the corresponding velocity for η = x 1 2 /t j . We use this velocity in the following computation for the metric components, Since the TOV metric is independent of the FRW metric matched with it, the TOV boundary condition remains the same, and we use the same matching as in the FRW-1 case (7.38).

FRW-2 and TOV Matched Simulation Results
Using the initial profiles and boundary conditions developed in the last section, we run the simulation of the FRW-2/TOV matched model, comparing it with the FRW-1/TOV   We continue on and test our same solution hypothesis numerically. We accomplish this task by using the FRW-1/TOV solution at a fine mesh refinement, n = 16, 384, as the model to compare against. We perform convergence calculations assuming the FRW-1 matched solution is the true solution, so we run the FRW-2/TOV matched model for 0.9579   Table 7.6, and the data shows the FRW-2/TOV model is converging to the FRW-1/TOV model at a linear rate, confirming our same solution at different times hypothesis numerically.
We find the discovery that the two FRW transformations produce the same result too uncanny to be just a random coincidence, and we search for a theoretical connection to justify this discovery. Since we are dealing with a solution to the Einstein equations of a spherically symmetric metric in standard Schwarzschild coordinates, we take another look at the general form for this metric (7.62) ds 2 = −B(t,r)dt 2 + 1 A(t,r) dr 2 +r 2 dΩ 2 .
Let us consider the coordinate freedoms that preserves this metric; it turns out the only freedom of this metric is an arbitrary nonlinear transformation of the time coordinatet. We cannot change the radial coordinater because the spheres of symmetry will no longer match, With this fact in mind, we look for the transformation between the two time coordinates t 1 andt 2 for the FRW-1/TOV and FRW-2/TOV models, respectively. Considering an arbitrary coordinate (t,r), the FRW-1 timet 1 is (7.63)t 1 = 1 +r 2 4t 2 t = 4t 2 +r 2 4t , and the FRW-2 timet 2 is where we set the integrating factor constant (7.65) Ψ 0 = 4t 2 +r 2 t , as before (7.48) to make Ψ = 1 and initially match the B metric component for the FRW-1/TOV and FRW-2/TOV models. Studying these equations, we notice a couple of things.
Looking at the specific coordinate (t 0 ,r 0 ) for the initial discontinuity, we combine the equations (7.66) and (7.67) to determine Consequently, running the FRW-2/TOV simulation is just an exercise in the time coordinate invariance of this solution of the Einstein equations in standard Schwarzschild coordinates.

Summary
We provide the specific details in building the initial profiles (shown in Figure 7.2) and the boundary data required to run the simulation using the locally inertial Godunov method.
We have the first glimpse of a shock wave in general relativity by providing snapshots of the solution in Figure 7.3. This solution results in two shock waves, a strong outgoing wave and a weak incoming wave, enveloping a region of higher density. Changing our one parameter In simulating the FRW-2/TOV model, we force the initial profiles for this model to closely resemble the FRW-1/TOV model, and after running the FRW-2/TOV model, we notice many similarities to its counterpart, providing evidence that these are the same solution at different times. We test and confirm this hypothesis numerically by proving convergence between the two solutions in Table 7.6. We proceed to find the theoretical jus- closer to satisfying the black hole criteria (8.1), we are content with passing the Buchdahl stability limit of 0.9, which loosely states, after reaching this limit, a star is unstable and is subject to it own gravitational collapse. That is, the theorem of Buchdahl is whenever a mass gets within 9/8ths of its Schwzarzschild radius, or in our language its black hole number gets within 0.9, there doesn't exist a static solution of the TOV equations with a finite pressure at the center capable of holding the star up [ST97b]. This fact gives us confidence that not only does black hole formation occur in the reverse FRW/TOV model, but our simulation gives us a first glimpse at it.
Plugging (8.3) into the Einstien field equations (2.1), along with the perfect and comoving fluid assumption, gives the following pair of constraint equations on the functions R(t), ρ(t), and p(t): There exists explicit solutions, ρ(t) and R(t), to these pair of ODEs for t > 0 [ST95].
We turn our attention to how this time reversal affects the FRW-1 metric. If we repeat the process in Chapter 3 to obtain the coordinate transformation, the time reversed

Reverse FRW and TOV Matched Simulation Results
To better compare and contrast the forward and reverse solutions, we use the same initial parameters as before (7.1). Figure 8.1 shows the initial profiles for the fluid variables, the metric components, and the mass function. Comparing these profiles to the corresponding ones for the forward time in Figure 7.2, we see they match except the fluid velocity along with the start timet 0 = −5.4554 has changed signs, as expected from the above analysis.
After setting up the initial profiles, we run the simulation for one unit of time, but since time is running in reverse,t 0 is getting closer to zero, evolving toward the big bang instead of away from it (i.e.t end =t 0 + 1 = −4.4554). Figure 8.2 gives a frame by frame view for the evolution of the fluid variables across this time frame, evenly distributed from the left frame att 0 to the right frame att end . From the initial discontinuity, two rarefaction waves are formed, a stronger wave moving in towards the FRW region and a weaker one moving out towards the TOV region, and between the two waves a pocket of lower density is formed. These results differ from the forward time case where there are two shock waves surrounding a region of higher density. As time progresses, the rarefaction waves and the pocket are both expanding. To the left of this pocket, there is a density spike growing and moving in towards the FRW region, meaning more matter is falling into the center of the universe as time progresses. One can notice that indeed the FRW region is accumulating more mass from the start time (Figure 8.1) to the end time (Figure 8.2). If more mass keeps coming into a finite radius, then µ will continue to grow and eventually the black hole criteria will be satisfied (8.1), leading us to believe, given enough time, a black hole will evolve out of this solution. We explore this idea in more detail in the next section.
We turn our attention to the end result at timet end , where Again, we use the successive mesh refinement technique to test for the convergence of this solution, recorded in   (Table 7.1), where the discontinuities of the shock waves should make these rates lower, but numerical convergence is obtained at an adequate rate.
We examine the preservation of the FRW and TOV metrics outside the region of interaction. Looking back at Figure 8.3, the ends of the rarefaction waves are outside the cone of sound. Before, there was numerical diffusion of shock waves bleeding outside the cone of sound, but we do not expect to see a similar phenomenon with a rarefaction wave. Since the speed at which matter travels is bound by the speed of sound, this fact is a cause for concern. Like the forward model, we need to determine where the wave ends, and the border criteria developed in the Chapter 7 is used here. On the FRW side, the border criteria is where the derivative of the fluid velocity changes sign, and instead of increasing then decreasing in the forward time model, the fluid velocity is decreasing until it hits the  Table 8.2, while the TOV side is recorded in Table  Number ρ 3. The convergence rates are significantly lower from before, but numerical convergence is still obtained, albeit a little slower. It is uncertain why both convergence rates are much slower, but these results give us confidence, as we mesh refine, that the outer edges to the rarefaction waves will converge to the edges of the cone of sound.
Using these border criteria, we test the numerical convergence of the FRW and TOV metrics outside the region of interaction. Recall, these borders act as a marker where we stop computing the error between the simulated and model metrics, which is discussed more thoroughly in Chapter 7. Using this procedure, the convergence results for the FRW and TOV region are recorded in Tables 8.4 and 8.5, respectively. Looking at both tables, both regions are converging to their respective models at the appropriate rate of 1 for a first order method. Again, we note how remarkable it is the metric components in the TOV region remain preserved after integration through the FRW metric and the region of interaction.

Black Hole Formation
In the previous section, as time progresses in the reversed time model more mass falls in towards the center of the universe. If this trend continues, then it is inevitable, given enough time, a black hole will form. In order to run the simulation long enough for a black hole to form, our initial parameters must be changed. Using the previous parameters (7.1), the region of interaction hits the boundaries before the black hole forms, and once this happens,  Figure 8.4, at the edge of the incoming rarefaction wave, the density has continued to increase while velocity has become more negative, so matter has been coming in towards the center at a faster rate. Another thing to notice is the metric B component has grown extremely large at the TOV edge of our simulation. This large metric B component causes a huge time dilation between both edges of our simulation, which can be seen more clearly in Figure 8.6 where we graph the (coordinate) speed of light √ AB in yellow and the black hole number µ in cyan. The black hole criteria is the most satisfied in a region of space around the pointr = 2.38 where µ = 0.9110. Now more extreme relativistic effects come into play because the significant time dilation, measured by the speed of light, between the two edges is around 11,267 to 1, which means light travels 11,267 times faster at one edge of the simulated region to the other. This is a problem because it causes our future time steps to be smaller, so the evolution of the system with greatest µ is slower than the development close to the TOV border. As a final note, there is a blip in the velocity profile as seen in Figure 8.5. Since there is much data to support the accuracy of our numerical simulation, we have no reason to believe this blip is due to numerical error. We interpret this blip in the velocity as an unknown phenomenon, and it is a future research project to determine the nature of it.
In order to get closer to the formation of the black hole, we could try to expand the region of simulated space again, but the huge time dilation causes the region of interaction to grow extremely fast, quickly reaching the new border while only a small amount of time passes near the region where the black hole appears to be forming (i.e. µ = 0.9110 at r = 2.38). Regardless of how far we put the maximum radius, we never get closer than µ = 0.9110 in the black hole criteria (8.1). Our alternative strategy is to realize we are only interested in the region of space where the black hole is forming, where µ is greatest. With this in mind, to continue the evolution of the black hole, we start chopping off the right side of our region of simulated space; even though the right most grid cell is no longer valid, the one next to it (on the left side) is still a valid solution. This process is like zooming into the region of the black hole formation. More precisely, the grid point x n+1 is used as the boundary data throughout the entire simulation. After the TOV edge has been reached, we use the gridpoint x n as the new boundary data for the next time step, discarding the data associated with the gridpoint x n+1 . During the subsequent time steps, we discard the current boundary data at the gridpoint x i and replace it with its predecessor x i−1 . Using this procedure and stopping when the radius of the right boundary matches the radius with the greatest value of the black hole number, we produce  formation. More interestingly, this radius for the greatest µ is moving farther out, from r = 2.38 tor = 4.5, meaning the region of black hole formation is expanding which is shown in Figure 8.8. In Figure 8.7, one notices the metric component B, along with the speed of light in Figure 8.8, is growing exponentially at this radius, causing the time dilation to be more extreme. Here, the time dilation from the center to this radius is 29,517 to 1, two and a half times greater than before.

Summary
We began the chapter by arguing for the validity of the existence of a reversed time solution, and then setup the initial profiles and boundary conditions needed to start the simulation. It turns out the only difference from the forward model is that positive time is Unlike the forward time model, the derivative of the metric is continuous in Figure 8.3, so this model is a strong solution to the Einstein equations. We mimic the convergence results obtained in the forward time model. In particular, in Table 8.1, we record numerical convergence of the entire solution taken as a whole, including the interaction and noninteraction regions, using successive mesh refinement. In Tables 8.2 and 8.3, we also record numerical convergence of the interaction region to the cone of sound, using border tracking of the edges to the simulated rarefaction waves. In Tables 8.4 and 8.5, we record the numerical convergence of the non-interaction regions, using the true solutions. Like the forward time model, these results affirm the region of interaction is precisely the region of the cone of sound (no formal mathematical proof of this is known).
As time progresses in this solution, more and more mass is falling in towards the center of the universe. This fact leads us to believe that if the solution were to progress in time beyond the simulation time limit, the black hole number (8.2) would continue to increase and eventually the black hole criteria (8.1) would be achieved exactly, resulting in black hole formation. Our method is to expand the region of simulation by running the simulation until the TOV region is completely gone, giving us a value of µ = 0.9110. We then continue to run the simulation by discarding the incorrect boundary data and zooming into the region of black hole formation, and by this refined method we obtain convergence all the way up to the black hole criteria value µ = 0.9218. These results taken together give us confidence in our solution as a whole and provide strong support for our claim that this solution has black hole formation contained within it. As a final note, with the successful results obtained here and Chapter 7, we have simulated general relativistic analogs of both the 1 and 2 family of shock and rarefaction waves, the elementary waves of conservation laws.
This fact gives us confidence in the ability and accuracy of our locally inertial Godunov method to simulate any solution to the Einstein equations for a perfect fluid in standard Schwarzschild coordinates.

Putting Units into the Simulation
We end our discussion of shock wave formation by putting back the units to give more To take this one step further, we also maket in units of mass related tot * in the units of length by (9.9)t * = G ct .
Therefore, by setting c = 1 and G = 1, we unify all the units (time, length, and mass) into one set of units. For our discussion, we choose the mass to be this unified set of units, where the units of length and time can be recovered with equations (9.8) and (9.9), respectively. To recover the density in units of M/L 3 , we compute (9.11) ρ * = 1 G 3 ρ.
To determine the value of Newton's gravitational constant, we perform dimensional analysis to recover the factors of c that have been suppressed by setting c = 1. Looking up the units, the Newton's true gravitational constant, denoted asĜ, has units of (9.12)Ĝ = L 3 M T 2 .

Programming Code for the Simulation
This appendix provides an overview of all the source code used to perform the simulations throughout this thesis. A copy of this source code is provided on the attached CD, which is approximately 8,000 lines contained in 22 files. This code is written in C++, with a Microsoft C++ 5.0 compiler, using the OpenGL graphics library on the Windows operating system. Our goal here is to give the reader a sense of the organization of the code to enable the reader to find specific areas of interest within the code.
To start, the main file (main.cpp) initializes the displayed window and the simulation itself, and contains the code specific to the Windows OS. This file also handles the interface between the operating system and the simulation. The rest of the code is segregated into different classes in accordance to an object-oriented programming paradigm. Each of these classes has two files associated with them, a header *.h file and an implementation file • The CRelEulerSolver class is the Riemann solver for the relative compressible Euler equations. This class is derived from the CRsolver class, and it contains the specific equations tied to the relative compressible Euler equations, scribed in Section 4.1.
• The CClasEulerSolver class is the Riemann solver for the classical compressible Euler equations [Smo83]. This class was used to test our Riemann solver on a more simple set of equations to find errors in the algorithms contained in the CRsovler class. This class has no affect on the simulations in this paper and is left over code from the testing stage of our simulations.
To handle the rendering of our simulations, we use a number of stand alone classes to serve specific graphical needs. We present a list of these classes along with a brief explanation of their utility: • The Color class is used to consolidate the 3 color components into one vector. This class allows for easy manipulation and use of colors in OpenGL.
• The Font class gives us a method to display all the text used in our simulations.
• The Interval class is a concise method of keeping track of intervals on the real line.