Wavefront Shaping with a Tunable Metasurface: Creating Coldspots and Coherent Perfect Absorption at Arbitrary Frequencies

Modern electronic systems operate in complex electromagnetic environments and must handle noise and unwanted coupling. The capability to isolate or reject unwanted signals for mitigating vulnerabilities is critical in any practical application. In this work, we describe the use of a binary programmable metasurface to (i) control the spatial degrees of freedom for waves propagating inside an electromagnetic cavity and demonstrate the ability to create nulls in the transmission coefficient between selected ports; and (ii) create the conditions for coherent perfect absorption. Both objectives are performed at arbitrary frequencies. In the first case a novel and effective optimization algorithm is presented that selectively generates coldspots over a single frequency band or simultaneously over multiple frequency bands. We show that this algorithm is successful with multiple input port configurations and varying optimization bandwidths. In the second case we establish how this technique can be used to establish a multi-port coherent perfect absorption state for the cavity.


I. INTRODUCTION
Extreme electromagnetic environments are prevalent in much of our day to day lives. While often unnoticed, this puts significant stress on the design of electronic systems that are expected to work under all conditions. Extraordinary care is taken to counter adverse effects whenever possible, but the operating environment is rarely known ahead of time. In addition, Electromagnetic Interference (EMI) takes the form of unwanted coupling between components. Some platforms, such as aircraft and spacecraft, can experience devastating consequences, resulting in severe mission degradation or even casualties [1]. Wave fields in electrically large enclosed areas such as passenger compartments in transportation systems show extreme sensitivity to frequency and geometric details even though these enclosures, termed chaotic cavities, are not intended to be reverberant. In such cavities the electromagnetic wave fields have specific statistical properties that depend upon a limited number of parameters [2]. Among these is the fact that the wave field is statistically equivalent to a random superposition of plane waves [3]. As such, we can leverage analytical tools from the active research area of quantum chaos [4,5] in the more generalized framework of wave chaos [6].
Our goal is to adaptively and intelligently control fields inside complex electromagnetic environments. We show that this can be achieved through programmable metasurfaces, which increase the available Degrees of Freedom (DOF) by manipulating boundary conditions [7][8][9]. In recent years, these devices have enabled novel concepts in * Benjamin.Frazier@jhuapl.edu a wide variety of applications throughout the microwave and optical domains [10][11][12][13][14][15][16]. The additional DOF in turn enhance diversity in space and time [17][18][19][20], and allow intricate control of the underlying scattering system; combining a programmable metasurface with a cavity unlocks applications not possible with a fixed system [21][22][23] and encourages research in new and under-explored domains [24].
One such unexplored area is Coherent Perfection Absorption (CPA), where coherent excitation of a lossy system can result in complete absorption of all incident waves [25,26]. It has applications in highly efficient notch filtering, energy conversion, and even detection; since the CPA state is extremely sensitive to parameter variation, it can be used to identify small changes in a complex scattering system [27]. The ability of a metasurface to manipulate additional DOF presents a novel capability for realizing CPA states.
In this article we describe the use of a binary programmable metasurface to create microwave coldspots at arbitrary frequencies, and to realize CPA states, both within the 1 GHz band of operation of the metasurface. The conceptual overview is shown in Figure 1 where the metasurface is installed in a complex reverberating cavity and controlled in a closed loop manner. Input directional diversity is introduced by simultaneously driving multiple ports with arbitrary relative phase shifts. An iterative optimization algorithm is used to generate coldspots at the output port, or to drive candidate scattering matrix eigenvalues towards the origin to achieve CPA. Figure 1. Conceptual overview of the metasurface enabled cavity as a closed loop system. The cavity S-parameters (scattering parameters) are measured with a network analyzer and passed to a controller that updates the metasurface elements with a new set of commands. The controller can generate coldspots at port 2 at an arbitrary set of frequencies, or drive candidate S-matrix eigenvalues towards the origin, and includes a stochastic iterative optimization algorithm. The three ports allow additional angular and spatial diversity to be added at the inputs. The inset shows a closeup view of one of the metasurface unit cells.

