---
abstract: |
    We prove the first explicit rate of convergence to the Tracy-Widom distribution for the fluctuation of the largest eigenvalue of sample covariance matrices that are not integrable. Our primary focus is matrices of type $ X^*X $ and the proof follows the Erdös-Schlein-Yau dynamical method. We use a recent approach to the analysis of the Dyson Brownian motion from [@bourgade2018extreme] to obtain a quantitative error estimate for the local relaxation flow at the edge. Together with a quantitative version of the Green function comparison theorem, this gives the rate of convergence.

    Combined with a result of Lee-Schnelli [@lee2016tracy], some quantitative estimates also hold for more general separable sample covariance matrices $ X^* \Sigma X $ with general diagonal population $ \Sigma $.
address: 'Courant Institute of Mathematical Sciences, New York University, 251 Mercer St., New York, NY 10012'
author:
- Haoyu Wang
bibliography:
- 'CovarianceTW.bib'
title: Quantitative Universality for the Largest Eigenvalue of Sample Covariance Matrices
---

Introduction
============

Overview and main results
-------------------------

Edge universality of sample covariance matrices has been a classical problem in random matrix theory. It is well known that the distribution of the largest eigenvalue (after appropriate rescaling) converges to the Tracy-Widom distribution. Early non-quantitative results were first proved in [@peche2009universality; @pillai2014universality; @Soshnikov]. Quantitative estimates, however, were only obtained for the Wishart ensemble (see [@el2006rate; @Johansson; @Johnstone; @ma2012accuracy]), which is essentially an integrable model. In this paper, we prove the explicit rate of convergence $ N^{-2/9} $ to the Tracy-Widom distribution for all sample covariance matrices of type $ X^*X $ with general distributed entries, and an analogous result with deterministic rate $ N^{-1/57} $ for all separable sample covariance matrices $ X^* \Sigma X $ with diagonal population $ \Sigma $. For simplicity and motivations from statistics, we only consider the real case, but the whole proof and the results are also true for complex sample covariance matrices.

Let $ X=(x_{ij}) $ be an $ M \times N $ data matrix with independent real valued entries with mean 0 and variance $ M^{-1} $, $$\label{e.Assumption1}
x_{ij} = M^{-1/2}q_{ij},\ \ \ {\mathbb{E}}q_{ij}=0,\ \ \ {\mathbb{E}}q_{ij}^2 =1.$$ Furthermore, we assume the entries $ q_{ij} $ have a sub-exponential decay, that is, there exists a constant $ \theta>0 $ such that for $ u>1 $, $$\label{e.Assumption2}
{\mathbb{P}}(|q_{ij}| > u) {\leqslant}\theta^{-1} \exp (-u^\theta).$$ This sub-exponential decay assumption is mainly for convenience, other conditions such as the finiteness of a sufficiently high moment would be enough. (For a necessary and sufficient condition for the edge universality we refer to [@ding2018necessary].)

The sample covariance matrix corresponding to data matrix $ X $ is defined by $ H := X^* X $. Throughout this paper, to avoid trivial eigenvalues, we will be working in the regime $$\xi=\xi(N) := N/M,\ \ \ \lim_{N \to \infty} \xi \in (0,1) \  \mbox{or}\ \xi \equiv 1.$$ We will mainly work with the rectangular case $ 0< \xi <1 $, but will also show how to adapt the arguments to the square case $ M \equiv N $. (The reason why we do not discuss the general case $ \lim \xi =1 $ is merely technical due to the lack of local laws at the hard edge. In particular, the rigidity estimate at the hard edge is only known for a fixed $ \xi \equiv 1 $ but not for $ \xi = \xi(N) \to 1 $.)

We order the eigenvalues of $ H $ as $ \lambda_1 {\leqslant}\cdots {\leqslant}\lambda_N $, and use $ \lambda_+ $ to denote the typical location of the largest eigenvalue (see for the definition). For the main result of this paper, we consider the Kolmogorov distance $${\mathsf{d_K}}(X,Y) := \sup_{x} \left| {\mathbb{P}}(X {\leqslant}x) - {\mathbb{P}}(Y {\leqslant}x) \right|.$$

\[t.Rate\] Let $ H_N $ be sample covariance matrices satisfying and . Let $ {\mathsf{TW}}$ be the Tracy-Widom distribution. For any $ {\varepsilon}>0 $, for large enough $ N $ we have $${\mathsf{d_K}}(N^{2/3}(\lambda_N - \lambda_+),{\mathsf{TW}}) {\leqslant}N^{-\frac{2}{9}+{\varepsilon}}.$$

The null case $ X^*X $ is our primary concern in this paper, but quantitative estimates are also valid for general diagonal population matrices $ X^* \Sigma X $ thanks to the comparison theorem for the Green function flow by Lee and Schnelli (see Section \[s.GeneralPopulation\] for more details). Combining our quantitative edge universality for the null case (Theorem \[t.Rate\]) with the Green function comparison by Lee-Schnelli (Proposition \[p.GreenFunctionFlow\]), we derive the rate of convergence to Tracy-Widom distribution for separable sample covariance matrices with general diagonal population.

\[t.RateGeneralPopulation\] Let $ Q := X^* \Sigma X $ be an $ N \times N $ separable sample covariance matrix, where $ X $ is an $ M \times N $ real random matrix satisfying and , and $ \Sigma $ is a real diagonal $ M \times M $ matrix satisfying . Let $ \mu_N $ be the largest eigenvalue of $ Q $. For any $ {\varepsilon}>0 $, for large enough $ N $ we have $$\label{e.RateGeneralPopulation}
{\mathsf{d_K}}\left( \gamma_0 N^{2/3} (\mu_N - E_+),{\mathsf{TW}}\right) {\leqslant}N^{-\frac{1}{57} + {\varepsilon}},$$ where $ E_+ $ defined in denotes the rightmost endpoint of the spectrum and $ \gamma_0 $ is a normalization constant defined in .

The method of the paper follows the three-step strategy of the Erdös-Schlein-Yau dynamical approach [@erdHos2011universality]: (i) a priori bounds on locations of eigenvalues; (ii) local relaxation of the eigenvalue dynamics; (iii) a density argument showing eigenvalues statistics have not changed after short time.

Specifically, in this paper, the three-step strategy is employed in the following way: (i) is the rigidity for singular values, which can be rephrased from classical results on eigenvalues (see [@bloemendal2014isotropic; @pillai2014universality]). For the particular square case $ M \equiv N $, we use a different rigidity estimate at the hard edge from [@alt2017local],which results in a slightly different proof; (ii) is the recent approach to the analysis of Dyson Brownian motion from [@bourgade2018extreme], which introduced an observable defined via interpolation with integrable models (see [@bourgade2016fixed; @landon2019fixed]). It describes the singular values evolution through a stochastic advection equation; (iii) is a quantitative version of the Green function comparison theorem, which is a slight extension of the classical result from [@GreenFunctionComparison].

As discussed in [@pillai2014universality Section 3], though we discuss the problems in the context of covariance matrices, the proof should also work for more generalized problems such as the quantitative edge universality of correlation matrices (see [@pillai2012edge]).

This paper is organized as follows. The main part of the paper (Section \[s.prelim\]-\[s.Comparison\]) is devoted to the null case $ X^*X $. In Section \[s.prelim\] we rephrase the classical results on the eigenvalues of sample covariance matrices to the version for singular values, including Dyson Brownian motion, local laws and rigidity estimates. In Section \[s.LocalRelaxationFlow\], we define an observable that describes the evolution of singular values, and then prove the error estimate for the local relaxation flow at the edge by studying the dynamics of the observable. In Section \[s.Comparison\] we prove the quantitative Green function comparison theorem and use it to derive the rate of convergence to the Tracy-Widom distribution. Finally, in Section \[s.GeneralPopulation\] we generalize our result to separable sample covariance matrices with diagonal population by studying the interpolation between general covariance matrices with the null case.

Notations
---------

Throughout this paper, we use the notation $ A \lesssim B $ if there exists a constant $ C $ which is independent of $ N $ such that $ A {\leqslant}CB $ holds. We also denote $ A \sim B $ if both $ A \lesssim B $ and $ B \lesssim A $ hold. If $ A $ and $ B $ are complex valued, $ A \sim B $ means $ {\text{\rm Re}\hspace{0.1cm}}A \sim {\text{\rm Re}\hspace{0.1cm}}B $ and $ {\text{\rm Im}\hspace{0.1cm}}A \sim {\text{\rm Im}\hspace{0.1cm}}B $. We also denote $ C $ a generic constant which does not depend on $ N $ but may vary form line to line. We use $ \llbracket A,B \rrbracket := [A,B] \cap {\mathbb{Z}}$ to denote the set of integers between $ A $ and $ B $.

We also denote $$\varphi = e^{C_0 (\log \log N)^2}$$ a subpolynomial error parameter, for some fixed $ C_0>0 $. This constant $ C_0 $ is chosen large enough so that the eigenvalues (and singular values) rigidity and the strong local Marchenko-Pastur law hold (see section \[s.LocalLaw\]).

