On Random Matrix Averages Involving Half-Integer Powers of GOE Characteristic Polynomials

Correlation functions involving products and ratios of half-integer powers of characteristic polynomials of random matrices from the Gaussian orthogonal ensemble (GOE) frequently arise in applications of random matrix theory (RMT) to physics of quantum chaotic systems, and beyond. We provide an explicit evaluation of the large-N\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N$$\end{document} limits of a few non-trivial objects of that sort within a variant of the supersymmetry formalism, and via a related but different method. As one of the applications we derive the distribution of an off-diagonal entry Kab\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K_{ab}$$\end{document} of the resolvent (or Wigner K\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K$$\end{document}-matrix) of GOE matrices which, among other things, is of relevance for experiments on chaotic wave scattering in electromagnetic resonators.


Introduction
The goal of the present article is to attract attention to the problem of systematic evaluation of the large-N asymptotics of random matrix averages of the form where μ Fi , i = 1, . . . , K and μ B j , j = 1, . . . , L are sets of complex parameters. The angular brackets here and henceforth denote the average over the ensemble of realsymmetric N × N matrices H with Gaussian entries characterised by the probability density P(H ) ∝ exp − N 4J 2 Tr H 2 and known as the Gaussian orthogonal ensemble (GOE). Note that the correlation functions involving products of square roots of the characteristic polynomials in the numerator can be always reduced to the above form by multiplying and dividing both the numerator and the denominator with the same corresponding factors. Although there are reasons to suspect that the correlation functions (1) may have a nice mathematical structure even for finite N , perhaps not unlike those determinantal or Pfaffian structures discovered in [1][2][3][4][5] for similar objects involving only integer powers (see also [6,7] for an alternative derivation) we were not yet able to reveal such structures beyond the simplest case K = 1, L = 1, see (12) below. Instead we are mainly concentrating on the large-N limit of a few simplest, yet nontrivial examples of the correlation function of the type (1). We start with considering correlation functions with two square roots in the denominator, and with one or two characteristic polynomials in the numerator, that is C 1,2 (μ F1 ; μ B1 , μ B2 ) and C 2,2 (μ F1 , μ F2 ; μ B1 , μ B2 ), and then treat a special case of the correlation function involving four square roots in the denominator, and two determinants in the numerator, that is C 2,4 in our notation. As it should be clear from the examples given below the most physically interesting (bulk) scaling regime in the large-N limit arises when all spectral parameters are close to some value E ∈ (−2J, 2J ) by a distance of the order of the mean spacing between neighbouring eigenvalues in the bulk, i.e. O(J/N ). Correspondingly we define the scaled version of the correlation function as and C (bulk) where the approximate equality sign above should be understood in the sense of extracting the leading asymptotic dependence on the parameters ω B and ω F when N → ∞. Our results for the above correlation functions are given in Eqs. (13) and (14) for C (bulk) 1,2 and in Eqs. (16) and (17) for C (bulk) 2,2 . In Eq. (22) we provide the result for a special limit (see Eq. (8)) of C (bulk) 2,4 . These objects are already rich enough to provide answers for quantities arising in applications of random matrices in the field of Quantum Chaos in closed and open (scattering) systems. We discuss such relations in much detail below.
Although our methods are specifically tailored for dealing with the GOE we expect our results in the bulk scaling limit to be universal and shared by a broad class of invariant measures on real symmetric matrices H [8] and by so-called Wigner ensembles of random real symmetric matrices with independent, identically distributed entries satisfying relevant moments conditions [9,10].