II. CAVITY CONFIGURATION
The metasurface used for this work is a reflectarray fabricated by the Johns Hopkins University Applied Physics Laboratory (JHU/APL) that is designed to operate from 3-3.75 GHz. It has a lattice of 10 × 24 squares occupied by unit cells with size < λ/6 [28], where λ is the wavelength. These 240 unit cells can be independently set to one of two states, which approximate electric or magnetic boundary conditions and provide a relative 180 • phase shift for waves reflected by the element. This results in the local surface impedance of the array varying from cell to cell and state to state. The array surface thus has 2 240 independent states, each of which reflects waves in a uniquely different set of directions. We refer to the settings of the 240 elements as a command.
The array was installed in a 0.76 m 3 cavity where it covers ∼ 1.5% of the total interior surface area. The cavity has 3 ports with one acting as a target for scoring and two used for signal injection; the input ports can be driven either individually or collectively with a relative phase shift. When all three ports are used, the underlying scattering system is represented by a 3 × 3 scattering matrix. The cavity has both low and high loss configurations to test how the behavior varies with the typical quality factor, Q, of the modes; here we consider the high Q case. Introducing the metasurface to the cavity reduced the average Q in the frequency band of operation of the metasurface by a factor of ∼2. However, once the metasurface was installed, the average Q was found to be independent of the number of active or inactive elements on the surface. The quality factor was determined to be roughly 5.5 × 10 3 , by measuring the power decay time (250 ns with the metasurface installed). Details of the metasurface, cavity construction, experimental setup, and impact of the metasurface on losses are provided in the supplemental material.
A measure of the loss in the cavity is the Q-width of a mode normalized to the mode spacing, α = f /(2∆f Q) = 3 for this cavity. For our cavity, the mean mode spacing is roughly 115 kHz at a 3.5 GHz center frequency. As discussed in the supplemental material, the mean spacing between nulls in the transmission coefficient, |S 21 |, was found to be ∼2 MHz and the average width of the nulls was found to be ∼200 kHz. This indicates a transmission coefficient null contains about 2 modes. Alternatively, it corresponds to a path difference of 750 m between two interfering signals.
A complex scattering system such as our cavity exhibits both universal fluctuations, which can be described by Random Matrix Theory (RMT) [4], and deterministic behavior arising from the system specific configuration of the ports and short orbits between the ports [30][31][32]. Due to the small relative size of the metasurface, the chaotic ray paths with many multipath bounces will experience the strongest influence. Since minimizing the power received at a port is accomplished by creating destructive interference of the ray paths, the relationship between commands and responses is quite complicated. This leads to utilizing stochastic iterative approaches, or machine learning, in place of linear deterministic methods for control. Here we consider iterative processes. A metasurface covering a larger fraction of the interior surface would likely produce stronger results [33] and allow us to use a transmission matrix based approach to determine the optimal metasurface commands [34][35][36]. For this reason most prior research utilizes metasurfaces that cover a significant portion of a wall (or multiple walls). However, using a relatively small metasurface coverage is better suited for real world applications where it is not practical to build or use a larger device.
A key step in evaluating system performance is to determine the range of possible responses of the scattering properties of the system so as to ensure that we have a sufficiently diverse command set. Unfortunately, with 2 240 possible commands (approximately 1.8 × 10 72 ), it is not feasible to test every one and we need to find a reduced number that produces the full range of outcomes. As discussed in the supplemental material, deterministic decomposition of commands into orthogonal basis functions, such as Hadamard bases [37,38], generated a very narrow range of system scattering responses. Diversity in the responses requires a distribution of commands with a variety of spatial frequencies, ratios of active to inactive elements, and localized groupings of active elements. Doubly random methods or compound distributions, such as a biased coin toss, or power law spectral density with the bias, or power exponent itself a random draw, were found to yield the widest range of responses. Details of our novel stochastic algorithm are discussed in the next section and the supplemental material.

