Central limit theorem for mesoscopic eigenvalue statistics of deformed Wigner matrices and sample covariance matrices
Yiting Li Kevin Schnelli Yuanyuan Xu






 0.85
 

 


Central limit theorem for mesoscopic eigenvalue statistics of deformed Wigner matrices and sample covariance matrices







	



 1.1

 0.33


Yiting Li  

		KTH Royal Institute of Technology

yitingl@kth.se



0.33
	

		Kevin Schnelli[3]

		KTH Royal Institute of Technology

		schnelli@kth.se
	


0.33


Yuanyuan Xu[2]

KTH Royal Institute of Technology

yuax@kth.se







[2]Supported by  the Göran Gustafsson Foundation and the Swedish Research Council Grant VR-2017-05195.
[3]Supported in parts by the Swedish Research Council Grant VR-2017-05195.


	






 0.83
 Abstract.
 We consider N by N deformed Wigner random matrices of the form X_N=H_N+A_N, where H_N is a real symmetric or complex Hermitian Wigner matrix and A_N is a deterministic real bounded diagonal matrix. We prove a universal Central Limit Theorem for the linear eigenvalue statistics of X_N for all mesoscopic scales both in the spectral bulk and at regular edges where the global eigenvalue density vanishes as a square root. The method relies on  studying the characteristic function of the linear statistics <cit.> by using the cumulant expansion method, along with local laws for the Green function of X_N <cit.> and analytic subordination properties of the free additive convolution <cit.>. We also prove the analogous results for high-dimensional sample covariance matrices.
 






 


 0.83
 Résumé.
 Nous  considérons des matrices aléatoires de Wigner déformées de taille N de la forme X_N = H_N + A_N,
où H_N est une matrice hermitienne de Wigner symétrique ou complexe réelle, et A_N est une 
matrice diagonale déterministe avec des entrées réelles et bornées. Nous prouvons un théorème de limite centrale universel pour les statistiques linéaires des valeurs propres de X_N pour toutes les échelles mésoscopiques à la fois dans le centre de spectre et aux bords réguliers où
la densité globale des valeurs propres disparait sous forme de racine carrée. La méthode repose sur l'étude de la fonction caractéristique des statistiques linéaires <cit.> en utilisant la méthode des cumulants, ainsi que les lois locales
pour la fonction de Green de X_N <cit.> et les propriétés de subordination analytique de la convolution libre additive <cit.>. Nous prouvons également les résultats analogues pour des matrices de corrélation empirique de haute dimension.
 



 
 

 

 
 Date: December 30, 2023

Keywords: linear eigenvalue statistics; deformed Wigner matrices; sample covariance matrices.


 
 


 


headings



§ INTRODUCTION

 



 §.§ Linear eigenvalue statistics of Wigner matrices



A Wigner matrix H_N is an N × N real symmetric or complex Hermitian random matrix with independent entries up to the constraint H_N=H_N^*. In the case the entries are Gaussian random variables, these matrices belong to the Gaussian Orthogonal Ensemble (GOE), Gaussian Unitary Ensemble (GUE), respectively. Wigner <cit.> proved the semicircle law stating that the empirical eigenvalue distribution of H_N converges to the semicircle distribution with density ρ_sc(x)=1/2 π√(4-x^2)1_[-2,2]. That is, for any test function f ∈ C_c(),

    1/N∑_i=1^Nf(λ_i) →∫_f(x) ρ_sc(x)  x  N→∞,

in probability, which can be understood as a Law of Large Numbers.


Johansson <cit.> obtained the corresponding Central Limit Theorem (CLT) for such linear eigenvalue statistics of the GUE, i.e.,

    ∑_i=1^Nf(λ_i)- N ∫_ f(x)ρ_sc(x) x

converges in distribution to a centered Gaussian random variable. Strikingly different from the classical CLT, the linear statistics need not be normalized by N^-1/2, which is a manifestation of the strong eigenvalue correlations.  Bai and Yao <cit.> used a martingale method to prove such CLTs for Wigner matrices and analytic test functions. Lytova and Pastur <cit.>, and Shcherbina <cit.> improved these results by weakening the regularity conditions on the test functions. More recently, Sosoe and Wong <cit.> obtained CLTs for Wigner matrices with H^1+ϵ test functions using Littlewood–Paley decompositions.

Boutet de Monvel and Khorunzhy initiated the study of mesoscopic linear eigenvalue statistics, i.e., the derivation of Gaussian fluctuations for the random variable 

    ∑_i=1^Nf(λ_i-E_0/η_0)-𝔼[∑_i=1^Nf(λ_i-E_0/η_0)],

with fixed E_0 ∈ (-2,2) on mesoscopic scales N^-1≪η_0 ≪ 1. In <cit.>, they obtained CLTs for the test function (x-)^-1 on all mesoscopic scales for the GOE, and  N^-1/8≪η_0 ≪ 1 for symmetric Wigner matrices, respectively. Lodhia and Simm <cit.> extended the CLT for complex Wigner matrices and general test functions on scales N^-1/3≪η_0≪1. He and Knowles <cit.> used moment estimates for Green functions to prove the CLTs for all arbitrary Wigner matrices on the optimal scales N^-1≪η_0≪1. More recently, Landon and Sosoe <cit.> obtained similar CLT by means of the characteristic function. 


Mesoscopic central limit theorems are important tools in the theory of homogenization of Dyson Brownian motion (DBM) introduced by Bourgade, Erdős, Yau and Yin <cit.> to prove fixed energy universality of local eigenvalue statistics of Wigner matrices.  Landon, Sosoe and Yau <cit.> subsequently derived a mesoscopic CLT to show fixed energy universality of the DBM. The dynamical approach using Dyson Brownian motion to prove the universality of the eigenvalue statistics on microscopic scale for all symmetry classes was initiated by Erdős, Schlein and Yau in <cit.>; we refer to the surveys <cit.> for further details.  Mesoscopic central limit theorems combined with DBM were used by Landon and Sosoe <cit.> and by Bourgade and Mody <cit.> to derive Gaussian fluctuations of single eigenvalues, and in <cit.> to show Gaussian fluctuations of the determinant of Wigner matrices.



Mesoscopic CLTs can also be studied at the spectral edges, where the mesoscopic scales are N^-2/3≪η_0 ≪ 1. For the GUE, Basor and Widom <cit.> used asymptotics of the Airy kernel to prove mesoscopic CLTs at the edges. Min and Chen <cit.> subsequently considered edge CLTs for the GOE. Recently, Adhikari and Huang <cit.> obtained the mesoscopic CLT at the edges down to the optimal scale η_0 ≫ N^-2/3 for the DBM. 





 §.§ Deformed Wigner matrices



In the present paper we are interested in deformed Wigner matrices. A deformed Wigner matrix is an N × N random matrix of the form

    X_N=H_N+A_N,

where H_N is a real symmetric or complex Hermitian Wigner matrix and A_N is a real deterministic diagonal matrix. It is also known as the Rosenzweig-Porter model in the physics literature. Suppose the empirical eigenvalue distribution of A_N has a deterministic limiting measure, denoted by μ_α. It was shown by Pastur <cit.> that the empirical eigenvalue distribution of X_N converges weakly in probability to the free additive convolution of μ_sc and μ_α, denoted by μ_fc=μ_sc⊞μ_α; see also <cit.>.

A CLT for the linear eigenvalue statistics with test functions in C^2_c() was obtained by Ji and Lee <cit.> under a one-cut assumption on μ_fc. They also computed the expectation and variance in terms of μ_α. Dallaporta and Fevrier <cit.> obtained the CLT for general μ_fc. Their results are summarized in Theorem <ref> below. 



In the present paper, we study the fluctuations of the linear eigenvalue statistics (<ref>) in the mesoscopic regime. We choose μ_α properly such that the free additive convolution μ_fc is supported on a single interval and vanishes as a square root at the end-points. This edge behavior of the limiting eigenvalue distribution is quite common in random matrix theory, and sometimes referred to as regular edge. Denoting κ_0=κ_0(E_0) the distance from E_0 to the closest edge of the free convolution measure, we derive a CLT at energy E_0 on scales η_0 with N^-1≪η_0√(η_0+κ_0)≤ 1; see Theorem <ref>. This range of η_0 covers the global scale as well as all mesoscopic scales up to the spectral edges. For energies E_0 in the bulk and at the edges respectively, we compute the variances and biases explicitly on the mesoscopic scales, where we recover the formulas for the Gaussian ensembles. This shows the expected universality of the linear eigenvalue fluctuations on mesoscopic scales. The universality of the eigenvalue statistics on the microscopic scale was derived for the deformed GUE in the bulk <cit.> and at the edge <cit.>. For deformed Wigner matrices, the local bulk universality was obtained in <cit.> for a special class of A_N using moment matching, under a one cut assumption in <cit.> and in <cit.> for the general case using the DBM methods. The edge universality was derived in <cit.> using a Green function comparison method, and in <cit.> using DBM. More recently, quantum unique ergodicity for deformed Wigner matrices was derived in <cit.>.

In the proof of the main results, we follow the idea of <cit.> to compute the characteristic function of (<ref>) in combination with the Helffer-Sjöstrand formula and cumulant expansions; see (<ref>) below. Cumulant expansions were used in e.g. <cit.> to study the linear eigenvalue statistics of random matrices. We also rely on local laws for Green functions <cit.> and analytic subordination for the free convolution measure, as used in <cit.>.  



On the global scale, the method used to derive the CLT for deformed Wigner matrices  <cit.>  is insensitive to the behavior of the free convolution measure μ_fc. An interesting aspect of the free additive convolution and deformed Wigner matrices is that the densities may show other edge behaviors than square roots. For such setups, one expects mesoscopic CLTs at the edges with different scalings, variances and biases. This is a main motivation for us to study linear eigenvalue statistics at spectral edges. The local eigenvalue statistics at such critical edges are only partly understood, see e.g.  <cit.> for some results. At cusp points in the interior of the bulk spectrum the universality of the local eigenvalue fluctuations was recently proved in <cit.>.





 §.§ Sample covariance matrix

Sample covariance matrices form another class of archetypal random matrix models, with applications in multivariate statistical analysis. We consider the separable sample covariance matrices of the form H=Σ^1/2 X X^* Σ^1/2, where X is a M × N matrix with independent random variables, and Σ^1/2 is the square root of the M × M diagonal and positive definite matrix Σ. The dimensionality M is chosen to be proportional to the sample size N. Assuming that the eigenvalue distribution of Σ has a deterministic limit μ_σ, it was proved by Marchenko and Pastur <cit.> that the spectral measure of H approaches a deterministic probability measure. In the null case Σ=I, the limiting measure is called the Marchenko-Pastur distribution, μ_MP. For the non-null case Σ≠ I, the limiting measure is given by the free multiplicative convolution of μ_MP and μ_σ, denoted by μ_MP⊠μ_σ, see <cit.>. A CLT for the fluctuations of the linear eigenvalue statistics was first studied by Jonsson <cit.> for Wishart matrices where X has Gaussian entries. CLTs for linear eigenvalue statistics with analytic test functions for general sample covariance matrices were then studied by Bai and Silverstein in <cit.>. The regularity condition on the test functions was weakened by <cit.> for the null case and <cit.> for the non-null case. In the second part of this paper, we extend the techniques to derive corresponding CLTs for the mesoscopic eigenvalue statistics of sample covariance matrices; see Theorem <ref>.



 



 §.§ Related models

Deformed Wigner matrices are closely related to Dyson Brownian motion, for which mesoscopic CLTs were obtained inside the bulk <cit.> and at the regular edges <cit.>. The mesoscopic linear statistics were also studied for random band matrices <cit.>, sparse Wigner matrices <cit.>, mesoscopic eigenvalue density correlations for Wigner matrices <cit.>, invariant β-ensembles <cit.> and orthogonal polynomial ensembles <cit.>. The global fluctuations of the deformed GOE/GUE can also be studied using the framework of second order freeness <cit.>.




 §.§ Structure of this paper


Section <ref> contains the precise definitions, assumptions and the main results. The proof the main theorem is carried out in Section <ref>-<ref>. In Section <ref> and <ref>, we compute the variances and the biases in the bulk and at the edges. In Section <ref>, we consider sample covariance matrices and obtain the corresponding results. Some auxiliary results are proved in the Appendices.







 §.§ Notation
  We denote the upper half-plane by ^+{z∈ :  z>0} and the positive real line by ^+{x∈ : x>0}. For any vector v ∈^N, we use v_2 to denote the Euclidean norm. For a matrix A ∈^N × N, we denote by A_op its operator norm induced by the Euclidean vector norm. We use c,k and C,K to denote strictly positive constants that are independent of N. Their values may change from line to line. We use standard big O and small o notations. For X, Y ∈, we write X ∼ Y  if there exist constants c, C>0 such that c |Y| ≤ |X| ≤ C |Y|. We write X ≪ Y if there exists a small τ>0 such that |X| ≤ N^-τ |Y| for large N. We will use the following definition on high-probability estimates from <cit.>. 

Let 𝒳≡𝒳^(N) and 𝒴≡𝒴^(N) be two sequences of nonnegative random variables. We say 𝒴 stochastically dominates 𝒳 if, for all (small) ϵ>0 and (large) D>0,

    (𝒳^(N)>N^ϵ𝒴^(N))≤ N^-D,

for sufficiently large N≥ N_0(ϵ,D), and we write 𝒳≺𝒴 or 𝒳=O_≺(𝒴).

We often use the notation ≺ also for deterministic quantities, then (<ref>) holds with probability one. Stochastic domination has the following properties.
(Proposition 6.5 in <cit.>)
	
		
  * X ≺ Y and Y ≺ Z imply X ≺ Z;
		
  * If X_1 ≺ Y_1 and X_2 ≺ Y_2, then X_1+X_2 ≺ Y_1+Y_2 and X_1X_2 ≺ Y_1Y_2;
		
  * If X ≺ Y, Y ≥ N^-c_1 and |X| ≤ N^c_2 almost surely with fixed constants c_1 and c_2, then we have X ≺ Y.
	
 







§ MODEL AND MAIN RESULTS





 §.§ Model and assumptions


Let H ≡ H_N be an N× N real or complex Wigner matrix satisfying the following assumption.

	For a real (β=1) symmetric Wigner matrix H we assume that:
	
		
  * {H_ij|i≤ j} are independent real-valued centered random variables with H_ij=H_ji.
		
  * For i ≠ j, [(√(N) H_ij)^2]=1; [(√(N) H_ii)^2]=m_2 for some constant m_2>0. In addition, [(√(N) H_ij)^4]=W_4 for some constant W_4>0.
		
  * All entries have uniform sub-exponential decay, that is, there exist C_0>0 and θ>1 such that
		
    ( |√(N) H_ij| ≥ x ) ≤ C_0 e^-x^1/θ,  ∀ i,j .

		In particular, we have 
		
    [|√(N) H_ij|^p] ≤ C (θ p)^θ p  (p ≥ 3) .

	
	For complex (β=2) Hermitian Wigner matrix we assume that:
	
		
  * {H_ij,H_ij|i≤ j} are independent centered real-valued random variables with H_ij=H_ji.
		
  * For i ≠ j, [H_ij^2]=0 and [(√(N) |H_ij|)^2]=1; [(√(N) |H_ii|)^2]=m_2 for some constant m_2>0. In addition, [(√(N) |H_ij|)^4]=W_4 for some constant W_4>0.
		
  * The sub-exponential tail assumption in (<ref>) holds.
	



Let {A_N}=Diag(a_i) be a sequence of real deterministic diagonal N × N matrices with A_op uniformly bounded in N. The empirical spectral measure of A_N is defined by μ_A1/N∑_i=1^Nδ_a_i. 


For a probability measure ν on  denote by m_ν its Stieltjes transform, i.e.

    m_ν(z)∫_ν(x)/x-z ,     z∈^+ .

Note that m_ν :^+→^+ is analytic and can be analytically continued to the real line outside the support of ν. Moreover, m_ν satisfies lim_η↗∞ηm_μ(η)=-1. Conversely, if m : ^+→^+ is an analytic function with lim_η↗∞η m(η)=-1, then m is the Stieltjes transform of a probability
measure ν, m(z) = m_ν(z), for all z∈^+; see  <cit.>.

The following assumption ensures the existence of the weak limiting measure of μ_A. 



	There exists a deterministic and compactly supported probability measure denoted as μ_α, such that μ_A converges weakly to μ_α. In addition, there exists α_0>0 such that for any fixed compact set D_0 ⊂^+∪ with D_0 ∩supp (μ_α) =∅, 
	
    max_z ∈ D_0| ∫_μ_A(x)/x-z -∫_μ_α(x)/x-z| =O(N^-α_0),

	for sufficiently large N.
	


Define the deformed Wigner matrix as

    X_N A_N+H_N.

The eigenvalues of X_N are denoted as λ_i ∈. The empirical spectral measure of X_N is defined by μ_N(x)=1/N∑_i=1^N δ_λ_i. For z∈^+, we introduce the Green function, G(z), and its normalized trace as

    G(z) (X_N-zI)^-1     m_N(z)1/N  G(z)=∫_μ_N(λ)/λ-z ,

i.e., m(z)≡ m_N(z) is the Stieltjes transform of μ_N.

The empirical spectral distribution μ_N converges as N tends to infinity to the free additive convolution of μ_α and the standard semicircle law, denoted by μ̃_fc:=μ_α⊞μ_sc. The free convolution measure can be described by analytic subordination <cit.>: Its Stieltjes transform, m̃_fc, is the unique solution to the Pastur equation
	
    m̃_fc(z)=∫_1/a-z-m̃_fc(z)μ_α(a) ,

subject to the constraint m̃_fc(z)>0, z∈^+. 

Since the convergence speed in (<ref>) can be very slow, we work with a finite N version of the free convolution measure. Let μ_fc:=μ_A⊞μ_sc. The Stieltjes transform of μ_fc, denoted by m_fc, is hence the unique solution to
 
    m_fc(z)=1/N∑_i=1^N 1/a_i-z-m_fc(z) ,

