What is the average value of $|c_i\bar c_j|$ for a random state $|\psi\rangle=\sum_i c_i|i\rangle$?

Question

Consider the density matrix $\rho=|\psi\rangle\!\langle\psi|$ of a random pure state in an $N$-dimensional space (in other words, an $N$-dimensional qudit, $|\psi\rangle\in\mathbb C^N$), $\rho_{ij}=c_i \bar c_j$.

The average value of the diagonal elements of such density matrix, $|c_i|^2$, can be seen to be $\mathbb E[\rho_{ii}]=1/N$, when sampling the states according to the uniform Haar random distribution (that is, for example, taking the states as a column of a random unitary matrix of the appropriate dimensions).

It seems also natural that the average off-diagonal terms vanish: $\mathbb E[\rho_{ij}]=\mathbb E[c_i \bar c_j]=\mathbb E[c_i c_j]=0$ for $i<j$.

What about $\mathbb E[\lvert c_i \bar c_j\rvert]$? Can anything be said in this case?

A quick numerical experiment, plotting this average for different values of $N$ in loglog scale, gives the following

which looks very much like a $C/N$ behaviour for some constant $C$.

---

The above plot can be obtained with the following python code

import numpy as np
import qutip
import matplotlib.pyplot as plt
N = 200
state_sizes_list = [2, 5, 10, 20, 30, 50, 70, 100, 200, 300]
averages = []
for state_size in state_sizes_list:
  cij_samples = np.zeros(shape=N)
  for idx in range(N):
    ket = qutip.rand_ket_haar(N=state_size).full()
    cij_samples[idx] = np.abs(ket[0] * ket[1])
  averages.append(cij_samples.mean())
plt.loglog(state_sizes_list, averages, 'ro-')
plt.xlabel('N');

Answer 1

The constant is $C = \frac{\pi}{4}$.

The $2N$-dimensional vector composed of the real and imaginary parts of the coefficients of a normalized state vector is uniformly distributed on an $S^{2N-1}$. Thus the needed expectation value is:

$$\frac{C}{N} = \frac{\int_{S^{2N-1}}|(c_{1x} + i c_{1y}) (c_{2x} - i c_{2y})| d\mu(S^{2N-1})}{\mathrm{Vol}(S^{2N-1})}$$

Integrals of homogeneous functions over spheres can be traded by Gaussian integrals over the ambient Euclidean space the with proper normalizations; since the angular integrations are the same and we need only to normalize the different radial integrals. In our case we get:

$$\frac{C}{N} = \frac{\int_{\mathbb{C}^N}|(c_{1x} + i c_{1y}) (c_{2x} - i c_{2y})| e^{-\sum_j|c_j|^2}\prod_j d\mathrm{Re}(c_j)d\mathrm{Im}(c_j)}{\int_{\mathbb{C}^N}(\sum_j|c_j|^2e^{-\sum_j|c_j|^2}\prod_j d\mathrm{Re}(c_j)d\mathrm{Im}(c_j)}$$

The denominator is the unique symmetric polynomial of the same homogeneity degree as the numerator, thus gives rise to the same radial integration value.

The Gaussian integral in the numerator is decomposable to two types of integrals and in the one in the denominator denominator is a sum of N similar complex gaussian integrals, thus we get:

$$\frac{C}{N} =\frac{\left(\int_{-\infty}^{\infty}\int_{-\infty}^{\infty} \sqrt{c_x^2 + c_y^2} e^{-c_x^2 - c_y^2}\right)^2 \left(\int_{-\infty}^{\infty}\int_{-\infty}^{\infty} e^{-c_x^2 - c_y^2}\right)^{N-2}}{N \left(\int_{-\infty}^{\infty}\int_{-\infty}^{\infty} (c_x^2 + c_y^2) e^{-c_x^2 - c_y^2}\right) \left(\int_{-\infty}^{\infty}\int_{-\infty}^{\infty} e^{-c_x^2 - c_y^2}\right)^{N-1}}= \frac{\pi}{4N}$$

(In the above, the integrals can be solved in polar and cartesian coordinates yielding Gamma functions of intgral or half integral arumements)

Remark:

The equivalence of the integration formulas is based on the following Dirac delta function limit representation:

$$ \delta(y) = \lim_{\epsilon \rightarrow 0} \frac{1}{\sqrt{4\pi \epsilon}} e^{-\frac{-y^2}{4\epsilon}}$$

For a homogeneous function $f(x): \mathbb{R}^{2N} \rightarrow \mathbb{R}$:

$$f(ax) = a^{\alpha} f(x)$$

Denoting: $ r^2 = \sum_{j=1}^{2N} x_j^2$, the ratio of integrals:

$$\frac{\int_{\mathbb{R}^{2N}} f(x) \frac{e^{-\frac{r^2-1}{4\epsilon}}}{\sqrt{4\pi \epsilon}}d\mu_L(\mathbb{R}^{2N})}{\int_{\mathbb{R}^{2N}} r^{\alpha} \frac{e^{-\frac{r^2-1}{4\epsilon}}}{\sqrt{4\pi \epsilon}}d\mu_L(\mathbb{R}^{2N})}$$

is invariant under the scaling $x_j \rightarrow \sqrt{\frac{\epsilon}{\epsilon’}}x_j$, and the integral conserves its form with $\epsilon$ replaced by $\epsilon’$.

Therefore, the ratio of integrals is independent of $\epsilon$. By taking the limit $\epsilon = \frac{1}{4}$, we get our Gaussian integral. By taking the limit $\epsilon \rightarrow 0$, we get a delta function concentrating the measure to the spherical shell:

$$r^2=1$$

in the Euclidean space: $\mathbb{R}^{2N}$, which the integral over the $2N-1$ dimensional sphere.