Level-Spacing Distributions and the Airy Kernel

이 논문은 random matrix models에서 edge of spectrum에 해당하는 level-spacing distribution functions을 분석하는 것을 목표로 한다. 이러한 분포는 Airy kernel K(x,y) = Ai(x)Ai'(y) - Ai'(x)Ai(y)/(x-y)에 의해 정의된다. 이 kernel은 sine kernel sin π(x-y)/π(x-y)의 scaling limit과 유사성이 있다.

이 논문에서는 다음 결과를 도출한다:

1. Airy kernel K(x,y)가 가진 완전한 적분 가능 시스템(P.D.E.)을 찾는다.
2. semi-infinite interval (s,∞)에서 Fredholm determinant D(J;λ)를 Painlevé transcendent과 연결시킨다.
3. level-spacing distribution function E(n;s)의 아сим포티카 확장식을 도출한다.
4. probability density for the nth largest eigenvalue F(n;s)도 도출한다.

이 논문은 random matrix models의 edge of spectrum에 해당하는 level-spacing distribution functions을 분석하여 새로운 결과를 도출하였다.

Level-Spacing Distributions and the Airy Kernel

arXiv:hep-th/9211141v1 1 Dec 1992ITD 92/93–9Level-Spacing Distributions and the Airy KernelCraig A. Tracy∗Department of Mathematics and Institute of Theoretical Dynamics,University of California, Davis, CA 95616, USAHarold Widom†Department of Mathematics,University of California, Santa Cruz, CA 95064, USAScaling level-spacing distribution functions in the “bulk of the spectrum” inrandom matrix models of N × N hermitian matrices and then going to the limitN →∞, leads to the Fredholm determinant of the sine kernel sin π(x −y)/π(x −y).Similarly a scaling limit at the “edge of the spectrum” leads to the Airy kernel[Ai(x)Ai′(y) −Ai′(x)Ai(y)] /(x −y). In this paper we derive analogues for this Airykernel of the following properties of the sine kernel: the completely integrable systemof P.D.E.’s found by Jimbo, Miwa, Mˆori and Sato; the expression, in the case of asingle interval, of the Fredholm determinant in terms of a Painlev´e transcendent; theexistence of a commuting differential operator; and the fact that this operator canbe used in the derivation of asymptotics, for general n, of the probability that aninterval contains precisely n eigenvalues.∗e-mail address: catracy@ucdavis.edu†e-mail address: widom@cats.ucsc.edu1

I. INTRODUCTION AND SUMMARY OF RESULTSA.

IntroductionIn this paper we present new results for the level spacing distribution functions obtainedfrom scaling random matrix models of N ×N hermitian matrices at the edge of the supportof the (tree-level) eigenvalue densities when the parameters of the potential V are not “finelytuned.” This universality class is already present in the Gaussian Unitary Ensemble. It isknown [3,10,20,23] that these distribution functions are expressible in terms of a Fredholmdeterminant of an integral operator K whose kernel involves Airy functions.It turns out that there are striking analogies between the properties of this Airy kernelK(x, y) = A(x)A′(y) −A′(x)A(y)x −ywhere A(x) =√λ Ai(x), and the sine kernelλπsin π(x −y)x −y,whose associated Fredholm determinant describes the classical level spacing distributionfunctions first studied by Wigner, Dyson, Mehta, and others [7,20,26].To be a little more explicit about how these kernels arise we remind the reader that inthe GUE the probability density that the eigenvalues of a random N × N hermitian matrixlie in infinitesimal intervals about x1, · · ·, xN is given by [20]PN(x1, · · ·, xN) = 1N!det (KN(xi, xj))i,j=1,···,NwhereKN(x, y) =N−1Xk=0ϕk(x)ϕk(y) .Here {ϕk(x)} is the sequence obtained by orthonormalizing the sequencenxk exp(−x2)o2

over (−∞, ∞). More generally the n-point correlation function Rn(x1, · · ·, xn), which isthe probability density that n of the eigenvalues (irrespective of order) lie in infinitesimalintervals about x1, · · ·, xn, is given byRn(x1, · · · , xn) = det (KN(xi, xj))i,j=1,···,nIn particular R1(x), the density of eigenvalues at x, equals KN(x, x).Replacing KN(x, y) by1R1(z)KN z +xR1(z), z +yR1(z)!has the effect of making the point z in the spectrum the new origin and rescaling so thatthe eigenvalue density at this point becomes equal to 1.

Now if z is fixed thenR1(z) = KN(z, z) ∼1π√2Nas N →∞, andlimN→∞π√2N KN z +πx√2N , z +πy√2N!= 1πsin π(x −y)(x −y),so we obtain the sine kernel in this way. (This limit is obtained by applying the Christoffel-Darboux formula to KN(x, y) and then using the known asymptotics of ϕN(x) as N →∞.

)The “edge of the spectrum” corresponds to z ∼±√2N and one has there the scaling limit[10,23]limN→∞121/2N1/6 KN√2N +x21/2N1/6 ,√2N +y21/2N1/6= Ai(x)Ai′(y) −Ai′(x)Ai(y)x −y.We now present some of the analogies which we have found. (That such analogiesbetween properties of the sine kernel and Airy kernel exist is perhaps not suprising sincethey are both scaling limits of the same family of kernels.

)The first is the analogue ofthe completely integrable system of P.D.E.’s of Jimbo, Miwa, Mˆori, and Sato [17] when theunderlying domain is a union of intervals. The second is the fact that in the case of thesemi-infinite interval (s, ∞) (the analogue of a single finite interval for the sine kernel) the3

Fredholm determinant is closely related to a Painlev´e transcendent of the second kind (thefifth transcendent arises for the sine kernel [17]). And the third is the existence of a secondorder differential operator commuting with the Airy operator K. (The existence of such adifferential operator in the sine kernel case has been known for some time, see e.g.

[14], page201.) This last fact leads to an explicit asymptotic formula, as the interval (s, ∞) expands,for the probability that it contains precisely n eigenvalues (n = 1, 2, .

. .) (the analogue ofresults in [2,31]).

Here are our results, stated in some detail. (A short announcement appearsin [29].)B.

