Wave diffusion and mesoscopic dynamics, towards a universal time-dependent random scattering matrix

We concern ourselves with the prediction of mesoscopic wave phenomena from statistical knowledge of classical trajectories. A diffusing particle picture for the flow of mean probability in chaotic systems is used to estimate dynamical features of mean square time-domain S matrices for waves coupled in and out through one perfectly open channel. A random process with that mean square, and with the additional constraint of unitarity, is then shown to lead to plausible S matrices with familiar mesoscopic wave dynamics. Features that are generated by this procedure include enhanced backscatter, quantum echo, power law tails, level repulsion and spectral rigidity. It is remarkable that such rich behaviours arise from such simple constraints. We conjecture that a generalization to n × n S matrices would exhibit behaviour identical to that of a Hamiltonian taken from the Gaussian Orthogonal or Unitary Ensembles (GOE or GUE) depending on its symmetries. Further constraining the S matrices to reproduce non universal aspects of classical dynamics, (known short time behaviours, periodic orbits, stable islands…) may generate mesoscopic wave features of such systems.


3
Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT of S matrices that includes such correlations, perhaps one that could incorporate knowledge of non-universal features where present.
The present study is inspired in part by a recent contribution from Lev Kaplan [8] which shows how incorporating exact short time wave behaviours into an otherwise ergodic ray dynamics, allows to derive nonuniversal features of the wave dynamics. It is also part of a project in structural acoustics in which cost-effective short-time direct numerical solutions (DNS) of many degree of freedom linear wave systems are to be used to make predictions for energy flow and field statistics over late times [6,11]. Here though, we eschew such information, and seek universal wave dynamics assuming only ergodic ray statistics. We make no attempt to incorporate deterministic information from short times and instead merely assume a simple diffusion of incoherent rays as would describe stochastic particles. It is thus suggested that the appropriate prior for the E dependence in an ensemble of S matrices can be taken from a diffuse thermodynamic picture of ergodic ray dynamics akin to that of room acoustics or statistical energy analysis [12].
The following section reviews the standard picture [3] for the relation between S matrices, scattering channels, Green's functions and Hamiltonians. It is followed by a section introducing an approximate scattering matrix s(t) that is motivated by a semi-classical picture of rapidly mixing rays. It is next shown how that S matrix may be made unitary in a minimally invasive manner. These ideas are then implemented numerically, and the resulting S matrix and Green's function are found to exhibit mesoscopic features familiar from Gaussian ensembles of Hamiltonians.

Green's functions, effective Hamiltonians and the S matrix
Most discussions of RMT wave dynamics begin with a single quantum particle in a closed system governed by the Schrödinger equation with an n × n Hermitian Hamiltonian matrix H with large n. It is associated with an n × n Green's function G(t) governed by with being the unit step function. In the frequency domain As H is Hermitian, G(ω) is Hermitian except at resonances. (Tildes are dropped except as needed for emphasis.) If the system is time-reversal invariant, such that H is symmetric and real, so is G(ω) is symmetric, and real except at resonances. These equations may also be understood as describing a classical wave system 1 , as shown in the appendix. In RMT, we typically take H to be a member of the GOE or GUE. 1 Classical wave systems are, in their time-independent form, essentially identical to their quantum counterparts; one need only map the eigenvalues; their eigenfunctions are identical. In the time-domain, differences can be more problematic. For narrow-band processes at a central frequency , the classical wave equation ∂ 2 t ψ = Oψ with ψ = (t) exp (i t) and slowly varying, reduces to a Schrödinger equation, 2i ∂ t = [O + 2 ] . A more exact mapping is discussed in the appendix.

4
Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT Following others [2,3] consider a set of m orthogonal functions W a a = 1, 2, . . . m representing source and receiver channels. Define a scattering matrix S ab , an m × m matrix function of frequency (m n) by H − i c W c W c * is called an effective Hamiltonian [2,3] equal to the bare Hamiltonian H plus an anti-Hermitian operator −i c W c W c * that represents losses into the channels. We define G to be the m × m measurable part of the Green's function so-called because this is what we would measure (except for two trivial factors of the infinitesimal coupling strength) as an S matrix, were the fully open channels to be made infinitesimally open. It is this quantity which is directly measured in an ideal ultrasonic experiment [2], but only inferred in microwave and nuclear physics. It may be recognized as identical to the Wigner reaction matrix, though here we emphasize instead its connection to the Hamiltonian H. It is concluded that and That S(ω) is unitary S † S = I follows from the Hermiticity of H. More abstractly, and without reference to a Hamiltonian, it follows from conservation of probability. Both G and S have timedomain versions given by inverse Fourier transforming: G(t) = G(ω) exp (−iωt) dω/2π. For a classical wave system G(t) is purely imaginary in the time domain (see appendix), and S(t) is purely real. For a quantum system we generally expect S(t) and G(t) to be complex.

