text:
# Statistics of transmitted power in multichannel dissipative ergodic structures

## I Introduction

Transport through open chaotic systems is often viewed as a scattering process. Standard examples of systems of that kind are compound-nuclei, mesoscopic quantum dots and wires, microwave cavities and acoustic or ultrasonic bodies \[1-12\]. Incident waves are introduced into a disordered or irregularly shaped part of the structure via channels, e.g. waveguides or infinitely long ideal leads. Assuming negligible dissipation, transport properties are obtained by relating incident and outgoing wave fields in terms of the unitary scattering matrix $`S`$. If internal dissipation is not negligible, it can be simulated by the action of additional open channels key-9 . This is the case, for example, with microwave cavities where non-perfectly reflecting walls cause loss of wave energy or with ultrasonic solids where internal friction acts in the bulk key-8 ; key-9 ; key-12 .

The scattering model thus applies to open systems both with and without internal dissipation. The scattering approach provides a useful tool for statistical characterization of chaotic transport. Assuming the wave dynamics inside the system to be ergodic so that the entire phase space of the system is explored \[2-6\], the scattering approach is combined with a statistical analysis based on Random Matrix Theory. For systems without losses, one can make assumptions on the statistics of the $`S`$-matrix key-2 ; key-11 . An alternative method uses a Random Matrix assumption on the level of the wave equation associated with the closed non-dissipative structure. Here the basic object is the Green function (resolvent) related to that wave equation, and the method works equally well in both open and closed, dissipative and non-dissipative complex structures \[1-9,10,12\]. Matrices from a Random Matrix Ensemble then replace the wave equation’s linear differential operator, and the problem of constructing various moments of the transport characteristics is expressed in terms of ensemble averages of the products of the resolvents. Transport characteristics calculated in that way are known, under certain conditions, to describe results of experimental measurements in systems with ergodic wave motion \[1-6\].

In this paper we are interested in characterizing the wave power $`T`$ transmitted between a source at site $`i`$ and a receiver at site $`j`$ in a closed system with internal losses. The statistics of $`T`$ are potentially useful for studies of power transmission in complex reverberant structures \[12-14\], where both mean power and the magnitude of its fluctuations away from the mean are important. In particular, we wish to calculate average $`T`$ and $`T^2`$ with the ultimate goal to compare with measurements such as these of Ref. . The complex amplitude of the transmitted wave is simply proportional to the off-diagonal matrix element of the resolvent: $`G(E)`$ $`\left[EI+i\epsilon IH+i\mathrm{\Gamma }\right]^1`$. Here $`H`$ corresponds to the Hamiltonian of the closed non-dissipative chaotic structure. The matrix $`\mathrm{\Gamma }`$ describes coupling to external channels or internal local-in-space losses, $`I`$ is the identity matrix, the parameter $`\epsilon >0`$ describes uniform dissipation and $`E`$ is the spectral variable. The quantity of prime interest is $`T=\left|G_{ij}\left(E\right)\right|^2,ij`$, i.e. the product of retarded and advanced Greens functions (propagators): $`G_{ji}^R(E)`$ $`\left[EI+i\epsilon IH+i\mathrm{\Gamma }\right]_{ji}^1`$ and $`G_{ij}^A=\left(G_{ji}^R(E)\right)^{}`$ respectively. Except for slowly varying factors of receiver gain and source strength, the quantity $`T`$ represents the ultrasonic power of Ref. ; see also .

The fluctuations in $`T`$, as measured in Ref.  were in only modest agreement with theoretical predictions based on a simplified version of the random matrix approach. The moments of $`T`$ were calculated there using a naive form of ensemble averaging, \[8,12-14\]. This relatively simple approach uses statistical assumptions for eigenfunctions and the real parts of the eigenvalues of the open (dissipative) system identical to those of the corresponding closed system. As will be seen later, such an assumption is strictly justified only for a special case of uniform dissipation. In a more general situation this approach fails. A proper treatment calls for a more elaborate technique, which we outline and present below.

When losses are negligible the systems discussed in Ref. \[12-14\] are invariant under time reversal. The appropriate choice for the corresponding random matrix $`H`$ should therefore be a real symmetric matrix taken e.g. from the Gaussian Orthogonal Ensemble (GOE). In principle, the powerful methods of ensemble averaging we employ here can be used for such an ensemble, but the calculations are technically involved and will be presented in a separate publication.

Here we address ourselves to the somewhat simpler case in which $`H`$ is complex Hermitian, generic for systems with broken time reversal invariance. Correspondingly, $`H`$ is treated as a Hermitian $`N\times N`$ matrix consisting of uncorrelated centered random complex numbers, their variances defined by: $`H_{ij}H_{lm}^{}=\left(\lambda ^2/N\right)\delta _{il}\delta _{jm}`$ with angular brackets indicating ensemble averaging. Such an $`H`$ is a member of the Gaussian Unitary Ensemble (GUE) of random matrices . Although our present results on the statistics of $`T`$ are not directly applicable to the time-reversal invariant systems discussed in \[12-14\], they may elucidate the discrepancies found in Ref.  between measurements and the predictions of the ’naive’ averaging. They also develop and illustrate the mathematical methods which will be used for a proper non-perturbative analysis of the time reversal-invariant problem. The present calculations are also relevant for scattering systems with broken time reversal invariance as exemplified in certain chaotic billiards , optical and semiconductor superlattices  and quantum graphs . In fact, our results on the distribution of the off-diagonal elements of the resolvent extend earlier studies concentrated on diagonal entries for the same quantity, see \[22-23\] and references therein. Let us finally mention that there exist clear analogy between our research and that presented in the paper , see also the review . However, the model considered in \[19-20\] did not take local dampers into account, but rather addressed effects of Anderson localization.

The damping matrix $`\mathrm{\Gamma }`$ is in general Hermitian positive semi-definite. In our model there is no loss of generality in assuming it to be diagonal. Indeed, in view of the rotational invariance of the Gaussian Unitary ensemble: $`HUHU^1`$ ($`U^1=U^{}`$) we always can select the basis which diagonalizes $`\mathrm{\Gamma }`$, bringing it to the form $`\mathrm{\Gamma }=diag\{\gamma ,\gamma ,\mathrm{}\gamma ,0,\mathrm{}0\}`$. The number $`M<N`$ of non-zero entries can be interpreted either as a number of equivalent open channels in the scattering system  or a number of equivalent localized ’dampers’ in a closed system with losses. While we take all the $`\gamma `$s to be equal, the expressions we develop are easily generalized to the case of varying damper strengths. It should be stressed that in general the matrices $`\mathrm{\Gamma }`$ and $`H`$ do not commute, and therefore the eigenvectors and eigenvalues of the ’effective non-Hermitian Hamiltonian’ $`Hi\mathrm{\Gamma }`$ are not trivially related to those of $`H`$. This very fact makes the naive averaging incorrect. In contrast, the term $`i\epsilon I`$ interpreted as the ’uniform damping’ preserves eigenvectors of $`H`$ and just adds a uniform shift $`i\epsilon `$ to all eigenvalues.

The presence of an $`N\times N`$ random Hamiltonian $`H`$ in the expression for the resolvent matrix $`G`$ enables us to carry out the ensemble averaging exactly using the supersymmetry method . Application of this non-perturbative technique leads to an expression for the entire probability distribution function of $`T`$. We present the corresponding derivation in Section II. In Section III, we compare the results for the first two moments of $`T`$ as obtained using the supersymmetry method and the methods of Ref. \[12-14\]. The results are then verified numerically by direct simulation of the model. Section IV contains conclusions.

## II Probability distribution function of the power. Supersymmetric calculation

