Gaussian Fluctuations in Complex Sample Covariance Matrices

Let X = (X i,j) m×n , m ≥ n, be a complex Gaussian random matrix with mean zero and variance 1 n , let S = X * X be a sample covariance matrix. In this paper we are mainly interested in the limiting behavior of eigenvalues when m n → γ ≥ 1 as n → ∞. Under certain conditions on k, we prove the central limit theorem holds true for the k-th largest eigenvalues λ (k) as k tends to infinity as n → ∞. The proof is largely based on the Costin-Lebowitz-Soshnikov argument and the asymptotic estimates for the expectation and variance of the number of eigenvalues in an interval. The standard technique for the RH problem is used to compute the exact formula and asymptotic properties for the mean density of eigenvalues. As a by-product, we obtain a convergence speed of the mean density of eigenvalues to the Marchenko-Pastur distribution density under the condition | m n − γ| = O(1 n)


Introduction and main results
Let X = (X i,j ) be a complex m×n, m ≥ n, random matrix, the entries of which are independent complex Gaussian random variables with mean zero and variance 1  n , namely Re(X i,j ), Im(X i,j ) form a family of independent real Gaussian random variables each with mean value 0 and variance 1 2n .Let S = X * X, then S can be viewed as a sample covariance matrix of m samples of n dimensional random vectors and it is of fundamental importance in multivariate statistical analysis.
The complex sample covariance matrices was first studied by Goodman (5) and Khatri (9).The distribution dµ(S) of S is given by The measure dµ on the matrices produces naturally a measure on the corresponding n real eigenvalues λ i .It turns out that this induced measure can be explicitly calculated and its joint density is given by where Z = Z n,m > 0 is again a normalization constant depending on m, n.
It is the explicit form (1.1) (sometimes called the Laguerre unitary ensemble) that makes possible the deep and thorough asymptotic analysis, both inside the bulk and at the edge of the spectrum.There are actually many other well-known ensembles with such a determinantal point process representation in random matrices and random growth models.See (10) for recent works.
Let V n,m (x) be such that e −nVn,m(x) = x m−n e −nx , x ∈ R + .Let p j (x), j ≥ 0, be a sequence of orthonormalized polynomials with γ j the highest coefficient with respect to the weight function e −nVn,m(x) .That is, R + p j (x)p k (x)e −nVn,m(x) dx = δ jk , j, k ≥ 0.
Define the kernel K n by Moreover, for 1 ≤ k ≤ n, the k-marginal dimensional density is given by In particular, ρ n,1 (x 1 ) describes the overall density of the eigenvalues.
The kernel K n (x, y) can also be represented by the so-called Christoffel-Darboux identity.For x = y it holds and for x = y one has m(x) . (1.5) The classic Marchenko-Pastur (11) theorem states that as n → ∞ such that m n → γ ≥ 1, almost surely in distribution, where µ γ is the probability density function of the M-P distribution with parameter γ, namely and α = ( √ γ − 1) 2 and β = ( √ γ + 1) 2 .
A recent remarkable work is on a limiting distribution of the largest eigenvalue.Let λ (1) > λ (2) > • • • > λ (n) be the ordered list of eigenvalues.While studying a random growth model of interest in probability, Johansson (8) derived the F 2 limit distribution for λ (1) .Specifically speaking, define in distribution, where F 2 is the Tracy-Widom distribution discovered in the Gaussian unitary ensemble (GUE).Analogs for λ (k) with k fixed has also been studied.
In this paper we deal with the distribution of λ (k) as n and k tend to infinity.Let (1.9) Our main results are as follows.
Theorem 1. (the bulk case).Set where k = k(n) is such that k n → a ∈ (0, 1) as n → ∞.Then as in distribution.
Remarks: (1) Analogous results have recently been established for the eigenvalues of the GUE by (7).Our proofs are again based on the Costin-Lebowitz-Soshnikov theorem.One should be able to apply the same methodology to many other ensembles.A work on the discrete Krawtchouk ensemble and Hahn ensemble is in progress.
(2) It is remarkable that with regard to the Plancherel measure on the set of partitions λ of n, the rows λ 1 , λ 2 , λ 3 , • • • of λ behaves, suitably scaled, like the 1st, 2nd, 3rd and so on eigenvalues of a random matrix from the GUE.On the other hand, the boundary of Young shape λ, a polygonal line, behaves like where Ω(x) is closely connected to Wigner's semicircle law and U (x) is a kind of generalized Gaussian random process.This is often referred to as Kerov's central limit theorem.
Borodin, Okounkov and Olshanski (1) established an exact determinantal formula for the correlation function of the Poissonized Plancherel measures.It would be expected that one could provide a proof of (1.12) using the Costin-Lebowitz-Soshnikov argument.
(3) It is well-known that µ γ (x) is also the asymptotic distribution for zeros of orthogonal polynomial p n (x).Therefore it is natural to ask whether the k-th eigenvalue λ (k) fluctuates around the corresponding zero of p n (x).The above theorems partly confirm the conjecture although a rigorous argument is still open.
We say as usual that λ (k) obeys the central limit theorem if there exist a n > 0 and b n such that Hence it suffices to prove where Φ(x) is the standard normal distribution function.
This in turn follows from the Costin-Lebowitz-Soshnikov theorem (2; 13) as long as (1.15) (1.15) will be used to determine the normalizing constants.
This paper is organized as follows.In Section 2 we shall give an exact formula and asymptotics for the mean density of eigenvalues λ 1 , • • • , λ n uniformly valid on the entire real axis based on the standard technique of asymptotic analysis of RH problems.Analogues for the GUE type ensemble have been studied and used to give a proof of a complete large N expansion for the partition function by Ercolani and McLaughlin (4).Just recently did Vanlessen (14) consider strong asymptotics of Laguerre type orthogonal polynomials on R + with respect to the weight function x α e −Q(x) where α > −1 and Q denotes a polynomial with positive leading coefficient.As a direct consequence, we obtain the convergence speed of ρ n,1 to µ γ when | m n − γ| = O( 1 n ).In the case when m = [γn], the speed is due to Götze and Tikhomirov (6) using recursive equations.
The proofs of Theorems 1 and 2 will be given in Section 3. Starting from (1.13)-(1.15),we need only compute the expectation and variance for the number of eigenvalues in an interval.The computation heavily depends on the exact formula and asymptotics for ρ n,1 (x) and K n (x, y) in the entire real axis.
Throughout the paper, there are lots of positive numerical constants, whose values are of no importance.We shall use c, C for simplicity, which may take different values in different places.

