Efficient semiclassical approach for time delays

The Wigner time delay, defined by the energy derivative of the total scattering phase shift, is an important spectral measure of an open quantum system characterising the duration of the scattering event. It is related to the trace of the Wigner-Smith matrix Q that also encodes other time-delay characteristics. For chaotic cavities, these exhibit universal fluctuations that are commonly described within random matrix theory. Here, we develop a new semiclassical approach to the time-delay matrix which is formulated in terms of the classical trajectories that connect the exterior and interior regions of the system. This approach is superior to previous treatments because it avoids the energy derivative. We demonstrate the method's efficiency by going beyond previous work in studying the time-delay statistics for chaotic cavities with perfectly connected leads. In particular, the universality for moment generating functions of the proper time-delays (eigenvalues of Q) is established up to third order in the inverse number of scattering channels for systems with and without time-reversal symmetry. Semiclassical results are then obtained for a further two orders. We also show the equivalence of random matrix and semiclassical results for the second moments and for the variance of the Wigner time delay at any channel number.


Introduction
The concept of time delay plays an important and special role in quantum collision theory. It was first introduced by Eisenbud and Wigner [1,2] in the context of elastic scattering, who related the time during which a monochromatic wave packet is delayed in the interaction region with the energy derivative of the scattering phase shift. Later on Smith [3] extended the concept to the case of inelastic scattering with many channels and introduced the time-delay (or Wigner-Smith) matrix where S(E) is the scattering matrix at energy E whose dimension is given by the number of open channels, M . For an energy and flux conserving scatterer (i.e. involving no absorption or dissipation), the S-matrix is unitary and Q is a Hermitian matrix, which is manifest from the second representation in (1). The real characteristics of Q (e.g. its diagonal elements and eigenvalues) provide us with various time-delay related observables [3]. Taking the trace, one arrives at the simple weighted-mean measure of the collision duration [4,5,6,7] This expression defines the Wigner time delay in multi-channel scattering. Further and more subtle discussions of the physical concepts of time delays and their applications to regular and chaotic scattering can be found in recent reviews [7,8,9]. As is clear from the definition, the dominant contributions to τ W come from the regions where the total scattering phase shift exhibits a sharp energy dependence. Away from the channel thresholds, this occurs in the vicinity of resonances that correspond to the meta-stable (decaying) states formed at the intermediate stage of the scattering event. Such resonances can be analytically viewed as the complex poles E n = E n − i 2 Γ n of the S-matrix in the complex energy plane, with E n and Γ n being the position and width of the n-th resonance, respectively. In the resonance approximation, neglecting the constant background connected to potential scattering and direct reactions, these poles are the only singularities of S. In view of the unitarity condition, their complex conjugates E * n serve as the S-matrix zeros, thus implying det S = n E−E * n E−En . This results in the important connection [4,6] between the Wigner time delay and the resonance spectrum of the intermediate open system. This expression is valid at arbitrary degree of the resonance overlap and thus the Winger time delay (2) can also be treated as the density of states of the open system, see [10,11] for relevant discussion. In particular, the spectral average of τ W over a narrow energy window isτ W = 2π ρ/M , whereρ is the mean density around E. The latter in turn determines the fundamental timescale in quantum systems, the Heisenberg time T H = 2π ρ, so thatτ W = T H /M . In many experimental situations, the intermediate system is characterised by very complicated internal motion, with representative examples being microwave cavities [12], quantum dots [13] and compound nuclei [14]. As a result, the time delay (as well as other scattering observables) reveals strong fluctuations around its mean value that need to be treated statistically. There are two main approaches to describe such fluctuations: random matrix theory (RMT) [15] and the semiclassical method [16].
The general RMT approach to the time delay, mostly developed in [6,7,17], is based on a representation of the Wigner-Smith matrix in terms of the effective non-Hermitian Hamiltonian H eff of the open system. The complex resonances E n are then given by the eigenvalues of the Hamiltonian [18]. The Hermitian part of H eff describes the internal Hamiltonian with a discrete spectrum and chaotic dynamics. It can therefore be modelled by a random N ×N matrix drawn from the Gaussian orthogonal (GOE) or unitary (GUE) ensemble, depending on the presence or absence of timereversal symmetry (TRS), respectively [15]. The anti-Hermitian part originates from coupling between the N internal and M channels states and has a specific algebraic structure due to the unitarity constraint [18]. The statistical averaging can then be performed by making use of the powerful supersymmetry technique, mapping the problem to a nonlinear supersymmetric σ-model of zero-dimensional field theory [19]. In the asymptotic limit N → ∞, the obtained results are universal in the sense that local fluctuations on the scale of mean level spacing 1/ρ ∼ N −1 become modelindependent provided that the appropriate natural units have been used [20]. The only parameters are then the number M of open channels and their strength of coupling to the continuum. The later is quantified by the average 'optical' S-matrix, with S = 0 (or S → 1) being the limiting case of a perfectly open (or almost closed) system.
The exact universal results were initially obtained for the autocorrelation function of the Wigner time delay, first for orthogonal [6] and unitary symmetry [7], and then for the whole crossover of partly broken TRS [21]. Along these lines, exact formulae were also derived for the distribution function of the partial time-delays [7,21] (defined by the energy derivatives of the S-matrix eigenphases and related to the diagonal elements of Q [22]) as well as that of the proper time-delays (given by the eigenvalues of Q) [17,23]. The developed method is actually flexible enough to also include the effects of finite absorption [23,24] and disorder [24,25]. The most recent overview of the relevant results can be found in [26].
Further progress is possible in the particular case of perfect coupling, when S is uniformly distributed in Dyson's circular ensemble of the appropriate symmetry [27]. Brouwer et al. [28] exploited the invariance properties of the energy-dependent S-matrix in this case to derive the distribution of the whole time-delay matrix, generalising the earlier one-channel result [29] to arbitrary M . The joint distribution function of its eigenvalues, i.e. the proper time-delays, turned out to be determined by the Laguerre ensemble from RMT [28]. This provided a route to apply orthogonal polynomials to compute marginal distributions [30] and various moments [31,32,33,34] or use a Coulomb gas method to study the total density [35]. Very recently, Mezzadri and Simm [36] applied the integrable theory of certain matrix integrals to the problem and developed an efficient method for computing cumulants of arbitrary order. In particular, the variance of the Wigner time delay was found in the simple form, where · · · stands for a statistical average and β = 1 (β = 2) indicates orthogonal (unitary) symmetry. This expression agrees with earlier results obtained at arbitrary coupling in [6,7]. We note, however, that the results for the distribution of Wigner time delay are still rather limited, with explicit expressions being available only at M = 1, 2 [22,37] or in the limit of M ≫ 1 [35]. Taking a semiclassical approach, an early expression for the Wigner time delay was derived through the connection to the density of states mentioned above. The fluctuating part of this quantity can be expressed, like Gutzwiller's trace formula for closed systems [38,16], as a sum over periodic orbits, but now only involving those orbits trapped inside the scattering system [39]. The trapped primitive orbits p and their repetitions m involve their actions S p and stability amplitudes A p,m which in turn depend on the period, monodromy matrix and Maslov index of the orbits. Since the dependence on the energy (and other parameters) is then encoded in the sum over the periodic orbits, this can be used to calculate correlation functions [40,41], in particular, the time-delay variance is expressed as When taking the semiclassical limit → 0, the rapid oscillations in the phase induced will be washed out by the overall spectral averaging unless there are systematic correlations between periodic orbits with action difference |S p − S p ′ | . The repetitions have been ignored in (6) since they are exponentially smaller than the m = 1 term. The semiclassical treatment of sums like in (6) then involves identifying and evaluating correlated pairs of orbits with a small action difference. The simplest pair is to set p = p ′ , known as the diagonal approximation [42], which can be evaluated using the open system version [43] of the Hannay-Ozorio de Almeida sum rule [44] Here the exponential weighting corresponds to the expected number of periodic orbits of the closed system which survive up to a time t, with µ = M/T H being the classical escape rate (see, however, the discussion in [6,45,46]). With TRS, one can also pair an orbit with its time reverse and obtain a further factor of 2, or var(τ W ) = 4 thus reproducing (4) to leading order. This is not surprising, of course, as the quasiclassical limit in the RMT treatment corresponds to the asymptotic case of M ≫ 1 [5,47]. The quantum effects are encoded in the higher order terms of the 1/M expansion, being in general responsible for slowing down the decay law in chaotic systems from the purely exponential one [45]. It is now well understood that the higher order terms can be obtained semiclassically through a systematic expansion of correlated periodic pairs which was first derived for the spectral statistics of closed systems [48,49,50]. The extension to open systems and the time delay was then developed in [51], providing the terms up to M −8 in full agreement with the RMT result. Further calculations along those lines become too involved because of the quickly growing complexity of relevant combinatorics. One possibility to overcome such difficulties might be in doing the semiclassical approximation at the level of the generating functions [52] used to derive the corresponding σ-model of RMT [6,7], the idea that has already proven to be success in establishing the full equivalence with RMT for the closed systems [53]. However, we will exploit another route here. An alternative starting point relies on the van Vleck approximation of the propagator that gives the semiclassical approximation for the scattering matrix elements [54,55,56,57] in terms of trajectories that connect the corresponding (say, input i and output o) channels The scattering trajectories now have different stability amplitudesÃ γ which do not explicitly depend on the duration T γ of the trajectories but still include a phase factor. Taking the energy derivative and noting that ∂S γ /∂E = T γ , one obtains the approximation for the Wigner time delay [46] where we have assumed that the amplitudes vary much more slowly than the phase. The diagonal approximation γ = γ ′ can again be evaluated with a sum rule [58] which directly gives the average time delay.
In the presence of TRS one can also partner a scattering trajectory with its time reverse if they start and end in the same channel. However there are further correlated pairs of trajectories, both with and without TRS. They can be related to correlated pairs of periodic orbits by formally cutting the periodic orbits open and deforming them into scattering trajectories [58,59,60]. It is worth stressing that the second form of the time-delay matrix in (1) turns out to be particularly useful for the semiclassical treatment. One can show in this way [61] that all further contributions to the average time delay cancel leaving only the diagonal terms intact. Furthermore, one can also derive (5) semiclassically from (10) by recreating the sum over periodic orbits in (5) from correlations of scattering trajectories that approach the trapped periodic orbits and follow them for many repetitions. This provides a duality between the two approaches and allows access to a range of statistics beyond the average time delay [61]. For example, the moments of the proper time-delays can be found to leading order in M −1 by using energy-dependent correlators of the S-matrix elements and performing a semiclassical expansion for the later [62]. The next two orders can similarly be obtained using a recursive graphical representation of the semiclassical diagrams of correlated trajectories [63], showing an agreement with RMT [32], though the energy dependence, which is later differentiated out and removed as in (1), complicates the treatment drastically.
In this paper, we derive yet another semiclassical approximation to the time delay which avoids using such an energy differentiation in the first place and significantly simplifies the semiclassical calculations. It builds upon the resonant representation [64] of the matrix elements Q cc ′ = b † c b c ′ as the overlap of the internal parts b of the scattering wave functions in the incident channels c and c ′ . We note that with such a factorised form, calculations of the moments of the Wigner time delay have some resemblance with those of the conductance and shot noise in the Landauer-Büttiker formalism [65]. The statistics of the latter in chaotic cavities is determined by the Jacobi ensembles of RMT [27], the exact RMT predictions for their first two cumulants being derived in [66]. These results agreed with those from the semiclassical approach to quantum transport problems developed in [59,60]. Here, we provide for the first time a similar semiclassical justification of the Laguerre ensembles of RMT. In particular, we derive the exact expression (4) for the variance of Wigner time delay semiclassically at any M both for unitary and orthogonal symmetry (along with the other second moments of time-delays). We also extend previous work [62,63] on establishing the universality for the moment generating functions.
In the next section, we formulate our starting point and develop a semiclassical approximation to the internal parts b c and establish diagrammatic sum rules. This representation is then applied in section 3 to derive the second moments of time-delays at arbitrary M . Section 4 deals with the moment generating function of the proper time-delays and works out the algorithmic approach for the first five terms in a 1/Mexpansion. Section 5 summarises our findings. Finally, we provide several Appendices with more technical details of our calculations, including sums for systems with TRS and a comparison to previous approaches, which we believe may be helpful for further development and applications of the method.

