## Abstract

We introduce the *locally inertial Godunov method with dynamical time dilation*, and use it to give a definitive numerical simulation of a point of shock wave interaction in general relativity starting from a new initial dataset. Prior work of Groah and Temple justifies meeting the Einstein constraint equations for the initial data only at the weak level of Lipshitz continuity in the metric. The forward time simulations, presented here, resolve the secondary 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 solution via collapse associated with an incoming rarefaction wave. A new feature is that space–time is approximated as locally flat in each grid cell so that Riemann problems and the Godunov method can be implemented. Clocks are then dynamically dilated to simulate effects of space–time curvature. Such points of shock wave interaction are more singular than points on single shock surfaces because the coordinate systems that make space–time locally flat on single shock surfaces (Gaussian normal coordinates), break down at points of shock wave interaction.

## 1. Introduction

We summarize the results in the thesis (Vogler 2010) in which Vogler introduces what we term the *locally inertial Godunov method with dynamic time dilation*, a fractional step method for simulating spherically symmetric shock wave solutions of the Einstein–Euler equations of general relativity (GR) in Standard Schwarzschild Coordinates (SSCs) (Groah & Temple 2004). The underlying issue is that the gravitational metric appears to be singular at shock waves in SSC coordinates—the coordinates in which the Einstein equations take the simplest form (Groah & Temple 2004). The simulations here give a definitive numerical demonstration that the locally inertial Godunov method is nevertheless a viable first-order numerical method for simulating shock waves in SSC. Numerical convergence of the method is demonstrated for a one parameter family of initial data obtained by matching a critically expanding Friedmann–Robertson–Walker (FRW) space–time Lipschitz continuously to the inside of a static Tolmann–Oppenheimer–Volkoff (TOV) solution, creating a point of shock wave interaction at the interface. (Shock surfaces generically intersect on a two-dimensional surface in space–time, but radial shocks with spherical symmetry intersect on a three-dimensional surface, and cross at a single point in the (*t*,*r*) plane.) The subsequent evolution provides a simple model for a GR explosion containing both an imploding and exploding shock wave. The most interesting new feature of the numerical method is that curved space–time is approximated by Minkowski flat space–time in each grid cell making approximation by Riemann problems (RPs) and the Godunov method viable, and this is compensated for by dilating the clocks in each grid cell to account for space–time curvature. In this paper, we record the elements of the numerical method at a level sufficient for replication, present a definitive case for convergence of the method, and establish the Lipschitz structure of the gravitational metric at points of shock wave interaction in the forward time simulation of a GR explosion. (See Vogler (2010) for black hole formation in the time-reversed problem.) Points of shock wave interaction are more singular than points on single shock surfaces because Gaussian normal coordinates, the coordinate systems that make space–time locally flat on single shock surfaces, break down at points of shock wave interaction (Israel 1966). In the follow-up paper, Reintjes & Temple (submitted) present Reintjes’ new proof that, in contrast to Israel’s theorem for single shock surfaces (Israel 1966), the gravitational metric is essentially only *C*^{0,1} (Lipschitz continuous) at points of shock wave interaction, and cannot be smoothed to *C*^{1,1} within the *C*^{1,1} coordinate atlas. As a consequence, points of shock wave interaction are a new kind of *regularity singularity* where space–time is not locally flat in the sense that unbounded second derivatives of the metric exist in every coordinate system (Reintjes 2011). (Such singularities are distinct from surface layer singularities by the absence of delta function sources in the curvature tensor (Israel 1966; Gaspar & Racz 2011; Reintjes & Temple submitted).)

The point of departure for this work is the existence theory (Groah & Temple 2004), in which Groah and Temple obtain existence of weak shock wave solutions of the Einstein equations in SSC for gravitational metrics that are only *C*^{0,1} at shock waves (Glimm 1965). Weak solutions are proven to exist for finite time starting from an initial *C*^{0,1} metric and initial density and velocity profiles of bounded total variation. The locally inertial formulation is amenable to approximation locally by flat Minkowski space–time, so exact RP solutions for the relativistic compressible Euler equations can be used in each grid cell, thereby making the analysis of a fractional step Glimm scheme feasible (Lax 1957; Glimm 1965; Smoller 1983). The conclusion is that Lipschitz continuous gravitational metrics make physical sense in SSC as weak solutions of the Einstein equations in the presence of arbitrary numbers of interacting shock waves of arbitrary strength. But this is a general theorem, and there is no known analysis that rigorously details the local structure of solutions at points where shock waves interact. In his thesis, Vogler (2010) developed these ideas into a viable numerical method, and used it to carefully simulate a point of GR shock wave interaction in forward time (and black hole formation in backward time, c.f. Neilsen & Choptuik (1999)). This is interesting from several standpoints.

