Let $\mathbf{X}$ follow a matrix variate F distribution with pdf $$ \begin{align} f\left(\mathbf{X} | \mathbf{\Sigma}, n, \nu\right) = \frac{\Gamma_p(\frac{n + \nu }{2})}{\Gamma_p(\frac{n}{2})\Gamma_p(\frac{\nu}{2})} |\mathbf{\Sigma}|^{\frac{\nu}{2}} |\mathbf{X}|^{\frac{n - p - 1}{2}} |\mathbf{\Sigma}+ \mathbf{X}|^{-\frac{n + \nu}{2}}, \end{align} $$ where $\mathbf{X}$ is the symmetric positive definite data matrix, $\mathbf{\Sigma}$ is the sym p.d. scale matrix and $n$ and $\nu$ are scalar parameters.
I am just interested in the Fisher Information Matrix w.r.t. the scale matrix $\mathbf{\Sigma}$.
Now $$ \begin{align} d \log(f\left(\mathbf{X} | \mathbf{\Sigma}, n, \nu\right)) &= \frac{\nu}{2} d \log|\mathbf{\Sigma}| - \frac{n + \nu}{2} d \log |\mathbf{\Sigma}+ \mathbf{X}| \\ &= \frac{\nu}{2} tr( \mathbf{\Sigma}^{-1} d \mathbf{\Sigma}) - \frac{n + \nu}{2} tr((\mathbf{\Sigma}+ \mathbf{X})^{-1}d\mathbf{\Sigma}), \end{align} $$ such that $$ \begin{align} d^2 \log(f\left(\mathbf{X} | \mathbf{\Sigma}, n, \nu\right)) &= \frac{\nu}{2} d tr( \mathbf{\Sigma}^{-1} d \mathbf{\Sigma}) - \frac{n + \nu}{2} d tr((\mathbf{\Sigma}+ \mathbf{X})^{-1}d\mathbf{\Sigma})\\ &= - \frac{\nu}{2} tr( \mathbf{\Sigma}^{-1} d \mathbf{\Sigma}\mathbf{\Sigma}^{-1} d \mathbf{\Sigma}) + \frac{n + \nu}{2} tr((\mathbf{\Sigma}+ \mathbf{X})^{-1} d\mathbf{\Sigma}(\mathbf{\Sigma}+ \mathbf{X})^{-1} d\mathbf{\Sigma})\\ \end{align} $$ so the Hessian of the log-likelihood obtains as, $$ \frac{\partial^2 \log(f\left(\mathbf{X} | \mathbf{\Sigma}, n, \nu\right))}{\partial vech(\mathbf{\Sigma}) \partial vech(\mathbf{\Sigma})^{\top} } = \mathbf{D}^{\top} \left( - \frac{\nu}{2} ( \mathbf{\Sigma}^{-1} \otimes \mathbf{\Sigma}^{-1} ) + \frac{n + \nu}{2} ((\mathbf{\Sigma}+ \mathbf{X})^{-1} \otimes (\mathbf{\Sigma}+ \mathbf{X})^{-1} ) \right) \mathbf{D}, $$ where $\mathbf{D}$ is the duplication matrix.
Now to get to the Fisher Information matrix w.r.t. $\mathbf{\mathbf{\Sigma}}$ I need the expectation $\mathbf{E}\left((\mathbf{\Sigma}+ \mathbf{X})^{-1} \otimes (\mathbf{\Sigma}+ \mathbf{X})^{-1}\right)$ or alternatively
$\mathbf{E}\left(tr((\mathbf{\Sigma}+ \mathbf{X})^{-1} d\mathbf{\Sigma}(\mathbf{\Sigma}+ \mathbf{X})^{-1} d\mathbf{\Sigma})\right)$
Unfortunately I couldn't find these expectations to have been derived anywhere in the literature. I tried to derive the distribution of $(\mathbf{\Sigma}+ \mathbf{X})^{-1}$, but it does not seem to be a well known one. I tried to compute the second expectation mentioned above via the usual formula as follows (ignoring the normalizing constant). $$ \begin{align} &\int_{\mathbf{X} > 0} tr((\mathbf{\Sigma}+ \mathbf{X})^{-1} d\mathbf{\Sigma}(\mathbf{\Sigma}+ \mathbf{X})^{-1} d\mathbf{\Sigma}) |\mathbf{X}|^{\frac{n-p-1}{2}} |\mathbf{\Sigma}+ \mathbf{X}|^{-\frac{n + \nu}{2}} d \mathbf{X} \\ =&\int_{\mathbf{X} > 0} tr((\mathbf{\Sigma}+ \mathbf{X})^{-1} d\mathbf{\Sigma}(\mathbf{\Sigma}+ \mathbf{X})^{-1} d\mathbf{\Sigma}) |\mathbf{X}|^{\frac{n}{2}} |\mathbf{\Sigma}+ \mathbf{X}|^{-\frac{n + \nu}{2}} d \left(\frac{\mathbf{X}}{|\mathbf{X}|^{\frac{p+1}{2}}} \right), \end{align} $$ now it is well known that the measure $\left(\frac{\mathbf{X}}{|\mathbf{X}|^{\frac{p+1}{2}}} \right)$ is invariant to the transformation $\mathbf{X} \rightarrow \mathbf{\Sigma}^{1/2} \mathbf{X} \mathbf{\Sigma}^{1/2}$, so that we can simplify the integral to \begin{align} & |\mathbf{\Sigma}|^{-\frac{\nu}{2}} \int_{\mathbf{X} > 0} tr( (\mathbf{I} + \mathbf{X})^{-1} \mathbf{\Sigma}^{1/2} d\mathbf{\Sigma}\mathbf{\Sigma}^{1/2} (\mathbf{I} + \mathbf{X})^{-1} \mathbf{\Sigma}^{1/2} d\mathbf{\Sigma}\mathbf{\Sigma}^{1/2}) |\mathbf{X}|^{\frac{n}{2}} |\mathbf{I} + \mathbf{X}|^{-\frac{n + \nu}{2}} d \left(\frac{\mathbf{X}}{|\mathbf{X}|^{\frac{p+1}{2}}} \right) \\ =& |\mathbf{\Sigma}|^{-\frac{\nu}{2}} \int_{\mathbf{X} > 0} tr( (\mathbf{I} + \mathbf{X})^{-1} \mathbf{\Sigma}^{1/2} d\mathbf{\Sigma}\mathbf{\Sigma}^{1/2} (\mathbf{I} + \mathbf{X})^{-1} \mathbf{\Sigma}^{1/2} d\mathbf{\Sigma}\mathbf{\Sigma}^{1/2}) |\mathbf{X}|^{\frac{n-p-1}{2}} |\mathbf{I} + \mathbf{X}|^{-\frac{n + \nu}{2}} d \mathbf{X}, \end{align} which means that all I need to get to the solution is to solve this integral: $$ \int_{\mathbf{X} > 0} tr( (\mathbf{I} + \mathbf{X})^{-1} \mathbf{A} (\mathbf{I} + \mathbf{X})^{-1} \mathbf{A}) |\mathbf{X}|^{\frac{n-p-1}{2}} |\mathbf{I} + \mathbf{X}|^{-\frac{n + \nu}{2}} d \mathbf{X}, $$ where $\mathbf{A}$ is just some arbitrary matrix. Here I am stuck again.
Any help would be much appreciated.