For completeness, let me just mention an alternative approach that makes it clearer that the energy is stored
locally in the field rather than in individual particles. The analogue in electromagnetism forms part of the foundation for the electromagnetic stress energy tensor.
Let us start from the expression for the total potential energy
$$
U = \frac{1}{2} \int \rho(\vec x) \phi(\vec x) d^3x
$$
where ##\rho(\vec x)## is the mass density and ##\phi(\vec x)## the gravitational potential such that ##\vec g = -\nabla \phi##. The Poisson equation for gravity ensures that the potential satisfies
$$
\nabla \cdot \vec g = - \nabla^2 \phi = - 4\pi G \rho
$$
or, equivalently,
$$
\rho = \frac{1}{4\pi G} \nabla^2 \phi.
$$
Inserting this into the expression for the potential energy leads to
$$
U = \frac{1}{8\pi G} \int [\nabla^2 \phi] \phi\, d^3x.
$$
Applying the divergence theorem and assuming the potential to vanish sufficiently fast as ##r \to \infty## now leads to
$$
U = - \frac{1}{8\pi G} \int (\nabla \phi)^2 d^3 x = -\frac{1}{8\pi G} \int \vec g(\vec x)^2 d^3x.
$$
We can thus associate a gravitational field ##\vec g(\vec x)## with a local
potential gravitational energy density*
$$
u_g(\vec x) = -\frac{1}{8\pi G} \vec g(\vec x)^2.
$$
Knowing the functional form of the gravitational field ##\vec g(r) = - \vec e_r GM(r)/r^2## with ##M(r) = 4\pi r^3 \rho/3## for ##r < R## and ##M(r) = 4\pi \rho R^3/3## for ##r > R## for the homogeneous sphere then leads to
\begin{align*}
U &= - \frac{1}{8\pi G} \int_0^\infty \frac{G^2 M(r)^2}{r^4} 4\pi r^2 dr \\
&= -\frac{G}{2} \int_0^\infty \frac{M(r)^2}{r^2} dr \\
&= - \frac{G \rho^2}{2} \left[\int_0^R \frac{16 \pi^2 r^4}{9} dr + \int_R^\infty \frac{16\pi^2 R^6}{9r^2} dr\right] \\
&= - \frac{8 \pi^2 G\rho^2}{9} \left(\left[\frac{r^5}{5}\right]_0^R - \left[\frac{R^6}{r}\right]_R^\infty \right) \\
&= - \frac{8\pi^2 G\rho^2 R^5}{9} \frac{6}{5} = - \frac{16 \pi^2 G \rho^2 R^5}{15}.
\end{align*}
This should now look familiar.
Note: The minus sign is essential. The gravitational potential energy of the assembled system is less than zero because you would have to do work to disassemble it.
Edit:
*: Compare this with the
electrostatic potential energy stored in an electric field
$$
u_e(\vec x) = \frac{1}{2} \varepsilon_0 \vec E(\vec x)^2.
$$
The constants are a bit different because of how ##\varepsilon_0## and constants appear in Gauss' law, but the essential argument is exactly the same. (Just replace ##\vec g \to \vec E## and ##-1/4\pi G \to \varepsilon_0##.)
Edit 2:
Including also the magnetic field, the EM energy density is given by
$$
u_{em} = \frac{1}{2}\left( \varepsilon_0 \vec E^2 + \frac{1}{\mu_0} \vec B^2\right).
$$
This is the time-time-component of the electromagnetic stress energy tensor. The time-space-components are given by the Poynting vector ##\vec S = \vec E \times \vec B/\mu_0## and the space-space-components are the components of the Maxwell stress tensor.