Are there cases where $\nabla\cdot\iiint\frac{\mathbf{J}(\mathbf{x}')}{\left|\mathbf{x}-\mathbf{x'}\right|}\mathrm{d}V' \neq 0$? In Jackson's Classical Electrodynamics, Section 5.4 (Vector Potential), the author seems to assume that because $\nabla\cdot\mathbf{J} = 0$, the following holds for the current density (where the integral is done over all space):
$$\nabla\cdot\iiint\frac{\mathbf{J}(\mathbf{x}')}{\left|\mathbf{x}-\mathbf{x'}\right|}\mathrm{d}V' = 0$$
However, in general we know that in a closed volume $V$, we have:
$$\nabla\cdot\iiint\limits_V\frac{\mathbf{J}(\mathbf{x}')}{\left|\mathbf{x}-\mathbf{x'}\right|}\mathrm{d}V' = \iiint\limits_V\nabla\cdot\frac{\mathbf{J}(\mathbf{x}')}{\left|\mathbf{x}-\mathbf{x'}\right|}\mathrm{d}V' = \iiint\limits_V\mathbf{J}(\mathbf{x}')\cdot\nabla\left(\frac{1}{\left|\mathbf{x}-\mathbf{x'}\right|}\right)\mathrm{d}V'$$
$$ = -\iiint\limits_V\mathbf{J}(\mathbf{x}')\cdot\nabla'\left(\frac{1}{\left|\mathbf{x}-\mathbf{x'}\right|}\right)\mathrm{d}V' = -\iiint\limits_V\nabla'\cdot\frac{\mathbf{J}(\mathbf{x}')}{\left|\mathbf{x}-\mathbf{x'}\right|}\mathrm{d}V' + \iiint\limits_V\frac{\nabla'\cdot\mathbf{J}(\mathbf{x}')}{\left|\mathbf{x}-\mathbf{x'}\right|}\mathrm{d}V'$$
Now in the magnetostatic case, we know from the continuity equation $\nabla\cdot\mathbf{J}+\frac{\partial\rho}{\partial t} = 0$ that in fact, $\nabla\cdot\mathbf{J} = 0$ over all space (because there are no local charge density fluctuations), and so the second term goes away and we are left with:
$$\nabla\cdot\iiint\limits_V\frac{\mathbf{J}(\mathbf{x}')}{\left|\mathbf{x}-\mathbf{x'}\right|}\mathrm{d}V' = -\iiint\limits_V\nabla'\cdot\frac{\mathbf{J}(\mathbf{x}')}{\left|\mathbf{x}-\mathbf{x'}\right|}\mathrm{d}V' = -\mathop{\LARGE\unicode{x222f}}\limits_{\partial V}\frac{\mathbf{J}(\mathbf{x}')}{\left|\mathbf{x}-\mathbf{x'}\right|}\cdot\mathrm{d}\mathbf{S}'$$
Now taking a limit as the volume becomes the entire space, we see that the surface integral is taken in regions increasingly far away from $\mathbf{x}$. If we make certain assumptions on the asymptotics of $\mathbf{J}$, such as $|\mathbf{J}(\mathbf{x}')| = o\left(\frac{1}{\left|\mathbf{x}-\mathbf{x}'\right|}\right)$, then we can bound the surface integral and show that it tends to zero.
Physically speaking, one could argue that this is enough because realistic systems will always have currents bound in a finite volume anyway. But sometimes we consider idealized scenarios such as infinite straight lines with a uniform current $I$. These cases might cause trouble. If you only have a single (or finitely many) such wires, I think you can still show the integral goes to zero because of the inverse dependence to the distance in the integrand (I am not certain). But even then, one could reasonably imagine idealized situations comprising infinitely many such wires that could make the surface integral not converge to zero.
Another difficulty is that the convergence ought to hold for any surface that encloses all of space eventually. If we restrict to spheres of radius $|\mathbf{x}-\mathbf{x}'| = R$ centered around $\mathbf{x}$, then the convergence is trivial because we get $$-\mathop{\LARGE\unicode{x222f}}\limits_{\partial V}\frac{\mathbf{J}(\mathbf{x}')}{\left|\mathbf{x}-\mathbf{x'}\right|}\cdot\mathrm{d}\mathbf{S}' = -\frac{1}{R}\mathop{\LARGE\unicode{x222f}}\limits_{\partial V}\mathbf{J}(\mathbf{x}')\cdot\mathrm{d}\mathbf{S}' = 0$$
(by the divergence theorem and continuity equation). But if your surfaces are non-spherical, this trick does not really work anymore. But maybe there is a way to avoid this issue.

I am interested in knowing the most general assumptions that can be made on $\mathbf{J}$ to satisfy the convergence, and also in knowing about situations that might be considered where this assumption is wrong.
 A: Why do you claim that $$ \int_V \nabla \cdot \frac{{\bf J}({\bf x}')}{|{\bf x}-{\bf x}'|} {\rm d}V' = -\int_V \nabla' \cdot \frac{{\bf J}({\bf x}')}{|{\bf x}-{\bf x}'|} {\rm d}V'$$? This doesn't seem to be true. You should have \begin{align} \int_V \nabla \cdot \frac{{\bf J}({\bf x}')}{|{\bf x}-{\bf x}'|} {\rm d}V' &= \int_V \Big(\nabla \frac{1}{|{\bf x}-{\bf x}'|}\Big) \cdot {\bf J}({\bf x}') {\rm d}V' =\\&= \int_V \Big(-\nabla' \frac{1}{|{\bf x}-{\bf x}'|}\Big) \cdot {\bf J}({\bf x}') {\rm d}V' = \\
&= \int_V \frac{1}{|{\bf x}-{\bf x}'|} \nabla' \cdot {\bf J}({\bf x}') {\rm d}V' - \int_V \nabla' \cdot \Big(\frac{1}{|{\bf x}-{\bf x}'|} {\bf J}({\bf x}') \Big){\rm d}V' =\\&= 0 - \int_{\partial V} \frac{1}{|{\bf x}-{\bf x}'|} {\bf J}({\bf x}') \cdot {\rm d}{\bf S}\end{align}
and this vanishes if ${\bf J}$ vanishes fast enough. Note that if it doesn't vanish fast enough, then you may have problem with the convergence of integral $\iiint \frac{{\bf J}({\bf x}')}{|{\bf x}-{\bf x}'|} {\rm d}V'$ taken over the whole space; to make sure that it is well-defined, you need ${\bf J}$ to vanish fast enough, and that makes the boundary term to vanish at infinity.
A: Your idea on spheres of radius $R$ gives the hint; if the current density behaves well, the large distance of the surface from $\mathbf x$ will make the integral, in the limit $R\to\infty$, go to zero.
The situation considered is electric current flows without accumulation of charges; what goes into a region, also goes out of the region. This is very common in practice, total current going in equals total current going out. In practice, this current going in is limited, no matter how the boundary surface is chosen, because the current in single wire is finite and the number of wires is finite.
The integral in question can be divided into two parts, one due to current going in to the region, and one due to current going out of the region:
$$
\int_{\partial V} \mathbf j/r\cdot d\mathbf S = C_{in} + C_{out}
$$
where
$$
C_{in} = \int_{\partial V} \chi_{in}(\mathbf x')\mathbf j(\mathbf x')/r\cdot d^2\mathbf x'
$$
and
$$
C_{out} = \int_{\partial V} \chi_{out}(\mathbf x')\mathbf j(\mathbf x')/r\cdot d^2\mathbf x'
$$
where $\chi_{in}$ is characteristic function of the part of the surface where the current goes in ($\mathbf j\cdot d\mathbf S < 0$).
Let us use the triangle inequality:
$$
\left| \int_{\partial V} \mathbf j/r\cdot d\mathbf S \right| \leq |C_{in}| +| C_{out}|.
$$
So the integral will obviously go to zero if both $C_{in}$ and $C_{out}$ go to zero. These two going to zero is a sufficient condition (it may not be a necessary one).
Let the smallest $r=|\mathbf x - \mathbf x'|$ for some stage of the limiting process be denoted $r_{min}$. Of course, in the limiting process, the whole boundary has to expand to infinity, so $r_{min}\to\infty$.
$$
|C_{in}| \leq \int_{\partial V} \chi_{in}(\mathbf x')|\mathbf j(\mathbf x')\cdot d^2\mathbf x'|/r_{min} = \frac{I_{in}}{r_{min}}
$$
where $I_{in}$ is (positive) value of current due to charges that come into the region. As long as this current does not grow too quickly with $r_{min}$, the contribution will go to zero as the boundary is expanded to infinity. Similarly for the other contribution.
So the sufficient condition for the surface integral to go to zero is that the electric current going in through the surface is not growing too fast as the surface is expanded. If the current is limited by a known maximum value no matter the boundary, as is the case if the system is made of finite number of (possibly infinitely long) wires of finite current, then the integral will go to zero. Thus one can consider arbitrary finite number of infinite wires, each carrying finite current. However, if the number of wires crossing the boundary increases as quick or more quickly than $r_{min}$, then there could be a problem and the integral may not have limit 0. This does not look like a common situation though.
A: I have few thoughts rather than an answer. Maybe it will help. Since your current density is diverge-less maybe it is a good idea to apply Helmholtz decomposition?
For a generic, well-behaved 3d vector field with $\boldsymbol{\nabla}.\mathbf{J}=0$ we have:
$\mathbf{J}\left(\mathbf{r}\right)=\boldsymbol{\nabla}\times \mathbf{M}$
For some vector field $\mathbf{M}$, but then (using Levi-Civita symbol $\epsilon_{\alpha\beta\gamma}$ to deal with curl):
$\int d^3 r' \mathbf{J}\left(\mathbf{r'}\right).\boldsymbol{\nabla}'\frac{1}{\left|\mathbf{r}-\mathbf{r}'\right|}=\epsilon_{\alpha\beta\gamma}\int d^3 r' \partial'_{\beta}M_{\gamma}\left(\mathbf{r'}\right).\partial'_{\alpha}\frac{1}{\left|\mathbf{r}-\mathbf{r}'\right|}=\epsilon_{\alpha\beta\gamma}\oint d^2 r' \hat{n}_{\beta} M_{\gamma}\left(\mathbf{r'}\right).\partial'_{\alpha}\frac{1}{\left|\mathbf{r}-\mathbf{r}'\right|}$
Where $\boldsymbol{\hat{n}}$ is the normal to the surface of the volume over which you are integrating, and I used the symmetry of the Levi-Civita to get rid of the other integral.
Did this help? Well we still have a surface integral, but now we are talking about magnetization ($\mathbf{M}$) which does not have to vanish that rapidly for the integral to converge to zero. Also, you can repeat the same trick, indeed no generality is lost by assuming that $\boldsymbol{\nabla}.\mathbf{M}=0$. I wonder whether you can build an iterative proof where you go deeper and deeper into the derivatives of the current density which then establishes that the initial integral can be made as close to zero as you want.