III. GENERATING COLDSPOTS
Our goal is to program the metasurface to minimize the transmission between two ports in a complex scattering system at an arbitrary frequency. Cases are scored by evaluating the difference in average power, ∆P 2 , in a specified frequency range at a given center frequency between the initial inactive (all 0s) state and the current state of the metasurface. To maximize this difference we take a directed random walk approach in which at each step a number of array element states are toggled (changed), ∆P 2 is evaluated, and the new state is accepted or rejected based on whether or not it decreases ∆P 2 . As discussed in the previous section, we need to have a mix of large and small spatial groupings of ele-ments and a varied number of active elements to ensure a diverse set of responses. To meet this requirement, our iterative algorithm operates in 2 distinct phases: multiple element toggling and individual element toggling.
In the multiple element toggling phase, we select M elements at random as a trial and toggle their state (0 → 1 and 1 → 0). If ∆P 2 is decreased, the trial set of commands becomes the new reference set and we repeat the process selecting another M elements at random and toggling their state. When T consecutive trials have been made without improving ∆P 2 we claim convergence and move to the next value of M . In a typical experimental run, M = [120, 48,24,12,6], and T = 30. After all values of M have been exhausted, we move to the individual element phase.
The individual element toggling phase has 3 cases associated with each trial. We select a single element at random and toggle it and then, in an adaptation of the neighbor toggling method of Ref. [39], we toggle the 4 nearest neighbors and the 4 diagonal neighbors. ∆P 2 is evaluated for each of these cases and the algorithm continues as in the 1st phase until T consecutive trials are performed without improving ∆P 2 .
The multiple element toggling phase tends to result in a local minimum which is difficult to escape when toggling only a single element. Adding neighbor toggling significantly improves the performance, as it provides larger localized changes in the command set and allows us to escape the local minimum. Even with the neighbor toggling, however, our stochastic approach does not guarantee that a global minimum is found. Increasing the convergence criteria, T , can increase the probability of finding the global minimum, but comes with the cost of increased time. The absolute minimum is not necessarily required, and our stochastic algorithm is able to provide substantially deep nulls at arbitrary frequencies in a reasonable amount of time. A typical experimental run will provide ∼350 trials, ∼25 iteration updates, and take ∼1.5 hours. The experimental setup is not optimized for run time. We use an ethernet connection to transfer 32001 frequency samples over the full 1 GHz band for each of the 4 complex S-parameter measurements using 64-bit precision. With the frequency values themselves included, this means 2.3 MB of data are transferred for each trial. In addition, the commands and measured |S 21 | are plotted at each trial for operator feedback, resulting in ∼15 seconds per trial. Disabling plotting and capturing only the processed frequency band could potentially reduce the time to 1-2 seconds per trial. Figure 2 shows the results obtained when minimizing the average power at the output port and compares the results of many different experiments and configurations. All the cases are scored by the change in average power, ∆P 2 , between the initial inactive (all 0s) state and the final state. The optimization algorithm was performed with ∆P 2 evaluated over a single frequency band as well as simultaneously over multiple separated frequency bands. As discussed previously, the widths of the nulls were observed to be ∼200 kHz. The initial bandwidth was selected to be 500 kHz in order to ensure that ∆P 2 was evaluated over an entire null. In addition, the cavity configuration was switched between driving a single input port and driving two input ports simultaneously with varying relative phase shifts. The achieved suppression ranges from 4-40 dB with most cases providing > 10 dB. The lower values of ∆P 2 arise in the following cases: working near the edges of the metasurface operational window, evaluating ∆P 2 over a large bandwidth, or evaluating ∆P 2 over multiple separated bands. This is not surprising as more bandwidth results in more features in the region where ∆P 2 is evaluated, which then means more degrees of freedom are required to be manipulated for destructive interference. The metasurface provides some benefit outside of the 3-3.75 GHz design window; the reflection phase change of the pixels is limited near the edges of the operational bandwidth, so performance is expected to be reduced under those conditions. More details on specific cases are provided in the supplemental material.
Since ∆P 2 is inherently a relative measurement, there is an implicit dependence on the initial state. Using the inactive (all 0s) state as the reference ensures the metasurface is always initialized with the same command even though the specific value is dependent on the selected frequency window. Starting with a condition where there was already a deep null would result in limited improvement; the average power in that case would already be quite low and there would not be a need for further reduction. Starting with a condition where there is a transmission peak however, would result in significant reduction. When using a single frequency band metric, we were able to drive deep nulls in each of the windows that were tested, as can be seen by the circles in Fig. 2a and the power at port 2 in panels b) through e). Panels b) and d) show moderate cases where there is not a clear peak in the initial P 2 measurement, while panels c) and e) show cases with a clear peak in the initial P 2 measurement and demonstrate significant improvement. This highlights the dependence of ∆P 2 on the initial state.
Deep transmission nulls were also observed when driving two input ports simultaneously with varying relative phase shifts, as shown by the diamonds in Fig. 2a. This indicates our approach is self-adaptive and can compensate for multiple input signals as well as signals coming in from different directions. With dual frequency bands however, we were generally unable to drive deep nulls in both bands simultaneously, which can be seen by the hexagrams in Fig. 2a and the power at port 2 over the frequency band in panels f) and g). This is because the metasurface frequency response in separated bands is correlated, as the metasurface induces wide bandwidth effects on the scattering properties of the enclosure. As shown in the supplemental material, different choices of metrics produce different out of band behavior.

