I am self-studying General Relativity.

Is there a method for obtaining the metric tensor exterior to a specified mass distribution numerically? In the simplest case of a spherical mass this should yield the Schwarzschild exterior geometry. I am primarily interested in such cases, without radiation fields (the simpler the better).

I realise that there is an ambiguity in the coordinate system chosen - presumably if such numerical methods exist they include a specification of the coordinate system.

My google-fu has been unable to find a simple answer to this question. The introductions to the topic of numerical GR that I found are dense, lacking in simple examples (if such things exist), and focus on gravitational waves. I realise that solving a system of non-linear coupled PDEs is not, in general, a simple task. I found a lot of literature talking about the so-called (3+1) method of foliating spacetime with space-like 3D hypersurfaces, but not much in the explicit sense of 'starting with these initial/boundary conditions and coordinate system, solve these equations using method x to obtain the metric tensor, here is some code for the simple example of a spherical mass'.

So basically: is it possible to start with a mass density function and obtain a numerical solution/approximation to the metric tensor exterior to this in some coordinate system, and if so, how?

If the answer to my question is 'no' or 'you are fundamentally misunderstanding something' I welcome being corrected.

  • $\begingroup$ Do you want a matrix at each gridpoint or do you have some other ansatz? $\endgroup$
    – Emil
    Commented Dec 26, 2018 at 16:04
  • 1
    $\begingroup$ I think you need to be more descriptive in your use of the term "numerically". Eistein's equations provide a second order non-linear PDE for the components of the metric for a given stress energy tensor. All you are really asking (imo) is how to solve a PDE using a numerical method. Am I correct in this? There are several techniques, relaxation, finite difference, finite element, just to name a few. $\endgroup$
    – user196418
    Commented Dec 26, 2018 at 16:09
  • 1
    $\begingroup$ @ggcg well, I don’t really know where to start, so a formulation of the EFE equations in a form that a solver could take as input for a simple case would suffice. But my confusion is deep: for the exterior solution the energy momentum tensor disappears, so how or what does one solve to find the metric there? Does one match at the boundary etc? I am primarily interested in understanding how to formulate the equations in a suitable manner, not in PDE solver techniques. $\endgroup$
    – Martin C.
    Commented Dec 26, 2018 at 18:10
  • $\begingroup$ @Emil that would certainly suffice, yes $\endgroup$
    – Martin C.
    Commented Dec 26, 2018 at 18:11
  • $\begingroup$ Related, if not dupe of, physics.stackexchange.com/q/71828/25301,https://…, physics.stackexchange.com/q/220562/25301, possibly some others $\endgroup$
    – Kyle Kanos
    Commented Dec 26, 2018 at 18:19

2 Answers 2


Such a simple question, but it opens up so many cans of worms. Here's my crack at a comprehensive answer simpler than what you'd find in a textbook. I'm sure a numerical relativist would have a much more inside scoop... this is coming from more of a mathematical perspective. First I'll try to clear up some general difficulties which I think were misunderstood in the question, then get to the example at the end. I was excited to see this question especially since, as you said, there's a big lack of simple examples here. Had a little too much fun, I hope the answer's not too long to be helpful!

Basically, yes, what you want to do can be done. But if you want to include a nontrivial matter distribution, it gets really hard really fast. To see why, we have to look at the overall problem we're trying to solve. It's not just "solve Einstein's field equations". It's "solve Einstein's field equations coupled to matter equations". We're trying to minimize the action $$ S = S_{gravity} + S_{matter}. $$ Varying the action with respect to the metric gives Einstein's equation, while varying with respect to whatever variables are in $S_{matter}$ gives the matter time evolution equations (this might be the curved-spacetime Klein-Gordon equation for a scalar field, for example).

Let's suppose, for definiteness, that the matter content is going to be modeled by some classical fields. Varying the total action gives a second-order coupled system of nonlinear PDE's, some coming from Einstein's equation, some coming from matter equations. Since the system is second order, we basically need to provide an initial configuration and initial time derivative for both the metric and each of the matter fields.

