Large-$N$ expansion for the time-delay matrix of ballistic chaotic cavities

We consider the $1/N$-expansion of the moments of the proper delay times for a ballistic chaotic cavity supporting $N$ scattering channels. In the random matrix approach, these moments correspond to traces of negative powers of Wishart matrices. For systems with and without broken time reversal symmetry (Dyson indices $\beta=1$ and $\beta=2$) we obtain a recursion relation, which efficiently generates the coefficients of the $1/N$-expansion of the moments. The integrality of these coefficients and their possible diagrammatic interpretation is discussed.


A. Background
The Wigner-Smith 11,37,40 time-delay matrix Q plays a central role in the theory of quantum transport. 14,38 It is defined in terms of the N-channel scattering matrix S via the relation where E is the energy of the incoming particle. If S is unitary, it is easily seen that Q is Hermitian. The eigenvalues τ 1 , . . . , τ N of Q are called proper delay times. Apart from the scalar case (N = 1), the individual proper delay times τ k have no immediate physical meaning. Physically relevant quantities are instead unitarily invariant functions of Q such as powers of traces TrQ k = τ k 1 + · · · + τ k N . In fact, several measurable observables are entirely determined by these powers of traces. For instance, the Wigner delay time TrQ is a bona fide measure of the time spent by an incident particle in the scattering region. Higher powers turn out to play a role in AC electronic transport, 5,6,15 e.g., in the low frequency expansion (ω → 0) of the AC dimensionless conductance G(ω) = [−iωTrQ + (1/2)ω 2 TrQ 2 + · · ·].
Using the random matrix approach to ballistic chaotic scattering, with the assumption that the internal Hamiltonian of the cavity belongs to a unitarily invariant ensemble with a large number of bound states, Brouwer, Frahm and Beenakker 4 showed (in quite an ingenious way) that the eigenvalues of (NQ) −1 are distributed as those of matrices in the Laguerre ensemble. In other words, denoting the latter rates by λ 1 , . . . , λ N , their joint probability density is supported on R N + and proportional to a) Electronic address: fabiodeelan.cunden@bristol.ac.uk where β ∈ {1, 2, 4} according to the physical symmetries of the cavity (hereafter the delay times τ k are measured in units of the Heisenberg time τ H ). The above joint distribution is valid in the presence of completely transparent contacts between the cavity and the external world. The situation of perfect coupling is indeed the most relevant, because it is possible to reduce arbitrary non-ideal coupling to the ideal case using the procedure in Refs. 33 and 35. Our main result gives an asymptotic expansion as N → ∞ for the average moments (k ≥ 0) where the expectation E[·] is taken with respect to (2). In particular, we derive explicit recurrence relations which efficiently provide the coefficients of the expansion to all orders in 1/N. 27). For all three symmetry classes β ∈ {1, 2, 4}, the following asymptotic expansion holds: It is tempting to call (4) a "genus" expansion. For complex Gaussian Hermitian matrices (Gaussian unitary ensemble (GUE)), the large-N expansion of the moments indeed enumerates maps of given genus. 12,41 Although we cannot prove that our expansion is related to the enumeration of maps, we have strong evidence of an underlying enumeration problem also for the moments of the ensemble modelling the Wigner-Smith time-delay matrix. More precisely, in this paper we extend a previous conjecture 8 on the integrality of the large-N expansion coefficients for the cumulants of TrQ k beyond the leading order.
Clearly, τ (β) 0 = 1. It is known that the average Wigner delay time is τ (β) 1 = 1 for every N and β. There is no general result for τ (β) k for k > 1 and generic β. For β = 1 (systems with time reversal symmetry) a finite-N formula is available, but is too lengthy to be reported here, see Ref. 26,Eq. (27). For β = 2 the problem somehow simplifies. First, the expansion (4) contains only powers of N −2 (that is, τ (2) k,g = 0 if g is odd); remarkably, two explicit finite-N formulae for τ (2) k are available in the literature. We report them here.
These two formulae have been derived independently with two different methods. Formula (5b) is quite hard to extract large-N asymptotics from (but this is possible, in principle, by using methods similar to those discussed in Refs. 9 and 22). Formula (5a) behaves better than (5b); although the number of terms in the sum is unbounded for large N, it is a sum of positive terms and the only asymptotic analysis one needs is the complete asymptotic series for the ratio of two Gamma functions. Along these lines of reasoning, the first three terms τ (2) k,g (g = 0, 2, 4) of the large-N expansion of τ (2) k have been obtained in Ref. 27. We also recall that the leading order coefficients τ (β) k,0 are independent of β and given by the large Schröder numbers, 34 where 2 F 1 is the classical Gauss hypergeometric function, and whose generating function has been computed explicitly in Refs. 3 and 7. A systematic study of the 1/N-expansion for the time-delay matrix moments for β = 1 and β = 2 is our objective.

