text:
# Calculations of radio pulses from High Energy Showers

## I INTRODUCTION

It was about 40 years ago that Askar’yan proposed the detection of very high energy particles through the coherent radio emission from all the particles in a high energy shower askaryan . Although radio emission from showers had been discussed before, Askar’yan noted that showers would develop an excess of electrons of order $`10\%`$ independently of shower energy and that such excess would lead to coherent radiation. Indeed if the emission is coherent the electric field scales as the excess charge and as the shower energy rises, the relative contribution of radio emission increases with respect to other wavelengths (say optical Čerenkov). Since the number of particles in a shower is proportional to shower energy, the power in the coherent signal scales with the square of the primary energy. This fact together with the technical simplicity of detectors makes the technique attractive for detecting very high energy showers. It is remarkable that many detection possibilities addressed in this meeting were already in Askar’yan original work.

Clearly coherence demands that the wavelength of the radiation exceeds the size of the radiating region. As electromagnetic showers develop through successive pair production and bremsstrahlung interactions, the particle distributions of showers of a given energy are, to a very rough approximation, similar when described in radiation length units 111Hadronic showers can be regarded as a jet of hadrons along shower axis which generates a chain of electromagnetic subshowers through neutral pion decays and have similar scales.. The radiation length of different materials in units of matter depth basically scales with the inverse of the atomic number, being $`37`$ g cm-2 for air and not too different for many other abundant materials that have been proposed for radio experiments such as water, salt and sand. Because of the medium density the scale of the showers can however be completely different. While air showers have sizes in the km range and coherence effects can be obtained at frequencies of order 10 MHz, for showers in ice the scale is a meter and coherence is kept to frequencies of order 10 GHz.

Since the power associated to the Čerenkov emission increases with frequency Askar’yan suggested detecting showers in dense medium, as those that could be induced by deeply penetrating particles underground or by cosmic rays on the Moon surface. Attempts were made to detect radio pulses in coincidence with other air shower detectors weekes  and, although the results were positive, it became clear that there were many difficulties. Calculations revealed that several additional mechanisms competed with excess charge in the generation of radio pulses, including dipole radiation and transverse currents because of the magnetic field of the Earth, and transition radiation when the shower reached the ground kahn . The interpretation of radio signals was obscured by these and other difficulties and the technique was almost completely abandoned but for a few isolated efforts that have been reviewed in this conference.

In the mid 1980’s the efforts to build high energy neutrino telescopes renewed the activity in the radio technique markov ; ralston  as a possible alternative to the projects to detect the Čerenkov light from muons under water or ice halzen . As a result early in the 1990’s a Monte Carlo program was developed for the calculation of distant radio signals from TeV electromagnetic showers in ice zhs91 ; zhs92 . These led to several efforts to assess the possibility to detect neutrino events with arrays of antennas in Antarctica provorov ; frichter ; jelley . Later on a series of approximations were made to extend the results of these calculations to both electromagnetic and hadronic showers of energies in the EeV and ZeV ranges alz97 ; alz98 ; alz99 ; alz00 . In this article we will concentrate on the efforts made towards understanding Čerenkov radio emission in dense media from the excess charge that develops in high energy showers. The recent experimental confirmation of this effect in a SLAC experiment gorham , in agreement with expectations, leaves now little doubt about the enormous potential of this technique. We believe that many developments will follow this conference and it is our hope that these calculational efforts will help in the task.

## II ZHS: A shower generator for radio emission.

The calculation of radio pulses from showers is a complex problem in which the emission from all particles in a high energy shower have to be carefully added in a coherent fashion, taking into account the velocity and direction, position and time of each of the particles involved. For a medium such as ice, having a refraction index $`n=1.78`$ for radio waves, electrons of energy exceeding 107 keV emit Čerenkov radiation. The simulation of the shower cascade from this perspective was tackled by Zas, Halzen and Stanev in zhs92  with a result that has become known as the ZHS Monte Carlo. Since the largest contribution to the excess charge in a shower is from electrons, only electromagnetic showers were considered. This program has been a reference point for much progress made during the 1990’s in radio detection and deserves a brief description. Efforts have been made by other groups to test the code against other shower codes besson .

The longitudinal development of an electromagnetic shower mostly follows from successive pair production and bremsstrahlung interactions with the electric field of atomic nuclei that are the heart of the showering process. For both these processes the effective interaction distance increases as the energy of the interacting photon or electron rises and atomic screening of the atoms becomes progressively more important. In this respect the Landau-Pomeranchuck-Migdal (LPM) effect can be regarded as the higher energy extreme of this situation in which the atomic potential of the individual atoms has to be modified because of collective effects of several nuclei. The lateral structure of the shower is mostly controlled by elastic scattering. Because of the long range of the Coulomb interaction it turns out that most angular deviations arise because of multiple elastic scattering which is treated in this Monte Carlo keeping the first two terms in Molière’s expansion moliere .

Besides these interactions that are common to all simulation programs, the ZHS code has special features for radio calculations. It is a three dimensional code that takes timing of the particles in great detail. Time delays are measured with respect to a plane perpendicular to shower axis that starts in phase with the primary particle and moves parallel to the axis at speed $`c`$. Most of such delays are due to angular deviations, to multiple elastic scattering and to the subluminal velocities of the shower particles all of which are taken into account. The program includes electron positron annihilation, Möller, Bhabha and Compton scattering, which only become relevant at low energies but are responsible for the excess charge in the shower. These routines were developed following the EGS4 code egs4 . In the first version zhs91  it did not include the LPM effect but this was soon remedied zhs92 . It can be run with threshold energies as low as 100 keV, but its accuracy is expected to decrease for thresholds below this value because it does not include the photoelectric effect. Finally it is designed to be fast using, when possible, numerical parameterizations instead of actual functions to minimize calculating time.

The results obtained for the longitudinal development of the showers are consistent with the Greisen parameterization. The lateral distributions are consistent with three dimensional theoretical shower results provided the Molière radius is divided by a factor of about 2, in agreement with other simulations hillas . The excess charge obtained is of order $`20\%`$ depending on particle threshold, it is only slightly dependent on shower depth, see Fig. 1a. The extra track due to electrons is $`21\%`$ of the total (the Čerenkov emission is proportional to the tracklengths of the charged particles) as shown in Fig. 1b. Compton scattering of atomic electrons that get accelerated into the shower account for most of the effect ($`12\%`$). It is remarkable that close to half ($`14\%`$) of the excess tracklength is due to particles below 5 MeV (1 MeV). The calculation of coherent signals requires a lower threshold than implied by these numbers because the lower energy particles are expected to interfere more destructively as they get longer delays and larger angular deviations. This is already important because it sets an upper bound to the energy threshold necessary for an accurate calculation.

## III Čerenkov radiation

The detection of radio pulses from air showers is clearly a very different problem from detecting them in a dense medium such as ice. Firstly in ice the coherence band extends to the GHz regime. Secondly most shower particles travel very small paths of tens of centimeters, and they are hardly deviated because of the Earth’s magnetic field. As a result the excess charge mechanism is expected to dominate the radiation and the problem becomes simpler from the calculation point of view. Moreover there is a substantial difference in the relation between the typical observation distances and the shower dimensions. While for air showers the observation distance is of the some order or even less than the shower dimensions, prospects of detecting radio pulses under the Earth surface or from the Moon are attractive precisely because they can be detected from very large distances compared with the shower dimensions which are of order one meter. For most interesting cases the Fraunhofer approximation turns out to be accurate what is useful for calculational purposes and simplifies both the interpretation and the presentation of results.

When a particle of charge $`z`$ moves through a medium of refractive index $`n`$ with velocity $`|\stackrel{}{v}|=\beta c>c/n`$ Čerenkov light is emitted at the Čerenkov angle $`\theta _C`$, verifying $`\mathrm{cos}\theta _C=(\beta n)^1`$, with a power spectrum given by the well known Frank-Tamm result fran  :

$$\frac{d^2W}{d\nu dl}=\left[\frac{4\pi ^2\mathrm{}}{c}\alpha \right]z^2\nu \left[1\frac{1}{\beta ^2n^2}\right],$$
(1)
with $`\nu `$ the frequency, $`c`$ the speed of light, $`dl=c\beta dt`$ a small element of particle track length, and $`\alpha `$ the fine structure constant. This is the standard approximation used for most Čerenkov applications for wavelengths orders of magnitude smaller than the tracks. This expression corresponds to an infinite track and in the case of a finite track there are diffraction effects.

Working with the time Fourier transform of the electric field, $`\stackrel{}{E}(\omega )`$, Maxwell’s equations can be solved in the transverse gauge for an arbitrary current zhs92 :

$$\stackrel{}{E}(\omega ,\stackrel{}{\mathrm{x}})=\frac{\mathrm{e}\mu _\mathrm{r}}{2\pi ϵ_0\mathrm{c}^2}i\omega \mathrm{dt}^{}\mathrm{d}^3\stackrel{}{\mathrm{x}^{}}\mathrm{e}^{i\omega \mathrm{t}^{}+i|\stackrel{}{k}||\stackrel{}{\mathrm{x}}\stackrel{}{\mathrm{x}^{}}|}\frac{\stackrel{}{\mathrm{J}}_{}(\mathrm{t}^{},\stackrel{}{\mathrm{x}^{}})}{|\stackrel{}{\mathrm{x}}\stackrel{}{\mathrm{x}^{}}|}$$
(2)
where $`|\stackrel{}{k}|=k=n\omega /\mathrm{c}`$ and $`\stackrel{}{\mathrm{J}}_{}`$ is a divergenceless component of the current transverse to the observation direction, $`\mu _r`$ is the relative permeability of the medium and $`ϵ_0`$ the permittivity of the vacuum. Eq. 1 can be obtained from this expression. For a single particle moving uniformly with velocity $`\stackrel{}{v}`$ between times $`t_1`$ and $`t_1+\delta t`$ in the Fraunhofer limit it simplifies to:

$$\stackrel{}{E}(\omega ,\stackrel{}{\mathrm{x}})=\frac{e\mu _\mathrm{r}i\omega }{2\pi ϵ_0\mathrm{c}^2}\frac{\mathrm{e}^{ikR}}{R}\mathrm{e}^{i(\omega \stackrel{}{k}\stackrel{}{v})\mathrm{t}_1}\stackrel{}{v}_{}\left[\frac{\mathrm{e}^{i(\omega \stackrel{}{k}\stackrel{}{v})\delta \mathrm{t}}1}{i(\omega \stackrel{}{k}\stackrel{}{v})}\right]$$
(3)
The Fraunhofer limit corresponds to an observation distance $`R>>[v\delta tsin\theta ]^2/\lambda `$. For a characteristic track of one radiation length ($`40`$ cm) and $`\lambda `$ of 30 cm corresponding to 1 GHz frequencies, the condition implies $`R>>50`$ cm which is well satisfied for detection distances of order one km. It should be remarked that the inequality can be always satisfied provided a sufficiently small time interval is chosen. The ZHS program uses this expression. For an infinitesimal track $`\delta l`$ it becomes:

$$R\stackrel{}{E}(\omega ,\stackrel{}{\mathrm{x}})=\frac{e\mu _\mathrm{r}i\omega }{2\pi ϵ_0\mathrm{c}^2}\stackrel{}{\delta }l_{}\mathrm{e}^{i(\omega \stackrel{}{k}\stackrel{}{v}_1)\mathrm{t}_1}\mathrm{e}^{ikR},$$
(4)
This expression shows that the electric field amplitude is proportional to the tracklength and that radiation is polarized in the direction of $`\stackrel{}{l}_{}`$, the apparent direction of the track as seen by the observer.

For a high energy shower the calculation of the radio pulse is now a matter of book-keeping. Results of the electric field amplitude in the Čerenkov direction display a characteristic frequency spectrum: a linear rise with frequency until the GHz region where destructive interference sets in and the spectrum gradually levels off as shown in Fig. 2a. Although in the Čerenkov direction all points of the trajectory of a particle moving at the speed of light along the shower axis are emitting in phase, lateral deviations of the particle trajectories and time delays destroy the coherent interference at sufficiently high frequencies. The effect of the lateral distribution is dominant. The amplitude at a given frequency is proportional to the excess tracklength which in turn accurately scales with shower energy as shown in Fig. 1b.

The polar diagram of the emitted radiation displays a characteristic peak in the Čerenkov direction of a width which is inversely proportional to the frequency. This is plotted in Fig. 2b and corresponds to a diffraction-like pattern. The width of the diffraction peak about the Čerenkov angle depends on the longitudinal distribution of the charge in the shower. As the observation angle gets away from the Čerenkov direction different longitudinal sections of the shower development start to interfere destructively and the amplitude of the electric field drops.

Long tracks in a shower can be subdivided in arbitrary subintervals. Provided the Fraunhofer approximation is valid for the shower as a whole, the algorithm of Eq. 3 for a particle that moves at constant speed displays cancellations between all the end points of the intermediate subintervals of the track. Only the end points of the total track need to be considered in this case what simplifies the computational task enormously. In approximation “a” (default) each electron gives rise to a single track to which Eq. 3 is applied. As particles do suffer some deceleration along the track the average velocity is used. Different criteria for this subdivision have given rise to approximations “b” and “c” rome . The approximation “a” has been shown to provide sufficiently good results when compared to longer calculations in which each electron track is subdivided into subtracks each of which contributing with Eq. 3 as shown in Fig. 3. Note that approximation “a” is conservative in the sense that it underpredicts the amplitude of the pulse at very high frequencies.

## IV 1-D Approximation

With the help of scaling it is possible to raise the shower simulation threshold and to generate showers of about 1 PeV but at higher energies the calculations become too lengthy to be handled. The calculation of radioemission from EeV showers is however of great interest because the technique is expected to be most competitive for these energies. Because of the LPM suppression of bremstrahlung and pair production, these showers have very different elongated depth distributions. As a result the angular width of the Čerenkov peaks is expected to be reduced with respect to showers of energy below 1 PeV. As the lateral distributions of these showers are not expected to be very different from the lower energy showers, not many effects can be anticipated in the frequency spectrum of the signal in the Čerenkov direction.

The one dimensional approximation has been developed to study the effects of the longitudinal elongation alz97 ; alz98 . In this approximation the particles in the shower are all assumed to be moving at the speed of light and the lateral distribution is neglected. The current density simply becomes alz00 :

$$\stackrel{}{J}_{}(\stackrel{}{\mathrm{x}^{}},\mathrm{t}^{})=Q(z^{})\stackrel{}{c}_{}\delta ^3(\stackrel{}{\mathrm{x}^{}}\widehat{n}_zct^{})$$
(5)
where the shower develops along the $`z`$ axis, $`\widehat{n}_z`$ is a unit vector in this direction and $`Q(z^{})`$ is the excess charge in the shower. The substitution of this current into Eq. 2 leads to:

$$\stackrel{}{E}(\omega ,\stackrel{}{\mathrm{x}})=\frac{e\mu _\mathrm{r}}{2\pi ϵ_0\mathrm{c}^2}i\frac{\omega }{c}\mathrm{sin}\theta \widehat{n}_{}𝑑z^{}Q(z^{})\frac{e^{i\frac{\omega }{c}z^{}+ik|\stackrel{}{\mathrm{x}}z^{}\widehat{n}_z|}}{|\stackrel{}{\mathrm{x}}z^{}\widehat{n}_z|}$$
(6)
where $`\theta `$ is the angle between the shower axis and the direction of observation $`\stackrel{}{\mathrm{x}}`$ and $`\widehat{n}_{}`$ is a unitary vector perpendicular to $`\stackrel{}{\mathrm{x}}`$. It should be stressed that no further approximations are made at this stage.

The phase factor in Eq. 6 can be approximated by $`ik|\stackrel{}{x}\stackrel{}{z}{}_{}{}^{}|ikRi\stackrel{}{k}\dot{\stackrel{}{z}}^{}`$ in the Fraunhofer limit. Here $`R=|\stackrel{}{\mathrm{x}}|`$ is the distance from the center of the shower to the observation point. The electric field amplitude simply becomes the Fourier transform of the longitudinal charge distribution:

$$\stackrel{}{E}(\omega ,\stackrel{}{\mathrm{x}})=\frac{e\mu _\mathrm{r}}{2\pi ϵ_0\mathrm{c}^2}i\omega \mathrm{sin}\theta \frac{\mathrm{e}^{ikR}}{R}\widehat{n}_{}𝑑z^{}Q(z^{})\mathrm{e}^{ipz^{}}$$
(7)
where we have introduced for convenience the parameter $`p(\theta ,\omega )=(1n\mathrm{cos}\theta )\omega /c`$ in Eq. 7. The angular pattern around the Čerenkov direction is precisely the analog of the classical diffraction pattern of an aperture function.

The accuracy of the approximation has been studied comparing it to complete simulations. Excellent agreement is obtained for frequencies below 100 MHz. For frequencies progressively increasing above this value, the electric field becomes overestimated in the Čerenkov direction, but the angular distribution of the pulse is otherwise preserved. This is not surprising since the lateral distribution is ignored and no levelling off in the frequency spectrum can be expected. Moreover the difference between the complete simulation and the approximation scales with shower energy as shown in Fig. 7. Ad hoc corrections can be implemented to make quantitative predictions in the GHz regime accounting in an effective way for the ignored effects alz98 ; alz00 . This is analogous to the form factor described in buniy . The approximation gives some insight into the complexity of the spectrum and angular distribution of the pulses relating them to features in shower development.

The development of EeV showers of electromagnetic alz97  and hadronic alz98  character has been studied using hybrid techniques that combine simulation and parameterizations of low energy showers. The LPM suppression leads to a very large elongation of electromagnetic showers which results in a reduction of the angular width of the diffracted pulses. The elongation has been calculated to scale as $`E_{sh}^{1/3}`$ for shower energies satisfying $`E_{sh}>20`$PeV alz97 . For hadronic showers the elongation does not affect shower development significantly until energies in the EeV range. This is because high energy hadronic interactions have large multiplicities and the particles that emerge carry small fractions of the primary energy. Moreover the neutral pions at very high energies interact before decaying into photons, and the transfer of energy to electromagnetic particles is postponed until the pions have energies in the few 10’s of PeV. The emerging photons induce ordinary showers. Photons or electrons from very short lived particles are the main source of the occasional shower elongations observed alz98 . This is illustrated in Fig. 4.

Eq. 6 can be also used for calculating the radio emission in an approximate fashion to establish Fresnel interference effects. Since this approximation neglects the lateral distribution the results are only expected to be valid for showers which are very elongated because of the LPM effect. For these showers the Fresnel interference effect becomes more important than the lateral distributions. As the observation distance decreases the angular distribution of the pulses becomes wider and the main peak drops. The results show that at the Fresnel distance defined as $`R_F=\pi n\nu (L_s\mathrm{sin}\theta /2)^2/c`$ ($`L_s`$ is the shower length) deviations of the electric field in the Čerenkov direction from the Fraunhofer approximation are below $`10\%`$.

For neutrino induced showers the two types of showers become relevant. Neutrinos always induce a hadronic character shower through the nucleon fragments. In the case of charged current electron neutrino interactions, the emerging shower has mixed hadronic and electromagnetic character and the relative weights of each component depend on $`y`$, the fraction of energy transferred to the nuclear fragments. For very high energy neutrinos this leads to an interesting combination of two showers of different lengths because of the relative importance of the LPM effect for both types of showers. The resulting angular distributions of the radio pulse displays a complex pattern of the two corresponding angular widths. It has been shown that it is in principle possible to use the angular information of the radio pulse to extract the corresponding value of $`y`$ for the interaction alz99 .

## V Summary

We have reviewed the development in calculations of the full diffraction patterns and frequency spectra of hadronic and electromagnetic showers developing in ice. The complex emerging patterns are well understood in terms of the shower properties. We have stressed the dependence on the angular width of the Čerenkov pulse on length and its implication for electromagnetic showers with strong LPM elongations.

## Acknowledgements

We thank D. Saltzberg for helpful comments after carefully reading the manuscript and the organizers P. Gorham and D. Saltzberg for providing this unique oportunity to give a good impulse to this field. This work was supported in part by the European Science Foundation (Neutrino Astrophysics Network N. 86), by the CICYT (AEN99-0589-C02-02) and by Xunta de Galicia (PGIDT00PXI20615PR). The research activities of J.A-M at Bartol Research Institute are supported by the NASA grant NAG5–7009.


format: text

text:
# 1 Introduction

## 1  Introduction