Motivations and Background
To explain the origin of interest in the correlation functions (1) we start with recalling that the phenomenon of Quantum Chaos attracted considerable theoretical and experimental interest for more than three decades and remains one of the areas where applications of Random Matrix Theory are most fruitful and successful [11]. The applications are based on the famous Bohigas-Giannoni-Schmit (BGS) [12] conjecture claiming that in appropriately chosen energy window sequences of highly excited discrete energy levels of generic quantum systems whose classical counterparts are chaotic are statistically indistinguishable from sequences of real eigenvalues of large random matrices of appropriate symmetry.
Although not yet fully rigorously proven, this conjecture has an overwhelming support in experimental, numerical and analytical work of the last decades [13]. Inspired by this analogy as well as by the fact of universality of many random matrix properties (i.e. insensitivity to the particular choice of the probability measure on the matrix space), see [9,10] and references therein, one of the common strategies for predicting universal observables of quantum chaotic systems has been expressing them in terms of resolvents of underlying Hamiltonians, then replacing the actual Hamiltonians by random matrices taken from analytically tractable (usually, Gaussian) ensembles of N × N random matrices. The characteristic functions of the probability densities of the observables under consideration can be frequently computed explicitly by appropriate ensemble averages. Note that the eigenvalues of the standard Gaussian Ensembles, Unitary (GUE, β = 2), Orthogonal (GOE, β = 1) or Symplectic (GSE, β = 4) are independent of the eigenvectors, with the matrix of N orthonormal eigenvectors being uniformly distributed over the Haar's measure of the Unitary U (N ), Orthogonal O(N ) or Symplectic Sp(2N ) group, correspondingly. To that end it is natural to evaluate the corresponding characteristic functions by performing first the ensemble average over the eigenvectors. For the β = 2 case the average can be frequently done exactly for any N by employing the so-called Itzykson-Zuber-Harish-Chandra [14,15] formula, which is not yet available for β = 1, 4 group averages. Nevertheless, one is able to perform the eigenvector averages in the limit N 1 by using a heuristic idea (going back to [16]) that the set of eigenvectors essentially behaves for N 1 as if their components were independent, identically distributed Gaussian variables with mean zero and variance 1/N . One can rigorously justify this procedure if only a number n N 1/2 of eigenvectors is involved in the set, see e.g. [17], but in general a rigorous justification of such a step requires some nontrivial estimates on the resolvents. The heuristic procedure is widely employed in Theoretical Physics for RMT applications to Quantum Chaos using the properties of the standard Gaussian integrals over complex or real variables. In this way the analysis of many distributions of practical interest is reduced to correlation functions of products and ratios involving integer (for β = 2, 4) or half-integer (for β = 1) powers of characteristic polynomials of random matrices. Similar averages arise if one is interested in statistics of the matrix elements of the resolvents computed in the basis of random Gaussian vectors, as it is frequently done in applications to scattering systems with Quantum Chaos, see e.g. the recent paper [18] for an example and further references. For those and other reasons averages of products and ratios of powers of characteristic polynomials of random matrices attracted much interest over the years. When only integer powers are involved in the average the corresponding theory was developed for β = 2 in [2][3][4] and extended to β = 1, 4 in [5]. The case of half-integer powers for β = 1 remains however outstanding, despite the fact that it is most relevant for an overwhelming majority of experiments in Quantum Chaos due to the preserved time-reversal invariance of the underlying Hamiltonians. Additional interest to this type of averages gives the fact that they are closely related to the problem of evaluating averages of quantities involving absolute values of characteristic polynomials due to the relation | det valid for matrices H with real eigenvalues. Such averages emerge, for example, when studying the statistics of the so-called "level curvatures" in quantum chaotic systems [19,20], see Eq. (5) below, as well as in the problem of counting the number of stationary points of random Gaussian surfaces, see [21,22].
To support the above picture we describe below explicitly a few examples of relations between the characteristic functions of the physical observables of interest in quantum chaotic systems which can be related to particular instances of the correlation function (1). The list is almost certainly not exhaustive (for example, when writing this article we have learned that the square roots of characteristic polynomials emerged very recently in [23]), but hopefully representative.
• LDoS distribution. One of the first examples of that sort which is worth mentioning is related to the statistics of the local density of states (LDoS) ρ(x; E, η) at a point x of a quantum system with energy levels broadening η due to a uniform absorption in the sample. Mathematically the LDoS is defined in terms of the diagonal matrix element of the resolvent as ρ(x; E, η) = 1 π Im x|(E − iη N − H ) −1 |x , and one is interested in understanding the statistics of the LDoS assuming a random matrix GOE Hamiltonian H of size N × N , with the parameter η being fixed when N → ∞. The Laplace transform for the probability density P(ρ) of the LDoS can be expressed in the large-N limit as [24] ∞ 0 e −sρ P(ρ) dρ = Evaluation of the above random matrix average (which in our notation is a particular case of C (bulk) 2,4 ) attempted in [24] resulted in a quite impractical 5-fold integral, and to this end remains an outstanding RMT problem. Note however that the density P(ρ) has been found via a different route avoiding (4) as a sum of two-fold integrals in [25,26]. • Probability distribution of "level curvatures". Consider a perturbation H = H + αV of the Hamiltonian H where α is a control parameter and V is a real symmetric matrix. "Level curvatures" are defined as second derivatives of the eigenvalues λ n (α) (interpreted as energy levels of a quantum-chaotic system) with respect to the external parameter α: C n = ∂ 2 λ n (α) Assuming the perturbation V to be taken as well from the GOE one can show that the probability density P E (c) = 1 ρ(E) N n=1 δ(c − C n )δ(E − λ n ) of the level curvatures for GOE matrices H with eigenvalues λ n and mean density of eigenvaluesρ(E) can be represented as [19,20] where the required random matrix average in the right-hand side was independently evaluated by several alternative methods in [19,20]. Note that heuristic arguments appealing to Gaussianity of GOE eigenvectors in the large-N limit suggest universality of the level curvature distribution for a "generic" choice of V , and a rigorous proof of this fact is under consideration [27]. • Statistics of S-matrix poles. Various questions related to the statistics of quantum chaotic resonances (poles of the scattering matrix in the complex energy plane [28]) in the regime of a weakly open scattering system can be related to evaluation of the averages where ω is considered as N -independent parameter. The first of these averages features in the statistics of resonance widths change under influence of a small perturbation of the Hamiltonian H → H + αV akin to that considered above for the level curvature case. Such change reflects the intrinsic non-orthogonality of the associated resonance eigenfunctions [29]. Another manifestation of the same non-orthogonality is the statistics of the so-called Petermann factor which again can be related to random matrix averages involving half-integer powers of characteristic polynomials, see [30]. The second average in (6) arose in a recent attempt of clarifying the statistics of resonance widths beyond the standard first-order perturbation theory, see [31]. Evaluating both averages featuring in (6) in a uniform way by a systematic procedure was one of our motivations behind writing the present paper. • Statistics of Wigner K -matrix. In the theory of quantum chaotic scattering the Wigner K -matrix is essentially defined as a certain projection of the resolvent of H . More precisely this is an M × M matrix with entries K ab = W T a (E − H ) −1 W b , with W a being an N -component vector of coupling amplitudes W ia between N energy levels of the closed system (modelled for a chaotic system by an N ×N random matrix Hamiltonian H ) and M scattering channels open at a given energy E of incoming waves. Note that the more standard M × M unitary S-matrix is related to K via a simple Cayley transform S = I −i K I +i K . In the random matrix approach one usually assumes for the amplitudes W ia either the model of fixed orthogonal channels with W T a W b = γ a δ ab [32] or independent Gaussian channels where the amplitudes are taken to be i.i.d. Gaussian variables with W T a W b = γ a δ ab [33]. The quantities K ab are of direct experimental relevance and can be measured in microwave experiments as they are related to the real part of the electromagnetic impedance [34,35]. For real E in the bulk of the spectrum the statistics of the diagonal entries K aa is long known to be given by the same Cauchy distribution for all β = 1, 2, 4, see e.g. [36,37], and very recently was actually shown to be very insensitive to spectral properties of H under rather general conditions [38]. Similarly, one can consider the probability density P(K ab ) of the individual off-diagonal entries K a =b for β = 1. For the model of Gaussian channels one arrives to the Fourier transformed P(K ab ) in the form: Note that the average featuring in the right-hand side does not follow from either C as a special case, but is rather a limiting case of the more general correlation as it can be seen from the following representation: To the best of our knowledge the probability density P(K ab ) for a = b (or its Fourier transform) was not yet given explicitly in the literature 1 and we will find it below for the center of the GOE spectrum, see Eq. (22). Note that it is expected that statistics of the K -matrix entries for a GOE Hamiltonian H is the same for the two choices of the coupling W as long as M stays finite for N → ∞.
As to the M × M matrix K as a whole, the probability density P(K ) for β = 1 and E in the bulk of the spectrum is expected to be given by a Cauchy-like expression: with E-dependent mean K and the width parameter λ. This distribution was conjectured in 1995 by P. Brouwer on the experience of working with H from the so-called Lorentzian ensemble, see [42]. A similar formula for invariant ensembles of complex Hermitian random matrices H ( i.e. β = 2) was proved rigorously very recently in [18], and in the same paper it was mentioned that for β = 1 and the case of random Gaussian coupling the following relation holds 2 : for negative x c and is zero otherwise. Although our attempts to verify Brouwer's conjecture for β = 1, M = 2 along these lines were not fully successful yet, we discuss partial results, see (24)-(26) below.
• A particular type of the correlation functions (1) was investigated in [43] where it has been shown that for any integer k > 0 and fixed real δ holds 3

