Tracy-Widom method for Janossy density and joint distribution of extremal eigenvalues of random matrices

The J\'{a}nossy density for a determinantal point process is the probability density that an interval $I$ contains exactly $p$ points except for those at $k$ designated loci. The J\'{a}nossy density associated with an integrable kernel $\mathbf{K}\doteq (\varphi(x)\psi(y)-\psi(x)\varphi(y))/(x-y)$ is shown to be expressed as a Fredholm determinant $\mathrm{Det}(\mathbb{I}-\tilde{\mathbf{K}}|_I)$ of a transformed kernel $\tilde{\mathbf{K}}\doteq (\tilde{\varphi}(x)\tilde{\psi}(y)-\tilde{\psi}(x)\tilde{\varphi}(y))/(x-y)$. We observe that $\tilde{\mathbf{K}}$ satisfies Tracy and Widom's criteria if $\mathbf{K}$ does, because of the structure that the map $(\varphi, \psi)\mapsto (\tilde{\varphi}, \tilde{\psi})$ is a meromorphic $\mathrm{SL}(2,\mathbb{R})$ gauge transformation between covariantly constant sections. This observation enables application of the Tracy--Widom method to J\'{a}nossy densities, expressed in terms of a solution to a system of differential equations in the endpoints of the interval. Our approach does not explicitly refer to isomonodromic systems associated with Painlev\'{e} equations employed in the preceding works. As illustrative examples we compute J\'{a}nossy densities with $k=1, p=0$ for Airy and Bessel kernels, related to the joint distributions of the two largest eigenvalues of random Hermitian matrices and of the two smallest singular values of random complex matrices.

The Jánossy density for a determinantal point process is the probability density that an interval I contains exactly p points except for those at k designated loci. The Jánossy density associated with an integrable kernel K . = (ϕ(x)ψ(y) − ψ(x)ϕ(y))/(x − y) is shown to be expressed as a Fredholm determinant Det(I −K| I ) of a transformed kernel K . = (φ(x)ψ(y) −ψ(x)φ(y))/(x − y). We observe thatK satisfies Tracy and Widom's criteria if K does, because of the structure that the map (ϕ, ψ) → (φ,ψ) is a meromorphic SL(2, R) gauge transformation between covariantly constant sections. This observation enables application of the Tracy-Widom method to Jánossy densities, expressed in terms of a solution to a system of differential equations in the endpoints of the interval. Our approach does not explicitly refer to isomonodromic systems associated with Painlevé equations employed in the preceding works. As illustrative examples we compute Jánossy densities with k = 1, p = 0 for Airy and Bessel kernels, related to the joint distributions of the two largest eigenvalues of random Hermitian matrices and of the two smallest singular values of random complex matrices.

