Royal Society Publishing

Simulation of general relativistic shock wave interactions by a locally inertial Godunov method featuring dynamical time dilation

Zeke Vogler, Blake Temple


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 C0,1 (Lipschitz continuous) at points of shock wave interaction, and cannot be smoothed to C1,1 within the C1,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 C0,1 at shock waves (Glimm 1965). Weak solutions are proven to exist for finite time starting from an initial C0,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 C0,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=σρ, σ=c2/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 C0,1 in SSC makes it feasible to create viable initial data by matching exact solutions continuously. Matching metrics at the stronger level of C1,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 Embedded Image2.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, Embedded Image is the coupling constant, Embedded Image is Newton’s gravitational constant, and we use the convention that c=1 and Embedded Image whenever convenient (Weinberg 1972). In component form, Embedded Image2.2where w=(w0,…,w3) 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 ds2=gij dxi dxj, and summation is assumed on repeated up-down indices.

We restrict to spherically symmetric gravitational metrics in SSC coordinates Embedded Image2.3where (t,r) are temporal and radial coordinates, Embedded Image is the line element of the unit 2-sphere, and x≡(x0,…,x3)≡(t,r,θ,ϕ) is the space–time coordinate system. It is well known that a general spherically symmetric gravitational metric of form ds2=−B(t,r) dt2+A(t,r)−1 dr2+2D(t,r) dt dr+C(t,r) 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 dt2 and dr2 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 ∇iTi0=0 and ∇iTi1=0, respectively, (cf. (2.1)), and these can be written in the locally inertial form Embedded Image2.4and Embedded Image2.5where Embedded Image, are the Minkowski stresses, related to Tij by Embedded Image2.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 Embedded Image2.7Using Embedded Image in (E1) gives the formula for the mass function M (Weinberg 1972), Embedded Image2.8

When p=σρ, σ=const., the components of TM are given by Embedded Image2.9where Embedded Imageis the velocity, cf. Groah & Temple (2004). The main point is that TM is independent of the metric, and unlike T, the equations (2.6) and (2.7) close when T0jM 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, Embedded Image2.10where Embedded Image are the Minkowski energy and momentum densities, A=(A,B) are the metric components, Embedded Image2.11is the flux, g(A,u,x)=(g0(A,u,x),g1(A,u,x)), with Embedded Image2.12and Embedded Image2.13gives the source term of the balance law, and h(A,u,x)=(h0(A,u,x),h1(A,u,x)), with Embedded Image

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 Embedded Image, a maximum radius Embedded Image, the number of spatial gridpoints n and a start time t0. In our simulations, the number of spatial grid points n is chosen to be a power of two (i.e. n=2k for some k). From these parameters, the mesh width Δx is determined to be Embedded Image3.1and is fixed throughout the scheme. Let (xi,tj) represent a mesh point in an unstaggered grid defined on the domain Embedded Image3.2The spatial points are defined as Embedded Image3.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 tj, a new time step is computed by Embedded Image3.4where the minimum is taken over all the spatial gridpoints at time tj of the metric Aij=(Aij,Bij), to be defined shortly. Starting at t0, the temporal mesh points are defined by Embedded Image3.5We assume at our current time tj for j≥0 there exists a solution u(tj,x) and A(tj,x) for (tj,x)∈D. This solution is either provided as the starting solution at t0 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(tj,x), let uΔx be given by piecewise constant states uij at time t=t+j as follows: Embedded Image3.6For notational convenience, we denote xi+xi+1/2 and xixi−1/2 throughout this paper. Define the grid rectangle Rij so the mesh point (xi,tj) is in the bottom centre, Embedded Image3.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 Aij, in each cell Rij throughout the simulation. These constant values are established by setting Embedded Image3.8This approximation makes AΔx discontinuous along each line x=xi,i=1,…,n, at each time step t=tj.

Figure 1.

Staggering of the metric A and the solution u.