The Results
• As it has been mentioned above, we were not yet able to reveal nice mathematical structures for (1) at finite values of the matrix size N beyond the simplest case K = 1, L = 1, where the methods outlined below yielded a determinantal structure which we give here for completeness: ] is a Hermite polynomial and the function may be associated with the Cauchy transforms of Hermite polynomials [2]. • The explicit forms for the "bulk" correlation functions C (3)) depend very essentially on the signs of ω B1 and ω B2 . In particular, if sgn ω B1 = sgn ω B2 the first correlation function is given by whereas for sgn ω B1 = − sgn ω B2 the same object takes instead the form with where we introduced ρ = 1 2π J 2 √ 4J 2 − E 2 for the mean eigenvalue density of large GOE matrices in the bulk of the spectrum and used the standard notation K m (z) for the modified Bessel (Macdonald) functions of second kind and index m. Note that the asymptotic expression (14) shows an interesting "parity effect": it behaves differently depending on whether N is even or odd for arbitrary large values of N . Similarly the second correlation function for sgn ω B1 = sgn ω B2 is given by Let us now discuss a few special cases motivated by applications mentioned above.
• The characteristic function of the "level curvatures", Eq. (5) can be represented as a special limit of C The Fourier transform of this result (for brevity we choose E = 0, J = 1) yields the curvature distribution, which coincides with the expression found in earlier works by alternative methods [19,20]. • The two averages featuring in Eq. (6) can be recovered as special cases from C (bulk) 2,2 and are for the choice J = 1 given by The above formulas have been already presented in [29,31], with derivation relegated to the present paper. We tested the validity of (21) by direct numerical simulations of GOE matrices of a moderate size, see Fig. 1. • For the characteristic function of an off-diagonal element K ab of the K -matrix, see Eq. (7), we choose to present the corresponding result only for the so-called "perfect coupling" case, i.e. E = 0 and γ a = γ b = 1, the case of general γ a = γ b following by a trivial rescaling. It is given by The ensuing distribution P(K ab ) is then consequently given by its Fourier transform, In the Appendix A we verify that this result is in complete agreement with Brouwer's conjecture claiming that K for the "perfect coupling" case is distributed as P(K ) ∝ det[1 + K 2 ] −(M+1)/2 . We also check these expressions against direct numerical simulations, see Fig. 2.
Assume that x 1 x 2 > 0 so that (−x 1 x 2 ) = 0 and the sign-factor is immaterial. The correlation function then takes the form of which simplifies even further to e −|x 1 |−|x 2 | 2J for the "perfect coupling" case E = 0, γ 1 = γ 2 = 1. In the opposite case x 1 x 2 < 0 on the other hand the correlation function takes the form which is again a special case of C (bulk) 2,4 . In the particular case γ 1 x 1 = −γ 2 x 2 ≡ γ x, the above expression assumes the same form as one needed for extracting the distribution of a single off-diagonal element K ab , see Eqs. (7) and (22). While a full proof that K is distributed according to the Cauchy distribution, Eq. (9), requires the knowledge of the above expression for arbitrary values of x 1 and x 2 , one can show that our partial results for γ 1 x 1 = −γ 2 x 2 ≡ γ x are indeed consistent with Eq. (9), see Appendix B.
• Finally we notice that an interesting special case of C (bulk) 1,2 is the average of the sign of the GOE characteristic polynomial given asymptotically by where A(E, N ) is defined in Eq. (15).

