Given any vector field $X$ in Minkowski space, we can use Stokes' theorem to derive $$ \int_{t_0}^{t_1}\text{div} X \text d V_g =\int_{\{t=t_0\}} g(X,\partial_t)-\int_{\{t=t_1\}} g(X,\partial_t), $$ where the first integral refers to integration over the slab between the two planes of constant time at $t_0$ and $t=1$. I have tripled check this calculation and convinced myself that the RHS is indeed correct; we subtract the flow through the bottom slice from the flow through the top slice. As a sanity check, if $X=T:=\partial_t$, then the flow through a section of the bottom plane is negative, as intuition would dictate, although if this is wrong, please let me know.
My problem comes when I plug in an energy current $X^a = Q^a{}_{b} T^b$ where $Q$ is the stress-energy tensor for some function $u$: $$ Q_{ab}=\nabla_a u\nabla_b u-\frac 12 g_{ab}\nabla^c u\nabla_c u. $$ In this case, say in spherical coordinates, I have $$ X =\begin{bmatrix} -\frac 12 |\partial u|^2 \\ u_t u_r \\ \frac 1{r^2}u_t u_\theta \\ \frac 1{r^2\sin^2 \theta}u_t u_\varphi \end{bmatrix}. $$ Using the coordinate formula for the divergence I am able to verify that $\text{div} X = u_t \square u = u_t(\Delta u - u_{tt})$, but $Q(X,\partial_t)=\frac 12 |\partial u|^2$, so that the RHS of $$ 2\int_{t_0}^{t_1}u_t\square u \text d V_g =\int_{\{t=t_0\}}|\partial u|^2 -\int_{\{t=t_1\}} |\partial u|^2, $$ has the opposite sign from the correct energy estimate. After going through some textbooks I am pretty sure the formula for the divergence theorem is correct, so I imagine the mistake must be somewhere in the calculations for the energy current $X$. Since I get the correct LHS, I think the expression for $X$ is correct, so maybe I am missing some detail with the inner product, but I do not see how. Where am I making a mistake?
EDIT: I am just now realizing that I do obtain the correct energy identity because my wave operator is the negative of the conventional one $u_{tt} - \Delta u$. However, I have seen in several notes (e.g. these ones on slide 78) that the regular energy identity (i.e. the one for the conventional wave operator) holds even when they have explicitly defined $$ \square u = g^{\mu \nu} \nabla_\mu \nabla_\nu u $$ do the authors actually mean $\square u= - g^{\mu \nu} \nabla_\mu \nabla_\nu u$?