Intensity statistics inside an open wave-chaotic cavity with broken time-reversal invariance

Using the supersymmetric method of random matrix theory within the Heidelberg approach framework we provide statistical description of stationary intensity sampled in locations inside an open wave-chaotic cavity, assuming that the time-reversal invariance inside the cavity is fully broken. In particular, we show that when incoming waves are fed via a finite number $M$ of open channels the probability density ${\cal P}(I)$ for the single-point intensity $I$ decays as a power law for large intensities: ${\cal P}(I)\sim I^{-(M+2)}$, provided there is no internal losses. This behaviour is in marked difference with the Rayleigh law ${\cal P}(I)\sim \exp(-I/\overline{I})$ which turns out to be valid only in the limit $M\to \infty$. We also find the joint probability density of intensities $I_1, \ldots, I_L$ in $L>1$ observation points, and then extract the corresponding statistics for the maximal intensity in the observation pattern. For $L\to \infty$ the resulting limiting extreme value statistics (EVS) turns out to be different from the classical EVS distributions.


I. INTRODUCTION
This work aims to contribute towards understanding the statistics of intensity of a monochromatic wave field inside an irregularly shaped enclosure (cavity) which could be fed with incoming waves through M open channels (antennae), see a sketch in the figure: Figure 1.A schematic sketch of a chaotic wave scattering in a cavity, with ψ l , l = 1, . . ., M being the wave in l th channel/antenna, with x l being the coordinate along the channel.An operator describing wave dynamics in the system decoupled from the channels/antennae is assumed to be effectively described by a large random matrix Ĥ.The M × M unitary scattering matrix Ŝ can be related to Ĥ in the framework of the Heidelberg approach.
According to the standard paradigm of Quantum Chaos, we assume that the shape of the enclosure ensures chaotic ergodization of a single classical particle motion in the same scattering domain.At this ergodic situation universal properties of closed wave-chaotic systems can be, following the famous Bohigas-Giannoni-Schmidt (BGS) conjecture [1], effectively modelled by replacing the microscopic system's Hamiltonian (or wave operator) by random matrices Ĥ of large dimension N ≫ 1.The standard choice is to use three ensembles with Gaussian-distributed entries, GOE, GUE and GSE, composed of real symmetric, complex Hermitian and real quanternionic matrices, respectively, and labelled by the Dyson parameter β = 1, 2, 4. The choice β = 1 is used to describe time-reversal invariant systems, β = 2 coresponds to broken time reversal symmetry, and β = 4 to systems with Kramers degeneracy of energy levels.
The ensuing statistical characteristics of closed quantum-chaotic systems turn out to be universal, i.e. independent of microscopic details, when studied in the energy/frequency intervals of the length comparable to the typical distance ∆ between neighbouring energy/frequency levels.It is expected that essentially the same statistics should be observed in regularly-shaped cavities with a finite density of randomly placed scatterers inside, provided one neglects the effects of Anderson localization.The scale ∆ is assumed to be much smaller than the energy scale of the order of inverse relaxation time t e ensuring full ergodization in the chaotic enclosure.In the context of systems with disorder such time is controlled by classical diffusion and the corresponding energy scale is known as the Thouless time.Although proving BGS conjecture remains one of the great challenges in Mathematics, see e.g.[2], its validity at the level of Theoretical Physics is beyond any reasonable doubt, being supported by extensive numerics, as well as by elaborate field theory [3,4] and semiclassical computations [5,6].
Chaotic wave scattering in enclosures is an object of intensive research effort extending over several decades, with application to studies in compound nucleus scattering [7], transport properties in mesoscopic electronic systems [8] and more recently in lasing [9] as well as in manipulating light in complex media for energy deposition and imaging purposes [10].One of the central objects in both theory and experiments is the energy-dependent (or, rather, in the classical wave scattering context, frequency-dependent) unitary scattering matrix (or simply Ŝ−matrix) Ŝ(E), the elements of which describe the relationship between the vector a = (a 1 , . . ., a M ) of amplitudes of M incoming waves in all open channels to the vector b = (b 1 , . . ., b M ) of the amplitudes of outgoing waves.Since the scattering process is essentially random, the properties of Ŝ-matrix must be described using statistical language, i.e. probability distributions and correlation functions.In developing this description the Random Matrix theory (RMT) plays the central role.
The modern use of RMT for describing chaotic wave scattering statistically goes back to the seminal work of the Heidelberg group [11] who suggested to model the scattering matrix elements in the form where Ĥeff = Ĥ − i c w c ⊗ w † c , with Ĥ being N × N random matrix replacing true Hamiltonian of the closed cavity, and the energy-independent vectors of coupling amplitudes w c = (W c1 , . . ., W cN ) relate N inner states in the chosen basis to M open channels.Without restricting generality, one can take vectors w c as fixed orthogonal satisfying The orthogonality condition ensures that the ensemble-averaged scattering matrix can be assumed to be diagonal where as N → ∞ at fixed M one finds that ⟨S c (E)⟩ = 1−iγc⟨G⟩ 1+iγc⟨G⟩ and we introduced the mean value ⟨G⟩ := lim η→0 ⟨G(r, r, E + iη)⟩ = ⟨Re G⟩ − iπρ(E) for the diagonal entry of the (retarded) Green's function of the underlying closed cavity: G(r, r ′ , E + iη) := ⟨r| (E + iη − Ĥ) −1 |r ′ ⟩.This implies that ρ(E) is the mean density of states in the cavity, which at the level of RMT is given as N → ∞ by the Wigner semicircle: implying The set of parameters g c , c = 1, . . ., M provides the complete description of coupling of the medium to scattering channels in the universal regime, with the "perfect coupling" value g c = 1 (happening when sin α = 1 and γc = 1) corresponding to |⟨S c ⟩| = 0.The latter condition physically implies absence of short-time (also known as "direct") scattering processes at the channel c entrance: all the incoming flux penetrates inside the medium and participates in formation of long-living resonant structures.This situation is thus most interesting from theoretical point of view, and is frequently described by most elegant formulas.In the RMT model the perfect coupling may occur only at the center of the spectrum E = 0, where ⟨G⟩ = −i, and we restrict our calculations henceforth to that point.Let us mention also the opposite limit g c → ∞ corresponding to the channel c closed for incoming and outgoing waves.The above-described choice for the model provides the most convenient framework for studying statistics of the scattering matrix on small energy/frequency scales, comparable with separation ∆ between neighbouring resonant frequencies/energy levels of a closed system by utilizing the powerful supersymmetry approach developed earlier by Efetov [12] in the context of disordered electronic systems.Over the years it allowed to compute explicitly many statistical characteristics of the Ŝ−matrix and other closely related objects, see e.g.[13][14][15][16][17] and references therein.
Nowadays the model experimental setups to test the theoretical predictions based on Random Matrix Theory (RMT) are mainly systems of classical waves (acoustic or electromagnetic) scattered from specially built resonators, shaped in the form of the so-called chaotic billiards or/and with added scatterers inside, see e.g.[18,19].Under appropriate conditions, the associated Helmholtz equation for the electric field strength is scalar and mathematically identical to the two-dimensional Schrödinger equation of a particle elastically reflected by the contour of the microwave resonator, i.e., of a quantum billiard.Alternatively, experiments on chaotic wave scattering are performed on systems built with microwave graphs, see e.g.[20].
Whereas a lot of efforts was devoted to study of transmission and reflection of waves, which pertains to measuring the wave field outside of the scattering medium (or at its boundary with external world), an interesting question is also to understand the statistics of wave patterns inside the chaotic enclosure.This question is especially natural in view of growing interest in various aspects of coherent manipulations of wave propagation in complex media for imaging, light storage, electromagnetic compatibility tests etc, see e.g.[10,[21][22][23][24][25] and references therein.The study of statistics of radiation intensity in random medium has a long history.In particular, it has been suggested to model the wave pattern as a random superposition of running plane waves with complex coefficients [26,27]: where all wavevectors k have the same length, while the amplitudes a(k) are chosen as random gaussian complex numbers.While in closed systems with preserved time-reversal invariance one has to assume a * (k) = a(−k), the correlations between a(k) and a(−k) gradually diminish with increased degree of openness of the scattering system.The simplest prediction of such a model was a oneparameter family of possible distributions for the point intensity I = |u(r)| 2 , reducing to the simple exponential/Rayleigh distribution P(I) ∝ e −I/I for completely uncorrelated a(k) and a(−k).Despite favourably agreeing with some experimental results [28], the use of simple Gaussian model Eq.( 6) looks largely phenomenological, and certainly calls for a proper microscopic justification.Motivated by this, the present paper aims to investigate statistics of the intensity of wave field at a given point r inside a chaotic cavity relying on the same assumptions as RMT-based model Eq.( 1) for the Ŝ−matrix.We will demonstrate that the framework of the Heidelberg approach gives a possibility to derive P(I) for any fixed number M of open channels without further assumptions, at least in the simplest case of chaotic systems with fully broken time-reversal invariance.The latter is described at the level of RMT by Ĥ taken from the Gaussian Unitary Ensemble.We note that for such systems the eigenfunctions of the closed cavity are already complex, with independent, identically distributed complex and imaginary parts.Hence when applying the Gaussian wave ansatz Eq.( 6) to the open system with broken time-reversal invariance one would naturally expect that a(k) and a(−k) are uncorrelated, implying the Rayleigh law as the reference for the intensity distribution.We will indeed see how such a law emerges in the limit of very open system, with the number of scattering channels tending to infinity.However for any finite number of channels the ensuing distribution P(I) of local intensity is found to be very different and shows a powerlaw rather than exponential decay.As scattering systems with broken time-reversal invariance are now routinely realized both in "billiard"-type scattering experiments [29][30][31][32][33] and in chaotic scattering in microwave graphs, see e.g [34,35], one may expect that the predicted behaviour may be eventually tested experimentally.

II. FORMULATION OF THE PROBLEM AND THE MAIN RESULTS
We recall that the incoming waves are fed into the cavity via M channels c = 1, . . ., M , with amplitudes given by the vector a = (a 1 , . . ., a M ).This creates a field inside the cavity which we think of as a vector u in the N −dimensional inner Hilbert space.In particular, for our purpose it is convenient to think of the position basis |r⟩, associated with an appropriate coordinate system inside the cavity domain, so that the quantities u(r) ≡ ⟨r|u⟩ give precisely the amplitude of the wave in a point r inside the cavity.The corresponding intensity is given by In the framework of the Heidelberg model one can relate the vector u at a given value of the energy/frequency to the scattering matrix as ( see e.g.[14] or Eq.( 27) in [13]) where 1M stands for the identity matrix and Ĥ is the random matrix representing the inner Hamiltonian, while Ŵ is the matrix whose M columns are channel vectors w c .Further using an equivalent form of the scattering matrix given by one can bring Eq.( 7) to another well-known form, cf.e.g.Eq.( 38) in [36], conveniently written in the bra-ket notations as and implying for the intensity a representation (10) This formula is the starting point of our calculation of the probability density P(I) for the single-point intensity Relegating the technical part of the calculation, largely inspired by the methods of the works [16,17], to the body of the paper, we start with presenting and discussing our main results below.

A. Single-point intensity distribution
Given the set of coupling parameters g c ≥ 1, c = 1, . . ., M , define for a given I > 0 the parameter λ 1 > 1 as the (unique) solution of the equation The existence and uniqueness of the solution follows from the fact that the right-hand side of Eq.( 11) is positive, monotonically increasing to infinity function of λ 1 in the whole range λ 1 ∈ [1, ∞), and is equal to zero when λ 1 = 1.The intensity distribution in a single spatial point is then characterized by the probability density given explicitly in terms of λ 1 as with F c (I) given by where we introduced the notation: There are two special cases when the solution to Eq.( 11) can be explicitly written.The first one pertains to the situation when the incoming wave is incident only via a single channel, which we can choose to correspond to the channels index c = 1, whereas all other channels with 2 ≤ c ≤ M may only support outgoing waves.Indeed, setting a 1 = 1 for simplicity, and a c = 0, ∀c = 2, . . ., M we see that Eq.( 11) becomes quadratic and one immediately finds that which after some manipulations allows to show that and which further implies (17) Correspondingly, Eq.( 13) takes very explicit and rather elegant form The remaining integral, hence the probability density for the intensity I, can be evaluated in a closed form for any coupling strengths but general results are quite cumbersome.In the simplest case of a single open channel one gets whereas for the two-channel case the cumulative distribution of intensities is given by with λ 1 defined in Eq.( 15).
The second special case corresponds to all scattering channels being of equal strength: g c = g ≥ 1 for all c = 1, . . ., M .Defining the total incoming flux in all channels as and further introducing the ratio J = I/I one finds that implying that again g = −J + 1 + 2gJ + J 2 , and further finding Evaluating the integral in the closed form, we may assume M ≥ 2 as we already considered M = 1 case above.
We then find where we defined with λ 1 defined in Eq.( 22).The probability density P M (I) is then obtained by substituting Eq.( 24) into Eq.(12).The most elegant result emerges if all channels are perfectly coupled, with g = 1 implying λ 1 = 2J + 1.After some algebra we get in that case and after substituting into Eq.(12) the probability density for the intensity I takes an especially simple form: In fact for any coupling the tail behaviour can be easily extracted from Eqs. ( 12) and ( 23) and has the same powerlaw form: setting I → ∞ at a fixed value of g one finds the tail P(I) ∼ I −(M +2) .We conclude that for any finite number of channels the ensuing powerlaw-tailed distribution is quite different from the Rayleigh law predictions of the "Gaussian random wave" model.Note however that setting in Eq.( 27) the number of channels to infinity in such a way that the incoming flux per channel remains finite: lim M →∞ I/M = I < ∞ restores the Rayleigh law: lim This fact supports the view that the Gaussian wave model is asymptotically accurate if scattering system is open in an essentially semiclassic way, with many incoming channels supporting finite flux per channel.
In the figure 2 we show the results for a direct numerical simulations of the Heidelberg model against our theoretical predictions, with a very satisfactory agreement between the two.A few remarks are now in order.
Remark 1.The one-point intensity distribution presented in Eqs.( 12)-( 13) has been obtained under a physical assumption of the observation point location r to be chosen "far enough" ( much further away than the wavelength at a given energy/frequency) from the point of attached antenna/channel.Mathematically this condition has been implemented by considering the value of all scalar products ⟨r|w c ⟩ to be negligible in comparison with the norms γ c = |w c | 2 for every c = 1, . . ., M .Remark 2. In the course of derivation it has been also assumed that no irreversible losses of flux occur inside the cavity domain.In real microwave experiments this is hardly a realistic assumption, unless resonator walls made of superconducting material, like e.g. in [33].It is however well-known how to account for the uniform  10) in a given realization of GUE matrix of size N = 100.Statistics is built using 10000 different GUE matrix realizations.The solid line is the theoretical prediction given by Eq.( 19) for M = 1 with the coupling constant g 1 = 2.067, by Eq.( 20) for M = 2 and the coupling constants g 1 = 2.067, g 2 = 1.551 and Eq.( 27) for M = 5 perfectly coupled channels with g i = 1.
absorption in cavity walls in the framework of the Heidelberg approach, see e.g.[15,37].The idea is that absorption can be treated as loss of flux in the multitude of unobserved open channels, very weakly coupled to the cavity.From this angle, one can add to M observed channels a big number M ≫ 1 of channels numbered by channel indices c = M + 1, . . ., M + M , all with the same coupling strength: g M +1 = . . .= g M +M := g a and consider the limits M → ∞ and g a → ∞ while keeping M /g a = ϵ fixed.It is easy to check the result of this procedure amounts to adding to the integrand, in expressions like Eq.( 13), an extra factor e −ϵ(λ1−λ2) .The dimensionless parameter ϵ should be then interpreted as the effective rate of absorption.An alternative procedure for arriving to the same result amounts to adding a small positive imaginary part iη to the energy E in the formulation of Heidelberg model, see Eq.( 1) or its equivalent formulation Eq. (8).In particular, it is evident from Eq.( 8) that the change E → E + iη entails the loss of S-matrix unitarity, indicating that such a procedure accounts for the irreversible loss of incoming flux in the cavity due to absorption.One then finds that at the level of final formulas the net result is exactly the same factor e −ϵ(λ1−λ2) in the integrands, with the parameter ϵ given by the ratio of the imaginary part η to the mean level spacing ∆.This fact shows the equivalence of the two methods.It is easy to see that the additional exponential factor in the integrand immediately converts the most distant tails of the intensity distribution P(I) from power law to exponential ones, so that the powerlaw behaviour can be observed only in a finite interval of intensities 1 ≪ I/I ≪ ϵ −1 .
Remark 3. If the incoming flux is nonvanishing only in a single channel c = 1, the density (Eq.19) coincides with M = 1 case of the distribution of the photodissociation cross-section σ(E) in quantum chaotic systems studied in [38], see eq. ( 18) in that paper.Such coincidence is not at all accidental and is to be expected.Namely, in the Heidelberg approach the cross-section can be represented as |m⟩, with Ĥeff defined after Eq.( 1) and |m⟩ being a fixed nonrandom vector, related to a dipole moment operator.On the other hand, for M = 1 one may use the identity which, when substituted to Eq. (10) shows that the local intensity in this case is proportional to the diagonal element of the resolvent: We see that indeed intensity for M = 1 is statistically equivalent to the photodissociation cross-section σ(E), up to a constant proportionality factor.If incoming fluxes are nonvanishing in more than one channel, then As the result, the correspondence between the diagonal part of the resolvent and the local intensity is lost, hence for M > 1 the distribution in Eq.( 27) is different from the corresponding distribution for the normalized crosssection q = σ ⟨σ⟩ at perfect coupling, derived originally in [39].The latter can be given for systems of any symmetry β = 1, 2, 4 as Note however that for β = 2 the same tail behaviour is shared by Eq.( 27) and by Eq.( 32).
Remark 4. Finally, in the case when the waves are fed via a single channel c = 1 one may consider the limit g 1 ≫ g c , ∀c = 2, . . ., N describing the case of extremely weak coupling of the feeding channel.In such a limit the point intensity studied in this paper should coincide, after appropriate normalization, with the so-called "transmitted power", whose distribution for β = 2 chaotic systems has been recovered in the Heidelberg approach framework in [40] by the method of moments.One can indeed check that our Eq.( 18) in this limiting case reproduces the distribution found in [40].Following the same logic, one should expect the Eq.( 18) to be itself deducible as a limiting case from the distribution of the modulus of the off-diagonal element of the scattering matrix found in [16,17].We check in the Appendix B that this is indeed the case.Let us however stress that (i) the distribution of intensity in the general case, Eq.( 13), can not be deduced in such a way and (ii) our computation despite being inspired by [16,17], was implemented somewhat differently which helped to arrive to the final results in a rather economic way.

