It's possible to define the energy stored in the gravitational field for Newtonian gravity using a similar approach to the way we define energy stored in the electric field. For the electric field, we start from Coulomb's Law (in CGS units) for the force on a point charge q_{1} which arises from another point charge q_{2}:
\vec{F}_{1} = \frac{q_{1} q_{2}}{{\Vert \vec{r}_1 - \vec{r}_{2} \Vert}^3} (\vec{r}_{1} - \vec{r}_{2})
For a bunch of point charges labeled q_{i}, where i goes from 1 to N, the force on q_{1} adds linearly, so we get
\vec{F}_{1} = \sum^{N}_{j=2} \frac{q_{1} q_{j}}{{\Vert \vec{r}_1 - \vec{r}_{j} \Vert}^3} (\vec{r}_{1} - \vec{r}_{j})
and in general, the force acting on the ith charge is given by
\vec{F}_{i} = \sum^{N}_{j \neq i} \frac{q_{i} q_{j}}{{\Vert \vec{r}_i - \vec{r}_{j} \Vert}^3} (\vec{r}_{i} - \vec{r}_{j})
We define a potential energy function for each charge q_{i}, using \vec{F}_{i} = - q_{i} \vec{\nabla} \Phi_{i} for some scalar potential functions \Phi_{i}, and a vector field \vec{E}_{i} = \frac{1}{q_{i}} \vec{F}_{i} corresponding to the electric field that each charge sees. Note that the potential energy function and electric fields are different for each charge using this method, since we assume that point charges do not see their own fields. This is not strictly true (if I remember correctly, point charges interacting with their own fields are what causes the radiation reaction force for accelerating charges), but the problem will go away when we take the continuum limit and radiation reaction can be ignored here anyway. We find that the potential energy functions have a nice expression
\Phi_{i} = \sum_{j \neq i}^{N} \frac{q_{j}}{\Vert \vec{r}_{i} - \vec{r}_{j} \Vert}
This allows us to calculate the amount of energy necessary to assemble a collection of point charges by bringing them in one at a time from infinity, that is, the potential energy of the system. This is given by
E_{potential} = \frac{1}{2} \sum^{N}_{i=1} q_{i} \Phi_{i} = \frac{1}{2} \sum^{N}_{i \neq j} \frac{q_{i} q_{j}}{\Vert \vec{r}_{i} - \vec{r}_{j} \Vert}
The real electric field is defined using the idea of a test charge. We imagine adding a small charge Q to the system at \vec{r} and define the electric field \vec{E} and electric potential \Phi as the electric field and scalar potential seen by our test charge. This gives the well-known formulae
\vec{E}(\vec{r}) = \sum^{N}_{j=1} \frac{q_{j}}{{\Vert \vec{r} - \vec{r}_{j} \Vert}^3} (\vec{r} - \vec{r}_{j})
and
\Phi(\vec{r}) = \sum_{j =1}^{N} \frac{q_{j}}{\Vert \vec{r} - \vec{r}_{j} \Vert}
With these formulae in hand, we can verify that Gauss' law holds,
- \nabla^2 \Phi = \vec{\nabla} \cdot \vec{E} = 4 \pi \rho(\vec{r})
where we have defined the charge density \rho(\vec{r}) = \sum^{N}_{j = 1} q_{j} \delta^{(3)}(\vec{r} - \vec{r}_{j}) and made use of the identity \nabla^2 \left( \frac{1}{r} \right) = - 4 \pi \delta^{(3)} (\vec{r}). When we take the continuum limit, we no longer have to worry about particles interacting with their own fields, which allows us to write
E_{potential} = \frac{1}{2} \int_{space} \rho(\vec{r}) \Phi(\vec{r}) d^3r = - \frac{1}{8 \pi} \int_{space} \Phi(\vec{r}) \nabla^2 \Phi(\vec{r}) d^3r
We can do the three-dimensional equivalent of integrating by parts using Green's theorem to get
E_{potential} = \frac{1}{8 \pi} \int_{space} \vec{\nabla} \Phi(\vec{r}) \cdot \vec{\nabla} \Phi(\vec{r}) d^3r + \text{boundary term}
and because we are integrating over all space, our boundary term vanishes for localized charge distributions. Using -\vec{\nabla} \Phi = \vec{E}, we discover
E_{potential} = \frac{1}{8 \pi} \int_{space} \vec{E} \cdot \vec{E} d^3r
which motivates the assertion that the density of energy stored in the electric field is given by \frac{1}{8 \pi} \vec{E} \cdot \vec{E}. The same derivation applies to the Newtonian gravitational field mutatis mutandis. In this case, Gauss' law must be changed to
- \nabla^2 \Phi_{g} = \vec{\nabla} \cdot \vec{E}_{g} = - 4 \pi G \rho_{g}
where \Phi_{g} is the gravitational potential, \vec{E}_{g} is the gravitational analogue of an electric field, and \rho_{g} is the mass density. There is an extra minus sign as compared to before because gravity is attractive between two positive masses, whereas two positive charges repel. Nonetheless, we can still define
E_{potential} = \frac{1}{2} \int_{space} \rho_{g}(\vec{r}) \Phi_{g}(\vec{r}) d^3r
This time, however, when we apply Green's theorem, we are short one minus sign as compared to last time, and we discover in the end
E_{potential} = - \frac{1}{8 \pi G} \int_{space} \vec{E}_{g} \cdot \vec{E}_{g} d^3r
motivating us to define gravitational energy density as -\frac{1}{8 \pi G} \vec{E}_{g} \cdot \vec{E}_{g}. This negative energy density is disturbing- what could it possibly mean? Fortunately for us, the question is only academic, because real gravity is described by general relativity. General relativity has its own problems with defining gravitational energy density since gravitational effects are contained in the spacetime metric instead of in a field living on the spacetime manifold, but there are clever constructions one can make to get around this difficulty. As far as Newtonian gravity goes, you really do get a negative energy density if you try to do things in the same way you would in electromagnetism. I've found a paper with more information about this, as well as a way to get around this by including a self-interaction of the gravitational field, here:
www.iop.org/EJ/article/0143-0807/28/6/016/ejp7_6_016.pdf[/URL]