In the previous Section we defined the power $`T`$ as a product of advanced and retarded Green functions $`G_{ij}`$. Our goal is to compute the statistics, i.e. ensemble averages: $`T_H`$, $`T^2_H`$, etc., where subscript $`H`$ designates averaging with the Gaussian weight $`\mathrm{exp}\left\{\frac{N}{2\lambda ^2}TrH^{}H\right\}`$. At the first stage of the supersymmetric calculation we make use of the following identities for the inverse propagator $`D_{ij}=\left[E+i\epsilon H+i\mathrm{\Gamma }\right]_{ij}`$ :

$$det𝔇_b^1=\left[dS^{}\right]\left[dS\right]\mathrm{exp}\left\{i𝔏_b(E,S)\right\},$$

$$det𝔇_f=\left(1\right)^N\left[d\chi ^{}\right]\left[d\chi \right]\mathrm{exp}\left\{i𝔏_f(E,\chi )\right\}.$$

Here we introduced $`2N`$dimensional vectors $`S^T=(S_1^T,S_2^T)`$ and $`\chi ^T=(\chi _1^T,\chi _2^T)`$, consisting of complex commuting or bosonic (b) variables and anticommuting or fermionic (f) variables respectively. $`𝔇_b=diag\{D,D^{}\}`$ and $`𝔇_f=diag\{D,D^{}\}`$ are $`2N\times 2N`$ block diagonal matrices, and $`𝔏_b(E,S)=S^{}𝔇_bS`$, $`𝔏_f(E,\chi )=\chi ^{}𝔇_f\chi `$. The negative sign of $`𝔇_b^{22}`$ is necessary for convergence of the integrals in what follows. Differentiating the first equality with respect to $`𝔇_{bji}^{11}`$ and $`𝔇_{bij}^{22}`$, and then combining the result with the second equality, we obtain:

$`T=`$ $`𝔇_{ij}^1𝔇_{ij}^1`$
$`=`$ $`{\displaystyle \left[d\mathrm{\Phi }^{}\right]\left[d\mathrm{\Phi }\right]S_{1j}^{}S_{1i}S_{2i}^{}S_{2j}\mathrm{exp}\left\{i𝔏(E,\mathrm{\Phi })\right\}},`$ (1)
where the integration involves four-component supervectors $`\mathrm{\Phi }^T=(S^T,\chi ^T)`$, and where $`𝔏(E,S)=𝔏_b(E,S)+𝔏_f(E,\chi )=\mathrm{\Phi }^{}𝔇\mathrm{\Phi }`$, $`𝔇=diag\{𝔇_b,𝔇_f\}`$ . Because the random matrix $`H`$ is in the exponent, $`T`$ is now suitable for ensemble averaging.

In a similar fashion, employing the Wick theorem one can verify the following formula necessary for the calculation of an arbitrary moment of transmitted power $`T^n_H`$:

$`T^n`$ $`=\left(n!\right)^2{\displaystyle \left[d\mathrm{\Phi }^{}\right]\left[d\mathrm{\Phi }\right]S_{1j}^nS_{1i}^nS_{2i}^nS_{2j}^n}`$
$`\times `$ $`\mathrm{exp}\left\{i𝔏(E,\mathrm{\Phi })\right\}=S_{1j}^nS_{1i}^nS_{2i}^nS_{2j}^n_\mathrm{\Phi }.`$ (2)
A shorthand notation $`\mathrm{}_\mathrm{\Phi }`$ has been introduced for the ’Gaussian’ integration over the supervector components. Hereafter we use the more convenient ’’ (’retarded-advanced’)-block notation for supervectors and supermatrices, see for example Ref. . With the supervector $`\mathrm{\Psi }^T=(S_1^T,\chi _1^T,S_2^T,\chi _2^T)`$ and the $`4\times 4`$ supermatrices $`L=diag\{1,1,1,1\}`$, $`\mathrm{\Lambda }=diag\{1,1,1,1\}`$ the exponent in the integrand reads:

$$𝔏(E,\mathrm{\Psi })=E\mathrm{\Psi }^{}\left(IL\right)\mathrm{\Psi }+i\mathrm{\Psi }^{}\left(\mathrm{\Gamma }\mathrm{\Lambda }L\right)\mathrm{\Psi }$$

$$\mathrm{\Psi }^{}\left(HL\right)\mathrm{\Psi }+i\epsilon \mathrm{\Psi }^{}\left(IL\right)\mathrm{\Psi },$$

and, Eq. (II) becomes:

$$T^n=\left(n!\right)^2F^n\left[\mathrm{\Psi }\right]_\mathrm{\Psi },F\left[\mathrm{\Psi }\right]=S_{1j}^{}S_{1i}S_{2i}^{}S_{2j}.$$
(3)
The ensemble averaging can now be performed with the aid of the identities :

$$\mathrm{exp}\left\{i\mathrm{\Psi }^{}\left(HL\right)\mathrm{\Psi }\right\}_H=\mathrm{exp}\left\{\frac{1}{2N}StrA^2\right\},$$

$$A_{pq}^{\left(kl\right)}=L_{kk,pp}^{1/2}\underset{i=1}{\overset{N}{}}\left(\mathrm{\Psi }_i\right)_k^p\left(\mathrm{\Psi }_i^{}\right)_l^qL_{ll,qq}^{1/2},$$

$$StrAL=\mathrm{\Psi }^{}\left(IL^{1/2}\mathrm{\Lambda }L^{1/2}\right)\mathrm{\Psi }=\mathrm{\Psi }^{}\left(I\mathrm{\Lambda }L\right)\mathrm{\Psi },$$

where we have set $`\lambda =1`$, and introduced a $`4\times 4`$ supermatrix $`A`$. Thus

$$T^n_H=\left(n!\right)^2\left[d\mathrm{\Psi }^{}\right]\left[d\mathrm{\Psi }\right]F^n\left[\mathrm{\Psi }\right]\mathrm{exp}\left\{iE\mathrm{\Psi }^{}\left(IL\right)\mathrm{\Psi }\mathrm{\Psi }^{}\left(\mathrm{\Gamma }\mathrm{\Lambda }L\right)\mathrm{\Psi }\frac{1}{2N}StrA^2\epsilon StrAL\right\},$$
(4)
where $`k`$, $`l`$ distinguish between retarded and advanced, ($`1`$ and $`2`$) supermatrix blocks indices and $`p`$, $`q`$ equal to $`b`$ or $`f`$ . The next stages of supersymmetric procedure include : 1) the Hubbard-Stratonovich transformation, that removes quartic (in $`\mathrm{\Psi }`$) term in the exponential; 2) $`\mathrm{\Psi }`$variables integration; 3) evaluation of the remaining integral using saddle point approximation in the limit $`N\mathrm{}`$. We have, after step 1):

$`T^n_H=\left(n!\right)^2`$ $`{\displaystyle \left[dR\right]\mathrm{exp}\left\{\frac{N}{2}StrR^2+i\epsilon NStrR\mathrm{\Lambda }+iStrRA\right\}}`$ (5)
$`\times `$ $`{\displaystyle \left[d\mathrm{\Psi }^{}\right]\left[d\mathrm{\Psi }\right]F^n\left[\mathrm{\Psi }\right]\mathrm{exp}\left\{i\left(E\mathrm{\Psi }^{}L\mathrm{\Psi }+i\mathrm{\Psi }^{}\left(\mathrm{\Gamma }\mathrm{\Lambda }L\right)\mathrm{\Psi }\right)\right\}}.`$
Since $`StrRA=\mathrm{\Psi }^{}L^{1/2}RL^{1/2}\mathrm{\Psi }`$, for an arbitrary $`4\times 4`$ supermatrix $`R`$,

$`T^n_H=\left(n!\right)^2`$ $`{\displaystyle \left[dR\right]\mathrm{exp}\left\{\frac{N}{2}StrR^2+i\epsilon NStrR\mathrm{\Lambda }\right\}\left[d\mathrm{\Psi }^{}\right]\left[d\mathrm{\Psi }\right]F^n\left[\mathrm{\Psi }\right]}`$
$`\times `$ $`\mathrm{exp}\left\{i\mathrm{\Psi }^{}L^{1/2}\left(EI_4I_NRI_Ni\mathrm{\Lambda }\mathrm{\Gamma }\right)L^{1/2}\mathrm{\Psi }\right\}.`$ (6)
Using the Gaussian nature of the integral:

