Let $\lambda_j$ be the eigenvalues of the Neumann Laplacian on a bounded domain $\Omega$ with eigenvectors $\varphi_j$, which we know are smooth.
Given a function $u \in H^{\frac 12}(\Omega)$ with $\int_\Omega u = 0$, define $$v(x,y) = \sum_{k=1}^\infty e^{-y\lambda_k^{\frac 12}}(u,\varphi_k)_{L^2}\varphi_k(x).$$
The claim is $v \in H^1(\Omega \times (0,\infty))$ and $v \in C^\infty((0,\infty);H^1(\Omega))$.
I don't see why $v(y) \in H^1(\Omega)$ for a.a. $y$. I know the basis functions $\varphi_k$ are smooth, but that can't be the full reason. If it were, then by this logic, can we not just say that every $u \in L^2(\Omega)$ is in $H^1(\Omega)$, since $u=\sum (u,\varphi_k)\varphi_k$, and $\nabla u = \sum (u,\varphi_k)\nabla \varphi_k$?