By symmetry of the integrand under the interchange $v_1 \leftrightarrow v_2$ we have
$$I = 2\int_0^\infty\int_0^{v_1}(v_2v_2)^2e^{-a(v_1^2+v_2^2)}(v_1-v_2)\:dv_2\:dv_1$$
Then in polar coordinates we have
$$I = 2\int_0^\infty dr \:r^6e^{-ar^2}\int_0^{\frac{\pi}{4}}d\theta\:\cos^2\theta\sin^2\theta(\cos\theta-\sin\theta)$$
The first integral evaluates to
$$2\int_0^\infty dr \:r^6e^{-ar^2} = a^{-\frac{7}{2}}\int_0^\infty u^{\frac{5}{2}}e^{-u}\:du$$
$$ = a^{-\frac{7}{2}}\Gamma\left(\frac{7}{2}\right) = a^{-\frac{7}{2}}\cdot\frac{5}{2}\cdot\frac{3}{2}\cdot\frac{1}{2}\cdot\Gamma\left(\frac{1}{2}\right) = \frac{15}{8}\sqrt{\frac{\pi}{a^7}}$$
The second integral evaluates to
$$\int_0^{\frac{1}{\sqrt{2}}} x^2-x^4\:dx - \int_{\frac{1}{\sqrt{2}}}^1x^2-x^4\:dx = \frac{7}{30}\frac{1}{\sqrt{2}}-\frac{2}{15}$$
which means the final answer is the product of the two results.