Register to reply

Light ray paths near schwarzschild blackhole

by nocks
Tags: blackhole, light, paths, schwarzschild
Share this thread:
Nabeshin
#19
Nov4-09, 04:11 PM
Sci Advisor
Nabeshin's Avatar
P: 2,193
W(r) should be the effective potential of the system.
nocks
#20
Nov20-09, 02:54 PM
P: 24
Could anyone expand on this please? I would appreciate the help.
I have the the effective potential in the schwarzschild metric as (L being angular momentum)
[tex]V_{eff} = ( 1 - \frac{r_{s}}{r})(mc^{2} + \frac{L^{2}}{mr^{2}})[/tex]

Would this be enough information to solve for r and [tex]\Phi[/tex] so that I could plot the trajectories.

Also could I use the same equation for plotting the course of an observer descending into the black hole?
bcrowell
#21
Nov20-09, 03:09 PM
Emeritus
Sci Advisor
PF Gold
bcrowell's Avatar
P: 5,597
Here is some python code I wrote that does pretty much what you're talking about. It's meant to be short and easy to understand, so it uses a pretty crude method of doing the numerical integration. If you want to do accurate numerical calculations of geodesics, you'd want to substitute a better integration method. There are various general-purpose subroutines out there, e.g., in the book Numerical Recipes in C. What my code does is to calculate the deflection of a light ray that grazes the sun. Actually, it calculates half he deflection for a ray that grazes the sun, with the mass of the sun scaled up by a factor of 1000 in order to keep the result from being overwhelmed by rounding errors in my el-cheapo integration method.

import math

# constants, in SI units:
G = 6.67e-11         # gravitational constant
c = 3.00e8           # speed of light
m_kg = 1.99e30       # mass of sun
r_m = 6.96e8         # radius of sun

# From now on, all calculations are in units of the
# radius of the sun.

# mass of sun, in units of the radius of the sun:
m_sun = (G/c**2)*(m_kg/r_m)
m = 1000.*m_sun

# Start at point of closest approach.
# initial position:
t=0
r=1 # closest approach, grazing the sun's surface
phi=-math.pi/2
# initial derivatives of coordinates w.r.t. lambda
vr = 0
vt = 1
vphi = math.sqrt((1.-2.*m/r)/r**2)*vt # gives ds=0, lightlike

l = 0    # affine parameter lambda
l_max = 20000.
epsilon = 1e-6 # controls how fast lambda varies
while l<l_max:
  dl = epsilon*(1.+r**2) # giant steps when farther out
  l = l+dl
  # Christoffel symbols:
  Gttr = m/(r**2-2*m*r)
  Grtt = m/r**2-2*m**2/r**3
  Grrr = -m/(r**2-2*m*r)
  Grphiphi = -r+2*m
  Gphirphi = 1/r
  # second derivatives:
  #  The factors of 2 are because we have, e.g., G^a_{bc}=G^a_{cb}
  at   = -2.*Gttr*vt*vr
  ar   = -(Grtt*vt*vt + Grrr*vr*vr + Grphiphi*vphi*vphi)
  aphi = -2.*Gphirphi*vr*vphi
  # update velocity:
  vt = vt + dl*at
  vr = vr + dl*ar
  vphi = vphi + dl*aphi
  # update position:
  r = r + vr*dl
  t = t + vt*dl
  phi = phi + vphi*dl

