Capture and Escape of a Photon at a Black Hole: as a trajectory r(t) and as a path r(φ)

By Alexander Gofen (2025)

Abstract

We study the Ordinary Differential Equations (ODEs) representing the path and trajectory of light approaching a (non-spinning) black hole. In the vicinity of a black hole, a light ray dramatically bends. To express the path of light via ODEs near a black hole is possible due to the Schwarzschild solution of the Einstein equations of the General Theory of Relativity.

We explore those ODEs and graph their solutions using the Taylor Center ODE solver. Some remarkable features of those ODEs present a particular interest for teaching – and a challenge for numeric integration.   

Preface

Since the beginning of science and mathematics, the light rays were viewed as the best physical realization of the mathematical abstraction of a straight line. Though a light beam can reflect in a mirror, or bend at the boundary of two mediums, in the vacuum, the light remained the ultimate standard of the very idea of straightness.

However, in 1916 Albert Einstein published his General Theory of Relativity. One of its consequences was that masses bend the curvature of the space, and the light bends too in the proximity of (big) masses. While sounding so bizarre in that time, the new theory predicted a small, but measurable effect of light bending near the Sun.

To measure this effect, astronomer needed a full solar eclipse. Such an eclipse happened in 1919, and the measurement of shifts of stars near Sun disk did demonstrate a small shift merely about one angular second predicted by the theory!   

Since then, the light bending has been observed many times while watching remote galaxes, and it was called gravitational lensing. Still, this effect measures only in parts of angular seconds.

In contrast, here we consider effects of light bending in a close proximity to black holes where light ray of a photon can make U-turn or even wind around the black hole before escaping it, or getting caught otherwise. It seems fascinating to study such paths of light – even if we will never be able to observe them directly.

Here we are to present two models of a photon motion approaching a non-spinning black hole based on the Schwarzschild solution of the Einstein equations of the General Theory of Relativity. Unlike the original equations of Einstein, the Schwarzschild model relies on Ordinary Differential Equations (ODEs) available for both static and kinematic models of a photon motion. In the static model the photon motion is given only geometrically as a path r(φ) in polar coordinates, while in the kinematic model it's a trajectory r(t), where time t corresponds to an infinitely remote inertial frame (the lab).  

The ODEs of both models are instructive in themselves for the following reasons.

· Their right-hand sides are square roots, and both branches are required to represent the solution. However…

· It's possible to rid of the square roots transforming the ODEs into the regular second order ODEs.

· At the critical value of the parameter b, numerical integration of the ODEs is a challenge because a catastrophic subtraction effect embedded in this problem.

Our setting is this:

Fig. 1

The origin of the polar coordinates is in the center of the black hole. The black circumference is the horizon whose radius is 2GM, while the blue circumference is the so-called circular orbit whose radius is 3GM (here G is the gravitational constant, and M – the mass of the black hole).

The model of photon motion for r(t) and that for  r(φ) both use a special impact parameter b conserving its value and defined in the following picture:

(1)

Fig. 2

For convenience, here we chose the point (r₁, 3π/2), where r₁  and α are arbitrary, however, formula (1) defines the parameter b for any point (r, φ) of the path. If α = 0 or α = π, the photon moves exactly along the line between itself and the center, and b = 0.

The fact that the value of b is conserved along the path once it was chosen means that the angle of direction α is actually a function α(r)  implicitly defined by formula (1).

The parameter b corresponds to the tangential component of the direction of motion of a photon given by the angle α. When α = π/2, b expresses purely tangential motion for any point (r, φ). The direction of the radial motion (inbound or outbound) is defined by the proper choice of sign of the square roots in the formula (2) for r(t), or (5) for u(φ) = 1/r(φ) below.

The meaning of b is the same in both models. In particular, we distinguish the critical value b_crit = 3GM√3 delivered by formula (1) for r₁ = 3GM and α = π/2 at the point (3GM, 3π/2) on the circular orbit. The solution corresponding to this point is the circular orbit. When b = b_crit, the circular orbit separates the inbound and outbound orbits corresponding to different initial values r and α, as we will see further.

Indeed, both models for r(t) or for r(φ) generate the same paths for the same initial values. However, in order to utilize the dynamic graphing of the Taylor Center, we are particularly interested to begin with visualizing the ODEs (2) modelling the time dilation of the trajectories r(t) when approaching the horizon.

In order to run all the graphs here as simulations under Windows, you have to download and install the application TCenter, and download and unzip the script files of all the simulations.

