Define the function $f_i:\mathbb{R}^3\rightarrow\mathbb{R}^3$, $i\in\{1,2,3\}$, by $f_i(\boldsymbol{x})=\delta(\boldsymbol{x-x_0})\boldsymbol{e}_i$ where $\delta$ is the Dirac Delta function and $\boldsymbol{e}_i$ is the $i^{th}$ Cartesian unit basis vector.
Does $f_i$ have a Helmholtz decomposition
$$f_i(\boldsymbol{x})=\nabla\psi+\boldsymbol{\zeta}$$ $$\nabla\cdot\boldsymbol{\zeta}=0$$
where $\nabla\psi$ and $\boldsymbol{\zeta}$ are of the same smoothness as $f_i$ (or smoother), and if so, what are $\psi$ and $\boldsymbol{\zeta}$?
I'm interested because I have a PDE of the form
$$L\boldsymbol{u}=\boldsymbol{g}, \boldsymbol{x}\in\Omega$$ $$\sigma(u)\cdot\boldsymbol{n}=0\in\partial\Omega$$
and I'd like to try to find a Green's function using a Helmholtz decomposition of $L\boldsymbol{u}$.
The answer is yes. It is known that $\frac{1}{4\pi|\boldsymbol{x-x_0}|}$ is the Green's function of Laplace's equation in 3-D. Thus $\Delta(\frac{1}{4\pi|\boldsymbol{x-x_0}|})=\delta(\boldsymbol{x-x_0})$. Due to the property that the vector Laplacian of a vector is equal to the vector whose components in Cartesian coordinates are the scalar Laplacians of the Cartesian coordinates of the original vector. Thus $\Delta(\frac{\boldsymbol{e_i}}{4\pi|\boldsymbol{x-x_0}|})=\delta{(\boldsymbol{x-x_0)}}\boldsymbol{e_i}=f_i$.
The vector Laplacian also has the property $\Delta y=\nabla(\nabla\cdot\boldsymbol{y})-\nabla\times(\nabla\times y)$, so we conclude that $f_i$ has the following Helmholtz decomposition:
$$f_i=\nabla\left(\partial_{x_i}\frac{1}{4\pi|{\boldsymbol{x-x_0}}|}\right)-\nabla\times\left(\nabla\times\frac{\boldsymbol{e_i}}{4\pi|\boldsymbol{x-x_0}|}\right)$$