# Direction of propagation, approximated in asymptotically flat coords.
# First, differentiate (x,y)=(r cos phi,r sin phi) to get vx and vy:
vx = vr*math.cos(phi)-r*math.sin(phi)*vphi
vy = vr*math.sin(phi)+r*math.cos(phi)*vphi
prop = math.atan2(vy,vx) # inverse tan of vy/vx, in the proper quadrant
prop_sec = prop*180.*3600/math.pi
print "final direction of propagation = %6.2f arc-seconds" % prop_sec
nocks
#22
Nov20-09, 03:15 PM
P: 24
I actually have the numerical recipes book next to me although I may avoid solving the elliptic integral I mentioned earlier, and just use the approximation for light deflection, i.e. 4GM/bc[tex]^{2}[/tex], to get the einstein ring effect, and focus on the trajectory of the observer descending into the black hole.
bcrowell
#23
Nov20-09, 04:44 PM
Emeritus
Sci Advisor
PF Gold
bcrowell's Avatar
P: 5,597
Quote Quote by nocks View Post
I actually have the numerical recipes book next to me although I may avoid solving the elliptic integral I mentioned earlier, and just use the approximation for light deflection, i.e. 4GM/bc[tex]^{2}[/tex], to get the einstein ring effect, and focus on the trajectory of the observer descending into the black hole.
Keep in mind that 4GM/bc2 is only a weak-field approximation. It won't give you the right answer if you're close to the black hole.
George Jones
#24
Nov21-09, 04:16 AM
Mentor
George Jones's Avatar
P: 6,239
Quote Quote by George Jones View Post
Yes, exactly. The geodesic equation when [itex]\theta = \pi /2[/itex] then is

[tex]\begin{equation*}
\begin{split}
\frac{d \phi}{d \lambda} &= \frac{L}{r^2} \\
\left( \frac{dr}{d \lambda} \right)^2 &= E^2 - L^2 W \left( r \left( \lambda \right) \right), \\
\end{split}
\end{equation*}
[/tex]

where [itex]W \left( r \right)[/itex] is function that I'll specify later, and [itex]E[/itex] and [itex]L[/itex] are constants of motion.

In its present form, the second equation is a little difficult to implement on a computer since sometimes a positive square root is needed (increasing [itex]r[/itex]) and sometimes a negative square root is needed (decreasing [itex]r[/itex]) in the same photon orbit. To get aorund this, differentiate the second equation with respect to the affine parameter [itex]\lambda[/itex] taking into account that the [itex]r[/itex] in [itex]W \left( r \right)[/itex] is itself a function of [itex]\lambda[/itex].

What do you get for the second equation after this differentiation?
Quote Quote by nocks View Post
Given it's been a while since i've done any differentiation, is it simply:

[tex]\frac{dr}{d \lambda} \right) &= \frac{ E - L}{W \left( r \right)}[/tex]
No, this isn't quite right. Let's start with the left side. What is

[tex]\frac{d}{d\lambda} \left[ \left( \frac{dr}{d \lambda} \right)^2 \right] ?[/tex]
nocks
#25
Nov21-09, 02:49 PM
P: 24
Quote Quote by George Jones View Post
...
Given

[tex]\left( \frac{dr}{d \lambda} \right)^2 &= E^2 - V^2(r) \right)[/tex]
and
[tex]V^2(r) = \left(1 - \frac{2M}{r} \right)\frac{L^2}{r^2}[/tex]

I have [tex] \frac{d^2r}{d\lambda^2} = -\frac{1}{2}\frac{d}{dr}V^2(r)[/tex]
nocks
#26
Nov22-09, 03:01 PM
P: 24
Quote Quote by nocks View Post
[tex]V^2(r) = \left(1 - \frac{2M}{r} \right)\frac{L^2}{r^2}[/tex]
Guess I should state that

[tex]\frac{d}{dr}V^2(r) = \frac{2L^2M}{r^4} - \frac{2L^2(\left 1 - \frac{2M}{r} \right)}{r^3}[/tex]
George Jones
#27
Nov23-09, 07:16 AM
Mentor
George Jones's Avatar
P: 6,239
After simplifying the right side of
Quote Quote by nocks View Post
[tex]\frac{d}{dr}V^2(r) = \frac{2L^2M}{r^4} - \frac{2L^2(\left 1 - \frac{2M}{r} \right)}{r^3}[/tex]
and plugging back into
Quote Quote by nocks View Post
[tex] \frac{d^2r}{d\lambda^2} = -\frac{1}{2}\frac{d}{dr}V^2(r)[/tex]
what do you get?
nocks
#28
Nov23-09, 11:47 AM
P: 24
That would be

