We know that given the divergence and curl of a vector field (and appropriate boundary conditions) it is possible to construct a unique vector field in $\mathbb R^3$. The specific problem I am thinking about is related to the PDE
$$\operatorname{div} F = g,$$
where $F \colon \mathbb R^n \to \mathbb R^n$ is a vector field and $g \colon \mathbb R^n \to \mathbb R$ is a scalar field, and $\operatorname{div}$ is the $n$-dimensional generalization of the divergence given by
$$\operatorname{div} F = \frac{dF_{i}}{dx_{i}}$$
(summation implied). What additional pieces of information are necessary to uniquely specify $F$ given the function $g$ (we know the answer is the curl of $F$ in 3D)?
[Math] Higher Dimensional Generalization of Helmholtz Theorem
partial differential equationsvector analysis
Related Solutions
Shuhao Cao has provided a strong answer in terms of the usual language of differential forms and tensor calculus. I present here an alternative in terms of geometric calculus which I hope will be enlightening.
As I related in my answer to your other question, geometric calculus can look somewhat different from differential forms, but most of the basic concepts are the same. For instance, there is an exterior derivative:
$$\mathrm d\omega \equiv \nabla \wedge \omega$$
And a coderivative, or interior derivative,
$$\delta \omega \equiv \nabla \cdot \omega$$
These are simply how these concepts are denoted in geometric calculus. The preference of using them over $d, \delta$ is that the geometric product makes meaningful the statement
$$\nabla \omega = \nabla \wedge \omega + \nabla \cdot \omega$$
which, written in terms of $\mathrm d\omega, \delta \omega$ might otherwise come off as patent gibberish. But the use of multivectors in geometric algebra and calculus makes this a perfectly valid mathematical statement. It allows us to separate what would otherwise be some tensor $\partial_i \omega_j$ into two distinct, meaningful parts. If $\omega$ is a $k$-vector, then $\nabla \cdot \omega$ is a $k-1$-vector and $\nabla \wedge \omega$ is a $k+1$-vector, and both pieces can be kept track of in the equation using the clifford algebra framework.
Again, treating $\nabla$ as an operator in and of itself is useful in geometric calculus because it has a Green's function:
$$\nabla G = \delta$$
where here, $\delta$ should be understood as the Dirac delta function. The form of the Green's function in 3d is $G(r) = \frac{r}{4\pi |r|^3}$, and you can ultimately relate it to having the same form in higher dimensions, just dependent on the hyperarea of the unit sphere and so on. Having a Green's function makes $\nabla$ invertible through integration in a way that the exterior and interior derivatives by themselves are not.
With that in mind, the Helmholtz decomposition can be expressed using geometric calculus quite simply. First, you ask about an identity that is used in the derivation. This is actually a general result:
$$\nabla(\nabla \omega) = \nabla^2 \omega = \nabla \wedge (\nabla \cdot\omega) + \nabla \cdot (\nabla \wedge \omega)$$
This is the geometric calculus analogue of your identity, and it unifies the Laplace-Beltrami and Laplace-de Rham operators that Shuhao Cao refers to into a single, scalar operator $\nabla^2$.
Now we can go about proving the Helmholtz decomposition. Let $E$ be a vector field. Let $M$ be some region with boundary $\partial M$. Let us presume that there exist a scalar function $\phi$ and a bivector field $A$ such that
$$E = \nabla \wedge \phi + \nabla \cdot A = \nabla \psi$$
such that $\nabla \wedge F = 0$. There is nothing particularly lacking in generality with these assumptions, as I will show. Consider the following integral:
$$\oint_{\partial M} G(r-r') dS' \, \psi(r') = (-1)^{n} \int_{M} \delta(r-r') \, dV' \psi(r') + \int_{M} G(r-r') \, dV' \, E(r')$$
This is where keeping track of grades, as they're used in clifford algebra, is important: even if $\nabla \wedge A$ were nonzero, we could separate it out of this equation by considering that 3-vector term alone. So I have done nothing that violates generality (yes, I've assumed some things about integrability and smoothness and such; we'll go with all those necessary things being applicable, for otherwise you wouldn't even be able to do the problem).
What this use of the fundamental theorem does is give us a formula for $\psi = \phi + A$, and since separate grades must always independently obey such an equation, we can write down the equation as
$$\begin{align*} \phi(r) &= \int_M G(r-r') \cdot E(r') \, |dV'|- \oint_{\partial M} |dS'| \left[ G(r-r') \cdot \hat n' \phi(r') - [G(r-r') \wedge \hat n'] \cdot A(r') \right] \\ A(r) &= \int_M G(r-r') \wedge E(r') \, |dV'| - \oint_{\partial M} |dS'| \left[G(r-r') \wedge \hat n' \phi(r') + G(r-r') \cdot \hat n' \, A(r')\right] \end{align*}$$
We're not actually going to use these formulas; I'm merely arguing that, given these potentials defined on some surface, the functions themselves can be constructed as well.
What we really want to do is look at the Green's function for $\nabla^2$ instead. There does exist such a function (which I will call $H$ such that $\nabla H = G$). In 3d, $H(r) = (4\pi |r|)^{-1}$, so this is exactly what you're used to seeing in expressions of the Helmholtz decomposition.
Again, our preferred hammer for deriving integral theorems is the fundamental theorem. Use $H$ and $E$ and observe that
$$\nabla \oint_{\partial M} H(r-r') \, dS' \, E(r') = (-1)^n (i E) + \nabla \int_M H(r-r') \, dV' \, \nabla' E(r')$$
(If you can't see how this is so, just pull the $\nabla$s into the integrals and remember $\nabla H = G$.)
But remember that $E = \nabla \psi$. The invertibility of $\nabla$ means that we can, through some integrals, essentially cancel $\nabla$ to get
$$(-1)^n \phi(r) = \oint_{\partial M} H(r-r') \, \hat n' \cdot E(r') |dS'| - \int_{M} H(r-r') \nabla' \cdot E(r') \, |dV'|$$
and the bivector potential $A$ follows by replacing dots with wedges. This is the Helmholtz decomposition. You need to know the Green's function for the Laplacian. You need to know the dimension of the space. But most importantly, since you seem hung up on this point, you need to know the tangential components of the vector field, not just the normal component.
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
The appropriate generalization you want is to consider the vector field $F$ and the scalar field $g$ as $(n-1)-$ and $n$-forms in the deRham complex. That is $$ F \in \Omega^{n-1}(\mathbb{R}^n) = \Omega^{n-1} \otimes C^\infty(\mathbb{R}^n) , $$ $$F = f_1 \ dx_2 \wedge ... \wedge dx_n + f_2 \ dx_1 \wedge dx_3 \wedge ... \wedge dx_n + $$ $$ f_3 \ dx_1 \wedge dx_2 \wedge dx_4 \wedge ... \wedge dx_n + ...+ f_n \ dx_1 \wedge ... \wedge dx_{n-1} \wedge dx_n $$ and $$ g \in \Omega^n(\mathbb{R}^n) $$ $$ g = G \ dx_1 \wedge dx_2 \wedge ... \wedge dx_n $$
Framing it in these terms the equation $$ div( F) = g $$ is $$dF = g $$ where d is the external derivative $$ d : \Omega^j(\mathbb{R}^n) \rightarrow \Omega^{j+1}(\mathbb{R}^n) $$
Now your question becomes: given $ g$, what $F$ satisfy $dF = g$. The answer is certainly there is no unique function/n-1 form $F$ which satisfies this equation. Given any such solution $F_0$, then $F = F_0 + dE$ is also a solution where $E$ is a $(n-2)-$form, because $d^2 = 0$. Further, given any such $g$ there must exist at least one $F$ because $g$ is a closed $n$-form, meaning $dg = 0$ by the Poincare theorem, every closed $n$-form is also exact, meaning there exists such a $(n-1)-$form $F$ such that $dF = g$. (In the language of deRham cohomology, $H^n_{DR}(\mathbb{R}^n) = 0$.)
The more interesting question is this: what happens if we look now not at the deRham complex on $\mathbb{R}^n$ but on subsets of $\mathbb{R}^n$ that are topologically non-trivial? Or on $n$-dimensional manifolds? Are there then $n$-forms $g$ for which there are no $F$ such that $dF = g$?
In other words, the existence of closed $n$-forms that are not exact will depend on the topology of the underlying space, and vice versa. It's a great topic.