Using Black Holes to Time Travel Into the Future

[PeterDonis]

It is also not clear from MTW whether that equation is precise only at the turning points, or everywhere in Schwarzschild coordinates.

It's exact everywhere. It's straightforward to derive that equation from the geodesic equation. I don't have my copy of MTW handy right now, but I'm pretty sure they discuss this.

Any geodesic that is even slightly 'in-going' at 4M should spiral in, because there is no maximum effective potential below that.

Huh? The maximum effective potential continues to increase all the way down to the (unstable) circular photon orbit at $r = 3M$.

in units where $c=1$, $2GM = 1$

I think it should be $GM = 1$. I don't like units where $M$ disappears from the equations for this reason; $M$ is not a physical constant like $G$ and $c$, it's a solution parameter, so to me it needs to remain explicit in the equation.

 

[Ibix]

I'll post the maths and some revised code later (the code is mathematically the same as above, but better laid out to allow multiple trajectories and plots). But here are some results - tangential launches in the equatorial plane from 3M and 4M. Since a trajectory grazing these radii must pass through this point, I think this points to the notion that 4M is escapable but 3M isn't:

 

[PeterDonis]

tangential launches in the equatorial plane from 3M and 4M

I would try values in between. I agree $3M$ shouldn't be escapable, since it's the photon sphere; but $3.001M$ should be, although with a value of $E$ much larger than $1$.

Also, what exactly does $E$ signify in your code? In the usual meaning of that symbol in the effective potential equation, a marginal escape trajectory (i.e., just coming to rest at infinity) should have $E = 1$. But it looks like your $E = 0.9$ trajectory at $r = 4M$ escapes.

 

[Ibix]

I would try values in between. I agree $3M$ shouldn't be escapable, since it's the photon sphere; but $3.001M$ should be, although with a value of $E$ much larger than $1$.

Here's a launch at 3.001M. E=100 does escape or, at least, passes 50M on a nearly-radial path:

Also, what exactly does $E$ signify in your code? In the usual meaning of that symbol in the effective potential equation, a marginal escape trajectory (i.e., just coming to rest at infinity) should have $E = 1$. But it looks like your $E = 0.9$ trajectory at $r = 4M$ escapes.

It does escape (rightly or wrongly). Assuming I didn't muck anything up, I'm following Sean Carroll's GR notes chapter 7, equations 7.43, 7.44, 7.47 and 7.48 with $\epsilon=c^2$ (=1 in his notation) for a timelike path. What I've actually implemented is$$\begin{eqnarray}
\frac{dp}{d\tau}&=&-\frac{GM}{r^2}+\frac{L^2c^2}{r^3}-\frac{3GML^2}{r^4}\\
p&=&\frac{dr}{d\tau}\\
\frac{d\phi}{d\tau}&=&\frac{L}{cr^2}\\
\frac{dt}{d\tau}&=&\frac{E}{c^2-2GM/r}
\end{eqnarray}$$
For initial values you set E, the radial coordinate at launch, R, and the angle $\psi$, which is what an observer hovering at R measures as the launch angle (0 is "straight up"). In this context L is a derived quantity:$$L=\frac{ER\sin\psi}{\sqrt{1-2GM/c^2R}}$$The initial value of $p=dr/d\tau=\left(\cos\psi\right)\sqrt{1-2GM/c^2R}$. See also George Jones' and your comments in an old thread of mine on Schwarzschild null geodesics.

 

[trainman2001]

How much energy would it take to break away from your orbit once you've reached your targeted future time? If you can't escape, what have you gained? I love how you guys can communicate using the math, but it's unintelligible to me. I do well thinking in terms of images and thought experiments, but don't have the math background to understand the language you use to relate your points. It's my problem not yours.

 

[PeterDonis]

What I've actually implemented is

Yes, this looks ok to me.

 

[Jorrie]

For initial values you set E, the radial coordinate at launch, R, and the angle ψ\psi, which is what an observer hovering at R measures as the launch angle (0 is "straight up"). In this context L is a derived quantity:

Looks good, but like Peter, I'm uncomfortable with your E=1 that plunges in. What E are you using for escape energy?
The large E looks like a particle that is very relativistic 'at infinity'. I would like to see what would happen to a particle that would be nearly at rest at infinity, i.e. just above escape energy.

 

[Jorrie]

How much energy would it take to break away from your orbit once you've reached your targeted future time? If you can't escape, what have you gained?

That's where the "zoom-whirl" scenario comes in. The breakaway is "free".

 

[stevendaryl]

Anyway, I'm pretty convinced that using black holes is a pretty wimpy sort of time travel, if you want to do it on the cheap (not spend much energy). It looks like the best you can do is get time dilation factors that are less than a factor of 2 or 3. So you're not going to get 1000 years into the future using relativity.

So broadening the question a little:

Suppose we have two events $e_1$ and $e_2$ that have a timelike separation, and we have two different timelike paths, $\mathcal{P}_1(s)$ and $\mathcal{P}_2(s)$. Let $\tau_1$ be the proper time for the first path and $\tau_2$ be the proper time for the second path.

Question for geodesics: If the two paths are geodesics, is there a limit to how large the ratio $\frac{\tau_1}{\tau_2}$ can be? My original thought was that if one of the geodesics got close enough to a black hole event horizon, while the other stayed away, then the ratio could become arbitrarily large. But it seems that that doesn't work. So is there a limit to the ratio imposed by GR? If there is no limit, then what kind of spacetime geometry could have arbitrarily large ratios?

Now that I think about it, the usual twin paradox, which is non-inertial, could probably be done with geodesics if you had stars arranged in the right pattern. You have the stay-at-home twin, who stays on Earth (or orbiting Earth, to make it geodesic), and you have the traveling twin who uses the gravitational slingshot effect again and again to get to higher and higher velocities, all the while remaining inertial.This would require a very specific arrangement of planets and stars, so it's not really what I was looking for.

Question for noninertial paths: Broadening the question a little: Let's let one of the geodesics be noninertial. Then the ratio can become arbitrarily large, if you're willing to expend energy to force an object to follow one of the paths. The sort of open-ended question is: Is there a "cheapest" way to achieve a given ratio, in terms of energy?

 

[checksix]

Slightly off topic for a physics forum, but an interesting sf novel written back in 1970 by Poul Anderson explored the issues of forward time travel using relativistic speeds: https://en.wikipedia.org/wiki/Tau_Zero

 

[Jorrie]

Maybe you could visit Sir Roger Penrose's "black hole energy extraction city" and hitch a ride on one of the garbage containers! See the figure below that I scanned from Gravitation, Misner, Thorne, Wheeler (MTW), fig.33.2.

It shows a mega-city built on a rigid structure around the equator of a massive, spinning black hole. The city is at just the right distance from the hole so that the local gravity is at a comfortable 1g and that tidal forces and inertial frame dragging do not affect the structures adversely.

The city's garbage is processed at the top left by dumping it into suitable containers and dropping the containers at a carefully chosen angle towards the black hole. Due to frame dragging, the containers swing around the black hole in its ergosphere, and so "steal" some of the black hole's angular momentum. Just before the container starts to leave the ergosphere, the garbage is ejected towards the hole, to be "swallowed" by it.

As I understand it, the garbage essentially attains negative energy as it enters the black hole, causing the hole's rest mass to decrease. The "lost" energy (or most of it) is added to the container as kinetic energy and it can become highly relativistic. I'm not sure what the time dilation ratio would be for a "realistic scenario", but at least it might solve your energy requirements problem. I suggest that you jettison your capsule from the container before the city catches it for converting the extra kinetic energy into electricity. The g-forces might be severe in the catching action!

Then all you need to do is find other black holes, stars, whatever, to shape your trajectory so that you can inertially return to the city. Perhaps then traveling around the city's black hole again, but in a retrograde sling-shot, can decelerate your capsule so that that you can land at the spaceport, much younger than your own generation.