So even before trying to actually solve the problem, it's not as simple as: "Let me just specify $T_{\mu\nu}$ and then solve Einstein's equation for the metric". First you have to choose how to model the matter, and hope it's going to evolve the way you want it to. Then you have to set the initial data for the metric, and choose initial data for the matter fields that give you the desired initial $T_{\mu\nu}$. There are probably multiple ambiguous ways to do this, so watch out. Now you have to check that all your initial data match any constraints on the system (e.g. gauge fixing constraints). If you managed all that, you can finally hit go and pray that the time evolution does something like what you expected. Clearly this is a pretty big undertaking, and we haven't even considered all the hairy details yet. So much for simple examples.

And, by the way, you can't just ignore the matter dynamics and treat $T_{\mu\nu}$ as a fixed source, since, as Elio Fabri has pointed out in his answer, you can't really specify $T_{\mu\nu}$ without knowing the spacetime structure.

Fortunately, however, you asked about the Schwarzschild solution, and the Schwarzschild solution doesn't have any matter in it! Well, actually, the matter is all relegated to "elsewhere", and manifests as a boundary condition on the Weyl tensor. But anyway, locally, no matter. That sounds a lot easier. So let's talk about vacuum solutions, where $T_{\mu\nu}=0$ and $S=S_{gravity}$ only.

Vacuum Solutions

First of all, we need to pose the question.

Let's work in $D=N+1$ spacetime dimensions.

We are given:

  • An initial N-dimensional manifold $\Sigma$.
  • An initial spatial (positive-definite) metric $h_{ij}$ on $\Sigma$.
  • An initial time derivative of the metric, usually specified by the extrinsic curvature tensor $K_{ij}$ on $\Sigma$.
  • Some kind of spatial boundary conditions.

We want to obtain a spacetime $M$ with metric $g_{\mu\nu}$ such that:

  • The spacetime obeys the vacuum Einstein equation $G_{\mu\nu}=0$.
  • You can slice up $M$ into spacelike hypersurfaces, each topologically equivalent to $\Sigma$, and each obeying the desired boundary conditions, such that one of the slices is equivalent to the initial data.
  • The spatial hypersurfaces foliating $M$ can be labeled smoothly by a parameter $t$, increasing toward the future.

To convert this into a differential equation for the spatial metric evolving on $\Sigma$ we have to write out

  • $G_{\mu\nu}=0$,
  • The "Gauss-Codazzi" constraints relating $h_{ij}$ to $K_{ij}$ on $\Sigma$ (these come from differential geometry of embedded surfaces),
  • Four extra (gauge) constraints accounting for the coordinate ambiguity,

in terms of $h_{ij}$, $K_{ij}$, and a parameter $\alpha$ called the lapse function (which is basically just $g_{00}$, and tells you how much proper time elapses while traveling normal to $\Sigma_t$).

After all that, you finally have a set of differential equations for the spatial metric, along with some constraints on the initial data.

Yikes, still got complicated pretty fast, huh?

If you want to do all this without a giant effort, it's only going to be in the absolute simplest cases, and you basically have to already know the answer. We'll take a look at that, but first we need to consider some more of the challenges.

EDIT: Thought of a better way to explain previous section, but erasing that much would be bad form, so here's a quick addendum to it.


Assume that the solution spacetime to be obtained has the topology $\mathbb{R} \times \Sigma$ and can be covered by coordinates $(t,x^i)$ with metric $$ ds^2 = - \alpha^2 \, dt^2 + h_{ij} \, (dx^i + \beta^i dt) (dx^j + \beta^j dt). $$ We call $\alpha$ and $\beta^i$ the lapse function and shift vector, respectively.

If we can get a spacetime metric of this form which

  • obeys $G_{\mu\nu}=0$,
  • matches our initial data on the slice $t=0$,
  • matches any other constraints or boundary conditions,

we'd be done.

The equations of motion (see example at the the end for $\beta^i=0$ case) are derived by calculating extrinsic curvature in the assumed metric, applying $G_{\mu\nu}=0$ in the assumed metric, and utilizing the Gauss-Codazzi relations. The timelike components of $G_{\mu\nu}=0$ are not dynamical, and are called the Hamiltonian and Vector constraints.

Solutions to the equations of motion are not unique unless you specify four additional "gauge" constraints on the spacetime metric. The simplest possible choice would be $\alpha=1$ and $\beta^i=0$. But often more complicated constraints, which may involve the metric and extrinsic curvature in addition to the lapse and shift, make the problem easier.

Slicing: A big can of worms

Fundamentally, the question of how spacetime is sliced into spatial parts is what makes the problem so hard.

