Solving Geodesics from Christoffel Symbols

Question

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):

Needs["OGRe`"]
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.

Answer 1

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.

Answer 2

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.

Answer 3

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.