Gauge invariance of pseudo stress-energy tensor of gravitational waves

The pseudo stress-energy tensor of gravitational waves is given by $T_{\mu\nu}^{(\mathrm{G}\mathrm{W})} = \frac{1}{32\pi}\left\langle \partial_{\mu} \bar{h}_{\alpha\beta} \partial_{\nu}\bar{h}^{\alpha\beta} - \frac{1}{2}\partial_{\mu} \bar{h} \partial_{\nu} \bar{h} -\partial_{\beta} \bar{h}^{\alpha\beta} \partial_{\nu} \bar{h}_{\alpha\mu} - \partial_{\beta} \bar{h}^{\alpha\beta} \partial_{\mu} \bar{h}_{\alpha\nu} \right\rangle.$ I want to prove that it is invariant under gauge transformation $h_{\mu\nu} \to h_{\mu\nu} - \partial_{\mu}\xi_{\nu} - \partial_{\nu}\xi_{\mu}.$ The bracket $\langle \rangle$ means taking average over the region of which the scale is several times that of gravitational wave's wavelength. So the terms of total divergence, e.g. $\langle \partial_{\mu} (\cdots) \rangle$ , will vanish. But I have met a problem. Under gauge transformation, $\partial_{\beta} \bar{h}^{\alpha\beta} \partial_{\nu} \bar{h}_{\alpha\mu}$ will produce a term $\partial_{\beta} \bar{h}^{\alpha\beta}\partial_{\nu}\partial_{\alpha}\xi_{\mu}$ . However, there is no other term of the order $O(h)$ would ever contain $\xi_{\mu}$ and so I find it impossible to eliminate it. I have check a few textbooks GR, such as MTW (page 969-970) and Sean Carroll (page 307-310). But they just state that the pseudo stress-energy tensor is gauge invariant but give no detailed proof. And I also failed to find a paper which cover the proof. Any help or recommendation for books/papers is appreciated! Some additional information: In Padmanabhan's book Gravitation , page 422, there is an example illustrating the gauge dependence of the energy of the gravitational waves. Consider a gravitational wave of the form $h_{ab}=f^{\prime\prime}(u)yzk_{a}k_{b}$ , where $u = t-x$ and $k_a = \partial_a u$ . Show that, while this is a valid gravitational wave solution, its energy-momentum tensor vanishes. Argue that it is, however, possible to change the gauge and transform $h_{ab}$ to the form $h_{ab}^{\prime}=f(u)\left(\delta_{a}^{y}\delta_{b}^{z}+\delta_{a}^{z}\delta_{b}^{y}\right)$ when the energy-momentum tensor is nonzero and has the form $t_{ab} \propto f^{\prime 2}k_a k_b$ . I have checked this exercise and find the demanding gauge transformation is $\xi_t = \frac{1}{2}f'(u)yz, \; \xi_x = -\frac{1}{2}f'(u)yz, \;\xi_y = -\frac{1}{2}f(u)z, \;\xi_z = -\frac{1}{2}f(u)y.$ So this example shows that $T_{\mu\nu}^{(\mathrm{G}\mathrm{W})}$ is not gauge independent strictly. Of course, I know the coordinate transformation generated by $\xi_a$ of the form above is not infinitesimal. I think that some additional constraints may be needed for the derivation of gauge invariance of $T_{\mu\nu}^{(\mathrm{G}\mathrm{W})}$ .