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) & xinbar{Omega}\
frac{partial}{partial vec{n}}varphi(x)=g(x)& xinpartialbar{Omega}
end{cases}label{np}tag{NP}}
$$
is solvable if and only if the following compatibility condition
$$
color{blue}{
intlimits_bar{Omega}f(x)mathrm{d}x=intlimits_{partialbar{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:
$$
intlimits_bar{Omega}Delta varphi(x)mathrm{d}x=intlimits_bar{Omega}f(x)mathrm{d}x. label{1}tag{1}
$$
Now, remembering that $Delta=nablacdotnabla=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}
intlimits_bar{Omega}Delta varphi(x)mathrm{d}x &=intlimits_bar{Omega}nablacdotnabla varphi(x)mathrm{d}x = intlimits_{partialbar{Omega}}nabla varphi(x)cdotvec{n}_xmathrm{d}sigma_x \
&= intlimits_{partialbar{Omega}}frac{partial}{partial vec{n}}varphi(x)mathrm{d}sigma_x = intlimits_{partialbar{Omega}}g(x)mathrm{d}sigma_x,
end{split}
$$
thus
$$
intlimits_bar{Omega}Delta varphi(x)mathrm{d}x=intlimits_{partialbar{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}{2pi}log|{x}-{y}| & n=2\
\
-dfrac{1}{(n-2)|sigma_n}|{x}-{y}|^{2-n} & nge 3
end{cases}
$$
be the fundamental solution supported in $yinOmegasubsetBbb R^n$ of the Laplace operator: if we formally define
$$
begin{align}
varphi^ast(x)&=varphi(x)-intlimits_bar{Omega}s(x-y)f(y)mathrm{d}y & xinOmega \
g^ast(x)&=g(x)-frac{partial}{partial vec{n}}intlimits_bar{Omega}s(x-y)f(y)mathrm{d}y & : xinpartialOmegalabel{3}tag{3}
end{align}
$$
we have that $varphi(x)$ solves eqref{np} if and only if $varphi(x)$ is harmonic and solves
$$
color{brown}{
begin{cases}
Delta varphi^ast(x)=0 & xinbar{Omega}\
frac{partial}{partial vec{n}}varphi^ast(x)=g^ast(x)& xinpartialbar{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)= intlimits_{partialbar{Omega}}s(x-y) mu(y)mathrm{d}sigma_y,label{4}tag{4}
$$
and searching for the unknown (charge) density $mu(x)$, $xinpartialbar{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{xto x_0 \ xin vec{n}_{x_0}^+}} frac{partial}{partial vec{n}_{x}} intlimits_{partialbar{Omega}}s(x-y)mu(y)mathrm{d}sigma_y = -frac{1}{2}mu(x_0) + frac{partial}{partial vec{n}_{x_0}} intlimits_{partialbar{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) + intlimits_{partialbar{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}{
intlimits_{partialbar{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}
intlimits_{partialbar{Omega}}g^ast(x)mathrm{d}sigma_x &= intlimits_{partialbar{Omega}}bigg[g(x) - frac{partial}{partial vec{n}}intlimits_bar{Omega}s(x-y)f(y)mathrm{d}ybigg]mathrm{d}sigma_x\
& = intlimits_{partialbar{Omega}} g(x) mathrm{d}sigma_x -intlimits_{partialbar{Omega}}frac{partial}{partial vec{n}}intlimits_bar{Omega}s(x-y)f(y)mathrm{d}y,mathrm{d}sigma_x\
& text{and applying the divergence theorem}\
& = intlimits_{partialbar{Omega}} g(x) mathrm{d}sigma_x -intlimits_bar{Omega} Deltaintlimits_bar{Omega}s(x-y)f(y)mathrm{d}y,mathrm{d}x\
& = intlimits_{partialbar{Omega}} g(x) mathrm{d}sigma_x -intlimits_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.