As John Hughes already mentioned, we require $\nabla \cdot \vec J=0$. Under that restriction, we proceed.
Since the curl of the gradient is zero ($\nabla \times \nabla \Phi=0$), then if
$$\nabla \times \vec B =\mu_0 \vec J$$
for the magnetic field $\vec B$, then we also have
$$\nabla \times (\vec B+\nabla \Phi) =\mu_0 \vec J$$
for any (smooth) scalar field $\Phi$. This means that there is not a unique solution to the problem since $\vec B +\nabla \Phi$ is also a solution for any (smooth) $\Phi$.
However, if we also specify the divergence of the magnetic field (we know that it is zero), then we can pursue a unique solution. For example,
$$\nabla \times \nabla \times \vec B =-\mu_0 \nabla \times \vec J$$
whereupon using the vector identity $\nabla \times \nabla \times \vec B= \nabla ( \nabla \cdot \vec B)-\nabla^2 \vec B$ and exploiting $\nabla \cdot \vec B=0$ gives
$$\nabla^2 \vec B=-\mu_0 \nabla \times \vec J$$
which has solution
$$\vec B(\vec r)=\mu_0 \int_V \frac{\nabla' \times \vec J(\vec r')}{4\pi |\vec r-\vec r'|}dV'$$
where the volume integral extends over all space where $\vec J=\ne 0$. We can integrate by parts in three dimensions by using the vector product rule identity $\nabla \times (\Phi \vec A) = \Phi \nabla \times \vec A+\nabla \Phi \times \vec A$ to write
$$\begin{align}
\vec B(\vec r)&=\mu_0 \int_V \frac{\nabla' \times \vec J(\vec r')}{4\pi |\vec r-\vec r'|}dV'\\\\
&=\mu_0 \int_V \left(\nabla' \times \left(\frac{\vec J(\vec r')}{4\pi |\vec r-\vec r'|}\right) -\nabla' \left(\frac{1}{4\pi |\vec r-\vec r'|}\right) \times \vec J(\vec r') \right)dV'\\\\
&=\mu_0 \oint_S \frac{\hat n' \times \vec J(\vec r')}{4\pi |\vec r-\vec r'|}dS'-
\mu_0 \int_V \nabla' \left(\frac{1}{4\pi |\vec r-\vec r'|}\right) \times \vec J(\vec r') dV'
\end{align}$$
Now, we may extend the integration region to all of space. Then, if $\vec J=0$ outside a finite region, then the surface integral vanishes and we have
$$\begin{align}
\vec B(\vec r)&= -\mu_0
\int_V \nabla' \left(\frac{1}{4\pi |\vec r-\vec r'|}\right) \times \vec J(\vec r') dV'\\\\
&=\mu_0 \int_V \nabla \left(\frac{1}{4\pi |\vec r-\vec r'|}\right) \times \vec J(\vec r') dV'\\\\
&=\nabla \times \left( \mu_0 \int_V \frac{\vec J(\vec r')}{4\pi |\vec r-\vec r'|}dV' \right)
\end{align}$$
where the first equality is effectively the Biot-Savart law and the last equality reveals that $\vec B =\nabla \times \vec A$ for the vector potential
$$\vec A(\vec r) = \mu_0 \int_V \frac{\vec J(\vec r')}{4\pi |\vec r-\vec r'|}dV'$$
Yes, there's a more elegant way! It uses the language of differential forms, which has replaced the 19th-century language of gradients, divergences, and curls in modern geometry. You can appreciate the simplicity of this language even before learning how to read it:
For any 1-form $A$,
$$\begin{align*}
(\star d)(\star d)A & = (\star d \star)dA \\
\operatorname{curl} \operatorname{curl} A & = d^\dagger dA
\end{align*}$$
Recalling that $\Delta = d d^\dagger + d^\dagger d$, we see that
$$\begin{align*}
\operatorname{curl} \operatorname{curl} A & = -d d^\dagger A + \Delta A \\
& = d(\star d \star)A + \Delta A \\
& = \operatorname{grad} \operatorname{div} A + \Delta A
\end{align*}$$
This is the identity you wanted to prove, where $-\Delta$ is the vector Laplacian.
My favorite place to learn about differential forms is in Chapters 4 and 5 of Gauge Fields, Knots, and Gravity by John Baez and Javier Muniain.
Here's a rough glossary that should help you move between the language of differential forms and the old language of vector calculus. I'll start by telling you the various kinds of differential forms, and the basic operations on them.
- In $n$-dimensional space, there are $n+1$ kinds of differential forms, from 0-forms to $n$-forms. You can think of a $k$-form as a $k$-dimensional density. A 0-form is a function, and a 1-form is a row-vector field (in coordinate-free language, a dual-vector field).
- The exterior derivative is an operation $d$ that turns $k$-forms into $(k + 1)$-forms. As its name suggests, it generalizes the operation of differentiating a function.
- The Hodge star is an operation $\star$ that turns $k$-forms in to $(n - k)$-forms. It comes from the dot product between column vectors. In fact, the Hodge star encodes the same geometric information as the dot product: if you know one, you can reconstruct the other.
- The codifferential is an operation $d^\dagger$ that turns $k$-forms into $(k - 1)$-forms. In an odd-dimensional space, like ordinary 3-dimensional space, applying $d^\dagger$ to a $k$-form is the same as applying $(-1)^k \star d \star$. In an even-dimensional space, $d^\dagger$ always acts like $-\star d \star$.
If you keep in mind that a 0-form is a function and a 1-form is a row-vector field, all the familiar operations of vector calculus can be written in terms of the ones above.
- The gradient of a function $f$ is the 1-form $df$.
- The curl of a 1-form $A$ is the 1-form $\star dA$.
- The divergence of a 1-form $A$ is the function $\star d \star A$.
- The Laplacian of a function or 1-form $\omega$ is $-\Delta \omega$, where $\Delta = dd^\dagger + d^\dagger d$. The operator $\Delta$ is often called the Laplace-Beltrami operator.
With this glossary in hand, you should be able to follow the steps of the calculation above, which is mostly just translating back and forth between languages. The only tricky bit is getting the sign right when you rewrite $d^\dagger$ as $\pm \star d \star$: you have to figure out what kind of form $d^\dagger$ is being applied to.
Best Answer
To make things easier, first change the $x,y,z$ coordinates to a new system of Cartesian coordinates $x',y',z'$ parallel to the original one but with the origin at $(0,0,-1)$, so $x'=x$, $y'=y$, $z'=z+1$. In the new coordinate system the vector $\vec g$ is:
$$\vec g = \frac{(x',y',z')}{(x^{\prime 2}+y^{\prime 2}+z^{\prime 2})^{3/2}} = \frac{\vec r}{r^3}\,,$$
where $\vec r=(x',y',z')$ is the radial vector field (from the new origin of coordinates) and $r$ its magnitude. This brings up the symmetry of the problem (and its physical connection; $\vec g$ is proportional to the electric field of a static electric charge at the new origin of coordinates, or to the the magnetic field produced by a Dirac's magnetic monopole at that point). Although you can work in Cartesian coordinates, to take advantage of the symmetry it is probably easier to work using standard spherical coordinates $r,\theta,\phi$, and that is what I am going to do, later converting back to Cartesian coordinates. Before we start it is a good idea to check that $\nabla\cdot\vec g=0$, otherwise the problem has no solution. But the divergence of this vector field is well known ($\vec g$ is essentially the gradient of the Green's function of the Poisson equation):
$$ \nabla\cdot\frac{\vec r}{r^3} = 4\pi\delta^3(\vec r)\,.$$
Although outside the origin the divergence is zero, because it is not zero everywhere in $R^3$, it is not possible to find a vector $\vec f$ that is non-singular everywhere $R^3$ minus the origin and such that $\nabla\times\vec f=\vec g$. However it is going to be possible to find $\vec f$ for the required region $z>0$. Note that the solution to this equation is far from unique, because adding to a solution any gradient also gives a solution. We need only to find one solution, all other solution are obtained form this one by adding an arbitrary gradient. To make the calculations easier we can impose some reasonable condition that will restrict the number of available solutions (physicsits call this to "fix the gauge"). For example, if you want to work in Cartesian coordinates you could require your solution to have $f_z=0$ (physicists call this choice the "axial gauge"). Because $\vec g$ is radial, and the radial component of $\vec f$ does not contribute to it, I am going to make the radial component of $\vec f$, $f_r$ equal to zero. With this choice, using the expression of the curl in spherical coordinates we have for $f_\theta$ and $f_\phi$:
$$ \begin{align} \frac{1}{r\sin{\theta}}\left[\frac{\partial}{\partial\theta}\left(\sin\theta f_\phi\right)-\frac{\partial f_\theta}{\partial\phi}\right] &= \frac{1}{r^2} \\ \frac{\partial}{\partial r}\left(rf_\phi\right) &= 0 \\ \frac{\partial}{\partial r}\left(rf_\theta\right) &= 0 \end{align} $$
The last two equations show that in this gauge ($f_r=0$), $f_\theta$ and $f_\phi$ have to be proportional to $1/r$, which makes also the radial part of the first equation correct. In fact, with this in mind, it is easily seen by inspection of the first equation that the following choice is a solution:
$$f_\phi = \frac{C-\cos\theta}{r\sin\theta},~~f_\theta = 0\,,$$
where $C$ is some constant. Because we want the solution to be valid if $z>0$ we must choose $C=1$. In other words, if $\hat u_\phi$ is the unit vector along the $\phi$ direction, then:
$$\vec f = \frac{1-\cos\theta}{r\sin\theta}\hat u_\phi$$
is a solution to the problem in spherical coordinates. The only thing that remains is to convert back to Cartesian coordinates. If $\hat u_x, \hat u_y, \hat u_z$ are the unit vectors along the Cartesian axes (old and new!), and using: $$\hat u_\phi = -\sin\phi\,\hat u_x+\cos\phi\,\hat u_y$$
and
$$\cos\phi = \frac{x'}{r\sin\theta},~~\sin\phi = \frac{y'}{r\sin\theta}, ~~\cos\theta = \frac{z'}{r},~~r^2\sin^2\theta = x'^2+y'^2$$
we get:
$$\vec f = \frac{y'(z'-r)}{r(x'^2+y'^2) }\hat u_x -\frac{x'(z'-r)}{r(x'^2+y'^2) }\hat u_y\,.$$
Reverting to the original $x,y,z$ coordinate system (just replace $x'~y'$ by $x~y$, and $z'$ by $z+1$) we get the final answer:
$$\vec f = \frac{y(z+1-r)\,\hat u_x-x(z+1-r)\,\hat u_y}{r(x^2+y^2)}$$
with
$$ r=\sqrt{x^2+y^2+(z+1)^2}\,.$$
Note that this solution satisfy the physicists "axial gauge" $f_z=0$, although I was not looking for it. Also note that this solution is regular for $z>-1$, as required, and becomes singular in the half line $x=y=0, z<=-1$.