First, the numerics and supporting analysis here establishes and clarifies the *C*^{0,1} structure of the gravitational metric at points of shock wave interaction, and thereby numerically verifies the starting assumptions in Reintes proof (Reintjes 2011). Also, we assume *p*=*σρ*, *σ*=*c*^{2}/3 (the equation of state for both pure radiation and the extreme relativistic limit of free particles, cf. Weinberg 1972), and in this case the relativistic Euler equations form a Nishida system (Nishida 1968; Smoller & Temple 1993). Vogler’s exact global Riemann solver takes advantage of the Nishida property of shock curves, namely, that shock curves are globally translationally invariant in the plane of Riemann invariants. This greatly simplifies the RP step of the locally inertial Godunov method (Groah & Temple 2004).

The initial dataset is also interesting in its own right. First, it agrees with exact solutions on either side of the FRW–TOV intersection, so finite speed of propagation implies that exact boundary conditions can be imposed and convergence can be tested on either side of a bounded region of interaction. In particular, the method entails starting with boundary data on the FRW side of the simulation, and integrating through the region of interaction to recover (up to transformation of the SSC time variable) the TOV metric on the other side. Thus, we have a GR framework tailored for a definitive test of convergence of numerical methods at shock waves. The initial data also provide a natural starting point for a rigorous proof of GR shock wave or black hole formation from smooth solutions. Keep in mind that the admissibility of metrics at the lower regularity *C*^{0,1} in SSC makes it feasible to create viable initial data by matching exact solutions continuously. Matching metrics at the stronger level of *C*^{1,1} would be problematic, c.f. Hawking & Ellis (1973). Finally, the simulation can be interpreted as resolving the secondary wave neglected in the Smoller Temple shock wave model (Smoller & Temple 1995)

The numerical demonstrations are backed up by theorem 6.1 below, which states that if a sequence of approximate solutions generated by the locally inertial Godunov method converge pointwise almost everywhere and the total variation of the fluid variables remains uniformly bounded in time, then the limit solution is an exact weak solution to the Einstein equations. By this general theorem, it suffices only to demonstrate numerical convergence to a limit, with bounded oscillations, in order to conclude the simulated solutions accurately represent exact (weak) solutions of the Einstein equations. (See Vogler 2010; Vogler & Temple 2011 for details of the proof.)

In §2, we introduce the Einstein–Euler equations, and the FRW and TOV exact solutions. In §3, we give the elements of the locally inertial Godunov method. In §4, we discuss the RP step of the method. In §5, we put in the time dilation. In §6, we state theorem 6.1 on convergence. In §7, we give the one parameter family of matched FRW–TOV initial data. In §8, we present the forward time simulations and give a definitive numerical demonstration of the convergence of the locally inertial Godunov method in the presence of shock waves. We conclude by putting in dimensions to give an indication of the physical regimes to which the simulations apply.

## 2. Preliminaries

The Einstein–Euler equations *G*=*κT* of GR are equations for the gravitational metric tensor *g*, coupled to the relativistic compressible Euler equations
2.1through the Bianchi identities Div *G*=0 (Weinberg 1972). Here, *G* is the Einstein curvature tensor, *T* is the stress energy tensor for a perfect fluid, is the coupling constant, is Newton’s gravitational constant, and we use the convention that *c*=1 and whenever convenient (Weinberg 1972). In component form,
2.2where **w**=(*w*^{0},…,*w*^{3}) is the unit 4-velocity, *ρ* the energy density, *p* the pressure, and we use the Einstein summation convention whereby indices *i*,*j*=0,…,3 are raised and lowered with the metric d*s*^{2}=*g*_{ij} d*x*^{i} d*x*^{j}, and summation is assumed on repeated up-down indices.