Because the Godunov step is a three-point method, we need boundary data u0,j, A0,j for a left ghost cell and un+1,j, An+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 Rij, and the discontinuity in uΔx at the bottom centre of Rij thus poses a classical RP for the compressible Euler equations in Minkowski space–time, Embedded Image3.9Let Embedded Image denote the solution of (3.9) within the Riemann cell Rij, and define the concatenation Embedded Image for (t,x)∈Rij as the RP step of the fractional step scheme. Given this, the Godunov method computes the average of uRPΔx across the intervals [xi,xi+] at the next time step tj+1. Since the metric A is different on both sides of xi, separate averages must be taken over the left and right half cells and combined to obtain the true average. We let Embedded Image and Embedded Image denote the averages on the left and right half cells, respectively, and Embedded Image and Embedded Image 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 Embedded Image3.10and do the same for the right half cell Embedded Image3.11where the 2 accounts for the half cell calculations. Taking the average of these results leads to Embedded Image3.12defining our Godunov step of the method.

To define the ODE step of our fractional step method, let Embedded Image denote the solution to Embedded Image3.13starting from initial data Embedded Image, where Embedded Image takes the form Embedded Image3.14and Embedded Image3.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 Embedded Image3.16Therefore, uΔx(t,x) is equal to uRPΔ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 Embedded Image3.17In terms of these equations, define the metric as Embedded Image3.18and Embedded Image3.19

Finally, in order to update the metric and conserved quantities, we use the RP averages Embedded Image to replace the RP solution uRPΔx(t,x) and perform numerical integration on the analytical equations (3.16)–(3.19). This process leads us to define Embedded Image3.20The mass is Embedded Image3.21with Embedded Image3.22and the metric becomes Embedded Image3.23and Embedded Image3.24where Embedded Image3.25with Embedded Image3.26Note that since the metric is staggered relative to the conserved quantities, we use the in-between values, like xk and uΔx(xk,tj+1) in the update step. Let Ai,j+1=(Ai,j+1,Bi,j+1) denote the constant value for AΔx on Ri,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 uRPΔ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 Embedded Image, p=σρ and we neglect the source terms on the right-hand side, equations (2.4) and (2.5) reduce to Embedded Image4.1where Embedded Image4.2and Embedded Image4.3The RP is the initial value problem when the initial data u0(x) consist of two constant states separated by a jump discontinuity at x=0, Embedded Image4.4System (4.1) with (4.4) define the RP step of the locally inertial method. Except for the factor Embedded Image, equations (4.1) agree with the special relativistic compressible Euler equations, Embedded Image4.5and we will account for the factor Embedded Image by time-dilation. When Embedded Image, 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=σρ, Embedded Image as long as uL and uR satisfy ρLR>0 and −c<vL,vR<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 (u0,u1) and the fluid variables (ρ,v) is given by Embedded Image4.6and Embedded Image4.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 Embedded Imageand the Riemann invariants r and s are Embedded Image4.8and Embedded Image4.9where Embedded Image4.10so that (ρ,v) are given by Embedded Image4.11

The following formulae are also useful: Embedded Image

The 1,2-rarefaction curves R1, R2 are the straight lines s=const., r=const. in the plane of Riemann invariants, respectively, and the 1,2-shock curves S1, S2 are given by the following parametrization with respect to β, Embedded Image: Embedded Image4.12and Embedded Image4.13where Embedded Image4.14and Embedded Image4.15with shock speeds Embedded Image4.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 uL and uR to Riemann invariants via (4.8) and (4.9), numerically computing UM in the rs-plane, then converting UM back to fluid variables using (4.11).

5. Time dilation between space–time cells

When Embedded Image, the factor can be accounted for in the solution of the RP for (4.1) in a grid cell by solving the RP for Embedded Image, and then rescaling the time in that grid cell, by a factor of Embedded Image. The factor Embedded Image 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.

Figure 2.

Effects of time dilation.

Theorem 5.1