The spectrum of eigenvalues of complex Hamiltonians are often modelled by a random matrix theory, in which the random matrices belong to various ensembles according to the symmetries of the physical problem. The most common space-time symmetries of the Hamiltonian lead to the consideration of ensembles of real, complex or quaternionic random matrices. In the simplest case one considers Gaussian probability distributions. This simple choice is in many cases sufficient since it is now understood that the local statistics of the eigenvalues are universal, i.e. largely independent of the probability distribution. The most commonly studied Gaussian ensembles , called GOE, GUE and GSE, are invariant under the orthogonal, unitary or symplectic groups, respectively, and they all have important applications .

In this article we follow our previous study for the GUE case of the characteristic polynomials of random matrices . If $`X`$ is an $`N\times N`$ random matrix , whose characteristic polynomial is $`det(\lambda X)`$, we consider the average of products of such characteristic polynomials defined by

$$F_k(\lambda _1,\mathrm{},\lambda _k)=<\underset{l=1}{\overset{k}{}}det(\lambda _lX)>.$$
(1)
In the GUE case we have derived in a previous article explicit formulae for those correlation functions, found then their asymptotic behavior for large N, and proved their universality in the short distance limit in which the differences $`\lambda _i\lambda _j`$ is the order of the mean spacing of the eigenvalues of $`X`$. As usual the orthogonal and symplectic cases are more difficult to handle. It turns out that there is a hidden duality in these problems between $`N`$ the size of the matrices, and $`k`$ the number of points in the correlation functions : we may turn the integrals over $`N\times N`$ matrices into integrals over $`k\times k`$ matrices and, since we are interested in large $`N`$ -finite $`k`$ limit, this is the required tool for obtaining the large $`N`$ limit by the saddle-point method. We return below to the GUE case and exhibit its self-duality. However the GOE case turns out to be dual to the GSE. This duality, in the simpler case of all $`\lambda _i`$’s equal, has been discussed recently within the orthogonal polynomial method . In both cases one may use geometry to reduce further the number of integrations. In the GUE case it relies on the Harish-Chandra-Itzykson-Zuber formula (HIZ) . For the GOE case it turns out that the $`Nk`$ duality maps the problem into the GSE case, and the use of the geometry of the symmetric space $`U(2k)/Sp(k)`$ leads to considerable simplifications. At the orders that we have considered one finds that there is a generalization of (HIZ) ; it is well-known that in that case the WKB approximation happens to be exact. For the GSE problem that one finds here, it is WKB, plus a finite number of corrections, which happens to be exact. Therefore we shall begin be re-exposing the unitary case at the light of this duality and of the HIZ formula. We shall then proceed to the GOE ensemble.

## 2 Survey of the unitary ensemble

For the Gaussian unitary ensemble (GUE), the random matrix $`X`$ is a complex $`N\times N`$ Hermitian matrix, with a probability distribution function

$$P(X)=\frac{1}{Z}\mathrm{exp}(\frac{N}{2}\mathrm{trX}^2).$$
(2)
The average $`<\mathrm{}>`$ means integration with the normalized weight $`P(X)`$ , with the Euclidean measure $`_idX_{ii}_{i<j}d\mathrm{}X_{ij}d\mathrm{}X_{ij}`$. It is easily shown that the $`F_2(\lambda _1,\lambda _2)`$ reduces, up to a trivial factor, to the the kernel  $`K_N(\lambda _1,\lambda _2)`$ which characterizes the correlation functions of the eigenvalues of $`X`$. When all the $`\lambda _j`$’s are nearby , i.e. in the short distance scaling region in which $`N`$ is large and the products $`N(\lambda _i\lambda _j)`$ are finite, $`F_{2k}(\lambda _1,\mathrm{},\lambda _{2k})`$ becomes, within an appropriate scaling, a universal function, i.e. independent of the specific distribution $`P(X)`$. When all the $`\lambda _j`$ are equal,

$$F_{2k}(\lambda )=F_{2k}(\lambda ,\mathrm{}.,\lambda )=<[det(\lambda X)]^{2k}>$$
(3)
the 2k-th moment of the characteristic polynomial. In the large N limit, we have derived earlier

$$F_{2k}(\lambda )=\gamma _k[2\pi N\rho (\lambda )]^{k^2},$$
(4)
in which $`\rho (\lambda )`$ is the density of eigenvalues, and $`\gamma _k`$ is a universal factor. This number had been first computed for the circular unitary ensemble by Keating and Snaith who used the Selberg integral formula .

There are several different derivations for those results. Let us here expose the duality which was mentioned in the introduction. We introduce Grassmann variables $`c`$ and $`\overline{c}`$, normalized to

$$𝑑c𝑑\overline{c}e^{iN\overline{c}c}=1.$$
(5)
Then the characteristic polynomial may be written as

$$det(\lambda X)=\underset{a=1}{\overset{N}{}}d\overline{c}_adc_ae^{iN\overline{c}_a(\lambda \delta _{ab}X_{ab})c_b}.$$
(6)
Repeating this k-times

$$\underset{\alpha =1}{\overset{k}{}}det(\lambda _\alpha X)=\underset{a=1}{\overset{N}{}}\underset{\alpha =1}{\overset{k}{}}d\overline{c}_{a\alpha }dc_{a\alpha }e^{iN_{a=1}^N_{\alpha =1}^k\overline{c}_{a\alpha }(\lambda _\alpha \delta _{ab}X_{ab})c_{b\alpha }}.$$
(7)
The (normalized) integration over $`X`$ , in presence of a matrix source $`Y`$, yields

$$𝑑Xe^{\frac{N}{2}\mathrm{TrX}^2+\mathrm{iNtrXY}}=e^{\frac{N}{2}\mathrm{TrY}^2}.$$
(8)
We may now apply this to the matrix $`Y_{ab}=_{\alpha =1}^k\overline{c}_{a\alpha }c_{b\alpha }`$, generated by (7). Then one finds easily that

$$\mathrm{Tr}(\mathrm{Y}^2)=\underset{\alpha \beta =1}{\overset{\mathrm{k}}{}}\underset{\mathrm{a}=1}{\overset{\mathrm{N}}{}}\overline{\mathrm{c}}_{\mathrm{a}\alpha }\mathrm{c}_{\mathrm{a}\beta }\underset{\mathrm{b}=1}{\overset{\mathrm{N}}{}}\overline{\mathrm{c}}_{\mathrm{b}\beta }\mathrm{c}_{\mathrm{b}\alpha }=\mathrm{tr}(\gamma ^2)$$
(9)
with

$$\gamma _{\alpha \beta }=\underset{a=1}{\overset{N}{}}\overline{c}_{a\alpha }c_{a\beta }.$$
(10)
Our notation are as follows : ”$`\mathrm{Tr}`$” refers here to N-dimensional space, whereas ”$`\mathrm{tr}`$” refers to matrices acting in the k-dimensional space. We introduce next an auxiliary matrix $`k\times k`$ hermitian matrix $`B`$, such that

$$e^{\frac{N}{2}\mathrm{tr}\gamma ^2}=𝑑B\mathrm{exp}(\frac{N}{2}\mathrm{trB}^2+\mathrm{Ntr}\gamma \mathrm{B}),$$
(11)
integrate over the Grassmann variables (which are now decoupled in the original N-dimensional space) and end up with

$`F_k(\lambda _1\mathrm{}\lambda _k)`$ $`=`$ $`{\displaystyle 𝑑Bdet(\mathrm{\Lambda }iB)^N\mathrm{exp}}({\displaystyle \frac{N}{2}}\mathrm{trB}^2)`$
$`=`$ $`e^{\frac{N}{2}\mathrm{tr}\mathrm{\Lambda }^2}{\displaystyle 𝑑B(detB)^N\mathrm{exp}}({\displaystyle \frac{N}{2}}\mathrm{trB}^2+\mathrm{iNtr}\mathrm{\Lambda }\mathrm{B})`$
in which $`\mathrm{\Lambda }`$ is the diagonal matrix $`(\lambda _1,\mathrm{},\lambda _k)`$. The problem is thus mapped into Gaussian integrals over $`k\times k`$ hermitian matrices as announced. This dual representation is of course well adapted to the $`k`$-fixed, $`N`$-large, limit that we are considering, since (2) contains $`k^2`$ variables instead of $`N^2`$ in our starting point. It is not difficult to proceed from (2) and derive the scaling results which were given in .

However it turns out that it is simpler , and necessary in view of what comes next for the GOE problem, to integrate out the unitary degrees of freedom in (2). This is done through the HIZ formula  which gives the integral over the unitary group $`U(k)`$ :

$$𝑑Ue^{iN\mathrm{trUXU}^{}\mathrm{Y}}=C_N\frac{det_{1i,jk}e^{iNx_iy_j}}{\mathrm{\Delta }(x_1,\mathrm{},x_k)\mathrm{\Delta }(y_1,\mathrm{},y_k)}$$
(13)
in which the $`x_i`$’s and $`y_j`$’s are the eigenvalues of the Hermitian $`X`$ and $`Y`$ respectively; $`\mathrm{\Delta }(x_1,\mathrm{},x_k)`$ is the Van der Monde determinant

$$\mathrm{\Delta }(x_1,\mathrm{},x_k)=\underset{i<j}{}(x_ix_j).$$
(14)
It is well-known that the formula (13) happens to be exact semi-classically, i.e. if one retains only the sum of the $`k!`$ stationary points in the space of unitary matrices, weighted by the Gaussian fluctuations around each of them. Higher corrections happen to cancel exactly. This leads immediately to an integral over the $`k`$ eigenvalues $`b_l`$ of $`B`$, rather than over the $`k^2`$ matrix elements :

$$F_k(\lambda _1,\mathrm{},\lambda _k)=Ce^{\frac{N}{2}\mathrm{tr}\mathrm{\Lambda }^2}\left(\underset{l=1}{\overset{k}{}}db_lb_l^N\right)e^{\frac{N}{2}_{l=1}^k\left(b_l^22ib_l\lambda _l\right)}\underset{l<l^{}}{\overset{k}{}}\frac{(b_lb_l^{})}{(\lambda _l\lambda _l^{})}$$
(15)
When we consider simply the $`k`$-th moment of the characteristic polynomials, namely the case in which ($`\lambda =\lambda _1=\mathrm{}=\lambda _k`$), the previous formula reduces to

