How do I get the Schwarzschild solution from the Einstein Equations?

Question

Consider the Schwartzchild-like metric $$ds^2=-A(r)dt^2+B(r)dr^2+r^2(d\theta^2+\sin^2\theta d\phi^2).$$ The Einstein field equations for this metric reduce to $$R_{\mu\nu}=0,$$ which is also known as the vaccum Einstein equations. Now consider the fact that $\mu,\nu=(0,1,2,3)$. This yields the fact that the EFEs are actually a set of 10 independent highly coupled non linear hyperbolic-eliptic second order partial differential equations. These equations are incredibly difficult to solve by hand, hence a need for technology.

If I know the Ricci tensor $R_{\mu\nu}$ e.g. computed with Mathematica, does there exist any Wolfram Language code that will allow me to solve for the unknown functions in the metric $g_{\mu\nu}$?

Can I use just DSolve only?

Answer 1

What is a solution to the Einstein equations?

Solutions to the Einstein equations describe structure of the spacetime which is encoded in gravitational potentials forming a tensor, while from mathematical point of view it is a four-dimensional Lorentzian manifold, or more appropriately an equivalence class of isometric Lorentzian manifolds with metric tensor satisfying the Einstein equations. A uniqe property of the Einstein equations among partial differential equations of mathematical physics is that the solution defines a manifold on which solution itself is satisfied.

Cauchy problem for the Einstein equations

Due to the general covariance principle the Einstein equations are hyperbolic in a specific sense quite different from other hyperbolic equations. We generally say that a partial differential system of equations is hyperbolic if initial value problem is well posed (i.e. for appropriately large class of initial data there exists a uniqe solution depending continuously on initial data) and its solution has the domain of dependence property i.e. at a time $t>0$ solution depends on compact subset of initial data, in other words this means that physical signals propagate with a finite speed. Cauchy problem is the problem of parametrizing the solution space by the space of initial data. Initial data for the Einstein equations is a three-dimensional Riemannian manifold with a metric tensor (Riemannian) and the second fundamental form. Appropriate compatibility conditions of embedding that Riemannian manifold in the Lorentzian one form the Einstein constraint equations which are of elliptic nature. Having solved the constraint equations (i.e. preparing the initial data) we have to solve a coupled system of nonlinear $2$-nd order wave equations. There is an analog with the Maxwell equations for which we have to prepare initial data satisfying the Gauss law (the constraint equations) and then we can solve the hyperbolic propagation equations (with partial time derivatives).

The Schwarzschild spacetime

Solving the Cauchy problem for the Einstein equations is not the only way of constracting a physical spacetime. We can get a uniqe solution imposing appropriate conditions, e.g. if we assume that we are looking for vacuum, static, spherically symmetric, asymptotically flat solution we can find a unique class of solutions parametrized by one parameter and this is the case here. Historically the Schwarzschild solution was the first solution of the Einstein equations (precisely the Minkowski spacetime had been discovered earlier is another solution found without solving equations) and may be regarded as the most important in Genral relativity.

How can we find it?

In order to calculate the Ricci tensor we evaluate functions in How to calculate scalar curvature, Ricci tensor and Christoffel symbols in Mathematica? from the snippet starting at InverseMetric and ending at RicciScalar. There is a brief summary related to solving the Einstein equations, nevertheless it appears that another post discussing in details how we can obtain the Schwarzschild solution would be welcome. For description of motion of massive particles see e.g. The time-like geodesics (orbits) in the Schwarzschild spacetime

We define coordinates and a metric tensor of $4$-dimensional static spherically symmetric Lorentzian spacetime assuming a specific ansatz for the metric:

xx = {t, r, θ, ϕ};

g  = { {-A[r],    0,         0,            0},
       { 0,     B[r],        0,            0},
       { 0,        0,      r^2,            0},
       { 0,        0,        0, r^2 Sin[θ]^2}};

Since there are two unknown functions we need choosing two independent equations from the tensor equation $R_{\mu \nu}=0$ (here on the RHS we have $4 \times 4$ matrix with $0$ entries). It appears that we can restrict to the following system of ordinary differential equations $ R_{0 0} = 0$ and $R_{1 1}=0$ (mind the specific convention related to enumeration of tensor indices):

eqs = { RicciTensor[g, xx][[1, 1]] == 0, RicciTensor[g, xx][[2, 2]] == 0}

the same output in copyable form:

{-((A'[r] B'[r])/(4 B[r]^2)) + (A'[r]/r - A'[r]^2/(4 A[r]) +
  (A'')[r]/2)/B[r] == 0,
  (A[r] (4 A[r] + r A'[r]) B'[r] + r B[r] (A'[r]^2 - 
   2 A[r] (A'')[r]))/(4 r A[r]^2 B[r]) == 0}

and the solution is simply:

sol = Flatten @
        Simplify[{A[r], B[r]} /. 
          DSolve[ eqs, { A[r], B[r]}, r, GeneratedParameters -> c]]

 {-(c[1]/r) + c[2], (r c[3])/(c[1] - r c[2])}

Now we should impose a physical assumption that the spacetime becomes asymptotically flat, that is $A(r) \to 1$ and $B(r) \to 1$ as $r \to \infty$, i.e.

parameters = 
  Flatten @ 
    Solve[
      Limit[ sol, r -> Infinity] == {1, 1}, 
              {c[1], c[2], c[3]}, Reals]

{c[2] -> 1, c[3] -> -1}

and finally

{A[r_], B[r_]} = sol /. parameters

{1 - c[1]/r, -(r/(-r + c[1]))}

we can see also that

1/A[r] == B[r] // Simplify

True

There is the only parameter c[1] which from the corespondence with the Newtonian gravity appears to be $c_1 = \frac{2G M}{c^2}$, where $G$ -gravitational constant, $c$- speed of light and $M$ is the total mass of a spherically symmetric distribution of matter, namely the solution is valid outside the matter distribution what says Birkhoff's theorem and there is an analog with the Newtonian gravity: How do I solve an integral equation related to the Newtonian gravity?

Answer 2

The original code related stuff that will be used in this post can be found here

Other relevant links from this site related to this discussion are the following: Artes' wonderful answer, another answer by myself, and another one.

pretty sure there are some more related material on this site...

See this recent answer of mine for a thorough explanation of the approach taken here. Here we are moving a bit faster to the end result.

If I understand the OP properly, we are interested in obtaining the functions $A(r)$ and $B(r)$ appearing in the metric.

We set the dimensionality of the spacetime

n = 4;

Now we write the coordinate system

coord = {t, r, θ, φ}

We insert the metric elements. Note that these are the ones with both indices down; $g_{\mu \nu}$. We do it in the following manner:

metric = {{-A[r], 0, 0, 0}, {0, B[r], 0, 0}, {0, 0, r^2, 0}, {0, 0, 0,
     r^2 Sin[θ]^2}};

Its inverse is of course $g^{\mu \nu}$

inversemetric = Simplify[Inverse[metric]];

Computation of the Christoffel symbols; indices are $\Gamma^{x}_{xx}$:

affine := affine = Simplify[Table[(1/2)*Sum[(inversemetric[[i, s]])*
       (D[metric[[s, j]], coord[[k]] ] +
         D[metric[[s, k]], coord[[j]] ] - 
         D[metric[[j, k]], coord[[s]] ]), {s, 1, n}],
    {i, 1, n}, {j, 1, n}, {k, 1, n}] ]

Computation of the Riemann tensor; indices are $R^{x}_{xxx}$

riemann := riemann = Simplify[Table[
    D[affine[[i, j, l]], coord[[k]] ] - 
     D[affine[[i, j, k]], coord[[l]] ] +
     Sum[affine[[s, j, l]] affine[[i, k, s]] - 
       affine[[s, j, k]] affine[[i, l, s]],
      {s, 1, n}],
    {i, 1, n}, {j, 1, n}, {k, 1, n}, {l, 1, n}] ]

Components of the Ricci tensor $R_{xx}$

ricci := ricci = 
  FullSimplify[
   Table[Sum[riemann[[i, j, i, l]], {i, 1, n}], {j, 1, n}, {l, 1, 
     n}] ]

Now, if you call ricci you get the output of the above which is a list of lists.

ricci

gives

{{-((Derivative[1][A][r] Derivative[1][B][r])/(4 B[r]^2)) + (
   Derivative[1][A][r]/r - Derivative[1][A][r]^2/(
    4 A[r]) + (A^\[Prime]\[Prime])[r]/2)/B[r], 0, 0, 0}, {0, (
  A[r] (4 A[r] + r Derivative[1][A][r]) Derivative[1][B][r] + 
   r B[r] (Derivative[1][A][r]^2 - 2 A[r] (A^\[Prime]\[Prime])[r]))/(
  4 r A[r]^2 B[r]), 0, 0}, {0, 0, 
  1/2 (2 - (2 + (r Derivative[1][A][r])/A[r])/B[r] + (
     r Derivative[1][B][r])/B[r]^2), 0}, {0, 0, 0, (
  Sin[\[Theta]]^2 (-r B[r] Derivative[1][A][r] + 
     A[r] (2 (-1 + B[r]) B[r] + r Derivative[1][B][r])))/(
  2 A[r] B[r]^2)}}

However, we can use the following

ricci[[1, 1]]

to get a single equation. This is the one related to the $\mu=\nu=0$ component of the tensor

-((Derivative[1][A][r] Derivative[1][B][r])/(4 B[r]^2)) + (
 Derivative[1][A][r]/r - Derivative[1][A][r]^2/(
  4 A[r]) + (A^\[Prime]\[Prime])[r]/2)/B[r]

This is where automatization ends to be honest. Then, one has to do some re-writing manually following the logic outlined here.