Density and dimensionality of zeros in inverse square force fields of randomly distributed sources in (at least) 1, 2 and 3 dimensions? Background: In this answer to Are there places in the Universe without gravity? in Astronomy SE I did a quick finite 2D calculation for 20 random sources to see if there was at least one zero, and without rigor convinced myself that there might always be some finite density of zeros. See image below which shows log10 normalize force magnitude, the script is in the linked answer.
Question: For a random distribution of discrete inverse square force sources (e.g. stars in space) of finite density, what is the density and dimensionality of zeros in the net force field relative to the density of sources? Please address the 1, 2 and 3 dimension cases at least.

Results of a very quick exploration, looking for zero-dimensional zeros in 2D.
The circa -14 values are simply a result of the cutoff in the minimization. The original script can be found in the linked answer. The labeled minima are the result of minimization from the lowest grid point value. I haven't made an attempt yet to look for all zeros.

 A: This is an almost complete answer. There is currently one small gap, with regards to non-degeneracy, but I am relatively certain that this can be fixed:
tl;dr If there are $n$ masses, then in general, there will be at least $n-1$ zeroes (in 1d precisely $n-1$ zeroes)  of the field, all of them isolated, within the convex hull of these masses and none of them stable.
Within the convex hull: This is the easy to prove independent of dimension. Assume that there is a hyperplane with all masses on one side and you on the other. Then  the force from each mass has a component pulling you towards that plane, so the total force will do the same and cannot be zero. Such a plane exists precisely for the points not in the convex hull.
For the rest, one should think of zeroes as critical points of the potential, which I am going to call $\Phi$. Then assuming they are not degenerate (which they generally aren't) the so called Morse theory tells us about their number. Also note that you can convert all masses into finite minima without creating additional critical points, since close to any given mass its $1/r$-potential dominates all other contributions, so you can replace it at that level with a finite potential well with a single minimum.
1 dimension: That case one can do without resorting to more advanced tools and it is somewhat illustrative of the idea: First note that $\frac{d^2}{dx^2}\frac{-1}{|x|} < 0$ for $x \neq 0$, so $\frac{d^2}{dx^2} \Phi < 0$ and thus every critical point is a maximum. Now its easy to see that on the real line minima and maxima need to alternate. So between the $n$ minima resulting from the masses there are precisely $n-1$ maxima, which are all the critical points.
2 dimensions: Assume for now that all critical points are non-degenerate, by which one means that the Hessian $D^2 \Phi$ has only non-zero eigenvalues. From this you automatically get that these critical points are isolated (apply the inverse function theorem to $D \Phi$), which by a compactness argument also implies that there are only finitely many of them.
Next we finally need some Morse theory: The index of a non-degenerate critical point is defined as the number of negative eigenvalues of $D^2\Phi$. Roughly speaking, this is the number of independent directions in which $\Phi$ is maximal. So a minimum has index $0$, a maximum has index equal to the dimension of the space and all the saddle-points have an index somewhere in between. Now denote by $C_\gamma$ the number of critical points of index $\gamma$.
The details are a bit complicated, but the fundamental theorem of Morse theory now tells us that $\sum_{\gamma} (-1)^\gamma C_\gamma = \xi(M)$, where $\xi(M)$ is the Euler-characteristic of $M$. In our case $M= \mathbb{R}^d$ which has characteristic $1$. Assuming also that the only minima are at the masses, we then get
$$C_2 - C_1 + n = 1 \Rightarrow C_1 \geq C_1 - C_2 = n-1$$
so there are at least $n-1$ critical points of index $1$ plus one more for each of index $2$. (Putting all masses in a line gives you none of index $2$, aranging them into a grid, gives you a lot of them, but there might be an upper bound)
3 dimensions: Before we start with Morse theory, lets do some partial differential equations. In 3 dimensions $\frac{1}{|r|}$ solves Laplace's equation, so in other words we have $\Delta \Phi = 0$ away from the masses. The maximum principle for Laplace's equation then tells us that $\Phi$ is either constant (which it isn't) or has no local maxima or minima (The masses don't count, as the equation no longer holds there). So if all critical points are non-degenerate, we get
$$ -C_3 + C_2 - C_1 + C_0 = C_2-C_1 + n = 1 $$
same as before.
Finally we need to clean up some assumptions:
3d Critical points are non-degenerate:
I am reasonable certain that this holds, by some PDE-argument that I am missing, but technically this is a gap.
2d Critical points are non-degenerate:
Embed $\mathbb{R}^2$ as a plane into $\mathbb{R}^3$ and put all masses on that plane. Then by symmetry, the normal derivative to that plane $\partial_n \Phi$ is zero, so critical points in 2d correspond to those in 3d and similarly, the normal direction is an eigenvector. Then the other eigenvectors are in the plane, so the eigenvalues of those in 3d correspond to eigenvalues in 2d, so if the 3d critical point is non-degenerate, the 2d critical point is as well.
There are no additional 2d minima: By the same embedding argument and the convex hull discussion at the top, critical points are minimal in normal direction to the plane. So if they would be minimal in the plane, they would be local minima in 3d, which contradicts the maximum principle.
