The vector field
$$
G(x,y,z)=\frac{1}{(x^2+y^2+z^2)^{3/2}}\begin{pmatrix}x\\y\\z\end{pmatrix}\,,\quad\text{ on }\mathbb R^3\setminus\{0\},
$$
describes Coulomb's law and
Newton's law of universal gravitation and is mathematically interesting in many ways: It is rotation free and divergence free. A scalar function
whose gradient is $G$ is
$$
H(x,y,z)=-\frac{1}{(x^2+y^2+z^2)^{1/2}}\,.
$$
The divergence freeness of $G$ suggests that there should be a vector field $F$ whose curl is $G\,.$ It turns out that
- $F$ exists but cannot be defined on all of $\mathbb R^3\setminus\{0\}\,.$
This was briefly remarked in this answer. After some confusion on my side an finally learning it from Ted Shifrin I'd like to sum it up a bit: If $F$ is defined on all of $\mathbb R^3\setminus\{0\}$ one can apply Stokes' theorem to each hemisphere separately to conclude that the flux of $G$ through the whole unit sphere must be zero but this cannot happen because the flux of $G$ is
\begin{align}
\oint_{S^2}\hat{\mathbf{n}}\cdot G\,dS&=\oint_{S^2}\frac{\left(\begin{smallmatrix}x\\y\\z\end{smallmatrix}\right)}{(x^2+y^2+z^2)^{1/2}}\cdot \frac{\left(\begin{smallmatrix}x\\y\\z\end{smallmatrix}\right)}{(x^2+y^2+z^2)^{3/2}}\,dS\\
&=\int_0^{2\pi}\int_0^\pi \sin\theta\,d\theta\,d\varphi=4\pi\,.
\end{align}
To see that on the contrary Stokes' theorem leads to a flux of zero (assuming that $F$ is defined on all of $\mathbb R^3\setminus\{0\}$) note that the line integral of $F$ along the boundary of each hemisphere equals the flux of $\nabla\times F=G$ through that hemisphere. But these line integrals cancel out because of their opposite orientation which leads to zero flux. A contradiction.
- Consequently, when $F$ whose curl is $G$ exists there must be at least one point on one of the two hemispheres on which $F$ is not defined. An example for such an $F$ is
$$
F(x,y,z)=\frac1{r^2+rz}\begin{pmatrix}-y\\x\\0\end{pmatrix}\,,\quad r:=(x^2+y^2+z^2)^{1/2}\,,
$$
which is a simpler notation for the solution found here. This $F$ is defined for all $x,y,z\in\mathbb R^3$ except on the ray $z\le 0,x=y=0$ (where $r+z$ becomes zero) and
describes a clockwise circular motion around the $z$-axis and parallel to the $xy$-plane. A way to visualise it is to look at the speed of this motion
which is
$$
\|F\|=\frac{\varrho}{\varrho^2+z^2+z\sqrt{\varrho^2+z^2}}\,,\quad\varrho:=\sqrt{x^2+y^2}\,.
$$
The following picture shows this speed as a function of the distance $\varrho$ to the $z$-axis for various levels of $z\,.$ For $z\le 0$ where $F$ is not defined on the $z$-axis the speed blows up when the $z$-axis is approached.

Another interesting picture shows the $[0,3\pi/2]$-sectors of the surfaces of constant speed $\|F\|\,.$
The blue surfaces has the slowest speed and the brown surface the highest.
These surfaces converge at the origin.

Further properties of $F\,:$
The one-form
$$
\omega=\frac{-y\,dx+x\,dy}{r^2+rz}
$$
is closed if and only if $z=0\,.$ This follows from $(\nabla\times F)_3=\partial_xF_2-\partial_yF_1=G_3=z/r^3\,.$ For $z=0$ this form has interesting properties that I summarized here.
On the unit sphere $S^2$ the speed of the integral curves can be written
more simply as
$$
\|F\|=\sqrt{\frac{1-z}{1+z}}\,.
$$
This is zero at the north pole $z=1$ and infinite at the south pole $z=-1\,.$ The vorticity of this vector field
is, as we know, $\nabla\times F=G$ which on $S^2$ is simply $\left(\begin{smallmatrix}x\\y\\z\end{smallmatrix}\right).$
Due to the circular motion in planes parallel to the $xy$--axis this seems odd at points away from the poles.
But this rotationally symmetric vorticity $G$ arises from different speeds of integral curves next to each other.
In particular, near the equator we have practically a
parallel flow with shear that has a nice animation here.