IV. GENERATING COHERENT PERFECT ABSORPTION
Coherent Perfect Absorption (CPA) is a situation in which all energy injected into a system is absorbed, no matter how small the losses are in the system. Creating CPA requires coherent excitation of all the ports in an eigenvector whose corresponding S-matrix eigenvalue is zero. Operationally, the first step in establishing CPA is to find an eigenvalue of the scattering matrix that is close to zero. For example, a 2 x 2 scattering matrix will have a pair of eigenvalues at each frequency. However, realizing CPA only requires driving 1 eigenvalue to zero, as the other eigenvalue corresponds to the anti-CPA state [27]. For the following discussion and experimental results, we only consider the smallest eigenvalue of each pair.
CPA has typically been investigated in simple, regular scattering scenarios and cavities but recently it has been demonstrated in more complex systems, specifically in the realm of wave chaos, and graphs [40][41][42][43]. These works analytically demonstrate the use of RMT to explore CPA states with semiclassical tools without relying on the limit of weak coupling. CPA states have also been experimentally investigated in multiple scattering environments [26], and in graphs that break time-reversal invariance [27]. The use of enhanced spatio-temporal diversity from a metasurface for realization of CPA has not yet been explored.
Recent research however has investigated the use of metasurfaces for Perfect Absorption (PA) inside a cavity and demonstrated a secure communication system as an application [44]. PA is a complementary idea to CPA for a single port system that relies only on the reflection coefficient, S 11 [45]. Coherent excitation of a single port with complete absorption has been demonstrated to enhance wireless power transfer [46]. Our work extends this to coherent operation with the full scattering matrix for a 2-port system, and can be generalized to an arbitrary number of ports.
Realizing a true absolute zero of the S-matrix eigenvalues is generally difficult, because the eigenvalues are complex numbers. Thus two parameters must be varied independently to drive an eigenvalue to zero. Further, a CPA state is highly dependent on the structure of the underlying scattering system. This is best understood in the framework of the Random Coupling Model (RCM) [6,47]. The eigenvalues accessible by means of the programmable metasurface tend to cluster around values determined by the coupling properties of the ports, which are characterized by the radiation S-matrix, S rad . We define S rad as the S-matrix corresponding to the freespace radiation condition with the cavity walls taken out to infinity such that no waves come back to the ports [48]. S rad can be determined by a number of means [49].
Here we employ the ensemble average of the time gated measured S-parameters in the cavity [29], as described in the supplemental material. Deviations of the scattering matrix from S rad have a number of causes. First there are deviations resulting from relatively direct ray paths between the ports [31]. These deviations are removed by averaging the S-matrix over a frequency window that is the reciprocal of the time of flight on the path. However, in finding the eigenvalues of the S-matrix in a narrow frequency range, these deviations are present. Second, there are deviations due the multitude of longer paths, and these are characterized statistically by RMT within the RCM. These fluctuations in S tend to be of the order of 1/(πα) 1/2 [30][31][32] where the loss parameter α = f /(2∆f Q) = 3 in the present experimental case. Finally, there are deviations dependent on the state of the metasurface. These deviations are constrained to be less than or equal to either the direct path or the statistical long path deviations.
Thus, to find a CPA state it is necessary for the ports to be sufficiently matched so that the statistical fluctuations can shift the eigenvalues to zero. If the ports are poorly matched and losses within the cavity are sufficiently high, the eigenvalues will naturally fall near values determined by the properties of the ports with statistical fluctuations around those values dictated by the amount of cavity loss. As such, it is generally not possible to realize a CPA state at arbitrary frequencies when limited to a single DOF [27]. The availability of additional DOF, such as those produced by the metasurface, allows greater control over the underlying scattering system and provides a greater likelihood of potential CPA states.
Characterization of the S-matrix eigenvalues from a distribution of 2000 command sets is presented in Figures 3 and 4. Figure 3 shows the probability distributions for all of the S-matrix eigenvalues over all frequencies and commands. Panel a) shows that the magnitude follows a Rician distribution as predicted by Ref. [50], which also tells us that the ν parameter of the Rician distribution is due to the presence of persistent short orbits [51]. Panel b) shows that the phase of the S-matrix eigenvalues is not truly uniformly distributed. The deviation of the eigenphase from uniformity indicates that the random distribution is not statistically independent and again tells us there are persistent short orbits present in the system. These short orbits are not captured explicitly in S rad , and will cause the eigenvalues of S rad to be offset from the center of the point cloud of S-matrix eigenvalues. Short orbits can be explicitly included analytically in the RCM [51,52] and do not prevent us from proceeding. Panel c) shows the CDF of the eigenvalue magnitudes and is useful in establishing thresholds for potential CPA candidates. Figure 4 shows point clouds of the S-matrix eigenvalues at selected frequencies and demonstrates that the eigenvalues can have very different behavior in how they approach the origin. The panels show the collection of eigenvalues of the 2000 S-matrices at selected frequencies along with the eigenvalues of S rad at that frequency. We can see that the eigenvalues of the distribution tend to cluster around the eigenvalues of S rad ; the offset from the center of the point cloud is due to the presence of short orbits, as discussed above. In panel a), the S-matrix eigenvalues are clustered in the upper right quadrant far from the origin and do not enter the inner rings. The S rad eigenvalue is in the upper right-hand quadrant outside of the plot area, at 0.1862+j0.2288. In panel b), the S-matrix eigenvalues are clustered in the upper half, with some getting close to the origin. In panel c), the S-matrix eigenvalues show a fairly uniform density throughout the full |λ s | < 0.15 range. In panel d), the S-matrix eigenvalues show a high density clustered around the origin. The results in these panels are from a random distribution of commands rather than a targeted search. During optimization, we will take smaller dithering steps for finer control as we approach the origin, and expect to see slightly different behavior.
The variance in eigenvalue magnitudes means we need to use a large threshold for identifying candidates because the overall global minimum S-matrix eigenvalue may not be identified as a candidate in every realization. In practice, we found that we were unable to realize CPA states when starting with a magnitude |λ s | ≥ 0.2 but were generally able to realize CPA states when starting with a magnitude |λ s | ≤ 0.15. Moving an eigenvalue far from the origin requires modifying the underlying scattering matrix more strongly than moving an eigenvalue that is already near the origin, so this behavior is ex- pected. Assessing the probability of finding a CPA state in a given frequency range a priori is difficult. The universal properties of a complex scattering system are not easily separated from the deterministic properties when working with S-parameters, as the statistics are dominated by S rad [53]. This means the existence of a CPA state is highly dependent on the coupling properties of the ports and therefore the specific antennas chosen. An analytical approach is possible through the framework of the RCM and will be left to future work.
An open question is how small do the eigenvalues need to be to realize CPA? This is dependent on the specific application and scattering system, as that determines how accurately the eigenvalues can be measured and maintained. For our experimentation, we set |λ s | ≤ 5 × 10 −3 as the upper bound and |λ s | ≤ 1 × 10 −3 as the goal for realizing CPA.
We adopt the same basic algorithm used for power minimization but initialize it differently. We apply a random set of commands to the metasurface and then select a candidate eigenvalue with a specified magnitude.  Fig. 5a shows the behavior of 4 selected cases away from the origin for |λ s | < 0.15, and Fig. 5b shows the behavior at the CPA condition for |λ s | < 5 × 10 −3 . Only 3 of the 4 selected cases reach the CPA threshold. Fig. 5c presents the collection of all 27 experiments, with the four shown in detail in panels a) and b) color coded. Each case was initialized with an eigenvalue magnitude chosen in the range 0.075 ≤ |λ s | ≤ 0.5. The case that started with |λ s | = 0.5 is enclosed by a triangle, the cases that started with |λ s | = 0.2 are enclosed by squares and the case that started with |λ s | = 0.175 is enclosed by a circle. All the rest started with |λ s | ≤ 0.15. Three cases initialized with |λ s | = 0.15 did not quite make the CPA threshold, |λ s | ≤ 5 × 10 −3 . Two cases were within a factor of 2, |λ s | ≤ 9 × 10 −3 , while the third was within ∼20%, |λ s | = 6 × 10 −3 .
Utilizing the iterative optimization algorithm to change the metasurface, we are able to drive eigenvalues towards the origin in all cases, but the algorithm stalls at different points. The closer we get to the origin, the more difficult it becomes to reduce the eigenvalue further. As with the coldspot optimization, the stochastic nature of the algorithm plays a role in where convergence is reached. The overall performance could be improved by increasing the convergence criteria or making the algorithm adaptive so that it tracks multiple candidates and switches to another candidate when the optimization stalls.
As a final step, we want to verify that the CPA state has been achieved. Because the CPA state is found by minimizing the eigenvalues of the scattering matrix, verification requires that we apply the corresponding Smatrix eigenvector. This can be done using a network analyzer with 2 independent sources and an external phase shifter [27]. After directing a particular eigenvalue to-wards the origin, the network analyzer was configured for independent source operation and the amplitude and phase were adjusted to generate the eigenvector, as described in the supplemental material. The presence of a CPA state is verified by looking at the ratio of all the power emerging from the cavity to all the power injected into the cavity, P out /P in . Sensitivity to changes in the eigenvector can be determined by making small deviations in the relative phase shift or amplitude between the two sources. Sensitivity to the eigenvalue can be determined by small changes in frequency.
A set of parameter sweeps that verify a CPA state was realized are presented in Figure 6. Fig. 6e shows the S-matrix eigenvalue magnitude trajectory during optimization prior to performing the verification sweeps. The overall experimental setup is shown in Fig. 6f, which shows that a 2-source network analyzer was configured with independent source operation and connected to the cavity with an external phase shifter on port 1. This allows us to produce the appropriate eigenvector by controlling the relative amplitude with the network analyzer and the relative phase with the phase shifter. The metric for the sweeps is the power ratio, P out /P in , of all the power emerging from the cavity to all the power injected into the cavity. At the CPA condition, all the energy should be absorbed. However, due to instrumentation limitations with the system noise floor, the smallest measurable power ratio is ∼10 −6 . Before performing the sweeps, the eigenvector was tweaked to provide the closest CPA state realization and then the parameters were varied to determine the sensitivity of the power ratio. Fig. 6a shows the results of the frequency sweep performed in a ± 10 MHz window around 3.6697 GHz, with the inset showing a closeup in a ± 200 kHz window. The width of the deep null is ∼200 kHz, which matches the null widths found during cold spot generation. Fig. 6b shows the results of the phase sweep, which was performed by adjusting the external phase shifter. Here, ∆φ represents the phase shift at port 1 away from the CPA eigenvector phase. Fig. 6c shows the results of the relative amplitude sweep. This was performed by sweeping the power injected into port 1 from -10 dBm to +10 dBm. The x-axis is then scaled to show the relative change in injected amplitude from the initial CPA state. In each of these cases, the minimium power ratio is ∼ 6 × 10 −6 and shows a steep cusp-like increase with the various parameters. Fig. 6d shows the results of the metasurface command sweep. In this case, the 240 individual metasurface elements were toggled to determine the impact of a single element on the CPA state. Several elements had negligible impact on the power ratio in comparison with the optimized value as seen in the dashed black line, but no toggles were found with clearly better performance. The elements in the center of the metasurface have a stronger impact than those at the edges of the metasurface, but the largest change from the CPA condition was observed by setting all the elements to 1s as shown in the dashed magenta line.