Evaluation of the Correlation Functions Eqs. (2) and (3)
At present the only systematic method for evaluating the ensemble averages C (bulk) 1,2 and C (bulk) 2,2 seems to be the so-called supersymmetric formalism, see [45] and references therein. Within RMT framework several variants of that method are by now well-developed and we will follow one of them proposed in [46]. We only outline the major steps of the procedure below referring the interested reader to the cited literature and leaving technical detail for [47].
To that end we start with replacing the square roots of determinants in the denominator by Gaussian integrals over N -component real vectors x i , and the determinants in the numerator by integrals over vectors ζ i whose N components are complex anticommuting (Grassmann) variables. In that way the correlation function C (bulk) 1,2 can be represented by and similarly for C where we have to introduce one more integration over a vector of N anticommuting components. Note that we have to introduce s i ≡ sgn ω Bi in order to render the integrals over the commuting variables convergent.
The ensemble average can now easily be performed and yields for C where we introduced the N × N matrix B = s 1 x 1 ⊗ x T 1 + s 2 x 2 ⊗ x T 2 as well as the 2 × 2 matrices A similar expression for C (bulk) 1,2 can be obtained from the above by replacing all terms containing ζ 2 with 0 so that Q F becomes a scalar in this case. At the next step we employ a Hubbard-Stratonovich transformation for the anticommuting variables only by exploiting the identity where Q F = q 11 q 12 q * 12 q 22 is a Hermitian 2 × 2 matrix of commuting variables for C we also need to bilinearise the term ζ T 1 ζ 2 ζ † 2 ζ * 1 which can be achieved by introducing an auxiliary Gaussian integral over a complex variable u, with u * standing for its conjugate: With the integrand being bilinear in the Grassmann vectors it is easy to perform the integration over the anticommuting variables explicitly. The resulting expression in both cases depends on the x-vectors only via the eigenvalues of the matrix Q B L. This allows us to follow the route explained in detail in [43,46] and to employ the identity from Appendix D of [48]: which helps one to replace the integration over n real vectors of dimension N by an integral over a positive definite real symmetric matrix Q B of dimension n × n, where in both our cases actually n = 2. In the first case this procedure leads us after a trivial rescaling of the integration variables to where for notational convenience we omitted the hats here and henceforth. Similarly, in the second case we arrive at Here we introduced the 2 × 2 matrices M B(F) = E1 2 + i N diag(ω B1(F1) , ω B2(F2) ) and used λ B1 and λ B2 for the real eigenvalues of the 2 × 2 non-selfadjoint matrix Q B L, see [43,46] for technical details.
Setting aside the issue of performing the integration over the matrix Q B for the time being, in the first case the procedure leaves us with a single q-integration whereas in the second case we have to deal with an integral over the 2 × 2 Hermitian matrix Q F which contains four independent variables, and in addition with integrals over the complex variable u. To simplify the integrand we then use that Q F can be diagonalised by a unitary transformation Q F = U diag(q F1 , q F2 )U † . The integration over the unitary group can then be performed using the Itzykson-Zuber-Harish-Chandra (IZHC) formula [14,15] which reduces the integration variables to the set q F1 , q F2 , u and u * . Next we note that by introducing a matrix R = q F1 u u * q F2 one can express the integrand in terms of R (note e.g. that det Q F − u * u = det R, tr Q 2 F +2u * u = tr R 2 etc.). This latter matrix is Hermitian as well, so can also be diagonalized by a unitary transformation R = U 2 diag(r 1 , r 2 )U † 2 . Although the group integral is not of the IZHC type in this case, it still can be performed explicitly. Following this procedure the correlation function simplifies to At the final step we aim at simplifying the integral over Q B , which in both cases is a 2 × 2 real symmetric positive definite matrix. As the integrands in (34) and (36) actually depend on the combination Q B L we change the integration from Q B to Q B L. Recall that the matrix L = diag(sgn ω B1 , sgn ω B2 ) reflects the signs of ω B1 and ω B2 and this fact will play now a crucial role. If ω B1 and ω B2 are of the same sign, L is proportional to the identity and hence Q B L is still positive definite real symmetric and can be diagonalized by an orthogonal transformation Q B L = ±O diag( p 1 , p 2 )O T . If, however, the signs are different (we may assume for definiteness ω B1 > 0 and ω B2 < 0), then the matrix Q B L will have an underlying hyperbolic symmetry and can be parametrised as [43,46] where p 1 , p 2 > 0 and θ ∈ (−∞, ∞). The only term in the integrands (34) and (36) which actually depends on θ is tr For the sgn ω B1 = sgn ω B2 case one obtains the same type of expression with p 2 → −p 2 and cosh(2θ) → cos(2θ). The θ -integration can be performed explicitly using where I 0 (x) and K 0 (x) stand for the modified Bessel function of the first and second kind, respectively. in this way we arrive at the final expression which is exact for arbitrary value of N , and The superscript +− is to remind us that the expression corresponds to the choice ω B1 > 0 and ω B2 < 0. The expression for equal signs can be obtained from the above by replacing p 1 → + sgn(ω B1 ) p 1 , p 2 → − sgn(ω B1 ) p 2 , p 1 + p 2 → |p 1 − p 2 | and K 0 → I 0 . So far our manipulations were exact and did not use any approximation. As was explained in the introduction we are mainly interested in extracting the "bulk" large-N asymptotic of these correlation functions. The most natural way to proceed from here is by performing a saddle-point analysis. We believe with due effort such analysis can be done with full mathematical rigor, see e.g. a recent paper [49], but we do not attempt it here concentrating on explaining the gross structures of the saddle-point analysis which yield the correct results.
For the case of different signs the saddle points of the integrand are given by For p 1 and p 2 only solutions with positive real parts are contributing to the asymptotics.
There is no such restriction for q or r 1 and r 2 , respectively, and we have two saddle points contributing in each of these variables. Hence for C there are in principle four different contributions. However, the contributions from the saddle points satisfying r S P 1 = r S P 2 are actually negligible due to the factor r 1 − r 2 in the integrand. Moreover the integrand is invariant under exchanging r 1 and r 2 , and hence the two remaining contributions are identical. It therefore suffices to choose for r S P 1 the solution with positive real part and for r S P 2 the one with negative real part. One may further notice that the integrand itself vanishes when evaluated at the saddle points due to the factors (q − p 1 )(q + p 2 ) and (r 1 + p 1 )(r 2 + p 1 )(r 1 − p 2 )(r 2 − p 2 ) . This fact makes it necessary to expand the integrand to a higher order around the saddle points. The corresponding calculation is rather tedious, but managable. We refrain from presenting it here and refer the interested reader to [47] for technical detail. The outcome of the analysis are precisely the formulae given in Eqs. (14) and (17). The case of same signs looks quite different. Here the saddle points are given by where s ≡ sgn ω B1 = sgn ω B2 . Again we must choose p S P 1 and p S P 2 to have a positive real part, so that two contributions arise for C the same arguments as before suggest to choose for r S P 1 the solution with positive real part and for r S P 2 with negative real part neglecting the other three contributions. While the integrand still vanishes at he saddle points due to the factor | p 1 − p 2 | and for C (bulk) 2,2 due to the factors (r 1 + sp 1 )(r 2 + sp 1 )(r 1 + sp 2 )(r 2 + sp 2 ), the saddle point analysis is now much simpler than in the previous case. Indeed, when extracting the leading-order contribution one has to replace p 1 = p S P 1 + ξ 1 (with ξ 1 parametrizing the integration around the relevant saddle point) and similarly for the other variables, and then expand the N -independent part of the integrand to zero-th order in ξ 1 etc. (apart from the factors which come naturally in first order like | p 1 − p 2 | = |ξ 1 − ξ 2 |). It is then readily seen that the corresponding integrals yield a nonvanishing contribution rather straightforwardly without need to expand the integrand to higher orders like it was necessary in the previous case of opposite signs. The results of such saddle-point analysis is then much simpler and is given in Eqs. (13) and (16).