$$F_k(\lambda )=F_k(\lambda ,\mathrm{},\lambda )=Ce^{\frac{Nk}{2}\lambda ^2}\left(\underset{l=1}{\overset{k}{}}db_lb_l^N\right)e^{\frac{N}{2}_{l=1}^k\left(b_l^22ib_l\lambda \right)}\underset{l<l^{}}{\overset{k}{}}(b_lb_l^{})^2$$
(16)
Note that the above representations of the correlation functions of the characteristic polynomials, in terms of integrals over $`k`$ variables, are exact for any size $`N\times N`$ of the random matrices. It is then simple to find the large-N limit of those functions by saddle-point integration. If we focus to even values of $`k`$, and substitute $`2k`$ to $`k`$ (the odd case is doable of course, but it leads to an oscillatory behavior) the saddle point equation for each $`b_l`$ is

$$b_l^2i\lambda _lb_l1=0$$
(17)
whose roots are $`b_l^\pm =(i\lambda _l\pm \sqrt{4\lambda _l^2})/2`$. In the scaling limit, in which the $`\lambda _l\lambda _l^{}`$ are of order $`1/N`$, the eigenvalues $`b_j`$ are very close and one must pay attention to the Vandermonde determinant in the integration measure. Finally the leading saddle-points correspond to equal numbers of $`b_l`$ close to either $`b^+`$ or $`b^{}`$, with $`b^\pm =(i\lambda \pm \sqrt{4\lambda ^2})/2`$. There are thus $`\left(\begin{array}{c}2k\\ k\end{array}\right)=2k!/k!k!`$ saddle-points of equal weight. The combinatorial factor $`\gamma _k`$ of (4) is then simply

$`\gamma _k`$ $`=`$ $`\left(\begin{array}{c}2k\\ k\end{array}\right){\displaystyle \frac{(h_k)^2}{h_{2k}}}`$ (18)
$`=`$ $`{\displaystyle \underset{l=0}{\overset{k1}{}}}{\displaystyle \frac{l!}{(k+l)!}}`$
where we have used

$`h_k`$ $`=`$ $`{\displaystyle \frac{1}{(2\pi )^{k/2}}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \underset{1}{\overset{k}{}}}dx_je^{\frac{1}{2}_{i=1}^kx_i^2}{\displaystyle \underset{l<l^{}}{\overset{k}{}}}(x_lx_l^{})^2`$ (19)
$`=`$ $`{\displaystyle \underset{l=0}{\overset{k}{}}}l!`$
(This formula is used when we consider the gaussian fluctuations near the saddle-point in which $`k`$ of the $`b_j`$’s are near $`b^+`$ and the other half are close to the saddle point $`b^{}`$.)

Finally this representation through an integral over a finite matrix matrix $`B`$, may be generalized to the case of an external matrix source $`A`$ coupled to the random matrix $`X`$ .

$`{\displaystyle 𝑑X\underset{l=1}{\overset{2k}{}}det(\lambda _lX)e^{\frac{N}{2}\mathrm{trX}^2+\mathrm{NtrAX}}}`$ (20)
$`=`$ $`{\displaystyle \frac{1}{\mathrm{\Delta }(\lambda )}}{\displaystyle 𝑑b_l\underset{i=1}{\overset{N}{}}\underset{l=1}{\overset{2k}{}}(a_ib_l)\underset{l<l^{}}{\overset{2k}{}}(b_lb_l^{})e^{\frac{N}{2}{\scriptscriptstyle b_l^2}+iN{\scriptscriptstyle \lambda _lb_l}}}`$
We shall now transpose these techniques to real symmetric random matrices.

## 3 Real symmetric matrices and characteristic polynomials

Again we consider

$$F_k(\lambda _1,\mathrm{},\lambda _k)=𝑑Xe^{\frac{N}{2}\mathrm{trX}^2}\underset{i=1}{\overset{k}{}}det(\lambda _iX),$$
(21)
in which $`\mathrm{X}`$ is a real symmetric $`N\times N`$ matrix, $`\mathrm{X}=\mathrm{X}^\mathrm{T}`$. It is worth remembering that real symmetric matrices form the Lie algebra of the symmetric space $`U(N)/O(N)`$ (the Lie algebra of $`U(N)`$ consists of $`N^2`$ complex hermitian matrices ; the imaginary part of those matrices are the $`N(N1)/2`$ antisymmetric generators of $`O(N)`$ ; the real parts are the $`N(N+1)/2`$ real symmetric generators of the coset)

Using again Grassmann variables, and the representation (6) of the characteristic determinants we are led again to an integration over real symmetric matrices in the presence of the matrix source $`Y=_{\alpha =1}^k\overline{c}_{a\alpha }c_{b\alpha }`$. This gives

$$𝑑Xe^{\frac{N}{2}\mathrm{TrX}^2+\mathrm{iNtrXY}}=e^{\frac{N}{4}\mathrm{Tr}(\mathrm{Y}^2+\mathrm{YY}^\mathrm{T})}.$$
(22)
We have dealt earlier with

$$\mathrm{Tr}(\mathrm{Y}^2)=\mathrm{tr}(\gamma ^2)$$
(23)
which led to the integral (11) over an hermitian $`k\times k`$ matrix. In addition we have here

$$\mathrm{Tr}(\mathrm{YY}^\mathrm{T})=\mathrm{tr}(\mathrm{UV}),$$
(24)
with the matrices $`U`$ and $`V`$ defined by

$$U_{\alpha \beta }=\underset{a=1}{\overset{N}{}}\overline{c}_{a\alpha }\overline{c}_{a\beta }$$
(25)
$$V_{\alpha \beta }=\underset{a=1}{\overset{N}{}}c_{a\alpha }c_{a\beta }$$
(26)
Defining the complex conjugation of Grassmann variables as $`(\overline{c}_1c_2)^{}=\overline{c}_2c_1`$, we have $`\gamma =\gamma ^{}`$, $`V^{}=U`$. Therefore, we may again decompose the remaining quartic terms in the $`c`$’s and $`\overline{c}`$’s as

$$e^{\frac{N}{4}\mathrm{tr}(\mathrm{UV})}=𝑑De^{N\mathrm{tr}(\mathrm{D}^{}\mathrm{D}+\frac{\mathrm{i}}{2}\mathrm{V}^{}\mathrm{D}+\frac{\mathrm{i}}{2}\mathrm{D}^{}\mathrm{V})}$$
(27)
where $`\mathrm{D}`$ is a complex $`k\times k`$ antisymmetrix matrix, $`\mathrm{D}=\mathrm{D}^\mathrm{T}`$. Then we have

$`F_k(\lambda _1,\mathrm{},\lambda _k)`$ $`=`$ $`{\displaystyle d\overline{c}_{a\alpha }dc_{a\alpha }e^{iN_{\alpha =1}^k_{a=1}^N\lambda _\alpha \overline{c}_{a\alpha }c_{a\alpha }}}`$ (28)
$`\times `$ $`{\displaystyle 𝑑B𝑑De^{N\mathrm{tr}(\mathrm{B}^2+\mathrm{D}^{}\mathrm{D})+\mathrm{N}{\scriptscriptstyle \overline{\mathrm{c}}_{\mathrm{a}\alpha }\mathrm{c}_{\mathrm{a}\beta }\mathrm{B}_{\beta \alpha }}}}`$
$`\times `$ $`e^{\frac{i}{2}D_{\alpha \beta }^{}c_{a\beta }c_{a\alpha }\frac{i}{2}D_{\beta \alpha }\overline{c}_{a\alpha }\overline{c}_{a\beta }}.`$
Those auxiliary matrices B and D allow us to integrate over each pair $`\overline{c}_a,c_a`$ independently of the other pairs. It is convenient to define

$$\psi _a=\left(\begin{array}{c}\overline{c}_a\\ c_a\end{array}\right).$$
(29)
For a given antisymmetric matrix M, ($`M=M^T`$), we then have the following formula :

$`{\displaystyle d\overline{c}_{a\alpha }dc_{a\alpha }e^{\frac{iN}{2}\psi _{\alpha a}M_{\alpha \beta }\psi _{\beta a}}}`$ (30)
$`=`$ $`[{\displaystyle 𝑑\overline{c}_\alpha 𝑑c_\alpha e^{\frac{iN}{2}\psi _\alpha M_{\alpha \beta }\psi _\beta }}]^N`$
$`=`$ $`[\mathrm{Pf}M]^N=[detM]^{N/2}`$
where $`\mathrm{Pf}`$ is the pfaffian of the antisymmetrix matrix $`M`$. Applying this to our problem, we deal here with

$$M=\left(\begin{array}{cc}D& \mathrm{\Lambda }iB^T\\ (\mathrm{\Lambda }iB)& D^{}\end{array}\right)$$
(31)
in which $`\mathrm{\Lambda }=\mathrm{diag}(\lambda _1,\mathrm{}\lambda _\mathrm{k})`$. Thus we finally obtain

$$F_k(\lambda _1,\mathrm{}.,\lambda _k)=dBdDe^{N\mathrm{tr}(\mathrm{B}^2+\mathrm{D}^{}\mathrm{D})}[\mathrm{Pf}M]^N$$
(32)
This representation (32) of $`F_k`$ in terms of a finite number of integrals, here $`2k^2k`$ integrals (one $`k\times k`$ hermitian matrix B, one complex antisymmetric $`k\times k`$ matrix D), again exact for any $`N`$, is a solution to the problem. However, contrary to the GUE case, it turns out out that a direct use of the saddle-point equations fail in the scaling limit. In other words every term of the perturbative expansion around the saddle-point turns out to be relevant in the regime in which the products $`N(\lambda _i\lambda _j)`$ are finite.

One may use the remaining invariances of this representation, to reduce further the number of integrations. A unitary transformation $`U`$, among the $`c_a`$, which diagonalizes the Hermitian matrix B, one of the block matrices of $`M`$, transforms $`D`$ into $`DU^{}DU^T`$; in other words one can diagonalize $`B`$ and keep for $`D`$ an antisymmetric matrix. Therefore, applying the Harish-Chandra-Itzykson-Zuber formula  for the integration over the unitary group, i.e. over the relative unitary transformation between the diagonal matrix $`\mathrm{\Lambda }`$ and the eigenbasis of $`B`$, we obtain