B. Main results
Our central result is a recursion relation for the coefficients τ (β) k,g of the asymptotic expansion (4). As discussed in Ref. 8, the Laguerre measure (2) defines a one-cut β-ensemble. In particular, as N → ∞, the one-point marginal of (2) concentrates on a single interval whose edges are 3 ± √ 8. In this paper, we denote the monic polynomial vanishing at the edges by y(z) = z 2 − 6z + 1. (This is the so-called spectral curve of the β-ensemble (2).) Theorem 1.3 (Finite-N recursion at β = 2). For every integer k ≥ 1, and every N ≥ 1, with τ (2) 0 = τ (2) 1 = 1. The recursion formula (7) is much more efficient than (5a) and (5b) to generate tables of the moments of Q. See Appendix B. From Theorem 1.3 one obtains the recursion for the large-N coefficients τ (2) k,g as a corollary.
Corollary 1.4 (Double-recursion for the large-N coefficients at β = 2). The coefficients τ (2) k,g satisfy the homogeneous linear recurrence for g ≥ 0 and k ≥ 1, with initial conditions In particular, all coefficients τ (2) k,g with odd g vanish identically.
Analogous results can be obtained for β = 1. In this case we found that τ (1) k (and therefore τ (1) k,g ) satisfy an inhomogeneous recursion. Theorem 1.5 (Finite-N recursion at β = 1). For every integer k ≥ 1, and every N ≥ 1, where the auxiliary sequence b k is the solution of Corollary 1.6 (Double-recursion for the large-N coefficients at β = 1). The coefficients τ (1) k,g satisfy the inhomogeneous linear recurrence where the auxiliary sequence b k,g is the solution of for g ≥ 0 and k ≥ 1. The initial conditions are Again, the recursive formulae (10a) and (10b) and (11a) and (11b) can efficiently generate tables of moments and their large-N expansion. See Appendix B.
The proof of Theorem 1.3 and Theorem 1.5 is given in Section III.

