Compute the equations for the velocities

The geodesic equation can be written either as a second-order differential equation for the spacetime position of the test particle as a function of proper time -- or as a first order differential equation for the spacetime velocity of the test particle as a function of proper time.

But remember -- right now we're looking for how light travels in this spacetime. A ray of light travels like a test particle of zero mass, and test particles of zero mass don't have a "proper time" defined
for them. (Because light goes at the speed of light, which is never zero, hence there is no "rest frame" for a massless particle where proper time can be measured.) Instead there is a proper "affine parameter" that measures the distance along a lightlike curve through spacetime.

The easiest way to solve the geodesic equations is to start by defining the spacetime velocity vector, or the geodesic tangent vector as mathematicians call it. The "dot" below refers to the derivative with respect to the proper affine parameter v, but it turns out we can solve for the geodesics without ever having to use the proper affine parameter explicitly. We do this through mathematical trickery explained below.

In[76]:=


  DefineTensor[U,"U",{{1},1}]

  PermWeight::def: Object U defined

In[77]:=


  U[1] = rdot[r]
  U[2] = thdot[r]
  U[3] = tdot[r]

Out[77]=

  rdot[r]

Out[78]=

  thdot[r]

Out[79]=

  tdot[r]

In[80]:=


  U[la_] = MakeSum[Metricg[la,lb] U[ub]]

Out[80]=

  g   rdot[r] + g   tdot[r] + g   thdot[r]
   1a            3a            2a

In[81]:=


  U[-1]

Out[81]=

  rdot[r]

In[82]:=


  U[-2]

Out[82]=

             2  2
  (1 - 4 G M)  r  thdot[r]

In[83]:=


  U[-3]

Out[83]=

  -tdot[r]


Now it's time for the mathematical trickery! We will take advantage of the symmetry of this problem to skip several steps. We wind up in the end with three equations to solve. The first equation we get from wanting the paths of photons. This means that the vector U[ua] squares to zero. We'll also absorb the factor of (1 - 4GM) into the single parameter c:

In[84]:=


  eqn1 = Simplify[MakeSum[U[ua] U[ub] Metricg[la,lb]]] == 0
/. (1 - 4 G M) -> c

Out[84]=

         2          2    2  2         2
  rdot[r]  - tdot[r]  + c  r  thdot[r]  == 0


The second and third equations give us our constants of motion for the paths we seek, from the metric "dot product" of each basis Killing vector with the tangent vector to our geodesic. The first Killing basis vector is in the t direction, and the constant of motion it determines is the energy of the photon.

In[85]:=

  DefineTensor[Kt,"Kt",{{1},1}]

  PermWeight::def: Object Kt defined

In[86]:=


  Kt[1] = 0
  Kt[2] = 0
  Kt[3] = 1

Out[86]=

  0

Out[87]=

  0

Out[88]=

  1

In[89]:=


  eqn2 = Simplify[MakeSum[Kt[ua] U[ub] Metricg[la,lb]]] == -E

Out[89]=

  -tdot[r] == -E


The third equation sets the second constant of motion, which is the angular momentum of the photon around the point mass M:

In[90]:=

  DefineTensor[Ktheta,"Ktheta",{{1},1}]

  PermWeight::def: Object Ktheta defined

In[91]:=


  Ktheta[1] = 0
  Ktheta[2] = 1
  Ktheta[3] = 0

Out[91]=

  0

Out[92]=

  1

Out[93]=

  0

In[94]:=


  eqn3 = Simplify[MakeSum[Ktheta[ua] U[ub] Metricg[la,lb]]] == J
/. (1 - 4 G M) -> c

Out[94]=

   2  2
  c  r  thdot[r] == J


Now we just solve the system of equations for rdot, thdot and tdot:

In[95]:=


  eqlist = {eqn1, eqn2, eqn3}

Out[95]=

          2          2    2  2         2
  {rdot[r]  - tdot[r]  + c  r  thdot[r]  == 0,

                     2  2
    -tdot[r] == -E, c  r  thdot[r] == J}

Up to

Compute and solve the equations for light propagation