$`F_k(\lambda _1,\mathrm{},\lambda _k)`$ $`=`$ $`{\displaystyle \frac{e^{N\mathrm{Tr}\mathrm{\Lambda }^2}}{\mathrm{\Delta }(\mathrm{\Lambda })}}{\displaystyle \underset{l=1}{\overset{k}{}}db_l\mathrm{\Delta }(b)𝑑De^{N{\scriptscriptstyle b_\alpha ^2}+2iN{\scriptscriptstyle \lambda _\alpha b_\alpha }N\mathrm{TrD}^{}\mathrm{D}}}`$ (33)
$`\times `$ $`(\mathrm{Pf}M)^N`$
where now the matrix $`B`$ in $`M`$ is diagonal ; we have reduced the integrations to $`k^2`$ variables, instead of $`2k^2k`$. However it turns out that this is still unsufficient : the inapplicability of the saddle-point method in the scaling limit is still a problem if we proceed from (33). It is thus necessary to return to the underlying geometry of the space of matrices $`M`$ in the representation (31). In order to make the quaternionic structure more apparent we return to (28) and define the spinor

$$\psi _{a\alpha }=\left(\begin{array}{c}\overline{c}_{a\alpha }\\ c_{a\alpha }\end{array}\right)$$
(34)
and the adjoint

$$\overline{\psi }_{a\alpha }=\left(\begin{array}{c}c_{a\alpha }\\ \overline{c}_{a\alpha }\end{array}\right)$$
(35)
Then the quadratic form in the Grassmann variables of (28) takes the form (repeated $`N`$ times for each index $`a`$ that we drop)

$$\underset{\alpha ,\beta =1}{\overset{k}{}}\left(\overline{c}_\alpha c_\beta B_{\beta \alpha }\frac{i}{2}D_{\alpha \beta }^{}c_\beta c_\alpha \frac{i}{2}D_{\beta \alpha }\overline{c}_\alpha \overline{c}_\beta +i\lambda _\alpha \overline{c}_\alpha c_\alpha \right)=\frac{1}{2}\underset{\alpha ,\beta =1}{\overset{k}{}}\overline{\psi }_\alpha \left(q_{\alpha ,\beta }+i\lambda _\alpha \delta _{\alpha \beta }\right)\psi _\beta $$
(36)
in which the $`q_{\alpha \beta }`$ are quaternionic matrix elements, i.e. linear combination of the Pauli matrices. The identification in terms of $`2\times 2`$ matrices is thus

$$q_{\alpha ,\beta }=\left(\begin{array}{cc}B_{\alpha ,\beta }& iD_{\alpha ,\beta }^{}\\ iD_{\alpha ,\beta }& B_{\alpha ,\beta }^{}\end{array}\right).$$
(37)
This defines a self-dual quaternion matrix  , i.e.

$$q_{\alpha \beta }^{}=q_{\beta \alpha }$$
(38)
and $`q_{\alpha \beta }q_{\beta \alpha }`$ is a multiple of identity. Let $`\stackrel{~}{M}_0`$ be the quaternionic matrix whose elements are the quaternion $`q_{\alpha \beta }`$ and $`\stackrel{~}{M}`$ the quaternionic matrix with elements.

$$\stackrel{~}{M}_{\alpha \beta }=q_{\alpha \beta }+i\lambda _\alpha \delta _{\alpha \beta }.$$
(39)
The Grassmannian integration leads to

$$dc_\alpha dc_\alpha \mathrm{exp}\frac{N}{2}\underset{\alpha ,\beta =1}{\overset{k}{}}\overline{\psi }_\alpha \left(q_{\alpha ,\beta }+i\lambda _\alpha \delta _{\alpha \beta }\right)\psi _\beta =Qdet\stackrel{~}{M}=\mathrm{PfM},$$
(40)
in which $`Qdet`$ denotes the quaternionic determinant . In addition

$$tr(B^2+DD^{})=\mathrm{tr}(\stackrel{~}{\mathrm{M}}_0^2)=\underset{\alpha \beta }{}\mathrm{q}_{\alpha \beta }\mathrm{q}_{\beta \alpha }.$$
(41)
It may be clarifying to show this quaternionic construction explicitely for the $`k=2`$ case. There one has

$$M=\left(\begin{array}{cccc}0& d& \lambda _1iB_{11}& iB_{21}\\ d& 0& iB_{12}& \lambda _2iB_{22}\\ \lambda _1+iB_{11}& iB_{12}& 0& d^{}\\ iB_{21}& \lambda _2+iB_{22}& d^{}& 0\end{array}\right),$$
(42)
and

$$\mathrm{PfM}=|\mathrm{d}|^2+(\lambda _1\mathrm{iB}_{11})(\lambda _2\mathrm{iB}_{22})+|\mathrm{B}_{12}|^2.$$
(43)
The equivalent quaternionic construction is

$$\stackrel{~}{M}=\left(\begin{array}{cc}q_{11}+i\lambda _1& q_{12}\\ q_{21}& q_{22}+i\lambda _2\end{array}\right)$$
(44)
with

$`q_{11}`$ $`=`$ $`B_{11}\mathrm{𝟏},q_{12}=(\mathrm{}B_{12})\mathrm{𝟏}+i(\mathrm{}B_{12})\sigma _3+i(\mathrm{}d)\sigma _1+i(\mathrm{}d)\sigma _2`$
$`q_{21}`$ $`=`$ $`q_{12}^{},q_{22}=B_{22}\mathrm{𝟏}`$ (45)
and

$`Qdet(\stackrel{~}{M})=(q_{11}+i\lambda _1)(q_{22}+i\lambda _2)q_{12}q_{12}^{}`$
$`=(B_{11}+i\lambda _1)(B_{22}+i\lambda _2)+|B_{12}|^2+|d|^2=\mathrm{PfM}.`$ (46)
Therefore we end up for the correlation functions with the following duality :

$$F_k(\lambda _1,\mathrm{},\lambda _k)=𝑑\stackrel{~}{M}(Qdet\stackrel{~}{M})^Ne^{N\mathrm{tr}(\stackrel{~}{\mathrm{M}_0}^2)}.$$
(47)
The original integral over real symmetric $`N\times N`$ matrices is replaced by integrals over quaternionic matrices which depend upon $`(2k^2k)`$ degrees of freedom. Those matrices are the generators for the symmetric space $`U(2k)/Sp(k)`$.

This representation (47) in terms of a finite number of integration variables is a priori well adapted to the large N-limit. However it turns out that in the scaling limit of interest, the contributions of the non-gaussian fluctuations around the saddle-points are all relevant. Therefore it is necessary to eliminate first the ”angular” degrees of freedom. When all the $`\lambda _i`$ are equal, namely if we consider the moments of the characteristic polynomials, one can simply diagonalize the symplectic matrices in terms of $`k`$ eigenvalues and then proceed to the large N limit. This is done in the next section. However if the $`\lambda _i`$’s are unequal we need some equivalent of the HIZ formalism, which will be described afterwards.

## 4 Moments of the characteristic polynomials

We first note the trivial $`k=1`$ case : $`F_1(\lambda )=det(\lambda X)`$ is simply

$$F_1(\lambda )=_{\mathrm{}}^{\mathrm{}}𝑑be^{Nb^2}(\lambda ib)^N,$$
(48)
which, up to a trivial factor, is the Hermite polynomial $`\mathrm{H}_\mathrm{N}(\sqrt{\mathrm{N}}\lambda )`$ which has an oscillatory behavior for large $`N`$ when $`\lambda `$ belongs to the support of Wigner’s semi-circle. Therefore we consider from now on the more interesting even correlation functions.

When all the $`\lambda _i`$’s are equal, the matrix $`\mathrm{\Lambda }`$ is a multiple of identity, and one can diagonalize the quaternionic matrix $`\stackrel{~}{M}`$ through a transformation belonging to the symplectic group $`Sp(2k)`$. The transformation of $`\stackrel{~}{M}`$ into the diagonal matrix $`T=\mathrm{diag}(\mathrm{t}_1,\mathrm{}\mathrm{t}_\mathrm{k})`$ yields the Jacobian $`J=[\mathrm{\Delta }(t)]^4`$, ($`\mathrm{\Delta }(t)`$ is the Vandermonde determinant $`_{i<j}(t_it_j)`$).

This gives simply

$`F_{2k}(\lambda )`$ $`=`$ $`det(\lambda X)^{2k}`$ (49)
$`=`$ $`C{\displaystyle \underset{l=1}{\overset{2k}{}}dt_l(\underset{l=1}{\overset{2k}{}}t_l)^N\underset{l<l^{}}{}(t_lt_l^{})^4e^{N{\scriptscriptstyle t_l^2}+i2N\lambda {\scriptscriptstyle t_l}}}`$
The integral representation (49) is well suited to the study of the large N limit. Exponentiating $`t_l^N`$ term as $`e^{N\mathrm{log}t_l}`$, the integrand is of the form $`\mathrm{exp}(N{\displaystyle \underset{1}{\overset{2k}{}}}f(t_l))`$ with

$$f(t)=t^22i\lambda t\mathrm{log}t$$
(50)
The saddle points for every $`t_l`$ are solutions of $`f^{}(t_c)=0`$,i.e.

$$2t^22i\lambda t1=0$$
(51)
The two solutions are given by

$$t^\pm =\frac{1}{2}(i\lambda \pm \sqrt{2\lambda ^2})$$
(52)
The difference $`(t^+t^{})`$ is proportional to the semi-circular density of eigenvalues of the GOE ensemble :

$$t^+t^{}=\pi \rho (\lambda )$$
(53)
where $`\rho (\lambda )=\sqrt{2\lambda ^2}/\pi `$. Expanding $`t_l`$ around either $`t^+`$ or $`t^{}`$, we find that the leading saddle-points are those in which half of the $`t_l(l=1,\mathrm{},2k)`$ are near $`t^+`$ and the remaining half near $`t^{}`$. (Other choices give oscillatory contributions in $`\mathrm{exp}(N{\displaystyle \underset{1}{\overset{2k}{}}}f(t_l))`$ , which damp the large $`N`$-limit). Therefore we have to add the $`(2k)!/(k!k!)`$ leading saddle-points corresponding to the distribution of half of the $`t_l`$’s near $`t^+`$, and the other half near $`t^{}`$. The measure term given by the 4-th power of the Vandermonde determinant, yields a factor $`(\pi \rho (\lambda ))^{4k^2}`$ from the $`k`$ variables near $`t^+`$ and the $`k`$ near $`t^{}`$. The exponent $`f`$ is then expanded around $`t^+`$ or $`t^{}`$, and the remaining integral factorizes into an integration around $`t^+`$ and another one around $`t^{}`$. The integration around $`t^+`$ is