Introduction
In the history of random matrix theory (RMT), which models the local fluctuation of energy levels of quantum chaotic and/or disordered Hamiltonians typified by the Sinai billiard, the Anderson tight-binding model and the QCD Dirac operator, Gaudin and Mehta's discovery that the distribution of ordered eigenvalues or their spacings is expressed in terms of the Fredholm determinant or Pfaffian of an integral kernel K I restricted on an interval I [1,2] has been known as long as the RMT itself. Specifically, the distribution P k (s) of the kth largest eigenvalue (centered and scaled) of random Hermitian matrices is given as with K being the Airy kernel [3], and that of the kth smallest singular values of random complex matrices is given by Eq. (1) with K being the Bessel kernel [4] (with replacements I = (s, ∞) → (0, s) and ∂ s → −∂ s ). These trains of peaks that gradually approach Gaussian in the spectral bulk [5] constitute the spectral densities ρ 1 (s) = ∞ k=1 P k (s), as plotted in Fig. 1. For the practical purpose of fitting some spectral data to the RMT to extract systemspecific constants (such as the chiral condensate and the pion decay constant in the case of QCD Dirac operators [6]), characteristic peaky shapes of the individual distributions are  Fig. 1 Distributions P k (s) of the scaled kth largest eigenvalues of random Hermitian matrices (Tracy-Widom distribution) (right) and of the kth smallest singular values of random complex square matrices (left) (red (k = 1) to blue (k = 8)), their sums 8 k=1 P k (s) (grey dotted), and the spectral densities ρ 1 (s) (black).
better suited than the spectral density, as the oscillation of the latter tends to smooth out in the bulk and the data-fitting would yield little more than the mean level density.
With this in view, the purpose of this article is to advance the formula (1) a step further and provide a "user-friendly" analytic method to compute the joint distribution P 1···k (s 1 , . . . , s k ) of the first to kth largest/smallest eigenvalues of unitary-invariant random matrices, which is a constituent of the k-point correlation function ρ k (s 1 , . . . , s k ). To this end we apply the strategy of Tracy and Widom [7] on the evaluation of Fredholm determinants of integrable integral kernels to the Jánossy density [8][9][10][11], i.e., the probability distribution that an interval contains no eigenvalue except for those at k designated loci. As the simplest examples we shall evaluate the joint distributions P 12 (t, s) of the first and second largest eigenvalues and smallest singular values (see Fig. 2 for their histograms), i.e., the first peak that constitutes the two-point correlation function ρ 2 (t, s) = k< P k (t, s), for the Airy and Bessel Fig. 2 Histograms of the first (t) and second (s) largest eigenvalues of random Hermitian matrices (left) and of the first (t) and second (s) smallest singular values of random complex square matrices (right). Matrix rank N = 128 and number of samples = 10 7 for each case, and eigen/singular values x are rescaled as: t or s = √ 2N 1/6 (x − √ 2N ) and √ 2N x, respectively.
2/18 kernels. Each case has been worked out previously in Refs. [12,13], which devised an elaborate analytic procedure involving the Painlevé II and III transcendents and the associated isomonodromic systems [14]. This approach is later simplified (for the Airy kernel) using a solution to the Lax pair associated with the Painlevé XXXIV system [15]. Our alternative method presented in this article, which does not explicitly refer to these systems and employs the familiar Tracy-Widom method, has a clear advantage of permitting straightforward generalizations to a general P 1···k and/or to various finite-N and large-N kernels (Hermite, Laguerre, and other hypergeometric; circular, beyond-Airy, q-orthogonal, etc.) appearing in the RMT. This article is composed of the following parts: In Sect. 2 we list known formulas on Jánossy densities of a determinantal point process, and then express them in terms of Fredholm determinants of the "transformed kernel"K. The latter is a novel presentation to the best of our knowledge, except for the simplest (k = 1) case of the sine kernel previously treated in Ref. [16]. In Sect. 3 we demonstrate thatK satisfies Tracy and Widom's criteria for their functional-analytic method to be applicable if the original K does. In Sect. 4 we evaluate Jánossy densities and joint distributions of the first and second extremal eigenvalues from the Airy and Bessel kernels by the Tracy-Widom method. In Sect. 5 we conclude with listing possible applications and extensions of our approach. Numerical data of Jánossy densities for the Airy kernel and the Bessel kernels at ν = 0, 1 are attached as supplementary material. Throughout this article we follow the notations of Ref. [7], hereafter denoted as TW.

