For those who use GRTensorII, here is a simple and effective method for(adsbygoogle = window.adsbygoogle || []).push({});

The Lie algebra of Killing vector fields "exponentiates" to give the Lie group of self-isometries on our spacetime; if this acts transitively on a family of three-dimensional submanifolds, we can obtain a complete list of first integrals.

- determining the Lie algebra of Killing vector fields
- using the result to find first integrals of the geodesic equations via Noether's theorem on symmetries

This method requires no particular knowledge of GRTensorII.

I'll explain using a simple example, the Melvin electrovacuum, which models a cylindrically symmetric magnetostatic field, in which the EM field energy is the sole source of the gravitational field. Here is an (extensively commented) GRTensorII file you can use to define the static coframe in this spacetime:

Now fire up Maple and in your session type the following commands:Code (Text):

Ndim_ := 4:

x1_ := t:

x2_ := z:

x3_ := r:

x4_ := phi:

eta11_ := -1:

eta22_ := 1:

eta33_ := 1:

eta44_ := 1:

bd11_ := -(1+q^2*r^2):

bd22_ := (1+q^2*r^2):

bd33_ := (1+q^2*r^2):

bd44_ := r/(1+q^2*r^2):

Info_ := `Melvin non-null electrovacuum (1964); cyl chart; static obsvrs`:

# Chart -infty < t,z < infty, 0 < r < infty, -Pi < u < Pi

# (Discovered by Bonnor 1954 and Witten 1962 and then Melvin 1964)

#

# Four dimensional Lie algebra of Killing vector fields:

# @_t time translation

# @_z spatial translation in z direction

# z @_t + t @_z boost in z direction

# @_phi axial rotation

# So this is boost/rotation symmetric as well as cyl sym static.

# Geodesic equations:

# From Lagrangian we obtain quadratic invariant

# epsilon = (1+q^2*r^2)^2*(-t1^2+z1^2+r1^2)+phi1^2*r^2/(1+q^2*r^2)^2

# From Killing vectors we obtain four linear quadratic invariants:

# B = t1*(1+q^2*r^2)^2

# P = z1*(1+q^2*r^2)^2

# K = (z*t1+t*z1)*(1+q^2*r^2)^2

# J = phi1*r^2/(1+q^2*r^2)^2

# Solving { epsilon=0,B,P,J} for {t1,z1,r1,phi1} we obtain

# t1 = B/(1+q^2*r^2)^2

# z1 = P/(1+q^2*r^2)^2

# r1 = sqrt((B^2-P^2)*r^2-(1+q^2*r^2)^4*J^2)/(1+q^2*r^2)^2/r

# phi1 = J/r^2*(1+q^2*r^2)^2

# Helical null congruence

# J/sqrt(B^2-P^2) = r/(1+q^2*r^2)^2

# 1/(1+q*r^2) (B e_1 + P e_2) + sqrt(B^2-P^2)/r e_4

#

# Physical experience of static observers:

# Acceleration vector

# -2*q^2*r/(1+q^2*r^2)^2 e_3

# Expansion and vorticity tensor vanish

# Vector Potential

# -q/2*(1+q^2*r^2) @_phi

# Axial regular magnetic field concentrated near z = 0

# q/(1+q^2*r^2)^2 e_2

# Einstein tensor

# G^(ab) = 4*q^2/(1+q^2*r^2) diag(1,-1,1,1)

# Principal gravitational field invariants

# RR = 64*q^4*(3*r^4*q^4-6*r^2*q^2+5)/(1+r^2*q^2)^8

# Rstar = 0

# diRiem = 256*r^2*q^8*(133-126*r^2*q^2+45*r^4*q^4)/(1+r^2*q^2)^12

# CC = 192*q^4*(r*q-1)^2*(r*q+1)^2/(1+r^2*q^2)^8

# Cstar = 0

# diWeyl = 768*q^8*r^2*(15*r^4*q^4-42*r^2*q^2+31)/(1+r^2*q^2)^12

# Curvature everywhere finite and concentrated near r = 0

# Tidal tensor or electro-Riemann shows less axial tension than electro-Weyl tensor

# E_(ab) = EW_(ab) + 4*q^2/(1+q^2*r^2)^4 diag(1,0,0)