We restrict to spherically symmetric gravitational metrics in SSC coordinates
2.3where (*t*,*r*) are temporal and radial coordinates, is the line element of the unit 2-sphere, and *x*≡(*x*^{0},…,*x*^{3})≡(*t*,*r*,*θ*,*ϕ*) is the space–time coordinate system. It is well known that a general spherically symmetric gravitational metric of form d*s*^{2}=−*B*(*t*,*r*) d*t*^{2}+*A*(*t*,*r*)^{−1} d*r*^{2}+2*D*(*t*,*r*) d*t* d*r*+*C*(*t*,*r*) *dΩ*^{2} can, under generic conditions, be transformed over to SSC (Weinberg 1972).

Putting the SSC metric ansatz (2.3) into MAPLE, the Einstein equations *G*=*κT* reduce to the four PDEs (1.2)–(1.5) of Groah & Temple (2004) referred to here as (E1)–(E4), respectively. (Beware that in Groah & Temple (2004), *A* and *B* stand for the d*t*^{2} and d*r*^{2} coefficients, respectively, cf. (2.3).)

In the presence of shock waves, the stress–energy tensor *T* is discontinuous, and thus *A* and *B* in (E1)–(E4) are Lipschitz continuous at best, and (E4) is satisfied only in the weak sense. In Groah & Temple (2004), it is shown that when the metric is Lipschitz and the stress–energy tensor is bounded in sup-norm, system (E1)–(E4) is weakly equivalent to the system obtained by replacing (E2) and (E4) with ∇_{i}*T*^{i0}=0 and ∇_{i}*T*^{i1}=0, respectively, (cf. (2.1)), and these can be written in the *locally inertial* form
2.4and
2.5where , are the Minkowski stresses, related to *T*^{ij} by
2.6Here, {}_{,0} and {}_{,1} denote derivatives with respect to *t* and *r*, respectively. The remaining two equations (E1) and (E3) can then be rearranged as
2.7Using in (E1) gives the formula for the mass function *M* (Weinberg 1972),
2.8

When *p*=*σρ*, *σ*=const., the components of *T*_{M} are given by
2.9where
is the velocity, cf. Groah & Temple (2004). The main point is that *T*_{M} is independent of the metric, and unlike *T*, the equations (2.6) and (2.7) close when *T*^{0j}_{M} are taken as the conserved quantities.

Using *x* in place of *r*, the equations (2.4), (2.5) and (E1), (E3) take the form of a system of conservation laws with sources,
2.10where are the Minkowski energy and momentum densities, **A**=(*A*,*B*) are the metric components,
2.11is the flux, *g*(**A**,*u*,*x*)=(*g*^{0}(**A**,*u*,*x*),*g*^{1}(**A**,*u*,*x*)), with
2.12and
2.13gives the source term of the balance law, and *h*(**A**,*u*,*x*)=(*h*^{0}(**A**,*u*,*x*),*h*^{1}(**A**,*u*,*x*)), with

Our purpose is to introduce an effective first-order method, the locally inertial Godunov method, and use it to compute a family of weak solutions of system (2.10). We next describe the method, then introduce an initial dataset, and end by simulating solutions starting from the initial data.

## 3. Locally inertial Godunov method

In this section, we define the algorithm for the locally inertial Godunov method starting from initial data with two boundary conditions. In §7, we define the initial data by matching an FRW to a TOV space–time at fixed time in SSC coordinates, and the boundary values will be determined by formulas for the FRW and TOV solutions on the left and right boundaries, respectively.