Jánossy density
First we collect some facts on determinantal point processes (DPPs). Let X be a countable set. Consider an ensemble of finite subsets of X consisting of N elements ("particles") (n 1 , . . . , n N ), and assign to them a joint probability in a determinantal form: Here K = [K(n, m)] n,m∈X is an operator in the Hilbert space L 2 (X), i.e., an infinitedimensional matrix indexed by the points of X. The operator K, which we also call a kernel, is required to be real, symmetric, projective, and normalized: These requirements lead to the k-point correlation function ρ k (n 1 , . . . , n k ), i.e., the joint probability that k particles occupy the points n 1 , . . . , n k , to be given in a determinantal form as well: e.g., page 341 of Ref. [9], where it is denoted as π(X)): We interchangeably use notations A(n, m) = n|A|m for the (n, m)-element of a matrix A.
Equation (6) indicates that the Jánossy density is a gap probability for a DPP with a transformed kernelK multiplied by the k-point correlation function det κ. The projectivity and the normalization conditions (3) with N → N − k can be verified forK in a straightforward manner. Note that ρ p (m 1 , . . . , m p ; n 1 , . . . , n k ) := ρ p+k (m 1 , . . . , m p , n 1 , . . . , n k ) ρ k (n 1 , . . . , n k ) represents the conditional joint probability that p particles occupy the points m 1 , . . . , m p under the presumption that k particles already occupy the points n 1 , . . . , n k . Obviously, this fact could as well be deduced from the very definition of the conditional joint probabilities; 4/18 e.g., for k = 1, the transformed kernelK(m, m ) = K(m, m ) − K(m, n)K(n, m ) K(n, n) satisfies , etc., without any recourse to Jánossy densities. We could have reversed the course of derivation of Jánossy densities (6) and started backward from Eq. (8).
Generalization to the probability J k,p (n 1 , . . . , n k ; I) that there are exactly p particles in I except for k particles, one at each of the k designated loci, is straightforward (see part (3) in Fig. 3); we introduce a parameter z so that J k,p (n 1 , . . . , n k ; I) is given by as the derivation of Eq. (1) carries over to this case. The product of two determinants in Eq. (10) can as well be written as Either by definition or from Eqs. (10) and (11), the Jánossy density J k,p (n 1 , . . . , n k ; I) reduces to: (i) for k = 0, the gap probability E 0 (I) = det(I − K I ) (or its generalization E p (I) = 1/p! (−∂ z ) p det(I − zK I )| z=1 ) of finding no (or exactly p) particles in I, and (ii) formally for I = ∅ and p = 0, to the k-point correlation function (4). Note also that the factor det κ = ρ k (n 1 , . . . , n k ) in Eqs. (6) and (11) is canceled when we consider the conditional probabilitỹ J k,p (n 1 , . . . , n k ; I) that there are exactly p particles in a subset I under the condition that k particles are already at each of the k designated loci, All the above formulas carry over to a continuous DPP on R and for a set of intervals I ⊂ R. Trivial modifications are needed to regard K and K I as integral operators acting on the Hilbert space of square-integrable functions L 2 (R) and L 2 (I), and to reinterpret joint probabilities ρ k ,ρ k , J k,p ,J k,p as joint probability distributions. For the case of continuous DPPs the expression x|A|y = A(x, y) (denoted as A . = A(x, y) in TW) means that the integral operator A has a kernel equal to A(x, y). Namely, the Jánossy density J k (x 1 , . . . , x k ; I) is defined as the probability density of finding exactly k particles in I and one at each of the k infinitesimal intervals (x i , x i + dx i ) ⊂ I, and is given by the Fredholm 5/18 determinant Det (I − K I ) times the (ordinary) determinant of the resolvent kernel of K I : Likewise its generalization J k,p (x 1 , . . . , x k ; I) is given by Finally we note that the joint probability distribution of k leftmost or rightmost particles is derived from the Jánossy density for a semi-finite interval I = (s k , ∞) or (−∞, s k ) as 3. Applicability of the Tracy-Widom method

Inheritance of the Tracy-Widom criteria
Consider a kernel of an integral operator K of the Christoffel-Darboux form with its component functions satisfying a pair of linear differential equations with some polynomials m, A, B, C.
The tracelessness of the 2 × 2 matrix on the right-hand side of Eq. (18) is essential. As a unifying approach to their preceding works on the sine [18], Airy [3] and Bessel kernels [4], Tracy and Widom have shown in TW that the Fredholm determinant Det(I − K I ) of an operator K satisfying the criteria (17), (18) is always determined through a closed system of PDEs in the boundary points {a i } ∈ ∂I. This involves the boundary values of the , and the inner products of Q j and P j with ϕ and ψ such as Now we present a theorem: Theorem 1. If the kernel of K satisfies the TW criteria (17), (18), so does the transformed kernel ofK.
Proof. Since the kernel ofK (k) for the Jánossy density J k (x 1 , . . . , x k ; I) is obtained from the kernel ofK (k−1) for the Jánossy density J k−1 (x 1 , . . . , x k−1 ; I) by adding an extra locus 6/18 of particle x k = t, by induction it is sufficient to prove Theorem 1 for k = 1. Then the transformed kernel is Here we assumed that the density of particles ρ 1 (t) = K(t, t) = ϕ (t)ψ(t) − ψ (t)ϕ(t) at the designated locus t is nonzero (otherwise the Jánossy density would vanish by definition).
The transformed kernel (20) is again of the Christoffel-Darboux form wherẽ Using Eq. (18), they are shown to satisfy a set of linear differential equations withÃ Since the coefficient functions m, A, B and C are polynomials in x, so are the new coefficient

