Ever since the initial development of general relativity in 1915, one nagging question remained open: How do we best define the energy-content of the gravitational field?

The gravitational field definitely has energy-content. When two objects accelerate toward each other due to gravity, the energy gain is balanced by a loss of energy in the gravitational field itself. This energy can be extracted. This happens in the case of binary pulsars, but also in the case of neutron star or black hole mergers: The energy extracted from the gravitational field is converted into gravitational radiation.

However, this gravitational field energy cannot be localized. Any attempt to pinpoint "where" the gravitational field energy is to be found, is doomed to fail. The principal reason is that locally we can always find a freely falling coordinate system in which gravitational effects vanish.

One of the earliest attempts to establish a quantity representing the energy-content of the gravitational field was the Einstein pseudotensor. As it is a pseudotensor, its value depends on the chosen coordinate reference frame. This has been demonstrated in 1919 when Schrödinger showed [see also my English translation] that it is possible to construct a coordinate system representing Schwarzschild spacetime, in which the Einstein pseudotensor vanishes everywhere. The same year, Bauer showed that the Einstein pseudotensor can be nonvanishing even in flat Minkowski spacetime, in a suitable reference frame.

Schrödinger's calculation can be validated by the following Maxima script:

load(ctensor)$
derivabbrev:true$
ct_coords:x$
depends(R,[x[1],x[2],x[3]])$
lg:zeromatrix(4,4)$
I:ident(4)$
for m thru 3 do for n thru 3 do lg[m,n]:I[m,n]+A*x[m]*x[n]/(R^3*(1-A/R));
lg[4,4]:1-A/R;
cmetric()$
christof(false)$
R:sqrt(x[1]^2+x[2]^2+x[3]^2)$
assume('R>0)$
t:zeromatrix(4,4)$
for m thru 4 do for n thru 4 do for a thru 4 do for b thru 4 do for s thru 4 do for l thru 4 do
  t[s,a]:t[s,a]+1/2*I[s,a]*ug[m,n]*mcs[m,b,l]*mcs[n,l,b]-ug[m,n]*mcs[m,b,a]*mcs[n,s,b];
t:ev(t,'R,diff,(x[1]^2+x[2]^2+x[3]^2)='R^2,eval)$
t:factor(t)$
t,'R,eval;

Similarly, Bauer's result can also be checked:

load(ctensor)$
derivabbrev:true$
ct_coords:x$
lg:diag_matrix(-1/(3*x[1])^(4/3),-(3*x[1])^(2/3)/(1-x[2]^2),-(3*x[1])^(2/3)*(1-x[2]^2),1)$
cmetric()$
christof(false)$
riemann(true)$
I:ident(4)$
t:zeromatrix(4,4)$
for m thru 4 do for n thru 4 do for a thru 4 do for b thru 4 do for s thru 4 do for l thru 4 do
  t[s,a]:t[s,a]+1/2*I[s,a]*ug[m,n]*mcs[m,b,l]*mcs[n,l,b]-ug[m,n]*mcs[m,b,a]*mcs[n,s,b];
factor(t);

As it turns out, Schrödinger and Bauer used a slightly different definition of the Einstein pseudotensor,

\begin{align}
\chi t^\alpha_\sigma = \frac{1}{2} \delta^\alpha_\sigma g^{\mu\nu} \Gamma^\lambda_{\mu\beta} \Gamma^\beta_{\nu\lambda} - g^{\mu\nu} \Gamma^\alpha_{\mu\beta} \Gamma^\beta_{\nu\sigma},
\end{align}

compared to what appears in some of the modern literature, and used on Wikipedia:

\begin{align}
{\displaystyle {t_{\mu }}^{\nu }={\frac {1}{2\kappa {\sqrt {-g}}}}\left(\left(g^{\alpha \beta }{\sqrt {-g}}\right)_{,\mu }\left(\Gamma _{\alpha \beta }^{\nu }-\delta _{\beta }^{\nu }\Gamma _{\alpha \sigma }^{\sigma }\right)-\delta _{\mu }^{\nu }g^{\alpha \beta }\left(\Gamma _{\alpha \beta }^{\sigma }\Gamma _{\sigma \rho }^{\rho }-\Gamma _{\alpha \sigma }^{\rho }\Gamma _{\beta \rho }^{\sigma }\right){\sqrt {-g}}\right)}.
\end{align}

Nonetheless, the same results apply. First, Schrödinger:

load(ctensor)$
derivabbrev:true$
ct_coords:x$
depends(R,[x[1],x[2],x[3]])$
lg:zeromatrix(4,4)$
I:ident(4)$
for m thru 3 do for n thru 3 do lg[m,n]:I[m,n]+A*x[m]*x[n]/(R^3*(1-A/R));
lg[4,4]:1-A/R;
cmetric()$
christof(false)$
R:sqrt(x[1]^2+x[2]^2+x[3]^2)$
assume('R>0)$
t:zeromatrix(4,4)$
for m thru 4 do for n thru 4 do for a thru 4 do for b thru 4 do for s thru 4 do for r thru 4 do
  t[m,n]:t[m,n]+diff(ug[a,b],x[m])*(mcs[a,b,n]-I[n,b]*mcs[a,s,s])-
     I[m,n]*ug[a,b]*(mcs[a,b,s]*mcs[s,r,r]-mcs[a,s,r]*mcs[b,r,s]);
t:factor(ev(t,'R,diff,(x[1]^2+x[2]^2+x[3]^2)='R^2,eval))$
t:ev(t,'R);

Next, Bauer:

load(ctensor)$
derivabbrev:true$
ct_coords:x$
lg:diag_matrix(-1/(3*x[1])^(4/3),-(3*x[1])^(2/3)/(1-x[2]^2),-(3*x[1])^(2/3)*(1-x[2]^2),1)$
cmetric()$
christof(false)$
riemann(true)$
I:ident(4)$
t:zeromatrix(4,4)$
for m thru 4 do for n thru 4 do for a thru 4 do for b thru 4 do for s thru 4 do for r thru 4 do
  t[m,n]:t[m,n]+diff(ug[a,b],x[m])*(mcs[a,b,n]-I[n,b]*mcs[a,s,s])-
     I[m,n]*ug[a,b]*(mcs[a,b,s]*mcs[s,r,r]-mcs[a,s,r]*mcs[b,r,s]);
factor(t);

In light of these results, one may wonder: What's the use of a pseudotensor if its value depends on the choice of coordinates? Well, there is an interesting twist: The pseudotensor by itself may not be a suitable measure of local energy-content, but its integral over all space in asymptotically flat spacetimes, as well as its flux across closed surfaces, can represent generally covariant physical quantities. In other words, the pseudotensor may not tell us how much gravitational field energy is at any particular place, but it can tell us how much gravitational field energy is there in total, and where this energy goes.