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

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?

• It isn't clear what your problem really is, but you can find a simple "code" computing Einstein equations here How to calculate scalar curvature, Ricci tensor and Christoffel symbols in Mathematica? and an example how to deal with the Schwarzschild type metric leading to a simple system of ODE's. In general even in simplest cases there is no automatic way to solve given reduced Einstein equations. Nonetheless DSolve appears to be helpful. Mar 1, 2022 at 1:16
• @Artes, My problem is in calculating A(r) and B(r). When I use the DSolve command to solve the results of computing the tensors, I get an error that says "Equation or list of equations expected instead of in the first argument". Problem is I am providing a list of equations, so I have no idea what the problem could be.
– aygx
Mar 1, 2022 at 1:52
• Part is also very useful. If you have a list, say efe={{},{},{},...} you can look up what efe[[stuff here]] is giving you. Also, with all due respect, the OP is lacking crucial details. No definition of the Ricci tensor, Ricci scalar, and the energy-momentum tensor. Not everyone has a physics degree on this site, and maybe you would find more help with more details :-)
– user49048
Mar 1, 2022 at 2:06
• @kcr, I apologize for the lack of clarity in this question. I have updated the question to include formulas for the tensors.
– aygx
Mar 1, 2022 at 3:11
• thanks a lot for that. as @Artes mentioned it is not possible to fully automate the process. recall the in pen-paper derivation there are some non-trivial improvised steps to compute $A(r)$ and $B(r)$.
– user49048
Mar 1, 2022 at 3:13

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]))}


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?

• Is there way to extend the code eqs = {RicciTensor[g, xx][[1, 1]] == 0, RicciTensor[g, xx][[2, 2]] == 0} to include the ricci scalar and corresponding metric element? For example eqs = {RicciTensor[g, xx][[1, 1]]-1/2*ricci scalar*metric element == some quanity, RicciTensor[g, xx][[2, 2]]-1/2*ricci scalar*metric element == some quanity}
– aygx
Mar 1, 2022 at 22:08
• You can do it successfully e.g. in case of the Reissner–Nordström spacetime, however in general you shouldn't expect that Mathematica will solve symbolically coupled systems of nonlinear ODE's. Mar 2, 2022 at 0:30
• If I input something and get the error "ReplaceAll::reps:term, is neither a list of replacement rules nor a valid dispatch table, \ and so cannot be used for replacing" or just get back what I typed in, does that indicate that mathematica cant solve or is it an indication that im doing something wrong?
– aygx
Mar 2, 2022 at 0:39
• In general it doesn't, however warning ReplaceAll::reps... may appear e.g. when you evaluate a function depending on another undefined function using Set (i.e. =) insted of SetDelayed (i.e. :=). On the other hand if DSolve returns back a correctly defined system it suggests that it cannot solve the given system, however this issue is quite subtle. Mar 2, 2022 at 0:55
• So for the error would I just replace the = with :=? Or is there another fix?
– aygx
Mar 2, 2022 at 1:11

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.

• Thank you for the answer, although I already have the code for the computation the tensors. Im looking for a code that can solve the resulting equations given from the tensors being computed.
– aygx
Mar 1, 2022 at 11:46