Black Holes and Fractal Basins

During one of my many spells stumbling around on the internet, I came across a very cool paper on modelling the trajectories of particles around black holes (bear with me…). It has all the required ingredients for this blog: maths, physics, nice pictures, and the chance to throw in a Muse reference. Let’s have a look.

Black holes, in case you aren’t aware, are the remnants of dead stars, stars so massive they have collapsed in upon themselves to an infinitely dense point. Their gravitational fields are too strong for Newton to keep up with, so we need to turn to General Relativity in order to examine their properties.

Einstein developed the Einstein Field Equations (the clue is in the name) which describe how the presence of mass distorts spacetime. More precisely the stress-energy tensor of an object, a mathematical structure encoding the mass and energy contained, defines the local structure of spacetime by modifying its metric. The metric in turn is another bit of maths which just tells you how to work out the length of lines. In the flat spacetime Newton was used to the metric is very simple, and in fact everybody knows it from the formula

\Delta s^2 = \Delta x^2 + \Delta y^2 + \Delta z^2

where \Delta s is the length of a line, and \Delta x,y,z are the coordinate displacements – this is just Pythagorus’ theorem. In the limit that \Delta x,y,z \rightarrow dx,y,z \rightarrow 0  this becomes a differential relation.  The metric defines a geodesic, which is the shortest path between any two points. In flat space this is just a straight line, on the Earth’s curved surface it is a ‘great circle’, which is the path that aeroplanes take to minimise fuel usage.

Around a black hole things become more complicated, but by calculating the geodesics it is possible to see which paths a particle will take. There is no real gravitational ‘force’ in general relativity then, you are simply free-falling all the time along the geodesics that matter and energy are imprinting on your local spacetime. The force you feel on your feet as you walk around is the resistance to your natural free-fall path to the centre of the Earth.

So, with that out of the way, the authors of the paper I mentioned are concerned with the paths particles take around a cluster of black holes. It turns out that there is a surprisingly simple metric for this, the Majumdar-Papapetrou metric

ds^2 = -U^{-2}dt^2 + U^2\left(dx^2 + dy^2 + dz^2\right)


U(x,y,z) = 1 + \sum_{i=1}^N\frac{M_i}{\sqrt{(x - x_i)^2 + (y - y_i)^2 + (z - z_i)^2}}

and we have N black holes at positions (x_i,y_i) with masses M_i. This metric is static, which means the black holes don’t move around. This is a good approximation if the particle moves over a timescale shorter than the black hole orbital timescale.

Letting the authors do the hard work for me, the equations of motion for a neutral particle in 2D with position \mathbf{x} and momenta \mathbf{u} are

\dot{\mathbf{u}} = U^{-2}\left((1+ 2u^2) \nabla U - \mathbf{u}\mathbf{u}\cdot \nabla U\right)

\dot{\mathbf{x}} = U^{-1}\mathbf{u}.

As everything in relativity is, well, relative, these momenta are the ones which would be observed by someone at rest with respect to the metric. The derivatives are taken with respect to the proper time of the particle. The actual velocity is \mathbf{v} = \mathbf{u}/\sqrt{1 + u^2}. Knowing the equations of motion, I can set up some black holes, stick in a particle and see what happens using whichever ODE integrator methods tickle our fancy. Let’s start the particle at rest between the two black holes and see where it goes. The starting position I’ll call (x_o, y_0). Red stars mark the black hole positions, and the triangle is the starting position of the particle.

Straight through

Naturally, due to the symmetry involved, it flies straight through and out the other side. If we move it up a bit:


Again, unsurprisingly, it is captured by the black hole nearest to it. Keeping it still in the top half of the plane, there are also trajectories like this one where the particle swoops round into the other black hole.


In fact, it can be difficult to predict which black hole the particle will end up in, if its initial position changes even by a small amount, it can switch between the two – see the animation below . This is exactly the hallmark of a chaotic system as the authors of the paper point out.