Resonance scattering approach
In the general scattering formalism, the resonance part of the scattering matrix can be represented in terms of the effective non-Hermitian Hamiltonian as follows [18,19]: Here, the Hermitian N × N matrix H corresponds to the Hamiltonian of the closed system, whereas the rectangular N ×M matrix V consists of the decay amplitudes that couple N discrete energy levels to M decay channels. These amplitudes are commonly treated as energy-independent quantities, which can be justified by considering resonance phenomena away from the open channel thresholds. Substituting (12) into (1), one can readily find the representation of the Wigner-Smith matrix [64], In such an approach, the energy dependence enters only via the resolvent (E − H eff ) −1 that describes the propagation in the open system governed by H eff . The N -component vector b c (i.e. the c-th column of b) can therefore be treated [64] as the intrinsic part of the scattering wave function initiated in the channel c at energy E. The diagonal elements q c = Q cc are then given by the norm of b c , providing the interpretation as the average time delay of a wave packet in a given channel [3]. Taking the sum over the diagonal elements, we find that the Wigner time delay is given by the total norm of the internal parts This norm is also known as the dwell time, thus the two time characteristics coincide in the resonance approximation considered. The factorized representation (13) already proved [17,23] to be successful for deriving the exact (RMT) distributions of the proper time-delays (the eigenvalues of Q) in chaotic cavities. On the other hand, the partial time-delays are defined by the energy derivative t c = dφ c /dE of the scattering eigenphases, φ c . The exact distribution of the partial time-delays was found in [7,21]. Since the order of the diagonalisation and energy derivative is reversed for the proper and partial time-delays, they follow different statistics [28,30]. In particular, for perfectly open cavities (the case of interest here) the density of rates t −1 c reduces to χ 2 Mβ distribution characterised by the meant = T H /M =τ W and the variance This should be compared with expression (4) for var(τ W ) that contains an extra factor 2 M+1 , thus reflecting the self-averaging property of the linear statistic (14) and diminishing correlations when M grows. It is worth noting that the covariance of two partial time-delays can also be computed exactly with the help of their joint distribution found in [22], with the explicit result being .
More importantly, it was also shown in [22] that at perfect coupling the statistical properties the partial time-delays become equivalent to those of the diagonal elements of the 'symmetrised' Wigner-Smith matrix, Q s = S −1/2 QS 1/2 , introduced and studied in [28]. For the unitary case of systems without TRS, S becomes statistically independent of Q [28] and the symmetrisation process does not change the statistics of the diagonal elements so that q c also follow the same distribution, in particular, var(q c ) = var(t c ). This is not the case for the other symmetry classes. However, traces of Q and its powers are insensitive to such a symmetrisation, giving the identity Equations (15) and (16) therefore readily yield the variance of the Wigner time delay in the form of (4). Likewise, we can arrive at the same result by using the variance [32] and covariance [33] of the proper time-delays. For later use we quote the relevant expression for the second moment [31,32] The semiclassical approximation for the internal parts b c developed below determines these time-delay quantities in terms of certain scattering trajectories and thus allows us to study the universality of RMT predictions for individual systems.

The semiclassical approximation
In this section we derive a semiclassical approximation for the Wigner time delay that is based on representation (13) of the Wigner-Smith matrix Q. It is the third semiclassical approximation after (5) and (10). The starting point is the Green function G(r, r ′ , E) of the open cavity with an arbitrary number of leads. Its semiclassical approximation is a sum over all trajectories from r ′ to r [38,16] Here, S γ = γ p dq is the action along the trajectory γ, ν γ the number of conjugate points, and v ′ γ (v γ ) is the speed at initial (final) point (for a cavity without potential v ′ γ = v γ ). Furthermore, M γ denotes the stability matrix that describes linearised motion near the trajectory. It connects perpendicular deviations from the trajectory at the end point to those at the initial point Formally, the effective Hamiltonian H eff of an open cavity is an infinite dimensional operator [67,68,69], corresponding to N → ∞ of the matrix truncation in (12) and (13). The position representation of the resolvent (E − H eff ) −1 can then be identified with the Green function G(r, r ′ , E), whereas V corresponds to a projection onto the transverse wavefunctions in the leads such that [70] The integration is over the cross section at the beginning of the lead that contains the c-th incoming mode, v ′ = m k = m k 2 − (nπ/W ) 2 is the longitudinal velocity (m is the mass) and W is the lead width. The corresponding transverse wavefunction is Note that 1 ≤ c ≤ M labels the modes in all leads, whereas n labels the modes in one particular lead, so the choice of the lead and n depend on c.
The semiclassical approximation for the internal part r| b c follows by evaluating the integral in (21) in stationary phase approximation. After writing the sine in (22) as sum of two complex exponentials, the stationary phase condition reads wheren = ±n. This fixes the starting angle of the trajectories (sin θ = p ′ y /p ′ ) entering the cavity. Performing the stationary phase approximation results in where the sum runs over all trajectories that enter the cavity with the angle fixed by (23) and end at r. The amplitudes are given by where µ γ is the number of conjugate points for neighbouring trajectories with the same entrance angle. The expressions (24) and (25) allow us to represent the elements of the time-delay matrix (13) in terms of the trajectories specified above. In particular, the semiclassical approximation for the Wigner time delay follows from (2) and (13) as where the integral is over the interior of the cavity. This is the new representation for τ W that serves as the starting point for the semiclassical calculations in this article. We first apply (26) to calculate the mean time delay. The approximation sums over pairs of trajectories that contribute with highly oscillatory terms. After spectral averaging most terms can be neglected and the only remaining terms are from pairs of trajectories that are correlated. These pairs will be discussed in the following.
The trajectories involved in (26) are similar to those that occur in the semiclassics of the current density [71] which in turn is related to the survival probability [72,73].

