Collapsing Dust

The Tolman-Bondi Solution to Einstein's Equations is a time-dependent solution describing the collapse of a pressureless perfect fluid, commonly termed "dust". We will use it to illustrate how one might go about solving Einstein's Equations from an ansatz, as well as to illustrate the use of curvature invariants as clues in uncovering relations between metrics. We will restrict our attention to 4-dimensional solutions in this section.

We begin with the most general form [Kramer] of a spherically symmetric, time-dependent spacetime, in our usual coordinate chart {r, c, φ, t}:

ds2 = B(r, t)2 dr2 + R(r, t)2 dc2 / (1 - c2) + (1 - c2) R(r, t)2 dφ2 - A(r, t)2 dt2
The velocity vector of an observer comoving with the source fluid is purely timelike. Choosing a unit comoving velocity vector
u = {0, 0, 0, 1 / A(r, t)}
and setting the pressure equal to zero, we find that the only nonzero component of the stress-energy tensor is
Ttt = - ρ(r, t)
where ρ(r, t) is the as-yet unknown fluid density.

We examine the conservation equation first, because it is first order in derivatives of the stress tensor. It is obviously going to be much simpler than Einstein's Equations, which of course are second order in derivatives of the metric. The nonzero components of the conservation equation are

Tar; a = ρ(r, t) A'(r, t) / A(r, t)
and
Tat; a = - (R(r, t) ρ(r, t) B. (r, t) + 2 B(r, t) ρ(r, t) R. (r, t) +
B(r, t) R(r, t) ρ. (r, t)) / (B(r, t) R(r, t))
where here and below we denote derivatives with respect to r with primes and derivatives with respect to t using dots. We see from the first equation that A(r, t) does not depend on r, and so with a simple rescaling of the time coordinate we can set gt t equal to -1. This is equivalent to the inverse coordinate transformation
tp -> Integrate[ A[t], t]
followed by a redefinition of the remaining metric functions in terms of the new time coordinate. We can use Mathematica to find a solution to the other conservation equation:
DSolve[ cons[[4]] == 0, B[r, t], {r, t}]
{{B[r, t] -> C[1][r] / (R[r, t]2 rho[r, t])
where cons[[4]] is the time component of the conservation equation and "C[1][r]" denotes an unknown function of r. We will call this function B[r] and work now from the metric
ds2 = B(r)2 dr2 / (R(r, t)4 ρ(r, t)2) + R(r, t)2 dc2 / (1 - c2) + (1 - c2) R(r, t)2 dφ2 - dt2
We now compute the Einstein Tensor and look at its simplest components first. By setting the "r r" component equal to zero, we find an equation for B(r) as a function of ρ(r, t) and R(r, t) and its derivatives, and by setting the "r t" component equal to zero we have a differential equation involving ρ(r, t) and R(r, t). Solving these equations we obtain
B(r)2 = R(r, t)4 ρ(r, t)2 R'(r, t)2 / (1 + R. (r, t)2 + 2 R(r, t) R.. (r, t))
and
ρ(r, t) = ρ(r) / (R(r, t)2 R'(r, t))
The former equation is accompanied by an implicit constraint: since B(r) is not a function of t, the time derivative of the right hand side must vanish. This gives us a constraint on R(r, t):
2 R. (r, t) R.. (r, t) + R(r, t) R... (r, t) = 0
= (R. (r, t)2 + 2 R(r, t) R.. (r, t)). / 2

== E(r).

(the justification for this nomenclature will appear below) and so
B(r) = R(r, t)2 ρ(r, t) R'(r, t) / (1 + 2 E(r))1/2
We now have the metric solely as a function of R:
ds2 = R'(r, t)2 dr2 / (1 + R. (r, t)2 + 2 R(r, t) R.. (r, t)) +
R(r, t)2 dc2 / (1 - c2) + (1 - c2) R(r, t)2 dφ2 - dt2
subject, of course to the constraint on R(r, t). The only remaining equation to solve is the "t t" component of Einstein's Equations, which has the solution
ρ(r) = - R(r, t) (2 R.. (r, t) R'(r, t) + R(r, t) R'.. (r, t)) / (4 π)
This solution too comes with an implicit constraint, which is satisfied so long as the previous one is. This means that
ρ(r, t) = - (2 R.. (r, t) R'(r, t) + R(r, t) R'.. (r, t)) / (4 π R(r, t) R'(r, t))
= - (R(r, t)2 R.. (r, t))' / (4 π R(r, t)2 R'(r, t))
so that
4 π R(r, t)2 ρ(r, t) dR(r, t) = - d(R(r, t)2 R.. (r, t))
Therefore, if we treat R(r, t) as the physical radial coordinate,
M(r, t) = - R(r, t)2 R.. (r, t)
is the mass contained in a ball of radius r at time t. r then becomes a label coordinate, labeling an infinitesimally thick spherical shell of dust. If we write the sum of the kinetic and potential energies of one such shell, we have
E = ε(r,t) R. (r, t)2 / 2 - ε(r,t) M(r, t) / R(r, t)
= ε(r,t) (R. (r, t)2 / 2 + R(r, t) R.. (r, t))

= ε(r,t) E(r)

where ε(r,t) is the mass of the shell. Hence E(r) is the energy per unit mass of the shell at coordinate label r.

We now have only to solve the constraint on R(r, t) in order to have a solution. We observe that both terms of the constraint have time derivatives of order 2 or larger. Hence one solution is any function linear in t. We therefore try a slightly more general solution of the form

R(r, t) = (f1(r) + f2(r) (a t + b))d
and substitute into the constraint to obtain
a (d - 1) d (3 d - 2) f2(r) (f1(r) + (a t + b) f2(r))d - 1 = 0
If you use the Mathematica ReplaceAll and Rule functions to naively substitute this ansatz for R(r, t) into the constraint equation, you will be disappointed: ReplaceAll is absolutely literal, and as a result, none of the derivative factors will be changed. The following code can help here:
rules[rlis_List] := Block[ {args}, (
args = Table[ rlis[[1, 1, i]], {i, Length[ rlis[[1, 1]]]}];

Return[ Flatten[ {rlis, Table[ Thread[ D[ rlis, args[[i]]]], {i, Length[args]}],

Table[ Table[ Thread[ D[ D[ rlis, args[[i]]], args[[j]]]],
{j, Length[args]}], {i, Length[args]}]}]])];
FullSimplify[ conseq /. rules[ rules[ {R[r,t] -> (f[r] + f2[r] (a t + b)) ^ d}]]]
The argument "rlis" to rules is itself a list of Mathematica rules. The variable "args" is a list of those coordinates which are arguments to the left hand side of the first Mathematica rule in rlis. The function returns a list of Mathematica rules which include rlis and its first and second derivatives with respect to all of the coordinates in args. As illustrated, you can apply the function to its own output to increase the order of the "derivative rules" it generates.

Returning to the constraint, we see that it is solved for d = 1 or d = 2/3. For d = 1 (the solution linear in t), we find that ρ is zero. But for d = 2/3, we have a nontrivial solution. The mass is then

M(r, t) = 2 a2 f2(r)2 / 9
= M(r) == f(r) / 2,
introducing the Tolman-Bondi mass function f(r). It remains to enforce the initial condition
R(r, 0) = r
which gives us
f1(r) = r3/2 - b f2(r)
Putting this all together we have the Tolman-Bondi metric
ds2 = (2 (r f(r))1/2 - t f '(r))2 dr2 / (4 f(r) (r3/2 - 3 t f(r)1/2 / 2)2/3) +
(r3/2 - 3 t f(r)1/2 / 2)4/3 dc2 / (1 - c2) +

(1 - c2) (r3/2 - 3 t f(r)1/2 / 2)4/3 dφ2 - dt2

with density function
ρ(r, t) = f(r)1/2 f '(r) / (r3/2 - 3 t f(r)1/2 / 2) (8 π ((r f(r))1/2 - t f '(r) / 2))
(We are ignoring the solution t -> -t, which represents expanding dust.) E(r) is zero for this solution, which indicates that its shells are marginally bound. If we note that for this solution,
R(r, t) = (r3/2 - 3 t f(r)1/2 / 2)2/3
we can put the metric in the form
ds2 = R'(r, t)2 dr2 + R(r, t)2 dc2 / (1 - c2) + (1 - c2) R(r, t)2 dφ2 - dt2
with density function
ρ(r, t) = f '(r) / (8 π R(r, t)2 R'(r, t))


Geometry

It is probably most useful to express geometric quantities in terms of the physical radius R(r, t). The nonzero Christoffel Symbols are
Γ rr r R''(r, t) / R'(r,t)
Γ rr t R'. (r, t) / R'(r,t)
Γ rc c - R(r, t) / ((1 - c2) R'(r,t))
Γ rφ φ - (1 - c2) R(r, t) / R'(r, t)
Γ cr c R'(r, t) / R(r, t)
Γ cc c c / (1 - c2)
Γ cc t R. (r, t) / R(r, t)
Γ cφ φ c (1 - c2)
Γ φr φ R'(r, t) / R(r, t)
Γ φc φ -c / (1 - c2)
Γ φ φ t R. (r, t) / R(r, t)
Γ tr r R'(r, t) R'. (r, t)
Γ tc c R(r, t) R. (r, t) / (1 - c2)
Γ tφ φ (1 - c2) R(r, t) R. (r, t)
It may be useful to note that
R'(r, t) = ((r f(r))1/2 - t f '(r) / 2) / (f(r) R(r, t))1/2

R. (r, t) = - (f(r) / R(r, t))1/2

R'. (r, t) = (r1/2 f(r) - (r3/2 - t f(r)1/2) f '(r)) / (2 R(r, t)2 f(r)1/2), and

