There's two ways of doing it, I implemented them both in my PhD and didn't have a ton of fun doing it.
(a) There's a method that works well for monopolar sources (gravitational + electrostatic particles) called the Barnes-Hut method. You effectively divide space up into a quadtree (2D) or octree (3D), and in each cell work out the center of mass / total charge. You make particles in "nearby" cells (using a distance criterion that can be adjusted to speed up/slow down the simulation in a trade off with accuracy) interact directly, and far away cells you just use the center of mass to work out the interaction between any given 'far' particle and the particles in that cell. The method is O(N log N) but in practice, this is 'good enough' for many applications.
(b) uses a more rigorous technique called the Fast Multipole Method which is O(N), where rather than just using the center of mass or sum of charges, you expand the potential from particles out into higher order components which captures the distribution of particles within each cell. This also means you can capture more complex potentials. The downside is that this is a nightmare to implement in comparison to the Barnes-Hut method. Each cell has it's own multipole expansion, and it is 'transferred' to work out the additive contribution to every 'far' cell, calculating a 'local' expansion. Typically people use the most compact representation of these potential expansions which uses Lagrange polynomials, but this is a pain.
(a) There's a method that works well for monopolar sources (gravitational + electrostatic particles) called the Barnes-Hut method. You effectively divide space up into a quadtree (2D) or octree (3D), and in each cell work out the center of mass / total charge. You make particles in "nearby" cells (using a distance criterion that can be adjusted to speed up/slow down the simulation in a trade off with accuracy) interact directly, and far away cells you just use the center of mass to work out the interaction between any given 'far' particle and the particles in that cell. The method is O(N log N) but in practice, this is 'good enough' for many applications.
(b) uses a more rigorous technique called the Fast Multipole Method which is O(N), where rather than just using the center of mass or sum of charges, you expand the potential from particles out into higher order components which captures the distribution of particles within each cell. This also means you can capture more complex potentials. The downside is that this is a nightmare to implement in comparison to the Barnes-Hut method. Each cell has it's own multipole expansion, and it is 'transferred' to work out the additive contribution to every 'far' cell, calculating a 'local' expansion. Typically people use the most compact representation of these potential expansions which uses Lagrange polynomials, but this is a pain.