Beautiful sci-fi, but all calculable with GR.

 

[PeterDonis]

I would like to see what would happen to a particle that would be nearly at rest at infinity, i.e. just above escape energy.

I think it's possible that a particle that has just the escape energy will fall in if its periapsis is less than $r = 4M$. For any $r$ less than that, the effective potential is positive for an unstable circular orbit; that translates to energy having to be greater than escape energy for a trajectory that is tangent at periapsis to the same circular orbit.

In other words, I think it's possible that your earlier argument does apply to the case of a particle with energy just equal to escape energy.

 

[PeterDonis]

like Peter, I'm uncomfortable with your E=1 that plunges in. What E are you using for escape energy?

I think the resolution of the apparent "paradox" here is that in the GR Schwarzschild geometry it is not case (as it is in Newtonian gravity) that a trajectory with $E = 1$ will escape regardless of its direction. Heuristically, the range of angles in which an $E = 1$ trajectory will escape gets narrower as the $r$ coordinate decreases. $r = 4M$ marks the point at which a tangential $E = 1$ trajectory no longer escapes, and as $r \rightarrow 2M$, the range of angles at which an $E = 1$ trajectory escapes gets narrower and narrower around the radial outward direction.

I have seen this presented in terms of outgoing light trajectories (for the lightlike case, $r = 3M$ represents the point at which a tangential trajectory no longer escapes), but it would have to work for timelike trajectories as well. I have not worked out the details of how the math shows this, though.

 

[Ibix]

Looks good, but like Peter, I'm uncomfortable with your E=1 that plunges in. What E are you using for escape energy?

I have not worked out the details of how the math shows this, though.

I'm slightly suspicious of my own results, I must say. For example, I agree $E <1$ shouldn't escape because as $r\rightarrow\infty$ $E\rightarrow dt/d\tau$, which in the nearly flat spacetime far from the hole should be $\gamma$. So an escape with $E<1$ implies $\gamma <1$. But my E=0.9 trajectory looks like it escapes. I think I need to do some checking.

Edit: I've either messed something up, the integrator isn't appropriate, or the turnaround on that E=0.9 trajectory is much further from the hole than I guessed (aware I can actually check, not guess!)

 

[Jorrie]

What Peter said makes sense, so your program is probably correct. You can set $\psi$ to zero (i.e. radially or near it) and it should escape with E=1.

 

[stevendaryl]

Slightly off topic for a physics forum, but an interesting sf novel written back in 1970 by Poul Anderson explored the issues of forward time travel using relativistic speeds: https://en.wikipedia.org/wiki/Tau_Zero

I think I would have called it "Beta 1" or "Gamma Infinity" instead of "Tau Zero".

 

[Ibix]

There was a bug - something to do with the way the energy and angular momentum were being set. I've modified the program to set them directly, rather than by specifying an energy and a launch angle. It's not quite so convenient, but it works. Here are launches at 4M, 3M, and 3.001M, which now behave as expected - it does look as if 3M-plus-a-tiny-bit is the closest free-fall approach. I'll have a play with zoom-and-whirl (if I can figure out the initial conditions!) but I can't see how it'll solve Steven's problem if you can't approach closer than 3M.

 

[PeterDonis]

I can't see how it'll solve Steven's problem if you can't approach closer than 3M.

Well, one thing to try would be to find a near-maximal Kerr black hole and look at prograde trajectories in the equatorial plane at marginal escape energy. For a near maximal Kerr hole, the photon sphere gets very close to the horizon, so there should be hyperbolic marginal escape trajectories with periapsis just above the photon sphere that are also very close to the horizon. The time dilation factor at periapsis should be much higher for that case.

(I still think, intuitively, that a single hyperbolic flyby won't accumulate enough time dilation to solve Steven's problem, because the trajectory won't spend enough coordinate time close to periapsis. But it's worth checking since math is better than intuitive guessing.)

 

