Is there a singularity theorem in higher-dimensional Newtonian gravity?

This should follow from the following Virial-type computation.

For convenience we assume all particles have the same mass; this is not essential.

Let $x_i$ denote the position vector of the $i$th particle, then Newton's law of universal gravitation, suitably normalized, gives

$$ \ddot{x_i} = (d-2) \sum_{j \neq i} \frac{ x_j - x_i}{|x_j - x_i|^d } $$

This gives

$$ x_i \cdot \ddot{x_i} = (d-2) \sum_{j \neq i} \frac{ x_j \cdot x_i - |x_i|^2 }{|x_j - x_i|^d} $$

summing over all particles we get

$$ \sum_{i} x_i \cdot \ddot{x_i} = (d-2) \sum_{i,j, i\neq j} \frac{ x_j \cdot x_i - |x_i|^2 }{|x_j - x_i|^d} = -\frac{d-2}{2} \sum_{i,j, i\neq j} |x_i - x_j|^{2-d} \tag{1}$$

Conservation of energy, on the other hand, states that

$$ - \sum_{i,j, i\neq j} |x_i - x_j|^{2-d} + \sum_i |\dot{x_i}|^2 = C < 0 \tag{2}$$

(by negative energy assumption). Summing (1) and (2) and using that $d \geq 4$ so that $(d-2)/2 \geq 1$ we get

$$ \frac{d^2}{dt^2} \sum_{i} |x_i|^2 < 0 $$

Which shows that $\sum |x_i|^2$ must impact $0$ at some finite time (either future or past).

(Note: by "dropping particles at infinity" and using time-symmetric property of Newtonian mechanics we see that there exists trajectory where there is only one time of impact in the past, and none in the future.)

Following the lines Willie followed, but allowing for unequal masses $m_i$, set $I(x) = \langle x, x \rangle $ where $\langle v, w \rangle = \Sigma m_i v_i \cdot w_i$ is the so-called mass metric built so that the kinetic energy is $K =\frac{1}{2}\langle v, v \rangle$ with $v = \dot x$. Here $x$ and $v$ lie in N copies of a d-dimensional Euclidean space. Then Newton's equations read $\ddot x = - \nabla V (x)$ where the gradient $\nabla$ is with respect to the mass metric. Suppose now that the potential $V$ is homogeneous of degree $-\alpha$. Using Euler's identity $-\alpha V (x) = \langle x, \nabla V (x) \rangle$ a fun elementary computation yields the virial identity $\ddot I = 4H +(2 \alpha - 4) V$, where $H = K +V$ is the total energy, which is conserved. (This identity is called by mathematical celestial mechanicians the Lagrange-Jacobi identity.) Now for $d$-dimensional ``Newtonian' gravity we have $V < 0$ and $\alpha = d-2$ so that Lagrange-Jacobi reads $\ddot I = 4H +(8-4d)U$, with $U = -V > 0$. Conclusion: for $d > 4$ and $H \le 0$ we have $\ddot I < 0$ and thus all such solutions either begin and end in total collision $I =0$, (or have singularities before these total collision times), as you conjectured. (Asides. If $d< 4$ then the virial identity implies that $\ddot I > 0$ for $H \ge 0$ and consequently bounded solutions have $H < 0$. The case $d =4$ is a wonderful bounding case between these extremes. At the critical dimension $4$ solutions must lie on spheres $I = const.$ and have energy $H =0$ in order to be bounded.)