$$\left[d\mathrm{\Psi }^{}\right]\left[d\mathrm{\Psi }\right]\mathrm{exp}\left\{i\mathrm{\Psi }^{}f\mathrm{\Psi }\right\}=Sdetf^1,$$

the following general relation, can be derived similarly to Eq. (3):

$`{\displaystyle \left[d\mathrm{\Psi }^{}\right]\left[d\mathrm{\Psi }\right]\left[\left(\mathrm{\Psi }_i\right)^b\left(\mathrm{\Psi }_j^{}\right)^b\left(\mathrm{\Psi }_i\right)^b\left(\mathrm{\Psi }_j^{}\right)^b\right]^n}`$ (7)
$`\times `$ $`\mathrm{exp}\left\{i\mathrm{\Psi }^{}f\mathrm{\Psi }\right\}=\left(n!\right)^2f_{12,bb_{ij}}^1f_{21,bb_{ji}}^1Sdetf^1`$
Setting $`f=L^{1/2}\left(EI_4IRIi\mathrm{\Lambda }\mathrm{\Gamma }\right)L^{1/2}`$, we integrate out the components of $`\mathrm{\Psi }`$ with the help of Eq. (7):

$$T^n_H=\left[dR\right]F^n\left[𝔊\right]\mathrm{exp}\left\{N\left[R\right]+\delta \right\},$$
(8)
where

$`F\left[𝔊\right]`$ $`=𝔊_{12,bb}^1𝔊_{21,bb}^1,𝔊=EI_4R,`$
$`\left[R\right]`$ $`={\displaystyle \frac{1}{2}}StrR^2+Str\mathrm{ln}\left(EI_4R\right),`$
$`\delta =`$ $`i\epsilon NStrR\mathrm{\Lambda }`$
$`MStr\mathrm{ln}\left[I_4i\gamma \mathrm{\Lambda }\left(EI_4R\right)^1\right].`$
See Appendix A for the details.

