I am currently trying to calculate deflection angles for the path of photons near Schwarzschild black holes, using the formula

$$\Delta\phi=2\int_{r_0}^\infty \frac{dr}{r^2\sqrt{1/b^2-1/r^2(1-1/r)}}$$

To do this, I solve the equation of motion of the photon numerically:


(in units where $2GM=1$, $c^2=1$ and hence $r_s=1$ (the Schwarzschild radius).

b = 2;
sol = NDSolveValue[{u''[p] + u[p] == 3/2*u[p]^2, u[0] == 0, 
u'[0] == -1/b}, u, {p, 0, 5}]

Then to find $r_0$, the radius of closest approach, I use NSolve[sol'[p] == 0, p], which returns {{p -> 1.26921}}. And finally, the integral is evaluated with

2*NIntegrate[ 1/( r^2*Sqrt[1/b^2 - 1/r^2 (1 - 1/r)]), {r,-1/sol[1.2692092946080926], Infinity}]

Which comes out to 1.8498 radians. This seems incorrect, because as shown in the screenshot below from these lecture notes (p.35, author mis-typed his lower integration bound on first integral), the angle should be at least larger than $\pi$. Where is the issue? enter image description here

My only suspicion is the initial condition $u'(0)=-1/b$ but it seems correct given the following argument:

$du/d\phi= (dr/dl)(du/dr)(dl/d\phi)$

Where $l$ is the affine parameter $dr/dl=1$, the speed of light (this is where I'm not too convinced as affine parameter is not proper time) then use E-L equations to get $dl/d\phi$ and $u=1/r$.

  • $\begingroup$ Equation is not needed in this problem. The angle is determined by the integral. $\endgroup$ Commented Dec 4, 2018 at 0:36
  • $\begingroup$ @AlexTrounev, completely agree that angle is determined by the integral, the equation (from lecture notes) I just included so that we can verify the answer Mathematica gives. My point is that the equation implies that the angle (given by the integral) must be greater than pi, but Mathematica gives a number smaller than pi. $\endgroup$
    – tbfr416
    Commented Dec 5, 2018 at 8:12
  • $\begingroup$ This should be exactly solvable in terms of Weierstrass elliptic functions; if you are interested in that, I can try to write something up. $\endgroup$ Commented Jan 5, 2019 at 8:00

1 Answer 1


The angle is calculated as an integral depending on the parameter b. Calculate the minimum distance

 r0 = r /. First[Solve[1/b^2 - 1/r^2 (1 - 1/r) == 0, r]]

(*Out[]= ((2/3)^(1/3) b^2)/(-9 b^2 + Sqrt[3] Sqrt[27 b^4 - 4 b^6])^(
 1/3) + (-9 b^2 + Sqrt[3] Sqrt[27 b^4 - 4 b^6])^(1/3)/(
 2^(1/3) 3^(2/3))*)

Calculate the angle for example when b = 4

  1/(r^2*Sqrt[1/b^2 - 1/r^2 (1 - 1/r)] /. b -> 4.), {r, 
   Re[r0 /. b -> 4.], Infinity}]

(* Out[]= 4.00032*)

In this problem there is a minimum value of b=b0=(3 Sqrt[3])/2. For b < b0 the minimum radius becomes complex, for example for b = 2

r0 /. b -> 2.
(*1.191487883953119` - 0.5088517788327379` I*)

This means that light is captured by a black hole.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.