To start, fix a minimum radius , a maximum radius , 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
3.1and is fixed throughout the scheme. Let (*x*_{i},*t*_{j}) represent a mesh point in an unstaggered grid defined on the domain
3.2The spatial points are defined as
3.3Unlike 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 CFL condition for every time step. So for every time *t*_{j}, a new time step is computed by
3.4where the minimum is taken over all the spatial gridpoints at time *t*_{j} of the metric **A**_{ij}=(*A*_{ij},*B*_{ij}), to be defined shortly. Starting at *t*_{0}, the temporal mesh points are defined by
3.5We assume at our current time *t*_{j} for *j*≥0 there exists a solution *u*(*t*_{j},*x*) and **A**(*t*_{j},*x*) for (*t*_{j},*x*)∈*D*. 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:
3.6For notational convenience, we denote *x*_{i+}≡*x*_{i+1/2} and *x*_{i−}≡*x*_{i−1/2} throughout this paper. Define the grid rectangle *R*_{ij} so the mesh point (*x*_{i−},*t*_{j}) is in the bottom centre,
3.7(cf. figure 1). Each grid rectangle is a Riemann cell, containing a solution to a distinct RP. We are limited to solving RPs within Riemann cells having a constant speed of light, so the metric source **A**=(*A*,*B*) must be approximated by a constant value, denoted **A**_{ij}, in each cell *R*_{ij} throughout the simulation. These constant values are established by setting
3.8This approximation makes **A**_{Δx} discontinuous along each line *x*=*x*_{i},*i*=1,…,*n*, at each time step *t*=*t*_{j}.

Because the Godunov step is a three-point method, we need boundary data *u*_{0,j}, **A**_{0,j} for a left ghost cell and *u*_{n+1,j}, **A**_{n+1,j} for a right ghost cell (cf. figure 1). In our simulation in §7, we use boundary data from the exact FRW and TOV solutions, this being valid until the interaction region reaches the boundary. Figure 1 displays the location of these ghost cells.

The metric discontinuities are staggered relative to the fluid variables as illustrated in figure 1. Thus, the metric constants **A**_{Δx} establish a locally inertial coordinate frame in each grid cell *R*_{ij}, and the discontinuity in *u*_{Δx} at the bottom centre of *R*_{ij} thus poses a classical RP for the compressible Euler equations in Minkowski space–time,
3.9Let denote the solution of (3.9) within the Riemann cell *R*_{ij}, and define the concatenation for (*t*,*x*)∈*R*_{ij} as the RP step of the fractional step scheme. Given this, the Godunov method computes the average of *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. We let and denote the averages on the left and right half cells, respectively, and and the zero speed states of the left and right, respectively, as shown in figure 1. To perform the Godunov step on the left half cell, we compute
3.10and do the same for the right half cell
3.11where the 2 accounts for the half cell calculations. Taking the average of these results leads to
3.12defining our Godunov step of the method.

To define the ODE step of our fractional step method, let denote the solution to 3.13starting from initial data , where takes the form 3.14and 3.15

We define the approximate solution *u*_{Δx}(*t*,*x*) and **A**_{Δx}(*t*,*x*) analytically to derive the piecewise formulae used to update the numerical scheme and to be used in the convergence proof of §6. The conserved quantities are defined by the formula
3.16Therefore, *u*_{Δx}(*t*,*x*) is equal to *u*^{RP}_{Δx}(*t*,*x*), the solution to the RPs, plus a correction term from the ODE step of the method. The metric is derived from the definition of the mass
3.17In terms of these equations, define the metric as
3.18and
3.19

Finally, in order to update the metric and conserved quantities, we use the RP averages to replace the RP solution *u*^{RP}_{Δx}(*t*,*x*) and perform numerical integration on the analytical equations (3.16)–(3.19). This process leads us to define
3.20The mass is
3.21with
3.22and the metric becomes
3.23and
3.24where
3.25with
3.26Note 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, the locally inertial Godunov method constructs the solution inductively with four major steps: a RP step, a Godunov averaging step (with time dilation), an ODE step for sources and a final step to update the metric. The RP step (3.9) generates the approximation *u*^{RP}_{Δx}(*t*,*x*). Formulae (3.10)–(3.12) denote the Godunov step. The ODE step is detailed in (3.13)–(3.16); and equations (3.20)–(3.26) update the metric.

## 4. The Riemann problem step