The statistical semi-classical S matrix
Kaplan [8] proposed to calculate S over short times by direct numerical simulation or by semiclassical ray methods, and where necessary extend S(t) to later times by appropriate Gaussian statistical assumptions. Then G, or other quantities of interest, could be constructed by an iteration procedure equivalent to the summation in equation (8). To explore that idea here we simplify the problem and neglect all knowledge of short time non-universal behaviour so as to seek a universal dynamics. Diffuse field arguments, i.e., stochastic ray arguments, regarding the time-dependent transport of incoherent energy (probability in quantum mechanics) motivate a statistical semi-classical matrix s(t). It may be, but it is not explored here, that the same s(t) may be motivated by other arguments, without reference to rays. In the case of m = 1 channel in a single rapidly mixing system (figure 1), the diffuse field arguments call for an S matrix that is the function s(t) whose expected envelope is Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT Figure 1. A single open channel into and out of a single rapidly mixing structure.
the overbar representing a short time or ensemble average. Exponential recapture of probability is a consequence of a hypothesis of well mixed rays, in which a ray revisits the source where it is then recaptured with a constant probability. Each ray has an equal chance to be recaptured, independent of its history. Thermodynamic arguments for a classical diffuse field without mesoscopic considerations establish that σ may be no greater than 1/T H [13], and equals

1/T H only if the channel is fully open.
A fully open channel [3] corresponds to 100% prompt transmission between channel and random structure. If the channel is not fully open, the recapture rate would be slower, σ would be smaller, and the process s(t) would have correlations related to multiple revisits to the source. By presuming fully open transmission we are better justified in neglecting correlations in s(t). This semi-classical diffuse scattering matrix s(t) is, in the time domain, a real (we presume a classical wave system with a real field; see appendix for one way to cast a time-domain classical wave equation into Schrödinger form. Preliminary tests have concluded that assuming s complex leads to the same results) centred Gaussian process r(t) times the envelope √ E(t) where r is taken to have unit mean square and Gaussian correlations consistent with an assumed frequency band of interest. The assumed lack of long-range correlations within s(t) is consistent with the statistical semi-classical diffuse ray assumption. It is inconsistent with known mesoscopic wave phenomena, but if the channel is fully open such that rays can visit the detector at most once, we imagine that such correlations are minimal. This is especially true if there are many open channels, as most rays will be recaptured by a channel other than that of their source. This construction lacks the unitarity that s ought to have, the unitarity that is needed before s can be incorporated successfully into the iteration (8). Thus the next section discusses how it may be repaired.

Unitarization
The above choice for s(t) is consistent with the diffuse ray assumption, but it is unlikely that its Fourier transforms(ω) will have unit modulus. A good choice for normalization of E(t) can make it unitary on average, but the condition |s(ω)| = 1 for all ω is in general not satisfied by the above construction. Thus we are led to seek a filter h(t) by which s can be rendered unitary, while 6 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT leaving its time-domain envelope √ E as little disturbed as possible. In particular we require that the filter be causal, as S should remain zero for negative times.
The filter h is to be chosen such that |S(ω)| = 1; i.e. |h(ω)| = 1/|s(ω)|. This leaves the phase of h unspecified. The zero-phase choice, h(ω) = 1/|s(ω)|, is not causal and must be discarded. Causality in h implies that it can have no poles in the upper half ω plane. If it is to be minimum phase as well [14] it must have no zeros in the upper half plane. This condition corresponds to choosing h to be as compact as possible in time so as to distort E(t) as little as possible, and retain (as much as possible) the time-domain information explicit in s(t). Under this condition, h's logarithm is also causal, and obeys a Kramers-Kronig [14,15] relation. We may write where the unknown phase of the filter φ is given by a Hilbert transform on log |s|. Equivalently and more explicitly, a(ω) is given by While the modulus of the filter is determined locally |h(ω)| = 1/|s(ω)|, the phase of the filter depends on |s(ω)| at distant frequencies. This recipe completely specifies the filter h; the resulting S = s h is unitary and causal by construction, and in a sense resembles s as much as is feasible.

