Currently, I compute the force between two gravitational interacting particles in a simulation with $n$ bodies according to
$$F = G\frac{m_1m_2}{r^2}.$$
Doing this, however, assumes that all bodies in the simulation interact instantly with each other meaning that I assume the speed of light to be infinity.
Let's assume I want so simulate the formation of galaxy clusters where the propagation velocity of light plays an essential role. How would one incorporate the fact that the propagation speed of light is finite into an $n$-body simulation? Any ideas?