I want to find a value for the expression:
$$\underset{U}{\mathbb{E}}[U^{\otimes 2} (P_1 \otimes P_2)^{\otimes 2} U^{*\otimes 2}],$$
where $U$ is a two-qubit unitary operator chosen Haar randomly, $P_1$ and $P_2$ are two single qubit Pauli operators, and the expectation is taken over choices of $U$.
When both Paulis are identity, the expression is trivially $1$. However, I wasn't entirely sure of the other cases.
We can prove a slightly more general statement, which works for any Hermitian "sandwiched" matrix $P$ in any dimension $d$... but only for $t=2$ in the number of copies :'(
Let us call $V$ this expectation, and let $W$ be a unitary matrix. We have: $$\begin{align*} W^{\otimes2}V &= \int_U\left(WUPU^\dagger\right)^{\otimes 2}\,\mathrm{d}\mu(U)\\ &= \int_U\left(UPU^\dagger W\right)^{\otimes 2}\,\mathrm{d}\mu(UW)\\ &= \int_U\left(UPU^\dagger W\right)^{\otimes 2}\,\mathrm{d}\mu(U)\\ &= VW^{\otimes 2}. \end{align*}$$ Thus, $V$ commutes with every $W^{\otimes 2}$ for $W\in \mathcal{U}(d)$.
We now want to use Schur's Lemma. We know that $U(d)\ni W\mapsto W^{\otimes2}$ is an irreducible representation of the unitary group in both the symmetric subspace $\vee^2\mathbb{C}^d$ and the antisymmetric one . Furthermore, the former has dimension $\binom{d+1}{2}$ and the latter $\binom{d}{2}$, which conveniently sums up to $d^2$. Thus, we know that $V$ can be written as: $$V=\lambda_1P_{\text{sym}}+\lambda_2P_{\text{antisym}}$$ where $P_{\text{sym}}$ and $P_{\text{antisym}}$ are the respective projectors onto these subspaces.
We can evaluate $\lambda_1$ by computing: $$\lambda_1=\langle00|V|00\rangle=\int_U(\langle0|UPU^\dagger|0\rangle)^2\,\mathrm{d}\mu(U).$$ Since $P$ is Hermitian, $\langle0|UPU^\dagger|0\rangle$ is real. Thus, we can write: $$\lambda_1=\int_U|\langle0|UPU^\dagger|0\rangle|^2\,\mathrm{d}\mu(U)$$ This is a known quantity and we have: $$\lambda_1=\frac{\mathrm{tr}^2(P)+d}{(d+1)d}$$ Note also that we have, once again using the fact that $P$ is Hermitian: $$\mathrm{tr}(V)=\lambda_1\binom{d+1}{2}+\lambda_2\binom{d}{2}=\mathrm{tr}^2(P).$$ Which gives us: $$\lambda_2=\frac{\mathrm{tr}^2(P)-d}{d(d-1)}.$$ If $P$ is a Pauli matrix, it simplifies a lot, since the identity is the only one with a non-nil trace. For $P=I$, we find $\lambda_1=\lambda_2=1$, which gives us $$V=P_{\text{sym}}+P_{\text{anti}}=I,$$ and for any other Pauli matrix this gives us $\lambda_1=\frac{1}{d+1}$ and $\lambda_2=-\frac{1}{d-1}$, which gives us $$V=\frac{1}{d+1}P_{\text{sym}}-\frac{1}{d-1}P_{\text{anti}}.$$
