Gravitational waves

SpaceSim uses the classical, non-relativistic theory for the gravity simulation, governed by the Newton's law of universal gravitation. The only exception is the relativistic orbital decay, implemented as a small correction on top of the Newtonian gravity simulation.

The reason why the orbit of massive objects (neutron stars, black holes, etc.) shrink over time, without violation of the law of energy conservation, is the emission of gravitational waves. These waves carry energy away from the orbiting bodies in form of space-time perturbations.

To visualize the propagation of gravitational waves, I added a simple two-dimensional simulation that shows the solution of a wave equation, using massive objects as sources.

Why waves?

The idea of gravitational waves comes directly from Einstein's field equations. These equations describe how the curvature of space-time relates to the sources of gravity, mathematically described by the stress-energy tensor. If you want to learn more about the theory of gravitational waves and the history of their research, have a look at this paper.

While the Einstein field equations are extremely complex and can be solved analytically only in very specific cases (e.g. Schwarzschild solution), we can simplify them substantially under these two assumptions:

1. Look for a vacuum solution, assuming an empty space without any matter.

2. Assume we are far away from the sources of gravity and their gravitational influence is very weak. The solution should therefore be very close to a flat space-time.

This reduces the field equations to:

□ℎᵤᵥ = 0

which is known as the wave equation. The solution is therefore a wave propagating through a vacuum at the speed of light.

The wave equation can be solved numerically, however, it requires a very different method compared to the rest of the simulation. While the hydrodynamics and Newtonian gravity can be handled by the particle-based solver (where we only know the physical quantities at the position of particles), the wave equation requires a grid-based method, such as finite differences, as we need to keep the state of the space-time perturbations everywhere in space, not only where the particles are.

This comes with some limitations:

• The simulation must have a finite computation domain. Although in reality the gravitational waves can propagate to unlimited distances, it's not possible to simulate infinite space for obvious reasons.

• The waves are simulated in two dimensions only. Three-dimensional simulation is not computationally feasible. If the simulation runs on a regular grid with 1024x1024 cells, extending it into the third dimension would require 1024³ ≈ 1 billion nodes. Memory requirements alone make it impractical.

Note that the gravity affecting particles is still Newtonian. The gravitational waves are only used for visual output and have no feedback to the rest of the simulation, making it possible to solve the wave equation efficiently on the GPU.

Boundary conditions

As the simulation domain is finite in size, we need to provide boundary conditions to make the solution unique. In other words, we have to define how the waves interact with the boundary. The usual choice is to either keep the value fixed (Dirichlet boundary condition) or keep the gradient fixed (Neumann boundary condition), however, in this case, neither gives us the desired result. In both cases, the wave is still reflected with the same amplitude, only the phase of the wave is changed.

Instead, the solver exponentially attenuates the wave when it gets close to the boundary. The domain has a boundary layer (20 cells wide) that reduces both the wave value and its velocity, with attenuation getting stronger the closer it is to the boundary. This causes the wave to slowly dissipate instead of reflecting.

That's it for now. Stay tuned for more development updates.