[tex]\frac{d^2r}{d\lambda^2} = - \frac{L^2(3M-r)}{r^4}[/tex]

Which gives me the closest radius for a stable orbit of a photon as 3M
nocks
#29
Nov23-09, 08:52 PM
P: 24
Quote Quote by nocks View Post
That would be

[tex]\frac{d^2r}{d\lambda^2} = - \frac{L^2(3M-r)}{r^4}[/tex]

Which gives me the closest radius for a stable orbit of a photon as 3M
Apologies for the lack of basic physics understanding but for a photon on approach what would it's angular momentum be?
George Jones
#30
Nov23-09, 09:07 PM
Mentor
George Jones's Avatar
P: 6,239
We now have a system of three first-order coupled differential equations that determine the worldlines for "photons". In order to solve these equations, we need three initial conditions. Imagine firing a photon from a laser. One possibility for the three initial conditions: the [itex]\left( r , \phi \right)[/itex] position from which the photon is fired and the direction in which the photon is fired. These initial conditions determine [itex]L[/itex]. Alternatively, [itex]L[/itex] could be used as one of the initial conditions.

Much more on all of this later.
George Jones
#31
Nov24-09, 08:44 AM
Mentor
George Jones's Avatar
P: 6,239
Quote Quote by nocks View Post
That would be

[tex]\frac{d^2r}{d\lambda^2} = - \frac{L^2(3M-r)}{r^4}[/tex]

Which gives me the closest radius for a stable orbit of a photon as 3M
Reduce this second-order equation to two first-order equations by setting [itex]p = dr/d\lambda[/itex], so that [itex]dp/dt = d^2r/d\lambda^2[/itex]. The set of equations that describes the worldline of a photon then is

[tex]
\begin{equation*}
\begin{split}
\frac{d \phi}{d \lambda} &= \frac{L}{r^2} \\
\frac{dr}{d\lambda} &= p \\
\frac{dp}{d \lambda} &= \frac{L^2(r - 3M)}{r^4}.\\
\end{split}
\end{equation*}
[/tex]

Assuming that all required values are given, write a few lines of (pseudo)code that uses the simplest, most intuitive method (Euler's method) to solve these equations.
nocks
#32
Nov26-09, 06:37 PM
P: 24
Having a bit of trouble with this but I have p as

p = [tex]\frac{L^2(2M-r)}{2r^3}[/tex]
nocks
#33
Mar9-10, 09:14 AM
P: 24
wow it's been a while since I looked at this project.
Thought i'd come back to it :)

Quote Quote by George Jones View Post
Reduce this second-order equation to two first-order equations by setting [itex]p = dr/d\lambda[/itex], so that [itex]dp/dt = d^2r/d\lambda^2[/itex]. The set of equations that describes the worldline of a photon then is

[tex]
\begin{equation*}
\begin{split}
\frac{d \phi}{d \lambda} &= \frac{L}{r^2} \\
\frac{dr}{d\lambda} &= p \\
\frac{dp}{d \lambda} &= \frac{L^2(r - 3M)}{r^4}.\\
\end{split}
\end{equation*}
[/tex]

Assuming that all required values are given, write a few lines of (pseudo)code that uses the simplest, most intuitive method (Euler's method) to solve these equations.

A few questions. Since I lack a good maths background, I would appreciate some advice on how I would implement these 3 equations into a numerical solver like Euler's or Runge-Kutta to get back some results that I could plot.
e.g. Do I pre-define L?, how do I know which direction the photon is travelling?

currently trying to use this : http://www.ee.ucl.ac.uk/~mflanaga/java/RungeKutta.html


Register to reply

Related Discussions
Speed of light in a blackhole Astronomy & Astrophysics 3
Light and blackhole Astronomy & Astrophysics 5
Blackhole Question (involves Light) Astronomy & Astrophysics 7
Light in blackhole Astronomy & Astrophysics 15