C. Generating functions
In this section, we derive explicit formulae for the generating functions of τ (β) k,g . Let us consider the formal power series Using the recursions of Corollary 1.4 and Corollary 1.6 it is possible to obtain a differential equation for ϕ (β) (z, ζ). For instance, for β = 2, the generating function ϕ (2) (z, ζ) satisfies but this third-order inhomogeneous differential equation (17) is not as tractable. (This was to be expected, since ϕ (β) (z, ζ) is only a formal power series.) To make some further progress in the problem we introduce the "partial" generating functions The series ϕ (β) (z, ζ) , F (β) g (z) and J (β) k (ζ) are of course related by Remark 1. Note that ϕ (β) (z, N −1 ) is the generating function of the finite-N moments τ (β) k and The partial generating functions F (β) g (z) are central objects in the perturbative semiclassical approach.
The partial generating functions have a remarkably simple (i.e., algebraic) structure. We discuss first the case β = 2.
and has the following functional form (g ≥ 1) where R g (z) is a polynomial of degree 2g − 2.
Corollary 1.8. The order k ≥ 2 generating function J (2) k (ζ) is a rational function of the form is a polynomial satisfying the three term recursion Outline of the proof. The proof of Corollary 1.8 goes as follows. First, one multiplies the recursion (8) by ζ g +1 and sums over g to obtain a recursion for J (2) k (ζ). Inserting the expression (23) in the obtained recursion, elementary steps provide (24). It remains to be proved that P k (ζ) is a polynomial. In fact P 0 (ζ) = P 1 (ζ) = 1 are polynomials. From (24) it also follows that if P k−1 (ζ) and P k−2 (ζ) are polynomials so is P k (ζ). This completes the proof. The proof of Corollary 1.7 is again a routine calculation (multiplication of (8) by z k and sum over k). However, in this case R g (z) satisfies a recursion relation too complicated to be reported here. Nevertheless, it is easy to see that R g +2 (z) is a polynomial if R g (z) is so, and to compute its degree. The details are omitted.
From the partial generating functions F (2) g and J (2) k , it is possible to extract estimates on τ (2) k,g as k → ∞ (respectively, g → ∞) with g (respectively, k) fixed. These asymptotic results are based on Darboux's method 10  Corollary 1.9. The following asymptotics hold: The constants A g and B k are given explicitly by It is possible to obtain similar recursions for generating functions in the orthogonal case β = 1. Below we write the relation satisfied by F (1) g (z) explicitly. In this case, the recurrence relation for τ (1) k,g is not homogeneous and involves the auxiliary sequence b k,g . Therefore, the recursion for F (1) g (z) is coupled to a (homogeneous) recursion for the auxiliary generating function This fact complicates the structure of the formulae; it turns out that the generating functions F (1) g (z) are algebraic functions but they do not have the simple functional form (22) of . From a purely algorithmic point of view this is not a problem.
 , where the functions f g (z) satisfy The generating functions can be computed systematically and efficiently on standard computer algebra packages. The first few functions F (β) g (z) for β = 1 and β = 2 are reported in Appendix B. We note that only the leading order (planar) g = 0 and first four corrections g = 1, . . . , 4 have appeared in the literature so far. (See Ref. 27 for the random matrix approach and Refs. 2, 3, and 23 for semiclassical techniques.)

II. INTEGRALITY CONJECTURE AND ITS HEURISTIC EXPLANATION
The inspection of the first values of τ (β) k,g for β = 1 and β = 2 (see Table I) suggests that they are positive integers. A similar fact has been recently observed 7,8,28 for the leading order in 1/N of higher order cumulants of TrQ k (covariance, third order cumulants, etc.). Conjecture 1. For β = 1 and 2, τ (β) k,g ∈ N for every k and g. This conjecture extends beyond the leading order a generic conjectural statement 8 for the cumulants of TrQ k at generic β. We have considerable evidence supporting the conjecture. In fact, for β = 2, using the functional form of the generating functions we can actually prove that infinitely many τ (2) k,g 's are positive integers. We proceed to prove the following: Proof. The proof is based on the partial generating functions J (2) k (ζ) and F (2) g (z). For a fixed k, in order to prove that τ (2) k,g ∈ N for all g, we consider the generating function J k (ζ). The representation (23) shows that a sufficient condition for τ (2) k,g to be non-negative integers is that the polynomial P k (ζ) has non-negative integer coefficients (the series expansion of  j 1 − j 2 ζ 2 −1 at ζ = 0 is a product of geometric series), as suggested by the inspection of the first few polynomials (see Appendix B). Therefore, the claim "τ (2) k,g ∈ N for all g" involving an infinite number of coefficients can be proved by exhaustion of a finite number of cases: first, one computes the polynomial P k (ζ) using the recursion (24); then, one checks that the finitely many coefficients of P k (ζ) are all non-negative integers. This can be easily done by using a symbolic algebra software. (We have run a Maple code to compute recursively P k (ζ) and verify that it has non-negative integer coefficients for all k ≤ k ⋆ .) The proof that if g ≤ g ⋆ then τ (2) k,g ∈ N for all k goes along similar lines. For g odd the proof is trivial. Let us consider, for g even, the partial generating functions Here we use the classical identity where p ℓ (t) is the Legendre polynomial of degree ℓ. For t = 3, Hence we can focus on the coefficients of the polynomial R g (z). This case, however, is complicated by the fact that, although these coefficients seem to be integers, they are not necessarily positive (see Appendix B). Nevertheless, we can take advantage of (31) and (32) as follows. We have Note that C ℓ+1 ≥ C ℓ . If we denote R g (z) =  2g −2 j=0 a g, j z j , then Since C ℓ is nondecreasing we have Therefore, we conclude that the two conditions (i) a g, j ∈ Z and (ii)  j a g, j ≥ 0 imply τ (2) k,g ∈ N. Again, these conditions can be verified case by case using symbolic algebra softwares.
Obviously, the values k ⋆ and g ⋆ in Theorem 2.1 are fixed by limited computational power.

