Note: this is a work-in-progress. Source code, when considered to be sufficiently stable, will be provided through the open-source Maxima project at maxima.sourceforge.net.

In 1919, Theodor Kaluza proposed an extension to general relativity: using an appropriately constructed fifth dimension, he was able to incorporate electromagnetism into Einstein's theory of gravity. Recently, I endeavored to replicate the most basic of Kaluza's results: the equation of motion for a particle in empty dive-dimensional space, as seen from a four-dimensional perspective.

Now that I am working with Maxima, the question arose: can the same result be reproduced using this computer algebra system?

The answer, surprisingly, is a yes. Well, sort of. Some minor modifications were needed to the itensor package, nothing Earth-shattering really, the changes are probably best described as bug fixes. (If you're a Maxima user, these changes are not yet in the CVS version. Soon...)

One of the key features of the itensor package is the ability to define components of a symbolic tensor in an algorithmic fashion. For instance, you could specify the covariant and contravariant components of the tensor G like this:

components(g([i,j],[]),e([i,j],[])*p([],[]));
                               done
ishow(g([i,j],[]))$
                              p e
                                 i j

components(g([],[i,j]),e([],[i,j])/p([],[]));
                               done
ishow(g([],[i,j]))$
                                i j
                               e
                               ----
                                p

While this formalism allows you to define separate equations for different index patterns, it does not make it possible to define components depending in the index value. "Fortunately," the components feature was, until recently, broken. Why do I think that that is good news? Because it taught people an alternate method for defining tensor components, using Maxima's function formalism:

k(l1,l2):=if l2=[] then e(l1,[])*p([],[]) else e([],l2)/p([],[])$
ishow(k([i,j],[]))$
                              p e
                                 i j

ishow(k([],[i,j]))$
                                i j
                               e
                               ----
                                p

So if this works, why could one not use a function definition that also takes into account the values of the tensor components? Why not indeed. Here is the result of my work, a definition for the five-dimensional Kaluza-Klein metric as a function of the four-dimensional metric, the electromagnetic vector potential, and the scalar component:

a(l1,l2,[l3]):=if member(5,l3) then 0 else funmake('a,append([l1,l2],l3))$
g4(l1,l2,[l3]):=if member(5,l3) then 0 else funmake('g4,append([l1,l2],l3))$
g5(l1,l2,[l3]):=
    if member(5,l3) then 0
    else if l1#[] then
    (
        if not (predval(l1[1]<=4) and predval(l1[2]<=4)) then
            funmake('g5,append([l1,l2],l3))
        else if l1[1]<=4 and l1[2]<=4 then
            apply('g4,append([l1,l2],l3))+
                      g55*difflist(a([l1[1]],[])*a([l1[2]],[]),l3)
        else if l1[1]<=4 then g55*apply('a,append([[l1[1]],[]],l3))
        else if l1[2]<=4 then g55*apply('a,append([[l1[2]],[]],l3))
        else if l3#[] then 0 else g55
    )
    else if l2#[] then
    (
        if not (predval(l2[1]<=4) and predval(l2[2]<=4)) then
            funmake('g5,append([l1,l2],l3))
        else if l2[1]<=4 and l2[2]<=4 then apply('g4,append([l1,l2],l3))
        else if l2[1]<=4 then -apply('a,append([[],[l2[1]]],l3))
        else if l2[2]<=4 then -apply('a,append([[],[l2[2]]],l3))
        else if l3#[] then sum(difflist(a([i],[])*a([],[i]),l3),i,1,4)
        else 1/g55+sum(a([i],[])*a([],[i]),i,1,4)
    )
    else funmake('g5,append([l1,l2],l3))$

These functions use a small set of helper functions. One is used to determine if a predicate can be evaluated (in order to avoid errors when it cannot); one constructs the argument list to a function; and one recursively computes the derivative of a tensor product without evaluation:

predval(prd):=block([retval,saved_prederror:prederror],
    prederror:false,
    retval:ev(prd,pred)=true or ev(prd,pred)=false,
    prederror:saved_prederror,
    retval
)$
difflist(exp,lst):=if length(lst)=0 then exp
                   else difflist(idiff(exp,lst[1]),rest(lst))$

One of the conditions on the Kaluza-Klein metric is that the four-dimensional metric components and the electromagnetic vector potential are not dependent on the fifth dimension; i.e., their derivative with respect to the fifth coordinate is zero. This can be expressed through the following definitions:

a(l1,l2,[l3]):=if member(5,l3) then 0 else funmake('a,append([l1,l2],l3))$
g4(l1,l2,[l3]):=if member(5,l3) then 0 else funmake('g4,append([l1,l2],l3))$

Equipped with these tools, we can proceed. The first step is to define the metric and the contraction rules for the two metric tensors:

(imetric:g5,dim:5,defcon(g4),defcon(g5),defcon(g4,g4,kdelta),defcon(g5,g5,kdelta))$

Next, we declare some indices to be 4-dimensional indices:

assume(k<=4,l<=4,m<=4)$

We may also want to do a DERIVABBREV:TRUE to display derivatives in a more compact format.

The equation of motion of a particle in empty five-dimensional space can be written up as:

ishow('diff(x([],[a]),t,2)+'ichr2([b,c],[a])*'diff(x([],[b]),t)*'diff(x([],[c]),t)=0)$
                                      b     c        a       a
                                    (x )  (x )  ichr2    + (x )    = 0
                                        t     t      b c       t t

The main "trick" is to break down the first component in this equation depending on whether the summation indices B and C are equal to five or are in the 1..4 range:

ishow(part(first(%),1))$
                                             b    c        a
                                           (x ) (x )  ichr2
                                               t    t      b c

ishow(subst(m,c,%)+subst(5,c,%))$
                                  b     m        a       5    b        a
                                (x )  (x )  ichr2    + (x ) (x )  ichr2
                                    t     t      b m       t    t      b 5

ishow(subst(l,b,%)+subst(5,b,%)+part(first(%th(3)),2)=last(%th(3)))$
       l     m        a       5     l        a          a     5     m       a           a      5   2
     (x )  (x )  ichr2    + (x )  (x )  ichr2    + ichr2    (x )  (x )  + (x )    + CHR2    ((x ) ) = 0
         t     t      l m       t     t      l 5        5 m     t     t       t t       5 5      t

The resulting equation breaks down into two separate equations depending on the value of A. We're only interested in the case where A is a four dimensional index:

ishow(subst(k,a,%))$
       l     m        k       5     l        k          k     5     m       k            k      5   2
     (x )  (x )  ichr2    + (x )  (x )  ichr2    + ichr2    (x )  (x )  + (x )    + ichr2    ((x ) ) = 0
         t     t      l m       t     t      l 5        5 m     t     t       t t        5 5      t

Before we evaluate the Christoffel symbols in the resulting equation, we need to do something to protect the one not containing any 5-indices, since we do not wish to have it expanded:

ishow(subst(chr2klm,'ichr2([l,m],[k]),%))$
       5     l        k       L     m                  k     5     m       k            k      5   2
     (x )  (x )  ichr2    + (x )  (x )  chr2klm + ichr2    (x )  (x )  + (x )    + ichr2    ((x ) ) = 0
         t     t      l 5       t     t                5 m     t     t       t t        5 5      t

Now we can evaluate the remaining Christoffel symbols:

%,ichr2$
ishow(rename(%))$
  l     m               k %1   5     %2                                 k
(x )  (x )  chr2klm + g5     (x )  (x  )  (g5         - g55 a     ) + (x )    = 0
    t     t                      t      t    5 %1,%2         %2,%1        t t

We now break this up into two parts depending on whether %1=5:

map(lambda([u],block(if freeof(%1,u) then u else u+subst(5,%1,u))),first(%th(2)))=last(%th(2))$
ishow(rename(%))$
assume(%1<=4)$
%th(2),g5$
ishow(%)$
  k %1   5     %2                                 l     m               k
g4     (x )  (x  )  (a      g55 - a      g55) + (x )  (x )  chr2klm + (x )    = 0
           t      t   %1,%2        %2,%1            t     t               t t

Now is the time to factor out G55 in order to be able to apply a substitution rule that will match the electromagnetic field tensor:

map(lambda([u],factorout(u,g55)),%)$
ishow(ratsubst(-f([%1,%2],[]),a([%1],[],%2)-a([%2],[],%1),%))$
    k %1   5     %2                 l     m               k
- g4     (x )  (x  )  f     g55 + (x )  (x )  chr2klm + (x )    = 0
             t      t  %1 %2          t     t               t t

We're almost done. What remains is contracting the equation, restoring the Christoffel symbol that we protected earlier, and moving one term to the right to recover the usual format:

contract(%)$
ishow(rename(%))$
    5     %1    k          l     m               k
- (x )  (x  )  f   g55 + (x )  (x )  chr2klm + (x )    = 0
      t      t  %1           t     t               t t

%-part(first(%),1)$
ishow(subst('ichr2([l,m],[k]),chr2klm,%))$
  l     m        k       k         5     %1    k
(x )  (x )  ichr2    + (x )    = (x )  (x  )  f   g55
    t     t      l m       t t       t      t  %1

which is the equation we sought. Except for a little bit of cheating, that is: The Christoffel-symbol in this apparently 4-dimensional equation is, in fact, the 5-dimensional Christoffel symbol. But is that really a surprise? The presence of an electromagnetic field does more than just produce an apparent force on a test particle. It also contains energy, which means it bends spacetime; which should manifest itself in our equation as a change in the Christoffel-symbol. Without further ado, here's the analysis:

ishow(ichr2([k,l],[m]))$
  m %4
g5     (g5       - g55 (a  a     + a     a ) - g4       + g5      )
          l %4,k         k  l,%4    k,%4  l      k l,%4     k %4,l
-------------------------------------------------------------------
                                 2

rename(%)$
forget(%1<=4)$
subst(5,%1,%th(2))$
%,g5,g4$
ishow(%)$
   m
  a  (g55 a    + g55 a   )
           l,k        k,l
- ------------------------
             2

assume(%1<=4)$
%th(6),g5$
ratsubst(-f([%1,k],[]),a([%1],[],k)-a([k],[],%1),%)$
ratsubst(-f([%1,l],[]),a([%1],[],l)-a([l],[],%1),%)$
ishow(factor(contract(expand(%))))$
  m %1                    m    m                 m        m %1              m %1             m
g4     g4       - g55 a  f  + a  g55 a    - g55 f  a  - g4     g4       + g4     g4       + a  g55 a
         l %1,k        k  l           l,k        k  l            k l,%1            k %1,l           k,l
-------------------------------------------------------------------------------------------------------
                                                   2

%th(6)+%$
%, nouns$
ratsubst(chr42([k,l],[m]),g4([],[m,%1])*(g4([l,%1],[],k)+g4([k,%1],[],l)-g4([k,l],[],%1))/2,%)$
ishow(%)$
          m        m             m
  g55 a  f  + g55 f  a  - 2 chr42
       k  l        k  l          k l
- ----------------------------------
                  2

contract(%)$
ishow(map(factor, combine(distrib(%))))$
                    m    m
           g55 (a  f  + f  a )
     m           k  l    k  l
chr42    - -------------------
     k l            2

I don't know for sure but I'm pretty certain that this is the correct value for the contribution of the electromagnetic field to the metric.