B. Joint probability distribution of intensities at several points
Consider now a finite number L ≪ N of the observation points at locations r 1 , . . .r L , each location being both far enough from each of the M antennae, as well as from each other.We have found that for the case of ergodic systems with broken time-reversal invariance the joint probability density of the corresponding intensities P With this relation it is then straightforward to calculate the probability density p M (I Σ ) for the sum of the intensities In particular, in the case of perfectly coupled channels the Eqs.(34) and Eq.( 27) imply together: Introducing the intensity "per point" i Σ = I Σ /L one finds that such object has the finite limiting probability density as L → ∞: (36) Remark 5. Summing up the intensities in Eq.( 10) over all N internal points in the cavity and using the completeness relation r |r⟩ ⟨r| = 1N one finds that r I r = a † Qa, where The M ×M matrix Q is one of most important objects in scattering theory known as the Wigner-Smith time delay matrix.Various aspects of its statistical properties in wave-chaotic systems enjoyed intensive research over several decades, both in the framework of RMT, see the review [41] and references therein, as well as by semiclassical methods [42][43][44].In particular, for the perfect coupling in all channels the distribution of Q is known explicitly for all β = 1, 2, 4 [45], see also [46,47] for non-perfect couplings.Combining that distribution with Eq.( 37) it is easy to verify that r I r /I, with I defined in Eq.( 21), is distributed in the same way as the diagonal entries Qcc of the matrix Q.In turn, for the perfect coupling the latter entries are known to be distributed in the same way as partial delay times τ [46] whose probability density is explicitly given by where ∆ ∼ N −1 is the mean level spacing in the closed cavity.We then see that the distribution Eq.( 36) of intensity per point i Σ /I considered in the limit of many observation points 1 ≪ L ≪ N → ∞ coincides for β = 2 with the distribution of the total scaled intensity ∆ r I r /I, i.e. sampled accross the whole cavity.This matching implies that the same result will be valid for the (properly scaled) sum of intensities over L ∼ N ϵ points, for any 0 < ϵ ≤ 1. Remark 6.For systems with preserved time reversal invariance with values β = 1, 4 the problem of finding the full distribution of the local intensity I r for arbitrary coupling constants g c , c = 1, . . ., M and its further L-point generalizations remains largely open, apart from M = 1 case and L = N case, where distributions of partial time delays are known even at the crossover between β = 1 and β = 2, see [48].However one may safely conjecture that the far tail for all these quantities should be universally given by P(I ≫ I) ∼ I − βM 2 +2 , as this behaviour in all cases is expected to be controlled by the density of narrow resonances, see discussion in p.1967 of [13].