A. Semiclassical explanation of the conjecture
Periodic orbit theory is a collection of diversified results in the semiclassical analysis of quantum systems. Since its creation, 16,18,36 the theory has played an important role in the mathematical investigations of quantum chaos. Later, the ideas of periodic orbit theory were adapted to study quantum transport in the chaotic regime. Not surprisingly, a scattering orbit approach to delay times has also been developed. 2,3,23,31 The semiclassical approach is formulated in terms of the classical trajectories connecting the exterior and interior regions of the cavity. Each observable (e.g., τ (β) k ) is written as a sum over classical trajectories of wave amplitudes. The semiclassical contribution of a trajectory is determined by its topological properties. Therefore, the set of trajectories is partitioned according to topological properties where each class of trajectories is represented by a diagram with a given number of incoming channels, links, and encounters. See the recent paper by Kuipers et al. 23 for details. It turns out that the semiclassical contribution of a class of trajectories represented by a diagram D is given by (−1) c 1 (D) N c 2 (D) , where c 1,2 (D) ∈ Z. Then, one should sum over all classes of admissible trajectories (a sum over diagrams). It turns out that the admissible diagrams depend on the presence or absence of time-reversal symmetry (the Dyson index β in random matrix theory). The sum over diagrams is usually the hard part of the semiclassical approach; however, since the contribution of each diagram is given by the number of scattering channels N to some integer power, it is clear that the coefficients in the 1/N-expansion of τ (β) k are integers. We also mention that Novaes 29 computed the leading order of τ (β) k by considering the asymptotics of Selberg-like integrals; his method consists in enumerating certain classes of lattice paths. As expected, those paths were found to be in bijection with Schröder paths.

III. INVERSE MOMENTS OF WISHART-LAGUERRE MATRICES
In this section, we present several new results on the moments of inverse Wishart matrices. According to (2) and (3), the moments of the Wigner-Smith time-delay matrix Q are related to the inverse moments of a set of random variables belonging to a specific Laguerre ensemble (see Remark 2).

