I'm trying to better understand the theory of second-order elliptic PDEs on a smooth compact region $\Omega \in \mathbb{R}^N$ with boundary conditions. For the Dirichlet case, there seems to be a rich literature on criteria for existence of a solution, culminating in a Fredholm alternative for the system
\begin{align} Lu - \lambda u &= f \qquad \text{in } \Omega \\ u &= 0 \qquad \text{on } \partial \Omega \end{align} where $L$ is a second-order elliptic differential operator and $f$ is smooth. This system is well-posed and has a (unique) solution $u$ for all $f$ if $\lambda$ is not an eigenvalue of $L$. If $\lambda$ is an eigenvalue of $L$, then the system has a solution if and only if $f$ is orthogonal to all eigenfunctions of $L$ with eigenvalue $\lambda$. See e.g. Sec. 6.2.3 in Evans (2nd ed., 2010).
I am interested in the analogous Neumann system: \begin{align} Lu - \lambda u &= f \qquad \text{in } \Omega \\ \hat{n} \cdot \nabla u &= g \qquad \text{on } \partial \Omega \end{align} where $\hat{n}$ is the unit normal perpendicular to the boundary and $f,g$ are smooth. I have found statements in some special cases. For example if $L = - \nabla^2$, $\lambda = 0$ then this system has a solution if and only if
\begin{equation} \int_\Omega f = \int_{\partial \Omega} g \end{equation}
but I don't see a way to generalize easily.
Is there some general-purpose theorem for existence of a solution to this system? Does a weak solution to this system have any guaranteed interior regularity properties? I am most interested in the case $L = -\nabla^2$, $\lambda \neq 0$, so I would be satisfied with a result that applies only to this case.