Conditioning particles' loci as gauge transformation
Below we unravel the origin of inheritance of the TW criteria (17) and (18) and Eq. (22) is an SL(2, R) gauge transformation, Then the gauge-transformed sectionΨ(x) must be covariantly constant for the gauge-transformed sl(2, R) connectionÃ(x) that remains traceless, trÃ(x) = 0. Repetition of gauge transformations of the form on Ψ(x) yields the kth-order Jánossy density J k (x 1 , . . . , x k ; I). Although the gauge transformation U(x) has poles at x = x 1 , . . . , x k , the transformed sectionΨ(x) is regular and vanishes there. (iii) Meromorphy of A(x) inherits down toÃ(x) by Eq. (27) (which is equivalent to Eq. (24)), as U(x) is meromorphic.
Accordingly, the TW method is applicable to the evaluation of Jánossy densities of any continuous DPP if it is applicable to the evaluation of its gap probability, and J k,p (x 1 , . . . , x k ; I) is expressible in terms of a solution to a system of partial differential equations (PDEs) (containing x 1 , . . . , x k parametrically) in the endpoints {a i } of I.
A few comments are in order: • By construction (20), the transformed kernelK(x, y) vanishes when one of the arguments is equal to t,K (x, t) =K(t, y) = 0.
This leads to, for ∀ f ∈ L 2 (I), (K I · f )(t) = 0 and thus (( 8/18 • As mentioned above, the functionsφ(x) andψ(x) also vanish at x = t. It leads to q j (t) := ((I −K I ) −1 · x jφ )(t) = t jφ (t) = 0, p j (t) := ((I −K I ) −1 · x jψ )(t) = t jψ (t) = 0 (31) for j ∈ N. These could serve as part of the boundary conditions for the TW system, but we later use them only for a consistency check of the solution q j (s) and p j (s) derived from a different boundary condition imposed either at s 1 or s 1.
• For the sine kernel K(x, y) = sin( x − y governing the spectral bulk of unitary ensembles, Forrester and Odlyzko [16] previously considered the Fredholm determinant of its transformed kernel [19] (which they denoted as K 1 instead of ourK) with k = 1 and t set to 0 without loss of generality, They did apply the TW method to K 1 (x, y) and expressed the Fredholm determinant on a symmetric interval Det(I − K 1 | (−s,s) ) in terms of a solution to the TW system of ordinary differential equations (ODEs). However, in order to invoke the TW method they paid attention to an apparent fact that K 1 is related to K by a unit shift of the indices of the Bessel functions, rather than by an SL(2, R) gauge transformation U(x) = 1 0 −x −1 1 that retains the tracelessness of A(x) = 0 −1 1 0 . Nor did they explicitly write K 1 in a form of the right-hand side of Eq. (32), which would have meant K(x, y) − K(x, 0)K(0, 0) −1 K(0, y). Our formulation is a systematic generalization of the spirit of their work to arbitrary kernels of the TW type, to any interval I, and to any number (k ≥ 2) of conditioned particles.

Applications to random matrix theory
In this section we consider a DPP of eigenvalues {x i } of an N × N unitary-invariant random matrix ensemble with measure There the functions ϕ(x) and ψ(x) are (asymptotic forms of) the N th and (N − 1)th of the polynomials orthogonal with respect to the weight w(x). In this case, the conditional probability distributionsρ p (x 1 , . . . , x p ; t 1 , . . . , t k ) (8) andJ k,p (t 1 , . . . , t k ; I) (12) described by the transformed kernelK are nothing other than unconditional probability distributions of eigenvalues of a random matrix ensemble with weight functionw(x) = w(x) k j=1 (x − t j ) 2 , andφ(x) andψ(x) are (asymptotic forms of) polynomials orthogonal with respect tow(x). 1 If the values of the resolvent kernel R(x, y) = x|K I (I − K I ) −1 |y for arbitrary x, y ∈ I (not just its boundary values R(a i , a j ) at a i , a j ∈ ∂I, as derived in TW) were analytically available for various kernels appearing in the RMT, the Jánossy densities would readily be computed from the first line of Eq. (14). Since this path is infeasible (despite the fact that numerical evaluation of x|K I (I − K I ) −1 |y or Det(I −K I ) by the quadrature approximation is always possible [22]), we choose the second line of Eq. (14) as an alternate route.
As illustrative examples, we compute Jánossy densities for the Airy and Bessel kernels by applying the TW method to their transformed kernels.

Jánossy density for the Airy kernel
The Airy kernel governs local fluctuation and correlation of scaled eigenvalues of random Hermitian N × N matrices at the soft edge where the global density descends to zero as from which it follows that As an example we concentrate on the simplest of Jánossy densities, J 1 (t; I) with I = (s, ∞) and already set z to unity. Note that P 12 (t, s) = Θ(t − s)∂ s J 1 (t; (s, ∞)) represents the joint distribution of the first and second largest eigenvalues (t, s) of unitary ensembles, previously derived in Ref. [13] via Ref. [12] using a much more elaborate method than this work. The coefficient functions (24), whose degrees are increased by two after the redefinition Equation (36) could be slightly simplified by using the relation a 2 − b 2 t = 1 but we refrain. By taking the right endpoint a 2 of I to +∞, all terms in the TW system that contain a 2 vanish because of the exponential decay of the Airy kernel. Then the quantities involved, R(s) = s|K I (I −K I ) −1 |s (abbreviated notation of R(s, s)), The exponential decay of the Airy kernel and thus the transformed kernel (logK(x, y) log K(x, y) − 2 3 x 3/2 for x 1 and y, t fixed) leads to the boundary conditions for s 1: The diagonal resolvent (Eq. (1.7b) of TW) and the Fredholm determinant of the transformed kernelK I are expressed in terms of the solution to the ODEs (39): For numerical evaluation of the solution, in practice we impose the boundary condition q 0 (Λ) = ϕ(Λ), etc., at a sufficiently large positive Λ (∼ 10). Since q 0 (s) and p 0 (s) are regular at s = t (they are actually zero by Eq. (31)), apparent "double poles" at s = t in the first two nonuniversal equations of Eq. (39) are guaranteed to be canceled by the double zeroes on the right-hand side. Nevertheless, this could potentially cause loss of numerical accuracy when solving the TW system of ODEs from s = Λ down to s < t, e.g., by the explicit Runge-Kutta method. We have verified that this apparent stiffness at s = t can be circumvented by adding 11/18 Fig. 4 The joint distribution of the first and second largest eigenvalues P 12 (t, s) (orange) and the two-point correlation function ρ 2 (t, s) for t > s (transparent blue) of random Hermitian matrices.
to t a tiny imaginary part of the order of O(10 −10 ). With appropriately chosen values of = m(t), the real parts of q 0 (s) and p 0 (s) are stable upon varying , and q 0 ( e(t)) and p 0 ( e(t)) vanish up to the accuracy of O( ) as they should. In Fig. 4 we display the joint distribution of the largest eigenvalue t and the second largest eigenvalue s obtained by this prescription. The two-point correlation function ρ 2 (t, s) = ρ 1 (t)ρ 1 (s) − K(t, s) 2 , which is composed of peaks of joint distributions P k (t, s) of the kth and th largest eigenvalues for t > s (k < ), is overlaid for comparison. We have checked that the Fredholm determinant Det(I −K (s,Λ) ) = exp − Λ s ds R(s ) obtained by the TW system is in perfect agreement with numerical values from the Nyström-type quadrature approximation [23,24] where (x 1 , . . . , x M ; w 1 , . . . , w M ) is the Gauss-Legendre quadrature of the interval (s, Λ). Specifically, for a cutoff value Λ = 10, relative deviations between Λ s ds R(s ) computed from the TW system (39), (40), (41) using Mathematica's NIntegrate (for the preparation of boundary values) and NDSolve (for solving coupled ODEs) with = 10 −12 and at quadruple WorkingPrecision, and − log Det(I −K (s,Λ) ) computed by the Nyström-type approximation (43) with quadrature order M = 200, are between 10 −11 and 10 −8 for a range of variables −7 ≤ s, t ≤ 5 (see Table 1 for t = −2 and s = −7, . . . , 5). The table of numerical data for the Fredholm determinant Det(I −K (s,∞) ) is attached as online supplementary material.

Jánossy density for the Bessel kernel
from which it follows that Again we concentrate on the simplest of Jánossy densities, J 1 (t; I) with I = (0, s), and already set z to unity. P 12 (t, s) = −Θ(t − s)∂ s J 1 (t; (s, ∞)) represents the joint distribution of the first and second smallest eigenvalues (t, s) of unitary ensembles, previously derived in Ref. [12] using more elaborate methods. The new coefficient functions in Eq. (24) are, after redefinition, 13/18 The TW system of ODEs again takes the form (39), with the first two nonuniversal equations replaced by and the next eight universal equations sign-flipped: Note that by setting the left endpoint a 1 of I to 0, all terms containing a 1 either vanish or decouple. Accordingly all quantities are treated as functions of the right endpoint s alone, and their parametric dependence on t is implicit. Boundary conditions for s 1 are: The diagonal resolvent and the Fredholm determinant of the transformed kernelK I are expressed in terms of the solution to the ODEs (47): R(s) = −∂ s log Det(I −K (0,s) ) = p 0 (s)q 0 (s) − q 0 (s)p 0 (s).
For numerical evaluation of the solution, in practice we impose the boundary condition q 0 (µ) =φ(µ), etc., at a sufficiently small positive µ ∼ 10 −10 . The apparent stiffness in the first two nonuniversal equations of Eq. (47) at s = t can be circumvented by adding to t a tiny imaginary part of the order of O(10 −10 ). With appropriately chosen values of = m(t), the real parts of q 0 (s) and p 0 (s) are stable upon varying , and at s = e(t) they vanish up to the accuracy of O( ). The joint distribution of the smallest eigenvalue t and the second smallest eigenvalue s of random positive-definite Hermitian matrices, obtained by this prescription for ν = 0 and 1, and the corresponding two-point correlation function ρ 2 (t, s) = ρ 1 (t)ρ 1 (s) − K(t, s) 2 for t < s are converted to those of the singular values of random complex matrices by the replacements t → t 2 , s → s 2 and P 12 → 4ts P 12 , ρ 2 → 4ts ρ 2 and are plotted in Fig. 5. Again we have confirmed that the Fredholm determinant Det(I −K (µ,s) ) = exp − s µ ds R(s ) obtained by the TW system is in perfect agreement with numerical values from the Nyström-type quadrature approximation. Specifically, for 14/18 a cutoff value µ = 10 −12 , relative deviations between s µ ds R(s ), computed from the TW system (47)-(50) using Mathematica with = 10 −10 and at quadruple WorkingPrecision, and − log Det(I −K (0,s) ) computed by the Nyström-type approximation with quadrature order M = 200, are between 10 −12 and 10 −9 for a range of (original) variables 0 ≤ s, t ≤ 81 (see Table 2 for ν = 0, t = 4 and s = 1, . . . , 13). The table of numerical data for the Fredholm determinant Det(I −K (0,s) ) (after replacements t → t 2 , s → s 2 ) is attached as online supplementary material.

Conclusion and perspectives
In this article we have shown that the TW method is applicable to the evaluation of Jánossy densities and joint eigenvalue distributions for a kernel K . = (ϕ(x)ψ(y) − ψ(x)ϕ(y))/(x − y) if it is applicable to the gap probability. Essential to the inheritance of the TW criteria from K to the transformed kernelK . = (φ(x)ψ(y) −ψ(x)φ(y))/(x − y) is the structure that the map between the component functions (ϕ, ψ) → (φ,ψ) is an SL(2, R) gauge transformation to a covariantly constant section of an R 2 -bundle with an sl(2, R) connection A(x) B(x) −C(x) −A(x) . Our formulation generalizes the spirit of Ref. [16], which computed a special case of Jánossy density for the sine kernel by the TW method. As the simplest examples we evaluated the joint distributions of the two extremal eigenvalues P 12 (t, s) for

Fig. 5
The joint distribution of the first and second largest singular values P 12 (t, s) (orange) and the two-point correlation function ρ 2 (t, s) for t > s (transparent blue) of random complex matrices with ν = 0 (left) and ν = 1 (right).