Acknowledgement {#acknowledgement .unnumbered}
---------------

The author would like to thank Prof. Paul Bourgade for suggesting this problem, helpful discussions, and useful comments on the early draft of the paper.

Preliminaries {#s.prelim}
=============

Dyson Brownian motion for covariance matrices
---------------------------------------------

Let $ B $ be an $ M \times N $ real matrix Brownian motion: $ B_{ij} $ are independent standard Brownian motions. We define the $ M \times N $ matrix $ M_t $ by $$M_t = M_0 + \dfrac{1}{\sqrt{N}}B_t.$$ The eigenvalues dynamics for the real Wishart process $ X_t := M_t^* M_t $ was first proved in [@bru1989diffusions]. Under our normalization convention, the equation is in the following form given in [@bourgade2017eigenvector Appendix C] $$\label{e.EigenvalueDBM}
{\mathrm{d}}\lambda_k = 2\sqrt{\lambda_k}\dfrac{{\mathrm{d}}B_{kk}}{\sqrt{N}} + \left( \dfrac{M}{N} + \dfrac{1}{N} \sum_{l \neq k}\dfrac{\lambda_k + \lambda_l}{\lambda_k - \lambda_l} \right){\mathrm{d}}t.$$ Due to technical issues, it is difficult to use the coupling method from [@bourgade2016fixed] to analyze in a direct way. This motivates us to consider the singular values instead.

Let $ s_k := \sqrt{\lambda_k} $ denote the singular values of $ X $. The Dyson Brownian motion for singular values dynamics of such sample covariance matrices is the following Ornstein-Uhlenbeck process [@erdHos2012local equation (5.8)]. $${\mathrm{d}}s_k = \dfrac{{\mathrm{d}}B_k}{\sqrt{N}} + \left[ -\dfrac{1}{2\xi}s_k + \dfrac{1}{2}\left( \dfrac{1}{\xi} -1 \right) \dfrac{1}{s_k} + \dfrac{1}{2N} \sum_{l \neq k} \left( \dfrac{1}{s_k - s_l} + \dfrac{1}{s_k + s_l} \right) \right]{\mathrm{d}}t, \ \ \ 1 {\leqslant}k {\leqslant}N.$$ An important idea in this paper is the following symmetrization trick (see [@che2019universality equation (3.9)]): $$s_{-i}(t)=-s_i(t),\ \ \ B_{-i}(t)=-B_i(t),\ \ \forall t {\geqslant}0,\ 1 {\leqslant}i {\leqslant}N.$$ From now we label the indices from $ -1 $ to $ -N $ and $ 1 $ to $ N $, so that the zero index is omitted. Unless otherwise stated, this will be the convention and we will not emphasize it explicitly. After symmetrization, the dynamics turns to the following form $$\label{e.SymmetrizedDynamics}
{\mathrm{d}}s_k = \dfrac{{\mathrm{d}}B_k}{\sqrt{N}} + \left[ -\dfrac{1}{2\xi}s_k + \dfrac{1}{2}\left( \dfrac{1}{\xi} -1 \right) \dfrac{1}{s_k} + \dfrac{1}{2N} \sum_{l \neq \pm k} \dfrac{1}{s_k - s_l} \right]{\mathrm{d}}t, \ \ \ -N {\leqslant}k {\leqslant}N, k \neq 0.$$

Local law and rigidity for singular values {#s.LocalLaw}
------------------------------------------

The local law and rigidity estimates are classical results for the eigenvalues of sample covariance matrices. In this section, for later use, we rephrase these results into the corresponding version in terms of singular values.

It is well known that the empirical measure of the eigenvalues converges to the Marchenko-Pastur distribution $$\rho_{{\mathsf{MP}}}(x) = \dfrac{1}{2\pi \xi}\sqrt{\dfrac{[(x-\lambda_-)(\lambda_+ -x)]_+}{x^2}},$$ where $$\label{e.EndPoints}
\lambda_{\pm} = (1 \pm \sqrt{\xi})^2.$$ Define the typical locations of the singular values: $$\gamma_k :=  \inf \left\{ E>0 : \int_{-\infty}^{E^2} \rho_{{\mathsf{MP}}}(x){\mathrm{d}}x {\geqslant}\dfrac{k}{N} \right\},\ \ \ 1 {\leqslant}k {\leqslant}N.$$ Following the symmetrization trick, we also define $ \gamma_{-k} = -\gamma_k $. By a change of variable, it is easy to check that $$\label{e.typical}
\int_{-\infty}^{\gamma_k} \rho(x){\mathrm{d}}x=\dfrac{N+k}{2N}, \ \ \ \int_{-\infty}^{\gamma_{-k}} \rho(x){\mathrm{d}}x = \dfrac{N-k}{2N},$$ where $ \rho(x) $ is the counterpart of Marchenko-Pastur law for singular values, defined by $$\label{e.measure}
\rho(x)=\dfrac{1}{2\pi \xi}\sqrt{\dfrac{[(x^2-\lambda_-)(\lambda_+ - x^2)]_+}{x^2}},\ \ \ \sqrt{\lambda_-} {\leqslant}|x| {\leqslant}\sqrt{\lambda_+}.$$

Denote $ s_1 {\leqslant}\cdots {\leqslant}s_N $ the singular values of the data matrix $ X $, and extend the singular values following the symmetrization trick by $ s_{-k}=-s_k $. For $ z=E+{\text{\rm i}}\eta \in \mathbb{C} $ with $ \eta>0 $, let $ m_N(z) $ and $ S_N(z) $ denote the Stieltjes transform of the empirical measure of the (symmetrized) singular values and eigenvalues, respectively: $$m_N(z):= \dfrac{1}{2N}\sum_{-N {\leqslant}k {\leqslant}N} \dfrac{1}{s_k -z},\ \ \ S_N(z) := \dfrac{1}{N}\sum_{k=1}^N \dfrac{1}{\lambda_k -z}.$$ As mentioned previously, in the summation from $ -N $ to $ N $ the $ 0 $ index is always excluded. Note that due to the symmetrization, this is equivalent to $$\label{e.DiscreteStieltjes}
 m_N(z) = \dfrac{1}{2N}\sum_{k=1}^N \left( \dfrac{1}{s_k -z} + \dfrac{1}{-s_k - z} \right) = \dfrac{1}{N} \sum_{k=1}^N \dfrac{z}{s_k^2 - z^2} = z S_N(z^2).$$ On the other hand, use $ m_{{\mathsf{MP}}}(z) $ to denote the Stieltjes transform of the Marchenko-Pastur law $$m_{{\mathsf{MP}}}(z) := \int_{{\mathbb{R}}} \dfrac{\rho_{{\mathsf{MP}}}(x)}{x-z}{\mathrm{d}}x=\dfrac{1-\xi-z+\sqrt{(z-\lambda_-)(z-\lambda_+)}}{2\xi z},$$ where $ \sqrt{\quad} $ denotes the square root on the complex plane whose branch cut is the negative real line. With this choice we always have $ {\text{\rm Im}\hspace{0.1cm}}m_{{\mathsf{MP}}}(z)>0 $ when $ {\text{\rm Im}\hspace{0.1cm}}z>0 $. For the singular values, recall the limit distribution $ \rho(x) $ for the empirical measure and use $ m(z) $ to denote its corresponding Stieltjes transform $$m(z) := \int_{{\mathbb{R}}} \dfrac{\rho(x)}{x-z}{\mathrm{d}}x = \int_{{\mathbb{R}}} \dfrac{1}{x-z}\dfrac{1}{2\pi \xi}\sqrt{\dfrac{[(x^2-\lambda_-)(\lambda_+ - x^2)]_+}{x^2}}{\mathrm{d}}x.$$ We have the following relation between $ m(z) $ and $ m_{{\mathsf{MP}}}(z) $ $$\begin{gathered}
\label{e.ContinuousStieltjes}
m(z) = \int_{\sqrt{\lambda_-}}^{\sqrt{\lambda_+}} \left( \dfrac{1}{x-z} - \dfrac{1}{x+z} \right) \dfrac{1}{2\pi \xi}\sqrt{\dfrac{[(x^2-\lambda_-)(\lambda_+ - x^2)]_+}{x^2}}{\mathrm{d}}x\\
= \int_{{\mathbb{R}}}\dfrac{z}{x-z^2}\rho_{{\mathsf{MP}}}(x){\mathrm{d}}x = zm_{{\mathsf{MP}}}(z^2).\end{gathered}$$ It is well known that we have the strong local Marchenko-Pastur law [@bloemendal2014isotropic; @pillai2014universality] for the estimate of $ S_N(z) $, i.e. for any $ D>0 $, there exists $ N_0(D)>0 $ such that for every $ N {\geqslant}N_0 $ we have $${\mathbb{P}}\left( \left| S_N(z) - m_{{\mathsf{MP}}}(z) \right| {\leqslant}\dfrac{\varphi}{N\eta} \right) > 1-N^{-D}.$$ By the relations and , we know that $$|m_N(z) - m(z)| = \left| z \left( S_N(z^2) - m_{{\mathsf{MP}}}(z^2) \right) \right| {\leqslant}|z| \left| S_N(z^2) - m_{{\mathsf{MP}}}(z^2) \right|.$$ Combining with the strong local Marchenko-Pastur law, this gives us $$\label{e.LocalLaw}
{\mathbb{P}}\left( \left| m_N(z) - m(z) \right| \lesssim \dfrac{\varphi}{N\eta} \right) > 1-N^{-D}.$$

![Graphs of the density $ \rho(x) $ for $ x{\geqslant}0 $ (i.e. the density for the actual singular values) with $ \xi=1,(3/5)^2,(3/10)^2 $ respectively. The curve for $ \xi=1 $ does not have a “square root”-type shape at the edge $ \sqrt{\lambda_-} =0 $ in this case due to the singularity of the Marchenko-Pastur distribution.](density.jpg){height="4cm"}

For the rigidity estimates, a key observation is that the critical case $ \xi=1 $ is significantly different from other cases. This is because the Marchenko-Pastur law $ \rho_{{\mathsf{MP}}} $ has a singularity at the point $ x=0 $ in this situation. When $ \xi < 1 $, the rigidity of singular values can be easily obtained from the analogous estimates for eigenvalues (see [@pillai2014universality]). Let $ {\widehat}{k}:=\min(k,N+1-k) $, for any $ D>0 $ there exists $ N_0(D) $ such that the following holds for any $ N {\geqslant}N_0 $, $$\label{e.RigidityNotOne}
{\mathbb{P}}\left(|s_k - \gamma_k| {\leqslant}\varphi^{\frac{1}{2}} N^{-\frac{2}{3}}({\widehat}{k})^{-\frac{1}{3}} \  \mbox{for all}\ k \in \llbracket 1,N \rrbracket \right) > 1-N^{-D}.$$

For the critical case $ \xi=1 $, now the Marchenko-Pastur distribution is supported on $ [0,4] $ and is given by $ \rho_{{\mathsf{MP}}}(x)= \frac{1}{2\pi}\sqrt{(4-x)/x} $. A key observation is that the scales of eigenvalue spacings are different at the two edges. Due to this phenomenon, we use the following two different results, depending on the location in the spectrum.

On the one hand, the Marchenko-Pastur distribution still behaves like a square root near the soft edge $ x=4 $, which implies that the result is the same as the rectangular case. The rigidity estimate near the soft edge can be easily adapted from the result for eigenvalues in [@bloemendal2014isotropic Theorem 2.10], i.e. for some (small) $ \omega>0 $ and any $ {\varepsilon}>0 $ we have $$\label{e.RigidityOneSoft}
{\mathbb{P}}\left( |s_k - \gamma_k| {\leqslant}N^{\varepsilon}(N-k+1)^{-\frac{1}{3}}N^{-\frac{2}{3}} \ \mbox{for all} \ k \in \llbracket (1-\omega)N,N \rrbracket \right) > 1-N^{-D}.$$

On the other hand, as explained in [@cacciapuoti2013local], at the hard edge $ x=0 $ the typical distance between eigenvalues and the edge is of order $ N^{-2} $, which is much smaller than the typical distance between neighbouring eigenvalues in the bulk (or at the soft edge). Note that in this situation, the measure for the symmetrized singular values coincides with the standard semicircle law, that is $ \rho(x) = \frac{1}{2\pi}\sqrt{(4-x^2)_+} $. By the relation , this means that the typical $ k $-th singular value $ \gamma_k $ of an $ N \times N $ data matrix shares the same position with the typical $ (N+k) $-th eigenvalue of a $ 2N \times 2N $ generalized Wigner matrix. The link between these two models can be illustrated by the symmetrization trick: Define the $ 2N \times 2N $ matrix $${\widetilde}{H} = \left(
\begin{matrix}
0 & X^*\\
X & 0
\end{matrix}
\right),$$ then we know that the eigenvalues of $ {\widetilde}{H} $ are precisely the symmetrized singular values of $ X^* $. Note that we have $ {\widetilde}{H}={\widetilde}{H}^* $, $ {\mathbb{E}}{\widetilde}{H}_{ij}=0 $ and $ \sum_{i=1}^{2N} {\mathbb{E}}{\widetilde}{H}_{ij}^2 = 1 $ for every $ j \in \llbracket 1,2N \rrbracket $. This shows that $ {\widetilde}{H} $ is indeed a Wigner-type matrix except the lack of nondegeneracy condition caused by the zero blocks. By considering the matrix of this type, the rigidity at the hard edge can be proved directly from [@alt2017local Theorem 2.7] $$\label{e.RigidityOneHard}
{\mathbb{P}}\left( |s_k - \gamma_k| {\leqslant}N^{-1+{\varepsilon}} \ \mbox{for all}\ k \in \llbracket 1,(1-\omega)N \rrbracket \right) > 1-N^{-D}.$$

Stochastic Advection Equation for Singular Values Dynamics {#s.LocalRelaxationFlow}
==========================================================

Stochastic advection equation
-----------------------------

We follow the comparison method via coupling in [@bourgade2016fixed]. As in [@landon2019fixed], consider the interpolation between a general sample covariance matrix and the Wishart ensemble for the initial data: for any $ \nu \in [0,1] $, let $$x_k^{(\nu)}(0)=\nu s_k(0)+(1-\nu)r_k(0),\ \ \ -N {\leqslant}k {\leqslant}N,\ k \neq 0.$$ where $ s_k(t) $ and $ r_k(t) $ satisfy the singular values dynamics , with respective initial conditions a general sample covariance matrix and the Wishart ensemble. Define the corresponding dynamics of $ x_k^{(\nu)} $ to be $$\label{e.DBM}
{\mathrm{d}}x_k^{(\nu)} = \dfrac{{\mathrm{d}}B_k}{\sqrt{N}} + \left[ -\dfrac{1}{2\xi}x_k^{(\nu)} + \dfrac{1}{2}\left( \dfrac{1}{\xi} -1 \right) \dfrac{1}{x_k^{(\nu)}} + \dfrac{1}{2N} \sum_{l \neq \pm k} \dfrac{1}{x_k^{(\nu)} - x_l^{(\nu)}} \right]{\mathrm{d}}t .$$ For this Dyson Brownian motion we consider the quantity $${\mathfrak{u}}_k^{(\nu)}(t):=e^{\frac{t}{2\xi}}\dfrac{{\mathrm{d}}}{{\mathrm{d}}\nu}x_k^{(\nu)}(t).$$ From now we set $ \nu \in (0,1) $ and omit it from the notation for simplicity. A significant property is that $ {\mathfrak{u}}_k $ satisfies a non-local parabolic differential equation $$\label{e.parabolic}
\dfrac{{\mathrm{d}}}{{\mathrm{d}}t}{\mathfrak{u}}_{k}=\dfrac{1}{2}\left( 1- \dfrac{1}{\xi} \right)\dfrac{{\mathfrak{u}}_k}{x_k^2} + \dfrac{1}{2N} \sum_{l \neq \pm k}\dfrac{{\mathfrak{u}}_l - {\mathfrak{u}}_k}{\left(x_l - x_k \right)^2}.$$ Let $ {\mathfrak{v}}_k={\mathfrak{v}}_k^{(\nu)} $ solve the same equation as $ {\mathfrak{u}}_k $ in but with initial condition $ {\mathfrak{v}}_k(0)=|{\mathfrak{u}}_k(0)|=|s_k(0)-r_k(0)| $. An important result is that this equation yields a maximum principle for $ {\mathfrak{v}}_k $, which will be useful in later analysis for the estimate of its growth.

\[l.MaxPrin\] For all $ t {\geqslant}0 $ and $ -N {\leqslant}k {\leqslant}N $, we have $${\mathfrak{v}}_k(t) {\geqslant}0, \ \ \ |{\mathfrak{v}}_k(t)| {\leqslant}\max_k |{\mathfrak{v}}_k(0)|,\ \ \ |{\mathfrak{u}}_k(t)| {\leqslant}{\mathfrak{v}}_k(t).$$

Note that for the coefficients in the summation part of equation we have $ \frac{1}{(x_l - x_k)^2} >0 $. Therefore for $ f(t) := \min_k {\mathfrak{v}}_k(t) $, we have $$f'(t) {\geqslant}\dfrac{1}{2} \left( 1-\dfrac{1}{\xi} \right) \dfrac{1}{x_1^2} f(t).$$ Combined with the fact $ {\mathfrak{v}}_k(0) {\geqslant}0 $, this gives us the first claim. For the second claim, since $ {\mathfrak{v}}_k $’s are nonnegative, we know that $ \frac{{\mathrm{d}}}{{\mathrm{d}}t}\max_{k} {\mathfrak{v}}_k  {\leqslant}0 $ and this yields the desired result. The third claim follows from linearity. We know that both $ {\mathfrak{v}}+ {\mathfrak{u}}$ and $ {\mathfrak{v}}- {\mathfrak{u}}$ satisfy the equation , and we also have $ ({\mathfrak{v}}+ {\mathfrak{u}})(0) {\geqslant}0 $ and $ ({\mathfrak{v}}- {\mathfrak{u}})(0) {\geqslant}0 $. Similarly to the first claim, this gives us $ {\mathfrak{v}}_k(t) + {\mathfrak{u}}_k(t) {\geqslant}0 $ and $ {\mathfrak{v}}_k(t) - {\mathfrak{u}}_k(t)  {\geqslant}0 $, which completes the proof.

Due to our choice for the initial value of $ {\mathfrak{v}}$, we will have that $ \{{\mathfrak{v}}_k\}_{-N {\leqslant}k {\leqslant}N} $ are symmetric with respect to the label $ k $. To see this, it is easy to check that $ {\widetilde}{{\mathfrak{v}}}_k := {\mathfrak{v}}_{-k} $ satisfy the same equation . Note that the equation is linear and $ {\widetilde}{{\mathfrak{v}}}_k(0)={\mathfrak{v}}_k(0) $ for all $ k $. Lemma \[l.MaxPrin\] then gives us $ {\mathfrak{v}}_k(t) = {\widetilde}{{\mathfrak{v}}}_k(t) = {\mathfrak{v}}_{-k}(t) $ for all $ k $ and $ t {\geqslant}0 $.

We now consider the observable $$f_t(z)=e^{-\frac{t}{2\xi}}\sum_{-N {\leqslant}k {\leqslant}N}\dfrac{{\mathfrak{v}}_k(t)}{x_k(t)-z}.$$ A key observation is that the quadratic singularities in will disappear when combined with the Dyson Brownian motion, so that the evolution of $ f_t $ has no shocks similarly to a result from [@bourgade2018extreme]. Denoting $$s_t(z)=\dfrac{1}{2N}\sum_{-N {\leqslant}k {\leqslant}N}\dfrac{1}{x_k(t)-z}$$ the Stieltjes transform of the (symmetrized) empirical spectral measure, then the observable $ f_t $ satisfies the following dynamics.

\[l.dynamics\] For any $  {\text{\rm Im}\hspace{0.1cm}}z \neq 0 $, we have $$\label{e.dynamics}
\begin{aligned}
{\mathrm{d}}f_t &=\left(s_t(z)+\dfrac{z}{2\xi}\right)(\partial_z f_t){\mathrm{d}}t+\dfrac{1}{4N}(\partial_{zz}f_t){\mathrm{d}}t+\left[\dfrac{e^{-\frac{t}{2\xi}}}{2N}\sum_{-N {\leqslant}k {\leqslant}N}\dfrac{{\mathfrak{v}}_k}{(x_k-z)^2(x_k+z)}\right]{\mathrm{d}}t\\
& \quad + \left[\left( 1- \dfrac{1}{\xi} \right) e^{-\frac{t}{2\xi}} \left( \sum_{-N {\leqslant}k {\leqslant}N} \dfrac{3z {\mathfrak{v}}_k}{2 x_k^2 (x_k -z)(x_k+z)} + \sum_{-N {\leqslant}k {\leqslant}N} \dfrac{z^3 {\mathfrak{v}}_k}{x_k^2 (x_k-z)^2 (x_k+z)^2} \right)\right] {\mathrm{d}}t\\
& \quad -\dfrac{e^{-\frac{t}{2\xi}}}{\sqrt{N}}\sum_{-N {\leqslant}k {\leqslant}N}\dfrac{{\mathfrak{v}}_k}{(x_k-z)^2}{\mathrm{d}}B_k.
\end{aligned}$$

This can be proved by direct computation via the Itô’s formula. First, we have $${\mathrm{d}}f=-\dfrac{f}{2\xi}{\mathrm{d}}t+e^{-t/2}\sum_{-N {\leqslant}k {\leqslant}N}\dfrac{{\mathrm{d}}{\mathfrak{v}}_k}{x_k-z}+e^{-t/2}\sum_{-N {\leqslant}k {\leqslant}N}{\mathfrak{v}}_k {\mathrm{d}}\dfrac{1}{x_k-z} =: A_1 + A_2 +A_3.$$ By using Itô’s formula again we have $ {\mathrm{d}}(x_k-z)^{-1}=-(x_k-z)^{-2}{\mathrm{d}}x_k+\frac{1}{N}(x_k-z)^{-3}{\mathrm{d}}t $. Thus, we can now decompose the term $ A_3 $ as $ I_1+[I_2+I_3+I_4+I_5]{\mathrm{d}}t $, where $$\begin{aligned}
I_1 &=
-\dfrac{e^{-\frac{t}{2\xi}}}{\sqrt{N}}\sum_{-N {\leqslant}k {\leqslant}N}\dfrac{{\mathfrak{v}}_k}{(x_k-z)^2}{\mathrm{d}}B_k, \\
I_2 &=
\dfrac{e^{-\frac{t}{2\xi}}}{2\xi}\sum_{-N {\leqslant}k {\leqslant}N}\dfrac{{\mathfrak{v}}_k x_k}{(x_k-z)^2}, \\
I_3 &=
\dfrac{1}{2} \left( 1-\dfrac{1}{\xi} \right) e^{-\frac{t}{2\xi}} \sum_{-N {\leqslant}k {\leqslant}N} \dfrac{{\mathfrak{v}}_k}{x_k(x_k-z)^2}, \\
I_4 &= e^{-\frac{t}{2\xi}}\sum_{-N {\leqslant}k {\leqslant}N}{\mathfrak{v}}_k\left(-\dfrac{1}{(x_k-z)^2}\right)\dfrac{1}{2N}\sum_{l \neq \pm k}\dfrac{1}{x_k-x_l}=\dfrac{e^{-\frac{t}{2\xi}}}{2N}\sum_{l \neq \pm k}\dfrac{{\mathfrak{v}}_k}{(x_k-z)^2(x_l-x_k)}, \\
I_5 &= \dfrac{e^{-\frac{t}{2\xi}}}{N}\sum_{-N {\leqslant}k {\leqslant}N}\dfrac{{\mathfrak{v}}_k}{(x_k-z)^3}.\end{aligned}$$ Note that $$\partial_z f=e^{-\frac{t}{2\xi}}\sum_{-N {\leqslant}k {\leqslant}N}\dfrac{{\mathfrak{v}}_k}{(x_k-z)^2},\ \ \ \partial_{zz}f=2e^{-\frac{t}{2\xi}}\sum_{-N {\leqslant}k {\leqslant}N}\dfrac{{\mathfrak{v}}_k}{(x_k-z)^3}.$$ For the term $ A_2 $, by the equation , we have $$A_2 = e^{-\frac{t}{2\xi}}\sum_{-N {\leqslant}k {\leqslant}N}\dfrac{1}{x_k-z} \left[ \dfrac{1}{2} \left( 1-\dfrac{1}{\xi} \right) \dfrac{{\mathfrak{v}}_k}{x_k^2}  + \dfrac{1}{2N}\sum_{l \neq \pm k}\dfrac{{\mathfrak{v}}_l-{\mathfrak{v}}_k}{(x_k-x_l)^2} \right]=: B_1+B_2.$$ Note that $$\begin{gathered}
B_2 = \dfrac{e^{-\frac{t}{2\xi}}}{4N}\sum_{l \neq \pm k}\dfrac{{\mathfrak{v}}_l-{\mathfrak{v}}_k}{(x_l-x_k)^2}\left(\dfrac{1}{x_k-z}-\dfrac{1}{x_l-z}\right)=\dfrac{e^{-\frac{t}{2\xi}}}{4N}\sum_{l \neq \pm k}\dfrac{{\mathfrak{v}}_l-{\mathfrak{v}}_k}{x_l-x_k}\dfrac{1}{x_k-z}\dfrac{1}{x_l-z}\\
=-\dfrac{e^{-\frac{t}{2\xi}}}{2N}\sum_{l \neq \pm k}\dfrac{{\mathfrak{v}}_k}{x_l-x_k}\dfrac{1}{x_k-z}\dfrac{1}{x_l-z}.\end{gathered}$$ Combining with $ I_4 $, we obtain $$B_2 + I_4 =\dfrac{e^{-\frac{t}{2\xi}}}{2N}\sum_{l \neq \pm k}\dfrac{{\mathfrak{v}}_k}{x_l-x_k}\dfrac{1}{x_k-z}\left(\dfrac{1}{x_k-z}-\dfrac{1}{x_l-z}\right)=\dfrac{e^{-\frac{t}{2\xi}}}{2N}\sum_{l \neq \pm k}\dfrac{{\mathfrak{v}}_k}{(x_k-z)^2}\dfrac{1}{x_l-z}.$$ Moreover, we have that $$B_2+I_4+I_5 = s(z)\partial_z f+\dfrac{e^{-\frac{t}{2\xi}}}{2N}\sum_{-N {\leqslant}k {\leqslant}N}\dfrac{{\mathfrak{v}}_k}{(x_k-z)^3}+\dfrac{e^{-\frac{t}{2\xi}}}{2N}\sum_{-N {\leqslant}k {\leqslant}N}\dfrac{{\mathfrak{v}}_k}{(x_k-z)^2(x_k+z)}.$$ Then it suffices to calculate $ B_1 + I_3 $, and note that $$\begin{aligned}
B_1 + I_3 &= \dfrac{1}{2} \left( 1- \dfrac{1}{\xi} \right) e^{-\frac{t}{2\xi}} \sum_{k=1}^N \left( \dfrac{{\mathfrak{v}}_k}{x_k -z}\dfrac{1}{x_k^2} - \dfrac{{\mathfrak{v}}_k}{x_k +z}\dfrac{1}{x_k^2} + \dfrac{{\mathfrak{v}}_k}{(x_k - z)^2}\dfrac{1}{x_k} - \dfrac{{\mathfrak{v}}_k}{(x_k +z)^2}\dfrac{1}{x_k} \right)\\
&= \left( 1- \dfrac{1}{\xi} \right) e^{-\frac{t}{2\xi}} \left( \sum_{-N {\leqslant}k {\leqslant}N} \dfrac{3z {\mathfrak{v}}_k}{2 x_k^2 (x_k -z)(x_k+z)} + \sum_{-N {\leqslant}k {\leqslant}N} \dfrac{z^3 {\mathfrak{v}}_k}{x_k^2 (x_k-z)^2 (x_k+z)^2} \right).\end{aligned}$$ Note that now all the singularities are removed.

The desired result then follows by combining the previous results and the term $ I_1 $.

Recall that in Section \[s.LocalLaw\] we have shown that the Stieltjes transform of the empirical measure for singular values satisfies the local law , so that the leading term of the stochastic differential equation satisfied by $ f_t $ is close to $$\dfrac{z}{2\xi} + z m_{{\mathsf{MP}}}(z^2) = \dfrac{z}{2\xi} + \dfrac{1-\xi-z^2 + \sqrt{(z^2 - \lambda_-)(z^2 - \lambda_+)}}{2\xi z} = \dfrac{(1-\xi) + \sqrt{(z^2 - \lambda_-)(z^2 - \lambda_+)}}{2\xi z}.$$ Thus, the dynamics of $ f_t $ can be approximated by the following advection equation $$\label{e.AdvectionPDE}
\partial_t r = \dfrac{(1-\xi) + \sqrt{(z^2 - \lambda_-)(z^2 - \lambda_+)}}{2\xi z} \partial_z r.$$

Geometric properties of the characteristics {#s.geometry}
-------------------------------------------

In order to estimate the evolution of the observable, we analyze its dynamics by studying the characteristics of the approximate advection PDE , similarly to [@huang2018rigidity; @bourgade2018extreme]. To do this, we first need some bounds on the shape of the characteristics $ (z_t)_{t {\geqslant}0} $, and some estimates for the initial value. As mentioned in Section \[s.LocalLaw\], discussions for the case $ \xi=1 $ are expected to be different due to the singularity of the Marchenko-Pastur distribution (which results in a distinct shape of the density $ \rho(x) $ for singular values). Therefore, in this and the subsequent section, we first discuss the case $ \xi \neq 1 $ and will show how to adapt the proof to the case for square data matrices in Section \[s.CaseOne\].

Denote $$\kappa(z) := \min\left\{\left|z-\sqrt{\lambda_-}\right|,\left|z-\sqrt{\lambda_+}\right|\right\},$$ and $$a(z) := \operatorname{dist}\left(z,\left[\sqrt{\lambda_-},\sqrt{\lambda_+}\right]\right),\ \ \ b(z):= \operatorname{dist}\left( z, \left[\sqrt{\lambda_-},\sqrt{\lambda_+}\right]^c \right)$$ We consider the curve $${\mathscr{S}}:= \left\{ z=E+{\text{\rm i}}y: \sqrt{\lambda_-} + \varphi^2 N^{-2/3} <E< \sqrt{\lambda_+} - \varphi^2 N^{-2/3},\ y=\varphi^2/\left( N \kappa(E)^{1/2} \right) \right\},$$ and the domain $ \mathscr{R} := \cup_{0<t<1}\{z_t:z \in {\mathscr{S}}\} $.

\[l.characteristics\] Uniformly in $ 0<t<1 $ and $ z=z_0 $ satisfying $ \eta := {\text{\rm Im}\hspace{0.1cm}}z>0 $ and $ |z-\sqrt{\lambda_+}|<\sqrt{\xi}/10 $, we have $${\text{\rm Re}\hspace{0.1cm}}(z_t -z_0) \sim t \dfrac{a(z)}{\kappa(z)^{1/2}}+t^2,\ \ \ \ {\text{\rm Im}\hspace{0.1cm}}(z_t-z_0) \sim t\dfrac{b(z)}{\kappa(z)^{1/2}}.$$ In particular, if in addition we have $ z \in {\mathscr{S}}$, then $$(z_t - z_0) \sim \left( t \dfrac{\varphi^2}{N \kappa(E)} +t^2 \right) + {\text{\rm i}}\kappa(E)^\frac{1}{2} t.$$ Moreover, for any $ \kappa >0 $, uniformly in $ 0<t<1 $ and $ z=E+{\text{\rm i}}\eta \in [\sqrt{\lambda_-} + \kappa,\sqrt{\lambda_+} -\kappa] \times [0,\kappa^{-1}] $, we have $ {\text{\rm Im}\hspace{0.1cm}}(z_t-z_0) \sim t $.

![Shape of the characteristics near the edge: the solid lines represent the density $ \rho(x) $ and the corresponding characteristics; the dashed lines represent the semicircle distribution and the related characteristics.](characteristics.jpg){height="4cm"}

It is too complicated to work with the ODE satisfied by $ z_t $ in a direct way. The main idea is to compare this characteristics with the corresponding curve for a semicircle distribution (which has a explicit and simple formula). Define the two functions $$g(z):=\dfrac{(1-\xi) + \sqrt{(z^2 - \lambda_-)(z^2 - \lambda_+)}}{2\xi z} , \ \ \ g_{sc}(z):=\dfrac{\sqrt{(z-1)^2 - \xi}}{\xi}.$$ For some $ |z_0-\sqrt{\lambda_+}| < \sqrt{\xi}/10,\ {\text{\rm Im}\hspace{0.1cm}}z_0>0 $, let $ z(t) := z_t $ and $ z_{sc}(t) $ solve the following two initial value problems $$\left\{
\begin{aligned}
& \dfrac{{\mathrm{d}}z}{{\mathrm{d}}t} = g(z)\\
& z(0)=z_0
\end{aligned}
\right.
,\ \ \ \ \ \ \
\left\{
\begin{aligned}
& \dfrac{{\mathrm{d}}z_{sc}}{{\mathrm{d}}t} = g_{sc}(z_{sc})\\
& z_{sc}(0)=z_0
\end{aligned}
\right.
.$$ Note that in $ \Omega=\{z_{sc}(t):0<t<1,|z_0-\sqrt{\lambda_+}| < \sqrt{\xi}/10,{\text{\rm Im}\hspace{0.1cm}}z_0>0 \} $ we have $ g(z) \sim g_{sc}(z) $. This shows that for $ 0<t<1 $, we have $ (z_t-z_0) \sim (z_{sc}(t) -z_0) $. The rest of the proof now follows from [@bourgade2018extreme Lemma 2.2].

Furthermore, we have the following lemma regarding the growth of the characteristics, which will be useful for the error estimates in the local relaxation.

\[l.ChracteristcsInt\] For any $ z=E+{\text{\rm i}}\eta \in {\mathscr{S}}$, we have $$\dfrac{\varphi^4}{N^2} \int_0^t ds \int \dfrac{d \rho(x)}{|z_{t-s} -x|^4 \max(\kappa(x),s^2)} \lesssim \dfrac{\kappa(E)}{\max(\kappa(E),t^2)}.$$

The proof is essentially the same as in [@bourgade2018extreme Lemma A.2], except now we need to consider the Stieltjes transform $ m(z) $ instead of the one for the semicircle law. Recall the relation $ m(z)=z m_{{\mathsf{MP}}}(z^2) $, and note that for the Stieltjes transform of the Marchenko-Pastur law we have the following results from [@bloemendal2016principal Lemma 3.6]: $$|m_{{\mathsf{MP}}}(z)| \sim 1,
\ \ \ \ \
{\text{\rm Im}\hspace{0.1cm}}m_{{\mathsf{MP}}}(z) \sim \left\{
\begin{aligned}
& \sqrt{\kappa(E)+\eta} & \mbox{if } & E \in [\lambda_-,\lambda_+],\\
&\frac{\eta}{\sqrt{\kappa(E)+\eta}} & \mbox{if } & E \notin [\lambda_-,\lambda_+].
\end{aligned}
\right.$$ The rest of the proof follows from the calculation in [@bourgade2018extreme Lemma A.2].

We now consider the initial value $ f_0 $ on the curve $ {\mathscr{S}}$. For this purpose we define the set of good trajectories such that the rigidity holds: $${\mathscr{A}}:= \left\{ |x_k^{(\nu)}(t) - \gamma_k|<\varphi^{\frac{1}{2}}N^{-\frac{2}{3}}({\widehat}{k})^{-\frac{1}{3}} \ \mbox{for all}\ 0 {\leqslant}t {\leqslant}1, k \in \llbracket 1,N \rrbracket,0 {\leqslant}\nu {\leqslant}1 \right\}.$$ We have the following important estimate for the probability of these events.

\[l.rigidity\] There exists a fixed $ C_0 $ large enough such that the following holds. For any $ D>0 $, there exists $ N_0(D)>0 $ such that for any $ N > N_0 $ we have $${\mathbb{P}}({\mathscr{A}}) > 1-N^{-D}.$$

As mentioned in the last section, the rigidity estimates are proved in [@pillai2014universality] for fixed $ t $ and $ \nu=0,1 $. The extension to all $ t $ and $ \nu $ is based on the arguments in [@bourgade2018extreme; @erdHos2015gap]: (1) discretize in $ t $ and $ \nu $; (2) use Weyl’s inequality to bound the increments in small time intervals; (3) use the maximum principle to bound the increment in small $ \nu $-intervals.

Conditioned on the rigidity phenomenon, we have the following estimate for the initial conditions. The proof is the same as [@bourgade2018extreme Lemma 2.4].

\[l.InitialValue\] In the set $ {\mathscr{A}}$, for any $ z=E+{\text{\rm i}}\eta \in \mathscr{R} $, we have $ {\text{\rm Im}\hspace{0.1cm}}f_0(z) \lesssim \varphi^{1/2} $ if $ \eta > \max(E-\sqrt{\lambda_+},-E+\sqrt{\lambda_-}) $, and $ {\text{\rm Im}\hspace{0.1cm}}f_0(z) \lesssim \varphi^{1/2} \frac{\eta}{\kappa(z)} $ otherwise. The same bound also holds for $ |{\text{\rm Im}\hspace{0.1cm}}f_0| $.

Quantitative relaxation at the edge {#s.relaxation}
-----------------------------------

To prove the edge universality, we first have the following estimate for the size of the observable $ f_t $.

\[p.estimate\] For any (large) $ D>0 $ there exists $ N_0(D) $ such that for any $ N {\geqslant}N_0 $ we have $${\mathbb{P}}\left( {\text{\rm Im}\hspace{0.1cm}}f_t(z) \lesssim \varphi \dfrac{\kappa(E)^{1/2}}{\max(\kappa(E)^{1/2},t)} \ \mbox{\rm for all} \ 0<t<1 \ \mbox{\rm and} \ z=E+{\text{\rm i}}\eta \in {\mathscr{S}}\right) > 1-N^{-D}$$

For any $ 1 {\leqslant}l,m {\leqslant}N^{10} $, we define $ t_l:=l N^{-10} $ and $$z^{(m)} := E_m + {\text{\rm i}}\eta_m=E_m+{\text{\rm i}}\dfrac{\varphi^2}{N \kappa(E_m)^{1/2}},$$ where $$E_m := \inf\left\{ E >0 : \int_{-\infty}^{E^2} \rho_{{\mathsf{MP}}}(x){\mathrm{d}}x {\geqslant}\left(m-\frac{1}{2}\right)N^{-10} \right\}.$$ We also define the following stopping times (with respect to $ \mathcal{F}_t = \sigma(B_k(s):0 {\leqslant}s {\leqslant}t,1 {\leqslant}k {\leqslant}N) $) which represent the bad events: $$\tau_{l,m} := \inf \left\{ 0 {\leqslant}s {\leqslant}t_l : {\text{\rm Im}\hspace{0.1cm}}f_s(z_{t-s}^{(m)}) \gtrsim \dfrac{\varphi}{2}\dfrac{\kappa(E_m)^{1/2}}{\max(\kappa(E_m)^{1/2},t_l)}  \right\},$$ $$\tau_0 := \inf \left\{ 0 {\leqslant}t {\leqslant}1: \exists k \in \llbracket -N,N \rrbracket \ \mbox{s.t.} \  |x_k(t)-\gamma_k|>\varphi^{\frac{1}{2}}N^{-\frac{2}{3}}({\widehat}{k})^{-\frac{1}{3}} \right\},$$ $$\tau := \min \left\{ \tau_0,\tau_{l,m}: 0 {\leqslant}l,m {\leqslant}N^{10}, \kappa(E_m) > \varphi^2 N^{-\frac{2}{3}} \right\}.$$ We also define the convention $ \inf \emptyset =1 $.

We claim that in order to prove the desired result, it suffices to show that for any $ D>0 $ there exists $ {\widetilde}{N}_0(D) $ such that for any $ N {\geqslant}{\widetilde}{N}_0(D) $, we have $$\label{e.ProbabilityTau}
{\mathbb{P}}(\tau=1) > 1-N^{-D}.$$

*Step 1.* To see the claim would be enough, we first show the following sets inclusion $$\label{e.SetsInclusion}
\{ \tau=1 \} \bigcap_{\substack{1{\leqslant}l,m {\leqslant}N^{10} \\ -N {\leqslant}k {\leqslant}N}} A_{l,m,k} {\subseteq}\bigcap_{z \in {\mathscr{S}}, 0<t<1} \left\{ {\text{\rm Im}\hspace{0.1cm}}f_t(z) \lesssim \varphi \dfrac{\kappa(E)^{1/2}}{\max(\kappa(E)^{1/2},t)} \right\},$$ where $$A_{l,m,k} := \left\{ \sup_{t_l {\leqslant}u {\leqslant}t_{l+1}} \left| \int_{t_l}^u \dfrac{e^{-\frac{s}{2\xi}} {\mathfrak{v}}_k(s) {\mathrm{d}}B_k(s)}{(z^{(m)} - x_k(s))^2} \right| <N^{-3} \right\}.$$ To prove this, for any given $ z $ and $ t $, choose $ t_l $ and $ z^{(m)} $ such that $ t_l {\leqslant}t <t_{l+1} $ and $ |z-z^{(m)}|<N^{-5} $. Note that by rigidity and the maximum principle (Lemma \[l.MaxPrin\]) we have $ |{\mathfrak{v}}_k(t)| \lesssim \varphi N^{-2/3} $. Combining this with the definition of $ f_t $, we have $ |f_t(z) - f_t(z^{(m)})| < N^{-2} $. Moreover, note that we have the following estimates $$\begin{gathered}
\sup_{t_l {\leqslant}u {\leqslant}t_{l+1}} \left| \int_{t_l}^u \dfrac{e^{-\frac{t}{2\xi}}}{2N}\sum_{-N {\leqslant}k {\leqslant}N} \dfrac{{\mathfrak{v}}_k}{(x_k-z)^2(x_k+z)} {\mathrm{d}}t \right|\\
\lesssim N^{-10}  \dfrac{1}{2N} \varphi N^{-2/3} \max_{t_l {\leqslant}t {\leqslant}t_{l+1}} \sum_{-N {\leqslant}k {\leqslant}N} \left| \dfrac{1}{(x_k-z)^2(x_k+z)} \right| < N^{-5},\end{gathered}$$ and $$\begin{gathered}
\sup_{t_l {\leqslant}u {\leqslant}t_{l+1}} \left| \int_{t_l}^u \left( 1 - \dfrac{1}{\xi} \right) e^{-\frac{t}{2\xi}} \sum_{-N {\leqslant}k {\leqslant}N} \dfrac{z {\mathfrak{v}}_k}{x_k^2(x_k -z) (x_k+z)} {\mathrm{d}}t \right|\\
\lesssim N^{-10} \varphi N^{-2/3} \max_{t_l {\leqslant}t {\leqslant}t_{l+1}} \sum_{-N {\leqslant}k {\leqslant}N} \left| \dfrac{1}{(x_k -z) (x_k+z)} \right| < N^{-5},\end{gathered}$$ $$\begin{gathered}
\sup_{t_l {\leqslant}u {\leqslant}t_{l+1}} \left| \int_{t_l}^u \left( 1 - \dfrac{1}{\xi} \right) e^{-\frac{t}{2\xi}} \sum_{-N {\leqslant}k {\leqslant}N} \dfrac{z^3 {\mathfrak{v}}_k}{ x_k^2 (x_k -z)^2 (x_k+z)^2} {\mathrm{d}}t \right|\\
\lesssim N^{-10} \varphi N^{-2/3} \max_{t_l {\leqslant}t {\leqslant}t_{l+1}} \sum_{-N {\leqslant}k {\leqslant}N} \left| \dfrac{1}{(x_k -z)^2 (x_k+z)^2} \right| < N^{-5}.\end{gathered}$$ Similarly, we also have $$\sup_{t_l {\leqslant}u {\leqslant}t_{l+1}} \left| \int_{t_l}^u  \left(s_t(z)+\dfrac{z}{2\xi}\right)(\partial_z f_t)  {\mathrm{d}}t \right| < N^{-5},\ \ \ \  \sup_{t_l {\leqslant}u {\leqslant}t_{l+1}} \left| \int_{t_l}^u  \dfrac{1}{4N}(\partial_{zz}f_t)  {\mathrm{d}}t \right| < N^{-5}.$$ Then based on the dynamics , under the event $ \cap_k A_{l,m,k} $, we have $ |f_t(z^{(m)}) - f_{t_l}(z^{(m)})| < N^{-2} $. This completes the proof for .

*Step 2.* To prove the claim, we also need to estimate the probability for the events $ A_{l,m,k} $. To do this, we use the following Burkholder-Davis-Gundy type inequality (see e.g. [@shorack2009empirical Appendix B.6]). For some fixed constants $ c>0 $ and $ \alpha>0 $, for any martingale $ M $ we have $$\label{e.BDG}
{\mathbb{P}}\left( \sup_{0 {\leqslant}u {\leqslant}t} |M_u| {\geqslant}\alpha {\langle}M {\rangle}_t^{\frac{1}{2}} \right) {\leqslant}e^{-c \alpha^2}.$$ Note that we have the deterministic bound $ \int_{t_l}^u \frac{|{\mathfrak{v}}_k (s)|^2 {\mathrm{d}}s}{|z^{(m)} - x_k(s)|^4} \ll N^{-6} $. By taking $ \alpha = \varphi^{1/10} $, this implies $ {\mathbb{P}}(A_{l,m,k}) {\geqslant}1 - e^{-c \varphi^{1/5}} $, and then a union bound yields $${\mathbb{P}}\left( \bigcap_{1{\leqslant}l,m {\leqslant}N^{10}, -N {\leqslant}k {\leqslant}N} A_{l,m,k} \right) > 1-N^{-D}.$$ Together with the sets inclusion , this concludes that implies the desired result.

*Step 3.* It remains to prove . For simplicity of notations, let $ t=t_l $, $ z=E+{\text{\rm i}}\eta=z^{(m)} $ for some arbitrary fixed $ 1 {\leqslant}l,m {\leqslant}N^{10} $. Consider the function $ g_u(z) := f_u(z_{t-u}) $. By Lemma \[l.characteristics\] and Lemma \[l.InitialValue\], we have $ {\text{\rm Im}\hspace{0.1cm}}g_0(z) \lesssim \frac{\varphi}{10}\frac{\kappa(E_m)^{1/2}}{\max(\kappa(E_m)^{1/2},t)} $. Therefore we only need to bound the increments of $ g $. Using Lemma \[l.dynamics\], by the Itô’s formula we know it satisfies the following stochastic differential equation $$d g_{u \wedge \tau}(z) = {\varepsilon}_u(z_{t-u}){\mathrm{d}}(u \wedge \tau) - \dfrac{e^{-\frac{u}{2\xi}}}{\sqrt{N}}\sum_{-N {\leqslant}k {\leqslant}N}\dfrac{{\mathfrak{v}}_k(u)}{(z_{t-u}-x_k(u))^2}{\mathrm{d}}B_k(u \wedge \tau),$$ where $$\begin{aligned}
{\varepsilon}_u(z) &:= (s_u(z)-m(z))\partial_z f_u + \dfrac{1}{4N} (\partial_{zz} f_u) + \dfrac{e^{-\frac{u}{2\xi}}}{2N} \sum_{-N {\leqslant}k {\leqslant}N} \dfrac{{\mathfrak{v}}_k(u)}{(x_k -z)^2 (x_k+z)}\\
&\quad  + \left( 1- \dfrac{1}{\xi} \right) e^{-\frac{t}{2\xi}} \left( \sum_{-N {\leqslant}k {\leqslant}N} \dfrac{3z {\mathfrak{v}}_k}{2 x_k^2 (x_k -z)(x_k+z)} + \sum_{-N {\leqslant}k {\leqslant}N} \dfrac{z^3 {\mathfrak{v}}_k}{x_k^2 (x_k-z)^2 (x_k+z)^2} \right).\end{aligned}$$ In this step, we aim to estimate the first term $ \sup_{0 {\leqslant}s {\leqslant}t}|\int_0^s {\varepsilon}_u(z_{t-u}) {\mathrm{d}}(u \wedge \tau)| $. First, we have $$\begin{gathered}
\int_0^t \left| \left( s_u(z_{t-u}) -m(z_{t-u}) \right) \partial_z f (z_{t-u}) \right|{\mathrm{d}}(u \wedge \tau)\\ \lesssim \int_0^t \dfrac{\varphi}{N {\text{\rm Im}\hspace{0.1cm}}(z_{t-u})} \sum_{-N {\leqslant}k {\leqslant}N} \dfrac{|{\mathfrak{v}}_k(u)|}{|z_{t-u} -x_k(u)|^2} {\mathrm{d}}(u \wedge \tau)
\lesssim \int_0^t \dfrac{\varphi {\text{\rm Im}\hspace{0.1cm}}f_u(z_{t-u})}{N \left( {\text{\rm Im}\hspace{0.1cm}}(z_{t-u}) \right)^2}{\mathrm{d}}(u \wedge \tau) \\
\lesssim \int_0^t \dfrac{\varphi^2 du}{N (\eta + (t-u)\kappa(z)^{1/2})^2}\dfrac{\kappa(E)^{1/2}}{\max(\kappa(E)^{1/2},t)} \lesssim \dfrac{\kappa(E)^{1/2}}{\max \left( \kappa(E)^{1/2},t \right)}.\end{gathered}$$ To bound the $ |s_u - m| $ term above, we have used the local law for singular values simultaneously for all $ 0 {\leqslant}u {\leqslant}t $ (which is similar to Lemma \[l.rigidity\]). The last two inequalities follow from Lemma \[l.characteristics\]. We also have $$\sup_{0 {\leqslant}s {\leqslant}t} \left| \int_0^s \dfrac{1}{4N}(\partial_{zz} f_u(z_{t-u})) {\mathrm{d}}(u \wedge \tau) \right| \lesssim \int_0^t \dfrac{{\text{\rm Im}\hspace{0.1cm}}f_u(z_{t-u})}{N \left( {\text{\rm Im}\hspace{0.1cm}}(z_{t-u}) \right)^2}{\mathrm{d}}(u \wedge \tau)  \lesssim \dfrac{\kappa(E)^{1/2}}{\varphi \max(\kappa(E)^{1/2},t)},$$ $$\begin{gathered}
\sup_{0 {\leqslant}s {\leqslant}t} \left| \int_0^s \dfrac{e^{-\frac{u}{2\xi}}}{2N} \sum_{-N {\leqslant}k {\leqslant}N} \dfrac{{\mathfrak{v}}_k(u)}{(x_k - z_{t-u})^2(x_k+z_{t-u})} {\mathrm{d}}(u \wedge \tau) \right|\\
\lesssim \int_0^t  \dfrac{{\text{\rm Im}\hspace{0.1cm}}f_u(z_{t-u})}{N \left( {\text{\rm Im}\hspace{0.1cm}}(z_{t-u}) \right)}{\mathrm{d}}(u \wedge \tau) \lesssim \dfrac{\kappa(E)^{1/2}}{\varphi \max(\kappa(E)^{1/2},t)},\end{gathered}$$ And similarly, $$\begin{gathered}
\sup_{0 {\leqslant}s {\leqslant}t} \left| \left( 1-\dfrac{1}{\xi} \right)e^{-\frac{u}{2\xi}} \int_0^s \sum_{-N {\leqslant}k {\leqslant}N} \dfrac{z_{t-u} {\mathfrak{v}}_k(u)}{x_k^2(x_k(u) - z_{t-u}) (x_k(u) + z_{t-u})} {\mathrm{d}}(u \wedge \tau) \right|\\
\lesssim \int_0^t {\text{\rm Im}\hspace{0.1cm}}f_u(z_{t-u}) {\mathrm{d}}(u \wedge \tau)  \lesssim \dfrac{\varphi}{2}\dfrac{\kappa(E)^{1/2}}{\max \left( \kappa(E)^{1/2},t \right)},\end{gathered}$$ $$\begin{gathered}
\sup_{0 {\leqslant}s {\leqslant}t} \left| \left( 1-\dfrac{1}{\xi} \right)e^{-\frac{u}{2\xi}} \int_0^s \sum_{-N {\leqslant}k {\leqslant}N} \dfrac{z_{t-u}^3 {\mathfrak{v}}_k(u)}{x_k^2(x_k(u) - z_{t-u})^2 (x_k(u) + z_{t-u})^2} {\mathrm{d}}(u \wedge \tau) \right|\\
\lesssim \dfrac{\varphi}{2}\dfrac{\kappa(E)^{1/2}}{\max \left( \kappa(E)^{1/2},t \right)}.\end{gathered}$$

*Step 4.* Finally we focus on the estimate for $ \sup_{0 {\leqslant}s {\leqslant}t}|M_s| $ where $$M_s := \int_0^s \dfrac{e^{-\frac{u}{2\xi}}}{\sqrt{N}} \sum_{-N {\leqslant}k {\leqslant}N} \dfrac{{\mathfrak{v}}_k(u)}{(z_{t-u} - x_k(u))^2} {\mathrm{d}}B_k(u \wedge \tau).$$ Note that for all $ k $ and $ u < \tau $ we have $ |z_{t-u} - \gamma_k| \lesssim |z_{t-u} - x_k(u)| $ due to the fact $ |x_k(u) - \gamma_k| \ll |z_{t-u} - \gamma_k| $. Using again we have $$\sup_{0 {\leqslant}s {\leqslant}t}|M_s|^2 \lesssim \varphi^{\frac{1}{10}} \int_0^t \dfrac{1}{N}\sum_{-N {\leqslant}k {\leqslant}N} \dfrac{{\mathfrak{v}}_k(u)^2}{|z_{t-u} - \gamma_k|^4} {\mathrm{d}}(u \wedge \tau)$$ with overwhelming probability. By a similar argument in [@bourgade2018extreme equation (2.17)] and Lemma \[l.ChracteristcsInt\], we conclude $$\sup_{0 {\leqslant}s {\leqslant}t} |M_s|^2 \lesssim \varphi^{1/5} \dfrac{\kappa(E)}{\max(\kappa(E),t^2)}.$$

Hence, based on the previous estimates and a union bound we have proved that for any $ D>0 $ there exists $ N_0 $ such that for every $ N > N_0 $ we have $${\mathbb{P}}\left( \sup_{0 {\leqslant}l,m {\leqslant}N^{10},\kappa(E_m)>\varphi^2 N^{-2/3},0 {\leqslant}s {\leqslant}t_l} {\text{\rm Im}\hspace{0.1cm}}f_{s \wedge \tau}(z_{t_l - s \wedge \tau}^{(m)}) \lesssim \dfrac{\varphi}{2}\dfrac{\kappa(E_m)^{1/2}}{\max(\kappa(E_m)^{1/2},t_l)} \right) {\geqslant}1- N^{-D}.$$ Together with Lemma \[l.InitialValue\], we have proved , which then completes the proof.

Based on the previous estimate on $ {\text{\rm Im}\hspace{0.1cm}}f_t $, we now state the quantitative relaxation of the singular values dynamics at the edge. Remember that $ \{s_k\} $ and $ \{r_k\} $ satisfy the same equation , with respective initial conditions that correspond to a general sample covariance matrix and the Wishart ensemble.

\[t.Relaxation\] For any $ D>0 $ and $ {\varepsilon}>0 $ there exists $ N_0>0 $ such that for any $ N > N_0 $ we have $${\mathbb{P}}\left( |s_k(t) - r_k(t)| \lesssim \dfrac{N^{\varepsilon}}{Nt} \ \mbox{for all} \ k \in \llbracket 1,N \rrbracket \ \mbox{in} \ t \in [0,1] \right) > 1-N^{-D}.$$

Let $ z=\gamma_k + {\text{\rm i}}\frac{\varphi^2}{N \sqrt{\kappa(\gamma_k)}} \in {\mathscr{S}}$, then conditioned on the rigidity phenomenon $ {\mathscr{A}}$, by the nonnegativity of $ {\mathfrak{v}}_k $’s shown in Lemma \[l.MaxPrin\] we have $$|{\mathfrak{v}}_k(t)| \lesssim \dfrac{\varphi^2}{N \sqrt{\kappa(\gamma_k)}} {\text{\rm Im}\hspace{0.1cm}}f_t(z).$$ By the arguments in [@bourgade2018extreme Corollary 2.7], based on Proposition \[p.estimate\] we conclude $${\mathbb{P}}\left( |{\mathfrak{v}}_k^{(\nu)}(t)| \lesssim \dfrac{\varphi^{10}}{N} \dfrac{1}{\max(({\widehat}{k}/N)^{1/3},t)} \ \mbox{for all} \ k \in \llbracket 1,N \rrbracket \ \mbox{in} \ t \in [0,1] \right) > 1-N^{-D}.$$ Again by Lemma \[l.MaxPrin\], we know $ -{\mathfrak{v}}_k {\leqslant}{\mathfrak{u}}_k {\leqslant}{\mathfrak{v}}_k $. Therefore we have $${\mathbb{P}}\left( |{\mathfrak{u}}_k^{(\nu)}(t)| \lesssim \dfrac{\varphi^{10}}{N} \dfrac{1}{\max(({\widehat}{k}/N)^{1/3},t)} \ \mbox{for all} \ k \in \llbracket 1,N \rrbracket \ \mbox{in} \ t \in [0,1] \right) > 1-N^{-D}.$$ The result then follows as the proof in [@bourgade2018extreme Theorem 2.8].

Proof for $ \xi=1 $ {#s.CaseOne}
-------------------

Due to the fact that $ \rho(x) $ is the semicircle law in this special case, now the advection equation is the same as [@bourgade2018extreme equation (1.12)], whose characteristics has an explicit formula. This coincidence makes it easy to adapt the previous proofs for $ \xi \neq 1 $ to this case. With a little abuse of notations, define the curve $${\mathscr{S}}:= \left\{ z=E+{\text{\rm i}}y: 0 <E< 2 - \varphi^2 N^{-2/3},\ y=\varphi^2/\left( N \kappa(E)^{1/2} \right) \right\},$$ where $ \kappa(z) := |z-2| $. Under the framework of such notations and the rigidity estimates and , by the arguments in [@bourgade2018extreme Section 2], all previous results in Section \[s.geometry\] still hold for $ \xi=1 $. Then using the same method we can prove Proposition \[p.estimate\] with few changes, and consequently the proof for Theorem \[t.Relaxation\] is completed.

Rate of Convergence to Tracy-Widom Law {#s.Comparison}
======================================

Quantitative Green function comparison
--------------------------------------

Following the general three-step strategy in the dynamical approach, the derivation of the rate of convergence relies on both the relaxation and the Green function comparison theorem from [@GreenFunctionComparison]. In the context of sample covariance matrices, this Lindeberg exchange strategy based on the fourth moments matching condition was first used by Tao and Vu in [@TaoVu]. To obtain an explicit convergence rate, we need a quantitative version of the comparison theorem.

For the statement, we consider a fixed $ |E-\lambda_+| < \varphi N^{-2/3} $, a scale $ \rho=\rho(N) \in [N^{-1},N^{-2/3}] $, and a function $ f=f(N):{\mathbb{R}}\to {\mathbb{R}}$ satisfying $$\|f^{(k)} \|_{L^\infty([E,E+\rho])} {\leqslant}C_k \rho^{-k},\ \ \ \|f^{(k)} \|_{L^{\infty}([E^+,E^++1])} = O(1),\ \ \ 0 {\leqslant}k {\leqslant}2.$$ where $ E^+ = E+\varphi N^{-2/3} $. We assume that $ f $ is non-decreasing on $ (-\infty,E^+] $, $ f(x) \equiv 0 $ for $ x<E $ and $ f(x) \equiv 1 $ for $ E+\rho<x{\leqslant}E^+ $; and also assume $ f $ is non-increasing on $ [E^+,\infty) $, $ f \equiv 0 $ for $ x>E^+ +1 $. Furthermore, let $ F $ be a fixed smooth non-increasing function such that $ F(x) \equiv 1 $ for $ x {\leqslant}0 $ and $ F(x) \equiv 0 $ for $ x {\geqslant}1 $.

\[t.GreenCompare\] There exists $ C>0 $ such that the following holds. Let $ {X^{\mathrm{v}}}, {X^{\mathrm{w}}}$ be data matrices satisfying assumptions and , and $ {H^{\mathrm{v}}},{H^{\mathrm{w}}}$ be the corresponding sample covariance matrices. Assume that the first three moments of the entries are the same, i.e. for all $ 1 {\leqslant}i {\leqslant}M $, $ 1 {\leqslant}j {\leqslant}N $ and $ 1 {\leqslant}k {\leqslant}3 $ we have $${{\mathbb{E}}^{\mathrm{v}}}(x_{ij}^k) = {{\mathbb{E}}^{\mathrm{w}}}(x_{ij}^k).$$ Assume also that for some parameter $ t=t(N) $ we have $$\left| {{\mathbb{E}}^{\mathrm{v}}}(\sqrt{M}x_{ij})^4 - {{\mathbb{E}}^{\mathrm{w}}}(\sqrt{M}x_{ij})^4 \right| {\leqslant}t.$$ With the above notations for the test functions $ f $ and $ F $, we have $$\left| ({{\mathbb{E}}^{\mathrm{v}}}- {{\mathbb{E}}^{\mathrm{w}}}) F \left( \operatorname{Tr}f(H) \right) \right| {\leqslant}\varphi^C \left( \dfrac{1}{N^{18} \rho^{20}} + \dfrac{t}{N \rho} + \dfrac{1}{(N \rho)^2} + \dfrac{1}{N^2} \right).$$

We follow the notations in [@pillai2014universality] and the reasoning from [@erdos2017dynamical Theorem 17.4]. Fix a bijective ordering map on the index set of the independent matrix elements, $ \phi:\{(i,j):1 {\leqslant}i {\leqslant}M, 1 {\leqslant}j {\leqslant}N\} \to \{1,\cdots,MN\} $ and define the family of random matrices $ X_\gamma $, $ 0 {\leqslant}\gamma {\leqslant}MN $ $$[X_\gamma]_{ij} = \left\{
\begin{aligned}
& [{X^{\mathrm{v}}}]_{ij} & \mbox{if} & \ \  \phi(i,j)>\gamma,\\
& [{X^{\mathrm{w}}}]_{ij} & \mbox{if} & \ \  \phi(i,j) {\leqslant}\gamma.
\end{aligned}
\right.$$ Note that in particular we have $ X_0 = {X^{\mathrm{v}}}$ and $ X_{MN} = {X^{\mathrm{w}}}$. Denote sample covariance matrices $ H_\gamma $ as $$H_\gamma := X_\gamma^* X_\gamma.$$

Let $ \chi $ be a fixed, smooth, symmetric cutoff function such that $ \chi(x)=1 $ if $ |x|<1 $ and $ \chi(x)=0 $ if $ |x|>2 $. By the Helffer-Sjöstrand formula, if $ \lambda_i $’s are the (real) eigenvalues of a matrix $ H $, we have $$\sum f(\lambda_i) = \int_{{\mathbb{C}}} g(z) \operatorname{Tr}\dfrac{1}{H-z} {\mathrm{d}}m(z),$$ where $ {\mathrm{d}}m $ is the Lebesgue measure on $ {\mathbb{C}}$, and the function $ g $ is defined as $$g(z) := \dfrac{1}{\pi} \left( {\text{\rm i}}y f''(y) \chi(y) + {\text{\rm i}}(f(x) + {\text{\rm i}}y f'(x))\chi'(y) \right),\ \ \ z=x+{\text{\rm i}}y.$$ Define $$\Xi^H := \int_{|y|>N^{-1}} g(z) \operatorname{Tr}(H-z)^{-1} {\mathrm{d}}m(z),$$ and we have the bound (see [@bourgade2018extreme Section 5.2]) $$\left|\sum f(\lambda_i) - \Xi^H \right| {\leqslant}O \left( \dfrac{\varphi^C}{(N \rho)^2} \right).$$ This shows that it suffices to show $$\label{e.TelescopicSum}
\left| {\mathbb{E}}F(\Xi^{H_\gamma}) - {\mathbb{E}}F(\Xi^{H_{\gamma-1}}) \right| {\leqslant}\dfrac{ \varphi^C }{N^2} \left( \dfrac{1}{N^{18} \rho^{20}} + \dfrac{t}{N \rho} + \dfrac{1}{(N \rho)^2} + \dfrac{1}{N^2} \right).$$ For an arbitrarily fixed $ \gamma $ corresponding to $ (i,j) $, we can write $$X_{\gamma-1}=Q+V,\ \ \ V:={X^{\mathrm{v}}}_{ij}E^{(ij)},\ \ \ \ X_\gamma = Q+W,\ \ \ W:= {X^{\mathrm{w}}}_{ij}E^{(ij)}.$$ where $ Q $ coincides with $ X_{\gamma-1} $ and $ X_\gamma $ except on the $ (i,j) $ position (where it is 0). We define the Green functions $$R:= (Q^*Q-z)^{-1},\ \ \ S:=(H_{\gamma-1}-z)^{-1}.$$ By Taylor expansion, for some fixed order $ m $, we have $$\begin{gathered}
\label{e.Taylor}
{\mathbb{E}}F(\Xi^{H_{\gamma}}) - {\mathbb{E}}F(\Xi^{H_{\gamma-1}}) = \sum_{l=1}^{m-1} {\mathbb{E}}\dfrac{F^{(l)}(\Xi^Q)}{l !}\left( (\Xi^{H_\gamma} - \Xi^Q)^l - (\Xi^{H_{\gamma-1}} - \Xi^Q)^l \right)\\
 + O \left( \|F^{(m)} \|_\infty \right) \left( {\mathbb{E}}\left( (\Xi^{H_\gamma} - \Xi^Q)^m + (\Xi^{H_{\gamma-1}} - \Xi^Q)^m \right) \right).\end{gathered}$$

First we estimate the $ m $-th order error term. By the first order resolvent expansion we have $$\begin{aligned}
|\Xi^{H_{\gamma}} - \Xi^Q| &{\leqslant}\int_{|y|>N^{-1},|x|<\lambda_+ +2} |g(z)|\left| \operatorname{Tr}R(z)(V^*Q+Q^*V+V^*V)S(z) \right| {\mathrm{d}}m(z) \\
&{\leqslant}\varphi^C N \int_{|y|>N^{-1},|x|<\lambda_+ +2} |g(z)| \|S(z) \|_\infty \|R(z) \|_\infty {\mathrm{d}}m(z)\end{aligned}$$ with overwhelming probability, where we use the fact that there are only $ O(N) $ nonzero entries with size $ O(N^{-1}) $ in the matrix $ V^*Q+Q^*V+V^*V $. By the strong local Marchenko-Pastur law ([@pillai2014universality Theorem 3.1]), for any $ D>0 $ we have $${\mathbb{P}}\left( \max_{j}|S_{jj}(z) - m_{{\mathsf{MP}}}(z)| + \max_{j \neq k}|S_{jk}(z)| {\leqslant}\varphi^C \left( \dfrac{1}{Ny}+\sqrt{\dfrac{{\text{\rm Im}\hspace{0.1cm}}m_{{\mathsf{MP}}}(y)}{Ny}} \right) \right) > 1- N^{-D}.$$ Same bound for $ \|R(z) \|_\infty $ also holds (see [@pillai2014universality Lemma 5.4]). This shows that $${\mathbb{E}}(\Xi^{H_\gamma} - \Xi^Q)^m = O \left( \varphi^C / (N^m \rho^m) \right),\ \ \ {\mathbb{E}}(\Xi^{H_{\gamma-1}} - \Xi^Q)^m = O \left( \varphi^C / (N^m \rho^m) \right).$$ Therefore the $ m $-th order term in can be bounded by $ \varphi^C N^{-2}(N^{-m+2} \rho^{-m}) $.

Next we consider the first order term in the Taylor expansion. By the resolvent expansion, we have $$S=R-R {A^{\mathrm{v}}}R+(R {A^{\mathrm{v}}})^2R-(R {A^{\mathrm{v}}})^3R+ \cdots - (R {A^{\mathrm{v}}})^{11}R + (R {A^{\mathrm{v}}})^{12}S,$$ where $${A^{\mathrm{v}}}= V^*Q + Q^*V + V^*V.$$ Denote $${{\widehat}{R}_{\mathrm{v}}}^{(n)} := (-1)^{n} \operatorname{Tr}(R {A^{\mathrm{v}}})^n R,\ \ \ \ {\Omega_{\mathrm{v}}}:= \operatorname{Tr}(R {A^{\mathrm{v}}})^{12} S.$$ Then we have $${\mathbb{E}}F'(\Xi^Q) \left( \Xi^{H_{\gamma-1}} - \Xi^{H_{\gamma}} \right) =  {\mathbb{E}}F'(\Xi^Q) \int g(z) \left( \sum_{n=1}^{11} \left( {{\widehat}{R}_{\mathrm{v}}}^{(n)} - {{\widehat}{R}_{\mathrm{w}}}^{(n)} \right) +({\Omega_{\mathrm{v}}}- {\Omega_{\mathrm{w}}}) \right) {\mathrm{d}}m(z).$$ Since the first three moments of the two matrices are identical, we know that the case $ n=1 $ gives null contribution.

For $ n=2 $, note that the entries of the matrix $ A $ satisfy the following relation $$A_{ab}=\left\{
\begin{aligned}
& x_{ij} x_{ib} & \quad & \mbox{if} \  a=j,b \neq j,\\
& x_{ij} x_{ia} & \quad & \mbox{if} \  a \neq j, b=j,\\
& 0 & \quad &\mbox{otherwise}.
\end{aligned}
\right.$$ This shows that $${\mathbb{E}}({{\widehat}{R}_{\mathrm{v}}}^{(2)} - {{\widehat}{R}_{\mathrm{w}}}^{(2)}) {\leqslant}N \left( \dfrac{t}{N^2} \right) \left(\max_{i \neq j} |R_{ij}|\right)^2 \left(\max_i |R_{ii}|\right),$$ where we used that in the expansion $$\begin{aligned}
\operatorname{Tr}(R {A^{\mathrm{v}}})^2R & = \sum_{k} \sum_{a_1,b_1,a_2,b_2} R_{k a_1}{A^{\mathrm{v}}}_{a_1 b_1}R_{b_1 a_2}{A^{\mathrm{v}}}_{a_2 b_2}R_{b_2 k}\\
& = \sum_{k} \sum_{(a_1,b_1,a_2,b_2) \neq (j,j,j,j)} R_{k a_1}{A^{\mathrm{v}}}_{a_1 b_1}R_{b_1 a_2}{A^{\mathrm{v}}}_{a_2 b_2}R_{b_2 k}  + \sum_k  R_{kj}{A^{\mathrm{v}}}_{jj}R_{jj}{A^{\mathrm{v}}}_{jj}R_{jk}\end{aligned}$$ due to the moment matching condition, the terms that make nontrivial contribution are only in the second summation, which is $$\sum_k R_{kj}R_{jj}R_{jk}({X^{\mathrm{v}}}_{ij})^4.$$ Here we also use the fact that the contribution for the terms with $ k $ equal $ i $ or $ j $ is combinatorially negligible. By the local law, we conclude $${\mathbb{E}}F'(\Xi^Q) \int g(z) \left( {{\widehat}{R}_{\mathrm{v}}}^{(2)} - {{\widehat}{R}_{\mathrm{w}}}^{(2)} \right) {\mathrm{d}}m(z) = O \left( \dfrac{\varphi^C t}{N} \right) \int \dfrac{|g(z)|}{(Ny)^2} {\mathrm{d}}m(z) = O \left( \dfrac{\varphi^C}{N^2} \dfrac{t}{N\rho}  \right).$$ For the terms $ n=3,\dots,11 $, as explained in [@pillai2014universality Lemma 5.4], their contributions are of smaller order. Similarly, for the term $ ({\Omega_{\mathrm{v}}}- {\Omega_{\mathrm{w}}}) $, as shown in [@pillai2014universality Lemma 5.4] we have $ {\Omega_{\mathrm{v}}}= O(N^{-4}) $. Therefore we have $${\mathbb{E}}F'(\Xi^Q) \int g(z) \left( ({\Omega_{\mathrm{v}}}- {\Omega_{\mathrm{w}}}) \right) {\mathrm{d}}m(z) = O \left( \dfrac{\varphi^C}{N^2} \dfrac{1}{N^2 \rho} \right)=O \left( \dfrac{\varphi^C}{N^2} \right) \left( \dfrac{1}{(N \rho)^2} + \dfrac{1}{N^2} \right).$$

Moreover, as explained in [@erdos2017dynamical Theorem 17.4], the contributions of higher order terms in Taylor expansion are of smaller order. Combining the estimates and taking $ m=20 $ (we will see the reason in the next section) gives us . Finally, a telescopic summation yields the desired result.

Proof of Theorem \[t.Rate\]
---------------------------

Let $ s \in {\mathbb{R}}$. If $ |s|>\varphi $, due to the rigidity we know that for any $ D>0 $ and large enough $ N $, we have $ {\mathbb{P}}\left( N^{2/3}(\lambda_N - \lambda_+) {\leqslant}s \right) = {\mathbb{P}}({\mathsf{TW}}{\leqslant}s) + O(N^{-D}) $. So in the following discussion we assume $ |s| {\leqslant}\varphi $.

Denoting a non-decreasing function $ f_1 $ such that $ f_1(x)=1 $ for $ x>\lambda_+ + sN^{-2/3} $ and $ f_1(x)=0 $ for $ x<\lambda_+ + sN^{-2/3} - \rho $. We also define $ f_2(x) := f_1(x-\rho) $. Then we have $$\label{e.ProbEst}
{\mathbb{E}}_H F \left( \sum_{i=1}^N f_1(\lambda_i)  \right) {\leqslant}{\mathbb{P}}_H \left( \lambda_N < \lambda_+ + s N^{-2/3} \right) {\leqslant}{\mathbb{E}}_H F \left( \sum_{i=1}^N f_2(\lambda_i)  \right).$$

Moreover, as discussed in [@erdos2011universality; @pillai2014universality], we can find an $ M \times N $ matrix $ {\widetilde}{X}_0 $ such that the Gaussian divisible ensemble $ {\widetilde}{X}_t := e^{-t/2}{\widetilde}{X}_0 + (1-e^{-t})^{1/2}X_G $, where $ X_G $ is a matrix whose entries are independent Gaussian random variables with mean 0 and variance 1, satisfies the following: for $ 1 {\leqslant}k {\leqslant}3 $, $${\mathbb{E}}(\sqrt{M}X_{ij})^k = {\mathbb{E}}[{\widetilde}{X}_t]_{ij}^k,\ \ \ \ |{\mathbb{E}}(\sqrt{M}X_{ij})^4 - {\mathbb{E}}[{\widetilde}{X}_t]_{ij}^4| \lesssim t.$$ By the quantitative Green function comparison theorem \[t.GreenCompare\], we obtain the following bound $$\begin{gathered}
{\mathbb{E}}_{{\widetilde}{X}_t} F \left( \sum_{i=1}^N f_1(\lambda_i) \right) - \varphi^C \left( \dfrac{1}{N^{18} \rho^{20}} + \dfrac{t}{N \rho} + \dfrac{1}{(N \rho)^2} + \dfrac{1}{N^2} \right)
{\leqslant}{\mathbb{P}}_H \left( \lambda_N < \lambda_+ + s N^{-2/3} \right) \\
{\leqslant}{\mathbb{E}}_{{\widetilde}{X}_t} F \left( \sum_{i=1}^N f_2(\lambda_i) \right) + \varphi^C \left( \dfrac{1}{N^{18} \rho^{20}} + \dfrac{t}{N \rho} + \dfrac{1}{(N \rho)^2} + \dfrac{1}{N^2} \right).\end{gathered}$$ Using for $ {\widetilde}{X}_t $, the estimate becomes $$\begin{gathered}
{\mathbb{P}}_{{\widetilde}{X}_t} \left( \lambda_N < \lambda_+ + s N^{-2/3} -\rho \right)  - \varphi^C \left( \dfrac{1}{N^{18} \rho^{20}} + \dfrac{t}{N \rho} + \dfrac{1}{(N \rho)^2} + \dfrac{1}{N^2} \right)\\
{\leqslant}{\mathbb{P}}_H \left( \lambda_N < \lambda_+ + s N^{-2/3} \right) {\leqslant}\\
{\mathbb{P}}_{{\widetilde}{X}_t} \left( \lambda_N < \lambda_+ + s N^{-2/3} +\rho \right)  + \varphi^C \left( \dfrac{1}{N^{18} \rho^{20}} + \dfrac{t}{N \rho} + \dfrac{1}{(N \rho)^2} + \dfrac{1}{N^2} \right).\end{gathered}$$ After combined with the edge relaxation Theorem \[t.Relaxation\], the estimate now gives us $$\begin{gathered}
{\mathbb{P}}_{{\mathsf{Wishart}}} \left( N^{2/3}(\lambda_N - \lambda_+) < s -N^{2/3}\rho - \dfrac{N^{\varepsilon}}{N^{1/3}t} \right) - \varphi^C \left( \dfrac{1}{N^{18} \rho^{20}} + \dfrac{t}{N \rho} + \dfrac{1}{(N \rho)^2} + \dfrac{1}{N^2} \right)\\
 {\leqslant}{\mathbb{P}}_H \left( N^{2/3}(\lambda_N - \lambda_+) < s \right) {\leqslant}\\
{\mathbb{P}}_{{\mathsf{Wishart}}} \left( N^{2/3}(\lambda_N - \lambda_+) < s +N^{2/3}\rho + \dfrac{N^{\varepsilon}}{N^{1/3}t} \right) + \varphi^C \left( \dfrac{1}{N^{18} \rho^{20}} + \dfrac{t}{N \rho} + \dfrac{1}{(N \rho)^2} + \dfrac{1}{N^2} \right).\end{gathered}$$ Moreover, as shown in [@el2006rate; @ma2012accuracy], we know $${\mathbb{P}}_{{\mathsf{Wishart}}} \left( N^{2/3}(\lambda_N - \lambda_+) <s \right)  = {\mathbb{P}}({\mathsf{TW}}<s)+ O(N^{-2/3}).$$ By using this Wishart result and the boundedness of the density for $ {\mathsf{TW}}$, we obtain $$\begin{gathered}
\label{e.RateNullCase}
{\mathbb{P}}_H \left(N^{2/3}(\lambda_N - \lambda_+) <s \right) - {\mathbb{P}}\left( {\mathsf{TW}}<s \right)\\
 = O \left( N^{\varepsilon}\right) \left( N^{2/3}\rho + \dfrac{1}{N^{1/3}t} + \dfrac{1}{N^{18} \rho^{20}} + \dfrac{t}{N \rho} + \dfrac{1}{(N \rho)^2} + \dfrac{1}{N^{2/3}} \right).\end{gathered}$$ The optimal bound $ N^{-2/9 + {\varepsilon}} $ is obtained for $ t=N^{-1/9} $ and $ \rho = N^{-8/9} $. This completes the whole proof for Theorem \[t.Rate\].

Generalization to General Population Matrices {#s.GeneralPopulation}
=============================================

In this section, we proceed to generalize our previous results for sample covariance matrices of type $ X^*X $ (which corresponds to the identity population) and aim to derive the rate of convergence to the Tracy-Widom distribution for the (rescaled) largest eigenvalue of separable sample covariance matrices with general population. Throughout this section, we will follow the notations and the setup in the the work by Lee and Schnelli [@lee2016tracy].

Let $ X=(x_{ij}) $ be defined as in and . For some deterministic $ M \times M $ matrix $ T $, the sample covariance matrices associated with data matrix $ X $ and population matrix $ \Sigma:=T^*T $ is defined as $ \mathcal{Q}:=(TX)(TX)^* $. Note that the $ M \times M $ matrix $ \mathcal{Q} $ and the matrix $$Q:=X^*\Sigma X$$ share the same non-trivial eigenvalues. Since we are studying the largest eigenvalue of the sample covariance matrix, it is more convenient to work with the matrix $ Q $ (called the separable sample covariance matrix) for some technical reasons. We denote the eigenvalues of $ Q $ in increasing order by $ \mu_1 {\leqslant}\cdots {\leqslant}\mu_N $.

As mentioned previously, for the null case (i.e. the population matrix is identity), it is well known that the empirical eigenvalue distribution of a sample covariance matrix converges weakly in probability to the Marchenko-Pastur law. Under the general setting, however, this results need to be modified and the limiting measure (called the deformed Marchenko-Pastur law) will depends on the spectrum of the population matrix. Let $ \sigma_1 {\leqslant}\cdots {\leqslant}\sigma_M $ be the eigenvalues of the population matrix $ \Sigma $, we denote by $ {\widehat}{\rho}={\widehat}{\rho}(M) $ the empirical eigenvalue distribution of $ \Sigma $, which is defined as $${\widehat}{\rho}:=\dfrac{1}{M}\sum_{j=1}^M \delta_{\sigma_j}.$$ The deformed Marchenko-Pastur law $ {\widehat}{\rho}_{{\mathsf{fc}}} $ is defined in the following way. The Stieltjes transform $ {\widehat}{m}_{{\mathsf{fc}}} $ of the probability measure is given by the unique solution of the equation $${\widehat}{m}_{{\mathsf{fc}}}(z)=\dfrac{1}{-z+\xi^{-1}\int \frac{1}{t  {\widehat}{m}_{{\mathsf{fc}}}(z) + 1}{\mathrm{d}}{\widehat}{\rho}(t)},\ \ \ {\text{\rm Im}\hspace{0.1cm}}{\widehat}{m}_{{\mathsf{fc}}}(z) {\geqslant}0,\ \ \ z \in {\mathbb{C}}^+.$$ It has been discussed in [@knowles2017anisotropic] that $  {\widehat}{m}_{{\mathsf{fc}}} $ is associated to a continuous probability density $  {\widehat}{\rho}_{{\mathsf{fc}}} $ with compact support in $ [0,\infty) $. Moreover, the density ${\widehat}{\rho}_{{\mathsf{fc}}}  $ can be obtained from $ {\widehat}{m}_{{\mathsf{fc}}} $ via the Stieltjes inversion formula $${\widehat}{\rho}_{{\mathsf{fc}}}(E) = \lim_{\eta \downarrow 0}\dfrac{1}{\pi}{\text{\rm Im}\hspace{0.1cm}}{\widehat}{m}_{{\mathsf{fc}}}(E+{\text{\rm i}}\eta).$$ The typical location of the largest eigenvalue, which is the rightmost endpoint of the support of the density $ {\widehat}{\rho}_{{\mathsf{fc}}} $ is determined in the following way. Recall that $ \xi:=N/M $, we define $ \xi_+ $ as the largest solution of the equation $$\int \left( \dfrac{t \xi_+}{1-t\xi_+} \right)^2 {\mathrm{d}}{\widehat}{\rho}_{{\mathsf{fc}}}(t) = \xi.$$ We remark that $ \xi_+ $ is unique and $ \xi_+ \in [0,\sigma_M^{-1}] $. We then introduce the typical location for the largest eigenvalue $ E_+ $ by $$\label{e.RightEndpoint}
E_+ := \dfrac{1}{\xi_+}\left( 1+\xi^{-1}\int \dfrac{t \xi_+}{1-t\xi_+} {\mathrm{d}}{\widehat}{\rho}_{{\mathsf{fc}}}(t) \right).$$

Now we state our assumptions on the population matrix $ \Sigma $ that are needed to prove the explicit rate of convergence. For general random matrices $ X $ we require $ \Sigma $ to be diagonal, and we will show later this diagonal condition can be removed if $ X $ is Gaussian. We further need the following assumption for the spectrum of the population matrix $ \Sigma $. Throughout this section, we assume the following: $$\label{e.AssumptionPopulation}
\liminf_M \sigma_1 >0,\ \ \ \limsup_M \sigma_M <\infty,\ \ \ \mbox{and}\ \ \ \limsup_M \sigma_M \xi_+ <1.$$

The assumption is the same as [@lee2016tracy Assumption 2.2]. It is first used in [@bao2015universality; @knowles2017anisotropic] to prove the local deformed Marchenko-Pastur law. In particular, the last inequality ensures that the density $ {\widehat}{\rho}_{{\mathsf{fc}}} $ exhibits a square-root behavior near the right edge of its support, which is crucial to derive the local law.

It is natural to note that with a general population matrix, the distribution of the largest eigenvalue should not behave exactly like the null case. Besides the typical location of the largest eigenvalue is changed, the normalization constant of the fluctuation is also different. Therefore, we introduce the following normalization constant $ \gamma_0 $ given by $$\label{e.ScalingConstant}
\dfrac{1}{\gamma_0^3} = \dfrac{1}{\xi}\int \left( \dfrac{t}{1-t\xi_+} \right){\mathrm{d}}{\widehat}{\rho}(t) + \dfrac{1}{\xi_+^3}.$$ Moreover, we remark that the Tracy-Widom limit for the general case is rescaled and it is not the same as the previous one we used for the null case. However they are just different by a simple scaling so that we do not emphasize this difference and still use the notation $ {\mathsf{TW}}$ to denote this distribution. Under this framework, our main result given in Corollary \[t.RateGeneralPopulation\].

Unlike the proof for the null case in the previous sections, we will not strictly follow the three-step strategy in the dynamical approach. Instead, we use the comparison theorem for the Green function flow, which is a method based on continuous interpolation, linearization and renormalization developed in [@lee2016tracy].

Local deformed Marchenko-Pastur law
-----------------------------------

For completeness, we will briefly introduce the local deformed Marchenko-Pastur law. Though we will not give a detailed proof, we emphasize that the local law is an indispensable part to prove the Green function comparison Proposition \[p.GreenFunctionFlow\], which further leads to the edge universality.

For small nonnegative $ c,\epsilon {\geqslant}0 $ and sufficiently large $ E_+<C<\infty $, we consider the domain $$\mathcal{D}(c,\epsilon) :=\left\{ z=E+{\text{\rm i}}\eta \in {\mathbb{C}}^+ : E_+-c {\leqslant}E {\leqslant}C,  N^{-1+\epsilon} {\leqslant}\eta {\leqslant}1 \right\}.$$ We also denote $ \kappa=\kappa(E):=|E-E_+| $. Then we have the following estimates for the density and the Stieltjes transform of the deformed Marchenko-Pastur law.

\[l.EstiamtesDeformedMP\] Under the assumption , there exists a constant $ c>0 $ such that $${\widehat}{\rho}_{\mathsf{fc}}(E) \sim \sqrt{E_+ - E}\ \ \ \ E \in [E_+-2c,E_+].$$ Moreover, the Stieltjes transform $ {\widehat}{m}_{{\mathsf{fc}}} $ satisfies the following: for $ z \in \mathcal{D}(c,0) $, we have $$|{\widehat}{m}_{{\mathsf{fc}}}(z)| \sim 1,\ \ \ \ {\text{\rm Im}\hspace{0.1cm}}{\widehat}{m}_{{\mathsf{fc}}}(z) \sim \left\{
\begin{aligned}
& \frac{\eta}{\sqrt{\kappa+\eta}} & \mbox{if} &\ E {\geqslant}E_++\eta,\\
& \sqrt{\kappa+\eta}, & \mbox{if} &\ E \in [E_+-c,E_++\eta).
\end{aligned}
\right.
.$$

The Green function and the Stieltjes transform are defined in the usual way: $$G_Q(z) := (Q-z)^{-1},\ \ \ m_Q(z) := \dfrac{1}{N}\operatorname{Tr}G_Q(z).$$ Then we have the following local law for the separable sample covariance matrix $ Q $.

\[l.LocalDeformedMPLaw\] Under the assumption , for any sufficiently small $ \epsilon>0 $, and for any (large) $ D>0 $, there exists $ N_0(D)>0 $ such that for any $ N {\geqslant}N_0(D) $ we have the following estimate uniformly in $ z \in \mathcal{D}(c,\epsilon) $: $${\mathbb{P}}\left( |m_Q(z) - {\widehat}{m}_{{\mathsf{fc}}}(z)| {\leqslant}\dfrac{N^\epsilon}{N\eta} \right) > 1-N^{-D},$$ and $${\mathbb{P}}\left( \max_{i,j}\left| (G_Q)_{ij}(z) - \delta_{ij}{\widehat}{m}_{{\mathsf{fc}}}(z) \right| {\leqslant}N^\epsilon \left( \sqrt{\dfrac{{\text{\rm Im}\hspace{0.1cm}}{\widehat}{m}_{{\mathsf{fc}}}(z)}{N\eta}} + \dfrac{1}{N\eta} \right) \right) > 1-N^{-D}.$$

It is clear to see that the estimates for the deformed Marchenko-Pastur law (Lemma \[l.EstiamtesDeformedMP\]) and the local law (Lemma \[l.LocalDeformedMPLaw\]) are greatly similar as the corresponding results for the null case (see e.g. [@pillai2014universality Theorem 3.1]). Heuristically, this implies the Tracy-Widom limit in the edge universality for the non-null case.

Interpolation and Green function comparison
-------------------------------------------

In classical theory of random matrix universality, the tool needed to prove the edge universality is the Green function comparison theorem. The usual approach is to compare two ensembles with some moments matching conditions, and then use the construction of Gaussian divisible ensembles together with estimates of the local relaxation flow to remove the moments matching requirement. In this section, however, we do not follow this traditional step. Instead, we compare the Green function of a general ensemble with its corresponding null sample covariance matrix. This argument was first introduced in [@lee2015edge] to handle the deformed Wigner matrices, and used in [@lee2016tracy] to identify the Tracy-Widom limit for general separable sample covariance matrices. The basic idea is to introduce a time evolution that deforms the population matrices continuously to the identity and offset the change of the Green function by a renormalization of the matrix.

Recall the scaling constant $ \gamma_0 $ defined in . We consider the following two rescaled matrices $${\widetilde}{\Sigma} := \gamma_0 \Sigma,\ \ \ {\widetilde}{Q} := X^* {\widetilde}{\Sigma}X.$$ We also denote the eigenvalues of $ {\widetilde}{Q} $ by $ {\widetilde}{\mu}_1 {\leqslant}\cdots {\leqslant}{\widetilde}{\mu}_N $, and let $ L_+ := \gamma_0 E_+ $. We remark that in the literature about sample covariance matrices with general population (e.g. [@bao2015universality; @el2007tracy; @lee2016tracy]), the scaling of the Tracy-Widom distribution is chosen in the way such that it is the limit for the distribution of the (rescaled) largest eigenvalue of the matrix $$W:=\sqrt{\xi}(1+\sqrt{\xi})^{-4/3}X^*X.$$ Specifically, we order the eigenvalues of the matrix $ W $ by $ \lambda_1 {\leqslant}\cdots {\leqslant}\lambda_N $, and let $ M_+ $ denote the rightmost endpoint of the rescaled Marchenko-Pastur law for $ W $.

It has been shown in [@bao2015universality equation (1.9)] that $$\gamma_0 = \sqrt{\xi}(1+\sqrt{\xi})^{-4/3} + o(1).$$ This can be regarded as a good motivation for considering the scaling constant $ \gamma_0 $.

For the diagonal population matrix $ \Sigma = \mbox{diag}(\sigma_j) $, we introduce the following time evolution $ t \mapsto (\sigma_j(t)) $ that deforms $ \Sigma $ to the identity matrix $ {\mathsf{Id}}$ and the Green function flow by $$\label{e.GreenFunctionFlow}
\dfrac{1}{\sigma_j(t)}=e^{-t}\dfrac{1}{\sigma_j(0)}+(1-e^{-t}),\ \ \ {\widetilde}{Q}(t)=\gamma_0X^*\Sigma(t)X,\ \ \ m_{{\widetilde}{Q}(t)}(z) := \dfrac{1}{N}\operatorname{Tr}({\widetilde}{Q}(t)-z)^{-1}.$$ Based on the local law Lemma \[l.EstiamtesDeformedMP\] and Lemma \[l.LocalDeformedMPLaw\], and a delicate analysis for the time derivative of the Green function for $ {\widetilde}{Q}(t) $, the Green function comparison theorem (see Proposition \[p.GreenFunctionFlow\]) is proved in [@lee2016tracy]. We note that though the original estimate in [@lee2016tracy] is not explicit, a careful examination of the proof will reveal that the result is actually quantitative.

\[p.GreenFunctionFlow\] Let $ {\varepsilon}>0 $ and set $ \eta = N^{-2/3-{\varepsilon}} $. Let $ E_1,E_2 \in {\mathbb{R}}$ satisfy $ E_1 < E_2 $ and $ |E_1|,|E_2| {\leqslant}N^{-2/3+{\varepsilon}} $. Let $ F: {\mathbb{R}}\to {\mathbb{R}}$ be a smooth function satisfying $$\max_x |F^{(l)}(x)|(|x|+1)^{-C} {\leqslant}C,\ \ \ \ l=1,2,3,4.$$ Then for any (small) $ {\delta}>0 $ and for sufficiently large $ N $ we have $$\begin{gathered}
\left| {\mathbb{E}}F \left( N \int_{E_1}^{E_2} {\text{\rm Im}\hspace{0.1cm}}m_{{\widetilde}{Q}}(x+L_++{\text{\rm i}}\eta) {\mathrm{d}}x \right) - {\mathbb{E}}F \left( N \int_{E_1}^{E_2} {\text{\rm Im}\hspace{0.1cm}}m_{W}(x+M_++{\text{\rm i}}\eta) {\mathrm{d}}x \right) \right|\\
{\leqslant}N^{-\frac{1}{3}+2{\varepsilon}+{\delta}}.\end{gathered}$$

Quantitative edge universality
------------------------------

In this section we can finally prove Corollary \[t.RateGeneralPopulation\]. The proof is based on our previous rate of convergence for the null case (Theorem \[t.Rate\]) and the estimate on the comparison theorem for the Green function flow (Proposition \[p.GreenFunctionFlow\]).

We first remark that we are not supposed to use the rate of convergence for the null case (Theorem \[t.Rate\]) and the triangle inequality in a naive way to derive the convergence rate for the general case. This is because Theorem \[t.Rate\] is obtain by choosing the optimal parameters in the estimate , and the scale parameter $ \rho $ is also related to the scale in the Green function comparison (Proposition \[p.GreenFunctionFlow\]). To illustrate this link more clearly, we first briefly review how Green function comparison is used to obtain the edge universality.

We introduce a smooth cutoff function $ K:{\mathbb{R}}\to {\mathbb{R}}$ satisfying $$K(x)=
\left\{
\begin{aligned}
& 1 &\ \mbox{if}\ & x {\leqslant}1/9,\\
& 0 &\ \mbox{if}\ & x {\geqslant}2/9.
\end{aligned}
\right.$$ and we also define the Poisson kernel $ \theta_\eta $, for $ \eta>0 $ $$\theta_\eta(x) := \dfrac{\eta}{\pi(x^2+\eta^2)}.$$ Let $ E_*:=L_+ +\varphi^C N^{-2/3} $, and denote $ \chi_E := 1_{[E,E_*]} $. For $ {\varepsilon}>0 $, let $ l:= \tfrac{1}{2}N^{-2/3-{\varepsilon}} $ and $ \eta:=N^{-2/3-9{\varepsilon}} $. Then for any (large) $ D>0 $, it is proved in [@lee2016tracy; @pillai2014universality] that for large enough $ N $ we have $${\mathbb{E}}K\left( \operatorname{Tr}(\chi_{E-l} * \theta_\eta({\widetilde}{Q}))  \right) {\leqslant}{\mathbb{P}}\left( {\widetilde}{\mu}_N {\leqslant}E \right) {\leqslant}{\mathbb{E}}K\left( \operatorname{Tr}(\chi_{E+l} * \theta_\eta({\widetilde}{Q}))  \right) + N^{-D}.$$ Here the parameter $ l $ plays the same role as the $ \rho $ in , and therefore we have $ N^{-{\varepsilon}}=N^{2/3}\rho $ and $ \eta=N^{-2/3}N^{6}\rho^{9} $. By the Green function comparison Proposition \[p.GreenFunctionFlow\] and the [@lee2016tracy Theorem 2.4], we have $$\begin{gathered}
{\mathbb{P}}\left( N^{2/3}(\lambda_N - M_+) {\leqslant}s \right) - N^{-\frac{1}{3}+18{\varepsilon}+{\delta}} {\leqslant}{\mathbb{P}}\left( N^{2/3}({\widetilde}{\mu}_N - L_+) {\leqslant}s \right)\\
{\leqslant}{\mathbb{P}}\left( N^{2/3}(\lambda_N - M_+) {\leqslant}s \right) + N^{-\frac{1}{3}+18{\varepsilon}+{\delta}}.\end{gathered}$$ This gives us $$\label{e.GeneralQuantEdgeUniversality}
{\mathsf{d_K}}\left( \gamma_0 N^{2/3}(\mu_N - E_+),N^{2/3}(\lambda_N-M_+) \right) {\leqslant}N^{{\delta}}N^{-1/3}N^{-12}\rho^{-18}.$$ Combined with Theorem , by triangle inequality we finally obtain $$\begin{aligned}
&\quad {\mathsf{d_K}}\left( \gamma_0 N^{2/3}(\mu_N - E_+),{\mathsf{TW}}\right)\\
& {\leqslant}{\mathsf{d_K}}\left( \gamma_0 N^{2/3}(\mu_N - E_+),N^{2/3}(\lambda_N-M_+) \right) + {\mathsf{d_K}}\left( N^{2/3}(\lambda_N-M_+),{\mathsf{TW}}\right)\\
& {\leqslant}N^{{\delta}}\left( N^{-1/3}N^{-12}\rho^{-18} + N^{2/3}\rho + \dfrac{1}{N^{1/3}t} + \dfrac{1}{N^{18} \rho^{20}} + \dfrac{t}{N \rho} + \dfrac{1}{(N \rho)^2} + \dfrac{1}{N^{2/3}} \right).\end{aligned}$$ The optimal result is obtained now by choosing $ \rho=N^{-13/19} $ and $ t=N^{-6/19} $, which gives us $${\mathsf{d_K}}\left( \gamma_0 N^{2/3}(\mu_N - E_+),{\mathsf{TW}}\right) {\leqslant}N^{-\frac{1}{57}+{\delta}}.$$ This completes the proof.

Based on the Corollary \[t.RateGeneralPopulation\] for diagonal population matrices $ \Sigma $ and general random matrices $ X $, we can easily obtain the following result for the case in which we can have a general population if the random matrix $ X $ is restricted to be Gaussian.

\[c.RateGaussian\] Let $ Q := X^* \Sigma X $ be an $ N \times N $ separable sample covariance matrix, where $ X $ is an $ M \times N $ real random matrix with independent Gaussian entries satisfying , and $ \Sigma $ is a real positive-definite deterministic $ M \times M $ matrix satisfying . For any $ {\varepsilon}>0 $ and large enough $ N $, we have $${\mathsf{d_K}}\left( \gamma_0 N^{2/3} (\mu_N - E_+),{\mathsf{TW}}\right) {\leqslant}N^{-\frac{1}{57} + {\varepsilon}}$$

Under these assumptions, we know that the population matrix $ \Sigma $ is diagonizable, i.e. there exists an $ M \times M $ real diagonal matrix $ D $ and an $ N \times N $ orthogonal matrix $ U $ such that $ \Sigma = U^*DU $. Since $ X $ is a matrix whose entries are independent Gaussian random variable, we know that $ UX $ is also a real random matrix with Gaussian entries satisfying the assumption . Therefore, by applying our Corollary \[t.RateGeneralPopulation\] to the matrix $ X^*\Sigma X = (UX)^*D(UX) $, we will get the desired result.