$$_{\mathrm{}}^{\mathrm{}}\underset{l=1}{\overset{k}{}}dt_le^{\frac{N}{2}f^{\prime \prime }(t^+)(tt^+)^2}\underset{i<j}{\overset{k}{}}(t_it_j)^4=(\frac{1}{\sqrt{Nf^{\prime \prime }(t^+)}})^{2k^2k}\underset{l=1}{\overset{k}{}}(2l)!.$$
(54)
Noting $`t^+t^{}=1`$, and $`t^+t^{}=\pi \rho (\lambda )`$, we find $`f^{\prime \prime }(t^+)f^{\prime \prime }(t^{})=(\pi \rho (\lambda ))^2`$. We need to fix the normalization constant C in (49). It is obtained from the integral,

$$dt_le^{\frac{N}{2}{\scriptscriptstyle t_l^2}}\underset{i<j}{\overset{2k}{}}(t_it_j)^4=(\frac{1}{N})^{4k^2k}\underset{l=1}{\overset{2k}{}}(2l)!$$
(55)
The constant $`C`$ in (49) is thus the inverse of this number. This constant appears as a normalization for the n-point correlation function of the Gaussian symplectic ensemble . Thus including this normalization constant, $`F_{2k}(\lambda )`$ becomes in the large N limit as

$$F_{2k}^{(GOE)}=\gamma _kN^{2k^2}(2\pi \rho (\lambda ))^{2k^2+k}$$
(56)
$`\gamma _k`$ $`=`$ $`{\displaystyle \frac{(2k)!}{k!k!}}{\displaystyle \frac{[\underset{l=1}{\overset{k}{}}(2l)!]^2}{_{l=1}^{2k}(2l)!}}`$ (57)
$`=`$ $`{\displaystyle \underset{l=1}{\overset{k}{}}}{\displaystyle \frac{(2l1)!}{(2k+2l1)!}}`$
For example, in the case of $`F_2(\lambda )=<(det(\lambda X))^2>`$, it gives

$$F_2(\lambda )^{(GOE)}=\frac{1}{6}N^2(2\pi \rho (\lambda ))^3$$
(58)
This result agrees with the result which one would deduce from $`\mathrm{lim}\rho (\lambda _1,\lambda _2)/(\lambda _1\lambda _2)`$, where the limit means $`\lambda _1\lambda `$ and $`\lambda _2\lambda `$, ( details are given in appendix A).

The value of $`\gamma _k`$ agrees with the result of COE (circular orthogonal ensemble) found by Keating and Snaith  through the Selberg integral formula. In the COE case, however, the density of state is a constant, and the factor $`\rho (\lambda )`$ is absent. This result is to be compared with the earlier result for the GUE,

$$F_{2k}^{(GUE)}=\gamma _k(2N\pi \rho (\lambda ))^{k^2}$$
(59)
where $`\gamma _k`$ is given by (18). This universal constant $`\gamma _k`$ appears also in the average of the moments of the Riemann $`\zeta `$-function, and it has a number theoretical meaning .

## 5 Correlations of characteristic polynomials

The integral representation (2) of the correlations functions $`F_k(\lambda _1,\mathrm{},\lambda _k)`$ is not unitary invariant unless the $`\lambda _i`$’s are all equal. Therefore if we parametrize the matrix $`B`$ as $`B=U^{}bU`$, in which $`b`$ is a diagonal matrix, we have to consider the HIZ integral

$$𝑑U\mathrm{exp}i\mathrm{trU}^{}\mathrm{bU}\mathrm{\Lambda }=\frac{det\mathrm{exp}i\lambda _ib_j}{\mathrm{\Delta }(b)\mathrm{\Delta }(\lambda )}$$
(60)
which is well-known to be WKB exact . This explains why , in the GUE case, it is equally possible, to apply the saddle-point method with or without integrating out the unitary group.

In the symplectic case we are not aware of any similar explicit result ; however it will be shown now, at least for the lowest values of $`k`$, that the integral over the symplectic group can be performed exactly. The result is remarkably that, in this case, WKB plus a finite number of corrections is exact.

In the $`k=2`$ case, we have

$$F_2(\lambda _1,\lambda _2)=<det(\lambda _1X)det(\lambda _2X)>$$
(61)
As shown in the previous section, it is given by the integral (49). We first evaluate explicitly the angular integral. We first diagonalize $`B`$ by a unitary transformation, and then write the eigenvalues $`b_1`$ and $`b_2`$ in terms of new parameters $`t`$ and $`c`$,

$`b_1`$ $`=`$ $`(1c)t_1+ct_2`$
$`b_2`$ $`=`$ $`ct_1+(1c)t_2`$ (62)
Then we have $`b_1+b_2=t_1+t_2`$, and $`b_1b_2=(12c)(t_1t_2)`$. Since the integrand is a function of $`|d|^2`$, we change variable to $`|d|^2=c(1c)(t_1t_2)^2`$ ; then we have $`b_1^2+b_2^2+2|d|^2=t_1^2+t_2^2`$, and $`b_1b_2|d|^2=t_1t_2`$. Finally since the parameter $`c`$ is restricted to the interval $`0<c<1`$, we replace it by $`c=\mathrm{sin}^2\theta `$.

This leads to

$`F_2(\lambda _1,\lambda _2)`$ $`=`$ $`{\displaystyle _0^1}𝑑c{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑t_1𝑑t_2(12c){\displaystyle \frac{(t_1t_2)^3}{(\lambda _1\lambda _2)}}(t_1t_2)^N`$ (63)
$`\times `$ $`e^{N(t_1^2+t_2^2)+2iN(t_1\lambda _1+t_2\lambda _2)2iNc(t_1t_2)(\lambda _1\lambda _2)}`$
The integration over $`c`$ yields

