Joint convergence of sample cross-covariance matrices

Suppose $X$ and $Y$ are $p\times n$ matrices each with mean $0$, variance $1$ and where all moments of any order are uniformly bounded as $p,n \to \infty$. Moreover, the entries $(X_{ij}, Y_{ij})$ are independent across $i,j$ with a common correlation $\rho$. Let $C=n^{-1}XY^*$ be the sample cross-covariance matrix. We show that if $n, p\to \infty, p/n\to y\neq 0$, then $C$ converges in the algebraic sense and the limit moments depend only on $\rho$. Independent copies of such matrices with same $p$ but different $n$, say $\{n_l\}$, different correlations $\{\rho_l\}$, and different non-zero $y$'s, say $\{y_l\}$ also converge jointly and are asymptotically free. When $y=0$, the matrix $\sqrt{np^{-1}}(C-\rho I_p)$ converges to an elliptic variable with parameter $\rho^2$. In particular, this elliptic variable is circular when $\rho=0$ and is semi-circular when $\rho=1$. If we take independent $C_l$, then the matrices $\{\sqrt{n_lp^{-1}}(C_l-\rho_l I_p)\}$ converge jointly and are also asymptotically free. As a consequence, the limiting spectral distribution of any symmetric matrix polynomial exists and has compact support.