$`T^n_H`$ is now written as an integral over $`4\times 4`$ supermatrix $`R`$. The stationarity condition for $`\left[R\right]`$, in the limit of large $`N`$, yields a stationary point $`R_s`$, satisfying: $`R_s=1/\left(EI_4R_s\right)`$. The solution is not unique, it is a saddle manifold in a space of $`4\times 4`$ supermatrices, spanned by ’pseudounitary’ supermatrices $`𝔗`$:$`R_s=EI_4/2+i\pi \nu 𝔗^1\mathrm{\Lambda }𝔗=EI_4/2\pi \nu Q`$, where $`\nu =\sqrt{4E^2}/\left(2\pi \right)`$ is Wigner’s semicircular mean density of eigenvalues (GUE, $`\lambda =1)`$. See Refs.  for the explicit form of supermatrix $`Q`$.

After integrating out local fluctuations over directions $`R`$ orthogonal to the manifold of stationary points (the procedure is asymptotically exact for large $`N`$) the remaining integration goes over the manifold parametrized by $`Q`$:

$`T^n_H`$ $`=\left(\pi \nu \right)^{2n}{\displaystyle \left[dQ\right]\left(Q_{12,bb}Q_{21,bb}\right)^n}`$
$`\times S\stackrel{M}{det}\left[I_4+i{\displaystyle \frac{E}{2}}\gamma \mathrm{\Lambda }+i\pi \nu \gamma Q\mathrm{\Lambda }\right]`$
$`\times \mathrm{exp}\left\{i\epsilon \pi \nu NStrQ\mathrm{\Lambda }\right\}.`$ (9)
The expression for the $`n`$-th moment of power allows one to find the entire distribution function $`P\left(T\right)`$, cf. :

$$𝒫\left(T\right)=\left[dQ\right]\delta \left(T\left(\pi \nu \right)^2Q_{12,bb}Q_{21,bb}\right)Y\left(Q\right),$$

or, for the ’scaled’ power $`y=T/\left(\pi \nu \right)^2`$:

$$𝒫\left(y\right)=\left[dQ\right]\delta \left(yQ_{12,bb}Q_{21,bb}\right)Y\left(Q\right),$$
(10)
where

$`Y\left(Q\right)`$ $`=S\stackrel{M}{det}\left[I_4+i{\displaystyle \frac{E}{2}}\gamma \mathrm{\Lambda }+i\pi \nu \gamma Q\mathrm{\Lambda }\right]`$
$`\times `$ $`\mathrm{exp}\left\{i\epsilon \pi \nu NStrQ\mathrm{\Lambda }\right\}.`$
Evaluation of the superintegral in Eq. (10) is presented in Appendix A. The result reduces to:

$$𝒫\left(y\right)=\delta \left(y\right)+\left(\frac{d}{dy}+y\frac{d^2}{dy^2}\right)_1^{\mathrm{}}𝑑\lambda _1_1^1𝑑\lambda _2\delta \left(y+1\lambda _1^2\right)\frac{\lambda _1^2\lambda _2^2}{\left(\lambda _1\lambda _2\right)^2}\mathrm{exp}\left\{ϵ\left(\lambda _1\lambda _2\right)\right\}\left(\frac{g+\lambda _2}{g+\lambda _1}\right)^M$$
(11)
$$=\left(\frac{d}{dy}+y\frac{d^2}{dy^2}\right)\frac{\mathrm{exp}\left\{ϵ\sqrt{1+y}\right\}}{2\sqrt{1+y}\left(g+\sqrt{1+y}\right)^M}_1^1𝑑\lambda _2\frac{\sqrt{1+y}+\lambda _2}{\sqrt{1+y}\lambda _2}\mathrm{exp}\left\{ϵ\lambda _2\right\}\left(g+\lambda _2\right)^M,$$

where $`g=\left(1/\gamma +\gamma \right)/\left(2\pi \nu \right)`$and $`ϵ=2\pi \nu N\epsilon `$.

Setting the ’uniform damping’ $`ϵ`$ to zero we were able to evaluate the remaining integral explicitly (see Appendix A for the details) and Eq. (11) yielded:

$`𝒫\left(y\right)`$ $`={\displaystyle \frac{\left\{\left(g1\right)^{M+1}\left(g+1\right)^{M+1}\right\}p_1+\left\{\left(g+1\right)^{M+1}+\left(g1\right)^{M+1}\right\}p_2}{8\sqrt{\left(1+y\right)^5}\left(g+\sqrt{1+y}\right)^{M+2}}},`$
$`p_1`$ $`={\displaystyle \frac{g\left(y2\right)}{\left(M+1\right)}}\left(g+\left(M+2\right)\sqrt{1+y}\right)\left(y+1\right)\left(2+\left(M+1\right)\left(y+2\right)\right),`$ (12)
$`p_2`$ $`=2\left(y+1\right)\left(g+\left(M+2\right)\sqrt{1+y}\right).`$
Eqs. (11) and (12) constitute the main result of this Section.

At this point it is interesting to observe that Eq. (12) can be in fact used to cover the case of uniform damping in closed system: $`ϵ>0,M=0`$. For this we note that in the limit $`\gamma 0`$ (i.e. $`g1/2\pi \nu \gamma \mathrm{}`$) and $`M\mathrm{}`$ but $`\gamma Mconst`$ the factor $`\left(g+\lambda _2\right)^M/\left(g+\lambda _1\right)^M`$ in the integrand of Eq. (11) is converted to $`\mathrm{exp}\left\{2\pi \nu \gamma M\left(\lambda _2\lambda _1\right)\right\}`$. Such a replacement is equivalent to generating an effective uniform damping $`ϵ=2\pi \nu \gamma M`$. The fact that the large number of weakly open channels (or weak local-in-space dampers) is essentially equivalent to uniform damping is well known, see e.g. . Performing the limit $`g\mathrm{},M\mathrm{}`$ when keeping the product $`2\pi \nu \gamma M=ϵ`$ finite we find that the distribution Eq. (12) is reduced to:

$`𝒫\left(y\right)`$ $`={\displaystyle \frac{\mathrm{exp}\left\{ϵ\sqrt{1+y}\right\}\mathrm{sinh}ϵ}{4ϵ\sqrt{\left(1+y\right)^5}}}`$
$`\times \left[ϵ^2\left(y+1\right)\left(y+2\right)\left(y2\right)\left(1+ϵ\sqrt{1+y}\right)\right]`$
$`+{\displaystyle \frac{\mathrm{exp}\left\{ϵ\sqrt{1+y}\right\}\mathrm{cosh}ϵ\left(1+ϵ\sqrt{1+y}\right)}{2\sqrt{\left(1+y\right)^3}}},`$ (13)
This distribution of transmitted power for systems with uniform damping is interesting and important on its own.

Let us consider now a few other regimes. For the weakly damped system ($`M`$ is fixed, $`g1`$), Eq. (12) can be approximated by:

$$𝒫\left(y\right)\frac{4+y}{4\sqrt{\left(1+y\right)^5}}+O\left(\frac{1}{g^2}\right).$$
(14)
The asymptotic behavior of $`𝒫\left(y\right)`$ in the limit $`y\mathrm{}`$ for any $`M`$ and $`g`$ is given by:

$`𝒫\left(y\right)`$ $``$ $`{\displaystyle \frac{\left(M+1\right)}{4y^{\left(M+3\right)/2}}}\left\{\left(g+1\right)^{M+1}\left(g1\right)^{M+1}\right\}`$ (15)
$`+O\left({\displaystyle \frac{1}{y^{\left(M+4\right)/2}}}\right),`$
which shows that moments $`y^n`$ exist only for $`n<\left(M+1\right)/2`$. At the same time, as it follows from Eq. (11), any non-zero $`ϵ`$ guarantees existence of all moments. For large $`y`$, the asymptotic forms of the probability distribution function at non-zero $`ϵ`$ are:

$$𝒫\left(y\right)\frac{ϵ\mathrm{sinh}ϵ}{\sqrt{y}}\mathrm{exp}\left\{ϵy^{1/2}\right\},$$
(16)
$$𝒫\left(y\right)\frac{g^M\mathrm{exp}\left\{ϵy^{1/2}\right\}}{y^{\frac{M+1}{2}}},$$
(17)
for $`M=0`$, and for finite $`M`$ respectively.

Finally we want to compare the results of this Section with numerical solution of the model RMT problem. For this goal we numerically generate an ensemble of $`N\times N`$ Hermitian random matrices $`[H]`$, typically choosing $`1500`$ samples from the ensemble and $`N=1000`$. The entries of the matrix $`H`$ are constructed using a random number generator, with $`H_{ij}H_{lm}^{}=\left(1/N\right)\delta _{il}\delta _{jm}`$. To simulate the uniform damping and the case of a finite number of local dampers we take $`\mathrm{\Gamma }=\epsilon I`$ and $`\mathrm{\Gamma }=diag\{\gamma ,\gamma ,\mathrm{},\gamma ,0,\mathrm{}0\}`$ (with $`MN`$ non-zero entries) respectively. Then, for all members of our ensemble we generate the off-diagonal elements of the resolvent matrix $`G_{ij}(E)=\left[EI+i\mathrm{\Gamma }H\right]^1`$ modeling the response at site $`i`$ due to excitation at site $`j`$, with $`E`$ being the spectral parameter.

We first consider the case of uniform damping: $`\mathrm{\Gamma }=\epsilon I`$. The modal density $`\nu `$ for such a system is approximated by Wigner’s semicircular law: $`\nu =\sqrt{4E^2}/\left(2\pi \right)`$. Therefore, for a fixed size of the matrices $`N`$ and spectral variable $`E`$, we can explore the range of $`ϵ`$ by changing $`\epsilon `$. For $`E=0`$ the modal density is $`\nu =1/\pi `$ and so we need not distinguish between $`T`$ and $`y`$. In Fig.1, the numerically obtained histograms are compared with $`𝒫\left(y\right)`$ (Eq. (11)) for several values of $`ϵ`$. We see that numerical results correspond well with the theoretical curves.

This procedure was repeated for the damping matrix $`[\mathrm{\Gamma }]=diag\{\gamma ,\gamma ,\mathrm{},\gamma ,0,\mathrm{}0\}`$ with $`M`$ non-zero entries, by computing $`G_{ij}(E=0)`$ for different combinations of parameters $`M`$ and $`g`$. We The results are presented in Fig. 2. Again, the predictions of the supersymmetry method agree well with the numerical results.

## III Moments of power

In this Section we analyze the first two moments of power $`y`$ using two different approaches, both based on the RMT. These moments as obtained using the supersymmetric calculation (Eqs. (11), (12)) will be compared to those obtained using the naive approach. We first consider uniform damping and then $`M0`$.

In the simpler, but inexact, approach $`G_{ij}`$ is constructed as a modal sum,

$$G_{ij}\left(E\right)=\underset{r}{}\frac{u_i^ru_j^r}{EE_ri\zeta _r},$$

and then averaged using the eigenstatistics of the undamped GUE system. Here $`u^r`$ is the $`r`$th eigenmode and we call the imaginary part $`\zeta _r`$ of the eigenenergy $`E_r`$ the resonance width.

We first consider the case when $`\zeta _r`$ is uniform: $`\zeta _r=\epsilon ,`$ for all $`r`$ :

$$y\left(\pi \nu \right)^2=\underset{r}{}\underset{m}{}\frac{u_i^ru_j^r}{EE_ri\epsilon }\frac{u_i^mu_j^m}{EE_m+i\epsilon }.$$
(18)
On averaging Eq. (18) over the eigenmodes $`u^r`$, and assuming they are uncorrelated, $`y`$ becomes:

$$y\left(\pi \nu \right)^2=\underset{r}{}\frac{\left|u\right|^2^2}{\left(EE_ri\epsilon \right)\left(EE_r+i\epsilon \right)}.$$

The summation over the eigenenergies $`E_r`$ is then replaced with an integral ($`_rN\nu 𝑑E_r`$)

$$y=\frac{\left|u\right|^2^2}{\left(\pi \nu \right)^2}_{\mathrm{}}^+\mathrm{}\frac{\left(N\nu \right)dE_r}{\left(EE_ri\epsilon \right)\left(EE_r+i\epsilon \right)},$$

where $`\nu =\sqrt{4E^2}/\left(2\pi \right)`$ is the GUE modal density. Therefore, in the uniform damping case, the naive procedure produces:

$`y`$ $`=\pi \nu {\displaystyle \frac{\left|u\right|^2^2}{\left(\pi \nu \right)^2}}{\displaystyle \frac{\pi }{\epsilon }}`$ (19)
$`=N^2\left|u\right|^2^2{\displaystyle \frac{2}{2\pi \nu N\epsilon }}={\displaystyle \frac{2}{ϵ}},`$
where $`\left|u\right|^2`$ has been set to $`1/N`$ (by normalization).

The second moment of power is calculated in Appendix B by means of the same approach, and is given by:

$`y^2`$ $`={\displaystyle \frac{1}{\pi ^3\nu ^3N^3\epsilon ^3}}+{\displaystyle \frac{1}{4\pi ^4\nu ^4N^4\epsilon ^4}}`$ (20)
$`\times \left[1+8\pi ^2\left(N\nu \right)^2\epsilon ^2\mathrm{exp}\left(4\pi N\nu \epsilon \right)\right]`$
For the uniform damping case ($`\zeta _r=\epsilon `$ for all modes), application of the results of Section II is especially straightforward. As already discussed, Eq. (11) shows that in our model this case is realized either by setting $`M=0`$ with finite $`ϵ`$ or by letting $`M`$ be large and $`\gamma `$ be small, such that $`\gamma M`$is finite. One can use Eq. (13) for this purpose, but it is more convenient to start with the first part of Eq. (11). Setting $`M=0`$ and introducing:

$`𝔣=`$ $`{\displaystyle _1^{\mathrm{}}}𝑑\lambda _1{\displaystyle _1^1}𝑑\lambda _2\delta \left(y+1\lambda _1^2\right)`$
$`\times {\displaystyle \frac{\lambda _1^2\lambda _2^2}{\left(\lambda _1\lambda _2\right)^2}}\mathrm{exp}\left\{ϵ\left(\lambda _1\lambda _2\right)\right\},`$
we integrate by parts in Eq. (11):

$$𝒫\left(y\right)=\delta \left(y\right)+\left(\frac{d}{dy}+y\frac{d^2}{dy^2}\right)𝔣,$$

$$y^n=_0^{\mathrm{}}y^n𝒫\left(y\right)𝑑y=y^n\left(n𝔣+y𝔣_y^{}\right)_0^{\mathrm{}}$$
(21)
$$+n^2_0^{\mathrm{}}y^{n1}𝔣𝑑y=n^2_0^{\mathrm{}}y^{n1}𝔣𝑑y.$$

Integration with respect to $`y`$ in Eq. (21) eliminates the delta function, and gives:

$`y`$ $`=`$ $`{\displaystyle _1^{\mathrm{}}}𝑑\lambda _1{\displaystyle _1^1}𝑑\lambda _2\mathrm{exp}\left\{ϵ\left(\lambda _1\lambda _2\right)\right\}`$ (22)
$`\times `$ $`\left({\displaystyle \frac{\lambda _1+\lambda _2}{\lambda _1\lambda _2}}\right)={\displaystyle \frac{2}{ϵ}},`$
$`y^2`$ $`=`$ $`4{\displaystyle _0^{\mathrm{}}}y𝔉𝑑y={\displaystyle \frac{4}{ϵ^4}}\left(1\mathrm{exp}\left\{2ϵ\right\}\right)`$ (23)
$`+{\displaystyle \frac{8}{ϵ^2}}\left(1+{\displaystyle \frac{1}{ϵ}}\right),`$
which takes the same form as Eqs. (19) and (20) upon substitution of $`2\pi \nu N\epsilon `$ for $`ϵ`$. Thus, for uniform damping, the ’naive’ and supersymmetric methods agree, for both $`y`$ and $`y^2`$. This is not unexpected, because uniform damping with $`M=0`$ leaves eigenstatistics identical to those of closed systems, merely shifting all eigenenergies by $`i\epsilon `$. The results (22) and (23) can readily be reproduced by using $`P(y)`$ as given by Eq. (13). These moments are plotted in Fig. 3 together with the results of numerical simulations. The first two moments of $`y`$ were obtained numerically, by inverting matrix $`EI+i\mathrm{\Gamma }H`$ for each member of the ensemble. More precisely - we computed the column vector $`G_{ij}(E)`$ ($`j`$ is fixed, $`i=1,..N`$) by solving the algebraic equations:$`[EI+i\mathrm{\Gamma }H]G_{ij}=\delta _{ij}`$ for a fixed value of $`E`$. Repeating this procedure $`1500`$ times and averaging over the ensemble of $`H`$ and over the $`N1`$ values of $`ij`$, we obtained $`y`$ and $`y^2`$ for $`ϵ=1,2,4,6`$. As seen in Fig. 3, the correspondence is excellent.

Our next goal is to calculate $`y`$ and $`y^2`$ for the system with $`M`$ equivalent dampers and without additional uniform damping (i.e. with $`ϵ=0`$). For this problem the probability distribution function of scaled power is given by Eq. (12) and closed form expressions for the mean value of $`y`$ and its variance are cumbersome. It is, therefore, reasonable to carry out the corresponding integrations numerically. In contrast, the naive averaging, which now also includes integration over the resonance widths $`\zeta _r`$, distributed  according to $`\chi ^2`$ distribution:

$$p\left(\frac{\zeta _r}{\overline{\mathrm{\Gamma }}}\right)=\frac{M^M}{\mathrm{\Gamma }\left(M\right)}\left(\frac{\zeta _r}{\overline{\mathrm{\Gamma }}}\right)^{M1}\mathrm{exp}\left\{M\frac{\zeta _r}{\overline{\mathrm{\Gamma }}}\right\},$$

where $`\overline{\mathrm{\Gamma }}`$ is mean resonance width, produces a relatively compact answer for the statistics of $`y`$ (see Appendix B):

$$y\left(\pi \nu \right)^2=\left(\pi \nu \right)^2\frac{2}{ϵ_M}\frac{M}{M1},$$
(24)
$`y^2`$ $`={\displaystyle \frac{4M^2}{\left(M2\right)\left(M1\right)ϵ_M^2}}f(ϵ_M,M)`$
$`f(ϵ_M,M)=2+{\displaystyle \frac{4M\left(M1\right)}{\left(M3\right)ϵ_M}}`$ (25)
$`{\displaystyle \frac{M\left(4Mϵ_M4ϵ_MM\right)}{\left(2M3\right)ϵ_M^2}}`$
$`{\displaystyle \frac{\left(ϵ_M+M\right)^4}{M^2\left(2M3\right)ϵ_M^2}}\left(1+{\displaystyle \frac{ϵ_M}{M}}\right)^{42M},`$
where $`ϵ_M=2\pi \nu N\overline{\mathrm{\Gamma }}`$. We note that the $`\chi ^2`$ distribution for $`\zeta _r`$, is strictly correct only for the case of $`\overline{\mathrm{\Gamma }}`$ much less than mean level distance. It does, however, correspond well with the actual distribution  for more arbitrary value of $`\overline{\mathrm{\Gamma }}`$, as long as $`M`$ is large.

In order to compare Eqs. (24) and (25) to corresponding results obtained by numerical integration of Eq. (12), it is necessary to establish how $`\overline{\mathrm{\Gamma }}`$ is related to the parameters $`N`$, $`M`$ and $`g`$ of the supersymmetric calculation. By definition, the mean scaled resonance width in open systems is proportional to the product of modal density $`\nu `$ and average resonance width $`\overline{\mathrm{\Gamma }}`$ . Under the condition of uniform damping, the eigenmodes are equally damped, and the latter quantity is just equal to the individual damper strength $`\epsilon `$. In general, the relationship is not that simple and is given by Moldauer-Simonius formula :

$`\overline{\gamma }`$ $`=`$ $`2\pi {\displaystyle \frac{\overline{\mathrm{\Gamma }}}{\mathrm{\Delta }}}=2N\nu \pi \overline{\mathrm{\Gamma }}={\displaystyle \frac{1}{2}}{\displaystyle \underset{a}{\overset{M}{}}}\mathrm{ln}{\displaystyle \frac{g_a+1}{g_a1}},`$ (26)
$`g_a={\displaystyle \frac{1}{2\pi \nu }}\left({\displaystyle \frac{1}{\gamma _a}}+\gamma _a\right),`$
for the mean scaled resonance width, where $`\mathrm{\Delta }=1/\nu N`$ is the mean eigenenergy spacing and $`\gamma _a`$ is the coupling constant of $`a`$th channel. Note that our definitions of $`\overline{\gamma }`$ and eigenwidth $`\zeta _r`$ are different by factor of two from notation of Ref. . For the uniformly damped system we find:

$$\overline{\gamma }\frac{1}{2}\underset{a}{\overset{N}{}}\mathrm{ln}\frac{\frac{1}{2\pi \nu \epsilon }+1}{\frac{1}{2\pi \nu \epsilon }1}2\pi \nu \epsilon Nϵ,$$
(27)
which was also shown at the end of Section II. Thus, we conclude that our parameter $`ϵ`$ coincides with the mean scaled resonance width in the limit of large number of equivalent weak channels. Moreover, parameter $`ϵ_M`$ in Eqs. (24), (25) is the same as $`\overline{\gamma }`$.

The averaging, performed in Ref.  for the uniformly damped GOE system led to the dependence of the first two moments of the transmitted ultrasonic power on a single structural parameter, the modal overlap $``$. $``$ was defined in Ref.  in terms of average imaginary part of the eigenfrequency $`\omega _r`$ and modal density: $`=2\pi Im\omega _r\left(N/\omega \right)`$. We see that modal overlap may be identified with $`\overline{\gamma }`$.

$`y`$ and $`y^2`$ as predicted by supersymmetric calculation (from Eq. (12)) and by ’naive’ averaging (Eqs. (24) and (25)) are compared in Fig. 4 and Fig. 5 with numerical results for several different values of $`M`$ and $`\overline{\gamma }`$. The prediction by the supersymmetry method agrees with numerical results. In contrast, the results of the ’naive’ averaging underestimate both first and second moments of the power, except for very large $`M`$, close to the uniform damping case.

Finally, in connection with discussion of Ref.  we present the comparison of the relative variance ($`Relvar=y^2/y^21`$) of power in Fig. 6, and compare the supersymmetric and naive predictions.

## IV Conclusions

We investigated the statistical behavior of the power transmitted in a closed RMT system with internal dissipation, or an open RMT system coupled to the exterior via a finite number of equally strong channels. Using the supersymmetry method for systems with broken time reversal invariance we derived an expression for the probability distribution function for this quantity and studied its first two moments. The theoretical predictions were compared to the results of numerical simulations on GUE systems with dissipation, and to the results of a ’naive’ theory based on the RMT eigenstatistics of a closed non-dissipative system. The results of the supersymmetric calculation agree with the numerical data for the full range of parameters studied.

The naive averaging predictions are in general inconsistent with numerical results, because its assumptions ($`\chi ^2`$ distribution of resonance widths and decoupled uncorrelated Gaussian eigenmode amplitudes) follow from the first order perturbation theory, valid for small scaled resonance width. However, because the $`\chi ^2`$ distribution reduces to the exact distribution for the case of uniform damping or in the limit of a large number of weak channels, the naive theory is accurate in this limit, for all values of scaled level width.

###### Acknowledgements.

This work was supported by grants from the National Science Foundation \[CMS 99-88645 and CMS 0201346\], by computational resources from the National Center of Supercomputing Applications, by EPSRC grant GR/R13838/01 "Random matrices close to unitary or Hermitian" and by Vice-Chancellor grant from Brunel University. IR would like to thank International Programs in Engineering in University of Illinois at Urbana-Champaign for financial support, the Department of Mathematical Sciences at Brunel University for financial support and hospitality during his visit and Antonio M. Garcia-Garcia for valuable suggestions.

## Appendix A Evaluation of the superintegral

The result of ’Gaussian’ integration in Eq. (II) has to be simplified. To bring the integrand into a form convenient for saddle point integration, we use the series of identities for the supermatrix $`f=L^{1/2}\stackrel{~}{f}L^{1/2}`$ ($`L^2=I_4`$)

$`\stackrel{~}{f}`$ $`=EI_4I_NRI_Ni\mathrm{\Lambda }\mathrm{\Gamma }`$
$`=\left(I_NI_4i\mathrm{\Gamma }\left(\mathrm{\Lambda }𝔊^1\right)\right)𝔊,`$ (28)
$`\stackrel{~}{f}^1`$ $`=𝔊^1\left(I_4I_N\left(\mathrm{\Lambda }𝔊^1\right)i\mathrm{\Gamma }\right)^1`$
$`=𝔊^1{\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}\left(i\right)^k\left[\left(\mathrm{\Lambda }𝔊^1\right)i\mathrm{\Gamma }\right]^k,`$
where $`𝔊=EI_4R`$ was introduced. Supermatrix $`\stackrel{~}{f}^1`$ is diagonal in $`i`$ and $`j`$, thus:

$$f^{bb}=L^{1/2}\stackrel{~}{f}^{bb}L^{1/2},f_{12,bb_{ii}}f_{21,bb_{jj}}=𝔊_{12,bb}𝔊_{21,bb}.$$

Substituting $`Sdetf^1=\mathrm{exp}\left\{\underset{𝑖}{Str}\mathrm{ln}\stackrel{~}{f}\right\}`$ into the result of Gaussian integration with respect to the supervector components, and considering

$`T^n_H`$ $`=\left(n!\right)^2{\displaystyle \left[dR\right]\left(𝔊_{12,bb}^1𝔊_{21,bb}^1\right)^nSdetf^1}`$
$`\times \mathrm{exp}\left\{{\displaystyle \frac{N}{2}}StrR^2+i\epsilon NStrR\mathrm{\Lambda }\right\},`$ (29)
we separated the terms in the exponent according to their order in $`N`$ and obtained Eq. (8):

$`\left[R\right]={\displaystyle \frac{1}{2}}StrR^2+Str\mathrm{ln}𝔊,`$
$`\delta =i\epsilon NStrR\mathrm{\Lambda }Str\mathrm{ln}\left[I_Ni\mathrm{\Gamma }\left(\mathrm{\Lambda }𝔊^1\right)\right]`$
$`=i\epsilon NStrR\mathrm{\Lambda }MStr\mathrm{ln}\left[I_4i\gamma \mathrm{\Lambda }𝔊^1\right].`$ (30)
The last identity was proved by expanding the logarithm into the series (see Ref. ).

After the Gaussian integration around the saddle point in Eq. (29) the probability distribution function for the scaled power is expressed as an integral over the manifold formed by supermarices $`Q=𝔗^1\mathrm{\Lambda }𝔗`$:

$$𝒫\left(y\right)=\left[dQ\right]\delta \left(yQ_{12,bb}Q_{21,bb}\right)Y\left(Q\right),$$

$`Q`$ is parametrized by four commuting variables $`\lambda _1`$, $`\lambda _2`$, $`\mu _1`$, $`\mu _2`$ and four anticommuting $`\alpha `$, $`\alpha ^{}`$, $`\beta `$, $`\beta ^{}`$ key-4 ; key-17 . Here $`\lambda _1(1,\mathrm{})`$, $`\lambda _2(1,1)`$, and $`\left|\mu _1\right|^2=\lambda _1^21`$, $`\left|\mu _2\right|^2=1\lambda _2^2`$. We can also introduce another set of variables, according to $`\lambda _1=\mathrm{cosh}\theta _1`$, $`\mu _1=\mathrm{sinh}\theta _1\mathrm{exp}\left\{i\varphi _1\right\}`$, $`\lambda _2=\mathrm{cos}\theta _2`$, $`\mu _2=\mathrm{sin}\theta _2\mathrm{exp}\left\{i\varphi _2\right\}`$, where $`\theta _1(0,\mathrm{})`$, $`\theta _2(0,\pi )`$, $`\varphi _1`$, $`\varphi _2(0,2\pi )`$. Next we observe

$`StrQ\mathrm{\Lambda }=2i\left(\lambda _1\lambda _2\right),Y\left(Q\right)=S\stackrel{M}{det}\left[I_4+i{\displaystyle \frac{E}{2}}\gamma \mathrm{\Lambda }+i\pi \nu \gamma Q\mathrm{\Lambda }\right]\mathrm{exp}\left\{i\epsilon \pi \nu NStrQ\mathrm{\Lambda }\right\}`$ (31)
$`=`$ $`\left({\displaystyle \frac{1+2\pi \nu \gamma \lambda _2+\gamma ^2}{1+2\pi \nu \gamma \lambda _1+\gamma ^2}}\right)^M\mathrm{exp}\left\{2\epsilon \pi \nu NStr\left(\lambda _1\lambda _2\right)\right\}=\left({\displaystyle \frac{g+\lambda _2}{g+\lambda _1}}\right)^M\mathrm{exp}\left\{2\epsilon \pi \nu NStr\left(\lambda _1\lambda _2\right)\right\},`$
where $`g=\left(1/\gamma +\gamma \right)/\left(2\pi \nu \right)`$,

$`Q_{12,bb}=\mu _1\left(1\alpha ^{}\alpha /2\right)\left(1+\beta ^{}\beta /2\right)\alpha ^{}\beta \mu _2^{},`$
$`Q_{21,bb}=\mu _1^{}\left(1\alpha ^{}\alpha /2\right)\left(1+\beta ^{}\beta /2\right)+\alpha \beta ^{}\mu _2,`$
$`Q_{12,bb}Q_{21,bb}=\left|\mu _1\right|^2+\left|\mu _1\right|^2\alpha ^{}\beta ^{}\alpha \beta `$
$`+\left|\mu _2\right|^2\alpha ^{}\beta ^{}\alpha \beta +\left|\mu _1\right|^2\left(\beta ^{}\beta \alpha ^{}\alpha \right)`$
$`+\left|\mu _1\right|\left|\mu _2\right|e^{i\left(\varphi _1+\varphi _2\right)}\alpha \beta ^{}`$ (32)
$`\left|\mu _1\right|\left|\mu _2\right|e^{i\left(\varphi _1+\varphi _2\right)}\alpha ^{}\beta ,`$
and the integration measure was defined as

$$dQ=\frac{d\alpha ^{}d\beta ^{}d\alpha d\beta d\lambda _1d\lambda _2d\varphi _1d\varphi _2}{\left(2\pi \left(\lambda _1\lambda _2\right)\right)^2}.$$

Substituting Eqs. (31), (32) we proceed with integration with respect to Grassman variables. First, we need to expand the delta function retaining only the terms of zero and maximum order in these variables key-4 ; key-17 ; key-18 . Setting $`Q_{12,bb}Q_{21,bb}=1+\lambda _1^2+z`$ we expand the delta function,

$$\delta \left(yQ_{12}^{bb}Q_{21}^{bb}\right)=\delta \left(y+1\lambda _1^2z\right)=\delta \left(y+1\lambda _1^2\right)$$

$$\left(\delta _z^{}\left(y+1\lambda _1^2\right)+\delta _{zz}^{\prime \prime }\left(y+1\lambda _1^2\right)\left(1\lambda _1^2\right)\right)$$

$$\times (\lambda _1^2\lambda _2^2)\alpha ^{}\beta ^{}\alpha \beta +\mathrm{}=\delta (y+1\lambda _1^2)$$

$$+\left(\frac{d}{dy}+y\frac{d^2}{dy^2}\right)\delta \left(y+1\lambda _1^2\right)\alpha ^{}\beta ^{}\alpha \beta +\mathrm{},$$

where we used the fact that the argument of delta function is linear in $`y`$ ($`z`$ and $`1\lambda _1^2`$), in order to be able to take the differential operator out of the integral. Then, we calculate the integral over $`\varphi _1,\varphi _2`$ and Grassman variables,

$`𝒫\left(y\right)`$ $`=\delta \left(y\right)+\left({\displaystyle \frac{d}{dy}}+y{\displaystyle \frac{d^2}{dy^2}}\right)`$
$`\times `$ $`{\displaystyle _1^{\mathrm{}}}𝑑\lambda _1{\displaystyle _1^1}𝑑\lambda _2\delta \left(y+1\lambda _1^2\right)`$ (33)
$`\times `$ $`{\displaystyle \frac{\lambda _1^2\lambda _2^2}{\left(\lambda _1\lambda _2\right)^2}}\mathrm{exp}\left\{ϵ\left(\lambda _1\lambda _2\right)\right\}\left({\displaystyle \frac{g+\lambda _2}{g+\lambda _1}}\right)^M,`$
where $`ϵ=2\epsilon \pi \nu N`$ and we have used

$`{\displaystyle }`$ $`\left[dQ\right]\mathrm{exp}\left\{ϵ\left(\lambda _1\lambda _2\right)\right\}`$
$`\times `$ $`\left({\displaystyle \frac{g+\lambda _2}{g+\lambda _1}}\right)^M\delta \left(y+1\lambda _1^2\right)=\delta \left(y\right),`$
for the ’Efetov-Wegner’ term . Integration with respect to $`\lambda _1`$ is simple because of the presence of the delta function.

$`𝒫\left(y\right)`$ $`=\delta \left(y\right)+\left({\displaystyle \frac{d}{dy}}+y{\displaystyle \frac{d^2}{dy^2}}\right){\displaystyle \frac{\mathrm{exp}\left\{ϵ\sqrt{1+y}\right\}\theta \left(y\right)}{2\sqrt{1+y}\left(g+\sqrt{1+y}\right)^M}}`$
$`\times {\displaystyle _1^1}d\lambda _2{\displaystyle \frac{\sqrt{1+y}+\lambda _2}{\sqrt{1+y}\lambda _2}}\mathrm{exp}\left\{ϵ\lambda _2\right\}(g+\lambda _2)^M`$
$`=\delta \left(y\right)+\left({\displaystyle \frac{d}{dy}}+y{\displaystyle \frac{d^2}{dy^2}}\right)\theta \left(y\right)𝔉(y).`$
Using property $`\delta \left(y\right)=y\delta ^{}\left(y\right)`$, we arrive at:

$`𝒫\left(y\right)`$ $`=\delta \left(y\right)+\delta \left(y\right)y{\displaystyle \frac{d}{dy}}𝔉(y)`$ (34)
$`+\left({\displaystyle \frac{d}{dy}}+y{\displaystyle \frac{d^2}{dy^2}}\right)𝔉(y),`$
Eq. (34) completes the calculation of probability distribution function for the scaled power. This equation yields Eq. (11) upon substitution of: $`lim_{y0}y\frac{d}{dy}𝔉(y)=1`$.

Finally, we set $`ϵ=0`$ and derive Eq. (12). We note that integral in $`𝔉(y)`$ is a table integral:

$`𝒫\left(y\right)`$ $`=\left({\displaystyle \frac{d}{dy}}+y{\displaystyle \frac{d^2}{dy^2}}\right)`$
$`\times {\displaystyle \frac{\left(g+1\right)^{M+1}f\left(g+1\right)\left(g1\right)^{M+1}f\left(g1\right)}{\left(M+1\right)\left(g+\sqrt{y+1}\right)^{M+1}}},`$
$`f\left(u\right)`$ $`=_2F_1(M+1,1,M+2,{\displaystyle \frac{u}{g+\sqrt{y+1}}}).`$
Thus, it is possible to apply differential operator to obtain the final form of $`𝒫\left(y\right)`$ for this case. However, we notice:

$`𝔉(y)=`$ $`{\displaystyle \frac{1}{\left(g+\sqrt{y+1}\right)^M}}{\displaystyle _1^1}{\displaystyle \frac{\left(g+\lambda _2\right)^M}{\lambda _2\sqrt{y+1}}}𝑑\lambda _2`$ (35)
$``$ $`{\displaystyle \frac{1}{2\sqrt{y+1}\left(g+\sqrt{y+1}\right)^M}}{\displaystyle _1^1}\left(g+\lambda _2\right)^M𝑑\lambda _2.`$
The second term in the above equation can be evaluated immediately, while for the first one we can use an identity:

$`\left(g+\lambda _2\right)^M=`$ $`{\displaystyle \underset{m=0}{\overset{M}{}}}\left(g+\sqrt{y+1}\right)^m`$
$`\times `$ $`\left(\lambda _2\sqrt{y+1}\right)^{Mm}{\displaystyle \frac{M!}{m!\left(Mm\right)!}}`$
We integrate each term in Eq. (35) separately, and after the first differentiation of the result with respect to $`y`$, the series can be summed back, so that the remaining procedure becomes straightforward and leads to Eq. (12).

## Appendix B Moments calculation for the case of $`M`$ equivalent dampers

In this Appendix we demonstrate the intermediate steps leading to the Eqs. (20), (24), (25). We start with the modal expansion for $`T^2`$ without making an assumption about uniform damping:

$`y^2\left(\pi \nu \right)^4`$ $`={\displaystyle \underset{r,m,l,k}{}}{\displaystyle \frac{u_i^ru_j^r}{EE_ri\zeta _r}}{\displaystyle \frac{u_i^mu_j^m}{EE_m+i\gamma _m}}`$ (36)
$`\times {\displaystyle \frac{u_i^lu_j^l}{EE_li\zeta _l}}{\displaystyle \frac{u_i^ku_j^k}{EE_k+i\gamma _k}}.`$
Absence of correlation between different eigenmodes produces the following result for the variance of $`y`$:

$`y^2\left(\pi \nu \right)^4`$ $`={\displaystyle \underset{r}{}}{\displaystyle \frac{\left|u_i^r\right|^4\left|u_j^r\right|^4}{\left(EE_ri\zeta _r\right)^2\left(EE_r+i\zeta _r\right)^2}}+{\displaystyle \underset{rl}{}}{\displaystyle \frac{\left|u_i^r\right|^2\left|u_j^r\right|^2\left|u_i^l\right|^2\left|u_j^l\right|^2}{\left(EE_ri\zeta _r\right)^2\left(EE_l+i\zeta _l\right)^2}}`$
$`+{\displaystyle \underset{rl}{}}{\displaystyle \frac{\left|u_i^r\right|^2\left|u_j^r\right|^2\left|u_i^l\right|^2\left|u_j^l\right|^2}{\left(EE_ri\zeta _r\right)\left(EE_r+i\zeta _r\right)\left(EE_li\zeta _l\right)\left(EE_l+i\zeta _l\right)}}`$ (37)
Next, we replace summation over $`E_r`$ and $`E_l`$ with integration ($`_rN\nu 𝑑E_r`$) and take into account the correlation between the GUE eigenvalues in Eq. (B) by introducing the factor $`1Y_2\left(\pi N\nu \left(E_rE_l\right)\right)`$, :

$`y^2`$ $`={\displaystyle \frac{N\nu \left|u\right|^4^2}{\left(\pi \nu \right)^4}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{dx}{\left(xi\zeta _r\right)^2\left(x+i\zeta _r\right)^2}}+{\displaystyle \frac{2\left(N\nu \right)^2\left|u\right|^2^4}{\left(\pi \nu \right)^4}}`$ (38)
$`\times {\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{\left(1Y_2\left(\pi N\nu z\right)\right)dxdz}{\left(x^2+\zeta _r\right)\left(\left(xz\right)^2+\zeta _l^2\right)}}+{\displaystyle \frac{\left(N\nu \right)^2\left|u\right|^2^4}{\left(\pi \nu \right)^4}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{\left(1Y_2\left(\pi N\nu z\right)\right)dxdz}{\left(xi\zeta _r\right)^2\left(xz+i\zeta _l\right)^2}},`$
where $`x=`$

$`EE_r`$, $`z=E_rE_l`$ and the Dyson two-level correlation function for the GUE is $`Y_2\left(\xi \right)=\left(\mathrm{sin}\xi /\xi \right)^2`$. Integration over $`x`$ and $`z`$ in Eq. (38) for the case of uniform damping $`\zeta _r=\zeta _l=\epsilon `$ yields Eq. (20):

$`y^2={\displaystyle \frac{N\nu \left|u\right|^4^2}{\left(\pi \nu \right)^4}}{\displaystyle \frac{\pi }{2\zeta _r^3}}+{\displaystyle \frac{2\left(N\nu \right)^2\left|u\right|^2^4}{\left(\pi \nu \right)^4}}`$
$`\times {\displaystyle \frac{\zeta _r+\zeta _l}{2\zeta _r\zeta _l}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{\left(1Y_2\left(\pi N\nu z\right)\right)dz}{z^2+\left(\zeta _r+\zeta _l\right)^2}}`$ (39)
$`=\left|u\right|^4^2{\displaystyle \frac{N}{2\left(\pi \nu \right)^3\epsilon ^3}}+\left|u\right|^2^4{\displaystyle \frac{1}{4\left(\pi \nu \right)^4\epsilon ^4}}`$
$`\times \left[1+8\left(N\nu \right)^2\pi ^2\epsilon ^2\mathrm{exp}\left\{4N\nu \pi \epsilon \right\}4N\nu \pi \epsilon \right],`$
upon substitution $`\left|u\right|^4/\left|u\right|^2^2=2`$, (as is the case for complex Gaussian random numbers) and $`\left|u\right|^2=1/N`$.

In the case of finite number $`M`$ of weak dampers ($`ϵ=0`$) the ensemble averaging includes an integration with over a distribution of widths, given by :

$$p\left(\frac{\zeta _r}{\overline{\mathrm{\Gamma }}}\right)=\frac{M^M}{\mathrm{\Gamma }\left(M\right)}\left(\frac{\zeta _r}{\overline{\mathrm{\Gamma }}}\right)^{M1}\mathrm{exp}\left\{M\frac{\zeta _r}{\overline{\mathrm{\Gamma }}}\right\},$$
(40)
where $`\overline{\mathrm{\Gamma }}`$ is average resonance width and $`\mathrm{\Gamma }\left(M\right)`$ is a Gamma function. Starting with Eq. (21) we average $`y`$ over the eigenmodes and eigenenergies to get:

$$y\left(\pi \nu \right)^2=N\nu \left|u\right|^2^2\frac{\pi }{\zeta _r}.$$
(41)
Which becomes Eq. (24) upon integration over $`p\left(\zeta _r\right)`$. We indicate this averaging by overbar:

$$y\left(\pi \nu \right)^2=\overline{N\nu \left|u\right|^2^2\frac{\pi }{\zeta _r}}=\left|u\right|^2^2\frac{\pi N\nu }{\overline{\mathrm{\Gamma }}}\frac{M}{M1}.$$

To evaluate the second moment of power we integrate over $`x`$, and then over $`\zeta _r`$ and $`\zeta _l`$ in Eq. (38):

$`y^2\left(\pi \nu \right)^4`$ $`=\overline{N\nu \left|u\right|^4^2{\displaystyle \frac{\pi }{\zeta _r^3}}}+\overline{2\left(N\nu \right)^2\left|u\right|^2^4I}`$
$`I=`$ $`{\displaystyle \frac{\zeta _r+\zeta _l}{2\zeta _r\zeta _l}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{\left(1Y_2\left(\pi N\nu z\right)\right)dz}{z^2+\left(\zeta _r+\zeta _l\right)^2}}.`$ (42)
The remaining integral with respect to $`z`$ in $`\overline{I}`$ it is convenient to use the Fourier transform of $`Y_2\left(\xi \right)`$, which has a form:

$`b\left(q\right)`$ $`={\displaystyle _{\mathrm{}}^{\mathrm{}}}Y_2\left(\xi \right)\mathrm{exp}\left\{2\pi i\xi q\right\}𝑑\xi =1\left|q\right|,\left|q\right|1,`$
$`b\left(q\right)`$ $`=0,\left|q\right|1,`$
$`\overline{I}=`$ (43)
$`=`$ $`\overline{{\displaystyle \frac{\zeta _r+\zeta _l}{4\zeta _r\zeta _l}}{\displaystyle _1^1}\left(1{\displaystyle \frac{\left|q\right|}{2\pi }}\right)\mathrm{exp}\left\{N\nu \left(\zeta _r+\zeta _l\right)\right\}𝑑q}.`$
The average over $`\zeta _r`$ and $`\zeta _l`$ is now straightforward. Finally, substituting the result into the Eq.(B), and taking into account $`\left|u\right|^4/\left|u\right|^2^2=2`$, we obtain the second moment in its closed form (Eq. (25)). The derivation presented above assumes that resonance widths $`\zeta _r`$ and eigenmodes $`u^r`$ are statistically independent.