I am somewhat new to using Mathematica and I am facing difficulties with a specific problem related to the geodesics of Einstein's field equation in a vacuum.

The metric I am working with is derived from "Introduction to General Relativity" (Ronald Adler, Maurice Bazin, and Menahem Schiffer). You can download the book here: https://www.jp-petit.org/books/asb.pdf

Specifically on page 402 (Equation 12.29):

TNewCoordinates["Spherical", {x0, r, θ, φ}]
TShow@TNewMetric["CustomMetric", "Spherical", 
  DiagonalMatrix[{1, -Exp[g[x0] + f[r]], -r^2 Exp[      
   g[x0] + f[r]], -r^2 Sin[θ]^2 Exp[g[x0] + f[r]]}]]

I have attempted to calculate the geodesics from the Christoffel symbols, and here are the equations I derived:

$$ 0_{x}^{x} = \ddot{x0} + \frac{1}{2} e^{f[r]+g[x0]} (\dot{r}^2 + r^2 (\dot{θ}^2 + \dot{φ}^2 Sin[θ]^2)) ∂_{x0}g[x0] 0_{r}^{r} = \ddot{r} - r (\dot{θ}^2 + \dot{φ}^2 Sin[θ]^2) + \frac{1}{2} (\dot{r}^2 - r^2 (\dot{θ}^2 + \dot{φ}^2 Sin[θ]^2)) ∂_r f[r] + \dot{r} \dot{x0} ∂_{x0}g[x0] 0_{θ}^{θ} = \ddot{θ} - Cos[θ] \dot{φ}^2 Sin[θ] + \frac{\dot{θ} (\dot{r} (2 + r ∂_r f[r]) + r \dot{x0} ∂_{x0}g[x0])}{r} 0_{φ}^{φ} = \ddot{φ} + \frac{\dot{φ} (\dot{r} (2 + r ∂_r f[r]) + r (2 Cot[θ] \dot{θ} + \dot{x0} ∂_{x0}g[x0]))}{r} $$

I would be grateful if anyone could assist me in solving these equations or providing insights.

  • 1
    $\begingroup$ Welcome to Mathematica S.E. To start: 1) take the intro tour now, 2) when you see good questions and answers, vote them up by clicking the tray triangles, since the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) consider accepting the answer, if any, that solves your problem, by clicking the checkmark sign 4) give help too, by answering questions in your areas of expertise. $\endgroup$
    – bmf
    Sep 12, 2023 at 14:34
  • 1
    $\begingroup$ Having said that I would like to point out that you have referenced a book, but no links to the book. Please, either add a link or consider using MathJax to write the equations to which you are referring $\endgroup$
    – bmf
    Sep 12, 2023 at 14:35
  • $\begingroup$ Are you interested in a symbolic or numerical solution? I would recommend you first do some research into how geodesic equations are approached, without the use of computer algebra. I.e. many approaches use variational methods, symmetries, etc. Read up on these, decide how you want to approach the problem, then ask for help once you've got something concrete for users on this site to look at. $\endgroup$
    – user87932
    Sep 12, 2023 at 18:33
  • $\begingroup$ Thank you for your advice. I am primarily looking for a numerical solution using Mathematica. However, I'm not yet certain about the specific initial conditions. Would you have any recommendations on how to proceed with Mathematica in this case? Any guidance or example would be greatly appreciated. $\endgroup$
    – HMZ
    Sep 12, 2023 at 21:17
  • $\begingroup$ One thing to consider is that energy/angular momentum are conserved quantities. What I meant in my earlier comment about symmetries is that you derive expressions for these terms and use them to simplify your equations - replace a group of terms with constants E and L. For a numerical solution, use NDSolve and related functions like ParametricNDSolve, once you've decided on initial conditions to use. You may be able to find other useful ideas tried in the past by typing "solve geodesic equations" in a search engine. I.e. useful variable substitutions, etc. $\endgroup$
    – user87932
    Sep 12, 2023 at 22:09