The mean density of eigenvalues: exact formula and asymptotics
This section is devoted to the asymptotic analysis of orthogonal polynomials p n (x) and p n−1 (x) with respect to the n dependent weight function e −nVn,m(x) , in which techniques are used for the asymptotic analysis of RH problems, first developed for singularity by Deift-Zhou.Deift (3) is a standard excellent reference in this area.For completeness of notations we summarize major steps of rigorous analysis, including a series of equivalent RH problems and explicit transformations below.
Let us start with the following RH problem: We remark that a similar analysis was given for x α e −nx with α > −1 by (14).
The unique solution of the RH problem for U is given by where C(•) denotes the Cauchy operator.In particular, We have by (1.4) and (1.5) and 3) The equilibrium measure plays an important role in the asymptotic analysis of RH problems.It turns out that the measure with density function µ n,m (x) is the equilibrium measure in the presence of the external field V n,m (x).It is actually easy to check that µ n,m satisfies the Euler-Lagrange variational conditions: there exists a real number l n,m such that and In order to normalize the RH problem for U at infinity, we use the log-transform of the equilibrium measure.Define where we take the principal branch of the logarithm, so that We now list several important properties of the function g n (z): (2.4) and g n,+ (x) − g n,− (x) = 0, for x ∈ (β n,m ∞).
Now we are ready to perform the transformation U → T .Define the matrix valued function T as where Note that, the function e ngn(z) has no jump across (−∞, α n,m ), so that T has an analytic continuation to C \ [α n,m , ∞).It is then straightforward to check, using the conditions of the RH problem for U , that T is the unique solution of the following equivalent RH problem: e n(g n,+ (x)−g n,− (x)) , for x ∈ (α n,m , β n,m ), ), as z → ∞.From (2.5) we see that the diagonal entries of ν T (x) on (α n,m , β n,m ) are rapidly oscillating for large n.From (2.4) the jump matrix ν T on (0, ∞) \ (α n,m , β n,m ) converges exponentially fast to the identity matrix as n → ∞.Now we will transform the oscillatory diagonal entries of the jump matrix ν T (x) on (α n,m , β n,m ) into exponentially decaying off-diagonal entries.This lies the heart of the Deift-Zhou steepest descent method.
In order to perform the transformation T → S we will introduce scalar functions ψ n and ξ n .Define with principal branches of powers.So, the + boundary value of ψ n (z) on [α n,m , β n,m ] is precisely µ n,m (z).In particular, we have where the path of integration does not cross the real axis.
The important feature of the function ξ n is that ξ n,+ and ξ n,− are purely imaginary on (α n,m , β n,m ) and satisfy for x ∈ (α n,m , β n,m ) So, 2ξ n,+ (z) and 2ξ n,− (z) provide analytic extensions for g n,+ (x) − g n,− (x) into the upper half plane and lower half plane, respectively.
Further one can prove the existence of a δ 1 > 0 such that (2.9) The jump matrix ν T for T can be written in terms of the scalar function ξ n as (2.10) A simple calculation, using the fact that ξ n,+ (x) + ξ n,− (x) = 0 for x ∈ (α n,m , β n,m ), then shows that ν T has on the interval (α n,m , β n,m ) the following factorization, Now we are ready to do the transformation T → S. Let Σ S = ∪ 4 j=1 Σ j be the oriented lens shaped contours as shown in Figure 2.1.The precise form of the lens (in fact of the lips Σ 1 and Σ 3 ) is not yet defined but for now we assume that it will contained in the region where (2.9) holds.
for z outside the lens, , for z in the upper part of the lens, , for z in the lower part of the lens, (2.12) with the upper part of the lens we mean the region between Σ 1 and Σ 2 , and with the lower part of the lens the region between Σ 2 and Σ 3 .
One can easily check, using (2.10), (2.11) and the conditions of the RH problem for T , that S satisfies the following RH problem: Note that the jump matrix ν S on Σ 1 , Σ 3 , and Σ 4 converges exponentially fast to the identity matrix.We shall see that the leading order asymptotics of U will be determined by a solution P ∞ , which is often referred to as the parametrix for the outside region, of the following RH problem: with

