I have two related questions concerning normalized complex Gaussian vectors. (Apologies in advance for the length; I wanted to include what related details I already have been able to determine.) All the vectors for this scenario are circularly symmetric and have zero mean, but have a general covariance matrix that is not necessarily a scaled identity matrix nor diagonal. In what follows, $[\cdot]^H$ denotes the complex (Hermitian) transpose, $[\cdot]^{1/2}$ is the matrix square root, and $\mathbf{I}_N$ is the $N{\times}N$ identity matrix.
Let $\mathbf{y} \in \mathbb{C}^N$ be distributed $\sim \mathcal{CN}(0,\mathbf{R}_y)$. Let $\mathbf{z} = \frac{\mathbf{y}}{\|\mathbf{y}\|^2}$. What is the covariance matrix of $\mathbf{z}$, i.e., $\mathbf{R}_z = \mathbb{E}\{\mathbf{zz}^H\} \in \mathbb{C}^{N{\times}N}$?
Let $\mathbf{x} \sim \mathcal{CN}(0,\mathbf{R}_x)$ and $\mathbf{y} \sim \mathcal{CN}(0,\mathbf{R}_y)$, both $\in \mathbb{C}^N$, be correlated with each other. The relation between $\mathbf{x}$ and $\mathbf{y}$ can be expressed as $\mathbf{y} = \mathbf{Ax}+\mathbf{b}$, where $\mathbf{A} \in \mathbb{C}^{N{\times}N}$ is a generic constant complex-valued matrix, and $\mathbf{b} \sim \mathcal{CN}(0,\mathbf{R}_b) \in \mathbb{C}^N$ is an independent complex Gaussian vector. Let $\mathbf{z} = \frac{\mathbf{y}}{\|\mathbf{y}\|^2}$ as before. What are the expectations related to the joint distribution of $\mathbf{x}$ and $\mathbf{z}$, specifically $\mathbb{E}\{\mathbf{x}^H \mathbf{z}\}$, $\mathbb{E}\{\mathbf{xz}^H\}$ (the cross-covariance matrix), and $\mathbb{E}\{\mathbf{x}^H \mathbf{zz}^H\mathbf{x}\}$?
It is the last expectation that I'm most interested in, for both cases when $\mathbf{x}$ and $\mathbf{y}$ are correlated and uncorrelated. In the uncorrelated case, I can at least express $\mathbb{E}\{\mathbf{x}^H \mathbf{zz}^H\mathbf{x}\}$ as $$ \mathbb{E}\{\operatorname{tr}(\mathbf{x}^H \mathbf{zz}^H\mathbf{x})\} = \mathbb{E}\{\operatorname{tr}(\mathbf{xx}^H \mathbf{zz}^H)\} = \operatorname{tr}(\mathbb{E}\{\mathbf{xx}^H \mathbf{zz}^H\}) = \operatorname{tr}(\mathbb{E}\{\mathbf{xx}^H\} \mathbb{E}\{\mathbf{zz}^H\}) = \operatorname{tr}(\mathbf{R}_x \mathbf{R}_z), $$ but that doesn't help a great deal if I don't know $\mathbf{R}_z$.
What I have been able to work out thus far:
The probability distribution for $\mathbf{z}$ for question 1. First convert $\mathbf{y}$ and $\mathbf{z}$ to their real-valued equivalents $\mathbf{u}$ and $\mathbf{v}$, respectively. (For example, $\mathbf{u} = \left[\begin{array}{c} \mathfrak{Re}(\mathbf{y}) \\ \mathfrak{Im}(\mathbf{y}) \end{array}\right]$ is distributed $\sim \mathcal{N}(0,\mathbf{R}_u)$ with length $M=2N$ and covariance $\mathbf{R}_u = \frac{1}{2}\left[\begin{array}{cc} \mathfrak{Re}(\mathbf{R}_y) & -\mathfrak{Im}(\mathbf{R}_y) \\ \mathfrak{Im}(\mathbf{R}_y) & \mathfrak{Re}(\mathbf{R}_y) \end{array}\right]$.) The elements of $\mathbf{u}$ and $\mathbf{v}$ have the forward-backwards relationship equations $$ v_i = \frac{u_i}{\sum\limits_{i=1}^M u_i^2}\, \longleftrightarrow \,\, u_i = \frac{v_i}{\sum\limits_{i=1}^M v_i^2}\quad \textrm{ or }\quad \mathbf{v} = \frac{\mathbf{u}}{\|\mathbf{u}\|^2}\, \longleftrightarrow \,\, \mathbf{u} = \frac{\mathbf{v}}{\|\mathbf{v}\|^2} $$ Taking the Jacobian of this transformation results in $$ |J|=\frac{1}{(v_1^2 + v_2^2 + \ldots + v_M^2)^M} = \frac{1}{\|\mathbf{v}\|^{2M}} $$ We would then have the probability distribution function (pdf) $f_{\mathbf{v}}(\mathbf{v}) = |J|\,f_{\mathbf{u}}\Big(\frac{\mathbf{v}}{\|\mathbf{v}\|^2}\Big)$. Noting $\det(\mathbf{R}_u) = \tfrac{1}{2^M}(\det(\mathbf{R}_y))^2$ and $\|\mathbf{v}\|^2 = \|\mathbf{z}\|^2$, converting from real values back to complex values results in the pdf $$ f_{\mathbf{z}}(\mathbf{z}) = \frac{\exp\left(-\tfrac{\mathbf{z}^H}{\|\mathbf{z}\|^2} \mathbf{R}_y^{-1} \tfrac{\mathbf{z}}{\|\mathbf{z}\|^2}\right)}{\pi^N \|\mathbf{z}\|^{4N} \det(\mathbf{R}_y)} $$
The joint pdf for $\mathbf{x}$ and $\mathbf{z}$ in question 2. The covariance matrices for $\mathbf{x}$ and $\mathbf{y}$ are related by $\mathbf{R}_y = \mathbf{AR}_x\mathbf{A}^H + \mathbf{R}_b$. $f_{\mathbf{x},\mathbf{y}}(\mathbf{x},\mathbf{y})$ can be found from $f(\mathbf{y}|\mathbf{x})\,f(\mathbf{x})$. Conditioned on knowing $\mathbf{x}$, $\mathbf{y}|\mathbf{x}$ then becomes a complex Gaussian vector distributed $\sim \mathcal{CN}(\mathbf{Ax},\mathbf{R}_b)$. Knowing $f(\mathbf{y}|\mathbf{x})$ and $f(\mathbf{x})$, $f_{\mathbf{x},\mathbf{z}}(\mathbf{x},\mathbf{z})$ can be found from $f_{\mathbf{x},\mathbf{y}}(\mathbf{x},\mathbf{y})$ with the same transformation as above. The joint pdf ends up being $$ f_{\mathbf{x},\mathbf{z}}(\mathbf{x},\mathbf{z}) = \frac{\exp\left[-\Big(\tfrac{\mathbf{z}}{\|\mathbf{z}\|^2} - \mathbf{Ax}\Big)^H \mathbf{R}_b^{-1} \Big(\tfrac{\mathbf{z}}{\|\mathbf{z}\|^2} - \mathbf{Ax}\Big) - \mathbf{x}^H\mathbf{R}_x^{-1} \mathbf{x} \right]}{\pi^{2N} \|\mathbf{z}\|^{4N} \det(\mathbf{R}_b \mathbf{R}_x)} $$
The expectation $\mathbb{E}\{\mathbf{z}\}$, which is somewhat trivial. If $\mathbf{y}$ has a mean of zero, then so will $\mathbf{z}$; scaling $\mathbf{y}$ by its squared norm won't change its mean. More formally, $f_{\mathbf{z}}(\mathbf{z})$ is an even function, so $\mathbf{z} f_{\mathbf{z}}(\mathbf{z})$ is an odd function, and thus $\int_{\mathbb{C}^N} \mathbf{z} f_{\mathbf{z}} (\mathbf{z}) d\mathbf{z}$ comes out to zero.
The expectation $\mathbb{E}\{\mathbf{z}^H\mathbf{z}\} = \mathbb{E}\{\|\mathbf{z}\|^2\} = \mathbb{E}\Big\{\frac{\mathbf{y}^H\mathbf{y}}{(\|\mathbf{y}\|^2)(\|\mathbf{y}\|^2)}\Big\} = \mathbb{E}\Big\{\frac{1}{\|\mathbf{y}\|^2}\Big\}$. Obtaining $\mathbf{y}$ with the desired covariance matrix $\mathbf{R}_y$ can be achieved by $\mathbf{y} = \mathbf{R}_y^{1/2} \mathbf{v}$, where $\mathbf{v} \sim \mathcal{CN}(0,\mathbf{I}_N)$. Let the eigen-decomposition of $\mathbf{R}_y$ be $\mathbf{U \Lambda U}^H$, with $\mathbf{\Lambda}$ containing the eigenvalues $e_1, e_2, \dots, e_N$. We then have $$ \|\mathbf{y}\|^2 = \mathbf{y}^H\mathbf{y} = \mathbf{v}^H \mathbf{R}_y \mathbf{v} = \mathbf{v}^H \mathbf{U \Lambda U}^H \mathbf{v} = \mathbf{a}^H \mathbf{\Lambda} \mathbf{a} $$ $\mathbf{a} = \mathbf{U}^H \mathbf{v}$ is complex Gaussian with a covariance matrix of $\mathbf{U}^H \mathbf{I}_N \mathbf{U} = \mathbf{U}^H \mathbf{U} = \mathbf{I}_N$. We therefore have that $\|\mathbf{y}\|^2 = \sum_{i=1}^N e_i |a_i|^2$, with the entries $a_i$ being independent complex Gaussian scalars with zero mean and unit variance. Hence, $|a_i|^2$ follows an exponential distribution with a unit rate parameter, and $e_i |a_i|^2$ is exponentially distributed with rate parameter $\lambda_i = 1/e_i$. It therefore follows that $\|\mathbf{y}\|^2$ has a hypoexponential distribution with parameters $\{\lambda_1, \lambda_2, \dots, \lambda_N\} = \{\tfrac{1}{e_1}, \tfrac{1}{e_2}, \dots, \tfrac{1}{e_N}\}$. I am primarily interested in the case where there are no repeated eigenvalues, in which case the pdf of $S = \|\mathbf{y}\|^2$ can be expressed as $f_S(s) = \sum\limits_{i=1}^N \lambda_i e^{-\lambda_i s} \Big(\prod\limits_{j=1, j \neq i}^N \tfrac{\lambda_j}{\lambda_j - \lambda_i}\Big)$. From this, by a standard random variable transformation, we get the pdf for $W = 1/S$ is $f_W(w) = \tfrac{f_S(1/w)}{w^2} = \tfrac{1}{w^2} \sum\limits_{i=1}^N \lambda_i e^{-\lambda_i/w} \Big(\prod\limits_{j=1, j \neq i}^N \tfrac{\lambda_j}{\lambda_j - \lambda_i}\Big)$. We then ultimately arrive at the result that $\mathbb{E}\{W\} = \mathbb{E}\Big\{\frac{1}{\|\mathbf{y}\|^2}\Big\} = \sum\limits_{i=1}^N \lambda_i \log(\lambda_i) \Big(\prod\limits_{j=1, j \neq i}^N \tfrac{\lambda_j}{\lambda_i - \lambda_j}\Big)$. (Note the reversal of order of $i$ and $j$ in the denominator compared to the pdf.) Alternatively, to put it directly in terms of the eigenvalues, substituting in $\lambda_i = 1/e_i$ and rearranging, we get $\mathbb{E}\Big\{\frac{1}{\|\mathbf{y}\|^2}\Big\} = - \sum\limits_{i=1}^N e_i^{N-2} \log(e_i) \Big(\prod\limits_{j=1, j \neq i}^N \tfrac{1}{e_j - e_i}\Big)$. (Now with the same order of $j$ and $i$ in the denominator as in the pdf.)
Corollary: If $\mathbf{R}_y = c\mathbf{I}_N$, $\|\mathbf{y}\|^2$ is a sum of exponentially distributed variables all with the same rate parameter $\lambda = 1/c$, or equivalently a sum of Gamma distributed variables, each $\sim \Gamma(1,c)$ with shape parameter 1 and scale parameter $c$. $\|\mathbf{y}\|^2$ is therefore distributed $\sim \Gamma(N,c)$, and $\tfrac{1}{\|\mathbf{y}\|^2}$ has an inverse Gamma distribution, with a mean of $\tfrac{1}{c(N-1)}$.Partially related: When $\mathbf{x}$ and $\mathbf{y}$ are independent with scaled identity covariance matrices, this post and this other post appear to cover the distributions and expectations of $\mathbf{x}^H \mathbf{y}$ and $\mathbf{x}^H\mathbf{yy}^H\mathbf{x}$, respectively, though without the desired normalization of one of the vectors.
I also know that $\frac{\mathbf{y}}{\|\mathbf{y}\|}$, without the square in the denominator, has what is known as a "projected normal distribution", with the vector points falling on the $2N$-dimensional unit hypersphere. However, results in that area don't seem to be of much help, largely because they are generally concerned with the distribution of the angle(s) of the vectors, and mostly are interested in at most 3-dimensional vectors. This post appears to cover the covariance matrix of $\frac{\mathbf{y}}{\|\mathbf{y}\|}$ for the case when $\mathbf{R}_y = \mathbf{I}_N$.
I am stuck at trying to obtain the expectations at this point. I have tried integrating the pdf (both real and complex Gaussian versions), but it quickly falls apart after getting about one or two integrals deep into the $2N$-fold multiple integral. I can't find anything useful in a table of integrals (Gradshteyn and Ryzhik), and Mathematica also gives up on it at that point. For the same reason, I'm also unable to derive a moment generating function for the pdf. Contour integration also hasn't helped, what little of it that I know. You end up with a square root of a sum of variables in the denominator after one integral, so the function no longer has simple poles and/or is no longer meromorphic. Thus, the method of residues no longer appears to work. (My apologies if I'm not using the terminology here properly; my knowledge of contour integration is extremely limited.) I've also tried working with the joint pdfs of the entries in the vector starting from the hypoexponential or gamma distribution functions, but again haven't gotten very far. (Though I did get derive a closed-form side result for $\int_{-\sqrt{S}}^{\sqrt{S}} e^{-\beta x^2} \operatorname{erf}(\sqrt{\beta(S-x^2)}) dx$, which doesn't seem to be in the tables.)
There seems to be some relation to the Beta and Dirichlet distributions, in that both result from calculating (random variable)$\div$(sum of random variables). In those cases, the summed random variables are i.i.d. and it results that the sum is statistically independent of its components. It's possible something similar could happen here, but notably the constituent variables are not identically distributed nor independent. However, perhaps they could be broken up at least into independent variables like in the eigen-decomposition trick to calculate $\mathbb{E}\{\|\mathbf{z}\|^2\}$. I haven't gotten far with this tactic either. It doesn't seem like the sum is statistically independent anymore, though I haven't been able to verify it either way.
I would be content to do the integration numerically with specific covariance matrices for $\mathbf{y}$, but Mathematica has trouble with that as well. For $f_{\mathbf{z}}(\mathbf{z})$, it takes a rather long time for even small dimensions, and complains that the global error increases too many times. (I suspect the reason stems from the singularity in the pdf at the origin of the multi-dimensional axes.) However, it does at least eventually spit out something reasonable; the resulting value matches what is obtained if you take an average over numerous realizations of $\mathbf{y}$. Unfortunately, it can't do even a numerical integration for $f_{\mathbf{x},\mathbf{z}}(\mathbf{x},\mathbf{z})$ properly. The end result is nonsensical or worthless, with something like "value ± error" and the error being larger than the value.
So anyway, that's where things stand. I would appreciate any pointers you may have in what direction to take to get to expressions for the expectations. It would not surprise me if they have to be expressed in terms of some higher-level expressions (e.g. erf, Marcum Q, Bessel functions, hypergeometric, etc.) I just hope it's something that some mathematical software like Matlab or Mathematica can calculate with greater ease, compared to the current trouble the latter is having with the integrals. In a pinch, I could simply take an average over numerous vector realizations to get the expectations, though I'd rather have mathematical expressions for the expectations instead.