Let $\Omega \subset \mathbb R^n$ be a smooth, bounded domain. Then there exist sequences $(\lambda_k)_{k \in \mathbb N} \subset [0, \infty)$ and $(\varphi_k)_{k \in \mathbb N} \subset C^\infty(\overline \Omega)$ satisfying $\int_\Omega \varphi_k \varphi_j = \delta_{kj}$, $$\begin{cases} -\Delta \varphi_k = \lambda_k \varphi_k, \\ \partial_\nu \varphi_k = 0 \end{cases}$$ and for each $f \in L^2(\Omega)$ we have $$ f = \sum_{k=1}^\infty \varphi_k \int_\Omega f \varphi_k \quad \text{in $L^2(\Omega)$.} $$
I am looking for reference for this fact (to include in my thesis). I know “the” proof is quite similar to the Dirichlet case (which is handeled for example in Evan's book), but I would like to quote a source dealing with Neumann boundary conditions. Any suggestions?
Jost's Partial Differential Equations has a full proof for both Dirichlet and Neumann conditions. In the first English edition you can find it in section 8.5, titled Eigenvalues of Elliptic Operators.