Integral is as follows:
$ I = \int d \Omega_x \ \int d \Omega_y \ \hat{a} \cdot \hat {x}\ \hat{x} \cdot \hat {y} \ \hat{y} \cdot \hat {b}\ \hat{b} \cdot \hat {x} \ \hat{x} \cdot \hat {y} \ \hat{y} \cdot \hat {a}$
0- I have tried some educated guesses such as $ I = C (\hat{a} \cdot \hat{b} )^2 $ and choose $\hat{a}$ and $\hat{b}$ to be unit vector $\hat{k}$. But since I don't have the answer I can't confirm if this solution makes any sense or not.
1- For the first step I assumed $\hat{x}$ to be constant while integrating with respect to $\hat{y}$ .Then I got:
$ \int d \Omega_y \ \hat{x} \cdot \hat {y} \ \hat{y} \cdot \hat {b} \ \hat{x} \cdot \hat {y} \ \hat{y} \cdot \hat {a}$ converted vectors to polar coordinates but since I cannot choose any direction for $\hat {x}$ and $\hat {y}$ the question has become much more complicated, and this is where I am stuck right now.
For those interested the question is from: An Introduction to Mathematical Methods of Physics (Lorella M. Jones), Section 6 and question 16.
I check the book and it seems my guess is correct, that is, $\mathrm{d}\Omega$ represents the surface measure for the unit sphere $\mathbb{S}^2$. Based on this interpretation, I will show that
$$ I = \left( \frac{4\pi}{15} \right)^2 (2 + 9 \langle \hat{\mathbf{a}}, \hat{\mathbf{b}} \rangle^2). $$
1. (Preliminary) To begin with, we first replace $\hat{\mathbf{a}}$ and $\hat{\mathbf{b}}$ by $\mathbf{a}$ and $\mathbf{b}$ without affecting the essence of computation, respectively, hence we consider the integral
$$ I(\mathbf{a}; \mathbf{b}) = \int_{\mathbb{S}^2} \mathrm{d}\Omega_{x} \int_{\mathbb{S}^2} \mathrm{d}\Omega_{y} \, \langle \mathbf{a}, \hat{\mathbf{x}} \rangle \langle \hat{\mathbf{x}}, \hat{\mathbf{y}} \rangle \langle \hat{\mathbf{y}}, \mathbf{b} \rangle \langle \mathbf{b}, \hat{\mathbf{x}} \rangle \langle \hat{\mathbf{x}}, \hat{\mathbf{y}} \rangle \langle \hat{\mathbf{y}}, \mathbf{a} \rangle. $$
(In this answer, we use $\langle \cdot, \cdot \rangle$ for inner product and reserve $\cdot$ for scalar multiplication.) Next, we introduce a de-symmetrized version of the integral
$$ J(\mathbf{a}; \mathbf{b}, \mathbf{c}) = \int_{\mathbb{S}^2} \mathrm{d}\Omega_{x} \int_{\mathbb{S}^2} \mathrm{d}\Omega_{y} \, \langle \mathbf{a}, \hat{\mathbf{x}} \rangle \langle \hat{\mathbf{x}}, \hat{\mathbf{y}} \rangle \langle \hat{\mathbf{y}}, \color{red}{\mathbf{b}} \rangle \langle \color{red}{\mathbf{c}}, \hat{\mathbf{x}} \rangle \langle \hat{\mathbf{x}}, \hat{\mathbf{y}} \rangle \langle \hat{\mathbf{y}}, \mathbf{a} \rangle. $$
It is easy to verify the following properties:
$I(\mathbf{b}; \mathbf{a}) = I(\mathbf{a}; \mathbf{b})$.
$ I(\mathbf{a}; \mathbf{b}) = J(\mathbf{a}; \mathbf{b}, \mathbf{b}) $ by definition.
$J(\mathbf{a}; \mathbf{c}, \mathbf{b}) = J(\mathbf{a}; \mathbf{b}, \mathbf{c})$. This is easily proved by swapping the role of $\hat{\mathbf{x}} $ and $\hat{\mathbf{y}}$.
$ I(\mathbf{a}; \mathbf{b}+\mathbf{c}) = I(\mathbf{a}; \mathbf{b}) + 2J(\mathbf{a}; \mathbf{b}, \mathbf{c}) + I(\mathbf{a}; \mathbf{c}) $.
2. (Simplification by Rotational Symmetry) Now we move on to simplifying $I(\hat{\mathbf{a}}; \hat{\mathbf{b}})$.
Note that rotating all the vectors $\mathbf{a}$, $\hat{\mathbf{x}}$, $\hat{\mathbf{y}}$, $\mathbf{b}$ in the integrand of $I(\mathbf{a}; \mathbf{b})$ by the same rotational transform leaves the integral unchanged. So by choosing a suitable rotation, we may assume that
$$ \hat{\mathbf{a}} = \hat{\mathbf{x}}, \qquad \hat{\mathbf{b}} = (\cos \alpha) \hat{\mathbf{x}} + (\sin\alpha) \hat{\mathbf{y}}, $$
where $\alpha = \arccos\langle \hat{\mathbf{a}}, \hat{\mathbf{b}} \rangle$ is the angle between $\hat{\mathbf{a}}$ and $\hat{\mathbf{b}}$. Plugging this back to the integral, we get
\begin{align*} I(\hat{\mathbf{a}}; \hat{\mathbf{b}}) &= (\cos^2\alpha) I(\hat{\mathbf{x}}; \hat{\mathbf{x}}) + 2(\cos\alpha\sin\alpha) J(\hat{\mathbf{x}}; \hat{\mathbf{x}}, \hat{\mathbf{y}}) + (\sin^2\alpha) I(\hat{\mathbf{x}}; \hat{\mathbf{y}}). \tag{1} \end{align*}
3. (Simplification by Radial Symmetry) Next, we observe that
$$ \frac{15}{4\pi} = \int_{0}^{\infty} r^2 \, \mathrm{d}r \cdot r^4 \frac{e^{-r^2/2}}{(2\pi)^{3/2}}. $$
What is the use of this seemingly arbitrary integral? It allows us to introduce radial direction back to the integral. Together with the relations
$$ r \cdot \hat{\mathbf{r}} = \mathbf{r} \qquad\text{and} \qquad r^2 \mathrm{d}r\mathrm{d}\Omega = \mathrm{d}^3\mathbf{r}, $$
we can convert the integral into cartesian form:
\begin{align*} \left( \frac{15}{4\pi} \right)^2 I(\hat{\mathbf{a}}, \hat{\mathbf{b}}) &= \int_{\mathbb{R}^3} \mathrm{d}^3\mathbf{r}_1 \int_{\mathbb{R}^3} \mathrm{d}^3\mathbf{r}_2 \, \frac{e^{-\|\mathbf{r}_1\|^2/2}}{(2\pi)^{3/2}} \frac{e^{-\|\mathbf{r}_2\|^2/2}}{(2\pi)^{3/2}} \\ &\hspace{9em}\cdot \langle \hat{\mathbf{a}}, \mathbf{r}_1 \rangle \langle \mathbf{r}_1, \mathbf{r}_2 \rangle \langle \mathbf{r}_2, \hat{\mathbf{b}} \rangle \langle \hat{\mathbf{b}}, \mathbf{r}_1 \rangle \langle \mathbf{r}_1, \mathbf{r}_2 \rangle \langle \mathbf{r}_2, \hat{\mathbf{a}} \rangle \\ &= \mathbf{E}[\langle\hat{\mathbf{a}},X\rangle\langle X,Y\rangle\langle Y,\hat{\mathbf{b}}\rangle\langle\hat{\mathbf{b}},X\rangle\langle X,Y\rangle\langle Y,\hat{\mathbf{a}}\rangle], \end{align*}
where $X = (X_1, X_2, X_3)$ and $Y = (Y_1, Y_2, Y_3)$ are independent standard gaussian vectors in $\mathbb{R}^3$, i.e., $X, Y \sim \mathcal{N}(\mathbf{0},\mathbf{I}_3)$. Expanding the expectation using linearity,
\begin{align*} I(\hat{\mathbf{a}}, \hat{\mathbf{b}}) &= \left( \frac{4\pi}{15} \right)^2 \sum_{i_1, i_2, i_3, i_4, i_5, i_6} \mathbf{E}[\hat{a}_{i_1}X_{i_1}X_{i_2}Y_{i_2}Y_{i_3}\hat{b}_{i_3}\hat{b}_{i_4}X_{i_4}X_{i_5}Y_{i_5}Y_{i_6}\hat{a}_{i_6}] \\ &= \left( \frac{4\pi}{15} \right)^2 \sum_{i_1, i_2, i_3, i_4, i_5, i_6} \hat{a}_{i_1}\hat{b}_{i_3}\hat{b}_{i_4} \hat{a}_{i_6}\mathbf{E}[X_{i_1}X_{i_2}X_{i_4}X_{i_5}]\mathbf{E}[Y_{i_2}Y_{i_3}Y_{i_5}Y_{i_6}]. \tag{2} \end{align*}
4. (Evaluation of the Integral) Now we evaluate each term in $\text{(1)}$ using $\text{(2)}$.
Case 1. We first evaluate $I(\hat{\mathbf{x}}; \hat{\mathbf{x}})$. The above formula $\text{(2)}$ gives
\begin{align*} I(\hat{\mathbf{x}}; \hat{\mathbf{x}}) &= \left( \frac{4\pi}{15} \right)^2 \sum_{i_2, i_5} \mathbf{E}[X_1^2X_{i_2}X_{i_5}]\mathbf{E}[Y_1^2 Y_{i_2}Y_{i_5}]. \end{align*}
Since $\mathbf{E}[X_i] = \mathbf{E}[X_i^3] = 0$, the sum only survives when $i_2 = i_5$. Hence,
\begin{align*} I(\hat{\mathbf{x}}; \hat{\mathbf{x}}) &= \left( \frac{4\pi}{15} \right)^2 \sum_{i} \mathbf{E}[X_1^2X_i^2]\mathbf{E}[Y_1^2 Y_i^2] \\ &= \left( \frac{4\pi}{15} \right)^2 (\mathbf{E}[X_1^4]\mathbf{E}[Y_1^4] + \mathbf{E}[X_1^2X_2^2]\mathbf{E}[Y_1^2Y_2^2] + \mathbf{E}[X_1^2X_3^2]\mathbf{E}[Y_1^2Y_3^2]) \\ &= \left( \frac{4\pi}{15} \right)^2 \cdot 11 \end{align*}
In the last line, we utilized the independence of $X_1, X_2, X_3, Y_1, Y_2, Y_3$ and the relations $\mathbf{E}[X_1^2] = 1$ and $\mathbf{E}[X_1^4] = 3$.
Case 2. Next, we evaluate $I(\hat{\mathbf{x}}; \hat{\mathbf{y}})$. Invoking the general formula $\text{(2)}$ again, we get
\begin{align*} I(\hat{\mathbf{x}}; \hat{\mathbf{x}}) &= \left( \frac{4\pi}{15} \right)^2 \sum_{i_2, i_5} \mathbf{E}[X_1 X_2 X_{i_2}X_{i_5}]\mathbf{E}[Y_1 Y_2 Y_{i_2}Y_{i_5}] \\ &= \left( \frac{4\pi}{15} \right)^2 2\mathbf{E}[X_1^2 X_2^2 ]\mathbf{E}[Y_1^2 Y_2^2] \\ &= \left( \frac{4\pi}{15} \right)^2 \cdot 2 \end{align*}
Case 3. Finally, we turn to the evaluatio nof $J(\hat{\mathbf{x}}; \hat{\mathbf{x}}, \hat{\mathbf{y}})$. By $\text{(2)}$,
\begin{align*} J(\hat{\mathbf{x}}; \hat{\mathbf{x}}, \hat{\mathbf{y}}) &= \left( \frac{4\pi}{15} \right)^2 \sum_{i_2, i_5} \mathbf{E}[X_1^2 X_{i_2}X_{i_5}]\mathbf{E}[Y_1 Y_2 Y_{i_2}Y_{i_5}]. \end{align*}
For the term $\mathbf{E}[X_1^2 X_{i_2}X_{i_5}]$ to survive, we must have $i_2 = i_5$. On the other hand, $\mathbf{E}[Y_1 Y_2 Y_{i_2}Y_{i_5}]$ is non-zero only when $\{i_2, i_5\} = \{1, 2\}$. Since both conditions cannot hold simultaneously, it follows that
$$ J(\hat{\mathbf{x}}; \hat{\mathbf{x}}, \hat{\mathbf{y}}) = 0. $$
(Alternatively, apply the reflection about the $xz$-plane to argue that $J(\hat{\mathbf{x}}; \hat{\mathbf{x}}, \hat{\mathbf{y}}) = J(\hat{\mathbf{x}}; \hat{\mathbf{x}}, -\hat{\mathbf{y}})$, which in turn implies the above equality.)
5. Conclusion. Combining all the computations together, we obtain
\begin{align*} I(\hat{\mathbf{a}}; \hat{\mathbf{b}}) &= \left( \frac{4\pi}{15} \right)^2 (11 \cos^2\alpha + 2 \sin^2\alpha) \\ &= \left( \frac{4\pi}{15} \right)^2 (2 + 9 \langle \hat{\mathbf{a}}, \hat{\mathbf{b}} \rangle^2). \end{align*}
Below is a Monte Carlo simulation of the integral for a randomly chosen $\hat{\mathbf{a}}$ and $\hat{\mathbf{b}}$, hinting that the above answer indeed agrees with the integral