Eigenvalue densities of sample covariance matrices when the population covariance matrix is a perturbed identity matrix TLDR: I'm looking for a random matrix theory reference for the eigenvalue densities of sample covariance matrices (both dimensions approaching infinity at the same rate) when the true (population) covariance matrix is a perturbed form of the identity matrix.

Consider an $N\times M$ matrix $X$, with elements in each column $X_{.\ j}\sim\mathcal{N}(\mathbf{0},\mathbf{I})$. Random matrix theory tells us that in the limit $N,M\to\infty,\ N/M\to\gamma$, the eigenvalue density of the sample covariance matrix 
$$R=\frac{1}{M}XX^T$$
is given by the Marchenko-Pastur density. 
Now consider an $N\times M$ matrix $Y$, with elements in each column $Y_{.\ j}\sim\mathcal{N}(\mathbf{0},\mathbf{C})$. Let the $N$ eigenvalues of $\mathbf{C}$ (assume diagonal) be $\lambda_i=1+v_i,\ i=\lbrace 1,\ldots,N\rbrace$ where $v_i$ is some zero-mean perturbation, i.e. $\sum_i v_i=0$. So $\mathbf{C}$ can be considered a perturbed identity matrix. As long as the perturbation is "not large", the eigenvalue density of the sample covariance matrix 
$$S=\frac{1}{M}YY^T$$
should be "loosely" described by the same Marchenko-Pastur distribution.
To illustrate, here's a quick MATLAB example for $N=1000$, $N/M=\gamma=1/4$:
X=randn(1e3,4e3);
C=eye(1000)+0.1*diag(randn(1000,1));
A=sqrt(C);
S=(A*X)*(A*X)'/4e3;
[g,z]=hist(eig(S),linspace(0,3,100));
plot(z,g/trapz(z,g),'r--');hold on

mp=@(x)real(sqrt((x-0.25).*(2.25-x))./(0.5*pi*x));
x=linspace(0,3,100);
plot(x,mp(x),'k');hold off

The plot on the left below shows the eigenvalues of $\mathbf{C}$ (black solid line shows eigenvalues of $\mathbf{I}$) and the plot on the right shows the eigenvalue density of $S$ (dashed red) and the asymptotic density of $R$ (solid black).

So you can fiddle around and see that as long as the perturbation is "small" and $\gamma$ is not close to zero, the density will still be given by the Marchenko-Pastur density. 
What are some canonical references that deal with this, so that I can properly reference it in my article?
One paper that kind of deals with a similar issue is 

J. Baik, G. Ben Arous and S. Péché, "Phase transition of the largest eigenvalue for nonnull complex sample covariance matrices", Ann. Probab. vol:33, no:5, 2005

The authors show there that if the identity matrix is "spiked" (i.e., one eigenvalue larger than the others), then depending on the magnitude of spiking, the largest eigenvalue of the sample covariance matrix may or may not fall out of the Marchenko-Pastur bounds. I guess some heuristic argument can be applied to my case here by considering all the perturbations as small positive and negative spikes, but I'm looking for existing literature that deals with it.
 A: You can rewrite 
$$
S=\frac{1}{M}C^{1/2}XX^T(C^{1/2})^T 
$$
where $X_{.j}\sim \mathcal{N}(\boldsymbol 0, \boldsymbol I)$. It is then classical that the limiting eigenvalue distribution of $S$ (in the almost sure weak convergence) is given by $MP_\gamma\boxtimes\nu$, the free multiplicative convolution of the Marchenko-Pastur distribution $MP_\gamma$ of paramter $\gamma$ with the limiting eigenvalue distribution $\nu$ of $C$. Thus, you indeed obtain $MP_\gamma$ if and only if $\nu=\delta_1$, namely if and only if
$$
\frac{1}{N}\sum_{i=1}^N\delta_{v_i} \longrightarrow \delta_0\qquad \mbox{weakly as $N\rightarrow\infty$ }.
$$
For a reference, you can look at the book of Anderson, Guionnet and Zeitouni, "Introduction to random matrices" (Chapter 5), but also any introduction to free probability I guess. 
As you can see, free probability provides a powerful language to describe such perturbed random matrix models!
NB : As an example, a spiked model is the situation where only a fixed finite number of $v_i$'s are non zero, which fit with the above characterization.
A: I am referencing "Spectral Analysis of Large Dimensional Random Matrices" - Bai and Silverstein.