To see why, let's consider the simplest possible example: flat, empty, Minkowski spacetime.

Instead of starting from scratch, work backwards. We'll start with the full spacetime, slice it up into time-dependent spatial slices, choose one slice to act as the initial data, integrate from that initial data, and see if we recover the initial slicing, and thereby recover the Minkowski spacetime (I'm not going to actually perform the integration now, but that's the idea we're considering).

Consider the following four slicings:

Why is the difference between these so important? Consider each one:

(a) If you slice like this, the spatial metric is always flat, and the integration is easy. You can only do this because you already knew the answer.

(b) If you slice like this, the spatial metric could be something horrible, and the integration could be difficult, even though the spacetime is trivial. You might do this by accident if you start from scratch.

For many reasons, it's usually impossible to have a single numerical calculation cover the "whole" spacetime (in Minkowski it is possible, of course). This might be because you only know initial data in a certain region of space, because the spacetime itself has no global Cauchy surface, or because the spatial or spacetime manifold has no global coordinate system. When you're only dealing with part of the spacetime, the solution is only valid within a diamond-shaped (in causal diagrams) region called the "domain of dependence" of the surface: basically the subregion of spacetime where no outside influences have time to come in. For the next two we look at slicings of the domain of dependence of an initial surface.

(c) Slicing like this could work, and would make the integration easy, but leaves an ambiguity in the boundary conditions. You could do it, but you have to make sure the numerical solution includes the domain of dependence of the initial surface, and discard parts of the solution outside.

(d) This slicing has obvious boundary conditions and a clear physical interpretation in terms of the domain of dependence, but the spatial metric and differential equations are much worse.

When you start from scratch, you might not even know which of these four cases you're dealing with.

More Worms: Specifying the initial manifold

There's a whole bunch of constraints on specifying initial data for the metric. They come from the coordinate gauge constraints, from the "Gauss-Codazzi" relations I mentioned earlier, and from the timelike components of Einstein's equation.

Let's forget about all that. Much more fundamentally, let's talk about just picking a manifold and throwing a metric on it, with the hope of it corresponding to some particular thing. You can't just do it willy-nilly.

You might think, let me just take the manifold $\mathbb{R}^2$, give it coordinates $(x,y)$, and put the metric $ds^2 = dx^2 + \sin^2 \! x \, \, dy^2$ on it everywhere. Nope! That's the metric of a sphere. It curves back on itself. The manifold you specified doesn't actually make sense, because the geometry and the topology of the manifold are mismatched. There's some fancy theorems about this (sphere theorem, Gauss-Bonnet), but in general the point is that you can't just assume a metric on some coordinate patch makes sense. Locally you're always fine, if you work in a small enough neighborhood. But if you extend the patch far enough, you can run into trouble. The same type of thing can sometimes happen when, for example, you try to define "regularized" versions of Schwarzchild spacetimes that pass through $r=0$.

I don't think this is usually a practical issue in numerical calculations, but I've run into it before on occasion and wanted to mention it because it's pretty neat.

Schwarzschild Example

There's pretty much just one way to do the Schwarzschild example that makes any sense to attempt as a simple case. That is, to cover the entire "exterior" Schwarzschild patch, parameterized by $2M < r < \infty$ and $ -\infty < t < \infty$, sliced on lines of $t=const$.

This is a good choice because it exactly covers the domain of dependence of the initial surface $t=0$, it has obvious spatial boundary conditions, and the metric and curvature are simple.

On the other hand, it's not very physical. For the vacuum exterior to either a black hole or star, we'd have to adjust the initial domain of $r$ in a way that would make the numerical calculation less trivial.

Let's go for it. The spacetime metric we're supposed to end up with is $$ds^2 = -f(r) \, dt^2 + f(r)^{-1} \, dr^2 + r^2 \, d\Omega^2,$$ where $f(r) = 1-2M/r$.