.15)
Before we can do the final transformation S → R we need to do a local analysis near α n,m and β n,m since the jump matrices for S and P ∞ are not uniformly close to each other in the neighborhood of these points.We will only construct the parametrix P n near the right endpoint β n,m below.As for the left endpoint, if m n → γ > 1, then α > 0 and a similar construction is valid in the left endpoint α n,m .While in the case m n → 1, α = 0 becomes a hard edge since all eigenvalues of covariance matrices are positive.The behavior of polynomials near the origin is described via Bessel function.We refer to (14) for modifications.
as n → ∞ uniformly for z on the boundary ∂U βn,m,δ of the disk U βn,m,δ , and for δ in a compact subsets of (0, δ 2 ).
The construction of P n is based on an auxiliary RH problem for Ψ in the ζ-plane with jumps on the following oriented contour γ σ consisting of four straight rays (see Figure 2.2 below) . These four rays divide the complex plane into four regions.The RH problem for Ψ: for ζ ∈ γ σ with ν 1 as shown in Figure 2.2, i.e., Ψ has the following asymptotic behavior at infinity, Here It is well-known that Ψ is defined by (2.16) The idea is now to construct Ψ(f n (z)) for appropriate biholomorphic maps f n : U βn,m,δ 2 → f n (U βn,m,δ 2 ) with f n (β n,m ) = 0. We will choose these biholomorphic maps to compensate for the factor e −nξn(z)σ 3 . Define with φ n (z) defined in the following proposition.
Proposition.There exists a δ 2 > 0 such that there are biholomorphic maps (1) There is a constant c 0 > 0 such that for all z ∈ U βn,m,δ 2 and all n ≥ 1 the derivative of φ n can be estimated by (2) (2.18) Note that the function φn (z) has no jump across (α n,m , β n,m ) so that φn (z) has analytic continuation to C \ ((−∞, α n,m ) ∪ {β n,m }).From (2.7 ) and (2.8 ) it follows that Using Cauchy's theorem there exists an n 2 such that for all n ≥ n 2 and |s − Inserting this into (2.19)we obtain (2.20) Therefore the isolated singularity of φn (z) at β n,m is removable so that φn (z) is analytic in C \ (−∞, α n,m ], and there exists δ > 0 such that Re φn (z) > 0 for all |z − β n,m | ≤ 1/4 and n ≥ n 2 .This yields Observe that φ n (z) is uniformly (in n and z) bounded in U βn,m,δ , i.e., sup This implies, by using Cauchy's theorem for derivatives, that φ ′′ n (z) is also uniformly (in n and z) bounded in U βn,m,δ for a smaller δ, i.e., sup , so that Therefore, there exists 0 < δ 2 < δ such that for all n ≥ n 2 the φ n (z) are injective and hence biholomorphic in U βn,m,δ .We conclude (1).Now (2) and (3) easily follows from (1), so the proof is complete.
For later reference, observe that the biholomorphic maps is given by where fn (z) = φn (z) 2/3 with φn given by (2.18).