C. Distribution of the maximal and minimal intensities in a multipoint observation
Having at our disposal the joint probability density P (L) M (I 1 , . . ., I L ) given by Eq.( 33) one can pose a natural question of the distribution of the maximal and minimal value in the observed pattern: I max = max(I 1 , . . ., I L ), I min = min(I 1 , . . ., I L ). (39) Note that extreme values of the intensity field in chaotic reverberation chambers were studied experimentally in [49].
The joint probability (33) implies that intensities in different spatial points are in general correlated, apart from the only case when the single-point intensity is given by the exponential Rayleigh law P M (I) ∝ e −I/I .Thus the posed questions belong to the domain of extreme value statistics of many correlated variables I k , which attracted a lot of attention in recent years, especially when L → ∞, see [50] for a review.One of the most studied cases in this area is one inspired by the pattern of repelling eigenvalues of large random matrices, with correlations induced via the presence of the Vandermonde factor k<l |I k − I l | in the associated joint probability density, ultimately leading to the famous Tracy-Widom distribution for the associated extreme values [50].To this end it is necessary to stress that the correlations in the pattern of intensities emerging in our problem are of a very different nature, and induced rather by the joint probability density depending on all individual intensities only via their sum k I k .Extreme value statistics for such case was not much studied, though a special case appeared in [51], which in our language would correspond to the particularly simple choice P (L) (I 1 , . . ., I L ) ∝ δ(I 1 + . . .+ I L ), and very recently also in the context of resetting problems in [52].This motivated us to perform the analysis in our case in some detail.
After some computations explained in detail in Section (III B) one finds the general relation in terms of the single-point density P M (I): whereas Prob (I min > Y ) = ∞ LY P M (I) dI.In particular, for the perfect coupling case one can use the Eq.( 27) and get and We see that in such a pattern of L observation points the typical minimal intensity scales as I typ min ∼ IL −1 and the limiting density of the variable σ min = L Imin I is given by ρ (σ min ) = (M + 1) (1 + σ min ) −(M +2) , thus of the same form as the density of the one-point intensity.
The statistics of I max is somewhat more interesting.To start with, consider the simplest case of the Rayleigh law P(I) = 1 I e −I/I obtained in the limit of many open channels, keeping the incoming flux per channel finite: lim M →∞ I/M = I < ∞, see Eq. (28).In this case it is easy to see that Setting Y /I = ln L+q we then recover in the limit L → ∞ the Gumbel distribution: smoothly interpolating between zero at q → −∞ and one for q → ∞.The Gumbel law is one of the classical Extreme Value Statistics (EVS), and is fully expected here as the intensities I l at different points are uncorrelated.Note also that the threshold of extreme values is located to the leading order sharply at I max /I = ln L + o(1).
Turning our attention now to the finite number of channels, the Eq.( 41) can be alternatively represented as In the figure 3 below we compare the probability density associated with the Eq.( 45) to the results of numerically generated intensity pattern in the Heidelberg model for the simplest case of a single-channel system.The results show a reasonable overall agreement.
Setting in Eq.( 45) Y = σ max I ln L and considering σ max > 0 fixed as L → ∞ one first notices that lim It is now straightforward to see that the typical maximal intensity in a pattern of many points is scaled logarithmically with the number L of observation points: max ∼ I ln L, and the limiting distribution for the properly rescaled maximum intensity is given by: implying that the probability density for the scaled maximal intensity σ max = Imax I ln L converges in the limit of many observation points to as long as cavity with broken time-reversal invariance is perfectly coupled to antennas.This type of extreme value statistics resembles the Frechet law density ρ (σ) = ασ −α−1 e −σ −α arising in random patterns of independent, identically distributed random variables I k , each distributed with a powerlaw density P (I) ∼ I −(α+1) .Comparing with Eq.( 27) we may identify α = M + 1, and see that had the intensities be independent, the associated Frechet density for the maximum would be different from Eq.( 47), replacing σ −1 with σ −(M +1) in the exponential factor.Eq.( 47) does not seem to appear in the literature on extreme values before.Note that it is the same law as the limiting intensity per point, eq.( 36), or partial delay times, eq.(38).An interesting open question is to verify if this property still holds for systems with preserved time-reversal invariance.
To make contact with the previously considered Rayleigh case, one may consider M → ∞ limit in Eq. (46).Recalling that in this limit we assume I = M I we further introduce σ max M = σ max and assume it to remain finite in the limit M → ∞.We also rescale t = τ /M which gives: which upon evaluating the integral by the Laplace method yields lim This agrees with the fact that the threshold of extreme values in this case is sharply at I max = I ln L, but such interchange of limits (first L → ∞, then M → ∞) misses the fine-scale Gumbel distribution, replacing it by the step function.To improve on that one has to consider the following double scaling limit in Eq.( 45): both M and L tending to infinity in such a way that ln L √ M = c, with c ∈ [0, ∞) kept constant, and also lim In such a limit one finds that lim providing a family of interpolating distributions and reducing to Gumbel for c = 0.

