Periodic-orbit theory of universal level correlations in quantum chaos

Using Gutzwiller's semiclassical periodic-orbit theory we demonstrate universal behaviour of the two-point correlator of the density of levels for quantum systems whose classical limit is fully chaotic. We go beyond previous work in establishing the full correlator such that its Fourier transform, the spectral form factor, is determined for all times, below and above the Heisenberg time. We cover dynamics with and without time reversal invariance (from the orthogonal and unitary symmetry classes). A key step in our reasoning is to sum the periodic-orbit expansion in terms of a matrix integral, like the one known from the sigma model of random-matrix theory.


I. INTRODUCTION
Classical chaos is characterized by sensitive dependence on initial conditions, and by the fact that long trajectories uniformly fill the available space. These classical features have profound consequences in quantum mechanics: Many quantum properties of chaotic systems become universal in the semiclassical limit, i.e., they no longer depend on the system in question but only on symmetries. For instance, as conjectured in [1] the statistics of energy eigenvalues is universal: Neighboring levels tend to repel each other, and the correlation function of the level density ρ(E), (where the brackets and overbars denote averages over the center energy E) has the same form for all dynamics belonging to the same symmetry class. The most important classes are those of systems without any symmetries (unitary class; complex Hermitian Hamiltonians) and systems whose sole symmetry is time-reversal invariance (orthogonal class; real symmetric Hamiltonians). If one accepts universal spectral correlations as a fact, a phenomenological prediction for R(ǫ) can be obtained by averaging over ensembles of systems sharing the same symmetries. Modelling the corresponding Hamiltonians through matrices, one arrives at random-matrix theory (RMT), as initially proposed by Wigner and Dyson in the fifties [2] in the context of atomic nuclei. RMT has since found broad applications in many fields of physics. For the correlation function R(ǫ), random-matrix averages yield the prediction [3,4,5] R(ǫ) = −s(ǫ) 2 unitary −s(ǫ) 2 + s ′ (ǫ) 1 π Si(ǫ) − 1 2 sgn(ǫ) orthogonal (2) with s(ǫ) = sin ǫ ǫ and Si(ǫ) = ǫ 0 dǫ ′ s(ǫ ′ ). It is convenient to access R(ǫ) through a complex correlation function C(ǫ) which is analytic in the energy upper half plane and connected with the real correlator as R(ǫ) = lim Im ǫ→+0 Re C(ǫ). The correlator C(ǫ) has the nice property of being retrievable by Borel summation [6] from its asymptotic expansion The foregoing power series contains both non-oscillatory terms and oscillatory terms proportional to e 2iǫ ; the latter are responsible for the singularity (discontinuity at τ = 1 of the first or third derivative in the unitary and orthogonal 2i n for n ≥ 4. This phenomenological RMT result leaves open the question why chaotic systems behave universally, and that question we want to address in the present paper. Taking a semiclassical approach one can follow Gutzwiller [7] and express the level density as a sum over contributions of classical periodic orbits. The correlation function then turns into a sum over pairs of orbits. Systematic contributions to that double sum are due to pairs of orbits whose actions are sufficiently close for the associated quantum amplitudes to interfere constructively. The task of finding correlations in quantum spectra is translated into a classical one, namely to understand the correlations between actions of periodic orbits [8]. The first step in this direction was taken by Berry [9] who showed that pairs of identical or mutually time-reversed orbits explain the leading coefficient c 2 . Higher-order contributions c n with n ≥ 3 are due to a still not widely known feature of chaos, a certain "bunching" of periodic orbits [10]. A long orbit has many close "self-encounters" (see Fig. 1) where it comes close to itself in phase space (possibly up to time reversal, for dynamics with time reversal invariance). Such an orbit is just one of a closely packed "bunch". All orbits in a bunch are nearly identical, except that the orbit pieces inside the encounters are "switched".
Based on insight from disordered systems [11], the study of bunches of orbits differing in encounters was pioneered by Sieber and Richter [12] who derived the next-to-leading coefficient c 3 for time-reversal invariant systems from orbit pairs differing in just one encounter where two stretches of an orbit are close (see Fig. 2). Full agreement with all coefficients c n from RMT was established by the present authors in [13,14].
In none of these works oscillatory contributions could be obtained (Note however courageous forays by Keating [16] and Bogomolny and Keating [17]), due to the fact that Gutzwiller's formula for the level density is divergent. To enforce convergence, one needs to allow for complex energies with imaginary parts large compared to the mean level spacing, Im ǫ ≫ 1. Oscillatory terms proportional to e 2iǫ then become exponentially small and cannot be resolved within the conventional semiclassical approach.
In [18] we proposed a way around this difficulty that we here want to elaborate in detail. The key idea is to represent the correlation function through derivatives of a generating function involving spectral determinants. Two such representations are available and entail different semiclassical periodic-orbit expansions. One of them recovers the non-oscillatory part of the asymptotic expansion (3), essentially in equivalence to [12,13,14]; the other representation breaks new ground by giving the oscillatory part of (3).
The full random-matrix result also, and in fact most naturally, arises within an alternative semiclassical approxima-tion scheme proposed by Berry and Keating in [19]. That scheme constrains the semiclassical periodic-orbit expansion of the spectral determinants det(E − H) to be real and to converge for real energy argument. Inserted into the generating function the resulting "Riemann-Siegel lookalike formula" for det(E − H) was shown in [20] to simultaneously produce both the non-oscillatory and oscillatory parts of the expansion (3). Interestingly, the Riemann-Siegel lookalike formula parallels an exact composition of the random-matrix averaged generating function originating from a Weyl group symmetry [21]. We want to stress that obtaining both the non-oscillatory and the oscillatory contributions to R(ǫ) also implies a derivation of level repulsion: By summing up the asymptotic series (3) we access the full R(ǫ) including its behavior for small ǫ, R(ǫ) ∝ |ǫ| β (with β = 2 for the unitary and β = 1 for the orthogonal class), indicating a suppression of small energy differences. To get access to small ǫ the oscillatory contributions (varying on the scale of the mean level spacing) are of crucial importance.
For our semiclassical approach to spectral statistics we have drawn inspiration from a field-theoretical implementation of RMT, the nonlinear sigma model. There, the ensemble averaged generating function is ultimately transformed into an integral over matrices. In a stationary-phase approximation to that integral, two different saddles are found to respectively contribute non-oscillatory and oscillatory terms. The relation between these two contributions bears close resemblance to semiclassics. Moreover the topologically different structures of orbit bunches in semiclassics turn out to be analogous to Feynman diagrams in a perturbative treatment of the sigma model, with encounters corresponding to vertices. That analogy helps to do the sums required in semiclassics and to compact the semiclassical series to a matrix integral like the one in the sigma model. This derivation has to be contrasted with a previous field-theoretical approach to spectral universality, the ballistic sigma model [22,23].
The present paper is organized as follows. We start by explaining the generating-function approach (Section II), and give the pertinent generalization of Berry's diagonal approximation (Section III). We then evaluate the contributions of bunches of periodic orbits to the generating function and the correlation function (Section IV) and sum them up (Section V). Our reasoning culminates in the semiclassical construction of a sigma model lookalike (Section VI). For pedagogical reasons, Sections III-VI are devoted to the unitary symmetry class, the simplest case. The generalizations necessary for time-reversal invariant dynamics are listed in Section VII.
The Appendix A contains a derivation of contraction rules needed for the semiclassical construction of the matrix integral. Further appendices contain (B) a brief introduction to the supersymmetric sigma model, (C, D) two replica versions of the sigma model and their relation to the semiclassical theory, (E) a characterization of orbit topologies in terms of permutations, and (F) the energy average employed in the definition of the correlator and the connection between the real and complex correlation functions.

A. Definition
All approaches to resolve non-oscillatory and oscillatory contributions start with representing the correlation function through derivatives of a generating function. To get there we express the complex correlator C(ǫ) in terms of Green functions (E − H) −1 of the Hamiltonian H, where the superscript + denotes a positive imaginary part iη. The average will be understood as over the real center energy E, with an averaging window large compared to the mean level spacing (see Appendix F). The first Green function in (4) is a retarded one due to the positive imaginary part of the energy argument while the second is an advanced one. A suitable generating function involves four spectral determinants ∆(E) = det(E − H), The four parameters ǫ A , ǫ B , ǫ C , ǫ D denote offsets from the center energy E, made dimensionless by referral to the mean level spacing 1/ρ; they are all taken to include a positive imaginary part iη 1 . It is worth noting the symmetry Z(ǫ A , ǫ B , ǫ C , ǫ D ) = Z(ǫ A , ǫ B , −ǫ D , −ǫ C ) (6) which in the framework of the sigma-model treatment of RMT reflects a Weyl symmetry [21]. The normalization Z(ǫ A , ǫ B , ǫ A , ǫ B ) = 1 also deserves being mentioned. The correlation function can be accessed from Z in two different ways. In the standard "columnwise" representation, we take two derivatives w.r.t. ǫ A and ǫ B , and subsequently identify all energy increments, Here we used the identity det = exp tr ln. The notation ( ) signals that when equating ǫ A , ǫ B , ǫ C and ǫ D and accounting for the signs in (5), we have identified the four energy arguments in (5) in a "columnwise" way. An alternative representation of the correlation function can be obtained if we identify the energy arguments in a "crosswise" way (in view of the signs in (5)), In (9), the ratio of determinants 2 becomes unity only in the limit η → 0, i.e., one can access R(ǫ) as B. Periodic-orbit representation We are now prepared to semiclassically approximate the generating function in terms of classical periodic orbits. Our starting point is Gutzwiller's formula for the trace of the Green function, Herein, the first term involves the average level density ρ given by Weyl's law as ρ = Ω (2πh) 2 (for systems with two degrees of freedom), with Ω the volume of the energy shell. The orbit a is represented by its period T a , action S a , and stability amplitude F a (defined to include the Maslov phase factor but not the period) 3 1 Due to the imposed positivity of the imaginary parts of the energy offsets the determinants involving ǫ A and ǫ C are associated with retarded Green functions. The minus signs in front of ǫ B and ǫ D in (5) let the respective spectral determinants be associated with advanced Green functions. Note that the present sign convention differs from the one in our letter [18]. 2 The above ratio averages to Z(ǫ + , ǫ + , −ǫ − , −ǫ − ) which can also be written as Z(ǫ + , ǫ + , −ǫ − , −ǫ − ) = exp[2πi∆Nη] with ∆Nη = Nη(E + ǫ 2πρ ) − Nη(E − ǫ 2πρ ). Here Nη denotes the level staircase smeared out over a range η, Nη(E) = P k θη(E − E k ) with θη(x) = π −1 Im ln(−x + iη) the similarly smeared step function. The function Z(ǫ + , ǫ + , −ǫ − , −ǫ − ) could be used as a starting point for calculating the oscillatory part of the correlator C(ǫ), and then some similarity with the courageous foray of Bogomolny and Keating [17] mentioned in the Introduction would arise. Indeed, replacing Nη by the sum of its energy average and a fluctuating part, Nη = Nη + N fluct η , and treating N fluct η in a Gaussian approximation we get exp[2πi∆Nη] ∼ e 2iǫ exp[−2π 2 (∆N fluct η ) 2 ]. The second factor contains the number variance in the exponent; the first factor signals that the oscillatory part of the two-point correlator becomes accessible. This rough argument anticipates the diagonal approximation below. 3 We disregard orbits which are multiple repetitions of a shorter orbit since these are exponentially suppressed in the limit of long periods.
As regards orbits incorporating multiple approximate traversals of short orbits, the respective corrections were shown to be immaterial for the non-oscillatory part of the spectral correlator (see [14]) and can similarly be shown to cancel in our four-determinant generating function Z. Interestingly, such orbits do contribute to the averaged product of two spectral determinants [15].
To proceed to a semiclassical approximation of the spectral determinant, we integrate and exponentiate in (11) as The resulting exponent inherits the smoothed number N (E) of eigenvalues below E and an orbit sum from the trace formula (11), The last member of the foregoing chain of equations has the exponentiated orbit sum expanded into a sum over pseudo-orbits, i. e., non-ordered sets of periodic orbits A [19]. In that pseudo-orbit sum, n A is the number of orbits inside A, S A the sum of their actions, and F A the product of their stability amplitudes F a . The sum over A includes the empty set, whose contribution is unity. (Note that in the foregoing expressions we suppressed a proportionality constant due to the lower limit of the E ′ -integral. That factor is irrelevant since it cancels in the generating function. Eq. (12) applies to energies with a positive imaginary part. For energies with a negative imaginary part complex conjugation yields For inverse spectral determinants, the minus sign in the periodic-orbit sum and hence the factor (−1) nA in the pseudo-orbit sum disappear, The four spectral determinants in Z now yield semiclassical asymptotics as a fourfold sum over pseudo-orbits, where all offset variables ǫ A , ǫ B , ǫ C , ǫ D have positive imaginary parts η. Here the superscript 1 signals that (15) is a semiclassical approximation based on the Gutzwiller trace formula which requires for convergence such that terms vanishing in that limit (such as terms of the order e −η ) will be lost. The complement missed by Z (1) becomes important for real energies (a manifestation of the Stokes phenomenon of asymptotic analysis) and will be introduced below.
2πh is the Heisenberg time, we simplify the quantity Z (1) to The foregoing sum over pseudo-orbit quadruplets combines classical chaos with quantum interference due to the phase factor e i∆S/h ≡ e i(SA(E)−SB(E)+SC (E)−SD(E))/h . For most quadruplets the interference is destructive since for h → 0 the phase factor e i∆S/h oscillates rapidly as the energy varies and therefore vanishes after averaging. For constructive interference, the action mismatch ∆S must be small, at most of the order ofh; to within such a quantum mismatch the cumulative action of the pseudo-orbits B ∪ D must coincide with the cumulative action of A ∪ C. Such action matching happens systematically only if the orbits in B ∪ D either coincide with those in A ∪ C (and thus constitute the "diagonal" contributions, see Sect. III) or differ from them only by their connections inside encounters ("off-diagonal" contributions) as in Figs. 1,2. Indeed, then, the phenomenon of orbit bunching is crucial for spectral universality. Now imagine Z (1) calculated from (17) and then analytically continued to real energy offsets, η → 0 + . If instead of Z (1) we had the full generating function we could obtain both the non-oscillatory and the oscillatory contributions to C(ǫ) by taking derivatives w.r.t. ǫ A and ǫ B and identifying the energy increments either in the columnwise way (7) or the crosswise one (9). However, for Z (1) the columnwise procedure (7) cannot yield oscillatory terms since it makes the (Weyl) exponential e i(ǫA+ǫB −ǫC −ǫD)/2 collapse to unity. In contrast, with the crosswise representation (9) the signs of ǫ C and ǫ D are flipped and the Weyl factor turns into e 2iǫ . Hence the crosswise representation recovers only the oscillatory Fourier component proportional to e 2iǫ .
The random-matrix result suggests that to get the full correlation function one has to add both results Here the first and the third summand are due to (7) while the second summand is due to (10) 4 . Of course it is vital to justify this additivity independently from RMT, a task to be attacked below.
The "crosswise" term ∂ 2 Z (1) ∂ǫA∂ǫB (×) differs from the "columnwise" term ∂ 2 Z (1) ∂ǫA∂ǫB ( ) because ǫ C and ǫ D are interchanged and flipped in sign. This interchange and sign flip may as well be performed on the level of the generating function, i.e., we can write Instead of adding the crosswise representation in (18) we can use Z (1) + Z (2) as a generating function, and get R(ǫ) in full through the columnwise (or crosswise) representation alone. Therefore no extra work is needed for the oscillatory part of the two-point correlation function, once the "standard" part Z (1) of the asymptotics of Z is found. The resulting asymptotics of the generating function respects the Weyl symmetry and remains valid when the arguments are real (η = 0). We note that if one wants to evaluate Z (2) using the Gutzwiller trace formula, the sign flip in ǫ C and ǫ D means that the corresponding arguments of Z (2) should be taken with large negative imaginary parts.

C. Riemann-Siegel lookalike
As pointed out by Keating and Müller in [20], the additive structure (21) of the generating function most naturally arises with the help of the Riemann-Siegel lookalike formula [19].
The essential idea is to eliminate a drawback of the semiclassical approximation (12,13) of the spectral determinant, namely its failure to manifestly preserve the reality of the energy eigenvalues, i. e., the unitarity of the quantum dynamics. In the limit η → 0, ∆(E + ) and ∆(E − ) must become real and identical to each other. That property is not manifest in but may be imposed on Eqs. (12,13). By imposing this property, it was shown in [19] that a duality arises between "long" pseudo-orbits (sum of periods of contributing orbits larger than half the Heisenberg time T H ) and "short" ones (cumulative period smaller than T H /2); the overall contribution of long orbits is the complex conjugate of the contribution arising from short orbits. The ensuing Riemann-Siegel lookalike expresses the real spectral determinant ∆(E) as a sum over pseudo-orbits with durations shorter than T H /2, To evaluate the generating function Z one now uses the Riemann-Siegel formula for the two spectral determinants in the numerator. For the determinants in the denominator we are not allowed to use Riemann-Siegel as the non-zero imaginary parts with appropriate signs in ǫ A , ǫ B are crucial for the definition of Z, and we just stick to Eq. (14). Since each determinant in the numerator becomes the sum of two mutually complex conjugate terms, Z has altogether four summands which involve the four phases ±iπ N (E + ǫ C /(2πρ)) − N (E − ǫ D /(2πρ)) and ±i π N (E + ǫ C /(2πρ)) + N (E − ǫ D /(2πρ)) . The two latter contributions oscillate rapidly as functions of E and thus vanish after averaging.
The two surviving additive terms yield the high-energy asymptotics of the generating function Z ∼ Z (1) + Z (2) found in the previous subsection, according to (17,20), except that the sums over C and D are now restricted to pseudo-orbits with durations shorter than T H /2. To make that difference irrelevant it is convenient to evaluate the pseudo-orbit sum for large imaginary parts η of the dimensionless energy offsets in the arguments of the action, S A (E) → S A (E + iη/2πρ), see (16). That artifice makes sure that in the relevant integrals over orbit periods to be encountered below the integrand falls to exponentially small magnitude as the integration variable rises to near half the Heisenberg time. If all four energy offsets ǫ A,B,C,D are thus provided with large positive imaginary parts, the two additive pieces obtained through the Riemann-Siegel lookalike become respectively Z (1) . Complex conjugation of the last two arguments of Z (2) is needed to enable its calculation in terms of periodic orbits according to (17). At the final stage we use analyticity to let η go to zero; the sum of Z (1) and Z (2) yields then the above asymptotic representation (21).
To obtain a quantitative understanding of spectral statistics, we now have to evaluate the quadruple sum over pseudo-orbits in Z (1) , taking into account orbits in B ∪ D that either (diagonal contributions) coincide with or (off-diagonal contributions) at most differ in encounters from the orbits in A ∪ C.

III. DIAGONAL APPROXIMATION
In the diagonal approximation we effectively neglect all correlations between orbits that are not identical. If we represent the generating functions as a sum over A, B, C, D as in (17) this means we should keep only contributions where the orbits in A ∪ C are repeated inside B ∪ D. However we can save labor if we write the four spectral determinants entering our generating function like the second member in (12), i. e., with the help of exponentiated orbit sums. We thus obtain the equivalent representation of the generating function wherein each periodic orbit participates with the factor If we assume that contributions of different periodic orbits are uncorrelated, Z (1) factorizes into a product of energyaveraged orbit factors z a , Inasmuch as the stability coefficient F a ∼ e −λTa/2 is small for long orbits we may expand the exponential in z a and limit ourselves to the first nontrivial term, z a ∼ 1 + |F a | 2 f a AC f a BD ≈ e |Fa| 2 f a AC f a BD ; here terms proportional to e ±i S(E)/h and e ±2i S(E)/h have been eliminated by the energy average. We so obtain This result allows for an intuitive interpretation: If we assume that the non-identical orbits are uncorrelated the exponent in the product of all orbit factors z a , namely Φ = a F a e iSa(E)/h f a AC + F * a e −iSa(E)/h f a BD , effectively becomes a sum of many independent random contributions; due to the central limit theorem Φ must then obey Gaussian statistics such that e Φ = e Φ 2 /2 which just yields (26).
To bring our result to a more explicit final form, we need to invoke ergodicity in the form of the sum rule of Hannay and Ozorio de Almeida [3,24]: To sum a smooth function of the period T over an ensemble of orbits weighed with |F a | 2 , we may as well compute an integral over periods weighed with 1 T , here T 0 is some minimal period starting from which orbits behave ergodically. Employing the sum rule (27) we get Since for T ≪ T H the integrand is ∝ (T /T H ) 2 we can set the lower integration limit to 0, accepting an error of the order (T 0 /T H ) 2 . Then using ∞ 0 dx x e iax − e ibx = ln b a we arrive at the diagonal part of the generating function The diagonal approximation yields an estimate for the correlation function C(ǫ). Taking into account both Z (1) diag and its counterpart Z This result already reproduces the full RMT prediction for the unitary symmetry class under discussion. It remains to explain, within our semiclassical approach, why all off-diagonal contributions do in fact vanish.

IV. OFF-DIAGONAL CONTRIBUTIONS
Off-diagonal contributions arise from quadruplets of pseudo-orbits where not all orbits of A ∪ C are simply repeated inside B ∪ D; some orbits of A ∪ C and B ∪ D differ in encounters. By studying these encounters in detail we shall show that for systems without time-reversal invariance the off-diagonal contributions cancel mutually in a nontrivial way. This also sets the stage for a generalization to time-reversal invariant dynamics where this cancellation does not occur, and off-diagonal quadruplets of pseudo-orbits give rise to terms of higher order in 1/ǫ.

A. Encounters and orbit bunches
Long periodic orbits often come close to themselves, and one another, in phase space. We speak of an l-encounter whenever l stretches of one or more periodic orbits in a chaotic system approach in this way. Each stretch begins and ends at phase-space points to be called "entrance ports" and "exit ports". These ports can be uniquely defined by setting an upper bound for the separation between the encountering stretches. The parts of an orbit in between encounter stretches will be called "links". Now the crucial point is that the orbits and pseudo-orbits in hyperbolic motion exist in bunches. The members of a bunch differ noticeably only in encounters where the ports are connected differently by encounter stretches. In contrast the links are very nearly identical. All members of such a bunch have almost identical actions. If two of them are taken as A ∪ C and B ∪ D in Eq. (17) their contributions interfere constructively and need to be taken into account for a full understanding of spectral statistics.
The phenomenon of bunching is a consequence of hyperbolicity. For a system with two degrees of freedom 5 , hyperbolicity means that deflections between phase-space points in the three-dimensional energy shell can be split into three components: one pointing along the flow, one unstable component u, and one stable component s. As two close-by phase-space points move, the unstable component of their separation grows exponentially whereas the stable component shrinks exponentially. In the limit of large times the rate of growth or shrinking is given by the Lyapunov exponent λ as As a consequence, given any two encounter stretches that are long compared with the characteristic size of the system and close enough to linearize the motion as in (30) we can draw a piece of a physical trajectory that starts close to the first stretch and ends close to the second stretch. Its separation from the first stretch must point along the unstable direction at the respective phase-space points (such that it approaches the first stretch for large negative times) whereas the separation from the second stretch must point along the stable direction (such that it approaches that stretch for large positive times). The new piece of trajectory is thus located where the unstable manifold of one encounter stretch of A ∪ C intersects the stable manifold of the other stretch [26]. As a consequence we can reconnect the ports of each encounter in any way we like while leaving the links in between the encounter stretches practically unchanged. (Note that subsequent to the construction just explained exponentially small corrections to the encounter stretches and links are needed to connect them to continuous orbits.)

B. Structures of pseudo-orbit quadruplets
The pseudo-orbit quadruplets can have many different topologies some of which are shown in Fig. 3. For their classification we introduce the notion of a structure by which we mean a diagram like the ones in Fig. 3 but with the encounters and their ports numbered and each periodic orbit assigned to a particular pseudo-orbit. All structures can be generated if we proceed as follows: S1: Encounters -We consider any number V of encounters. For each of the encounters σ = 1, 2, . . . , V , we take into account all possible numbers of stretches l(σ) = 2, 3, 4, . . .. The l(σ) entrance and exit ports of each encounter are numbered from 1 to l(σ) in such a way that the i-th entrance port is connected to the i-th exit port before reconnection, and to the (i − 1)-st exit port after reconnection; the first entrance reconnects to the last exit. The encounters can thus be depicted as in Fig. 4 where full and dashed lines denote the encounter stretches before and after the reconnection.
S2: Links -We connect the ports by links in all possible ways. Each link must lead from an exit port to an entrance port of either the same or a different encounter.
S3: Orbits -We include all possible distributions of the orbits among the pseudo-orbits A, B, C, D. The graph formed by all links and encounter stretches before reconnection falls into a certain number of disjoint parts each standing for a periodic orbit. Each of these can be included in either A or C. After changing connections inside the encounters the graph consists of a (possibly different) number of disjoint orbits which have to be distributed between B and D.
Thus equipped we shall split the sum over quadruplets in (17) into a sum over structures and a sum over quadruplets associated to each structure. The sum over quadruplets associated to each structure will be done using ergodicity. We shall determine the likelihood for having encounters inside and between periodic orbits; this likelihood depends on the length of the links and the phase-space separation between the encounter stretches. The sum can then be replaced by suitable integrals. The subsequent summation over structures becomes a purely combinatorial problem.
When doing the sums and integrals mentioned, we have to be aware of a certain overcounting. Since our definition of a structure involves an ordering of encounters and ports, each physical quadruplet can be associated with a structure in several different ways. To begin with, the encounters of a given quadruplet can be numbered in any one of V ! different orders. Then, in each encounter any of its l(σ) entrance ports can be chosen as the first one. Note that as soon as the first entrance in an encounter is chosen, the numbers of all other encounter ports of the quadruplet become uniquely defined by our numbering scheme. In all, we get V ! σ l(σ) alternatives either leading to different structures or to the same structure but with different numberings and thus to a different set of the integration parameters. In both cases after summation over structures and integration over the continuous parameters the same quadruplet is taken into account several times. In order to compensate this overcounting, the contribution of each structure needs to be divided by V ! σ l(σ). Details can be found in Appendix E where a formal definition of structures in terms of permutations is given.

C. Classical properties of orbit bunches
We now follow [13,14] to investigate the classical properties associated to each structure, and then evaluate their contribution to the generating function. We start with the example of orbits differing in a 2-encounter (see Fig. 5). First we need coordinates measuring the phase-space separations between the two encounter stretches. We therefore place a transversal Poincaré surface of section at an arbitrary position inside the encounter and consider the two points where the encounter stretches pierce through that section. The separation of the piercings can be decomposed into a stable component s and an unstable component u [26,27,28]. These coordinates determine the difference between the action of the original orbit and the cumulative action of the two partner orbits as ∆S = su (see also [12,29]). Moreover, if we define the encounter as the region with phase-space separations |s|, |u| < c (where c is arbitrary but classically small) the time until the end of the encounter can be estimated as 1 λ ln c |u| and the time since the beginning of the encounter is obtained as 1 λ ln c |s| . Summation thus yields the overall duration 6 of the encounter t enc ∼ 1 λ ln c 2 |su| [14]. For the semiclassically relevant encounters we need ∆S = su ∼h, i.e., the encounter duration is of the order of the Ehrenfest time T E = 1 λ ln c 2 h . Encounters are thus considerably shorter than the orbits, whose periods must be of the order of the Heisenberg time T H ∝h −1 .
The likelihood of finding encounters in a periodic orbit can be determined using ergodicity, i.e., the fact that long trajectories and (ensembles of) long periodic orbits uniformly fill the energy shell [30]. Ergodicity implies that the probability for finding two piercings through a Poincaré section with separations t, s, u in differential intervals dt, ds, du is uniform and reads dtdsdu Ω , with Ω the volume of the energy shell. Statistical independence of the two piercings is assumed here; that assumption is an allowable one since for orbits with a period of the order T H the links between the encounter stretches typically have durations much longer than the classical timescale needed for the decay of correlations. Now integration of the probability density 1/Ω yields the number density w T (s, u) for finding, inside an orbit of period T , 2-encounters with stable and unstable separations s and u. That density was determined in [14,31] as Here the factor T is the period of the orbit (originating from integration over the time of the first piercing). dt is an integral over the time of the second piercing reckonned from the first one, with the following restriction: The time separation between the two piercings must be at least t enc so that in between the piercings the orbit can go through the rest of the first encounter stretch, a link of positive duration, and the beginning of the second stretch. This reasoning leads to dt = T − 2t enc , but we prefer to just write dt. Now let us consider more complicated structures with, say, V different encounters inside or between, say, n orbits of periods T 1 , T 2 , . . . , T n . The encounters are labelled by an index σ running from 1 to V and the σth encounter may have l(σ) stretches. Drawing a Poincaré section in each encounter we face altogether L = σ l(σ) piercings, l(σ) ones for the σth encounter. The separations between the piercings belonging to one encounter are characterized by l(σ) − 1 stable and unstable coordinates s σm , u σm with m = 1, 2, . . . , l(σ) − 1 (see [14] for the precise definition). These coordinates determine the action difference as ∆S = σ,m s σm u σm and the duration of the σ-th encounter as [14]. By following the same steps as for (31) the number density of encounters is generalized to with s = {s σm }, u = {u σm }. Here the orbit periods T i arise from integrating over the time of the first piercing in every orbit, and the integral d L−n t goes over the possible times of the remaining piercings reckoned from the first one.

D. Contribution of orbit bunches to Z
We can now investigate the impact orbit bunches have on spectral statistics. We thus evaluate the generating function Z (1) as given by the fourfold pseudo-orbit sum in (17), taking into account quadruplets with A ∪ C and B ∪ D differing in any number of encounters. As before A ∪ C and B ∪ D may also coincide in any number of component orbits, and therefore Z (1) diag still arises as a factor. We can thus rewrite (17) as with Z off picking up non-empty pseudo-orbits A ∪ C and B ∪ D that differ in encounters while no longer comprising identically repeated orbits, The off-diagonal part Z off of the generating function can now be written as a sum over structures, and then over all quadruplets of pseudo-orbits pertaining to each structure. We thus obtain where as explained above we divided out the number of possible orderings of encounters and stretches, V ! σ l(σ), since otherwise each quadruplet would be counted multiple times. The sum over pseudo-orbits permitted by a structure can be done by summing over all choices for the n component orbits p 1 , p 2 , . . . , p n of A ∪ C and integrating over the L − V = σ (l(σ) − 1) pairs of stable and unstable coordinates characterizing the V encounters, with the number density (32) as weight, × w T1,T2,...Tn (s, u) e i P σ,m sσmuσm/h e i(TAǫA+TB ǫB +TC ǫC +TD ǫD)/TH .
Each of the n orbit sums can now be done using the sum rule of Hannay and Ozorio de Almeida, pi |F pi | 2 → ∞ T0 dT T . Importing the density w from (32) we see the orbit periods in w cancel against those brought in by the sum rule and face an L-fold time integral, over the n periods and L − n piercing times.
To simplify the foregoing expression for Z (1) off we change the L-fold time integral just mentioned to an integral over the L link durations. The Jacobian of that transformation is unity. We thus obtain where we inserted compensating powers of the Heisenberg time. The summand coming from a structure in (37) then becomes a product over all links times a product over all encounters. It is easy to evaluate the pertinent factors. Each link appears precisely once in A ∪ C as well as in B ∪ D. The integral over its duration picks up the integrand 1 TH e i(ǫA or C +ǫB or D )t link /TH and contributes i ǫ A or C + ǫ B or D link factor ; (38) clearly, the index alternatives have to be decided according to which pseudo-orbit the link pertains to before reconnection (A or C) and after reconnection (B or D).
We now turn to the factor coming from an l-encounter. In the original pseudo-orbit A ∪ C, each of the l stretches may belong to either A or C. We denote by l A the number of stretches belonging to A and by l C the number of stretches belonging to C, with l A + l C = l. When switching connections inside the encounter these stretches are replaced by l new stretches differently connecting the same ports as the old stretches. These new stretches belong to either B or D, and the corresponding numbers are denoted by l B and l D with l B + l D = l. The encounter now yields a contribution l A t enc to the duration of pseudo-orbit A; analogously for B, C and D. Hence the encounter draws the part e i(lAǫA+lB ǫB +lC ǫC +lD ǫD)tenc/TH from the phase factor e i(TAǫA+TB ǫB +TC ǫC +TD ǫD)/TH in (35). Collecting everything from (37) pertaining to the encounter under consideration we get the dimensionless term For the calculation of the remaining integral we refer to [14,32] and here but briefly sketch the argument. The key is to Taylor expand the second exponential and show that only the term linear in t enc survives in the limith → 0. In that surviving term the encounter duration t enc cancels and the factor i(l A ǫ A + l B ǫ B + l C ǫ C + l D ǫ D )/T H can be pulled out from the integral. Using T H = Ω 2πh and TH The sum over structures (37) thus simplifies to A further simplification can be achieved after inspecting more closely the encounter contributions. For an l-encounter of some structure we include the overcounting factor 1 l and face the contribution 1 . Each of the l structures obtained from the given one by cyclic renumbering of the stretches of the selected encounter contributes identically, and altogether these l structures contribute i( . We now claim that the same result is obtained if we choose to assign to each structure the encounter factor i(ǫ A or C + ǫ B or D ) where the index alternatives must be decided according to the placement of the first entrance port; e.g., for a structure where the entrance port of the first encounter stretch belongs to the pseudo-orbit A before reconnection and to D afterwards we would write i(ǫ A + ǫ D ). Taking a sum over our group of l equivalent structures we then obtain ǫ A or C = l A ǫ A + l C ǫ C since the first entrance belongs to A and C in l A and l C members of the group, respectively. Similarly we have can be replaced by i(ǫ A or C + ǫ B or D ) since after summation over the group of equivalent structures associated with the encounter the two expressions give the same result. Eq. (41) can thus be replaced by the simpler formula in which the numerator refers to the first entrance port of each encounter, A symmetry between link and encounter factors thus results and entails an interesting property of the foregoing "master formula" for Z (1) off : Apart from a factor −1, each encounter factor is cancelled by the link factor for the link ending at the first entrance port of that encounter. Indeed, the first entrance port and the link ending there belong to the same pseudo-orbit before the reconnection, either A or C, and likewise after the reconnection, either B or D. After that cancellation, the numerator of the ǫ-dependent fraction becomes (−1) V and only the denominator remains ǫ-dependent. That property will be important for proving the cancellation of all encounter contributions in Sect. V.
According to (42) the contribution of a structure to the generating function is of the order 1 ǫ L−V . The sum over structures therefore generates a power series in 1 ǫ . Low-order diagrams can easily be drawn; for instance, Fig. 3 depicts all diagrams contributing in the leading two orders, 1 ǫ and 1 ǫ 2 , both for the unitary and orthogonal class (apart from a further diagram containing two copies of Fig. 5 and a diagram containing one copy of Fig. 5 and one Sieber-Richter pair).

E. Example: Structures for single 2-encounter
It is well to illustrate the formula (42) for the simple example of a single 2-encounter with two entrance and exit ports, depicted in Fig. 6 (a); the full and dashed lines respectively show the encounter stretches before and after reconnection. We can connect the ports by links (stage S2) in L! = 2 ways; the two resulting diagrams are shown in parts (b) and (c) of Fig. 6 .
We first deal with Fig. 6 (b) where we face a single pre-reconnection orbit γ 0 incorporating all ports and links. After the reconnection the graph disconnects into two orbits γ 1 and γ 2 each of which contains one entrance and one exit port as well as one link; the labels 1 and 2 signal the number of the respective entrance port.
At the stage S3 we generate structures by dividing the orbits among the pseudo-orbits. Since we deal with three orbits there are 2 3 distinct such allotments; four of these are and four similar structures have A interchanged with C. Consider the first structure in the foregoing list. The entrance port 1 belongs to γ 0 ∈ A before, and to γ 1 ∈ B after reconnection. Therefore the encounter factor in the numerator reads i (ǫ A + ǫ B ). Both links belong to A before and B after reconnection; hence the denominator is The factors in (42) are calculated using V = 1, n C = n D = 0; the structure thus yields Z 1 = i ǫA+ǫB . The contribution of the second structure differs by the replacement B → D; there are two orbits in D and none in C such that the sign is unchanged and we get Z 2 = i ǫA+ǫD . In the third structure the first entrance port belongs to A before and B after reconnection. The two links are in A before reconnection; after reconnection one of them goes to B, another one to D. The sign factor is negative since n C = 0, n D = 1. Hence we get Obviously, Z 3 cancels with Z 2 and Z 4 with Z 1 . The reader is invited to repeat the calculations for the remaining four structures and again find full cancellation. Finally, a patient reader will check that Fig. 6 (c) leads to 8 more structures whose contributions coincide with those of Fig. 6 (b) and also cancel pairwise.
We now show that for systems without time-reversal invariance all off-diagonal contributions mutually cancel. Our proof is based on two cancellation mechanisms: Each structure is cancelled either by a structure in which the orbits are partitioned differently into pseudo-orbits or in which one link has been removed. Different partition: To start we consider structures where one of orbits (say, in A) only contains one encounter stretch with two ports and a link, as in Fig. 7, and the stretch is the first in its encounter. The ports and the link then also have to belong to the same partner orbit (included, say, in B). Then according to Eq. (42) the contribution of the structure contains a factor i(ǫ A + ǫ B ) from the entrance port of the stretch and a factor 1 −i(ǫA+ǫB) from the link; these link and stretch factors multiply to −1.
Now we compare with a structure where the orbit from Fig. 7 is included in C rather than A. Then the above factors are replaced by i(ǫ C + ǫ B ) and 1 −i(ǫC +ǫB) , and the product is again −1. However, there is an extra factor −1 because the number of orbits in C is increased by one. Hence the two structures considered have cancelling contributions.
The same argument goes through if B is replaced by D, or if the orbit in Fig. 7 is a partner orbit; in the latter case its inclusion in B and D yields mutually cancelling contributions. In all cases we are free to drop from Eq. (42) all structures where an orbit includes just one link and one encounter stretch which is the first in its encounter. Shrunken links: To get rid of the remaining contributions we consider an arbitrary structure involving a 2encounter, and there focus specifically on the link ending at the first entrance port of the 2-encounter. Then where does that link start? The only possibility that we still have to consider is that it starts at the exit port of a different encounter (with an arbitrary number l of stretches). All other possibilities have been dealt with above: if the link would start at the exit port of the same stretch the corresponding original orbit would just include one stretch and one link; if the link would start at the exit of the second stretch, there would be a partner orbit with only one stretch and one link. Hence we must have a situation as depicted schematically in Fig. 8a for l = 4. Now let us consider a structure that looks topologically the same as Fig. 8a but has the link just considered shrunk away. In that structure the two ports connected by the link in question (the black ports in the picture) are ignored, leaving altogether l + 1 entrance ports and l + 1 exit ports. These ports have to be connected as in Fig. 8a, with connections in the original orbits following the full lines and connections in the partner orbits following the dashed lines. If we leave out the black ports, this leads to connections as in Fig. 8b. We conclude that when the link between the 2-encounter and the l-encounter is removed these two encounters merge into a single (l + 1)-encounter. The contribution of the new structure is proportional to the old one. Just the factor − 1 i(ǫA or C +ǫB or D ) from the link and the corresponding factor i(ǫ A or C + ǫ B or D ) from the removed entrance port are gone, and a change of sign is again incurred. Hence the two structures depicted have a chance to cancel.
What remains to be done is a little bit of book-keeping (e.g., to account for the numbering of encounters, and for the possibility that removing links from different structures yields the same result). To do so, we single out one stretch in each structure as a reference; that stretch must not be the first stretch in its encounter. The number of eligible stretches (i.e. the number of all stretches L minus the number of encounters V ) will be denoted by n. We now sum over all structures and over all possibilities of choosing references stretches inside a structure; to avoid an obvious overcounting we modify the summand of Eq. (42) to include an additional divisor n. Now for each contribution where the reference stretch belongs to a 2-encounter we can find exactly one contribution with opposite sign where the reference stretch belongs to an l-encounter with l > 2. This is done by (i) merging the 2-encounter with a different encounter as described above, (ii) renumbering the encounters to make up for the loss of one encounter (i.e. decrease by 1 the labels of all encounters which are larger than the label of the removed encounter), and (iii) choosing the encounter stretch following the merged one (e.g., the fourth encounter stretch in Fig. 8b) as our new reference stretch. The last step guarantees that the reference stretch ends up in any position but the first in its encounter, as demanded above. It is clear that the new contribution has the same n as the initial one (as both L and V have been reduced by one), hence the additional divisor 1/n weighing every contribution remains unchanged.
Each contribution where the reference stretch belongs to an l-encounter can be accessed by following this procedure. This becomes clear if we take a contribution as in Fig. 8b and then undo the steps described above to arrive at Fig.  8a. The same contribution as in Fig. 8b results from V different contributions as in Fig. 8a since the removed 2-encounter may have an arbitrary index between 1 and V . However this degeneracy is compensated by the change of the factorial in the denominator of Eq. (42). Since the remaining factors are the same but the sign is different, we thus see that every contribution where the reference stretch belongs to an l-encounter exactly cancels with V "parent" contributions with the reference stretch in the 2-encounter.
To summarize, we have seen that all off-diagonal contributions cancel due to two mechanisms: For quadruplets where at least one orbit is limited to just one encounter stretch (which is the first in its encounter) and one link, the structures where this orbit is assigned to different pseudo-orbits mutually cancel. For quadruplets where no orbit is of this type we singled out a reference stretch. We then distinguished between contributions where this reference stretch is part of a 2-encounter or an l-encounter with l > 2. By the construction of Fig. 8 we were able to show that these two types of contributions mutually cancel. Hence for systems without time-reversal invariance we have Z (1) off = 0 and only the diagonal term remains.

VI. SEMICLASSICAL CONSTRUCTION OF A SIGMA MODEL
An alternative proof for the cancellation of all off-diagonal contributions can be based on field-theoretical techniques. We shall show that the whole series (42) can be summed to become an integral over a certain matrix manifold. That matrix integral turns out to coincide with the one known from the σ-model in RMT, and gives the same result as the diagonal approximation.

A. Matrix elements for ports and contraction lines for links
We recall that structures of pseudo-orbit quadruplets were defined by (i) fixing the numbers of encounters and stretches involved (the stretches being ordered according to the way they are connected in A ∪ C and B ∪ D), (ii) connecting the encounter ports by links, and (iii) distributing the resulting pre-reconnection orbits over the pseudoorbits A and C, and the post-reconnection orbits between B and D.
Now we encode each structure in (42) by a sequence of alternating symbols B kj andB jk respectively standing for entrance and exit ports. We start with the first entrance port of the first encounter, continue with the first exit port of the first encounter, then list all further ports of that encounter, then of all further encounters. The indices of B,B indicate the affiliation to the pseudo-orbits. Namely, the second index of B and the first index ofB determine whether before reconnection the port belongs to A (case j = 1) or to C (j = 2). Similarly k shows whether the port belongs after reconnection to B (if k = 1) or to D (k = 2). Due to the two possible values for both k and j we altogether confront four B kj 's and fourB jk 's.
Our conventions for the numbering of ports imply that the ports corresponding to a B and the subsequentB (entrance and exit ports with the same number) are connected by an encounter stretch of an original orbit. Hence they must belong to the same original pseudo-orbit and have the same index j, as in B k ′ jBjk . Similarly the ports associated to aB and the subsequent B (representing an exit port and the entrance port with the following number) must be connected by an encounter stretch of a partner orbit. Consequently they belong to the same partner pseudoorbit and share the same index k, as inB jk B kj ′′ . The same applies for the lastB and the first B in each encounter. In any case neighboring indices of successive port symbols must coincide. For a 2-encounter with two entrance and two exit ports as in Fig. 6 (a) we thus get the four-symbol sequence just like in a trace of a matrix product; note that the sub-indices of the indices j, k signal port or encounter stretch labels within the encounter. For a structure with arbitrarily many encounters indexed by σ = 1, 2, . . . V and stretches inside each encounter indexed by i = 1, 2, . . . l(σ) the sequence of ports gives rise to the product Note that we do not imply summation over repeated indices either here or anywhere else in the paper. The indices j and k now carry two subindices each, σ = 1, 2, . . . , V to label the encounter and i = 1, 2, . . . , l(σ) to label the stretches within an encounter. Such "proliferation" of subindices notwithstanding there are still only four different B kj 's and fourB jk 's. Whenever possible below we shall suppress the "address labels" on k and j. If we take B kj ,B jk to be and similarly forB we may consider the expression above as one of the summands in the expansion of the trace product V σ=1 tr BB l(σ) .
The links of a structure will be indicated by drawing "contraction lines" between the respective exit and entrance ports. Each such line connects one entrance and one exit port of some orbit and its partner (recall that links do not change when shifting connections in encounters). The indices j, k of B andB indicating accommodation in original and partner pseudo-orbits therefore have to coincide.
Example: In our above example of a product of four symbols the ports can be connected by links in two ways: note the agreement of indices of the connected ports. We can even complement these diagrams to topological pictures of the orbits. If we draw lower lines to indicate how the ports are connected inside the encounter we obtain Here the upper pair represents the structures in the list (43) and pertains to Fig. 6. The first picture shows the single orbit before reconnection (with encounter stretches connecting the ports as BB). The second picture depicts the two orbits after reconnection (with encounter stretches connecting the ports asBB). Analogously the lower pair represents the structures pertaining to Fig. 6 (c).
Diagrammatic rules: We can now translate the previous diagrammatic rules for the contribution of a structure in (42) to the new language of contracted sequences of matrix elements: (i) The ports B kj andB j ′ k ′ connected by a link must have the same subscripts j = j ′ and k = k ′ , and the contribution from that link in (42) We can thus write the factor provided by a contraction between B kj ,B j ′ k ′ as [−i(ǫ j + ǫ ′ k )] −1 δ kk ′ δ jj ′ where ǫ j stands for ǫ A (j = 1) or ǫ C (j = 2); similarly ǫ ′ k equals ǫ B for k = 1 and ǫ D for k = 2. (The primed energy will always refer to the post-reconnection pseudo-orbits B or D.) (ii) The indices in an encounter factor i (ǫ A or C + ǫ B or D ) in (42) are determined by the pseudo-orbits containing the first entrance port of the encounter. In (45) the first entrance port of the σth encounter is denoted by B kσ1,jσ1 such that the respective encounter factor can be rewritten as i(ǫ jσ1 + ǫ ′ kσ1 ). (iii) Finally Eq. (42) implies that there should be a factor (−1) nC +nD depending on the numbers of orbits included in the pseudo-orbits C and D.

B. Wick's theorem and link summation
We now want to sum over all structures, which in particular involves summation over all possible ways to draw links alias all contractions. A useful tool for that summation is provided by Wick's theorem. Consider the integral over the complex plane with dbdb * ≡ dReb dImb π and f (b, b * ) a product involving an equal number of b's and b * 's, and Imǫ > 0 to ensure convergence. Wick's theorem states (i) that such Gaussian integrals can be written as sums over "diagrams" where each b in f (b, b * ) is connected to one b * by a "contraction line", and (ii) that for complex b, b * these diagrams can be evaluated by performing the Gaussian integral over contracted elements as if the rest of the integrand were absent. The contraction lines can then be eliminated as in Using this contraction rule one can stepwise remove contraction lines and the b's and b * 's connected by these lines and thereby in each step earn a factor − 1 iǫ . For complex b and b * Wick's theorem is just a fancy way of writing down the integral dbdb * e iǫbb * (bb * ) n = n! (−iǫ) n+1 ; the factorial gives the number of possible ways to draw contraction lines. But remarkably, Wick's theorem also holds, up to a sign in Eq. (49) , if we replace b and b * by Grassmann variables η and η * . Grassmann variables are formally defined as anticommuting, ηη * = −η * η, and with integrals dη = dη * = 0, dηη = dη * η * = 1 [3]. While η and η * should in general be seen as independent variables, for the purposes of this paper it is convenient to imagine η * to be the complex conjugate of η and define the complex conjugate of η * to be (η * ) * = −η. Imagining exp (iǫηη * ) Taylor-expanded we obtain This is the analogue of (49) with g = 1 which is the only case of interest since all powers of the Grassmann variables higher than 1 vanish. Below we shall apply contraction rules of the foregoing types to integrals over an arbitrary number of pairs of commuting and anti-commuting variables. Following traditions of the physical literature we shall refer to the anti-commuting variables as "Fermionic" as opposed to the commuting "Bosonic" ones. We now want to use Wick's theorem to represent sums over all possible ways of drawing contraction lines between B kj 's andB jk 's, for general products of the form (45). This can be done if we declare these symbols to be either complex Bosonic or Fermionic variables and make the variables B kj andB jk between which we want to draw contraction lines mutually conjugate up to a sign, i.e., we setB jk ≡ ±B * kj . (The choice of the Bosonic or Fermionic character of the variables and of the sign will be made later.) Then all possible contractions are generated by taking the corresponding expression without contraction lines, and integrating with a Gaussian weight that has a term proportional to B kj B * kj in the exponent. The prefactor of B kj B * kj should be chosen as i(ǫ j + ǫ ′ k ). Then Wick's theorem yields the desired factor [−i(ǫ j + ǫ ′ k )] −1 for each link belonging to the pseudo-orbits j and k (up to a sign factor arising for Fermionic variables). If we also include the encounter factor i(ǫ jσ1 + ǫ ′ kσ1 ) we are led to expressions of the type where the integration measure involves the product of all independent differentials, All (convergent) expressions of this type contain all relevant structures with the appropriate link and encounter factors.

C. Signs
The most delicate task is to fix the signs in (51) and the Bosonic vs. Fermionic nature of the variables so as to produce the sign factor (−1) nC +nD in the sum over structures (42). We shall make the following choices which are inspired by the supersymmetric sigma model of RMT (see [3] and Appendix B ) but will here be justified purely on semiclassical grounds: We take B 11  . Then we pick the signs in the Gaussian and the product of B's andB's of (51) such that these terms can be written in terms of supertraces Str X ≡ tr σ z X.
We also add an overall minus sign.
For the Gaussian integral we thus write here the curly brackets {. . .} have the same content as in (51); to compact the notation we employ the matriceŝ moreover, in the last member of the foregoing chain of equations we abbreviate the B-integral with the Gaussian weight (including the negative sign factor) by double angular brackets . . . . Since the Gaussian is not normalized, . . . is not an average, i.e., 1 = 1. Instead we have reminiscent of the diagonal part Z diag of Z (1) . The factor 1 remains after all contractions are performed and thus all B,B removed. Reassuringly it accounts, once multiplied with e i(ǫA+ǫB −ǫC−ǫD )/2 , for all diagonal quadruplets.
When writing the product of B's andB's in (51) in terms of a supertrace, we still need to keep the indices fixed (note that no summation over port labels is implied!). We thus employ projection matrices P 1 = diag(1, 0), P 2 = diag(0, 1) and write (51) as Note that just like for a trace, it is possible to cyclically permute the supermatrices appearing in this supertrace. If we also install the weight 1 V ! of a structure and the Weyl factor, we obtain which gives the cumulative contribution to Z (1) = Z off ) of all structures with V encounters, each having a fixed number of stretches l(σ), σ = 1 . . . V and a fixed distribution of ports {j σi , k σi } among the pseudo-orbits; each allowed way of drawing contraction lines makes for one structure. The purely diagonal additive term is due to V = 0.
The only thing still to be shown for a proof of this statement is that our conventions incorporate the correct sign (−1) nC +nD for each structure depending on the number of orbits included in C and D. To fill this gap we reexpress the contraction rules (49) and (50) in terms of B andB. Let us consider one way of drawing contraction lines in Eq. (58), and moreover select a particular contraction, either between ports within the same encounters (inside the same supertrace) or from different encounters. Then we can prove two identities connecting the contribution of a structure with that of another structure containing two ports and one link less (see A for details of the derivation), In these formulae, contraction lines other than the selected one are not drawn and are assumed to coincide in the left and right hand sides. In terms of periodic orbits, the rules (59, 60) respectively account for the two ports connected by a link belonging to either different encounters or to the same one. The symbols X, Y, U, V represent matrix products as in (58)  are consistent with the diagrammatic rules stated above: Links can only connect ports associated to the same original pseudo-orbit j = j ′ and the same partner pseudo-orbit k = k ′ , and contribute factors − 1 i(ǫj +ǫ ′ k ) . The contraction rules can be used to stepwise simplify our expressions for the structure contributions. Each step removes two matrices or, in other words, a link connecting two ports from a periodic-orbit structure, without changing other links.
At first sight, the rules (59,60) do not seem to yield any sign factors. But we must pay special attention to steps where the contraction line pertaining to the last link of an orbit is removed, after removal of all other contraction lines of some "parent" orbit in previous steps. The B andB so connected represent the only ports of the remnant orbit. That orbit is still periodic, and therefore the two ports are connected not only by the link but also by an encounter stretch. The two matrices B's andB's then have just one projection matrix in between. If the orbit in question is an original one, the matrices follow each other like BP jB . We thus have to use contraction rule (60) with U = 1. Then the first supertrace on the r.h.s. turns into Str P j which equals 1 if j = 1 (the orbit belongs to pseudo-orbit A) and −1 if j = 2 (the orbit belongs to C). As desired a sign factor −1 thus arises for each orbit included in C.
An analogous result holds for partner orbits. Imagine that we remove the last contraction line associated to a partner orbit. Then the ports connected by this line must also be connected by a partner encounter stretch. The rule (60) then applies with V = 1. The second supertrace on the r.h.s. thus equals Str P k and yields −1 if k = 2, i.e., if the orbit belongs to D.
To conclude, by turning B andB into supermatrices, we have successfully incorporated the sign factor (−1) nC +nD . Other ways to deal with that sign factor (the replica trick) are discussed in Appendices C and D.

D. Towards a sigma model
To determine Z (1) we sum over all structures of pseudo-orbit quadruplets. In Eq. (58) we had already collected all ways of connecting ports with links. Now we sum over all possibilities of assigning ports to pseudo-orbits, i. e., over the indices j σi , k σi in (58). This means that we have to drop the projection matrices. To keep the factors i(ǫ σ1 + ǫ ′ σ1 ) connected to the indices in the supertrace we use the diagonal offset matricesǫ,ǫ ′ and write Here the term Strǫ ′ (BB) l(σ) contains the factor ǫ ′ kσ1 from (58) while Strǫ(BB) l(σ) includes ǫ jσ1 . It remains to sum over the number V of encounters and over their sizes l(σ), σ = 1, 2, . . . , V . We let V run from 0 to ∞, the term V = 0 taking care of the purely diagonal contribution, and obtain Now recall that (62) was derived using a semiclassical approximation, which converges only with large imaginary parts η for all energy offsets. In the limit η ≫ 1 the B-integral draws overwhelmingly dominant contributions from near B = 0 whereas other B only make for exponentially small corrections of order e −η . Therefore in the important integration region we have ||BB|| ≪ 1 and the geometric series in the foregoing exponent can be replaced by its sum. Upon explicitly writing the Gaussian integral we get the compact result The integration domain L must include the vicinity of the point B =B = 0 and obey 1 − BB > 0, to exclude regions where the integrand is exponentially large; otherwise it is arbitrary. Under these limitations and η ≫ 1 the asymptotic 1/ǫ-expansion of the integral (63) coinciding with the semiclassical series (61) is created by the saddle point of the exponent at B =B = 0. If we evaluate (63) in the limit η ≫ 1 (see Appendix B ) we obtain the same result Z (1) as in the diagonal approximation, Eq (28). Then the same Z (2) as in the diagonal approximation is obtained using the Riemann-Siegel relation Z (2) . We see once again that all off-diagonal terms in the generating function cancel.

E. Comparison with the RMT sigma model
In the so called rational parametrization of the zero-dimensional sigma model of RMT (see, e.g., Appendix B ) the exact generating function Z is given by the integral identical in appearance with (63) and with the integration measure and the structure of the supermatrices the same as there. However there are two important differences: First, the offset variables ǫ A,B,C,D are all allowed to be real. Second, the complex Bosonic integration variables are assigned specific integration domains, the open unit disk |B 11 | < 1 and the full complex plane |B 22 | < ∞. It is easy to show that one of the eigenvalues of BB is then between 0 and 1 while the other one can take any negative value. Consequently 1 − BB > 0 and the integration domain L 0 is in fact a special case of L in (63). Now suppose that instead of calculating the integral exactly we consider its high-energy asymptotics. We then get a sum of contributions of the two existing stationary points of the exponent. The first stationary point is at B = 0; its contribution coincides with our semiclassical result Z (1) .
The second stationary "point" associated with the so called Andreev-Altshuler saddle corresponds to B 11 = 0, B 22 → ∞ (see Appendix B ). For η ≫ 1 its contribution is exponentially small compared to the first one and can be dropped, but it is of the same order for small or zero η. Now importantly this contribution is related to . This is the same relation as between the two summands in our semiclassical approximation of the generating function, based on the Riemann-Siegel lookalike formula and thus unitarity. As anticipated the generating function and the two-point correlation function of fully chaotic dynamics without time-reversal invariance thus agree with those derived from the Gaussian Unitary Ensemble of RMT.
The following extension of the comparison with the RMT sigma model requires some more knowledge of that model. Previously uninitiated readers may want to first consult with introductions provided in Appendix B and in [3] or with the pioneering work of K. Efetov [4] .
Q-integral and geometrical interpretation: We first compact the RMT matrix integral (64) to Z = dQ e 1 . The matrices Q and Λ (1) obey the nonlinear constraints −Q 2 = Λ (1) 2 = 1. The integration domain for the Bosonic elements in Q can be shown to be the two-sphere S 2 and (one of the disjoint halves of) the two dimensional hyperboloid H 2 . In particular, the complex variable B 22 in the foregoing rational parametrization can be identified as a stereographic projection variable for the sphere S 2 . The point B 22 = 0 corresponds to the "south pole" while the antipodal "north pole" is approached as B 22 → ∞. The complex variable B 11 , on the other hand, is associated with H 2 .
The foregoing rational parametrization of the manifold of matrices Q is just one of many. In fact, the manifold in question can be characterized without specifying any "coordinates", and we would like to just state that characterization here, referring to Appendix B of the present paper for some details. For a geometrical interpretation of the matrix integral it is well to reorder the lines and columns in the 4 × 4 matrices Q, T, Λ as 1234 → 1324, so as to be able to employ 2 × 2 blocks as Q = QBB QF B QBF QF F ; the Bose-Bose block Q BB and the Fermi-Fermi Block Q F F contain only commuting entries while all anticommuting entries are assembled in the off-diagonal blocks Q BF , Q F B . In that "Bose-Fermi notation" (the previously used ordering is known as the advanced-retarded notation) [3] the integration domain for the blocks Q BB , Q F F in the sigma model integral become where U(2) is the 2-dimensional unitary group, and U(1, 1) is the group of 2-dimensional pseudo-unitary matrices, i.e. the group of matrices obeying U † σ z U = σ z , where σ z = diag(1, −1). Geometrically, H 2 and S 2 define the two-dimensional hyperboloid and the two-sphere mentioned before.
Andreev-Altshuler saddle point [33,34,35]: We now comment on another specific parametrization of the RMT sigma model connected to the contributions Z (2) in our semiclassical approach. This parametrization is obtained from the above rational one by the transformation Q = RQ ′ R with R a permutation matrix interchanging rows and columns as 1234 → 1432 and obeying R 2 = 1. The pertinent Jacobian is unity and the integration manifold is unchanged. The exact generating function can thus be written as the "Q ′ -integral" and the rational parametrization chosen for T in Q ′ = T Λ (1) T −1 . We have thus recovered the Weyl symmetry, , since conjugation of the matrix of energy offsets with R just amounts to ǫ C ↔ −ǫ D . On the other hand, it is most illuminating to check that the "point" B =B = 0 for Q ′ corresponds to Q = RΛ (1) R ≡ Λ (2) , the so called Andreev-Altshuler saddle. The two diagonal "saddles" are not within reach of one another through the rational parametrization with any finite B,B. Rather, Λ (2) can only be approached as T Λ (1) T −1 for B,B → ∞ and vice versa. In particular, south and north pole of the sphere S 2 are swapped in the Q ′ -integral relative to the Q-integral. In that sense one may say that the Q-integral (65) rationally parametrizes Q w.r.t. the standard saddle while the Q ′ -integral (67) rationally parametrizes Q w.r.t. the Andreev-Altshuler saddle.
We would like to emphasize that either variant of the rational parametrization yields, upon exact evaluation of the superintegral, the full generating function; the inaccessibility of the antipodal saddle does not matter since all of its neighborhood is accounted for. As already stated, our semiclassical construction of the [B,B]-integral corresponds to choosing the rational parametrization according to the Q-integral in (65), since that choice for the RMT sigma model integral does yield Z (1) in a Gaussian approximation around B =B = 0. Had we done our semiclassical analysis for Z (2) (ǫ A , ǫ B , ǫ C , ǫ D ), we would have ended up with a [B,B]-integral corresponding, in the RMT context, to the alternative rational parametrization in the Q ′ -integral in (67): The Gaussian approximation about B =B = 0 in (67) does yield Z (2) (ǫ A , ǫ B , ǫ C , ǫ D ), as is easily checked. In other words, the asymptotic contribution which stems from the infinitely remote integration point in the parametrization using the standard saddle, becomes the contribution of the stationary point B =B = 0 in the Andreev-Altshuler parametrization and vice versa.

VII. TIME-REVERSAL INVARIANT SYSTEMS
We now propose to generalize the above to time-reversal invariant systems. For them, the possibility to revert orbits (and their pieces) in time makes room for more constructive interference. We shall see that the contributions of orbits differing in encounters no longer cancel but rather generate terms of arbitrarily high orders in 1/ǫ. The summation of all relevant pseudo-orbit quadruplets can again be done with the help of a sigma-model type matrix integral.

A. Diagonal approximation and off-diagonal corrections
For time-reversal invariant systems the orbits of (A, C) can not only be repeated in (B, D) but may also be reversed in time and only then included in (B, D). Thus if we again formulate the diagonal approximation in terms of an exponentiated double sum over orbits (26), the relevant contributions to the exponent involve pairs of orbits which are either identical or mutually time-reversed. This doubling of orbit pairs leads to an additional factor 2 in the exponent, and the resulting term (ǫA+ǫD )(ǫC+ǫB ) (ǫA+ǫB )(ǫC +ǫD) in Eq. (28) needs to be squared. Moreover, we now have to consider encounters where some of the stretches are almost mutually time-reversed rather than close in phase space. For instance, Sieber-Richter pairs as depicted in Fig. 2 contain one encounter with two almost time-reversed stretches. If we change connections inside this encounter, the two links outside the encounter still look almost the same in configuration space, but one of them must be traversed with opposite sense of motion, and this is possible only if the dynamics is time reversal invariant.
Since time reversal changes the sides where a stretch enters or leaves an encounter, it becomes cumbersome to work with the notion of entrance and exit ports. Instead we arbitrarily refer to the ports on one side of each encounter as "left" ports and those on the opposite side as "right" ports. Encounter stretches may lead either from left to right or from right to left. The links may connect left and left, right and right, or left and right ports.
Our formal definition of structures is then changed as follows: S1: When numbering ports we now make sure that in A ∪ C the i-th left port is connected to the i-th right port, whereas in B ∪ D it is connected to the (i − 1)-st right port. S2: When connecting ports through links we have to allow for arbitrary connections between ports on either side of the encounters. S3: Apart from dividing the orbits before reconnection between the pseudo-orbits A and C and the orbits after reconnection among B and D we have to choose the sense of traversal of every orbit.
Again, several structures can describe the same quadruplet. In particular, we decide arbitrarily which side of each encounter is chosen as left or right. For V encounters there are 2 V equivalent possibilities of labelling the sides. To avoid overcounting, this factor needs to be divided out.
The rest of our argument proceeds in analogy to Section IV. Our semiclassical generating function turns into where the subscripts in the encounter factor depend on the first left port.
We note that some of these structures already appear for systems without time-reversal invariance, or coincide with such structures up to the directions of motion in some of the orbits. One can show that the proof of section V generalizes accordingly, i.e., all structures trivially related to those in the unitary case mutually cancel.
Let us use the semiclassical expansion (68) to calculate the leading-order oscillatory contribution to the correlation function. The oscillatory term is obtained by application of ∂ 2 /∂ǫ A ∂ǫ B to and then taking the limit ǫ A , ǫ B , ǫ C , ǫ D → ǫ. The factor [(ǫ A − ǫ C ) (ǫ B − ǫ D )] 2 tends to zero in this limit even after taking the derivatives. Therefore the oscillatory contribution to the correlator is purely off-diagonal. It must be connected to summands in Z Since inverse powers of the energy offsets can only originate from the link factors, such terms can only be due to quadruplets that contain a link belonging to the pseudo-orbit A before and to D after reconnection, and another link belonging to C before and to B after the reconnection. In particular this means that the corresponding quadruplets must contain at least two orbits, both before and after reconnection.
In the lowest order these are the quadruplets depicted in column 2:2 of Fig. 3. (We omitted quadruplets containing two copies of Fig. 5 or one copy of Fig. 5 and one Sieber-Richter pair. These quadruplets are irrelevant since the contribution of a disconnected diagram can be written as the product of its components, and the factor due to Fig.  5 vanishes by the argument in section IV E.) Among them the only diagram that explicitly requires time-reversal invariance is the one consisting of two Sieber-Richter pairs (the middle diagram in column 2:2 of Fig. 3). In order to satisfy the above conditions on membership in pseudo-orbits we only need to consider the situation where one of these Sieber-Richter pairs involves one orbit included in A and one orbit included in D, whereas the other pair involves one orbit included in B and one orbit included in C. There are 2 4 = 16 ways to choose directions of motion in each of the four orbits involved, and two choices for considering either of the two encounters as the first, leaving altogether 32 equivalent structures. Their total contribution to Z The leading oscillatory contribution to the correlator thus reads which coincides with the predictions of RMT.

B. Semiclassical construction of a sigma model
To go beyond the leading order, we again write our semiclassical results in terms of a matrix integral. The derivation of this matrix integral is very similar to the unitary case. First we encode every structure by a sequence of symbols, each symbol standing for a port, and links depicted by contraction lines. Second, we replace every symbol by a Bosonic or Fermionic variable and integrate the product with a suitably chosen Gaussian weight. The contribution of every structure will then coincide with the corresponding contractions in the resulting Gaussian integral; the sum over all possible link connections will be provided by the integral itself. Finally, summing over all encounter sets we get an integral representation of the generating function known from the zero dimensional sigma model of RMT.

Matrix elements for ports and contraction lines for links
We again write down the sequence of ports starting with the first left port of the first encounter, continuing with the first right port, etc. The left ports will be denoted by B νk,µj and the right ports byB µj,νk . The first and the second pair of indices of B respectively refer to the properties of the port before and after reconnection, and vice versa forB. As in the unitary case the Latin index j = 1, 2 indicates whether the port belongs to the original pseudo-orbit A (j = 1), or C (j = 2). The index k reveals whether after reconnection the port belongs to the partner pseudo-orbit B (k = 1) or D (k = 2).
The Greek indices µ, ν account for the directions of motion through the port and the attached encounter stretches, both in the pertinent original and partner orbits: If µ = 1 the direction is from left to right in the original orbit; if µ = 2 it is from right to left. The index ν = 1, 2 indicates the same for the partner orbits. The symbols B νk,µj ,B µj,νk can be considered as elements of 4 × 4 matrices B,B; to simplify the notation we shall often substitute a capital Latin letter for the composite index like J = (µ, j).
As for systems without time-reversal invariance the ports are ordered such that an element of B and the immediately followingB represent ports that are connected by an encounter stretch of an original orbit. This means that they have to belong to the same original pseudo-orbit, and that the direction of motion (left to right or right to left) must be the same as well. Hence their indices J = (µ, j) have to coincide. Similarly eachB and the immediately following B represent ports that are connected by a stretch of a partner orbit, and the corresponding subscripts K = (ν, k) coincide. (When aB represents the last right port in its encounter, its subscripts ν, k have to coincide with those of the first B.) Thus the arrangement of indices J and K is the same as in a product of matrices and we obtain a product of the same form as in Eq. (45) but with j and k replaced by J, K.
We now build in the links and depict them by "contraction lines" above the symbol sequences. Since for timereversal invariant systems links can be drawn between ports on either side of the encounters, the contraction lines can now connect not only B toB but also B to B andB toB. Any ports connected by a link must belong to the same pseudo-orbit, before and after reconnection and hence their Latin indices must coincide, as in the unitary case. The Greek direction indices require new reasoning. A contraction between aB and a B stands for a link connecting a right port to a left port. In this case, if an orbit leaves one encounter at the right port it must enter the other encounter stretch at the left port. The direction of motion (left to right) is thus the same for both encounter stretches. The same applies for a direction of motion from right to left, and for original and partner orbits alike. Hence for contractions betweenB and B the Greek subscripts indicating the directions of motion must coincide, µ 1 = µ 2 , ν 1 = ν 2 . On the other hand, in the contractions B . . . BB . . .B the connected ports lie on the same side, and hence their directions of motion must be opposite. These possibilities are demonstrated in Fig. 9. We shall use a bar over ν, µ to indicate that the port direction of the motion is flipped like 1 ⇀ ↽ 2 such thatμ = 3 − µ. Again denoting the energy offsets as The contribution of an encounter is defined by its first left port and can be written as i(ǫ jσ1 + ǫ ′ kσ1 ). And again, there will be a factor (−1) nC +nD depending on the numbers of orbits included in C and D.

From symbol sequence to Gaussian integral
In order to represent additive terms in (68) by contractions in a Gaussian integral we interpret the symbols B νk,µj ,B µ ′ j ′ ,ν ′ k ′ as (Bosonic or Fermionic) variables. Again any two ports which may be connected by a link must be represented by a pair of variables which are mutually complex conjugate, up to a sign. Since rule (72) allows for contractions between two B's differing in both their subscripts µ and ν we thus need to have To write the matrices B in a more elegant way, it is helpful to divide B into four 2 × 2 blocks, each block having fixed direction indices µ, ν and matrix elements indexed by k, j. The diagonally opposite blocks must be mutually complex conjugate up to a sign (which may even vary with j and k). We may write with the question mark on the equality a reminder of the as yet unspecified signs. The index p on b p , b * p signals µ = ν, i.e., the ports associated with these subblocks are traversed in parallel directions by the original and the partner orbits. In contrast, the ports associated with b a , b * a have µ = ν and their traversal directions by original and partner orbits are antiparallel.
In addition, rule (71) allows contractions between matrix elements of B andB with identical subscripts. Hence these matrix elements, too, must be chosen as mutually complex conjugate (up to a sign). Since the subscripts in B andB are ordered differently, this means thatB must coincide with the adjoint of B (again up to signs), i.e., A Gaussian-weight integral of a product of mutually conjugate matrix elements will yield a sum over all permissible ways of drawing contraction lines. To obtain the correct contribution from each link the prefactor of B µk,νj B * µk,νj in the exponent of the Gaussian again has to be chosen as i(ǫ j + ǫ ′ k ). Furthermore, due to the internal structure of B the sum over all B µk,νj B * µk,νj will include each pair of complex conjugate elements twice, which means that we have to divide the exponent by 2. Using the composite indices J, K and defining ǫ J ≡ ǫ j and ǫ ′ K ≡ ǫ ′ k independent of the Greek part of J, K we can write the exponent in the Gaussian as JK ± i 2 (ǫ J + ǫ ′ K )B KJ B * KJ . Incorporating also the factor i(ǫ Jσ1 + ǫ ′ Kσ1 ) for each encounter we obtain each (convergent) integral of this type includes the correct encounter and link factors for all structures with V encounters involving l(1), l(2) . . . l(V ) encounter stretches.

Signs
In order to incorporate the sign factor (−1) nC +nD and the correct prefactor of the diagonal approximation, we need to fully specify the matrices B andB: The off-diagonal entries in all four blocks of B are chosen Fermionic and the diagonal entries Bosonic, and the signs in (75) are picked as here the multiplication with σ z inB means that each 2 × 2 block of B † is multiplied with σ z ; taking the adjoint of our 4 × 4 supermatrices involves, apart from complex conjugation and transposition in the sense of ordinary matrices, a sign flip in the lower left (Fermi-Bose) entries of all 2 × 2 blocks. The signs in the Gaussian weight and the matrix product are again chosen such that both quantities can be written in terms of supertraces. The supertrace is now defined as the sum of the two upper left (Bose-Bose) entries minus the sum of the lower right (Fermi-Fermi) entries of the two diagonal blocks; in other words, the supertrace includes a negative sign for all diagonal elements associated to a Latin index 2 (regardless of the Greek index). The integral (76) then acquires the form withǫ,ǫ ′ 4×4 diagonal matrices whose diagonal elements areǫ J = ǫ j ,ǫ ′ K = ǫ ′ k . After eliminating all contraction lines from any integral of this type we are left with the Gaussian integral as in the diagonal factor for time-reversal invariant systems. Therefore, upon adding to (78) the Weyl factor and the factor 1/(2 V V !) from (68) we obtain the contribution to the generating function (incorporating the diagonal term) as where P J (P K ) is a projection matrix whose diagonal element associated to J (K) equals 1 while all other elements vanish.
We still have to show that with the choices made above the sign factor (−1) nC +nD for each structure is obtained correctly. The contraction rules (59) and (60) we used to prove the corresponding result in the unitary case carry over. In addition it is now possible to draw contraction lines between two B's or twoB's. These lines represent links connecting two left ports or two right ports. The corresponding contraction rules are established in Appendix A; they read analogously forB. Here Y and Z are matrix products of the type required in (80) starting and ending withB. The sequencesỸ ,Z are obtained from Y , Z by (i) reverting the order of matrices; (ii) replacing all B's byB's and vice versa; (iii) replacing all projectors P K by PK and P J by PJ . Eqs. (81) and (82) again contain the familiar factor for each link. We just chose the more compact notation The Kronecker deltas in the contraction rules (81,82) imply that the connected ports must belong to the same pseudo-orbits and have opposite direction of motion.
We can now use the contraction rules (59), (60), (81) and (82) and their counterparts for contractions of twoB's to step by step remove the contraction lines corresponding to all links of a structure. In doing so we do not meet any sign factors apart from the steps where the final link belonging to an orbit is removed. This final link must always connect a left to a right port, since otherwise it would not be possible to close it to an orbit just by drawing an encounter stretch. Hence in this step we only have to use the rules (59) and (60) for removing links between left ports (denoted by B) and right ports (denoted byB). By the same argument as in Section VI we can show that the final step yields a factor Str (P J ) for each orbit in A ∪ C and a factor Str (P K ) for each orbit in B ∪ D. According to our definition of the supertrace these factors are equal to −1 if the corresponding small indices are j = 2 (corresponding to the pseudo-orbit C) or k = 2 (corresponding to the pseudo-orbit D). We thus obtain a negative sign factor for every orbit included in C or D, and these factors combine to the desired term (−1) nC +nD .

Towards a sigma model
Summing over structures as in Section VI and again using the assumed largeness of the imaginary part η of all energy offsets, see (16), we get the generating function To see that this asymptotic series agrees with the Gaussian Orthogonal Ensemble (GOE) of RMT, we can argue like in the unitary case: The sigma model for the GOE leads to a generating function as (84) but with a specified integration domain and without the restriction η ≫ 1. The large-ǫ asymptotics of the integral is determined by two saddle points of the exponent both of which have to be taken into account unless η ≫ 1. The so called standard saddle at B = 0 leads to a contribution dominated by small B and thus coinciding with (84). The contribution of the Andreev-Altshuler saddle is related to Z (1) by Z (2)

again as in semiclassics.
Our semiclassical results thus recover the asymptotics of the generating function in agreement with the sigma model.
The same must then be true for C(ǫ). The two-point correlation function of fully chaotic time-reversal invariant systems is thus asymptotically given by the sum of two divergent series in powers of ǫ −1 in Eq. (3). Now we proceed to the full C(ǫ). As is well known, an analytic function can, under rather general conditions, be uniquely restored from its asymptotic series by Borel summation [6]. That method involves term-by-term Fourier transform of the asymptotic expansion leading to a converging series, followed by the inverse Fourier transform of the resulting analytic function. The first stage of the Borel summation applied to the two components of the correlation function gives the spectral form factor K(τ ) both for τ < 1 and τ > 1. The inverse Fourier transform recovers the closed RMT expression for the correlation function in Eq. (2), for all energies ǫ, thus completing the semiclassical explanation of universality.

VIII. CONCLUSIONS
We have shown that chaotic systems without any symmetries (beyond time reversal invariance) display universal spectral statistics. To obtain both the non-oscillatory and oscillatory contributions to the spectral correlator R(ǫ) we had to (i) represent the correlator through derivatives of a generating function Z, (ii) semiclassically expand Z as a sum over quadruplets of pseudo-orbits (sets of classical periodic orbits), and (iii) uncover a Weyl symmetry of Z. The non-oscillatory and the oscillatory part of the correlation function are due to two summands of Z related by the Weyl symmetry. The oscillatory one could only be recovered in a semiclassical approximation that explicitly preserves unitarity; in this approximation the oscillations are connected to a phase factor which accounts for the average level density.
Constructively interfering contributions to the sum over pseudo-orbit quadruplets originate from orbits that are either identical or differ in close encounters in phase space. They define asymptotic 1 ǫ -series. For dynamics without time reversal invariance those series terminate after the first term, both for the non-oscillatory and oscillatory part. For time-reversal invariant systems, the series are infinite ones but can be summed to give R(ǫ) in closed form.
To implement the summation we write the semiclassical series for the generating function Z as a matrix integral. The sum over all possible topologies then turns into an integral well known from the zero dimensional sigma model of random-matrix theory.
For both the unitary and the orthogonal symmetry class, the correlator R(ǫ) comes out as the closed-form expression familiar from random-matrix theory. In contrast to the primary 1 ǫ -series which makes sense only for |ǫ| ≫ 1, the final result holds for any real value of the energy offset and yields, for ǫ ≪ 1, the familiar power law R(ǫ) ∝ |ǫ| β , with the exponent β revealing the degree of level repulsion characterizing each symmetry class. Our results thus imply the semiclassical explanation of universal level repulsion. Likewise, all other characteristics of universal spectral fluctuations based on the two-point correlator R(ǫ), like the so called spectral rigidity [36,37], are now semiclassically understood.
The techniques here elaborated (in particular the analogy between orbit-based semiclassics and the sigma model) will be, and in fact are already helpful in tackling further challenges in the field. Examples are (i) localization effects in long thin wires [38], (ii) non-universal effects related to the Ehrenfest time (which has no counterpart in RMT) [39], (iii) a semiclassical treatment of arithmetic billiards [40] (where intrinsically quantum Hecke symmetries have strong action degeneracies as classical counterparts), (iv) spectral correlations of higher order [41] and the density of nearest-neighbor spacings, (v) transport through ballistic chaotic conductors [32,42] including the "full counting statistics" [43], (vi) further symmetry classes, (vii) transitions between symmetry classes [44], and (viii) last but not least, the relative status of the ballistic sigma model and periodic-orbit theory [23].
Using Eq. (A1) we obtain the rule (59) for contractions between B andB in different supertraces (inter-encounter contractions) which comes about as Here we used the cyclic invariance of the supertrace, then wrote out the supertrace in components, with the help of Str (P k Z) = s k Z kk for any supermatrix Z, and finally used Eq. (A1). Similarly Eq. (60) for intra-encounter contractions between B andB follows from here in the second line we could interchange U jj ′ andB j ′ k ′ , since a nonzero result arises only for j = j ′ in which case U jj is Bosonic and thus commutes with all other variables. We now turn to time-reversal invariant systems. For these systems the derivation of rules (59) and (60) carries over directly, the only difference being that J, K are now double indices consisting of the Greek and Latin parts, J = (µ, j), K = (ν, k); the sign factor s J has the values 1 for j = 1 and −1 for j = 2 regardless of the Greek index, and s K is defined analogously. In addition the fact that the matrix elements of B andB are pairwise mutually conjugate allows to draw contraction lines between two B's or twoB's. We thus obtain further contraction rules (81) and (82) for contractions between two B's, and analogous rules for contractions between twoB's. To derive these contraction rules we use the important symmetry of B as defined in (77), where B t is the transpose of B, and we have introduced Σ = 0 Now let us consider the matrix products Y and Z on the l.h.s. of (81) and (82), involving alternating sequences of B's andB's and interspersed projection matrices. Under transposition and conjugation with Σ these products have (i) the order of matrices inverted (due to the transposition), (ii) B replaced byB and vice versa (due to (A4)), and (iii) the subscripts J and K of all projection matrices replaced byJ andK (due to (A5)). Since the numbers of B's andB's are equal the sign factors from (A4) all mutually compensate. Thus equipped we can prove rule (81) for inter-encounter contractions between two B's. By using the invariance of the supertrace both under transposition of its argument and under conjugation with Σ we write Here in the third lineZ denotes the product obtained from Z by steps (i)-(iii) above. We have thus proven Eq. (81). We now move on to rule (82) for intra-encounter contractions between two B's, as in First of all note that the contracted matrix elements are B KJ and B K ′ J ′ , the latter coinciding up to a sign with B * K,J (cf. Eq. (77)). Hence the only relevant case is that J ′ =J, K ′ =K. For this case we now write and evaluate the bracket in (A8) as Here we first replaced the matrix sequence by itself double transposed; note that double transposition of a supermatrix (not an identity operation!) leaves its Bosonic elements like the one with the subscripts JJ unchanged. Then we invoked Eqs. (A4,A5) and noted that for the matrix elements at hand transposition just amounts to interchanging the two subscripts. In the third line we used that only those elements of Σ, Σ t are nonzero for which one subscript is equal to another one barred, i.e., the port directions of motion are opposite. Finally, we used Σ tJ ,J ΣJ ,J = s J Inserting (A9) into (A8) we now obtain which is (82). Here we used (A1), definedỸ in analogy toZ, and in the end rearranged the result as a supertrace. One easily shows that rules analogous to (81) and (82) also hold for contractions between twoB's.

APPENDIX B: SUPERSYMMETRIC SIGMA MODEL
We propose to briefly review the basics of the sigma model; details can be found in [3,4,14,31,33,34]. The sigma model is a way to implement averages over random matrices. The appropriate random-matrix ensemble for systems without time-reversal invariance is the Gaussian Unitary Ensemble (GUE), consisting of complex Hermitian N × N matrices H with a Gaussian weight proportional to e − N 2 tr H 2 ; the matrix dimension N is sent to infinity. The pertinent generating function Z looks like (5), but with the energy average . . . replaced by the GUE average. To perform that average, three variants of the sigma model are available, all based on representing determinants by suitable Gaussian integrals. The supersymmetric variant deals with superintegrals over both ordinary and Grassmannian variables while two replica variants exclusively use either Grassmann or ordinary Gaussian integrals.

From GUE to supermatrix integral over saddle-point manifold
In the supersymmetric approach each of the determinants in the generating function (5) is represented by a Gaussian integral, Fermionic for the numerator and Bosonic for the denominator, at the expense of the overall factor (−1) N the signs of the "Bosonic" quadratic forms z * A (. . .)z A and z * B (. . .)z B have been chosen so as to make the Bosonic Gaussian integrals convergent, in consistency with Im ǫ A,B,C,D > 0. η C , η * C , η D , η * D are anticommuting Grassmann variables. As usual, the asterisk on commuting variables means complex conjugation.
On anticommuting variables the asterisk may, for the purposes of the present paper, be given the same interpretation (complemented by (η * ) * = −η, see [3]). At this point hooking onto the presentation in Chapt. 10 of Ref.
[3] 7 , we introduce the supervector where each of the four Φ α 's has N components Φ αm ; the ordering of the Φ α in (B2) is known as the advancedretarded notation. To accommodate the relative sign of the Bosonic quadratic forms mentioned and the different energy arguments above we employ the (4 × 4) matrices, whose block form will be convenient in the future. Our generating function may then be succinctly written as with dΦ * dΦ denoting the integration measure written explicitly in (B1). The Gaussian average over Hamiltonian matrices H with a normalized Gaussian weight ∝ e − N 2 trH 2 leads to Z = (−1) N exp − 1 2N allows to further compact Z as where the symbol Str denotes the supertrace. The integrand would be invariant under the transformationQ → TQT † with pseudounitary matrices T defined by if the (real part of) the energy matrix were proportional to the unit matrix. Due tō the energy offsets ǫ A,B,C,D /2πρ are of only order 1/N , which means that in any case this symmetry persists to leading order in 1/N . The exponent in (B6) is quartic in the integration variables Φ, Φ * . Upon employing a suitable Hubbard-Stratonovich transformation [3] we render the Φ-integral Gaussian and end up with an integral over 4 × 4 supermatrices Q, with a flat integration measure dQ. For a details of the transformation, in particular the overall minus sign and the integration range in (B9), we refer the reader to Ref. [3]. Important for what follows is that the integration range is restricted to matrices Q diagonalizable by pseudounitary transformations T . The large-N limit allows for a saddle-point approximation. The saddle manifold is determined by We may confine the discussion of this equation to E = 0, which means that we evaluate Z in the middle of Wigner's semicircle. As seen the energy offsetsÊ − E are of order 1/N and thus irrelevant for saddles, we may even dropÊ entirely. The saddle-point equation then simplifies to Q 2 = −1.
The eigenvalues of the matrices Q must be equal to ±i. This means that the saddle-point manifold can be divided into 5 conjugacy classes differing by the numbers n + , n − of the eigenvalues +i and −i. In fact, only one of them is of relevance [3], with n + = n − = 2. Since the matrices Q must be diagonalizable by pseudo-unitary transformations the relevant Q's can thus be represented as with the diagonal matrix

Rational parametrization
There is an ambiguity in the choice of the matrices T = Taa   Tra   Tar Trr . Due to the degeneracy of the eigenvalues of Λ (1) , Q remains unchanged if we replace T by the product T K with K any invertible block diagonal matrix. Choosing with B = T ar T −1 rr andB = T ra T −1 aa . The pseudounitarity of T entails the matrix (T K) † L(T K) = K † LK = to be block diagonal. On the other hand such that block diagonality yieldsB = σ z B † . With the latter relation between B andB established we can write these 2 × 2 supermatrices as The final integral over the manifold becomes an integral over the 2 × 2 matrices B,B, π dB 12 * dB 12 dB * 21 dB 21 . The integration ranges for the Bosonic integrals over B 11 and B 22 are respectively |B 11 | 2 < 1 and |B 22 | 2 < ∞. Note that the foregoing integral representation of Z (1) equals the one established semiclassically, see (63).
We now proceed to the asymptotics of the remaining integral (B18) in the limit of large energy offsets, |ǫ A,B,C,D | ≫ 1. To that end we employ one more saddle-point approximation 8 . The only finite stationary point of the exponent in (B18) occurs at B =B = 0. Let us assume, momentarily, that the energy offsets are provided with a large positive imaginary part, η ≫ 1. Then this will be the only critical point of the integrand responsible for the high energy asymptotics of the integral (B18). To get an asymptotic series of Z in powers of ǫ −1 we expand the matrix (1 − BB) −1 and (1 −BB) −1 in the exponent in powers of B,B about the unit matrix and further expand the resulting exponential about its Gaussian part. Scaling the integration variables as B = B ′ / |ǫ| and similarly forB with a typical dimensionless energy offset ǫ close to ǫ A,B,C,D we see that the higher-order terms in the expansion indeed become small perturbations in the high-energy limit. The Gaussian superintegral giving the leading-order contribution reads The integration domain for the Bosonic variables can be extended to the full complex plane; in view of η ≫ 1 the error thus introduced is exponentially small. Upon doing the remaining elementary integrals we find

Evaluation of the matrix integral
In fact the result (B21) fully exhausts Z (1) : In agreement with the Duistermaat-Heckmann theorem [45], all further terms of the perturbation series vanish for the unitary symmetry class; in the framework of the semiclassical theory this was already demonstrated in Section V.
To prove the absence of perturbative corrections in the present RMT context we transform integration variables such that the integral in (B18) becomes Gaussian and similar to (B20). To that end we introduce an analogue of the singular-value decomposition of the matrix B, with unitary U and pseudo-unitary V , and parametrize the latter matrices as The matrix B is diagonal with positive (numerical parts of the) diagonal entries, with the Bose-Bose entry obeying with the Fermi-Fermi entry negative. It further follows that the matrix products BB andBB are respectively diagonalized by U and V , with the eigenvalues 0 ≤ l B < 1, l F ≤ 0. The matrices U, V also diagonalize the rational functions of BB andBB appearing in the sigma-model integral (63), The representation (B22) can be inverted to express the variables on the right-hand side in terms of the original elements of B,B. The uniqueness of the decomposition of our B is thus ascertained. Moreover, the eight new variables l B , l F , φ B , φ F , η * , η, τ * , τ can be employed as integration variables instead of the original elements of B,B. A somewhat tedious calculation of the Berezinian yields the measure Our next step is yet another change of integration variables which will be elements of matrices analogous to B,B but with l B,F replaced by m B,F , We obviously have such that the integrand in the sigma-model integral becomes Gaussian in terms of the matrices C andC. Moreover, in complete analogy with (B26) the integration measure can be written as Comparing the measures (B26) and (B29) we get Therefore, in terms of the integration variables C ik , C * ik the sigma-model integral becomes Gaussian. Again assuming that the energy offsets are all provided with the large positive imaginary part η we can extend the integration domain for the Bosonic variables to the full complex plane, at a risk only of corrections of the order exp(−η). Then the resulting C-integral will become identical with the leading-order integral (B20). We have thus proven that contribution of all higher-order terms in the expansion of the integrand about the saddle point B = 0 cancels, and the asymptotics of the generating function, in all orders of 1 ǫ , is exhausted (for η ≫ 1) by (B21).

Andreev-Altshuler saddle point
If the imaginary part η of the energy arguments is small or zero Z (1) does not capture the full high-energy asymptotics of the generating function: an additive component Z (2) appears due to another, the Andreev-Altshuler or "anomalous" stationary point B AA . To show it let us use the parametrization (B22)) introduced in the previous subsection; note that the elements of the diagonal matrix B satisfy 0 ≤ x B < 1, 0 ≤ x F < ∞ . Using (B24),(B28) we can rewrite the supertrace in the exponent of the sigma-model integral like At a stationary point the derivatives of this expression by x B and x F must be zero, from which it follows that x B must be zero, whereas x F can be zero or infinity. Consequently there exist two stationary points in the matrix space, (a) the standard saddle B = 0 giving birth to Z (1) , and (b) the so called Andreev-Altshuler saddle B AA = diag (0, ∞ C ) where ∞ C stands for the infinitely remote point of the complex plane. The anomalous saddle is responsible for the part Z (2) of the generating function exponentially small when η ≫ 1 but of the same order as Z (1) if the energy offsets are real. An elegant way to find the contribution of the anomalous stationary point involves retracing the steps leading to the rational parametrization and finding the point Q AA of the manifold (B12) associated with B AA . Let us represent B AA andB AA as limits B AA = lim a→0 P K −1 a ,B AA = lim a→0 P K −1 −a with the diagonal 2 × 2 matrices for simplicity we take a real, i.e., we approach the infinite point of the complex plane along the real axis. When a is finite the respective 4 × 4 matrix T = 1 factorizes, (B33) Since the second factor in T a commutes with Λ (1) we find The matrix R describes a permutation transposing the second and fourth entries, i. e., the Fermi-Fermi ones. We thus have To get the part of the high-energy asymptotics connected with the Andreev-Altshuler stationary point we reparameterize the integration variables in (B11) as Q ′ = RQR. The Jacobian of the transformation Q → Q ′ is unity. Denoting byÊ ′ the offset matrix with the Fermi entries transposed,Ê ′ ≡ RÊR = diag(ǫ B , −ǫ C , −ǫ A , ǫ D ) we bring the Q-integral to an equivalent form Z = − Q ′2 =−1 dQ ′ e i 2 Str Q ′Ê′ . Now we can set Q ′ = T Λ (1) T −1 and rationally parameterize T which leads to a superintegral differing from (B18) by the interchange ǫ D ↔ −ǫ C . Expanding the exponent around the stationary point B =B = 0 corresponding to Q ′ = iΛ (1) , or Q = Q AA = iΛ (2) we obtain an asymptotic series in ǫ −1 which truncates on its first term, Z (2) = Z (1) | ǫD ↔−ǫC . Adding the contributions of both stationary points we get the full GUE generating function Z = Z (1) + Z (2) .

APPENDIX C: BOSONIC REPLICAS
The main technical difficulty in establishing a one-to-one correspondence between semiclassics and the sigmal model of RMT are the spectral determinants in the numerator of the generating functions. These determinants lead to unwelcome sign factors (−1) nC +nD in the semiclassical treatment, and in RMT they cannot be expressed through Gaussian integrals over complex numbers. In the supersymmetric treatment the way around these difficulties is to introduce anticommuting Fermionic variables. These variables incorporate the appropriate signs in semiclassics, and allow for an integral representation of non-inverted determinants in RMT. In the present appendix we want to present an alternative solution based on the so-called replica trick. (In Appendix D we will consider a combination of the replica trick with Fermionic variables.) While the replica trick offers considerable technical convenience, there is a price to pay: First, the replica trick is hard to justify mathematically; second, at least in the present (Bosonic) variant only the part Z (1) is accessible. That price has let replicas fall into some disrepute. From a semiclassical point of view, the price looks more affordable. We have seen Z (1) to yield the full two-point correlator with the help of the "crosswise" procedure of Sect. II A or, equivalently, Z (2) to be accessible from Z (1) by the Riemann-Siegel lookalike of Sect. II C.

Semiclassical construction of a sigma model
To construct a replica sigma model from semiclassics, we start from formula (37). We the assume that the pseudoorbits C and D are each divided into r − 1 subsets yielding altogether 2r pseudo-orbits. We denote the pseudo-orbit A by j = 1, and the r − 1 parts of C by j = 2, . . . , r. Similarly B is indicated by k = 1, whereas the r − 1 parts of D are labelled by k = 2, . . . , r. Each of the n C orbits in C can be included in any of the r − 1 subsets of C, and similarly for D. This leaves altogether (r − 1) nC +nD ways for distributing the orbits of C and D among subsets. If we now formally take the limit r → 0, this number of possibilities turns into the factor (−1) nC +nD needed in our semiclassical treatment. We can thus drop the sign factor if we write Z . (C1) where we sum over structures of 2r-tuples rather than quadruplets of pseudo-orbits. This means that now every way of distributing the orbits in C and D among the subsets is considered to give rise to a separate structure. The indices of ǫ j , ǫ ′ k in each link factor are chosen according to the pseudo-orbits j, k that the link belongs to; the corresponding indices in the encounter factor are chosen according to the pseudo-orbits the first entrance port belongs to. In line with the above definitions of j, k we take For the diagonal approximation we write where the factors −i are just inserted for later convenience; they cancel mutually in the replica limit r → 0. Now can now construct a sigma model similarly as in the main part: Each entrance port belonging to an original pseudo-orbit j (either A or a part of C) and a partner pseudo-orbit k (either B or a part of D) is indicated by a symbol B kj . Correspondingly each exit port is denoted by a symbolB jk . The indices of B's andB's belonging to the same encounter have to coincide like the subscripts in the trace of a matrix product. We can again sum over all ways of connecting these ports by links with the help of a Gaussian integral over B kj ,B jk . This is done in the same way as in the supersymmetric approach, with only one difference: Since the sign factor (−1) nC +nD has been eliminated by the replica trick, all sign factors in the previous treatment can be omitted including those in the commutation relation of B,B. The coefficients B jk ,B kj can thus be taken as ordinary complex variables withB jk = B * kj (i.e.B = B † ). All supermatrices are thus replaced by ordinary matrices which are now larger due to the r different values of j and k.
To go through the argument in detail, we split the sum over structures into sums over the possible numbers and sizes of encounters, assignments of ports to pseudo-orbits, and ways of connecting the ports through links: The three subsums are done step by step. Link sum: We start by summing over link connections for fixed encounters and assignment of ports to pseudoorbits. We denote by L in kj , L out kj the numbers of entrance and exit ports forming part of the original pseudo-orbits j = 1, 2, . . . r and the partner pseudo-orbits k = 1, 2, . . . r. The links now have to connect each of the L out kj exit ports of j, k to one of the L in kj corresponding entrance ports, and we can conclude L in kj = L out kj for all j and k. With that condition satisfied, all connections between entrance and exit ports associated with the same j and k are permissible, i.e., there are jk L in kj ! = jk L out kj ! possibilities. For each j, k we obtain a combinatorial factor δ L in kj ,L out kj L in kj !. We now combine the latter combinatorial factor with the L kj +1 contributions 1 −i(ǫj+ǫ ′ k ) , one from each link within j and k and one from the diagonal approximation. We can represent the product as an integral over complex parameters B kj , where we letB jk = B * kj , , for each pair j and k. The integral converges due to Im ǫ j , Im ǫ ′ k > 0. The only terms not yet taken into account are the encounter contributions i(ǫ j + ǫ ′ k ), and the factor e i( V ! from Eqs. (C1) and (C3). With these terms included the link sum yields the generating function We have thus re-expressed Z (1) r as an integral over arbitrary complex r × r matrices B. Port sum: Next, we sum over all possibilities to assign entrance and exit ports of encounters to pseudo-orbits j and k. These assignments will determine L in kj and L out kj and thus the numbers of factors B kj ,B jk for the summands in Eq. (C6). To perform the sum in a compact manner, we arrange the factors B kj ,B jk such that all factors originating from ports of the same encounter are grouped together. Each of these groups starts with the factor associated to the first entrance port, followed by the first exit port, the second entrance port, etc., up to the last exit port. This leads to an alternating sequence of B kj 's andB * jk 's. As we had seen, the subscripts in this sequence must coincide as for the trace of a matrix product. We can also absorb the encounter factor from (37) into the latter sequence. The indices in that factor represent the original and partner pseudo-orbits to which the first entrance port belongs and thus must coincide with the subscripts of the first B in the sequence above. If we sum over all assignments of ports to pseudo-orbits by summing over all independent indices, each l-encounter then gives rise to a factor j1,...,j l ,k1,...,k l i(ǫ j1 + ǫ ′ k1 )B k1,j1Bj1,k2 B k2,j2 . . . B k l ,j lB j l ,k1 = i trǫ(BB) 2 + i trǫ ′ (BB) 2 (C7) withǫ = diag(ǫ 1 , ǫ 2 , . . . , ǫ r ) = diag(ǫ A , ǫ C , . . . , ǫ C ) If we proceed in this way, and also write the exponentials of Eq.
Encounter sum: The remaining summation over the number V of encounters and over their sizes l(σ), σ = 1, 2, . . . , V , is trivial. We get 2. The sigma model of RMT Eq. (C10) agrees with the random-matrix result, which we shall derive in the following. In random-matrix theory the replica trick amounts to replacing the spectral determinants in the numerator of the generating function by powers ∆(E) −(r−1) with the number of replicas r an integer larger than unity. One evaluates the thus generalized generating function for r ≥ 1, but in the very end takes the limit r → 0. This leads to the following representation of the generating function The advantage of the replica trick is that now only inverted spectral determinants show up, which can be represented by Gaussian integrals over ordinary complex variables. If we go through the same steps as for the supersymmetric sigma model (average over matrices representing the Hamiltonian, Hubbart-Stratonivich transformation, stationaryphase approximation) we obtain the following integral representation for Z(r), where we dropped a proportionality factor converging to 1 in the replica limit r → 0. In contrast to the supersymmetric case, the matrices Q are now complex 2r × 2r matrices. They have to be diagonalizable bt pseudo-unitary transformations T which are now defined by where Λ (1) = (1, 1, . . . , −1, −1, . . .). The saddle point manifold Q 2 = −1 has the shape of a hyperboloid. In the replica case the relevant solutions of Q 2 = −1 are matrices that can be obtained through conjugation from iΛ (1) where , Multiplication with a block diagonal matrix as in the supersymmetric now leads to the rational parametrization where B is an arbitrary r × r matrix andB = B † . The fact that the off-diagonal blocks of this matrix are Hermitianconjugate to one another follows from the pseudo-unitarity condition.
In the limit of large ǫ the integral (C12) can be calculated in a saddle point approximation (not to be mixed with the already implemented limit N → ∞). Using the parametrization (C15) in Q and expanding the exponent in the vicinity of the saddle point B = 0 we obtain the integral agreeing with our semiclassical formula. Both the semiclassical and the RMT results can be generalized to the orthogonal case. One then obtains a formula as in (C16) but with B as a 2r × 2r matrix subject to symmetries analogous to the supersymmetric treatment. In fact, even the crossover between the unitary and the orthogonal symmetry class is accessible in this way [44].
A further variant of the sigma model makes use both replicas and Fermionic variables (the latter appearing only in the intermediate stages of the derivation from RMT). To derive this variant from semiclassics, we write the sign factor in Eq. (37) as Here the factor (−1) nA+nB can be dealt with in a similar way as the factor (−1) nC +nD in the Bosonic approach. This time we split the pseudo-orbits A and B into r − 1 parts. There are (r − 1) nA+nB ways to do this. If we take the limit r → 0 this number of possibilities converges to (−1) nA+nB . The remaining sign factor (−1) (nB +nD )−(nA+nC ) depends on the difference between the overall numbers of partner orbits and original orbits. To keep track of this factor we envisage the reconnections inside encounters as a step-bystep process. In each step the connections of just two stretches are changed, as in Fig. 5. This changes the number of orbits by one and thus flips the above sign. If we only have one 2-encounter we need only one such step and the sign factor is negative. For the example of a 4-encounter, Fig. 10 shows how we can get from connections in (a) to the final connections in (d) in three steps, each time only reshuffling two stretches. In general for each l-encounter we need l − 1 steps. Let us denote by v l the number of l-encounters. Then the overall number of steps, and thus sign flips, is given by the number L = l v l l of encounter stretches or links minus the number V = l v l of encounters. We are thus free to replace the sign factor by (−1) L−V .
Altogether we then obtain .

(D2)
Here the definition of structures is adapted to consider separately situations where the orbits in A or B are distributed differently among the r − 1 subsets of A and B. Furthermore ǫ j , ǫ ′ k (and the corresponing diagonal matrices) are now defined asǫ = diag(ǫ 1 , ǫ 2 , . . . , ǫ r ) = diag(ǫ A , . . . , ǫ A , ǫ C ) For the diagonal term we have If we use (D2) as a starting point we can repeat the semiclassical construction of the sigma model from the previous appendix. We just have to use the new definitions ofǫ,ǫ ′ and insert sign factors in the oscillatory prefactor, for each of the V encounters, and for each of the L encounter stretch represented byBB, BB. This yields

The sigma model of RMT
In RMT, the Fermionic replica sigma model is obtained if we represent all spectral determinants the generating function (5) by Gaussian Grassmann integrals. To that end one brings the determinants of the denominator "upstairs" as [34] det −1 = det r−1 , with the replica parameter r first taken as an integer larger than unity; only in the end, r is assumed real and sent to zero. The spectral determinants can the be represented through Gaussian integrals over Fermionic variables. By going through the same steps as in the supersymmetric and the Bosonic replica variant, these integrals can be traded for and integral over complex 2r × 2r matrices which are diagonalizable by unitary transformations and satisfy Q 2 = −1, withǫ,ǫ ′ as in (D3). The matrices Q relevant for the non-oscillatory terms in the correlator are again of the type The rational parametrization then amounts to which inserted into the integral (D6) yields just as in semiclassics. Again both the semiclassical and the RMT calculation can be generalized to the orthogonal case. The contribution of the saddle B =B = 0 does not exhaust the high-ǫ asymptotics of the generating function. Another contribution is connected with the vicinity of the so called anomalous or Andreev-Altshuler saddle, Λ (2) = diag(−1, 1, . . . , 1; 1, −1, . . . , −1) which can be obtained from Λ (1) by swapping the first and the (r + 1)-th entry, i. e., where all missing entries are zero. Λ (2) cannot be reached using the rational parametrization (D8); its vicinity corresponds to B,B tending to infinity.
To appreciate the contribution of the Andreev-Altshuler saddle, we substitute Q = RQ ′ R † in the integral (D6) and adopt Q ′ as the new integration variable. The Jacobian of this transformation is unity such that we get Z r = Q ′2 =1 dQ ′ e 1 2 tr Q ′ R † XR an expression still equivalent to (D6), the replacement of X by R † XR = diag (−ǫ ′ 1 , ǫ 2 , . . . , ǫ r ; ǫ 1 , −ǫ ′ 2 , . . . , −ǫ ′ r ) notwithstanding. We then employ the rational parametrization for Q ′ and obtain the asymptotic series in ǫ −1 expanding the integrand around the stationary point B =B = 0 corresponding to Q ′ = Λ (1) or Q = Λ (2) . Again, the Duistermaat-Heckmann theorem guarantees that the series breaks up after the leading term. The final result for Z (2) differs from Z (1) only by the appearance of the permuted offset matrix R † XR instead of X, i. e. Z (2) (ǫ A , ǫ B , ǫ C , ǫ D ) = Z (1) (ǫ A , ǫ B , −ǫ D , −ǫ C ) .

APPENDIX E: STRUCTURES, PERMUTATIONS AND SYMMETRY
Here we formalize the definition of a structure in terms of permutations, in order to clarify the relation of structures to diagrams and investigate the symmetry aspects of that relation. A "diagram" is understood as the set of all physically distinct pseudo-orbit quadruplets with a given topology, i.e, with a fixed set of encounters whose ports are connected by links in a specified way. Such a diagram can have many equivalent structures associated. For simplicity we limit the following to the unitary symmetry class.
We imagine that the encounters are numbered by σ = 1, . . . , V and the stretches and ports in the σ-th encounter by j = 1, . . . , l(σ). Each port is thus denoted by a label σ.j. According to the convention of section IV B in the original pseudo-orbits A, C each entrance port is connected with the exit port of the same number. After the reconnection leading to the pseudo-orbits B, D the entrance port σ.j will be connected to the exit port σ.(j − 1), with the exception of the first entrance port of each encounter which will be connected to the last exit port. The port reconnections can be depicted by a table in which the lower and the upper row list the entrance and exit ports, respectively, and each encounter is represented by a block of width l(σ), The connection of exit to entrance ports by links can also be represented by a similar permutation P link where the upper row refers to exit ports and the lower one to entrance ports. Unlike P enc the link permutation is totally arbitrary.
The permutation products P link I = P link and P link P enc define connections between the consecutive entrance ports before and after reconnection, respectively. Each cycle of P link represents the sequence of the entrance ports of a pre-reconnection orbit in A, C; similarly, each cycle of P link P enc defines such a sequence for a post-reconnection orbit in B or D. The numbers n, n ′ of cycles in these permutations give the numbers of pre-and post-reconnection orbits.
Our definition of a structure can thus be formalized as a particular choice of (S1) the numbers V and l(σ) determining the permutation P enc , (S2) a link permutation P link , and (S3) a division of the n cycles of P link between the pseudoorbits A and C, and of the n ′ cycles of P link P enc between the pseudo-orbits B and D.
The order in which the encounters are numbered and the choice of the first stretch in an encounter is irrelevant. Their change amounts to port renumbering, i.e., also to a permutation (an arbitrary permutation of the encounter numbers, and a cyclic permutation of the labels for each encounter). The possible renumberings corresponding to a given V , l(σ) form a group C with N (C) = V ! V σ=1 l(σ) elements. The renumbering C ∈ C transforms the encounter and link permutations as P → CP C −1 . For a structure devoid of any symmetry, none of the non-trivial C ∈ C commute simultaneously with the link and encounter permutation, hence renumberings will produce V ! V σ=1 l σ structures described by P ′ enc = CP enc C −1 , P ′ link = CP link C −1 , all of them corresponding to the same diagram and making the same contribution to the generating function. This overcounting is cancelled by the denominator 1/N (C) in (41) such that we are left with the contribution of the respective diagram in the form of the ratio of the encounter and link factors.
In a symmetric topology C may simultaneously commute with P link and P enc , and hence the number of equivalent structures is smaller than N (C). However renumbering the encounters and stretches will still change the variables s σj , u σj parametrizing the separations inside the encounters; the indexing of the link durations is also changed. Since we integrate over all values of s σj , u σj as well as the link durations, relabelling thus leads to a separate contribution even in the symmetric case. To count each quadruplet only once we thus always have to divide by the number of possible relabellings N (C).
Alternatively, we could formulate our approach in terms of a to sum over diagrams rather than structures. However, then the diagrams with symmetries would require special treatment. This was done in [46] for the diagrams with L − V = 2, but it becomes too cumbersome for a generalization to arbitrary orders.
For an example we consider a quadruplet containing, before and after the reconnection, one periodic orbit with two 2-encounters connected as shown in Fig. 11. The encounter and link permutations with the port numbering shown in Fig. 11  where we put primes over the indices of the exit ports, to distinguish them from the entrance ports. We have the following possibilities to relabel encounters and stretches (i) interchange all first indices 1 and 2, (ii) interchange 1. 1  FIG. 11: A symmetric orbit pair with two encounters, already shown in the gallery of Fig. 3 as the uppermost one in the left column for L − V = 2; here, one of the possible port numberings is indicated. and 1.2, and/or (iii) interchange 2.1 and 2.2. As (i) and simultaneous application of (ii) and (iii) both leave P enc and P link invariant there are only two equivalent structures associated with the diagram. However application of (i) and simultaneous application of (ii) and (iii) still changes the parameters s σj and u σj . Hence when calculating the contribution of each structure by integrating over all parameters each quadruplet will be counted four times per structure. To avoid avoid overcounting we thus have to divide out N (C) = 8, as we would without any symmetries.

APPENDIX F: COMPLEX VS REAL CORRELATOR
Let us prove the relation R(ǫ) = lim Im ǫ→+0 Re C(ǫ) with the complex correlator defined by (4), in particular the presence of the constant 1/2 in (4). Averaging in the definition of the spectral correlator, be it real or complex, is carried out over an interval of its both arguments, namely the central energy E and the energy offset ǫ/πρ; to obtain a smooth function the product of the two averaging intervals must be large compared with the squared mean level spacing, δE δǫ/πρ ≫ (1/ρ) 2 .
If we want oscillations of the correlator on the scale of the mean level spacing to be resolved the averaging interval of ǫ has to be taken small, δǫ ≪ 1. In view of (F1) this means that the averaging interval [a, b] of the central energy with δE = b − a must be large compared with the mean level spacing multiplied by 1/δǫ. On the other hand the averaging interval must be small compared with the range populated by the consecutive levels E 1 < . . . < E N . (The generalization to unbounded spectra or bounded spectra with N = ∞ is straightforward.) We define the average over the central energy as . . . E = (b − a) −1 b a dE (. . .). In the energy range E 1 ≤ E ≤ E N the mean level spacing 1/ρ can be taken as constant. Without loss of generality we choose to discuss spectral correlations in the middle of the range considered, i. e. at E = EN +E1 2 and can therefore take that energy as the center of the averaging interval [a, b]. By explicitly doing the average over E in the definition of the correlator (4) we obtain where the remaining average is over ǫ only. Using (x + iη) −1 ↔ P x −1 − i πδ (x) for η → 0 + and ln x = ln |x| + i Arg x, we can write the real part of C as Re C (ǫ) = − 1 2 + R I,I + R R,R where R I,I (R R,R ) stem from products of the imaginary (real) parts of the fraction and the square bracket.
We begin with R I,I assuming Im ǫ + = η → +0. In the summand (i, k) the phase of the expression in the square brackets can be (i) 2π (both E i and E k lie within the averaging interval of the central energy [a, b]); (ii) π (only one eigenvalue lies within [a, b]); (iii) zero (both eigenvalues are outside [a, b]). Let us assume that the energy offset ǫ/πρ is small compared with b − a; then considering the factor δ ǫ πρ + E i − E k , we may neglect the possibility (ii). Up to an additive constant, R I,I must then coincide with the real correlation function, In the sum of products of real parts are cancelled by the zero of the respective logarithms; the mild logarithmic singularities of R R,R at E i = a, b ± ǫ 2πρ are unconnected with the level correlation and thus of no physical relevance. (They are caused by the rectangular cut-off of the integrand at the borders of the averaging interval in . . . E and disappear if this cut-off is smoothed out). Hence the summand has are no relevant singularities or other fluctuations on the scale of ǫ. We are thus free to replace the double sum over a discrete eigenvalue spectrum by a double integral over a uniformly distributed density, as in Note that the smoothing w.r.t. the offset variable is thus rendered superfluous. Consider now the semiclassical limith → ∞ which means thatρ → ∞; we also assume the spectral span constant (in classical units) which means that the number N of levels taken into account goes to infinity. Let the averaging interval [a, b] simultaneously shrink compared with the spectral span in such a way that the hierarchy E N −E 1 ≫ b−a ≫ ǫ/πρ is observed. Neglecting ǫ πρ and letting the parameter A ≡ EN −E1 b−a go to infinity we get Collecting (F2), (F3), and (F4) we come to the desired relation, lim η→0 Re C (ǫ + ) = R (ǫ).