Numerical implementation
This procedure has been implemented numerically using discrete Fourier transforms on a long, but finite, time interval. An N (equal to 2 18 ) element array of Gaussian white noise s t (t = 1, 2, . . . , N) is generated by assigning a real centred Gaussian random number to each element of the array. It is multiplied by an envelope function √ E = exp (−t/2T H ), with Heisenberg time T H equal to 5000 independent of frequency. An example of the resulting s t is plotted in figure 2. The envelope has been normalized to s 2 t = 1. The absolute value of its spectrum |s f | is plotted in figure 3, where it may be seen that the random process is not unitary.
A complex filter h f (real in time domain h t ) is then constructed by means of a discrete Fourier transform version of equation (13). Causality is a subtler issue for finite length arrays and discrete Fourier transforms than it is on an unbounded time domain and for the ordinary Fourier transform. A discrete Fourier transform pair is of the form with integer t and f . The integer index f may be identified with frequency ω = 2π(f − 1)/N. Frequency ω/2π The analogue of equation (13) is where t is the analogue of the unit step function (t) and is defined by This choice for t permits us to establish Re(ã f ) = −log |s f |; i.e., |exp (ã f )| = −1/|s f |.
So thatS f =s f exp (ã f ) is unitary: |S f | = 1. Frequency ω/2π Real and Imaginary parts This operation is implemented numerically and applied to the s t as generated by the random number generator and the envelope √ E. The modulus |S f | of the spectrum of the resulting process S t was found to be unity to within a part in 10 7 . Its real and imaginary parts, and its phase, are shown in figures 4 and 5. Figure 6 demonstrates the uniform distribution of the phase of S f . Were the histogram to be taken from an ensemble of S matrices, such uniform phase would be an indication that the ensemble is invariant under unitary transformations, the prime requirement of ensembles of S at fixed frequency [10].
S t is shown in the time domain in figure 7, where it may be compared to the original signal s t of figure 2. As desired, S t has an envelope similar to that of s t . Figures 8 and 9 compare the smoothed (over a time interval of 100) squared S t to the smoothed squared s. We see that S 2 has an initial value approximately twice that of s 2 , and an initial exponential decay rate twice that of s 2 . A factor of two is precisely the elastic enhancement (also called coherent backscatter or weak localization) effect for time-reversal invariant systems. This factor is predicted by RMT, and also by simple arguments about rays and their time-reversed images that apply even without RMT [16].
The envelope of S t has a longer tail than does the simpler s t . Its intensity decays like a power law (figure 10) with an exponent close to the -5/2 it would take in the GOE case [2,3,17] and the −3 that it would take in the GUE case. The S(t) 2 generated here is compared to the integral form expressions for S(t) 2 as calculated by averaging over a Gaussian orthogonal ensemble of Hamiltonians (figure 10 solid curve) [3,17]. There are small differences, especially at late times, but they may perhaps be ascribed to errors related to our use of a finite time interval. S 2 (t) is the inverse Fourier transform (with respect to ) of the ω averaged S(ω)S * (ω + ) , and so this plot also describes the frequency-frequency correlations in the repaired S.
Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT Figure 6. A histogram of the phase of the repaired S matrix. Each count is the number of discrete frequencies (out of a total of about 25 000) at which the phase is within the given bin. It is striking that a simple and naive diffuse ray picture, as modified by the concepts of unitarity and minimum phase, imply an S matrix that agrees so well with predictions from the GOE, and in particular manifests both enhanced backscatter and power law tail.
We furthermore note that many of the detailed fluctuations in s 2 t are reproduced in S 2 t ; thus the minimum phase repair of s has honored even minor features in the guessed dynamics. This is encouraging for the conjectured extension to s t 's which include some non-universal information, based perhaps on direct reactions or short-time direct numerical simulations. It may be noted that the repaired S has twice the energy at early times, and twice the decay rate, of the original signal. We also note that the repaired S reproduces many of the fluctuations of the original signal. It may be seen that the exponent in the power law tail takes a value close to that predicted for GUE and GOE systems. The smooth solid curve is from an evaluation of the Verbaarschot, Weidenmuller and Zirnbauer (VWZ) integral [3] for the case of the GOE (courtesy T Gorin) with an assumed Heisenberg time of 5000, and a vertical shift for fit.
Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT Figure 11. The constructed function iG (t).
We also examine features of the G consequent to our S, and compare them with those that are expected for random Hamiltonians. The associated Green's function G is given by the power series, equation (8). In the time domain, this corresponds to successive convolutions with S t . In the frequency domain, it involves simple multiplication. Slightly simpler than summing the power series is to directly evaluate the ratio G = S/(I + S) = S − S 2 + S 3 − S 4 + · · · = (1/2)(I − iG). G = i(2G − 1). The presence of discrete frequencies at which S = −1, can complicate evaluation of the ratio S/I + S. One method for coping with the complication is to multiply S t by an artificial absorption exp (−ηt) with a judiciously chosen value of η. A Fourier transform is then taken and a sub-unitary quantity S D f is obtained for the absorbing system. The ratio G D f = S D f /I + S D f , may be inverse Fourier transformed to reveal G D t in the time domain. The artificial absorption is then removed by multiplying G D t by exp (+ηt) to construct the correct time-domain G t . As we use discrete Fourier transforms here, there is a danger when multiplying frequency-domain functions of incurring 'wrap-around'. Time domain elements S t for t > N/2 can behave as if they correspond to negative times. Thus η must be chosen sufficiently large. On the other hand, too large an η may lead to numerical imprecision. G t as constructed by the above procedure with η chosen at 10/N, is shown in figure 11. Figure 12 depresses the random character by squaring and smoothing and shows a remarkable evolution over times less than T H , in which the smoothed energy increases by a factor of about 3/2. That factor is a well-known property of random matrix systems, as the GOE enhanced backscatter factor of two at short times becomes three at times of the order of T H or greater [18]. (In the GUE, the quantum echo corresponds to a transition from an elastic enhancement of one at short times, to a late time value of two.) It is remarkable that a procedure based on such simple assumptions leads to the quantum echo, a phenomenon that is often considered a consequence of spectral rigidity [18].
The present procedure has made no reference to modes, but has nevertheless led to a spectrum, figure 14, with resonances. Mathematically it is no mystery; these resonances are merely those discrete and generic points where S(ω) = −1. The inverse of the mean spacing between these resonances ( f = 1/4701, see the staircase function in figure 15) is close to the initially specified Heisenberg time T H = 5000. This also is no mystery, as the group delay (which corresponds to the rate at which the phase of s f advanced and thus the approximate rate at which S f successively achieves a value of −1) of the random process is simply T H . Frequency ω/2π Figure 14. The imaginary part of −G f shows distinct resonances, at the frequencies where the phase of S passes through an odd multiples of π. Side lobes, and resonance widths, are related to the time-windowing used on G t , a cosine bell tapering to zero for times greater than 2 17 . Figure 15 shows what is commonly termed the 'staircase function'. It shows that the levels (as determined from those points at which the phase of S f is π) are distributed uniformly and with some degree of rigidity. The nearest neighbour level spacing histogram is shown in figure 17, where it may be seen that this statistic corresponds neither to the GOE nor to the GUE. Figure 18 shows a histogram of the strengths of the resonances. These quantities would be, were the system a reverberant cavity, the squares of the mode amplitudes at the antenna. Frequency ω/2π The staircase function Were the reverberant cavity to be time-reversal invariant, we would expect GOE to apply and the modal amplitudes to be distributed as real Gaussian random variables. Were it to be the GUE that applied, the modes at the antenna would be complex Gaussian random variables, and the expected distribution of its absolute value squared would be that of the sum of the squares of two real Gaussian variables. It may be seen that the observed resonance strength distribution corresponds well to neither that of the GUE nor the GOE. Unfortunately G 11 (ω), and S 11 (ω) are ill-suited for unambiguously determining the levels, because weak resonances are difficult to figures 15-18 are based on informal identification of the levels and are uncertain; some levels may have been missed, noise may have been erroneously identified as levels. We therefore do not attempt a detailed analysis of the observed eigen-statistics. While there are unmistakable differences between the level statistics seen here and those of the GOE and GUE, we nevertheless do observe level repulsion and spectral rigidity. It would perhaps be surprising if S did show perfect GOE or GUE statistics, as the weaker resonances are of little importance in G 11 and S 11 and so can be only weakly dictated by the considerations of unitarity and group delay that informed the present construction.

