Is Helmholtz decomposition inherently a non-local operation? Helmholtz decomposition, the process for splitting a vector field into parts which have vanishing divergence and curl, plays a central role in our ability to quantize the electromagnetic field because it permits us to separate the gauge invariant solenoidal part of the vector potential, $\mathbf{A}_{\mathrm{sol}}$, from its gauge dependent irrotational part, $\mathbf{A}_{\mathrm{irrot}}$. The formulae for these components of the field can be defined from the Wikipedia article as
\begin{align}
   \mathbf{A}_{\mathrm{irrot}}(\mathbf{x}) &= -\frac{\nabla}{4\pi} \int_V \frac{\nabla' \cdot \mathbf{A}(\mathbf{x}')}{\left|\mathbf{x}-\mathbf{x}'\right|} \operatorname{d}V' + \frac{\nabla}{4\pi} \oint_{\partial V} \hat{n}'\cdot \frac{\mathbf{A}(\mathbf{x}')}{\left|\mathbf{x}-\mathbf{x}'\right|} \operatorname{d}S' \\
     &= -\frac{\nabla}{4\pi} \int_V \mathbf{A}(\mathbf{x}') \cdot \nabla \frac{1}{\left|\mathbf{x}-\mathbf{x}'\right|} \operatorname{d}V',\ \mathrm{and}\\
   \mathbf{A}_{\mathrm{sol}}(\mathbf{x}) &= \mathbf{A}(\mathbf{x}) - \mathbf{A}_{\mathrm{irrot}}(\mathbf{x}).
\end{align}
Now, for a field that already satisfies $\nabla\times\mathbf{A}=0$ the first operator is the identity and the second is $0$, which both act locally, and the same is true for the swapped versions. Equivalently, if a field is either solenoidal or irrotational then that can be determined locally by application of the divergence and curl and seeing which one is zero. Given a mixed field, though, the operator that constructs its irrotational component looks highly non-local, especially if $V$ is expanded to cover all of space.
Can this determination be made locally? That is, does going the other direction with $V$, towards zero size by, say, partitioning $V$ into compact sub-volumes, work? My instinct is to say that the partition is inherently non-local because of the existence of magnetic scalar potential techniques.
 A: Yes, it is non-local. 
A useful way to think about the decomposition is that $\mathbf{A}_{\text{irrot}}$ is the gradient vector field closest to $\mathbf{A}$, with the closeness of two fields measured by the average norm-squared difference of the vectors in the fields. This averaging operation makes the decomposition non-local.
For example, if we have a vector field 
$$
\mathbf{A}(x,y) = \begin{bmatrix} 0 \\ x \end{bmatrix},
$$
then the closest gradient is the constant field equal to $\mathbf{A}$ at the center of the domain, and thus is not locally calcuable from $\mathbf{A}$.
A: Yes, Helmholtz decomposition is inherently a non-local process because it is only unambiguous when the integrals are performed over all of a simply connected space that contains all relevant sources of the field (where a "source" is defined as any region in the field where the divergence or curl is non-zero). As the usage of the magnetic scalar potential demonstrates, it fails on sub-spaces because of the ambiguity inherent in the classification. In standard Helmholtz decomposition, divergenceless and solenoidal are taken as synonymous, likewise with divergent and irrotational (curl free). Fields that have both zero divergence and zero curl are assumed to be excluded by the boundary condition at infinity that the fields vanish there.
When you limit your considerations to finite sized regions, though, the assumption breaks down because a source that is outside of your bounded region can produce a field that is both divergenceless and irrotational everywhere in the region of interest. For that reason, local categorization is only unambiguous if it is done via positive properties. In other words, a vector field can be locally described unambiguously by breaking it down into three components:
\begin{align}
    \mathbf{F}(\mathbf{x}) & = \mathbf{F}_{\mathrm{div}}(\mathbf{x}) + \mathbf{F}_{\mathrm{curl}}(\mathbf{x}) + \mathbf{F}_{\mathrm{har}}(\mathbf{x}),\ \mathrm{where} \\
    \nabla\cdot\mathbf{F}_{\mathrm{div}}(\mathbf{x}) &\neq 0\ \quad\mathrm{somewhere}, \\
    \nabla\times\mathbf{F}_{\mathrm{curl}}(\mathbf{x}) &\neq 0\ \quad\mathrm{somewhere}, \\
    \nabla\cdot \mathbf{F}_{\mathrm{har}}(\mathbf{x}) & = 0\ \quad\mathrm{everywhere},\ \mathrm{and} \\
    \nabla\times \mathbf{F}_{\mathrm{har}}(\mathbf{x}) & = 0\ \quad\mathrm{everywhere}.
\end{align}
The subscript of $\mathbf{F}_{\mathrm{har}}(\mathbf{x})$ is meant to be short for "harmonic" in analogy to harmonic functions because it can always be expressed in the region of interest as the negative gradient of a harmonic function. You could also call the harmonic term $\mathbf{F}_{\mathrm{ext}}(\mathbf{x})$, since it can also be modeled as being produced by sources external to the region of interest (when that region has finite size).
In this categorization, the reconstruction would work like this:
\begin{align}
    \mathbf{F}_{\mathrm{div}}(\mathbf{x}) & = \frac{1}{4\pi} \int_V \frac{\mathbf{x}-\mathbf{x}'}{\left|\mathbf{x}-\mathbf{x}'\right|^3} \nabla'\cdot \mathbf{F}(\mathbf{x}') \operatorname{d}V' \\
    \mathbf{F}_{\mathrm{curl}}(\mathbf{x}) & = \frac{1}{4\pi} \int_V \nabla\times\frac{\nabla'\times \mathbf{F}(\mathbf{x}')}{|\mathbf{x}-\mathbf{x}'|} \operatorname{d}V' \\
    \mathbf{F}_{\mathrm{har}}(\mathbf{x}) & = \mathbf{F}(\mathbf{x}) - \mathbf{F}_{\mathrm{div}}(\mathbf{x}) - \mathbf{F}_{\mathrm{curl}}(\mathbf{x}).
\end{align}
What this exercise makes clear is if you try to take the limit as $V\rightarrow 0$, the region will either have some kind of delta function source that leaves the field ill defined there (e.g. point charge, line charge, or sheet charge) or the integrals that define the divergent and solenoidal components vanish, leaving only a locally harmonic field. Thus, the classification is only useful when done in regions of finite size, and only when the region of interest is the entire possible space. 
That is why the Wikipedia formulae quoted in the question do not have an ambiguity in their derivation, they have an unstated assumption that no sources beyond the boundary contribute to the local field. As an example, consider a uniform field in the region of interest. Was it produced by one (or more) infinite charged sheet(s) (divergent source), or infinite current sheet(s) (solenoidal source)? Assuming the region of interest is a cube with faces aligned parallel and perpendicular to the field, the Wikipedia sources would have a combination of charged sheets on the faces the field pierces, and a current sheet circulating around the other four.
Point being, the process is only unambiguous when sources are only allowed in a closed region of interest.