[Jorrie]

I'll have a play with zoom-and-whirl (if I can figure out the initial conditions!) but I can't see how it'll solve Steven's problem if you can't approach closer than 3M.

You can get them from https://arxiv.org/abs/0802.0459 "A Periodic Table for Black Hole Orbits" by Janna Levin, et al. But whirl-zoom can only give you a factor 2 at best, because for Schwarzschild they need a periapsis of 4M (+ a little). I'm not sure how close one can 'whirl' with a Kerr solution, but they do discuss that, I think.

The other obvious solution is to have some ultra-advanced rocket attached (for orbit changes) and combining that with a Kerr hole. And for fun, add sir Roger's "free energy" city to the mix.

 

[PeroK]

Slightly belatedly, here's a plan with some benefit, I think.

The idea is to use an unstable circular orbit just above $r = 3M$. Initially, you must accelerate up to some $\gamma$ factor well away from the black hole. But, with the correct combination of energy and angular momentum, you should be able to plunge into the unstable orbit, where you will need to use some energy to avoid falling into the hole, or prematurely exiting.

The greater your original energy, the closer to $r=3M$ you can get. You then orbit there for a time and with a small thrust to exit the orbit when you want to come home.

Using conservation of energy, I calculate that the time dilation in the orbit will be approximately 3 times that of the original $\gamma$ factor.

This seems worth doing. E.g. if you consider the relativistic rocket equation. To achieve a $\gamma$ factor of $5$, say, then the rocket would need to be approximately 99% fuel, for the acceleration and deceleration phases. Using the black hole orbit, I believe you could potentially get a $\gamma$ factor of $15$ with the same fuel consumption.

 

[Ibix]

Well, one thing to try would be to find a near-maximal Kerr black hole and look at prograde trajectories in the equatorial plane at marginal escape energy. For a near maximal Kerr hole, the photon sphere gets very close to the horizon, so there should be hyperbolic marginal escape trajectories with periapsis just above the photon sphere that are also very close to the horizon. The time dilation factor at periapsis should be much higher for that case.

(I still think, intuitively, that a single hyperbolic flyby won't accumulate enough time dilation to solve Steven's problem, because the trajectory won't spend enough coordinate time close to periapsis. But it's worth checking since math is better than intuitive guessing.)

It's taken me a while to finish this, for a couple of reasons. First of all, I tried using a straightforward extension of what I was doing for Schwarzschild spacetime (taken from chapter 7 of Carroll's lecture notes), which is to use the conserved Killing quantities to generate equations for $d\phi/d\tau$ and $dt/d\tau$ and then the normalisation of the four-velocity to generate $dr/d\tau$.

In principle, this works for Kerr spacetime. However, the more complicated metric and Killing equations yield a ratio of high-order polynomials for $d^2r/d\tau^2$, which is a recipe for rounding-error induced disaster with a numerical integrator. But I remembered that @m4r35n357 had done this before, and his/her thread cited a paper by Kraniotis that does something clever with Hamiltonians to produce a much friendlier set of equations (equation 13 in that paper). Then the debugging took a while.

I was going to attach a zip file, but apparently you can't do that any more. Instead I've just pasted in (in spoiler tags at the bottom) a Maxima batch file to do some algebraic manipulation of Kraniotis' maths (for interest, really), and a python script that uses the numerical integrator from the scipy library to generate orbits. It also uses matplotlib to display orbits in the Boyer-Lindquist $r$-$\phi$ plane, but that's optional (you'll need to comment out the import, the body of the plot function, and the final call to matplotlib.pyplot.show() if you don't have this library). The code can only handle orbits in the equatorial plane of a Kerr black hole, but extending it to handle out-of-plane orbits wouldn't be hard. The differential equations are defined in the __de() method, and the orbit() method calculates the values of the coefficients as well as initialisation and managing the integrator.

