Closed expression for $\mathbb{E} \left\lbrace \Re \frac{(\textbf{x} + \textbf{y})^{H}\textbf{x}}{\| \textbf{x} + \textbf{y}\|^{2}} \right\rbrace$? Given the random variables $\textbf{x} \sim \mathcal{CN}(\textbf{0}_{M},\sigma_{x}^{2}\textbf{I}_{M})$ and $\textbf{y} \sim \mathcal{CN}(\textbf{0}_{M},\sigma_{y}^{2}\textbf{I}_{M})$ are independent, therefore, as I have calculated,
\begin{equation}\label{eq:lemma7}
\mathbb{E} \left\lbrace \Re \left[ \frac{(\textbf{x} + \textbf{y})^{H}\textbf{x}}{\| \textbf{x} + \textbf{y}\|^{2}} \right] \right\rbrace = \int_{0}^{\infty} \int_{-1}^{1} \frac{(ku + w)f_{U}(u)f_{W}(w)}{ku+\frac{1}{ku}+2w} dwdu,
\end{equation}
where $\mathcal{CN}(.,.)$ is the complex normal random variable, $\Re(.)$ is the real part, $k \triangleq \sqrt{\frac{\sigma_{x}^{2}}{\sigma_{y}^{2}}}$, $f_{W}(w)$ and $f_{U}(u)$ are defined below.
\begin{equation}
f_{W}(w) = \frac{M}{\pi} {B}\left( \frac{1}{2}, M \right) \left( 1 -w^{2} \right)^{M - \frac{1}{2}}, \;\; |w| < 1.
\end{equation}
\begin{equation}\label{eq:lemma6}
f_{U}(u) = \frac{2\Gamma(2M)}{(\Gamma(M))^{2}} \frac{u^{2M-1}}{(u^{2} + 1)^{2M}}, \;\; u > 0.
\end{equation}
where ${B}(.,.)$ is the Beta function and $M > 1$.
I have been using the integral in my calculations, however, I'd like to know if a closed form expression is possible for this double integral. I've tried Mathematica, however, it never finishes running. Follows below the Mathematica code I used.
fw = (M/Pi)*Beta[1/2, M]*((1 - (w^2))^(M - (1/2)));
fu = (2*Gamma[2*M])/((Gamma[M])^2)*(u^(2*M - 1))/(((u^2) + 1)^(2*M));
f = ((a*u + w)*fu*fw)/(a*u + (1/(a*u)) + 2*w);
ii = Integrate[f, {w, -1, 1}, Assumptions -> (a > 0 && M > 2 && 0 < u < \[Infinity])]
Integrate[ii, {u, 0, \[Infinity]}, Assumptions -> (a > 0 && M > 2)]

Follows a link to a matlab Monte Carlo simulation comparison with the proposed closed-form solution.
Matlab script
UPDATE 11/02/2019
By applying the following approximation to the expectation above:
$$\mathbb{E} \left[ \frac{\textbf{w}}{\textbf{z}} \right] \approx \frac{\mathbb{E}[\textbf{w}]}{\mathbb{E}[\textbf{z}]} - \frac{\text{cov}(\textbf{w},\textbf{z})}{\mathbb{E}[\textbf{z}]^{2}} + \frac{\mathbb{E}[\textbf{w}]}{\mathbb{E}[\textbf{z}]^{3}}\text{var}(\mathbb{E}[\textbf{z}]).$$
I found
$$\mathbb{E} \left\lbrace \Re \left[ \frac{(\textbf{x} + \textbf{y})^{H}\textbf{x}}{\| \textbf{x} + \textbf{y}\|^{2}} \right] \right\rbrace \approx \frac{\sigma^2_x}{\sigma_x^2+\sigma_y^2}$$.
 A: My guess in the comments was based on the fact that the expectation depends on $\sigma_x^2/\sigma_y^2$ only, and must be $1/2$ for $\sigma_x^2 = \sigma_y^2$, with limits of $0$ for $\sigma_x^2 \ll \sigma_y^2$ and $1$ for $\sigma_x^2 \gg \sigma_y^2$.
However, I believe this can be shown with an old trick by Marsaglia, although you should make sure there are no problems with the vectors being complex normal.
Write $\Sigma_x = E[x x^H]$ for the covariance matrix.
Let $z = x + y$, and $S = \Sigma_x \, (\Sigma_x + \Sigma_y)^{-1}$. We have $$E[(x - S z) \, z^H] = 0 = E[x - S z] \, E[z]^H \;,$$ so that $x - S z$ and $z$ are uncorrelated and jointly normal, hence independent. Then $$E\left[\frac{z^H x}{\Vert z \Vert^2}\right] = E\left[\frac{z^H (x - S z)}{\Vert z \Vert^2}\right] + E\left[\frac{z^H S z}{\Vert z \Vert^2}\right] \\ = E\left[\frac{z^H}{\Vert z \Vert^2}\right] E\left[x - S z\right] + E\left[\frac{z^H S z}{\Vert z \Vert^2}\right] \\ = E\left[\frac{z^H S z}{\Vert z \Vert^2}\right] \;.$$
Since $S = \frac{\sigma_x^2}{\sigma_x^2 + \sigma_y^2} I$, the result follows.
A: Mathematica evaluates the integral over $w$ in terms of a hypergeometric function,
$$  \int_{-1}^{1} \frac{(ku + w)f_{U}(u)f_{W}(w)}{ku+\frac{1}{ku}+2w} dw=$$
$$\frac{k^2 4^M u^{2 M+1}\Gamma(M+\tfrac{1}{2})}{\sqrt{\pi } \left(k^2 u^2-1\right)^3\left(u^2+1\right)^{2 M}}  \left[\frac{1}{(M-1)!} \left(k^4 u^4+4 M k^2 u^2+4M+3\right)-\frac{4 M}{(M+1)!} \left(k^2 u^2+1\right) (\tfrac{3}{4}+m^2+2m)\, \, _2{F}_1\left(-\tfrac{1}{2},1;M+2;\frac{4 k^2 u^2}{\left(k^2 u^2+1\right)^2}\right)\right].$$
The integral $\int_0^\infty du$ of the first term between square brackets has a closed form expression (again involving a hypergeometric function), but the integral of the second term does not.
In the large-$M$ limit we may average numerator and denominator separately, so
$$\mathbb{E} \left\lbrace \frac{(\textbf{x} + \textbf{y})^{H}\textbf{x}}{\| \textbf{x} + \textbf{y}\|^{2}} \right\rbrace \rightarrow\frac{\sigma^2_x}{\sigma_x^2+\sigma_y^2}, \;\;M\rightarrow\infty.$$
In the accepted answer @student has now elegantly shown this holds actually for all $M$ (quite an impressive "student" !)