Diagonal approximation for the mean time delay
The leading contribution to the average time delay comes from the diagonal approximation where γ ′ = γ: The evaluation of this expression requires a sum rule for the type of trajectories in (27); see also [74] for related sum rules. To obtain the sum rule, we fix one of the leads and consider the probability density that trajectories starting in the opening with angle θ and energy E will arrive after time T at point r, where the initial conditions of r(T ) are determined by y ′ , θ and E. The integral can be evaluated in local coordinates that are parallel and perpendicular to the trajectory, This sum runs over all trajectories and has wild oscillations which can be damped by a conventional smoothing of the delta-function δ → δ ε .
In an open chaotic cavity with area A the asymptotic form of the probability density is P ε (r, T, θ, E) ∼ e −µT /A as T → ∞. Here the exponential term describes the asymptotic escape of trajectories and the 1/A reflects the fact that each end point in the cavity is equally likely. We hence obtain the following sum rule γ(c→r) Note that channel c corresponds to two angles θ which cancels a factor 1/2 coming from (25). With this sum rule we can evaluate the diagonal approximation and obtain This is already the correct expression for the mean time delay. We will now show that off-diagonal contributions leave this result intact.

First off-diagonal corrections for the mean time delay
Off-diagonal contributions come from trajectories that have close self-encounters in which two or more stretches of the trajectory are almost parallel or anti-parallel [48,49,50,58,59,60]. These trajectories have close neighbours that differ in the way in which the remaining longer parts of the trajectory, the links, are connected in the encounter regions. The simplest example is a trajectory with one encounter as shown in figure 1(a). The neighbouring trajectory (dashed line) starts with the same angle θ and arrives at the same end point r, but it traverses the loop in the opposite direction. For this reason these pairs exist only in systems with TRS. For the evaluation of the semiclassical contribution of these orbits we follow [59,60].
The encounter is described in a Poincaré surface of section in the encounter region. The relative distance of the two piercings of the original trajectory through the Poincaré surface is specified by coordinates s and u along the stable and unstable manifolds. This information is sufficient to determine the neighbouring trajectory and the action difference that is given by S γ − S γ ′ = su in the linearised approximation. The duration of the encounter region is specified by requiring that the distance of the Figure 1. two stretches of the trajectory along the stable and unstable directions remain smaller than some small constant c, leading to where λ is the Lyapunov exponent. The summation over the trajectory pairs is done by applying probabilistic arguments. Let w T (s, u) be the probability density in a chaotic system that a trajectory of long duration T has a close self-encounter that is specified by coordinates s and u. The contribution of the trajectory pairs to the average time delay can then be expressed as The density w T (s, u) is given by an integral over the first two link durations where Ω is the volume of the surface of constant energy in phase space. Applying the sum rule (30) and changing the integral over the orbit length T to an integral over the third link duration results in Equation (35) contains a small correction to the sum rule (30) that is necessary when applied to trajectories with self-encounters. Namely, the trajectory time T in the exponent has to be replaced by the exposure time that counts each encounter duration only once. The reason is that if a trajectory does not escape during the first traversal of an encounter region, it will not escape during subsequent traversals of this region. The integral over s and u can now be evaluated by noting that the only contribution that survives in the semiclassical limit is one where the encounter duration t enc in the denominator is exactly cancelled by an encounter duration in the numerator.
In other words, after expanding exp(−µt enc ) in a Taylor series the only term that survives is the linear term in t enc . With ds du exp(isu/ ) = 2π , we finally obtain where we have used T H = Ω/2π and µ = M/T H . The contribution (36) would lead to a deviation from the correct result (27). There is, however, a further contribution of the same order. It arises from the socalled one-leg loops, which are correlated trajectories that have an end point in an encounter region as in figure 1(b). These type of correlations do not occur in transport problems, but they arise when trajectories have one or both their end points in the cavity as for example in problems involving the survival probability or the current density [72,73,71].
The encounter regions of one-leg loops require a different treatment than the usual encounters. It can be shown, however, that this difference can effectively be taken into account by adding a further factor of t enc in the integral over s and u [71]. So the contribution of the one-leg loop trajectories differs from (36) by a missing integration over the third link time t 3 and an additional factor of t enc , The integrals are then evaluated similarly as before and result in This contribution cancels exactly the contribution (36).

Higher off-diagonal corrections for the mean time delay
For higher-order corrections one considers trajectories with arbitrarily many selfencounters. These self-encounters can involve two or more stretches of a trajectory that are almost parallel or anti-parallel, where the latter case requires TRS. One speaks of an l-encounter if it involves l stretches of a trajectory. The types of a trajectory's encounters are detailed in a vector v whose l-th component v l specifies the number of l-encounters. The total number of encounters is thus V = l v l , and the total number of stretches in all encounter regions is L = l l v l . Trajectories with self-encounters have close neighbouring trajectories that differ in the way in which the links are connected in the encounter regions. For a given vector v there are many different configurations in which the encounters and the reconnections can be arranged along a trajectory pair. The number of these structures or families is denoted by N (v). The action difference of a trajectory pair can again be determined in terms of the separation of the trajectory stretches in the encounter regions. There are now altogether L − V pairs of coordinates in the stable and unstable unstable directions (l − 1 pairs for each l-encounter). These coordinates are combined into vectors s and u, and in the linearised approximation one has S γ − S γ ′ = su.
The definition of the encounter duration (32) is generalised to arbitrary lencounters by requiring that the separations of the l trajectory stretches remain smaller than some constant c in the stable and unstable directions, where α labels the encounters, 1 ≤ α ≤ V . One applies again probabilistic arguments to replace the summation over trajectory pairs in (26) by one over self-encounters, where w T (s, u) is the probability density that a trajectory of long duration T has self-encounters that are specified by the vector v and the separations s and u. This density can be expressed by an integral over the first L link durations At a final step, the sum rule (30) is applied. As earlier in (35), it requires replacing the time in the exponent by the exposure time counting all encounter times only once. After replacing the integral over the trajectory time T by an integral over the final link duration, one obtains . (42) In the integrals over the s and u coordinates the only terms that survive the semiclassical limit are those where all the encounter times in the denominator are exactly cancelled by corresponding encounter times in the numerator. The leading contribution thus again arises from the linear terms in the Taylor expansion of exp(−µ α t α enc (s, u)). Using ds du exp(isu/ ) = 2π , the final result reads As in the transport problem [59,60], one can identify simple diagrammatic rules from this result. Each link contributes by a factor 1/M , and each encounter contributes a factor (−M ). The Heisenberg times cancel up to one. Note further that the sum over the channels gives a factor of M that cancels the prefactor 1/M in (2).
As we have seen in section 2.3, there are additional trajectory correlations for the new semiclassical representation (26) of the Wigner time delay due to the one-leg loops that do not occur in the transport problem. In fact, for every trajectory configuration or structure in the above calculation there is a corresponding configuration where the end point is now inside the last encounter. The semiclassical calculation can be easily modified to obtain these additional contributions. The difference to (42) is that the integral over the final link duration is missing, and the integrals over the s and u coordinates contain an additional factor of the last encounter time [71]. These modifications lead to As expected all off-diagonal terms cancel. This calculation allows us to extend the diagrammatic rules: Encounters that include an end point simply contribute a factor of 1.