A. Technical results
Let W N be a N × N random matrix distributed according to the Wishart-Laguerre density where α is a generic parameter satisfying ℜ(α) > 0. This density is defined on the space of N × N positive definite real symmetric and complex Hermitian matrices for β = 1 and β = 2, respectively. We are interested in computing the moments (k ∈ Z) and their large-N expansions where from now on E[·] denotes averaging with respect to (36), and α is large enough to ensure that (37) is finite. Let us define the generating function where in what follows s ≤ 0. For simplicity, the dependence of (38) on α is omitted. First we present a lemma.
The next ingredient is that M (β) N (s) satisfies a differential equation. N satisfies the following homogeneous second-order differential equation: Theorem 3.3 (Generating function for β = 1). For all N ≥ 1, M (1) N satisfies the following inhomogeneous second-order differential equation: From these differential equations, it is easy to get a recurrence relation on the moments D (β) N (k, α) using Lemma 3.1. Indeed, for k < 0, multiplying (40) and (41) by s k and integrating from −∞ to 0, one removes all mention of M (β)′′ N (s) and M (β)′ N (s) by successive integrations by parts and then applies identity (39) to compute all integrals. Similarly, for k > 0, one differentiates (40) and (41) k times and then applies (39). Several cancellations simplify the outcome considerably. The final results are the following recurrence relations. N (k, α) satisfy Theorem 3.5 (Moments of real Wishart matrices, β = 1). The moments D (1) N (k, α) satisfy  (40) and then deduced the finite-N recurrence (42) for positive moments k > 0, thus generalizing the Harer-Zagier recursion formula for GUE matrices 19 to the complex Wishart ensemble. Here we show that the same recursion holds for negative moments (k ≤ 0). The differential equation (41) and the recursion (43) for β = 1 are new results. The inhomogeneous term in the differential equation (41) for β = 1 is related to a β = 2 ensemble; the reason for this will become clear in the proof below (see Eq. (47)-(51)). Theorem 1.3 is a specialization of (42) renaming τ (2) . Theorem 1.5 is a specialization of (43) renaming τ (1) In both real and complex cases, traces of powers of Wishart matrices were considered for generic covariance matrix Σ using orthogonal and unitary Weingarten functions. 24,25 However, even in the simplest case Σ = I (the one considered here), the Weingarten functions are not easy to compute. Furthermore, it is not clear how to obtain an asymptotic 1/N-expansion from such formulae.
Proof of Lemma 3.1. For k > 0, Eq. (39) is a classical formula. Let us consider the case of negative moments. By a unitary transformation we can write W N = UΛU † where Λ is diagonal, so that W N e sW N = UΛe sΛ U † and the left-hand side of (39) is where the integral acts on each diagonal entry via Then the trace is the sum over j and one obtains the right-hand side of (39).