Since we were interested in the best time dilation, I wrote an optimisation function that accepts an energy-per-unit-mass (E) as an input and finds the minimum angular momentum that doesn't fall in. It calculates the elapsed coordinate time (at 50M) divided by the proper time along the path, which is pretty close to the time dilation factor (it ought to be calculated at infinity, but the integrator takes too long to get there ).

Note that E, far from the black hole, can be interpreted as the Lorentz gamma factor relative to the hole (more precisely, relative to an observer hovering at constant Boyer-Lindquist coordinates in the $r\rightarrow\infty$ limit). Given that @stevendaryl was trying to avoid high delta-v engines, this has to be near 1. Unfortunately, that results in fairly modest time dilation factors. Using E=1 the two closest orbits (co- and counter-rotating) around an M=1, a=0.99 black hole give me 1.93 and 1.25 coordinate time at 50M versus the proper time along the path. You can do better by upping E. With E=4, you get $t/\tau$ 15.24 or 6.15 - obviously an improvement over the SR $\gamma$=4, but it requires you to accelerate to around 0.97c. Here are some orbits, co-rotating with different initial energies:

So, unpowered doesn't work. I haven't tried any powered or semi-powered (e.g. "parabolic" infall kicked over to circular orbit then kicked back) orbits, but you could use this integrator to explore circular orbits and connect them to escaping orbits. Perhaps you could change orbits by ejecting a slug of some kind, rather than rocket exhaust, and pick it up again later to recover the braking energy?

Spoiler: Maxima code

/* Kerr geodesics, following "Precise relativistic orbits in Kerr and */
/* Kerr-(anti) de Sitter spacetimes", G V Kraniotis                   */
/* arXiv:gr-qc/0405095v5 13 Oct 2004                                  */

/* Solve equation 2 (the metric) for dt in the case Lambda=0. This */
/* is useful for setting initial values - just choose the spatial  */
/* components of the initial four velocity and the vector type.    */
dssq=Delta/rhosq*(cdt-a*(sin(theta))^2*dphi)^2
     -(rhosq/Delta)*dr^2
     -rhosq*dtheta^2
     -((sin(theta))^2/rhosq)*(a*cdt-(r^2+a^2)*dphi)^2;
ratsimp(solve(%,cdt));
tex(%);

/* Rearrange the cdt and dphi expressions from equation 13 (P is */
/* defined in equation 14) to get E and L. This will give us the */
/* E and L values (and Q, but that's trivial from the dtheta and */
/* equation 11) corresponding to a given initial r, theta and    */
/* four velocity.                                                */
P:E*(r^2+a^2)-a*L;
dteq:rhosq*cdt=((r^2+a^2)/Delta)*P-a*(a*E*(sin(theta))^2-L);
dphieq: rhosq*dphi=(a/Delta)*P-a*E+L/((sin(theta))**2);
collectterms(expand(dteq),E,L);
collectterms(expand(dphieq),E,L);

/* Various definitions - equations 3, 11, 12 and in text below 12. */
/* Note P (equation 14) was defined above).                        */
rhosq:r^2+a^2*(cos(theta))^2;
Delta:r^2+a^2-2*G*M*r/c^2;
R:((r^2+a^2)*E-a*L)^2-Delta*(musq*r^2+(L-a*E)^2+Q);
Theta:Q-(a^2*(musq-E^2)+L^2/(sin(theta))^2)*(cos(theta))^2;

/* Equation 13 - the geodesic equations */
dtdlambda:((r^2+a^2)*P/Delta-a*(a*E*(sin(theta))^2-L))/(c*rhosq);
drdlambda:sqrt(R)/rhosq;
dthetadlambda:sqrt(Theta)/rhosq;
dphidlambda:(a*P/Delta-a*E+L/(sin(theta))^2)/rhosq;

