Spherical integration of vector function $\frac{\vec{r}\times(\vec{\mathbb{e}_z}\times\vec{r})}{r^3} \exp^{-2 r}\mathrm{d}{\vec{r}}$ In the context of a problem about the magnetic response of the electron cloud in the hydrogen atom I need the solution of the following volume integral of a vector function:
\begin{align} I &= \frac{1}{4\pi}\int_{\mathbb{R}^3} \frac{\vec{r}\times(\vec{\mathbb{e}_z}\times\vec{r})}{r^3} \exp^{-2 r}\mathrm{d}{\vec{r}} \\ &= \frac{1}{4\pi}\int_{\mathbb{R}^3} \begin{pmatrix}-xz \\ -yz \\ x^2+y^2 \end{pmatrix} \frac{\exp^{-2 r}}{r^3}\mathrm{d}{x}\mathrm{d}{y}\mathrm{d}{z}\end{align}
with $\vec{r}=\begin{pmatrix}x \\ y \\ z \end{pmatrix},\;\vec{\mathbb{e}}_z=\begin{pmatrix}0 \\ 0 \\ 1 \end{pmatrix},\; r=|\vec{r}|=\sqrt{x^2+y^2+z^2}.$
I have tried solving it with Mathematica, but it seems to evaluate for ever. Can anyone help me?

By symmetry the $x$ and the $y$ component is $= 0$. For the $z$ component we obtain
\begin{align}
 I_z & = \frac{1}{4\pi} \int_0^\infty\int_0^\pi \int_0^{2\pi} (r^2-z^2) \frac{\exp^{-2r}}{r}   \sin{(\vartheta)}    dr d\vartheta  d\varphi \\
& = \frac{1}{4\pi} \int_0^\infty r \exp^{-2r} dr\int_0^\pi \sin^3(\vartheta)d\vartheta \int_0^{2\pi} d\varphi \\
& = \frac{1}{4\pi}\frac{1}{4} \frac{4}{3} 2\pi \\
& = \frac{1}{6}
\end{align}
Hence \begin{align}
 I & = \begin{pmatrix} 0\\ 0\\ \frac{1}{6}\end{pmatrix} 
\end{align}
