I am not a mathematician nor physicist. I just know the basics of the representation theory. In my research, I realized that there is an orthogonality relation between the unitary group matrix elements as follows:

$$I_1 = \int \mathrm{D}\mathbf{U} \; U_{i j}^{(\mathbf{r})} U_{k l }^{*(\mathbf{r}^{\prime})} = \frac{1}{ d_{ \mathbf{r} } } \delta_{\mathbf{r} \mathbf{r}^{\prime} }  \delta_{i k} \delta_{j l}
$$

where $\mathbf{U} \in \mathcal{U}(N)$, $\mathrm{D}\mathbf{U}$ is the standard Haar measure, $U_{ij}^{(\mathbf{r})}$ denotes the $(i,j)$-th element of the representation matrix of $\mathbf{U}$, and $d_{ \mathbf{r} }$ is the dimension of the irreducible representation $\mathbf{r}$. 

Now, I need to know the answer for this integral:

$$I_2 = \int \mathrm{D} \mathbf{U} \; U_{i_1 j_1}^{(\mathbf{r})} U_{ k_1 l_1 }^{ * ( \mathbf{r} ) } U_{ i_2 j_2 }^{(\mathbf{r}^{\prime})} U_{ k_2 l_2 }^{* ( \mathbf{r}^{ \prime \prime } ) }
$$

I appreciate any help.


p.s. Here is my conjecture for the answer:

$$
I_2 = \delta_{ \mathbf{r}^{\prime} \mathbf{r}^{\prime \prime} } \times \left\{
 \eqalign{
\frac{1}{ d_{ \mathbf{r} } d_{ \mathbf{r}^{ \prime } } -1 } \delta_{ i_1 k_1 } \delta_{ j_1 l_1 } \delta_{ i_2 k_2 } \delta_{ j_2 l_2 } ( 1- \delta_{ \mathbf{r} \mathbf{r}^{\prime} } ) \\
 + \delta_{ \mathbf{r} \mathbf{r}^{\prime} } \left[
     \eqalign{
      \frac{ 1 }{ d_{ \mathbf{r} }^2 -1 }
 ( \delta_{ i_1 k_1 } \delta_{ j_1 l_1 } \delta_{ i_2 k_2 } \delta_{ j_2 l_2 }
 + \delta_{ i_1 k_2 } \delta_{ j_1 l_2 } \delta_{ i_2 k_1 } \delta_{ j_2 l_1 } ) \\
 - \frac{ 1 }{ d_{ \mathbf{r} } ( d_{ \mathbf{r} }^2 -1 ) } ( \delta_{ i_1 k_1 } \delta_{ j_1 l_2 } \delta_{ i_2 k_2 } \delta_{ j_2 l_1 } + \delta_{ i_1 k_2 } \delta_{ j_1 l_1 }  \delta_{ i_2 k_1 } \delta_{ j_2 l_2 } )  }
    \right] 
} \right\} $$

**UPDATE:**

I have been advised that it might be helpful if I can find the tensor product of two irreducible representations, $ \mathbf{s} = \mathbf{r} \otimes \mathbf{r}^{\prime}$, which most likely leads to a reducible representation, and then I need to decompose $\mathbf{s}$ into its irreducible components (by using the Clebsch–Gordan coefficients, according to wikipedia), to be able to use the Schur's lemma to get the answer!!!

However, it is hard for me to do this, and needs awful background.