Diagrammatic rules for higher moments
The previous section has shown one main advantage of the new semiclassical approach. It results in the simple diagrammatic rules that are similar to the ones in transport problems [60], thus strongly simplifying the calculations. As in transport, one can generalise the diagrammatic rules to higher moments. We sketch this by considering the moments of the proper time-delays defined as The representation (13) and the semiclassical approximation (24) lead to Here i 1 , . . . , i n label the n incoming channels and r 1 , . . . , r n denote the n final points. The symbols γ = {γ 1 , . . . , γ n } and γ ′ = {γ ′ 1 , . . . , γ ′ n } stand for two sets of n trajectories, where γ j goes from channel i j to the point r j , and γ ′ j from i j+1 to r j (we identify i n+1 with i 1 ). The total actions of the two sets are S γ = j S γj and S γ ′ = j S γ ′ j , respectively. Dealing with the correlated trajectories that survive the spectral averaging in (46) now involves two trajectory sets, γ and γ ′ . The trajectories in the set γ have encounters in which two or more trajectory stretches are almost parallel or antiparallel. The set γ ′ follows the set γ very closely along the links, but differs in the way those are connected in the encounter regions. One can then replace the double sum over γ and γ ′ by a single sum over γ plus an integral over a probability density for the self-encounters. The result can be split into contributions from links and encounters according to the following diagrammatic rules: • The summation over each incoming channel gives a factor of M .
• Each link contributes a factor of 1/M . • Each encounter gives a factor of (−M ), unless it contains an end point.
• An encounters contributes a factor of 1 if it contains one end point, and a factor of 0 if it contains more than one end point.
There is furthermore an overall factor of T n H , and a factor of 1/M from (45). The rules for the encounters follow from the fact that each end point inside an encounter provides an additional factor of the encounter time.
As an example, we discuss the leading order contribution to the second moment m 2 . In analogy to transport problems we denote the j-th end point by o j instead of r j . The simplest trajectory configuration with an encounter is the one that is schematically shown in figure 2(a). The two trajectories belonging to γ are shown by the full lines and go from channel i 1 to end point o 1 and from i 2 to o 2 , respectively. They have one encounter that is indicated by the circle. The neighbouring trajectories belonging to γ ′ (dashed lines) go from i 1 to o 2 and from i 2 to o 1 , respectively. According to the above rules we obtain the following contribution of this configuration to m 2 This contains (−M ) from the encounter, M 2 from the incoming channels, 1/M 4 from the links, and the overall factor T 2 H /M . The simpler configuration in figure 2(b) does not usually contribute to m 2 since the trajectories of γ ′ (dashed lines) don't connect the correct initial and final points. Note, however, that this configuration is possible if the incoming channels coincide. Then one obtains the configuration in figure 3(a). It contributes at the same order as (47) since it has two links, one encounter and one incoming channel less. There is the further possibility that one of the end points is in the encounter as in figures 3(b) and (c), which again contribute at the same order. All in all the result is This is indeed the leading order term of m 2 in systems with or without TRS. We have obtained it by considering figures 2(a) and 3(a)-(c) as different trajectory configurations that all contribute to the moment. There is even simpler and alternative point of view in which one considers all the contributions in (48) to come from figure 2(a). The diagrams in figures 3(a)-(c) are considered to be limiting cases of figure 2(a) where either the encounter moves into the incoming lead or one of the end points moves into the encounter. These limiting cases can be included by changing the contribution of the encounter.
For the description we adopt the language from transport and call the initial points of γ i-leaves and the final points o-leaves. For the calculation of the moments m n we then have to take into account the following cases: • Encounters of size l can move into the incoming lead when connected directly to l i-leaves. • The o-leaves can be moved into the encounter they are connected to, but each encounter can only take one at a time.
In terms of the semiclassical contributions, in both of these situations we change the rule for the affected encounter to include these cases. In the first case we lose l links, (l − 1) channel summations and the (−M ) from the encounter, and in the second case we lose one link and the (−M ) from the encounter. The contribution of encounter α is hence changed to M (−1 + δ + s α ) with δ being 1 when the encounter can move into the lead (and 0 otherwise) and s α the number of o-leaves attached.
In our example, if we change the encounter contribution for figure 2(a) in (47) from (−M ) to M (−1 + 1 + 2) = 2M then we obtain the same result as in (48). One can also easily see that the off-diagonal contributions for the average time delay all vanish, because the last encounter must have δ = 0 (l is at least 2 and there is only 1 i-leaf) and s v = 1 since it is connected to the outgoing channel. The product of contributions is then 0.
We will show below that these rules can be employed systematically to calculate the second moments of the proper time-delays and the variance of the Wigner time delay, and they can as well be incorporated into the graphical framework [75,62,76,63,77,78,79] to give moment generating functions. Previous approaches to the time delay involved including an energy dependence so that calculating the semiclassical contribution of any diagram required knowledge of its complete structure.
Here instead the semiclassical contributions are much closer to standard transport moments where only the difference in the number of links and encounters matters. The only additional complication is that now some information about the position of the o-leaves is necessary. Determining this is however much less demanding than treating full energy dependence.

Second moments
The diagrammatic rules turn semiclassical calculations into a combinatorial problem of counting trajectory configurations (also called structures or diagrams). Before calculating the second moments of various time-delay quantities, we recall that the structures that are relevant for transport moments [59,80,60] (and hence also for the time delay) can be related to the structures of periodic orbit pairs that contribute to the two-point correlators in spectral statistics [48,49,50].
Correlated periodic orbit pairs are also described by a vector v whose elements v l count the number of l-encounters of the orbits. The total number of encounters is V = l v l while the number of links is L = l lv l . The number of periodic orbit structures with a vector v and a labelled first link is N (v). These numbers can be determined recursively [50].
In the following we deal with systems without TRS. The case of preserved TRS is briefly discussed in section 3.7, being mostly deferred to the Appendices.

Counting diagrams in transport problems
If a periodic orbit pair is cut along one of its links, it can be deformed into a pair of scattering trajectories which contributes to the first transport moment, the average conductance [59]. The cut link becomes two so that conductance diagrams have L + 1 links. The number of structures N (v) of the scattering trajectories is the same as the one for periodic orbits. Since encounters contribute with a minus sign, terms in a M −1 expansion of the conductance can be related to the sum [59] This is the relevant quantity to evaluate. Without TRS, one can show C K = 0 for K ≥ 1 so only the diagonal pair is important [59]. It can be included in the formalism by defining L = V = 0 and C 0 = 1 for it. For the second transport moments (the shot noise and the conductance variance), the semiclassical diagrams involve four trajectories. These trajectory quadruplets can be divided into two groups: d-quadruplets and x-quadruplets. Consider a pair of trajectories γ 1 and γ 2 connecting channel i 1 to o 1 and i 2 to o 2 respectively. If the partner trajectory γ ′ 1 also connects channel i 1 to o 1 then we have a d-quadruplet. Otherwise if γ ′ 1 connects i 1 to o 2 we have an x-quadruplet. The remaining partner trajectory γ ′ 2 must connect i 2 to the other outgoing channel. For example, the diagram in figure 2(a) is an x-quadruplet of trajectories meeting at a single 2-encounter while the diagram in figure 2(b) is the simplest d-quadruplet. It is made up of two independent links.
Denote the number of x-quadruplet diagrams corresponding to a vector v by N x (v), and the number of d-quadruplet diagrams by N d (v), then the conductance variance and the shot noise can be related to the sums These sums can again be related to periodic orbit pairs. Without TRS, as detailed in [60], cutting periodic orbit pairs twice along links leads to a d-type quadruplet. In general this can be done in L(L+1) ways since we may cut the same link twice (and the order matters). Note that the final d-quadruplet will have L + 2 links. On the other hand, an x-quadruplet can be created by cutting out an entire 2-encounter from a periodic orbit pair. In general this can be done in 2v 2 ways and the final x-quadruplet will have one 2-encounter fewer and the same number of links as the periodic orbit pair.
With TRS, the cutting is more complicated, but in both symmetry cases the quantities D K and X K can be related to the following sums Without TRS, both are equal to 1 for even K and 0 otherwise; and we have D K = A K while X K = −B K+1 [60]. We further have D 0 = 1 and X 0 = 0.

Grouping diagrams
The sums X K are known, but for the time delays we need to make a further distinction because of the different rules for the o-leaves. We divide the x-quadruplets into two groups: we denote those where both final points are connected to the same encounter byx, and those where they are connected to different encounters by x ′ . For each vector v we count the structures in each group as Nx(v) and N x ′ (v), with and we define the sums Only theX K part will contribute to the second moment of the time delays. This follows from the rules in section 2.5. If the two final points are connected to different encounters then at least one of these encounters cannot be moved into the incoming lead and contributes with a factor of M (−1 + δ + s α ) = 0. We likewise partition the d-quadruplets and define the sums where again we are only interested in theD K part. However, along with Nd(v) and N d ′ (v), we also need to take into account a third group where one or both end points are not connected to an encounter at all. These diagrams involve a direct link between an incoming channel and an end point, while the other trajectory pair can be any arbitrary conductance (first moment) diagram. We then have where the last term is a correction to avoid overcounting the diagram in figure 2(b) made up of two direct links.

