Hi Daniel,
Let me first note that, since the Boyer-Lindquist coordinates are adapted to the spherical symmetry of Schwarzschild spacetime, there are at least as many Boyer-Lindquist charts on Schwarzschild spacetime as there are spatial rotations: performing a spatial rotation on the coordinates of a given Boyer-Lindquist chart preserves the form of the metric, so that the new system of coordinates should also be called a Boyer-Lindquist chart.
Now, as you say, a system of geodesic equations requires eight initial conditions to be well-posed. Four of them are the coordinates (t_0, r_0, th_0, ph_0) of the starting point P of the geodesic in a given Boyer-Lindquist chart, and the other four are the components (Dt_0, Dr_0, Dth_0, Dph_0) (with respect to the coordinate basis associated with our Boyer-Lindquist chart) of the vector V tangent to the geodesic at P with respect to the parameter s that is used to parametrise the geodesic. (More precisely, in SageManifolds, the geodesic equations that are solved numerically when you call method « solve » correspond to the « canonical » form of the geodesic equation that we are the most used to, i.e. without any right-hand side that would reveal a non-affine parametrization of the geodesic. Therefore, the parameter s is automatically enforced to be affine.)
Now, on one hand, one can actually prove that the geodesics in Schwarzschild are planar, i.e. the spatial coordinates of our geodesic (which is the unique affinely parametrized geodesic starting at P with initial tangent vector V) eternally remain in the unique plane of R^3 containing the vectors (r_0, th_0, ph_0) and (Dr_0, Dth_0, Dph_0).
On the other hand, the coordinates (t_0, r_0, th_0, ph_0) of P and the components (Dt_0, Dr_0, Dth_0, Dph_0) of V could be anything so far. Yet, one can always perform a spatial rotation from our original Boyer-Lindquist chart to another Boyer-Lindquist chart in which the spatial parts of P and V are contained in the equatorial plane of the new Boyer-Lindquist chart, i.e. th_0’ = pi/2 and Dth_0’ = 0. Then, let us use this new Boyer-Lindquist chart (and immediately abandon the primes although we have changed coordinates), since it is adapted to the initial conditions, at least, of our geodesic. Actually, using the previous paragraph, the spatial part of our geodesic will necessarily remain in the equatorial plane of the new Boyer-Lindquist chart, so that this chart is well adapted to the whole geodesic: for any value of the affine parameter s, theta is pi/2 and the theta-component of the tangent vector is 0.
This answers your last question, if I understood it well: there can be no intrisically 3D figure in this notebook since spatial parts of the Schwarzschild geodesics are planar. Setting the initial coordinates out of the equatorial plane would render the same results as in the whole notebook, only tilted.
Regarding your first question, the goal in this part of the notebook is to compute a spatially bounded, timelike geodesic. It is thus necessary and sufficient for the initial tangent vector V to be timelike, since it is a general property of geodesics that the squared norm of their tangent vectors is a constant of motion, so that the causal type of the geodesic is preserved throughout motion. Then, we want g(V,V) < 0 at P. In addition, we can actually try to get g(V,V) = -1 at P rather that any strictly negative value, since this would make the affine parameter s to be the proper time along the geodesic.
Now, the actual computation of g(V,V) at P requires to know the values of the metric components at P (with respect to the tensorial basis associated with the Boyer-Lindquist chart…). To do so, one first needs to prescribe the actual coordinates (t_0, r_0, th_0, ph_0) of P, before injecting them into the expression of the Schwarzschild metric. Since Schwarzschild is static and spherically symmetric, only the value of r_0 matters: as the horizon sits at r = 2 (since m is set to 1), r_0 = 8 seems fair.
But the important part is actually the prescription of the components of V : these play a more important role than r_0 in determining the aspects of the geodesic. These are chosen in order for the relation g(V,V) = -1 to hold. Now that the coordinates of P are set, we know that the (diagonal) metric components at P are the following:
g_00 = - (1 - 2m/r_0) = - (1 - 2/8) = -3/4
g_11 = (1 - 2m/r_0)^(-1) = 4/3
g_22 = r_0^2 = 8^2 = 64
g_33 = r_0^2 sin^2(th_0) = 8^2 sin^2(pi/2) = 64
Then :
g(V,V) = -3/4 Dt_0^2 + 4/3 Dr_0^2 + 64 Dth_0^2 + 64 Dph_0^2
For some reason (which is probably making the computations to come easier), I chose not to give any initial radial kick to the particle, i.e. Dr_0 = 0. Since we already have Dth_0 = 0, and we want g(V,V) = -1, we are left with:
-1 = -3/4 Dt_0^2 + 64 Dph_0^2
One is then free to chose any Dt_0 (or any Dph_0) and compute Dph_0 (or Dt_0) accordingly for the relation above to hold. For another reason, it seems that, after a few tests, I picked Dph_0 = 4,1/64, which does impose (requiring Dt_0 to be positive) Dt_0 = sqrt(80,81)/(4*sqrt(3)) as it appears in the notebook. The reason for choosing Dph_0 = 4,1/64 is probably that it seemed to produce nice visual results, allowing to see very clearly the advance of the periastron when one plots the spatial projection of the geodesic (figure below instruction [19] in the notebook).
Finally, Dt_0 is indeed the velocity of the Boyer-Lindquist time coordinate with respect to the affine parameter s/proper time along the geodesic. In physical words, consider a Boyer-Lindquist observer, i.e. someone who affects four numbers to any event he witnesses, those four numbers being the coordinates of the Boyer-Lindquist chart. If this observer sits at infinity, the Boyer-Lindquist time coordinate t corresponds to his own proper time.
If the particle sends a pulse every microsecond of its proper time s, the Boyer-Lindquist time coordinates t of the events corresponding to the emissions of the pulses will not be separated by a microsecond, due to relative motion and non-uniform gravitational field. Therefore, if you plot the Boyer-Lindquist time coordinate of the emission with respect to the corresponding proper time of the particle for each pulse, you get a function s -> t(s) different from the identity function s -> s, so that its derivative is different from the constant function equal to 1. Dt_0 is the value of this derivative when s is 0.
I might have gone into more details than you needed, but I really hope this answers your questions.
Karim
The 3d plots display correctly for me. Which browser are you using?
Microsoft Edge.
I am using Windows 10 on a 64bit machine with SageMath native.
Just to clarify, I can't see the 3D plots in the Jupyter nbviewer rendition of the file. If I download the ipynb file to my machine, change the 3D viewer to 'tachyon' then I can see the 3D plots.
I'm using Microsoft Edge to open Jupyter notebooks with no problems.
As I reported a while ago, the only 3D viewer that works on Windows 10 is 'tachyon'.
Daniel