# 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 * ket)
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! – Norbert Schuch Aug 15 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$.) – Norbert Schuch Aug 16 at 10:27
• Much better notation now! – Norbert Schuch Aug 17 at 22:42

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 Aug 16 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)... – David Bar Moshe Aug 18 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/… – David Bar Moshe Aug 18 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. ... – David Bar Moshe Aug 21 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. – David Bar Moshe Aug 21 at 12:38