Relativistic Astrophysics of Black Holes and Neutron Stars

Introduction

This is a course about black holes and neutron stars. My main goal is to give you facility in computing physical processes in black hole and neutron star spacetimes. I will emphasize general relativistic processes that are relevant to directly observable effects. I will de-emphasize observations and the associated astrophysical details, giving, only enough of the former to motivate a calculation and of the latter to carry one out. I hope that we will build a toolkit useful for exploring the theoretical properties of compact objects and, perhaps, for providing a framework for interpreting the observations. This is not a course about the foundations of relativity or about many of the fundamental theoretical physics questions that have been impacted by carrying out gedanken experiments in a black hole context. I will eschew charged holes, Hawking radiation, singularities, string entropy, two dimensional black holes and so on. It is part of the current excitement of this field that we will have plenty to occupy ourselves dealing with astrophysics-inspired matters. 

In order accomplish my goal without imposing an unreasonable calculational burden on you, I shall make heavy use of Mathematica. I have chosen this program simply because I am most familiar with it. There are alternative programs with which you may be familiar and which are often superior. The course will be structured so that they can be substituted if you want to make the investment. What I want is for you to finish the course with a collection of routines that you understand and which you can use with confidence.The homeworks will often entail performing calculations in this manner and producing suitable graphical output. They will be structured in such a way that you can pass the course with modest investment of effort learning. 

I will start by discussing the Schwarzschild spacetime dealing with familiar material and allowing you to establish your own framework for performing calculations. I shall then go on to discuss the Kerr spacetime, exploring the geometry using alternative coordinate systems and navigating them with null geodesics. I will next generalize these to discuss particle orbits and gravitational lensing. In the fourth week, I plan to discuss thin accretion disks using relatively unsophisticated prescriptions for their angular momentum transport. I hope that we will be able to compute their optical appearance. Thick accretion disks introduce us to relativistic fluid dynamics. I will show you how to compute their structure, again making simple assumptions about their internal thermodynamics etc. This will then allow us to consider two additional fluid problems, wave modes in disks and tidal effects on passing stars. In week 7, I will explain how Maxwell's equations are solved in a Kerr spacetime under the force-free approximation.This will allow us to see how the rotational energy stored in the spacetime can be extracted and put to use. I will next discuss the structure and appearance of spherical neutron stars and outline how to tackle the exterior spacetime of rapidly spinning stars.The course will finish with a relatively unsophisticated discussion of the calculation of gravitational radiation by orbiting stars and compact objects.

Prerequisites for the course are relatively modest. An understanding of Hartle's book will probably suffice. I am more concerned that you have a practical understanding of general relativity, that you can perform simple calculations with speed and confidence - more engineering than theoretical physics!   

Some more texts that should be useful for the course are:
Hartle - Gravity
Misner Thorne Wheeler - Gravitation
Frolov Novikov - Black Hole Physics
Caroll - Spacetime geometry
De Felice  Clarke - Relativity on Curved Manifolds
Plebanski Krasinski - General Relativity and Cosmology

1. Spherical Holes

Schwarzschild Metric

Metric tensor

Use units in which GM = c = 1.
(ds^2)= -(1-2/r)(dt^2)+ (((1-2/r))^(-1))(dr^2)+(r^2) (((dθ)^2) + (sin^2)θ ((dϕ)^2)) 
from which we read off the metric tensor

Clear[coord,metric,inversemetric,affine,riemann,co2on,on2co];
coord={t,r,θ,ϕ};
metric={{-1+2/r,0,0,0},{0,1/(1-2/r),0,0},{0,0,(r^2),0},{0,0,0,(r^2)  ((Sin[θ])^2)}};
Print[MatrixForm[metric]];

Note that time has index 1.

Vectors

Vectors live in the tangent space at a point and can be expanded in components
A = (A^α)eα = (Subscript[A, α])e αetc
which transform according to
(A^(β'))=  (Subscript[(L^(β')), α])(A^α);  (Subscript[A, α])=  (Subscript[A, β'])(Subscript[(L^(β')), α]) etc
The components are variously called vectors and dual vectors or contravariant and covariant components.
Tensors are real functions of vectors 
T(eα,eβ) ≡ (Subscript[(T^α), β]) etc
and their components transform according to the natural generalization of these rules.

Differentiation

The directional derivative of a vector along eαis
(Subscript[∇, (Subscript[e, α])])A≡(Subscript[∇, α]) A =((Subscript[(A^β), ,α])+(A^γ)(Subscript[(Γ^β), γα])) eβ = (Subscript[(A^β), ;α])eβ
         =((Subscript[(Subscript[A, β]), ,α])-(Subscript[A, γ])(Subscript[(Γ^γ), βα])) eβ = (Subscript[(Subscript[A, β]), ;α])eβ
where (Subscript[(A^β), ,α])=(∂_(Subscript[e, α]))(A^β). 
The connection for a coordinate basis is 
(Subscript[(Γ^α), βγ])= ((e )^α)· (Subscript[∇, (Subscript[e, β])])(Subscript[e , γ])= (1/2) (g^(αδ))((Subscript[g, δβ,γ])+(Subscript[g, δγ,β]) - (Subscript[g, βγ,δ]))
and we have used the inverse of the metric tensor. The derivative for a non-coordinate basis contains extra terms coming from the non-zero commutator.

inversemetric:=inversemetric=Simplify[Inverse[metric]];
Print[MatrixForm[inversemetric]];
affine:=affine=Simplify[Table[(1/2)*Sum[(inversemetric[[α,δ]])*
(D[metric[[δ,β]],coord[[γ]] ]+
D[metric[[δ,γ]],coord[[β]] ]-D[metric[[β,γ]],coord[[δ]] ]),{δ,1,4}],
{α,1,4},{β,1,4},{γ,1,4}] ];
Print[MatrixForm[affine]]

We can check that the covariant derivative of the metric tensor vanishes, as it should.
(Subscript[g, αβ;γ])=(Subscript[(Subscript[g, αβ]), ,γ])-(Subscript[g, αδ])(Subscript[(Γ^δ), βγ])-(Subscript[g, δβ])(Subscript[(Γ^δ), αγ])

delmetric:=delmetric=
Simplify[Table[D[metric[[α,β]],coord[[γ]]]-
Sum[metric[[α,δ]]affine[[δ,β,γ]]+metric[[δ,β]]affine[[δ,α,γ]],{δ,1,4}],{α,1,4},{β,1,4},{γ,1,4}]]
Print[MatrixForm[delmetric]];

Geodesic Motion

Stationarity of the proper time implies 
δ ∫ ds L = δ ∫ ds[-(Subscript[g, αβ])((dx^α)/ds)((((dx^β)/ds)])^(1/2))=0
where L is the Lagrangian. Performing the variation, we obtain the equation of motion of a freely-falling particle moving with velocity u 
(Subscript[∇, u])u = 0
where (u^α)= (dx^α)/ dτ with u · u = -1 and so we get the geodesic equation
(((d^2)(x^α))/((dτ)^2))+(Subscript[(Γ^α), βγ])((dx^β)/(dτ))((dx^γ)/(dτ))=0
For photons τ is replace by an affine parameter λ, chosen so that u=((dx)/(dλ))and  u · u = 0; we can think of λ as (lim_(m->0))τ/m though only defined up to multiplicative and additive constants.

Killing's Equation

If there is a symmetry, a direction ξ along which the metric does not change, then it satisfies (Subscript[ξ, α;β]) + (Subscript[ξ, β;α]) = 0.
Lagrange's equations then imply that ξ · u is conserved along the geodesic. This gives us constants of the motion. For motion in the equatorial plane, t and ϕ symmetries imply E ≡ -(Subscript[u, t]) and  L ≡ (Subscript[u, ϕ]) are constant.

Freefall

As an example, consider the geodesic motion of a particle falling from rest at infinity. We have (1-2/r)dt/dτ ≡E=1. We also know that 
d(τ^2)= -(ds^2)= (1-2/r)(dt^2)-(((1-2/r))^(-1))(dr^2)
from which we obtain (u^r)≡ dr/dτ = -(((2/r))^(-1/2)). The full velocity is then given by ensuring that (u^α)(Subscript[u, α]) = -1.
(u^α) = [(((1-2/r))^(-1)),-(((2/r))^(1/2)), 0, 0]
The speed is β = (((2/r))^(1/2)) similar to the Newtonian expression and the associated Lorentz factor is γ = (1-2 / r)-1/2.

Riemann and Ricci tensors

The law of geodesic deviation for a connection vector n is
(Subscript[∇, u])(Subscript[∇, u])n +R[,u,n,u]=0
or in components,
((((Subscript[∇, u])(Subscript[∇, u]) n))^α) + (Subscript[(R^α), βγδ])(u^β)(n^γ)(u^δ)=0
The components of the Riemann tensor R are defined in a coordinate basis by
(Subscript[(R^α), βγδ])=(Subscript[(Γ^α), βδ,γ]) - (Subscript[(Γ^α), βγ,δ]) + (Subscript[(Γ^α), ϵγ])(Subscript[(Γ^ϵ), βδ]) - (Subscript[(Γ^α), ϵδ])(Subscript[(Γ^ϵ), βγ]).

riemann:=riemann=
Simplify[Table[D[affine[[α,β,δ]],coord[[γ]]]-D[affine[[α,β,γ]],coord[[δ]]]+Sum[affine[[α,ϵ,γ]]*affine[[ϵ,β,δ]]-affine[[α,ϵ,δ]]affine[[ϵ,β,γ]],{ϵ,1,4}],{α,1,4},{β,1,4},{γ,1,4},{δ,1,4}]]
Print[MatrixForm[riemann]];

The Ricci tensor is defined by
(Subscript[R, αβ]) = (Subscript[(R^γ), αγβ]) 
All elements vanish as they must following the field equations.

ricci:=ricci=Simplify[Table[Sum[riemann[[γ,α,γ,β]],{γ,1,4}],{α,1,4},{β,1,4}]]
Print[MatrixForm[ricci]];

Riemann tensor in freefalling frame

First transform the Riemann tensor into an orthonormal basis, at rest in the Schwarzschild frame
 (Subscript[(R^(OverHat[α])), (OverHat[β])(OverHat[γ])(OverHat[δ])])=(Subscript[(L^(OverHat[α])), ϵ])(Subscript[(R^ϵ), ζηθ])(Subscript[(L^ζ), (OverHat[β])])(Subscript[(L^η), (OverHat[γ])])(Subscript[(L^θ), (OverHat[δ])])

onriemann:=onriemann=
Table[Sum[on2co[[α,ϵ]]riemann[[ϵ,ζ,η,θ]]co2on[[β,ζ]]co2on[[γ,η]]co2on[[δ,θ]],{ϵ,1,4},{ζ,1,4},{η,1,4},{θ,1,4}],{α,1,4},{β,1,4},{γ,1,4},{δ,1,4}]
Print[MatrixForm[onriemann]]

Next, boost into the freefalling frame with infall speed β by introducing the transformation matrices (Subscript[(L^(OverTilde[α])), (OverHat[β])]) and its inverse (Subscript[(L^(OverHat[α])), (OverTilde[β])])  where the ~ denotes the freefall coordinates.

on2ff={{gammaff,-gammaff betaff,0,0},
       {-gammaff betaff,gammaff,0,0},
       {0,0,1,0},
       {0,0,0,1}}/.{betaff->(((2/r))^(1/2)),gammaff->(((1-2/r))^(-1/2))};
ff2on=Simplify[Inverse[on2ff]];
Print[MatrixForm[on2ff],MatrixForm[ff2on]];

Finally, we boost from the rest orthonormal basis to the freefall orthonormal basis
 (Subscript[(R^(OverTilde[α])), (OverTilde[β])(OverTilde[γ])(OverTilde[δ])])=(Subscript[(L^(OverHat[α])), (OverHat[ϵ])])(Subscript[(R^(OverHat[ϵ])), (OverHat[ζ])(OverHat[η])(OverHat[θ])])(Subscript[(L^(OverHat[ζ])), (OverTilde[β])])(Subscript[(L^(OverHat[η])), (OverTilde[γ])])(Subscript[(L^(OverHat[θ])), (OverTilde[δ])])

ffriemann:=ffriemann=Simplify[Table[Sum[
ff2on[[α,ϵ]]onriemann[[ϵ,ζ,η,θ]]on2ff[[β,ζ]]on2ff[[γ,η]]on2ff[[δ,θ]],{ϵ,1,4},{ζ,1,4},{η,1,4},{θ,1,4}],{α,1,4},{β,1,4},{γ,1,4},{δ,1,4}]];
Print[MatrixForm[ffriemann]]

This is non-singular at the horizon.

Photon Orbits

Motivation

The flux that we measure from a neutron star or a disk will be influenced by the path that the photons have followed from the sources to us. Let us first consider this problem in a Schwarzschild spacetime as a prelude to repeating the exercise in a Kerr spacetime. We have not yet been able to image the spacetime directly, though there are futuristic plans to do so using submm VLBI and X-ray interferometry.

Gravitational Redshift

From symmetry, orbits will be confined to a plane and there will be two nontrivial constants of the motion associated with the t and ϕ symmetries. (For photons, we use the momentum p instead of the velocity with no change in the argument.)
Now the photon energy measured by an observer moving with velocity u is ω = -u·p. For an observer at rest on the surface of a star (u^α)= [(((1-2/r))^(-1/2)),0,0,0]. Therefore, u = (((1-2/r))^(-1/2)) (ξ^((t))) where (ξ^((t))) is the time Killing vector. The energy of a photon on the surface of a star is (SuperStar[ω])= -(((1-2/(SuperStar[r])))^(-1/2)) (ξ^((t)))·p. However, along a geodesic, E ≡ -(ξ^((t)))·p is conserved. Therefore 
(ω^∞)= (((1-2/(SuperStar[r])))^(1/2))(SuperStar[ω]).  
We can also use this relation to transform the intensity using the fact that the occupation number is invariant
 (Subscript[I, ω]^∞) =(((1-2/(SuperStar[r])))^(3/2)) (Subscript[I, ω]^( * )). 
 Note that this intensity must be computed at the transformed frequency.

Null Geodesic Equation

For non-radial motion, there is a Killing vector  (ξ^((ϕ))) and a second constant L = (ξ^((ϕ)))·p .  We can therefore write down two equations  of motion assuming θ = π/2:
(dt/(dλ)) = (((1-2/r))^(-1))E
((dϕ)/(dλ)) = (r^(-2)) L
where λ is the affine parameter. We obtain the radial equation from the interval (Subscript[g, αβ])d(x^α)d(x^β) = 0. This gives:
((((dr/(dλ))))^2)=(E^2) - (1-2/r)(L^2) /(r^2)
Now change Lλ->λ; L/E->b to get
(dt/(dλ)) = (b^(-1)) (((1 - 2/r))^(-1))
((dϕ)/(dλ)) = (r^(-2))
((((dr/(dλ))))^2) = 1/(b^2) - V(r)
where V(r) = (r^(-2)) - 2(r^(-3)). Note that far from the star (r^2)dϕ/dt -> b, which identifies b as the impact parameter.
The effective potential V(r) has the form:

veff[r_]:=(1/(r^2))-(2/(r^3));
Plot[veff[r],{r,2,10},Frame->True,PlotRange->{{0,10},{0,0.05}},FrameLabel->{"r","V(r)"},RotateLabel->False]

From this we can deduce that the cross section for capture of an infalling photon is π/V((SuperStar[r])). (Note that the capture cross-section for a hole is 27π.) 
We want to map the stellar surface onto the observer's sky. To do this introduce u=1/r and observe that
((((du/(dϕ))))^2) = 1/(b^2)- (u^2)(1-2u)
and we can obtain the longitude of emission by quadrature