Distribution of K ab via Eq.(8)
For the correlation function (8) associated with the distribution of an individual off-diagonal K -matrix element we consider for simplicity only the perfect coupling case E = 0 and γ a = γ b = 1, see Eq. (22)). For evaluating the ensemble average we first tried to follow the same method as described in the previous section. In this way we started with writing det( and then replaced the square roots of characteristic polynomials in the denominator by four Gaussian integrals over real commuting vectors and those in the numerator by Gaussian integrals over two vectors with anticommuting components. The ensemble averaging then yields a 4 × 4 Q B -matrix, but we found no efficient ways of evaluating the ensuing group integral over the diagonalizing matrices. We also attempted a direct saddle-point analysis for large N along the same lines as before, and found it to become very tedious as not only the zero-th and first, but also the second order of the integrand expansion in fluctuations around the relevant saddle points turned out to be vanishing at the saddle points. Expanding to an even higher order with the group integrals still present did not seem to us as a viable option. Confronted with those difficulties we followed a different method (inspired by the insights from [30]) which avoids introducing anticommuting variables altogether. We demonstrate it first for the correlation function C (bulk) 1,2 . For brevity we will consider only the simplest case E = 0 where such object can be written as We start with representing only the denominator by a Gaussian integral over a real N -component vector S and hence get where (46) Note that the above integral is well-defined only for ω B1 and ω B2 having different signs, otherwise the term ω B1 ω B2 /N 2 > 0 would render the integral divergent.
Let us further assume that ω B1 = −ω B2 ≡ ω B , such that the linear term −i H ω B1 +ω B2 N vanishes. Such assumption is not necessary to make the method functional but helps to simplify the presentation considerably. Next we parametrize the vector S of integration variables as where h is a real N − 1-component vector, H N −1 is the (N − 1) × (N − 1) subblock of H and H 11 is the first element of H . With such a decomposition one is able to integrate out the variable H 11 as well as the vector h, which leads to where we have introduced the short-hand notations where the ensemble average should be performed over the (N − 1) × (N − 1) GOE matrix H N −1 . Moreover, it actually suffices to know only I 1 since I 2 = −i N d I 1 dω F . As is well-known I 1 is proportional to the Hermite polynomial: , so that asymptotically we have I 1 ∝ e ω F /J +(−1) N e −ω F /J . It remains to perform the S-integration for which it is advantageous to introduce rescaled polar coordinates, such that |S| 2 = N 2 R. The problem then reduces to performing the single integral For large N 1 it is easy to verify that the leading contribution to the integral can be written as which indeed coincides with the earlier derived expression for C (bulk) (14). Now we follow the same route for evaluation of the correlation function (8). We will only outline the key steps and differences from the previous case, but refrain from presenting intermediate results relegating them to [47]. One starts with replacing only the square roots of the characteristic polynomials in the denominator by Gaussian integrals, which leads us to where S 1 and S 2 are two real N -component vectors, and In contrast to a single vector S in the previous case we now have to deal with two real vectors S 1 and S 2 , which we can conveniently combine into the matrix Q. Such a ranktwo N × N matrix has two nonzero eigenvalues which we call q 1 and q 2 , all other N − 2 eigenvalues being identically zero. Being real symmetric Q can be diagonalised by an orthogonal transformation: Q = O diag(q 1 , q 2 , 0, . . . , 0)O T and the orthogonal matrices can be omitted from the integrand by the same invariance reasons as before. Owing to this structure we can conveniently decompose H into its upper left 2 × 2 block, its lower right (N −2)×(N −2) block H N −2 and the two ensuing off-diagonal blocks. It is easy to integrate out all variables apart from those entering H N −2 and get, with a slight abuse of notations: where we used the notations The result then reduces to performing ensemble averages over expressions det H 2 N −2 multiplied with various powers of traces of the inverse matrices H −k N −2 for a few instances of positive integers k. One may notice that all the required averages can be represented as derivatives of the correlation function of two GOE characteristic polynomials, using e.g. the identities and similarly for the higher powers. As a result for the object featuring in (53) we have: where the differential operator D ξ 1 ,ξ 2 (q 1 , q 2 ) is explicitly given by The ensemble average of the product of two GOE characteristic polynomials is known and for large N is given asymptotically by (see e.g. [50]) Using this result, and taking the necessary derivatives and the limits ξ 1 , ξ 2 → 0, we finally get an explicit expression for (q 1 , q 2 ). The last step is to perform the integrals over S 1 and S 2 , see Eq. (51). In the previous case we could reduce integration over S 1 to a single integration in polar coordinates. Similarly we can now exploit the invariance of the integrand and exploit the identity (33). In this way we can restrict the integration to the manifold of positive definite real symmetric 2 × 2 matrices with eigenvalues q 1 and q 2 . Extracting the leading large-N -asymptotics is then a straightforward exercise and we finally end up with the integral representation Note that here the limit → 0 is implied, which can now trivially be performed. It turns out that this rather complicated-looking integral is actually proportional to A way to verify this claim is to differentiate both equations (assuming for definiteness x > 0, J = 1) with respect to x. The derivative of Eq. (58) is x K 1 (x), and the derivative of (57) is x times a certain two-fold integral which with some efforts can be shown to be proportional to K 1 (x). The details of this calculation are relegated to [47].

Conclusions and Open Problems
In this paper we have started the program of systematic evaluation of correlation functions (1) involving half-integer powers of the characteristic polynomials of N × N GOE matrices. Motivated by diverse applications outlined in the introductory section we mainly concentrated on extracting the asymptotic behaviour of several objects of that type as N → ∞. Our calculations were based on variants of the supersymmetry method or related techniques. The method in a nutshell amounts to replacing the initial average involving the product of K characteristic polynomials divided by L square roots of characteristic polynomials of N × N GOE matrices H with an average over the sets of K × K matrices Q F and L × L matrices Q B > 0 with Gaussian weights augmented essentially with the factors det Q B and det Q F raised to powers of order N , see e.g. (35). As we are eventually mostly interested in K , L fixed but N → ∞ this replacement is very helpful as it allows to employ saddlepoint approximations. In this paper we managed to perform all steps of such a procedure successfully only for relatively small values of K and L, but we hope that the general case can eventually be treated along similar lines. One reason and guiding principle for a moderate optimism is as follows. An inspection of a somewhat simpler example of β = 2 shows, see in particular [2], that the success of our method is deeply connected to the existence of the so-called duality relations for Gaussian ensembles, see [51] for a better understanding of such dualities. In particular, the Proposition 7 of the latter paper shows that one of such duality relations exists for general Gaussian β-ensembles with β > 0 for an object involving the ensemble average of the product of the corresponding characteristic polynomials raised to the power −β/2. For the GOE with β = 1 that object (see Proposition 2 in [51]) is exactly the particular case of (1) with K = 0 and arbitrary integer L which makes a contact to the present context; e.g. one can employ such a duality to reproduce the relation (11) in an alternative way. A deeper understanding of connections between the supersymmetric approach and the duality relations for Gaussian ensembles will certainly be helpful in dealing efficiently with asymptotics of (1) for arbitrary integer values K and L. The problem of revealing possible Pfaffian-determinant structures behind (1) for finite matrix size N remains at the moment completely outstanding. It may well be that the methods of [6,7] or relations to generalized hypergeometric functions noticed for some particular instances in [44] could be useful for clarifying that issue.
In order to obtain the probability density for K 12 we need to integrate out the other two variables. We start with integrating out the variable K 22 The integration on the right-hand side can be easily performed as with the last integral on the right yielding artanh a. In this way we arrive at the probability density for K 12 in the form

Appendix B: Consistency Between Eq. (26) and Brouwer's Conjecture
We show that the characteristic function of the probability density P(K ) in the case M = 2 given in Eq. (26) is fully consistent with the claim that P(K ) ∝ det[1 + K 2 ] −3/2 . For the particular choice γ 1 x 1 = −γ 2 x 2 ≡ γ x the expression Eq. (26) is equivalent to Eq. (22) (for brevity we choose γ = 1). Our task then amounts to demonstrating that where X can be chosen diagonal, X = diag(x, −x). Since K is symmetric we can diagonalise it by an orthogonal transformation, The integral over the angle yields a Bessel function, and can also be rewritten in the form 2π 0 dφ e i 2 x(k 1 −k 2 ) sin(2φ) . Now note that 1 2 (k 1 −k 2 ) sin(2φ) ≡ −K 12 , which allows to present Eq. (65) in the form This is precisely the Fourier transform of P(K 12 ), which due to Appendix A is proportional to x K 0 (x) + ∞ x dy K 0 (y). This shows the validity of the claim (64).