Manipulating diagrams in systems without TRS
To obtain some information aboutD K andX K we need to build a recursion relation between them. We start by considering ways in which and-quadruplet can be obtained from other diagrams. For ad-quadruplet both end points are connected to the same encounter (of size l say), and we can add a 2-encounter at the end. Now we have an xquadruplet whose final 2-encounter has two links connected to the same l-encounter.
Next we shrink the connecting links and merge the two encounters. If l = 2 both 2-encounters vanish and we have an arbitrary x-quadruplet with one encounter and 2 links fewer than the originald-quadruplet. This process is depicted in figure 4. If l > 2 there are two possibilities: A link could separate from the encounter which becomes one smaller (l → l − 1) to leave an arbitrary x-quadruplet with the same number of encounters but one link fewer than before. An example is given in figure 5.
Alternatively the encounter could break into two separate encounters each connected to an outgoing channel hence giving a x ′ -quadruplet. This has one more encounter than the originald-quadruplet and the same number of links.
Reversing the three processes, we can obtain anyd-quadruplet in exactly one of the following ways. We may add a 2-encounter to an arbitrary x-quadruplet ( figure 4). Alternatively for any x-quadruplet, we may join the link before an end point into the last encounter before the other end point ( figure 5). We may also join the final distinct encounters of any x ′ -quadruplet by pulling out a 2-encounter to create ad-quadruplet. Accounting for the minus sign from each encounter and the change in L and V we arrive at the relatioñ since X ′ K = X K −X K . Another way to look at this is the following: one can describe an l-encounter as a cyclic permutation (a 1 . . . a l ) where the a j are labels corresponding to the order encounter stretches are visited in the entire diagram [50]. If stretches a i and a j are connected to the outgoing channel, then adding a 2-encounter and shrinking the intervening links corresponds to multiplying (on the left) by (a i a j ). This breaks the l-cycle into a k and l − k cycle. If k or l − k are 1 then a link is separated from the encounter (leaving only a link if l = 2) otherwise the encounter breaks into two.
We can repeat the same process of adding a 2-encounter and shrinking the adjoining links but starting with anx-quadruplet. This again leads to three cases of d-quadruplets with a lower value of L − V : one with a 2-encounter removed, one with an l > 2 encounter reduced by 1 and one where the l encounter breaks into two separate ones. Reversing the steps, for any d ′ -quadruplet we can connect the two distinct encounters before the two end points by pulling out an 2-encounter (which we then remove) to give a contribution of −D ′ K−1 toX K . The minus sign derives from the final diagram having one encounter less. Then, for any d-quadruplet where o 2 is connected to an encounter, we can join the link to o 1 into the encounter giving a contribution of D K−1 − C K−1 . Here we remove the cases where o 2 connects directly to i 2 . Swapping the roles of o 2 and o 1 gives a factor of 2. Finally we may connect the final links of any d-quadruplet with a 2-encounter creating anx-quadruplet with an additional encounter (and 2 extra links) and hence a contribution of −D K−1 toX K . Putting it all together, using (59). Since the only diagram for K = 1 is in figure 2(a), we haveX 1 = −1 andD 1 = 0, so thatX K = −1 for odd K whileD K = −1 for even K and both are 0 otherwise. (Both are also 0 for K = 0).
We also need to consider the cases when an encounter of a diagram can be moved into an incoming lead. This is only possible if both incoming channels are connected to the same 2-encounter. If, for example, anx-quadruplet starts with such a 2encounter, then this encounter can be cut out (by moving the incoming channels to after the encounter) leaving ad-quadruplet with a value of L − V which is one smaller.
[The only exception is figure 2(a) where the removal of the 2-encounter leads to the d-quadruplet in figure 2(b)]. Reversing the process, all suchx-quadruplets can be built by adding a 2-encounter to the front of anyd-quadruplet [or figure 2(b)]. Similarly, one can interchange the roles of thex-andd-quadruplets and create anyd-quadruplet starting with a 2-encounter by adding a 2-encounter to the front of anyx-quadruplet with a value of L−V which is smaller by one. If we denote by an undertilde quadruplets with an encounter that can be moved into a lead then we havẽ

The second moment of the proper time-delays
We now return to the calculation of m 2 . It is based on thex-quadruplets. For each quadruplet we have L + 2 links, V encounters and a factor of M for each channel. Both end points can also be moved into the last encounter, and one encounter can possibly be moved into an incoming lead, giving a combined contribution of where δ = 1 if an encounter can move into a lead, and 0 otherwise. The case L−V = 1 is different, because then we have the diagram in figure 2(a) where the same encounter can receive the end points and move into a lead. Then the brackets (−1 + δ)(−1 + 2) are replaced by (−2). When we sum over all diagrams using the results of section 3.3, using (62) for the δ = 1 case (and further multiply by T 2 H /M = M/µ 2 ), we have where the term −2X 1 gives the leading order result in (48). The final expression in (64) is exactly the RMT result (18) at β = 2.

Variance of the Wigner time delay
This task requires the second moment of a different type, The trajectories belonging to γ ′ connect now i 1 to o 1 and i 2 to o 2 and hence the dquadruplets are the relevant diagrams. We note that for computing the variance we only need to consider connected diagrams, since the remaining diagrams merely cancel the average time delay squared (in fact, the first d-quadruplet in figure 2(b) does this). Compared to the calculation for the moment m 2 in the previous section, the role of thex-andd-quadruplets are simply reversed, the case K = 1 does not contribute, and due to the normalisation of the variance we don't have to multiply by M/µ 2 . The variance of the Wigner time delay is therefore given by the sum This result fully agrees with RMT, as can be seen by setting β = 2 in (4).

The second moments of the partial time-delays
As already discussed in the beginning of section 2, for systems without TRS the statistics of the partial time-delays turns out to be equivalent (at perfect coupling) to those of the diagonal elements, q c . Hence, we can consider var(q c ) which involves pairs of trajectories going to different end points but all starting in the same channel.
Since the incoming channels coincide, thed-andx-quadruplets now all have the same number of free channels and contribute at the same order. The leading order d-quadruplet cancels the mean part squared [as for var(τ W )], leaving This is in line with the RMT result (15) at β = 2, see also Appendix A for an alternative calculation of var(t c ) explicitly.

Time-reversal symmetry
By creating the relations above between thex-andd-diagrams, we could derive the second moments without any recourse to cutting periodic orbit diagrams. This is notably simpler than the results for the standard transport second moments (e.g. shot noise and the conductance variance) [60]. Treating the case with TRS requires, however, a more involved calculation which we pursue in the Appendices. First, we map the problem without TRS to periodic orbits in Appendix B and then evaluate the semiclassical sums which arise in Appendix C. The corresponding sums for systems invariant under time-reversal are worked out in Appendix D. This allows us to finally obtain, at any M , semiclassical results for the second moments of time-delays in the case of orthogonal symmetry in Appendix E. We find identical results to RMT for the variance of the Wigner time delay and the second moment of the proper time-delays, and derive a new result for the variance of the diagonal elements (which in the case of orthogonal symmetry is not the same as the variance of the partial time-delays). We also demonstrate the equivalence of the semiclassical approach for the time-delay problem developed here to previous treatments in Appendix F.

Moment generating functions
We now apply the diagrammatic rules established in section 2.5 to derive an expansion in inverse channel number of the moment generating function of proper time-delays, As shown in [75], and further developed in [62,76], the contribution of leading order in M −1 comes from diagrams that can be represented as rooted plane trees and which can therefore be constructed recursively. For example, the diagram in figure 2(a) can be redrawn as the boundary walk around the tree in figure 6(a). Rooted trees can of course be considered as unrooted trees, which we will call subtrees below, appended to a single point. Going beyond the leading order amounts instead to adding subtrees to increasingly intricate base structures, following the formalism of [63]. Leaving the details of the diagrammatic approach to [63], we highlight below how it can be adapted to the semiclassical method presented here.

