What is the expectation value of $|\langle \psi|U|\psi \rangle|$ over Haar random states $|\psi\rangle$?

Question

We know the average unitary fidelity, $\int |\langle \psi|U|\psi \rangle|^2 d\psi$, has a nice closed-form solution: $\frac{1+\frac{1}{d}|Tr (U)|^2}{1+d}$, thanks to Horodecki and Nielsen.

However, I am looking for an analytical expression for the Haar average $\int |\langle \psi|U|\psi \rangle| d\psi$, where $|\psi \rangle\in \mathbb{C}^d$. Is it computable? I am particularly interested in the case where the unitary $U$ is traceless.

Additionally, it would be helpful to have some useful references on this.

Remarks:

1. Without the absolute, one can calculate the integral $\int \langle \psi|U|\psi \rangle d\psi$: $U\int|\psi\rangle\langle\psi|d\psi{=}U.Id/d$, hence taking trace on both sides, we have the integral reduced to $Tr(U)/d$. This is not what I am looking for.
2. One can use the Cauchy-Schwartz inequality on the average gate fidelity to derive an upper bound of the integral in question. However, a closed-form solution would have been helpful.

Answer 1

This partial answer calculates the integral for $d=2$. In this case, every traceless unitary $U$ is equivalent to the Pauli $Z$ up to similarity and global phase, so, by rotational invariance of the Haar measure, we have $$ \begin{align} \int|\langle\psi|U|\psi\rangle|\,d\psi&=\int |\langle \psi|Z|\psi \rangle| \,d\psi\tag1\\ &=\int|\psi_0\overline{\psi_0}-\psi_1\overline\psi_1|\,d\psi\tag2\\ &=\frac{1}{4\pi}\int_0^\pi\int_0^{2\pi}\left|\cos^2\frac{\theta}{2}-\sin^2\frac{\theta}{2}\right|\sin\theta\,d\phi\,d\theta\tag3\\ &=\frac{1}{2}\int_0^\pi\left|\cos\theta\right|\sin\theta\,d\theta\tag4\\ &=\int_0^{\pi/2}\sin\theta\cos\theta\,d\theta\tag5\\ &=\frac{1}{2}\int_0^{\pi/2}\sin 2\theta\,d\theta\tag6\\ &=\frac{1}{2}.\tag7 \end{align} $$ This agrees with the bound $$ \int|\langle\psi|U|\psi\rangle|\,d\psi\leqslant\frac{1}{\sqrt{d+1}}\tag8 $$ obtained from Jensen's inequality and Horodecki's formula. It should be possible to generalize the calculation above. For example, for $d=3$, we have $$ \begin{align} \int|\langle\psi|U|\psi\rangle|\,d\psi&=\int |\langle \psi|Z|\psi \rangle| \,d\psi\tag9\\ &=\int|\psi_0\overline{\psi_0}+\omega\psi_1\overline\psi_1+\omega^2\psi_2\overline\psi_2|\,d\psi\tag{10} \end{align} $$ where $Z=\mathrm{diag}(1,\omega,\omega^2)$ and $\omega=e^{2\pi i/3}$. One could use the parametrization of $\psi$ from this answer. Note that for larger $d$ the spectrum of $U$ has continuous degrees of freedom. For example, for $d=4$ it can take the form $\mathrm{diag}(1, -1, z, -z)$ for any $z\in\mathbb{C}$ with $|z|=1$.

Answer 2

As a follow-up discussion with Adam, I evaluated the integral till $d=4$ for high-dimensional unitary $Z$, with ${Z}_{kl}{=}\exp(\frac{i2\pi k}{d})\delta_{kl}$. For $d{=}4$, I haven't considered the other possibility that Adam mentioned, i.e., $\tilde{U}{=}V_1\sigma_zV_1^\dagger{\oplus}e^{i\alpha}V_2\sigma_zV_2^\dagger$, $\sigma_z$ is Pauli-z.

I referred to this paper to parameterise the normalized, unitary invariant measure $d\psi$ on the manifold of $|\psi\rangle{=}\sum_{j=1}^d \sqrt{r_j}e^{i\theta_j}|j\rangle$ ($r_j{\geq}0$, $\theta_j{\in}[0,2\pi]$). The parameterisation is introduced after Section V, Lemma 4. For completeness, the parameterisation for $d\psi$ is given by the following Dirac-delta representation:

$$ \begin{align} d{\psi}\equiv\frac{\Gamma(d)}{2\pi^d}\delta\Big(1-\sum_{j=1}^d r_j\Big)\ \prod_{j=1}^d dr_jd\theta_j. \end{align} $$ Here $\Gamma(d){=}(d-1)!$ is the Gamma function. First, let us check if the above parameterisation agrees with the $d{=}2$ case. In this case, $|\langle \psi |Z|\psi\rangle|{=}|r_1-r_2|$, and the integral $I_{d=2}$ becomes $$ \begin{align} I_{d=2}&= \frac{1}{4\pi^2}\int_{r_1{=}0}^1\int_{r_2{=}0}^1 |r_1-r_2|\delta\Big(1-r_1-r_2\Big)dr_1dr_2 \int_{0}^{2\pi}\int_{0}^{2\pi} d\theta_1d\theta_2 \\ &=\int_{r_1,r_2{=}0}^1 |r_1-r_2|\delta\Big(1-r_1-r_2\Big)dr_1dr_2\\ &{=}\frac{1}{2}. \end{align} $$ Then I moved on to $d{=}3$, in which case, $|\langle \psi |Z|\psi\rangle|{=}|r_1+\omega r_2 + \omega^2 r_3|$, with $\omega{=}\exp(\frac{i2\pi}{3})$. With little intention to work it out myself, I fed it to Mathematica and found $I_{d=3}=\frac{1}{3}+\frac{\ln(2+\sqrt{3})}{6\sqrt{3}}{\approx}0.460058$, not conforming to my guess, $\frac{1}{\sqrt{d+2}}$. Similarly, I calculated for $I_{d=4}{\approx}0.405806$ with a clumsy analytical expression as shown in the attached screenshot.

I cross-verified the evaluation with numerical estimations with $10^6$ samples. Here are the list of numerical results for $d{\in}[2,6]$, $[0.4998, 0.4600, 0.4058, 0.3695, 0.3409]$.

So at this stage, I don't have much intuition about the general analytical expression. However the upper bound of $\frac{1}{\sqrt{d+1}}$ is helpful. I will be open to further interesting perspectives on this.