Let Embedded Image5.1represent the maximum time in a Godunov cell before the CFL condition is violated. If the change in time is shortened from Embedded Image to Embedded Image in a Godunov cell containing the solution to the RP’s u(t,x), the average across that grid cell at time Embedded Image Embedded Image and at time Δt, Embedded Image are related by Embedded Image5.2where Embedded Image is the ratio between the two times and uC 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(t0,x) and AΔx(t0,x) for t0>0. Assume these approximate solutions exist up to some time tend>t0 and converge to a solution Embedded Image as Embedded Image along with a total variation bound at each time step tjEmbedded Image6.1where Embedded Image represents the total variation of the function uΔx(tj,x) on the interval Embedded Image. Assume the total variation is independent of the time step tj 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 Embedded Image7.1in co-moving coordinates are Embedded Image7.2When p=ρc2/3, the solution is Embedded Image7.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 Embedded Image and k=0. Then the FRW metric (7.1) under the coordinate transformation Embedded Image7.4goes over to the following metric in SSC Embedded Image7.5where the fluid velocity v is related to Embedded Image by Embedded Image7.6

A direct consequence is the following corollary.

Corollary 7.1

The fluid variables (ρ,v) corresponding to (7.5) satisfy Embedded Image7.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 Embedded Image7.8with Embedded Image7.9and Embedded Image7.10where Embedded Image7.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 Embedded Image 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 Embedded Image comes out of the simulation, cannot be imposed ahead of time, and as a result, the TOV metric must be re-matched by adjusting Embedded Image at each time step.

To construct the initial dataset, assume an initial time Embedded Image and radius Embedded Image 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 Embedded Image, thereby posing an initial discontinuity in the fluid variables. By (7.5) and (7.9), we have Embedded Image7.12

Let Embedded Image represent the fluid velocity on the FRW side of the discontinuity so (7.12) implies Embedded Image7.13Note that v0 is independent of the free parameter r0. Using (7.6), we find the unknown starting time Embedded Image as Embedded Image7.14The independence of v0 from Embedded Image along with (7.14) implies the initial start time is proportional to the initial radius of the discontinuity. Finding Embedded Image enables us to build the initial profile of the FRW metric for any radial coordinate Embedded Image by computing Embedded Image 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 Embedded Image7.15forcing the constant B0 to take the form Embedded Image7.16

Combining (7.12)–(7.16) in the case σ=1/3, and letting Embedded Image, we define the initial data Embedded Image, Embedded Image, Embedded Image and Embedded Image posed at time Embedded Image, depending on the free parameter Embedded Image, as follows: Embedded Image7.17and Embedded Image7.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 x0, use (7.7) to express the fluid velocity and density as a function of time Embedded Image Embedded Image7.19Since the metric components are staggered relative to the fluid variables, we need to compute the half gridpoint x1/2=x0x/2, and use it to find the corresponding velocity Embedded Image7.20for Embedded Image. We use this velocity to compute the metric components, Embedded Image7.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 Bt in (7.9), and must be rematched during each time step. Using the above criteria for the TOV ghost cell, let xn be the gridpoint position of this border. We rematch the time scale at time Embedded Image by the following formula Embedded Image7.22where Embedded Image is the simulated solution at the coordinate Embedded Image.

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 Embedded Image, equation (7.14) gives the initial start time of Embedded Image. Putting σ=1/3 into the initial density (7.17) gives Embedded Image, 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. Embedded Image). Figure 3a 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 Embedded Image to the right frame at Embedded Image. After the initial time Embedded Image, 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.

Figure 3.

Solution after one time step. (a) Evolution of the fluid variables. (b) Derivatives of the metric. (Online version in colour.)

Next, we focus attention on the resulting solution at the end time, Embedded Image. Figure 3b 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 3b 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 3a 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, Embedded Image and Embedded Image, 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×1011M, 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 C1,1, even though there are no delta function sources in the curvature tensor G.


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.


View Abstract