such that m_fc(z)>0, z∈^+. Note that μ_fc depends on N, but is deterministic.

 Biane <cit.> proved that μ_fc and μ_fc are absolutely continuous probability measures whose densities, are analytic wherever positive. We denote the density functions by ρ_fc and ρ_fc. In general the measures ρ_fc and ρ_fc are supported on several disjoint intervals and may have irregular edges where the densities do not vanish as a square root or have cusp points inside the support. The following assumption will rule out such scenarios.



	Let ℐ be the smallest interval that contains the support of μ_α, and assume that
	
    inf_x ∈ℐ∫_μ_α(a)/(a-x)^2≥ 1+ w,

	for some constant w>0 (the left side may be infinite). Similarly, let ℐ̂ be the smallest interval that contains the support of μ_A, and assume that
	
    inf_x ∈ℐ̂∫_μ_A(a)/(a-x)^2≥ 1+ w,

	for sufficiently large N.

The above assumption ensures that the density functions ρ_fc and ρ̃_fc  are supported on a single interval (for N sufficiently large) and vanish as square roots at the endpoints of the support. 



(Lemma 2.4, 3.2 and 3.5 in <cit.>)
	Under Assumption <ref>, there exists L̃_- and L̃_+ ∈, such that supp ρ̃_fc =[L̃_-,L̃_+], and ρ̃_fc is strictly positive in (L̃_-,L̃_+). Moreover, there exists C>1 such that
	
    C^-1√(κ̃)≤ρ̃_fc(E) ≤ C √(κ̃) ,         E∈[L̃_-,L̃_+] ,

	where κ̃min{ |E-L̃_-|, |E-L̃_+|}. The endpoints L̃_± are the two real solutions to the equation 
	
    ∫_μ_α(a)/(a-L_±-m_fc(L_±))^2 =1 .

	The  same holds true, for sufficiently large N, if we replace μ_α, ρ̃_fc, L̃_± and κ̃ by μ_A, ρ_fc, L_± and κ, respectively. Here [L_-,L_+] is the support of ρ_fc and κmin{ |E-L_-|, |E- L_+|  }. 
	



 §.§ Local law for the deformed Wigner matrices


We introduce the spectral domain, 

    D':={z=E+η:  |E| ≤ M, N^-1+c≤η≤  3 },

where M>1+max{ |L̃_-|, |L̃_+|} and c>0 is small. Define the deterministic control parameters


    Ψ(z):=√( m_fc(z)/N |η|) +1/N |η| ,    Θ(z):=1/N |η| ,      z=E+η∈∖.

Using (<ref>), (<ref>) in Lemma <ref> below, we have 

    C N^-1/2≤Ψ(z) ≪ 1 ,        z ∈ D' .

The following local law for the Green function was proved in <cit.>.
(Local law for the deformed Wigner matrix, Theorem 2.10 in <cit.>)
	Under the Assumptions 2.1-2.3, the following holds

    max_ij| G_ij(z)- δ_ij1/a_i-z-m_fc(z)| ≺Ψ(z),   | N^-1 G(z)-m_fc(z) | ≺Θ(z),

uniformly for z∈ D'.


 The local law gives  strong rigidity estimates for the eigenvalues of X_N. It also gives an upper bound, up to factors of N^ϵ, on the size of the fluctuations G(z) - G(z). It is hence natural to study the fluctuations of G(z) - G(z). The CLT for the linear eigenvalue statistics for general test functions is proved in <cit.> and <cit.> on global scale when z is order one. Via the Helffer-Sjöstrand functional calculus, a CLT for the resolvent can be translated to a CLT for the linear statistics.



	
Under the Assumptions <ref>-<ref>, for any φ∈ C_c() which is analytic on a neighborhood of [L̃_-, L̃_+], the random variable ∑_i=1^N φ(λ_i)-N ∫_φ(x) ρ_fc(x)  x converges in distribution to the Gaussian random variable with mean M(φ)=-1/2 π∫_Γφ(z) b̃(z)  z, and variance
 V(φ)=1/(2 π)^2∫_Γ∫_Γφ(z_1) φ(z_2) K̃(z_1,z_2)  z_1  z_2, where 
 
    b̃(z):=m̃”_fc(z)/2(1+m̃'_fc(z))^2( (m_2-1)+m̃'_fc(z)+(W_4-3)m̃'_fc(z)/1+m̃'_fc(z)),
 
 and K̃(z_1,z_2):=(m_2-2) ∂^2 Ĩ/∂ z_1 ∂ z_2+(W_4-3) ( Ĩ∂^2 Ĩ/∂ z_1 ∂ z_2 +∂Ĩ/∂ z_1∂Ĩ/∂ z_2)+2/(1-Ĩ)^2( ∂Ĩ/∂ z_1∂Ĩ/∂ z_2 +(1-Ĩ) ∂^2 Ĩ/∂ z_1 ∂ z_2), with 
 
    Ĩ(z_1,z_2)∫_1/(x-z_1-m̃_fc(z_1))(x-z_2-m̃_fc(z_2))μ_α(x).

	Here Γ is a rectangular contour with vertices (a_±± v_0) so that ±(a_±-L̃_±)>0 and Γ lies within the analytic domain of φ.

Using ideas of M. Shcherbina <cit.>, the above result can be extended to C^2_c() test functions.  In <cit.>, the corresponding result was obtained for the multi-cut regime. 



 §.§ Main results


Choose E_0 ∈ [-1+L̃_-, 1+L̃_+] and N^-1≪η_0 ≪ 1. Consider a test function g ∈ C_c^2() and set

    f_N(x):=g( x-E_0/η_0).

We will write f_N as f for notational simplicity. Define

    κ_0:=dist((f),{L_+,L_-}).


Following <cit.>, we study the characteristic function 

    ϕ(λ):=[e(λ)],   e(λ):=exp{ıλ( f(X_N)- f(X_N)) },     λ∈.

Let τ>0 be an arbitrary small constant and define

    Ω_0 := { x + ı y ∈ :  |y| ≥ N^-τη_0  } .

A key observation in <cit.> is that working on Ω_0 instead of all , effectively removes the ultra-local scales without affecting the mesoscopic linear statistics.




	Let X_N be a deformed Wigner matrix satisfying Assumptions <ref>-<ref>. Let η_0 √(κ_0+η_0)≥ N^-1+c_0 for some c_0>0. Then there exists a small 0 < τ< c_0/16, such that the characteristic function (<ref>) satisfies
	
    ϕ'(λ)=-λϕ(λ) V(f)+ℰ̃,      V(f):=1/π^2∫_Ω_0∫_Ω_0(z_1) (z_2) K(z_1,z_2) ^2z_1 ^2z_2,

	where  is an almost analytic extension of f given in (<ref>) below and β=1,2 is the symmetry parameter. 
	The kernel K is given by
	
    K(z_1,z_2):= ( (m_2-2/β) I +(W_4-1-2/β)/2I^2)+2/β∂/∂ z_1(1/1- I∂  I/∂ z_2),

	with
	
    I(z_1,z_2):=∫_1/(x-z_1-m_fc(z_1))(x-z_2-m_fc(z_2))μ_A(x),

	and the error ℰ̃ is bounded by
	
    |ℰ̃|=O_≺( |λ| log N N^- τ)+O_≺( (1+|λ|^4)N^3 τ/N η_0 √(κ_0+ η_0))+O_≺( (1+|λ|^4)N^2 τ/√(N η_0 √(κ_0+η_0))),

	provided that V(f)=O(1). 
	



Proposition <ref> implies the following result.

	Under the same assumptions as in Proposition <ref>, if we further assume that there exist c,C>0 such that c ≤ V(f) ≤ C for sufficiently large N, then f(X_N)- f(X_N)/√(V(f)) converges in distribution to a standard Gaussian random variable.


We remark that the Theorem <ref> applies to the global scale as well as the mesoscopic scales. The expectation of f(X_N) has the following asymptotic expansion, which matches the result in <cit.> on the global scale.
 Under the same assumptions as in Proposition <ref>, the so-called bias is given by
	
    f(X_N)-N ∫_ f(x) ρ_fc(x)  x=1/2 π∫_Ω_0(z) b(z)  ^2z +O(N^-τ)+O_≺( N^2 τ/√(N η_0 √(κ_0+η_0))),

	where  is given in (<ref>) below, and
	
    b(z):=(2/β-1)1/1-I_s(z) I_s(z)/ z +(m_2-2/β)  I_s(z)/ z+(W_4-1-2/β) I_s(z) I_s(z)/ z,

	with
	
    I_s(z):=∫_1/(x-z-m_fc(z))^2μ_A(x).



 Note that the variance V(f) in (<ref>) and the bias in (<ref>) are N-dependent and their formulas depend explicitly on the free convolution measures. We will compute their limits on the mesoscopic scales as N goes to infinity in order to obtain the following universal CLTs in the bulk and at the regular edges respectively.
(Mesoscopic CLT in the bulk)
	Let X_N be a deformed Wigner matrix satisfying Assumptions <ref>-<ref>. Let N^-1+c≤η_0 ≤ N^-c with some small c>0, fix E_0 ∈ (L̃_-,L̃_+) such that κ_0 > c_0, for some c_0>0 and large N. Then, for any test function g ∈ C^2_c(), the linear statistics 
	
    ∑_i=1^N g ( λ_i-E_0/η_0)-N ∫_ g (x-E_0/η_0) ρ_fc(x)  x

	converges in distribution to a Gaussian random variable with mean zero and variance
	
    1/2 βπ^2∫_∫_(g(x_1)-g(x_2))^2/(x_1-x_2)^2 x_1  x_2=1/βπ∫_ |ξ| |ĝ(ξ)|^2 ξ,

	where ĝ(ξ):=(2 π)^-1/2∫_ g(x) e^-ıξ x x. In particular, the bias vanishes in the bulk regime.



(Mesoscopic CLT at the edge)
	Let X_N be a deformed Wigner matrix satisfying Assumptions <ref>-<ref>. Let N^-2/3+c≤η_0 ≤ N^-c with some small c>0.  For any function g ∈ C^2_c(),
	the linear statistics (<ref>) with E_0=L_+
	converges in distribution to a Gaussian random variable with mean (2/β-1) g(0)/4 and variance
	
    1/4 βπ^2∫_∫_( g(-x^2)-g(-y^2)/x-y)^2  x  y=1/2 βπ∫_ |ξ| |ĥ(ξ)|^2 ξ,

	where h(x)=g(-x^2) and ĥ(ξ):=(2 π)^-1/2∫_ h(x) e^-ıξ x x. At the left edge E_0=L_-, we obtain a similar CLT with h(x)=g(x^2).




Remark: The bulk variance (<ref>) agrees with the GOE/GUE. For the edges, the bias and variance in (<ref>) coincide with those of the GUE/GOE obtained in <cit.> and the Dyson Brownian motion in <cit.>.


Remark: We remark that our assumption that the fourth moments of the off-diagonal entries are identical can easily  be relaxed in the above theorems. The regularity condition we impose on the test function g is clearly not optimal, and we expect results can be extended to C^1,r,s() functions; see <cit.>. The CLTs also hold true if we consider the resolvent test function g(x)=1/x-.

Finally, for test functions in C^2_c(), we can relax the single support condition for μ_fc by assuming instead that the cuts of the support of μ_fc are separated by order one and the density ρ_fc has square-root decay near the edges. 





§ PROOF OF PROPOSITION <REF>

In this section, we prove Proposition <ref> by reducing it to the main technical result Lemma <ref>. Recall the scaled test function f on scale η_0 from (<ref>). There are constants such that

    f_1 ≤ C η_0;     f'_1 ≤ C';     f”_1 ≤C”/η_0.


