[proofplan]
We differentiate the energy with respect to time and use the regularity of $u$ to justify passing the derivative under the spatial integral. The spatial-gradient term is converted by [integration by parts](/theorems/210) into a Laplacian term plus a boundary flux term. The boundary flux vanishes because the Dirichlet condition holds for all times, so its time derivative is zero on $\partial U$. The remaining volume integral is exactly $\int_U \partial_t u(\partial_t^2u-\Delta u)\,d\mathcal{L}^n$, which vanishes by the [wave equation](/page/Wave%20Equation).
[/proofplan]
[step:Differentiate the energy under the spatial integral]
Fix $t \in (0,T)$. Since $u \in C^2(\overline{U}\times[0,T])$ and $\overline{U}\times[0,T]$ is compact, the functions $\partial_t u$, $\partial_t^2u$, $\nabla u$, and $\nabla\partial_tu$ are continuous and bounded on $\overline{U}\times[0,T]$. Hence the classical differentiation-under-the-integral rule applies to the energy integral. We obtain
\begin{align*}
E'(t) = \int_U \partial_t u(x,t)\,\partial_t^2u(x,t)\, d\mathcal{L}^n(x) + \sum_{i=1}^n \int_U \partial_{x_i}u(x,t)\,\partial_{x_i}\partial_tu(x,t)\, d\mathcal{L}^n(x).
\end{align*}
[guided]
Fix $t \in (0,T)$. The first task is to justify differentiating $E(t)$, not merely to perform the formal computation. Because $u \in C^2(\overline{U}\times[0,T])$, the functions $\partial_t u$, $\partial_t^2u$, $\partial_{x_i}u$, and $\partial_{x_i}\partial_tu$ are continuous on the compact set $\overline{U}\times[0,T]$ for every $i \in \{1,\dots,n\}$. In particular, they are bounded and integrable over $U$ with respect to $\mathcal{L}^n$.
The integrand of the energy is
\begin{align*}
\frac{1}{2}\left(|\partial_t u(x,t)|^2 + |\nabla u(x,t)|^2\right).
\end{align*}
For each fixed $x \in U$, differentiating this expression with respect to $t$ gives
\begin{align*}
\partial_t u(x,t)\,\partial_t^2u(x,t) + \sum_{i=1}^n \partial_{x_i}u(x,t)\,\partial_{x_i}\partial_tu(x,t).
\end{align*}
The boundedness just noted provides the domination needed for the classical differentiation-under-the-integral rule on the bounded domain $U$. Therefore
\begin{align*}
E'(t) = \int_U \partial_t u(x,t)\,\partial_t^2u(x,t)\, d\mathcal{L}^n(x) + \sum_{i=1}^n \int_U \partial_{x_i}u(x,t)\,\partial_{x_i}\partial_tu(x,t)\, d\mathcal{L}^n(x).
\end{align*}
[/guided]
[/step]
[step:Integrate the spatial-gradient term by parts]
Let $\nu:\partial U\to\mathbb{R}^n$ denote the outward unit normal field on $\partial U$, and let $\nu_i$ denote its $i$th component. Let $\mathcal{H}^{n-1}$ denote the $(n-1)$-dimensional [Hausdorff measure](/page/Hausdorff%20Measure) restricted to $\partial U$, which is the surface measure appearing in the boundary term. For each $i \in \{1,\dots,n\}$, the functions $x\mapsto \partial_{x_i}u(x,t)$ and $x\mapsto \partial_tu(x,t)$ are $C^1(\overline{U})$. Applying the classical integration-by-parts formula on the smooth bounded domain $U$ gives
\begin{align*}
\int_U \partial_{x_i}u(x,t)\,\partial_{x_i}\partial_tu(x,t)\, d\mathcal{L}^n(x) = \int_{\partial U} \partial_tu(x,t)\,\partial_{x_i}u(x,t)\,\nu_i(x)\, d\mathcal{H}^{n-1}(x) - \int_U \partial_tu(x,t)\,\partial_{x_i}^2u(x,t)\, d\mathcal{L}^n(x).
\end{align*}
Summing over $i$ yields
\begin{align*}
\sum_{i=1}^n \int_U \partial_{x_i}u(x,t)\,\partial_{x_i}\partial_tu(x,t)\, d\mathcal{L}^n(x) = \int_{\partial U} \partial_tu(x,t)\,\nabla u(x,t)\cdot\nu(x)\, d\mathcal{H}^{n-1}(x) - \int_U \partial_tu(x,t)\,\Delta u(x,t)\, d\mathcal{L}^n(x).
\end{align*}
[/step]
[step:Use the Dirichlet condition to eliminate the boundary flux]
For each $x \in \partial U$, define the boundary time trace $\gamma_x:[0,T]\to\mathbb{R}$ by $\gamma_x(s)=u(x,s)$. The boundary condition says $\gamma_x(s)=0$ for every $s\in[0,T]$, so $\gamma_x'(t)=\partial_tu(x,t)=0$ for every $t\in(0,T)$. Therefore the boundary integrand vanishes pointwise on $\partial U$:
\begin{align*}
\partial_tu(x,t)\,\nabla u(x,t)\cdot\nu(x) = 0.
\end{align*}
Hence
\begin{align*}
\int_{\partial U} \partial_tu(x,t)\,\nabla u(x,t)\cdot\nu(x)\, d\mathcal{H}^{n-1}(x) = 0.
\end{align*}
[/step]
[step:Apply the wave equation to show the energy derivative is zero]
Substituting the integration-by-parts identity and the vanishing boundary term into the expression for $E'(t)$ gives
\begin{align*}
E'(t) = \int_U \partial_tu(x,t)\,\partial_t^2u(x,t)\, d\mathcal{L}^n(x) - \int_U \partial_tu(x,t)\,\Delta u(x,t)\, d\mathcal{L}^n(x).
\end{align*}
Combining the two volume integrals,
\begin{align*}
E'(t) = \int_U \partial_tu(x,t)\left(\partial_t^2u(x,t)-\Delta u(x,t)\right)\, d\mathcal{L}^n(x).
\end{align*}
Since $u$ solves $\partial_t^2u-\Delta u=0$ in $U\times(0,T)$, the integrand is identically zero on $U$. Thus
\begin{align*}
E'(t)=0
\end{align*}
for every $t\in(0,T)$.
[/step]
[step:Conclude that the energy is constant on the whole time interval]
The preceding step shows that $E$ is differentiable on $(0,T)$ and satisfies $E'(t)=0$ for every $t\in(0,T)$. Hence $E$ is constant on $(0,T)$. Since the integrand defining $E(t)$ depends continuously on $t$ and is bounded on the compact set $\overline{U}\times[0,T]$, the map $E:[0,T]\to\mathbb{R}$ is continuous. Therefore the constant value on $(0,T)$ extends to the endpoints, and
\begin{align*}
E(t)=E(0)
\end{align*}
for every $t\in[0,T]$.
[/step]