Distribution of dot product of two unit random vectors Consider $\mathbf{u}, \mathbf{v}\in \mathcal{C}^M$ to be two independent unit norm random vectors on the $M-1$ dimensional complex sphere $\mathcal{S}^{M-1}$. In addition, $\mathbf{u}$ follows an isotropic distribution (i.e., $\mathbf{u}$ is uniformly distributed on the complex sphere $\mathcal{S}^{M-1}$. What is the distribution of $Z=|\mathbf{u}\cdot\mathbf{v}|^2$?
I know that if $\mathbf{v}$ is also uniformly distributed on the complex sphere $\mathcal{S}^{M-1}$, then $Z$ follows Beta$(1, M-2)$ distribution. (I don't know how to prove this!) Does the same result hold if $\mathbf{v}$ follows an arbitrary distribution?
 A: If $u$ is uniformly distributed over the sphere, we can write it as $u=Uv$, where $U$ is a unitary transformation uniformly distributed over the unitary group. Then the quantity $|u\cdot v|^2$ is just the modulus square of a matrix element of $U$. So your question is: what is the distribution of the modulus square of a matrix element of a random unitary matrix?
Consider a column of the matrix as a random vector. The only constraint is must satisfy is that its norm must be unit, $\sum_{i=1}^N|z_i|^2=1$. In other words, the joint distribution of the entries is simply 
$$\delta(1-\sum_{i=1}^N|z_i|^2)=\int ds e^{is(1-\sum_{i=1}^N|z_i|^2)}$$
The distribution of a single element, $z_1$ say, is obtained integrating over the remaining ones. Indeed, you get something proportional to $(1-|z_1|^2)^{N-2}$ in dimension $N$, which is your Beta distribution.
This is independent of $v$, so the distribution of $v$ is irrelevant.
A very similar argument applies to the real case and random orthogonal matrices.
This is discussed in several articles, such as
Kus, Mostowski and Haake, J. Phys. A 21, L1073 (1988)
Haake and Zyczkowski, Phys. Rev. A 42, 1013 (1990)
K. Zyczkowski, H.-J. Sommers, J. Phys. A 33, 2045 (2000) 
A: $\newcommand{\bu}{\boldsymbol{u}}$ $\newcommand{\bv}{\boldsymbol{v}}$ $\newcommand{\bP}{\mathbb{P}}$ $\newcommand{\bR}{\mathbb{R}}$
Here are the details in Andrea Blass' comment. We define
$$F:S^{m-1}\times S^{m-1} \to \bR, \;\;F(\bu,\bv):=|\bu\cdot\bv|^2. $$ Note that the range of $F$ is $[0,1]$. Denote by $p(d\bv)$ the probability distribution of $\bv$. For every interval  $[a,b]\subset (0,1)$ we have
$$ \bP[a\leq f\leq b] =\int_{S^{m-1}} \bP[a\leq F\leq b| \bv=\bv_0] p(d\bv_0) $$
(use the independence of $\bu$ and $\bv$)
$$= \int_{S^{m-1}} \Bigl(\;{\rm Area}\;\{ \sqrt{a}\leq |\bu\cdot\bv_0|\leq \sqrt{b}\}\;\Bigr)\; p(d \bv_0) $$
$$=2 \int_{S^{m-1}} \underbrace{ \Bigl(\;{\rm Area}\;\{ \sqrt{a}\leq \bu\cdot\bv_0\leq \sqrt{b}\}\;\Bigr)}_{=:I(a,b,\bv_0)}\; p(d \bv_0) $$
Due to the rotational symmetry, the integrand $I(a,b,\bv_0)$  is independent of $\bv_0$ so  I will denote it by $I(a,b)$.   Hence
$$\bP[a\leq f\leq b]= 2 I(a,b). $$
To compute $I(a,b)$  use   the coarea formula exactly as in Example 9.1.10 of these notes.  We have $\newcommand{\bsi}{\boldsymbol{\sigma}}$
$$ I(a,b) = \bsi_{m-2}\int_{\sqrt{a}}^{\sqrt{b}}(1-t^2)^{\frac{m-3}{2}} dt, $$
where $\bsi_k$ denotes the  area of the unit $k$-dimensional sphere
$$\bsi_k=\frac{2\pi^{\frac{k+1}{2}}}{\Gamma\left(\frac{k+1}{2}\right)}. $$
If $\rho_F(x)$ denotes the probability density of $F$, then we deduce that
$$\rho_F(x)=2\frac{d}{dh}\Bigl|_{h=0} I(x,x+h)=\bsi_{m-2} x^{-1/2}(1-x^2)^{\frac{m-3}{2}}. $$
Thus the answer is independent of the distribution   of $\bv$.
