Boundary values of boundary value problems Let $M$ be a manifold with smooth boundary. We can consider the Dirichlet or the Neumann problem on $M$. Let $(\phi_k)$ be an orthonormal basis of eigenfunctions to the Dirichlet problem and let $(\psi_k)$ and orthonormal basis of eigenfunctions to the Neumann problem.
What can be said about the functions $\partial_{\mathbf{n}}\phi_k|_{\partial M}$ and $\psi_k|_{\partial M}$? For example, does a subsequence of them constitute an orthonormal basis of $L^2(\partial M)$? Are they solutions to some PDE?
 A: Consider the Dirichlet eigenfunctions. Let $\Delta \phi_k + \lambda_k^2\phi_k=0$, where $0<\lambda_1\leq \lambda_2\leq \lambda_3\dots$, and $\displaystyle u_k = \frac{\partial \phi_k}{\partial n}\bigr|_{\,\partial M}$. Here is a quotation from  Bäcker, Fürstberger, Schubert, and Steiner when $\dim M=2$:

Our results show that the mean behaviour of the normalized boundary functions is very similar to the mean behaviour of eigenfunctions. The crucial difference between the two sequences of functions $\{\phi_k\}_{k\in\mathbb N}$, $\{u_k\}_{k\in\mathbb N}$ is that the eigenfunctions live on a two–dimensional space whereas the boundary functions live on a one–dimensional space. Since both $u_k$ and $\phi_k$ oscillate roughly with the same de Broglie wave length $2π/\lambda_k$, this leads to an overcompleteness of the set $\{u_k\}_{k\in\mathbb N}$. This statement can be made more explicit $\ldots$
More precisely, this means that for every $w \in C^\infty(\partial M)$ 
  $$ 
  w(s) = \sum_{k\in\mathbb N} \rho(\lambda-\lambda_k) w_k u_k(s) +
  O\left(1/\lambda\right) \tag{54}
$$
  holds with coefficients 
  $$
  w_k := \frac{\pi}{2\lambda_k^2}\int_{\partial M} w(s)u_k(s)\,ds. \tag{55}
$$ 
  This follows from ($\ldots$) by the method of stationary phase (see $\ldots$). Since $\rho$ is a rapidly decreasing function, this means that the boundary functions with spectral parameter $\lambda_k$ in an interval of fixed width around $\lambda$ form a complete set in the limit $\lambda\to\infty$. The number of these states grows like
  $\lambda$, in contrast to the number of all states up to energy $\lambda^2$, which, according to the Weyl formula, grows like $\lambda^2$. Therefore this result gives a quantitative measure of the overcompleteness of the set $\{u_k\}$.

In formulas (54), (55), $s$ is the generic point in the boundary $\partial M$,
while $\rho\in S(\mathbb R)$ has Fourier transform which is even and compactly supported in a small neighborhood of $0$. Also the notation has been slightly adjusted (e.g., $\lambda_k$ in place of $k_n$ which is used by BFSS)
Similar results hold for the Neumann eigenfunctions and in higher dimensions, see e.g. Hassell and Barnett.
