This is not an unswer, just a hint and too long for a comment: The integral above may be transformed with $\theta = \frac{\arccos(y)}{2}$
$$2\, \left(1+2\, r^2\right)^{1/3} \int_{-1}^1 \left(1-r^2 y\right)^{-1/3} \left(1-y^2\right)^{-\frac{1}{2}} \, dy$$
Next we can split the integral in two parts with $p=-y$ in the first one:
$$2\,\left(1+2\, r^2\right)^{1/3} \left(\int_0^1 \left(1+r^2 p\right)^{-1/3} \left(1-p^2\right)^{-\frac{1}{2}} \, dp + \int_0^1 \left(1-r^2 y\right)^{-1/3} \left(1-y^2\right)^{-\frac{1}{2}} \, dp \right) $$
and evaluate them in form of hypergeometric functions. For this step e.g. we can express the integrand in form of MeijerG - Function. The first integral will be transformed with $q = p^2$ apply the formular Mathematical Functions Site and simplify:
$$I_1 = \frac{\left(1+2\, r^2\right)^{1/3}}{\gamma\left[\frac{1}{3}\right]} \int_0^1 (1-q)^{-\frac{1}{2}} q^{-\frac{1}{2}} \,G_{1,1}^{1,1}\left( \sqrt{q}\, r^2\left\vert
\begin{array}{c}
\frac{2}{3}\\
0%
\end{array}%
\right. \right)\,
\, dq$$
$$ = \left(1+2 \,r^2\right)^{1/3} \left(\pi \,_{2}F_{1}\left[\frac{1}{6},\frac{2}{3},1,r^4\right]-\frac{2}{3} r^2 \,_{3}F_{2}\left[\left\{\frac{2}{3},1,\frac{7}{6}\right\},\
\left\{\frac{3}{2},\frac{3}{2}\right\},r^4\right]\right)$$
The second integral $I_2$ is performed in the same way. The limit $r \to 1$ Mathematical Functions Site and Mathematical Functions Site exists e.g. for $I_1$
$$\lim_{r \to 1}\,I_1=\frac{\frac{3 \pi\, \gamma\
\left[\frac{1}{6}\right]}{\gamma\left[\frac{1}{3}\right]\, \gamma\left[\frac{5}{6}\right]}-2\,_{3}F_{2}\left[\left\{\frac{2}{3},1,\frac{7}{6}\right\},\left\{\frac{3}{2},\frac{3}{2}\right\},1\right]}{3^{2/3}}$$
and for $I_2$ in the same way. The sum of $I = I_1 + I_2$
delivers just the half of your answer:
$$I= \frac{2^{1/3}\, 3^{5/6}\, \gamma\left[\frac{1}{3}\right]^2}{\gamma\left[\frac{2}{3}\right]} $$.
Therefore, inspite of a calculation error of factor 2, @user1055 is right, the limit may be also calculated inside the integral.