B. Proof of the technical results
Proof of Theorems 3.2 and 3.3. We introduce the finite-N average density of eigenvalues of W N where the expectation is taken with respect to P β,α in (36). The generating function M (β) N (s) can be written as The fundamental ingredient is that for β = 2 and β = 1 the finite-N eigenvalue density ρ (β) N (x) is given explicitly in terms of Laguerre polynomials. We denote the standard Laguerre polynomial of degree N and parameter α by They satisfy the second order differential equation The following formulae for the mean eigenvalue densities are well-known. For β = 2, the Christoffel-Darboux formula leads to the expression For β = 1, explicit forms for the mean eigenvalue density for the Laguerre ensembles were given in Refs. 13 and 39. The result is that ρ (1) N (x) can be expressed in terms of ρ (2) N −1 (x) plus a correction term. The correction term was understood in a more general context in Ref. 1, where it was formulated in a slightly different way that turns out to be useful here. Moreover, the correction term depends on the parity of N. Nevertheless, the moments τ (1) k are always rational functions of N, and hence they are uniquely determined by subsequences like even N's. For simplicity, in the case β = 1 we will perform our computations for N even, but our final results do not depend on the parity of N. From Ref. 1, for N even we have where the constant d N is given by and The plan is now to insert these expressions into (47) and derive differential equations for the generating functions. β = 2: Derivation of equation (40) Equation (40) has been derived in at least two ways in the literature, first by Haagerup and Thorbjørnsen 17 and then rediscovered in a more general setting by Ledoux. 20 We will sketch below the proof given in Ref. 17, where the idea is to prove that It satisfies a classical second-order ordinary differential equation, which after some lengthy algebraic manipulations yields (40). To prove (54), a crucial role is played by the classical second-order differential equation (49) satisfied by the Laguerre polynomials. Combining the differential equation with the expression (50) leads to d dx Then integration by parts in (47) results in To compute the integral in (56), substitute u = x(1 − s) and make use of the scaling identity with c = (1 − s) −1 , which reduces the integral to the orthogonality relation of Laguerre polynomials. The single summation that remains is recognised as the series definition of the hypergeometric function, and hence (54). β = 1: Derivation of equation (41) The proof of Theorem 3.3 for β = 1 is given below and is based on Ref. 21 where the analogous computation was done for the Gaussian orthogonal ensemble. It was suggested in Ref. 21 that the computation done there could in principle be carried out for other classical ensembles, but to our knowledge this has not been done before.
The starting point is the finite-N formula (51) which we insert into (47), where From Theorem 3.2, we may treat M (2) N −1 (s) in (58) as a known quantity. Hence we seek a differential equation for the second addendum in (58). We write the differential equation (49) as the eigenfunction relation where Next, for any sufficiently smooth f and g, the following identity is a direct consequence of integration by parts: We now derive a differential equation for Y (s). Taking a derivative with respect to s and integrating by parts yields In the last integral we used that ψ ′ (x) = 2x (α−1)/2 e −x/2 L (α) N −2 (x). We can rewrite this as where Note that u N −1 (s) is closely related to the Laplace transform of ρ (2) N −1 (x) (cf. identity (56)), Differentiating (63) one more time we arrive at On the other hand, we can use (60) and (61) to show that Differentiating (68) gives the identity where Solving (63) and (67) for χ N (s) and χ ′ N (s) and inserting the result into (69) gives Our aim is now to express K N in terms of known quantities. Integrating by parts in (70), we find Adding the two representations (70) and (72) shows that where The difference of Laguerre polynomials in (76) is nothing but the Christoffel-Darboux form of the eigenvalue density ρ (2) N −1 (x), cf. formula (50). We deduce that and hence K N can be expressed in terms of explicitly known quantities. Using (66) and solving (74) and (75) for 2K N , we insert the results into (71) and obtain the closed equation All that remains is to rewrite this in terms of M (1) N (s) and its derivatives using (58) The result (41) now follows after using (40) of Theorem 3.2 to eliminate the variable M (2)′′ N −1 (s) appearing in (79).

IV. CONCLUSIONS AND REMARKS
We considered the average of power traces E TrQ k of the time-delay matrix for ballistic chaotic cavities. The sample-to-sample average can be computed using a random matrix theory ansatz. The large-N expansion of these averages has been computed (recursively) for systems with and without broken time reversal symmetry ( β = 2 and β = 1, respectively); we suggest that the coefficients of the large-N expansion are whole numbers (Conjecture 1), thus extending a previous conjecture for the leading order of higher cumulants. A heuristic explanation of the conjecture comes from the semiclassical approach to chaotic scattering.
We conclude with a last remark. We recall that, up to a scaling factor, the moments of the time-delay matrix are statistically identical to the inverse moments Tr(W −k N ) of a specific Wishart ensemble. More precisely, for β = 2 for concreteness, W N = X X † where X = (x i j ) is a N × 2N matrix with independent standard complex Gaussian entries. It is known that the large-N expansions of moments of Gaussian (GUE) matrices are integers and they are related to a very precise enumeration problem. Similar enumeration problems emerge for the moments of Wishart matrices. The salient feature of Gaussian and Wishart matrices, and the one which is fundamental for applications to graphical enumeration, is that expectations of general polynomial functions can be reduced to a counting of Wick pairings. When considering the average of traces of powers of inverse Wishart matrices (our situation), Wick's calculus no longer applies, but the integer nature of the coefficients nevertheless suggests an underlying combinatorial structure yet to be unveiled. The authors are grateful to a referee for useful comments. All underlying data are included in full within this paper.