Summary
It is concluded that the single channel case reproduces many of the dynamical predictions of RMT. That it does so suggests that RMT dynamics follows from the simple constraint that the otherwise incoherent exponentially decaying S matrix must be unitary. It is surprising that such a simple constraint leads with such seeming accuracy to familiar mesoscopic probability flows. The repaired unitary S shows elastic enhanced backscatter, the quantum echo, the power law tail, and at least some of the level statistics of RMT. Reproducing the dynamical behaviours of RMT without reference to an internal Hamiltonian is of considerable interest [10,19]. The 1 × 1 case treated here is, as discussed above, ill suited for fully dictating level statistics and it is therefore conjectured that the multi-channel case may lead to more precise correspondence with predictions based on RMT Hamiltonians [3,17]. There are systems for which RMT does not apply, and for which non-universal mesoscopic predictions beyond those of RMT are needed. These include those with short time correlations (e.g. short periodic orbits), direct reactions, and those with substructures that mix only slowly. The initial model of purely incoherent decay used to motivate equation (9) could be modified to include such information. That the procedures used here may lend themselves to such generalization is exciting.
(A being real and skew, B being real and symmetric, H is Hermitian) which may be combined to derive a classical wave equation for ψ 1 ∂ 2 t ψ 1 = A∂ t ψ 1 − B 2 ψ 1 allowing us to identify A with −m −1/2 c m −1/2 , and B 2 with m −1/2 k m −1/2 and ψ 1 with m −1/2 q. The Schrödinger Green's function is a 2 × 2 matrix of operators The (11) element of this is i times the time derivative (∂ t → −iω) of the classical wave Green's function. Thus by taking W a = {w a (x), 0} T , we find