R''(r, t) = (-((r f(r))1/2 - t f '(r) / 2)2 +

R(r, t)3/2 (2 (f(r) / r1/2) + t (f '(r)2 - 2 f(r) f ''(r)) / f(r)1/2) / 4) /

(2 f(r) R(r, t)2)

As expected from the spherical symmetry of the solution, the r-t hypersurface is geodesic, as is the r-c-t hypersurface. The structure of the Christoffel Symbols leads us to expect that the coordinate chart breaks down at the zeroes of both R(r, t) and R'(r, t). Zeroes of R'(r, t) are unphysical: if the physical radius does not depend on the shell label, then multiple shells are at the same physical radius.

To investigate the zeroes of R(r, t), we will examine the equation for null geodesics

ds2 = 0
When applied to radial geodesics (dc = dφ = 0) it implies
dr / dt = 1 / R'(r, t)
For such geodesics, the change in physical radius with time is
dR(r, t) / dt = R'(r, t) dr / dt + R. (r, t)
= 1 + R. (r, t)

= 1 - (f(r) / R(r, t))1/2

Hence if f(r) approaches zero at the same rate as R(r, t) approaches zero, an asymptotic observer will not see the singularity ar R(r, t) = 0. Otherwise it will be visible as a naked singularity.

The nonzero components of the Riemann Tensor are

R r c r c R(r, t) R'(r, t) R. (r, t) R'. (r, t) / (1 - c2)
R r φ r φ (1 - c2) R(r, t) R'(r, t) R. (r, t) R'. (r, t)
R r t r t - R'(r, t) R'.. (r, t)
R c φ c φ R(r, t)2 R. (r, t)2
R c t c t - R(r, t) R.. (r, t) / (1 - c2)
R φ t φ t - (1 - c2) R(r, t) R.. (r, t)
and the nonzero components of the Ricci Tensor are
R r r R'(r, t) (2 R. (r, t) R'. (r, t) + R(r, t) R'.. (r, t)) / R(r, t)
R c c (R'(r, t) R. (r, t)2 + R(r, t) R'(r, t) R.. (r, t) + R(r, t) R. (r, t) R'. (r, t)) / ((1 - c2) R'(r, t))
R φ φ (1 - c2) (R'(r, t) R. (r, t)2 + R(r, t) R'(r, t) R.. (r, t) + R(r, t) R. (r, t) R'. (r, t)) / R'(r, t)
R t t - (2 R'(r, t) R.. (r, t) + R(r, t) R'.. (r, t)) / (R(r, t) R'(r, t))
Many of the invariants are less than meaningful, but those which are, are given below, along with the denominators of the others:
R 8 π ρ(r, t)
Ra b Ra b 64 π2 ρ(r, t)2
Ra b c d Ra b c d ~ 1 / (R(r, t)4 R'(r, t)2)
Ra b Rca Rb c 128 π3 ρ(r, t)3
Ra b c d Ra c Rb d 384 π3 ρ(r, t)3
Ra b c d Rea Rb c d e ~ 1 / (R(r, t)6 R'(r, t)3)
Ra b c d Re fa b Rc d e f ~ 1 / (R(r, t)6 R'(r, t)3)
Ra b c d Reafc Rb e d f ~ 1 / (R(r, t)4 R'(r, t)2)
Ra b c d Reafc Rb f d e ~ 1 / (R(r, t)6 R'(r, t)3)
Ra b; c Ra b; c ~ 1 / (R(r, t)6 R'(r, t)6)
Ra b; c Ra c; b ~ 1 / (R(r, t)6 R'(r, t)6)
Ra b; a Rcb; c ~ 1 / (R(r, t)6 R'(r, t)6)
Ra b c d; e Ra b c d; e ~ 1 / (R(r, t)6 R'(r, t)6)
Ra b c d; a Reb c d; e ~ 1 / (R(r, t)6 R'(r, t)6)
Euler class 32 R. (r, t) (2 R.. (r, t) R'. (r, t) + R. (r, t) R'.. (r, t)) / (R(r, t)2 R'(r, t))
εa b c i ... εe f gi ... Rb c eh Rf g a h ~ 1 / (R(r, t)3 R'(r, t)2)
It is obvious from the above that the singularity at R(r, t) = 0 is essential.


Schwarzschild Limit

Since f(r) is twice the mass contained within a ball whose surface is labeled by the coordinate r, there is a simple limit to this metric: when f(r) = 2M, a constant, all of the mass is contained within all of the "shells", and so it is concentrated at r = 0. We should have then the Schwarzschild Spacetime. In this limit, we have
ds2 = 8 M r dr2 / (8 M (r3/2 - 3 t (M / 2)1/2)2/3) +
(r3/2 - 3 t (M / 2)1/2)4/3 dc2 / (1 - c2) +

(1 - c2) (r3/2 - 3 t (M / 2)1/2)4/3 dφ2 - dt2

This is not obviously related to the Schwarzschild metric, and so we need to find a coordinate transformation which makes explicit the relationship between the two. To do this, we examine the Kretschmann Invariant for both metrics. For the Tolman-Bondi metric in the Schwarzschild limit,
Ra b c d Ra b c d = 48 M2 / (r3/2 - 3 t (M / 2)4
while for the Schwarzschild metric it is
Ra b c d Ra b c d = 48 M2 / rS6
Equating these and solving for rS(r, t), we are led to the coordinate transformation
coordxform[ gs, {rs, c, ph, ts}, {(r3/2 - 3 t (M / 2)1/2)2/3, c, ph, h(r, t)}, {r, c, ph, t}]
from the Schwarzschild metric to the Tolman-Bondi metric. If we set the latter's t t component equal to -1, Mathematica is unable to solve the resulting differential equation for h(r, t) But if we apply the rule
(2 r3/2 - 3 21/2 M1/2 t) -> u(r, t),
the t t component becomes
(4 M + (8 M - 8 M2 / (u(r, t) / 2)2/3 - 21/3 u(r, t)2/3) h. (r, t)2 / (21/3 u(r, t)2/3 - 4 M)
Setting this equal to -1 and using DSolve we obtain the solution
h(r, t) = Integrate[ u(r, t)2/3 / (25/3 M - u(r, t)2/3), t] + h(r)
Once again, Mathematica is unable to perform this integration after the substitution for u(r, t). But we can cheat a little: we do not need h(r, t), but only some function whose time derivative is the same in order to effect the coordinate transformation from Schwarzschild time to Tolman-Bondi time. After a little experimentation, we notice that the time derivatives of
f ( u(r, t) ) dt
and
( f ( u(r, t) ) u. (r, t) dt ) / u. (r, t)
are identical if u(r, t) is linear in t. And of course it is. The latter integral, after substitution of our u(r, t), is
r3/2 / (3 (M / 2)1/2) - t + 27/6 M1/2 (2 r3/2 - 3 (2 M)1/2 t)1/3 -
4 M Tanh-1((2 r3/2 - 3 (2 M)1/2 t)1/3 / (25/6 M1/2))
Retaining the terms in t, we have the desired coordinate transformation relating the Tolman-Bondi metric in the Schwarzschild limit to the Schwarzschild metric itself. Hence we see that the Tolman-Bondi metric in the Schwarzschild limit is simply the Schwarzschild metric in comoving coordinates.

The next section discusses several computations involving BPS membranes in D=10 and D=11 supergravity.

©2005, Kenneth R. Koehler. All Rights Reserved. This document may be freely reproduced provided that this copyright notice is included.

Please send comments or suggestions to the author.