In the context of a seminar lecture on some physico-chemical subject (magnetic shielding by electrons on atomic scale) I wanted to present a simple analytical example for the solution of an equation that is usually solved numerically. For that I have chosen a simple but still relevant case. Finally the thing boils down to the solution of this integral:
$$ \begin{align} \vec{\sigma}(\mathbf{r}_0) = & \int_{\mathbb{R}^3} \mathbf{r}\times (\vec{e}_z \times \mathbf{r})% \frac{\exp(-2|\mathbf{r}+\mathbf{r}_0|)}{r^3}\mathrm{d}\mathbf{r} \\ \tag{1}\end{align} $$
here $\vec{\sigma},\mathbf{r}, \vec{e}_z$ and $\mathbf{r}_0\in\mathbb{R}^3$ in particular the "position" vector $\mathbf{r}=\begin{pmatrix}x\\ y\\ z\end{pmatrix}$, $\mathbf{r}_0=\begin{pmatrix}x_0\\ y_0\\ z_0\end{pmatrix}$ and $\vec{e}_z=\begin{pmatrix} 0\\ 0\\ 1\end{pmatrix}$. $r:=\sqrt{\mathbf{r}\cdot\mathbf{r}}$ (i.e. the norm of $\mathbf{r}$). Accordingly for the resulting function we have $\vec{\sigma}: \mathbb{R}^3\to\mathbb{R}^3$. I use the integration symbol $\int_{\mathbb{R}^3} \mathrm{d}\mathbf{r}$ for the volume integral $\int_{-\infty}^\infty\int_{-\infty}^\infty\int_{-\infty}^\infty \mathrm{d}x\mathrm{d}y\mathrm{d}z$, which is sometimes also denoted $\int_{-\infty}^\infty\int_{-\infty}^\infty\int_{-\infty}^\infty \mathrm{d}V$
I have tried various strategies, namely all 4 combinations of A) with and without transformation of the variable ($\mathbf{r}':=\mathbf{r}+\mathbf{r}_0$) and B) using cartesian or spherical coordinates ($\int_{-\infty}^\infty\int_{-\infty}^\infty\int_{-\infty}^\infty \cdot \mathrm{d}x\mathrm{d}y\mathrm{d}z$=$\int_{0}^\pi\int_{0}^{2\pi}\int_{0}^\infty \cdot r^2\sin\theta\mathrm{d}r\mathrm{d}\phi\mathrm{d}r$), but in each case I obatin finally expressions which I fail to solve. For example in cartesian coordinates and without transformation I finally obtain the expression $$ \int_{\mathbb{R}^3} \begin{pmatrix}-xz\\-yz\\x^2+y^2\end{pmatrix} \frac{\exp(-2\sqrt{x^2 +y^2 +z^2 + 2(xx_0 + yy_0 + zz_0) + r_0^2})}{(x^2+y^2+z^2)^\frac{3}{2}} \tag{2}$$
I believe that a anaylitcal (symbolic) integration shall be feasible, since by some physical arguments one can "guess" a solution that fulfills certain constraints one would expect the solution to fulfill. This solutions is
$$ \vec{\sigma}^{guess}(\mathbf{r}_0)=\begin{pmatrix}0\\0\\1\end{pmatrix}\frac{1+2r}{\pi}e^{-2r_0}$$
However the knowledge of this didin't really help me to solve $(1)$ except to find some symmetry arguments on why the $x$ and $y$ components are $=0$ (For example this can be seen from $(2)$ with 'ungerade' vector components in $x,y$ and 'gerade' ones $z$).
Who can help me to solve $(1)$ "systematically"?
(Note: This question is an extension/generalzation of this one I have asked previously.)
Not a solution, but too long for a comment
This is a more complicated task then the previous one, because you have two constant vectors $\big(\vec r_0$ and $\vec a=\vec e_z\,\big)$. First of all, you decide how to choose the system of coordinates for $\vec r$ (the origin and direction). Making the change $\vec\rho=\vec r_0+\vec r_0$ and simplifying the double vector multiplication, you can present the integral in the form $$\vec\sigma(\vec r_0, \vec a)=\int_{R^3}\bigg((\vec\rho-\vec r_0)\;\vec a\cdot(\vec\rho-\vec r_0)-\vec a(\vec\rho-\vec r_0)^2\bigg)\frac{e^{-2\rho}}{|\vec\rho-\vec r_0|^3}d^3\rho$$ The natural choice is to direct the axis $\vec e_z$ along the vector $\vec r_0$, and the integral takes the form $$\vec\sigma(\vec r_0, \vec a)=\int_{R^3}\bigg((\vec\rho-\vec r_0)\;\vec a\cdot(\vec\rho-\vec r_0\big)-\vec a(\vec\rho-\vec r_0)^2\bigg)\frac{e^{-2\rho}}{(\rho^2-2r_0\rho\cos\theta+r_0^2)^{3/2}}\sin\theta d\theta d \phi \rho^2 d\rho$$ Generally speaking, this integral as a vector can be expressed via two constant vectors we have. So, the most general form is $$\vec\sigma(\vec r_0, \vec a)=\vec a f_1(\vec r_0, \vec a)+\vec r_0 f_2(\vec r_0, \vec a)+[\vec r_0\times\vec a]f_3(\vec r_0, \vec a)$$ where all functions $f_i(\vec r_0, \vec a)$ are scalars, depending on $|\vec a|^2,|\vec r_0|^2, (\vec a,\vec r_0) $.
Multiplying consequentially $\vec\sigma(\vec r_0, \vec a)$ by $\vec a, \vec r_0,[\vec r_0\times\vec a]$, we get a system of integrals (scalars), from which we get $f_i(\vec r_0, \vec a)$.
For example, $$\vec a\cdot\vec\sigma(\vec r_0, \vec a)=\int(\vec a\cdot(\vec\rho-\vec r_0))^2-(\vec a)^2(\vec\rho-\vec r_0)^2\bigg)\frac{e^{-2\rho}}{(\rho^2-2r_0\rho\cos\theta+r_0^2)^{3/2}}\sin\theta d\theta d \phi \rho^2 d\rho$$ where $\vec\rho\cdot\vec r_0=\rho r_0\cos\theta$ and $\vec a\cdot\vec \rho=a\rho\big(\sin\theta_a\cos\phi_a\sin\theta\cos\phi+\sin\theta_a\sin\phi_a\sin\theta\sin\phi+\cos\theta_a \cos\theta\big)$ (the vector $\vec a$ has the polar coordinates $a, \theta_a, \phi_a$ in the chosen system).
I'm not sure whether we can get a closed form in general case, but in some specific cases (for example, for $[\vec r_0\times\vec a]=0$) it can be done.
$\mathbf{Addendum}$ The case $\vec a=\lambda\vec r_0$. In this case we have only one constant vector, and $$\vec\sigma(r_0)=\lambda\int_{R^3}\Big(\vec r_0 r^2-\vec r(\vec r,\vec r_0)+\vec r r_0^2-\vec r_0(\vec r,\vec r_0)\Big)\frac{e^{-2r}}{|\vec r-\vec r_0|^3}d^3r=2\pi\lambda\vec r_0\frac{f(r_0)}{r^2_0}\quad(1)$$
Multiplying both sides by $\vec r_0$, integrating with respect to the angle $\phi$, and denoting $\cos\theta =x$ $$f(r_0)=r_0^2\int_0^\infty e^{-2r}r^4dr\int_{-1}^1\frac{1-x^2}{(r^2-2rr_0x+r^2_0)^{3/2}}dx\quad(2)$$ First of all, to check our answer afterwards, we note that at $r_0\to 0$ $$f(r_0)\to r_0^2\int_0^\infty e^{-2r}rdr\int_{-1}^1(1-x^2)dx=\frac{r_0^2}{3}\qquad(3)$$ To evaluate the integral we can use the decomposition into the Legendre' polynomials, or simply to integrate by part several times. Using $$\frac{1}{rr_0}\frac{d}{dx}\frac{1}{\sqrt{r^2-2r_0rx+r_0^2}}=\frac{1}{(r^2-2r_0rx+r_0^2)^{3/2}}$$ we integrate by part with respect to $x$ several times: $$f(r_0)=r^2_0\int_0^\infty e^{-2r}r^4dr\int_{-1}^1\frac{1-x^2}{r_0r}d\Big(\frac{1}{\sqrt{r^2-2r_0rx+r_0^2}}\Big)$$ $$=-2r_0\int_0^\infty e^{-2r}r^3dr\int_{-1}^1\frac{x}{r_0r}d\big(\sqrt{r^2-2r_0rx+r_0^2}\big)$$ $$=-2\int_0^\infty e^{-2r}r^2dr\Big(\sqrt{(r+r_0)^2}+\sqrt{(r-r_0)^2}\Big)-\frac{2}{3r_0}\int_0^\infty e^{-2r}rdr\int_{-1}^1d\Big((r^2-2r_0rx+r_0^2)^{3/2}\Big)$$ $$=-2\int_0^\infty e^{-2r}r^2(r+r_0)dr-2\int_0^{r_0} e^{-2r}r^2(r_0-r)dr - 2\int_{r_0}^\infty e^{-2r}r^2(r-r_0)dr$$ $$+\frac{2}{3r_0}\int_0^\infty e^{-2r}r(r+r_0)^3dr-\frac{2}{3r_0}\int_0^{r_0} e^{-2r}r(r_0-r)^3dr-\frac{2}{3r_0}\int_{r_0}^\infty e^{-2r}r(r-r_0)^3dr$$ $$=-4r_0\int_0^\infty e^{-2r}r^2dr-4\int_0^{r_0} e^{-2r}r^2(r_0-r)dr$$ $$+\frac{4}{3}\int_0^\infty e^{-2r}r(r^2_0+3r^2)dr-\frac{4}{3r_0}\int_0^{r_0} e^{-2r}r(r_0-r)^3dr$$ Performing integration, we get: $$f(r_0)=\frac{1}{r_0}\Big(1-e^{-2r_0}(1+2r_0+2r_0^2+r_0^3)\Big)\qquad(4)$$ At $r_0\to 0$, decomposing the exponent into the Taylor's series $$f(r_0)=\frac{1}{r_0}\Big(1-\big(1-2r_0+2r_0^2-\frac{4}{3}r_0^3+O(r_0^4)\big)\big(1+2r_0+2r_0^2+r_0^3\big)\Big)$$ $$f(r_0)=\frac{r_0^2}{3}+O(r_0^3)\,\,\text{at}\,\,r_0\to0$$ in accordance with our expectations (see (3)).
Finally, for this specific case we get the answer: $$\boxed{\,\,\vec\sigma(r_0)=2\pi\lambda\,\vec r_0\frac{f(r_0)}{r^2_0}=2\pi\lambda\,\frac{\vec r_0}{r_0^3}\Big(1-e^{-2r_0}(1+2r_0+2r_0^2+r_0^3)\Big)\,\,}$$