rstar=5;
long[b_]:=NIntegrate[(1/(Sqrt[(1/(b^2))-(x^2) (1-2 x)])),{x,0,(1/rstar)}]
Plot[long[b],{b,0.001`,(1/(Sqrt[veff[rstar]]))},Frame->True,PlotRange->{{0,8},{0,3}},FrameLabel->{"b","ψ"},RotateLabel->False]

Observations

Doeleman et al 2008 Nature measure an angular diameter at 1.3mm of 37μas.

Ingoing Eddington-Finkelstein Coordinates

This coordinate system is generally introduced to clarify the nature of the Schwarzschild non-singularity. We will use it as a convenience for numerical computations which we don't want to blow up near the horizon. The coordinates are Schwarzschild except that t->v = t+r+2 log(r-2). The metric becomes 
(ds^2)= -(1-2/r)(dv^2)+ 2dvdr +(r^2) (((dθ)^2) + (sin^2)θ ((dϕ)^2)) 
Radial light rays satisfy
-(1-2/r)(dv^2) + 2dv dr=0
Ingoing photons follow v= v(0); outgoing photons v = 2[r - (Subscript[r, 0])+ 2log {(r-2)/((Subscript[r, 0])-2)}]+v(0). 
In order to see its significance, let us first consider the trajectories of two photons, one moving radially inwards from r=3 at t=0, the other moving outward in Schwarzschild coordinates.

tin=NDSolve[{t'[r]==-(r/(r-2)),t[3]==0},t,{r,3,2.01}];
tout=NDSolve[{t'[r]==(r/(r-2)),t[3]==0},t,{r,3,10}];
inplot=Plot[t[r]/.tin,{r,3,2.01},Frame->True,AspectRatio->Automatic,PlotRange->{{0,10},{-1,10}},FrameLabel->{"r","t","Schwarzschild",""},RotateLabel->False,DisplayFunction->Identity];
outplot=Plot[t[r]/.tout,{r,3,10},DisplayFunction->Identity];
Show[inplot,outplot,DisplayFunction->$DisplayFunction]

The ingoing photon never reaches the horizon.
Now repeat the exercise in Eddington-Finkelstein coordinates.

vin=NDSolve[{v'[r]==0,v[3]==0},v,{r,3,2}];
vout=NDSolve[{v'[r]==((2 r)/(r-2)),v[3]==0},v,{r,3,10}];
inplot=Plot[v[r]/.vin,{r,3,2},Frame->True,AspectRatio->Automatic,PlotRange->{{0,10},{-1,10}},FrameLabel->{"r","v","Eddington-Finkelstein",""},RotateLabel->False,DisplayFunction->Identity];
outplot=Plot[v[r]/.vout,{r,3,10},DisplayFunction->Identity];
Show[inplot,outplot,DisplayFunction->$DisplayFunction]

The ingoing photon passes straight through the horizon. It is conventional to replot this using (OverTilde[v]) = v-r as the ordinate.

inplot=Plot[v[r]+3-r/.vin,{r,3,2},Frame->True,AspectRatio->Automatic,PlotRange->{{0,10},{-1,10}},FrameLabel->{"r","\!\(\*OverscriptBox[\(v\), \(~\)]\)","Eddington-Finkelstein",""},RotateLabel->False,DisplayFunction->Identity];
outplot=Plot[v[r]+3-r/.vout,{r,3,10},DisplayFunction->Identity];
Show[inplot,outplot,DisplayFunction->$DisplayFunction]

Problems 1

1. Killing Equation

Explain Killing's equation and show that ξ·u is conserved along a geodesic.

2. Fun? with Riemann tensors

This is an exercise to coerce you to start writing your own software and develop your own conventions for calculating. It also brings out an important point.  We computed the Riemann tensor in the freefall frame by a two stage process - first transform to the orthonormal basis at rest in the Schwarzschild coordinates, then Lorentz boost to the orthonormal freefall frame. We can actually do it in one step.
a. Write your own software for computing R in Schwarzschild coordinates starting from the metric tensor
b. Observe that the velocity u is (1,0,0,0) in the free fall orthonormal basis. It is therefore the time component of the freefall basis vector set, (Subscript[e,  (OverTilde[0])]).What are the components of the vector (Subscript[e,  (OverTilde[0])])in Schwarzschild coordinates.
c. Use the fact that the other three freefall basis vectors must be normalized and orthogonal to (Subscript[e,  (OverTilde[0])])to determine their components in Schwarzschild coordinates. Hint: Use symmetries to construct this basis
d. Now construct the freefall components of the Riemann tensor directly by evaluating 
(Subscript[R, (OverTilde[α])(OverTilde[β])(OverTilde[γ])(OverTilde[δ])]) = ((((Subscript[e, (OverTilde[α])])))^α) ((((Subscript[e, (OverTilde[β])])))^β) ((((Subscript[e, (OverTilde[γ])])))^γ) ((((Subscript[e, (OverTilde[δ])])))^δ)(Subscript[R, αβγδ]).

3. Appearance of a Neutron Star

EXO 0748-676 is a neutron star. The Fe line normally seen at 1.27keV is now observed at 0.94keV. The observed continuum spectrum may be crudely modeled as a black body with temperature 1.6 keV. (See astro-ph 0303147)
a. Compute the radius, assuming that the neutron star has a "standard" mass of 1.4 times that of the sun.
b. Compute the intensity measured on the stellar surface assuming that it is isotropic and uniform. In other words there is no limb-brightening or darkening and there are no "hot spots" on the surface.
c. Compute the total luminosity emitted by this star as measured at large distance
d. Assume that a local geographer has painted lines of latitude and longitude on the surface of the star and 
a future X-ray interferometer can resolve the surface. Compute the image you would see.

2. Spinning Holes

Kerr metric

Boyer-Lindquist coordinates

The Kerr metric in Boyer-Lindquist coordinates is
(ds^2)= -(1-((2r)/(ρ^2)))(dt^2) - ((4 a r (sin^2)θ )/(ρ^2))  dϕ dt +((ρ^2)/Δ) (dr^2)+ ρ2 ((dθ)^2) + (((Σ^2)(sin^2) θ)/((ρ^2)))((dϕ)^2)
     = -(λ^2)(dt^2)  +(ρ^2)( ((dr^2)/Δ)+ ((dθ)^2)) + (((ϖ^2)(dϕ-ωdt))^2)
where (ρ^2)= (r^2)+(a^2)(cos^2)θ,  Δ = (r^2)-2r + (a^2), (Σ^2)=((((r^2)+(a^2)))^2)-Δ (a^2)(sin^2)θ  and the lapse, shift and cylindrical radius are  λ=(((ρΔ)^(1/2))/Σ), ω=((2 a r)/(Σ^2)),ϖ=((Σ sinθ)/ρ).
The parameter a measures the specific angular momentum of the hole and lies in [0,1] assuming cosmic censorship.
For derivations and discussions of uniqueness, see Chandrasekhar or Isham. For lapse and shift, see MacDonald and Thorne 1982  MNRAS 198 345.

Clear[coord,metric,inversemetric,affine,riemann,ricci,a];
coord={t,r,θ,ϕ};
ρ=((((r^2)+(a^2)((Cos[θ])^2)))^(1/2));
Δ=(r^2)-2r+(a^2);
Σ = (((((((r^2)+(a^2)))^2)-Δ (a^2)((Sin[θ])^2)))^(1/2));
λ=((ρ (Δ^(1/2)))/Σ);
ω=((2a r)/(Σ^2));
ϖ=((Σ Sin[θ])/ρ);
metric={{-(1-((2r)/(ρ^2))),0,0,((-2a r ((Sin[θ])^2))/(ρ^2))},{0,((ρ^2)/Δ), 0,0},{0,0,(ρ^2),0},{((-2a r ((Sin[θ])^2))/(ρ^2)),0,0,((Σ^2)/(ρ^2))((Sin[θ])^2)}};
Print[MatrixForm[metric]];

The metric is asymptotically flat and becomes Schwarzschild when a=0. 
Note that there are two Killing vectors in the t, ϕ directions. We say that the metric is stationary and axisymmetric.
The event horizon is a 3D (null) surface. It is located where Δ=0 or r=1+(((1-(a^2)))^(1/2)).

a=.95;
rh=1+(((1-(a^2)))^(1/2));Ωh=(a/(2rh));
holeplt=Graphics[{Black,Disk[{0,0},1.02rh,{0,π/2}]}];
Δplt=ContourPlot[Tooltip[Δ/.{r->((((x^2)+(y^2)))^(1/2)),θ->ArcTan[y,x]}],{x,0,4},{y,0,4},PlotRange->Full,Contours->{.5,1,2,4,8,16},ContourShading->Automatic,PlotLabel->Style["Delta",Black]];
aplt=Graphics[{Text[Style["a=",Black],{2,3}],Text[Style[a,Black],{3,3}]}];
Δplt=Show[Δplt,holeplt,aplt];ϖplt=ContourPlot[Tooltip[ϖ/.{r->((((x^2)+(y^2)))^(1/2)),θ->ArcTan[y,x]}],{x,.01,4},{y,0,4},PlotRange->Full,Contours->{1,2,3,4,5},ContourShading->Automatic,PlotLabel->Cyl. radius];
ϖplt=Show[ϖplt,holeplt];
ωplt=ContourPlot[Tooltip[If[r>rh,ω,Ωh]/.{r->((((x^2)+(y^2)))^(1/2)),θ->ArcTan[y,x]}],{x,0,4},{y,0,4},PlotRange->Full,Contours->{.025,.05,.1,.2,.3,.4,.5},ContourShading->Automatic,PlotLabel->Shift];
ωplt=Show[ωplt,holeplt];
λplt=ContourPlot[Tooltip[If[r>rh,λ,0]/.{r->((((x^2)+(y^2)))^(1/2)),θ->ArcTan[y,x]}],{x,0,6},{y,0,6},PlotRange->Full,Contours->{.1,.2,.3,.4,.5,.6,.7,.8,.9},ContourShading->Automatic,PlotLabel->Lapse];
λplt=Show[λplt,holeplt];
Show[GraphicsGrid[{{Δplt,ϖplt},{ωplt,λplt}}]]
Clear[a];

Riemann and Ricci

We will need the inverse metric coefficients, 
(g^tt)=-((Σ^2)/((ρ^2)Δ)),  (g^(tϕ))=-((2ar)/((Δρ)^2)), (g^rr)=(Δ/(ρ^2)), (g^(θθ))= (1/(ρ^2)), (g^(ϕϕ))=(((ρ^2)-2r)/(((Δρ)^2)(sin^2)θ))
the connection and the Riemann tensor.

inversemetric=FullSimplify[Inverse[metric]]

affine=FullSimplify[Table[Sum[(inversemetric[[α,δ]])
(D[metric[[δ,β]],coord[[γ]]]+D[metric[[δ,γ]],coord[[β]]]-D[metric[[β,γ]],coord[[δ]]]),{δ,1,4}]/2,{α,1,4},{β,1,4},{γ,1,4}]];
Print[MatrixForm[affine]];

riemann=FullSimplify[Table[D[affine[[α,β,δ]],coord[[γ]]]-D[affine[[α,β,γ]],coord[[δ]]]+Sum[affine[[α,ϵ,γ]]affine[[ϵ,β,δ]]-affine[[α,ϵ,δ]]affine[[ϵ,β,γ]],{ϵ,1,4}],{α,1,4},{β,1,4},{γ,1,4},{δ,1,4}]];
Print[MatrixForm[riemann]];

Again the vanishing of the Ricci tensor is a sanity check.

ricci=Simplify[Table[Sum[riemann[[γ,α,γ,β]],{γ,1,4}],{α,1,4},{β,1,4}]];
Print[MatrixForm[ricci]];

Why not Kerr-Newman?

Gravitationally significant charges discharge immediately.

Event Horizon

Location

The event horizon, located where (Subscript[g, rr])→ ∞,

f[x_]:=Δ/.r->x
soln=Solve[f[x]==0,x]

Clear[rh,a];rh[a_]:=1+(Sqrt[1-(a^2)])

at r = (SubPlus[r])= 1 + (1 - (((a^2)))^(1/2)),

Geometry of the horizon

The event horizon is a null, 3-surface, separating the spacetime region from which there is a ray to infinity from one where there is no such ray. It is generated by a null vector field lying in the t-ϕ plane. We will be interested in three subspaces spanned by t-θ-ϕ, t-ϕ, θ-ϕ coordinates with metrics:

tθϕmetric=Simplify[
{{metric[[1,1]],metric[[1,3]],metric[[1,4]]},{metric[[3,1]],metric[[3,3]],metric[[3,4]]},{metric[[4,1]],metric[[4,3]],metric[[4,4]]}}
/.r->1+ (Sqrt[1-(a^2)])];
Print[MatrixForm[tθϕmetric]];
tϕmetric=Simplify[
{{metric[[1,1]],metric[[1,4]]},{metric[[4,1]],metric[[4,4]]}}/.r->1+ (Sqrt[1-(a^2)])];
Print[MatrixForm[tϕmetric]];
θϕmetric =Simplify[
{{metric[[3,3]],metric[[3,4]]},{metric[[4,3]],metric[[4,4]]}}
/.r->1+ (Sqrt[1-(a^2)])];
Print[MatrixForm[θϕmetric]]

Angular Velocity

Let the null vector field in the t-ϕ space have components [1,(Subscript[Ω, H])]. Solve for (Subscript[Ω, H])=a/2(Subscript[r, H]) <1/2, the hole's angular velocity.

f[x_]:=Simplify[{1,Ω}.tϕmetric.{1,Ω}]/.Ω->x
soln=FullSimplify[Solve[f[x]==0,x]]

Ωh[a_]:=(a/(2rh[a]))
Plot[Ωh[a],{a,0,1}, Frame->True,FrameLabel->{"a","Ω"},RotateLabel->False]
Export

Photons in the horizon propagate with this angular velocity independent of θ.

Curvature

Now take a 2D slice at constant t. Let us look at the metric coefficients.

a=1;
gθθ[θ1_,a1_]:=θϕmetric[[1,1]]/.{θ->θ1,a->a1}
gϕϕ[θ1_,a1_]:=θϕmetric[[2,2]]/.{θ->θ1,a->a1}
Plot[{gθθ[θ,a],gϕϕ[θ,a]((Csc[θ])^2)},{θ,0,0.5 π},Frame->True,PlotRange->{{0,0.5 π},{0,5}},FrameLabel->{θ,"\!\(\*SubscriptBox[\(g\), \(θθ\)]\),\!\(\*SubscriptBox[\(g\), \(ϕϕ\)]\)\!\(\*SuperscriptBox[\(csc\), \(2\)]\)θ"},RotateLabel->False,DisplayFunction->Identity]
Clear[a];

The internal geometry of the timeslice is anisotropic and not embedable in E3 for a>(Sqrt[3])/2 when the polar curvature becomes negative.  To see this, compute the two dimensional Riemann and Ricci tensors

θϕcoord={θ,ϕ};
inverseθϕmetric=Simplify[Inverse[θϕmetric]];
θϕaffine=Simplify[Table[Sum[(inverseθϕmetric[[i,s]])
(D[θϕmetric[[s,j]],θϕcoord[[k]]]+D[θϕmetric[[s,k]],θϕcoord[[j]]]-D[θϕmetric[[j,k]],θϕcoord[[s]]]),{s,1,2}]/2,{i,1,2},{j,1,2},{k,1,2}]];
θϕriemann=Simplify[Table[D[θϕaffine[[i,j,l]],θϕcoord[[k]]]-D[θϕaffine[[i,j,k]],θϕcoord[[l]]]+Sum[θϕaffine[[i,s,k]]θϕaffine[[s,j,l]]-θϕaffine[[i,s,l]]θϕaffine[[s,j,k]],{s,1,2}],{i,1,2},{j,1,2},{k,1,2},{l,1,2}]];
θϕricci=Simplify[Table[Sum[θϕriemann[[k,i,k,j]],{k,1,2}],{i,1,2},{j,1,2}]];
Print[MatrixForm[θϕricci/.{a->.99,θ->0}]];

Area

We can compute the area (Subscript[A, H])= 8(Subscript[πr, H])from

FullSimplify[4Pi Integrate[Sqrt[θϕmetric[[1,1]]θϕmetric[[2,2]]],{θ,0,Pi/2}]]
ah[a_]:=8 Pi rh[a]

Ergosphere

Photon Orbits

Consider photons moving instantaneously azimuthally outside the horizonon circular, equatorial orbits. Let Ω=dϕ/dt. The interval gives us (Subscript[g, ϕϕ])(Ω^2)+2(Subscript[g, ϕ0])Ω+(Subscript[g, 00])=0. Solving for the two angular velocities,

Ω1=Simplify[-(metric[[1,4]]+Sqrt[metric[[1,4]]^2-metric[[1,1]]metric[[4,4]]])/metric[[4,4]]];

Ω2=Simplify[-(metric[[1,4]]-Sqrt[metric[[1,4]]^2-metric[[1,1]]metric[[4,4]]])/metric[[4,4]]];
Ωmin[r1_,a1_]:=Ω1/.{r->r1,a->a1,θ->(π/2)}
Ωmax[r1_,a1_]:=Ω2/.{r->r1,a->a1,θ->(π/2)}
Ωplt[a_]:=Plot[{Ωmin[r,a],Ωmax[r,a],Ωh[a]},{r,rh[a],6},Frame->True,PlotRange->{{rh[a],6},{-.2,.5}},FrameLabel->{"r","Ω"}]
Ωplt[.6]

These two angular frequencies define a range in which azimuthally moving timelike geodesics exist. Note that (Subscript[Ω, min])=(Subscript[Ω, max])=(Subscript[Ω, H]) at the horizon.

Static Limit

The boundary of the region within which there are no stationary, azimuthal timelike orbits is defined by (Subscript[g, 00])(r,θ;a)=0 or (r^2)=(ρ^2).  Sometimes this is called the static limit.

rergo[θ_,a_]:=1+(Sqrt[1-(a^2) ((Cos[θ])^2)])

Plotting these up,

a=1;
ergoplt[a_]:=PolarPlot[rergo[(π/2)-θ,a],{θ,0,2 π},Frame->True,AspectRatio->Automatic,DisplayFunction->Identity]
hplt[a_]:=Graphics[Disk[{0,0},rh[a]]]
Show[{ergoplt[a],hplt[a]},DisplayFunction->$DisplayFunction]
Clear[a];

The region within the static limit is called the "ergosphere".

Hawking Area Theorem

The theorem states that the area of the horizon must increase. (Recall that the horizon depends upon the future spacetime which must be asymptotically flat.) Introduce irreducible mass  (Subscript[m, i])=(Sqrt[(Subscript[A, H])/16π])=(Sqrt[a/4(Subscript[Ω, H])]). Reinstating the gravitational mass m we have that ((Subscript[m, i])^2)=m ( m+(Sqrt[(m^2)-(a^2)]))/2. Equivalently, (((m^2)=((Subscript[m, i])^2)+S)^2)/4((Subscript[m, i])^2), where S=m a is the spin
or ((m=(Subscript[m, i ])(1-4((Subscript[Ω, H])^2)((Subscript[m, i])^2)))^(-1/2)). (Note the whimsical similarity to special relativity!) The extractable mass is given by (Subscript[m, r])=m-(Subscript[m, i]).

ParametricPlot[{Ωh[a],1-(Sqrt[(a/(4 Ωh[a]))])},{a,0,1},Frame->True,PlotRange->{{0,0.5},{0,0.3}},FrameLabel->{(Subscript[Ω, H]),((Subscript[m, r])/m)},RotateLabel->False]

Up to 0.29 of the gravitational mass is extractable.

ZAMO Frame

Angular Velocity

Physical observers must rotate with respect to infinity close to the hole. A particular family of rotating frames at constatn r, θ, have  Zero Angular MOmentum (also known as FIDucial Observers). 
Consider non-geodesic circular motion with (Subscript[u, ϕ])=0. This is the angular momentum conserved along a geodesic and associated with the ϕ Killing vector. Using the length of the four velocity we 
(Subscript[u, α])= [- (((-(g^tt)))^(-1/2)),0,0,0]
Compute its angular velocity 
(Subscript[Ω, Z]) =  ((u^ϕ)/(u^t))=((2ar)/(((((r^(2 ))+(a^2)))^2)-((Δa)^2)(sin^2)θ))

Clear[a];ublcov=FullSimplify[{-(((Sqrt[Δ]) ρ)/Σ),0,0,0}];
ublcon=FullSimplify[{(Σ/((Sqrt[Δ]) ρ)),0,0,((2 a r)/(Σ (Sqrt[Δ]) ρ))}];
FullSimplify[inversemetric.ublcov-ublcon]
FullSimplify[(∑_(α=1)^4)(∑_(β=1)^4)ublcon[[α]] metric[[α,β]] ublcon[[β]]];

These are checks on the algebra. Now compute the angular velocity.

FullSimplify[((ublcon[[4]])/(ublcon[[1]]))]

Therefore,

Ωz[r_,θ_,a_]:=((2 a r)/(((((r^2)+(a^2)))^2)-((r^2)-2 r+(a^2)) (a^2) ((Sin[θ])^2)))

It turns out that (Subscript[Ω, Z])= (1/2)((Subscript[Ω, 1])+(Subscript[Ω, 2])). This is therefore a useful frame to carry out physics calculations.

FullSimplify[(1/2)(Ω1+Ω2)-Ωz[r,θ,a]]

Next plot (Subscript[Ω, Z]) for different values of r, θ, a.

a=.5;θ=(π/2);
plt1=Plot[Ωz[r,θ,a],{r,rh[a],6},Frame->True,FrameLabel->{"r","\!\(\*SubscriptBox[\(Ω\), \(z\)]\)"},RotateLabel->False,DisplayFunction->Identity]
Clear[a,θ];

a=.95;r=3;
plt2=Plot[Ωz[r,θ,1],{θ,0,(π/2)},Frame->True,FrameLabel->{"θ","\!\(\*SubscriptBox[\(Ω\), \(z\)]\)"},RotateLabel->False,DisplayFunction->Identity]
Clear[a,r];

r=3;θ=(π/2);plt3=Plot[Ωz[3,(π/2),a],{a,0,1},Frame->True,FrameLabel->{"a","\!\(\*SubscriptBox[\(Ω\), \(z\)]\)"},RotateLabel->False,DisplayFunction->Identity]
Clear[r,θ];

Another useful property is that if we imagine a set of ZAMOs stringing a loop of optical fiber (along which photons travel at the speed of light) around the hole at fixed r, θ, then the pro- and retrograde circuit times are equal. For this reason, ZAMOs are sometimes also called locally non-rotating observers.

η={{-1,0,0,0},{0,1,0,0},{0,0,1,0},{0,0,0,1}};td={{4,1,0,0},{1,4/3,0,0},{0,0,4/3,0},{0,0,0,4/3}};td.η

Transformation to ZAMO

We transform the components using, 
(u^α) = (Subscript[(L^α), (OverHat[β])])(u^(OverHat[β])),
(u^(OverHat[α])) = (Subscript[(L^(OverHat[α])), β])(u^β), 
(Subscript[u, (OverHat[α])]) = (Subscript[u, β])(Subscript[(L^β), (OverHat[α])]),
(Subscript[u, α]) = (Subscript[u, (OverHat[β])])(Subscript[(L^(OverHat[β])), α]),
 (Subscript[(L^α), (OverHat[γ])])(Subscript[(L^(OverHat[γ])), β]) = (Subscript[(δ^α), β])
If the velocity of the ZAMO in B-L coordinates is (u^α), then the timelike orthonormal basis vector in the ZAMO rest frame satisfies (Subscript[(L^α), (OverHat[t])]) = ((((Subscript[e, (OverHat[t])]) ))^α)=(u^α)and we can construct the remaining basis vectors by imposing orthonormality.

bl2z={{Σ/(ρ  (Δ^(1/2))),0,0,0},{0,(Δ^(1/2))/ρ,0,0},{0,0,1/ρ,0},{2 a r/(Σ  ρ (Δ^(1/2))),0,0,ρ/(Σ  Sin[θ])}};
z2bl=FullSimplify[Inverse[bl2z]];
Print[MatrixForm[bl2z]]
Print[MatrixForm[z2bl]]
FullSimplify[bl2z.z2bl]

The sanity check comprises computing the components of the velocity in the ZAMO rest frame.

uzcov=FullSimplify[Table[Sum[ublcov[[β]]bl2z[[β,α]],{β,1,4}],{α,1,4}]]
uzcon=FullSimplify[Table[Sum[z2bl[[α,β]]ublcon[[β]],{β,1,4}],{α,1,4}]]

3+1 Split

It is possible to rewrite the metric in the form
(ds^2)= -(λ ((dt))^2)+ (Subscript[g, jk])((dx^j)+ (β^j)dt)((dx^k)+(β^k)dt)
where the lapse λ = dτ / dt = ρ (Δ^(1/2))/Σ; (x^j)= (r, θ, ϕ) and the shift (β^j)= (0, 0, -(Subscript[Ω, Z])). The spatial metric is diagonal with (Subscript[g, rr])=(ρ^2)/Δ, (Subscript[g, θθ])= (ρ^2), (Subscript[g, ϕϕ])= (Σ sinθ ((/ρ))^2). (Subscript[g, jk])is the metric of a surface of constant t inhabited by ZAMOs, moving on worldlines that are orthogonal to it.  This is often a convenient way to develop physics so as to make look as much as possible like flat space.

Gravity

ZAMOs have to accelerate and their acceleration can be shown to be a = du / dτ = -∇ln λ along (Subscript[e, (OverHat[r])]). This diverges as the 
horizon is approached. A finite surface gravity can be defined for the horizon by taking the limit of g=du/dt = λ a as r->(Subscript[r, H]).
Calculating this gives (Subscript[g, H])= (1/2)-(1/(2(Subscript[r, H]))).

gh[a_]:=(1/(2 ((1/(Sqrt[1-(a^2)]))+1)))
Plot[gh[a],{a,0,0.99999},Frame->True,FrameLabel->{"a","\!\(\*SubscriptBox[\(g\), \(H\)]\)"},RotateLabel->False]

Mass Formulae

Returning to our  formula for the irreducible mass, we can write down a "first law" with clear interpretation
dm = ((((Subscript[g, H]) )/(8π))) d (Subscript[A, H])+ (Subscript[Ω, H])d S
Note that reducing the mass of the hole requires reducing the spin at constant area. 
We can also derive the Smarr formula (using Euler's theorem on homogeneous functions)
m = 2[(((Subscript[g, H])/(8π)))(Subscript[A, H])+ (Subscript[Ω, H])S]

Kerr-Schild Coordinates

Metric

The natural generalization of infalling Eddington-Finkelstein coordinates, in which (Subscript[g, rr])=0, is Kerr coordinates defined by the transformation
t->v = t+r+ log Δ+(1/(((1-(a^2)))^(1/2))) log [((r-1-(((1-(a^2)))^(1/2)))/(r-1+(((1-(a^2)))^(1/2))))];   ϕ->ψ = ϕ+(a/(2(((1-(a^2)))^(1/2))))log [((r-1-(((1-(a^2)))^(1/2)))/(r-1+(((1-(a^2)))^(1/2))))]	
so that 
dt = dv-(((r^2)+(a^2))/Δ)dr;  dϕ = dψ -(a/Δ)dr
The metric tensor is easily computed.

kcoord={v,r,θ,ψ};
bl2k={{1,0,0,0},{-((r^2)+(a^2))/Δ,1,0,-a/Δ},{0,0,1,0},{0,0,0,1}};
kmetric=Simplify[Table[Sum[bl2k[[α,γ]]bl2k[[β,δ]]metric[[γ,δ]],{γ,1,4},{δ,1,4}],{α,1,4},{β,1,4}]];
Print[MatrixForm[kmetric]]

In order to make the light cone be oriented conventionally far from the hole, we change to Kerr-Schild coordinates, where we repeat the same simple transformation (OverTilde[v]) = v-r and suppress the tildes

kscoord={v,r,θ,ψ};
bl2ks={{1,0,0,0},{-2r/Δ,1,0,-a/Δ},{0,0,1,0},{0,0,0,1}};
ksmetric=FullSimplify[Table[Sum[bl2ks[[α,γ]]bl2ks[[β,δ]]metric[[γ,δ]],{γ,1,4},{δ,1,4}],{α,1,4},{β,1,4}]];
Print[MatrixForm[ksmetric]];

Note that (Subscript[g, rr]) no longer vanishes and that we have not introduced any off diagonal elements involving θ.

Ricci Tensor

As a sanity check, compute the Ricci tensor.

inverseksmetric=FullSimplify[Inverse[ksmetric]];
ksaffine=FullSimplify[Table[Sum[(inverseksmetric[[α,ϵ]])
(D[ksmetric[[ϵ,β]],kscoord[[γ]]]+D[ksmetric[[ϵ,γ]],kscoord[[β]]]-D[ksmetric[[β,γ]],kscoord[[ϵ]]]),{ϵ,1,4}]/2,{α,1,4},{β,1,4},{γ,1,4}]];
Print[MatrixForm[ksaffine]];
ksriemann=FullSimplify[Table[D[ksaffine[[α,β,δ]],kscoord[[γ]]]-D[ksaffine[[α,β,γ]],kscoord[[δ]]]+Sum[ksaffine[[α,ϵ,γ]]ksaffine[[ϵ,β,δ]]-ksaffine[[α,ϵ,δ]]ksaffine[[ϵ,β,γ]],{ϵ,1,4}],{α,1,4},{β,1,4},{γ,1,4},{δ,1,4}]];
Print[MatrixForm[ksriemann]];
ksricci=Simplify[Table[Sum[ksriemann[[γ,α,γ,β]],{γ,1,4}],{α,1,4},{β,1,4}]];
Print[MatrixForm[ksricci]];

Homework 2 (due before class Lecture 6)

1. Meaning of Coordinates.

In one concession to thought experiments, imagine that you command a platoon of trained and personally brave surveyors, equipped with clocks, synchronized far from a hole (with specific angular momentum a = 0.95), short strands of optical fiber and maneuverable vehicles. 
a. Explain how  you could position them at rest on a  grid of colatitude and longitude at the same r =6, beyond the static limit.
b. Compute the time it takes light to complete a polar orbit
c. Plot the time it takes light to traverse an orbit with fixed colatitude as a function of θ. 
e. Explain how you could get them to synchronize their clocks among themselves and with you, (who have wisely 
chosen to remain far from the horizon).
f. After your time has advanced by 10, resynchronise. Plot the associated change in proper time of your associates as a function of θ
g. Now have a very few of your associates venture inside the ergosphere to r = 1.5, θ = π/2. Have them become ZAMO
observers by calculating what their proper orbital periods should be. Compute the time it will take light to complete an equatorial orbit and return to them and show that it is the same for forward and backward directed orbits

2. Black Hole Thermodynamics

a. Use the area theorem to verify that ((((m=(Subscript[m, i ])(1-β)^2)))^(-1/2)) where β = 2(Subscript[Ω, H])(Subscript[m, i]).
b. Plot the reducible mass m - (Subscript[m, i]) as a function of a and (Subscript[Ω, H]) on the same graph for m = 1.
c. Verify the computation of surface gravity and the Smarr formula.

3. Particle Orbit.

Compute the trajectory of a particle falling into a spinning hole with (a=0.9) in Kerr-Schild coordinates assuming it starts from infinity at rest with zero angular momentum.

4. Volume of Ergosphere

Compute the 3 volume of the ergosphere as a function of a for a surface t=constant.

3. Orbits and Rays

Motivation

Stellar orbits

Stars are observed in orbit around the 4 million solar mass hole at our Galactic center (Sgr A*). 
http://en.wikipedia.org/wiki/Sagittarius_A*
http://www.youtube.com/watch?v=EvuV3GdVaY4
They approach as close as 60 AU ~ 1500m and follow Keplerian orbits. There is every reason to believe (and some evidence to show from X-ray flares) that most massive black holes are orbited by star clusters extending in close to the horizon. What is remarkable about the Sgr A* cluster is that the stars are young. They were probably formed in an accretion disk.

Particle Orbits

BL Coordinates

Clear[coord,metric,inversemetric,affine];
coord={t,r,θ,ϕ};
metric={{-(((r^2)-2r+(a^2)((Cos[θ])^2))/((r^2)+(a^2)((Cos[θ])^2))),0,0,((-2a r ((Sin[θ])^2))/((r^2)+(a^2)((Cos[θ])^2)))},
{0,(((r^2)+(a^2)((Cos[θ])^2))/((r^2)-2r+(a^2))), 0,0},{0,0,(r^2)+(a^2)((Cos[θ])^2),0},{((-2a r ((Sin[θ])^2))/((r^2)+(a^2)((Cos[θ])^2))),0,0,((((r^2) +(a^2))((r^2)+(a^2)((Cos[θ])^2))+2r (a^2)((Sin[θ])^2))/((r^2)+(a^2)((Cos[θ])^2)))((Sin[θ])^2)}};
ρ=(Sqrt[(r^2)+(a^2)((Cos[θ])^2)]);
Δ=(r^2)-2r+(a^2);
Σ = (Sqrt[((((r^2)+(a^2)))^2)-Δ (a^2)((Sin[θ])^2)]);
inversemetric=FullSimplify[Inverse[metric]];
affine=FullSimplify[Table[Sum[(inversemetric[[α,δ]])
(D[metric[[δ,β]],coord[[γ]]]+D[metric[[δ,γ]],coord[[β]]]-D[metric[[β,γ]],coord[[δ]]]),{δ,1,4}]/2,{α,1,4},{β,1,4},{γ,1,4}]];

Horizon and Ergosphere

rh[a_]:=1+(Sqrt[1-(a^2)])
rergo[θ_,a_]:=1+(Sqrt[1-(a^2)((Cos[θ])^2)])

Killing Vectors

Killing's equation assures us that
(Subscript[ξ, α;β])+(Subscript[ξ, β;α])= ((Subscript[g, αδ])(Subscript[(Γ^δ), βγ])+(Subscript[g, βδ])(Subscript[(Γ^δ), αγ]))(ξ^γ) = 0

in BL coordinates. ξ·u is conserved along a geodesic. Verifying for the temporal and azimuthal Killing vectors,

γ=4;
Print[MatrixForm[FullSimplify[Table[Sum[metric[[α,δ]]affine[[δ,β,γ]]+metric[[β,δ]]affine[[δ,α,γ]],{δ,1,4}],{α,1,4},{β,1,4}]]]];

Killing Tensor

There is also a Killing tensor that we construct from the ingoing and outgoing principal null vectors
(l^α) = [ (r^2)+ (a^2), Δ, 0, a ] / Δ
(n^α)= [ (r^2)+ (a^2), -Δ, 0, a ] / 2(ρ^2)
where l.l = n.n = 0; l.n = -1. These geodesics are shear free.

lnull={((r^2)+(a^2))/Δ,1,0,a/Δ};
nnull={((r^2)+(a^2))/(2(ρ^2)),-Δ/(2(ρ^2)),0,a/(2(ρ^2))};
FullSimplify[lnull.metric.lnull]
FullSimplify[nnull.metric.nnull]
FullSimplify[lnull.metric.nnull]

The covariant form of the Killing tensor is given by
(Subscript[ξ, αβ]) = 2 (Subscript[l, (α])(Subscript[n, β)])(ρ^2) + (r^2)(Subscript[g, αβ])

lnullcov=FullSimplify[metric.lnull];
nnullcov=FullSimplify[metric.nnull];
ktcov=FullSimplify[Table[(lnullcov[[α]]nnullcov[[β]]+nnullcov[[α]]lnullcov[[β]])(ρ^2)
+(r^2)metric[[α,β]],{α,1,4},{β,1,4}]];
Print[MatrixForm[ktcov]];

Killing's equation generalizes to 
(Subscript[ξ, (αβ;γ)])= 0

term1=FullSimplify[Table[D[ktcov[[α,β]],coord[[γ]]]-
2 Sum[affine[[δ,α,γ]]ktcov[[δ,β]],{δ,1,4}],
{α,1,4},{β,1,4},{γ,1,4}]];
term2=FullSimplify[Table[D[ktcov[[α,β]],coord[[γ]]]-
2 Sum[affine[[δ,α,γ]]ktcov[[δ,β]],{δ,1,4}],
{γ,1,4},{α,1,4},{β,1,4}]];
term3=FullSimplify[Table[D[ktcov[[α,β]],coord[[γ]]]-
2 Sum[affine[[δ,α,γ]]ktcov[[δ,β]],{δ,1,4}],
{β,1,4},{γ,1,4},{α,1,4}]];
FullSimplify[term1+term2+term3]

Carter Constant

This provides us with a constant of the motion, 
K= (Subscript[ξ, αβ]) (u^α)(u^β) = 2 (l^α)(Subscript[u, α])(n^β)(Subscript[u, β])(ρ^2)-(r^2)
where
(Subscript[u, α])= [ -E, (Subscript[u, r]), (Subscript[u, θ]), L]

ucov={-e,(Subscript[u, r]),(Subscript[u, θ]),l};

fac1=FullSimplify[Sum[lnull[[α]]ucov[[α]],{α,1,4}]]
fac2=FullSimplify[Sum[nnull[[α]]ucov[[α]],{α,1,4}]]

We use the invariant length of the velocity to eliminate (Subscript[u, r]), choosing the positive value

ur=FullSimplify[Solve[ucov.inversemetric.ucov==-1,(Subscript[u, r])][[2]]]
fac3=FullSimplify[fac1/.ur];
fac4=FullSimplify[fac2/.ur];
fac5=Simplify[fac3 fac4];

K is then given by

int3=FullSimplify[2fac5 (ρ^2)-(r^2)]

This is good for both signs of (Subscript[u, r]). It turns out to be more convenient to use the "Carter constant"
Q = K-(L - a E )^(2 )= ((Subscript[u, θ])^2)+ (a^2)(1 - (E^2)) (cos^2)θ +(L^2)(cot^2)θ
Q provides an energy equation for vertical epicyclic motion about the equatorial plane. An alternative approach follows Hamilton Jacobi-equations.

Equations of Motion

We can now solve for the equations of motion directly, without having to appeal to Lagrange's equations. Taking them in turn,
(dt/(dτ))=(u^t)= (g^(tα))(Subscript[u, α])= ((((r^2)+(a^2)){((r^2)+(a^(2 )))E-a L})/((Δρ)^2))-((a(a E (sin^2)θ -L))/(ρ^2)) 
It turns out to be convenient to replace the proper time with a new parameter x defined by
((dτ)/dx)=(ρ^(2 ))
We end up with 
t'≡(dt/dx)= ((((r^2)+(a^2)){((r^2)+(a^(2 )))E-a L})/Δ)-a(a E (sin^2)θ -L)

ucov={-e,(Subscript[u, r]),(Subscript[u, θ]),l};
FullSimplify[(ρ^2)(inversemetric .ucov)[[1]]];

ft[r_,θ_]:=((((r^2)+(a^2))(((r^2)+(a^2))e-a l))/(((r^2)-2r+(a^2))))-a(a e ((Sin[θ])^2)-l)
fτ[r_,θ_]:=(r^2)+(a^2)((Cos[θ])^2)

Next, turn to the radial equation
((((r')^2)=((((ρ^2)(u^r)))^2)=((ρ^2)(g^rr)(Subscript[u, r])))^2)=R(r)= (({((r^2)+(a^2))E-aL})^2)-Δ{(r^2)+(((aE-L))^2)+Q}

eq1=FullSimplify[ucov.inversemetric.ucov==-1];
eq2=FullSimplify[q==((Subscript[u, θ])^2)+(a^2)(1-(e^2))((Cos[θ])^2) + (l^2) ((Cot[θ])^2)];
FullSimplify[(ρ^4)((inversemetric[[2,2]])^2)((Subscript[u, r])^2)/.FullSimplify[Solve[{eq1,eq2},(Subscript[u, r]),(Subscript[u, θ])]][[2]]];

fr[r_]:=(((((r^2)+(a^2))e-a l))^2)-((r^2)-2r+(a^2))((r^2)+(((a e-l))^2)+q)

The radial equation has turning points and it is simpler, though slower, to differentiate and deal with the second order equation
r''= R'(r)/2 = 2r E {((r^2)+(a^2))E-aL} - (r-1){(r^2)+(((aE-L))^2)+Q} - r Δ

Simplify[D[fr[r],r]/2]

dfr[r_]:=2r e(((r^2)+(a^2))e-a l)-(r-1)((r^2)+(((a e-l))^2)+q)-r((r^2)-2r+(a^2))

Treat the θ equation similarly.
((((θ')^2)=((((ρ^2)(u^θ)))^2)=Θ(θ)=((ρ^2)(g^(θθ))(Subscript[u, θ])))^2)= Q - (a^2)(1-(E^2))(cos^2)θ -(L^2)(cot^2)θ

FullSimplify[(ρ^4)((inversemetric[[3,3]])^2)((Subscript[u, θ])^2)/.FullSimplify[Solve[eq2,(Subscript[u, θ])][[2]]]];

fθ[θ_]:=q-(a^2)(1-(e^2))((Cos[θ])^2)-(l^2)((Cot[θ])^(2 ));

Differentiating,
θ''' = (a^2) (1-(E^2))cosθ sinθ +(L^2)cotθ (csc^2)θ

Simplify[D[fθ[θ],θ]/2]

dfθ[θ_]:=(a^2)(1-(e^2))Sin[θ]Cos[θ]+(l^2)Cot[θ]((Csc[θ])^2)

The azimuthal equation becomes 
ϕ' = (ρ^2)(g^(ϕα))(Subscript[u, α])= ((a{((r^2)+(a^2))E-aL})/Δ)-aE + L (csc^2)θ

fϕ[r_,θ_]:=((a(((r^2)+(a^2))e-a l))/(((r^2)-2r+(a^2))))-a e+l ((Csc[θ])^2);

Numerical Integration

Although some of these equations are expressible in terms of elliptic functions and integrals, we will simply use a brute force orbit integrator.

orbit:=NDSolve[{
t'[x]==ft[r[x],θ[x]],
r''[x]==dfr[r[x]],
θ''[x]==dfθ[θ[x]],
ϕ'[x]==fϕ[r[x],θ[x]],
τ'[x]==fτ[r[x],θ[x]],
t[0]==0,
r[0]==(Subscript[r, 0]),
r'[0]==Sqrt[fr[(Subscript[r, 0])]],
θ[0]==(Subscript[θ, 0]),
θ'[0]==Sqrt[fθ[(Subscript[θ, 0])]],
ϕ[0]==0,
τ[0]==0
},{t,r,θ,ϕ,τ},{x,0,xend }]

Initial Conditions, parameters and orbit calculation

We set the initial conditions for the orbit and check that they are consistent

a=ap=.5;
(Subscript[r, 0])=10.;
(Subscript[θ, 0])=1.6;
e=ep=.98;
l=lp=4.22;
q=qp=1;
Sqrt[fr[(Subscript[r, 0])]];
Sqrt[fθ[(Subscript[θ, 0])]];
xend=3.;
run1=orbit;
Clear[a,e,l,q];

Orbit Plots

orbitplt:=Show[GraphicsGrid[{{ParametricPlot[Evaluate[{τ[x],t[x]}/.run1],{x,0,xend},Frame->True,AspectRatio->1,PlotRange->All,FrameLabel->{"τ","t",ap},DisplayFunction->Identity],ParametricPlot[Evaluate[{τ[x],r[x]}/.run1],{x,0,xend},Frame->True,AspectRatio->1,PlotRange->All,FrameLabel->{"τ","r",ep},DisplayFunction->Identity]},{ParametricPlot[Evaluate[{τ[x],θ[x]}/.run1],{x,0,xend},Frame->True,AspectRatio->1,PlotRange->All,FrameLabel->{"τ","θ",qp},DisplayFunction->Identity],ParametricPlot[Evaluate[{τ[x],ϕ[x]}/.run1],{x,0,xend},Frame->True,AspectRatio->1,PlotRange->All,FrameLabel->{"τ","ϕ",lp},DisplayFunction->Identity]}}]]
orbitplt

2D Representation

scale=40;
hole2D=Graphics[{PointSize[((rh[ap])/scale)],Point[{0,0}]},AspectRatio->1];
xyplt=Show[ParametricPlot[Evaluate[{r[x] Sin[θ[x]] Cos[ϕ[x]],r[x] Sin[θ[x]] Sin[ϕ[x]]}/.run1],{x,0,xend},Frame->True,FrameLabel->{"x","y"},PlotRange->{{-scale,scale},{-scale,scale}},AspectRatio->Automatic,DisplayFunction->Identity],hole2D,DisplayFunction->Identity];
xzplt=Show[ParametricPlot[Evaluate[{r[x] Sin[θ[x]] Cos[ϕ[x]],r[x] Cos[θ[x]]}/.run1],{x,0,xend},Frame->True,FrameLabel->{"x","z"},PlotRange->{{-scale,scale},{-scale,scale}},AspectRatio->Automatic,DisplayFunction->Identity],hole2D,DisplayFunction->Identity];
yzplt=Show[ParametricPlot[Evaluate[{r[x] Sin[θ[x]] Sin[ϕ[x]],r[x] Cos[θ[x]]}/.run1],{x,0,xend},Frame->True,FrameLabel->{"y","z"},PlotRange->{{-scale,scale},{-scale,scale}},AspectRatio->Automatic,DisplayFunction->Identity],hole2D,DisplayFunction->Identity];
Show[GraphicsGrid[{{xyplt,xzplt},{yzplt}}]]

3D Representation

hole3D=Graphics3D[{PointSize[((rh[a])/40)],Point[{0,0,0}]},AspectRatio->Automatic]
xyzplt=Show[ParametricPlot3D[Evaluate[{r[x] Sin[θ[x]] Cos[ϕ[x]],r[x] Sin[θ[x]] Sin[ϕ[x]],r[x] Cos[θ[x]]}/.run1],{x,0,3},Boxed->False,AspectRatio->Automatic,PlotRange->{{-40,40},{-40,40},{-40,40}},DisplayFunction->Identity],hole3D,DisplayFunction->$DisplayFunction]

Circular Equatorial Orbits

Energy

We require that r'=r''=0 and note that Q=0. Solving for E and choosing the relevant root, adopting the convention that a>0 is a prograde orbit, we obtain

FullSimplify[Solve[{fr[(x^2)]==0,dfr[(x^2)]==0}/.q->0,e,l][[1]]]/.x->(r^(1/2))

We find
E=((1-2(r^(-1))+a (r^(-3/2)))/(((1-3(r^(-1))+2a (r^(-3/2))))^(1/2)))for Keplerian orbits. Specializing to a=0, E=((1-2(r^(-1)))/(((1-3(r^(-1))))^(1/2))).

ek[r_,a_]:=((1-2(r^(-1))+a (r^(-3/2)))/(((1-3(r^(-1))+2a (r^(-3/2))))^(1/2)))

Plot up the energy for different spins. It comprises the rest mass minus the binding energy.

a=.99;
Plot[ek[r,a],{r,rh[a],20},Frame->True,PlotRange->{{0,20},{0,2}},FrameLabel->{"r","(Subscript[E, K])"},RotateLabel->False]
Clear[a];

Photon Orbit

The energy becomes infinite at the photon orbit where the denominator of E->0 or a=(3/2)(r^(1/2))-(1/2)(r^(3/2)) .

rph[ap_]:=r/.NSolve[3(r^(1/2))-(r^(3/2))-2a==0,r][[1]]
Plot[rph[a],{a,-1,1},Frame->True,FrameLabel->{"a","(Subscript[r, ph])"},RotateLabel->False]

Marginally Bound Orbit

The marginally bound circular orbit occurs when E=1

FullSimplify[Solve[(((1-2(r^(-1))+a (r^(-3/2))))^2)==1-3(r^(-1))+2a (r^(-3/2)),r][[2]]]

or (Subscript[r, mb])= 2-a + 2(((1-a))^(1/2))

rmb[a_]:=2-a+2(((1-a))^(1/2));
Plot[{rmb[a],rph[a]},{a,-1,1},Frame->True,FrameLabel->{"a","(Subscript[r, mb])"},RotateLabel->False]

Angular Momentum

Similarly, solving for L

FullSimplify[Solve[{fr[(x^2)]==0,dfr[(x^2)]==0}/.q->0,l,e][[2]]]/.x->(r^(1/2))

We find
L=(((1-2a (r^(-3/2))+(a^2)(r^(-2))) (r^(1/2)))/(((1-3(r^(-1))+2a (r^(-3/2))))^(1/2)))

lk[r_,a_]:=(((1-2a (r^(-3/2))+(a^2)(r^(-2))) (r^(1/2)))/(((1-3(r^(-1))+2a (r^(-3/2))))^(1/2)))

Stability

Look at the "effective energy" equation effective potential, ((r')^2)-R(r)=0. The particle oscillates between turning points on the x-axis.

a=.99;q=0.;rk=4.;
e=ek[rk,a]+.001
l=lk[rk,a]
Plot[-fr[r],{r,1,10},Frame->True,PlotRange->{{0,10},{-.5,.5}},FrameLabel->{r,"-R(r;a,e,l)",a,""},RotateLabel->False]
Clear[a,e,l,q,r];

Stability requires that the radial potential have positive curvature or 
R''(r) < 0

FullSimplify[D[2dfr[r]/.q->0,r]/.{e->ek[r,a],l->lk[r,a]}]

FullSimplify[Solve[(((r(r-6)-3 (a^2)))^2)==64(a^2)r,r][[3]]]

which leads to 
(Subscript[r, ms])=3∓(Sqrt[3+(a^2)+(-3+a) (((-1+a))^(1/3)) (((1+a))^(2/3))+(((-1+a))^(2/3)) (((1+a))^(1/3)) (3+a)])+√(18-(-3+a) (((-1+a))^(1/3)) (((1+a))^(2/3))-(((-1+a))^(2/3)) (((1+a))^(1/3)) (3+a)+2 (-6+(a^2))∓((16 (a^2))/(Sqrt[3+(a^2)+(-3+a) (((-1+a))^(1/3)) (((1+a))^(2/3))+(((-1+a))^(2/3)) (((1+a))^(1/3)) (3+a)])))
after choosing the relevant signs.

FullSimplify[D[dfr[r]/.q->0,r]/.{e->ekep[r],l->lkep[r]}];
FullSimplify[Solve[(((r(r-6)-3 (a^2)))^2)==64(a^2)r,r]];

rms[a_]:=Chop[
3-Sign[a]√(3+(a^2)+(-3+a) (((-1+a))^(1/3)) (((1+a))^(2/3))+(((-1+a))^(2/3)) (((1+a))^(1/3)) (3+a))+√(18-(-3+a) (((-1+a))^(1/3)) (((1+a))^(2/3))-(((-1+a))^(2/3)) (((1+a))^(1/3)) (3+a)+2 (-6+(a^2))-Sign[a](16 (a^2))/(√(3+(a^2)+(-3+a) (((-1+a))^(1/3)) (((1+a))^(2/3))+(((-1+a))^(2/3)) (((1+a))^(1/3)) (3+a))))]

Plot of Equatorial Radii

Let us now plot these various radii on a single graph.

rhplt=Plot[rh[a],{a,-1,1},Frame->True,PlotRange->{{-1,1},{0,9}},FrameLabel->{"a","r"}];
rergoplt=Plot[rergo[(π/2),a],{a,-1,1}];
rphplt=Plot[rph[a],{a,-1,1}];
rmbplt=Plot[rmb[a],{a,-1,1}];
rmsplt=Plot[rms[a],{a,-1,1}];
Show[rhplt,rergoplt,rphplt,rmbplt,rmsplt,DisplayFunction->$DisplayFunction]

Energy and Angular Momentum Plots

The energy of a circular orbit of radius r with a range of values of a.

ek1plt=Plot[ek[r,1],{r,rms[1],10},Frame->True,PlotRange->{{0,10},{0.5,1}},FrameLabel->{"r","E"},RotateLabel->False];
ek2plt=Plot[ek[r,.99],{r,rms[.99],10}];
ek3plt=Plot[ek[r,0.],{r,rms[0.],10}];
ek4plt=Plot[ek[r,-1.],{r,rms[-1],10}];
Show[ek1plt,ek2plt,ek3plt,ek4plt,DisplayFunction->$DisplayFunction]

It will be useful to plot the binding energy of the marginally stable circular orbit (ISCO).

beplt1=Plot[1-ek[rms[a],a],{a,-1,1},Frame->True,PlotRange->{{-1,1},{0.,0.5}},FrameLabel->{"a","1-E"},RotateLabel->False]

Note that it is 1-(3^(-3/2))5=0.038 for a=-1,1-(2^(3/2))/3=0.057 for a=0 and 1-(3^(-1/2))=0.42 for a=1.

The corresponding angular momenta, shown as positive functions.

lk1plt=Plot[lk[r,1],{r,rms[1],9},Frame->True,PlotRange->{{0,10},{0.,5}},FrameLabel->{"r","L"},RotateLabel->False];
lk2plt=Plot[lk[r,0.99],{r,rms[0.99],10}];
lk3plt=Plot[lk[r,0.],{r,rms[0.],10}];
lk4plt=Plot[lk[r,-1.],{r,rms[-1.],10}];
Show[lk1plt,lk2plt,lk3plt,lk4plt,DisplayFunction->$DisplayFunction]

Angular Frequency

We next turn to the two equations of motion, reverting to τ as the independent variable
(dt/(dτ))=(((r^3) (-2 a L+E (r^3)+(a^2) E (2+r)))/((a^2)+(-2+r) r))

((dϕ)/(dτ))=(((2 a E+L (-2+r)) (r^3))/((a^2)+(-2+r) r))
from which we obtain 
Ω = (1/((r^(3/2))+a))

ut=FullSimplify[(r^2)ft[r,π/2]];
uϕ=FullSimplify[(r^2)fϕ[r,π/2]];
FullSimplify[uϕ/ut/.{e->ek[r,a],l->lk[r,a]}]

Evaluating

Ωk[r_,a_]:=((((r^(3/2))+a))^(-1))

A useful relation for circular geodesics

The acceleration vanishes. Therefore 
(Subscript[a, α])=(u^β)(Subscript[u, α;β]) = (u^β)((Subscript[u, α;β])-(Subscript[u, β;α])) = (u^β)((Subscript[u, α,β])-(Subscript[u, β,α])) = 0 
Hence 
(Subscript[a, r])=-(u^t)(Subscript[u, t,r]) - (u^ϕ)(Subscript[u, ϕ,r])=0
As we change the radius of the orbit, 
(u^t)(Subscript[du, t])+(u^ϕ)(Subscript[du, ϕ])=0
and so
(dE/dL)=Ω
Another way to see this is to note that, if there are no other forces involved, the work done by a tangential force equals the change in the energy. 
To verify,

FullSimplify[D[ek[r,a],r]-Ωk[r,a]D[lk[r,a],r]]

Rays

Equations of Motion

The equations for null geodesics can be obtained simply from the those for timelike geodesics by taking the limit E, L ->∞ (recall that these are quantities per unit rest mass m)
and choosing a suitable affine parameter, essentially the limit of dτ/m. As we will mainly be interested in what a distant observer along direction (Subscript[θ, 0]), might see, we introduce an impact vector b with components (Subscript[b, 1]), (Subscript[b, 2]) resolved perpendicular and parallel to the projected hole spin axis. We let Ex->x, L/E≡l->-(Subscript[b, 1])(Subscript[sinθ, 0]), (Subscript[u, θ])/E->(Subscript[b, 2]) at large distance and (Q/(E^2))->(q^2)= (Subscript[b, 2]^2)+(cos^2)(Subscript[θ, 0])((Subscript[b, 1]^2)-(a^2)).
The orbits just depend upon l and q which are, in turn, fixed by b. As in the Newtonian case we substitute u=1/r and use look up tables. The equations become:
t'= ((((r^2)+(a^2)){((r^2)+(a^(2 )))-a l})/Δ)-a(a (sin^2)θ -l) 
((r')^2)=(({((r^2)+(a^2))-a l})^2)-Δ{(((a-l))^2)+(q^2)}
r''= 2r{((r^2)+(a^2))-a l}-(r-1){(((a-l))^2)+(q^2)} 
((θ')^2)=(q^2)+(a^2)(cos^2)θ -(l^2)(cot^2)θ
θ'' = -(a^2) cosθ sinθ +(l^2)cotθ (csc^2)θ 
ϕ'= a{((((r^2)+(a^2))-a l)/Δ)}+(l(csc^2)θ-a)

gt[u_,θ_]:=((((r^2)+(a^2))(((r^2)+(a^2))-a l)-a((r^2)-2r+(a^2))(a  ((Sin[θ])^2)-l))/(((r^2)-2r+(a^2))))/.r->(1/u)
gu[u_]:=((((((r^2)+(a^2))-a l))^2)-((r^2)-2r+(a^2))((((a -l))^2)+(q^2)))/.r->(1/u)
dgu[u_]:=2r (((r^2)+(a^2))-a l)-(r-1)((((a -l))^2)+(q^2))/.r->(1/u)
gθ[θ_]:=(q^2)+(a^2)((Cos[θ])^2)-(l^2)((Cot[θ])^(2 ))

dgθ[θ_]:=-(a^2)Sin[θ]Cos[θ]+(l^2)Cot[θ]((Csc[θ])^2)
gϕ[u_,θ_]:=((a(((r^2)+(a^2))-a l)+((r^2)-2r+(a^2))(l ((Csc[θ])^2)-a))/(((r^2)-2r+(a^2))))/.r->(1/u)

Numerical Integration

We integrate backwards along a ray from the observer plane located at a radius = ((Subscript[u, 0])^(-1)).

ray:=NDSolve[{
t'[x]==gt[u[x],θ[x]],
u''[x]-(2/(u[x]))((u'[x])^2)==-((u[x])^2)dgu[u[x]],
θ''[x]==dgθ[θ[x]],
ϕ'[x]==gϕ[u[x],θ[x]],
t[0]==0,
u[0]==(Subscript[u, 0]),
u'[0]==-((Subscript[u, 0])^2)Sqrt[gu[(Subscript[u, 0])]],
θ[0]==(Subscript[θ, 0]),
θ'[0]==Sqrt[gθ[(Subscript[θ, 0])]],
ϕ[0]==0
},
{t,u,θ,ϕ},
{x,0,xstart }]

Initial Conditions and Parameters

a=.5;
(Subscript[u, 0])=.01;
(Subscript[θ, 0])=1.;
b={5,5};
l=-b[[1]]Sin[(Subscript[θ, 0])]
q=Sqrt[((b[[2]])^2)+((Cos[(Subscript[θ, 0])])^2)(((b[[1]])^2)-(a^2))]
Sqrt[gu[(Subscript[u, 0])]]
Sqrt[gθ[(Subscript[θ, 0])]]
xstart=-.62

Ray Calculation

run2=ray

Ray Plots

rayplt:=Show[GraphicsGrid[{{Plot[Evaluate[t[x]/.run2],{x,0,xstart},Frame->True,PlotRange->All,FrameLabel->{"x","t"},DisplayFunction->Identity],Plot[Evaluate[(1/(u[x]))/.run2],{x,0,xstart},Frame->True,PlotRange->All,FrameLabel->{"x","r"},DisplayFunction->Identity]},{Plot[Evaluate[θ[x]/.run2],{x,0,xstart},Frame->True,PlotRange->All,FrameLabel->{"x","θ"},DisplayFunction->Identity],Plot[Evaluate[ϕ[x]/.run2],{x,0,xstart},Frame->True,PlotRange->All,FrameLabel->{"x","ϕ"},DisplayFunction->Identity]}}]]
rayplt

2D Representation

scale=100;
rh[a_]:=1+(Sqrt[1-(a^2)])
xyrayplt=Show[ParametricPlot[Evaluate[{Sin[θ[x]] Cos[((ϕ[x])/(u[x]))],((Sin[θ[x]] Sin[ϕ[x]])/(u[x]))}/.run2],{x,0,xstart},Frame->True,FrameLabel->{"x","y"},PlotRange->{{-scale,scale},{-scale,scale}},AspectRatio->Automatic,DisplayFunction->Identity],hole2D,DisplayFunction->Identity];
xzrayplt=Show[ParametricPlot[Evaluate[{((Sin[θ[x]] Cos[ϕ[x]])/(u[x])),((Cos[θ[x]])/(u[x]))}/.run2],{x,0,xstart},Frame->True,FrameLabel->{"x","z"},PlotRange->{{-scale,scale},{-scale,scale}},AspectRatio->Automatic,DisplayFunction->Identity],hole2D,DisplayFunction->Identity];
yzrayplt=Show[ParametricPlot[Evaluate[{((Sin[θ[x]] Sin[ϕ[x]])/(u[x])),((Cos[θ[x]])/(u[x]))}/.run2],{x,0,xstart},Frame->True,FrameLabel->{"y","z"},PlotRange->{{-scale,scale},{-scale,scale}},AspectRatio->Automatic,DisplayFunction->Identity],hole2D,DisplayFunction->Identity];
Show[GraphicsGrid[{{xyrayplt,xzrayplt},{yzrayplt}}]]

3D Representation

hole3D=Graphics3D[{PointSize[((rh[a])/scale)],Point[{0,0,0}]},AspectRatio->Automatic];
xyzrayplt=Show[ParametricPlot3D[Evaluate[{((Sin[θ[x]] Cos[ϕ[x]])/(u[x])),((Sin[θ[x]] Sin[ϕ[x]])/(u[x])),((Cos[θ[x]])/(u[x]))}/.run2],{x,0,xstart},Boxed->False,AspectRatio->Automatic,PlotRange->{{-scale,scale},{-scale,scale},{-scale,scale}},AxesLabel->{"x","y","z"},DisplayFunction->Identity],hole3D,DisplayFunction->$DisplayFunction]

Penrose Process

Negative Energy States

The key to this process is the realization that (Subscript[ξ, (t)])≡(∂_t) becomes spacelike inside the ergosphere as (Subscript[g, tt])= (Subscript[ξ, (t)]). (Subscript[ξ, (t)])< 0.
This means that the energy E=-p.(Subscript[ξ, (t)]) can be negative as it is now a component of 3-momentum.  Orbits will be most negative with vanishing rest mass, ie photons. Wlog set the absolute value of the angular momentum as unity. Also, keep things simple by restricting to the equatorial plane. We examine the orbits at maximum radius when there is no radial motion. There are two orbits, both forwardly progressive within the ergosphere one with angular momentum parallel to the spin of the hole, the other anti-parallel.  The interval vanishes.

energyplus[rmax_,a1_]:=e/.FullSimplify[Solve[{-e,0,0,1}.inversemetric.{-e,0,0,1}==0/.θ->(π/2),e]][[2]]/.{r->rmax,a->a1}
energyminus[rmax_,a1_]:=e/.FullSimplify[Solve[{-e,0,0,-1}.inversemetric.{-e,0,0,-1}==0/.θ->(π/2),e]][[2]]/.{r->rmax,a->a1}
Plot[{energyplus[rmax,1],energyminus[rmax,1]},{rmax,rh[1],rergo[(π/2),1]},Frame->True,FrameLabel->{"(Subscript[r, max])","E/|L|"},RotateLabel->False]

As (Subscript[r, max]) approaches the horizon, E ->(Subscript[Ω, H])L

So there exist isolated negative energy orbits in the ergosphere.

Energy Extraction

Suppose that the particle has negative radial motion so that it approaches the horizon (in BL) and actually crosses it in KS. It will change the mass by E and the spin by L.  If we recall the mass formula
dm = ((((Subscript[g, H]) )/(8π))) d (Subscript[A, H])+ (Subscript[Ω, H])d S
we see the the area(entropy) is destined not to decrease and the purely retrograde photon hovering just above the horizon corresponds to a reversible change.

Penrose Orbits

Energy extraction is possible in principle because we can scatter a photon (or massive particle) onto a negative energy orbit within the ergosphere. Momentum must be conserved in this process and the remaining outgoing energy must exceed the ingoing energy. However, the phase space for this to happen is small and this is not an important process naturally.

Weak Field Limit

Gravitomagnetism

In the gravitomagnetic approach the field equations and geodesic equation are taken to the weak field, slow motion limit and expressed in the language of electromagnetic theory. A gravitoelectric and gravitomagnetic field are introduced so that
E=-∇ϕ; ϕ = -(1+(Subscript[g, 00]))/2
B = ∇ × A ; (Subscript[A, i]) = (Subscript[g, 0i])
∇·E = -4πρ;  ∇·B = 0; ∇ × E =0; ∇ × B = -4(4πρv-((∂E)/(∂t)))
(dv/dt)= E+v × B
Note the sign changes and absence of a term in the third equation and th presence of a factor 4 in the fourth equation.
Proceeding by direct analogy with electromagnetic theory, the gravitomagnetic field for a body with spin S is
S = ((3(S·r)r-((Sr)^2))/(r^5))
A slightly inclined orbit will then precess at a rate given by integrating the torque around the orbit
(Subscript[Ω, p])=((2S)/(r^3))
as may be verified directly by expanding the equations for θ',ϕ'.
This formalism can also be used to discuss the precession of spinning black holes.  Specifically the gravitomagnetic 
field associated with motion through the gravitoelectric field of a non-spinning mass, - v×E induces a spin precession
at a rate -H / 2 = Ω / 2r. This must be added to the geometrical precession induced by space curvature to give a combined geodetic precession of 3Ω/2r. 
Finally, there is Lense Thirring precession which arises when an orbiting spin precesses in the gravitomagnetic field. This is retrograde and of magnitude -S/2(r^3)when the orbit is equatorial and prograde
with twice the magnitude when the orbit is polar. We will return to these matters.

Homework 3 (due Tuesday May 6)

1. Penrose Orbits

We discussed the nature of Penrose orbits and showed how they could not exist outside the ergosphere. Consider a single, horizon-crossing, negative energy orbit. Where does it come from?

2. Quasi-periodic X-ray Signal from a Star Orbiting a Massiove Black Hole

A phenomenon which has been sought in X-ray observations (and which may well have been detected in recent observations of Sgr A*) is a  burst of hard X-rays produced every time a star crashes through an equatorial accretion disk. Consider a 3 million solar mass, spinning black hole with a=0.9.  Suppose that there is a half solar mass star on a prograde orbit with r = 6, inclined at (30^o) to the equator and crossing it first time at the Eobserver's longitude ϕ. 
a. Determine the BL coordinates of successive crossings of the equatorial plane.
b. Compute rays which connect these events to the observer's sky plane and hence determine the time interval
between successive, observed X-ray pulses.  
c. Perform some sanity checks of your routines.
(If you have difficulty debugging your software, generous partial credit will be given for a good description of what you are trying to do.)

4. Thin Disks

Motivation

X-ray Binary and AGN Accretion Disks

We regularly observe accretion disks around stellar and massive black holes and neutron stars.It appears that AGN holes create radiant energy per unit mass with efficientcy~0.2, suggesting that the holes are spinning rapidly as this exceeds the binding energy of the most bound,stable circular orbits. In addition, broad Fe lines are observed from disks around both stellar and massive holes, suggesting that the disks exist very close to the horizon, again indicating rotation. In fact the power in the line is large that it has been suggested that it must also derive from the ergosphere.

Simplifying Assumptions

Observed disks

Accretion disks are complex. In reality, they are non-stationary, non-axisymmetric, thick, non-planar, inhomogeneous and the (magnetic) torque is non-Newtonian. The emission is nonthermal and anisotropic. The relativity is the best understood feature of their description although numerical simulation is leading to a much more sophisticated description. In order to bring out their relativistic properties I shall suppose that disks are 2D, equatorial, cold, prograde, circular stationary and axisymmetric while the emission is thermal and isotropic. A good general reference is Krolik: Active Galactic Nuclei, Princeton.

Geometry and Kinematics

Metric

The 3D metric is given by:

dcoord={t,r,ϕ};
dmetric={{-(1-(2/r)),0,((-2a)/r)},{0,((r^2)/((r^2)-2r+(a^2))),0},{((-2a)/r),0,(((r^2)((r^2)+(a^2))+2(a^2)r)/(r^2))}};
Print[MatrixForm[dmetric]]
Δ=(r^2)-2r+(a^2);
inversedmetric=FullSimplify[Inverse[dmetric]];
daffine=FullSimplify[Table[Sum[(inversedmetric[[i,s]])
(D[dmetric[[s,j]],dcoord[[k]]]+D[dmetric[[s,k]],dcoord[[j]]]-D[dmetric[[j,k]],dcoord[[s]]]),{s,1,3}]/2,{i,1,3},{j,1,3},{k,1,3}]];

and the determinant is:

ddetg=FullSimplify[Det[dmetric]]

g=-(r^2).

Angular Frequency, Energy, Angular Momentum and Marginally Stable Radius

From previous sections we define

Ωk[r_,a_]:=((((r^(3/2))+a))^(-1))
ek[r_,a_]:=((1-2(r^(-1))+a (r^(-3/2)))/(((1-3(r^(-1))+2a (r^(-3/2))))^(1/2)))
lk[r_,a_]:=(((1-2a (r^(-3/2))+(a^2)(r^(-2))) (r^(1/2)))/(((1-3(r^(-1))+2a (r^(-3/2))))^(1/2)))
rms[a_]:=Chop[
3-√(3+(a^2)+(-3+a) (((-1+a))^(1/3)) (((1+a))^(2/3))+(((-1+a))^(2/3)) (((1+a))^(1/3)) (3+a))+√(18-(-3+a) (((-1+a))^(1/3)) (((1+a))^(2/3))-(((-1+a))^(2/3)) (((1+a))^(1/3)) (3+a)+2 (-6+(a^2))-(16 (a^2))/(√(3+(a^2)+(-3+a) (((-1+a))^(1/3)) (((1+a))^(2/3))+(((-1+a))^(2/3)) (((1+a))^(1/3)) (3+a))))]

Angular Frequency Plot

For a=0, 1.

Ωkplt0=Plot[Ωk[r,0.],{r,rms[.00001],10},Frame->True,FrameLabel->{"r","(Subscript[Ω, K])"},RotateLabel->False,PlotRange->{{0,10},{0,0.5}}];
Ωkplt1=Plot[Ωk[r,1.],{r,rms[1.],10}];
Show[Ωkplt0,Ωkplt1]

Disk Velocity

We need the velocity (u^α) of the disk in the BL frame.

ucov={-ek[r,a],0,lk[r,a]};
ucon=FullSimplify[Table[Sum[inversedmetric[[i,j]]ucov[[j]],{j,1,3}],{i,1,3}]]

Coordinate Transformation

We generalize the approach we used for the ZAMO frame.
(Subscript[(L^α), (OverHat[t])])= ((((Subscript[e, (OverHat[t])])))^α)= (u^α) 
e(OverHat[α]) · e(OverHat[β]) =(Subscript[η, αβ])
and we find that the transformation from BL to disk coordinates is given by

bl2d={{ucon[[1]],0,((ucov[[3]])/(Δ^(1/2)))},{0,((Δ^(1/2))/r),0},{ucon[[3]],0,-((ucov[[1]])/(Δ^(1/2)))}};
Print[MatrixForm[bl2d]]
Print[MatrixForm[FullSimplify[Table[Sum[dmetric[[γ,δ]]
bl2d[[γ,α]]bl2d[[δ,β]],{γ,1,3},{δ,1,3}],{α,1,3},{β,1,3}]]]]
d2bl=FullSimplify[Inverse[bl2d]];
Print[MatrixForm[d2bl]]

Mass conservation

Mass flux vector

Introduce a proper mass surface density Σ. The mass flux vector  J=Σu has vanishing divergence. Therefore, as (Sqrt[-g])=r, 2πrJ = 2πr (u^r)Σ= -(OverDot[M])=const. Note that J < 0.

Stress energy tensor

Conservation Laws

The divergence of the full stress energy tensor is symmetric and must vanish. If we have a Killing vector ξ, satisfying Killing's equation, then
(Subscript[((Subscript[ξ, α])(T^(αβ))), ;β])=0.
Using the temporal and axial Killing vectors and imposing stationarity and axisymmetry, we immediately define a conserved energy flux vector
(ℰ^α)=-(Subscript[(T^α), t])
and a conserved angular momentum flux vector
(ℒ^α)=(Subscript[(T^α), ϕ]) .
Again, we are only interested in the radial components, ℰ = (ℰ^r), ℒ = (ℒ^r).

Mass transfer

In absence of pressure, the matter contribution to the stress tensor is given by (Subscript[(T^α), β])=((Σu)^α)(Subscript[u, β]). The contribution to the energy flux vector is ℰ = -(Subscript[Ju, t]) = JE and to the angular momentum flux vector is ℒ = (Subscript[Ju, ϕ])=JL

Torque

We do not have a good understanding of  the torque, (Subscript[(T^r), ϕ])= G. (Note that it might even be non-local.) However we do know that it will do work at a rate GΩ and there will be a corresponding energy flux. Therefore the contribution to the stress energy tensor in BL coordinates must have the form
 (Subscript[(T^r), t])=-(Subscript[((ΩT)^r), ϕ])
As a sanity check, assume that we have a Newtonian shear viscosity η. The rate of shear tensor is given by the symmetric 
trace free part of (Subscript[u, α;β]) projected perpendicular to u. (The acceleration (Subscript[a, α])= (u^β)(Subscript[u, α;β]) and the expansion θ=(Subscript[(u^α), ;α])vanish.)
(Subscript[σ, αβ]) = (Subscript[u, (α;β)])+(Subscript[u, (α])(Subscript[u, β);γ])(u^γ)
The viscous part of the stress energy tensor is then given by (Subscript[(T^α), β])=-2η (Subscript[(σ^α), β]). 
(eg MTW pp566-568 ). evaluating,

ducov=FullSimplify[Table[D[ucov[[α]],dcoord[[β]]]-Sum[daffine[[γ,α,β]]ucov[[γ]],{γ,1,3}],{α,1,3},{β,1,3}]];
thet=FullSimplify[Sum[inversedmetric[[α,β]]ducov[[α,β]],{α,1,3},{β,1,3}]];
pducov=FullSimplify[Table[ducov[[α,β]]+Sum[ucov[[α]]ducov[[β,γ]]ucon[[γ]],{γ,1,3}],{α,1,3},{β,1,3}]];
sigcov=FullSimplify[(1/2)Table[pducov[[α,β]]+pducov[[β,α]],{α,1,3},{β,1,3}]];
newtorque=-2η Table[Sum[inversedmetric[[α,γ]]sigcov[[γ,β]],{γ,1,3}],{α,1,3},{β,1,3}];
Print[MatrixForm[newtorque]]

The contribution to the energy flux automatically has the form described above.

FullSimplify[newtorque[[2,1]]+newtorque[[2,3]]Ωk[r,a]]

We can also transform this tensor into the disk frame

Print[MatrixForm[FullSimplify[Table[Sum[d2bl[[α,γ]]newtorque[[γ,δ]]bl2d[[δ,β]],{γ,1,3},{δ,1,3}],{α,1,3},{β,1,3}]]]]

Note that there is no energy flux in this frame as should be the case.

Radiation

We must next consider the contribution of the escaping radiation to the stress-energy budget. We assume that the disk is efficiently radiative so that the energy that is dissipated is immediately radiated isotropically in the disk frame. If we denote the radiation flux leaving one surface by F, the four dimensional stress energy tensor has non-zero components 
(T^((OverHat[θ])  (OverHat[t])))= (T^((OverHat[t])(OverHat[θ])))= F; (T^((OverHat[t])  (OverHat[t])))= 4F; (T^((OverHat[r])  (OverHat[r]))) = (T^((OverHat[θ])  (OverHat[θ]))) = (T^((OverHat[ϕ])  (OverHat[ϕ])))= ((4F)/3).
Note that this is trace-free, as it should be. 
In BL coordinates, -(Subscript[(T^θ), t])= (FE/r), while the vertical angular momentum flux is  (Subscript[(T^θ), ϕ])= (FL/r).

Torque and Flux

Energy and Angular Momentum Conservation

The fluxes of radiative energy and angular momentum escaping from the disk can be equated to the minus the divergence of the 3D stress tensor (((-g))^(-1/2)) (Subscript[((g^(1/2))(Subscript[(T^β), α])), ,β]).  In BL coordinates, this implies that
(1/r)(d/dr)r(G + JL) + 2FL=0
(1/r)(d/dr)r(GΩ + JE) + 2FE=0

Solution

Solving these two equations and mass conservation, using dE=ΩdL and assuming that G and F vanish at (Subscript[r, min]) gives
G  = ((OverDot[M])/(2πr))(((Subscript[(∫^r), (Subscript[r, min])])dL(E-ΩL))/(E-ΩL))
and 
F = -((OverDot[M])/(2πr))((dΩ)/dr)(((Subscript[(∫^r), (Subscript[r, min])])dL(E-ΩL))/(((E-ΩL))^2))
which reduce to the nonrelativistic expressions. Note that these are independent of the nature of the torque excepting that it and the associated dissipation be local.

a=.95;
eoml[r_,a_]:=FullSimplify[ek[r,a]-Ωk[r,a]lk[r,a]]
dldr[r_,a_]:=FullSimplify[D[lk[r1,a],r1]]/.r1->r
domdr[r_,a_]:=FullSimplify[D[Ωk[r1,a],r1]]/.r1->r
torque[r_]:=((NIntegrate[eoml[rp,a]dldr[rp,a],{rp,rms[a],r}])/(2π r eoml[r,a]))
flux[r_]:=(-NIntegrate[eoml[rp,a]dldr[rp,a],{rp,rms[a],r}]domdr[r,a])/(2π r ((eoml[r,a])^2))

Analytic expressions can be given for these functions (eg Krolik pp150-152), but they are not very illuminating and it is easier to compute them numerically.

Plot[torque[r],{r,rms[a],10},Frame->True,FrameLabel->{"r","G"}]
Plot[flux[r],{r,rms[a],10},Frame->True,FrameLabel->{"r","F"}]
Clear[a]

Homework 4 (due Thursday May15)

Accretion Disk Image

1. Use the prescription given in class to compute the image of an accretion disk, viewed from a specific inclination.
a) First make a disk model assuming zero torque at the marginally stable circular orbit. Evaluate the rate of energy dissipation in the disk frame in a steady state and the temperature of the radiation in the disk frame as a function of radius.
Ignore returning radiation.
b) Compute the Doppler-boosted temperature as a function of direction and radius in the BL frame.
c) Introduce a grid of points on the observer's sky and for each point, trace a ray back onto the disk or the hole/gap and assign the appropriate intensity at a given observer frequency.
d) Hence exhibit the image of the disk.
e) Compute the observed flux and explore how it varies with inclination and frequency.

5. Thick Disks

Motivation

Low Mass Accretion Rate

When the gas accretes at a low rate, the electron-ion equipartition is too slow to keep the electrons in equilibrium 
with the ions, assuming that the ions acquire all the dissipated energy. We can ignore radiative loss. This is the ion-supported disk or torus.

High Mass Accretion Rate

At high mass accretion rate the escaping radiation may be trapped and advected into the hole. We can also ignore radiative loss. This is the radiation-supported disk or torus.

Equilibrium Flow

Metric  etc

Clear[coord,metric,inversemetric];
coord={t,r,θ,ϕ};
metric={{-(((r^2)-2r+(a^2)((Cos[θ])^2))/((r^2)+(a^2)((Cos[θ])^2))),0,0,((-2a r ((Sin[θ])^2))/((r^2)+(a^2)((Cos[θ])^2)))},
{0,(((r^2)+(a^2)((Cos[θ])^2))/((r^2)-2r+(a^2))), 0,0},{0,0,(r^2)+(a^2)((Cos[θ])^2),0},{((-2a r ((Sin[θ])^2))/((r^2)+(a^2)((Cos[θ])^2))),0,0,((((r^2) +(a^2))((r^2)+(a^2)((Cos[θ])^2))+2r (a^2)((Sin[θ])^2))/((r^2)+(a^2)((Cos[θ])^2)))((Sin[θ])^2)}};
ρ=(Sqrt[(r^2)+(a^2)((Cos[θ])^2)]);
Δ=(r^2)-2r+(a^2);
Σ = (Sqrt[((((r^2)+(a^2)))^2)-Δ (a^2)((Sin[θ])^2)]);
inversemetric=FullSimplify[Inverse[metric]];
ek[r_,a_]:=((1-2(r^(-1))+a (r^(-3/2)))/(((1-3(r^(-1))+2a (r^(-3/2))))^(1/2)))
lk[r_,a_]:=(((1-2a (r^(-3/2))+(a^2)(r^(-2))) (r^(1/2)))/(((1-3(r^(-1))+2a (r^(-3/2))))^(1/2)))
Ωk[r_,a_]:=((((r^(3/2))+a))^(-1))
rh[a_]:=1+(Sqrt[1-(a^2)])
rmb[a_]:=2-a+2(((1-a))^(1/2))
rms[a_]:=Chop[
3-√(3+(a^2)+(-3+a) (((-1+a))^(1/3)) (((1+a))^(2/3))+(((-1+a))^(2/3)) (((1+a))^(1/3)) (3+a))+√(18-(-3+a) (((-1+a))^(1/3)) (((1+a))^(2/3))-(((-1+a))^(2/3)) (((1+a))^(1/3)) (3+a)+2 (-6+(a^2))-(16 (a^2))/(√(3+(a^2)+(-3+a) (((-1+a))^(1/3)) (((1+a))^(2/3))+(((-1+a))^(2/3)) (((1+a))^(1/3)) (3+a))))]

Fluid

Write the enthalpy as 
w = nm +((γP)/((γ-1)))
where n is the baryon density and m is the rest mass per baryon. This is good for ion tori with γ=5/3 and radiation tori with γ=4/3. The fluid stress tensor is 
 (T^(αβ))= (wu^α)(u^β)+ (Pg^(αβ))
 It will also be useful to define an entropy function, which is monotonically related to the true entropy 
 S=((P^(1/γ))/nm).

Euler equation

Ignore inflow and assume axisymmetry. Project the stress energy conservation law perpendicular to u=[(u^t),0,0,(u^ϕ)]. Introduce a, b, c for r,θ, i, j, k, for t, ϕ.
(Subscript[P, αβ])(Subscript[(T^(βγ)), ;γ])
=> (Subscript[wa, a])+(Subscript[P, ,a])=0
where (Subscript[a, α])= (u^β)(Subscript[u, α;β]) is the acceleration.
Using (u^β)(Subscript[u, β;α])=0, it is helpful to write 
(Subscript[a, α])= (u^β)((Subscript[u, α,β])-(Subscript[u, β,α]))
Hence 
((Subscript[P, ,a])/w)= (u^t)(Subscript[u, t,a])+(u^ϕ)(Subscript[u, ϕ,a])
which implies
((∇P)/w)= (u^t)∇(Subscript[u, t])+(u^ϕ)∇(Subscript[u, ϕ]), 
where ∇ here represents the two dimensional derivative [(∂/(∂r)),(1/r)(∂/(∂θ))] and we can map the meridional structure onto the Euclidean plane.

Fluid Angular Momentum

It is convenient to replace the dynamical angular momentum L with the fluid angular momentum
H =(L/( E))=((-(Subscript[u, ϕ]))/(Subscript[u, t]))
Using 
(u^t)=(1/(E(1-ΩH)))
we obtain 
-((∇P)/w)=∇ln E -((Ω∇H)/(1-ΩH))
For a Keplerian orbit, the fluid angular momentum is:

hk[r_,a_]:=((lk[r,a])/(ek[r,a]))

Barytropic Disks

von Zeipel's theorem

To proceed we need to make an assumption about the equation of state. A convenient (and possibly relevant) choice is that it is barytropic, ie w=w(P). This implies that Ω=Ω(H). This is the relativistic von Zeipel theorem. (The non-relativistic theorem is that 'isogyres' are cylinders.)
If we use the relation,
Ω=(((g^(tϕ))-(g^(ϕϕ))H)/((g^tt)-(g^(tϕ))H)) 
then we see that this defines isorotational surfaces on which Ω and H are simultaneously constant, independent of the equation of state.

Isobars

Given a relation Ω(H), we can define
Φ= ln E -∫((dH Ω)/(1-ΩH))  
where we compute E using
E=(((-(g^tt)+2(g^(tϕ))H-(g^(ϕϕ))(H^2)))^(-1/2))
The isobars are then surfaces of constant Φ obtained by solving
(dP/(dΦ))=-w(P)

Equatorial Angular Momentum

Consider the plot of the two definitions of angular momentum with radius for equatorial Keplerian orbits.

a=0.95;
Plot[{lk[r,a],hk[r,a]},{r,rmb[a],6},Frame->True,FrameLabel->{"r","L,H"},RotateLabel->False]

Clearly (Subscript[H, mb])=(Subscript[L, mb]). The minima also coincide in radius (because dE = Ω dL for a circular geodesic). Now this radius is simply (Subscript[r, ms]). To see this recall that the stable circular orbits for given L are located at the minima of R(r)≡(((r'))^2); the unstable orbits at the maxima and as we reduce L, the two solutions join at a point of inflexion in R(r) and then disappear. Therefore, for (Subscript[H, ms ])< H < (Subscript[H, mb]) (or (Subscript[L, ms]) < L < (Subscript[L, mb])), there are two bound (E<1) geodesic orbits - the inner one unstable, the outer one stable. 
When we introduce pressure, there will be an additional, radial force which will modify the orbits except at the pressure maximum (Subscript[r, max]) and the inner radius (Subscript[r, cusp]) where the pressure vanishes. (We assume that the gas continues to flow slowly inwards until it reaches a point where the pressure vanishes at the disk surface and it is on an unstable circular geodesic from which it quickly plunges into the hole.) For (Subscript[r, cusp])<r<(Subscript[r, max]), the angular velocity is super-Keplerian; for (Subscript[r, max ])<r, it is sub-Keplerian.

Constant Angular Momentum Disks

The simplest case is constant H. We choose a value in the middle of the allowed range.

a=.95;h0=.5(hk[rmb[a],a]+hk[rms[a],a]);
rcusp=r/.FindRoot[hk[r,a]==h0,{r,rms[a]-.1}];ecusp=ek[rcusp,a];Ωcusp=Ωk[rcusp,a];
rmax=r/.FindRoot[hk[r,a]==h0,{r,rms[a]+.1}];emax=ek[rmax,a];Ωmax=Ωk[rmax,a];
Print["(Subscript[r, mb])=",rmb[a],"   (Subscript[r, ms])=",rms[a],"   H=",h0]
Print["(Subscript[r, cusp])=",rcusp,"  (Subscript[E, cusp])=",ecusp,"  (Subscript[Ω, cusp])=",Ωcusp]
Print["(Subscript[r, max])=",rmax,"  (Subscript[E, max])=",emax,"  (Subscript[Ω, max])=",Ωmax]

Isorotational surfaces

The isorotational surfaces, Ω = constant or, equivalently, H=constant.

hxy[x_,y_]:=h0;rout=2rmax;Ωxy[x_,y_]:=((inversemetric[[1,4]]-h inversemetric[[4,4]])/(inversemetric [[1,1]]-h inversemetric[[1,4]]))/.{r->Sqrt[(x^2)+(y^2)],θ->ArcTan[y,x],h->hxy[x,y]}
Ωcusp=Ωxy[rcusp,0];Ωout=Ωxy[rout,0];ΔΩ=.2(Ωcusp-Ωout);
Ωcon=ContourPlot[Ωxy[x,y],{x,.99rcusp,rout},{y,-rout,rout},Contours->Table[Ωout+i ΔΩ,{i,1,5}],ContourShading->False,AspectRatio->Automatic,PlotPoints->100];

Isobars

The isobars are surfaces of constant Φ=lnE. As H is increased, and the cusp approaches (Subscript[r, mb]),  and the binding energy (1-H) of the gas that flows over into the hole decreases. If (Subscript[r, cusp])=(Subscript[r, mb]), it is possible to accrete without liberating any binding energy.

en[x_,y_,h_]:=(((-inversemetric[[1,1]]+2inversemetric[[1,4]]h-inversemetric[[4,4]](h^2)))^(-.5))/.{r->Sqrt[(x^2)+(y^2)],θ->ArcTan[y,x]}
Φxy[x_,y_]:=Log[en[x,y,h]]/.h->hxy[x,y];
isobars=ContourPlot[Φxy[x,y],{x,rcusp,rout},{y,-rout,rout},Contours->Table[Log[emax+.2i (ecusp-emax)] ,{i,1,5}],AspectRatio->Automatic,ContourShading->False,ContourStyle->Red,PlotPoints->100];
Show[Graphics[{Disk[{0,0},rh[a],{-.5π,.5π}],Point[{rmb[a],0}],Point[{rcusp,0}],Point[{rms[a],0}],Point[{rmax,0}],Line[{{0,0},{2rmax,0}}]}],Ωcon,isobars, PlotRange->{{0,2rmax},{-rmax,rmax}},Frame->True]

Pressure Variation

The density and entropy are also constant on the isobars. The simplest choice is to suppose that the entropy function, S, is constant in which case,  the equation of hydrostatic equilibrium is integrable.

Simplify[Integrate[(((((γ p)/(γ-1))+((p^(1/γ))/S)))^(-1)),p]]

1+(γ/(γ-1))(P^((γ-1)/γ))S=((Subscript[E, cusp])/E) 
Set S=1 wlog and specialize to a radiation-dominated torus.We can now label the isobars.

p[en_]:=((((((γ-1))/γ) ((ecusp/en)-1)))^(γ/(γ-1)))/.γ->(4/3)
Plot[p[en],{en,emax,ecusp},PlotRange->All,Frame->True,RotateLabel->False,FrameLabel->{"E","P"}]

General angular momentum variation

The angular velocity must have the Keplerian value at the cusp and the pressure maximum, (Subscript[r, max]). It must also rise eventually as H=(Subscript[H, 0]) (r^(1/2)), (Subscript[H, 0]) <1. We can interpolate a simple, radial variation H(r) with these properties.

Clear[a,h,rmax];a=0.9;h0=0.9;
rcusp=rmb[a]+0.1;hcusp=hk[rcusp,a];ecusp=ek[rcusp,a];
rmax=8;hmax=hk[rmax,a];emax=ek[rmax,a];
coeff=Solve[{(c1/(rcusp^0.5))+c2+h0 (rcusp^0.5)==hcusp,(c1/(rmax^0.5))+c2+h0 (rmax^0.5)==hmax},{c1,c2}][[1]];
heq[r_]:=If[r<rcusp,hcusp,(c1/(r^0.5))+c2+h0 (r^0.5)/.coeff]
rout=30;hout=heq[rout];
Plot[{hk[r,a],heq[r]},{r,rcusp,rout},Frame->True,FrameLabel->{"r","H"}]

We next invert this relation to obtain r(H)

req[h_]:=(((c2^2)-2c2 h+(h^2)-2c1 h0+ (Sqrt[-4 (c1^2) (h0^2)+(((- (c2^2)+2c2 h-(h^2)+2 c1 h0))^2)]))/(2(h0^2)))/.coeff

and compute the function Ω(H) in the equatorial plane.

gtt[r1_]:=inversemetric[[1,1]]/.{r->r1,θ->π/2}
gtϕ[r1_]:=inversemetric[[1,4]]/.{r->r1,θ->π/2}
gϕϕ[r1_]:=inversemetric[[4,4]]/.{r->r1,θ->π/2}
Ωh[h_]:=If[h>hcusp,((gtϕ[r]-h gϕϕ[r])/(gtt[r]-h gtϕ[r])),10]/.r->req[h]

The same relation Ω(H) must be satisfied everywhere giving a second relation for Ω(H,x,y). Eliminating Ω, between these two relations gives the required function H(x,y). Surfaces of constant H are the gyrentropes.

Ωxyh[x_,y_,h_]:=If[h>hcusp,FullSimplify[((inversemetric[[1,4]]-h inversemetric[[4,4]])/(inversemetric [[1,1]]-h inversemetric[[1,4]]))/.{r->Sqrt[(x^2)+(y^2)],θ->ArcTan[y,x]}],10]
hxy[x_,y_]:=If[x>rcusp,h/.FindRoot[Ωxyh[x,y,h]==Ωh[h],{h,{heq[x],heq[x]+.1}}],0]
vonzeipel=ContourPlot[hxy[x,y],{x,.001,rout},{y,-.5rout,.5rout},Contours->{hcusp,heq[5],heq[10],heq[20]},ContourShading->False]

We next compute the isobars Φ(x,y)= ln E(x,y,H) -(Subscript[(∫^H), (Subscript[H, cusp])])((dHΩ)/(1-ΩH))=constant.

enxyh[x_,y_,h_]:=(((Max[-inversemetric[[1,1]]+2inversemetric[[1,4]]h-inversemetric[[4,4]](h^2),.01]))^(-.5))/.{r->Sqrt[(x^2)+(y^2)],θ->ArcTan[y,x]}
ndsg=NDSolve[{g'[h]==((Ωh[h])/(1-h Ωh[h])),g[hcusp]==0},g,{h,hcusp,hout}];
g1[h1_]:=g[h]/.ndsg[[1]]/.h->h1
phixy[x_,y_]:=If[x>rcusp&&hxy[x,y]>hcusp,(Log[((enxyh[x,y,h])/ecusp)]-g1[h])/.h->hxy[x,y],1]
isobars=ContourPlot[phixy[x,y],{x,0.001,rout},{y,-.5rout,.5rout},Contours->{0,-.04,-.07,-.09,-.1,-.105},ContourShading->False,PlotPoints->100]

Isobars and von Zeipel surfaces

The plot of the isobars exhibits a very steep funnel with this choice of parameters.

holeplt=ContourPlot[(x^2)+(y^2)-(((1+(Sqrt[1-(a^2)])))^2),{x,0.001,rout},{y,-0.5 rout,0.5 rout},Contours->{0},PlotPoints->100,DisplayFunction->Identity];Show[holeplt,vonzeipel,isobars,DisplayFunction->$DisplayFunction]

Gyrentropic Disks

An alternative prescription is that the disks are gyrentropic. This means that L = L(S). This may be more 
physical for fluid disks. The reason is that this is the condition for marginal stability to axisymmetric convection.
It is convenient to introduce the relativistic Bernoulli function
B=(wE/nm)
and it is then straightforward to show that  
∇ln B=((w-nm)/w)∇ln S+((Ω∇H)/(1-ΩH))
so that the "gyrentropes" are also "isoberns".  We use a similar angular momentum variation to before

Clear[a,h,rmax];a=0.9;h0=0.75;γ=1.67;n=0.58;
agam=((3 γ-5-2 n γ+2 n)/(γ-5-2 n γ+2 n));
rcusp=0.5 (rmb[a]+rms[a]);
hcusp=hk[rcusp,a];ecusp=ek[rcusp,a];
rmax=4;hmax=hk[rmax,a];emax=ek[rmax,a];
coeff=Solve[{(c1/(rcusp^0.5))+c2+h0 (rcusp^0.5)==hcusp,(c1/(rmax^(0.5`)))+c2+h0 (rmax^0.5)==hmax},{c1,c2}][[1]];
hr[r_]:=(c1/(r^0.5))+c2+h0 (r^0.5)/.coeff
dhr[r_]:=-((0.5 c1)/(r^1.5))+((0.5 h0)/(r^0.5))/.coeff
rin=rcusp-0.03;hin=hr[rin];
rh[h_]:=If[h<hin,rin,(1/(2 (h0^2)))((c2^2)-2 c2 h+(h^2)-2 c1 h0+(Sqrt[-4 (c1^2) (h0^2)+(((-(c2^2)+2 c2 h-(h^2)+2 c1 h0))^2)]))]/.coeff
rout=10;hout=hr[rout];
sr[r_]:=Min[(r^(-(agam/(1-agam)))),10];
lnsr[r_]:=Log[sr[r]]
dlnsr[r_]:=10000 (lnsr[r+0.0001]-lnsr[r])
gtt[r1_]:=inversemetric[[1,1]]/.{r->r1,θ->(π/2)}
gtϕ[r1_]:=inversemetric[[1,4]]/.{r->r1,θ->(π/2)}
gϕϕ[r1_]:=inversemetric[[4,4]]/.{r->r1,θ->(π/2)}
energyr[r_]:=(1/((Max[-gtt[r]+2 gtϕ[r] hr[r]-gϕϕ[r] ((hr[r])^2),0.01])^0.5))
Ωr[r_]:=((gtϕ[r]-hr[r] gϕϕ[r])/(gtt[r]-hr[r] gtϕ[r]))

