Is there an algorithm for N body simulations in General Relativity

I am new to general relativity but have a background in computer science. Why is it so hard to do n-body simulations in GR? For example, there could be a program which takes the properties (mass, position, velocity, etc.) of each particle as input and numerically integrates the evolution of the system using discrete but very small time steps. If we throw enough memory and computing power at it, we can get arbitrarily close to the 'real' solution. I assume the entire histories of the worldlines would need to be known, in order to calculate things like gravitational waves. I have downloaded programs which show the evolution of the Schrodinger equation, so I am surprised they are so hard to find for GR. Some responses I've heard elsewhere are 'the equations are nonlinear,' but from what I hear numerical algorithms are great for these types of problems.

What types of numerical algorithms are out there for GR?

Answer

N-body simulations in full general relativity are difficult because gravity is a field theory and because it is non-linear.

Let's deal with the field theory part first. In Newtonian mechanics gravity is static. The field itself has no energy or momentum, no degrees of freedom at all. It is simply an instantaneous force law between all matter. Remove the matter, no gravity. This is the context where most N-body simulations are done: your system is defined completely by the masses, positions, and velocities of your particles.

Now, consider electromagnetism. Electromagnetism is a full dynamical field theory. The electromagnetic field carries energy and momentum independent of whatever charges happen to be present, like electromagnetic waves. For charged particles interacting electromagnetically, you cannot describe your system by their instantaneous masses/charges, positions, and velocities alone. You must know the electromagnetic field as well. For instance, the evolution of the same initial set of particles will be very different if the initial field contains strong electromagnetic waves. There are codes/algorithms that simulate this system, typically called Particle-In-Cell (PIC) codes. They use Newton's laws to evolve the particles (as in an N-body) and also evolve the electromagnetic field dynamically by solving Maxwell's equations.

General relativity (GR) is the field theory extension of gravity, and is exactly like electromagnetism except the equations are much more complicated. Newton's Laws for particle motion are replaced by the geodesic equation, and Maxwell's equations are replaced by the Einstein Field Equations. I'm not sure if this has been done, but I think you could in principle write a GR PIC simulation as with electromagnetism. This strictly speaking would not be an N-body simulation because you'd be simultaneously evolving the gravitational field, but under a certain approximation it could probably be done. The difficulty is that the Einstein Field Equations are MUCH harder to solve than Maxwell because of the non-linearity.

However, this normally is not done because it doesn't gain you much. For systems that you'd want to approach with an N-body simulation (like stars in a galaxy, or galaxies in a cluster, for instance), the static Newtonian approximation to gravity is extremely good and you would gain nothing but a headache trying to approach it relativistically. In electromagnetism, the dynamics of the field (its own energy and momentum) are important when it has waves that can interact and affect the matter. In GR this only happens in the "strong field" limit, typically very close to a black hole.

Simulations of one or several black holes are an entirely different ball game than N-body simulations. In these situations any matter in the system is probably an astrophysical plasma that you would model with the hydro equations, not as a bunch of particles. The black holes form massive singularities on your grid that you excise as they move, orbit, and merge with each other. The first successful calculation of a black hole binary, 2 black holes orbiting and merging, was only done in 2005. We've gotten better at it, but it is still a very hard problem.

tl;dr The equations governing the true evolution of a multiple black hole system are incredibly different from a simple N-body system with a force law.