{'fraction_non_alphanumeric': 0.03293655872690814, 'fraction_numerical': 0.01296983659712445, 'mean_word_length': 4.6332612352799805, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Implementing Non-Universal Features with a Random Matrix Theory Approach: Application to Space-to-Configuration Multiplexing

We consider the efficiency of multiplexing spatially encoded information across random configurations of a metasurface-programmable chaotic cavity in the microwave domain. The distribution of the effective rank of the channel matrix is studied to quantify the channel diversity and to assess a specific system's performance. System-specific features such as unstirred field components give rise to nontrivial inter-channel correlations and need to be properly accounted for in modelling based on random matrix theory. To address this challenge, we propose a two-step hybrid approach. Based on an ensemble of experimentally measured scattering matrices for different random metasurface configurations, we first learn a system-specific pair of coupling matrix and unstirred contribution to the Hamiltonian, and then add an appropriately weighted stirred contribution. We verify that our method is capable of reproducing the experimentally found distribution of the effective rank with good accuracy. The approach can also be applied to other wave phenomena in complex media.

We consider the efficiency of multiplexing spatially encoded information across random configurations of a metasurface-programmable chaotic cavity in the microwave domain. The distribution of the effective rank of the channel matrix is studied to quantify the channel diversity and to assess a specific systems performance. System-specific features such as unstirred field components give rise to nontrivial inter-channel correlations and need to be properly accounted for in modelling based on random matrix theory. To address this challenge, we propose a two-step hybrid approach. Based on an ensemble of experimentally measured scattering matrices for different random metasurface configurations, we first learn a system-specific pair of coupling matrix and unstirred contribution to the Hamiltonian, and then add an appropriately weighted stirred contribution. We verify that our method is capable of reproducing the experimentally found distribution of the effective rank with good accuracy. The approach can also be applied to other wave phenomena in complex media.
The phenomenon of wave multiplexing underpins a wide range of applications in wave engineering [1]. Sometimes, multiplexing devices are carefully engineered [2][3][4]; often, however, multiplexing is a result of the complete scrambling of wave fronts propagating through a complex medium (e.g., multiply scattering layer, chaotic cavity, multi-mode fiber). The most common example in the microwave domain is spatial multiplexing in multi-channel wireless communication due to an irregular propagation environment [5,6]. More recently, in order to circumvent the need for costly coherent receiver networks, the idea of multiplexing spatial information across different frequencies [7][8][9][10] or configurations [11][12][13][14][15] of a complex medium gained traction in electromagnetic imaging and sensing.
Common to all multiplexing schemes is the description with the matrix formalism Y = HX, linking the input vector X to the output vector Y via the channel matrix H. The amount of correlations between different channels determines the quality of the multiplexing: the higher the correlations are, the more redundant information is acquired. Incidentally, this insight recently motivated efforts to tweak the disorder of complex media in order to reduce channel correlations via "disorder-engineering", both in space-to-space [16] and space-to-configuration [15] multiplexing. A convenient metric to quantify channel correlations is the effective rank, defined as R eff (H) = exp (− n i=1σ i ln(σ i )), wherẽ are the singular values of H and n is the smaller one of the two dimensions of H [17]. Note that only perfectly orthogonal channels yield R eff (H) = n.
The intuition used by the wave engineering community to interpret multiplexing phenomena and to conceive applications building upon them is largely based on the assumption that H resembles a random matrix with i.i.d. Gaussian entries. Consider for concreteness space-to-configuration multiplexing with a reconfigurable chaotic cavity, as depicted in Fig. 1 and discussed in Refs. [11][12][13][14][15]. Ideally, such a system with considerable losses is indeed characterized by a perfectly stirred field distribution following the Rayleigh model [18,19]. In practice, however, a substantial unstirred field component persists, resulting in additional channel correlations and a lower-than-expected R eff [15]. As a result, the wide-spread notion of degrees of freedom as corresponding to the dimensions of H, rather than being related to its singular value spectrum [15,20,21], needs to be revisited [22,23].
The unstirred field component also presents a significant challenge from a modelling point of view. Over the past decades, random matrix theory (RMT) has seen a large success in predicting universal features of a wide range of wave-chaotic systems [24][25][26][27]. In other words, non-universal features are usually removed before comparing experimental data to RMT predictions [28][29][30]. The crux of the space-to-configuration multiplexing system we consider lies, however, precisely in the nonuniversal features which cause additional correlations between channels. Our goal here is to model the statistical properties of our specific system's scattering matrix and to thereby reproduce the experimentally observed distribution of R eff .
System-specific features arise from reflections at the ports due to impedance mismatches [31][32][33], from rays connecting two ports without ergodically exploring the cavity [34,35] and from ergodic rays that are not affected by the stirring mechanism. Refs. [36][37][38] report efforts to capture system-specific properties via semiclassical calculations of the average impedance matrix in terms of ray trajectories between ports. Such an analysis is not feasible for a geometrically complex three-dimensional cavity like the one we consider, see Fig. 1. An approach recently proposed in Refs. [39,40] addresses yet another systemspecific mechanism: the presence of an established direct transmission mediated by a resonant mode coupled to the (isolated) complex environment. This approach is not applicable in our case, where scattering has a multimode nature with both modes and channels being cou-pled mainly to the same chaotic environment.
To overcome the above limitations, we propose a twostep approach to incorporate non-universal features in an RMT framework relying on an effective non-Hermitian Hamiltonian to describe the open wave system. First, we learn a system-specific coupling matrix. Second, we determine an appropriate weighting between the unstirred (deterministic) and stirred (chaotic) components of the Hamiltonian by employing a variant of parametric motion modelling [41]. Rather than relying on knowing exact geometrical details of our system, we extract the scattering distributions of our system from experimental measurements of the scattering matrix for an ensemble of configurations. We find a good agreement between the distributions of R eff observed in our experiments and RMT simulations.
Our experiment, depicted in Fig. 1, consists of an irregularly shaped electrically large metallic cavity whose walls are partially covered with a reconfigurable reflectarray metasurface [15]. At the working frequency of 5.6 GHz, each of its 304 elements can be configured individually via the bias voltage of a PIN diode to mimic Dirichlet or Neumann boundary conditions [42]. By choosing different random configurations, the field inside the chaotic cavity can therefore be stirred in an allelectronic manner [11,[13][14][15]43]. We modulate in situ eight electromagnetic signals in amplitude and phase and inject them into the cavity via eight randomly located antennas. These eight pieces of information are multiplexed across P configurations of the system by probing the field inside the cavity with a single antenna for P random metasurface configurations. In total, M = 9 antennas are thus connected to the cavity. For concreteness, we consider P = 8 in this work. As illustrated in Fig. 1, each row of the resulting 8 × 8 channel matrix H is part of a different system's 9 × 9 scattering matrix S. In our experiment, considering the ninth port as receiver, based on 500 realizations we find the ensemble-averaged R eff = 4.11 ± 0.20 which is well below the value of 6.48 ± 0.23 expected for an 8 × 8 matrix whose random complex entries are i.i.d. Gaussian distributed [15].
Following the RMT scattering approach [24,[44][45][46], we model our system as consisting of N levels (modes) which are coupled to the environment via M scattering channels. The coupling is described by Hamiltonian H of the closed system and the coupling matrix V as  Ensemble (GOE) and the entries of V are mutually independent zero-mean real Gaussian random variables. This ensures that S RM T is a symmetric matrix. We evaluate S RM T at E = 0 and follow the normalization conventions used in Ref. [47].
To account for system-specific non-universal features, we alter two details in the above-described conventional RMT approach. First, we compose the Hamiltonian of an unstirred component H 0 and a stirred component H s , weighted by a parameter λ: Both H 0 and H s are drawn from the GOE, but H 0 is kept fixed for different system realizations. The cosine and sine terms ensure that the Hamiltonian's eigenvalue probability density function (PDF) is not altered through our modification (see Sec. 11.10 in Ref. [41]). Second, while H 0 is chosen randomly, we optimize V and λ such that the statistical properties of the resulting ensemble of S RM T match as closely as possible those of the experimental ensemble of S EXP . We assume that the nondiagonal nature of the ensemble-averaged scattering matrix S , i.e. the unstirred field components, are related to correlations between V and H 0 in light of Eq. (1); in principle, one could thus also keep V fixed and optimize H 0 [48]. Our two-step optimization is based on the intuition that V and λ determine center and diameter, respec-tively, of the clouds of S i,j in the complex Argand diagram. We hypothesize that using H = H 0 or H = cos(λ)H 0 + sin(λ)H s in Eq. (1) will approximately yield the same S RM T -an intuitive assumption that we observe to be heuristically verified in Fig. 2(a). Therefore, using H = H 0 (and hence independent of λ), we can first optimize V to ensure that S RM T Once V is fixed, we identify a suitable value of λ to ensure that the standard deviations (SDs) of the S i,j , a measure of the cloud sizes in the Argand diagram, match the experimentally observed ones.
We choose N = 250 which is sufficiently large to ensure the local character of fluctuations but at the same time not too large in terms of the computational cost of optimizing V . Assuming that global absorption effects dominate, we extract Γ a = 7.4 (in units of the mean level spacing) from the average decay rate of inverse Fourier transforms of experimentally measured transmission spectra. The optimization of V is based on sequential least squares programming as originally introduced in Ref. [49] and achieves equality of S RM T 0 and S EXP within large precision, as evidenced in Fig. 2(a). The conventional RMT approach would assume that S is diagonal and norm For our optimized coupling matrix, this is no longer the case in order to account for the abovementioned correlations between V and H 0 . Moreover, we observe that for different optimization runs (for the same objective ensemble of S EXP but using a different H 0 ), completely different norm(V i ) are obtained. This is a typical observation in inverse problems, where usually a large number of different configurations (local optima) yields optimization outcomes of comparable quality.
Having optimized V , we now identify the optimal value of λ based on the SDs of S i,j across 500 realizations. For SD[S EXP ], shown in Fig. 2(c), we notice that the average of the diagonal SDs (0.0673) is roughly double that of the off-diagonal SDs (0.0302). Moreover, the distribution of off-diagonal SDs shown in Fig. 2(c) is narrowly peaked close to the average value but a weaker second peak at roughly double the value indicates that for a few entries of S the SD is much stronger than on average. The fluctuations of the diagonal SDs for the experimental data are much larger than for the RMT data. To optimize λ, we search for the value that yields the lowest average error of the SDs of S i,j , SD . As shown in Fig. 2(b), λ ≈ 0.05 is optimal if the error of the SD is averaged over all entries of S or only over the off-diagonal entries of S. For the diagonal entries of S, however, λ ≈ 0.09 would be ideal. We use the former since our channel matrix is exclusively built from off-diagonal entries of S.
Using the optimized V and λ in combination with an ensemble of 500 realizations of H s drawn from the GOE, we obtain an ensemble of 500 realizations of S RM T via Eq. (1). As seen in Fig. 2(a), its average value S RM T is still extremely similar to S EXP , heuristically confirming our hypothesis that S RM T is approximately inde- pendent of λ. The SDs obtained for the RMT ensemble, shown in Fig. 2(c), nicely match the experimental ones for most off-diagonal entries. We notice that the RMT SDs are very uniform without any outliers, as evidenced by the narrow single-peaked PDF of the off-diagonal SDs. Consequently, our RMT model essentially appears to assume that all off-diagonal entries have the same SD of 0.03 and that all diagonal entries have a SD of 0.04. Different optimization runs starting with a different random H 0 yield very similar results.
Given our system-specific RMT model, we proceed with constructing the space-to-configuration multiplexing channel matrix H by picking eight random realizations out of the 500 available ones. We repeat this 10 5 times to compute the PDF of R eff . We perform this analysis for all nine possible choices of receiving antenna and for five optimization runs starting with a different random H 0 . The results are summarized in Fig. 3. To evidence how closely our proposed RMT model matches the experimental distribution in contrast to conventional RMT, we include the latter as benchmark in Fig. 3, too. Specifically, for each of the five runs with conventional RMT, we generate 500 scattering matrices by drawing 500 Hamiltonians from the GOE but keeping the random coupling vector V fixed for all 500 realizations. We also plot the PDF of R eff obtained for an 8 × 8 matrix for which the real and imaginary components of the entries are simply drawn from a normal distribution. Fig. 3 that the PDFs of R eff for conventional RMT and for such a Gaussian matrix are identical and independent of the choice of receiving port (i.e. of antenna coupling etc.). They clearly fail completely to predict the experimentally observed PDF, which motivates this paper. Moreover, it can be seen that the experimentally observed PDF is very well fitted with a Gaussian normal distribution. The PDFs resulting from five runs with our proposed RMT model differ slightly from run to run. Overall, the agreement of our proposed RMT model's PDF with the experimentally observed one is occasionally perfect (receiving port #7) but always decent such that our proposed model serves its purpose. Fig. 2(c) already hints at the reason why the agreement with the experiments is not perfect: the centered distributions of the S i,j are not identical for all i, j. To shed some light on this limitation of our model, we plot in Fig. 4 the distribution of selected S i,j in the complex Argand diagram as well as PDFs of real and imaginary component, using the data from the optimization run already considered in Fig. 2. The first example, S 4,7 , is representative of the behavior seen for most off-diagonal entries of S: the cloud in the complex plane is circular and the PDFs from our proposed RMT model match the experimentally obtained ones. One exception, illustrated by S 2,8 , is due to clouds with a diameter that is significantly larger than the average. Another exception, namely S 1,4 , is a case in which the cloud is not circular but elliptical. In both cases, our model fails to reproduce the PDFs correctly. For the off-diagonal entries of S, as already obvious in Fig. 2(c), our model's prediction is not very good. In general the diameter of the cloud is much larger in the experimental data, e.g. for S 1,1 , and in some cases the cloud is not circular, as for S 8,8 .

It is evident in
The S i,j clouds thus occasionally deviate from the expected behavior, either having a different radius or being deformed. Deformed (non-circular) clouds of S i,j were not seen in Ref. [43], possibly by chance because only a single S i,j was considered. Such exceptions imply that simply using centered quantities S ij − S ij does not always guarantee Rayleigh statistics in a tunablemetasurface-stirred chaotic cavity. Indeed, it was recently shown that in the presence of a deterministic scattering component the transmission amplitude and phase develop nontrivial statistical correlations even at strong absorption [50]. Such correlations may also impact techniques using such centered quantities, for instance, for antenna characterization in reverberation chambers [51].
Finally, we mention an alternative parametric model in which H = H 0 + λCC † , where C is a N × N pix matrix whose entries are zero-mean real Gaussian random variables and N pix is the (effective) number of metasurface pixels. This model could be motivated by a treatment of each metasurface element as an active scattering center or "channel" in one of two possible distinct reflecting states [52], which can be described by a low-rank perturbation of that type [53]. The weighted sum of a Gaussian Wigner and a Wishart matrix [54] brings about complications related to the distortion of the Hamiltonian's eigenvalue PDF. Our findings with this model, summarized in Ref. [55], show that overall it is not as compatible with our optimization protocol as is the model that we pre-sented above. In particular, we saw that using H = H 0 or H = H 0 + λCC † does not yield the same S (within reasonable precision), a crucial assumption for the applicability of our two-step optimization procedure.
To summarize, we have introduced a modified RMT framework to capture system-specific features based on the measurements of an ensemble of realizations of the experimental scattering matrix. Our approach, first, optimizes the coupling matrix to ensure the correct S is obtained and, second, adds an appropriately weighted stirred contribution to the Hamiltonian to adjust the fluctuations of the scattering parameters. We found a good agreement with the experimentally obtained distribution of R eff that characterizes our space-to-configuration multiplexing system. Our modified RMT scheme can also be applied to other scenarios where an ensemble of realizations of a chaotic system contains non-universal features, for instance, in the context of antenna characterization in stirred reverberation chambers [51,56] For the interested reader, here we provide additional details on the use of an alternative parametric model with H = H 0 + λCC † . This model could be motivated by a treatment of each metasurface element as an active scattering center or "channel" in one of two possible distinct reflecting states [1], which can be described by a low-rank perturbation of that type [2]. However, as summarized in the main text, this model has two important drawbacks compared to the H = cos(λ)H 0 + sin(λ)H s model used in the main text: (i) it is significantly more complicated to evaluate, and (ii) its performance in terms of reproducing the distributions of the experimentally measured scattering parameters and effective ranks is inferior. Nonetheless, for completeness, we summarize the results that we obtained with this alternative model in this Supplemental Material.

A. Overview of alternative parametric model
The alternative parametric model that we investigate in this Supplemental Material is where I α denotes the α × α identity matrix, H 0 is a N × N matrix drawn from the GOE, Γ a represents the average mode width, V is the N × M coupling matrix, and C is a N × N pix matrix. C is generated in the same way that the coupling matrix V of a traditional RMT model would be generated for a system with N pix ports with all κ i = 1. As in the main text, we choose N = 250 in sight of the computational cost of optimizing V . N pix is the (effective) number of independent metasurface pixels and we define the variable β = N pix /N .
The variable λ weights stirred and unstirred component of the Hamiltonian.
B. Impact of λ and β on the eigenvalue PDF of the Hamiltonian Since the diagonal entries of C are real and positive, they alter the eigenvalue PDF of the Hamiltonian H = H 0 + λCC † . We begin by plotting the eigenvalue PDF for a range of different values of λ for four values of β in Fig. S1. The results suggest that the position of the eigenvalue PDF's peak does not appear to be impacted by the addition of the stirred Hamiltonian. However, the maximum value of the eigenvalue PDF is reduced.
We continue to work at E = 0 as in the main text. To counteract the lower value of the PDF of eigenvalues at E = 0, we consider replacing Γ a by Γ a = αΓ a , where α is the ratio of the slope of the curve of sorted eigenvalues around zero for H and H 0 . However, since we find that accounting for the change in slope or not does not have an appreciable impact on the final results, we neglect the correction in the following. We choose N pix = 50 in the following which ensures N pix < N to preserve the GOE shape of the level density near the origin and N pix 1 to reflect the fact that we have a multitude of pixels rather than a single or a few pixels. Since the optimization of V is assumed to only dependent on H 0 in our two-step optimization procedure, it is the same as for the parametric model discussed in the main text; therefore, we only discuss the optimization of λ for the alternative parametric model here.
We present in Fig. S2 the counterpart to Fig. 2 from the main text. The optimal value of λ ≈ 0.06 is comparable to that found with the parametric model in the main text. However, S RM T does not resemble S EXP as much as for the parametric model presented in the main text. In Fig. S2(a), the differences of the magnitude can be spotted by the naked eye and the magnitude of the difference between S RM T and S EXP is an order of magnitude larger. The alternative model is thus not as suitable for our two-step optimization procedure  Experiment and simulation are presented in blue and red, respectively.

E. S i,j distribution
For completeness, we finally provide the counterpart to Fig. 4 from the main text in Fig. S4. We note that even for scattering matrix entries like S 1,4 for which the SD is correctly captured and which are correctly reproduced by the model discussed in the main text, the center of the cloud generated by the alternative model is dislocated. This is essentially due to the inferior match in terms of S already seen in Fig. S2(a).