Understading geodesic equation/code for (+1,-1,-1,-1) Sch. Metric

[Arman777]

Hey guys, as you may remember, I have posted a question about plotting the orbit of timelike and null-like particles for a given metric.

I think this discussion might be helpful for me and some other people in future studies. I have found an article, and in that article, the authors are using Sagemath to run the code. Luckily to run the Sagemath, you can use an online editor called CoCalc (https://cocalc.com/).

To run a Sagemath code, open a notebook and choose kernel as SageMath 9.3

Here is the code given by the authors.

reset()

# Define 4-dim. the manifold "Man":
Man = Manifold(4, 'Man', r'\mathcal{M}')

# Define the parameter "M" (mass):
M = var('M')

# Define the coordinates {t=0, r=1, theta(=th)=2, phi=3} with ranges
# (BL = Boyer-Lidquist)
BL.<t,r,th,ph> = Man.chart(r't r:(0,+oo) th:(0,pi):\theta ph:(0,2*pi):\phi')

# Define the metric "g" on manifold "Man":
g = Man.lorentzian_metric('g')

# Enter the Schwarzschild metric components:
g[0,0] = (1-(2*M)/r)
g[1,1] = -1/(1-(2*M)/r)
g[2,2] = -r^2
g[3,3] = -(r*sin(th))^2

# Display the metric
show('The Schwarzschild metric:')
show(g.display())

# Define variables, functions and calculate inverse metric
var('eta,m,E,L,S,HJfull')
F=function('F')(r)
G=function('G')(th)
ginv = g.inverse()

# Define the principal function Ansatz
S=((eta*m^2)/2)-E*t+L*ph+F+G

# Calculate the Hamilton-Jacobi equation
HJfull=0
for i in range(len(BL[:])):
    for j in range(len(BL[:])):
        HJfull=HJfull+ginv[i,j].expr()*diff(S,BL[i])*diff(S,BL[j])
HJfull=(diff(S,eta)-(1/2)*HJfull)
show('The Full Hamilton-Jacobi equation (variable name is HJfull):')
show(HJfull)

show('The geodesic equations in LaTeX form (variable name is geodeqnrhs):')
geodeqnrhs=zero_vector(SR, len(BL[:]))

for mu in range(len(BL[:])):
    for nu in range(len(BL[:])):
        geodeqnrhs[mu]=geodeqnrhs[mu]-(ginv[mu,nu].expr())*diff(S,BL[nu])
    writeresult='D[0](%s)($\eta$) = $%s$' %(BL[mu],latex(geodeqnrhs[mu]))
    show(writeresult)

show('Right hand sides of the geodesic equations as a vector')
show(geodeqnrhs)

# Take theta = pi/2:
var('HJst')

# Hamilton-Jacobi equation for theta=pi/2 (variable name is HJst)
HJst=(HJfull.subs(diff(G)==0)).subs(th=pi/2)
show(HJst)

# Right hand sides of the geodesic equations for theta=pi/2
# (variable name is geodeqnrhsst)
geodeqnrhsst=(geodeqnrhs.subs(diff(G)==0)).subs(th=pi/2)
show(geodeqnrhsst)

# Importing the solver
from sage.calculus.desolvers import desolve_odeint
#Importing circle for visualizing the black hole
from sage.plot.circle import Circle

var('m_aux,L_aux,E_aux,M_aux,r_initial,ph_initial,eta_end,step_size')

##########################
# Variables
m_aux=1 # Either 1 (timelike) or 0 (null).
M_aux=1
L_aux=4
E_aux=1.0
step_size=0.1
eta_end=500
r_initial=2.1*M_aux
ph_initial=0.3
##########################

# Derivative of F (the radial function)
# (variable name is derofradfun)
# We will use the first root
derofradfun=solve(HJst,diff(F,r))

# Define equations to solve
# dr/dt = geodeqn1
# dphi/dt = geodeqn2
geodeqn1=((geodeqnrhsst[1]/geodeqnrhsst[0]).subs(diff(F,r)==derofradfun[1].rhs())).subs(E=E_aux,L=L_aux,m=m_aux,M=M_aux)
geodeqn2=(geodeqnrhsst[3]/geodeqnrhsst[0]).subs(E=E_aux,L=L_aux,m=m_aux,M=M_aux)

# Solve the equations
sol=desolve_odeint([geodeqn1,geodeqn2],[r_initial,ph_initial],srange(0,eta_end,step_size),[r,ph])
p=line(zip(sol[:,0]*cos(sol[:,1]),sol[:,0]*sin(sol[:,1])))

# Plot the black hole as a circle
# Show the geodesics and the circle on the same plot
C=circle((0,0),2*M_aux,fill=True,rgbcolor='black')
show(C+p)

If you be careful, you'll see that the metric defined as $diag(1,-1,-1,-1)$, and the variables are given in lines 75-85. I have tried to plot the images that I have provided in this post (https://www.physicsforums.com/threa...ulate-orbits-in-schwarzschild-metric.1003484/)
but for some reason, I could not produce the pictures.

In the code given above the $L_{aux}$ (the angular momentum) is given by $4$ but in my code it will be $L_{aux} = 4.3$.

So the problem is, when I change the energy levels to see how the code/program reacts, I am getting either strange results and sometimes errors. Do you guys have an idea how can make the code work to produce those images?

In general in many books the metric is given by $diag(-1,1,1,1)$. When you change the metric in lines $16-23$, you guys should also change the equation in $33$ (I guess).

Here is the article: https://arxiv.org/abs/1703.09738

Look Section 4.3 in page 15

Thanks

 

[Ambitwistor]

for sharing the code and the article. It seems like you are trying to plot the geodesics of particles in the Schwarzschild spacetime. The code you have shared looks correct, and it is using the Schwarzschild metric in the Boyer-Lindquist coordinates. However, I noticed a few things that might be causing issues with your results.

Firstly, in line 35, the variable "M" is not defined. It should be defined as the parameter "M" which represents the mass of the black hole. Also, in line 36, the variable "eta" is not defined. It should be defined as the affine parameter along the geodesic.

Secondly, in line 84, the variable "L_aux" is defined as 4.3, but in the article, it is mentioned as 4. It is possible that this discrepancy is causing issues with your results.

Thirdly, in line 87, the variable "E_aux" is defined as 1.0, but in the article, it is mentioned as 1. It is possible that this discrepancy is also causing issues with your results.

Lastly, in line 98, the variable "M_aux" is defined as 1, but in the article, it is mentioned as 1.5. It is possible that this discrepancy is also causing issues with your results.

I would suggest checking these variables and values to see if they match with the ones mentioned in the article. Also, make sure to define all the necessary variables before running the code. If you are still having issues, please provide more details about the errors or strange results you are getting. I would be happy to help further.