III. OUTLINE OF THE METHOD AND DERIVATIONS OF THE MAIN RESULTS.
A

. Distribution of one-point intensity
To characterize the distribution of one-point intensity I r we use the method of Laplace transform generating functions, and aim to calculate for p > 0 the function where u(r) := ⟨r|u⟩ is the amplitude of the wave in a point r inside the cavity, and angular brackets stand for averaging performed over the random matrix Hamiltonian Ĥ, assumed to be represented by a GUE matrix.
As the first step of the evaluation we find it to be expedient to use a variant of the Gaussian (Hubbard-Stratonovich) transformation, representing L(p) as e −pu * (r)u(r) = dq * dq π e −q * q R(q, q * ), (52) where we defined which by using Eq.( 56) below can be equivalently written as We also recall a more general Gaussian identity valid, as long as the integral over z ∈ C N is convergent, for any N × N matrix Â and any complex-valued N −component vectors a, b.
Let us recall the definition: where we defined with η > 0 being a regularization parameter, physically chosen to be proportional to the uniform absorption in the sample.Now we use Eq.( 55) for Â = −i E − Ĥ + i Γη to observe that Similarly, for Â = i E − Ĥ − i Γη we may see that Note that here and below we systematically disregard the proportionality constants, restoring them in final expressions by normalization conditions, and also find it convenient intermittently use the bra-ket notations for the scalar products, e.g.⟨z|w a ⟩ ≡ z † w a .Another useful remark is that there is a certain freedom in choosing the arrangement of the variables q, q * in front of scalar products in the exponents, and we exploited it in two different ways in Eq.(58 and Eq.( 59).This choice will be aposteriori justified by very essential simplification of the forthcoming calculations.
On the other hand, the determinant factors entering Eq.(58 and Eq.( 59) can be represented as Gaussian integrals over anticommuting N − vectors χ σ and χ * σ with σ = 1, 2: with no issues of convergence arising in this case by definition.
It is convenient to combine vectors with commuting and anticommuting components in a single 4Ndimensional supervector Φ defined as (60) and also introduce supermatrices L = diag (1, 1, −1, 1) , Λ = diag (1, 1, −1, −1).To shorten the notation we in most cases do not distinguish between the number 1 and identity matrix 1 N when the dimensions are evident from the context.
As the result, we can rewrite the function R(q, q * ) in Eq.( 54) as where the supervectors ξ σ are given by Closely following a variant of the supersymmetry approach as exposed e.g. in [13] (one may consult also the lectures [53] for the detail of similar procedures) one can perform the average over GUE matrices Ĥ and, after exploiting a supermatrix version of the Hubbard-Stratonovich transformation and peforming the Gaussian integrals over the supervectors arrive at the following representation in terms of a 4 × 4 supermatrix R (see the above references for its structure motivated by convergence arguments): where we introduced the 4N component supermatrix In what follows we assume the scaling η = ϵ/2N , with ϵ fixed as N → ∞, and in the limit N ≫ 1 perform the R− integral in Eq.( 63) by the saddle-point method, assuming the number of channels M being fixed.Repeating the same steps as in [13], the R− integral is reduced to one over a saddle-point manifold parametrized by a 4 × 4 supermatrix Q = T Λ T −1 where supermatrices T satisfy T † L T = L.In the Appendix A we give an explicit parametrization of these matrices for convenience of the reader.
To simplify the presentation we also assume for simplicity E = 0, the results for general E are obtained via the well-known rescaling using the semicircular density of GUE eigenvalues as in [13].After all these steps one arrives at the following representation: where we introduced the short-hand notation τk ≡ L−1/2 Q Λ Q k L−1/2 and used the parameters γ c as defined in Eq.( 2).
To evaluate the expression in the exponent of Eq.( 64) we use the definition Eq.(62) of supervectors ξ 1,2 to write the k th term in the sum as Due to the condition of orthogonality of channels, see Eq.( 2), we have Γk Next important assumption is to consider only the observation of intensity in points far from the channel entrances.Such a condition is taken into account assuming that ⟨w c |r⟩ = 0, ∀c = 1, . . ., M , implying that ⟨r| Γk |r⟩ = ⟨r|r⟩ δ k,0 .In principle, here one may put ⟨r|r⟩ = 1, but we leave it in this form as it will help to understand some arising structures later on.Further using the identity (67) where we introduced the supermatrices (68) Substituting Eq.(67) into Eq.(64)gives Substituting further such a R(q, q * ) into Eq.(53) one may notice that the integral over q, q * can be easily performed resulting in the following Laplace transformed density of the single-point intensity: We see that the dependence on the Laplace parameter p in Eq.( 70) is extremely simple, which is a direct consequence of the specific choice made by us in Eqs(58)-(59).This fact allows us to invert the Laplace transform immediately, getting the probability density for the singlepoint intensity via where we eventually replaced ⟨r|r⟩ = 1.Explicit evaluation of such an integral is sketched in the Appendix A, and leads to the form featuring in Eq.( 12).

B. Joint probability of intensities at L observation points and extreme value statistics
Let us now consider the computation of the joint probability density P (L) M (I 1 , . . ., I L ) of wave intensities , where u(r l ) := ⟨r l |u⟩ is the amplitude of the wave in a point r l inside the cavity, l = 1, . . ., L. To start with, we define for the parameters p 1 > 0, . . ., p L > 0 the joint Laplace transform: which after applying Gaussian (Hubbard-Stratonovich) transformations L times takes the form q * l q l R({q l , q * l }), (73) where R(q 1 , q * 1 , . . ., q L , q * L ) := e −i L l=1 √ p l (q * l u(r l )+q l u * (r l )) .
(74) Now one may notice that Eq.( 9) implies where we defined and similarly Using the above one can see that we need to evaluate R(q 1 , q * 1 , . . ., q L , q * Now, comparing Eq.(78)with Eq.( 54) one may notice that putting p = q = 1 in the later, and replacing also |r⟩ → |X⟩ makes the two expressions identical.Moreover, assuming that all observation points to be located far from every channel entrance implies that the vector |X⟩, being a linear combination of |r l ⟩, will be orthogonal to all channel vectors |w c ⟩. Therefore the evaluation of ensemble average in Eq.( 78) should be simply read off from the expression Eqs.(69) for R(q 1 , q * 1 ) implying that R(q 1 , q * 1 , . . ., q L , q * L ) := R (⟨X|X⟩) , where we made explicit the fact that R depends on the variables q l , q * l for all l = 1, . . ., L (as well as on the Laplace parameters p l ) only via the norm: where in the above we exploited the inner basis orthogonality: ⟨r l1 |r l2 ⟩ = δ l1l2 .Substituting such dependence back into Eq.(73) , passing to polar coordinates: q l = √ R l e iθ l , and finally rescaling R l → p −1 l R l leads to: In such a form the joint Laplace transform can be easily inverted due to the well-known identity involving the Bessel function J 0 (z): yielding the joint probability density of L intensities in the form: (82) At the next step we use the following chain of identities: Substituting this back to Eq.( 82), changing the order of integrations and using that one arrives to the following representation for the joint probability density: where for the function Φ L (I; t) one easily finds that Here in the last step we used the inversion of Eq.( 84).
To reflect that this joint probability density depends on individual intensities only via their sum We will concentrate on the former as the most interesting.Using the same type representation as in Eq.(83): 90) one easily finds: where we defined Expanding the binomial and using the identity: where θ(t) = 1 for t > 0 and zero otherwise, one finds In particular, one can see that and T L (t; Y ) = 0 for t > LY .This fact, together with the relation Eq.( 33) allows to integrate by parts in Eq.( 91), which eventually results in the first of relations Eq.( 40).

IV. CONCLUSION
With this work we obtained a pretty complete description of intensity statistics inside irregularly shaped microwave resonator in the quantum chaos regime with broken time-reversal invariance, including multipoint distributions and extreme value statistics.In case of finite number of open channels and no absorption inside all expressions can be, in principle, reduced to elementary functions.In such a case the one-point intensity is generically powerlaw-distributed, in clear difference with the well-known random Gaussian wave conjecture, cf.Eq.( 6), predicting the exponential Rayleigh law.The latter is only recovered in the very open system limit, while keeping the incoming flux per channel constant.If however uniform losses in the cavity (modelled e.g. by infinite number of weakly coupled hidden channels) are taken into account, the power law remains only valid in a restricted range of intensities, being cut exponentially at larger values.Interestingly, we demonstrated that the joint probability density of intensities sampled at many points depends only on the sum of individual intensities.We do not have a transparent explanation of such a pattern, though it may be traced back to the statistical independence of the real and imaginary parts of the complex wavefunctions in a closed cavity, so is definitely expected to hold only for systems with fully broken time reversal invariance.Even with such a simple dependence, the intensities at different spatial points are clearly correlated, unless the system is in the Rayleigh regime.In particular, by extracting the statistics of the highest intensity in an observation pattern of L points explicitly in the perfect coupling regime we were able to demonstrate that the ensuing extreme values distribution for fixed M and L → ∞, Eq.( 47), differs from the classical extreme value statistics.This provides an example of nontrivial EVS which is potentially accessible in experiments, provided the losses due to absorption can be effectively controlled.The problem of characterizing multipoint and extreme value statistics in systems with preserved time reversal invariance remains currently open.Calculations in the supersymmetry approach for that case are expected to yield much more cumbersome structures, cf.[16,17] for the statistics of the modulus of off-diagonal entries of the S−matrix.In particular, we expect that the property of the joint distribution of intensities depending only on the sum of individual intensities will be lost in the systems with preserved time-reversal invariance.Such a study presents therefore an interesting challenge for the future research.
Another possible extension is to consider modifications of intensity statistics by the effects of Anderson localization, which are operative in disordered systems of finite spatial dimension.To this end it is worth mentioning that some aspects of wave intensity statistics inside quasi-1D disordered samples have been recently under experimental investigation, see e.g.[54,55].In the framework of the supersymmetry formalism exploited in the present work this would require to go beyond the effectively zero-dimensional limit and combine the 1D nonlinear σ− model description of interacting diffusive modes, see [56] and references therein, with the Heidelberg model formalism.For a few examples of recent studies of not dissimilar problems see e.g.[55] and [57].
It is immediate to check that in this parameterization Str Λ Q and Sdet 1 4 + γ c Λ Q take the form where F M (I) will be defined in (100) below.Here we note that as explained in the Appendix of the paper [40] the so-called "Efetov-Wegner" term δ (I) in Eq.( 99) gets eventually cancelled and can be omitted.The function F M (I) is given explicitly by After further manipulations using 2 and noticing that the δ−functional constraint implies and U(r) is given by with λ 1 for a given r being defined via λ 1 = (g 1 + g M )r 2 + (g 1 − g M ) 2 r 4 + 4r 2 (g 1 g M − 1) + 4 2(1 − r 2 ) .
(104) On the other hand recall that according to Eq.( 1) we have where Ŵ Ŵ † = M c=1 w c ⊗ w † c .Consider now the limit γ M = |w M | 2 → 0 while keeping |w c | 2 = γ c of the order of unity for all c ̸ = M .Physically this corresponds to almost closing the channel with c = M , with the effective coupling g M ≈ 1 2|w M | 2 ≫ g c , ∀c < M .It is then easy to see that in such a limit |S 1M | 2 → 0, whereas the product |S 1M | 2 g M /2 remains finite and simply proportional to the intensity I at a single point inside the cavity given by Eq.( 10), provided we reduce the number of open channels by one and consider the incoming wave amplitudes a c to be nonzero only for c = 1.We therefore can extract the probability distribution of the point intensity in such a case by performing the limit g M → ∞ and r 2 → 0 in Eq.(102) while keeping r 2 g M = 2I and g c , ∀c < M finite.In such a limiting procedure we get: Further we have (107) where we used the definitions Eq.( 14) and Eq.( 16):

1 √
k a l e −ikx l + b l e ikx l .

Figure 2 .
Figure 2. The histogram shows numerically evaluated probability density P (I) of the single-point intensity I in the Heidelberg model for an M-channel system, with the incoming flux present only in the first channel.The intensity data are generated according to Eq.(10) in a given realization of GUE matrix of size N = 100.Statistics is built using 10000 different GUE matrix realizations.The solid line is the theoretical prediction given by Eq.(19) for M = 1 with the coupling constant g 1 = 2.067, by Eq.(20) for M = 2 and the coupling constants g 1 = 2.067, g 2 = 1.551 and Eq.(27) for M = 5 perfectly coupled channels with g i = 1.

M
(I 1 , . . ., I L ) is very simply related to the previously studied one-point density P (1) M (I) := P M (I) via:

Figure 3 .
Figure 3.The histogram shows numerically evaluated probability density of the maximum intensity in a pattern of L = 2, L = 6 and L = 25 internal points in the Heidelberg model for a single perfectly open channel system.The intensity data are first generated according to Eq.(10) at 200 randomly chosen inner points in a given realization of GUE matrix of size N = 200.Statistics is built using 100 random subsets of L points per every realization, and the displayed data correspond to 1000 different GUE matrix realizations.The solid line for a fixed L is the theoretical prediction given by Eq.(45).
k .In particular, for finding the probability density for the sum of all intensities, Eq.(34), we use the identity: next step is to consider the simplest extreme value statistics, the distributions of the maximal and the minimal value in the pattern, defined as Prob (I max < Y ) =