Subtrees
Starting with unrooted trees, there are two types: one with an excess of outgoing leaves whose generating function we call f and one with a excess of incoming leaves (a) counted inf . For example, removing the top i 1 channel from the tree in figure 6(a) to obtain an unrooted tree, an excess of o-leaves remains forming a subtree included in f . The generating variable will be r whose power counts the number of leaves, which is related to the order of the moment. Breaking the trees at the top encounter node gives the recursion where the first term is an empty tree going straight to an outgoing leaf, the next term are trees with encounter nodes of size l with l ≥ 2 and the last term is the correction for allowing one outgoing leaf to move into the encounter (and replacing the corresponding general tree f ) where we have the factor of r to account for the lost leaf. Summing we get with h = ff . This also gives us the useful relation The trees of typef with an excess of incoming leaves can actually also move into the incoming lead if all thef type trees after the top encounter end directly in an incoming channel. Then we have the recursion where the first three terms again correspond to an empty tree, trees with top encounter of size l and moving outgoing leaves into the encounter while the last term is the new possibility of moving the encounter into the incoming lead. Summing giveŝ Substituting for r from (70) in just the first term on the right gives the simplification from which we can get a quadratic for f ,f and more importantly h with s = r 2 .

Leading order
To get the full leading order moments with generating function G 0 , we need to root our trees by adding an incoming channel to the f trees (and divide by M ) or adding an outgoing leaf to thef trees at the top. The first option then allows the trees to also move into the incoming lead to give Using (73) this is just Alternatively, and more simply, we can place an outgoing leaf on top of thef type trees where the sum is over the additional possibility of placing the new outgoing leaf into the encounter. Either way, the end result is that from the correct solution of (74). This is exactly what was in [62] and equivalent to the previous RMT result [30]. Compared to the energy dependent cubic equations of [62] with corresponding energy derivatives and identity matrix corrections, the result (78) can now be obtained much more simply and directly with our new semiclassical approach for the proper time-delays.

Subleading order
Moving to the subleading order, we can continue looking for the dominant diagrams. The non-vanishing contribution exists only for systems with TRS, whereas for those without TRS the first correction occurs in the second subleading order (see below). As shown in [63], the relevant diagram in the former case has the topology of a Möbius strip that arises after merging the quadruplets with time-reversal partners. Around the Möbius strip there are two types of nodes, those with an even number of subtrees on each side and those with an odd number. We shall denote these as an even node and odd node respectively. For an l encounter there are (l − 1) subtrees of each type (2 stretches are the Möbius loop itself) which can be arranged in l ways with an even number on each side and (l − 1) ways with an odd number on each side. These nodes cannot move into the incoming lead, but directly connecting odd leaves can be moved into the encounter node (one at a time). Using slightly different notation than [63] which is more useful for the higher orders in section 4.4 and Appendix G, the contribution of an even node is while an odd node contributes Both include the link leading up to the encounter node (but not the one leaving as this is included in the next node).
Around the Möbius loop we have an arbitrary number of encounter nodes, but we need to divide by their number because of the rotational symmetry. As such we define a generating function of an arbitrary arrangement of nodes around a loop where we divide by 2 to account for the inside/outside symmetry of the loop. A factor p is included with the odd nodes since we actually need to have an odd number of odd nodes to close the loop properly. This is then achieved by comparing the values ofK 1 at p = ±1, giving which is the integrated moment generating function. For the generating function itself, we differentiate so that by using the solution for h from (74) we get the result which is the same as in [63] but obtained in much more straightforward way, without needing to add and then remove an energy dependence.

Algorithmic approach
To treat higher-order corrections, we can use the algorithmic approach of [79] which makes use of combinatorial base structures built on permutations describing the diagram's edges and vertices. The possible permutations are generated via a computer search, while the permissible edge and vertex components are matched up according to the prescription of the permutation. All that is then needed by the algorithm are the general semiclassical contribution of the possible edge and vertex components which we list in Appendix G.

Conclusions and discussion
An efficient method for the semiclassical calculation of the statistics of time delays in chaotic cavities is developed. The method relies on the resonant representation of the Wigner-Smith time-delay matrix that does not involve taking the energy derivative and can be expressed in terms of semiclassical trajectories that enter the system and terminate inside. Under spectral averaging, the results for time-delay moments are produced by sums over sets of classical trajectories which can be evaluated using simple diagrammatic rules. In particular, for individual systems we establish the universality of the RMT predictions for the second moments, including the variance of the Wigner time delay, at arbitrary number of open channels.
We also significantly advance the computation of the moment generating function of the proper time-delays, for which the first five orders are found in section 4.4. This has been achieved by incorporating the present approach into the algorithmic formalism of [63], leading to much simpler diagrammatic rules for the trees which arise. Although previous semiclassical approaches involving energy correlations can be used to obtain these results, the quadratic subtree equations become cubic making the derivations and semiclassical contributions notably more involved. For the second moment of the proper time-delays and the variance of the Wigner time delay, there are tricks to reduce the difficulty of the previous semiclassical treatment, as discussed in Appendix F. However, the results can only be obtained using the sums evaluated here in Appendix C and Appendix D while the variance of the diagonal elements of the time-delay matrix cannot be treated. The new approach developed here presents a simpler and, more importantly, a unified approach to all the second moments.
But the main advantage of our approach is that the diagrams and their contributions are now more similar to those of quantum transport problems, as developed in [50]. Notably for the unitary case, the second time-delay moments can be calculated without any recursive sums, which is not even possible for their transport analogues (like the conductance variance and shot noise). Transport moments are expressed in terms of the transmission eigenvalues that follow the Jacobi ensemble of RMT [81]. Their lower-order cumulants can be obtained most efficiently by exploiting the connection with the Selberg integral, as first recognized in [66] and further developed in [82,83,84]. Combined with the theory of symmetric functions [85], this led to a powerful method for computing transport moments of arbitrary order [86]. Note that arbitrary moments can also be expressed in terms of recursively generated functions called Weingarten coefficients. Moreover, the corresponding semiclassical diagrams can be mapped to certain types of primitive factorisations [77,78] which in turn match the Weingarten coefficients so that semiclassics and the Jacobi ensemble can be proven to be identical [78]. The time-delay matrix (in its 'symmetrised' form [28]) follows, however, an inverse Wishart distribution and its eigenvalues form the Laguerre ensemble [30]. Intriguingly, the moments of inverse real Wishart matrices, corresponding to time delays with TRS, were recently expressed in terms of deformed Weingarten coefficients [87]. Indeed, appropriately substituting γ = M/2 and σ −1 = T H I/2 into the example formulae in [87], one quickly obtains the second moments in (4), (15) and (18) with β = 1. Therefore, the new semiclassical approach, with its simplified diagrammatic rules, opens the door for a dynamical justification of the use of RMT and the Laguerre ensemble in the time-delay problem.
To this end, we mention that another method was recently put forward by Novaes [88] who suggested to use a matrix model which generates the same diagrams as in semiclassics but can be calculated exactly in RMT. For systems with broken TRS, the method was originally implemented to the transport problem, yielding successfully the exact RMT results for general counting statistics, and then further applied to time delays [34], producing the exact moments of proper time-delays up to 8-th order semiclassically. However, showing the full equivalence for arbitrary moments was not yet possible, mainly due to the diagrammatic complexity induced by energy correlations. We believe that the incorporation of the semiclassical approach developed here into the formalism of [34] is a promising way to go forward in establishing the full RMT-semiclassics equivalence. Further development includes generalizations to the case of preserved TRS as well as other symmetry classes, e.g. Andreev billiards for which the statistics of the time-delay matrix was recently derived in [89].
Throughout this article, we have considered the universal regime described by RMT, which neglects the effects due to a finite Ehrenfest time [90]. The latter is responsible for system-specific corrections which can only be obtained semiclassically [91]; various applications to transport were already discussed in [92,93,94] and more recently in [95,96]. The semiclassical representation developed here is already tailored for taking into account such effects in the time-delay problem, calling for further study in this direction. We also mention the challenge of generalising the approach to treat other 'real-world' effects, e.g. absorption [23] and non-ideal coupling [97].
Appendix A. Variance of the partial time-delays without time-reversal symmetry A partial time-delay may be related to the matrix Q as follows [22] where U is the matrix which diagonalises the scattering matrix S. First, we note that the averages over Q and U can be taken independently for the unitary case [28]. Then, averages over U from the CUE are known both from RMT [81,98] and semiclassically [78]. Combined together, this readily yields the mean time delay For a partial time-delay squared, we use the known average to obtain Substituting here the expressions from (64) and (65), we find which gives a (rescaled) variance of var(t c ) = 1/(M − 1), i.e. identical to (66).

Appendix B. Relation of the second moment diagrams to periodic orbit structures for the unitary case
For the calculation of the second moments without TRS in section 3 we did not need to resort to using periodic orbit structures. Since the numbers N (v) of periodic orbit pairs are known however we may explore the relation between orbits and quadruplets with both outgoing leaves attached to the same encounter.

