Integrating exponential function over all 3D rotation matrices (SO(3)) with Euler angle parameterization. I have a random variable $\mathbf{X}$ following an isotropic matrix normal distribution with mean $\mathbf{M}$, i.e. $\mathrm{vec}(\mathbf{X}) \sim \mathcal{N}(\mathbf{\mathrm{vec}({\mathbf{M})}}, \sigma \mathbf{I})$.
I am, essentially, interested in figuring out how the equivalence classes / orbits induced by the $\mathcal{SO}(3)$ group (the set of all 3D rotation matrices) is distributed.
That is, how the random variable $[\mathbf{X}] = \{\mathbf{X} \mathbf{R}^T \mid \mathbf{R} \in \mathcal{SO}(3)\}$ is distributed.
Working with ZYX Euler angles (and omitting some constants) led me to an integral of the following form, which I would like to solve:
$$
\int_{0}^{2\pi}\int_{-\pi}^\pi \int_{0}^{2\pi}
\cos(\beta)
\exp\left(
\sum_{n=1}^N \left(\mathbf{R}(\alpha,\beta,\gamma)\mathbf{A}_n\right)^T \mathbf{B}_n
\right)
\ d \alpha \ d\beta \ d\gamma
\tag{1}
$$
with $\mathbf{A}, \mathbf{B} \in \mathbb{R}^{N \times 3}$
and
$$\mathbf{R}(\alpha,\beta,\gamma)
=
\begin{bmatrix}
\cos(\gamma) & -\sin(\gamma) & 0 \\
\sin(\gamma) & \cos(\gamma) & 0\\
0 & 0 & 1
\end{bmatrix}
\cdot
\begin{bmatrix}
\cos(\beta) & 0 & \sin(\beta) \\
0 & 1 & 0 \\
-\sin(\beta) & 0 & \sin(\beta)
\end{bmatrix}
\cdot
\begin{bmatrix}
1 & 0 & 0\\
0 & \cos(\alpha) & -\sin(\alpha) \\
0 & \sin(\alpha) & \cos(\alpha)
\end{bmatrix}.
$$

I am aware that Eq. 1 is (up to some normalizing constant) equivalent to a Haar integral of the form
$$
\int_{\mathcal{SO}(3)} f(\mathbf{R}) \mu(\mathbf{R}).
$$
But I do not know how to use that information.
There exist some posts (e.g. 1, 2, 3) about Haar integration for $\mathcal{SO}(3)$, but they don't discuss a specific function or how to solve the integral.
There is also a question about numerically integrating over Euler angles, but I am interested in an analytical solution.

I have the following two questions
1.) Is there some way of deriving a closed-form analytic expression for Eq. 1?
2.) Or does the integral in Eq. 1 correspond to some special function (e.g. Bessel function?)

I think you are over thinking this: The original distribution is the same as $\bf{X}\sim \bf{M} + \bf{Y}$ where $\bf M$ is the deterministic mean, and $\bf Y$ is distributed as an isotropic Gaussian. These two terms are idependent. that applying a rotation $\bf R$ to $\bf X$ will yield  $ \bf{X}' = \bf{RM} + \bf{RY}$. Since you are asking about what happens when $R$ is isotropic, then $ \bf{X}'$ is rotationally invariant, hence $P({\bf X}')\equiv p(|{\bf X}'|)\equiv p(|{\bf X}|) $ where $p$ is a scalar function. To compute $p$ we are left with computing the distribution of $|{\bf X}| = |\bf{M} + \bf{Y}|.$