V. CONCLUSIONS
We have demonstrated the ability of a programmable metasurface to generate microwave coldspots in a chaotic cavity at arbitrary frequencies and showed this capability exists even when applied over multiple frequency bands simultaneously. The coldspots can be generated for different bandwidths and mulitple input port configurations that induce additional angular and spatial diversity. We have also utilized the programmable metasurface to control the eigenvalues of the scattering matrix and direct them towards the origin to realize a CPA state for the cavity. Finally, we verified the existence of a CPA state and demonstrated the sensitivity to parameter sweeps in frequency, phase, amplitude, and metasurface configura-tion. All of this is accomplished with a metasurface that covers only 1.5% of the interior surface area of the cavity and a unique and effective stochastic algorithm to find desired outcomes despite the enormous space of possible metasurface commands.
Future research directions include quantitatively analyzing CPA in the framework of the Random Coupling Model [6,30,32,47], and using deep learning to facilitate generating optimal metasurface commands to minimize power and/or realize CPA states. Of the four cases shown, only 3 were able to get inside the inner rings near the origin where |λs| < 5 × 10 −3 . c) Minimum achieved eigenvalue magnitude for each performed experiment. The circles indicate data that is shown in the upper plots and are color coded to match. The gray squares indicate an experiment that was performed but whose detailed trajectory is not shown in the upper plots. The dashed black lines indicate the cross over points of 5 × 10 −3 and 1 × 10 −3 . The enclosed squares indicate cases where the initial eigenvalue magnitude |λs| > 0.15. |λs| = 0.5 for the triangle, |λs| = 0.2 for the squares, and |λs| = 0.175 for the circle.