B.1. d-quadruplets
To obtain an arbitrary d-quadruplet, one can cut any pair of links of any correlated orbit pair (including the same link twice) [60]. If we cut links that leave different encounters however, then the resulting quadruplet will have one o-leaf on each encounter and the diagram's contribution will be 0.
To remain attached to the same l-encounter, we simply need to cut different links leaving that encounter (we also cannot cut the same link twice) which can de done in l(l − 1) ordered ways. The total number ofd-quadruplets for each vector v is then and soD Figure B1. A periodic orbit pair with a 2-encounter connected to two different encounters. Shrinking the links between them and the 2-encounter leads to a single larger encounter.

B.2. x-quadruplets
To obtain an arbitrary x-quadruplet, one can cut out any 2-encounter [60]. But again we only need to consider the case when both outgoing leaves are connected to the same encounter as all other cases cancel. Let us first consider instead the opposite case where the outgoing leaves are connected to different encounters and we have an x ′ -quadruplet. The parts of the periodic orbit must be arranged as in figure B1. Say that an l 1 -encounter and an l 2 -encounter are connected to the 2-encounter with encounter links numbered by a i and b i respectively. When the links between the encounters and the 2-encounter are shrunk, a single l-encounter with l = l 1 + l 2 is created. If the reconnection of the original l 1 -encounter is represented by the permutation (a 1 . . . a l1 ) and of the l 2 encounter by (b 1 . . . b l2 ) and the 2-encounter swaps links (a i b j ) then the l-encounter has the permutation Reversing the process, we may pull a 2-encounter out of an l-encounter (and break it into an l 1 -and l 2 -encounter) if we can turn the cycle (1 . . . l) into a 2-cycle followed by a l 1 and l 2 cycle with l 1 , l 2 ≥ 2. For this we can easily check that multiplying (1 . . . l) on the left by (ij) leads to cycles of the correct size as long as i and j are not equal or adjacent (cyclically). There are then l(l − 3) ways of pulling a 2-encounter out of an l-encounter for l > 3. Note that this process adds two links and two encounters so the value of L − V remains constant.
For any periodic orbit pair described by a vector v, an x ′ -quadruplet can therefore be created in l>3 l(l − 3)v l ways while x-quadruplets could be created in the 2v 2 ways of cutting out any 2-encounter. Taking the difference to obtainx-quadruplets leaves possibilities.
Removing a 2-encounter from a periodic orbit leaves a x-quadruplet with L links, V − 1 encounters and a vector with a 2-encounter removed. Recalling the minus contribution of encounters leads tõ

Appendix C. Evaluating unitary sums
Now we turn to evaluating such sums. Since without TRS, (for the general case of K > 0) the sums for both thed-andx-quadruplets reduce to evaluating In fact from the results in section 3.3, H K must be -1 for even K and 0 otherwise (for K > 0). Here though instead we wish to evaluate the sum directly since this will be useful for systems with TRS.
In terms of N (v, l) = lv l N (v)/L we wish to know while these numbers satisfy the recursion relation [50] N where v [a1,...,am→b1,...,bn] denotes a vector that is obtained from v by decreasing all components v ai by one, and increasing all v bj by one. The recursions are derived by shrinking each link of the periodic orbits. Either a k-and l-encounter merge to form a (k + l − 1)-encounter, giving the first term, or an l-encounter splits into a k-and (l − k − 1)-encounter, giving the second term. Applying (C.4) to l = 2, one finds by relabelling the sums. Substituting this into (C.3) yields Since l N (v, l) = N (v) this is just an example of using (C.1).
Substituting next the result for l = 3 leads to Substituting the result for increasing values of l increases the lower limit of the first sum. Since N (v, l) = 0 for l > K + 1, finally we have Now we consider particular terms in this sum. The case where l = 3 (and hence k = 1) can be simplified to in terms of vectors with a lower value of L − V . See [50] for details of the 1-cycles. Likewise the cases where l > 3 and k = 1 or k = l − 2 (whose results must be identical) boil down to (C.10) Finally the more interesting case of when k and l − k − 1 are both at least 2. This gives where we set (l − k − 1) = j. The Kronecker delta arises since the (v l−k−1 + 1) in the original equation (C.8) is only the number of (l − k − 1)-encounters in v [l→k,l−k−1] when l − k − 1 = k. If they are equal it is 1 fewer. Performing the sums over j and k, this result simplifies to Combining the results from the three cases, Since C K−2 = 0 for (K − 2) > 0 then this is simply H K = H K−2 for K > 2. Furthermore, H 2 can be easily checked to be -1 while H 1 is 0, so that H K is always -1 for even K. For K = 0 with C 0 = 1 we could set H 0 = 0 in line with a sum over a 0 vector.

C.1. Orbit interpretation
The terms in the recursion relations above have a simple interpretation in terms of orbits. If we take an arbitrary orbit with L − V = K − 2, we can take any point in any of the L links, combine it with any point in the (L − 1) remaining links or any point either side of the original point and join them together to make a new 3-encounter. This adds 3 links and one encounter. Including the (-1) for each encounter gives (C.14) Alternatively one may take a point from any of the L links and move it to (before or after) any of the L encounter stretches to increase the size of that encounter by 2. This adds two more links and no new encounters Finally one can pick stretches from different encounters, a k-and l-encounter say and combine them into a k + l − 1 encounter. This adds one link and reduces the number of encounters by 1. To count them we can pick any two (ordered) encounter stretches in L(L − 1) ways and remove the l l(l − 1)v l which come from the same encounter

C.2. A further sum
Later we will also need the following sum which we now treat in the same way as H K for systems without TRS. First for l = 2 we have following the steps in (C.12) and since merging the encounters reduces the number of links by 1. Relabelling the sums and substituting into (C.17) leads to The term for l = 3 is since when two encounters are merged the number of links decreases by 1 while when a 3-encounter breaks into links, 3 links are lost in the end. Substituting into (C.17) and relabelling the sum For l = 4 we have (L − 1)N (v, 4) inside the first sum. We can continue to replace terms using (C.4) and the same steps we used for H K as long as we keep track of how the number of links changes for different terms in the recursions. We find that as l increases we continue to have this factor of (L − 1) and the sum reduces to This result can be read off directly from the orbit interpretation if we remember to multiply by the new number of links minus two when making a new 3-encounter and instead by the new number of links minus one in the other cases. We can now simplify the results from the three cases since the L 2 factors cancel. This leaves With C K−2 = 0 for (K − 2) > 0 then J K = J K−2 for K > 2. Since J 2 is -1 and J 1 is 0, then J K is always -1 for even K. For K = 0 with C 0 = 1 we could also set J 0 = 0 in line with a sum over a 0 vector.

Appendix D. Evaluating the orthogonal sums
With TRS, the numbers N (v, l) satisfy slightly different recursion relations [50] N with an extra factor of 2 and a new term for when a link returns to the same encounter in the opposite direction. Performing the same steps as for the unitary case, one finds The first term on the second line derives from the l = 2 case in (D.1) with the 1-cycle treated as in [50]. Since the 2-encounter is removed entirely in the end, a minus sign appears. Otherwise for encounters with l > 2, the encounter is merely made one smaller providing the remaining term on the second line, which is simply H K−1 .
With TRS, we further have so that those terms in (D.2) cancel. All told we have the recursion relation we then get Including the factor of L for J K and performing the same steps we have The delta function arises since the l = 3 case in the recursions starts with a value of (L − 2) as opposed to the (L − 1) for l > 3. This result can be rewritten in terms of known sums To proceed further we use the relation [60] and (D.3) to obtain Using also that we can rearrange the result to and obtain a recursion relation for I K = J K −A K . With the starting values of I 1 = −1 and I 2 = −5 we have a general result of from which we can find J K by adding A K from (D.11). For K = 0 with A 0 = C 0 = 1 we could also set J 0 = 0 in line with a sum over a 0 vector and hence I 0 = −1 in agreement with (D.13).

Appendix E. Second moments with time-reversal symmetry
When we turn to systems with TRS we have the complication that the encounter stretches can be traversed in either direction. As such we can now divide the xquadruplets into three groups:x where the outgoing leaves are connected to the same encounter and those encounter stretches travel in the same direction,x where the outgoing leaves are connected to the same encounter but those encounter stretches travel in opposite directions, and x ′ containing the remaining diagrams with the outgoing leaves connected to different encounters. For each vector v we count the number in the new group as Nx(v) and definê Since all quadruplets belong to one group Using the same notation for the d-quadruplets, the new groupd contains those where the outgoing leaves connect to the same encounter from encounter stretches travelling in opposite directions with number Nd(v). Likewisê

