Edit: I have changed the proof of the \eqref{cc}$ \implies $\eqref{np} to correct an error in the reasoning pointed out by the asker. The development is now necessarily more complicated but entirely correct. I'd like to thank Prof. Alberto Cialdea for the useful discussion on the topic and the suggestion to use Fredholm theory and the equivalent Neumann problem for the Laplace equation.
What we want to prove is that the following Neumann problem
$$
\color{green}{
\begin{cases}
\Delta \varphi(x)=f(x) & x\in\bar{\Omega}\\
\frac{\partial}{\partial \vec{n}}\varphi(x)=g(x)& x\in\partial\bar{\Omega}
\end{cases}\label{np}\tag{NP}}
$$
is solvable if and only if the following compatibility condition
$$
\color{blue}{
\int\limits_\bar{\Omega}f(x)\mathrm{d}x=\int\limits_{\partial\bar{\Omega}}g(x)\mathrm{d}\sigma_x.
\label{cc}\tag{CC}}
$$
holds (with obvious meaning of the symbols), i.e. \eqref{np}$ \iff $\eqref{cc}. Let's proceed with proving the two opposite implications.
- \eqref{np}$ \implies $\eqref{cc}. Integrating on the domain $\bar{\Omega}$ both members of the Poisson equation $\Delta \varphi=f$ from \eqref{np} we get:
$$
\int\limits_\bar{\Omega}\Delta \varphi(x)\mathrm{d}x=\int\limits_\bar{\Omega}f(x)\mathrm{d}x. \label{1}\tag{1}
$$
Now, remembering that $\Delta=\nabla\cdot\nabla=\operatorname{div}\operatorname{grad}$ and applying the Gauss-Green (divergence) theorem and the Neumann boundary condition from \eqref{np} to the first side of equation \eqref{1}, we have
$$
\begin{split}
\int\limits_\bar{\Omega}\Delta \varphi(x)\mathrm{d}x &=\int\limits_\bar{\Omega}\nabla\cdot\nabla \varphi(x)\mathrm{d}x = \int\limits_{\partial\bar{\Omega}}\nabla \varphi(x)\cdot\vec{n}_x\mathrm{d}\sigma_x \\
&= \int\limits_{\partial\bar{\Omega}}\frac{\partial}{\partial \vec{n}}\varphi(x)\mathrm{d}\sigma_x = \int\limits_{\partial\bar{\Omega}}g(x)\mathrm{d}\sigma_x,
\end{split}
$$
thus
$$
\int\limits_\bar{\Omega}\Delta \varphi(x)\mathrm{d}x=\int\limits_{\partial\bar{\Omega}}g(x)\mathrm{d}\sigma_x, \label{2}\tag{2}
$$
and this implies condition \eqref{cc}.
- \eqref{cc}$ \implies $\eqref{np}. Perhaps the easiest (and classical) way to prove this step is to transform problem \eqref{np} into a Fredholm integral equation of the second kind, and verify that the necessary (and sufficient, see the notes below) condition for its unique solvability is \eqref{cc}. To proceed in this way, we'll reduce problem \eqref{np} to the corresponding one for the Laplace equation and then apply the result known for the Neumann problem for harmonic functions. Let
$$
s(x-y)=
\begin{cases}
\dfrac{1}{2\pi}\log|{x}-{y}| & n=2\\
\\
-\dfrac{1}{(n-2)|\sigma_n}|{x}-{y}|^{2-n} & n\ge 3
\end{cases}
$$
be the fundamental solution supported in $y\in\Omega\subset\Bbb R^n$ of the Laplace operator: if we formally define
$$
\begin{align}
\varphi^\ast(x)&=\varphi(x)-\int\limits_\bar{\Omega}s(x-y)f(y)\mathrm{d}y & x\in\Omega \\
g^\ast(x)&=g(x)-\frac{\partial}{\partial \vec{n}}\int\limits_\bar{\Omega}s(x-y)f(y)\mathrm{d}y & \: x\in\partial\Omega\label{3}\tag{3}
\end{align}
$$
we have that $\varphi(x)$ solves \eqref{np} if and only if $\varphi^\ast(x)$ is harmonic and solves
$$
\color{brown}{
\begin{cases}
\Delta \varphi^\ast(x)=0 & x\in\bar{\Omega}\\
\frac{\partial}{\partial \vec{n}}\varphi^\ast(x)=g^\ast(x)& x\in\partial\bar{\Omega}
\end{cases}\label{hnp}\tag{HNP}}
$$
Now let's search for $\varphi^\ast(x)$ by expressing it as a single layer potential, i.e.
$$
\varphi^\ast(x)= \int\limits_{\partial\bar{\Omega}}s(x-y) \mu(y)\mathrm{d}\sigma_y,\label{4}\tag{4}
$$
and searching for the unknown (charge) density $\mu(x)$, $x\in\partial\bar{\Omega}$.
By applying the classical jump formula for the first derivatives of harmonic functions to \eqref{4} (see for example [2], §22.7, pp. 304-306), i.e.
$$
\lim_{\substack{x\to x_0 \\ x\in \vec{n}_{x_0}^+}} \frac{\partial}{\partial \vec{n}_{x}} \int\limits_{\partial\bar{\Omega}}s(x-y)\mu(y)\mathrm{d}\sigma_y = -\frac{1}{2}\mu(x_0) + \frac{\partial}{\partial \vec{n}_{x_0}} \int\limits_{\partial\bar{\Omega}}s(x_0-y)\mu(y)\mathrm{d}\sigma_y
$$
we get the following Fredholm integral equation for the unknown density $\mu(x)$
$$
g^\ast(x)= - \frac{1}{2}\mu(x) + \int\limits_{\partial\bar{\Omega}} \frac{\partial}{\partial \vec{n}_{x}} s(x-y)\mu(y)\mathrm{d}\sigma_y \label{5}\tag{5}
$$
whose sufficient condition for solvability, by Fredholm's alternative theorem (see for example [2], §16.3, pp. 225-227), is
$$
\color{purple}{
\int\limits_{\partial\bar{\Omega}}g^\ast(x)\mathrm{d}\sigma_x=0
\label{hcc}\tag{HCC}}
$$
Now let's analyze condition \eqref{hcc}: substituting \eqref{3} in it, we have
$$
\begin{split}
\int\limits_{\partial\bar{\Omega}}g^\ast(x)\mathrm{d}\sigma_x &= \int\limits_{\partial\bar{\Omega}}\bigg[g(x) - \frac{\partial}{\partial \vec{n}}\int\limits_\bar{\Omega}s(x-y)f(y)\mathrm{d}y\bigg]\mathrm{d}\sigma_x\\
& = \int\limits_{\partial\bar{\Omega}} g(x) \mathrm{d}\sigma_x -\int\limits_{\partial\bar{\Omega}}\frac{\partial}{\partial \vec{n}}\int\limits_\bar{\Omega}s(x-y)f(y)\mathrm{d}y\,\mathrm{d}\sigma_x\\
& \text{and applying the divergence theorem}\\
& = \int\limits_{\partial\bar{\Omega}} g(x) \mathrm{d}\sigma_x -\int\limits_\bar{\Omega} \Delta\int\limits_\bar{\Omega}s(x-y)f(y)\mathrm{d}y\,\mathrm{d}x\\
& = \int\limits_{\partial\bar{\Omega}} g(x) \mathrm{d}\sigma_x -\int\limits_\bar{\Omega} f(x) \mathrm{d}x
\end{split}\label{6}\tag{6}
$$
and the proof is then finished.
Final notes on the method of proof of the implication \eqref{cc}$ \implies $\eqref{np}.
- The same very same method can be used to prove directly the equivalence \eqref{cc}$ \iff $\eqref{np}: as alluded above, condition \eqref{hcc} (and his equivalent condition \eqref{cc} for Poisson's equation), is de facto a necessary and sufficient condition for the solvability of integral equation \eqref{5}. This is the approach followed for example by V.S. Vladimirov for the Laplace equation ([2], §23.5, pp. 315-318) which, however, does not explicitly deals with Poisson's equation, i.e. does not provide the ansatz \eqref{3} nor the development \eqref{6}.
- The same method can also be used to prove the existence and uniqueness of the weak solution (and thus the weak formulation) of the Neumann problem \eqref{np} even for the more general divergence form equation,
$$
\operatorname{div}\big(k(x)\nabla \varphi(x)\big)-a(x)\varphi(x)=f(x)
$$
obviously without any reference nor use of potential theory. To see why, it is sufficient to remember that the weak formulation of problem \eqref{np} is simply a set of integral identities to be satisfied by the solution $\varphi$: the details can be found in reference [1], chapter IV, §1.2, 1.6.
[1] V. P. Mikhailov (1978), Partial differential equations, Translated from the Russian by P.C. Sinha. Revised from the 1976 Russian ed., Moscow: Mir Publishers, p. 396 MR0601389, Zbl 0388.3500.
[2] V. S. Vladimirov (1971)[1967], Equations of mathematical physics, Translated from the Russian original (1967) by Audrey Littlewood. Edited by Alan Jeffrey, (English), Pure and Applied Mathematics, Vol. 3, New York: Marcel Dekker, Inc., pp. vi+418, MR0268497, Zbl 0207.09101.