VI. ACKNOWLEDGEMENTS
We thank Joseph Miragliotta, David Shrekenhamer, and Robert Schmidt for providing the metasurface and supporting detailed discussions on operational use. We also thank Tsampikos Kottos for insights into CPA, and Lei Chen for instructions on performing the various CPA verification sweeps. Funding for this work was provided through AFOSR COE Grant FA9550-15-1-0171 and ONR Grant N000141912481.

VIII. COMPETING INTERESTS
The authors declare no competing interests.

IX. MATERIALS AND CORRESPONDENCE
Correspondence and requests for materials should be addressed to B. Frazier (email: benjamin.frazier@jhuapl.edu).

X. DATA AVAILABILITY
The data and analysis scripts that support results presented within this paper and other findings of this study are available from the corresponding author upon reasonable request.

I. EXPERIMENTAL SETUP
The metasurface used for the experimentation is a reflectarray fabricated by the Johns Hopkins University Applied Physics Laboratory (JHU/APL) and is shown in Figure 1. The individual LC resonator unit cells can be seen on the front side and the GaAs FET amplifiers can be seen on the back side. It is designed to operate in the S-band from 3-3.75 GHz and provides 240 binary unit cells arranged in a rectangular grid of 10 rows and 24 columns. The unit cells are electric LC resonators driven by GaAs FET amplifiers with characteristic length < λ/6 [1].
A test cavity was constructed using aluminum angle brackets for the frame and 0.019 inch thick aluminum sheets for the sides. Each side length is 3ft (0.9144 m), so the total cavity volume is 0.76 m 3 and the total surface area is 5.02 m 2 . This geometry ensures the cavity is overmoded with at least 9 wavelengths per side, but unfortunately means that the active area of the metasurface only covers a small portion (∼ 1.5%) of the total surface area. All interior joints were lined with copper tape to minimize losses and hemispherical scatterers were installed in the corners of the cavity to force an irregular shape, after which the effective volume was reduced to 0.74 m 3 . A higher loss version of the cavity was realized by distributing absorbing material in the bottom of the cavity. The metasurface was installed on a wall of the cavity as shown in Figure 2, with a 1/4 inch gap between the metasurface and the wall. This figure also shows the power and 3 USB cables running out through the top of the cavity.
An overview of the experimental setup is shown in Figure 3 and shows the cavity configuration, the locations of the 3 ports relative to the metasurface, and the overall connectivity. The cavity has 3 ports with port 2 used for scoring and ports 1 and 3 used for signal injection. The input ports can be driven either individually or collectively with a relative phase shift. When they are driven collectively, the underlying scattering matrix is 3 x 3; * Benjamin.Frazier@jhuapl.edu this extra dimensionality is hidden when using a 2-port network analyzer. To ensure that wave interaction with the metasurface dominates the results, a sheet of foil is used to block the direct line of sight path from the input ports to port 2. Ultra Wideband (UWB) antennas designed for operation from 3-10 GHz are connected to each port. The antennas connected to ports 1 and 3 are Taoglas FXUWB10 antennas and are mounted so they radiate outward into the cavity from the walls in a vertically polarized configuration. The antenna connected to port 3 is a PCB module mounted orthogonally so that it is horizontally polarized.
Finally, the metasurface is controlled through 3 USB interfaces using custom C++ software with a Python wrapper on a MacBook Pro laptop, which also controls the Agilent network analyzer through an ethernet connection. In order to prevent corruption from noise, multiple measurements are averaged.

II. CAVITY LOSSES
The cavity time constant, τ , is an intrinsic aspect of the cavity that is dependent on losses but independent of the specifics of the underlying scattering system. This means τ is not dependent on the positioning or characteristics of the ports or antennas used to obtain it [2]. τ is an important characteristic of the cavity as it is related to the quality factor, Q, through Q = ωτ . One way to estimate τ is through Power Delay Profile (PDP), which is the ensemble average of the magnitude squared of the Inverse Fourier Transform (IFT) of the transmission coefficient, PDP = |IFT {S 21 }| 2 [3,4]. Since the power in the cavity decays exponentially, we can perform a linear fit on the logarithmic magnitude and estimate τ as 4.34/ν where ν is the slope of the PDP in dB/s [5]. Figure 4 shows the time constant estimated for the cavity under various configurations. There were 3 antenna pairs used in the PDP measurement: dipoles with both horizontal and vertical polarization, loops, and Ultra Wide Band (UWB). Measurements were taken with the cavity empty before installing the metasurface and in low and high loss cases with the metasurface installed.  Each data point in Figure 4 came from a 100 MHz window centered at the corresponding frequency. Installing the metasurface in the cavity had a significant impact on the time constant, reducing it by a factor of 2. Adding absorbing material for the high loss configuration reduced the time constant by another factor of 2. Powering on and off the metasurface however, had little impact on the time constant in either configuration. The PDP was measured between ports 1 and 2 with 3 different antenna pairs: dipole (in both horizontally and vertically polar- ized orientations), loop and Ultra Wide Band (UWB). The On and Off curves indicate whether the metasurface was powered on or off.