/* Special case: Schwarzschild spacetime, equatorial orbits. */
/* a=0, theta=pi/2 (=> Q=0). We know what the answers are,   */
/* from Carroll 7.43, 7.44, 7.47 and 7.48 - compare. NB:     */
/* Kraniotis' mu^2 (musq, here) is Carroll's epsilon.        */
V:epsilon/2-epsilon*G*M/(c^2*r)+L^2/(2*r^2)-G*M*L^2/(c^2*r^3);
scSubst(f):=block(
    ratsimp(substitute(0,Q,substitute(%pi/2,theta,substitute(0,a,f))))
);
ratsimp(scSubst(dtdlambda)-E/(c-2*G*M/(c*r)));
radcan(substitute(epsilon,musq,scSubst(drdlambda))-sqrt(E^2-2*V));
ratsimp(scSubst(dthetadlambda)-0);
ratsimp(scSubst(dphidlambda)-L/r^2);

/* We are actually interested in equatorial orbits (theta=pi/2, */
/* dtheta/dlambda => Q=0) for general a. Since the dr/dlambda   */
/* is a square root, we will want to know:                      */
/* dt/dlambda                                                   */
/* dphi/dlambda                                                 */
/* (d/dr)((dr/dlambda)^2)/2 (=d^2r/dlambda^2)                   */
eqSubst(f):=block(
    ratsimp(substitute(0,Q,substitute(%pi/2,theta,f)))
);
phide: dphi/dlambda=eqSubst(dphidlambda);
tde: dt/dlambda=eqSubst(dtdlambda);
pde: dp/dlambda=ratsimp(diff(eqSubst(R/(rhosq**2)),r)/2);

/* Rearrange dp/dlambda into a more helpful form */
part(rhs(pde),1,1);                       /* Numerator...               */
collectterms(expand(%),r,r^2);            /* ...grouped by order of r...*/
pde: dp/dlambda=ratsimp(rhs(%th(3))/%)*%; /* ...and back in its place.  */

/* In tex for easy exporting */
tex(phide);
tex(tde);
tex(pde);

Spoiler: python script

"""
KerrOrbits.py

Solves differential equations describing orbits of test particles (i.e.,
particles of negligible mass) in the equatorial plane of a Kerr (uncharged,
rotating) black hole using Boyer-Lindquist coordinates. Covers Schwarzschild
black holes in Schwarzschild coordinates if you set a=0.

Cannot handle orbits outside the equatorial plane.

Based on equation 13 in "Precise relativistic orbits in Kerr and
Kerr-(anti) de Sitter spacetimes", G V Kraniotis, arXiv:gr-qc/0405095v5
13 Oct 2004  

Uses scipy.integrate to solve the equations and matplotlib.pyplot to
display paths.
"""

from __future__ import division

import math,scipy.integrate,matplotlib.pyplot

class EnumError(Exception):
    def __init__(self,varname,varlist):
        Exception.__init__(self,varname+" must be one of "+", ".join(varlist))

# Coding constants (no physics, just labels)
ELQ_INDICES="IE","IL","IQ"
for i in range(len(ELQ_INDICES)): globals()[ELQ_INDICES[i]]=i
TERMINATE_TYPES="FAILED","ESCAPED","CRASHED","TIMED_OUT","UNKNOWN"
for i in range(len(TERMINATE_TYPES)): globals()[TERMINATE_TYPES[i]]=i

# Physical constants - Newton's G and the speed of light, c
PATH_TYPES="SPACELIKE","LIGHTLIKE","TIMELIKE"
for i in range(len(PATH_TYPES)): globals()[PATH_TYPES[i]]=i-1
PATH_DIRECTIONS="INWARD","OUTWARD"
for i in range(len(PATH_DIRECTIONS)): globals()[PATH_DIRECTIONS[i]]=2*i-1
G_SI=6.67e-11
c_SI=3e8