Foliate on slices $t=const$. The spatial metric on each slice is $ds^2 = f(r)^{-1} dr^2 + r^2 d\Omega^2$. The intrinsic Ricci curvature on each slice is $R_{ij} = \frac{2M}{r^3} \, \textrm{diag}(-f(r)^{-1}, \frac{r^2}{2}, \frac{r^2}{2} \sin^2 \theta)$. The extrinsic curvature on each slice is $K_{ij}=0$ (this isn't trivial).

We're integrating to find the spacetime metric [1] $$ ds^2 = - \alpha^2 \, dt^2 + h_{ij} \, (dx^i + \beta^i dt) (dx^j + \beta^j dt). $$

Call the spatial coordinates $(x^1,x^2,x^3) = (r, \theta, \phi)$.

As our coordinate gauge conditions choose $\alpha = f(r)^{1/2}$ and $\beta^i=0$ (the $\beta^i$ are "shift vectors", which I left out earlier for simplicity).

For spatial boundary conditions, $h_{ij}$ is periodic in $\theta$ and $\phi$ in the usual spherical way, and $h_{ij}$ is equal to the initial value at $r = 2M$ and $r \to \infty$.

The initial data is $h_{ij} = \textrm{diag}(f(r)^{-1},r^2,r^2 \sin^2 \theta)$ and $K_{ij}=0$.

The initial data already automatically obey all other constraints, since we initally got them by slicing a spacetime.

From [1], the equations of motion with the above setup are $$ \partial_t h_{ij} = -2 \alpha K_{ij},$$ $$ \partial_t K_{ij} = \alpha [ R_{ij} - 2 K_{il} {K^{l}}_{j} + K K_{ij}] - D_i D_j \alpha, $$ where $K = {K^i}_i$, and $D_{i}$ and $R_{ij}$ are the spatial covariant derivative and intrinsic Ricci curvature on the slice.

At the first timestep, the right hand side of both equations of motion are zero (it's a bit tedious and you need to calculate the Christoffel symbols for the spatial metric on the slice, but the two nonzero terms in the bottom equation eventually cancel out). Since nothing changes, at each subsequent timestep the right hand sides remain zero, and the metric and extrinsic curvature remain constant. We can do this for all $t$.

So the answer we get by integrating is the above spacetime metric, with $h_{ij}(t) = h_{ij}(0)$.

That's the right answer! Exterior Schwarzschild is static! We're done! Pretty trivial answer, but that's how it goes. Clearly, a lot of the steps we took were only simple because we already knew the answer.

Thanks for the opportunity to write this out!

Here's the references I consulted:

[1] Luis Lehner. Numerical Relativity: A review. 2001. https://arxiv.org/abs/gr-qc/0106072.

[2] Wald. General Relativity. Chapter 10.

[3] Hawking and Ellis. Large Scale Structure of Spacetime. Beginning of chapter 7 (it gets ridiculous after the nice intro bit in 7.1).

  • 1
    $\begingroup$ This is a great answer, thanks! I'm going to have to take some time to digest it, but +1. $\endgroup$
    – Martin C.
    Commented Jan 8, 2019 at 9:10

$\def\bg{\mathbf g} \def\bu{\mathbf u} \def\bT{\mathbf T}$ Mine isn't properly an answer. In a sense it's another question, but I also believe I'm going to pinpoint a deep difficulty previous comments (and also the other links given by @KyleKanos) didn't address.

Actually I never attempted numerical GR. I taught GR for several years, but beyond the usual exact solutions I only showed problems with spherical symmetry (e.g. dust star collapse) where the difficulty I've in mind had a clear solution.

Here is what I mean. Suppose I want to study a general collapse, with no symmetry assumed a priori. I have some matter, and have to characterize it some way. Think of the simplest case, a perfect fluid. Then $\bT$ is known: $$\bT = (\rho + p)\,\bu \otimes \bu - p\,\bg \tag1$$ ($\bu$ is 4-velocity). Of course eq. (1) holds where matter is present, whereas $\bT=0$ in the empty space around. I may also assume that spacetime is asymptotically Minkowskian.

And now? Let me compare the present problem with a classical e.m. problem. There the unknowns may be (e.g.) $A^\mu$, four scalar functions of four coordinates. Frequently density and current of charge are taken as known (otherwise things get considerably complicated.). But above all spacetime structure is given and the choice of spacetime coordinates is at my discretion.

In case of GR it's not so, spacetime itself is the unknown. The same is true for $\rho$ and $p$ as functions of spacetime points. So I don't know how to choose coordinates. There are other difficulties I leave aside for brevity.

Thats'all. Maybe someone will tell me that the problem is well known (no doubt on that) and the solution is ... I look forward to learn about it.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.