Galaxies would require coupling of the particles, interactions through gravity... this would kill the efficiency of the solver. N-Particle simulations are always O(N^2) (or more).
It's not cutting off. It's replacing many long-range effects with a single one weighted more (and done it carefully so that the floating-point error is literally smaller than you could tell).
There's even better methods that take O(n) time, like the fast multipole expansion ones.
Btw this goes into numerical analysis, and O(-) is important but so is error convergence. I imagine these big-O figures are for a constant error, but I figure the asymptotic may change behavior as the other parameter changes.
One could also use FFTs and the particle-mesh method or a TreePM method. Basically, the particle masses are deposited on a fine mesh. Poisson's equation gets solved using FFTs and then the particle accelerations are computed by interpolating on the gradient of the potential field.
You could potentially do better than N^2 (ish) by approximating. Quantize the world into a grid, add up mass at each grid point. Then each particle is affected by this weighted grid, which is made up of a relatively small number of pseudo-particles.
Of course, you'd still need to fall back to N^2 for the nearest particles, or the accuracy would be unusably bad. But it's a much smaller N.