nds1=NDSolve[{bernr'[r]==(bernr[r]-energyr[r]) dlnsr[r]+((Ωr[r] bernr[r] dhr[r])/(1-Ωr[r] hr[r])),bernr[rcusp]==energyr[rcusp]},bernr,{r,rcusp,rout}]

fig1=Plot[{hr[r],hkep[r,a]},{r,rcusp,rout},Frame->True,PlotRange->{{0,rout},{2.4,3.6}},PlotStyle->{{Thickness[0.005]},{Thickness[0.005],Dashing[{0.02,0.02}]}},FrameLabel->{"r","\!\(\*OverscriptBox[\(L\), \(~\)]\)"},RotateLabel->False];
fig2=Graphics[{Line[{{rmb[a],2.4},{rmb[a],3}}],Line[{{rcusp,2.4},{rcusp,3.3}}],Line[{{rms[a],2.4},{rms[a],3.3}}],Line[{{rmax,2.4},{rmax,3.3}}],Text["\!\(\*SubscriptBox[\(r\), \(mb\)]\)",{1.2,3}],Text["\!\(\*SubscriptBox[\(r\), \(c\)]\)",{1.7,3.3}],Text["\!\(\*SubscriptBox[\(r\), \(ms\)]\)",{2.8,3.3}],Text["\!\(\*SubscriptBox[\(r\), \(\(m\)\( \)\)]\)",{4.5,3.3}]}];
fig3=Show[fig1,fig2];
fig4=Plot[{sr[r],bernr[r]/.nds1[[1]]},{r,rcusp,10},Frame->True,PlotRange->{{0,rout},{0.6,1}},PlotStyle->{{Thickness[0.01]},{Thickness[0.01],Dashing[{0.02,0.02}]}},FrameLabel->{"r","S,\!\(\*OverscriptBox[\(B\), \(~\)]\)"},RotateLabel->False,DisplayFunction->Identity];
fig5=Graphics[{Line[{{rmb[a],0.6},{rmb[a],0.8}}],Line[{{rcusp,0.6},{rcusp,0.9}}],Line[{{rms[a],0.6},{rms[a],0.9}}],Line[{{rmax,0.6},{rmax,0.9}}]}];
fig6=Show[fig4,fig5];
gyr=Show[GraphicsRow[{fig3,fig6}]]
bernh[h_]:=bernr[r]/.nds1[[1]]/.r->rh[h]
dbernh[h_]:=10000 (bernh[h+0.0001]-bernh[h])
lnsh[h_]:=Log[sr[rh[h]]]
dlnsh[h_]:=10000 (lnsh[h+0.0001]-lnsh[h])
Ωxyh[x_,y_,h_]:=If[x<rin,omegar[rin],((inversemetric[[1,4]]-h inversemetric[[4,4]])/(inversemetric[[1,1]]-h inversemetric[[1,4]]))/.{r->(Sqrt[(x^2)+(y^2)]),θ->ArcTan[y,x]}]
energyxyh[x_,y_,h_]:=If[x<rin,energyr[rin],1/((Max[-inversemetric[[1,1]]+2 inversemetric[[1,4]] h-inversemetric[[4,4]] (h^2),0.01])^0.5)/.{r->(Sqrt[(x^2)+(y^2)]),θ->ArcTan[y,x]}]
hxy[x_,y_]:=If[x<rin,hin,h/.FindRoot[dbernh[h]==(bernh[h]-energyxyh[x,y,h]) dlnsh[h]+((Ωxyh[x,y,h] bernh[h])/(1-Ωxyh[x,y,h] h)),{h,hr[x],hr[x]+0.01}]]
Ωxy[x_,y_]:=((inversemetric[[1,4]]-h inversemetric[[4,4]])/(inversemetric[[1,1]]-h inversemetric[[1,4]]))/.{r->(Sqrt[(x^2)+(y^2)]),θ->ArcTan[y,x]}/.h->hxy[x,y]
energyxy[x_,y_]:=1/((Max[-inversemetric[[1,1]]+2 inversemetric[[1,4]] h-inversemetric[[4,4]] (h^2),0.01])^0.5)/.{r->(Sqrt[(x^2)+(y^2)]),θ->ArcTan[y,x]}/.h->hxy[x,y]
bernxy[x_,y_]:=bernh[hxy[x,y]]
sxy[x_,y_]:=(E^(lnsh[hxy[x,y]]))
pxy[x_,y_]:=If[x<rcusp+Abs[y]||bernxy[x,y]<energyxy[x,y],0,((((((γ-1) (((bernxy[x,y])/(energyxy[x,y]))-1))/(γ sxy[x,y]))))^(γ/(γ-1)))]
holeplt=ContourPlot[(x^2)+(y^2)-(((1+(Sqrt[1-(a^2)])))^2),{x,0.001`,rout},{y,-0.5 rout,0.5 rout},Contours->{0},PlotPoints->100,AspectRatio->Automatic];gyrentropes=ContourPlot[hxy[x,y],{x,0.0001,rout},{y,-0.5 rout,0.5 rout},Contours->{hxy[3,0],hxy[4,0],hxy[5,0],hxy[6,0],hxy[7,0],hxy[8,0],hxy[9,0]},ContourShading->False,PlotPoints->50,DisplayFunction->Identity];
pcontours={pxy[2.15,0],pxy[2.5,0],pxy[3,0],pxy[3.5,0]};
tor=Show[holeplt,gyrentropes]

Stability

denergyr[r_]:=10000 (energyr[r+0.0001]-energyr[r])
dpxy[r_]:=((denergyr[r])/(energyr[r]))-((Ωr[r] dhr[r])/(1-Ωr[r] hr[r]))
bernour[r_]:=bernr[r]/.nds1
dsr[r_]:=10000 (sr[r+0.0001]-sr[r])
dΩr[r_]:=10000 (Ωr[r+0.00001]-Ωr[r])
detg[r1_]:=((metric[[1,4]])^2)-metric[[1,1]] metric[[4,4]]/.{r->r1,θ->(π/2)}
gam[r_]:=((((energyr   r])^4) dhr[r])/(detg[r]))-((dΩr[r])/(((1-Ωr[r] hr[r]))^2))
hoilandone[r_]:=(((((energyr[r])/(bernour[r]-energyr[r]))))^2) ((pxy[r,0])^(1/γ)) dpxy[r] dsr[r]-gam[r] dhr[r]
Plot[hoilandone[r],{r,2.5,10}]

6. Tides and Oscillations

Tidal Effects

Motivation

The discovery of large numbers of young stars in close orbit around our central black hole suggests strongly that the supply of stars in galactic nuclei is much greater than suspected.  It further raises the possibility that stars can be captured onto relativistic orbits quite readily and that tidal disruption events might be quite common in nearby galaxies.

Metric etc

Clear[coord,metric,inversemetric,affine,a];
coord={t,r,θ,ϕ};
metric={{-(((r^2)-2r+(a^2)((Cos[θ])^2))/((r^2)+(a^2)((Cos[θ])^2))),0,0,((-2a r ((Sin[θ])^2))/((r^2)+(a^2)((Cos[θ])^2)))},
{0,(((r^2)+(a^2)((Cos[θ])^2))/((r^2)-2r+(a^2))), 0,0},{0,0,(r^2)+(a^2)((Cos[θ])^2),0},{((-2a r ((Sin[θ])^2))/((r^2)+(a^2)((Cos[θ])^2))),0,0,((((r^2) +(a^2))((r^2)+(a^2)((Cos[θ])^2))+2r (a^2)((Sin[θ])^2))/((r^2)+(a^2)((Cos[θ])^2)))((Sin[θ])^2)}};
ρ=(Sqrt[(r^2)+(a^2)((Cos[θ])^2)]);
Δ=(r^2)-2r+(a^2);
Σ = (Sqrt[((((r^2)+(a^2)))^2)-Δ (a^2)((Sin[θ])^2)]);
inversemetric=FullSimplify[Inverse[metric]];
affine=FullSimplify[Table[(1/2)Sum[(inversemetric[[α,δ]])
(D[metric[[δ,β]],coord[[γ]]]+D[metric[[δ,γ]],coord[[β]]]-D[metric[[β,γ]],coord[[δ]]]),{δ,1,4}],{α,1,4},{β,1,4},{γ,1,4}]];
riemann=FullSimplify[Table[D[affine[[α,β,δ]],coord[[γ]]]-D[affine[[α,β,γ]],coord[[δ]]]+Sum[affine[[α,ϵ,γ]]affine[[ϵ,β,δ]]-affine[[α,ϵ,δ]]affine[[ϵ,β,γ]],{ϵ,1,4}],{α,1,4},{β,1,4},{γ,1,4},{δ,1,4}]];
ricci=FullSimplify[Table[Sum[riemann[[γ,α,γ,β]],{γ,1,4}],{α,1,4},{β,1,4}]];
Print[MatrixForm[ricci]]
Ωk[r_,a_]:=((((r^(3/2))+a))^(-1))
ek[r_,a_]:=((1-2(r^(-1))+a (r^(-3/2)))/(((1-3(r^(-1))+2a (r^(-3/2))))^(1/2)))
lk[r_,a_]:=(((1-2a (r^(-3/2))+(a^2)(r^(-2))) (r^(1/2)))/(((1-3(r^(-1))+2a (r^(-3/2))))^(1/2)))
rms[a_]:=Chop[
3-Sign[a]√(3+(a^2)+(-3+a) (((-1+a))^(1/3)) (((1+a))^(2/3))+(((-1+a))^(2/3)) (((1+a))^(1/3)) (3+a))+√(18-(-3+a) (((-1+a))^(1/3)) (((1+a))^(2/3))-(((-1+a))^(2/3)) (((1+a))^(1/3)) (3+a)+2 (-6+(a^2))-Sign[a](16 (a^2))/(√(3+(a^2)+(-3+a) (((-1+a))^(1/3)) (((1+a))^(2/3))+(((-1+a))^(2/3)) (((1+a))^(1/3)) (3+a))))]

Equatorial Riemann tensor

We confine our attention to stars on circular orbits. These might be formed in an accretion disk and migrate slowly inwards or be captured and dragged into the equatorial plane by interaction with a disk. The Riemann tensor becomes

eriemann=riemann/.{θ->π/2};

Transformation to the star frame

Transform from BL to star coordinates

ucov={-ek[r,a],0,0,lk[r,a]};
ucon=FullSimplify[Table[Sum[inversemetric[[i,j]]ucov[[j]],{j,1,4}],{i,1,4}]]/.θ->π/2
bl2d={{ucon[[1]],0,0,((ucov[[4]])/(Δ^(1/2)))},{0,((Δ^(1/2))/r),0,0},{0,0,(1/r),0},{ucon[[4]],0,0,-((ucov[[1]])/(Δ^(1/2)))}};
d2bl=FullSimplify[Inverse[bl2d]];

Newtonian Tidal Tensor

Transform the Riemann tensor to a coordinate basis free-falling with the star using the transformation matrices
(Subscript[(R^(OverHat[α])), (OverHat[β])(OverHat[γ])(OverHat[δ])])=(Subscript[(L^(OverHat[α])), ϵ])(Subscript[(L^ζ), (OverHat[β])])(Subscript[(L^η), (OverHat[γ])])(Subscript[(L^θ), (OverHat[δ])])(Subscript[(R^ϵ), ζηθ])
As 
(Subscript[∇, u])(Subscript[∇, u])n +R[,u,n,u]=0
the co-moving, Newtonian tidal tensor can be evaluated as
(Subscript[T, (OverHat[j])]^î) ≡ (Subscript[Φ, ,(OverHat[j])]^(OverHat[i]))=(Subscript[(R^(OverHat[i])), (OverHat[t])(OverHat[j])(OverHat[t])])
It should have a vanishing trace according to Laplace's equation.

oriemann=FullSimplify[Table[Sum[d2bl[[i,m]]bl2d[[n,j]]bl2d[[o,k]]bl2d[[p,l]]eriemann[[m,n,o,p]],{m,1,4},{n,1,4},{o,1,4},{p,1,4}],{i,1,4},{j,1,4},{k,1,4},{l,1,4}]];
tidal=Table[oriemann[[i,1,j,1]],{i,2,4},{j,2,4}];
FullSimplify[Tr[tidal]]
Print[MatrixForm[tidal]]

Plots of tidal forces

The radial element is negative as the force is expansive; the two tangential elements are positive and compressive.
As  r->∞, (Subscript[T, (OverHat[r])]^(OverHat[r]))->-2(Subscript[T, (OverHat[θ])]^(OverHat[θ])),-2(Subscript[T, (OverHat[ϕ])]^(OverHat[ϕ])), -(2/(r^3))as it should.   
As r->(Subscript[r, ms]),(Subscript[T, (OverHat[r])]^(OverHat[r]))->-(3/2)(Subscript[T, (OverHat[θ])]^(OverHat[θ])),-3(Subscript[T, (OverHat[ϕ])]^(OverHat[ϕ]))->-3(a->1),=-(1/72)(a=0),=-(1/243)(a=-1).

tidalplt=GraphicsColumn[Table[Plot[{tidal[[1,1]]/.a->0.999 i,tidal[[2,2]]/.a->0.999 i,tidal[[3,3]]/.a->0.999 i},{r,rms[0.999 i],10},Frame->True,FrameLabel->{"r","T"},RotateLabel->False],{i,-1,1}]]

Numerical Values

A star will be ripped apart outside (Subscript[r, ms])if its mean density < 70(Subscript[T, 1]^1)(Subscript[M, 8]^(-2))g (cm^(-3)). This is the Roche limit.

Rotation of Inertial Frame

The basis vectors associated with an inertial frame tied to a free-falling star will be parallel propagated and will rotate with respect to our BL basis and our choice of local orthonormal basis used above. Consider a 2D unit vector, i chosen to lie in the equatorial plane. Write
e = [-(Subscript[Ωi, (OverHat[ϕ])]),(Subscript[i, (OverHat[r])]),0,(Subscript[i, (OverHat[ϕ])])]
in BL, ensuring that u·e = 0, and ignoring normalization. The law of parallel propagation is (Subscript[∇, u])e=0 or
((Subscript[de, α])/(dτ))=(Subscript[e, β])(Subscript[(Γ^β), αγ])(u^γ) 
from which we obtain 
(((d^2)(Subscript[e, α]))/((dτ)^2))=(Subscript[(Γ^δ), αγ])(Subscript[(Γ^β), δϵ])(u^γ)(u^ϵ)(Subscript[e, β])
as u has only t,ϕ components and ((Subscript[de, α])/(dτ)) is a function of only r,θ. Evaluating the rhs,

Clear[a];
evec={-Ωk[r,a] iϕ,ir,0,iϕ};
rhs=FullSimplify[Table[Sum[eqaffine[[β,δ,ϵ]]eqaffine[[δ,α,γ]]evec[[β]]ucon[[γ]]ucon[[ϵ]],{β,1,4},{γ,1,4},{δ,1,4},{ϵ,1,4}],{α,1,4}]]

so that 
(((d^2)i)/((dτ)^2))=-(i/(r^3))
independent of a! The dragging of inertial frames is really a consequence of the choice of coordinates .
The combined rate of so-called "geodetic" plus "Lense-Thirring" precession (in the BL frame) is given by: 
(Subscript[Ω, gLT])=Ω-(1/((r^(3/2))(u^t)))=((1-(((1-(3/r)+((2a)/(r^(3/2)))))^(1/2)))/((r^(3/2))+a))=(3/(2(r^(5/2))))-(a/(r^3))+⋯
in agreement with the weak field result.

Ωgeo[r1_,a1_]:=Ωk[r1,a1]-((r1^(-3/2))/(ucon[[1]]/.{r->r1,a->a1}))
geoplt1=Plot[Ωgeo[r,-1],{r,rms[-1],10},Frame->True,FrameLabel->{"r","\!\(\*SubscriptBox[\(Ω\), \(geo\)]\)"},RotateLabel->False]
geoplt2=Plot[Ωgeo[r,0],{r,rms[0],10},Frame->True,FrameLabel->{"r","\!\(\*SubscriptBox[\(Ω\), \(geo\)]\)"},RotateLabel->False]
geoplt3=Plot[Ωgeo[r,1],{r,rms[1],10},Frame->True,FrameLabel->{"r","\!\(\*SubscriptBox[\(Ω\), \(geo\)]\)"},RotateLabel->False]

This might be observable if there is a transfer of energy from the orbit to the star as it passes through a resonance with an internal oscillation frequency. The generalization to an  inclined, eccentric orbit is given by Diener et al ApJ 479 164.

QPOs and Diskoseismology

Motivation

Many accretion disks exhibit quasi-perioidc oscillations - modulations of the hard X-ray emission with amplitude of a few percent and Q of order  10-100. Both BH and NS show these oscillations which are mostly supposed to be disk modes. Specific frequencies are associated with individual sources. NS often show pairs of frequencies which drift together with a spacing related to the spin frequency, or half this value. They do not drift much in frequency as the luminosity changes. If these are disk modes then the disk is only the clock as it radiates soft X-rays. Presumably the emission comes from a hot, Comptonizing corona which is modulated by a variable soft photon spectrum from the disk or a variable heating

Small perturbations for cold disks

When  the disk is cold and thin we can ignore pressure. The equations of motion are: 
(dt/dx)=t'= ((((r^2)+(a^2)){((r^2)+(a^(2 )))E-a L})/Δ)-a(a E (sin^2)θ -L) 
r''= 2r E {((r^2)+(a^2))E-aL} - (r-1){(r^2)+(((aE-L))^2)+Q} - r Δ
θ''' = (a^2) (1-(E^2))cosθ sinθ +(L^2)cotθ (csc^2)θ

Specializing to small departures from circular equatorial orbits with θ=π/2-δψ gives

ft[r_]:=((((r^2)+(a^2))(((r^2)+(a^2))e-a l))/(((r^2)-2r+(a^2))))-a(a e -l)
dfr[r_]:=2r e(((r^2)+(a^2))e-a l)-(r-1)((r^2)+(((a e-l))^2)+q)-r((r^2)-2r+(a^2))
dfθ[δψ_]:=((a^2)(1-(e^2))+(l^2))δψ

Vertical epicyclic motion

We see that the vertical equation of motion is 
(((d^(2 ))δψ)/(dx^2))=-((a^2)(1-(e^2))+(l^2))δψ
Evaluating the rhs for circular orbits, we find

kψ=FullSimplify[-((a^2)(1-(e^2))+(l^2))/.{e->ek[r,a],l->lk[r,a]}]

We also need(dt/( dx))

dtdx=FullSimplify[((((r^2)+(a^2))(((r^2)+(a^2))e-a l))/(((r^2)-2r+(a^2))))-a(a e -l)/.{e->ek[r,a],l->lk[r,a]}]

so that  the square of the vertical epicyclic frequency is

FullSimplify[kψ/(dtdx^2)]

((Subscript[κ, v])^2)=(Ω^2)(1-((4a)/(r^(3/2)))+3((a^2)/(r^2)))

(We are interested in frequencies measured from infinity.)
Expanding at large radius, we obtain
(Subscript[κ, v])≃Ω-((2a)/(r^3))

This corresponds to prograde, orbital Lense-Thirring precession as measured from infinity.  
(Subscript[Ω, LT])=((2a)/(r^3))

Radial epicyclic motion

Setting Q=0 in the linear approximation and using 
δr''=((R(r)'')/2)δr

kr=FullSimplify[D[
2r e(((r^2)+(a^2))e-a l)-(r-1)((r^2)+(((a e-l))^2))-r((r^2)-2r+(a^2)),r]/.{e->ek[r,a],l->lk[r,a]}]

and so the square of the radial epicyclic frequency is given by

-FullSimplify[(kr/(dtdx^2))]

((Subscript[κ, r])^2)=(Ω^2)(1-(6/r)+((8a)/(r^(3/2)))-((3(a^2))/(r^2)))

Again, expanding at large radius,
(Subscript[κ, r])≃Ω-(3/(r^(5/2)))
corresponding to prograde apsidal motion at a rate
(Subscript[Ω, aps])=(3/(r^(5/2)))

Frequency Plots

κv[r_,a_]:=(Sqrt[((3 (a^2)-4 a (Sqrt[r])+(r^2))/((r^2) (((a+(r^(3/2))))^2)))])
κr[r_,a_]:=(Sqrt[((-3 (a^2)+8 a (Sqrt[r])+(-6+r) r)/((r^2) (((a+(r^(3/2))))^2)))])
GraphicsColumn[Table[Plot[{Ωk[r,a],κv[r,a],κr[r,a]},{r,rms[a],3 rms[a]},Frame->True,PlotRange->{{rms[a],3 rms[a]},All},FrameLabel->{"r","Ω,\!\(\*SubscriptBox[\(κ\), \(r\)]\),\!\(\*SubscriptBox[\(κ\), \(v\)]\)"},RotateLabel->False],{a,-1,1}]]

Diskoseismology

There are three types of high frequency mode present in simple disk models. g-modes, p-modes and c-modes.  We must solve the equation
of continuity:
(Subscript[(((ρu)^α)), ;α])=0
energy
(Subscript[u, α])(Subscript[(T^(αβ)), ;α])=0
and momentum
(Subscript[P, αβ])(Subscript[(T^(βγ)), ;γ]) =0
In the simplest, limiting case, we ignore thermal modes and treat the gas as isentropic, and assume axisymmetry.
This leads eventually to a dispersion relation for g- and p-modes in the WKB limit 
((ω^2)-((Subscript[κ, r])^2))((ω^2)-((ΥΩ)^2))=(ω^2)(k^2)(s^2)
where k is the wave number and s is the sound speed and Υ is roughly 1. There are two types of mode, g-modes with vertical motions that have ω<(Subscript[κ, r]), Ω and are trapped below the maximum in the function (Subscript[Ω, r])(r), and p-modes wsith radial motion and ω>(Subscript[κ, r]), Ω that propagate outside a turing point.

7. Neutron Stars

Consider a non-rotating, spherically symmetric neutron star. The exterior metric must be reducible to the Schwarzschild metric through Birkhoff's theorem. Let us continue to define r in the same way as for black holes. The metric can then be written in a standard form:
(ds^2)=-(e^(2Φ))(dt^2)+((dr^2)/(1-2m/r))+(r^2)(((dθ)^2)+(sin^2)θ ((dϕ)^2)) 
without loss of generality, where Φ=Φ(r) and m=m(r). (Note absence of dtdr term as a photon should take the same Killing time t to traverse a radial interval forwards and backwards.) 
Now (u^α) =((e^(-Φ)),0,0,0), (Subscript[u, α]) =(-(e^Φ),0,0,0). Hence, 
(Subscript[P, ,r])=(wu^t)(Subscript[u, t,r]) = -(Subscript[wΦ, ,r]) or (dP/(dΦ))=-w=-(P+ρ)  
Next, evaluate the mixed Einstein tensor.

Clear[coord,metric,inversemetric,affine,ρ,m];
coord={t,r,θ,ϕ};
metric={{-Exp[2Φ[r]],0,0,0},{0,(((1-2m[r]/r))^(-1)), 0,0},{0,0,(r^2),0},{0,0,0,(r^2)((Sin[θ])^2)}};
inversemetric=FullSimplify[Inverse[metric]];
affine=Simplify[Table[(1/2)Sum[(inversemetric[[α,δ]])
(D[metric[[δ,β]],coord[[γ]] ]+
D[metric[[δ,γ]],coord[[β]] ]-D[metric[[β,γ]],coord[[δ]] ]),{δ,1,4}],
{α,1,4},{β,1,4},{γ,1,4}] ];
riemann=
Simplify[Table[D[affine[[α,β,δ]],coord[[γ]]]-D[affine[[α,β,γ]],coord[[δ]]]+Sum[affine[[α,ϵ,γ]]affine[[ϵ,β,δ]]-affine[[α,ϵ,δ]]affine[[ϵ,β,γ]],{ϵ,1,4}],{α,1,4},{β,1,4},{γ,1,4},{δ,1,4}]];
ricci=FullSimplify[Table[Sum[inversemetric[[α,δ]]riemann[[γ,δ,γ,β]],{γ,1,4},{δ,1,4}],{α,1,4},{β,1,4}]];
einstein=FullSimplify[ricci-(1/2)Tr[ricci]{{1,0,0,0},{0,1,0,0},{0,0,1,0},{0,0,0,1}}]

so that 
(Subscript[(G^t), t]) = -2((Subscript[m, ,r])/(r^2)) =8π (Subscript[(T^t), t])=-8πρ => m = 4 π ((∫_0)^r)dr'((r')^2)ρ(r')
and 
(Subscript[(G^r), r]) =2(1-((2m)/r))(Subscript[Φ, 'r])-((2 m)/(r^3)) =8(Subscript[((πT)^r), r])=8πP  
so that
(Subscript[P, ,r])=-(Subscript[wΦ, ,r])=-(((m+4((πr)^3)P)(ρ+P))/(r(r-2m))).
This is the Tolman-Oppenheimer-Volkoff equation. m is the interior mass energy and equals the m determined by eg orbits at large r at the surface of the star. m+4((πr)^3)P is the active gravitational mass and w the passive. The (1-2m/r) is kinematic. If the rest mass per baryon is μ and the proper density n, then 
(Subscript[m, rest])=4π∫dr (((r^2)(1-2m/r))^(-1/2))μ n. We can also define quasi-Newtonian internal energy density ρ-μn and gravitational  energy density -ρ[1-(((1-2m/r))^(1/2))] which agree with standard non-relativistic relations in the limit.
Note that at the surface where P=0, r=R, m=M, g=(1/ρ)(dP/dr)=(M/(((R^2)(1=2M/R))^(1/2))) and that 
One simple solution is ρ=3P = 3/56((πr)^2). More interesting is the incompressible star ρ const. so that m=4((πρr)^3)/3. The observed temperature is lower than the surface temperture by a factor (e^(Φ(R)))=(((1-2M/R))^(1/2)).

Simplify[Exp[Integrate[((4ρ)/((ρ+3P)(ρ+P))),P]+Integrate[((16π ρ r)/(3-8π ρ (r^2))),r]]];

and ((ρ+3P)/(ρ+P))=(({((3-8((πρr)^2))/(3-8((πρR)^2)))))^(1/2)).  The ccntral pressure is finite so long as R>9M/4.
Real neutron stars are commonly of M=1.4 solar masses and R=8-16km so that R/2M ~2-4.

8. Electromagnetic Effects

Motivation

There is good, circumstantial evidence that AGN jets, pulsar wind nebulae, gamma ray bursts (GRB) and now Galactic superluminal sources countain ultrarelativistic outflows emanating from black holes. The inferred outflow Lorentz factors are typically 10, (10^6), 300 and 2-30 respectively. It is commonly thought that the origin of these outflows is electromagnetic although the details remain controversial. Most models rely on some variant of unipolar induction - orbiting a flux Φ with angular velocity Ω generates an EMF ℰ. Under the prevailing electromagnetic conditions, the relevant impedance is that of free space, (Subscript[Z, 0])~100Ω to order of magnitude. Therefore the current is I ~ ℰ / (Subscript[Z, 0])and the electromagnetic (Poynting) power is L ~ ℰ I. The electromagnetic torque is ~ L/Ω.
The most powerful source we know is a GRB with L ~ (10^43)W. This implies ℰ ~ 3 ZV, I~30EA, Ω~(10^4)rad (s^(-1))and Φ~ 300 PWb with less powerful objects scaling accordingly.
Under these powerful conditions, the vacuum becomes an extremely powerful conductor - it only takes a MV in principle to break it down and less than 1 TV in practice. We can therefore suppose that electromagnetic field dominates the stress tensor around the prime mover. If the associated current closes nearby, the power will be transferred into a fluid form and we will be doing gas dynamics; if it closes in the emission region, we are dealing with Poynting jets.

Basic Equations

Metric etc

Clear[coord,metric,inversemetric,affine,a];
coord={t,r,θ,ϕ};
metric={{-(((r^2)-2r+(a^2)((Cos[θ])^2))/((r^2)+(a^2)((Cos[θ])^2))),0,0,((-2a r ((Sin[θ])^2))/((r^2)+(a^2)((Cos[θ])^2)))},
{0,(((r^2)+(a^2)((Cos[θ])^2))/((r^2)-2r+(a^2))), 0,0},{0,0,(r^2)+(a^2)((Cos[θ])^2),0},{((-2a r ((Sin[θ])^2))/((r^2)+(a^2)((Cos[θ])^2))),0,0,((((r^2) +(a^2))((r^2)+(a^2)((Cos[θ])^2))+2r (a^2)((Sin[θ])^2))/((r^2)+(a^2)((Cos[θ])^2)))((Sin[θ])^2)}};
ρ=(Sqrt[(r^2)+(a^2)((Cos[θ])^2)]);
Δ=(r^2)-2r+(a^2);
Σ = (Sqrt[((((r^2)+(a^2)))^2)-Δ (a^2)((Sin[θ])^2)]);
inversemetric=FullSimplify[Inverse[metric]];
affine=FullSimplify[Table[Sum[(inversemetric[[α,δ]])
(D[metric[[δ,β]],coord[[γ]]]+D[metric[[δ,γ]],coord[[β]]]-D[metric[[β,γ]],coord[[δ]]]),{δ,1,4}]/2,{α,1,4},{β,1,4},{γ,1,4}]];
rh[a_]:=1+(Sqrt[1-(a^2)])
Ωh[a_]:=(a/(2rh[a]))
Ωk[r_,a_]:=((((r^(3/2))+a))^(-1))
ek[r_,a_]:=((1-2(r^(-1))+a (r^(-3/2)))/(((1-3(r^(-1))+2a (r^(-3/2))))^(1/2)))
lk[r_,a_]:=(((1-2a (r^(-3/2))+(a^2)(r^(-2))) (r^(1/2)))/(((1-3(r^(-1))+2a (r^(-3/2))))^(1/2)))
rms[a_]:=Chop[
3-Sign[a]√(3+(a^2)+(-3+a) (((-1+a))^(1/3)) (((1+a))^(2/3))+(((-1+a))^(2/3)) (((1+a))^(1/3)) (3+a))+√(18-(-3+a) (((-1+a))^(1/3)) (((1+a))^(2/3))-(((-1+a))^(2/3)) (((1+a))^(1/3)) (3+a)+2 (-6+(a^2))-Sign[a](16 (a^2))/(√(3+(a^2)+(-3+a) (((-1+a))^(1/3)) (((1+a))^(2/3))+(((-1+a))^(2/3)) (((1+a))^(1/3)) (3+a))))]

Field Tensor

The field tensor can be expressed in terms of the vector potential
(Subscript[F, αβ])=2(Subscript[A, [β,α]])
which guarantees that the homogeneous field equations, 
(Subscript[F, [αβ,γ]])=0
are satisfied. The relevant gauge condition is Lorentz
(Subscript[(A^α), ;α])=0.

Field Equations

The four field equations can be written (in Heaviside units) as
(Subscript[(F^(αβ)), ;β])=(J^α)
(Subscript[F, (αβ,γ)])=0

Vacuum Solutions

Axisymmetric, Schwarzschild electrostatics

Vacuum solutions are instructive though unrealistic as they will quickly discharge. Restrict attention to a charge at rest on the pole of a Schwarzschild hole will have a field described by -(ξ^((t)))·A=-(Subscript[A, t])=V(r,θ), the electrostatic potential.   This trivially satisfies the gauge condition and the homogeneous Maxwell equations. The one inhomogeneous equation away from the charge is (Subscript[[(((-g))^(1/2))(F^(tα))], ,α])=0, as (F^(αβ)) is antisymmetric.

coord={t,r,θ,ϕ};metric={{-1+2/r,0,0,0},{0,1/(1-2/r),0,0},{0,0,(r^2),0},{0,0,0,(r^2)  ((Sin[θ])^2)}};
inversemetric:=inversemetric=Simplify[Inverse[metric]];
fFcov={{0,-D[fV[r,θ],r],-D[fV[r,θ],θ],0},{D[fV[r,θ],r],0,0,0},{D[fV[r,θ],θ],0,0,0},{0,0,0,0}};
fFcon=Table[Sum[inversemetric[[α,γ]]fFcov[[γ,δ]]inversemetric[[δ,β]],{γ,1,4},{δ,1,4}],{α,1,4},{β,1,4}];
MatrixForm[fFcov]
MatrixForm[fFcon]
Simplify[Table[(1-2/r)Csc[θ]Sum[D[(r^(2 ))Sin[θ] fFcon[[α,β]],coord[[β]]], {β,1,4}],{α,1,4}]]

From this, we can read off r(r-2)(Subscript[V, ,rr])+2(r-2)(Subscript[V, ,r])+(Subscript[V, ,θθ])+cotθ (Subscript[V, ,θ])=0. Separating variables as usual gives V(r,θ)=R(r)(Subscript[P, l])(cosθ), where r(r-2)(Subscript[R, ,rr])+2(r-2)(Subscript[R, ,r])-l(l+1)R=0. There is a polynomial series solution (Subscript[R, 1l]) = (∑_(n=0)^l)(Subscript[c, n/])(r^n) where (Subscript[c, n+1l])=[((n(n+1)-l(l+1))/(2(n+1)(n+2)))](Subscript[c, nl]). If we normalize to (Subscript[R, 1l])(0)=(((-1))^l), then (Subscript[R, 10]) = 1, (Subscript[R, 11])=-1+(r/2), (Subscript[R, 12])=1-((3 r)/2)+((r^2)/2)... Note that (Subscript[R, 1l])(2)=0. l≥1.
The second solution can be generated in the usual way as (Subscript[R, 2l])=(Subscript[R, 1l])∫(dr/((r^2)((Subscript[R, 1l])^2))). The first three solutions normalized so that (Subscript[R, 2l])(2)=1 are (Subscript[R, 20])=(2/r),(Subscript[R, 21])=2 (1-(r/2)) ((1/(2-r))-(1/r)-Log[-2+r]+Log[r]),(Subscript[R, 22])=-2 (1-((3 r)/2)+((r^2)/2)) (-(1/r)+((9-5 r)/(2-3 r+(r^2)))-3 Log[-2+r]+3 Log[r]). We also need the Wronskians (Subscript[W, l])=(Subscript[R, 1l])(Subscript[R, 2l]^('))-(Subscript[R, 1l]^('))(Subscript[R, 2l])=-(2/(r^2)).

lmax=6;fR1=Table[1,{l,0,lmax}];
Do[tm=(((-1))^l);fR1[[l+1]]=tm;Do[tm=r((n(n+1)-l(l+1))/(2(n+1)(n+2)))tm ;fR1[[l+1]]=fR1[[l+1]]+tm,{n,0,l}],{l,0,lmax}];
fR2=Table[fR1[[i]] Integrate[((-2 )/((r^2)((fR1[[i]])^2))),r],{i,1,lmax+1}];
Limit[fR2,r->2];Simplify[Table[fR1[[i]]D[fR2[[i]],r]-D[fR1[[i]],r]fR2[[i]],{i,1,lmax+1}]]

Place a unit charge on the axis at height s and assume no charge inside the hole. For 2≤r≤s, V=(∑_l)(Subscript[A, l])(Subscript[R, 1l])(Subscript[P, l]) and for s≤r, V=(∑_l)(Subscript[B, l])(Subscript[R, 2l])(Subscript[P, l]). Clearly (Subscript[B, 0])=(1/2). Continuity implies that (Subscript[A, l])(Subscript[R, 1l])(s)=(Subscript[B, l])(Subscript[R, 2l])(s) and (∑_(l=0)^∞)((Subscript[A, l])(Subscript[R, 1l]^('))(s)-(Subscript[B, l])(Subscript[R, 2l]^('))(s))(Subscript[P, l])∝δ(1-μ). Hence (Subscript[A, l])(Subscript[R, 1l]^('))(s)-(Subscript[B, l])(Subscript[R, 2l]^('))(s)∝2l+1 and (Subscript[A, l])=(((Subscript[R, 2l])(s))/((Subscript[R, 1l])(s)))(Subscript[B, l])=(l+(1/2))(Subscript[R, 2l])(s). 
This approach leads to a Greens function for any axisymmetric distribution of charge including inside the hole.

s=6.;
cA=Chop[Table[(l+(1/2))fR2[[l+1]],{l,0,lmax}]/.r->s];cB=Chop[Table[(l+(1/2))fR1[[l+1]],{l,0,lmax}]/.r->s];
pltin=Plot[Table[cA[[l+1]]fR1[[l+1]],{l,0,lmax}],{r,2,s}];
pltout=Plot[Table[cB[[l+1]]fR2[[l+1]],{l,0,lmax}],{r,s,10}];
Show[pltin, pltout,Frame->True,PlotRange->{{2,10},{0,1.2/s}}]
fVin[r1_,θ_]:=Sum[cA[[i]] fR1[[i]] LegendreP[i-1,Cos[θ]],{i,1,lmax+1}]/.r->r1
fVout[r1_,θ_]:=Sum[cB[[i]] fR2[[i]] LegendreP[i-1,Cos[θ]],{i,1,lmax+1}]/.r->r1
Plot[If[r<s,fVin[r,θ],fVout[r,θ]]/.θ->1,{r,2,8}];
ContourPlot[{fVin[((((x^2)+(y^2)))^(1/2)),ArcTan[y,x]],fVout[((((x^2)+(y^2)))^(1/2)),ArcTan[y,x]]},{x,.01,10.},{y,.01,10.},Contours->{.002},AspectRatio->Automatic]

Following Copson(1928), there is actually a closed form version of this expansion V=(1/(s r))(1+(((r-1)(s-1)-μ)/((((((r-1))^2)+(((s-1))^2)-1-2(r-1)(s-1)μ+(μ^2)))^(1/2)))) which satisfies the inhomogeneous equation, has the constant value of (1/s)on the horizon and ~ (1/r) as r->∞. We see the horizon behaving as a conductor.

fV[r_,θ_]:=(1/(s r))(1+(((r-1)(s-1)-Cos[θ])/( ((((((r-1))^2)+(((s-1))^2)-1-2(r-1)(s-1)Cos[θ]+((Cos[θ])^2)))^(1/2)))))
Simplify[r(r-2)D[fV[r,θ],{r,2}]+2(r-2)D[fV[r,θ],r]+D[fV[r,θ],{θ,2}]+Cot[θ] D[fV[r,θ],θ]];
s=6.;
Show[ContourPlot[fV[((((x^2)+(y^2)))^(1/2)),ArcTan[y,x]],{x,.01,10.},{y,.01,10.},AspectRatio->Automatic,Contours->(1/s){.6,.8,1,1.5,2}],Graphics[Disk[{0,0},2,{0,.5π}]]]
Clear[s];

Axisymmetric Schwarzschild  magnetostatics

The second scalar is 2π(ξ^((ϕ)))·A=2(Subscript[πA, ϕ])≡W which can be identified with the enclosed magnetic flux.

coord={t,r,θ,ϕ};metric={{-1+2/r,0,0,0},{0,1/(1-2/r),0,0},{0,0,(r^2),0},{0,0,0,(r^2)  ((Sin[θ])^2)}};
inversemetric:=inversemetric=Simplify[Inverse[metric]];
fFcov={{0,0,0,0},{0,0,0,D[fΦ[r,θ],r]},{0,0,0,D[fΦ[r,θ],θ]},{0,-D[fΦ[r,θ],r],-D[fΦ[r,θ],θ],0}};
fFcon=Table[Sum[inversemetric[[α,γ]]fFcov[[γ,δ]]inversemetric[[δ,β]],{γ,1,4},{δ,1,4}],{α,1,4},{β,1,4}];
MatrixForm[fFcov]
MatrixForm[fFcon]
Simplify[-(r^2)Sin[θ]Table[Sum[D[(r^(2 ))Sin[θ] fFcon[[α,β]],coord[[β]]], {β,1,4}],{α,1,4}]]

From this, we can read off r(r-2)(Subscript[W, ,rr])+2(Subscript[W, ,r])+(Subscript[W, ,θθ])-cot θ (Subscript[W, ,θ])=0. Again we can separate variables and solve for arbitrary axisymmetric current loops. Alternatively we can obtain simple analytical solutions: 
W=(Subscript[πB, ∞])(r^2)(sin^2)θ, corresponding to a uniform field (Subscript[B, ∞]) as r->∞
W = (Subscript[W, 0])(1-cosθ) corresponding to a split monopole field with total flux (Subscript[W, 0])threading the hole and toroidal current in the equatorial plane.
W = ((Subscript[W, 0])/(4ln2))[(r-2)(1-cosθ)+4ln2-2(1+cosθ)ln(1+cosθ)] corresponding to a paraboloidal field with total flux (Subscript[W, 0])threading the hole and toridal current in the equatorial plane.

fW[r_,θ_]:=(((r-2)(1-Cos[θ])+4Log[2]-2(1+Cos[θ])Log[1+Cos[θ]])/(4Log[2]))
Simplify[r(r-2)D[fW[r,θ],{r,2}]+2D[fW[r,θ],r]+D[fW[r,θ],{θ,2}]-Cot[θ] D[fW[r,θ],θ]];
Show[ContourPlot[fW[((((x^2)+(y^2)))^(1/2)),ArcTan[y,x]],{x,.001,4.},{y,.001,6.},AspectRatio->Automatic,Contours->{.2,.4,.6,.8,1.,1.2,1.4,1.6}],Graphics[Disk[{0,0},2,{0,.5π}]]]

Force Free Electrodynamics

Introduction

In order to demonstrate how power may be extracted from a black hole, we allow for current and spin.

Force-Free Equation

In the simplest treatment, the plasma which carries the current is assumed to have negligible inertia so that
(Subscript[F, αβ])(J^β)=0. Equivalently, the (Maxwell) electromagnetic stress tensor (T^EM)has vanishing divergence.

Stationary, Axisymmetric Fields

We assume that fields as well as the spacetime are stationary and axisymmetric so that the enclosed magnetic flux W is constant on a set of nested surfaces generated by the local magnetic field. The t and ϕ components of the force-free equation give (J^i)(Subscript[V, ,i])=(J^i)(Subscript[W, ,i])=0, where i=r, θ, implying that the flux surfaces are also equipotentials, V=V(W).

Introduce the quantity 
(Subscript[B, ϕ ])≡Ψ=(Subscript[A, θ,r])-(Subscript[A, r,θ])
The field tensor is gauge invariant

Flux Surfaces

The t,ϕ force-free conditions imply
(Overscript[J, ->])·∇(Subscript[A, t])=(Overscript[J, ->])·∇V=(Overscript[J, ->])·∇(Subscript[A, ϕ])=(((Overscript[J, ->])·∇Φ)/(2π))=0⇒∇V=-(Ω/(2π))∇Φ⇒V=V(Φ), Ω=Ω(Φ).
The field equations give
(Subscript[[(((-g))^(1/2))(g^rr)(g^(θθ))Ψ], ,θ])=(Subscript[[(Subscript[ΔsinθB, ϕ])/(ρ^2)], ,θ])=(((-g))^(1/2))(J^r)=(Subscript[I, ,θ])/2π
(Subscript[[(((-g))^(1/2))(g^rr)(g^(θθ))Ψ], ,r])=(Subscript[[(Subscript[ΔsinθB, ϕ])/(ρ^2)], ,r])=-(((-g))^(1/2))(J^θ)=(Subscript[I, ,r])/2π
Hence 
I=I(Φ). 
Flux surfaces are equipotential, isorotational and contain the conserved current I measured at infinity.

fieldtensor={{0,Ω[Φ[r,θ]] D[Φ[r,θ],r],Ω[Φ[r,θ]] D[Φ[r,θ],θ],0},{-Ω[Φ[r,θ]] D[Φ[r,θ],r],0,(((ρ^2)Ψ[r,θ])/(Δ Sin[θ])), D[Φ[r,θ],r]},{-Ω[Φ[r,θ]] D[Φ[r,θ],θ],-(((ρ^2)Ψ[r,θ])/(Δ Sin[θ])),0, D[Φ[r,θ],θ]},{0,- D[Φ[r,θ],r],-D[Φ[r,θ],θ],0}};
mixedfieldtensor=FullSimplify[Table[Sum[inversemetric[[α,γ]]fieldtensor[[γ,β]],{γ,1,4}],{α,1,4},{β,1,4}]]

Kinematic Condition

The r or θ force-free conditions imply
(J^ϕ)=((ΩJ)^t)+(((J^r)Ψ)/(Subscript[A, ϕ,θ]))=((ΩJ)^t)+(1/(((-g))^(1/2)))(dI/(dΦ))Ψ=((ΩJ)^t)+(I/(2π Δ (sin^2)θ))(dI/(dΦ))
This identifies Ω as the kinematic angular velocity of the field lines.

Stress Balance Equation

The two remaining field equations furnish a nonlinear PDE in (Subscript[A, ϕ]) which imposes  a stress balance condition across the flux surfaces and fixes their shape.

Boundary conditions

In order to avoid the singularity at the horizon, we need to evaluate quantities in the free-falling frame. One way to impose this is to transform to Kerr-Schild coordinates. A more physical approach is to consider the situation from the from the vantage point of the ZAMO frame suspended just above the horizon. This observer will see a highly blueshifted transverse electromagnetic wave traveling radially inward which diverges as the horizon is approached and a diminishing outgoing wave which tends to zero in this limit. The radial fields are unaffected. In other words, (Subscript[B, r]),(Subscript[E, r ])are finite,  and (E^(OverHat[θ]))-(B^(OverHat[ϕ])),(E^(OverHat[ϕ]))+(B^(OverHat[θ]))diverge and (E^(OverHat[θ]))+(B^(OverHat[ϕ])),(E^(OverHat[ϕ]))-(B^(OverHat[θ]))→0. The boundary condition for our problem is that (E^(OverHat[θ]))→-(B^(OverHat[ϕ])). In BL coordinates this becomes 
Δ Ψ→(Ω-(Subscript[Ω, H]))(((SubPlus[r])^2)+(a^2))(Subscript[A, ϕ,θ]).
The finite transverse fields measured in an infalling frame will be equal to (Subscript[E, H]). 
The toroidal field diverges but in just such a way as to keep the current measured at infinity 
I=((2π sinθ Δ Ψ)/(ρ^2))
finite. It is helpful to think of this condition between the either the toroidal magnetic field and the poloidal electric field,

or the current that must ultimately cross the field lines beyond the horizon and the poloidal electric field as an Ohm's law
where the horizon is given a surface resistivity E / J of 377Ω. This is not a true resistance of the horizon - in an infalling frame, currents cross it rather than flow along it. However, it is a useful way to consider the problem just like the 
ascribing an impedance to free space is useful to a radio engineer, 
The boundary conditions at large radius depend upon the electromagnetic environment. If the electromagnetic field matches onto a relativistic outflow, as we hope, the same impedance condition will apply; if it matches onto a conducting plasma - eg an accretion disk - then the effective resistance will be small. We call this region the load.

Circuit analysis

These considerations motivate a circuit analysis. Consider two neighboring axisymmetric flux surfaces trapping a flux 
ΔΦ. We can decompose the magnetosphere into many nested surfaces.  The current flowing out and back along these
surfaces is simply I.  At the horizon, the renormalized potential difference between the two surfaces evaluated in the 
ZAMO frame is 
(Subscript[ΔV, H])=((((Subscript[Ω, H])-Ω)ΔΦ)/(2π))=(Subscript[IΔR, H])
Evaluating,
(Subscript[ΔR, H])=((60((SubPlus[ρ])^2)Δθ)/(Σ sin θ))Ω
A similar analysis at the load gives
(Subscript[ΔV, L])=(((Ω-(Subscript[Ω, L]))ΔΦ)/(2π))=(Subscript[IΔR, L])
Hence, the field line rotates with the compromise angular velocity
Ω=(((Subscript[Ω, H])(Subscript[ΔR, L])+(Subscript[Ω, L])(Subscript[ΔR, H]))/((Subscript[ΔR, L])+(Subscript[ΔR, H])))
Typically, the load and the horizon have matched impedance and as much power is dissipated behind the horizon, increasing its irreducible mass as is dissipated in the load. Note that the direction of the electric field for the ZAMO observer at the horizon is reversed from that for an observer at the load. There appears to be energy going into the horizon  and out to large radius. In order to understand what is happening, it is necessary to transform into BL coordinates where there are conserved fluxes of energy and angular momentum associated with the two Killing vectors. In this frame the direction of the electric field is reversed and energy appears to flow from the horizon. Of course the "source" of the power is really the spinning spacetime in the ergosphere and these frame dependent associations can mislead.
The conserved electromagnetic torque in the BL frame is 
ΔG=((IΔΦ)/(2π))=(((I^(2 ))((Subscript[ΔR, H])+(Subscript[ΔR, L])))/((Subscript[Ω, H])-(Subscript[Ω, L])))
Rearranging gives the energy budget 
(Subscript[ΔGΩ, H])=(Subscript[ΔGΩ, L])+(I^(2 ))(Subscript[ΔR, H])+(I^2)(Subscript[ΔR, L])

MHD vs Force-Free

This formalism has ignored the inertia of the plasma. A more general and more complex formalism includes the inertial terms in the relativistic MHD approximation.

Homework 5

1. Electrostatic Field in Schwarzschild Spacetime

Consider a spherically symmetric electric field in a Schwarzschild spacetime described by the single vector potential component (Subscript[A, 0])(r).Write down the source-free field equation and solve it for a central charge Q located on the black hole.  Interpret your solution carefully in terms of the electric field strength as measured in a locally orthonormal frame by an observer at rest with respect to Schwarzschild coordinates. The charge never crosses the horizon according to a distant observer. Explain what is going on.

2. Magnetostatic Field in Schwarzschild Spacetime

Consider a magnetostatic field described by (Subscript[A, ϕ])=(r/2)(1-cosθ)+(1+cosθ)(1-ln(1+cosθ)). Show that this satisfies the source-free field equations and plot the magnetic field lines in some representation of your choosing, again interpreting your diagram carefully.

3. Hole-Disk Connection

A black hole is surrounded by a thin accretion disk extending down to its marginally stable circular orbit.  Magnetic flux 
connects the disk surface to the event horizon. Show that the angular momentum transfer is only from the hole to the disk if the hole angular velocity exceeds 0.093.

9. Gravitational Radiation

Motivation

The discovery of the binary pulsar in 1974 led to a test of the theory of gravitational radiation for two neutron stars in a eight hour orbit with an eccentricity ~0.6. The theory for rate at which the period decreases is now corroborated to a few parts per thousand and several more such systems are known. However, we have still not yet succeeded in detecting gravitational radiation directly. On class of detectors, epitomized by LIGO, uses  Michelson interferometers to seek kHz sources. Advanced LIGO should see several sources beginning in 2014. Another proposed project is LISA, which is sensitive to mHz frequencies, is designed to be able to see merging modestly massive black holes at large redshift though a foregound of white dwarf binaries. A third approach for the μHz range is to use pulsar timing.

Summary of slow motion, weak field formalism

Waves are weak and we linearize. They are described with a metric perturbation in the far field
(Subscript[g, αβ])=(Subscript[η, αβ])+(Subscript[h, αβ])
in the Lorentz gauge, (Subscript[(h^(αβ)), ;β])=0, which is automatically satisfied as it is a contraction of the convariant derivative of the metric tensor. In a suitably oriented coordinate system, we   have (Subscript[h, xx])=-(Subscript[h, yy])=(SubPlus[h]); (Subscript[h, xy])=(Subscript[h, yx])=(Subscript[h, X])
If a wave passes a test mass at rest, (((d^2)x )/((dτ)^2)) =-(Subscript[(Γ^x), tt])=0 etc and x is const. However the proper distances between two test masses increases by a fraction h/2 - the strain.
The wave satisfies the linearized field equations (Subscript[h, xx,tt])-(Subscript[h, xx,zz])=16π ((Subscript[T, xx])^((TT))) etc.

Clear[coord,metric,inversemetric,affine,h];
coord={t,x,y,z};
metric={{-1,0,0,0},{0,1+h[t,z], 0,0},{0,0,1-h[t,z],0},{0,0,0,1}};
inversemetric=Inverse[metric];
affine=FullSimplify[Table[Sum[(inversemetric[[α,δ]])
(D[metric[[δ,β]],coord[[γ]]]+D[metric[[δ,γ]],coord[[β]]]-D[metric[[β,γ]],coord[[δ]]]),{δ,1,4}]/2,{α,1,4},{β,1,4},{γ,1,4}]];
affine[[2,1,1]]
riemann=FullSimplify[Table[D[affine[[α,β,δ]],coord[[γ]]]-D[affine[[α,β,γ]],coord[[δ]]]+Sum[affine[[α,ϵ,γ]]affine[[ϵ,β,δ]]-affine[[α,ϵ,δ]]affine[[ϵ,β,γ]],{ϵ,1,4}],{α,1,4},{β,1,4},{γ,1,4},{δ,1,4}]];
ricci=Simplify[Table[Sum[inversemetric[[α,δ]]riemann[[γ,δ,γ,β]],{γ,1,4},{δ,1,4}],{α,1,4},{β,1,4}]];
einstein=FullSimplify[ricci-(1/2)Tr[ricci]{{1,0,0,0},{0,1,0,0},{0,0,1,0},{0,0,0,1}}];
lineinstein=einstein/.{h[t,z]->0,Derivative[0,1][h][t,z]->0,Derivative[1,0][h][t,z]->0}

The solution of the wave equation is (h^ij)=(4/r)∫dV (T^ij)at the retarded time. Next combine the gauge with stress-energy conservation to obtain (Subscript[(T^tt), ,tt])=(Subscript[(T^kl), ,kl]), multiply by (x^i)(x^j), integrate over the volume of the source and twice by parts to obtain:
(T^ij)= (2/r) ∫(dVx^i)(x^j)(Subscript[(T^tt), ,tt])
Taking the transverse trace-free part of the leading term in the stress tensor gives
(h^ij)=(2/r)((Subscript[(δ^i), k])-(n^i)(Subscript[n, k]))((Subscript[(δ^j), l])-(n^j)(Subscript[n, l]))((Overscript[I, ..])^kl)
where the trace free inertia tensor is
(I^ij)=∫dm((x^i)(x^j)-((x^2)/3)(δ^ij)) 
The flux of gravitational radiation is
F=(1/(32π))<((OverDot[h])^ij)(Subscript[(OverDot[h]), ij])>=(1/(8((πr)^2)))<((Overscript[I, ...])^ij)(Subscript[(Overscript[I, ...]), ij])>
and the power is 
P=∫((dΩr)^2)F=(1/5)<((Overscript[I, ...])^ij)(Subscript[(Overscript[I, ...]), ij])>.
The torque is
G=∫((dΩr)^2)r xF=(2/5)(Subscript[ϵ, ijk])<((Overscript[I, ...])^jl)(Subscript[((Overscript[I, ..])^k), l])>.
It is possible to derive these using a radiation reaction potential 
ϕ=(1/5)(Subscript[(Overscript[I, .....]), ij])(x^i)(x^j).

Binary Stars

Circular Orbits

The power is computed to be 
P=((32(μ^2)(m^3))/(5(r^5)))=ΩG=(32/5)(μ^2)(m^(4/3))(Ω^(10/3))
nb the observed frequency is 2Ω. The multiplying constant is (c^5)/G=3.6 x (10^59)erg (s^(-1))or 2 x (10^5) solar masses per second.
A measure of the decay rate is the number of orbits left: N=((3(Ω^2))/(10π(OverDot[Ω])))=(1/(64π))(m/μ)(((Ωm))^(-5/3))

energy=-((μ m)/(2r));power=((32(μ^2)(m^(3 )))/(5(r^5)));
((-3(Ω^(2 ))D[energy/.r->(m^(1/3))(Ω^(-2/3)),Ω])/(10π power/.r->(m^(1/3))(Ω^(-2/3))))

Elliptical Orbits

Computing the orbit-averaged power with a the semi major axis, P=-(d/dt)(((-mμ)/(2a)))=((32(μ^2)(m^3))/(5(a^5)))(((1+(73/24)(e^2)+(37/96)(e^4)))/(((1-(e^2)))^(7/2))) and the torque is G=-(d/dt)[((μ(ma(1-(e^2))))^(1/2))] =((32(μ^2)(m^(5/2)))/(5(a^(7/2))))(((1+(7/8)(e^2)))/(((1-(e^2)))^2)).  
We can then obtain expressions for the rate of change of a, period T and e. 
 ((dln a)/dt)=(2/3)(dlnT/dt)= -((2 aP)/(μm)) = -((64((μm)^2))/(5(a^4)))(((1+(73/24)(e^2)+(37/96)(e^4)))/(((1-(e^2)))^(7/2))). 
 (dlne/dt)=(G/((μe)^2))(((((1-(e^2))/ma)))^(1/2))+((1-(e^2))/(2(e^2)))((dln a)/dt)= -((304((μm)^2)(1+(121/304)(e^2)))/(15(((a^4)(1-(e^2)))^(5/2))))

Parabolic Orbits

Take the limit with s=a(1-e) fixed, ΔE=((85((πμ)^2)(m^(5/2)))/((2^(5/2))3 (s^(7/2)))); ΔL=((6((πμ)^2)(m^2))/(s^2)). If we introduce  (Subscript[Ω, max])=(((2^(1/2))(m^(1/2)))/(s^(3/2))), ΔE=(85/144)(Subscript[ΔLΩ, max]). If there is enough energy radiated a star can become bound.

Extreme Mass Ratio Inspirals

Motivation

Stars are observed to orbit massive black holes including the one in the nucleus of our Galaxy and they can spiral in under gravitational radiation friction. We may have sees stars on unbound orbits being tidally disrupted. Here we will consider tidal disruption of bound stars on circular equatorial orbits. RX J1034+396 with a one hour period could be an example.

Angular Momentum

The relativistic corrections to circular orbits are significant. Relative to Kerr,

lgrplt[a_]:=Plot[(((1-2a (r^(-3/2))+(a^2)(r^(-2))) )/(((1-3(r^(-1))+2a (r^(-3/2))))^(1/2))),{r,rms[a],30},Frame->True,PlotRange->{{1,30},{1,1.5}},FrameLabel->{"r","L/(Subscript[L, K])"},RotateLabel->False]
Show[lgrplt[-1],lgrplt[0],lgrplt[.9],lgrplt[.99],lgrplt[.999]]

Tidal Potential

In units where lengths are measured in (μ^(1/3))r and the potential in units of (μ^(2/3))/r, we can write ϕ[x_,y_,z_]:=-((((x^2)+(y^2)+(z^2)))^(-1/2))-(((k (x^2)-(k-1) (y^2)- (z^2)))/2);where x,y,z  are in the r,θ,ϕ directions. Note that the tidal potential is Laplacian.

r=6.;a=.95;Print[MatrixForm[tidal]]
t11=-tidal[[1,1]];t22=tidal[[2,2]];t33=tidal[[3,3]];
ContourPlot[-((((x^2)+(y^2)))^(-1/2))-.5(t11  (x^2)-t22 (y^2)),{x,-2(t11^(-1/3)),2(t11^(-1/3))},{y,-2(t11^(-1/3)),2(t11^(-1/3))},Contours->{-2(t11^(1/3)),-1.5(t11^(1/3)),-(t11^(1/3))}]
ContourPlot[-((((x^2)+(z^2)))^(-1/2))-.5(t11  (x^2)-t33 (z^2)),{x,-2(t11^(-1/3)),2(t11^(-1/3))},{z,-2(t11^(-1/3)),2(t11^(-1/3))},Contours->{-2(t11^(1/3)),-1.5(t11^(1/3)),-(t11^(1/3))}]

Roche Volume

Clear[r,a];
k[r_,a_]:=((3 (a^2)-4 a (r^(1/2))+(r^2))/(2 a (r^(1/2))+(r-3) r))
(*xm=(((k[r,a]+1))^(-1/3));
ym=((-k[r,a] (((k[r,a]+1))^(1/3))+(((((k[r,a])^2)+(Sqrt[((k[r,a])^3) (2 k[r,a]+1)])))^(2/3)))/(k[r,a] (((((k[r,a])^2)+(Sqrt[((k[r,a])^3) (2 k[r,a]+1)])))^(1/3))));
zm=((-(((k[r,a]+1))^(1/3))+(((1+(Sqrt[(k[r,a]+2 )])))^(2/3)))/( (((1+(Sqrt[(k[r,a]+2 )])))^(1/3))));ϕm=-1.5(((k[r,a]+1))^(1/3));*)
v[r_,a_]:=Chop[8NIntegrate[((((s^2)-(x^2)-(y^2)))^(1/2))/.Solve[-1/s+.5( (s^2)-(k[r,a]+2)(x^2)+(k[r,a]-1)(y^2))==-1.5(((k[r,a]+1))^(1/3)),s][[3]],{x,0,(((k[r,a]+1))^(-1/3))},{y,0,((((s^2)-(x^2)))^(1/2))/.Solve[-1/s+.5(k[r,a] (s^2)-(2k[r,a]+1)(x^2))==-1.5(((k[r,a]+1))^(1/3)),s][[3]]}]]
v[6,.95]

Gravitational Radiation

alst={-.99,-.9,-.5, 0.0,0.2,0.5,0.8,0.9,0.99,0.999};Ωlst=((1-(((1-(a1^2)))^(1/2)))/(2a1))/.a1->alst+.000001;
gwlst={{1.000,1.240,1.233,1.197,1.143,1.114,1.053,0.9144,0.7895,0.4148,0.2022},
{1.001,1.239,1.232,1.196,1.142,1.114,1.053,0.9140,0.7894,0.4154,0.2032},
{1.002,1.238,1.231,1.196,1.141,1.113,1.052,0.9137,0.7894,0.4160,0.2041},
{1.005,1.235,1.228,1.193,1.139,1.111,1.050,0.9126,0.7891,0.4177,0.2069},
{1.01,1.231,1.224,1.189,1.135,1.107,1.047,0.9109,0.7887,0.4207,0.2116},
{1.02,1.222,1.215,1.181,1.127,1.100,1.041,0.9076,0.7880,0.4263,0.2208},
{1.05,1.198,1.192,1.159,1.108,1.081,1.025,0.8988,0.7867,0.4434,0.2473},
{1.1,1.165,1.159,1.128,1.080,1.055,1.002,0.8876,0.7859,0.4701,0.2881},
{1.2,1.115,1.110,1.082,1.039,1.017,0.9706,0.8726,0.7882,0.5182,0.3581},
{1.3,1.081,1.075,1.051,1.012,0.9913,0.9493,0.8638,0.7920,0.5587,0.4160},
{1.4,1.055,1.051,1.028,0.9919,0.9733,0.9348,0.8583,0.7960,0.5930,0.4648},
{1.7,1.011,1.007,0.9888,0.9591,0.9435,0.9119,0.8524,0.8075,0.6665,0.5723},
{2.0,0.9893,0.9862,0.9705,0.9448,0.9312,0.9034,0.8530,0.8171,0.7117,0.6411},
{2.5,0.9734,0.9709,0.9580,0.9363,0.9248,0.9012,0.8589,0.8302,0.7556,0.7089},
{3.0,0.9674,0.9653,0.9542,0.9352,0.9250,0.9040,0.8662,0.8415,0.7813,0.7469},
{4.0,0.9651,0.9634,0.9546,0.9391,0.9306,0.9129,0.8807,0.8597,0.8121,0.7882},
{5.0,0.9665,0.9651,0.9577,0.9448,0.9371,0.9216,0.8930,0.8742,0.8320,0.8118},
{6.0,0.9687,0.9675,0.9611,0.9490,0.9430,0.9291,0.9031,0.8858,0.8469,0.8286},
{7.0,0.9709,0.9699,0.9641,0.9533,0.9480,0.9354,0.9116,0.8955,0.8589,0.8416},
{8.0,0.9730,0.9720,0.9669,0.9588,0.9522,0.9407,0.9186,0.9036,0.8689,0.8524},
{9.0,0.9749,0.9740,0.9693,0.9607,0.9558,0.9452,0.9246,0.9105,0.8774,0.8616},
{10.0,0.9765,0.9757,0.9714,0.9616,0.9589,0.9491,0.9298,0.9164,0.8847,0.8695}};
gwdat=Table[{alst[[i]],gwlst[[1,i+1]],gwlst[[2,i+1]],gwlst[[3,i+1]],gwlst[[4,i+1]],gwlst[[5,i+1]],gwlst[[6,i+1]],gwlst[[7,i+1]],gwlst[[8,i+1]],gwlst[[9,i+1]],gwlst[[10,i+1]],gwlst[[11,i+1]],gwlst[[12,i+1]],gwlst[[13,i+1]],gwlst[[14,i+1]],gwlst[[15,i+1]],gwlst[[16,i+1]],gwlst[[17,i+1]],gwlst[[18,i+1]],gwlst[[19,i+1]],gwlst[[20,i+1]],gwlst[[21,i+1]],gwlst[[22,i+1]]},{i,1,10}];
raddat=Table[gwlst[[i,1]],{i,1,22}];
gwtab=Table[{raddat[[i]],gwdat[[j,1+i]]},{j,1,10},{i,1,22}];
gwplt=ListPlot[Table[gwtab[[i]],{i,1,10}],Frame->True,Joined->True,PlotRange->{{1,10},{0,1.5}},ImageSize->600];
Show[gwplt,Graphics[Line[{{1,1},{10,1}}]]]

Gravitational Radiation from Major Coalescence

Motivation

Perhaps the most exciting prospect of all is to measure gravitational radiation from comparable mass black hole coalescence. Here the expansion procedures break down and teh resulting wave forms depend upon the masses and spin vectors (8 numbers) as well as the observer inclination (two numbers). Assume an initially circular orbit - not necessarily true. The calculation requires numerical solution of the full dynamical field equations which is a formidable task. However, the comparison with observation, when it comes will provide the first comprehensive test of (classical) GTR. Computing the recoil is also important.

Numerical Evolution

It has taken 35 years to solce this problem excpet in limiting cases. Features of the solution include :
-choice of coordinates for stability (harmonic) 
-economical imposition of constraint conditions
-compactify r->∞
-AMR
-excision of holes as they move through the grid
-active damping of high frequency instability in evolution and initial value formulation 
-convergence tests