E.1. First relations
As for the unitary case, we can first consider adding a 2-encounter to the end of andquadruplet. Since the relevant encounter stretches are traversed in the same direction, when the links to the new 2-encounter are shrunk, the steps are identical to those in section 3.3 giving the same relatioñ albeit with an extra term arising when we replace X ′ using (E.2). Also nothing changes when we start with anx-quadruplet giving using (E.4).

E.2. Second relation
Now we consider adding the same 2-encounter to the end of and-quadruplet. Since the relevant encounter stretches are traversed in opposite directions, the l-encounter at the end before the new 2-encounter must have size l > 2. To see what happens, it is simplest to express the l-encounter as a permutation, and we need to record that of both the encounter stretches and their time reversals. We use a bar to represent time reversal and hence stretches travelling in the opposite direction. Say that stretches a . We start with ad-quadruplet with both end links connected to the same 3-encounter, but travelling in opposite directions through the encounter.
Appending an additional 2-encounter in (a) is akin to reconnecting the dashed trajectories as in (b). Untwisting the lower solid trajectory leads to (c) while we must also untwist the top loop so that all the encounter stretches end up parallel or anti-parallel as in (d). This diagram is now ax-quadruplet still with a 3-encounter traversed in opposite directions by the end links. Appending another 2-encounter to (d) would effectively reverse the process and simply cancel the one added in (a).
andb were originally connected to the outgoing leaves so that the l-encounter has the permutation (ax . . . ybz . . . w)(w . . .zbȳ . . .xā) , (E.7) consisting of a single l cycle and its time reversal. The entries w . . . z could contain bars while the actual numbering of the elements would correspond to the order of traversal. Adding the 2-encounter between the stretches from a and b corresponds to multiplying later by the permutation (ab), and for the time reversal before by the permutation (āb) leaving us with again a single l-encounter with both outgoing leaves still attached to stretches travelling in opposite directions. Adding a 2-encounter however changes a d-type quadruplet into an x-type and once the connecting links are shrunk we still have the same number of links and encounters. We illustrate this process in figure E1. Performing the same steps to ax-quadruplet we just go back to ad-one so this is an involution andD

E.3. Third relation
So far these relations are not sufficient to determine the four quantitiesD,D,X,X but to proceed we can use TRS to our advantage. With TRS, we can connect the outgoing leaves o 1 to o 2 together, and also the incoming channels i 2 to i 1 and create a periodic orbit. If we start with ad-quadruplet the link formed from connecting the outgoing leaves must return to the same encounter, but now in the opposite direction.
Starting the partner orbit in the same direction it will also follow the link made from connecting the incoming channels in the same direction. Performing the same steps to ax-quadruplet we again have a link connecting an encounter to itself (but returning to the encounter in the opposite direction), as well as one where the orbit and its partner travel in different directions. Given all the periodic orbits with a vector v we can cut any link which connects an encounter to itself (traversed in the opposite direction) and place the outgoing leaves there. Then we can cut any of the remaining L − 1 links, placing the incoming channels appropriately, to create Nx(v) + Nd(v). Fortunately the number of links connecting encounters to themselves (in the opposite direction) corresponds to the third line in (D.1). When we multiply by the required (L − 1) (which cancels with the denominator) this leads tõ by relabelling sums and treating the 1-cycles that arise from the l = 2 term as in [50]. The l = 2 term in the first line gives the first sum on the second line which has a simple orbit interpretation. Given an orbit with a vector v ′ with L ′ links, we can replace any of those links by a 2-encounter that returns to itself, adding two more links and one more encounter. Cutting the new link returning to the 2-encounter, L ′ + 1 other links can be cut to create ax-ord-quadruplet.

E.4. Fourth relation
Likewise, if we connect the outgoing leaves of ax-ord-quadruplet we arrive at a periodic orbit where a link connects an encounter to itself, but now in the same direction. From the periodic orbits with a vector v we therefore cut any link which connects an encounter to itself in the same direction and then cut any of the remaining L − 1 links to create Nx(v) + Nd(v). The number of links connecting encounters to themselves (in the same direction) corresponds to the second line in (D.1). Again the factor of (L − 1) cancels and following the same steps as in Appendix D.

E.5. Final results
SinceD =X from (E.9) we therefore havê an explicit result for these quantities using (D.13). ForD andX we first consider the differenceX using (E.5) and (E.6) for K > 1. Again sinceD =X, this simplifies tõ . (E.14) Next (X 1 −D 1 ) can be checked to equal −1 so that for K > 0. Note that for K = 0 bothX 0 andD 0 are 0. Substituting into (E.10) gives the explicit formulae 2D K = I K−1 − (−1) K 2X K = I K−1 + (−1) K (E. 16) in terms of (D.13). Useful for the second moments are the sums Finally, we can again consider the case where both incoming channels are connected to the same 2-encounter which can then be moved into the lead so that both incoming channels coincide. As for systems without TRS in section 3.3 such diagrams can be generated by simply appending a 2-encounter to the start of appropriate quadruplets. We therefore again have the relation (62) while for thex andd cases we analogously haveX K = −D K−1 ,D K = −X K−1 , K > 1 . (E.18)

E.6. The second moments
For the calculation of m 2 the contribution of each diagram is as in section 3.4. When we sum over all diagrams (and further divide by µ 2 M ) we have to evaluate , (E. 19) where the leading order term from (48) is the same as for systems without TRS, and is included as the K = 0 term in the final sum. As discussed in section 2.5 this may be viewed as either a combination of D 0 −X 1 −X 1 or simply −2X 1 . The sum over d-quadruplets in (E.19) is again derived from x-quadruplets where the initial 2encounter can be moved into the incoming lead as in (62) , (E.20) sinceD 1 =D 1 = 0. This result again agrees with RMT; set β = 1 in (4).

E.7. Spin-orbit interaction
For the symplectic symmetry class of spin 1 2 particles, the semiclassical diagrams are identical to those for the orthogonal TRS case. Spin orbit interactions are instead included as additional spin propagators along the classical trajectories [99,100]. Each channel or leaf is also split into a spin up and spin down version though observables are appropriately rescaled so that the semiclassical contribution of the diagonal pair to the average time delay remains unchanged. At subleading order, each of the diagrams in figure 1 gain an additional factor of − 1 2 [101] but they still cancel. Continuing to all orders [73] all off-diagonal contributions similarly cancel.
Treating the second moment of the proper delay times, the spin semiclassical contributions of each diagram become more complicated [101] but effectively the leading order term remains unchanged while each higher order gains a factor of − 1 2 . This can be simply incorporated by substituting M → −2M in (E.19) to give µ 2 m 2 = 4M 2 (M + 1)(2M − 1) , (E. 21) and enormously reduce the complexity of the problem. We state the semiclassical result for the correlator [61] C(ǫ) = 1 + 2 − β βM in terms of a sum over periodic orbits (since these are cut once to obtain the first moment diagrams) and we include the diagonal term as a vector with 0 entries. The prefactor accounts for the coherent backscattering diagrams: With TRS, when the incoming and outgoing channel coincide, the time reversal of the partner trajectory can also be paired with the original trajectory. Differentiating (F.6) Since σ l σ = l lv l = L the second term is simply 1 and so These are sums we have already evaluated and hence Without TRS, β = 2, the sum A K + C K − J K = 2 for all even K (including K = 0) and hence in agreement with (64). With TRS, β = 1, we have the sum in agreement with (E. 19). Even for the unitary case we are forced to use the semiclassical sums from Appendix C and the calculation is only tractable in this way because we used (F.3). Nonetheless, this shows the agreement between the new semiclassical method presented here and previous approaches.

F.2. The variance of the Wigner time delay
For the variance of the Wigner time delay, to avoid treating energy dependent correlations between quadruplets of scattering trajectories, we turn to the expression in terms of periodic orbit correlations like (6) as discussed in section 1. In particular we may derive an expression for the two point correlator of the time delay at different energies as a second differential (c.f. the appendix in the preprint version of [60]). This For edges without TRS and traversed in the same direction on either side we also need to consider odd nodes with an excess of f orf type subtrees which can also touch the incoming lead.
The edge contributions are then [79] E For the vertices of degree k, we label the edge stumps by the components of a vector b. If adjoining components are identical, we need an even number of subtrees in that sector, otherwise an odd number with an appropriate excess of one type of subtree f orf . The normal contribution is where q is the number of times i follows i in the sequence b (taken cyclically) and p is the number of times o follows o. Next, any of the f type trees can connect directly to an outgoing leaf, which we can move into the encounter at the vertex. This give a further contribution of Finally, if all the sectors are odd with an excess off types subtrees so that p = k, the encounter can move into the lead giving an extra contribution of (G.10) All combined we have These contributions can be plugged into the algorithm of [79] to give the moment generating functions detailed in section 4.4.