# Evaluate the system of differential equations
class KerrEquatorialSpacetime:
    def __init__(self,M,a,G=1,c=1):
        # Store black hole parameters
        self.__M,self.__a=M,a
        # Store unit constants
        self.__G,self.__c=G,c
    def radiusM(self):
        return self.__G*self.__M/(self.__c*self.__c)
    def __minmaxL(self,E,R):
        M,a,G,c=self.__M,self.__a,self.__G,self.__c
        discrim=math.sqrt(-4*c**2*R**2*G**2*M**2 \
                            +2*c**2*(-R**3*E**2+2*c**2*R**3+a**2*c**2*R)*G*M \
                            +c**4*R**2*(R**2+a**2)*E**2 \
                            -c**4*c**2*R**4 \
                            -a**2*c**4*c**2*R**2)
        L1=(2*a*E*G*M+discrim)/(2*G*M-c**2*R)
        L2=(2*a*E*G*M-discrim)/(2*G*M-c**2*R)
        return (L1,L2)
    def minL(self,E,R):
        return min(self.__minmaxL(E,R))
    def maxL(self,E,R):
        return max(self.__minmaxL(E,R))
    def __de(self,t,y,args=None):
        r=y[1]
        r2=r*r
        r3=r2*r
        rhosq=self.__rho1*r+self.__rho2*r2+self.__rho3*r3
        return [self.__p0/(r2*r2)+self.__p1/r3+self.__p2/r2,
                y[0], \
                (self.__phi0+self.__phi1*r)/rhosq, \
                (self.__t0+self.__t1*r+self.__t3*r3)/rhosq]
    def orbit(self,pathType,pathDirection,E,L,r0,rmax,dtau,taumax):
        # Convert the orbit parameters into the prefactors for the
        # differential equations
        if not pathType in (globals()[pt] for pt in PATH_TYPES):
            raise EnumError("pathType",PATH_TYPES)
        if not pathDirection in (globals()[pd] for pd in PATH_DIRECTIONS):
            raise EnumError("pathDirection",PATH_DIRECTIONS)
        musq=pathType*self.__c*self.__c
        GM=self.__G*self.__M
        self.__rho1=-(self.__a**2)*(self.__c**2)
        self.__rho2=2*GM
        self.__rho3=-(self.__c**2)
        self.__p0=-3*GM*((L-self.__a*E)**2)/(self.__c**2)
        self.__p1=L**2+(self.__a**2)*musq-(self.__a*E)**2
        self.__p2=-musq*GM/(self.__c**2)
        self.__phi0=2*GM*(L-self.__a*E)
        self.__phi1=-(self.__c**2)*L
        self.__t0=self.__phi0*self.__a
        self.__t1=-(self.__a**2)*(self.__c**2)*E
        self.__t3=-(self.__c**2)*E
        # Initialise track with proper time zero and initial
        # values for dr/dtau,r,phi,t in a list
        tau=[0]
        drdtau2=((r0**2+self.__a**2)*E-self.__a*L)**2 \
                -(r0**2+self.__a**2-2*GM*r0/(self.__c**2))*(musq*r0**2+(L-self.__a*E)**2)
        if drdtau2<0:
            print "Negative (dr/dtau)^2",drdtau2," - setting to zero"
            drdtau2=0
        y=[[pathDirection*math.sqrt(drdtau2)/(r0*r0),r0,0,0]]
        # Set up the integrator, and feed it the initial values
        r=scipy.integrate.ode(self.__de).set_integrator("dopri5")
        r.set_initial_value(y[0],tau[0])
        # Run the orbit
        rmin=1.01*(GM/(self.__c**2)+math.sqrt(GM/(self.__c**2)-self.__a**2))
        while r.successful() and rmin<=r.y[1] and r.y[1]<=rmax and r.t<=taumax:
            r.integrate(r.t+dtau)
            tau.append(r.t)
            y.append(r.y)
        # Return a list of proper times and corresponding coordinates
        self.exitReason=UNKNOWN
        if not r.successful():
            self.exitReason=FAILED
        elif rmin>r.y[1]:
            self.exitReason=CRASHED
        elif r.y[1]>rmax:
            self.exitReason=ESCAPED
        elif r.t>taumax:
            self.exitReason=TIMED_OUT
        return [tau,y]