# This describes tidal accelerations of clouds of electrically neutral test particles.

# Weyl tensor is "pure electric"

# EW_(ab) = 2*q^2*(1-r^2*q^2)/(1+q^2*r^2)^4 diag(-2,1,1)

# Orthogonal hyperslices have three-dimensional Riemann tensor

# rv^2_(323) = rv^2_(424) = 2*q^2*(q^2*r^2-1)/(q*2*r^2+1)^4

# rv^3_(434) = 4*q^2*(q^2*r^2-2)/(q^2*r^2+1)^4

#

# For cylindrically symmetric functions, the Laplace-Beltrami equation

# becomes the usual wave equation

# h_(tt) = h_(rr) + (h_r)/r

# A simple exact solution is

# h = cos(a*t)*BesselJ(0,a*r)

# The corresponding vector field is

# sin(a*t)*BesselJ(0,a*r)/(1+q^2*r^2) e_1

# -a*cos(a*t)*BesselJ(1,a*r)/(1+q^2*r^2) e_3

#

# References:

# Bonnor, prsl A67 (1954): 225.

# Melvin, Phys. Lett. 8 (1964): 65.

Here is what we did:Code (Text):

qload(Melvin_nnevac_mg_stataxs_reg_cyl);

grcalcalter(x(up),basisv(dn),basisv(up),g(dn,dn),g(up,up),1):grdisplay(_);

grcalcalter(LieD[kv0,g(dn,dn)],1):grdisplay(_);

keqs := [seq(seq(grcomponent(LieD[kv0,g(dn,dn)], [grcomponent(x(up),[j]),grcomponent(x(up),[k])]),j=1..4),k=1..4)]:

# Enter

# T(t,z,r,phi)

# Z(t,z,r,phi)

# R(t,z,r,phi)

# Phi(t,z,r,phi)

op(%)[1];

pdsolve(%,{T,Z,R,Phi});

subs(%,[T,Z,R,Phi](t,z,r,phi));

sol := subs(_C1=c[1],_C2=c[2],_C3=c[3],_C4=c[4],%);

# Output

# [c[1]*z+c[4], c[1]*t+c[2], 0, c[3]]

subs(c[1]=1,c[2]=0,c[3]=0,c[4]=0,sol);

subs(c[1]=0,c[2]=1,c[3]=0,c[4]=0,sol);

subs(c[1]=0,c[2]=0,c[3]=1,c[4]=0,sol);

subs(c[1]=0,c[2]=0,c[3]=0,c[4]=1,sol);

# Output shows four dimensional Lie algebra of Killing vectors

# z @_t + t @_z

# @_t, @_z, @_phi

# Transitive on cylinders r=r0

#

multiply(matrix([[1,0,0,0]]), matrix( [[-(1+q^2*r^2)^2,0,0,0],[0,(1+q^2*r^2)^2,0,0],[0,0,(1+q^2*r^2)^2,0],[0,0,0,r^2/(1+q^2*r^2)^2]]), matrix( [[t1],[z1],[r1],[phi1]]))[1,1] = -EE; t1=solve(%,t1)

# Output shows t1 = EE/(1+q^2*r^2)^2

multiply(matrix([[0,1,0,0]]), matrix( [[-(1+q^2*r^2)^2,0,0,0],[0,(1+q^2*r^2)^2,0,0],[0,0,(1+q^2*r^2)^2,0],[0,0,0,r^2/(1+q^2*r^2)^2]]), matrix( [[t1],[z1],[r1],[phi1]]))[1,1] = A; z1=solve(%,z1);

# Output shows z1 = A/(1+q^2*r^2)^2

multiply(matrix([[0,0,0,1]]), matrix( [[-(1+q^2*r^2)^2,0,0,0],[0,(1+q^2*r^2)^2,0,0],[0,0,(1+q^2*r^2)^2,0],[0,0,0,r^2/(1+q^2*r^2)^2]]), matrix( [[t1],[z1],[r1],[phi1]]))[1,1] = L; phi1=solve(%,phi1); lprint(%);

# Output shows phi1 = L*(1+q^2*r^2)^2/r^2

