The existence and uniqueness of the pure Neumann boundary value problem for smooth data can be proved using the Green representation formula, explicitly.
First notice that for $-\Delta \Phi = \delta(x)$, by Green's second identity and convolution formula: for $x\in \Omega$
$$
\int_{\Omega} \Big(u(y) \Delta \Phi(y-x) - \Phi(y-x)\Delta u(y)\Big) dy = \int_{\partial \Omega}\left(u(y)\frac{\partial \Phi}{\partial n}(y-x) - \Phi(y-x)\frac{\partial u}{\partial n}(y) \right) dS(y),
$$
we have
$$
u(x) = -\int_{\Omega} \Phi(y-x)\Delta u(y) \, dy + \int_{\partial \Omega}\left(\color{red}{\Phi(y-x)}\frac{\partial u}{\partial n}(y) - u(y)\color{blue}{\frac{\partial \Phi}{\partial n}(y-x)} \right) dS(y).
$$
For Dirichlet problem, a correction function is constructed so that red term vanishes.
For Neumann problems, Green function is constructed so that $u(x)$ on each point is unique up to a constant. Some people prefer to set:
$$
\left\{\begin{aligned}
-\Delta \Phi &= \delta(x) \quad \text{in }\Omega,
\\
\frac{\partial \Phi}{\partial n} &= \frac{1}{|\partial \Omega|} \quad \text{on }
\partial\Omega,
\end{aligned}\right.
$$
where absolute value just means the Lebesgue measure of codimension $1$, i.e., surface area, so that the constant difference is the average of $u$ on $\partial \Omega$. I myself prefer to set something like:
$$
\left\{\begin{aligned}
-\Delta \Phi &= \delta(x) - \frac{1}{|\Omega|}\quad \text{in }\Omega,
\\
\frac{\partial \Phi}{\partial n} &=0 \quad \text{on }
\partial\Omega,
\end{aligned}\right.\tag{$\star$}
$$
so that you have:
$$
u(x) = \frac{1}{|\Omega|}\int_{\Omega}u(y)dy -\int_{\Omega} \Phi(y-x)\Delta u(y) \, dy + \int_{\partial \Omega}
\Phi(y-x)\frac{\partial u}{\partial n}(y)dS(y).
$$
This means the solution is also unique up to a constant, that constant is the average of $u$ in $\Omega$, in other words, once we fix this average, $u$ is pinned down to just one function. The existence of the Green function is another story, which relies on functional analysis, e.g. see here.
The unique of $u$ up to a constant can be either proved by maximum principle or energy method, but for Helmholtz decomposition does not address the uniqueness, we don't have to address this issue here.
Now for your equation, we know that the compatibility condition is satisfied because:
$$
\int_{\Omega}\Delta u = \int_{\Omega}\nabla \cdot u = \int_{\partial \Omega} \frac{\partial u}{\partial n} dS = \int_{\partial \Omega} E\cdot n\,dS.
$$
Consider $u$ is smooth and $\int_{\Omega} u = 0$, then
$$
u(x) = -\int_{\Omega} \Phi(y-x) (\nabla_y \cdot E(y)) \, dy + \int_{\partial \Omega}
\Phi(y-x) (E(y)\cdot n) dS(y).
$$
Here are another two proofs for smooth vector field. The Neumann boundary approach is good for less regular vector field (say the one has component in some Sobolev space), for smooth vector field, using vector potentials would be more preferable (at least for me).
The first one is using a pointwisely hold formula if $\Omega$ is convex with respect to a point $x_0 \in \mathbb{R}^3$:
$$
E = \nabla \phi + A,
$$
where
$$
A = -(x - x_0)\times \int^1_0 t \nabla \times E\big(x_0 + t(x- x_0)\big)\,dt ,
$$
and
$$
\phi = (x - x_0)\cdot \int^1_0 tE\big(x_0 + t(x- x_0)\big).
$$
Though this does not bear the exact same form, but the formula gives you more physical intuition. For reference please see here.
The second one is using the Dirichlet problem's green function. Now we want to solve a vector Poisson equation:
$$
\left\{\begin{aligned}
-\Delta F &= E\quad \text{in }\Omega,
\\
F &=0 \quad \text{on }
\partial\Omega,
\end{aligned}\right.
$$
Then for the Dirichlet Green function:
$$
\left\{\begin{aligned}
-\Delta \Phi &= \delta(x) \quad \text{in }\Omega,
\\
\Phi &=0 \quad \text{on } \partial\Omega.
\end{aligned}\right.
$$
Use the Green representation formula component-wise, we have
$$
F_i(x) = -\int_{\Omega} \Phi(y-x)\Delta F_i(y) \, dy - \int_{\partial \Omega} F_i(y) \frac{\partial \Phi}{\partial n}(y-x) dS(y),
$$
and this is
$$
F_i(x) = \int_{\Omega} \Phi(y-x)E_i(y) \, dy.
$$
Now use the vector Laplacian's identity for smooth vector fields
$$
E(x) = -\Delta F(x) = \nabla \times \nabla \times F - \nabla \nabla \cdot F
\\
=\nabla_x \times \left(\nabla_x \times \int_{\Omega} \Phi(y-x)E (y) \, dy\right) - \nabla _x\left(\nabla_x \cdot \int_{\Omega} \Phi(y-x)E(y) \, dy\right),
$$
where the subscript $x$ means taking gradient/div/curl w.r.t. $x$. Further moving the derivative into the integral:
$$
E(x) = \nabla_x \times \left(\int_{\Omega} \nabla_x \times \Big(\Phi(y-x)E (y)\Big) \, dy\right) - \nabla_x\left(\int_{\Omega} \nabla_x \cdot \Big(\Phi(y-x)E(y)\Big) \, dy\right)
\\
= \nabla_x \times \left(\int_{\Omega} \nabla_x \Phi(y-x)\times E (y) \, dy\right) - \nabla_x\left(\int_{\Omega} \nabla_x \Phi(y-x)\cdot E(y) \, dy\right)
\\
= - \nabla_x \times \left(\int_{\Omega} \nabla_y \Phi(y-x)\times E (y) \, dy\right) + \nabla_x\left(\int_{\Omega} \nabla_y \Phi(y-x)\cdot E(y) \, dy\right)
\\
= \nabla_x \times \left(\int_{\Omega} \Big( \Phi(y-x) \nabla_y \times E (y) - \nabla_y \times (\Phi(y-x)E(y))\Big) \, dy\right)
\\
+ \nabla_x\left(\int_{\Omega} \Big( -\Phi(y-x) \nabla_y \cdot E (y) + \nabla_y \cdot (\Phi(y-x)E(y))\Big) \, dy\right)
\\
= \nabla_x \times \left(\int_{\Omega} \Phi(y-x) \nabla_y \times E (y)\,dy - \int_{\partial \Omega} \Phi(y-x)n\times E(y)\,dS(y) \right)
\\
+ \nabla_x\left(-\int_{\Omega} \Phi(y-x) \nabla_y \cdot E (y) \, dy +\int_{\partial \Omega}\Phi(y-x)E(y)\cdot n\,dS(y) \right).
$$
So that: in $\Omega$
$$
E = \nabla \phi + \nabla \times A,
$$
where
$$
\phi(x) = -\int_{\Omega} \Phi(y-x) \nabla_y \cdot E (y) \, dy +\int_{\partial \Omega}\Phi(y-x)E(y)\cdot n\,dS(y),
\\
A(x) = \int_{\Omega} \Phi(y-x) \nabla_y \times E (y)\,dy - \int_{\partial \Omega} \Phi(y-x)n\times E(y)\,dS(y).
$$
Best Answer
This question has already a bounty winning answer: however, since I found it interesting, upvoted it the very first day it was posted, and found unusual difficulties in finding a nice answer, I nevertheless want to share my own researches and thoughts.
The only reference I am aware of which gives a proof of this result, is the textbook by Mikhailov [1]. Precisely, he deals with the regularity problem for the first boundary problem (Dirichlet problem) and the second boundary problem (Neumann problem) for the Poisson equation in ยง2.3, pp. 216-226. Mikhailov proves the result directly: he first solves the problem in the case of homogeneous boundary conditions, proving the following theorem:
Theorem 4 ([1], p. 217). If $f\in H^k(\Omega)$ and $\partial\Omega\in C^{k+2}$ for certain $k\ge 0$, then the generalized solutions $u(x)$ of the first and second boundary-value problems with homogeneous boundary condition for the Poisson equation belong to $H^{k+2}(\Omega)$ and satisfy (in the case of the second boundary problem it is assumed that $\int_\Omega u \mathrm{d}x=0$) the inequality $$ \Vert u\Vert_{ H^{k+2}(\Omega)}\le C\Vert f\Vert_{ H^{k}(\Omega)} $$ where the constant $C>0$ does not depend on f.
By using this result, he extends the regularity theorem 4 above to the case of non homogeneous boundary conditions ([2], p. 226) by reducing non homogeneous boundary conditions to homogeneous ones. He explicitly does it only for the Dirichlet problem but the same method nevertheless works for the Neumann problem: let's see this. Consider a solution $u(x)$ of the problem \eqref{eqlocalvar-Neumann} above and a function $\Phi\in H^{k+2}(\Omega)$ whose normal derivative on $\partial\Omega$ is $g\in H^{k+3/2}(\partial\Omega)$, i.e. $$ \frac{\partial\Phi}{\partial\nu}=g\;\text{ on }\;\partial\Omega\label{1}\tag{1} $$ Define $w=u-\Phi$: then, for every $v\in H^{1}(\Omega)$, $$\label{GeneralizedNeumann}\tag{GN} \begin{split} \int_{\Omega} \nabla w \cdot \nabla v \, \mathrm{d}x&=\int_{\Omega} \nabla u \cdot \nabla v \, \mathrm{d}x - \int_{\Omega} \nabla \Phi \cdot \nabla v \, \mathrm{d}x \\ &= \int_\Omega f v \, \mathrm{d}x + \int_{\partial \Omega}gv \, \mathrm{d}\sigma(x) - \int_{\Omega} \nabla \Phi \cdot \nabla v \, \mathrm{d}x \\ &=\int_\Omega f v \, \mathrm{d}x + \int_{\partial \Omega}gv \, \mathrm{d}\sigma(x) - \int_{\Omega} \nabla\cdot(v\nabla\Phi) \, \mathrm{d}x + \int_{\Omega} v\Delta\Phi \, \mathrm{d}x \\ &=\int_\Omega \big[\,f + \Delta\Phi\big]v\, \mathrm{d}x + \int_{\partial \Omega}gv \, \mathrm{d}\sigma(x) - \int_{\partial\Omega} \nu \frac{\partial\Phi}{\partial\nu} \, \mathrm{d}\sigma(x), \end{split}\label{2}\tag{2} $$ and substituting \eqref{1} in \eqref{2} we get $$ \int_{\Omega} \nabla w \cdot \nabla v \, \mathrm{d}x=\int_\Omega f_1 v \, \mathrm{d}x $$ where $f_1=f+\Delta\Phi\in H^k(\Omega)$. This means that $w$ solves the Neumann problem \eqref{eqlocalvar-Neumann} with homogeneous boundary conditions (i.e. $g\equiv 0$), and by applying theorem 4 we have $$ w=u-\Phi\in H^{k+2}(\Omega) \iff u=w+\Phi\in H^{k+2}(\Omega) $$
Notes
[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.