The simulations for the section below are in the subdirectory \r(t).   

The ODE for u(φ) is simpler and more elegant. We will explore it afterwards.

The motion as a trajectory r(t)

Unlike the static model in the next section, this kinematic model for the trajectory r(t) makes sense only where the time t exists, i.e. only in the space outside of the horizon. Consequently, the solution r(t) is not continuable inside the radius r = 2GM, while the solution of ODEs  (5, 6), being a path, extends also inside the horizon.

The original ODEs of a photon trajectory r(t) are

(2)

Remark 1: Two ODEs for the same r(t). The ODE (2) for r contains the square root in the right-hand side with double signs representing the two branches of the complex square root. Despite that, the solution r(t) is a single holomorphic solution of these ODEs.

How do we know that r(t) is a single holomorphic solution of these singular ODEs? – Because r(t) is also a solution of holomorphic ODEs (3), just as u(φ) is a solution of holomorphic ODE (6) so that both  r(t) and  u(φ) are holomorphic even at the point where the square roots turn zero.

     The two branches in the right-hand side represent the decreasing and increasing segments of the same solution r(t). During attempt of escape from inside out, r(t) increases satisfying ODE (2) with "+"sign until possibly reaching the maximum, after which r(t) decreases satisfying ODE (2) with "–"sign. During approach from infinity, r(t) decreases satisfying ODE (2) with "–"sign, until possibly reaching the minimum, after which r(t) increases satisfying ODE (2) with "+"sign.

Remark 2: "Removable singularity" of the ODE. ODE (2) has a branching singularity when the under-the-root expression becomes zero – and this may happen in reality. However, as we will demonstrate later, the ODEs for r(t) or r(φ) belong to those remarkable ODEs, which, despite being singular, have regular solutions. An example of such an ODE is  x' = nx/t  which is singular at t = 0. However, for natural n, it has a holomorphic solution x = tⁿ. Besides, there exists also a regular ODE x' = nt^{n – 1} having the same solution. We will see, that for the holomorphic solutions r(t) of (2) and u(φ) of (5), regular ODEs without a square root also exists.