3 Answers 3


Euler-Lagrange equation for geodesics produce the Christoffel symbols as coefficients of the products of velocities in the equation of motion, solved for the second order derivatives of the coordinates

$ds^2 = \sum_{i,k} \ g_{ik} \ dx^i \ dx^k \ \to \ L = \frac{1}{2} \ \sum_{i,k} \ g_{ik}\ x'[s]^i \ x'[s]^k $

where s is the arc length as invariant parameter.

The Euler-Lagrange equation are

$$\frac{d}{ds} \ \frac {d L}{d x'[s]^i } - \frac {d L}{d x[s]^i } = 0$$

The Christoffel symbols have two sources, the time derivative of the metric from the first term and the gradient of the second, both producing quadratic expressions in the velocity components.

In Mathematica code

  tangentVectors = Through[{vt, vr, \[Omega]\[Theta], \[Omega]\[Phi]}[\[Tau]]]

$$\left( \text{vt}(\tau ),\text{vr}(\tau ),\omega \theta (\tau ),\omega \phi (\tau )\right)$$

    lagrangian = 1/2 tangentVectors . 
           DiagonalMatrix[{1, -Exp[g[x0] + f[r]], 
                    -r^2 Exp[ g[x0] + f[r]], 
                    -r^2 Sin[\[Theta]]^2 Exp[g[x0] + f[r]]}] . tangentVectors

$$\frac{1}{2} \left(-r^2 \sin ^2(\theta ) \omega \phi (\tau )^2 e^{f(r)+g(\text{x0})}-r^2 \omega \theta (\tau )^2 e^{f(r)+g(\text{x0})}+\text{vr}(\tau )^2 \left(-e^{f(r)+g(\text{x0})}\right)+\text{vt}(\tau )^2\right)$$

   cotangentVectors =( D[lagrangian, #] & /@ 
         tangentVectors) /. {r -> r[\[Tau]], \[Theta] -> \[Theta][\[Tau]]}

$$\left\{\text{vt}(\tau ),\text{vr}(\tau ) \left(-e^{f(r(\tau ))+g(\text{x0})}\right),r(\tau )^2 \omega \theta (\tau ) \left(-e^{f(r(\tau ))+g(\text{x0})}\right),r(\tau )^2 \omega \phi (\tau ) \sin ^2(\theta (\tau )) \left(-e^{f(r(\tau ))+g(\text{x0})}\right)\right\}$$

The Langrangian equations now yield

     ({0 == #}) & /@ ((D[ cotangentVectors, \[Tau]] - (D[lagrangian, #] &) /@ {t[\[Tau]], r[\[Tau]], \[Theta][\[Tau]], \[Phi][\[Tau]]}) //. 
 {vt' :> t'',  vr' -> r'', v\[Theta]' -> \[Theta]'', v\[Phi] -> \[Phi]'', 
 vt :> t',  vr -> r', \[Omega]\[Theta] -> \[Theta]', \[Omega]\[Phi] -> \[Phi]',
f_[\[Tau]] :> f}) // TableForm

$$\left( \begin{array}{l} 0=t'' \\ 0=\left(r'\right)^2 f'(r) \left(-e^{f(r)+g(\text{x0})}\right)-r'' e^{f(r)+g(\text{x0})} \\ 0=-r^2 \theta ' r' f'(r) e^{f(r)+g(\text{x0})}+r^2 \theta '' \left(-e^{f(r)+g(\text{x0})}\right)-2 r \theta ' r' e^{f(r)+g(\text{x0})} \\ 0=-r^2 \sin ^2(\theta ) r' \phi ' f'(r) e^{f(r)+g(\text{x0})}-2 r^2 \theta ' \sin (\theta ) \cos (\theta ) \phi ' e^{f(r)+g(\text{x0})}-r^2 \sin ^2(\theta ) \phi '' e^{f(r)+g(\text{x0})}-2 r \sin ^2(\theta ) r' \phi ' e^{f(r)+g(\text{x0})} \\ \end{array} \right)$$

Further simplifications now have to be made eg for pure radial lightlike rays.


On my first post, I found the geodesics equations found from Christoffel symbols:

0 = \ddot{x0} + \frac{1}{2} e^{f[r]+g[x0]} (\dot{r}^2 + r^2 (\dot{θ}^2 + \dot{φ}^2 Sin[θ]^2)) ∂_{x0}g[x0]

0 = \ddot{r} - r (\dot{θ}^2 + \dot{φ}^2 Sin[θ]^2) + \frac{1}{2} (\dot{r}^2 - r^2 (\dot{θ}^2 + \dot{φ}^2 Sin[θ]^2)) ∂_r f[r] + \dot{r} \dot{x0} ∂_{x0}g[x0]

0 = \ddot{θ} - Cos[θ] \dot{φ}^2 Sin[θ] + \frac{\dot{θ} (\dot{r} (2 + r ∂_r f[r]) + r \dot{x0} ∂_{x0}g[x0])}{r}

0 = \ddot{φ} + \frac{\dot{φ} (\dot{r} (2 + r ∂_r f[r]) + r (2 Cot[θ] \dot{θ} + \dot{x0} ∂_{x0}g[x0]))}{r}

It seems to be different from yours. Right?

Would it be possible for you to help me modify the given Mathematica code to integrate these suggestions? In particular, I'm interested in:

Incorporating the proposed metric changes and simplifications. Factoring in the conservation of energy and angular momentum, if they can be identified.

Moreover, given the equations of motion, what initial conditions would you suggest, especially for illustrative or insightful scenarios? For instance, initial conditions for radial lightlike rays or other types of trajectories in this spacetime? Do you mean $\omega_{\theta}(\tau) = 0, \quad \omega_{\phi}(\tau) = 0, \quad \theta'(\tau) = 0, \quad \phi'(\tau) = 0$?

Thank you once again for your invaluable help.

  • 1
    $\begingroup$ This answer would better as a new question, or you should edit your question to add this information. $\endgroup$
    – creidhne
    Sep 13, 2023 at 19:49

In an earlier comment I made, I mentioned that the term Exp[g[t] + f[r]] was common to all the spatial terms in your metric, and that you might want to consider factoring the exponential and redefining your radial coordinate to absorb the Exp[f[r]] term. When I did so, the resulting form of the metric reminded me of the Friedmann-Lemaître-Robertson-Walker (FLRW) metric for a homogeneous, isotropic and uniformly-expanding/contracting universe. I checked the book you linked to in your post (Adler et. al.) to confirm my suspicions, and in fact, the title of the section you reference is "Robertson-Walker Metric". The authors are developing their material using minimal assumptions and work their way toward some of the more common forms you find in the literature. I notice that section 12.4 in Adler et. al. addresses the problem of computing geodesics for this metric, which you're asking about, using more common forms of the metric. I'll present a somewhat different approach here.

First let me transform your metric to one which uses "hyperspherical coordinates". Your metric was

DiagonalMatrix[{1, -Exp[g[t] + f[r]], -r^2 Exp[g[t] + f[r]], -r^2 Sin[th]^2 Exp[g[t] + f[r]]}

Let Exp[g[t]] = a[t]^2, and let Exp[f[r]] = h[r]^2. The line element of your metric after these substitutions is

$ds^2 = dt^2 - a[t]^2 (h[r] ^2 dr^2 + h[r]^2 r^2 (d\theta^2 + \sin[\theta]^2 d\phi^2))$

Now transform to a new radial coordinate, $\chi$, using $d\chi = h[r] dr$. You can get r by integrating both sides. So $\chi = \chi[r]$, and the inverse gives $r = r[\chi]$. So you can go back and forth if you later wish to return to your original coordinate system. Define $h[r[\chi]] r[\chi]) = S_k(\chi)$, where k = 1, 0, or -1. (Adler et. al. introduce k in their derivation, so you can consult the text for an explanation.) Your metric now takes the form

$ds^2 = dt^2 - a[t]^2 (d\chi^2 + S[\chi]_k^2 (d\theta^2 + Sin[\theta]^2 d\phi^2)$

This metric was independently developed by FLRW in the 1920's/1930's, and is now the basis of the current standard model of cosmology. It's reasonable to assume that someone has worked out procedures for computing geodesics. The more you can simplify your expressions, the easier it's likely to be to find a solution. Rather than reinvent the wheel, I searched the web, and found a set of lecture notes (http://www.ncra.tifr.res.in:8081/~tirth/Teaching/Cosmology-2021/Lectures/lecture-04-web.pdf) (Tirthankar Roy Choudhury, NCRA/TIFR, 2021), which illustrates a method of obtaining the solutions for the FLRW geodesics which relies on symmetries of the equations. To make this answer self contained, I'll put most of the steps here, though I'm pulling most of this from his slides. (He uses R(t) instead of a(t).)

Start with the standard geodesic equation, using $\lambda$ as an affine parameter

$\frac{d^2 x^i}{d\lambda^2} + \Gamma^i_{jk} \frac{d x^j}{d\lambda}\frac{dx^k}{d\lambda} = 0$

Rewrite this in terms of derivatives of the coordinates; define $u^i = \frac{dx^i}{d\lambda}$. Then

$\frac{du^i}{d\lambda} + \Gamma^i_{jk} u^j u^k = \frac{du^i}{dx^k}\frac{dx^k}{d\lambda} + \Gamma^i_{jk} u^j u^k = 0$

Factor out the $u^k$ term. The other factor is a covariant derivative, so we get

$u^i_{;k}u^k = 0 $.

Use the metric to lower the index 'i'. Since the covariant derivative of the metric is zero, we get $u_{i;k} u^k = 0$.

Now expand this to obtain a different form for the geodesic equation. The point of these index gymnastics is that at the end, he obtains an expression without Christoffel symbols. Expand the covariant derivative (in parenthesis), rewrite $u_j$ as $g_{jl} u^l$, then bring the $u^k$ term at the end inside the parenthesis. $(\frac{\partial u_i}{\partial x^k} - \Gamma^j_{ik} u_j)u^k = (\frac{\partial u_i}{\partial x^k} \frac{dx^k}{d\lambda} - \Gamma^j_{ik} (g_{jl}u^l)u^k = 0$

Use the $g_{jl}$ term to lower the index on the Christoffel symbol. The Christoffel symbol can then be expanded into a sum of partial derivatives of the metric. (Remember that the metric is also symmetric in its indices.)

$0=\frac{du_i}{d\lambda} - g_{jl}\Gamma^j_{ik}u^lu^k = \frac{du_i}{d\lambda} -\frac12 (\frac{\partial g_il}{\partial x^k} + \frac{\partial g_ik}{\partial x^l} -\frac{\partial g_ik}{\partial x^l}) u^l u^k$

Now we come to the point. The $u^l u^k$ terms are symmetric under interchange of indices, and the first and third terms in the expanded Christoffel symbols are antisymmetric under this interchange. When summed over repeated indices, these products drop cancel out, leaving us with the following alternate form of the geodesic equation

$\frac{du_i}{d\lambda} = \frac12 \frac{\partial g_{kl}}{\partial x^i} u^l u^k$

We begin to make use of symmetries (closely associated with conserved quantities). The space described by this metric is homogenous and isotropic. This means we can have the geodesic pass through the spatial origin at $\chi=r=0$ without loss of generality.

Note that none of the terms in the metric depend on $\phi$. This implies that there is a symmetry in $\phi$. Since all of the derivatives of the metric wrt $\phi$ vanish, the rhs vanishes, implying that $u_\phi$ is a constant along the geodesic. We now use the relation $u_{\phi} = g_{\phi\phi} u^{\phi} = -a(t)^2 S(\chi)^2 \sin(\theta)^2 u^{\phi}$. Using the definitions for $\chi$ and $S_k(\chi)$ a few paragraphs up, at r=0, $\chi=S_k(\chi) = 0$. Since the geodesic passes through the origin where r = 0, we have $u_{\phi} = 0$ here. It's a constant along the geodesic, so it's zero all along the geodesic, as is $u^{\phi}$. Therefore $\phi$ is a constant. By symmetry, this holds for geodesics passing through all points in the space.

Now consider $\frac{du_{\theta}}{d\lambda}$. The only term in the metric which depends on $\theta$ is the $g_{\phi\phi}$ term, but this is multiplied by $u^{\phi}$ terms, which are zero, so $u_{\theta}$ is constant along the geodesic. Using reasoning similar to that for the $\phi$ case, we conclude that $u_{\theta} = 0$ and that $\theta$ is also a constant along the geodesic. So the geodesics are radial, and we're free to choose $(\theta,\phi)$ anywhere on the sphere. Not surprising, since we started with spherical symmetry.

Since $u^{\phi} = u^{\theta} = 0$, these terms drop out for the radial/time equations regardless of what the metric coordinate dependencies are. For the radial equation, note that $g_{tt} = 1$ and $g_{\chi\chi} = a(t)^2$, neither of which depend on $\chi$. So $u_{\chi} = g_{\chi\chi} u^r = -a(t)^2 \frac{d\chi}{d\lambda} = K$, where K is a constant. For this case, the constant need not be zero. It's value will depend on the initial conditions chosen.

That leaves the equation for the coordinate time, t. Since the goal is to derive an equation for $\frac{dt}{d\lambda}$, Choudhury considered two cases, massive and massless particles. For a massive particle, it's common to use "proper time" ($\tau$) as the affine parameter. In units where the speed of light, c, is unity, proper time is just 's' in the metric. So we have

$d\tau^2 = g_{\mu\nu} dx^\mu dx^\nu$, or dividing through by $(d\tau)^2$, $1 = g_{\mu\nu} \frac{dx^\mu}{d\tau} \frac{dx^nu}{d\tau}$.

For a massless particle (i.e. a photon/graviton/etc.), the proper time is zero all along the length of the geodesic, so the affine parameter is used to characterize the trajectory instead. However, for this case, we have $ds^2=0$. We can use this to determine the final equation we're after for radial geodesics

$((\frac{dt}{d\lambda})^2 - a[t]^2 (\frac{d\chi}{d\lambda})^2 =(\frac{dt}{d\lambda})^2 - \frac{K^2}{a[t]^2}= T$, where T = 0 or 1 for the massless/massive particles cases.

We're left with two equations to solve, for $\chi$ and t as functions of affine parameter $\lambda$. These are a bit simpler than the ones you started with.

A side note - I'm not sure whether you used coordinate time (x0 or t) as the denominator for your expressions, but if you did, that's a problem, because the geodesics are generally defined in terms of either proper time or an affine parameter, not coordinate time. It's possible to derive an expression using coordinate time, but it contains extra terms.

At this point, I'm going to leave the rest to you. Since a(t) depends on coordinate time, you'll need to solve for that and eliminate your affine parameter, but that should be straightforward. You can use Mathematica, but these equations should be simple enough to do by hand.

Of course this approach works because your metric possesses a high degree of symmetry. In more realistic cases, these techniques won't work and you may be left with no choice but to brute-force the problem. You can look at the link posted in another comment to see an example for a different metric to get an idea of how to proceed in those cases. Still, it's often useful to look at simplified cases to get a feel for the problem before diving in.


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.