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

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.

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');

• When you write "N modes", do you mean an N-level system? "modes" sounds like fermions or bosons. (Especially with c_i next to it!) ---- Also, rather than doing a log-plot and saying "this is close to C/N", you should do a loglog vs. log plot, to see that there is a line! Commented Aug 15, 2019 at 16:48
• N-level system? qu-N-it, if you wish? $|\psi\rangle\in\mathbb C^d$? -- But modes is implies bosons/fermions. (Not to mention that the "Hilbert space" of N fermionic modes is $2^N$, not $N$.) Commented Aug 16, 2019 at 10:27
• Much better notation now! Commented Aug 17, 2019 at 22:42
• Very stupid question, but why is $\mathbb E[c_i \bar c_j]=\mathbb E[c_i c_j]$ for your case? Commented Oct 17, 2020 at 15:53
• @BlackHat18 tbh in the question I was just stating what my intuition was, I don't know that I had a real proof of the fact. Still, I think you can see it more formally writing down what the expectation value amounts to. You get an integral of the form $\int (c_{ix}+ic_{iy})(c_{j x}\pm ic_{j y}) d c_{ix} d c_{iy}d c_{jx}d c_{jy}$. The integral over $c_{jy}$ is symmetric, being the measure uniform, thus the sign shouldn't matter
– glS
Commented Oct 17, 2020 at 20:15

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.

• this is extremely interesting, but I would love some more detail about the details. Do you have some references to point to that explain how this sort of calculations can be done? Also can you add some more detail about the fact that "Integrals of homogeneous functions over spheres can be traded by Gaussian integrals over the ambient Euclidean space the with proper normalizations"?
– glS
Commented Aug 16, 2019 at 9:09
• @NorbertSchuch $U(N)$ is a subgroup of SO(2N), e.g., through the following inclusion: $$U(N) \ni A+iB \rightarrow \begin{pmatrix} A & B \\ -B & A \end{pmatrix}\in SO(2N)$$ Thus, the probability measure of $2N$ dimensional real random vectors obtained by the action of Haar distributed SO(2N) matrices on a fixed vector, will be automatically invariant under U(N)... Commented Aug 18, 2019 at 6:33
• ...The representation of $S^{2N-1}$ as a homogeneous space of both $SO(2)$ and $U(N)$ can be seen in equations (3.148), (with $N$ to be replaced by $2N$) and (3.149) of Bengtsson and Życzkowski researchgate.net/profile/Karol_Zyczkowski/publication/… Commented Aug 18, 2019 at 6:34
• @gls Yes, the result can be used for any monomial in the moduli of amplitudes. (For more than the second degree in the amplitudes, the Gaussian integral in the denominator is a little bit harder but it can be solved exactly). Please, see the following question in math.stackexchange, where two references are given in the answer. The first reference by Fulton is not open access, but its result is given explicitly in section 5 of arxiv.org/abs/0709.1999v1. ... Commented Aug 21, 2019 at 12:38
• @gls … The second reference is open access, but from a brief reading, it covers the same material as the first with respect to integration over spheres. There is a second technique for integration over spheres: The Pizzetti formula (Please see arxiv.org/abs/1409.8207v1 section 5) which gives the integral of a not-necessarily homogeneous function as an infinite series. The problem with Fulton's result is that it is valid for monomial functions only and can't be used when you have an amplitude modulus. I couldn't find a reference for this exact case. Commented Aug 21, 2019 at 12:38