In a particle simulation, it's pretty easy to compute gravitational interactions directly from Newton's Law of Universal Gravitation: f = g m1 m2 / r

g is the gravitational constant

m1 and m2 are the masses of the objects

r is the distance between the centers of the objects

f is the magnitude of the force

Since you want a vector force, you need to scale this magnitude f by a vector, which can be the normalized direction vector R between the object centers. Or you can save one vector scale factor by using the non-normalize vector between centers R, and calculating F = R * g m1 m2 / r

The other problem in a simulation is particles will occasionally pass very close to one another, resulting in a distance r near zero, which gives a force near infinity. To avoid this, the force is typically "softened" by dividing by 1/(a+r

One fairly simple way to get much better performance is to only
consider nearby neighbors. You can find your nearby
neighbors by storing lists of particles by their integer
coordinates in a hash table, and looking through your list to find
your closest neighbors. But there's a problem: the galaxy
can't hold together anymore, even if we artificially inflate the
gravitational constant, because the long-range forces that bind
one side of the galaxy to the other are never evaluated.

The hashtable version gives the wrong answer because it ignores
long-range interactions. There are several neat ways to
approximate these interactions.

A very simple approach is to sum all the particles into one
"megaparticle", which will end up at the center of the
galaxy. We can then compute each particle's interaction with
its local neighbors, and then with the megaparticle. This
doesn't even conserve momentum, because the influence of a
particle on the megaparticle isn't symmetric, but at least the
galaxy holds together!

A more accurate version of this is called Barnes-Hut:
for neighboring particles, we compute full particle-particle
interactions. To interact with a group of distant particles,
we lump the group into a megaparticle. To accurately define
which particles form a "group", we impose a spatial tree across
the simulation, using the coarse levels of the tree to approximate
distant regions, and increasingly fine levels for more nearby
regions. See my Barnes-Hut
Gravity Tree code.

Here's some measured performance of gravity via Barnes-Hut:

n Particles | Tree Build | Lump Tree | Gravity | Export |

20 | 0.3 | 0.015 | 0.31 | 3.2 |

200 | 1.2 | 0.03 | 5.4 | 23.5 |

2000 | 23.6 | 0.3 | 307.7 | 227.3 |

Scaling |
O(n) + big constant? |
O(n) |
O(n log n) |
O(n) with huge constant |