We use the Helffer-Sjöstrand formula to link f(X_N) to the Green function of X_N.
(Helffer-Sjöstrand formula) 
	Let f ∈ C_c^2() and χ(y) be a smooth cutoff function with support in [-2,2], with χ(y)=1 for |y| ≤ 1. Define its almost-analytic extension 
	
    f̃(x+ y):=(f(x)+ y f'(x)) χ(y).

	Then we have
	
    f(λ)=1/π∫_f̃(z)/λ-z^2z=1/2 π∫_^2 y f”(x) χ(y)+( f(x)+ y f'(x) ) χ'(y)/λ-x- y x  y,

	where z=x+ y, =1/2(+), and ^2z is the Lebesgue measure on .


Therefore, we write

    f(X_N) - f(X_N) =1/π∫_(z) ((G(z))- G(z))^2z.

Plugging the above equation in e(λ) given by (<ref>), we have

    e(λ)=exp{ıλ/π∫_(z) ((G(z))- G(z)) ^2 z}.

Taking the derivative of the characteristic function given in (<ref>), and applying (<ref>), we get

    ϕ'(λ)=ı/π∫_(z)   [ e(λ)((G(z))- G(z)) ]  ^2z.


Following <cit.>, we restrict the domain of the spectral parameter to Ω_0, as the very local scales do not contribute to ϕ(λ).  We write

    f(X_N) - f(X_N)=1/π( ∫_Ω_0 +∫_Ω^c_0) ∂f̃(z)/∂z̅((G(z))- G(z)) ^2z.

Recall f̃ in (<ref>) and the definition of Ω_0 in (<ref>). Since χ(y)=1 for |y| ≤ 1, we can write the second integral in (<ref>) with z=x+ y as

    N/π∫_∫_0^η_0/N^τ y f”(x) (m_N(z) - m_N(z)  )  x  y   =-N/π∫_∫_0^η_0/N^τ   y f”(x)  (m_N(z)- m_N(z) )  x  y,

where we used the fact that m_N(z)=m_N(z). We now choose a small τ>0 such that N^-1≪ y_0:=√(η_0/N^1+τ)≤ N^-τη_0. In the regime y ∈ [y_0, N^-τη_0], the integral can be estimated using the local law (<ref>), (<ref>) and Lemma <ref>, i.e.,

    |N/π∫_∫_y_0^N^-τη_0   y f”(x) (m_N(z)- m_N(z) )  x  y| ≺| ∫_ f”(x) x ∫_y_0^N^-τη_0 y|=O_≺(N^-τ).

In the regime y ∈ [0, y_0], the local law is not sharp but instead we use the fact that y → m_N(x+ y)y is increasing. That is,

    |N/π∫_∫_0^y_0 y f”(x)  (m_N(z)- m_N(z) )  x  y| =O_≺(N y_0^2/η_0)=O_≺(N^-τ).

Therefore, we have from (<ref>) that

    f(X_N) - f(X_N)=1/π∫_Ω_0f̃(z)  ((G(z))- G(z)) ^2z+O_≺(N^-τ).



Using the same argument, since |e(λ)|=1, we have

    ϕ'(λ)=ı/π∫_Ω_0(z)   [ e(λ)((G(z))- G(z)) ]  ^2z+O_≺(N^-τ).

Similarly, we restrict the integration domain of e(λ) in (<ref>) to Ω_0. Let

    :=exp{ıλ/π∫_Ω_0(z) ((G(z))- G(z)) ^2 z  }.

In addition, (<ref>) implies that |e(λ)-| = O_≺( |λ| N^ -τ). We also have ||=1, using |e(λ)|=1 and G(z)= G(z). If we further replace e(λ) by  in (<ref>), then we get

    ϕ'(λ)=ı/π∫_Ω_0(z)  [ ( ((G(z))- G(z)) ] ^2 z+O_≺( |λ| log N N^ -τ).


The last error term on the right side, and many error terms below, are estimated using the following lemma, which is a variant of Lemma 4.4 in <cit.>. The proof is provided in Appendix B.

		Suppose h(z) is a holomorphic function on Ω_0 and |h(z)| ≤K/| z|^s for some constants s, K ≥ 0, then there exists some constant C such that
	
    |∫_Ω_0(z) h(z) ^2z | ≤ CK N^τ sη_0^1-s.

	For 1 ≤ s ≤ 2, the bound is sharpened to  CK log (N) η_0^1-s.


Thus, in order to study ϕ'(λ), it is sufficient to estimate [  ((G(z))- G(z)) ]. The key input is the following cumulant expansion formula.
(Cumulant expansion formula)
	Let h be a real-valued random variable with finite moments, and f is a complex-valued smooth function on  with bounded derivatives. Let c_k be the k-th cumulant of h, given by c_k(h):=(-)^k/ tlog e^ t h|_t=0. Then for any fixed l ∈, we have
	
    [h f(h)]=∑_k=0^l 1/k! c_k+1(h)[ f^(k)(h) ]+R_l+1,

	where the error term satisfies
	
    |R_l+1| ≤ C_l  |h|^l+2sup_|x| ≤ M |f^(l+1)(x)| +C_l [ |h|^l+2 1_|h|>M] f^(l+1)_∞,

	and M>0 is an arbitrary fixed cutoff.

For reference, we refer e.g. to Lemma 3.1 in <cit.>. We give the  proof the following lemma in Section 5. 


	For any z:=E+η∈Ω_0∩ D', see (<ref>), and κ:=min{ |E-L_-|, |E- L_+|  }. we have
	
    [ ( G(z)-  G(z))] =λ/π [] ∫_Ω_0(z') K(z,z')  ^2 z'+ℰ(z),

	where K is given in (<ref>)
	and ℰ(z) is analytic in Ω_0 and satisfies
	
    ℰ(z)=O_≺( 1+|λ|^4/√(κ+η)) ( (κ+η)^1/4/√(N η^3)+1/√(N η^2)+1/√(N η_0 η)+ 1/N η_0 η +1/Nη^2).



Admitting Lemma <ref> and plugging in (<ref>),  we have

    ϕ'(λ)=- λ [] V(f)+O_≺( |λ| log N N^- τ)+ℰ̃,

where 

    V(f)=1/π^2∫_Ω_0∫_Ω_0(z)  (z') K(z,z') ^2 z ^2 z',     ℰ̃=ı/π∫_Ω_0(z)  ℰ(z) ^2 z.

By the definition of κ_0 in (<ref>), κ≥κ_0. Moreover |η| ≥ N^-τη_0, for z∈Ω_0. Using Lemma <ref>, we hence obtain the estimate 

    ℰ̃=O_≺( (1+|λ|^4)N^3 τ/N η_0 √(κ_0+ η_0))+O_≺( (1+|λ|^4)N^2 τ/√(N η_0 √(κ_0+η_0))).
 
Assuming V(f) ≺ O(1), we can replace  by e(λ) with error O_≺( |λ| N^ -τ). Thus we have completed the proof of Proposition <ref>.







§ PROPERTIES OF THE FREE CONVOLUTION




 §.§ Properties of m_fc and m̃_fc


In this subsection, we recall some properties of the Stieltjes transforms m_fc and m̃_fc of the free convolution measures.
Let κ=κ(E) be the distance from E to the closest spectral edge, i.e.,

    κ:=min{ |E-L_-|, |E-L_+|  }.

Similarly define κ̃:=min{ |E-L̃_-|, |E-L̃_+|  }. 
Define the spectral domain

    D:={z=E+ıη:  |E|<M, 0< η≤ 3 }.

(Lemma 3.5, Lemma A.1 in <cit.>)
	
		
  * For all z ∈ D, there exists C>1 such that 
		
    C^-1√(κ̃ +η)≤|m̃_fc(z)| ≤ C √(κ̃ +η),

		if E ∈ [L̃_-,L̃_+]. If E ∈ [L̃_-,L̃_+]^c, then 
		
    C^-1η/√(κ̃+η)≤|m̃_fc(z)| ≤ C η/√(κ̃+η).

                
		
  * (Stability bound) There exists C>1, such that 
		
    C^-1≤  |a-z-m̃_fc(z)| ≤ C,

		uniformly for z ∈ D and a ∈supp(μ_α) .
		
  * For all z ∈ D, there exist k,K>0 such that
		
    k √(κ̃ +η)≤| 1-∫_1/(x-z-m̃_fc(z))^2μ_α(x) | ≤ K √(κ̃ +η).

		
		
  * There exist C>0 and c_0>0 such that for all z ∈ D satisfying κ̃+η≤ c_0,
		
    C^-1≤| ∫_1/(x-z-m̃_fc(z))^3μ_α(x) | ≤ C;

		moreover, there exists C>1 such that for all z ∈ D,
		
    | ∫_1/(x-z-m̃_fc(z))^3μ_α(x) | ≤ C.

	



The following lemma implies that m_fc behaves similarly as m̃_fc, for sufficiently large N.
 (Lemma 3.6 in <cit.>)
	Under Assumptions 2.2 and 2.3, for sufficiently large N, statements 1-4 in Lemma <ref> hold true with m̃_fc, κ̃, μ_α and L̃_± replaced by m_fc, κ, μ_A and L_± respectively. Moreover, the constants in these inequalities can be chosen uniformly in N for sufficiently large N. Furthermore, there exists c>0 such that
	
    max_z ∈ D| m̃_fc(z)- m_fc(z) | ≤  N^-cα_0/2,     |L̃_±-L_±|≤  N^-cα_0,

	for sufficiently large N.





Recall the function I(z_1,z_2) given in (<ref>) and I_s(z) in (<ref>). By direct computation, one proves the following lemma.

	For z_1 ≠ z_2, we have
	
    I(z_1,z_2)=m_fc(z_1)-m_fc(z_2)/z_1+m_fc(z_1)-z_2-m_fc(z_2);     I_s(z) =  m'_fc(z)/1+m'_fc(z).



As a result of Lemmas <ref>, <ref> and <ref>, we have the following lemma. 



	There exists C>1 such that
	
    |I(z_1,z_2)| ≤ C;      |I_s(z)| ≤ 1;      C^-1√(κ+η)≤ |1-I_s(z)| ≤ C √(κ+η);

	
    |m_fc(z)| ≤ C;      |m_fc'(z)| ≤C/√(κ+η);     |m_fc”(z)| ≤C/√(( κ+η)^3),

	uniformly for z, z_1, z_2 ∈ D.


The proof of the above two lemmas can be found in Appendix B.



 §.§ Properties of the Green function

As a more general version of the local law in Theorem <ref>, we introduce the anisotropic local law. Recall the control parameters Ψ and Θ from (<ref>).
(Theorem 12.2, 12.4 in <cit.>; Theorem 2.1, 2.2 in <cit.>; Theorem 2.6 in <cit.>)
	For any deterministic vector v,w ∈^N and matrix B ∈^N × N, we have
	
    | ⟨ v, G(z) w ⟩  - ⟨ v, G(z) w ⟩| ≺v_2 w_2 Ψ(z),    | N^-1 (B (G(z)-G(z)))  | ≺B_opΘ(z),

	uniformly in z ∈{z=E+η:  |E| ≤ρ^-1, N^-1+ρ≤η≤ρ^-1}, where ρ is small so that ρ^-1≥A_op, and G=Diag(1/a_i-z-m_fc(z)).
	









§ PROOF OF LEMMA <REF>


For the simplicity of the presentation, we consider only the real symmetric case here. The complex case being similar is proved in Appendix A.  For notational simplicity, let 

    g_i(z):=1/a_i-z-m_fc(z),      z ∈∖supp(μ_fc).


Before we proceed the proof of Lemma <ref>, we state a useful lemma. 



For any i,j, we have	

    ∂/∂ H_ij   =- (2-δ_ij)λ/π∫_Ω_0(z) / zG_ji^2 z;
    ∂^2 /∂^2 H_ij   = (2-δ_ij) λ/π∫_Ω_0(z) / z (g_i(z) g_j(z)) ^2 z+O_≺( (1+|λ|)^2/√(N η_0)).

	In general, for any integer k ∈, we have
	
    | ∂^k G_ij/∂ H^k_ij| ≺ O(1);   | ∂^k  /∂^k H_ij| ≺ O((1+|λ|)^k).


The above lemma follows from the relation
	
    ∂ G_ij/∂ H_ab=-G_ia G_bj+G_ibG_aj/1+δ_ab.

The details are provided in Appendix <ref>. Now we are ready to prove Lemma <ref>. 





	By the definition of the resolvent function, we have
	
    (z-a_i) G_ii=(HG)_ii-1.

	Thus we obtain that
	
    (z-a_i)  [(G_ii- G_ii)]=∑_j=1^N (  [ H_ijG_ji ]-  [H_ijG_ji]  [] ).

Using the cumulant expansion Theorem <ref>, we obtain

    (z-a_i)  [(G_ii- G_ii)]=I_1+I_2+I_3+O_≺(N^-3/2(1+|λ|^4)),
	
where

    I_1:   =1/N∑_j=1^N c^(2)_ij( [∂/∂ H_ij G_ji]+[ ( ∂ G_ji/∂ H_ij - [∂ G_ji/∂ H_ij] ) ] );
    
    I_2:   =1/2! N^3/2∑_j=1^N c^(3)_ij( [∂^2 /∂^2 H_ij G_ji]+2[∂/∂ H_ij∂ G_ji/∂ H_ij]+[(1-)(∂^2 G_ji/∂^2 H_ij) ] );
    
    I_3:   =1/3! N^2∑_j=1^N c^(4)_ij( [∂^3  /∂^3 H_ij G_ji]+3 [∂^2  /∂^2 H_ij∂ G_ji/∂ H_ij]+3[∂/∂ H_ij∂^2 G_ji/∂^2 H_ij]
           +[ ( ∂^3 G_ji/∂^3 H_ij-[∂^3 G_ji/∂^3 H_ij]) ]  ).



Here c^(k)_ij denotes the k-th cumulant of √(N) H_ij. In particular, 
	
    c^(1)_ij=0;      c^(2)_ij=1+(m_2-1)δ_ij;      c^(4)_ij=W_4-3   (i ≠ j).

The last term on the right side of (<ref>) is estimated by (<ref>), (<ref>) and Lemma <ref>. Note that for z ∈Ω_0 ∩ D', we have the deterministic bound |G_ij| ≤G_op≤ (z)^-1=O(N^c). Combining with ||=1, we can use the last statement of Lemma <ref>. We will use this argument throughout the proof. The error terms in this section are all uniform in z ∈Ω_0∩ D'. In the following, we estimate I_1, I_2 and I_3 respectively.
	
	
	
	

 §.§ Estimate on I_1

	Using (<ref>), we have for each i, 
	
    I_1=   -1/N [ ((G^2)_ii - (G^2)_ii)] -1/N [ ( G G_ii - G G_ii) ]
       -m_2-2/N [ (G_ii G_ii- G_ii G_ii )]
    	+1/N∑_j=1^N (1+(m_2-1) δ_ij) [ ∂/∂ H_ij G_ji]
    =:   A_1(i)+A_2(i)+A_3(i)+A_4(i).

	The first term can be written as
	
    A_1(i)=-1/N[   (1-) / z   (G(z))_ii] ≺Ψ(z)/Nz,

	with Ψ(z) as in (<ref>). The last step follows from the local law and the Cauchy integral formula.
	Similarly using the local law, the second term A_2(i) can be written as
	
    A_2(i)   =-1/N [ ( G (G_ii -  G_ii) +  G_ii( G- G) + G  G_ii- G G_ii)  ] 
       =- m_fc(z)  [ (G_ii- G_ii)]-1/Ng_i(z)  [ ( G-  G)]+O_≺( Θ(z) Ψ(z)),

	with Θ(z) as in (<ref>) and g_i(z) in (<ref>). Here the first term of A_2 will be moved to the left side of the equation (<ref>).
	In addition, the local law also implies that A_3(i)=O_≺(Ψ(z)/N).
	
	Note that A_4 is a leading term of I_1. Using the local law, (<ref>) and Lemma <ref>, we write
	
    A_4(i)=   A_41(i)+A_42(i)+O_≺((1+|λ|)N^-1Ψ(z) ),

	where
	
    A_41(i)=1/N∑_j=1^N  [ ∂/∂ H_ij (1+ δ_ij)G_ji],     A_42(i)=m_2-2/N[ ∂/∂ H_ii g_i(z)].

	
	We compute these two terms below in the Section 5.4.
	
	
	
	

 §.§ Estimate on I_2

	In this subsection, we will show that I_2 is negligible, which can be written as
	
    1/2 N^3/2∑_j=1^N c^(3)_ij(  [∂^2  /∂^2 H_ij G_ji] +2[∂/∂ H_ij∂  G_ji/∂ H_ij]  +[ ( ∂^2  G_ji/∂^2 H_ij -∂^2  G_ji/∂^2 H_ij) ] )
       =:B_1(i)+B_2(i)+B_3(i).

	
	
	First, we study the last term B_3(i). Using (<ref>) and the local law, we have
	
    B_3(i)   =-1/2 N^3/2∑_j=1^N c_ij^(3)[ (  6 G_iiG_jj G_ij +2 (G_ij)^3-6  [G_iiG_jj G_ij] -2 (G_ij)^3) ]
       =-3/N^3/2∑_j=1^N [  c_ij^(3) g_i(z) g_j(z) ( G_ij - G_ij) ] +O_≺( N^-1/2Ψ^2(z) ).

	Next, we estimate 1/√(N)∑_j=1^N c^(3)_ij g_j(z) G_ij, using the anisotropic local law Theorem <ref>. Let v_j=δ_ij and w_j=1/√(N)c^(3)_ij g_j(z). And w_2 is bounded because of the stability bound (<ref>) and the moment condition (<ref>). 	
	
	Note that Theorem <ref> holds for vector entries w_j and v_j that are deterministic constants. As in our setting w_j depend on z, we use a continuity argument to show that 
		
    | 1/√(N)∑_j=1^N c^(3)_ij g_j(z) G_ij(z) -1/√(N) c^(3)_ii (g_i(z))^2 | ≺Ψ(z) ,

uniformly in z∈ D'. Indeed, choose a lattice Δ of the domain D' in (<ref>), with |Δ|=N^100. Then for any z ∈ D', there exists some point p ∈Δ, such that |z-p| ≤ N^-10. The anisotropic local law (<ref>) combined with a union bound implies 

    ℙ( ∃ p ∈Δ: | 1/√(N)∑_j=1^N c^(3)_ij g_j(p) G_ij(p) -1/√(N) c^(3)_ii (g_i(p))^2 | ≥ N^ϵΨ(p)) ≤ N^-D+100.

Recall that g_j(z)=1/a_j-z-m_fc(z). Using (<ref>), Lemma <ref> and the fact that |G_ij(z)| ≤1/η, the function 1/√(N)∑_j=1^N c^(3)_ij g_j(z) G_ij(z)-1/√(N) c^(3)_ii (g_i(z))^2 as well as Ψ(z) are Lipschitz continuous on D' with Lipschitz constant at most N^3. Thus we obtain from (<ref>) that

    ℙ( ∃ z∈ D': | 1/√(N)∑_j=1^N c^(3)_ij g_j(z) G_ij(z) -1/√(N) c^(3)_ii (g_i(z))^2 | ≥ 2N^ϵΨ(z)) ≤ N^-D+100 ,
		
		which implies (<ref>).

		
		
	Next, using (<ref>) and Ψ(z) ≥ C N^-1/2, we have | 1/√(N)∑_j=1^N c^(3)_ij g_j(z) G_ij(z) | ≺Ψ(z). Therefore, we obtain the upper bound
	
    B_3(i)=O_≺( N^-1Ψ(z) )+O_≺(N^-1/2Ψ^2(z) )=O_≺(N^-1/2Ψ^2(z) ).

	
	For the second term, by (<ref>), (<ref>), and the local law we have
	
    B_2(i)=   -1/N^3/2∑_j=1^N c^(3)_ij[∂/∂ H_ij (G_jiG_ji +G_iiG_jj) ]
    
    	=   -1/N^3/2∑_j=1^N c^(3)_ij[∂/∂ H_ij g_i(z) g_j(z)]   +O_≺( (1+|λ|) Ψ(z)/N √(η_0))
    
    	=   2 λ/π N^3/2[ ∑_j=1^N c^(3)_ij∫_Ω_0(z') / z'(G(z'))_ji^2 z'  g_i(z) g_j(z)  ]+O_≺( (1+|λ|) Ψ(z)/N √(η_0)).

	By the same argument as in (<ref>) and the Cauchy integral formula, we have
	
    | / z'1/√(N)( ∑_j=1^N c^(3)_ij g_j(z) (G(z'))_ji) |=O_≺( Ψ(z')/|z'|).

	Using the stability bound (<ref>) and Lemma <ref>, we have
	
    |B_2(i)| ≺1+|λ|/N √(N η_0)+(1+|λ|) Ψ(z)/N √(η_0)=O_≺((1+|λ|) Ψ(z)/N √(η_0)).

	
	Similarly, by plugging (<ref>) in the expression of B_1, we have
	
    B_1(i)=   λ/π N^3/2∑_j=1^N c^(3)_ij[  ( ∫_Ω_0(z') / z'  ( g_i(z') g_j(z'))  ^2 z') G_ji]
    	+O_≺(  (1+|λ|^2) Ψ(z)/N √(η_0)).

	Using the anisotropic local law, we have	
	
    B_1(i) =O_≺(  (1+|λ|^2)N^-1Ψ(z))+O_≺(  (1+|λ|^2) Ψ(z)/N √(η_0))=O_≺(  (1+|λ|^2) Ψ(z)/N √(η_0)).

	
	
	
	
	

 §.§ Estimate on I_3

	It is not hard to show that the diagonal terms for i=j are negligible. Thus we can replace the fourth cumulants by W_4-3.  There are four terms in I_3 and we denote them as D_1(i), D_2(i), D_3(i) and D_4(i) respectively.
	
	First, we look at D_1. By the local law and (<ref>), we have |D_1(i)| ≺ (1+|λ|^3) N^-1Ψ(z). Similarly, using (<ref>), (<ref>) and the local law, we have |D_3(i)| ≺(1+|λ|) Ψ(z) / N √(N η_0).
	For the last term D_4, using (<ref>) and the local law, we obtain that
	
    D_4(i)=1/6 N^2∑_j=1^N[  (1-) ( 36 G_iiG_jj(G_ij)^2+6 (G_ii)^2 (G_jj)^2+6 (G_ij)^4 ) ] ≺ N^-1Ψ(z).

	Finally, we look at the leading term D_2(i). Using the local law and (<ref>), we have
	
    D_2(i)=   -W_4-3/2 N^2∑_j=1^N [∂^2  /∂^2 H_ij( (G_ji)^2+G_ii G_jj) ]
    
    	=   -W_4-3/2 N^2∑_j=1^N [∂^2  /∂^2 H_ij g_i(z) g_j(z)  ]+O_≺( (1+|λ|^2)N^-1Ψ(z) ).

	
	
	
	

 §.§ Adding up the contributions to (<ref>)

	Summing up the contributions from the previous subsections, we write (<ref>) as
	
    (z-a_i+m_fc)  [(G_ii- G_ii)]=-1/N g_i(z) [ ( G-  G)]+A_41(i)+A_42(i)+D_2(i)+ϵ(i),

	where D_2 is given in (<ref>) and A_41, A_42 in (<ref>), and ϵ(i) is the error term obtained in the previous subsections. Thanks to the stability bound (<ref>), we can divide  both sides by z-a_i+m_fc to get
	
    [(G_ii- G_ii)]=   1/N (g_i(z))^2  [ ( G-  G)]-g_i(z)( A_41(i)+A_42(i)+D_2(i) +ϵ(i) ).

	Summing over i and rearranging, we find 
	
    (1-I_s(z))  [ ( G-  G)] =-∑_i=1^N g_i(z) ( A_41(i)+A_42(i)+D_2(i) )+ℰ_1,

	where ℰ_1 is the linear statistics of ϵ(i). By the argument in Section 5.1-5.3, we get 
	
    ℰ_1=O_≺( (1+|λ|^4) N Ψ(z) Θ(z) )+O_≺( (1+|λ|^4) √(N)Ψ^2(z) )+O_≺( (1+|λ|^4) Ψ(z)/√(η_0)).
 
	
	Next, we study the leading terms of the right side of (<ref>). Plugging (<ref>) in (<ref>), we have
	
    ∑_i=1^N A_41(i)/z-a_i+m_fc= 2 λ/π N∑_i=1^N g_i(z) [ ∫_Ω_0(z')  (G(z')G(z))_ii^2 z' ].

	By the resolvent identity,
	
    G(z) G(z')=G(z)-G(z')/z-z',    z ≠ z',

	 we can write
        
    F(z,z'):=1/N∑_i=1^N g_i(z) (G(z')G(z))_ii=1/N∑_i=1^N 1/z' -zg_i(z) (G_ii(z')-G_ii(z)).

 We separate into two cases:

	Case 1: If z and z' belong to different half-planes, then we have 1/|z-z'|≤1/| z |.
	Thus by the anisotropic local law, we have
	
    | F(z,z')- 1/z' -z1/N∑_i=1^N g_i(z) (g_i(z')-g_i(z)) |    =1/|z'-z||1/N∑_i=1^N g_i(z) (G_ii(z')-g_i(z')) |  
    
    	+1/|z'-z|   |1/N∑_i=1^N g_i(z) (G_ii(z)-g_i(z)) | =O_≺( Θ(z)+Θ(z')/|z|).

	
	Case 2: If z and z' are in the same half-plane, without loss of generality, we can assume they both belong to the upper half plane. If |z-z'| ≥1/2z, then we can use the same argument as in Case 1. Thus it is sufficient to study when |z-z'| ≤1/2z, which means 2/3z' ≤z ≤ 2 z'. Note that
	
    | F(z,z')-1/z' -z1/N∑_i=1^N  g_i(z) (g_i(z')-g_i(z)) |  ≤|1/N∑_i=1^N (g_i(z)-g_i(z')) (G_ii(z')-g_i(z'))/z-z'| 
    
    	+| 1/N∑_i=1^N g_i(z') (G_ii(z')-g_i(z'))-1/N∑_i=1^N g_i(z) (G_ii(z)-g_i(z)) /z-z'|.

	Next, we look at the first term on the right side. By direct computation, we get
	
    | g_i(z)-g_i(z')/z-z'| ≤ |g_i(z) || g_i(z') | ( 1+| m_fc(z)-m_fc(z')/z-z'| ).

	When z,z' are in the same half plane, m_fc is analytic in the neighborhood of the segment connecting z and z', denoted as L(z,z'). Thus
	
    | m_fc(z)- m_fc(z')/z-z'| ≤sup_ω∈ L(z,z')|m'_fc(ω)| ≤C/z.

	Combining with (<ref>), we have
	
    | g_i(z)-g_i(z')/z-z'| ≤C'/z.

	Using the second statement of the anisotropic local law by letting B=(g_i(z)-g_i(z')/z-z') and the continuity argument as in (<ref>), we obtain that the first term is bounded as O_≺( Θ(z')/z). 
	
	For the second term, we write it as h(z)-h(z')/z-z', where 
    h(z):=1/N∑_i=1^N g_i(z)(G_ii(z)-g_i(z)).
 Since h is analytic in the neighborhood of L(z,z'), we have
	
    | h(z)- h(z')/z-z'| ≤sup_ω∈ L(z,z')|/ω h(ω)|.

	The anisotropic local law implies that sup_w ∈ L(z,z')|h(w)| ≺Θ(z). Using the Cauchy integral formula, the second term is O_≺( Θ(z)/z). Then we obtain the same upper bound as in Case 1.
	
	Therefore, in both cases, we have
	
    F(z,z')=1/z'-z(1/N∑_i=1^N g_i(z) g_i(z')-1/N∑_i=1^N g_i^2(z) ) +O_≺( Θ(z)/z)+O_≺( Θ(z')/z).

	
	Taking the derivative and using the Cauchy integral formula, we have
	
    F(z,z')
    	=( 1/1-I(z,z')∂ I(z,z')/∂ z) (1-I_s(z))+O_≺( Θ(z)/z|z'|)+O_≺( Θ(z')/z|z'|).

	Then by using Lemma <ref>, we have
	
    ∑_i=1^N A_41(i)/z-a_i+m_fc=    2 λ/π [] ∫_Ω_0(z') ( 1/1-I(z,z')∂ I(z,z')/∂ z (1-I_s(z)) ) ^2 z' 
       +O_≺( Θ(z)/z)+O_≺( 1/N η_0z).

	
	Similarly, plugging (<ref>) in (<ref>), we have
	
    ∑_i=1^N A_42(i)/z-a_i+m_fc= (m_2-2) λ/π [] ∫_Ω_0(z') ( ∂ I(z,z')/∂ z (1-I_s(z)) ) ^2 z' +O_≺( 1/√(N η_0)).

	
	Finally, plugging (<ref>) in the leading term of (<ref>) we have
	
    ∑_i=1^N D_2(i)/z-a_i+m_fc(W_4-3)  λ/π[] ∫_Ω_0(z') (  ∂ I(z,z')/∂ z (1-I_s(z))  I(z,z') )^2 z'+O_≺( 1+|λ|^2/√(N η_0)).

	
	
	Therefore, we have
	
    (1-I_s(z))  [ ( G-  G)] =(1-I_s(z)) ıλ/π [] ·∫_Ω_0(z')  K(z,z') ^2 z' +ℰ_2,

	where K(z,z') is given by (<ref>) and 
	
    ℰ_2=(1+|λ|^4) [ O_≺(  N Ψ(z)Θ(z) )+O_≺( ( √(N)Ψ^2(z) )+O_≺( Ψ(z)/√(η_0))+O_≺( Θ(z)/η_0) +O_≺( Θ(z)/η)].
 
Dividing both sides by 1-I_s(z), recalling from Lemma <ref> that |1/1-I_s(z)| ∼1/√(κ+η), and using (<ref>), (<ref>), we have completed the proof of Lemma <ref>.	






§ PROOF OF THEOREM <REF> AND THEOREM <REF>



In this section, we compute the variances of the mesoscopic CLT in the bulk and at the edges.



 §.§ In the bulk

We compute the variance V(f) defined in (<ref>) with f given in (<ref>). 

Under the assumptions and notations of Theorem <ref>,  we have
	
    lim_N →∞ V(f)=1/2 βπ^2∫_∫_(g(x_1)-g(x_2))^2/(x_1-x_2)^2 x_1  x_2.



Assuming that we have proved the above lemma, V(f) converges to some positive constant since g ∈ C^2_c(). Theorem <ref> is a direct result of Proposition <ref> after integrating ϕ'(λ) and using the Arzelá-Ascoli theorem and the Lévy continuity theorem.


	Recall that
	
    V(f)=1/π^2∫_Ω_0∫_Ω_0(z_1) (z_2) (K_1+K_2+K_3)^2z_1 ^2z_2,

	where 
	
    K_1   =(m_2-2/β) I;   K_2=(W_4-1-2/β)  ( II +II);
    
    	K_3   =2/β∂/∂ z_1(1/1- I∂  I/∂ z_2)=2/β( 1/1-I(z_1,z_2)I+ 1/(1-I(z_1,z_2))^2II).

	Using Lemma <ref> and the stability bound (<ref>), we have
	
    I(z_1,z_2)   =1/N∑_i=1^N 1+m'_fc(z_1)/(a_i-z_1-m_fc(z_1))^2(a_i-z_2-m_fc(z_2)) =O( 1/√(κ_1+η_1));
     I(z_1,z_2)   =1/N∑_i=1^N 1+m'_fc(z_2)/(a_i-z_1-m_fc(z_1))(a_i-z_2-m_fc(z_2))^2 =O ( 1/√(κ_2+η_2));
     I(z_1,z_2)   =1/N∑_i=1^N (1+m'_fc(z_1))(1+m'_fc(z_2))/(a_i-z_1-m_fc(z_1))^2(a_i-z_2-m_fc(z_2))^2 =O ( 1/√((κ_1+η_1)(κ_2+η_2))).

	In addition, recalling (<ref>), for z_1 ≠ z_2, we have
	
    1/1-I(z_1,z_2)  = 1+m_fc(z_1)-m_fc(z_2)/z_1-z_2.

	If z_1 and z_2 are in the same half plane, m_fc is analytic in a neighborhood of the segment connecting z_1 and z_2, denoted as L(z_1,z_2).  By Lemma <ref>, then we have
	
    | 1/1-I(z_1,z_2)| ≤ 1+| m_fc(z_1)-m_fc(z_2)/z_1-z_2| ≤sup_z ∈ L(z_1,z_2) |m'_fc(z)| ≤ C sup_z ∈ L(z_1,z_2)(1/√(κ+η)).

	If z_1, z_2 belong to different half planes, using Lemma <ref>, then we have
	
    | 1/1-I(z_1,z_2)| ≤ 1+| m_fc(z_1)-m_fc(z_2)/z_1-z_2|  ≤C/|z_1-z_2|≤C/|η_1|+|η_2|.

	Now, we are ready to compute V(f). Since K_i(z,z')=  K_i(z,z')=0, (i=1,2,3), and by Stokes' formula, we have
	
    V(f)=-1/4 π^2∫_Γ_1∫_Γ_2f̃(z_1) f̃(z_2) (K_1+K_2+K_3)  z_1  z_2:=V_1+V_2+V_3,

	where Γ_1={x_1+ı y_1: |y_1|= N^-τη_0 } and Γ_2={x_2+ı y_2: |y_2|= 1/2 N^-τη_0 }. We choose the orientation of both contours to be counterclockwise. The parts on the upper half plane are denoted as Γ_1^+, Γ_2^+, while the parts on the lower half plane are Γ_1^-, Γ_2^-.
	
	Using (<ref>)-(<ref>), since κ≥κ_0 ≥ c_0 for some positive constant c_0>0, we have |K_1+K_2| = O(1).
	Combining with (<ref>), by direct computation, we have |V_1+V_2| =O(η_0^2). It then suffices to estimate V_3. We consider two cases.
	
	Case 1: If z_1, z_2 are in the same half plane, by (<ref>) and (<ref>)-(<ref>), we have  |K_3| = O(1). Therefore, 
	
    ( ∫_Γ_1^+∫_Γ_2^+ +∫_Γ_1^-∫_Γ_2^-) (z_1) (z_2) K_3(z_1,z_2)  z_1  z_2=O(η^2_0).

	
	Case 2: Consider z_1, z_2 are in different half planes. For notational simplicity, we define m_1=m_fc(z_1) and m_2=m_fc(z_2).
	Differentiating I given in (<ref>), we have 
	
    I=  (z_1-z_2) m'_1-m_1+m_2/(z_1+m_1-z_2-m_2)^2;     I=(z_2-z_1) m'_2+m_1-m_2/(z_1+m_1-z_2-m_2)^2.

	Using (<ref>) (<ref>), (<ref>)-(<ref>) and Lemma <ref>, we have
	
    K_3 =2/β1/(z_1-z_2)^2((z_1-z_2) m'_1-m_1+m_2)((z_2-z_1) m'_2+m_1-m_2)/(z_1+m_1-z_2-m_2)^2+O( η^-1_0N^τ).

	Note that if z ∈^+ and in the bulk,  then there exists k, K >0 such that k ≤ m_fc(z) ≤ K. 
	If z_1, z_2 are in different half planes, there exists some constant c>0 such that |z_1+m_1-z_2-m_2|>c. Combining with Lemma <ref>, we have
	
    K_3=-2/β(m_1-m_2)^2/(z_1-z_2)^2(z_1+m_1-z_2-m_2)^2 +O_≺( η^-1_0N^τ)+O(1)=-2/β1/(z_1-z_2)^2 +O( η^-1_0N^τ).

	
Therefore, recalling the definition of f̃ in (<ref>), by symmetry and (<ref>), we have
	
    V_3=1/βπ^2∫_Γ_1^+∫_Γ_2^-f̃(z_1) f̃(z_2) /(z_1-z_2)^2 z_1  z_2+O(η_0N^τ),

	with opposite integral directions on the contours. Since Γ_1^+ and Γ_2^- are disjoint and  has compact support, we obtain from Cauchy's integral theorem that
	
    1/βπ^2∫_Γ_1^+∫_Γ_2^-f̃(z_2)^2/(z_1-z_2)^2 z_1  z_2= 1/βπ^2∫_Γ_2^-f̃(z_2)^2  ( ∫_Γ_1^+1/(z_1-z_2)^2 z_1)  z_2=0.

	The integral of f̃(z_1)^2/(z_1-z_2)^2 vanishes similarly. Thus, we have from (<ref>) and (<ref>) that
	
    V_3=   -1/2 βπ^2∫_Γ_1^+∫_Γ_2^-(f̃(z_1)- f̃(z_2)^2 /(z_1-z_2)^2 z_1  z_2+O(η_0N^τ)
    
    	=   1/2 βπ^2∫_∫_(f(x_1) -f(x_2)+ı N^-τη_0 (f'(x_1)- f'(x_2)))^2/(x_1-x_2+3 /2 N^-τη_0)^2 x_1  x_2+O(η_0N^τ).

	Changing the variable 
	
    x̃_1=x_1-E_0/η_0;     x̃_2=x_2-E_0/η_0,

	we hence obtain from (<ref>) that
	
    V_3   =1/2 βπ^2∫_∫_(g(x̃_1)- g(x̃_2) +ı N^-τ (g'(x̃_1)-g'(x̃_2)))^2/(x̃_1-x̃_2+ 3 ı/2 N^-τ)^2 x̃_1  x̃_2 +O(η_0N^τ).

Note that the integrand can be bounded uniformly by

    |  (g(x̃_1)- g(x̃_2) +ı N^-τ (g'(x̃_1)-g'(x̃_2)))^2/(x̃_1-x̃_2+3 ı/2 N^-τ)^2 | ≤(g(x̃_1)- g(x̃_2) )^2/(x̃_1-x̃_2)^2+(g'(x̃_1)- g'(x̃_2) )^2/(x̃_1-x̃_2)^2.

	Since g ∈ C_c^2(), we then conclude the proof using dominated convergence.





 §.§ Near the edge

Theorem <ref> is a result of the following lemma and Proposition <ref>.

Under the assumptions and notations of Theorem <ref>,  we have
	
    lim_N →∞V(f)=1/ 2 βπ^2∫_∫_(g(-x^2)-g(-y^2)/x-y)^2  x  y.





	Similarly as in the bulk, using Stokes' formula, we have
	
    V(f)=-1/4 π^2∫_Γ_1∫_Γ_2(z_1) (z_2) (K_1+K_2+K_3)  z_1  z_2:=V_1+V_2+V_3,

	where K_1,K_2,K_3 are given by (<ref>) and (<ref>), and we use the same notations and definitions as in the previous subsection. Using (<ref>)-(<ref>) and Lemma <ref>, we have 
	
    |K_1+K_2| =O( 1/√(|z_1 z_2)|)=O(N^τη_0^-1).

	Using (<ref>), one can show |V_1+V_2|=O(η_0N^τ). Thus it is sufficient to study the integral involving K_3. Using (<ref>), (<ref>), and 
	
    ∂^2/∂ z ∂ z' I(z,z')=∂^2/∂ z' ∂ z I(z,z')=(m_1'+m_2'+2m_1'm_2')(z_1-z_2)-(m_1'+m_2'+2)(m_1-m_2)/(z_1+m_1-z_2-m_2)^3,

	by direct computation, we have
	
    K_3=2/β( (1+m_1')(1+m_2')/(z_1+m_1-z_2-m_2)^2  -1/(z_1-z_2)^2).

	
	For the second integrand, using the similar arguments as in the previous subsection, we have
	
    lim_N →∞∫_Γ_1^±∫_Γ_2^±(z_1) (z_2)/(z_1-z_2)^2 z_1  z_2=   -1/2lim_N →∞∫_Γ_1^±∫_Γ_2^±((z_1) -(z_2))^2/(z_1-z_2)^2 z_1  z_2
    
    	=   -1/2∫_∫_(g(x_1)-g(x_2))^2/(x_1-x_2)^2 x_1  x_2.

	Due to the opposite integral directions, we have
	
    lim_N →∞∫_Γ_1^±∫_Γ_2^∓(z_1) (z_2)/(z_1-z_2)^2 z_1  z_2=1/2∫_∫_(g(x_1)-g(x_2))^2/(x_1-x_2)^2 x_1  x_2.

	The whole integral with respect to the second term of K_3 will hence vanish when N →∞. Thus it is sufficient to study the integral of the first term (1+m_1')(1+m_2')/(z_1+m_1-z_2-m_2)^2, that is,
	
    V_3(f)   =-1/2 βπ^2( ∫_Γ_1^+∫_Γ_2^++∫_Γ_1^-∫_Γ_2^-+∫_Γ_1^+∫_Γ_2^-+∫_Γ_1^-∫_Γ_2^+) (z_1) (z_2) (1+m_1')(1+m_2')/(z_1+m_1-z_2-m_2)^2 z_1  z_2
       :=V_3^+++V_3^–+V_3^+-+V_3^-+.

	
	
	
	
	
	
	Let ζ=z+m_fc(z) and ζ_±=L_±+m_fc(L_±) ∈. Define F(ζ):=ζ-1/N∑_i=1^N 1/a_i-ζ so that (<ref>) is equivalent to z=F(ζ). Assumptions 2.2 and 2.3 imply that there exists some constant c_0 independent of N such that
	
    ( {ζ_±,ℐ̂} ) ≥ c_0,

	for all sufficiently large N, where ℐ̂ is the smallest interval that contains the support of μ_A; see (<ref>). Hence, F(ζ) is analytic in a neighborhood of ζ_+, where we write
	
    F(ζ)=F(ζ_+)+F'(ζ_+)(ζ-ζ_+) +F”(ζ_+)/2 (ζ-ζ_+)^2 +O(|ζ-ζ_+|^3).

	By (<ref>), F'(ζ_+)=1-1/N∑_i=1^N 1/(a_i-ζ_+)^2=0. Moreover, F”(ζ_+)=-2/N∑_i=1^N 1/(a_i-ζ_+)^3, and by (<ref>) it is bounded uniformly from below. In general, we have |F^(k)(ζ_+)|= | k!/N∑_i=1^N 1/(a_i-ζ_+)^(k+1)| = O(1) because of (<ref>). Inverting F(ζ)=z in the neighborhood of ζ_+, we have the expansion
	
    ζ=z+m_fc(z)=ζ_++c_+ √(z-L_+)(1+A_+(√(z-L_+)) ),

	where A_+ is an analytic function depending on N with A_+(0)=0. This has been shown in the proof of Lemma 3.6 and  Lemma A.1 in <cit.>.
	Note that c_+= (- 1/N∑_i=1^N 1/(a_i-ζ_+)^3)^-1/2 is some positive number depending on N but is uniformly bounded. Furthermore, the coefficients of the expansion of A_+ are also uniformly bounded. Thus
	
    z+m_fc(z)=ζ_++c_+ √(z-L_+) +O(|z-L_+|),

	where the square root is taken in a branch cut such that √(z-L_+)>0 as z>0. Similarly, we have
	
    1+m'_fc(z)=c_+ /2√(z-L_+)+d_++O(√(|z-L_+|)),

	where d_+ is some number which depends on N but is uniformly bounded.
	Let z=L_++η_0 x+ı N^-τη_0. Then
	
    z+m_fc(z)=ζ_++ c_+√(η_0 (x + ı N^-τ ) ) +O(η_0);  1+m'_fc(z)=c_+/2√(η_0 (x+ ı N^-τ))+d_++O(√(η_0 )).

Therefore, after changing the variable as in (<ref>), we have 
	
    V_3^++   =-1/8 βπ^2∫_∫_(x_1) (x_2) (1/√(x_1+ı N^-τ)+O(√(η_0)) )(1/√(x_2+ı/2 N^-τ)+O(√(η_0)) )/(√( x_1+ı N^-τ)-√( x_2+ı/2 N^-τ)+O( √(η_0)))^2 x_1  x_2
       =-1/8 βπ^2∫_∫_(x_1) (x_2)/√(x_1+ı N^-τ)√(x_2+ı/2 N^-τ)(√( x_1+ı N^-τ)-√( x_2+ı/2 N^-τ))^2 x_1  x_2+O(√(η_0)N^3τ),

	where (x)=g(x)+ı N^-τ g'(x). The last step follows from the fact that | √(x_1+ı N^-τ)-√(x_2+ı/2 N^-τ)| ≥ C N^- τ, when x_1,x_2 belong to some compact set.
	Let γ_1^±:={ x_1 ±ı N^-τ: x_1 ∈} and γ_2^±:={ x_2 ±ı/2 N^-τ: x_2 ∈}.
	Then we obtain
	
    V_3^++=-1/8 βπ^2∫_γ_1^+∫_γ_2^+(z_1) (z_2)/√(z_1)√(z_2) (√( z_1)-√( z_2))^2 z_1  z_2+O(√(η_0 ) N^3τ ),

	where (x + ı y)=g(x)+ı y g'(x)χ(y). Since γ_1^+ and γ_2^+ are disjoint and  has compact support, changing the variable w=√(z) and using Cauchy's integral theorem, we have
	
    ∫_γ_1^+∫_γ_2^+(z_1)^2 /√(z_1)√(z_2) (√( z_1)-√( z_2))^2 z_1  z_2=0=∫_γ_1^+∫_γ_2^+(z_2)^2 /√(z_1)√(z_2) (√( z_1)-√( z_2))^2 z_1  z_2,

	and thus
	
    V_3^++=1/16 βπ^2∫_γ_1^+∫_γ_2^+ ((z_1) -(z_2))^2/√(z_1)√(z_2) (√( z_1)-√( z_2))^2 z_1  z_2+O(√(η_0 ) N^3 τ).

	Therefore, we get
	
    lim_N →∞ V_3^++ =1/16 βπ^2lim_N →∞∫_∫_ (g(x_1)-g(x_2) +ı N^-τ g'(x_1)- ı/2 N^-τ g'(x_2))^2 /√(x_1+ı N^-τ)√(x_2+ı/2 N^-τ)(√( x_1+ı N^-τ)-√( x_2+ı/2 N^-τ))^2 x_1  x_2.

	We denote the integrand as h_N(x_1,x_2).
	Next, we interchange the limit and the integral . One shows that there exists C>0 such that
	
    |√( x_1+ı N^-τ)-√( x_2+ı/2 N^-τ)| ≥ C | √(x_1)-√(x_2)|.

	Set
	
    h(x_1,x_2):= C^-2(g(x_1)-g(x_2))^2+(g'(x_1)-g'(x_2))^2/√(|x_1|)√(|x_2|)|√(x_1)-√(x_2)|^2,

	and observe that |h_N(x_1,x_2)| ≤ h(x_1,x_2). Next, we will show that h(x_1,x_2) is integrable.
	
	Suppose supp(g)⊂[-M,M] for some M>0. 
	Then if x_1 and x_2 are both in [-2M,2M] then we have the following estimation.

Case 1: If x_1, x_2 have the same sign, then
		
    h(x_1,x_2)=   (g(x_1)-g(x_2))^2+(g'(x_1)-g'(x_2))^2/√(x_1)√(x_2) ( √(x_1) -√(x_2))^2
    
    		=   1/√(|x_1|)√(|x_2|)((g(x_1)-g(x_2)/x_1-x_2)^2+(g'(x_1)-g'(x_2)/x_1-x_2)^2)(√(|x_1|)+√(|x_2|))^2
    ≤   8M/√(|x_1|)√(|x_2|)(g'_∞^2+g”_∞^2).

		
Case 2: If x_1 and x_2 are of opposite signs, using |x_1-x_2|=(√(|x_1|) -ı√(|x_2|))(√(|x_1|) +ı√(|x_2|)), we have
		
    h(x_1,x_2)=   (g(x_1)-g(x_2))^2+(g'(x_1)-g'(x_2))^2/√(|x_1|)√(|x_2|)| √(|x_1|) -ı√(|x_2|)|^2
    
    		=   1/√(|x_1|)√(|x_2|)((g(x_1)-g(x_2)/x_1-x_2)^2+(g'(x_1)-g'(x_2)/x_1-x_2)^2)|√(|x_1|)+ı√(|x_2|)|^2
    ≤   8M/√(|x_1|)√(|x_2|)(g'_∞^2+g”_∞^2).

If x_1∉[-2M,2M], then x_2∈[-M,M] otherwise h(x_1,x_2)=0. So for (x_1,x_2)∈[-2M,2M]^c×[-M,M],
	
    h(x_1,x_2) ≤4g_∞^2+4g'_∞^2/√(|x_1|)√(|x_2|)(√(|x_1|)-√(|x_2|))^2≤4g_∞^2+4g'_∞^2/√(|x_1|)√(|x_2|)(√(|x_1|)-1/√(2)√(|x_1|))^2=C/|x_1|^3/2|x_2|^1/2.

	Therefore, h(x_1,x_2) is integrable. Thus by dominated convergence, 
	
    lim_N →∞ V_3^++   =1/16 βπ^2∫_∫_ (g(x_1)-g(x_2))^2 /√(x_1+ı 0)√(x_2+ı 0) (√( x_1+ı 0)-√( x_2+ı 0))^2 x_1  x_2
       =1/4 βπ^2∫_ψ(+ı 0)∫_ψ(+ı 0)(g(w^2_1)-g(w^2_2))^2/(w_1-w_2)^2 w_1  w_2,

	where we change the variable ψ(z):=√(z); with branch cut such that ψ: ^+ →^+.
	
	Similarly, we have
	
    lim_N →∞ V_3^–=1/4 βπ^2∫_ψ(-ı 0)∫_ψ(-ı 0)(g(w^2_1)-g(w^2_2))^2/(w_1-w_2)^2 w_1  w_2;

	
    lim_N →∞ V_3^+-=1/4 βπ^2∫_ψ(+ı 0)∫_ψ(-ı 0)(g(w^2_1)-g(w^2_2))^2/(w_1-w_2)^2 w_1  w_2;

	
    lim_N →∞ V_3^-+=1/4 βπ^2∫_ψ(-ı 0)∫_ψ(+ı 0)(g(w^2_1)-g(w^2_2))^2/(w_1-w_2)^2 w_1  w_2.

	
	The contours are shown in Figure <ref>. 
	
	Note that the horizontal parts of the blue and the red lines of above graph will cancel because of the opposite integral direction. To sum up, we have
	
	
    lim_N →∞ V_3=1/4 βπ^2∫_-ı∞^ı∞∫_-ı∞^ı∞(g(w^2_1)-g(w^2_2))^2/(w_1-w_2)^2 w_1  w_2=1/4 βπ^2∫_∫_(g(-x_1^2)-g(-x_2^2)/x_1-x_2)^2  x_1  x_2.
 
	This concludes the proof of Theorem <ref>.
	




§ PROOF OF PROPOSITION <REF> AND COMPUTATION OF THE BIAS

In this section, we first prove Proposition <ref>, using the same technique as in Proposition <ref>. After this, we compute the bias on mesoscopic scales inside the bulk and at the edges.

We treat the expectation similarly using the cumulant expansion and (<ref>):

    (z-a_i)  G_ii = (HG)_ii-1=1/N∑_j=1^N c^(2)_ij∂ G_ji/∂ H_ij-1 +1/2! N^3/2∑_j=1^N c^(3)_ij∂^2 G_ji/∂^2 H_ij +1/3! N^2∑_j=1^N c^(4)_ij∂^3 G_ji/∂^3 H_ij +O_≺(N^-3/2)


    =   -1/N∑_j=1^N  G_ii G_jj -1/N (G^2)_ii -m_2-2/N (G_ii)^2-1+1/2 N^3/2∑_j=1^N c^(3)_ij(6 G_ii G_ij G_jj+2 G^3_ij)
       +1/6 N^2∑_j=1^N (W_4-3) (-36 G_ii G_jj G^2_ij-6 G^2_iiG^2_jj-6 G^4_ij)+O_≺(N^-3/2).

Combining with the local law, we have

    (z-a_i)  G_ii=-1/N G_ii G -1/Nd/dz1/a_i-z-m_fc -m_2-2/N1/(a_i-z-m_fc)^2-1


    +3/ N^3/2∑_j=1^N   c^(3)_ij/(a_i-z-m_fc)(a_j-z-m_fc) G_ij-1/ N^2∑_j=1^N (W_4-3)  1/(a_i-z-m_fc)^2(a_j-z-m_fc)^2+O_≺( Ψ(z)/N η).



Using the anisotropic local law and the argument as in (<ref>), one can show that the second term of the last line of above equation is O_≺( N^-1Ψ(z)).
Therefore, we have

    (z-a_i)  G_ii=-1/N( G_ii-1/a_i-z-m_fc)  G -1/N1/a_i-z-m_fc G-1-1/N1+m'_fc(z)/(a_i-z-m_fc)^2


    -m_2-2/N1/(a_i-z-m_fc)^2-1/ N  I_s(z) W_4-3/(a_i-z-m_fc)^2+O_≺( Ψ(z)/N η),

and thus

    (z-a_i+m_fc)(  G_ii-1/a_i-z-m_fc)   =-1/N1/a_i-z-m_fc(  G-N m_fc)
    
    -1/N1+m'_fc(z)/(a_i-z-m_fc)^2   -m_2-2/N1/(a_i-z-m_fc)^2-1/ N  I_s(z) W_4-3/(a_i-z-m_fc)^2+O_≺( Ψ(z)/N η).

Dividing both sides by a_i-z-m_fc∼ O(1) and summing over i, we obtain

    (1-I_s(z)) ( G -N m_fc)=   1/N∑_i=1^N 1+m'_fc(z)/(a_i-z-m_fc)^3 +m_2-2/N∑_i=1^N 1/(a_i-z-m_fc)^3
       +W_4-3/ N∑_i=1^N I_s(z)/(a_i-z-m_fc)^3+O_≺( Ψ(z)/η).

Dividing both sides by 1-I_s(z)  and using the relation 1-I_s(z)=1/1+m'_fc(z)∼√(κ+η), we obtain

    ( G-N m_fc)=1/1-I_s(z)1/2 I_s(z)/ z +m_2-2/2 I_s(z)/ z +W_4-3/2 I_s(z) I_s(z)/ z+O_≺( 1/η√(N η√(κ+η))).

Plugging into (<ref>) (here we replace μ_N by μ_fc), using Lemma <ref> and Stokes' formula, we have

    f(X_N)-N ∫_ f(x) ρ_fc(x)  x=1/4 πı∫_∂Ω_0(z) b(z)  z+O_≺( N^2 τ/√(N η_0 √(κ_0+η_0)))+O_≺(N^-τ),

where b(z) is given by (<ref>).
Using the relation I_s=m'_fc/1+m'_fc, it coincides with the expectation that obtained in the global CLT given in Theorem <ref>. Thus we conclude the proof of Proposition <ref>.



Next, we explicitly compute the bias in the bulk and at the edges, for the scaled test function in (<ref>). 



 §.§ Bias in the mesoscopic bulk

Note that

    dI_s/dz=2/N∑_i=1^N 1+m'_fc/(a_i-z-m_fc)^3= O( 1/√(κ+η));      1-I_s(z) ∼1/√(κ+η);      |I_s(z)| = O(1).


If κ≥κ_0> c>0, then |b(z)|=O(1). In combination with (<ref>), we have 

    f(X_N)-N ∫_ f(x) ρ_fc(x)  x =   O_≺(η_0+N^2τ/√(N η_0 √(κ_0+η_0))+N^-τ),

hence we see that the bias vanishes as N goes to infinity.



 §.§ Bias at the mesoscopic edge

Similarly, using (<ref>) and (<ref>), the last two terms of b(z) will contribute O_≺(√(N^τη_0)). We have

    f(X_N)-N ∫_ f(x) ρ_fc(x)  x=1/4 πı∫_∂Ω_0(z) m_fc”/1+m_fc' z+   O_≺(N^-τ+ N^2τ/√(N η_0 √(κ_0+η_0))+√(N^τη_0)).

Using (<ref>), we obtain the following expansions:

    1+m'_fc(z)=c_+ /2√(z-L_+)+O(1),   m”_fc(z)=-c_+ /4(√(z-L_+))^3+ O ( 1/√(|z-L_+|)),

and then

    m_fc”/1+m_fc' =-1/2(z-L_+)+O( 1/√(|z-L_+|)).

Changing variables and using (<ref>), we have

    f(X_N)-   N ∫_ f(x) ρ_fc(x)  x=-1/8πı∫_ (g(x)+ı N^-τg'(x)) 1/x+ı N^-τ x
       +1/8πı∫_ (g(x)-ı N^-τg'(x)) 1/x-ı N^-τ x+O_≺(N^-τ+ N^2τ/√(N η_0 √(κ_0+η_0))+√(N^τη_0))
    
    =   -1/8πı∫_g(x)/x+ı N^-τ x+1/8πı∫_g(x)/x-ı N^-τ x+O_≺(N^-τ+ N^2τ/√(N η_0 √(κ_0+η_0))+√(N^τη_0)).

Using the Sokhotski-Plemelj lemma, we have

    f(X_N)-N ∫_ f(x) ρ_fc(x)  x=g(0)/4+O_≺(N^-τ+ N^2τ/√(N η_0 √(κ_0+η_0))+√(N^τη_0)),

where we used the regularity g ∈ C_c^2(). This finishes the computation of mesoscopic bias.




§ SAMPLE COVARIANCE MATRIX


In this section, we use the previous arguments to derive the mesoscopic eigenvalue statistics of sample covariance matrix and prove similar CLTs in the bulk and at the regular edges. We start by introducing the model in detail.



 §.§ Setup, assumptions and main results

 Let X_N=(X_ij) be an M × N matrix satisfying the following assumption.



  * {X_ij| 1 ≤ i ≤ M, 1 ≤ j ≤ N} are independent real-valued centered random variables.

  * For all i,j, we have |√(N)X_ij|^2=1. In addition, √(N)X_ij has uniformly bounded moments, that is, there exists C_p>0 independent of N such that for all i,j,

    |√(N) X_ij|^p≤ C_p.


  * To simplify the statement,  we also assume that there exists a constant K_4 such that

    K_4:=1/N∑_j=1^N c^(4)_ij,



Note that M depends on N and set 

    γ≡γ_N:=M/N→γ_0,       0< γ_0 <∞.

We study the M × M sample covariance matrices 

    H_M:=Y_NY^*_N,     Y_N:=Σ^1/2 X_N,     Σ:=(σ_i),

where Σ is an M × M positive definite, deterministic and diagonal matrix with 

    ∞ > σ_1 ≥σ_2 ≥⋯≥σ_M > 0,     lim sup_N →∞σ_1 <∞,     lim inf_N →∞σ_M>0.

We denote by μ_Σ the empirical eigenvalue distribution of Σ, i.e. μ_Σ:=1/M∑_j=1^M δ_σ_j. The following assumption ensures that the limit of μ_Σ exists.

Together with (<ref>), μ_Σ  converges weakly to a deterministic measure μ_σ as N →∞ such that μ_σ is compactly supported in (0, ∞).

 The eigenvalues of H ≡ H_M are denoted as λ_i ∈, 1 ≤ i ≤ M. The empirical spectral measure of H_M is defined by μ_M=1/M∑_i=1^M δ_λ_i. 
 The Stieltjes transform of μ_M is then given by

    m_M(z):= M^-1 G(z),    G(z):=(H_M-zI)^-1,    z ∈^+.

We further introduce the N × N matrices

    ℋ_N:=Y_N^* Y_N,     𝒢:=(ℋ-z)^-1,     𝔪_N(z):=N^-1𝒢,    z ∈^+.

The eigenvalues of ℋ≡ℋ_N are denoted by {μ_i}_i=1^N. It is straightforward that {λ_i}_i=1^M differs from {μ_i}_i=1^N by |N-M| zeros, hence we have the relation

    𝔪_N(z)=γ m_M+γ-1/z.



In the null case Σ=I, the Marchenko-Pastur law states that the empirical eigenvalue distribution of H=XX^* converges weakly to the Marchenko-Pastur distribution with aspect ratio γ_0, whose density is given by μ_MP, γ_0:=1/2 πγ_0√([(x-γ_-)(γ_+-x)]_+/x^2) x+(1-γ_0^-1)_+ δ_0 with γ_±=(1 ±√(γ_0))^2. Its Stieltjes transform m_MP, γ_0, or denoted by m_γ_0 for short, is characterized as the unique solution of 

    1+ (z-1+γ_0) m(z) +γ_0 zm^2(z)=0,      m(z)=1/1-γ_0 -γ_0 z m(z)-z,

such that m(z)>0, z ∈^+. Because of (<ref>), the Stieltjes transform of the limiting spectral measure of ℋ=X^*X, denoted by _γ_0^-1, is then given by

    𝔪_γ_0^-1(z)=γ_0 m_γ_0(z)+γ_0-1/z.


In the non-null case Σ≠ I, under Assumption <ref>, the limiting spectral measure of H=Σ^1/2 X X^* Σ^1/2 exists, henceforth referred to as the deformed Marchenko-Pastur law. Its Stieltjes transform, denoted by m_fc, γ_0, or m_fc for short, is the unique solution of

    m(z)=∫_1/t(1-γ_0 -γ_0 z m(z))-zμ_σ(t),

such that m(z)>0, z ∈^+. The corresponding limiting measure, denoted by μ_fc,γ_0 or μ_fc for short, is the free multiplicative convolution of μ_σ and the standard Marchenko-Pastur law with ratio γ_0, i.e., μ_fc, γ_0=μ_σ⊠μ_MP, γ_0; see <cit.>. It was proved in <cit.> that the free multiplicative convolution measure is absolutely continuous and its density function is analytic whenever positive in (0, ∞).

According to (<ref>), the Stieltjes transform of the limiting spectral measure of ℋ=X^* Σ X is the unique solution to

    (z)=1/-z+γ_0 ∫_t/t(z)+1μ_σ(t),     γ_0-1 -z (z)= ∫_γ_0/1+ t (z)μ_σ(t),

such that (z)>0, z ∈^+.


 For general Σ, μ_fc could be supported on several disjoint intervals; we refer to <cit.> for discussions on the support of the density. The following assumption ensures that the free multiplicative convolution measure is supported on a single interval and the edges behave like square roots. It also rules out the possibilities of outliers. 

	Let [σ_-, σ_+] ⊂^+ be the smallest interval that contains the support of μ_σ and set ℐ:=[σ_+^-1, σ_-^-1]. Assume that
	
    inf_x ∈ℐ∫_( t x/1-t x)^2 μ_σ(t)≥γ_0^-1+w,

	for some constant w>0 (the left side may be infinite). Similarly, set ℐ̂:=[ σ_1^-1, σ_M^-1]. Assume that
	
    inf_x ∈ℐ̂∫_( t x/1-t x)^2 μ_Σ(t)≥γ^-1+w,

	for sufficiently large N.


Let ξ=-(z) so that (<ref>) is equivalent to 

    z=F(ξ),      F(ξ):=1/ξ + γ_0 ∫_t /1-t ξμ_σ(t).
 
As an analogue of (<ref>), it was argued in <cit.> that the edges of the support of μ_fc, denoted as E_±, are given by
E_±=F(ξ_±), where ξ_±∈ are the solutions to

    H(ξ):=∫_( t ξ/1-t ξ)^2 μ_σ(t)=γ_0^-1.

Under Assumption <ref>, we have at most two solutions of (<ref>), since H(ξ) is monotone outside [σ_+^-1, σ_-^-1]. Let ξ_+ be the unique solution of (<ref>) in (0, σ^-1_+). The right boundary of the spectrum is given by E_+=F(ξ_+). As for the left edge, we split into three cases. If 0<γ_0<1, there is a unique solution of (<ref>) in the interval (σ_-^-1, ∞), denoted by ξ_-, and the corresponding left edge is given by E_-=F(ξ_-). If γ_0>1, similarly, there is a unique solution of (<ref>) in the interval (-∞,0). These edges are referred to as soft edges. For γ_0=1, the solution does not exist (or say ξ_-=∞), corresponding to E_-=0. This scenario is referred to as the hard edge and the density there goes to infinity at rate of κ^-1/2. In this paper, we only consider the right edge for all 0<γ_0 < ∞. Same discussion easily extends to the soft left edges when γ_0 ≠ 1. 

In addition, Assumption <ref> implies that  

    ( {σ^-1_+, σ^-1_-}, ξ_±) ≥ c_0>0,

provided that ξ_- is finite (γ_0 ≠ 1). The above condition is crucial for the density of the limiting measure to have the square root behavior at the soft edges. It will be proved later in Lemma <ref>. Similar assumptions also appeared in <cit.> for the rightmost edge and <cit.> for all edges in multi-cuts. 


Since the convergence rate of μ_Σ and γ could be very slow, from now on, we work with the finite-N version μ_fc=μ_MP, γ⊠μ_Σ. The corresponding Stieltjes transforms are given by (<ref>) and (<ref>) replacing the limiting measure μ_σ by μ_Σ and γ_0 by γ. The following notations, e.g., m_fc, , E_±, ξ_± and κ are corresponding to μ_MP, γ⊠μ_Σ and are N-dependent.  To be consistent with the previous sections, we will use the tilde sign to denote the ones with respect to μ_MP, γ_0⊠μ_σ. 
The second condition of Assumption <ref> ensures the same properties for the support of μ_MP, γ_0⊠μ_Σ. In particular, for sufficiently large N we have

    min_i{|1- ξ_±σ_i|}≥ c_0>0,

if γ_0≠ 1. If γ_0=1, it only holds true with respect to ξ_+.

Next, we state the local law for the Green function of sample covariance matrix, which is an essential tool in our proof. Let ≡(z) be the unique solution of finite-N version of (<ref>), i.e.,

    γ-1 -z (z)=1/M∑_i=1^M γ/1+ σ_i (z),

such that (z)>0, z∈^+. Define the deterministic control parameters

    Ψ(z):=√(𝔪(z)/N |η|) +1/N |η| ,    Θ(z):=1/N |η| ,      z=E+η∈∖.

We also introduce the spectral domain, for some small c>0,

    S':={z=E+η:  |E| ≤ c^-1, N^-1+c≤η≤  c^-1,|z| ≥ c }.

We further introduce the N+M by N+M matrices

    R:=[ -Σ^-1     X;   X^*    -z ]^-1;     Π:=[ -Σ (1+𝔪Σ)^-1            0;            0         𝔪(z) ];     Σ':=[ Σ 0; 0 I ],    z ∈^+.

Using the Schur decomposition/Feshbach formula, we see that

    R=[ z Σ^1/2GΣ^1/2         Σ X 𝒢;       𝒢 X^* Σ             𝒢 ]=
    [ z Σ^1/2GΣ^1/2       Σ^1/2GY;    Y^* GΣ^1/2             𝒢 ].



We are ready to state the (anisotropic) local law for such random matrix.
(Theorem 2.4 in <cit.>, Theorem 3.6 in <cit.>)
	For any deterministic unit vector v,w ∈^N, we have
	
    | ⟨ v, Σ'^-1( R(z)-Π(z) )Σ'^-1 w ⟩| ≺Ψ(z),
 
	uniformly in z  ∈ S'. It also implies that
	
    |(𝒢(z))_ij-𝔪(z) δ_ij| ≺Ψ(z);     |  (G(z))_ij+1/z(1+𝔪σ_i)δ_ij| ≺Ψ(z).

	In addition, we have the averaged result
	
    | N^-1𝒢(z)-𝔪  (z)| ≺Θ(z),     | M^-1 G(z)+  ∫_1/z(1+𝔪t)μ_Σ(t) | ≺Θ(z).



	
In the following, we state some properties of the Stieltjes transform 𝔪 in (<ref>), whose proofs are given in Appendix B. Define the spectral domain, for some small c>0,

    S:={z=E+η:  |E| ≤ c^-1, 0 < η≤  c^-1,|z| ≥ c }.

And set κ≡κ(E):=min{ |E_+-E|, |E_–E| }.




  * For z ∈ S and sufficiently large N, we have 
	
    |𝔪(z)| ∼ 1;     min_i |1+ σ_i 𝔪(z)|>c_0.
	

  * For z ∈ S and sufficiently large N, we have 
	
    |(z)|  ∼√(κ+η),     E ∈ [E_-,E_+], 
    η/√(κ+η),    .


  * For z ∈ S with E_–c<E< E_++c and η≤ c for some small c>0, we have
	
    1-1/M∑_i=1^M γσ_i/z(1+𝔪(z) σ_i)^2=-(z)/z '(z)∼√(κ+η).


  * Under the same condition as in (3), we have
         
    '(z)=-(z)/z-1/M∑_i=1^M γσ_i/(1+𝔪(z) σ_i)^2∼1/√(κ+η).




We are now prepared to state our main results for the sample covariance matrix.


Consider a sample covariance matrix satisfying Assumption <ref>, <ref> and <ref> and E_0 is chosen to be away from zero, then Proposition <ref> and <ref> hold true with 

    K(z_1,z_2)=2 ( _1' _2'/(_1-_2)^2-1/(z_1-z_2)^2)+K_4 γ( 1/M∑_i=1^M 1/(1+𝔪_1 σ_i)(1+𝔪_2 σ_i)),

where we use _1 and _2 to denote (z_1) and (z_2), and 

    b(z)=(  ('(z))^2/(z) +K_4 (z) '(z)) 1/M∑_i=1^Mγσ_i^2 /(1+(z) σ_i)^3.




Proposition <ref> implies that Theorem <ref> and <ref> hold true for sample covariance matrix. More specifically, we have the following theorem.

	Let H_N be a sample covariance matrix of the form (<ref>) satisfying Assumptions <ref>-<ref>. Let N^-1+c_1≤η_0 ≤ N^-c_1 with some c_1>0 and fix E_0 ∈ (E_-,E_+), such that κ_0:=dist(supp(f_N),{E_±}) > c_0, for some c_0>0 and sufficiently large N. Then, for any function g ∈ C^2_c(), the linear eigenvalue statistics (<ref>)
	converges in distribution to the Gaussian random variable 𝒩(0, 1/π∫_ |ξ| |ĝ(ξ)|^2 ξ),
	where ĝ(ξ):=(2 π)^-1/2∫_ g(x) e^-ıξ x x. 

	In addition, the linear statistics (<ref>) with E_0=E_+ and N^-2/3+c_2≤η_0 ≤ N^-c_2 for some c_2>0, converges in distribution to a Gaussian random variable 𝒩( g(0)/4, 1/2 π∫_ |ξ| |ĥ(ξ)|^2 ξ),
	where h(x)=g(- x^2) and ĥ(ξ):=(2 π)^-1/2∫_ h(x) e^-ıξ x x. 
	Furthermore, if γ_0 ≠ 1, a similar CLT for E_0=E_- can be obtained with h(x)=g(x^2).


Remark: We remark that (<ref>) in Assumption <ref>  can  be removed. In addition, we can relax the single support condition for μ_fc by assuming instead that the cuts of the support of μ_fc are separated by order one and the density has square root behaviors at the edges away from zero.





 §.§ Proof of the CLT and variance computation


	From the definition of the Green function and (<ref>), we get
	
    z G_ii=(HG)_ii-1=(Σ^1/2 X X^* Σ^1/2G)_ii-1=√(σ_i)∑_j=1^N  X_ij (G Y)_ij -1.

	Similarly as (<ref>), by the cumulant expansion formula, we have

    z [(G_ii-  G_ii)]=I_1+I_2+I_3+O_≺(N^-3/2 (1+|λ|^4)),

where

    I_1:   =√(σ_i)/N∑_j=1^N  c^(2)_ij( [∂/∂ X_ij ( G Y)_ij]+[ ( ∂ (G Y)_ij/∂ X_ij - [∂ (G Y)_ij/∂ X_ij] ) ] ),
    
    I_2:   =√(σ_i)/2! N^3/2∑_j=1^N  c^(3)_ij( [∂^2 /∂^2 X_ij (G Y)_ij]+2[∂/∂ X_ij∂ (G Y)_ij/∂ X_ij]+[(∂^2 (G Y)_ij/∂^2 X_ij -[∂^2 (GY)_ij/∂^2 X_ij] ) ] ),
    
    I_3:   =√(σ_i)/3! N^2∑_j=1^N c^(4)_ij( [∂^3  /∂^3 X_ij (G Y)_ij]+3 [∂^2  /∂^2 X_ij∂ (G Y)_ij/∂ X_ij]+3[∂/∂ X_ij∂^2 (G Y)_ij/∂^2 X_ij]
            +[(1-) ( ∂^3 (GY)_ij/∂^3 X_ij)]  ).

The last term on the right side of (<ref>) is estimated by (<ref>), (<ref>) and Lemma <ref>. The argument is similar as in (<ref>). The only things to check are the deterministic bounds of (Y^*GY)_ii and (GY)_ij. Note that from (<ref>) and (<ref>), YG=𝒢Y and  |G_ij| =O(N^c), for z ∈Ω_0 ∩ S', thus we have

    (Y^*GY)_ii=(Y^* Y 𝒢 )_ii=(ℋ𝒢 )_ii=(1+z 𝒢)_ii=O(N^c_1);
       |(GY)_ij| ≤√(N) (G Y Y^* G^*)_ii^1/2=√(N) (z(G G^*)_ii+G_ii^*)^1/2=√(N)(z/2 Imz(G_ii- G^*_ii)+G_ii^*)^1/2=O(N^c_2),

where we use Cauchy-Schwarz inequality and the resolvent identity (<ref>). 

Using the formulas,

    ∂ G_ab/∂ Y_jk=-G_aj (Y^*G)_kb-(GY)_ak G_jb,     ∂ G_ab/∂ X_jk=∂ G_ab/∂ Y_jk√(σ_j),

we obtain the analogue of Lemma <ref>:


For any i,j, we have	

    ∂/∂ X_ij   =- 2 √(σ_i)λ/π∫_Ω_0(z) / z (GY)_ij^2 z;
    ∂^2 /∂^2 X_ij   =- 2 σ_i λ/π∫_Ω_0(z) / z(𝔪/1+𝔪σ_i) ^2 z+O_≺( (1+|λ|)^2/√(N η_0)).

	In general, for any integer k ∈, we have
	
    | ∂^k (GY)_ij/∂ X^k_ij| ≺ O(1);   | ∂^k  /∂^k X_ij| ≺ O((1+|λ|)^k).




We first look at I_1. Using (<ref>) and (<ref>), we have
	
    I_1=   σ_i/N∑_j=1^N   [ (1-) (   G_ii)  ] - σ_i/N∑_j=1^N  [ (1-) ( (Y^*GY)_jj G_ii)  ]
       - σ_i/N∑_j=1^N  [ (1-) ( (GY)_ji(GY)_ji)  ]+√(σ_i)/N∑_j=1^N  [∂/∂ X_ji (GY)_ji] =: A_1+A_2+A_3+A_4.

	Note that
	
    A_1=σ_i [ (1-) G_ii ]=σ_i [ (1-) G_ii ].

 In addition, using the definition of the resolvent and the local law in Theorem <ref>, we have
 
    A_2=-σ_i/N [ (1-)  (   (G H) G_ii)]=-γσ_i  [ (1-)  ( z M^-1 G G_ii)]-γσ_i  [ (1-)  (G_ii)]

 
    =(-γσ_i  z m_fc-γσ_i)  [ (1-)  ( G_ii)] +γσ_i/1+𝔪σ_i [ (1-)  ( M^-1 G  )]+O_≺( Ψ(z) Θ(z) ).

         Next, we use the local law to estimate the third term,
         
    A_3=- 1/N[ (1-) ( Y^* G^2 Y)_ii]=- 1/N[ (1-) / z ( Y^*GY)_ii]=O_≺( Ψ(z)/ N η).

	Finally, we study the last term A_4 using (<ref>), which can be written as

    A_4   =-ı 2 σ_i λ/π N[  ∫_Ω_0∂/∂z_2(z_2) ∂/∂ z_2 ( G(z_2) H G(z_1) )_ii^2 z_2  ]
       =-ı 2 σ_i λ/π N[  ∫_Ω_0∂/∂z_2(z_2) ∂/∂ z_2 (z_2 G(z_2) G(z_1))_ii^2 z_2  ].

Using the resolvent identity (<ref>), the local law Theorem <ref> and Lemma <ref>,
we obtain that
	
    A_4 = -ı 2 λ/π N[  ∫_Ω_0∂/∂z_2(z_2) ∂/∂ z_2( σ_i z_2 (_i(z_1)-_i(z_2))/z_1-z_2)^2 z_2  ]+e(i),
	
	where _i(z):=-1/z(1+𝔪σ_i) for simplicity, and e(i) is the error term.
	If we consider the linear statistics of the error term e(i), using the same argument as for deformed Wigner matrix in Section <ref>, we have
	
    | ∑_i=1^N _i(z) e_i(z) | = O_≺( 1/N η^2_1)+O_≺( 1/N η_0η_1).

	
	
As for the second cumulant expansion term I_2, we apply the same argument as in (<ref>) and the anisotropic law Theorem <ref> to find that
 
    I_2 =O_≺( (1+|λ|^2) Ψ(z)/N √(η_0))+O_≺(  Ψ^2(z)/√(N)).

	
	
We compute the third cumulant expansion term I_3 similarly using (<ref>), the local law Theorem <ref> and (<ref>). The leading term comes from the second term of I_3, denoted by D_2, i.e.,

    I_3=σ_i/2 N^2∑_j=1^N [∂^2  /∂^2 X_ij c^(4)_ij( G_ii-(Y^*GY)_jj G_ii-(GY)_ij (GY)_ij) ]+O_≺ ((1+|λ|^3) N^-1Ψ(z)).

Using (<ref>) and (<ref>), D_2 can be written as

    D_2=   σ_i/2 N^2∑_j=1^N [∂^2  /∂^2 X_ij c^(4)_ij𝔪/1+𝔪σ_i]+O_≺( (1+|λ|^2)N^-1Ψ(z) )
     
    =   -ı  K_4 σ_i^2 λ/π N[  ∫_Ω_0(z') ( 𝔪(z)/1+𝔪(z) σ_i𝔪(z')/1+𝔪(z') σ_i) ^2 z' ] 
       +O_≺( 1+|λ|^2/ N √(N η_0))+O_≺( (1+|λ|^2) N^-1Ψ(z) ),

where K_4 is given in (<ref>) and is independent of the index i.
	

	Summing up and rearranging terms in (<ref>), the coefficient in front of [ (1-) G_ii] is given by
	
    (z-σ_i+γσ_i +γσ_i z m_fc )=z(1+ 𝔪(z) σ_i),

	which is away from zero for z ∈ S'. Dividing both sides of (<ref>) by this coefficient, we have
	
    [ ( G_ii -  G_ii)]=γσ_i/M z(1+𝔪(z) σ_i)^2 [ (1-)  (   G  )]+ A_4(i)/z(1+ 𝔪(z) σ_i)+D_2(i)/z(1+ 𝔪(z) σ_i) +ℰ(i),

	where A_4 and D_2 are given in (<ref>) and (<ref>), and the error terms ℰ(i) is analytic in Ω_0 and estimated as before. Summing over i and rearranging,  we have
	
    ( 1-1/M∑_i=1^M γσ_i/z(1+𝔪(z) σ_i)^2)  [ (1- )  G]=∑_i=1^M  A_4(i)/z(1+ 𝔪(z) σ_i)+∑_i=1^M D_2(i)/z(1+ 𝔪(z) σ_i) +ℰ,

	where the error term ℰ has an upper bound as in (<ref>).	
	Dividing by the coefficient 1-1/M∑_i=1^M γσ_i/z(1+𝔪(z) σ_i)^2 and recalling (<ref>), the first two terms on the right side of (<ref>) are denoted as A and D respectively, and the error term is bounded by (<ref>). Using (<ref>) and (<ref>), we have
	
    A   =-z_1 '_1/_1ı 2 λγ/π[  ∫_Ω_0∂/∂z_2(z_2) ∂/∂ z_2( z_2/z_1-z_21/M∑_i=1^M σ_i _i(z_1)(_i(z_1)-_i(z_2))  )^2 z_2  ]
       =ı 2 λ/π[  ∫_Ω_0∂/∂z_2(z_2)  ( _1' _2'/(_1-_2)^2-1/(z_1-z_2)^2)^2 z_2  ],

	where we set _1:=(z_1) and _2:=(z_2) for short.
	This follows from
	
    z_1z_2 1/M∑_i=1^M σ_i  _i(z_1)  _i(z_2)    =1/M∑_i=1^N σ_i/(1+_1 σ_i) (1+_2 σ_i)=z_1_1-z_2 _2/γ(_1-_2);
    
    	z_1^2 1/M∑_i=1^M σ_i  ^2_i(z_1)     =1/M∑_i=1^N σ_i/(1+_1 σ_i)^2=(z_1 _1 )'/γ_1'.

	Similarly, recalling (<ref>), we obtain that
	
    D   =-z_1 '_1/_1∑_i=1^N  D_2(i)/z(1+ 𝔪(z) σ_i)=ı  K_4 γλ/π[  ∫_Ω_0(z_2)  1/M∑_i=1^M ( σ_i^2 𝔪'_1/(1+𝔪_1 σ_i)^2𝔪_2/1+𝔪_2 σ_i) ^2 z_2 ]
       =ı  K_4 γλ/π[  ∫_Ω_0(z_2)  ( [ 1/M∑_i=1^M 1/(1+𝔪_1 σ_i)(1+𝔪_2 σ_i)] ) ^2 z_2 ].

Therefore, we obtain an analogue of Lemma <ref>, where the kernel is given  instead by (<ref>).


Next, we compute the explicit formula for the variance V(f) given by (<ref>) with K in (<ref>) and test function f in (<ref>). Both in the bulk and at the edge, the second term of (<ref>) will only contribute O(η_0N^τ), because of (<ref>), (<ref>) and (<ref>). It is sufficient to look at the first term. In the bulk, the main contribution of V(f) comes from the term -2/(z_1-z_2)^2 with z_1, z_2 in different half planes. Hence by similar arguments as in Lemma <ref>, we obtain that

    lim _N →∞ V(f)=1/2  π^2∫_∫_(g(x_1)-g(x_2))^2/(x_1-x_2)^2 x_1  x_2.

At the edge, we recall the expansion of (z) near z=E_+ from (<ref>), 

    (z)=(E_+)+c_+ √(z-E_+)(1+A_+(√(z-E_+)))=ξ_++c_+√(z-E_+)+O(|z-E_+|),

where the square root is taken in a branch cut such that >0 when z>0, and c_+>c_0>0 for sufficiently large N. Differentiating it, we have

    '(z)=c_+/2√(z-E_+)+d_++O(√(|z-E_+|)).

Therefore, repeating the arguments in the proof of Lemma <ref>, we get the variance at the edge

    lim _N →∞ V(f)=1/4 π^2∫_∫_(g(-x_1^2)-g(-x_2^2)/x_1-x_2)^2  x_1  x_2.







 §.§ Expectation and bias computation


Starting from (<ref>), using the cumulant expansion and (<ref>), we obtain that
	
    z G_ii=   √(σ_i)∑_j=1^N   X_ij (G Y)_ij -1=√(σ_i)/N∑_j=1^N c^(2)_ij∂ (GY)_ij/∂ X_ij-1 
       +√(σ_i)/2! N^3/2∑_j=1^N c^(3)_ij∂^2 (GY)_ij/∂^2 X_ij +√(σ_i)/3! N^2∑_j=1^N c^(4)_ij∂^3 (GY)_ij/∂^3 X_ij +O_≺(N^-3/2)
    
    =   σ_i/N∑_j=1^N  G_ii-σ_i/N∑_j=1^N  (Y^*GY)_jj G_ii-σ_i/N∑_j=1^N  ((GY)_ij)^2-1
       +σ_i^3/2/2 N^3/2∑_j=1^N c^(3)_ij(-6 G_ii (GY)_ij+6 G_ii(GY)_ij(Y^* G Y)_jj+2 ((GY)_ij)^3)
       +σ_i^2/6 N^2∑_j=1^N c^(4)_ij(-6(G_ii)^2+12 (G_ii)^2 (Y^* G Y)_jj-6 (G_ii)^2((Y^*GY)_jj)^2 )+O_≺(N^-3/2).

The first line (<ref>) can be written as

    σ_i  G_ii-σ_i/N [ (Y^*GY) G_ii]-σ_i/N  (GH G)_ii-1
    
    =   σ_i (1-γ) G_ii-zσ_iγ [M^-1 G G_ii]+σ_i/N( 1/1+σ_i )'-1+O_≺(N^-3/2)
    
    =   σ_i (1-γ)  G_ii-zσ_iγ m_fc( G_ii+1/z(1+ σ_i))+σ_iγ/M(1+σ_i) G-σ^2_i/N'/(1+σ_i )^2-1+O_≺((N η)^-3/2).

Using the anisotropic local law and the same arguments as in (<ref>), one shows that the second line (<ref>) is O_≺( N^-1Ψ(z)). The local law implies that the last line (<ref>) becomes

    σ_i^2/ N^2∑_j=1^N c^(4)_ij(-(G_ii)^2+2 (G_ii)^2 (Y^* G Y)_jj- (G_ii)^2((Y^*GY)_jj)^2 )+O_≺(N^-3/2)


    =-σ_i^2/ N^2∑_j=1^N c^(4)_ij(^2/(1+σ_i)^2)+O_≺(N^-3/2)=-K_4σ_i^2 / N(^2/(1+σ_i)^2)+O_≺(N^-3/2).

Therefore, we have

    z(1+ 𝔪σ_i) (  G_ii+1/z(1+σ_i))=σ_iγ/1+σ_i [M^-1 G-m_fc]-1/Mγσ^2_i '/(1+σ_i )^2
    
    -K_41 / Mγσ_i^2 ^2/(1+σ_i)^2+O_≺((N η)^-3/2).

Dividing both sides by z(1+ 𝔪σ_i) ∼ O(1) from (<ref>) and summing over i, we obtain

    (1-1/M∑_i=1^M γσ_i/z(1+𝔪(z) σ_i)^2) ( G -M m_fc)=  -  1/M∑_i=1^Mγσ_i^2 ' /z(1+σ_i)^3-K_4  1/M∑_i=1^M γσ_i^2 ^2/z(1+σ_i)^3 +O_≺(1/√(N η^3)).

Dividing both sides by the coefficient of ( G -M m_fc) and using (<ref>), the error becomes O_≺( 1/η√(N η√(κ+η))) and the leading term on the right side becomes (<ref>).
Therefore, we obtained an analogue of Proposition <ref> with the integral kernel b(z) given instead by (<ref>). 

Finally, we compute the explicit formula for the bias. Using (<ref>), (<ref>) and (<ref>), the second term of b(z) will contribute O(√(η_0 N^τ)) both in the bulk and at the edge. In addition, (<ref>) implies that the first term of b(z) can be written as

    (')^2/1/M∑_i=1^M γσ_i^2 /(1+σ_i )^3=”/2 '-'/.

The second term on the right side contributes O(√(η_0 N^τ)). The first term vanishes if κ_0>c>0 and thus the bias in the bulk vanishes. At the edge, we use the expansion of (z) around E_+ from (<ref>), 

    '(z)=c_+/2√(z-E_+)+d_++O(√(|z-E_+|)),     ”(z)=-c_+/4(√(z-E_+))^3+O(1/√(|z-E_+|)).

Hence we have

    ”/2 '=- 1/4(z-E_+)+O(1/√(|z-E_+|)),

and using the Sokhotski-Plemelj lemma, the bias at the edge becomes g(0)/4.








§ COMPLEX CASE


In this appendix, we extend previous results from real symmetric to complex Hermitian matrices. We will use the complex analogue of Lemma <ref>. 
(Complex cumulant expansion)
	Let h be a complex-valued random variable with finite moments, and f is a complex-valued smooth function on  with bounded derivatives. Let c_p,q be the (p,q) cumulant of h, which is defined as 
	
    c_p,q:=(-)^p+q( ∂^p+q/∂ s^p ∂ t^qlog e^ı s h+ı t h) |_s,t=0.
 
	Then for any fixed l ∈, we have
	
    [ h̅ f(h, h̅)]=∑_p+q=0^l 1/p!q! c_p,q+1(h)[ f^(p,q)(h)] +R_l+1,

	where the error term satisfies
	
    |R_l+1| ≤ C_l  |h|^l+2max_p+q=l+1{sup_|x| ≤ M |f^(p,q)(z,z̅)| }+C_l [ |h|^l+2 1_|h|>M] max_p+q=l+1f^(p,q)(z,z̅)_∞,

	and M>0 is an arbitrary fixed cutoff.

Instead of (<ref>), we have

    ∂ G_ij/∂ H_ab=-G_iaG_bj,

from which we obtain the analogue of Lemma <ref>.

The assumption H_ij^2=0 implies that c_ij^(1,1)=1, c_ij^(2,2)=W_4-2 for i ≠ j.  Using the anisotropic law and (<ref>), one shows similarly that the expansion terms corresponding to p+q=3 are negligible. Using (<ref>) and the analogue of Lemma <ref>, we obtain that
         
    (z-a_i)  [(G_ii- G_ii)]=   1/N∑_j=1^N c^(1,1)_ij[∂/∂ H_ji(  (G_ji- G_ji) ) ] 
       +1/2! N^2∑_j=1^N c^(2,2)_ij[∂^3 /∂^2 H_ji∂ H_ij(  (G_ji- G_ji) ) ] +⋯
    
    =   1/N[( ∂ G_ji/∂ H_ji- ∂ G_ji/∂ H_ji) ]+1/N∑_j=1^N  [ ∂/∂ H_ji G_ji]+m_2-1/N[ ∂/∂ H_ii G_ii]
       +1/N^2∑_j=1^N (W_4-2) [ ∂/∂ H_ji∂ H_ij∂ G_ji/∂ H_ji]+⋯.

Thus Proposition <ref> holds with modified variance, i.e. m_2-2 be replaced by m_2-1, W_4-3 be replaced by W_4-2, and the coefficient of the remaining term be 1 instead of 2. 
Similarly, as for the expectation, 

    (z-a_i)  G_ii = (HG)_ii-1=1/N∑_j=1^N c^(1,1)_ij∂ G_ji/∂ H_ji-1 +1/2! N^2∑_j=1^N c^(2,2)_ij∂^3 G_ji/∂^2 H_ji∂ H_ij +⋯.

Thus the first term of b(z) given in (<ref>) vanishes, m_2-2 is replaced by m_2-1 and W_4-3 is replaced by W_4-2.



§ PROOFS OF AUXILIARY LEMMAS




         For 1 ≤ s ≤ 2, the proof is given in Lemma 4.4 in <cit.>. Since h(z) is holomorphic on Ω_0, (z) h(z)=((z) h(z)). Using Stokes' formula, we have
	
    ∫_Ω_0(z) h(z) ^2z =-ı/2∫_∂Ω_0(z) h(z)  z.

	Since g is compactly support, (z)=0 on ∂Ω_0 except 
    Γ_0:={ x+ı y: x ∈(f),  |y|=N^-τη_0 }.
 Using (<ref>) we have
	
    | ∫_Ω_0(z) h(z) ^2z | ≤ C K ∫_Γ_0(  |y |^-s |f(x)|+ |y|^1-s |f'(x)|  )  z ≤ C' K N^τ sη^1-s_0 .
	



	Using the self-consistent equation of m_fc in (<ref>), we have
	         
    m_fc(z_1)-m_fc(z_2)   =1/N∑_i=1^N ( 1/a_i-z_1-m_fc(z_1)-1/a_i-z_2-m_fc(z_2))
       =1/N∑_i=1^N ( z_1+m_fc(z_1)-z_2-m_fc(z_2)/(a_i-z_1-m_fc(z_1))(a_i-z_2-m_fc(z_2))).

	Dividing z_1+m_fc(z_1)-z_2-m_fc(z_2) on both sides and we get the first identity.
	Taking the derivative of (<ref>), we have
	
    1/N∑_i=1^N 1+m'_fc(z)/(a_i-z-m_fc)^2=m'_fc(z).

	 We treat Ĩ and Ĩ_s similarly. Thus we complete the proof.




Note that
	
    |I_s(z)| ≤1/N∑_i=1^N 1/|a_i-z-m_fc(z)|^2 =m_fc(z)/m_fc(z)+η < 1.

	
	By (<ref>), we have m'_fc/1+m'_fc=I_s(z) and thus m'_fc(z)=I_s(z)/1-I_s(z). Using  Lemma <ref>, we have
	
    |m'_fc(z)| ≤1/|1-I_s(z)|∼1/√(κ+η).

	Differentiating (<ref>) again, we obtain that
	
    m”_fc/(1+m'_fc)^3=2/N∑_i=1^N 1/(a_i-z-m_fc)^3.

	Combining (<ref>) and (<ref>),  we get the upper bound of m”_fc. The rest inequalities follow directly from Lemma <ref>




	Using (<ref>), we have
	
    ∂/∂ H_ij= ıλ/π∫_Ω_0(z) (∑_l=1^N ∂ G_ll/∂ H_ij) ^2 z=-ı (2-δ_ij)λ/π∫_Ω_0(z) (G^2)_ji^2 z.

	Note that (G^2)_ji=d/dzG_ji(z). Since G_ij is analytic in D', using the Cauchy integral formula and the local law, we have that for i ≠ j,  (G^2)_ji≺Ψ(z)/z. Combining with Lemma <ref>, we obtain that, for i ≠ j,
	
    |∂/∂ H_ij|=O_≺( 1+|λ|/√(N η_0)).

	Similarly, if i=j, we have 
	
    ∂/∂ H_ii  = -ıλ/π∫_Ω_0(z) 1/a_i-z-m_fc(z)^2 z+O_≺(1+|λ|/√(N η_0)).

	Furthermore, we compute that
	         
    ∂^2 /∂^2 H_ij=   -λ^2 (2-δ_ij)^2/π^2( ∫_Ω_0(z) (G^2)_ji^2 z  )^2
       +ı (2-δ_ij) λ/π∫_Ω_0(z) (2 (G^2)_ji G_ij+ (1-δ_ij)(G^2)_ii G_jj+(1-δ_ij)(G^2)_jj G_ii) ^2 z.

	For i ≠ j, combining the local law and Lemma <ref>, we have
	
    ∂^2 /∂^2 H_ij   =ı 2 λ/π∫_Ω_0(z) / z(G_ji G_ij+ G_ii G_jj) ^2 z+O_≺( (1+|λ|)^2/N η_0)
       =ı 2 λ/π∫_Ω_0(z) 1/(a_i-z-m_fc)(a_j-z-m_fc)^2 z+O_≺( (1+|λ|)^2/√(N η_0)).

	Similarly, for i=j, we have
	
    ∂^2 /∂^2 H_ii=ıλ/π∫_Ω_0(z) 1/(a_i-z-m_fc)^2^2 z+O_≺( (1+|λ|)^2/√(N η_0)).

	In general, using the local law, (<ref>) and Lemma <ref>, we complete the proof of (<ref>).	




We start by proving the first two statements, using which we will show (<ref>). The last statement then follows directly from (<ref>). Note that the first equation in (<ref>) implies the first inequality in (<ref>). To prove the rest, we divide the spectral domain S into three regimes, corresponding to the bulk, the edge and the outside.
First, we consider z near the edge E_+, i.e., z ∈ S^e:={ z=E+η∈ S: E ∈ [E_+-τ', E_+ +τ']} for some small τ'>0. Let ξ=-(z) so that (<ref>) is equivalent to 

    z=F(ξ),      F(ξ):=1/ξ + γ∫_t /1-t ξμ_Σ(t).
 
Due to (<ref>), F(ξ) is analytic around ξ_+, and we have

    F(ξ)=F(ξ_+)+F'(ξ_+)(ξ-ξ_+)+1/2!F”(ξ_+)(ξ-ξ_+)^2+O(|ξ-ξ_+|^3),

where the linear term vanishes because of (<ref>). It also yields

    F”(ξ_+)=2/ξ_+^3+ 2 γ∫t^3/(1-t ξ_+)^3μ_Σ(t)=2 γ∫t^2/ξ_+(1-t ξ_+)^3μ_Σ(t) ≥ c>0.

The last step follows from the fact that ξ_+ ≥ c>0 and 1-tξ_+ ≥ c>0 for sufficiently large N. Thus we have the expansion of  near the edge z=E_+, 

    (z)=(E_+)+c_+ √(z-E_+)(1+A_+(√(z-E_+)))=ξ_++c_+√(z-E_+)+O(|z-E_+|),

where c_+>c_0>0 for sufficiently large N, A_+ is an analytic function with A_+(0)=0, and the square root is taken with the branch cut such that (z)>0 when z>0. Hence the corresponding density has the square root behavior at the right edge. The left edge can be treated similarly when γ_0 ≠ 1. Using the definition of Stieltjes transform, one shows (<ref>). Similar arguments can be found in Lemma A.5 <cit.>. 

 If E is inside the bulk, i.e., z ∈ S^b:={ z ∈ S: E ∈ [E_-+τ', E_+ -τ']}, then ≥ c>0. Thus (<ref>) and (<ref>) follows.  Finally, for the outside spectral domain, if  z ∈ S^o:={ z ∈ S: (E, [E_-, E_+]) ≥τ'}, it follows from ∼η and (<ref>). Similar arguments can be found in Appendix A in <cit.>. 

Next, we will prove (<ref>). We first prove the upper bound. Taking the real and imaginary part of (<ref>), we have
         
    E + η=γ/M∑_i=1^M σ_i 𝔪/|1+𝔪σ_i|^2;      E - η=-γ/M∑_i=1^M  1+σ_i 𝔪/|1+𝔪σ_i|^2+γ-1,

        with ≡(z). Then we have
         
    | z-γ/M∑_i=1^M σ_i/|1+𝔪(z) σ_i|^2|=| ıη-/η| =|| η/.

         Using |(z)| ∼ 1 and (<ref>). we obtain an upper bound of the right side as C√(κ+η).
        In addition,
        
    | 1/M∑_i=1^M σ_i/|1+𝔪(z) σ_i|^2-1/M∑_i=1^M σ_i/(1+𝔪(z) σ_i)^2|
       =2 |1/M∑_i=1^M σ_i ( (1+𝔪(z) σ_i))^2+σ_i  (1+𝔪(z) σ_i)  (1+𝔪(z) σ_i)/|1+𝔪(z) σ_i|^4|  ≤  C √(κ+η),

        hence we obtain an upper bound for the left side of (<ref>). Next, it is sufficient to show the lower bound. 
         If z ∈ S^e, (<ref>) implies that (1+σ_i )>c if we choose τ' sufficiently small. We split into two cases. If E ∈ [E_-, E_+], we have
         
    |z-1/M∑_i=1^M γσ_i/(1+𝔪(z) σ_i)^2| ≥| ( z- 1/M∑_i=1^M γσ_i/(1+𝔪(z) σ_i)^2) | ≥| η+ 2γ/M∑_i=1^M σ^2_i   (1+σ_i )/|(1+𝔪(z) σ_i)|^4|

         
    ≥2γ/M∑_i=1^M σ_i    (1+σ_i )/|(1+𝔪(z) σ_i)|^4≥ C √(κ+η).

         Otherwise, E ∈ [E_-, E_+]^c, we have
         
    | z-γ/M∑_i=1^M σ_i/(1+𝔪(z) σ_i)^2|    ≥| |z|-γ/M∑_i=1^M σ_i/|1+𝔪(z) σ_i|^2|=| √(E^2+η^2)-E -/η|
       ≥||/η -η≥ C √(κ+η)-η≥ C' √(κ+η),

         if we choose τ' sufficiently small. The last second step follows from the fact that near the edge, || ≥ C>0, because of the second equation of (<ref>).
Next, if z ∈ S^b, then we have >c. We also split into two cases. If > 0, we repeat (<ref>) to get (<ref>). If ≤ 0, from (<ref>) we have

    γ/M∑_i=1^M σ_i /|1+𝔪σ_i|^2=E + η/≤ E.
 
In addition, we have

    1/M∑_i=1^M γσ_i/(1+𝔪(z) σ_i)^2 =γ/M∑_i=1^M σ_i [ (1+𝔪(z) σ_i)^2- (1+𝔪(z) σ_i)]^2] /|1+𝔪(z) σ_i|^4≤ E -C .

Therefore, we have
          
    |z-1/M∑_i=1^M γσ_i/(1+𝔪(z) σ_i)^2| ≥   E- ( γ/M∑_i=1^M σ_i/(1+𝔪(z) σ_i)^2) ≥ C ≥ C√(κ+η).

Finally, taking the derivative of (<ref>), we obtain that
	
    (z )'=+z '=1/M∑_i=1^M γσ_i '/(1+ σ_i )^2.

	Hence, we finish the proof of (<ref>), which directly implies (<ref>).






00

 	huang  A. Adhikari, J. Huang. Dyson Brownian Motion for General β and Potential at the Edge, preprint, arXiv:1810.08308 (2018).

Ak N. I. Akhiezer. The classical moment problem: and some related questions in analysis, Hafner Publishing Co., New York (1965).

 	isotropic3  J. Alt, L. Erdős, T. Krüger and D. Schröder. Correlated Random Matrices: Band Rigidity and Edge Universality, Annals of Probability 48(2), 963-1001 (2020).
	
	baisample Z. D. Bai, J. W. Silverstein. CLT for linear spectral statistics of large-dimensional sample covariance matrices, Ann. Probab. 32(1A), pp 553-605 (2004).
	
	baibook Z. D. Bai, J. W. Silverstein. Spectral Analysis of Large Dimensional Random Matrices. Mathematics Monograph Series 2, Beijing: Science Press (2006).

	baisample2 Z. D. Bai, X. Wang and W. Zhou. Functional CLT for sample covariance matrices, Bernoulli 16(4), pp 1086-1113 (2010).
				
 	baiwigner Z. D. Bai, J. F. Yao. On the convergence of the spectral empirical process of Wigner matrices, Bernoulli 11(6), 1059-1092 (2005).	
		
	Bao+Pan+Zhou Z. Bao, G. Pan and W. Zhou. Universality for the Largest Eigenvalue of Sample Covariance Matrices with General Population, Ann. Statist. 43(1), 382-421 (2015).
 	
 	basor+widom E. L. Basor, H. Widom. Determinants of Airy Operators and Applications to Random Matrices, J. Stat. Phys.  96, 1-20 (1999).
		
 	
 	Bekerman+lodhia F. Bekerman, A. Lodhia. Mesoscopic central limit theorem for general β-ensembles, Ann. Inst. H. Poincare Probab. Statist.  54, 1917-1938 (2018).	
 	
	 quantum  L. Benigni. Eigenvectors distribution and quantum unique ergodicity for deformed Wigner matrices, preprint, arXiv:1711.07103 (2018).
		
 	free2 P. Biane. On the Free Convolution with a Semi-circular Distribution, Indiana Univ. Math. J. 46, 705-718 (1997).
	
	Bia98 P. Biane. Processes with free increments, Math. Z., 227(1), 143-174 (1998).
	
        Alex+Erdos+Knowles+Yau+Yin  A. Bloemendal, L. Erdős, A. Knowles, H.T. Yau and J. Yin. Isotropic local laws for sample covariance and generalized Wigner matrices, Electron. J. Probab. 19(33), 53pp (2014).	
        
        Alex+Knowles+Yau+Yin  A. Bloemendal,  A. Knowles, H.T. Yau and J. Yin. On the principal components of sample covariance matrices, Probab. Theo. and Relat. Field. 164(1-2), 459-551 (2016).       
	
 	Bourgadegap P. Bourgade. Extreme gaps between eigenvalues of Wigner matrices, preprint arXiv:1812.10376, (2018).
 	
 	BEYY P. Bourgade, L. Erdős, H.-T. Yau and J. Yin. Fixed energy universality for generalized Wigner matrices, Comm. Pure Appl. Math. 69(10), 1815-1881 (2016).
		
          logarithm P. Bourgade, K. Mody. Gaussian fluctuations of the determinant of Wigner matrices, Electron. J. Probab. 24(96), 28 pp (2019).
		
		meso1 A. Boutet de Monvel, A. Khorunzhy. Asymptotic distribution of smoothed eigenvalue density. I. Gaussian random matrices, Random Oper. and Stoch. Equ. 7(1), pp 1-22 (1999).
	
 	meso2 A. Boutet de Monvel, A. Khorunzhy. Asymptotic distribution of smoothed eigenvalue density. II. Wigner random matrices, Random Oper. and Stoch. Equ. 7(2), pp 149-168 (1999).

 	
 	Breuer+Duits J. Breuer, M. Duits. Universality of mesoscopic fluctuations for orthogonal polynomial ensembles, Comm. Math. Phys. 342 (2), 491-531 (2016).
 	
 	free M. Capitaine, C. Donati-Martin, D. Feral and M. Fevrier. Free Convolution with a Semicircular Distribution and Eigenvalues of Spiked Deformations of Wigner Matrices, Electron. J. Probab. 16(64), 1750-1792 (2011).

	
		Cusp2 G. Cipolloni, L. Erdős, T. Krüger and D. Schröder. Cusp universality for random matrices II: the real symmetric case, Pure and App. Analy. 1(4), 615-707 (2019).
		
			global2 S. Dallaporta, M. Fevrier. Fluctuations Of Linear Spectral Statistics Of Deformed Wigner Matrices. Preprint, arXiv:1903.11324 (2019).
 	

	 	Duits+Johansson  M. Duits, K. Johansson. On mesoscopic equilibrium for linear statistics in Dyson's Brownian Motion, Mem. Amer. Math. Soc. 255 (1222) (2018).
 	
 	erdos+knowles L. Erdős, A. Knowles. The Altshuler-Shklovskii formulas for random band matrices I: the
 	unimodular case, Comm. Math. Phys. 333(3), 1365-1416 (2015).
 	
 	erdos+knowles2 L. Erdős, A. Knowles. The Altshuler-Shklovskii formulas for random band matrices II: the
 	general case, Ann. H. Poincaré 16(3), 709-799 (2015).
 	
 	
EKY L. Erdős, A. Knowles and H.-T. Yau. Averaging fluctuations in resolvents of random band matrices, Ann. Henri Poincaré 14(8), 1837-1926 (2013). 
	
 	
 	isotropic2  L. Erdős, T. Krüger and D. Schröder. Random Matrices with Slow Correlation Decay, Forum of Mathematics, Sigma 7(8) (2019).
	
	 Cusp1 L. Erdős, T. Krüger and D. Schröder. Cusp universality for random matrices I: Local law and the complex Hermitian case. Preprint arXiv:1809.03971 (2018).
	 
	 ESY L. Erdős, B. Schlein and  H.-T. Yau.  Universality of random matrices and local relaxation flow, Invent. Math. 185(1), 75-119 (2011). 
	
	 SchnelliErdoes L. Erdős, K. Schnelli.  Universality for random matrix flows with time-dependent density, Ann. Inst. Henri Poincaré Probab. Stat. 53(4), 1606-1656 (2017). 
	
	 bull_american_math_society L. Erdős, H.-T. Yau.  Universality of local spectral statistics of random matrices, Bull. Amer. Math. Soc. 49(3), 377-414 (2012). 
		
 	book  L. Erdős, H.-T. Yau. A dynamical approach to random matrix theory, volume 28 of Courant Lecture Notes. American Mathematical Soc. (2017).
 	
	Hachem+Hardy+Najim_1 W. Hachem, A. Hardy, and J. Najim. Large complex correlated Wishart matrices: Fluctuations and asymptotic independence at the edges, Ann. Probab. 44(3), 2264-2348 (2016).
	
	Hachem+Hardy+Najim_2 W. Hachem, A. Hardy, and J. Najim. Large complex correlated Wishart matrices: the Pearcey kernel and expansion at the hard edge, Electron. J. Probab. 21(1), 36pp. (2016).
	
  	sparse Y. He. Bulk eigenvalue fluctuations of sparse random matrices, Preprint, arXiv:1904.07140 (2019).	
 	
 	HeKnowlesDensityCorrelation Y. He,  A. Knowles. Mesoscopic Eigenvalue Density Correlations of Wigner Matrices, Probability Theory and Related Fields 177, 147-216 (2020).

 	moment Y. He,  A. Knowles. Mesoscopic eigenvalue statistics of Wigner matrices, Ann. Appl. Probab. 27(3), 1510-1550 (2017).
 	 	
 Huang+Landon  J. Huang, B. Landon. Rigidity and a mesoscopic central limit theorem for Dyson Brownian
 motion for general beta and potentials, Prob. Theory and Related Fields 175(1-2), 209-253 (2019).
	
 	global H. C. Ji, J. O. Lee. Gaussian fluctuations for linear spectral statistics of deformed Wigner matrices, Random Matrices: Theory and Applications, 2050011 (2019).
 	
 	Johansson K. Johansson. On fluctuations of eigenvalues of random Hermitian matrices, Duke Math. J. 91(1), 151-204 (1998).
	
	 gumbel K. Johansson. From Gumbel to Tracy-Widom, Prob. Theory and Related Fields 138(1-2), 75-112 (2007).
	 
	Jonsson D. Jonsson. Some limit theorems for the eigenvalues of a sample covariance matrix, J. Multiva. Analy. 12(1), 1-38 (1982).	 
	
	edgecondition N. El Karoui. Tracy-Widom limit for the largest eigenvalue of a large class of complex
sample covariance matrices, Ann. Probab. 35(2), 663-714 (2007).	
		
	 isotropic  A. Knowles, J. Yin. Anisotropic local laws for random matrices, Prob. Theory and Related Fields 169(1-2), 257-352 (2017).
	 
	  character B. Landon, P. Sosoe. Applications of mesoscopic CLTs in random matrix theory, Preprint, arXiv:1811.05915 (2018).
 	
 	character2 B. Landon, P. Sosoe and H.-T. Yau. Fixed energy universality of Dyson Brownian motion, Advances in Math, 346, pp. 1137-1332 (2019).
	
 	LandonYau B. Landon, H.-T. Yau. Convergence of Local Statistics of Dyson Brownian Motion, Comm. Math. Phys., 355(3), pp. 949-1000 (2017).

 	LandonYau_edge B. Landon, H.-T. Yau. Edge statistics of Dyson Brownian motion, preprint, arXiv 1712.03881 (2017).

 	locallaw J. O. Lee, K. Schnelli. Local deformed semicircle law and complete delocalization for Wigner matrices with random potential, J. Math. Phys 54(10), 103504 (2013).
	
	edge J. Oon Lee, K. Schnelli. Edge universality for deformed Wigner matrices, Rev. Math. Phys 27(8), 1550018 (2015).	

 	eigenvector J. O. Lee, K. Schnelli. Extremal eigenvalues and eigenvectors of deformed Wigner matrices, Probab. Theory Related Fields 164(1), 165-241 (2016).
	
	sample_lee_kevin J. Oon Lee, K. Schnelli. Tracy-Widom distribution for the largest eigenvalue of real sample covariance matrices with general population, Ann. Appl. Probab. 26(6), 3786-3839 (2016).	
	
 	bulk J. O. Lee, K. Schnelli, B. Stetler and H.-T. Yau. Bulk universality for deformed Wigner matrices, Ann. Probab. 44(3), 2349-2425 (2016).
	
 	mesowigner A. Lodhia, N. J. Simm. Mesoscopic linear statistics of Wigner matrices, Preprint, arXiv:1503.03533 (2015).
 	
 	lytova+pastur A. Lytova, L. Pastur. Central limit theorem for linear eigenvalue statistics of random
 	matrices with independent entries, Ann. Prob. 37(5), 1778-1840 (2009).
	
	M+Pastur V.A. Marchenko, L. Pastur. Distribution of Eigenvalues for Some Sets of Random Matrices, Math. USSR Sb. 1(4), 457-483 (1967).		
 	
 	Min+Chen C. Min, Y. Chen. Linear Statistics of Random Matrix Ensembles at the
 	Spectrum Edge Associated with the Airy Kernel, Nuclear Physics B 950, 114836 (2020).
		
	secondfreeness J. A. Mingo, R. Speicher. Second order freeness and fluctuations of random matrices: I. Gaussian and Wishart matrices and cyclic Fock spaces, J. Func. Analy. 235(1), 226-270 (2006).
	
         Najim+Yao J. Najim, J. Yao. Gaussian fluctuations for linear spectral statistics of large random covariance matrices, Ann. Appl. Probab. 26(3), 1837-1887 (2016).
         
           O;rourke_Vu S. O'Rourke, V. Vu. Universality of local eigenvalue statistics in random matrices with external source, Random Matrices: Theory and Applications 3(2), 1450005 (2014).

 	solution L. Pastur. The spectrum of random matrices, Teor. Math. Phys 10, 64-74 (1972).
	
	Silverstein+Choi J. W. Silverstein, S. I. Choi. Analysis of the Limiting Spectral Distribution of Large Dimensional Random Matrices, J. Multiv. Anal. 54(2), 295-309 (1995).	
 	
 	M.Shcherbina M. Shcherbina. Central limit theorem for linear eigenvalue statistics of the Wigner and sample covariance random matrices, Zh. Mat. Fiz. Anal. Geom. 7(2), 176-192 (2011).
 	
 	edgesolution  T. Shcherbina. On universality of local bulk regime for the deformed Gaussian unitary ensemble, Zh. Mat. Fiz. Anal. Geom. 5(4), 396-433 (2009).
	
 	edgesolution2  T. Shcherbina. On universality of local edge regime for the deformed Gaussian unitary ensemble, J. Stat. Phys. 143, 455-481 (2011).

 	SosoeWongP. Sosoe, P. Wong. Regularity conditions in the CLT for linear eigenvalue statistics of Wigner matrices. Advances in Mathematics 249(20), 37-87 (2013)
 	
          Voi86 D. Voiculescu. Addition of certain non-commuting random variables, J. Funct. Anal. 66 (3), 323-346 (1986).

           freetimes2 D. Voiculescu. Multiplication of certain non-commuting random variables, J. Operator Theory 18(2), 223-235 (1987).
                      
           Voi93 D. Voiculescu. The analogues of entropy and of Fisher's information theory in free probability theory, I,  Comm. Math. Phys. 155, 71-92 (1993).

	freetimes D. Voiculescu, K. J.  Dykema and A. Nica. Free Random Variables: A Noncommutative Probability Approach to Free Products with Applications to Random Matrices, Operator Algebras and Harmonic Analysis on Free Groups, American Mathematical Society (1992).
	
 	wigner_annals  E. P. Wigner. Characteristic vectors of bordered matrices with infinite dimensions, Ann. of Math.  62 (3), 548-564 (1955).
	





