Solution (source code)

= Solution

Expand the <Einstein tensor> about the <Minkowski metric>, keeping a second-order <metric perturbation> as well:
$$
g=\eta+\varepsilon h+\varepsilon^2j+O(\varepsilon^3),\qquad G(g)=\varepsilon G^{(1)}[h]+\varepsilon^2\bigl(G^{(1)}[j]+G^{(2)}[h,h]\bigr)+O(\varepsilon^3).
$$
Here $G^{(2)}[h,h]$ includes all quadratic terms from the <Einstein tensor>, including the inverse-metric corrections in its contractions. In vacuum, the first-order <Linearized Einstein equations> give $G^{(1)}[h]=0$, while the next order gives $G^{(1)}[j]=-G^{(2)}[h,h]$. Define the <quadratic gravitational effective stress-energy tensor> by
$$
\boxed{t_{\mu\nu}=-\frac1{8\pi}G^{(2)}_{\mu\nu}[h,h].}
$$
It is symmetric because the <Einstein tensor> is symmetric, and is quadratic in $h$ and its derivatives. The linearized <contracted Bianchi identity> is the off-shell identity $\partial^\mu G^{(1)}_{\mu\nu}[j]=0$. The second-order vacuum equation therefore implies
$$
\boxed{\partial^\mu t_{\mu\nu}=0.}
$$
One can obtain the same conclusion directly by expanding the full <contracted Bianchi identity> to second order: every correction to the flat divergence multiplies $G^{(1)}[h]$ and vanishes on the first-order vacuum solution. Thus conservation holds for vacuum solutions of <linearized gravity>, rather than for an arbitrary off-shell $h$. Relative to the fixed Minkowski background, $t$ is a tensor under background changes of coordinates, but this does not make it invariant under the perturbative gauge freedom.