multiply(matrix([[z,t,0,0]]), matrix( [[-(1+q^2*r^2)^2,0,0,0],[0,(1+q^2*r^2)^2,0,0],[0,0,(1+q^2*r^2)^2,0],[0,0,0,r^2/(1+q^2*r^2)^2]]), matrix( [[t1],[z1],[r1],[phi1]]))[1,1] = B; factor(%); -t*z1+z*t1=solve(%,-z*t1+t*z1): simplify(%,factor);

# Output shows z*t1-t*z1 = B/(1+q^2*r^2)^2

[t1 = EE/(1+q^2*r^2)^2, z1 = A/(1+q^2*r^2)^2, phi1 = L*(1+q^2*r^2)^2/r^2];

grcomponent(g(dn,dn),[t,t])*t1^2 + grcomponent(g(dn,dn),[z,z])*z1^2 + grcomponent(g(dn,dn),[r,r])*r1^2 + grcomponent(g(dn,dn),[phi,phi])*phi1^2;

subs(%%,%): simplify(%,factor): collect(%,EE): collect(%,A): collect(%,L);

lprint(%);

# Output shows

# r1^2 = (EE^2-A^2)/(1+q^2*r^2)^3 - L^2*(1+q^2*r^2)/r^2

The grand result is that we have the following list of first integrals for null geodesics:

- We computed the Killing equations for a generic vector field

[tex]

T \, \partial_t + Z \, \partial_z + R \, \partial_r + \Phi \, \partial_\phi

[/tex]

You can call the undetermined functions [itex]T,Z,R,\Phi[/itex] anything you like, of course.- We put these PDEs in a list.
- The key step is feeding this list of first order PDEs to Maple's powerful "casesplit" commend, which "triangularizes" them. (The form of the Killing equations guarantees than in principle, this can be done.)
- We use Maple's "pdesolve" command solves these system.
- We extract a list of four vector fields which generate the Lie algebra,
which confirms that the solution is static and cylindrically symmetric, and what it is more, it is a boost-rotation symmetric solution (see the review paper by Pravda and Pravdova in the "Useful Links" thread).

- the boost [itex]z \, \partial_z + t \, \partial_z[/itex]
- time translation [itex]\partial_t[/itex]; this turns out to be a vorticity free timelike vector field, so the constant t slices form a family of spatial hyperslices everywhere orthogonal to the world lines of the static observers,
- translation along axis of cylindrical symmetry [itex]\partial_z[/itex], a spacelike Killing vector field,
- rotation about the axis of cylindrical symmetry [itex]\partial_\phi[/itex], a spacelike cyclic Killing vector field
- Courtesy of Noether's theorem on variational symmetries of a system of PDEs determined by a Lagrangian, we immediately obtain several first integrals (more than we need, actually). You can use any notation you like for the first derivatives of a geodesic wrt an arc length parameter you like, but I find t1, z1, r1, phi1 particularly convenient for working with Lagrangians (when at various moments we need to consider derivatives as independent variables, or to work purely algebraically, as here). I used EE for energy in the Maple session so that there is no possibility of confusion with the defined constant [itex]\exp(1)[/itex]. In particular, we obtain first integrals for [itex]\dot{t}, \; \dot{z}, \; \dot{\phi}[/itex], but we still need a first integral for [itex]\dot{r}[/itex].
- We plug these into the line element [itex] 0 = ds^2[/itex] to obtain the first integral for [itex]\dot{r}[/itex], in the case of null geodesics.

[tex]

\begin{array}{rcl}

\dot{t} & = & \frac{E}{(1+q^2 \, r^2)^2} \\

\dot{z} & = & \frac{A}{(1+q^2 \, r^2)^2} \\

\dot{r} & = & \pm \sqrt{\frac{E^2-A^2}{(1+q^2 \, r^2)^3} - \frac{L^2 \, (1+q^2 \, r^2)}{r^2}} \\

\dot{\phi} & = & \frac{L \, (1+q^2 \,r^2)^2}{r^2}

\end{array}

[/tex]

where [itex]E, \, A, \, L[/itex] are the invariants of null geodesic motion, and where the choice of sign determines an ingoing or outgoing motion.

We can immediately write down a null geodesic vector field whose integral curves form a null geodesic congruence:

[tex]

\vec{k} =

\frac{1}{(1+q^2 \, r^2)^2} \; \left( E \, \partial_t + A \, \partial_z \right)

\pm \sqrt{ \left( \frac{E^2-A^2}{(1+q^2 \, r^2)^3} - \frac{L^2 \, (1+q^2 \, r^2)}{r^2} \right)} \; \partial_r

+ \frac{L \, (1+q^2 \,r^2)^2}{r^2} \; \partial_\phi

[/tex]

Far from the axis of cylindrical symmetry, this reduces to

[tex]

\vec{k}_{\rm far} = \left( E \, \partial_t + A \, \partial_z \right)

\pm \sqrt{E^2-A^2} \; \partial_r + \frac{L}{r^2} \; \partial_\phi

[/tex]

which corresponds to a plane wave in cylindrical coordinates in flat spacetime, with a constant wave vector pointing in an arbitrary direction (determined by E,A,L). This also shows that for our null geodesics, we can loosely interpret E as "the energy of a photon measured at spatial infinity", and likewise L as "linear momentum parallel to the axis of cylindrical symmetry" and "angular momentum about the axis of cylindrical symmetry". This "photon" doesn't remain "at spatial infinity", but because E, A, L are constants of geodesic motion, these parameters still have meaning as it approaches its point of closest approach to the axis of cylindrical symmetry, even though we interpret them in terms of measurements made when the "photon" was "at spatial infinity".

By using [itex]-1 = ds^2[/itex] instead, you can find first integrals for an arbitrarytimelike geodesicin the Melvin electrovacuum. Using this you can study such interesting phenomena as the light bending and geodesic precession due to the gravitational attraction of the magnetic field, which is strongly concentrated near the symmetry axis.

In spacetimes with smaller Lie groups of self-isometries, it may not be possibly to obtain in such trivial fashion a complete list of first integrals for null, timelike, or spacelike geodesics, but one almost always obtains some valuable information. The case of the Kerr vacuum is remarkable because, although we have only a two dimensional Lie group of self-isometries, and thus we can expect only three first integrals for (null, timelike or spacelike) geodesicsl. In this case, it turns out that there also exists a more subtle symmetry, given by a nontrivial Killing tensor field. Brandon Carter was able to exploit this to solve the geodesic equations for the Kerr vacuum. Trivial Killing tensor fields are common and useless; nontrivial Killing tensor fields seem to be quite rare, in fact, IIRC no example is known other than the Kerr vacuum!

At this point, if you want to examine the Killing equations with which we started, you can try this

The first command instructs Maple to abbreviate partials by subscripts. The second prints out each equation, one line per equation. From this you can see that as I said, this is a first order system of PDEs. Compare with the "triangularized" version output by casesplit!Code (Text):

declare(T(t,z,r,phi), Z(t,z,r,phi), R(t,z,r,phi), Phi(t,z,r,phi);

for guy in keqs do guy od;

I should say that casesplit uses nontrivial Groebner-basis like methods in differential rings, which we can think of as a sophisticated analog of using Gaussian reduction to "triangularize" a system of linear equations. In both cases, the point is that we can systematically solve a triangularized system by solving for the variable for which we have obtained equations involving only that variable, then plugging the result into the equations involving that variable and one other, solving, and so forth until we are done.

There is a general principle in computational algebra which says that the fastest methods tend to be non-deterministic. Thus, casesplit may sometimes present results in a different order, so the directions I gave above must be interpreted with this in mind. Also, casesplit is not perfect and sometimes the above method will fail to find all the solutions it should, so be on the lookout for such problems.

**Physics Forums - The Fusion of Science and Community**

Join Physics Forums Today!

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

# BRS: Solving Killing/Geodesic Eqns with Maple's casesplit

Loading...

Similar Threads - Solving Killing Geodesic | Date |
---|---|

A Solving Schwarzschild Field Equations in this Form | Feb 8, 2018 |

I Solving General Relativity Equations | Jun 24, 2017 |

A Flat s-t 4d killing vectors via solving killing equation | Feb 12, 2017 |

B Solving field equations and schwarzschild metric | Dec 10, 2016 |

**Physics Forums - The Fusion of Science and Community**