In fact we came across this type of behaviour before on this blog, when I was talking about Newton fractals. There we picked some initial conditions and iterated a complex mapping, and our point moved around until it settled at a particular place. An equivalent thing is happening here, but instead of a complex map there is a system of differential equations. The concept of a basin of attraction is therefore useful again here – regions of phase space which send a particle to one black hole or another. In the Newton fractal it turns out that the basins of attraction have fractal boundaries for most polynomials. The same is true here, I’ve plotted below a small region of phase space where black means the particle ended up at one black hole, and white the other.


This shape has some strange properties, however far you ‘zoom in’ you can always find some holes or gaps. For this reason it isn’t quite two dimensional, but then again neither is it one dimensional as it can have some thickness. Rather, it has fractal dimension 1.43 as calculated in the paper, something in between a line and a plane. In fact, if you haven’t already you really should have a flick through, the authors manage to find some beautiful structures. In a very real sense though this object is one defined purely mathematically and deterministically, you really have to go in and explore it to see what surprises are in store. I encourage all of you with the means to simulate this to have a go – my code is attached at the bottom of the post.

I’ve also plotted below a selection of actual paths from this calculation. At this scale it looks like they all come from the same point, but actually they’re displaced very slightly. A very small displacement then leads to a wildly different path, chaos again! In fact the authors predict just how sensitive this system is to displacements in initial conditions, in terms of something called the Lyapunov exponent. This measures how quickly two similar trajectories diverge, and you may be interested to know for this system it is 18\cdot 5^{-3/2}, something the authors verified numerically.


That takes a long time to compute though, for every point one must integrate a set of differential equations for many time steps. Instead let’s have a look at some nice trajectories:


And some with lots of black holes:


If you’d like to try this for yourself in Matlab, try the code here, unfortunately attached as a Word document as that’s all I’m allowed. You’ll need to supply

  • xbh – vector of black hole x positions. I used (0,0)
  • ybh – vector of black hole y positions. I used (-1 1)
  • M – vector of black hole masses. I used (1/3,1/3).
  • x0 – initial x position of particle
  • y0 – iniial y position of particle
  • ux0 – initial x velocity of particle
  • uy0 – initial y velocity of particle
  • plotflag – whether to plot or not, 1 or 0

Let me know if you get any nice trajectories!


11 thoughts on “Black Holes and Fractal Basins

  1. Hi, do you have any further idea about the nature of the trajectories that you display right after the video? Do you know if for any k, there is always a trajectory that loops k times around one of the holes, or both? Thanks!


    1. Good question, I was thinking along similar lines. I found these trajectories by the following process:
      1. Start somewhere with a trajectory which collides with a black hole
      2. Increase a component of initial velocity until particle escapes
      3. Drop that velocity below escape threshold and switch to other velocity
      4. Repeat
      I found that this produced trajectories which were gradually more ‘loopy’. I was thinking about running some search algorithm but needed a good metric for ‘loopiness’ which I could measure. Any thoughts?


  2. The high-level answer would be that you want some kind of equivalent of the winding number (or index of a curve) used in Cauchy’s integral formula.

    I think that the notion of fundamental group of the plane with two points removed could help you with this definition. The problem I would see is that these are defined for closed curves, not for the trajectories that you have here. Looks very interesting though.


  3. Hi Jason,

    My PhD student and I came across your interesting blog article some months ago. At that time, we had just started thinking about the shadow cast by a binary black hole (essentially, the basin of attraction). Your article, among others, encouraged us to start investigating the Majumdar–Papapetrou case. This has led to a preprint, “Binary black holes shadows, chaotic scattering and the Cantor set”, available here:

    Please let me know if there is some appropriate way to acknowledge your blog post in our article.



    1. Hi Sam,

      Thanks for getting in touch, and for the link to the interesting paper. I’m perfectly happy to be acknowledged, though I only re-implemented the ideas in the Dettmann preprint. As such a citation probably isn’t appropriate, perhaps in an acknowledgements section? Good luck with publication either way!

      Cheers, Jason


Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s