Remark 3: Specialty of the regularity of this solution. We know, that some functions expressed as a superposition of two functions, are holomorphic despite that one of the two is singular. For example, f (z) = cos(√z) is holomorphic at z = 0 despite that √z is singular at zero. The external holomorphic function cosine "kills" the singularity of the internal one. The situation with the holomorphic functions r(t) or r(φ) is the opposite. Their derivatives r'(t) and r'(φ) (or u'(φ)) are holomorphic too. When an approaching photon escapes, both radii reach their minimum – the minimal distance to the black hole. Indeed, at the point of minimum, r'(t) = 0 and u'(φ) = 0 both expressed as square roots in formulas (2) and (4) below: r'(t) = √g(t)  and  u'(φ) = √f(φ). Here the under-the-root functions g(t)  and  f(φ) are holomorphic. However, despite that the square root is singular at zero, the superposition is holomorphic, i.e. now the internal holomorphic function somehow kills the singularity of the external one! Wow! What kills the singularity in this case? – The special properties of the under-the-root functions g(t)  and  f(φ). Not only don't they ever get negative under the root, but even when they do reach zero (at points of the nearest approach to the black hole), their minimums are so "flat" at the zero, that their Taylor expansion somehow "regularizes" the singularity of the square root so that analytic continuation of r'(t) and u'(φ) smoothly transitions from one branch of the square root to the other! For example, in the superposition f(z) = ±√(z³ – z²) the internal function is polynomial, therefore holomorphic. It turns into zero at z = 0. However, both branches of f(z) are also holomorphic, which gets obvious after this simple transform f(z) = ±z√(z –1). So too, the properties of the solutions r'(t) and r'(φ) in the under-the-root functions g(t)  and  f(φ) kill the singularity of the square root.    

Remark 4: When to switch the sign. The solution r(t) may have an extremal point where r'(t) = 0 because the under-the-root expression turns into zero. If this is the case, r'(t) in (2) must change the sign going over this point. This translates into the necessity to switch the sign at the square root, both branches of which represent the same holomorphic solution r(t). An attempt to integrate ODE (2) with the Taylo method over this extremal point would fail, but we will integrate a regular ODE (3) instead.

Remark 5: A quadrature. Here the ODE (2) in r is uncoupled, and the variables t and φ are absent in the right-hand side of it. As such, it has a quadrature – the integral expressing t(r), though not in the Liouville elementary functions.    

Theorem 1: The ODE (2) has two types of special stationary solutions:

1. For any point (r = 2GM, φ) on the horizon with φ and b arbitrary, the exceptional motionless solution r, φ = const and r', φ' = 0 satisfies the ODE so that the entire orbit reduces to just one point of no motion.

2. For any point (r = 3GM, φ) with b = b_crit and φ arbitrary, r = const = 3GM, r' = 0, satisfies the ODE presenting the circular orbit.

Proof: It's easy to see that this system has a stationary motionless solution r ≡ 2GM= const and φ = const  – a point on the horizon. At this point, for an observer in a remote inertial lab, the photon will be perceived as motionless, and the time stopped. Any solution touching the horizon would hit into a point of no motion. If the numeric integration were infinitely accurate, we would see just that r(t) → 2GM, and (1 – 2GM/r) → 0 always remaining positive. However, due to inevitable rounding errors and finite number of digits, at certain step the parenthesis gets negative so that the further integration evolves into a numerical artefact with radius r increasing (making no physical sense). This is the first type of special solution. 

           Now consider the under-the-root expression in the formula (1). For b = b_crit = 3GM√3 and r = 3GM the under-the-root expression turns into zero, proving the case (2). ■

See here for more theorems concerning these ODEs and their solutions.  

In order to rid of the square root, let's transform the ODE (2) into a second order ODE by differentiating both sides of (2). Such a differentiation is particularly simple when applied to an ODE of the format r' = √f(r). Then

r" =	f'(r)r'	=	f'(r)√f(r)	=	f'(r)
	2√f(r)		2√f(r)		2

In our case, f(r) is the result of placing the external parenthesis of formula (2) under the square root. With that in mind, we will get a second order ODE 

(3)

where

The required initial value for dr/dt is obtainable from the respective ODE (2) at any regular initial point with the appropriate choice of the sign to specify the inbound or outbound direction, and this is the only usage of ODE (2).

In terms of r, φ, the speed c of photon is

containing time derivatives of r and φ so that the speed of a photon (in the lab frame) varies in the vicinity of the black hole (rather than being the known constant c of the speed of light in the vacuum afar from any masses).  

Below are results of integration of these ODEs for several special settings.

Inbound motion: from infinity toward the black hole.

Here are a few pictures of photon trajectory r(t) from infinity toward the black hole meaning that in ODE (2) for initial velocity we chose the sign "+".  (In order to watch the animated motion, you need to download and run the TCenter application).  

Fig. 3. A typical trajectory of a capture. The parameter b < b_crit forcing approach of a photon (in red) towards the horizon (in black) in a finite spiral fall.

Fig. 4. The graph of the speed c(r) of a photon on Fig. 3, which moves towards the horzion along a finite spiral.

The vertical red line denotes the radius of the circular orbit. The blue line denotes the horizon.

Fig. 5. A typical trajectory of an escape enforced by the parameter b > b_crit

Fig. 5. An atypical trajectory of an escape enforced by the parameter b > b_crit yet almost equal to b_crit.

The photon (in red) makes several lapses around the circular orbit.

The lapses are so tight that they are undiscernible on the graph.

Fig. 6. A capture by the circular orbit rather than by the horizon.

The parameter b = b_crit so that the photon (in red) approaches and is caught by the circular orbit along an infinite spiral

whose lapse are so tight that they are undiscernible on the graph. The horizon is in black.

The integration was conducted with the highest accuracy in order that the trajectory do not go astray.

Fig. 7. The graph of the speed c(r) of a photon on Fig. 6. The photon starts at the big radius (at the right)

approaching and being caught by the circular orbit in infinite spiral so that radius diminishes.

The vertical red line denotes the radius of the circular orbit. The blue line denotes the horizon.

Fig. 8. The graph of the speed c(r) of a photon which moves towards the center along a straight line (b = 0).

The blue line denotes the horizon. The vertical red line denotes the radius of the circular orbit.

Here the c(r) decreases to zero steeper than in other cases above.

A puzzle of the c(t)

Unlike in the previous subsection, here we consider a dependency of the speed of a photon c(t) on time rather than on the distance to the center.

All three graphs of the varying speed c(r) of photon above are displayed as functions of radius-vector r, i.e. the distance from the center of the black hole. In all those cases we see that the speed of the photon monotonously decreases as the r decreases getting closer to the horizon. At that, all the curves c(r) remain clearly convex.

It's not so for graphs c(t) having an inflection point: at least one.

Fig. 9. The graph of the speed c(t) of a photon for the circular approach Fig. 6.

Fig. 10. The graph of the speed c(t) of a photon for the strait line approach Fig. 8.

These two graphs above have one inflection point.

Now take the look at this combined graph of photon speed c(t) (green), corresponding to the case Fig. 3 when a photon moves along a finite spiral towards the horizon.  

Fig. 11. Here φ(t) is blue, r(t) is red, and photon speed c(t) (displayed magnified as 25c(t)) is green.

All three graphs are still monotonous, however, unlike before,

now c(t) switch convexity/concavity in the process twice.

Why is the case of finite spiraling different?  

The expression for the speed of photon is

We can obtain either ∂²/ ∂r², or d²/dt² of the expression for c (the latter expression being more cumbersome than the former), and ask where these second-order derivatives turn zero. We can also observe that c²(t) depends also on the parameter b which may explain the difference in inflection points taking place on Fig. 8.

Outbound motion from near the horizon: capture or escape.

Here we will take a look at a few examples how photon can escape away from a very close vicinity of a black hole slightly outside the horizon while inside the circular orbit.

If a photon was fired away slightly outside of the horizon, its chances to escape depend on the direction of firing controlled by the parameter b. It surely escapes if fired away exactly along the radius so that b = 0. Otherwise, the outcome depends on the angle α and b.  

Escape from a point near the horizon: b < b_crit.

Capture from a point near the horizon: b > b_crit

Capture of a photon fired from a point near the horizon with b = b_crit.

The trajectory is an infinite tight spiral approaching the circular orbit infinitely close from inside.

Capture of two photons: one (in blue) fired outbound from a point near the horizon,

the other (red) fired from outside inbound: both with b = b_crit.

The two trajectories are infinite tight spirals approaching the circular orbit infinitely close from inside and outside.

This is a zoomed-in image demonstrating tightness of the infinite spiral

winding around the circular orbit.

The approach as a path r(φ)

The simulations for this section are in the subdirectory \r(fi).

In this model the ODE is simpler. The ODE, as originally written for u(φ) = 1/r(φ) is:

(5)

Remark 6: Here the variable φ is absent in the right-hand side so that this ODE has a quadrature as an integral for φ(u), though not in the Liouville elementary functions.  

Remark 7: Also here, at the point when the distance r_min of a photon to the origin gets minimal (and  u_max= 1/r_min gets maximal), it must be that du/dφ = 0, possible only if the under-the-root expression becomes zero – the singularity of the root. An attempt to integrate this ODE until point r_min would inevitably fail. We must try to get rid of the square root also here. 

As before, we transform (5) into an ODE of order two eliminating the radical:

u" = 3GMu^2 – u

(6)

Not only is this ODE simpler than those for the temporal model, but it also uncovers a deeper reality: literally.

 This picture replicates the case of a spiral fall Fig. 3 considered above, but unlike there…

Now the red path of a photon penetrates inside the horizon until reaching the center of the black hole.

Wow! That's the new reality!

This is possible because, unlike the temporal model above, where the lab time stops at the horizon, this ODE describes merely the path of the photon motion which surely moves on also after crossing the horizon.

Here a photon makes a U-turn

And here the path replicates letter γ upside-down.

The accuracy challenge

When b = b_crit, the right-hand side  3GMu² – u  of the ODE (6) approaches zero when u approaches 1/3GM – the circular orbit. Both terms in the right-hand side are finite values, but their difference gets closer and closer to zero, leading to the situation of the so-called catastrophic subtraction error in the derivative u" = 3GMu² – u. With every next lap of the spiral, the accuracy of u" drops more and more, until this error causes a numeric deviation from the circular limit, so that both the trajectory r(t) or the path r(φ) go astray.

Here is the case with b = b_crit = 3GM√3, and the initial r > 3GM outside the circular orbit. 

We know that this solution must wind around the circular orbit closer and closer never escaping.

However, due to integration errors in the mode of default error tolerance,

the path deviates into an erroneous orbit of escape.

This is the same case with b = b_crit = 3GM√3, and the initial r > 3GM outside the circular orbit.

However here it was integrated with the smallest relative error tolerance

meaning that all 63 binary digits in the mantissa are enforced to be correct.

Now this correctly integrated solution very tightly winds around the circular orbit,

closer and closer, never escaping, as it should.

Approach to and capture by the circular orbit