Suppose I have a state $|\psi\rangle = U|0\rangle$ where $U$ is a $d$-dimensional unitary sampled uniformly with respect the Haar measure. I'm interested in computing or approximating analytically an average quantity of the form
$$ \mathbb{E}_{U \sim U(d)} \left[ \frac{|\psi_i|^2}{|\psi_i|^2 + |\psi_j|^2} \right] = \int_{U(d)} \frac{|U_{i0}|^2}{|U_{i0}|^2 + |U_{j0}|^2} \mu(dU) \tag{1} $$ for some $0 \leq i,j \leq d-1$ and $\mu(dU)$ denotes the Haar measure. What are possible ways to compute or approximate this?
I figured one possibility is to take the limit $d \rightarrow \infty$ and then use a Gaussian approximation for $\psi_j$ (as explained in the answer to this question, for example). However I am a bit concerned about singularities in this polynomial making it an ill-defined computation somehow.
Are there other approaches? Alternatively a reference for explicit computation of rational polynomials of the form of (1) would be appreciated.