III. EXTRACTING RADIATION S-PARAMETERS
In many cases, it is not feasible to directly measure S rad , but it can extracted from time gating the S cav measurements [5]. To ensure that the important features of S rad are maintained, the optimal gate width is determined by examining the S-parameters in the time domain and selecting the point in time where the individual ray trajectories start to diverge from the average.
The conventional time-gating approach is to convert the signal into the time domain with an Inverse Fast Fourier Transform (IFFT), multiply the signal with a rectangular gate, and then convert the signal back to the frequency domain with a Fast Fourier Transform (FFT). This approach compounds truncation effects through the IFFT/FFT sequence and induces band edge roll off effects due to the fact that we are using a finite, single-sided spectrum [6]. An alternative approach is to perform gat- ing in the frequency domain through a convolution operation and use the concept of renormalization to remove band edge artifacts. In order to optimize frequency domain gating, the gate needs to be centered at the time of the response being gated [6]. Because we are interested in gating the initial response, we will center the time window at 0 seconds; accordingly, the overall width of the gate will be double the desired end time to include both positive and negative time extents.
The gate can be designed in 2 segments: a rectangular gate in time transformed to the frequency domain, and a window to reduce side lobes and ringing artifacts due to the sharp transitions of the rectangular gate. The generalized gate function in the frequency domain, G(f ), for a rectangular time domain gate defined between times t 1 and t 2 is given by a sinc function, where we define sinc(x) = sin(πx)/(πx).
For a gate centered at t = 0 with end time t g (t 1 = −t g , t 2 = t g ), this expression is simplified.
The next step is to design a window function in the frequency domain to suppress side lobes. A common window is the Kaiser-Bessel window, which is defined by a shape parameter, β and the window length, M [7]. For the analysis described here, the Agilent PNA provided 32,001 points of data and the window was defined with β = 6.5 and M = 23, 897. We can then apply the gate by convolving the product of G(f ) and W (f ) with the measured S-parameter for a given cavity realization. Renormalization is done by dividing the gated S-parameter by the convolution of the product of G(f ) and W (f ) with a constant unit frequency response, which removes band edge roll off effects [6]. S rad is then found by taking the average of the results over the ensemble of cavity configurations as shown in Equation 3.

IV. DIVERSE CAVITY REALIZATIONS
Attempts were made to decompose the input commands into a deterministic set of orthogonal basis functions. This included driving single elements, columns of elements and even Hadamard matrices. Unfortunately, these all produced a narrow range of responses that did not cover the full extent of possibilities. A Hadamard matrix provides an orthonormal basis and decomposes sequences similarly to spatial frequencies; it is less computationally intensive than 2-D Fourier transforms and has many applications in multi-input multi-output communication theory and synthetic aperture imaging [8,9]. Figure 5 presents the resulting ensembles from driving the metasurface with Hadamard basis functions and shows that the responses are narrow and do not cover the full extent. While the spatial frequency content and grouping of elements in the command sets changed, the number of active elements did not. An unbiased random coin toss approach likewise yielded a narrow distribution as roughly 1/2 the elements were active in each command set. The ensembles do not span the range covered by the inactive (all 0s) and active (all 1s) cases and do not cover the full extent. More variation is present in the low loss case because the ray trajectories survive longer and have more opportunities to interact with one another The overall best approach to generating diverse ensembles was to use doubly random methods, in which compound probability distributions are used. This implies the statistics follow a Cox process, which is a generalization of a Poisson process with the underlying intensity or local mean itself a random process [10]. A random biased coin toss, where the bias itself was selected from a random draw for each command set worked well but only excited high spatial frequencies on the metasurface. To include lower spatial frequencies, we added a doubly random inverse power spectrum approach, where the power spectrum is just a power law on spatial frequencies with the exponent a random draw.
A small exponent will excite lower spatial frequencies, while a larger exponent will excite higher spatial frequencies. the number of active elements was allowed to change with each trial. We generally used a combination, with half the ensemble generated through an inverse power  spectrum and the other half generated through a random biased coin toss. The ensemble for a set of 4000 realizations is shown in Figure 6. The doubly random methods allow the number of active elements to change and provide a wider range of responses than the deterministic methods. This can be seen as the distribution from the biased coin toss covers the entire area between the inactive (all 0s) and active (all 1s).