In this section, we discuss the RP step of the locally inertial method. When , *p*=*σρ* and we neglect the source terms on the right-hand side, equations (2.4) and (2.5) reduce to
4.1where
4.2and
4.3The RP is the initial value problem when the initial data *u*_{0}(*x*) consist of two constant states separated by a jump discontinuity at *x*=0,
4.4System (4.1) with (4.4) define the RP step of the locally inertial method. Except for the factor , equations (4.1) agree with the special relativistic compressible Euler equations,
4.5and we will account for the factor by time-dilation. When , *p*=*σρ*, the RP for the compressible Euler equations was given in closed form in Smoller & Temple (1993). We summarize the results in the following theorem:

### Theorem 4.1

*There exists a solution of the RP for system (4.5) with an equation of state p=σρ,* *as long as u*_{L} *and u*_{R} *satisfy ρ*_{L}*,ρ*_{R}*>0 and −c<v*_{L}*,v*_{R}*<c. 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 waves and admissible shock waves.*

We now record the exact formulae required in the construction of the RP solutions (cf. Smoller & Temple (1993)). (These formulae correct typos in Smoller & Temple (1993), specifically in (2.5.73), (2.5.74), (4.2.12) and (4.2.13), cf. Groah & Temple (2004).)

To start, the mapping between the conserved variables (*u*^{0},*u*^{1}) and the fluid variables (*ρ*,*v*) is given by
4.6and
4.7which is 1−1 and non-singular in *ρ*>0, |*v*|<*c*. The eigenvalues (wave speeds) *λ*_{1,2}=*λ*_{−,+} of *DF* in (4.3) are given by
and the Riemann invariants *r* and *s* are
4.8and
4.9where
4.10so that (*ρ*,*v*) are given by
4.11

The following formulae are also useful:

The 1,2-rarefaction curves *R*_{1}, *R*_{2} are the straight lines *s*=const., *r*=const. in the plane of Riemann invariants, respectively, and the 1,2-shock curves *S*_{1}, *S*_{2} are given by the following parametrization with respect to *β*, :
4.12and
4.13where
4.14and
4.15with shock speeds
4.16The formulae (4.12) and (4.13) display the so-called Nishida property of this system, that *i*-shock curves based at different points are rigid translations of one another, and 1-shock curves are reflections of 2-shock curves (cf. Smoller & Temple (1993)). Constructing the Godunov averages requires computing the zero speed middle state of the RP, which is greatly simplified in plane of Riemann invariants by the Nishida property. Vogler’s RP algorithm solves for the middle state *u*_{*}≡*u*(*t*,*x*)=(*ρ*_{*},*v*_{*}) by first converting *u*_{L} and *u*_{R} to Riemann invariants via (4.8) and (4.9), numerically computing *U*_{M} in the *rs*-plane, then converting *U*_{M} back to fluid variables using (4.11).

## 5. Time dilation between space–time cells

When , the factor can be accounted for in the solution of the RP for (4.1) in a grid cell by solving the RP for , and then rescaling the time in that grid cell, by a factor of . The factor has the effect of dilating the physical (geodesic) time in a grid cell relative to the coordinate time (cf. figure 2). Thus, our GR Godunov method has the nice property that it allows for the use of exact RP solutions in each grid cell, and the consequent errors due to neglecting the local curvature are accounted for by simply rescaling the local time relative to coordinate time. The following theorem, based on this, expresses that if the time in a Godunov cell is shortened, the resulting average is an affine combination of the original average and the centre state, based on the ratio of the original and new time change.

### Theorem 5.1

*Let
*5.1*represent the maximum time in a Godunov cell before the CFL condition is violated. If the change in time is shortened from* *to* *in a Godunov cell containing the solution to the RP’s u(t,x), the average across that grid cell at time* *and at time Δt,* *are related by
*5.2*where* *is the ratio between the two times and u*_{C} *is the intermediate state between adjacent RPs.*

## 6. Convergence to a weak solution

Our main convergence theorem regarding the locally inertial Godunov method with dynamic time dilation is the following.

### Theorem 6.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* *as* *along with a total variation bound at each time step t*_{j}6.1*where* *represents the total variation of the function u*_{Δx}*(t*_{j}*,x) on the interval* *. Assume the total variation is independent of the time step t*_{j} *and the mesh length Δx. Then the solution (u,**A**) is a weak solution to the Einstein equations (E1)–(E4).*

