As noted by Christian Remling, the Fourier transform of $|x|^{-1}$ is $|\xi|^{-2}$, so the relevant integral you wish to bound is $$\iiint |\eta - \xi|^{-2} \bar{\hat{u}}(\xi) |\eta|^{2\alpha} \hat{u}(\zeta) |\eta - \zeta|^{-2} ~d\eta ~d\xi ~d\zeta $$ Do the $\eta$ integral first. This converges as long as $2\alpha < 1$. (At issue is the decay at infinity, if $2\alpha \geq 1$ the integrand is bounded below by $|\eta|^{-3}$ for large $|\eta|$; this issue with decay at infinity is directly tied to the lack of regularity of $|x|^{-1}$ at the origin.) Note that by the symmetries of $\mathbb{R}^3$ the resulting quantity is only a function of $|\xi - \zeta|$. So a scaling argument shows that $$ \int \frac{|\eta|^{2\alpha}}{|\eta - \xi|^2 |\eta - \zeta|^2} ~d\eta = C |\xi - \zeta|^{2\alpha - 1}$$ So your integral boils down to $$ \iint \bar{\hat{u}}(\xi) \hat{u}(\zeta) |\xi - \zeta|^{2\alpha - 1} ~d\xi ~d\zeta \lesssim \|u\|_{H^2}^2 \left(\iint \frac{|\xi - \zeta|^{4\alpha - 2}}{(1 + |\xi|^2)^2(1 + |\zeta|^2)^2} ~ d\xi~d\zeta \right)^{\frac12} $$ where the estimate is by Cauchy-Schwarz (once in $\xi$ and once in $\zeta$), and the initial $H^2$ bound on $u$ is assumed.

Since the set $|\xi - \zeta| = 0$ is three dimensional, and we know that $4\alpha - 2 \geq -2 > -3$ when $\alpha \in [0,\frac12)$, the integral in the parentheses converges.

The above sketch shows that $r^{-1} u$ is in $H^\alpha$ for $\alpha \in [0,\frac12)$. My comment shows that this range is sharp.

As noted by Christian Remling, the Fourier transform of $|x|^{-1}$ is $|\xi|^{-2}$, so the relevant integral you wish to bound is $$\iiint |\eta - \xi|^{-2} \bar{\hat{u}}(\xi) |\eta|^{2\alpha} \hat{u}(\zeta) |\eta - \zeta|^{-2} ~d\eta ~d\xi ~d\zeta $$ Do the $\eta$ integral first. This converges as long as $2\alpha < 1$. (At issue is the decay at infinity, if $2\alpha \geq 1$ the integrand is bounded below by $|\eta|^{-3}$ for large $|\eta|$; this issue with decay at infinity is directly tied to the lack of regularity of $|x|^{-1}$ at the origin.)

Writing $\mu = \frac12( \xi + \zeta)$ and $\rho \epsilon = \frac12 ( \xi - \zeta)$, where $\epsilon$ is a unit vector and $\rho$ a scalar, we see, with the change of variables $\eta = \mu + \rho \psi$ that we are down to estimating $$ \int \frac{|\eta|^{2\alpha}}{|\eta - \xi|^2 |\eta - \zeta|^2} ~d\eta = \int \frac{|\mu + \rho \psi|^{2\alpha}}{\rho |\psi -\epsilon|^2|\psi + \epsilon|^2} ~d\psi \leq C_1 \frac{|\mu|^{2\alpha}}{\rho} + C_2 \rho^{2\alpha - 1}. $$
Inserting it into the original integral (and taking moduli of the $u$ terms) we need to bound $$ \iint |\hat{u}(\xi)| |\hat{u}(\zeta)| \Big(\frac{|\xi + \zeta|^{2\alpha}}{|\xi - \zeta|} + |\xi - \zeta|^{2\alpha - 1} \Big) ~d\xi ~d\zeta \lesssim \\ \|u\|_{H^2}^2 \left(\iint \frac{\big(|\xi + \zeta|^{4\alpha} + |\xi - \zeta|^{4\alpha}\big) |\xi - \zeta|^{- 2}}{(1 + |\xi|^2)^2(1 + |\zeta|^2)^2} ~ d\xi~d\zeta \right)^{\frac12} $$ where the estimate is by Cauchy-Schwarz (once in $\xi$ and once in $\zeta$), and the initial $H^2$ bound on $u$ is assumed.

The integrability of the term inside the parentheses is not too hard to check. Split the integral into dyadic rectangles with $|\xi| \approx 2^k$ and $|\zeta|\approx 2^\ell$. When both $k,\ell \leq 0$ we are on a compact region, and the set $|\zeta - \xi| = 0$ has codimension 3 so the integral converges. When $k > \ell$ we have that $|\xi - \zeta| \approx |\xi + \zeta| \approx |\xi|$ and so the integral is bounded by $$ 2^{(4\alpha - 3)k} / (2^{-3\ell} + 2^{\ell}) $$ . When $k = \ell > 0$ we can replace the denominator by $|\xi|^{4}|\zeta|^4$. On the compact region the integral converges again due to $|\xi-\zeta| = 0$ having codimension 3, and the integrand scales like $2^{(4\alpha - 2 - 8)k}$. So the integral is bounded by $2^{(4\alpha - 4)k}$. The double summation in $k$ and $\ell$ therefore converges obviously.

The above sketch shows that $r^{-1} u$ is in $H^\alpha$ for $\alpha \in [0,\frac12)$. My comment shows that this range is sharp.