It is classical that if $\Omega \subset \mathbb R^n$ is a bounded domain with smooth boundary, then the equation
\begin{equation}\label{pf2}
\begin{aligned}
\begin{cases}
\partial^2_{t}u- \Delta u=0\,\quad &\text{on $(0,T)\times \Omega$},
\\
u=f\,\quad &\text{on $(0,T)\times \partial \Omega$,}\\
\end{cases}
    \end{aligned}
\end{equation} 
admits a unique solution $u\in C(0,T;L^2(\Omega))\cap C^1(0,T;H^{-1}(\Omega))$ with $\partial_\nu u \in H^{-1}((0,T)\times \partial \Omega)$.

I want to derive these estimates via a slightly unorthodox method, namely by considering spectral decomposition in the spatial variables. Indeed, let us write $\{\phi_k\}_{k=1}^{\infty}$ and $\{\lambda_k\}_{k=1}^{\infty}$ for the Dirichlet eigenfunctions and eigenvalues of $-\Delta$ respectively.

Writing $u(t,x)= \sum_{k=1}^{\infty} u_k(t)\phi_k(x)$, it is trivial to see that
$$\partial^2_t u_k(t) + \lambda_k u_k(t) = -\int_{(0,T)\times \partial \Omega} f(t,x)\,\partial_\nu \phi_k(x)\,dx$$

How can I conclude the classical properties of $u$ above from this representation.