Representation theorem for matrices (reference request)

My comments converted to an answer:

1st comment. I know no reference but the proof need not be [as per OP] lengthy — one just computes the two scalars by which the map $\mathscr I:$ $$ \textstyle A\mapsto\mathscr I(A)=\int_S h\langle h,Ah\rangle\langle h,\cdot\rangle d\lambda(h) \tag1 $$ acts on the irreducible components of $\mathfrak{gl}(n,\mathbf C)=\mathfrak{sl}(n,\mathbf C)\oplus\mathbf C I$ (Schur’s lemma).

(In detail: write $r$ and $s$ for the scalars in question. Taking the trace of $sI=\mathscr I(I)$ in (1) gives $s=1$. It follows that we have $ \mathscr I(A)= r\left(A-\tfrac{\operatorname{tr}A}{n}I\right)+\tfrac{\operatorname{tr}A}{n}I $ and hence $$ \operatorname{tr}(\mathscr I(A)B)=r\operatorname{tr}(AB)+\tfrac{1-r}n\operatorname{tr}(A)\operatorname{tr}(B)\rlap{\qquad\quad\forall A, B.} \tag2 $$ Writing this out for $(A,B)=(E_{12},E_{21})$, resp. $(E_{11},E_{22})$ where $E_{ij}=e_i\langle e_j,\cdot\rangle$, one obtains that $\int_S|\langle e_1,h\rangle|^2|\langle e_2,h\rangle|^2d\lambda(h)$ equals both $r$ and $\frac{1-r}n$. Therefore $r=\frac1{n+1}$, and with that (2) becomes (3) below. QED)

2nd comment. For a reference: your desired formula is equivalent to $$ \mathrm{tr}(AB)+\mathrm{tr}(A)\mathrm{tr}(B)=(n+1)\int_S\langle h,Ah\rangle\langle h,Bh\rangle\,d\lambda(h) \tag3 $$ which is e.g. (3.8) of Gibbons (1992), taking differences of notation into account. (The integrand descends to $P^{n-1}(\mathbf C)$ where he uses the measure of total volume $\pi^{n-1}/(n-1)!$) I think it should also follow from Archimedes-Duistermaat-Heckman (1982, Prop. 3.2).

More context: Generally I think you’ll find many similar formulas in the (somewhat repetitive) literature on “coherent states” or “quantum mechanics as classical mechanics on $P\mathscr H$”. Also compare Schur’s proof of his orthogonality relations (1924, p.199 or Bröcker-tom Dieck 4.5i), specialized to the adjoint representation.

This is an easy consequence of the $k=2$ case of the complex version of the Isserlis-Wick theorem for moments of Gaussian measures, i.e., the identity $$ \int_{\mathbb{C}^n} z_{i_1}\cdots z_{i_k}\ {\bar{z}}_{j_1}\cdots {\bar{z}}_{j_k} \ e^{-|z|^2} \prod_{a=1}^{n}\frac{d(\Re z_a) d(\Im z_a)}{\pi}\ =\ \sum_{\sigma\in\mathfrak{S}_k} \delta_{i_1 j_{\sigma(1)}}\cdots \delta_{i_k j_{\sigma(k)}}\ . $$

Going to spherical coordinates produces the integral $\int_{S}\cdots d\lambda(h)$, while the sum over the permutation $\sigma$ gives the other two terms of the wanted identity.

Indeed, it is good to see this really as an integral over $\mathbb{C}\mathbb{P}^{n-1}$ for the Fubini-Study metric, but appealing to Duistermaat-Heckman is not necessary (as known to François).