Distribution of zeros of the S-matrix of chaotic cavities with localized losses and Coherent Perfect Absorption: non-perturbative results

We employ the Random Matrix Theory framework to calculate the density of zeroes of an $M$-channel scattering matrix describing a chaotic cavity with a single localized absorber embedded in it. Our approach extends beyond the weak-coupling limit of the cavity with the channels and applies for any absorption strength. Importantly it provides an insight for the optimal amount of loss needed to realize a chaotic coherent perfect absorbing (CPA) trap. Our predictions are tested against simulations for two types of traps: a complex network of resonators and quantum graphs.

One of the undisputed successes of the Random Matrix Theory (RMT) in last decades is its faithful description of scattering properties of disordered or/and chaotic cavities ‡ verified in a number of experiments: distribution of conductances, delay times, resonance widths, the consequences of preserved/violated time reversal symmetry, the presence of uniform absorption etc, see e.g. recent review articles [1,2,3,4,5].
Inspired by all these successes it is natural to expect that RMT should play a role among essential tools employed for the understanding of a new category of scattering problems, associated with the design of chaotic coherent perfect absorbers (CPA). CPAs are typically weakly lossy cavities which act as perfect constructive interference traps for incident coherent radiation [6]. CPA protocols [6,7] have been found increasingly relevant because of our capabilities in manipulating incoming wavefront shapes in frameworks ranging from acoustics, to microwaves and optics. Technological applications of CPAs might cover areas as diverse as optical light-light control and switching, plasmonic sensors, lasing/anti-lasing cavities [6,8,9,10,11,12] to THz and radiofrequency wave trapping and stabilization of wireless communications [13,14] and sound absorption [15]. Despite the broad range of technological applications the study of CPA was confined up to now to the investigation of simple (integrable) cavities. Only recently researchers started investigating the effects of complexity (chaos) in the realization of CPAs [16]. This is quite surprising since at the heart of CPA protocols are wave interference phenomena which are abundantly present in wave chaotic cavities. In fact one can further capitalize on these complex (chaotic) interferences in order to provide universal predictions for properties of chaotic CPA cavities.
In this paper we are utilizing the RMT toolbox in order to shed new light on the realization of chaotic CPAs with M open channels and spatially non-uniform losses. For simplicity we shall assume that the losses are described by a single point-like absorber with loss-strength γ 0 embedded inside the cavity. Specifically we provide a detailed statistical description of the density of complex zeros of the scattering matrix S ( i.e. solutions of det S = 0 in the complex energy plane) which describes such chaotic CPA cavities, for a given strength of absorber γ 0 . In particular, the number of those Smatrix zeroes which with increased γ 0 move away from the upper-half plane and cross the real energy axis can be interpreted as the number of "perfect absorbing states" supported by such cavities. We provide analytical expressions for the density of zeros as a function of γ 0 for a fixed M and given universality class (violated/preserved timereversal invariance). We point out that CPAs associated with violated time-reversal symmetry are novel and do not fall in the original category of CPAs as defined in [6,7,8] the latter are typically recognized as the time-reversal of an outgoing lasing ‡ In the standard terminology of the field of quantum chaos a scattering domain is called chaotic (respectively, integrable) if after detaching all external leads or waveguides the dynamics of a classical particle bouncing inside the associated closed domain is chaotic (resp. integrable). The generally accepted conjecture is that highly excited eigenstates of the Laplacian operator on the domains corresponding to chaotic cavities are faithfully described by eigenvalues of large random matrices, whereas eigenspectra of integrable cavities follow very different (Poisson) statistics.
field. Nevertheless, this new family of incident coherent fields lead to perfect absorption as well.
In contrast to [16] we do not assume a weak-coupling limit, thus providing nonperturbative insights into the optimal amount of loss (perfect tuning) needed to realize a crossover of a zero from positive to negative complex energy semi-plane. Our results are tested against detailed numerical simulations using both RMT modeling and an actual wave chaotic network (quantum graphs).
The main object of interest is a (non-unitary) M -channel scattering matrix S (E, γ 0 ) which, in the presence of localized absorbers of strength γ 0 , connects incoming |I to outgoing |O wave amplitudes at a fixed real energy E via the relation |O = S|I . We start by a brief demonstration of the adjustments needed to be done in the approach of [17,18,19] in order to make it capable of describing the CPA problems.
In fact, one may conveniently think of a scattering system with M ≥ 1 open scattering channels and L ≥ 1 local absorbers of equal strength ( generalization to different strengths is obvious) as represented by a fictitious system of M = M + L open channels, such that its corresponding M × M scattering matrix S is unitary due to the flux conservation. The building blocks for the associated S-matrix are (i) an N × N (N 1) random Gaussian matrix H (real symmetric GOE, β = 1 or Hermitian GUE, β = 2) used to model the Hamiltonian (or, in context of classical waves, the wave operator) of an isolated chaotic cavity, (ii) an N × M matrix W (containing the coupling amplitudes to M open scattering channels) satisfying, for simplicity, the relation N l=1 W * cl W bl = γδ cb , ∀c, b = 1, . . . , M ('equivalent orthogonal channels assumption') and finally (iii) N × L matrix A of coupling amplitudes A li describing L absorptive channels and taken as another set of fixed orthogonal vectors satisfying N i=1 A li A mi = γ 0 δ cd , ∀(l, m) = 1, . . . , L. We also assume that the vector space spanned by coupling vectors A l , l = 1, . . . , L for the absorptive part is orthogonal to one spanned by vectors W b , b = 1, . . . M . In such an approach the (sub-unitary) scattering matrix S(E, γ 0 ) defined above is simply a M × M diagonal sub-block of the whole M × M scattering matrix S and can be written in the standard form, see [17,18,19]: where K A = W † 1 E1−H+iΓ A W and we defined an effective non-Hermitian Hamiltonian as By using the identity det (1 − P Q) = det (1 − QP ) one can further find from (1) that The real values E for which det S(E, γ 0 ) vanishes define the energies of incident travelling waves for which a CPA can be realized [6]. In fact, to have CPA one has to impose a condition that the incident waveform |I CPA must be realized as a linear combination of channel modes with amplitudes given by the components of the eigenvector of the scattering matrix S(E CPA , γ 0 ) associated with a zero eigenvalue. Indeed, in this case S(E CPA , γ 0 )|I CPA = 0 = |O implying that there is no outgoing wave, hence CPA. The condition of zero eigenvalue in turn implies det S(E, γ 0 ) = 0.
The equation (3) makes it clear that the complex zeroes of det S(E, γ 0 ) are eigenvalues of the non-Hermitian Hamiltonian H (a) As long as absorption is non-vanishing, i.e. Γ A = 0, those eigenvalues are situated in both positive and negative half-planes of the complex energy z = E + iY . When studying their distribution in the complex plane one can exploit the statistical rotational invariance of the Hamiltonian matrix H and safely replace the matrix Γ = W W † by a diagonal one γ M l=1 |l l| , where |l is assumed to be the eigenbasis of W W † . In addition, the physical assumption of no direct overlap between the scattering channels and the position of the absorbers allows us to consider the matrix Γ A to be diagonal as well. In the simplest representative example case L = 1 the matrix A is rank-one, and denoting the eigenvector corresponding to its only non-vanishing eigenvalue as |0 we finally see that zeroes z n = E n + iY n of det S(E, γ 0 ) are given by the N eigenvalues of the effective non-Hermitian Hamiltonian [16] H (a) The local point absorber is modeled by the last term in (4), with γ 0 > 0 standing for the loss-strength and the vector |0 (normalized as 0|0 = 1 and characterizing the position of the absorber) is assumed to be orthogonal to all channel vectors |l .
In what follows we will describe the statistics of the rescaled imaginary parts of complex zeros y n = βπY n /∆ = πβN ρ(E)Y n as a function of the relative strength of the absorber a ≡ γ 0 /γ. Here ∆ = (ρ(E)N ) −1 is the mean spacing between neighboring real eigenvalues for H close to the energy E and ρ(E) = 1 π 1 − E 2 /4 is the standard RMT density of eigenvalues of H. Our main object of interest is the probability density of the imaginary parts for complex zeroes (whose real parts are close to the energy E) defined as Perturbative treatment in the weak-coupling regime. Before presenting the non-perturbative results, valid for arbitrary coupling and absorption strength, we briefly consider the limit of weak coupling γ, γ 0 1 where the form of P M (y) can be found via a standard perturbative analysis, cf. [20,21]. In this regime one expects the zeroes to be located in the vicinity of the real axis in such a way that their 'height' in the complex plane satisfy Y n ∆. First order perturbation theory results in replacing E n → z n = E n + iY n with where |n stand for the eigenvectors of H. In the following we assume that E n is close to E = 0, so that the appropriately rescaled variables are y n = βN Y n . The probability density P M (y) = δ (y − y n ) H in this case can be further evaluated by using the fact that the projections n|l = u l + i(β − 1)v l of the eigenvectors of H in any arbitrary basis, for large N 1, can be effectively replaced with a set of independent, identically distributed mean-zero Gaussian-distributed variables u l , v l with variance 1/βN .
Performing the averaging over the eigenvector components, the distribution of the rescaled imaginary parts can be after some manipulations represented as where C , the relative strength of the absorber is conveniently characterized by the ratio a ≡ γ 0 /γ, and θ(x) stands here and henceforth for the Heaviside step function. Eq. (6) is a generalization of the famous Porter-Thomas distribution for M open channels, see [21] and references therein, and is reduced to it in the limit a → 0.
Equation (6) allow us to evaluate the probabilityΠ M (a) ≡ 0 −∞ P(y) dy that, for given a and M , a zero of the S-matrix crosses the real axis. One finds (7) for β = 2 and β = 1, correspondingly. In the small a-limit the latter becomes Equation (8) is relevant for CPA protocols, which usually are concerned with situations where the absorber is weak enough to absorb by itself the incident waves. Next let us consider a spectral strip in the complex energy plane parallel to the imaginary axis, with a width W around E = Rez = 0, so that it contains on average L = W/∆ complex zeroes.
The number of "perfect absorbing" states for a given set of parameters is given by the number L a of those zeroes which are below the real axis.
This number has the binomial distribution  M (a). Specifically for β = 1, we can use the above to calculate the probability P (β=1) (L a = 0) of having no "perfect absorbers" at all in the a 1 limit: which allow us to identify the scaling of the minimal a = γ 0 /γ associated with at least one CPA state. Setting L ∼ N and requesting that P(L a = 0) = O(1) we get that a min ∼ N −2/M . The large value of a min for M 1 is a direct consequence of the fact that, in this limit, the density of zeroes for γ 0 = 0 is strongly depleted. Hence one needs larger values of γ 0 in order to drive such zeroes through the real axis.
Non-perturbative treatment of the problem.
For a general coupling/absorption strengths beyond perturbation regime the density P M (y) can be inferred by a judicious re-interpretation of the RMT results for resonance statistics obtained originally in [22], and further elaborated in [19,23,24]. A fully controlled ab initio derivation is relatively involved even for the simplest β = 2 case, and will be published elsewhere [25].
As is well-known [17,18,19], as long as N M the main control parameters of the theory are the non-perturbative coupling constants 1 ≤ g < ∞ and 1 ≤ g 0 < ∞ defined in terms of 'bare' values 0 ≤ γ < ∞ and 0 ≤ γ 0 < ∞ as Correspondingly, for β = 2 the density of complex zeroes turns out to be given by One can easily calculate the mean imaginary part y for the zeros: where in the latter form the expression is also valid for β = 1 case considered later on. This formula provides a generalization of the well-known Moldauer-Simonius relation, see [22,19] and refs. therein. The logarithmic divergence of y for g → 1 or g 0 → 1 reflects the characteristic power-law tails P (β=2) M (y) ∼ y −2 in the upper/lower half plane typical for the 'perfect coupling' case.
The following remark is due here. Looking at the invariance of coupling constants Eq. (10) with respect to changing γ → 1/γ and γ 0 → 1/γ 0 one may naively think that there is a complete equivalence of the regimes of small and large γ's. The situation is however more subtle. If both parameters γ and γ 0 take values in the interval 0 ≤ γ, γ 0 < 1 the distribution (11) is asymptotically exact for all N zeroes. At the same time, for γ > 1 exactly M 'broad zeroes' separate and move far from the real axis in the upper half-plane, with their imaginary parts being non-random and given by Y upper = γ − γ −1 = O(1), see [19]. The overwhelming majority of N − M zeroes stay however close to the real axis at a distance of order of 1/N and their density is still faithfully described by the same distribution (11). Such a restructuring is natural to call the zeroes self-trapping phenomenon [16], and its analogue for the case of resonances is well-known and was observed experimentally [26]. Similarly, for γ 0 > 1 a single 'broad' zero with imaginary part Y lower = − γ 0 − γ −1 0 emerges in the lower half plane. Use of Eq. (11) allow us to calculate the probabilityΠ ≡ 0 −∞ P M (y) dy that, for given values of M, g and g 0 , a zero of the S-matrix crosses the real axis and hence be located in the lower half-energy plane: From Eq. (13) we immediately infer that whenever the absorber is not 'perfectly tuned', so that g 0 < 1, the probabilityΠ (β=2) decays exponentially with the number of open channels M 1 -in a way similar to the perturbative regime. This is a consequence of the 'zeroes depletion', or 'gap', arising for many open scattering channels in the upper half-plane close to real axis. A strikingly different behavior is observed in the case of "perfectly tuned" absorber g 0 = 1 when the probability to have zeroes below the real axis takes the formΠ (β=2) (g, and decays only as 1/(M + 1), irrespective of the coupling strength γ of the cavity with the scattering channels. In particular, a spectral window around E = 0 with a large enough number of levels (> M ) will on average contain a zero below the real axis. Hence, a single 'perfectly tuned' absorber dramatically increases the probability of realizing a CPA state. In fact, even in the case of slight detuning g 0 = 1 + δ M (where δ ∼ O(1) and M → ∞) we get from (13)Π(g, g 0 ) ∼ g+1 2(M +1) e −δ/(g+1) . Re-interpreting the results of [23] one may also provide exact non-perturbative expressions for the distribution of zeroes for β = 1. In this case, however, the formulas are quite cumbersome, even for the simplest M = 1 case. In the latter case the density of zeroes in the upper half-plane takes the form with F(λ) = ∞ g dp 1 e −yp 1 The density of zeroes for the lower half plane y < 0 is obtained by replacing y → −y in the density Eq. (14) and exchanging g ↔ g 0 everywhere.
In the limit of a single perfect channel g = 1 we can further proceed and evaluate the probabilityΠ 1 (g = 1, g 0 ) = 1− ∞ 0 P(y) dy that for a fixed absorptive coupling value 1 ≤ g 0 ≤ ∞ a given zero of the scattering matrix becomes negative: Finally, adapting the methods of [21] one can show that the exact non-perturbative distribution (14) takes the form (6) in the weak coupling limit g 1, g 0 1.
Simulations. We have tested the validity of the RMT predictions via direct diagonalization of the effective RMT Hamiltonian Eq. (4) and by an exact search of zeros of the scattering matrix describing an actual wave chaos system i.e. a complex network of one-dimensional waveguides coupled together via splitters. The latter system has been established during the last years as a prototype model for wave chaos studies [27,28] while its experimental realization in the microwave domain has been demonstrated by a number of groups [29,30,31,32].
The RMT model can be used to describe (in the coupled mode theory approximation [33]) a network of coupled cavities (acoustic/microwave or optical) or LC circuits which are randomly coupled with one another. The random coupling can be introduced by the randomness in the distances between the cavities (or by a random capacitive coupling in the case of LC circuits). The time-reversal invariance can be violated via magnetooptical effects in the case of electromagnetic network resonators [34] or via gyrators in the case of LC circuits [35]. The breaking of time-reversal symmetry is much more challenging in the acoustic domain but recent developments have demonstrated that it can be achieved by incorporating circulating fluid elements [36,37].
We start with the presentation of the RMT simulations. We have considered a small window around E ≈ 0 and we have collected at least 5000 S-matrix zeroes for statistical processing. In Figs. 1a,b,c we show the distribution of zeros P(y) for a GUE (β = 2) case with M = 2 number of channels and matrices of the size N = 600. The agreement between the numerical data and the theory is excellent when the zeroes selftrapping effect is absent (see Figs. 1a,b), for 0 < γ, γ 0 ≤ 1. For γ = 2, a = 2 Eq. (11) the predictions of RMT are still describing our data quite well on the scale of 1/N , see Fig. 1c, in agreement with the discussion of the resonance trapping phenomenon. In the insets of Figs. 1a,b,c we summarize our results for the probabilityΠ β=2 (γ, a) for fixed γ-values. The numerical analysis ofΠ β=1 (a) for the GOE case is shown in Figs. 2a,b and compared with Eqs. (S6) and (14) for weak and strong coupling respectively. Representative probability densities P(y) are shown in the insets of Fig. 2.
Next we proceed with the numerical analysis of zeroes for quantum graphs. The system consists of n = 1, · · · , V vertices connected by B bonds. The number of bonds emanating from a vertex n is the valency v n of the vertex and the total number of bonds can be expressed as B = 1 2 n v n . The length of the bonds l n,m are taken from a uniform box distribution l n,m ∈ [l 0 − w/2, l 0 + w/2] where l 0 = 1 is its mean and w = 1. On each bond, the component Ψ nm of the total wavefunction is a solution of the free Helmholtz equation i d dx + A 2 Ψ nm (x) = k 2 Ψ nm (x) where we have included a "magnetic field vector" A = A nm = −A mn which breaks the time-reversal symmetry.
At the vertices (where scattering events occur), we can introduce a δ-potential with a strength λ n . Furthermore the wavefunction Ψ n,m at the vertices must be continuous and must satisfy appropriate current conservation relations. We will be assuming that losses are concentrated in only one vertex which will have a potential strength with a negative imaginary part −iλ 0 . The system is turned to a scattering set-up once leads are attached at some of the vertices. In this case the system is described by a scattering matrix S(k, λ 0 ) where k is the wavenumber of an incident wave (for further details see [16,27]). The zeros of the S(k, λ 0 ) matrix are then evaluated numerically. Our numerical data for P(y) in two typical cases for M = 3 of a GUE tetrahedron graph is shown in the inset of Fig. 3a,b. The value of the scattering parameter g (associated with the real part of λ) has been extracted in two different ways from our data for graphs. The first approach invokes Eq. (12) in the absence of any loss. In this case the last term is zero and from the numerical evaluation of y we extract the value of g describing the coupling strength between the leads and the graph. The second approach utilizes the relation | S | 2 = g−1 g+1 which by using the expression (1) can be easily shown to be insensitive to any degree of loss incorporated in the network. Both approaches gave us the same value of g ≈ 1.65. Then the absorption coupling g 0 (describing the loss-strength λ 0 ) has been extracted using Eq. (12) -this time in the presence of absorption. In Fig. 3a we report the case with g 0 ≈ 1.07 corresponding to (almost) perfectly tuned absorber while in Fig. 3b we report the data associated with g 0 ≈ 2.86. From the figures it is obvious that in the former case we have a proliferation of negative zeroes as predicted by Eqs. (13). The effect is amplified even more in the case of chaotic graphs where the numerical data for y < 0 are consistently above the RMT prediction Eq. (11), see Fig. 3a. We speculate that the origin of this discrepancy is associated with the absence of a statistical gap in the resonance width distribution due to the presence of scar states, see Ref. [27,28]. A summary of our results forΠ (β=2) and various g 0 values are shown in Fig. 3c. Conclusions. We have investigated the statistics of complex zeroes of a scattering matrix describing a chaotic cavity with a single point-like lossy defect. Our approach assumes that the cavity is coupled with M leads and it is described by a RMT model. Using non-perturbative calculations we were able to evaluate the density of imaginary parts of S−matrix zeroes. This allowed us to calculate the probability of such a zero to move to the negative part of the complex energy plane as the loss-strength at the defect increases. We have tested our predictions both against direct numerical calculations within the RMT model and with a dynamical system described by a complex networks of coupled waveguides. Our results are "universal" and provide new insights into the problem of chaotic coherent perfect absorbers. Specifically they are directly addressing the question of the optimal amount of loss needed to realize a chaotic coherent perfect absorber. It will be interesting to extend these studies beyond universality and identify dynamical effects (like scars etc) that can further enhance the efficiency of such chaotic cavities to act as CPAs. Finally, we note that perhaps the easiest way to experimentally test our predictions for S−matrix complex zeroes is to use that in the framework of RMT such zeroes are statistically identical to zeroes of the diagonal (M − 1) × (M − 1) sub-block of the M × M unitary scattering matrix without absorbers. Thus we expect our results for M = 1 should describe complex zeroes of the S 11 element of the twoantenna device, which may be accessible by the same harmonic inversion technique used earlier in [38] for extracting the resonance poles.