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
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.