The System of Partial Differential EquationsWe setJ =m[j=1(a2j−1, a2j)and write D(J; λ) for the Fredholm determinant of K acting on J. We think of this as afunction of a = (a1, .

. .

, a2m). Thenda log D(J; λ) = −2mXj=1(−1)jR(aj, aj) daj(1.1)where R(x, y) is the kernel of the operator K(I−K)−1 and da denotes exterior differentiationwith respect to the aj’s.

We introduce quantitiesqj = (I −K)−1A(aj) ,pj = (I −K)−1A′(aj) ,(1.2)(which are the analogue of the quantities r±j of [17]; see also [9,15,21,22]) as well as twofurther quantitiesu =A, (I −K)−1A,v =A, (I −K)−1A′(1.3)where the inner products refer to the domain J. Then the equations read4

∂qj∂ak= (−1)k qjpk −pjqkaj −akqk(j ̸= k),(1.4)∂pj∂ak= (−1)k qjpk −pjqkaj −akpk(j ̸= k),(1.5)∂qj∂aj= −Xk̸=j(−1)k qjpk −pjqkaj −akqk + pj −qju ,(1.6)∂pj∂aj= −Xk̸=j(−1)k qjpk −pjqkaj −akpk + ajqj + pju −2qjv ,(1.7)∂u∂aj= (−1)jq2j ,(1.8)∂v∂aj= (−1)jpjqj . (1.9)Moreover the quantities R(aj, aj) appearing in (1.1) are given byR(aj, aj) =Xk̸=j(−1)k (qjpk −pjqk)2aj −ak+ p2j −ajq2j −2pjqju + 2q2jv .

(1.10)These equations are derived very much in the spirit of [28]; see also [9,15,21].C. The Ordinary Differential EquationsFor the special case J = (s, ∞) the above equations can be used to show that q(s; λ)(the quantity q of the last section corresponding to the end-point s) satisfiesq′′ = s q + 2q3 ,(′= dds)(1.11)withq(s; λ) ∼√λ Ai(s)as s →∞.

(1.12)This equation is a special case of the PII differential equation [1,16,18,25]. One can similarlyderive for R(s) := R(s, s), which in view of (1.1) equalsdds log D(J; λ) ,the third-order equation12 R′′R′!′−RR′ + R′ = 0 .

(1.13)5

It is also the case thatR′(s) = −q(s; λ)2(1.14)and this gives the following simple formula for D(J; λ) in terms of a PII transcendent:D(J; λ) = exp−Z ∞s (x −s)q(x; λ)2 dx. (1.15)This is much simpler than the corresponding representation of D(J; λ) for the sine kernelin terms of a PV transcendent.

The fact that q(s; λ) satisfies (1.11) can also be obtained bycombining some results in [1,5,13]. Thus in this case of a semi-infinite interval our resultsconnect with inverse scattering.D.

Asymptotics and the Commuting Differential OperatorAgain we take J = (s, ∞) and consider asymptotics as s →−∞. (Asymptotics ass →∞can be obtained trivially from the Neumann series for (I −K)−1.) From the randommatrix point of view the interesting quantities areE(n; s) := (−1)nn!∂n∂λn D(I; λ)λ=1 .This is the probability that exactly n eigenvalues lie in J.

Using our differential equationsplus the known asymptotics of q(s; 1) as s →−∞[13] we obtain the asymptotic expansionsR(s; 1) ∼14s2 −18 s +964 s4 −189128 s7 + 21663512 s10 + · · ·(s →−∞). (1.16)(where R(s; 1) denotes R(s) corresponding to λ = 1) andE(0; s) = D(J; 1) ∼τ0(−s)1/8 exps3/12×1 −326 s3 + 2025213 s6 −2470825219 s9+ · · ·(s →−∞)(1.17)where τ0 is an undetermined constant.The analogue of this formula for the sine kernel was obtained by Dyson [8].

The ana-logue of our constant τ0 in the expansion was found to be equal to 21/12e3ζ′(−1) (ζ(s) is the6

Riemann zeta function). This constant was obtained by scaling a result of the second author[30] on Toeplitz matrices.

Such a result is not available to us now. Using (1.16) and the factthat R(s) is the logarithmic derivative of E(0; s), it is straightforward to derive an integralrepresentation for log τ0.

This integral can be evaluated approximately by numerically inte-grating (1.11) and using (1.14). This leads to the value log τ0 ≃−0.136540.

We conjecturethat in factτ0 = eζ′(−1) 2124 = 0.87237 14149 54127 . .

. .This agrees with the numerically computed value, up to its degree of accuracy.For asymptotics of E(n; s) for general n we introducer(n; s) := E(n; s)E(0; s) .Successive differentiation of (1.11) with respect to λ, plus the known asymptotics of q(s; 1),allows us to find asymptotic expansions for the quantitiesqn(s) := ∂nq∂λnλ=1(1.18)(for the analogue in the sine kernel case see [2]); and these in turn can be used to findexpansions for r(n; s).

One drawback of this approach is that yet another undeterminedconstant factor enters the picture. (In [2,31] Toeplitz and Wiener-Hopf techniques, notavailable for the Airy kernel, fixed this constant.

)Another drawback is that when oneexpresses the r(n; s) in terms of the qn(s) a large amount of cancellation takes place, withthe result that even the first-order asymptotics of r(n; s) are out of reach by this methodwhen n is large.There is, however, another approach (briefly indicated in [2], and with details in [28], forthe sine kernel case). We haver(n; s) =Xi1<··· λ1 > · · · are the eigenvalues of the integral operator K (with λ = 1).

Now justas the operator with the sine kernel commutes with the differential operator for the prolate7

spheroidal wave functions [14], so does the Airy operator commutes with the differentialoperator L given byLf(x) = ((x −s) f ′(x))′ −x(x −s)f(x) . (1.20)An application of the WKB method, plus a trick allows us to derive the following asymptoticformula for λi with i fixed:1 −λi ∼√πi!

25i+3 t3i/2+3/4 exp−83t3/2(s = −2t →−∞). (1.21)(The analogue of this for the sine kernel is in [11].) From this and (1.19) we deducer(n, −2t) ∼1!2!

· · ·(n −1)!πn/2 2(5n2+n)/2 t−3n2/4 exp83nt3/2. (1.22)(Note that this can be used to fix the constant mentioned in the last paragraph.)E.

Probability Density for the nth Largest EigenvalueAn elementary probability argument shows that if F(n; s) (n = 0, 1, · · ·) denotes theprobability density for the (scaled) position of the nth largest eigenvalue (n = 0 is thelargest), thenF(0; s) = dE(0; s)dsand for n ≥1F(n; s) −F(n −1, s) = dE(n; s)ds.From the preceding expressions for E(n; s) in terms of q(s; 1) and qn(s) we can deriveformulas for F(n; s) which can in turn be numerically evaluated to produce tables of valuesfor these probability densities. We have done this for n = 0 and n = 1 and in Fig.

1 we plotE(0; s) and E(1; s) and in Fig. 2 we plot F(0; s) and F(1; s).8

II. THE SYSTEM OF PARTIAL DIFFERENTIAL EQUATIONSA.

Derivation of the EquationsA central role in the derivation of our equations is played by the commutator identityhL, (I −K)−1i= (I −K)−1 [L, K] (I −K)−1(2.1)for any operators K and L. Taking L to be D (= d/dx) or M (=multiplication by theindependent variable) give two of the three basic ingredients from which the equations arederived. The third is the equally simple fact that if K depends upon a parameter α thenddα (I −K)−1 = −(I −K)−1 dKdα (I −K)−1 .

(2.2)It will be very useful for us to think of K as acting not on J but on (−∞, ∞) and tohave kernelK(x, y)χJ(y)(2.3)where χJ is the characteristic function of J. We continue to denote the resolvent kernel ofK by R(x, y) but note that while it is smooth in x it is discontinuous at y = aj.

Thus thequantity R(aj, aj) appearing in (1.1) must be interpreted to meanlimy→ajy∈JR(aj, y);similarly for pj and qj in (1.2) and like quantities appearing later. The definitions of u andv must be modified to readu =AχJ, (I −K)−1 A,v =AχJ, (I −K)−1 A′.

(2.4)Notice that since(I −K)−1 A = (I −K)−1 AχJin Jthis agrees with the original definitions of u and v given by (1.3).9

We can think of K as an operator taking smooth functions to smooth functions andso the operators KD and DK make sense, the former having distributional kernel. Thetranspose operator Kt takes distributions to distributions.With these preliminaries out of the way we begin with the derivation of the representa-tions of the quantities R(aj, ak) in terms of the qj, pj, u and v. We introduce the notationQ = (I −K)−1 A,P = (I −K)−1 A′ ,(2.5)so that qj = Q(aj), pj = P(aj).

(For the moment we ignore the dependence of Q and P onthe parameter a = (a1, · · · , a2m); this will come later. )The first commutator relation we use is the simple[M, K] .= (A(x)A′(y) −A′(x)A(y)) χJ(y)(where .= means “has kernel equal to”).

Hence from (2.1) with L = MhM, (I −K)−1i .= Q(x)I −Kt−1 A′χJ(y) −P(x) (I −K)−1 AχJ(y)(2.6)(The transpose here arises from the general fact that if L .= U(x)V (y), then T1LT2.=T1U(x) T t2V (y).) If we use the fact thatI −Kt−1 AχJ = (I −K)−1 Aon J,(2.7)we deduce from (2.5) and (2.6) that [15]R(x, y) = Q(x)P(y) −P(x)Q(y)x −y(x, y ∈J, x ̸= y) .

(2.8)In particular, therefore,R(aj, ak) = qjpk −pjqkaj −ak(j ̸= k),(2.9)R(aj, aj) = Q′(aj)pj −P ′(aj)qj ′ = ddx!. (2.10)To compute Q′(x) and P ′(x) we consider first the commutator [D, K].

If an operator Lhas distributional kernel L(x, y) then10

[D, L] .= ∂∂x + ∂∂y!L(x, y) . (2.11)It is an easy verification, using the equation A′′(x) = xA(x), that ∂∂x + ∂∂y!K(x, y) = −A(x)A(y)(2.12)and so (recall the definition (2.3) of the kernel of K)[D, K] .= −A(x)A(y)χJ(y) −Xk(−1)kK(x, ak)δ(y −ak) .Hence, applying (2.1) with L = D (and recalling definitions (2.5)),hD, (I −K)−1i= −Q(x)I −KtAχJ(y) −Xk(−1)kR(x, ak)ρ(ak, y)(2.13)whereρ(x, y) = δ(x −y) + R(x, y)is the distributional kernel of (I −K)−1.We proceed with the computations of Q′(x) and P ′(x).

Differentiating the first relationin (2.5) and using (2.13) giveQ′ = D(I −K)−1A = P +hD, (I −K)−1iA= P −QI −Kt−1 AχJ, A−Xk(−1)kR(·, ak)qk .Thus (recall (2.4))Q′(x) = P(x) −Q(x)u −Xk(−1)kR(x, ak)qk . (2.14)SimilarlyP ′ = (I −K)−1A′′ +hD, (I −K)−1iA′ .Since A′′(x) = xA(x) this is equal toM (I −K)−1 A −hM, (I −K)−1iA +hD, (I −K)−1iA′ .11

We now use formula (2.6) and (2.13) to deduce thatP ′(x) = xQ(x) −2Q(x)v + P(x)u −Xk(−1)kR(x, ak)pk . (2.15)Setting x = aj in (2.14) and (2.15) and using (2.10) giveR(aj, aj) = p2j −ajq2j + 2q2jv −2pjqju +Xk(−1)kR(aj, ak)(qjpk −pjqk) .In view of (2.9) this is precisely equation (1.10).Most of the work is already done and we can derive the equations (1.4)–(1.9) very quickly.First, we have the easy fact∂∂akK .= (−1)kK(x, ak)δ(y −ak)and so by (2.2)∂∂ak(I −K)−1 .= (−1)kR(x, ak)ρ(y −ak) .

(2.16)At this point we must keep in mind that Q and P are functions of a as well, and so wedenote them now by Q(x, a) and P(x, a), respectively. We deduce immediately from (2.16)and the definitions (2.5) that∂∂akQ(x, a) = (−1)kR(x, ak)qk,∂∂akP(x, a) = (−1)kR(x, ak)pk .

(2.17)Since qj = Q(aj, a) and pj = P(aj, a) this gives∂qj∂ak= (−1)kR(aj, ak)qk,∂pj∂ak= (−1)kR(aj, ak)pk(j ̸= k).In view of (2.9) these are equations (1.4) and (1.5). Moreover∂qj∂aj= ∂∂x + ∂∂aj!Q(x, a)x=aj,∂pj∂aj= ∂∂x + ∂∂aj!P(x, a)x=ajand (2.17), (2.14) and (2.15) give∂qj∂aj= pj −qju −Xk̸=j(−1)kR(aj, ak)qk ,∂pj∂aj= ajqj −2qjv + pju −Xk̸=j(−1)kR(aj, ak)pk .12

In view of (2.9) again, these are equations (1.6) and (1.7).Finally, using the definition of u in (2.4), the fact∂∂ajχJ(y) = (−1)jδ(y −aj) ,and (2.16) we find that∂u∂aj= (−1)jA(aj)qj + (−1)j AχJ, R(·, aj)qj .ButAχJ, R(·, aj)=ZJ R(x, aj)A(x) dx =ZJ R(aj, x)A(x) dxsince R(x, y) = R(y, x) for x, y ∈J. Since R(y, x) = 0 for x ̸∈J the last integral equalsZ ∞−∞R(aj, x)A(x) dx = A(aj) −qj .Thus∂u∂aj= (−1)jq2j ,which is (1.8).

Equation (1.9) is derived analogously.We end this section with two relations between u, v, the qj and the pj which in fact canbe used to represent u and v directly in terms of the qj and pj. These relations are2v −u2 =Xj(−1)jq2j ,(2.18)u = −Xjp2j −ajq2j −2pjqju + 2q2jv.

(2.19)To obtain the first of these observe that (1.4) and (1.6) imply Xk∂∂ak!qj = pj −qju(2.20)while from (1.8) and (1.9)∂∂aj2v −u2= 2(−1)jqj (pj −qju) .13

Hence if we multiply both sides of (2.20) by 2(−1)jqj and sum over j we obtain Xk∂∂ak! Xj(−1)jq2j= Xk∂∂ak!

2v −u2.It follows that the two sides of (2.18) differ by a function of (a1, · · · , a2m) which is invariantunder translation by any vector (s, · · · , s). Since, clearly, both sides tend to zero as allai →∞, their difference must be identically zero.To deduce (2.19) we add (2.13) and (2.16) to obtain ∂∂x + ∂∂y +Xk∂∂ak!R(x, y) = −Q(x)Q(y)(x, y ∈J).

(Observe that the kernel of [D, (I −K)−1] is (∂/∂x + ∂/∂y)R(x, y).) Hence Xk∂∂ak!R(aj, aj) = −q2j = −(−1)j ∂u∂aj,the last by (1.8).

Multiplying by (−1)j and summing over j gives Xk∂∂ak! Xj(−1)jR(aj, aj)= − Xk∂∂ak!u .From this we deduce, by an argument similar to the one at the end of the last paragraph,thatXj(−1)jR(aj, aj) = −u .

(2.21)If we substitute formula (1.10) into this we see that the resulting double sum vanishes andwe are left with (2.19).B. Hamiltonian StructureThe partial differential equations (1.4)–(1.9) can be rewritten in a compact form thatreveals a symplectic structure with (1.10) defining a family of commuting Hamiltonians.

Tosee this we first change the notation slightlyq2j = −i2x2j,p2j = −iy2j,q2j−1 = 12x2j−1,p2j−1 = y2j−1,14

for j = 1, 2, . .

., m, andv = 12x0,u = y0,and introduce the canonical symplectic structure{xj, xk} = {yj, yk} = 0,{xj, yk} = δjk,for j, k = 0, 1, . .

., m. Then equations (1.4)–(1.9) becomedaxj = {xj, ω(a)}anddayj = {yj, ω(a)}(j = 0, 1, . .

., m)(2.22)whereω(a) = da log det (I −K)=2mXj=1Gj(x, y) daj(2.23)andGj(x, y) = y2j −14ajx2j −xjyjy0 + 14x2jx0 −142mXk=1k̸=j(xjyk −xkyj)2aj −ak.Furthermore, the Gj’s are in involution{Gj, Gk} = 0 .This last result can be verified by a direct calculation, but as we discuss below, there is abetter way to understand why such Hamiltonians are in involution. We can summarize thesystem of equations by saying, in words, that to find the ak-flow of the coordinates xj andyj, one flows according to Hamilton’s equation with Hamiltonian Gk.

The consistency ofthese equations (i.e. mixed partials are equal) follows immediately from this Hamiltonianformulation (see, e.g.

[28]).When equations (1.4)–(1.9) are expressed in the Hamiltonian form (2.22), the equationsare exactly analogous to the case of the sine kernel [17] and in both cases the Hamiltoniansare defined via (2.23). For the JMMS equations (see, e.g.

[28]) this Hamiltonian system is15

closely related to the Hamiltonian system studied by Moser [24]. In fact in both cases thesystems of partial differential equations can be formulated in terms of the classical R-matrixapproach to completely integrable systems.

Such an approach identitfies the Lax operatorand hence the commuting flow result above follows without any calculations. We refer thereader to [12] for a discussion of this classical R-matrix formulation of these systems ofpartial differential equations.III.

THE ORDINARY DIFFERENTIAL EQUATIONSIn this section we specialize to the case J = (s, ∞) and derive the differential equations(1.11), (1.13) and (1.14). In the notation of the last section m = 1, a1 = s, a2 = ∞.

Allterms in the formulas we derived which contain a2 vanish because of the decay of the Airykernel at +∞. Recall that we now denote the quantity q corresponding to the endpointa1 = s by q(s), and write R(s) for R(s, s).First, equation (1.14) is immediate from (2.21), which now says that R = u, and (1.8).For (1.13) we observe that (1.6) and (1.7) in the present case readq′ = p −qu ,(3.1)p′ = sq + pu −2qv ,(3.2)whereas (1.10) readsR = pq′ −qp′ .Thus, by (1.14),RR′ = pq!′.

(3.3)Taking the logarithmic derivative of both sides of (1.14) gives12R′′R′ = q′qand by (3.1) this equals16

pq −u = pq −R .Thus12 R′′R′!′= pq!′−R′ .By (3.3) this is just (1.13).Finally, to obtain (1.11) we differentiate (3.1) and use (3.1), (3.2), (1.8) and (1.9). Wefind by an easy computation thatq′′ = sq + q3 + q(u2 −2v) .Conveniently, we have (2.18), which says in this case that u2 −2v = q2.

Thus (1.11) isestablished.IV. ASYMPTOTICSA.

Asymptotics via Painlev´eHastings and McLeod [13] (see also [5,6]) have shown that q(s; 1), the unique solution toequation (1.11) which is asymptotic to Ai(s) as s →+∞, is asymptotically equal toq−s/2as s →−∞. In fact there is a complete asymptotic expansion in decreasing powers of −s,beginnning with this term [19].

Successive terms are easily computed. This result and laterones, become a little easier to state if we write s = −t/2:q(−t/2; 1) = 12√t1 −1t3 −732t6 −106572t9−139122778t12+ O( 1t15)(t →+∞)(4.1)Squaring and integrating and using (1.14) we deduceR(s) = 14s2 + c −18 s +964 s4 −189128 s7 + 21663512 s10 + O( 1s13)(s →−∞)where c is a constant.

Substituting this into (1.13) shows that c = 0. Thus we obtain (1.16)and then, after integrating and exponentiating, (1.17).17

To obtain the asymptotic expansions of qn(s) in (1.18), we differentiate (1.11) withrespect to λ and set λ = 1. We obtainq′′1 −s + 6q20q1 = 0(where q0(s) = q(s; 1)).

Using (4.1) we find two linearly independent solutions of this equa-tion, one exponentially large at −∞and one exponentially small. Assuming the exponen-tially large one actually appears with a nonzero factor we obtain the asymptotic expansionq1(−t/2) = c1exp( 13t3/2)t1/41 +1724 t3/2 + 15132732 t3 + 85019321034 t9/2 + 40711752121535 t6+ · · ·(4.2)where c1 is a constant factor to be determined.Successive differentiations of (1.11) with respect to λ, followed by setting λ = 1, yields asequence of linear differential equations for the qn,q′′n −s + 6q20qn = ϕn(q0, · · · , qn−1),where ϕn is a polynomial in the qm with m < n. For n > 1 a particular solution dominatesall solutions of the homogeneous equation.

Therefore, in the asymptotics, no new unknownconstant factors are introduced (unless we want to go “beyond all orders”).To deduce expansions for the r(n; s) we must differentiate (1.15) n times with respect toλ and set λ = 1. We find expressions of the formr(n; s) =Z ∞s (x −s)ψ(q0, · · ·, qn) dx(4.3)where ψn is a polynomial in q0, · · ·, qn.

In particular we obtain (after setting s = −t/2)r(1; −t/2) = 12Z t−∞(t −x)q0(−x/2)q1(−x/2) dx .Substituting (4.1) and (4.2) into this integral gives the asymptotic expansion for r(1; −t/2).The first order result isr(1; −t/2) ∼c1exp(t3/2/3)t3/418

and comparing this with (1.22) reveals thatc1 =12√2π .We now state the expansions we have obtained in this way with the help of Mathematica.For qn with n > 2, as t →+∞qn(−t/2) =n!23n/2 πn/2exp( n3t3/2)t3n/4−1/2 1 + 59n −3624t−3/2 +3481n2 + 8244n −84961152t−3 + O(t−9/2)!. (4.4)For n = 1, 2 the leading order term is as above.

For n = 1 the correction terms are givenabove in (4.2). For n = 2 the coefficient of t−3/2 is as given above but the coefficient of t−3is 5461/288.

(Similar exceptional values of n occurred in analogous expansions for the sinekernel [2].) These expansions have been computed for n ≤10 and are undoubtedly correctfor all n. (In [2] the analogous expansions were proved for all n by a backwards inductionargument—such arguments should extend to the present case of qn.

)For r(n; s) we have obtained the following sharpening of (1.22): as t →+∞r(n; −t/2) = 1!2! · · ·(n −1)!2n2+n/2πn/2exp( n3t3/2)t3n2/4 1 + n(34n2 + 31)24t−3/2 +n2(1156n4 + 6608n2 + 11509)1152t−3 + O(t−9/2)!.

(4.5)This was established for n ≤6. As was mentioned in the Introduction, when we substituteexpansions (4.4) into (4.3) a large amount of cancellation takes place, so that r(n; s) is muchsmaller than might be anticipated.This also means that to deduce even the first-orderasymptotics in (4.5) for moderate values of n one has to go very far out in the expansions(4.4).

Thus it is important that we have an alternative approach which yields asymptoticresults for all values of n.B. Asymptotics via the Commuting Differential OperatorWe begin with the integral identity19

K(x, y) =Z ∞0A(x + z)A(y + z) dz . (4.6)This appears in [5] but we give a simple proof here.

Think of (2.12) as a differential equationfor the “unknown” function K(x, y). Since ∂∂x + ∂∂y!A(x + z)A(y + z) = ∂∂z A(x + z)A(y + z)one solution is the right side of (4.6) and, of course, another solution is the Airy kernelK(x, y).

The difference is therefore of the form ϕ(x −y) for some function ϕ. Since bothsides of (4.6) tend to zero as x and y tend to +∞independently, we must have ϕ ≡0.Rewriting (4.6) asK(x, y) =Z ∞sA(x + z −s)A(z + y −s) dzshows that K as an operator on (s, ∞) is the square of the operator with kernel A(x+y−s).Integration by parts shows that an integral operator with kernel L(x, y) on (s, ∞) com-mutes with a differential operatorddxα(x) ddx + β(x)if α(s) = 0 and ifα(y)∂2L∂y2 + α′(y)∂L∂y + β(y)L = α(x)∂2L∂x2 + α′(x)∂L∂x + β(x)L .

(Of course we also require, in the end, appropriate vanishing at ∞.) If we set L(x, y) =A(x + y −s) and use the fact A′′(x) = xA(x) we can easily check that a solution is given byα(x) = x −s,β(x) = −x(x −s)and so we obtain the commuting differential operator (1.20).Here is an outline of how (1.21) is found.

If we let f1, f2, · · · denote the eigenfunctionsof L then the logarithmic derivative of each λi (thought of as a function of s) is expressiblevery simply in terms of the corresponding fi. The asymptotics of these eigenfunctions, andso of the logarithmic derivative, is established by WKB techniques.

Integrating, and using20

the fact that each λi →1 as s →−∞, gives (1.21). The WKB argument we give is quiteheuristic although it could very likely be made rigorous.We make a preliminary change of variables, replacing x by x + s (so that L acts on(0, ∞)) and s by −2t.

The eigenvalue problem can then be written as(xf ′(x))′ +µ −(x −t)2f(x) = 0(4.7)or in the equivalent formxg′′(x) +µ −(x −t)2 + 14xg(x) = 0 ,g(x) = √xf(x)(4.8)The eigenvalues satisfy 0 < µ0 < µ1 < · · · and oscillation considerations applied to (4.8)show that µi = O(t1/2), for each i as t →∞. So we think of µ in our equations as O(t1/2),we normalize our eigenfunctions so that f(0) = 1, and proceed to find asymptotics for larget.1.

The region x << t−1/2If we make the substitutionsf(x) = h(t2x) = h(y)then (4.7) becomes(yh′(y))′ −h(y) = t−2 y2t4 −2yt −µ!h(y) .The solution of the equation with the right hand side replaced by 0 which satisfies h(0) = 1is I0(2√y) where I0 is the usual modified Bessel function. The actual solution h(y) will beasymptotic to this in any interval over which the integral of the factor of h(y) on the righthand side of the equation is o(1).

Recalling that µ = O(t1/2) we see that y << t3/2 suffices.Hencef(x) ∼I02t√x,(x << t−1/2)(4.9)21

2.The region x >> t−2, t −x >> t1/4If in equation (4.8) we setα(x) = (t −x)2 −µx,y =Z x0qα(z) dz ,h(y) = α(x)1/4g(x)then the equation for h(y) ish′′(y) −h(y) ="14α′′α2 −516α′2α3 −14x2α#h(y) ,where on the right side ′ =ddx. Since dy = α(x)1/2 dx, the integral with respect to y (of thefactor of h(y) on the right side of the equation) over any interval equalsZ 14α′′α3/2 −516α′2α5/2 −14x2α1/2!dxover the corresponding x interval.

This will be o(1) as long as, in this interval, x >> t−2and t −x >> t1/4.In this same region y >> 1 and soh(y) ∼aey(4.10)for some constant a = a(t). In case x << 1 we have y = 2t√x + o(1) and sof(x) = x−1/2g(x) ∼at−1/2x−1/4e2t√x ,(t−2 << x << 1).Comparing this result with (4.9) in the overlapping region t−2 << x << t−1/2 and using theknown asymptotics of I0(z) as z →∞shows that a ∼1/2√π.

Hence from (4.10)f(x) ∼12√π(t −x)−1/2x1/4expZ x0s(t −z)2 −µzdz,(x >> t−2, t −x >> t1/4) .We also used here once again the estimate µ = O(t1/2). Because of this same estimate wecan write the integral in the exponential asZ x0t −z√z(1 −12µ(t −z)2 + O µ2(t −z)4!

)dz .22

Because t −x >> t1/4 the last term in brackets contributes o(1) to the integral, and so weobtain2tx1/2 −23x3/2 −µ2Z x0dz√z(t −z) .Replacing z by tz2 shows that the integral here equalst−1/2 log√t + √x2t −x.Thusf(x) ∼12√π(t −x)12 t−1/2µ−12x1/4√t + √xt−1/2 µ exp2tx1/2 −23x3/2(4.11)(x >> t−2, t −x >> t1/4).Let us see what happens if x is not far from t, ift1/4 << t −x << t1/2 .Then since2tx1/2 −23x3/2 = 43t4/3 −12(x −t)2 + o(1)in this region, we findf(x) ∼12√π2−t−1/2µt−3/4tt −x−12 t−1/2µ+ 12 exp43t4/3 −12(x −t)2(4.12)(t1/4 << t −x << t1/2)3. The region x −t >> t1/4This is quite similar to the preceding.

The main difference is that instead of making thevariable change y =R α(x)1/2 dx in (4.8) we now sety = 23x3/2 −2tx1/2 −Z ∞xqα(z) −z1/2 + tz−1/2dz .23

Since an eigenfunction cannot grow exponentially at ∞only the negative exponential e−ycan appear in the asymptotics and we find thatf(x) ∼b(x −t)−1/2x1/4exp2tx1/2 −23x3/2 +Z ∞xs(z −t)2 −µz−z1/2 + tz−1/2dzfor some constant b = b(t). Now the integral in the exponential is−µ2 t−1/2 log√x +√t2x −t+ o(1)and sof(x) ∼b(x −t)12t−1/2µ−12x1/4√x +√tt−1/2µ exp2tx1/2 −23x3/2(x −t >> t1/4)(4.13)and specializing to the region t1/4 << x −t << t1/2 givesf(x) ∼b 2−t−1/2µt−3/4tx −t−12 t−1/2µ+ 12 exp43t3/2 −12(x −t)2(4.14)(t1/4 << x −t << t1/2)4.

The region |x −t| << tIf in equation (4.8) we setg(x) = h√2t−1/4(t −x)= h(y)and write ε = (√2t3/4)−1 then the equation becomes(1 −εy) h′′(y) +" µ2t1/2 −y24 +ε24(1 −εy)#h(y) = 0 .Now we expect that for bounded y, and therefore also if |y| →∞at some rate (dependingupon ε), the solution of this will be asymptotic to a solution ofh′′(y) + µ2t1/2 −y24!h(y) = 0 ,Weber’s equation. Now unless24

µ2t1/2 = i + 12(i = 0, 1, · · ·)there is no solution which vanishes at both ±∞. But (4.12) and (4.14) show that in theirregions of validity (which overlap with our present region) f(x) is a factor depending ont times a factor tending exponentially to 0 as x −t →±∞.

Hence our eigenvalues mustsatisfyµi2t1/2 = i + 12 + o(1)and h(y) is asymptotically a constant times the parabolic cylinder function Di(y) (see, e.g. [4], Chp.

8). Thus (since x ∼t in the present situation)f(x) ∼c Di√2t−1/4(t −x)(4.15)for some constant c = c(t).

Using the known asymptotics of Di we findf(x) ∼c√2t−1/4(t −x)i exp−12t−1/2(t −x)2(4.16)Comparing this with (4.12) givesc ∼2−52i−2√π t−34 i−34 exp43t3/2. (4.17)We also find, using (4.14), that b ∼1/2√π.5.

The asymptotics of λiIf we look at the asymptotics of fi(x) (the eigenfunction associated with the eigenvalueµi of L, normalized so that fi(0) = 1) given by (4.9), (4.11), (4.13), and (4.16) we see thatthe main contribution toR ∞0 fi(x)2 dx comes from any regiont−1/4|t −x| < η(t)as long as η(t) →∞. For some such region we have, by (4.15) and (4.17),fi(x) ∼2−52 i−2√π t−34i−34Di√2t−1/4(t −x).25

Using this, and the factR ∞−∞Di(x) dx =√2πi!, we deduceZ ∞0fi(x)2 dx ∼i!√π2−5i−4t−32 i−54 exp83t3/2. (4.18)This was the goal of all that went before.Since L has simple eigenvalues and commutes with our integral operator K (translatedto (0, ∞)) the set {fi} is the set of eigenfunctions of K corresponding to its eigenvalues {λi}.However fi corresponds to µi and there is no assurance that the corresponding eigenvalueof K is its ith largest, which we have denoted by λi.

We shall eventually show that this isactually the case. The first lemma of this section shows that the eigenvalues of K are simple,and this is crucial for the argument we shall later use.

(The proof of the corresponding resultfor the sine kernel is in [27]. )Lemma 1 If f1 and f2 are eigenfunctions of L corresponding to distinct eigenvalues, thenthey are eigenfunctions of K corresponding to distinct eigenvalues.Proof.

The eigenvalues of K are the squares of the eigenvalues of the operator on (0, ∞)with kernel A(x + y −2t). So we have to show that if for some ν we have eitherZ ∞0A(x + y −2t)fi(y) dy = νfi(x)(i = 1, 2)(4.19)orZ ∞0A(x + y −2t)f1(y) = νf1(x) ,Z ∞0A(x + y −2t)f2(y) dy = −νf2(x)(4.20)then f1 = f2.

(Recall that all eigenfunctions are normalized to satisfy f(0) = 1. )Assuming first that (4.19) holds, apply d2/dx2 to both sides of the identity with i = 1and integrate the resulting integral by parts twice.

What results isνf ′′1 (x) = −A′(x −2t) + A(x −2t)f ′1(0) +Z ∞0A(x + y −2t)f ′′1 (y) dy .If we multiply both sides by f2(x) and integrate we obtain, using (4.19) with i = 2 (and itsdifferentiated version),26

νZ ∞0f ′′1 (x)f2(x) dx = −νf ′2(0) + νf ′1(0) + νZ ∞0f2(y)f ′′1 (y) dy .Thus, since ν ̸= 0 (this follows easily from the fact that the Fourier transform of A(x) isnonzero), we have f ′1(0) = f ′2(0). But it is clear from (4.7) that for any eigenfunction f wehave f ′(0) = t2 −µ.

Thus f1 and f2 correspond to the same eigenvalue of L and so are equal.The case when (4.20) holds is even easier. Differentiating both sides of the first relationonce and integrating by parts givesνf ′1(x) = −A(x −2t) −Z ∞0A(x + y −2t)f ′1(y) dy .Multiplying both sides of this by f2(x) and integrating, using the second relation of (4.20),we obtainνZ ∞0f ′1(x)f2(x) dx = −ν + νZ ∞0f2(y)f ′1(y) dy ,contradicting ν ̸= 0.Since, as we now know, the eigenvalues of K are simple, there is a permutation σ ofN = {0, 1, 2, · · ·} such that fi is the eigenfunction of K corresponding to its σ(i)th largesteigenvalue λσ(i); this permutation is independent of t or else, by the continuity of the eigen-values in t, there would be a multiple eigenvalue for some t.Up to now our kernel K corresponded to an arbitrary factor λ.

Of course formula (1.21)refers to the case λ = 1.Lemma 2 For each i we have λi →1 as s →−∞.Proof. In this case λ = 1 formula (4.6) readsK(x, y) =Z ∞0Ai(x + z)Ai(z + y) dz .

(4.21)For the purposes of this proof we denote by A (resp. K) the operator on (−∞, ∞) withkernel Ai(x + y) (resp.

K(x, y)) and by Ps the projection from L2(−∞, ∞) to L2(s, ∞). Sothe λi are the eigenvalues of PsKPs.

Now it is known [13] that A2 = I, and (4.21) saysthat K = AP0A. Thus K2 = AP0A2P0A = AP0A = K. In other words, K is a projection27

operator. (This is not surprising since its kernel is a scaling limit of the kernels KN(x, y) ofprojection operators.) Clearly it has infinite-dimensional range.

It now follows easily fromthe minimax characterization of the eigenvalues that, for each i, the ith largest eigenvalue ofPsKPs tends to 1 as s →−∞.We now state a lemma which is the basis of a trick used by Fuchs [11] in his derivationof the asymptotics of the eigenvalues of the sine kernel.Lemma 3 Let J(x, y) be a symmetric kernel, on a fixed interval, depending smoothly on aparameter t. Let g be an eigenfunction normalized so thatR g(x)2 dx = 1 and let λ be thecorresponding eigenvalue. Then, with the subscript t denoting ∂/∂t, we haveλt =Z ZJt(x, y)g(y)g(x) dydx .Proof.

FromZJ(x, y)g(y) dy = λg(x)(4.22)we getZJt(x, y)g(y) dy +ZJ(x, y)gt(y) dy = λtg(x) + λgt(x) .Multiplying both sides by g(x) and integrating over x, we obtain (using the normalizationof g and (4.22))Z ZJt(x, y)g(y)g(x) dydx +Z ZJ(x, y)gt(y)g(x) dxdy = λt + λZgt(x)g(x) dx= λt +Z ZJ(x, y)g(y)gt(x) dydx .By symmetry of J the two double integrals involving it cancel, and we obtain the statementof the lemma.We apply the lemma to the opertor with kernel A(x + y −2t) on (0, ∞) whose squareis our Airy operator (translated to act on (0, ∞)). Thus we may write its eigenvalues asqλσ(i), where the square roots might have either sign.

The lemma gives28

∂∂tqλσ(i) = −2Z ∞0Z ∞0A′(x + y −2t)gi(y)gi(x) dydx(4.23)wheregi(x) = fi(x)Z ∞0fi(y)2 dy1/2.But fromZ ∞0A(x + y −2t)gi(y) dy =qλσ(i)gi(x)we obtain upon differentiation thatZ ∞0A′(x + y −2t)gi(y) dy =qλσ(i)g′i(x)so (4.23) may be rewritten∂∂tqλσ(i) = −2qλσ(i)Z ∞0g′i(x)gi(x) dx .The integral is of course −12gi(0)2. Recalling that fi(0) = 1 we see that we have establishedthe relation∂∂tqλσ(i) =qλσ(i)Z ∞0fi(x)2 dx,or∂∂t log λσ(i) = 2Z ∞0fi(x)2 dxWe now use (4.18) and integrate the resulting relation with respect to t from t to ∞.Recalling Lemma 2, we deduce−log λσ(i) ∼√πi!

25i+3t32i+ 34 exp−83t3/2,precisely the right side of (1.21).All that remains now is to show that σ(i) = i for all i. But it is clear from this asymptoticrelation that i < j implies λσ(i) > λσ(j) for large t (and so for all t), so that σ(i) < σ(j).Thus σ is order preserving, and since σ : N →N is onto we must have σ(i) = i for all i.This concludes the derivation of (1.21).

The deduction of (1.22) from (1.21) is completelyanalogous to the corresponding deduction for the sine kernel in ref. [28], Sec.

VIIB.29

ACKNOWLEDGMENTSWe wish to thank E. Br´ezin, P. J. Forrester, F. A. Gr¨unbaum and J. Harnad for helpfulcomments. This work was supported in part by the National Science Foundation, DMS–9001794 and DMS–9216203, and this support is gratefully acknowledged.30

REFERENCES[1] M. J. Ablowitz and H. Segur, Exact linearization of a Painlev´e transcendent, Phys. Rev.Letts.

38 (1977) 1103–1106. [2] E. L. Basor, C. A. Tracy, and H. Widom, Asymptotics of level spacing distributions forrandom matrices, Phys.

Rev. Letts.

69 (1992) 5–8. [3] M. J. Bowick and E. Br´ezin, Universal scaling of the tail of the density of eigenvaluesin random matrix models, Phys.

Letts. B268 (1991) 21–28.

[4] A. Erd´elyi (ed. ), Higher Transcendental Functions, Vol.

II (McGraw-Hill, New York,1953). [5] P. A. Clarkson and J.

B. McLeod, A connection formula for the second Painlev´e tran-scendent, Arch. Rat.

Mech. Anal.

103 (1988) 97–138. [6] P. A. Clarkson and J.

B. McLeod, Integral equations and connection formulae for thePainlev´e equations, in Painlev´e Transcendents: Their Asymptotics and Physical Appli-cations, eds. D. Levi and P. Winternitz (Plenum Press, New York, 1992), pgs.

1–31. [7] F. J. Dyson, Statistical theory of energy levels of complex systems, I, II, and III, J.Math.

Phys. 3 (1962) 140–156; 157–165; 166–175.

[8] F. J. Dyson, Fredholm determinants and inverse scattering problems, Commun. Math.Phys.

47 (1976) 171–183. [9] F. J. Dyson, The Coulomb fluid and the fifth Painlev´e transcendent, IASSNSS-HEP-92/43 preprint, to appear in the proceedings of a conference in honor of C. N. Yang, ed.S.-T.

Yau. [10] P. J. Forrester, The spectrum edge of random matrix ensembles, preprint.

[11] W. H. J. Fuchs, On the eigenvalues of an integral equation arising in the theory ofband-limited signals, J. Math.

Anal. and Applic.

9 (1964) 317–330.31

[12] J. Harnad, C. A. Tracy, and H. Widom, Hamiltonian structure of equations appearingin random matrices, in preparation. [13] S. P. Hastings and J.

B. McLeod, A boundary value problem associated with the secondPainlev´e transcendent and the Korteweg-de Vries equation, Arch. Rat.

Mech. Anal.

73(1980) 31–51. [14] E. L. Ince, Ordinary Differential Equations, (Dover, New York, 1956).

[15] A. R. Its, A. G. Izergin, V. E. Korepin, and N. A. Slavnov, Differential equations forquantum correlation functions, Int. J. Mod.

Physics B 4 (1990) 1003–1037. [16] K. Iwasaki, H. Kimura, S. Shimomura, and M. Yoshida, From Gauss to Painlev´e: AModern Theory of Special Functions (Vieweg, Braunschweig, 1991).

[17] M. Jimbo, T. Miwa, Y. Mˆori, and M. Sato, Density matrix of an impenetrable Bose gasand the fifth Painlev´e transcendent, Physica 1D (1980) 80–158. [18] B. M. McCoy, C. A. Tracy, and T. T. Wu, Connection between the KdV equation andthe two-dimensional Ising model, Phys.

Letts. 61A (1977) 283–284.

[19] J. B. McLeod, private communication.

[20] M. L. Mehta, Random Matrices, 2nd edition (Academic, San Diego, 1991). [21] M. L. Mehta, A non-linear differential equation and a Fredholm determinant, J. dePhys.

I France, 2 (1992) 1721–1729. [22] M. L. Mehta and G. Mahoux, Level spacing functions and non-linear differential equa-tions, preprint.

[23] G. Moore, Matrix models of 2D gravity and isomonodromic deformation, Prog. Theor.Physics Suppl.

No. 102 (1990) 255–285.

[24] J. Moser, Geometry of quadrics and spectral theory, in Chern Symposium 197932

(Springer, Berlin, 1980), 147–188. [25] P. Painlev´e, Sur les ´equations diff´erentielles du second ordre et d’ordre sup´erieur dontl’int´egrale g´en´erale est uniforme, Acta Math.

25 (1902) 1–85. [26] C. E. Porter, Statistical Theory of Spectra: Fluctuations (Academic, New York, 1965).

[27] D. Slepian and H. O. Pollak, Prolate spheroidal wave functions, Fourier analysis anduncertainty-I, Bell Systems Tech. J.

40 (1961) 43–64. [28] C. A. Tracy and H. Widom, Introduction to random matrices, to appear in the pro-ceedings of the 8th Scheveningen Conference, Springer Lecture Notes in Physics.

[29] C. A. Tracy and H. Widom, Level spacing distributions and the Airy kernel, submittedto Phys. Letts.

B. [30] H. Widom, The strong Szeg¨o limit theorem for circular arcs, Indiana Univ.

Math. J.

21(1971) 277–283. [31] H. Widom, The asymptotics of a continuous analogue of orthogonal polynomials, toappear in J. Approx.

Th.33

FIGURESFIG. 1.

The probabilities E(0; s) and E(1; s). Of course, E(0; s) →1 as s →+∞.FIG.

2. The probability densities F(0; s) and F(1; s).

The expected position of the largesteigenvalue is approximately −1.7711 with variance 0.8132. The expected position of the next-largesteigenvalue is approximately −3.6754 with variance 0.5405.34


출처: arXiv:9211.141원문 보기

Subscribe to koineu.com

Don’t miss out on the latest issues. Sign up now to get access to the library of members-only issues.
jamie@example.com
Subscribe