The integral of interest is:
$$I = \int_{\mathbb{R}^3}\int_{\mathbb{R}^3} \boldsymbol{1}\left(\frac{1}{2}\frac{(x_2^2 - x_1^2) + (y_2^2 - y_1^2) + (z_2^2 -z_1^2)}{x_2-x_1} \in [0,1]\right) \nonumber \\ \times \boldsymbol{1}\left(2\mathrm{arcsin}\left(\frac{\sqrt{(x_2-x_1)^2 + (y_2-y_1)^2 + (z_2-z_1)^2}}{2\sqrt{z_1^2 + y_1^2 + \left(x_1 - \frac{1}{2}\frac{x_2^2 - x_1^2 + y_2^2 - y_1^2 + z_2^2 - z_1^2}{x_2-x_1}\right)^2}}\right) > \tau\right) \\ \times \exp\left(-C\left(z_1^2 + y_1^2 + \left(x_1 - \frac{1}{2}\frac{x_2^2 - x_1^2 + y_2^2 - y_1^2 +z_2^2 -z_1^2}{x_2-x_1}\right)^{2}\right)^{3/2}\right)\mathrm{d}Q_1\mathrm{d}Q_2,$$ where $C > 0$, $\tau \in (0, \pi]$, $\boldsymbol{1}(\cdot)$ is the indicator function, $Q_1 = (x_1, y_1, z_1)$ and $Q_2 = (x_2, y_2, z_2)$.
How does one solve this integral or at least simplify to have the expression with the least possible number of integrals?
My attempt: I tried the following transformations:
Substituting $t_1 = x_2 + x_1$, $t_2 = x_2 - x_1$, $t_3 = y_2 + y_1$, $t_4 = y_2 - y_1$, and $t_5 = z_2 + z_1$, $t_6 = z_2 - z_1$, and then eliminating the first indicator function, we get
$$I = \frac{1}{4}\int_{\mathbb{R}^5}\boldsymbol{1}\left(\mathrm{arcsin}\frac{\rho}{2 D} > \frac{\tau}{2}\right)\exp\left(-CD^3\right)\mathrm{d}t_2\mathrm{d}t_3\mathrm{d}t_4\mathrm{d}t_5\mathrm{d}t_6,$$ $\rho = \sqrt{t_2^2 + t_4^2 + t_6^2}$ and $$D = \frac{1}{2}\sqrt{t_3^2 + t_5^2 + \frac{t_3^2t_4^2}{t_2^2}+\frac{t_5^2t_6^2}{t_2^2} + \frac{2t_3t_4t_5t_6}{t_2^2} + \rho^2}.$$
Now, we do the transformation to the spherical coordinates as $t_2 = \rho \sin \theta \cos \varphi$, $t_4 =\rho \sin \theta \sin \varphi$, and $t_6 = \rho\cos \theta$. Then, we have
$$D = \frac{1}{2}\sqrt{t_3^2\left(1 + \tan^2 \varphi\right) + t_5^2\left(1 + \frac{\cot^2 \theta}{\cos^2 \varphi}\right) + \frac{2t_3t_5 \cot \theta \tan \varphi}{\cos \varphi} + \rho^2}$$ and $$\mathrm{d}t_2 \mathrm{d}t_4\mathrm{d}t_6 \mathrm{d}t_3 \mathrm{d}t_5 = \rho^2 \sin\theta \mathrm{d}\theta \mathrm{d}\varphi \mathrm{d}r\mathrm{d}t_3 \mathrm{d}t_5.$$
After doing the substitution, the resulting integration looks still difficult to further simplify.
You already made it to spherical coordinates. Using $r$ instead of $\rho$ and introducing polar coordinates for $t_3=\rho \cos \phi$ and $t_5=\rho \sin \phi$, the integral becomes $$I = \frac{1}{4} \int_0^{2\pi}\mathrm{d}\varphi\int_0^\pi\mathrm{d}\theta \int_0^\infty \mathrm{d}r\, r^2\sin \theta \int_0^{2\pi} \mathrm{d}\phi \int_0^\infty \rho \mathrm{d}\rho \boldsymbol{1}\left(\mathrm{arcsin}\frac{r}{2 D} > \frac{\tau}{2}\right)\exp\left(-CD^3\right)$$ where $$D=\frac{1}{2}\sqrt{\frac{1-\left(\sin\theta \sin\phi \sin\varphi - \cos\theta \cos\phi\right)^2}{\sin^2\theta \cos^2\varphi} \, \rho^2 + r^2} \equiv \frac{1}{2} \sqrt{a^2\rho^2+r^2} \, .$$ In order to carry out the $\rho$ integral, we replace the indicator function by bounds for $\rho$: $\rho=0$ leads to ${\rm arcsin}(1)=\frac{\pi}{2}>\frac{\tau}{2}$. The other bound follows from $$\frac{r^2}{a^2\rho^2+r^2}=\sin^2\frac{\tau}{2} \qquad \Rightarrow \qquad a\rho=r\cot\frac{\tau}{2} \, .$$ Hence the inner $\rho$ integral together with the $r$ integral becomes $$\int_0^\infty {\rm d}r \, r^2 \int_0^{\frac{r}{a}\cot \frac{\tau}{2}} {\rm d}\rho \, \rho \, e^{-\frac{C}{8}(a^2\rho^2+r^2)^{3/2}} \stackrel{u=\frac{C}{8}(a^2\rho^2+r^2)^{3/2}}{=} \frac{4}{3a^2C^{2/3}} \int_0^\infty {\rm d}r \, r^2 \int_{\frac{Cr^3}{8}}^{\frac{Cr^3}{8\sin^3\frac{\tau}{2}}} u^{-1/3} e^{-u} \, {\rm d}u \\ \stackrel{v=\frac{Cr^3}{8}}{=} \frac{32}{9a^2C^{5/3}} \int_0^\infty {\rm d}v \underbrace{\int_v^{\frac{v}{\sin^3\frac{\tau}{2}}} u^{-1/3} e^{-u} \, {\rm d}u}_{=f(v)} = \frac{32}{9a^2C^{5/3}} \left\{ v f(v) \Big|_0^\infty - \int_0^\infty v f'(v) \, {\rm d}v \right\} \\ =\frac{32}{9a^2C^{5/3}} \int_0^\infty {\rm d}v \left( v^{2/3} e^{-v} - \frac{v^{2/3}}{\sin^2 \frac{\tau}{2}} \, e^{-\frac{v}{\sin^3\frac{\tau}{2}}}\right) \\ \stackrel{\text{second term: }v\rightarrow v\sin^3\frac{\tau}{2}}{=} \frac{32}{9a^2C^{5/3}} \left( 1 - \sin^3 \frac{\tau}{2}\right) \int_0^\infty {\rm d}v \, v^{2/3} e^{-v} = \frac{32\,\Gamma(5/3)}{9a^2C^{5/3}} \left( 1 - \sin^3 \frac{\tau}{2}\right) \, .$$
It remains to calculate the angular integrals $$I=\frac{8\,\Gamma(5/3)}{9\,C^{5/3}} \left(1-\sin^3\frac{\tau}{2}\right) \int_0^{2\pi} {\rm d}\varphi \int_0^\pi {\rm d}\theta \, \sin\theta \int_0^{2\pi} {\rm d}\phi \, \frac{1}{a^2} \\ =\frac{8\,\Gamma(5/3)}{9\,C^{5/3}} \left(1-\sin^3\frac{\tau}{2}\right) \int_0^{2\pi} {\rm d}\varphi \int_0^\pi {\rm d}\theta \, \sin\theta \int_0^{2\pi} {\rm d}\phi \, \frac{\sin^2\theta \cos^2\varphi}{1-\left(\sin\theta \sin\phi \sin\varphi - \cos\theta \cos\phi\right)^2} \\ =\frac{16\pi\,\Gamma(5/3)}{9\,C^{5/3}} \left(1-\sin^3\frac{\tau}{2}\right) \int_0^{2\pi} {\rm d}\varphi \int_0^\pi {\rm d}\theta \, \sin^2\theta \, |\cos \varphi| \\ =\frac{32\pi^2\,\Gamma(5/3)}{9\,C^{5/3}} \left(1-\sin^3\frac{\tau}{2}\right) \, .$$