Let $(M^2,g)$ be a compact Riemannian surface with boundary and let $L = \Delta_g + q$ be a Schrödinger operator, where $\Delta_g = -\operatorname{div} \nabla$ is the Laplacian with respect to the metric $g$ and $q \in C^\infty(M)$. Also let $\nu$ denote the outward unit normal to $\partial M$ and $f \in C^\infty(\partial M)$.

I am looking for either a reference or a proof that the eigenvalue problem

$$\cases{Lu = \lambda u, \quad \text{on } M \\ \frac{\partial u}{\partial \nu} = fu, \quad \text{on } \partial M}$$

has a discrete real spectrum of the form

$$\lambda_1 < \lambda_2 \leq \cdots \leq \lambda_k \leq \cdots \nearrow \infty.$$

Moreover, is there a variational (minmax) characterization of these eigenvalues?