$`F_2(\lambda _1,\lambda _2)`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑t_1𝑑t_2(t_1t_2)^Ne^{N(t_1^2+t_2^2)+2iN(t_1\lambda _1+t_2\lambda _2)}`$ (64)
$`\times `$ $`[{\displaystyle \frac{1}{N}}({\displaystyle \frac{t_1t_2}{\lambda _1\lambda _2}})^2+{\displaystyle \frac{i}{N^2}}{\displaystyle \frac{t_1t_2}{(\lambda _1\lambda _2)^3}}]`$
When $`\lambda _1=\lambda _2=\lambda `$, it reduces as expected to (49). This formula may be easily checked for finite values of $`N`$ , since it reduces to Gaussian integrals over $`t_1`$ and $`t_2`$ ; for instance in the simplest case $`N=1`$, it gives $`F_2(\lambda _1,\lambda _2)\lambda _1\lambda _2+1`$, which agrees with the direct calculation (in this case the trivial integral over the real axis $`{\displaystyle 𝑑x(\lambda _1x)(\lambda _2x)e^{\frac{1}{2}x^2}}`$).

However the representation (64), which is exact for any $`N`$, makes it clear
i) that the large $`N`$-limit may be found through a saddle-point integration over $`t_1`$ and $`t_2`$ ; this will be done below.
ii) that in the universal local limit of interest, in which $`N`$ goes to infinity, $`\lambda _1\lambda _2`$ goes to zero and $`N(\lambda _1\lambda _2)`$ remains finite, the large $`N`$-limit could not have been taken earlier. If, for instance, we had used the saddle-point method at the level of (33), we would have missed the second term in the bracket of (64). If, at the early level of (33), we had recognized that the regime of interest requires to expand beyond the Gaussian approximation to the saddle-point, it would have appeared unexpectedly that the expansion stops after the first correction. Therefore (64) could have been obtained by a semi-classical approximation with a finite number of corrections, here just one. This is analogous, although not as simple, to the Harish-Chandra-Itzykson-Zuber formula for the GUE case, which is semi-classically exact, without any correction term . The above integration over $`c`$ is therefore , for $`k=2`$, the corresponding HIZ formula for the symplectic group.

For higher values of $`k`$ we need a more elaborate strategy. The HIZ formula may be easily derived by considering the Laplacian operator ,

$$L=\frac{^2}{X_{ij}^2}.$$
(65)
Its eigenfunctions are plane waves

$$Le^{iN\mathrm{tr}\mathrm{\Lambda }\mathrm{X}}=(N^2\mathrm{tr}\mathrm{\Lambda }^2)\mathrm{e}^{\mathrm{iNtr}\mathrm{\Lambda }\mathrm{X}}.$$
(66)
One can construct a unitary invariant eigenfunction of $`L`$, for the same energy $`N^2\mathrm{tr}\mathrm{\Lambda }^2`$, by the superposition

$$I=𝑑Ue^{iN\mathrm{tr}\mathrm{\Lambda }\mathrm{UXU}^{}},$$
(67)
which is nothing but the HIZ integral. The integral beeing unitary invariant, it is a function of the $`k`$ eigenvalues $`t_i`$ of $`X`$. The same considerations hold for the three ensembles $`\beta =1,2`$ and 4, corresponding to the orthogonal, unitary and symplectic ensemble, with

$$I=e^{N\mathrm{tr}\mathrm{\Lambda }\mathrm{gXg}^1}𝑑g.$$
(68)
The Laplacian, expressed in terms of a differential operator on the eigenvalues $`t_i`$ reads

$$[\underset{i=1}{\overset{k}{}}\frac{^2}{t_i^2}+\beta \underset{i=1,(ij)}{\overset{k}{}}\frac{1}{t_it_j}\frac{}{t_i}]I=ϵI,$$
(69)
with the eigenvalue $`ϵ`$

$$ϵ=N^2\underset{i=1}{\overset{k}{}}\lambda _i^2$$
(70)
The $`t`$-dependent eigenfunctions of this Schrödinger operator have a scalar product given by the measure

$$\phi _1|\phi _2=𝑑t_1\mathrm{}𝑑t_k|\mathrm{\Delta }(t_1\mathrm{}t_k)|^\beta \phi _1^{}(t_1\mathrm{}t_k)\phi _2(t_1\mathrm{}t_k)$$
(71)
The measure becomes trivial if one multiplies the wave function by $`|\mathrm{\Delta }|^{\beta /2}`$ . Thus if one changes $`I(t)`$ to

$$\psi (t_1\mathrm{}t_k)=|\mathrm{\Delta }(t_1\mathrm{}t_k)|^{\beta /2}I(t_1\mathrm{}t_k),$$
(72)
one obtains the Hamiltonian,

$$[\underset{i=1}{\overset{k}{}}\frac{^2}{t_i^2}\beta (\frac{\beta }{2}1)\underset{i<j}{}\frac{1}{(t_it_j)^2}]\psi =ϵ\psi .$$
(73)
For $`\beta =2`$, the solution is again given by plane waves in the $`t_i`$ and (taking into account the symmetry under permutations of $`I`$), one obtains the HIZ formula.

In the $`\beta =4`$ case, the problem is less trivial, but simple for finite values of $`k`$. For $`k=2`$, a solution of this equation is

$$\psi _0=e^{iN(\lambda _1t_1+\lambda _2t_2)}(1+\frac{2i}{N(t_1t_2)(\lambda _1\lambda _2)})$$
(74)
The symmetry of $`I`$ under permutation of thed $`t_i`$’s leads then to the solution

$$\psi =e^{iN(\lambda _1t_1+\lambda _2t_2)}(1+\frac{2i}{N(t_1t_2)(\lambda _1\lambda _2)})+e^{iN(\lambda _1t_2+\lambda _2t_1)}(1+\frac{2i}{N(t_2t_1)(\lambda _1\lambda _2)})$$
(75)
Then, after multiplication by the Vandermonde factor, we obtain the required symplectic HIZ formula (for $`k=2`$),

$$I=\frac{1}{(t_1t_2)^2(\lambda _1\lambda _2)^2}\psi $$
(76)
For general $`k`$ ($`\beta `$ = 4), the solution of (73) is of the form

$$\psi _0=e^{iN(\lambda _1t_1+\mathrm{}+\lambda _kt_k)}\chi $$
(77)
where $`\chi `$ satifies

$$[\underset{i=1}{\overset{k}{}}\frac{^2}{t_i^2}+2iN\underset{i=1}{\overset{k}{}}\lambda _i\frac{}{t_i}\underset{i<j}{}\frac{4}{(t_it_j)^2}]\chi =0$$
(78)
The operator $`{\displaystyle \underset{i=1}{\overset{k}{}}}{\displaystyle \frac{^2}{t_i^2}}{\displaystyle \underset{i<j}{}}{\displaystyle \frac{4}{(t_it_j)^2}}`$ annihilates the function $`\mathrm{\Delta }^1(t_1,\mathrm{},t_k)`$ . Consequently the solution of (78) may be written

$$\chi (t_1\mathrm{}t_k)=\frac{f(t_1\mathrm{}t_k)}{\mathrm{\Delta }(t_1\mathrm{}t_k)}$$
(79)
in which $`f(t_1\mathrm{}t_k)`$ is a polynomial of degree $`k(k1)/2`$ in the $`t_a`$’s. Defining

$$\tau _{ij}=N(\lambda _i\lambda _j)(t_it_j)$$
(80)
one finds for $`k=3`$

$`\chi =[1{\displaystyle \frac{2}{i}}({\displaystyle \frac{1}{\tau _{12}}}+{\displaystyle \frac{1}{\tau _{23}}}+{\displaystyle \frac{1}{\tau _{31}}})`$
$`4({\displaystyle \frac{1}{\tau _{12}\tau _{23}}}+{\displaystyle \frac{1}{\tau _{23}\tau _{31}}}+{\displaystyle \frac{1}{\tau _{31}\tau _{12}}})12i{\displaystyle \frac{1}{\tau _{12}\tau _{23}\tau _{31}}}].`$ (81)
Again, as for $`k=2`$, one sees that the successive terms in the r.h.s. of (5) are of same order in the limit of interest, and again they could have been obtained through a finite number of corrections to a semi-classical calculation. It is remarkable that the series of $`\chi `$ stops at the order the inverse of the Vandermonde ; thus the symplectic HIZ integral is expressed as the sum of a finite number of terms. The successive coefficients of each term are determined by the equation (73).

Using this modified HIZ formula for the symplectic case, we obtain for the k=3 case, $`F_3(\lambda _1,\lambda _2,\lambda _3)`$ which is expressed by

$`F_3(\lambda _1,\lambda _2,\lambda _3)={\displaystyle 𝑑t_1𝑑t_2𝑑t_3e^{N(t_1^2+t_2^2+t_3^2)+2iN(\lambda _1t_1+\lambda _2t_2+\lambda _3t_3)}(t_1t_2t_3)^N(\frac{\mathrm{\Delta }(t)}{\mathrm{\Delta }(\lambda )})^2}`$ (82)
$`\times `$ $`[1+i({\displaystyle \frac{1}{N(\lambda _1\lambda _2)(t_1t_2)}}+{\displaystyle \frac{1}{N(\lambda _2\lambda _3)(t_2t_3)}}+{\displaystyle \frac{1}{N(\lambda _3\lambda _1)(t_3t_1)}})`$
$``$ $`{\displaystyle \frac{1}{N^2(\lambda _1\lambda _2)(\lambda _2\lambda _3)(t_1t_2)(t_2t_3)}}{\displaystyle \frac{1}{N^2(\lambda _2\lambda _3)(\lambda _3\lambda _1)(t_2t_3)(t_3t_1)}}`$
$``$ $`{\displaystyle \frac{1}{N^2(\lambda _3\lambda _1)(\lambda _1\lambda _2)(t_3t_1)(t_1t_2)}}i{\displaystyle \frac{3}{2}}{\displaystyle \frac{1}{N^3\mathrm{\Delta }(\lambda )\mathrm{\Delta }(t)}}]`$
where $`\mathrm{\Delta }(t)=(t_1t_2)(t_2t_3)(t_3t_1)`$. (Note that the symmetries of the integrand allowed us to keep only the single solution (5), without adding permutations). For low values of $`N`$ (i.e. N=1 or 2, one verifies easily this result by a direct integration over $`1\times 1`$ or $`2\times 2`$ matrices.

Higher values of $`k`$ may be handled in a similar way, but the combinatorics become quite heavy. For instance in an appendix the solution of the case $`k=4`$ is given explicitly and, although again it consists of a finite number of terms, it is quite cumbersome.

As is now clear , those integral representations make it easy to find the scaling limit (large N, finite $`N(\lambda _i\lambda _j)`$). For instance for $`k=2`$ one finds the saddle point values of $`t_1`$ and $`t_2`$ from (64) in the large N limit,

$$t_i=\frac{i}{2}\lambda _i+\frac{ϵ}{2}\sqrt{2\lambda _i^2}$$
(83)
where $`i=1,2`$ and $`ϵ=\pm 1`$. We use the parametrization

$`\lambda _i=\sqrt{2}\mathrm{cos}\theta _i`$. There are a priori four saddle-points given by (83), but the two dominant ones are $`t_1=\frac{i}{\sqrt{2}}e^{i\theta _1ϵ}`$, $`t_2=\frac{1}{\sqrt{2}}e^{i\theta _2ϵ}`$ for $`ϵ=\pm 1`$. In the short distance limit, $`N`$ large and $`N(\theta _1\theta _2)=N\theta _{12}`$ finite , we obtain

$$F_2(\lambda _1,\lambda _2)=\underset{ϵ=\pm 1}{}e^{2iϵN\theta _{12}\mathrm{sin}^2\theta }[\frac{1}{(N\theta _{12})^2}+\frac{iϵ}{2(N\theta _{12})^3\mathrm{sin}^2\theta }]$$
(84)
where $`\theta =(\theta _1+\theta _2)/2`$. The semi-circle density of states is given by $`\rho (\lambda )=\frac{\sqrt{2}}{\pi }\mathrm{sin}\theta `$, and $`\theta _{12}=\frac{1}{\sqrt{2}}\frac{\lambda _1\lambda _2}{\mathrm{sin}\theta }`$. Thus we obtain in the scaling short distant limit,

$$F_2(\lambda _1,\lambda _2)=C[\frac{\mathrm{cos}x}{x^2}\frac{\mathrm{sin}x}{x^3}]$$
(85)
where $`x=\pi N(\lambda _1\lambda _2)\rho (\lambda )`$. It is interesting to note that this function may be expressed as a half-integer Bessel function, since $`({\displaystyle \frac{\mathrm{cos}x}{x^2}}{\displaystyle \frac{\mathrm{sin}x}{x^3}})=\sqrt{{\displaystyle \frac{\pi }{2x^3}}}J_{3/2}(x)`$. In the unitary case, the sine kernel is similarly a half integer Bessel function since $`{\displaystyle \frac{\mathrm{sin}x}{x}}=\sqrt{{\displaystyle \frac{\pi }{2x}}}J_{1/2}(x)`$.

## 6 Extension to an external matrix source

In the GUE case, when an external matrix source $`A`$ is coupled to an Hermitian random matrix $`X`$, as we have discussed earlier, $`F_{2k}(\lambda _1,\mathrm{},\lambda _{2k})`$ is given by (20). The degrees of freedom provided by the eigenvalues of $`A`$ are useful to study a number of new universality classes . For instance by tuning the eigenvalues $`a_i`$ of the external source matrix $`A`$, we can study the problem of a closing gap in the spectrum of random hermitian matrices . Thus it is interesting to consider this external source problem for real symmetric matrices as well .

One can always assume that the external source matrix $`A`$ is diagonal. In the method of integration over Grassmann variables used in section three, it is simple to include the external matrix $`A`$ :

$$e^{\frac{N}{2}\mathrm{trX}^2+\mathrm{NtrAX}+\mathrm{iNtrXY}}𝑑X=e^{\frac{N}{2}\mathrm{tr}[(\mathrm{Y}\mathrm{iA})^2+(\mathrm{Y}\mathrm{iA})(\mathrm{Y}^\mathrm{T}\mathrm{iA})}$$
(86)
where $`Y=_\alpha \overline{c}_{a\alpha }c_{b\alpha }`$. Since $`A`$ is diagonal, the term $`\mathrm{trAY}`$ gives simply the extra term $`\mathrm{exp}[iNa_n\overline{c}_{n\alpha }c_{n\alpha }]`$ in the integrand of $`F_k`$ in (28), . Therefore, repeating the calculations of section 2, we find that eq.(30) is modified as follows :

$`{\displaystyle d\overline{c}_{j\alpha }dc_{j\alpha }e^{\frac{iN}{2}\psi _{\alpha j}M_{\alpha \beta }^{(j)}\psi _{\beta j}}}`$ (87)
$`=`$ $`{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle 𝑑\overline{c}_\alpha 𝑑c_\alpha e^{\frac{iN}{2}\psi _\alpha M_{\alpha \beta }^{(j)}\psi _\beta }}`$
$`=`$ $`{\displaystyle \underset{j}{}}[\mathrm{Pf}M^{(j)}]={\displaystyle \underset{j}{}}[detM^{(j)}]^{/2}`$
where $`\mathrm{Pf}M^{(a)}`$ is the pfaffian of the antisymmetrix matrix $`M^{(a)}`$ given by

$$M^{(j)}=\left(\begin{array}{cc}D& \mathrm{\Lambda }+a_j\mathrm{𝟏}iB^T\\ (\mathrm{\Lambda }+a_j\mathrm{𝟏}iB)& D^{}\end{array}\right).$$
(88)
in which $`\mathrm{\Lambda }=\mathrm{diag}(\lambda _1,\mathrm{}\lambda _\mathrm{k})`$ and $`a_j\mathrm{𝟏}=diag(a_j,\mathrm{},a_j)`$. Thus we finally obtain

$$F_k(\lambda _1,\mathrm{}.,\lambda _k)=dBdDe^{N\mathrm{tr}(\mathrm{B}^2+\mathrm{D}^{}\mathrm{D})}\underset{j=1}{\overset{N}{}}[\mathrm{Pf}M^{(j)}]$$
(89)
This integral can be expressed in terms of a quarternion matrix $`Q`$, which can be diagonalized by the symplectic group $`Sp(k)`$. When all the $`\lambda _i`$’s are equal to a single $`\lambda `$, we get

$$F_k(\lambda ,\mathrm{},\lambda )=\underset{l=1}{\overset{k}{}}\underset{j=1}{\overset{N}{}}(t_lia_j)\underset{l<l^{}}{}(t_lt_l^{})^4e^{N{\scriptscriptstyle t_l^2}+2iN\lambda {\scriptscriptstyle t_l}}\underset{l=1}{\overset{k}{}}dt_l.$$
(90)
For the case of different $`\lambda _i`$’s , this formula is modified by an extra factor as in the previous section.

As an example of the usefulness of the above representation, we choose an external source with only two opposite eigenvalues $`\pm c`$, with half of the eigenvalues equal to $`+c`$ and the other half to $`c`$. This gives a factor $`(t_l^2+c^2)^{N/2}=\mathrm{exp}\frac{N}{2}\mathrm{log}(t^2+c^2)`$ in the integrand. Expanding it in powers of $`t^2`$, the total coefficient of $`t^2`$ in the exponent vanishes for $`c=1/\sqrt{2}`$. Therefore in the large N limit, we obtain at this new critical point

$$<[det(X)]^k>=e^{N\mathrm{trQ}^4}𝑑Q$$
(91)
where $`Q`$ is a $`k\times k`$ symmetric quaternionic matrix, and the the average $`<\mathrm{}>`$ is evaluated with the distribution in the presence of the external source whose eigenvalues are $`\pm c=\pm 1/\sqrt{2}`$.

One could make other choices for the eigenvalues of the external source matrix $`A`$, and obtain thereby higher multicritical points with terms such as $`\mathrm{trQ}^\mathrm{n}`$ in the exponent, in analogy with the GUE case in an external matrix source .

## 7 Summary

In this article, an exact representation of the k-point functions $`_{a=1}^kdet(\lambda _aX)`$, averaged over $`N\times N`$ real symmetric random matrices, has been derived in terms of an integral over quaternionic $`k\times k`$ matrices, invariant under the unitary symplectic group. This representation leads to an easy calculation of the moments of the characteristic polynomials ($`\lambda _1=\mathrm{}=\lambda _k`$). In the large $`N`$-limit one finds

$$F_{2k}^{(GOE)}=\underset{l=1}{\overset{k}{}}\frac{(2l1)!}{(2k+2l1)!}N^{2k^2}(2\pi \rho (\lambda ))^{2k^2+k},$$
(92)
to be compared to the earlier result for the GUE,

$$F_{2k}^{(GUE)}=\underset{l=0}{\overset{k1}{}}\frac{l!}{(k+l)!}(2N\pi \rho (\lambda ))^{k^2}.$$
(93)
For unequal $`\lambda _a`$’s, in spite of the fact that the integral representation involves a finite number of variables, in the large $`N`$-limit the corrections to the saddle-point, in the scaling regime $`N(\lambda _i\lambda _j)`$ finite, are not negligible. A generalization of the HarishChandra-Itzyson-Zuber formula is shown to solve the problem. Remarkably this formula is”nearly” semi-classical, in the sense that it happens that the semi-classical expansion terminates after a few terms, a number of terms which increases with $`k`$ but not with $`N`$. Then the saddle-point method may easily be applied for large $`N`$, and this leads to explicit asymptotic formulae for the correlation functions of the characteristic polynomials. Finally this may be generalized to include an external matrix source in the probability measure. Real symmetric random matrices appear as models of numerous physical time-reversal invariant Hamiltonians. For instance the orthogonal matrix model with an external source has been investigated as a model of glassy behavior . The results of the present work for the moments and for the correlation functions in an external source may be of interest for such problems.

Appendix A: The solution for $`k=4`$

From (78) and (79), the polynomial $`f`$ satisfies

$$\underset{i=1}{\overset{k}{}}\frac{^2}{t_i^2}f+2\underset{i=1}{\overset{k}{}}(\frac{f}{t_i}+iN\lambda _if)(\mathrm{\Delta }\frac{}{t_i}\frac{1}{\mathrm{\Delta }})+2iN\underset{i=1}{\overset{k}{}}\lambda _i\frac{f}{t_i}=0$$
(A.1)
The solution of this equation is obtained by a perturbation expansion in powers of the $`\lambda _i`$’s, but it ends at the level of the Vandermonde $`\mathrm{\Delta }(\lambda _1\mathrm{}\lambda _4)`$. Using the notation of (80), $`\tau _{ij}=N(t_it_j)(\lambda _i\lambda _j)`$, we obtain

$`f`$ $`=`$ $`C[1{\displaystyle \frac{i}{4}}(\tau _{12}+\tau _{13}+\tau _{14}+\tau _{23}+\tau _{24}+\tau _{34})`$ (A.2)
$``$ $`{\displaystyle \frac{1}{12}}(\tau _{12}\tau _{13}+\tau _{12}\tau _{14}+\tau _{13}\tau _{14}+\tau _{12}\tau _{23}+\tau _{23}\tau _{24}+\tau _{12}\tau _{24}`$
$`+`$ $`\tau _{14}\tau _{34}+\tau _{14}\tau _{24}+\tau _{24}\tau _{34}+\tau _{23}\tau _{34}+\tau _{13}\tau _{34}+\tau _{13}\tau _{23})`$
$``$ $`{\displaystyle \frac{1}{18}}(\tau _{12}\tau _{34}+\tau _{13}\tau _{24}+\tau _{14}\tau _{23})`$
$`+`$ $`{\displaystyle \frac{i}{24}}(\tau _{12}\tau _{13}\tau _{14}+\tau _{12}\tau _{23}\tau _{24}+\tau _{13}\tau _{23}\tau _{34}+\tau _{14}\tau _{24}\tau _{34})`$
$`+`$ $`{\displaystyle \frac{i}{36}}(\tau _{12}\tau _{13}\tau _{23}+\tau _{12}\tau _{14}\tau _{24}+\tau _{13}\tau _{14}\tau _{34}+\tau _{23}\tau _{24}\tau _{34}`$
$`+`$ $`\tau _{14}\tau _{34}\tau _{23}+\tau _{14}\tau _{24}\tau _{23}+\tau _{12}\tau _{24}\tau _{34}+\tau _{12}\tau _{23}\tau _{34}`$
$`+`$ $`\tau _{12}\tau _{14}\tau _{23}+\tau _{13}\tau _{14}\tau _{23}+\tau _{12}\tau _{13}\tau _{34}+\tau _{12}\tau _{14}\tau _{34}`$
$`+`$ $`\tau _{13}\tau _{34}\tau _{24}+\tau _{13}\tau _{24}\tau _{23}+\tau _{14}\tau _{24}\tau _{13}+\tau _{13}\tau _{12}\tau _{24})`$
$`+`$ $`{\displaystyle \frac{1}{72}}(\tau _{12}\tau _{23}\tau _{34}\tau _{14}+\tau _{12}\tau _{13}\tau _{24}\tau _{34}+\tau _{13}\tau _{14}\tau _{24}\tau _{23}+\tau _{12}\tau _{14}\tau _{24}\tau _{34}+\tau _{12}\tau _{14}\tau _{24}\tau _{23}`$
$`+`$ $`\tau _{12}\tau _{14}\tau _{24}\tau _{13}+\tau _{12}\tau _{13}\tau _{23}\tau _{34}+\tau _{12}\tau _{13}\tau _{23}\tau _{14}+\tau _{12}\tau _{13}\tau _{23}\tau _{24}+\tau _{12}\tau _{24}\tau _{23}\tau _{34}`$
$`+`$ $`\tau _{14}\tau _{24}\tau _{23}\tau _{34}+\tau _{13}\tau _{24}\tau _{23}\tau _{34}+\tau _{12}\tau _{14}\tau _{13}\tau _{34}+\tau _{14}\tau _{13}\tau _{23}\tau _{34}+\tau _{14}\tau _{13}\tau _{34}\tau _{24})`$
$``$ $`{\displaystyle \frac{i}{144}}(\tau _{12}\tau _{13}\tau _{24}\tau _{23}\tau _{34}+\tau _{12}\tau _{14}\tau _{24}\tau _{13}\tau _{23}+\tau _{12}\tau _{14}\tau _{24}\tau _{13}\tau _{34}`$
$`+`$ $`\tau _{14}\tau _{13}\tau _{24}\tau _{23}\tau _{34}+\tau _{12}\tau _{14}\tau _{13}\tau _{34}\tau _{23}+\tau _{12}\tau _{14}\tau _{24}\tau _{23}\tau _{34})`$
$``$ $`{\displaystyle \frac{1}{288}}\tau _{12}\tau _{13}\tau _{14}\tau _{23}\tau _{24}\tau _{34}].`$
The HIZ integral is obtained by requiring the symmetry under permutation of the $`t_i`$’s in the final expression for $`I`$, thus from (79) (and a replacement of $`\lambda `$ by $`2\lambda `$),

$$I=C\frac{1}{[\mathrm{\Delta }(t)\mathrm{\Delta }(\lambda )]^3}(f(t_1,\mathrm{},t_k)+\mathrm{perm}.\mathrm{of}f)$$
(A.3)
where the last term means that one adds the terms in which one permutes the $`t_i`$’s for fixed $`\lambda `$’s .