# Create a plot from a list of trajectories
def plot(trajectories,title,rmax,rstep,filename=None):
    fig=matplotlib.pyplot.figure()
    ax=fig.add_subplot(111, projection='polar')
    for trajectory in trajectories:
        tau,y,colour,linewidth,label=trajectory
        phi_path=[p[2] for p in y]
        r_path=[p[1] for p in y]
        ax.plot(phi_path,r_path,colour,linewidth=linewidth,label=label)
    ax.set_rmax(rmax)
    ax.set_rticks([r*rstep for r in range(int(math.ceil(rmax/rstep)))])
    ax.set_title(title)
    ax.set_xticklabels([])
    ax.legend(loc=2,bbox_to_anchor=(0.8,0.15))
    if not filename is None:
        fig.savefig(filename)

def optimiseLaunch(st,E,Lcra,Lesc,froot):
    while True:
        Ltry=(Lesc+Lcra)/2
        RM=st.radiusM()
        orbitData=st.orbit(TIMELIKE,INWARD,E,Ltry,10*RM,15*RM,0.001,1000)
        if st.exitReason==ESCAPED:
            Lesc=Ltry
            ttau=orbitData[1][-1][-1]/orbitData[0][-1]
            print "Escaped - t/tau="+str(ttau)+" L="+str(Ltry)
        elif st.exitReason==CRASHED or st.exitReason==FAILED:
            Lcra=Ltry
            print "Crashed - L="+str(Ltry)
        else:
            print "Failed at L="+str(Ltry)+ \
                        "\nExit reason: "+TERMINATE_TYPES[st.exitReason]
            break
        if abs(Lesc-Lcra)<0.00000001:
            print "Found optimum - building plots"
            orbitData=st.orbit(TIMELIKE,INWARD,E,Lesc,50*RM,55*RM,0.001,1000)
            ttau=orbitData[1][-1][-1]/orbitData[0][-1]
            plot([orbitData+["#000000",1,"L="+str(Lesc)]], \
                "Best: E="+str(E)+", t/tau="+str(ttau),60,10,froot+"1.png")
            plot([orbitData+["#000000",1,"L="+str(Lesc)]], \
                "Best: E="+str(E)+", t/tau="+str(ttau),10,1,froot+"2.png")
            break
def optimiseLaunches(st,E):
    RM=st.radiusM()
    Lcra=0
    for i in range(11):
        Lesc=st.maxL(E,50*RM)*i/10
        print "Trying escape trajectory - L="+str(Lesc)
        orbitData=st.orbit(TIMELIKE,INWARD,E,Lesc,10*RM,15*RM,0.001,1000)
        if st.exitReason==ESCAPED:
            print "Found - beginning optimisation"
            break
    else:
        raise ValueError("Couldn't find non-crashing orbit")
    optimiseLaunch(st,E,Lcra,Lesc,"Best_"+str(E)+"_Max_")
    Lcra=0
    for i in range(11):
        Lesc=st.minL(E,50*RM)*i/10
        print "Trying escape trajectory - L="+str(Lesc)
        orbitData=st.orbit(TIMELIKE,INWARD,E,Lesc,10*RM,15*RM,0.001,1000)
        if st.exitReason==ESCAPED:
            print "Found - beginning optimisation"
            break
    else:
        raise ValueError("Couldn't find non-crashing orbit")
    optimiseLaunch(st,E,Lcra,Lesc,"Best_"+str(E)+"_Min_")
st=KerrEquatorialSpacetime(1,0.99)
optimiseLaunches(KerrEquatorialSpacetime(1,0.99),4)

matplotlib.pyplot.show()

 

[m4r35n357]

Just re-checked my references and, if anyone is interested in the principle (Hamilton-Jacobi analysis), I think I used the equations (14-17) in this link before I even found Kraniotis' paper.

Of course, the whole method is an extension of Wilkins' Kerr spacetime approach from 1972.

Anyway, please pardon the brief hijacking!