V. POWER OPTIMIZATION RESULTS
The optimization algorithm was successfully performed with different cavity configurations and different performance metrics. From Figure 6, we can see that the bandwidth of the narrowest null is ∼ 500 kHz and the closest spacing between nulls is ∼ 1 MHz. This matches the observed trends over the full 1 GHz measurement window, with typical bandwidths of 0.5-1 MHz and spacings of 1-2 MHz. The metric was chosen to be the average power measured at port 2, P 2 , when driving an input from some combination of ports 1 and 3. To isolate narrow band features, our initial metric bandwidth was chosen to be 500 kHz.
For initial validation of the algorithm, we used single frequency bands. To ensure that the band selection was not biased to only bands where the metasurface performed well, the center frequency for the metric was chosen as a random draw during each experimental run. Figure 7 shows the results for optimizing at 3 different frequencies when driving port 1. This figure shows that we were able to realize deep nulls in many arbitrary locations and were able to reduce the power by at least 4 dB in each band attempted. The next step was to determine how the optimization algorithm responded when the metric was defined with different bandwidths. Figure 8 shows the optimization results with 500 kHz, 5 MHz, and 10 MHz bandwidths centered at 3.652 GHz when driving port 3. This figure shows that suppression of P 2 decreases as the bandwidth increases. This is not surprising as a wider bandwidth means more modes and therefore more Degrees Of Freedom (DOF) that need to be matched. We can also see how the shape of the spectrum changes as features outside the narrow bandwidth metric get included in the The optimal set of commands is shown in the middle and the evolution of the metric is shown in the right hand side. The worst case improvement was 7 dB and the best case improvement was 40 dB.
wider bandwidth metrics.
Next, we evaluated the performance when the metric was defined over 2 separated frequency bands. Figure 9 shows the optimization results when driving port 1 and defining the metric as the average power contained in two separated frequency bands. The suppression capability when using multiple frequency bans is reduced; as in the case with a wider bandwidth metric, we are including more features and need to match more DOF.
The last step for the iterative algorithm was to observe the performance when ports 1 and 3 were driven simultaneously with a relative phase shift. Figure 10 shows the results with phase shifts of 0 deg/GHz, 25 deg/GHz, and 50 deg/GHz. This shows we are able to generate coldspots in complex cases where the inputs themselves add spatial and angular diversity.
Finally, we look at the repeatability of the optimization process and determine if we can recover the same final state each time. Figure 11 shows the results from repeating an optimization when driving ports 1 and 3 collectively with a 50 deg/GHz phase shift. The final states are different; however, the metrics followed the same path and the differences are due to the algorithm declaring convergence and stopping at different times. Increasing the convergence criteria would likely result in more closely related final states.

VI. COHERENT PERFECT ABSORPTION STATE GENERATION AND VERIFICATION
A Coherent Perfect Absorption (CPA) state is not guaranteed at any specific frequency, so the approach needs to identify candidates. We repeated the iterative optimization algorithm from the coldspot generation but initialized it by applying a random command set to the metasurface and then finding the eigenvalue with magnitude closest to a pre-selected value. Figure 12 shows the results for 7 experiments with the eigenvalue initialized so that 0.075 ≤ |λ s | ≤ 0.15. This figure clearly shows the directed random walk nature of the iterative algorithm and demonstrates that the eigenvalues can in general move significantly.
Verifying the CPA state requires exciting the scattering system with the corresponding eigenvalue. We used an Agilent N5242A PNA-X Network Analyzer configured for 2-independent source operation. This configuration requires selecting the appropriate signal path on the setup menu and making the jumper connections on the back of the instrument as per the network analyzer documentation. An external phase shifter was connected between port 1 and the cavity. The CPA state condition was tuned by adjusting the relative amplitude between the elements with the power control on the network analyzer and the relative phase between the elements with the external phase shifter. We can perform parameter sweeps The metric was centered at 3.652 GHz and bandwidths of 500 kHz, 5 MHz, and 10 MHz were used. The measured P2 responses are given in the left hand side and show the initial all 0s state as the dashed black line, the final state as the bold red line and the incremental states as the remaining lines. The region where the metric was evaluated is shown as the shaded gray region. The optimal set of commands is shown in the middle and the evolution of the metric is shown in the right hand side. The worst case improvement was 6 dB and the best case improvement was 40 dB. Figure 9. Optimization results when driving port 1 with the metric defined over 2 separated frequency bands. The bandwidths were 500 kHz and the center frequencies were 3.652 GHz and 3.826 GHz. The measured P2 responses are given in the left hand side and show the initial all 0s state as the dashed black line, the final state as the bold red line and the incremental states as the remaining lines. The region where the metric was evaluated is shown as the shaded gray region. The optimal set of commands is shown in the middle and the evolution of the metric is shown in the right hand side. The total reduction in average power was 10 dB.
to determine the sensitivity of the CPA to frequency, relative phase, relative amplitude, and metasurface commands.
Because there are 2 independent sources, S-parameter measurements are not available in this configuration and we need to make use of the network analyzer receivers.
The network analyzer has reference and test port receivers to measure incoming and outgoing signals. R 1 measures the reference signal out of port 1 and R 2 measures the reference signal out of port 2. A measures the signal into port 1 and B measures the signal into port 2. The metric we are interested in is the power ratio Figure 10. Optimization results when driving ports 1 and 3 collectively with a relative phase shift. The metric was defined over a 500 kHz bandwidth centered at 3.33 GHz. The measured P2 responses are given in the left hand side and show the initial all 0s state as the dashed black line, the final state as the bold red line and the incremental states as the remaining lines. The region where the metric was evaluated is shown as the shaded gray region. The optimal set of commands is shown in the middle and the evolution of the metric is shown in the right hand side. The impact of the phase shift can be seen in the difference of the initial P2 response between the various cases and the total power reduction was ∼ 20 dB in each case. Figure 11.
Results of repeating an optimization experiment. Optimization results when driving ports 1 and 3 collectively with a 50 deg/GHz relative phase shift with the metric defined over a 500 kHz bandwidth centered at 3.33 GHz. The measured P2 responses are given in the left hand side and show the initial all 0s state as the dashed black line, the final state as the bold red line and the incremental states as the remaining lines. The region where the metric was evaluated is shown as the shaded gray region. The optimal set of commands is shown in the middle and the evolution of the metric is shown in the right hand side. This shows the results trended in the same direction, the final achieved state depends on of the total outgoing power from the cavity to the total incoming power from the network analyzer, P out /P in = (A + B)/(R 1 + R 2 ).