.35)
Taking derivatives at both sides of (2.35) yields 1 0 e −2ξn(x) 1 e n(gn(x)− (2.37) Also, P ∞ (x) can be written as so we have where and in the sequel .
(3).We need only to give a proof in the case of x ∈ (β n,m − δ, β n,m + δ) since the other case is similar.
then it holds true . (2.54) The proof of Theorem 3 is now complete.
Recall that the Airy function is bounded on the real line, and for r > 0 (see ( 7),( 13)) (2.60) and

A simple calculation shows
In this section we shall give the proof of Theorems 1 and 2. By virtue of (1.13), (1.14) and (1.15), we shall basically focus on the computation of expectation and variance, which are given in the following lemmas.
Proof.For some δ > 0, we have By virtue of (2.32), for x ≥ β n,m + δ, ρ n,1 (x) is exponentially small in n and exponentially decaying in x.So it holds For Here and We now look at the different terms in the asymptotic expression for nρ n,1 (x) above.Observe that (2.60) and (2.61) give the asymptotic expansion of Airy function for large negative values.
We need also the following asymptotic expansion of Airy function for large positive values (see (12), ( 13)): for r > 0 and Hence it follows for any x ∈ (−∞, ∞) and Thus we have Also, integrating the third term of (3.4)only gives a contribution of order n −1 .The main contribution comes from the second term.In fact, a primitive function can be found for this expression: where c > 0.
One can now use (2.60) and (2.61) to get the desired result.Indeed, |f n (t n )| ≥ C under the hypothesis, so we have and Thus it follows where in the last equation we use the condition β n,m − t n → 0. We now conclude the proof.
Lemma 2. Assume t is defined as in Theorem 1.Let Proof.We use the fact again The above integral (3.16) can be written as where δ > 0 is a constant.
We now look at the integrals at the right hand side of (3.17).First, one easily sees from (2.56) that In turn, by a primitive function argument used in (3.11), we need to prove Now combining (3.17 Observe that by (2.1) (3.25) We will compute the above determinant in two basic cases separately.The first case is when t ≤ β n,m − δ for some δ > 0, i.e., in the bulk.The second case is when β n,m − t → 0 as n → ∞, i.e., near the right spectrum edge.
Proof.First consider the sub-domain where both variables are in the bulk: For z ∈ (α n,m + δ, β n,m − δ), define θ(z) to be such that Then a simple computation shows So it follows from (2.14) that .
(3.30) Define (3.32) Thus using (2.29) and (3.32) we have for (x, y) ∈ Γ, then it follows from (2.7), (2.8) and (3.31) which in turn implies A similar argument to (3.41) gives In order to calculate the integral of K n (x, y) 2 over Γ 1 1 , we shall use an elementary fact

.45)
It immediately follows from (3.45) that Also, one easily sees (3.47) Both the integrals are easy to estimate: and both F ′ (y) and F ′′ (y) are bounded in (α n,m + δ, β n,m − δ).We easily have (3.49) Combining (3.48) and (3.49) yields Also, it is easy to see then using (3.45) again, we have Similarly, it follows (3.52) Finally, we need to look at the sub-domain Γ 3 =: I n × I c n − Γ.We divide Γ 3 into the following four sub-domians: 3 but there are no difficulties.One can just take the absolute values in the integral since |x − y| ≥ β n,m − α n,m and the result Proof.The argument is similar to that of Lemma 3 but requires a finer partition.Now we divide I n × I c n into the following nine sub-domains: where δ > 0, C > 0, r(n) and ε are defined by In the corner Ω 0 we shall use the representation (1.2).By use of the Cauchy-Schwarz inequality we have K 2 n (x, y) ≤ K n (x, x)K n (y, y).Having separated the variables one can now use the calculations of the expected value giving Now let us look at the sub-domain Ω 1 .We shall prove the following estimate for K n (x, y): To see this, we use for z It is now sufficient to calculate the integrals of K n (x, y) 2 in Ω i , i = 3, 4, 5 are O(log(r(n))).The calculation is straightforward and almost same as in the case of GUE (see (7)), so some details will be skipped.
In the sub-domain Ω 6 one can perform the same calculations as in In Ω 7 one can use the fact x − y ≥ δ to show that the contribution form this domain is O(1).In Ω 8 one can easily get K n (x, y) is exponentially small in n and exponentially decaying in x (or y).Thus the contribution from this domain is o(1).The proof of Lemma 4 is complete.
dµ(S) = 1 Z (det S) m−n e −trS dS, for S ∈ M n (C) + , where Z = Z n,m > 0 is a normalization constant depending on m, n and dS = ( n j=1 dS jj ) j<k d(ReS jk )d(ImS jk ).

Figure 1 :
Figure 1: Σ S = ∪ 4 j=1 Σ j Define an analytic matrix valued function S on C \ Σ S as
then by (2.60) and (2.61) again we have