The proof, adaptation of the proof in Groah & Temple (2004), involves demonstrating that the terms in the residual of numerical approximations of the locally inertial Godunov method that do not converge first order in Δ*x*, come back to cancel other terms of the same order, producing only terms first order in Δ*x*. The result implies that only convergence and stability need be verified numerically in order to conclude the limit is a weak solution of the Einstein equations. The proof is omitted. See Vogler (2010) and Vogler & Temple (2011) for details.

## 7. One parameter family of shock wave initial data

We obtain a one parameter family of initial data by matching FRW metrics Lipschitz continuously to TOV metrics in SSC. The Einstein equations for the *k*=0 FRW metric
7.1in co-moving coordinates are
7.2When *p*=*ρc*^{2}/3, the solution is
7.3where we have assumed *t* is time since the big bang (Temple & Smoller 2009).

To match FRW to TOV in SSC, we use the following representation of FRW found in Temple & Smoller (2009). (We use bar to distinguish SSC coordinates from unbarred FRW coordinates.)

### Theorem 7.1

*Assume* *and k=0. Then the FRW metric (7.1) under the coordinate transformation
*7.4*goes over to the following metric in SSC
*7.5*where the fluid velocity v is related to* *by
*7.6

A direct consequence is the following corollary.

### Corollary 7.1

*The fluid variables* (*ρ*,*v*) *corresponding to (7.5) satisfy*
7.7

The general relativistic version of TOV metrics that model static singular isothermal spheres was described in Smoller & Temple (1993). Assuming *p*=*σρ*, these are given by
7.8with
7.9and
7.10where
7.11depends on *σ*. The velocity is zero because the TOV metric is a time-independent metric in SSC coordinates (cf. eqn (3.4) of Temple & Smoller (2009)). The arbitrary function is included to account for the time scale freedom in (7.8), a freedom required to match the simulated FRW–TOV solution at a later time. That is, our simulation involves integrating the metric starting from FRW boundary data on the left-hand side of the simulation, and the TOV metric comes out of the simulation to the right of a region of interaction. But the time-scaling function comes out of the simulation, cannot be imposed ahead of time, and as a result, the TOV metric must be re-matched by adjusting at each time step.

To construct the initial dataset, assume an initial time and radius to be determined later. We match the metric components (*A*,*B*) of the FRW and TOV metrics (7.1) and (7.8) Lipschitz continuously at , thereby posing an initial discontinuity in the fluid variables. By (7.5) and (7.9), we have
7.12

Let represent the fluid velocity on the FRW side of the discontinuity so (7.12) implies
7.13Note that *v*_{0} is independent of the free parameter *r*_{0}. Using (7.6), we find the unknown starting time as
7.14The independence of *v*_{0} from along with (7.14) implies the initial start time is proportional to the initial radius of the discontinuity. Finding enables us to build the initial profile of the FRW metric for any radial coordinate by computing and using equations (7.5)–(7.7).

On the TOV side, *A* is already determined as the constant (7.9). To find *B*, use (7.5) and (7.9) to get
7.15forcing the constant *B*_{0} to take the form
7.16

Combining (7.12)–(7.16) in the case *σ*=1/3, and letting , we define the initial data , , and posed at time , depending on the free parameter , as follows:
7.17and
7.18

Finally, we give the boundary data for the ghost cells using the exact FRW and TOV solutions at their respective boundaries. For the ghost cell on the FRW side at the gridpoint *x*_{0}, use (7.7) to express the fluid velocity and density as a function of time
7.19Since the metric components are staggered relative to the fluid variables, we need to compute the half gridpoint *x*_{1/2}=*x*_{0}+Δ*x*/2, and use it to find the corresponding velocity
7.20for . We use this velocity to compute the metric components,
7.21The boundary condition for the TOV is static, so values of the fluid variables and the metric component *A* are constant in time, but the function *B* changes by the time scale factor *B*^{t} in (7.9), and must be rematched during each time step. Using the above criteria for the TOV ghost cell, let *x*_{n} be the gridpoint position of this border. We rematch the time scale at time by the following formula
7.22where is the simulated solution at the coordinate .

In the final section to follow, we present our simulation of solutions of (2.10) starting from initial boundary data given in (7.17)–(7.18) and (7.19)–(7.22), using the locally inertial Godunov method developed above.

