I Using Black Holes to Time Travel Into the Future

[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!