As is well-known, the vector field $F = (-y,x)/(x^2+y^2)$ on the punctured plane is irrotational yet not conservative, i.e., there is no scalar potential $f$ such that $F = \nabla f$, because for instance $\oint_{S^1} F = 2\pi$. From a more abstract point of view, vector fields on the punctured plane correspond to 1-forms, irrotational fields correspond to closed forms, conservative fields correspond to exact forms, and $H^1_{dR}(\mathbb R^2 \setminus 0) \cong H^1_{dR}(S^1)$ is 1-dimensional and spanned by $d\theta$, so that any irrotational and non-conservative field on the punctured plane is a scalar multiple of $F$ plus a conservative field.
I am looking for the 3D analogue of the above vector field. Since $H^2_{dR}(\mathbb R^3 \setminus 0) \cong H^2_{dR}(S^2)$ is 1-dimensional, there must exist an incompressible vector field $F$ (meaning that $\nabla \cdot F = 0$, or that the corresponding 2-form is closed) on $\mathbb R^3 \setminus 0$ that does not admit a vector potential (meaning that there is no field $G$ such that $F = \nabla \times G$, or that the corresponding 2-form is not exact). How could one come up with a formula for such a field? I tried to understand how one could get the formula for $F$ in the 2D case from the 1-form $d\theta$ on the circle by pulling it back through the deformation retract $\mathbb R^2 \setminus 0 \to S^1$ given by $(x,y) \mapsto (x,y)/\sqrt{x^2 + y^2}$, but I was not able to fill in the details.
As suggested in the comments, the "point mass field" (better known to PDE people as the gradient of the fundamental solution of Laplace's equation) $$\def\x{\mathbf{x}}F(\x)=|\x|^{-3}\x$$ is such an example. A routine calculation shows that $\nabla \cdot F = 0.$ In order to show $F$ is not a curl, we can use the same idea as in your 2-dimensional example: test its integral over a closed chain. In this case the unit sphere $S^2 \subset \mathbb R^3 \setminus \{0\}$ is the obvious choice, since $F$ coincides with the outwards normal and thus $$\int_{S^2}F\cdot\nu\;dA=4 \pi.$$
Since $S^2$ is closed (i.e. has no boundary), Stokes' theorem tells us that $$\int_{S^2} (\nabla \times G) \cdot \nu \, dA = 0$$ for any $G;$ so we conclude that $F$ is not a curl.
As you seem to have guessed, this is really a very direct generalization of the two-dimensional version: in terms of differential forms, they're both constructed from the volume form of the sphere, which is a natural choice because we can view $\mathbb R^{n+1}\setminus \{0\}$ as the cylinder $S^n \times \mathbb R.$
In both cases (and in higher dimensions too), $F$ is a vector field representation of the differential $n$-form $\omega = \phi^* dA$ where $\phi : \mathbb R^{n+1} \setminus\{0\}\to S^n$ is the radial projection map and $dA$ is the $n$-volume form of $S^n$. By pullback invariance we can conclude $d\omega = \phi^*d(dA) = 0,$ but integration over the sphere gives $\int_{S^n} \omega = \int_{S^n} dA > 0.$
The correspondence between the vector field and form perspectives comes from the Hodge dual and the musical isomorphism. In spherical coordinates $(r,\theta,\varphi)$ on $\mathbb R^3\setminus \{0\}$ and $(\theta,\varphi)$ on $S^2,$ the projection map is simply $\phi(r,\theta,\varphi)=(\theta,\varphi).$ Thus $$\omega = \phi^*(\sin \theta\, d \varphi\wedge d \theta) = \sin \theta \, d \varphi \wedge d \theta.$$ If we let $\hat\cdot$ denote metric normalization, then since $|d\varphi|=1/(r \sin \theta), |d \theta| = 1/r$ we have $$\omega = \frac{1}{r^2} \widehat{d \varphi} \wedge \widehat{d \theta}.$$ Thus $$\star \omega = \frac{1}{r^2}\widehat{dr} = \frac1{r^2}dr,$$ which has corresponding vector field $$(\star \omega)^\sharp = \frac1{r^2}\partial_r = \frac\x{|\x|^3}.$$