How to Code a Blackhole in Python

Let's use python and matplotlib to gain an intuition for coding difficult physics problems! Thank you to my advisor Dominic Chang at Harvard BHI for making and teaching me how to do this.
Source Material »

Table of Contents

1. About The Project
   • Built With
2. Getting Started
   • Set-up
3. Pseudocode
4. Breakdown of Code
5. Results
6. Contact

About The Project

This mini lecture is made to guide you through the process of coding a black hole in python but also to gain a more fundamental way of thinking when trying to code a physics problem.

My goal is that you are able to comment the code provided on your own, knowing what the functions do conceptually, along with feeling more confident with coding physical systems.

The project simulates gravitational lensing around a Schwarzschild black hole by:

• Calculating light ray trajectories using geodesic equations
• Visualizing how gravity bends light paths
• Creating an accretion disk with multiple light-path intersections

(back to top)

Built With

• 
• 

(back to top)

Getting Started

Go into 'src', there you'll have the choice of a jupyter notebook or regular python file. I suggest clicking the badge above for Jupyter notebook and importing the .ipynb file into that, since WPI gives us JupyterHub! I've included the notes below in the .ipynb file, so if you want everything in one place I highly suggest that.

Set-up

Either use the Jupyter Hub link for WPI, or make sure you have python downloaded. You can easily make sure you have python, or any other program downloaded by typing:

python -V

This should return the latest version of python.

Once you have python, run command prompt as administrator and install the required packages:

• scipy

  pip install scipy

• matplotlib

  pip install matplotlib

• numpy

  pip install numpy

Pseudocode

So we want to trace rays from uniform rings in the bulk space of a black hole to our observer.

Here's the step-by-step approach:

1. Set up the geodesic equations

   • Define the radial potential function for light rays in Schwarzschild geometry
   • This determines how light moves in the curved spacetime around the black hole

2. Solve for light trajectories

   • Use numerical integration to calculate φ (azimuthal angle) as a function of r (radius)
   • Handle the turning points where light rays reach minimum radius
   • Implement alternative parameterization to avoid numerical issues

3. Generate multiple trajectories

   • Loop through different impact parameters (b values)
   • Each b value represents light coming from different distances from the black hole
   • Plot trajectories to visualize gravitational lensing

4. Create the accretion disk

   • Define observer angle θ and viewing geometry
   • Calculate where light rays from the disk intersect with different radii
   • Handle multiple intersections (n=0, n=1, n=2) representing light that wraps around the black hole
   • Map to screen coordinates (α, β) for visualization

(back to top)

Breakdown of Code

Core Physics Functions

radial_potential(r, b)

• Calculates the effective potential for light rays
• Formula: r⁴ - b²r² + 2b²r
• This determines whether a light ray at radius r with impact parameter b will continue outward or fall into the black hole

minrs(b)

• Finds the minimum radius a light ray can reach for a given impact parameter
• Uses the cubic formula to solve for turning points
• Returns 0 if the ray falls into the black hole (complex result)

d_phi_dr(r, b) and phi(rs, b)

• Calculates how the azimuthal angle φ changes with radius
• Uses numerical integration (scipy's quad function)
• Integrates from the screen radius to infinity

Trajectory Generation

generate_trajectory(ximin, ximax, b, steps)

• Creates a complete light ray path for visualization
• Uses alternative parameterization to handle turning points
• Returns x,y coordinates for plotting

Accretion Disk Visualization

The final section creates the iconic black hole image you see in movies:

1. n=0 (blue): Direct light from the disk
2. n=1 (green): Light that wraps halfway around the black hole
3. n=2 (red): Light that makes a full orbit before reaching us

Each color represents light from the same accretion disk but following different paths around the black hole!

Key Parameters

• theta: Observer viewing angle (85° gives a nice tilted view)
• rmin, rmax: Inner and outer radii of the accretion disk
• steps: Resolution of the calculation (higher = smoother but slower)
• M=1: Black hole mass in geometric units (sets the scale)

The Math Behind It

The code solves the geodesic equation in Schwarzschild geometry:

dφ/dr = ±b/√(r⁴ - b²(r² - 2Mr))

Where:

• φ is the azimuthal angle
• r is the radial coordinate
• b is the impact parameter (related to angular momentum)
• M is the black hole mass (set to 1)

The shadow you see has a radius of √27M ≈ 5.2M, which is the photon sphere radius where light can orbit the black hole.

(back to top)

Results

Running the complete code produces:

• Individual trajectory plots showing gravitational lensing
• Multiple trajectories at different impact parameters
• Complete accretion disk visualization with color-coded light paths
• Black photon sphere at the center (radius √27 ≈ 5.2)

The resulting image shows how a black hole warps spacetime, bending light around it in beautiful and complex ways!

(back to top)

Contact

Adam Field - adfield@wpi.edu

Project Link: https://github.com/AdamField118/Blackhole-in-Python

(back to top)

---