## 8. FRW–TOV simulation in forward time

We now present the results of the forward time simulation of solutions of the Einstein equations (2.10) starting from the matched FRW–TOV initial data given in (7.17)–(7.18), together with the boundary conditions (7.19)–(7.22), in the special case *σ*=1/3. By selecting the initial discontinuity at , equation (7.14) gives the initial start time of . Putting *σ*=1/3 into the initial density (7.17) gives , so the density jumps down by one-third from the FRW side to the TOV side of the initial discontinuity.

With these initial profiles, we run the simulation for one unit of time (i.e. ). Figure 3*a* depicts the evolution of the fluid variables (*ρ*,*v*), giving us a frame-by-frame view for the evolution of the fluid variables, from the left frame at to the right frame at . After the initial time , two shock waves form, the stronger shock moving out towards the TOV side and the weaker shock moving in towards the FRW side, creating an intermediate pocket of higher density expanding and interacting with the FRW and TOV metrics on either side. We conclude that the incoming wave, viewed as secondary to the strong outgoing shock wave, is another shock wave, reflected back in.

Next, we focus attention on the resulting solution at the end time, . Figure 3*b* highlights where the two shock positions are relative to the cone of sound and the cone of light. The cone of light is represented by the white region while the cone of sound, embedded in the cone of light, is represented by the grey region. Note that the edges of the cone of sound align with the shock waves on either side, confirming that the interaction region between the two metrics lies completely within the cone of sound. In particular, the edges of the cone of sound move at the local sound speed, so if one of the edges of the cone of sound were to get slightly ahead or behind the shock position, then that edge would get pushed back into the shock as the characteristics impinge on the shock, thereby explaining the alignment. Figure 3*b* also displays the spatial derivatives *A*′ and *B*′ in the metric components *A* and *B*, respectively. The derivatives (*A*′,*B*′), found using numerical differentiation, have discontinuities aligned with the discontinuities associated with the fluid variables at the edges of the cone of sound. Figure 3*a* shows the profiles for the metric (*A*,*B*) are no better than Lipschitz continuous at shocks, reinforcing the fact that we have a weak solution to the Einstein equations.

Convergence was tested by successive mesh refinements confirming a first-order convergence rate. Convergence to FRW and TOV was confirmed on each side of the interaction region. A range of the initial data parameters were explored resulting in two conclusions. Solutions always exhibit a region of higher density between two shock waves, a strong shock on the TOV side and a weak shock on the FRW side, and shock speeds change with parameters, resulting in quantitatively different solutions. We conclude that the initial dataset produces a one parameter family of distinct but quantitatively similar shock wave solutions of the Einstein equations, each representing a point of shock wave interaction emanating from an initial discontinuity in the fluids.

To get a physical sense, the range of values in our simulation on a solar scale is 0.09 *M*_{⊙}<*M*<1.5 *M*_{⊙}, and , or about 1.4 times solar mass across a distance of about 5.94 km for a time interval of about 4.9 ms. Setting one unit of mass equal to a galactic mass of about 1.8×10^{11} *M*_{⊙}, the simulation corresponds to a mass of about 1.4 times galactic mass across a distance of about 0.12 light-years in a time interval of about 10 days.

## 9. Conclusion

In summary, this is a successful demonstration of convergence of an effective numerical method, at a point of GR shock wave interaction, in a coordinate system where the gravitational metric appears singular at shocks, in a framework tailored for a definitive test of numerical convergence. A proof by Reintjes presented in a follow-up paper (Reintjes & Temple (submitted)) shows that points of shock wave interaction exhibiting the structure simulated here are a new kind of singularity in GR where the gravitational metric tensor cannot be smoothed to *C*^{1,1}, even though there are no delta function sources in the curvature tensor *G*.

## Acknowledgements

This work summarizes results credited to Vogler’s doctoral dissertation, Vogler (2010), which was supervised by Blake Temple. Both authors were partially supported by the second author’s NSF Grant, where the problem was first proposed.

- Received June 17, 2011.
- Accepted March 6, 2012.

- This journal is © 2012 The Royal Society