ABSTRACT
Title of Dissertation: WAVE CHAOS IN MICROWAVE COUPLED
ENCLOSURES, RESERVOIR COMPUTING,
AND PHOTONIC TOPOLOGICAL INSULATOR
GRAPHS: THEORY AND EXPERIMENT
Shukai Ma
Doctor of Philosophy, 2022
Dissertation Directed by: Professor Steven M. Anlage
Department of Physics
Complex scattering exists in many diverse physical and real-life scenarios. Examples
include reactions of atomic nuclei, transport through quantum dots, and the propagation
of electromagnetic (EM) waves in over-moded resonant systems. These systems are
Wave Chaotic, meaning that minute perturbations will lead to a drastic change in the
wave properties of the system. The underlying chaotic property in the short wavelength
limit makes deterministic modeling of wave properties vulnerable to small perturbations.
Because of this, statistical methods play a central role in wave chaotic system studies. The
Random Coupling Model (RCM) has been successfully applied to predict the statistics of
single chaotic EM enclosures. We here expand RCM to systems consisting of multiple
volumes that are coupled together, and do so with highly reduced computational complexity.
Going beyond knowledge-based modeling, we employ machine-learning techniques to
identify hidden information embedded in the scattering properties of wave chaotic systems.
Reservoir computing (RC) is a genre of neural networks employed in machine
learning studies. Its training is radically simplified (compared to a full back-propagation
process in neural networks) because the input and reservoir layers remain unchanged
during the process. Recent work shows that RC can reasonably predict the future evolution
of spatio-temporal chaotic systems. We aim to reverse the thinking: to emulate a software
RC using the spatio-temporal chaotic wave fields in physical EM enclosures. A proof-
of-principle hardware RC is demonstrated experimentally, and tested through a series of
complex tasks carried out at ns-time scales.
The concept of photonic topological insulator (PTI) is translated from the study of
topological insulators (TIs) in condensed matter physics. For a TI material, the charge
will only flow in topologically-protected states on the boundary surrounding the material,
rather than in the bulk. We experimentally demonstrated a novel composite PTI system
with quantum Hall (QH) and quantum spin Hall (QSH) topological phases. The TI
effect also introduces a synthetic spin-1 degree-of-freedom to the guided waves. The
2
Fermionic two-state property, absent in the bosonic photon world, is now accessible using
the PTI system. Using this PTI system, I realize a chaotic graph system that falls in the
Gaussian Symplectic Ensemble (GSE) universality class, which in principle only exists
in the Fermionic world. We use simulations to show that GSE statistics will emerge in an
appropriately designed PTI graph obeying anti-unitary symmetry.
WAVE CHAOS IN MICROWAVE COUPLED
ENCLOSURES, RESERVOIR COMPUTING,
AND PHOTONIC TOPOLOGICAL INSULATOR
GRAPHS: THEORY AND EXPERIMENT
by
Shukai Ma
Dissertation submitted to the Faculty of the Graduate School of the
University of Maryland, College Park in partial fulfillment
of the requirements for the degree of
Doctor of Philosophy
2022
Advisory Committee:
Professor Steven M. Anlage, Chair/Advisor
Professor Yanne K. Chembo, Dean?s Representative
Professor Thomas M. Antonsen
Professor Edward Ott
Professor Rajarshi Roy
? Copyright by
Shukai Ma
2022
Acknowledgments
This thesis would not have been possible without the inspiration and support of
many remarkable individuals ? my thanks and appreciation to all of them for being part
of this journey and making this thesis possible. I owe my deepest gratitude to my principal
supervisor Professor Steven M. Anlage. Without his enthusiasm, encouragement, support,
and endless optimism, this thesis would hardly have been completed. He had faith in my
skills, allowed me to explore seemingly parallel subjects, and helped me grow as an early-
career professional.
Besides my advisor, I express my warmest gratitude to my associate supervisor,
Prof. Thomas M. Antonsen, and Prof. Ed Ott. Their guidance into the domain of
theoretical modeling, statistical analysis, and machine learning has been a valuable input
for this thesis. I want to thank Prof. Yanne K. Chembo and Prof. Rajarshi Roy for
agreeing to be the committee member. I also want to express my gratitude to my collaborators,
Prof. Gregor Tanner, Prof. Gabriele Gradoni, and Prof. Sendy Phang at the University
of Nottingham, Prof. Zhen Peng at UIUC, Prof. Gennady Shvets and Dr. Yang Yu at
Cornell, and Dr. Zachary Drikas and Dr. Bisrat Addissie at Naval Research Laboratory.
Their unwavering help and immense knowledge have significantly influenced many of
the concepts presented in this thesis. Their work has significantly influenced many of the
concepts presented in this thesis.
ii
It is a pleasure to thank my friends at the basement/sub-basement: Bo Xiao, Seokjin
Bae, Bakhrom Oripov, Min Zhou, Jingnan Cai, and Lei Chen, for the beautiful times we
shared, especially the Tuesday group lunches. I am also thankful for the help from Tamin
Tai, Ben Frazier, Arthur Carlton-Jones, and Jared Erb. In addition, I would like to thank
all my friends in College Park who gave me the necessary distractions from my research
and made my stay in Maryland memorable. I would also like to thank Quantum Material
Center for hosting my experiments and taking care of all the administrative matters. It is
impossible to remember all, and I apologize to those I have inadvertently left out.
Finally, my deep and sincere gratitude to my family for their continuous and unparalleled
love, help, and support. I am grateful to my partner, Dr. Qixin Yang, for always being
there for me. If I have ever faltered, she has always stood by to pick me up. Saying
?Thank you? is hardly enough. I am forever indebted to my parents, Mr. Qi Ma and Mrs.
Xiangyun Shen, for giving me the opportunities and experiences that have made me who
I am. They selflessly encouraged me to explore new directions in life and seek my own
destiny. This journey would not have been possible if not for them, and I dedicate this
milestone to them.
iii
Table of Contents
Acknowledgements ii
Table of Contents iv
List of Tables vii
List of Figures viii
List of Abbreviations xviii
Chapter 1: Introduction 1
1.1 Introduction to the Phenomenon of Chaos . . . . . . . . . . . . . . . . . 1
1.2 Introduction to Random Matrix Theory . . . . . . . . . . . . . . . . . . . 4
1.3 Introduction to the Random Coupling Model . . . . . . . . . . . . . . . . 7
1.4 Outline of the Thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
Chapter 2: Random Coupling Model for Multi-Enclosure Wave Chaotic Systems 12
2.1 Overview of the Multi-enclosure RCM . . . . . . . . . . . . . . . . . . . 13
2.2 Multi-enclosure RCM Formalism . . . . . . . . . . . . . . . . . . . . . . 17
2.3 Scaled-structure Experimental Setup . . . . . . . . . . . . . . . . . . . . 18
2.4 Experimental and RCM Results . . . . . . . . . . . . . . . . . . . . . . . 20
Chapter 3: Hybrid Power Balance-RCM Method for Multi-enclosure Wave Chaotic
Systems 26
3.1 Power Balance Method . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.2 Hybrid Method Overview . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.3 Hybrid Method Details Part I: RCM Treatment . . . . . . . . . . . . . . 30
3.4 Hybrid Method Details Part II: Detailed Single Cavity Treatments . . . . 32
3.5 Hybrid Method Details Part III: Finalizing the Hybrid Model . . . . . . . 34
3.6 Experimental Results and Discussion . . . . . . . . . . . . . . . . . . . . 35
3.7 Hybrid Model Limitations . . . . . . . . . . . . . . . . . . . . . . . . . 39
Chapter 4: Application of Machine Learning to Wave Chaotic Systems 42
4.1 Machine Learning Applications, Overview . . . . . . . . . . . . . . . . . 43
4.2 Multi-enclosure Experimental Setup . . . . . . . . . . . . . . . . . . . . 43
4.3 ML Model I: Neural Network . . . . . . . . . . . . . . . . . . . . . . . . 48
4.4 Task I: Multi-cavity System Classification Task . . . . . . . . . . . . . . 50
iv
4.5 Task II: Chaotic Cavity S-parameter Series Prediction . . . . . . . . . . . 54
4.6 Discussion and Summary . . . . . . . . . . . . . . . . . . . . . . . . . . 55
Chapter 5: Late-Time Time-Domain Random Coupling Model 58
5.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
5.2 Details of TD-RCM Formulation . . . . . . . . . . . . . . . . . . . . . . 59
5.2.1 TD-RCM Equation of Motion . . . . . . . . . . . . . . . . . . . 60
5.2.2 System Port Treatments . . . . . . . . . . . . . . . . . . . . . . . 62
5.2.3 Generating of System Eigenmodes with RMT . . . . . . . . . . . 64
5.3 Gigabox Short-pulse Experiment . . . . . . . . . . . . . . . . . . . . . . 65
5.4 Statistical Analysis and Discussion . . . . . . . . . . . . . . . . . . . . . 67
5.5 Discussion and Summary . . . . . . . . . . . . . . . . . . . . . . . . . . 70
Chapter 6: Wave-Based Reservoir Computing 73
6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
6.2 Software Origin of RC . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
6.3 Hardware Realization of RC . . . . . . . . . . . . . . . . . . . . . . . . 77
6.4 Chaotic Wave-based RC . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
6.4.1 Operating Principles . . . . . . . . . . . . . . . . . . . . . . . . 78
6.4.2 Reservoir Expansion Techniques . . . . . . . . . . . . . . . . . . 80
6.4.3 Experimental Set-up . . . . . . . . . . . . . . . . . . . . . . . . 82
6.4.4 Benchmark Test Results . . . . . . . . . . . . . . . . . . . . . . 83
6.4.5 Discussions and Summary . . . . . . . . . . . . . . . . . . . . . 87
Chapter 7: Chaotic 1D Graphs with Photonic Crystal Waveguides 92
7.1 Introduction of PC Defect Waveguide . . . . . . . . . . . . . . . . . . . 93
7.2 Photonic Crystal Graph . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
7.3 Eigenmode Spacing Statistical Analysis . . . . . . . . . . . . . . . . . . 101
7.4 Eigenfunction Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
Chapter 8: Quantum Hall - Quantum Spin Hall Composite Photonic Topological
Insulator 110
8.1 Introduction to PTI Technology . . . . . . . . . . . . . . . . . . . . . . . 111
8.2 Design of QH-QSH BMW Structure . . . . . . . . . . . . . . . . . . . . 113
8.3 Composite QH-QSH Systems . . . . . . . . . . . . . . . . . . . . . . . . 115
8.4 Topological Y-junction . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
8.5 Topological 4-port Circulator . . . . . . . . . . . . . . . . . . . . . . . . 121
Chapter 9: Chaotic Gaussian Symplecticc Ensemble Structures with Photonic
Topological Insulator Graphs 124
9.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
9.2 Design of QVH-QSH Graphs . . . . . . . . . . . . . . . . . . . . . . . . 127
9.3 Symplectic Group Analysis for BMW Systems . . . . . . . . . . . . . . 130
9.4 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134
9.5 Room-temperature Experimental Results . . . . . . . . . . . . . . . . . . 136
v
9.6 Low-temperature Experimental Results . . . . . . . . . . . . . . . . . . . 143
9.7 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147
Chapter 10: Conclusions 148
Appendix A: Multi-cavity RCM Modeling 151
A.1 Single Cavity Treatment . . . . . . . . . . . . . . . . . . . . . . . . . . . 151
A.2 Total Admittance of an N-Cavity Chain . . . . . . . . . . . . . . . . . . 153
A.3 Input and Trans-Impedances . . . . . . . . . . . . . . . . . . . . . . . . 154
Appendix B: Multi-cavity RCM Experimental Details 156
B.1 Aperture and port radiation admittance studies . . . . . . . . . . . . . . . 156
B.2 Convergence of the RCM formalism with inclusion of aperture modes . . 157
B.3 Additional cavity cascade experimental studies . . . . . . . . . . . . . . 159
Appendix C: Multi-channel analysis for the cavity cascade system 160
Appendix D: Applying the Hybrid Model to Generic Cavity Systems 167
Appendix E: Simplified PWB-RCM Hybrid Method 169
Appendix F: Wave-based Reservoir Computing Details 173
F.1 Simulation setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173
F.2 Observer tasks for dynamical systems . . . . . . . . . . . . . . . . . . . 173
F.3 Benchmark tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 175
F.4 RC operation speed analysis . . . . . . . . . . . . . . . . . . . . . . . . 176
Appendix G: Heterogeneous BMW Structure Details 178
G.1 Experimental Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . 178
G.2 Ferrite Material . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 179
G.3 Kink mode Back-scattering Studies . . . . . . . . . . . . . . . . . . . . . 180
G.4 Coupled waveguide issue . . . . . . . . . . . . . . . . . . . . . . . . . . 183
G.5 QH-QSH structure with no external magnetic field . . . . . . . . . . . . . 185
G.6 Topological Y-junction . . . . . . . . . . . . . . . . . . . . . . . . . . . 187
Appendix H: Extended Symplectic Group Analysis of BMW Hamiltonian 188
Bibliography 190
vi
List of Tables
7.1 The summarized consecutive mode spacing ratios of the photonic crystal
graph (PC-Graph) system, the graph system with random ratio order, and
the theoretical predictions for the GOE, GUE and GSE systems [1]. . . . 103
vii
List of Figures
1.1 Figures (a) and (b) show the trajectories of two billiard balls inside a non-
chaotic square billiard, and a chaotic sinai billiard (square billiard with a
center circular block), respectively. Figure (c) shows the quarter bowtie
billiard (open plate view). Figure (d) shows the eigenfunction profile of
one system eigenmode at 12.57 GHz. Figures (a) and (b) are from the
internet. Figure (d) is based on Ref. [2]. . . . . . . . . . . . . . . . . . . 5
2.1 Schematic of the cavity cascade experimental set-up. N identical cavities
are connected in a linear chain with aperture connections. Single-mode
ports are installed at the first and last cavities in the chain. The scattering
properties are measured with a VNA and frequency extenders. Independently
controlled mode stirrers are located in each cavity. a and b refer to the two
sides of a cavity in the cascade. More details on the experimental setup
can be found in Section 2.3. . . . . . . . . . . . . . . . . . . . . . . . . 14
2.2 (a) The scaled 3 cavity cascade case (where the second one is hidden
behind the vertical copper supporting plate) connected to single-mode
WR10 waveguides. The inset shows how the second cavity is arranged
and the location of the two circular apertures. (b) Picture of the full scale
cavity cascade set-up with one wall of cavity 2 removed. The cavities
are connected through circular apertures. Individual perturbers and RF
absorber cones are employed inside the enclosures (not shown). . . . . . . 18
2.3 Comparison of distributions of the imaginary part of trans-impedance
(Im(Zt)) between the experimental curves and the theoretical predictions
for 2 and 3 cavity cascade system. The red curves are for two-cavities and
blue are for three-cavities. Figure (a) and (b) are from scaled measurements
with rectangular and circular shaped apertures, respectively. Figure (c) is
from the full-scale system. The inset shows schematically the shape of
the aperture that is applied in the experiment. . . . . . . . . . . . . . . . 20
2.4 The PDFs of induced voltage on a 50 Ohm load attached to the last cavity
(cav) of the full scale cavity cascade systems (3.95-5.85GHz) and its
scaled counterparts (75-110GHz), assuming 1W input on the single-mode
port of the first cavity. The curves corresponding to the 1, 2 and 3 cavity
system are color-coded in blue, yellow and green, respectively. The full
scale system experimental (theoretical) data are shown at dashed (solid)
lines, and the scaled experiment data are shown as dotted lines. . . . . . . 25
viii
3.1 Schematic illustration of the hybrid model applied to a 3-cavity cascade.
In the hybrid model, we use the PWB method to characterize the power
flow from the first cavity to the next to last cavity. The fluctuations in the
final cavity are described using the RCM method using the mean power
flow values obtained from PWB as an input. . . . . . . . . . . . . . . . . 27
3.2 The PDFs of load induced voltage |UL| of 2- and 3-cavity experiments
(solid) and hybrid model calculations (dashed). The single cavity loss
parameter is varied from 9.7, 7.5, 5.7 and 1.7 from (a-d), respectively.
The inset in (b) shows the 2- and 3-cavity experiment averaged induced
voltage values ?|UL|? with respect to different loss parameters ?. Multi-
mode (? 100 modes) circular apertures are employed between the cavities. 36
3.3 The 3-cavity case load induced voltage |UL| values computed with the
PWB-only model. The single cavity loss parameter is varied from 9.7,
7.5, 5.7 and 1.7 from (a-d), respectively. Multi-mode (? 100 modes)
circular apertures are employed between the cavities. . . . . . . . . . . . 38
3.4 Comparison of induced load voltage statistics for the hybrid model (dotted),
RCM predictions (dashed) and experimental results (solid) for the case of
2 and 3 cavity cascades with single-cavity loss parameter ? = 9.7, and
circular multi-mode apertures between the cavities. The frequency range
is from 3.95 to 5.85 GHz. . . . . . . . . . . . . . . . . . . . . . . . . . . 40
4.1 Schematic of the experimental set-up. We measure the S-parameter of a
3-cavity cascade system with a VNA. The cavities are connected through
circular apertures. Rotatable mode stirrers are employed in each cavity
to generate different system configurations. The scale of the operating
wavelength ?op is shown as the bar. . . . . . . . . . . . . . . . . . . . . 44
4.2 Multi-cavity diagonal impedances Z11 from 2-port measurements on either
1-cavity (blue), 2-cavity (red) or 3-cavity (green) cascade. The inset in (c)
is the schematic diagram of the Z-parameter measurement of 1-, 2- and 3-
cavity cascade structures. (a) and (d): the real and imaginary parts of Z11
from 3.95-5.85GHz for a single realization of the system, whose detailed
views are shown in (b) and (e), and the corresponding statistical analysis
over a large ensemble of Re(Z11) and Im(Z11) are plotted in (c) and (f).
The single cavity loss parameter is ? = 9.7. . . . . . . . . . . . . . . . . 46
4.3 Generalized (recurrent) neural network architecture with two hidden layers.
The NN consists of the one input layer, hidden layers and one output layer.
The RNN is built upon the NN structure with the addition of ?memory
effects? at each hidden layer (shown as the dashed arrows). . . . . . . . . 47
4.4 The NN classification algorithm performance for 1-, 2-, and 3-cavity Re(Z11)
data. (a-c): The cost function fcost (blue, left logarithmic axis) and test
accuracy (red, right linear axis) results are evaluated to 500 iterations of
the algorithm. The dimension of the input vector is 2000, 5000 and 16001
in (a), (b) and (c), respectively. (d): The test accuracy results with and
without additional Gaussian noise. The dimension of the input vector is
5000. The signal to noise ratio values are -5, 5, 15 dB and the original data. 52
ix
4.5 The |S21| prediction of a wave chaotic electromagnetic system future states
using measured |S11|. The measured |S21| data (colored solid lines) and
algorithm predictions (colored dashed lines) from 5 frequency points are
compared for cavity realizations 900 to 1000. The black lines show the
absolute magnitude of the difference between the measured data and the
algorithm prediction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
5.1 Diagrams of the TD-RCM port load modeling. a. The linear load at port j
is represented with an impedance Zj . b. The schematic of a diode-loaded
port. c. The schematic of a multi-mode aperture. Note that the subscript
?in? corresponds to waves following the arrow, going into the cavity. . . 62
5.2 The TD-RCM simulated RX port voltage signal from 0 to 1000ns. Five
realizations with distinct cavity eigenmode spectra and coupling are shown.
The injected waveform is a 5GHz center frequency and 5ns long sine
pulse. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
5.3 Schematic of the short pulse time-domain experiment setup. The injected
waveform (left inset) is generated through IQ modulation using the arbitrary
waveform generator (AWG) and the programmable signal generator (PSG).
The short pulse is broadcast into the enclosure through the transmitting
(TX) port (star). The data ensemble is created by rotating a motorized
mode-stirrer inside the enclosure. The receiver (RX) port (star) induced
voltage signal (right inset) is measured by the oscilloscope. A lab computer
is used for instrument control and data transmission/collection. . . . . . . 66
5.4 We plot the measured gigabox RX port voltage signals in (a), and the
corresponding TD-RCM simulation in (b). In each figure, the port signals
from 20 realizations are overlaid to show the diversity of system response. 67
5.5 The statistics of the maximum RX port voltage studied at different frequencies
merged in the gigabox. The input pulse is a single-frequency 20ns sine
pulse. The center frequency fc changes from 5GHz to 10GHz. Note that
the input pulse has a finite rise/fall time (1 ? 2ns). Data histograms (blue)
are compared to the predictions of the TD-RCM model (red). . . . . . . 68
5.6 Mean values of the receiver port maximum voltage PDFs shown in Fig.
5.5 as a function of the center frequency of the 20ns pulse in the Gigabox.
The error bar is one standard deviation. . . . . . . . . . . . . . . . . . . 69
5.7 The root-mean-square ensemble averaged RX port signal comparison.
The injected pulse is a 5GHz single frequency sine wave with a duration
sweeping from 5 to 30ns. Experimental and TD-RCM simulation results
are shown in a and b, respectively. The insets are the blow-up of the
early-time (before 100ns) curves. . . . . . . . . . . . . . . . . . . . . . 71
x
6.1 Wave chaotic systems for reservoir computing. a, Schematic of the
experimental setup. The input information is transferred from a lab computer
to the AWG and broadcast into the reservoir. The amplified waveform
is injected into the chaotic enclosure through an electric dipole antenna.
Several diode-loaded antennas are used to probe the EM field, whose
voltage signals are measured by the oscilloscope and further transferred
to the lab computer for training. b, Measured port signals (solid) under
random input stream (dashed). c, Schematics of the diode-loaded port. d,
The dynamics of diode-port voltages (red) under single frequency input
at 4 GHz (black). e, The FT of the diode-port signal shown in d. . . . . . 78
6.2 Reservoir Expansion Techniques (RETs). a and b show the schematics
of the boundary condition perturbation method and the frequency-stirring
method, respectively. The first method may be realized by the translation
of a metallic perturber shown as the cylindrical rod. Under the same input
waveform u0, uniquely different evolutions of the EM fields inside the
enclosure area are created by means of translating the perturber to new
locations. In b, the frequency stirring technique utilizes small changes
of center frequency of a given input waveform to create new wave field
configurations. In the experiment, the input wavelet is stretched/shrunken
(shown as u1 ? u3) in order to shift the center frequency. c, Schematics
of the RET concept. The combined reservoir consists of the observable
nodes from various single realizations of reservoirs. . . . . . . . . . . . . 81
6.3 Benchmark tests and the RC performance map on system parameters.
a and c, Ro?ssler observer task testing set performance (blue) from a combined
RC of 630 nodes. The true states are shown as red dashed lines in a, c
and e. b, The deviation of normalized mean square error for the Ro?ssler
observer task: contour as a function of the input duration Tosc and the
system decay time Tdecay. The dashed boundary outlines the optimal
parameter island (error deviation < 20%) from where the combined 630-
node RC is built. e, NARMA10 task testing set performance at the optimal
system parameters from a 90-node RC, shown as the cross in f. f, The
deviation of normalized mean square error for the NARMA10 task. . . . . 84
6.4 Experimental demonstration of the task-versatility. a illustrates that
a single wave chaotic RC is able to execute different tasks with a simple
switch of output couplers, and thus serves as an ideal platform an edge
device. The RC performances (blue) for the function simulator, He?non-
map observer, and the nonlinear channel equalizer tasks are shown in b-d,
respectively. The targets are shown as red dashed lines in b-d. In e-g, we
show the simulated (dashed) and experimental (solid) RC performance as
a function of the dimension of the reservoir. . . . . . . . . . . . . . . . . 88
xi
7.1 a. The open-plate view of the photonic crystal lattice with an L-shaped
defect waveguide region. b. The side view of the PC unit cell. A dielectric
rod is sandwiched between two PEC (perfect electric conductor) surfaces.
The quantities h0 and d0 are the height and diameter of the rod. c. The
photonic band structure from a supercell defect waveguide simulation.
The waveguide modes appear in the bulk bandgap region (delineated by
the blue dashed lines). d. The E-field profile of one of the waveguide
mode solutions (the red circle in c). The black dashed line marks the
center of the waveguide. . . . . . . . . . . . . . . . . . . . . . . . . . . 95
7.2 a. The simulated S-parameters of the PC-graph right-angle junction (solid),
as well as the S-parameter from a 3.5mm uniform coaxial cable 2-port
right-angle connector (dashed), which serves as a reference. Inset: the
schematics of the PC right-angle junction and the 2-port coaxial cable
junction. b. The simulated S-parameters of the PC-graph T-junction
(solid) and the S-parameters from a coaxial cable 3-port ?tee? connector
(dashed). Inset: the schematics of the PC T-junction and the 3-port coaxial
cable junction. All PC-junction S-parameters are frequency averaged with
a 500MHz window to remove the effect of spurious reflections. . . . . . 98
7.3 a. Schematic diagram of the PC waveguide graph system. The graph
bonds are shown as the clear gray channels and the rectangular lattice is
shown as a lattice of dielectric rods (black circles). b. The simulated |Ez|
profile of the graph system eigenmode at 2.8GHz. c shows the statistics
of the normalized mode spacing from graph simulation (histogram). The
theoretical predictions for the GOE, GUE and GSE systems are shown
in red, yellow and purple curves. d shows the spectral rigidity ?(L)
of the normalized spectrum from graph simulation (blue circles). The
theoretical predictions for the GOE and GUE systems are shown in red
and yellow curves. e shows the statistics of the consecutive mode spacing
ratio r (histogram) and the theoretical predictions for the GOE, GUE and
GSE systems. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
7.4 a. Probability amplitude distribution of photonic crystal defect waveguide
graph eigenmodes values (v) obtained from the grid-wise method (symbols
connected by blue line). The theoretical prediction for systems with GOE
statistics are in red. b. shows the statistics of the normalized bond amplitude
value |a| (symbols connected by blue line) and the Gaussian distribution
(red dashed). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
7.5 a. The inverse participation ratio (IPR) of graph eigenmodes obtained
from the bond-value method, computed based on the method IIa. The
simulated |Ez| profile of the red (black) circled mode in a is shown in
c (d). b. The inverse participation ratio (IPR) of graph eigenmodes
computed based on the method IIb. The simulated |Ez| profile of the
red (black) circled mode in b is shown in e (f). . . . . . . . . . . . . . . 106
xii
8.1 (a) Quantum Hall (QH) unit cell model with h0 = 36.8mm, h1 = 1.45mm,
d0 = 12.7mm and d1 = 8.8mm. A QH rod consists of two ferrite
disks (blue) straddling a middle metallic rod, sandwiched by two metallic
plates. (b) Photonic band structure for a QH bulk region. The horizontal-
axis is the kx vector mapping an excursion in the 1st Brillouin zone and the
vertical-axis is the eigenmode frequency. Under Hz = 1350Oe the QH
bulk bandgap (shaded area from 5.9 to 6.3 GHz) matches the quantum
spin Hall bandgap [3]. The inset offers an oblique view of a typical open-
plate BMW hexagonal lattice with lattice constant a0 = 36.8mm. . . . . . 113
8.2 |S21| transmission results through a bulk QH structure in the ? ?K(K ?)
direction for both magnetized and un-magnetized ferrites. With Hz on
(magnetized ferrites, red curve), we observe a decrease in transmission S-
parameter at 5.9-6.3 GHz, which is caused by the absence of bulk modes
inside the bandgap region. Inset: schematic of the experiment. . . . . . . 116
8.3 (a) and (b) are the measured edgemode transmission along the QH-QSH
interface with opposite H-field directions applied to the QH region. The
black-dot(red-line) curves represent the right(left)-going waves. The insets
show the locations of port A and B along with the direction of the applied
H-field. The white arrows mark the expected edgemode propagating directions.
We find that the edgemode propagation direction is reversed when the
biasing magnetic field direction is inverted. . . . . . . . . . . . . . . . . 118
8.4 (a): The Y-junction design consist of three different topological domains:
two QSH domains (green and blue) and a QH region (red), creating three
interfaces: ab, bd and bc. The ??/? represents an edgemode with spin
up/down. Purple/yellow arrows: propagation directions of ??/? modes.
Plot (b-d) are the transmission measurements of the Y-junction with ?Hz
in the QH region. The edgemode will follow the 2 ? 1 and the 1 ? 3
paths strictly without backscattering. Port 3 becomes a forbidden source
as shown in (b). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
8.5 Experimental results of edge waves propagating in a topologically protected
4-port circulator. The QH (QSH with ?SOC > 0 and ?SOC < 0) domains
are color-coded as red (green and blue) for all insets. The applied H-field
direction is inverted between plots (a-b) and (c-d), creating a clockwise
and counter-clockwise circulation around the center QH island, respectively.121
xiii
9.1 a. The schematic of the unit cell structure that define the QVH and QSH
regions. The quantities lv, gv, hv are the prism side length, the air gap
distance, and the prism height of one QVH rod. The quantities ds, hs are
the diameter and height of the QSH cylinder rod. b. The photonic band
structure of the QSH-QVH interface waveguide. Inset shows a schematic
of the supercell simulation model. The red-dashed line marks the interface
between the QSH (left) and the QVH (right) region. c. The side-view
of the QSH-QVH structure (same structure as the b inset). The top and
bottom metallic plates are separated by h0. All QSH rods are attached to
one of the metallic plates, and the QVH rods are mounted on both plates,
leaving an air gap. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126
9.2 The open-plate view of the PTI-graph structure. Shown are the top and
bottom metallic plates and the rods/prisms attached to each one of them.
The top plate has the top part of the QVH rods. On the bottom plate, we
have installed the QSH rods as well as the bottom part of the QVH rods.
Several supporting rods are installed which will spread out the weight of
the top-plate evenly, and provide a fixed and uniform separation between
the plates, h0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
9.3 The simulated eigenmode results of a PTI-graph realization. It is clearly
shown that the two-fold degeneracy is present at the bulk bandgap frequencies.
At the edge of the bandgap (? 23.5GHz), we observe a hybridization of
the edgemode and bulk modes. Inset shows the graph structure used in
this numerical study, consisting of a QSH region show by blue circles
and a QVH region shown by grey triangles. The second inset shows an
enlarged region of the QVH-QSH boundary. . . . . . . . . . . . . . . . 135
9.4 a The histogram approximation probability distribution of the PTI graph
nearest-neighbor mode-spacing s. The prediction of the GOE, GUE, and
GSE groups are plotted as solid lines for reference. b. The integrated
probability function of s along with the theory predictions of the three
groups. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
9.5 Top view of the graph bottom plate, showing all QSH rods, and showing
the bottom half of the QVH rods. The white dashed lines are drawn to
outline the interface between the QVH and QSH domains, which represents
the graph that hosts the topologically-protected modes. The ports locations
are marked by the black stars, along with their numbering. . . . . . . . . 137
9.6 a. The measured reflection S-parameters as a function of frequency in
the bulk bandgap at Ports 3 and 4 of the composite QSH-QVH shown in
Fig. 9.5 at room temperature. The measured complex Wigner time delay
curves near 24.20 and 24.46 GHz, shown in b and c, respectively. The
fitting results are shown in solid red curves. . . . . . . . . . . . . . . . . 140
9.7 a and b show the real and imaginary parts of the ?W as a function of
frequency computed using different values of ?, respectively. The Re(?W )
and Im(?W ) show two divergences as ? is varied (marked out by the blue
arrows) and match the values of ?1 and ?2. . . . . . . . . . . . . . . . . 141
xiv
9.8 a and b show the front and top views of the PTI-graph structure mounted
below the mixing chamber plate, respectively. c. The thermal connection
(copper connectors and heat straps) and cable connection on the PTI-
graph structure. The inset shows the heat-strap fixations on the mixing
chamber plate. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142
9.9 a and b show the real and imaginary part of the measured PTI-graph ?W
versus temperature and frequency, respectively. c and d are the top view
of a and b. To guide the eye, white dashed lines are added to mark the
temperatures at which Re(?W ) changes sign and Im(?W ) diverges. . . . 144
9.10 a shows the 4 cable connections at the top of the fridge. These cables
connect to the VNA in room temperature, and go to the PTI-graphs in low
temperature environment. b and c show the measured S-parameter and
the real part of ?W , respectively. Different colors represent data measured
at 130 different temperatures (from 80K to 4K). . . . . . . . . . . . . . . 146
10.1 Summary of the projects. There are three themes in the thesis: wave
chaos, microwave topological photonics, and machine learning. The corresponding
chapter indices for the projects are also shown in this diagram. Wave
chaos takes the center stage. The topological photonics and the machine
learning projects also intersects with the wave chaos themes. . . . . . . . 149
A.1 Schematic diagram of the ith cavity in a cascade with port a on the left
and port b on the right. Here U and I refer to the voltage and current at
each port (a vector quantity in general). . . . . . . . . . . . . . . . . . . 152
B.1 The model set-up for aperture admittance matrix calculations in the CST
simulations. A PEC plate is carved with the dimension of the aperture,
and attached to a vacuum box. The port is assigned at the 2D surface of
the aperture. The wave enters the aperture (blue arrow) and flows into the
vacuum space before absorbed by the surrounding radiation boundaries
(orange arrows). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157
B.2 Statistics of the imaginary part of the trans-impedance for a scaled 3-
cavity cascade with circular aperture connections, illustrating the convergence
of the theoretical model with increasing mode number included in the
aperture admittance matrix. The circular aperture allows 100 propagating
modes. The model PDFs curves are nearly unchanged when the non-
propagating modes are included. . . . . . . . . . . . . . . . . . . . . . . 158
B.3 (a) and (b): Statistics of the imaginary part of the trans-impedance of 2-
and 3-cavity cascades when the loss of the system changes. (c): Comparison
of the Re(Zt) statistics experimental and model results for the scaled 2
cavity set-up, with circular aperture connection. In (a), (b) and (c), the
single cavity loss parameter ? is measured to be ? = 5.7, 7.5 and 9.1,
respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 159
xv
C.1 Schematic of a cavity cascade array with many-channel connecting apertures
between all cavities except that between cavity i and cavity i + 1. Note
the definition of input and output powers at each cavity aperture. . . . . . 160
E.1 The PDFs of load induced voltage |UL| of 2- and 3-cavity experiments
(solid), hybrid model (dashed), and the simplified hybrid model (dash-
dotted). The single cavity loss parameter is varied from 9.7, 7.5, 5.7 and
1.7 from (a-d), respectively. . . . . . . . . . . . . . . . . . . . . . . . . 170
G.1 An open-plate view of the experimental set-up of the 4-port BMW circulator
structure. Each region is color coded to match the insets in Fig. 5. . . . . 179
G.2 The photonic band structure of a 15 + 15 period QH-QSH Supercell
structure from COMSOL simulations. The insets are the supercell structures
near the QH-QSH interface. The dashed lines represent the location of the
QH-QSH interface. Under applied magnetic field in opposite directions
(?Hz/ + Hz in (a)/(b)), the supported kink mode will have opposite
synthetic spin (up/down-spin shown as ??/?) and group velocity. Dirac
points are at kxa0 = ?2?/3. Boxed tags: the Chern numbers of the bulk
modes below the bandgap for QSH and QH domains. . . . . . . . . . . . 181
G.3 An open-plate view of the QH-QSH interface 2-port isolator experimental
set-up. Note that the biasing magnets on the QH region are not shown. . . 182
G.4 Plots of kink mode transmission measurement in (a) propagating and (b)
non-propagating directions along a QH-QSH interface. In both plots,
the S-parameters represent the transmission between two ports: S2 is
measured at a port further from the source (port A) than S1. The plots
show that the amplitude of the kink mode along the propagating direction
does not attenuate significantly, while the amplitude of kink mode in the
backward direction is always low. . . . . . . . . . . . . . . . . . . . . . . 182
G.5 The 3-QSH supercell PBS near the Dirac point. The positive/negative
group velocity parts (blue/red) are mixture of the up/down-spin mode
localized at the top interface and down/up-spin mode at the bottom interface.
Near the secondary bandgap at about 6.06-6.09 GHz, there is a four-mode
mixing. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 184
G.6 Supercell simulations results of the kink mode at the Dirac point. Color:
|E| of the eigenmodes near 6.1 GHz. The 5 and 7 periods of unit cell
separation cases cause E-field spreading in the entire middle QSH region,
while the 9-period separation case makes the kink modes better localized
at the interfaces. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 185
G.7 Zero magnetic field transmission measurement between port A and B.
Two ports are in the same locations as in the measurement shown in Fig.
3, located on the QH-QSH interface. We find that the transmission profile
between two ports are very similar, indicating the time-reversal symmetry
in not broken in the un-magnetized ferrite case. . . . . . . . . . . . . . . 186
xvi
G.8 The open-plate view of the topological Y-junction. The magnets for QH
region and the upper QSH rods and are not shown in this view. The three
topological regions are color-coded as in Fig. 4(a). The dashed lines
mark the interface between different topological phases. The stars are the
locations of the ports. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 187
xvii
List of Abbreviations
AWG Arbitrary waveform generator
BMW Bi-anisotropic meta-waveguide
CPA Coherent perfect absorption
COMSOL COMSOL Multiphysics
CST Computer simulation technology
DL Deep learning
EM Electromagnetic
EMC Electromagnetic compatibility
GHz Gigahertz
GOE Gaussian orthogonal ensemble
GUE Gaussian unitary ensemble
GSE Gaussian symplectic ensemble
HFSS High-frequency structure simulator
LSTM Long-short term memory
ML Machine learning
MXC Mixing chamber plate
NC Neuromorphic computing
NN Neural network
PC Photonic crystal
PDF Probability density function
PEC Perfect electrical conductor
PTI Photonic topological insulator
PSG Programmable signal generator
PWB Power balance method
xviii
QH Quantum Hall
QVH Quantum valley-Hall
QSH Quantum spin-Hall
RC Reservoir computing
RCM Random coupling model
RET Reservoir expansion technique
RMT Random matrix theory
RNN Recurrent neural network
RIS Re-configurable intelligent surface
RX Receiver
SO Short Orbit
TD Time domain
TI Topological insulator
TX Transmitter
VNA Vector network analyzer
xix
Chapter 1: Introduction
1.1 Introduction to the Phenomenon of Chaos
?Chaos? is a common concept in life. A good hint for ?chaoticity? is to examine
the complex behavior of the system and see if it is sensitive to the initial conditions [4].
One may consider two states of the system with very similar initial conditions. In a
chaotic system, these two states can evolve quickly into starkly different states under
time evolution. Such a characteristic of chaotic systems makes the precise analysis of
any long-term behaviors challenging, for example, the weather prediction for an extended
period of time. A great example of the phenomenon of chaos is the famous butterfly
effect proposed by Lorenz in 1972 [5]. The butterfly effect presents a poetic picture of
the sensitivity of the outcome to minute changes, expressed as the exact formation time
of a tornado in Texas and a butterfly flapping its wings in the Amazonian forest. Such
an effect is very different from the linear system where the superposition principle will
hold. However, linear systems are usually the idealized simplifications of many real-
life phenomena. Many real-life systems can be considered nonlinear, meaning that the
outcome of two closely-spaced initial conditions is not proportional to the initial condition
difference. Examples include the swinging of a double pendulum, the trajectory of billiard
balls in a chaotic billiard, and the mechanism underlying self-oscillating circuits.
1
Suppose we go back to the butterfly effect as an introduction to chaotic phenomena.
In that case, one may argue that the above example is more attributed to the lack of precise
modeling of the weather system, and/or the unforeseeable events that may perturb the
weather condition, rather than the ?chaotic? nature of the system. Even though weather
forecasting is challenging to model, we mean to say that even with the precise knowledge
of a particular system, for example, knowing all equations of a dynamical system, the
initial condition would need to be precisely known to produce a future state prediction of
good accuracy. Here we use the double pendulum system to illustrate the idea better. The
governing equations can be written out explicitly, and the computation of the pendulum
motion can be quickly conducted on a personal computer with reasonably small numerical
error. However, the system is still chaotic as a slight initial condition difference in the
pendulum angle will lead to very different future evolution. This is a good example of
a deterministic system with chaotic dynamics. One may produce a good approximation
of the system?s future evolution within a specific error tolerance range. However, this
prediction is only sound in the near-term future.
Besides the qualitative description of chaos, here we introduce a more mathematical
way of characterizing the chaoticity of a system. For a dynamical system, a Lyapunov
exponent can be calculated for each degree of freedom of the system. The Lyapunov
exponent is defined as the separation rate, with respect to time, between two close trajectories.
A chaotic system must have at least one positive Lyapunov exponent. For example, the
Rossler and Lorenz system will each have 3 Lyapunov exponents corresponding to the
2
x, y, z variables [6, 7]. The chaotic Rossler system is governed by
dx/dt = ?y ? z,
dy/dt = x+ ay, (1.1)
dz/dt = b+ z(x? c),
where typical parameter values are a = 0.5, b = 2.1, c = 3.5. The Lorenz attractor is a
set of equations that models fluid motion. The chaotic Lorenz system is expressed as
dx/dt = ?ax+ ay,
dy/dt = bx? y ? xz, (1.2)
dz/dt = cz + xy,
where typical parameter values are a = 10, b = 28, c = 8/3. The Lyapunov exponent
of the dynamical variable z for the two above systems will be positive. Depending on
the choice of system parameters, the system can either fall into the chaotic regime or
not. The specific value of the Lyapunov exponent is linked to the divergence rate of
system states. One may evaluate the Lyapunov exponent values of a particular system
through numerical calculation, and the process of ensemble averaging over different initial
conditions is usually applied. Note that we will re-encounter these equations in Chapter
6 in the context of reservoir computing with wave chaotic systems.
To sum up the section, chaotic phenomena are commonly found in real-life scenarios
and also scientific studies. The detailed, deterministic study of chaotic systems are usually
computationally intensive. Moreover, the computed results are tightly linked to the initial
3
conditions. Any uncertainty in the initial conditions will eventually invalidate the computed
deterministic solutions. I will introduce a better way to characterize chaotic systems,
which is to look at the statistical descriptions of the quantity of interest. In the next two
sections, I will introduce the Random Matrix Theory (RMT) and the Random Coupling
Model (RCM). Random Matrix Theory is a well-studied method of characterizing the
spectral properties of wave chaotic systems. Here at Maryland, our Wave Chaos group
has invented the Random Coupling Model for wave chaotic systems. It will be clear later
that RMT is an important part in the RCM.
1.2 Introduction to Random Matrix Theory
The Random Matrix Theory (RMT) is proposed by Wigner and later developed
by Dyson and Mehta [8, 9, 10]. RMT is originally developed to study the energy level
spacing statistics in nuclear physics. The study of spectral properties investigates the
correlation between system energy levels, such as the level spacing between nearest-
neighbor levels and the total number of levels within a given frequency range. Differing
from the level-repulsion phenomenon in chaotic systems (i.e., sinai billiard in Fig. 1.1
(b)), the nearest-neighbor level spacing of systems with integrable dynamics (e.g., square
billiard in Fig. 1.1 (a), and circular billiard) shows Poissonian statistical behavior. Note
that one usually studies these statistical properties of a large number of quantum systems
rather than a particular realization. Besides studying the eigenvalue spacing statistics,
there is also studies for the eigenfunction profiles for chaotic system states. The eigenfunction
study is related to the distribution of the probability density of system eigenstates.
4
Figure 1.1: Figures (a) and (b) show the trajectories of two billiard balls inside a non-
chaotic square billiard, and a chaotic sinai billiard (square billiard with a center circular
block), respectively. Figure (c) shows the quarter bowtie billiard (open plate view). Figure
(d) shows the eigenfunction profile of one system eigenmode at 12.57 GHz. Figures (a)
and (b) are from the internet. Figure (d) is based on Ref. [2].
5
Depending on the system?s symmetry, the spectral properties of a certain chaotic
system can be categorized into one of 3 groups. The first two groups are the Gaussian
Orthogonal Ensemble (GOE) and Gaussian Unitary Ensemble (GUE) of random matrices
for a system with and without time-reversal symmetry, respectively. The third group
is the so-called Gaussian Symplectic Ensemble (GSE) of random matrices for a system
with half-integer spins and time-reversal symmetry. It is often the case that the spectral
statistical properties of systems from the same group show a universal behavior. It was
later conjectured by Bohigas-Giannoni-Schmit (BGS conjecture) that the spectral statistical
properties of systems that show chaos in the classical limit are well described by the
statistics of random matrices [11]. Based on the BGS conjecture, RMT has found a
wider range of application in many other chaotic physical systems (atoms, quantum dots,
quantum wires, microwave enclosures, and others) [12, 13].
The field of microwave wave chaos aims to study the EM field propagation in
electrically large, chaotic shaped enclosures/resonators. In the regime where the microwave
wavelength is much smaller than the dimension of the system, the propagation of waves
can be thought of as particles traveling in straight trajectories and having specular reflection
at the system boundaries. Because the systems have a chaotic shape, the future propagation
of the ray trajectories are extremely sensitive to the boundary condition details. Fig. 1.1
illustrates the behavior of wave chaotic enclosures. Fig. 1.1 (a) shows the ray orbits inside
a non-chaotic square billiard. After several bounces, the two billiard balls still keep the
same displacements as the initial value. However, in a chaotic billiard shown in Fig.
1.1 (b), the displacement between the two billiard balls will increase exponentially with
time. Thus chaotic billiard systems are sometimes said to be ray-chaotic. Many numerical
6
and experimental studies have shown that the properties of chaotic microwave billiards
and graphs agreed nicely with the RMT predictions [14, 15]. In microwave experiments,
the chaotic dynamics of the system can be nicely manifested in measurable quantum
quantities (i.e., scattering parameters). Because of this, many experimental studies of
quantum chaos are conducted in microwave wave chaotic systems. In our lab, we studied
different types of microwave wave chaotic systems, ranging from 1D graphs, 2D billiard,
and 3D enclosures. Fig. 1.1 (c) shows the bowtie-billiard used in our lab. Because of
its chaotic shape, the eigenfunction for a system eigenmode at 12.57GHz shows highly
complicated patterns (Fig. 1.1 (d)).
1.3 Introduction to the Random Coupling Model
It is of interest to study the scattering properties of complex ray-chaotic systems in
the semi-classical limit [16]. Examples include reactions of atomic nuclei [12], quantum
dot transport [13], and the flow of electromagnetic (EM) waves through electrically-large
resonant systems (defined as those for which ? ? V 1/3, where ? is the wavelength of light
and V is the system volume) [17, 18, 19, 20, 21]. These settings are typically found to be
ray-chaotic and highly over-moded (i.e., many resonant modes at and below the frequency
of interest), posing challenges to both numerical and experimental analysis means. A
ray-chaotic enclosure has the property that two rays launched with slightly different
initial conditions will separate exponentially with time as they continue bouncing from
boundary walls and obstructions inside the structure [4, 22]. On the other hand, a minute
change of the structure boundary condition can drastically affect the pre-established field
7
profile inside the system [23, 24, 25]. Though deterministic approaches are available
[26], the chaotic properties make the numerical solutions vulnerable to small changes and
uncertainties of interior structure details.
As introduced above, the statistical and/or approximate approaches in these wave
systems can be more beneficial than deterministic methods (e.g., detailed numerical computations
for a specific configuration). The Random Coupling Model (RCM) determines the statistical
properties of the impedance (Z) and scattering (S) parameters for complex enclosures
[3, 8, 27, 28, 29, 30, 31]. In contrast to the other methods mentioned above, the RCM
generates both deterministic and statistical predictions, treats interference, and utilizes
a minimum of information, namely, the system coupling details, and the enclosure loss
parameter [32, 33, 34, 35, 36]. The scattering parameters describe the steady-state electrical
transmission and reflection properties of a system. For an N-port system, the N-by-N S-
matrix connects the incoming and outgoing voltage waves by [V o] = [S][V i]. Here the N-
by-1 vector [V o,i] represents the outgoing/incoming complex voltage waves at all N ports.
The impedance parameters connects the total voltages and currents at the terminals of the
port transmission lines in an electrical system by [V ] = [Z][I]. The Z-parameters and
S-parameters can be computed by knowing the other one. Both parameters are of great
importance in terms of characterizing the reflection an?d tran?smission properties of an?0 1?
electrical system. The S-matrix can be written as |S| = ?? ?? for a 2-port reflectionless
1 0
and lossless transmission line. The S-matrix and Z-matrix are linked by the bi-linear
transformation S = Z1/2(Z + Z )?1(Z ? Z )Z?1/2, where Z is a diagonal matrix
0 0 0 0 0
whose elements are the transmission line characteristic impedance (usually 50 Ohm).
8
The RCM cavity impedance matrix of an individual cavity is created via a normalized
impedance matrix ? derived from a random matrix ensemble [27, 28, 37, 38]. The
RCM
distribution function of the complex fluctuating ? is governed by a single parameter,
RCM
called the loss parameter ?. In terms of the normalized impedance matrix, the fluctuating
cavity impedance matrix is written as
Z = iIm[Z ] +Re[Z ]1/2 ? ? ?Re[Z ]1/2,
cav avg avg RCM avg
where the quantity Z is the ensemble averaged system impedance (discussed more in
avg
the next two chapters). Here, ? is defined as
RCM
i ?N? wnwT? = n ,
RCM ? (k20 ? k2n)/?k2n + i?n
where the sum over n represents a sum over the modes inside the closed cavity (n =
1, 2, . . . , N). The vector wn, whose number of elements equals the number of ports,
consists of independent, zero mean, unit variance random Gaussian variables, which
represent the coupling between each port and the nth cavity mode. This random choice
of mode-port coupling originates from the so-called Berry hypothesis, where the cavity
modes can be modeled as a superposition of plane waves with random direction and
random phase [39]. The quantities k0 and kn are wave-numbers corresponding to the
operating frequency ?0 = k0 c and the resonant frequencies of the cavity modes, ?n =
kn c. Rather than use the accurate resonant frequencies of the cavity, a representative set
of frequencies is generated from a set of eigenvalues of a random matrix selected from the
9
relevant Random Matrix Theory ensemble [18, 37, 40, 41]. These RMT eigenvalues are
appropriately normalized to give the correct spectral density via the parameter ?k2n, the
density of states near the operating wave-number. The loss parameter ? = k2/(Q?k2n)
where k is the wave vector of interest and Q is the quality factor [31]. The loss parameter
varies from 0 (no loss) to ? (high loss). The RCM loss parameter is also related to the
finesse parameter F that is widely used in the optical cavities. We have ? = F?1 = k
2Q?k
where the ?k2 = 2k?k approximation is used [42].
1.4 Outline of the Thesis
The next chapter will introduce our effort to extend the Random Coupling Model
(RCM) from a single chaotic enclosure to more complex systems with multiple coupled
chaotic sub-systems. Although the full RCM description of multiple enclosure systems
coupled with electrically-large apertures is successful, it is computationally intensive
when the coupling channel (aperture size) grows. In chapter 3, I further discuss a hybrid
multi-enclosure method formulated with both RCM and the power balance model. This
new method is much simpler and faster than the full RCM treatment. To better utilize
the large ensemble of measured data collected for the multi-enclosure projects, we have
applied machine learning techniques to distinguish subtle details in the experiment that
are not visible to the eye (chapter 4). I further study the time-domain modeling of wave
chaotic systems in chapter 5. In chapter 6, I employed the rich dynamics of the chaotic
EM waves to execute several different machine learning tasks. Unlike the 3D cavity wave
chaos studies in chapters 2 to 5, in chapter 7, I studied the properties of quantum graphs
10
based on defect photonic crystal waveguides. The theme of chapter 8 is not centered
around wave chaos. Here we discuss our experimental realization of a composite photonic
topological insulator (PTI) system with QH and QSH domains. Chapter 9 turns back to
the chaos-related study, where I utilize the spin-1/2 property of PTI system modes to
realize GSE statistics using PTI graphs. I will summarize the thesis and discuss future
directions in chapter 10.
11
Chapter 2: Random Coupling Model for Multi-Enclosure Wave Chaotic
Systems
This work was a collaborative effort with the Naval Research Laboratory group and
is published as Ref. [16]. In the previous chapter, we introduced the study of single wave
chaotic systems using the Random Coupling Model. We think that concatenating two or
more such systems is also of great interest, but not as extensively studied. Such coupled
cavity systems can be realized in a wide range of physical platforms from inter-connected
photonic crystal cavities [43], to Cooper pair boxes in superconducting resonators [44], to
microwave (MW) billiards [45], and the complex acoustic and electromagnetic environment
in ships and aircraft containing multiple connected compartments [46]. It has proven
possible to perform experiments for such inter-connected systems and to measure transmission
as a function of coupling. Examples include measurements of conductance of quantum
dots systems with coupled electron billiards [47, 48, 49], resistance of disordered nanowires
modeled by a cascade of quantum dots [50, 51, 52, 53], and simulating resonance strength
of coupled quantum mechanical systems with superconducting MW billiards [54], etc.
Likewise, the EM wave properties of inter-connected electrically large enclosures, like
the power flow and the impedance or scattering parameters, are also widely studied in
engineering [55, 56] and realistic situations ranging from computer enclosures to rooms
12
or buildings [46, 57, 58, 59]. It was recently demonstrated that two single enclosures with
a specific scaling relationship concerning size, frequency, and the RCM loss parameter
? share the same normalized impedance statistics [60]. In this work, the authors study
a scaled cavity (75-110 GHz, cavity longest dimension = 6.375 cm) and its ?20 full
scale version (3.7-5.5 GHz, cavity longest dimension = 127.5 cm). The degrees of
lossyness for the two systems are matched by increasing the scaled cavity Q-factor in
a low-temperature environment. The loss parameter ? describes the degree of lossyness
of the system. With ? known, one is able to generate the statistical ensemble of system
impedance by generating new system modes and mode-port coupling coefficients using
RCM. In the experiment, an ensemble of system impedance can be created by rotating
a mode-stirrer inside the cavity. The impedance statistics of the full scale and scaled
cavities were found to be the same, and in agreement with RCM predicted statistics for
the same ?. This work paves the way for experimentally testing the wave properties of
large, coupled complex systems in a typical lab environment by studying their scaled-
down-in-size counterparts [61, 62].
2.1 Overview of the Multi-enclosure RCM
The statistics of the scattering of waves inside single ray-chaotic enclosures have
been successfully described by the Random Coupling Model (RCM). We expand the
RCM to systems consisting of multiple complex ray-chaotic enclosures with various
coupling scenarios. The statistical properties of the model-generated quantities are tested
against measured data of electrically large multi-cavity systems of various designs. The
13
Figure 2.1: Schematic of the cavity cascade experimental set-up. N identical cavities
are connected in a linear chain with aperture connections. Single-mode ports are installed
at the first and last cavities in the chain. The scattering properties are measured with
a VNA and frequency extenders. Independently controlled mode stirrers are located in
each cavity. a and b refer to the two sides of a cavity in the cascade. More details on the
experimental setup can be found in Section 2.3.
14
statistics of model-generated trans-impedance and induced voltages on a load impedance
agree well with the experimental results. The RCM coupled chaotic enclosure model is
general and can be applied to other physical systems including coupled quantum dots,
disordered nanowires, and short-wavelength electromagnetic and acoustic propagation
through rooms in buildings, aircraft and ships.
Here, we apply an RCM-based model which can be used to make statistical predictions
of impedance values in inter-connected systems of chaotic cavities. The multi-enclosure
RCM method builds on an earlier work for single-mode connection between cavities
[33]. As shown schematically in Fig. 2.1, our approach first adopts RCM to describe
each individual chaotic enclosure [31, 60]. The system-specific details are identified in
orange in Fig. 2.1. These consist of the single-mode waveguide ports and the multi-
mode apertures between cavities, and are described by the radiation impedance Zport,
and the radiation admittance matrix Yaper, respectively. With known geometry, these
frequency dependent complex coupling quantities can be determined through either full-
wave numerical simulations or direct radiation measurements. For an N-mode aperture
we utilize the aperture admittance matrix Yaper (an N?N matrix) to describe its deterministic
properties as a function of frequency. We use Zport to represent the deterministic properties
of the single-mode ports in a manner identical to previous treatments of the port radiation
impedance [30, 37, 60, 63].
An important consideration is the number of aperture modes (both propagating and
evanescent) to include in Yaper. The single cavity radiation admittance matrix can be
15
written as a block-diagonal complex and frequency-dependent matrix
?? ?Yrad,a 0 ?
Yrad = ?? ?? (2.1)
0 Yrad,b
As shown in Fig. 2.1, the Yrad subscripts ?a? and ?b? refer to the connection on the left
and right sides of the cavity in the linear chain. The choice of Yrad solely depends on
the specific cavity connection (i.e., a port or aperture that allows Na modes in ?a? and
Nb modes in ?b?). The off-diagonal zeros reflect the assumption that the apertures and
ports are sufficiently separated such that no direct connection exists between them. This
assumption can be lifted if direct illumination exists between apertures and ports [42].
The RCM single cavity admittance matrix Ycav is then constructed as shown in Eqn. 2.2
from Yrad and the (Na + Nb) ? (Na + Nb) dimensionless universal fluctuation matrix
?RCM whose statistics are determined solely by the loss parameter ? of the cavity,
( ) ( )0.5 ( )0.5
Ycav = i ? Im Yrad + Re Yrad ? ?RCM ?Re Yrad (2.2)
The matrix ?RCM can be calculated using a Monte-Carlo technique [27, 28, 58, 60]. The
single cavity admittance matrix Ycav reflects the chaotic universal fluctuations from ?RCM
[31], ?dressed? by the system-specific properties of the ports and apertures described by
Yrad.
16
2.2 Multi-enclosure RCM Formalism
For a description of the statistical properties of the entire cavity cascade, cavities
can be connected together by enforcing continuity of voltages and currents at the cavity
coupling planes. ?????? ? ? ????
???
??Ia??? ?Ua?= Y (i) ?cav ? ? ??
????
I U (2.3)b b
? ?Ib = Y (i+1)L ? Ub
In Eqn. 2.3, the superscript i refers to the index of the cavity running from 1 to N as
shown in Fig. 2.1. Y (i+1)L is the overall load admittance including the (i+1)st cavity and
everything beyond it. By solving Eqn. 2.3 in matrix form, we have:
( )?1
Ub = ? Y (i) + Y (i+1) Y (i)bb L ba Ua (2.4)
( )?1
Y (i) = Y (i) ? Y (i) Y (i) + Y (i+1) (i)L aa ab bb L Yba (2.5)
Eqn. 2.4 connects the voltages on two sides of a cascade unit, and Eqn. 2.5 gives the
YL recursion relationship which calculates the load admittance Y (i)L of the ith cavity given
the knowledge of Y (i+1)L . With this recipe, the complete RCM cavity cascade model is
created as follows. Starting by computing the load admittance Y (i)L presented by the load
impedance at the N th cavity ( (N)YL ) and working back to the first cavity using Eqn. 2.5.
Combined with the single cavity admittance matrices and the input voltage at the first
cavity, we use Eqn. 2.4 to calculate the voltages at each cavity connection, and finally at
the output port of the N th cavity. (Detailed information of the Yab matrices can be found
17
Figure 2.2: (a) The scaled 3 cavity cascade case (where the second one is hidden
behind the vertical copper supporting plate) connected to single-mode WR10 waveguides.
The inset shows how the second cavity is arranged and the location of the two circular
apertures. (b) Picture of the full scale cavity cascade set-up with one wall of cavity 2
removed. The cavities are connected through circular apertures. Individual perturbers
and RF absorber cones are employed inside the enclosures (not shown).
in Appendix A.) One can then make predictions for the statistics of Zin and Zt of the
entire system (see Eqns. A.8, A.9 in Appendix A) based on the minimal information of
cavity loss and system coupling.
2.3 Scaled-structure Experimental Setup
The details of the cascade cavity experimental set-up are shown in Fig. 2.2. The
scaled and full scale 3-cavity cascade structures are shown in Fig. 2.2 (a) and (b),
respectively. The scaled experiments are conducted at the University of Maryland, and
the full scale version at the Naval Research Laboratory. My predecessor Bo Xiao acquired
the equipment and set up many parts of the scaled structure facility. I took over in June
of 2017. In our experiments, we use the Vector Network Analyzer (VNA) to conduct
S-parameter measurements. The S-parameter (S-matrix) is introduced in the previous
chapter. The VNA can produce faithful measurement of the full S-matrix for an electrical
18
system. Here we will briefly describe the experimental method. Before the measurement,
we will first calibrate the VNA so that the calibration plane is right at the end of the
VNA cable. A good quality calibration will essentially render the VNA as a matched
load throughout the entire bandwidth. During the measurement, the VNA will first inject
a signal into the system from one of its ports, and measure the reflection/transmission
signals at the other ports. In our cascade cavity setup, energy will be injected into the first
cavity through the waveguide port, then transmitted into the rest of the cavities through
the apertures, and finally leak out through both the transmitting port and the receiving
ports. The full S-matrix is measured by computing the ratio between the measured
reflected/transmitted signals to the injected signal. The resulting S-parameters are complex
functions of frequency. In all our experiments, we measure the S-matrix with the VNA
and then store the data on a workstation. Other quantities (Z-matrix, induced voltages
values at the ports, etc.) can be computed from the S-matrix. For our setup, the characteristic
length of a single cavity is about 9 cm. We conduct measurements from 75 to 110 GHz
with the help of frequency extenders made by Virginia Diodes (Tx/Rx WR10 module).
At room temperature and mmWave frequencies, the single cavity Q-factor is in the range
of 5000 to 6000 (which results in RCM loss parameter values ? ? (7.5, 9.5)). Note that
we connected the frequency extenders to the VNA to boost the bandwidth from 26GHz to
mmWave frequencies. More details of the experimental setup are included in Appendix
B.
19
Figure 2.3: Comparison of distributions of the imaginary part of trans-impedance
(Im(Zt)) between the experimental curves and the theoretical predictions for 2 and
3 cavity cascade system. The red curves are for two-cavities and blue are for three-
cavities. Figure (a) and (b) are from scaled measurements with rectangular and circular
shaped apertures, respectively. Figure (c) is from the full-scale system. The inset shows
schematically the shape of the aperture that is applied in the experiment.
2.4 Experimental and RCM Results
To kick-off the build-up of multi-cavity RCM, the loss parameter of each single
cavity is required. The loss parameters, obtained from single cavity measurements, are
? = 9.1 for the scaled cavity and ? = 9.7 for the full scale cavity with 6 absorber cones
present in the cavity [16]. With the single cavity losses made equal and the physical
dimensions carefully scaled, nominally identical electromagnetic conditions are achieved
between the two set-ups [60]. We have put the scaled cavity system into the dilution
refrigerator to decrease the loss of the system. Since all the cavities inside the linear array
are electrically identical (nominally), the same cavity loss parameter ? will be assigned to
each cavity in the cascade. The port and the aperture radiation admittances are obtained
by direct measurements and numerical simulations in Computer Simulation Technology
(CST), respectively (see Appendix B. (a)). Combining this information, we generate an
ensemble of 2 and 3-cavity cascade system impedances, and compare its statistics with
20
those of the measured data. Note that the multi-enclosure RCM model is a numerical
method rather than an analytical method. We conduct Monte-Carlo simulations of random
matrices so that a statistical ensemble of data can be generated by the RCM program. In
Fig. 2.3, we compare the statistics of Im(Zt) for cavity cascade systems with various
cavity numbers, coupling strengths and physical dimensions. Fig. 2.3 (a) show PDF?s
of Im(Zt) for 2 and 3 cavity systems when the cavities are connected by rectangular
shaped apertures having 5 propagating modes (10?4 coupling strength). The statistics
of the 2 and 3 cavity cascade theory generated Im(Zt) (solid lines) match well with the
measured data (dashed lines). A minor mismatch between the 2-cavity Im(Zt) statistical
comparison is observed. Compared to the measurements, the theory generated Im(Zt)
PDFs have higher peak values near zero, which indicates that the theory predictions tend
to have smaller magnitude values. The inter-cavity coupling strength is increased to 10?2
by employing the circular shaped aperture, and as shown in Fig. 2.3 (b), this increases
the magnitude of fluctuations of the trans-impedance for both the 2 and 3 cavity systems.
Results for the full-scale measurements with circular aperture connections are shown in
Fig. 2.3 (c). Good agreement between model and measurements are obtained in all cases.
Additional cavity cascade experimental results are shown in Appendix B (c).
Aside from validating the prediction performance of the RCM cavity cascade model,
another key aspect of our experiment is to study the miniature cavity technique for the
multi-cavity systems. As introduced in section I, the full scale cavity is built with linear
dimensions 20 times larger compared with its scaled counterpart. With the operating
frequency properly scaled and loss parameter made equal by adjusting the wall conductivity,
the statistical wave properties of the two set-ups are expected to be identical [60]. The
21
direct comparison of system trans-impedance statistics can be examined by comparing
the Im(Zt) PDFs shown in Fig. 2.3 (b) and (c). The peak values and spreads of the
scaled and full scale 2 and 3 cavity Im(Zt) PDFs are in similar ranges. For the scaled
2 and 3 cavity Im(Zt) PDFs, the peak values are 0.17 and 0.28 and the FWHMs are
5.1 and 3?, respectively. For the full scale 2 and 3 cavity Im(Zt) PDFs, the peak
values are 0.15 and 0.21 and the FWHMs are 5.8 and 4?, respectively. We believe
that this imperfect agreement between scaled and full scale Im(Zt) PDFs is caused by a
difference in the aperture thickness. The circular aperture thickness is about 1?op in the
scaled cavities, but only 0.04?op in the full scale set-up (?op represents the characteristic
operating wavelength). Note that the finite thickness of the apertures are included in the
full-wave Yrad simulations, resulting in good agreement between model and measurements.
By scaling the thickness of the full scale aperture to ? 1?op, we calculate the frequency
dependent Yrad with CST. The corresponding full scale Yrad is identical to that of the
scaled aperture. If one conducts full scale multi-cavity RCM calculations using this ?thick
aperture? Yrad, the obtained impedance statistics match well with the scaled cases (solid
lines in Fig. 2.3 (b)).
We are also able to calculate the statistics of the magnitude of the induced voltage
delivered to a 50? load impedance attached to the last cavity in the 1D chain due to a
given input to the port on the first cavity. On the experimental side, the load induced
voltage VL is calculated from the measured impedance based on the analysis presented
in Refs. [58, 62]. M?ore specifically, the induced voltage signal at the load can be
expressed as |V | = | 2Pin|Z 2 2
2
p| |Z11|
L |, where ZZ = Z12ZL/Zeq and Z = Z ? 12 .Re(Z p11) Z eq 1122+ZL Z22+ZL
The quantities Pin is the source power, ZL is the load impedance, and Z11, Z12, Z22 are
22
the measured system impedance matrix elements. On the RCM theory side, the model-
generated induced voltage is calculated using Eqns. A.8 and A.9 in appendix A. The
input powers used in the experimental and model-generated VL are set to be 1W , and the
statistical analysis of VL are conducted throughout the entire frequency range. Thus, with
the input power known, one is able to compute the fluctuating induced voltage signal at
the RX port through the fluctuating quantities Zin and Zt (Eqns. A.8 and A.9). Despite
the differences in aperture thickness, good agreement between these two set-ups is found
for the induced voltage statistics shown in Fig. 2.4. The experimental results of the load
induced voltage PDFs for scaled and full scale cavity systems are represented by dotted
and dashed lines, respectively. The results show that such a scaling technique can be
very conveniently extended from single to multi-cavity systems, allowing investigation of
systems with a large number of cavities and more sophisticated connection topology. The
experimental results in Fig. 2.4 are also in good agreement with the model predictions
(solid lines).
The proposed theoretical formulation is not expected to work at the extreme high-
loss limits (? ? ?) due to the failure of the random plane wave hypothesis, which is
a prerequisite of the RCM. This breakdown can be expected when the estimated 3dB
width of a mode becomes comparable to the operating frequency (Q ? 1). However,
the model is valid for moderately large loss (? ? 1), and the impedance statistics
simplify to Gaussian distributions in this limit [27, 28]. The RCM theoretical formulation
can be applied to lower loss systems (? ? 1), but the stronger impedance fluctuations
of very low-loss systems poses great challenges for the acquisition of good statistical
ensemble data by either numerical or experimental methods [64, 65]. The formulation
23
will also require modifications of the cavity total admittance matrix when the inter-cavity
coupling strength is increased substantially. Non-zero off-diagonal components of the
Yrad matrix (Eqn. 2.1) must be determined when direct coupling between the input
and output channels of the cavity is prominent. Failure to include these off-diagonal
terms may contribute to the lack of detailed agreement between the model-generated and
experimental results in Fig. 2.3.
To sum-up, we have established the multi-enclosure RCM model to treat interconnected
single chaotic systems. The proposed model is able to cover different coupling details
between each sub-enclosure. We have conducted experimental tests and found good
agreement between the measured data and RCM simulation. Limitations of the RCM
formalism is also discussed.
24
Figure 2.4: The PDFs of induced voltage on a 50 Ohm load attached to the last cavity
(cav) of the full scale cavity cascade systems (3.95-5.85GHz) and its scaled counterparts
(75-110GHz), assuming 1W input on the single-mode port of the first cavity. The curves
corresponding to the 1, 2 and 3 cavity system are color-coded in blue, yellow and green,
respectively. The full scale system experimental (theoretical) data are shown at dashed
(solid) lines, and the scaled experiment data are shown as dotted lines.
25
Chapter 3: Hybrid Power Balance-RCM Method for Multi-enclosure Wave
Chaotic Systems
This work was a collaborative effort with the University of Nottingham group and
the Naval Research Laboratory group, and is published as Ref. [66]. The Random
Coupling Model (RCM) has been successfully applied to predicting the statistics of currents
and voltages at ports in complex electromagnetic (EM) enclosures operating in the short
wavelength limit [27, 28, 37, 58]. Recent studies have extended the RCM to systems of
multi-mode aperture-coupled enclosures. However, as the size (as measured in wavelengths)
of a coupling aperture grows, the coupling matrix used in the RCM increases as well, and
the computation becomes more complex and time consuming. A simple Power Balance
Model (PWB) can provide fast predictions for the averaged power density of waves inside
electrically-large systems for a wide range of cavity and coupling scenarios. However,
the important interference induced fluctuations of the wave field retained in the RCM
are absent in PWB. Here we aim to combine the best aspects of each model to create
a hybrid treatment and study the EM fields in coupled enclosure systems [66]. The
proposed hybrid approach provides both mean and fluctuation information of the EM
fields without the full computational complexity of coupled-cavity RCM. We compare
the hybrid model predictions with experiments on linear cascades of over-moded cavities.
26
Figure 3.1: Schematic illustration of the hybrid model applied to a 3-cavity cascade. In
the hybrid model, we use the PWB method to characterize the power flow from the first
cavity to the next to last cavity. The fluctuations in the final cavity are described using the
RCM method using the mean power flow values obtained from PWB as an input.
We find good agreement over a set of different loss parameters and for different coupling
strengths between cavities. The range of validity and applicability of the hybrid method
are tested and discussed.
3.1 Power Balance Method
The PWB method can be used to determine mean values of EM power flow and
energy in systems of coupled cavities [56, 57, 67, 68]. For a multi-enclosure problem,
the PWB method assumes uniform energy density in each enclosure, and solves for
the mean power density Si (i is the cavity index) by balancing the powers entering
and leaving each cavity. These power transfer rates are characterized in terms of area
cross-sections (?), such that the power transferred is ?Si. Various loss channels, such
as aperture/port leakage, cavity wall absorption and lossy objects inside the enclosure
are characterised through the corresponding cross sections ?o, ?w and ?obj , respectively
27
[57]. Constant power is injected into the coupled systems through sources in some or all
of the enclosures. The PWB method solves for a steady state solution when the time-
independent inputs and losses are made equal for each individual cavity in the system.
This requires reaching the so-called power balanced condition. The PWB method does
not contain phase information of the EM fields and thus does not describe fluctuations
due to interference. This can lead to an inaccurate prediction of enclosure power flow in
the case of small apertures.
3.2 Hybrid Method Overview
As stated in previous sections, RCM is able to produce accurate prediction of multi-
enclosure systems. The computational complexity of the RCM grows, however, with
the addition of large apertures connecting together multiple enclosures. In the RCM
an aperture is treated as a set of M correlated ports, the number of which scales with
the area of the aperture as measured in wavelengths squared. For example, a circular
shaped aperture whose diameter corresponds to four operating wavelengths allows ? 100
propagating modes, leading to M ? 100 ports in the RCM modelling of the inter-
connected cavities. A cavity with M ports is described by an (M ? M) matrix [16, 28,
37, 42]. When large apertures are present, connecting multiple cavities, the RCM model
can become cumbersome. First there is the need to calculate the matrix elements for each
aperture that describe the passage of waves through an aperture radiating into free space.
Second these matrices mush be combined with random matrices that give the statistical
fluctuations. Third, the RCM is a Monte-Carlo method in which the matrices simulating
28
the cavities are constructed for each realization, and many realizations may be required
to get accurate statistical results. Finally, the matrices must be connected together which
involves inverting the sub-matrices representing the sub-volumes inside a complex system
for each realization [16]. There is thus a need to develop a simple statistical method that
applies in cases where the apertures are larger than a wavelength, but small enough so that
the two enclosures connected by the aperture can be considered as two separate volumes.
Here, we introduce a hybrid approach that combines the PWB and RCM to generate
statistical predictions of the EM field for multi-cavity systems without the computational
complexity of a full RCM treatment. The hybrid approach is valid in cases where the
coupling between adjacent cavities is carried by many channels due to, for example, large
apertures as described above. Using the hybrid method, we apply the PWB method for
computing the average EM field intensities in each cavity and use RCM to predict the
fluctuations in the cavity of interest only. The modeling of multiple channels between
adjacent cavities is thus reduced to computing a scalar coupling coefficient in the PWB
model, often given through simple expressions involving the area of the aperture.
A coupled RCM model still needs to be applied where the number of connecting
channels between enclosures is small. Such small apertures act as ?bottlenecks? for the
wave dynamics and a full RCM treatment is necessary to characterize the fluctuations
correctly, see Appendix C and D. One may further reduce the computational cost of the
hybrid model using a simplified treatment as proposed in Appendix E. This simplified
version of the hybrid model does not require additional knowledge of the frequency-
dependent aperture admittance, which is usually obtained through full-wave simulations.
29
3.3 Hybrid Method Details Part I: RCM Treatment
Random Coupling Model provides an alternative method to describe the statistics
of the EM fields in a wide variety of complex systems. In contrast with PWB, the RCM
deals with both the mean and fluctuations of the cavity fields.
For coupled-cavity systems, the RCM multi-cavity treatment begins with the modeling
of the fluctuating impedance matrix of each individual cavity [16]. These matrices relate
the voltages and currents at the ports of a cavity. When cavities are connected the voltages
at the connecting ports are made equal and the connecting currents sum to zero. The input
port on the first cavity is excited with the known signal. This leads to a linear system of
equations that can be solved for all the voltages and currents. This system is resolved for
each realization of the cavity impedance matrices.
As introduced earlier, the RCM normalized impedance matrix can be written as
Z = iIm[Z ] +Re[Z ]1/2 ? ? ?Re[Z ]1/2.
cav avg avg RCM avg
System-specific information Zavg can be calculated by assuming the walls of the cavity
have been moved far from each port, and each port responds as if there were only outgoing
waves from the port. In the case of apertures as ports, the transverse electric and magnetic
fields in the aperture opening are expanded in a set of basis modes with amplitudes that
are treated as port voltages in the case of electric field, and port currents in the case of
magnetic field. The linear relation between these amplitudes is calculated for the case of
radiation into free space, and this becomes the average aperture admittance/impedance
30
matrix. Each mode in the aperture field representation is treated as a separate port in the
cavity matrix. Thus, the dimensions of the matrix grows rapidly with the addition of a
large aperture [42].
In the cavity cascade system, the above mentioned apertures (with M propagating
modes) are adopted as the connecting channel between neighbouring cavities. With
M connecting channels between cavities, the dimension of the above defined cavity
impedance matrix becomes M ? M . The matrix multiplications and inversions in the
calculation of RCM multi-cavity formulations [16] thus have complexity which grows as
O(M2.4) using common algorithms [69]. Thus for large M , the computational cost of
the RCM scales roughly as N ?M2.4, where N represents the number of cavities in the
system.
Here we propose a hybrid method for multi-cavity problems that combines both
PWB and RCM as shown schematically in Fig. 3.1. In an N?cavity cascade system
with multiple channel connections between adjacent cavities, we utilize the PWB to
characterize the mean flow of EM waves from the input port of the first cavity to the input
aperture of the last (N th) cavity. The RCM method is now applied to the last cavity and
the connected load at the single-mode output port of that cavity. Thus the hybrid method
combines the strengths of both methods: the fluctuations of the EM field will be captured
with RCM, and the computational cost is greatly reduced using the PWB method.
In the following, we discuss the hybrid model formulation in detail based on the
3-cavity system example shown in Fig. 3.1. We will first introduce the PWB treatment
to the first two cavities in the chain, followed by the modeling of the last cavity using
RCM. We will then discuss how to connect the PWB and RCM models at the aperture
31
plane between the last two cavities. With the model formulation introduced, we will look
into the validity of the hybrid model in section. A step-by-step protocol to apply the
PWB-RCM fusion to generic cavity systems is detailed in Appendix D.
3.4 Hybrid Method Details Part II: Detailed Single Cavity Treatments
PWB characterizes the flow of high-frequency EM waves inside a complex inter-
connected system based on the physical dimensions of the cavities, the cavity quality
factors Q, and the coupling cross sections ?o as well as the incident power Pin driving
the system [57, 67]. PWB then calculates the power densities of each individual cavity in
steady state. For the 3-cavity cascade system in Fig. 3.1, the PWB equations are
?? ?? ? ? ????? w1
+ ?o1 + ?o2 ??o2 0 ??? S????? 1???? ???
P
? in????? ??o2 ?w2 + ?o2 + ?o3 ??o3 ??????S2???
?
= ??? 0 ???? (3.1)
0 ??o3 ?w3 + ?o3 + ?o4 S3 0
where the ?wi?s refer to the wall loss cross section and Si is the power density of the ith
cavity, ?o1 and ?o4 are the cross section of the input and output ports, while ?o2 and ?o3
represent the aperture cross sections, see Fig. 3.1. The cross sections can be expressed
explicitly from known physical dimensions and cavity wall properties [57]. Pin is the
(assumed steady) incident power flow into the first cavity. In this example, we assume
that only the first cavity receives EM power from external sources. The balance between
the input and loss is achieved by solving Eq. (3.1) for the steady state power densities Si.
For example, the power balance condition of the first enclosure is expressed as (?w1 +
32
?o1 + ?o2) ? S1 = Pin + ?o2 ? S2. The LHS of this equation represent the loss channels
of the cavity, including the cavity wall loss and the leakage through the input port and
the aperture. The RHS describes the power fed into the cavity, consisting of the external
incident power and the power flow from the second cavity. The net power that flows into
the last cavity in the cascade is expressed as P2?3 = ?o3(S2 ? S3).
The last cavity is characterized by the RCM method. With the knowledge of the
cavity loss parameter ? and the port coupling details, the full cavity admittance matrix of
the last cavity can be expressed as [42]
( ) ( )0.5 ( )0.5
Y = i ? Im Y + Re Y ? ? ?Re Y .
cav rad rad rad
The quantity Y is a frequency-dependent block-diagonal matrix whose components are
rad
the radiation admittance matrices of all the ports and apertures of that cavity. We assume
no direct couplings between apertures because the direct line-of-sight effect is small in
the experimental set-up. Consider a cavity with two M?mode aperture connections,
the dimension of the corresponding matrix Y is 2M ? 2M . The matrix elements
rad
are complex functions of frequency in general and can be calculated using numerical
simulation tools. We use here the software package CST Studio to calculate the aperture
radiation admittance. The RCM normalized impedance ? is a detail-independent fluctuating
?kernel? of the total cavity admittance Y . The loss parameter ? describes the degree
cav
of lossyness of the system. With known ?, an ensemble of the normalized admittance ?
can be generated through random matrix Monte Carlo approaches [58]. Combining the
fluctuating ? and Y , an ensemble of ?dressed? single cavity admittance matrices for the
rad
33
final cavity can be generated. It is later shown in Appendix E that a substantial reduction
of the hybrid model computational cost is made possible using an aperture-admittance-
free treatment, at the price of reduced accuracy for longer cascade chains.
3.5 Hybrid Method Details Part III: Finalizing the Hybrid Model
We next connect the PWB and RCM treatments at the interface between the second
and the third cavity. As discussed in the previous section, the power flow into the third
cavity, P2?3, is calculated from the 3-cavity PWB calculation. Identical system set-ups
are utilized in the PWB and RCM treatments, including the operating frequency range,
the dimensions of the cavities, ports and apertures, and the loss of the cavity (achieved by
a simple analytical relationship between the RCM ? parameter and ?w). To transfer the
scalar power values P2?3 generated by PWB into an aperture voltage vector required for
RCM, we assign random voltages U o3 drawn from a zero mean, unit variance Gaussian
distribution for the M -mode aperture and calculate the random aperture power using
1
Po3 = Re(U
?
o3 ? Y ? U o3).2 rad
These randomly assigned aperture voltages U o3 are then normalized by the ratio P2?3/Po3
to match with the value calculated from PWB. Combined with the RCM generated cavity
admittance matrix Y , an ensemble of induced voltage values UL at the load on the lastcav
cavity is computed utilizing Eq. (A.9) in Ref. [16]. The power delivered to the load is
obtained using
1
P ?L = Re(UL ? YL ? UL),
2
34
where YL is the load admittance, taken to be 1/(50?) here.
With the formulation of the hybrid model now explained, here we discuss the
improvement in computational time and memory usage by replacing the full RCM multi-
cavity method with the hybrid model. For an N?cavity system connected through M?mode
apertures, the hybrid model requires only a fraction of 1/[(N ? 1)M2] of the memory
consumption as compared to the full RCM method, enabled by the reduced cavity impedance
matrix storage for the first N?1 cavities. Moreover, the computation time is also reduced
by eliminating the RCM modeling of the prior cavities. For example, the computation/CPU
time for the two-cavity cascade system with circular aperture connections reduces from
about 120s to 40s with application of the hybrid method (tests run on a typical workstation).
In addition, the computation time of the full RCM method scales with the total number of
cavities in the system N, while the hybrid method is insensitive to the further addition of
cavities. Aside from the difference in the CPU time, there is a finite loading time to move
the data (e.g. random matrices to model single cavities) from the hard drive to RAM.
Compared to the full RCM method, the hybrid model reduces the number of cavities
that need RCM modeling from N to 1, and thus decreases the data loading time. The
advantage of the method becomes more prominent for larger lengths of the cavity chain
N and larger apertures (having M modes).
3.6 Experimental Results and Discussion
We now conduct the PWB-RCM hybrid analysis for the multi-cavity experiment
and compare with measurements. The experiment setup is shown in Fig. 2.1, and photos
35
Figure 3.2: The PDFs of load induced voltage |UL| of 2- and 3-cavity experiments
(solid) and hybrid model calculations (dashed). The single cavity loss parameter is varied
from 9.7, 7.5, 5.7 and 1.7 from (a-d), respectively. The inset in (b) shows the 2- and
3-cavity experiment averaged induced voltage values ?|UL|? with respect to different loss
parameters ?. Multi-mode (? 100 modes) circular apertures are employed between the
cavities.
36
of the actual measurement setup can be found in Fig. 2.2. We consider 2- and 3-cavity
cascades with large (on the scale of the wavelength) apertures between the cavities and
single-mode connections to the load in the last cavity. The induced voltages at the load
|UL| are calculated from data using the methods reported in Refs. [58, 62], and these
experimental results are shown as solid lines in Fig. 3.2. An ensemble of induced load
voltages for the multi-cavity system is created by moving the mode stirrers in all cavities
between each measurement. The losses in the single cavities is altered by inserting equal
amounts of RF absorbers in each cavity. In addition, the hybrid PWB-RCM method is
used to calculate |UL| and the resulting distributions are shown as dotted lines in Fig. 3.2.
Good agreement between the measured and model generated results are observed over a
range of different total cavity numbers and single cavity loss values. Under varying cavity
loss conditions, the probability density function (PDF) of the induced load voltage |UL|
of the 3-cavity system has a lower mean value and smaller fluctuations compared to the
2-cavity system results. Going from two to three cavities will decrease the energy density
in the last cavity and thus the power delivered to the load. This difference between the
2- and 3-cavity |UL| becomes smaller when the single cavity loss decrease as can be seen
following Fig. 3.2 (a) to Fig. 3.2 (d), see also the inset in Fig. 3.2 (b).
For comparison purposes, we computed the induced load voltage UL value of the 3-
cavity system using just PWB for the entire 3-cavity system (Fig. 3.3). The computation
covers the entire operating bandwidth (3.75 to 5.50 GHz). The ensemble averaged values
of UL obtained from the PWB-RCM hybrid model computations (Fig. 3.2) are 0.19, 0.24,
0.32, 0.82 V for cases (a-d), respectively. With window averaging, the corresponding
averaged UL values from the PWB-only method are 0.17, 0.22, 0.30, 0.85 V. Thus, we
37
Figure 3.3: The 3-cavity case load induced voltage |UL| values computed with the PWB-
only model. The single cavity loss parameter is varied from 9.7, 7.5, 5.7 and 1.7 from
(a-d), respectively. Multi-mode (? 100 modes) circular apertures are employed between
the cavities.
38
find relatively good agreement between the mean induced voltage from the hybrid model
and the PWB-only model result.
The induced load voltage PDFs of the multi-cavity system can be generated solely
with the RCM formulation [16]. A comparison between the |UL| PDFs generated with
full RCM method, the hybrid method, and the experiments are shown in Fig. 3.4. Both
theoretical approaches are able to generate statistical ensembles which agree well with
the experimental results. We find that the RCM results (dashed lines) slightly outperform
the hybrid method (dotted lines) for the two-cavity case. However, the computation time
and storage cost of the full RCM method is N times larger than the hybrid method, where
N refers to the total number of connected cavities in the cascade.
3.7 Hybrid Model Limitations
The hybrid model is based on the assumption that the fluctuations in a given cavity
are independent of the fluctuations in adjacent cavities and thus of the fluctuations in the
power flowing between cavities (as a function of frequency, for example). We assess
the validity of these assumptions by analysing a multi-channel cascaded cavity system
in Appendix C. We study in particular the effect of the total number of effective cavity-
cavity coupling channels Mn between the nth and (n + 1)st cavity on the fluctuation
levels of the load-induced power PL connected to the final cavity. Since PL ? |UL|2,
the conclusions drawn from the power flow studies can also be applied to voltage-related
results as presented below in Section IV. Defining the load power fluctuation levels as the
ratio ? = ?PL?2 / ?P 2L?, where ?? ? ? ? represents averaging over an ensemble, we treat ? as
39
Figure 3.4: Comparison of induced load voltage statistics for the hybrid model (dotted),
RCM predictions (dashed) and experimental results (solid) for the case of 2 and 3 cavity
cascades with single-cavity loss parameter ? = 9.7, and circular multi-mode apertures
between the cavities. The frequency range is from 3.95 to 5.85 GHz.
40
a measure characterising the level of fluctuations of the power. Here, ? ? 1 and ? ? 1
refer to low and high fluctuations of the power values, respectively. In Appendix C, it is
shown that ?N
? ? (1 +M?1n ), (3.2)
n=1
where the product is over all the cavities in the cascade. If cavities n = 1 to n = N ? 1
have strong multi-channel connections due to, for example, large apertures with Mn ? 1
propagating modes at the operating frequency, then M?1? n ? 0 and the contributions to
(3.2) not including the coupling to the load, that is, N?1 ?1n=1 (1 +Mn ) ? 1; this holds for
the experiments described in Section IV with Mn ? 100 at the frequencies considered.
The quantity Mn is small when a single-mode waveguide connects the last cavity (N )
to the load. At the last cavity (N ) there is a single mode output port, Mn = 1 and
the 1 + M?1n = 2 which induces higher fluctuations at the load compared to the case
where all apertures are large. Similar small Mn situations appear when a ?bottleneck? is
introduces between cavities. It is therefore sensible to adopt RCM for just the last cavity
to capture the power fluctuations at the load, while it is sufficient to include the influence
of the intervening cavities with PWB only giving the required information about the mean
power flow.
41
Chapter 4: Application of Machine Learning to Wave Chaotic Systems
This work was a collaborative effort with the Naval Research Laboratory group and
is published as Ref. [70]. The wave properties of complex scattering systems that are
large compared to the wavelength, and show chaos in the classical limit, are extremely
sensitive to system details. A solution to the wave equation for a specific configuration
can change substantially under small perturbations. Due to this extreme sensitivity, it is
difficult to discern basic information about a complex system simply from scattering data
as a function of energy or frequency, at least by eye. In this work, we employ supervised
machine learning algorithms to reveal and classify hidden information about the complex
scattering system presented in the data [70]. As an example we are able to distinguish the
total number of connected cavities in a linear chain of weakly coupled lossy enclosures
from measured reflection data. A predictive machine learning algorithm for the future
states of a perturbed complex scattering system is also trained with a recurrent neural
network. Given a finite training data series, the reflection/transmission properties can be
forecast by the proposed algorithm.
42
4.1 Machine Learning Applications, Overview
Machine learning (ML) techniques have enjoyed intense research interest in recent
years [71, 72]. A ML algorithm treats any given task as a mathematical problem and
does not utilize knowledge of the specific physics underlying the data. Benefiting from
this ?model-free? nature, the ML algorithms find broad application in various sub-fields
of physics, such as the identification of phase transitions in condensed matter studies
[73, 74], the classification of multi-qubit states of trapped-ion experiments [75], the auto-
tuning of gate voltages in quantum dots system [76, 77], and the future state predictions
of spatio-temporal chaotic systems [78, 79, 80]. Although successfully applied in various
studies, one crucial drawback of the ML techniques is the trade-off between a successfully
trained program and the amount of data required during its training phase. However, the
problem of training data acquisition does not pose a great challenge to the analysis of
wave chaotic systems. In order to compare experimental data to theoretical predictions
based on the statistical methods mentioned in the previous paragraph, a large statistical
ensemble of measured data is required. This feature of wave chaotic system analysis is
suitable for building a successfully trained machine learning algorithm.
4.2 Multi-enclosure Experimental Setup
We study the transmission and reflection of electro-magnetic (EM) waves in wave
chaotic enclosures. Short-wavelength EM waves from 3.95-5.85 GHz are injected into
cavities of dimension 0.762?0.762?1.778m3 through WR178 single-mode waveguides.
43
Figure 4.1: Schematic of the experimental set-up. We measure the S-parameter of
a 3-cavity cascade system with a VNA. The cavities are connected through circular
apertures. Rotatable mode stirrers are employed in each cavity to generate different
system configurations. The scale of the operating wavelength ?op is shown as the bar.
The cavities are electrically large (? 104 modes at the operating frequency range) in order
to simulate real-life examples of wave chaotic settings such as rooms in buildings and
cabins in a ship. A series of individual cavities can be connected into a linear cascade
chain through apertures as shown in Fig. 4.1. The total number of connected cavities is
varied from 1 to 3. Independent mode stirrers are employed inside each cavity to create a
large ensemble of statistically distinct realizations of the system [38, 81, 82, 83]. Single-
mode ports are mounted on the first and last cavity in the cascade, shown as T(R)X in Fig.
4.1. The 2? 2 Scattering(S)-parameters of the entire cavity cascade system are measured
with a Vector Network Analyzer (VNA), and the 2 ? 2 impedance(Z)-parameters are
calculated. The S and Z parameters are connected through the bilinear transformation,
S = Z 0.5(Z + Z )(Z ? Z )?1Z ?0.50 0 0 0 , where Z0 is a diagonal matrix whose elements are
44
the characteristic impedance of the waveguide channels leading to the ports. All mode
stirrers are rotated simultaneously to ensure a low correlation between each transmission
measurement, and the S-parameter measurements are carried out when the mode stirrers
are stationary. A total number of 200 distinct realizations of the cavity cascade are created.
The degree of loss in each cavity is altered by placing RF absorber cones in each cavity.
The single cavity ?lossyness? is described by the RCM loss parameter ? which is defined
as the ratio of the 3-dB bandwidth of a resonance mode to the mean spacing between the
modes [31, 60], and is basically the number of overlapping modes at a given frequency.
The loss parameter can have values ranging from 0 (no loss) to infinity (extremely lossy).
With ? known, one is able to generate the statistical ensemble of system impedance using
RCM.
The real and imaginary parts of diagonal impedance Z11 measurements for a single
realization is shown in Fig. 4.2. The curves correspond to systems with a different number
of connected cavities, varying from 1 to 3. The cavities are nominally identical, in that
all cavities share the same physical dimension and a uniform single cavity loss parameter
value: ? = 9.7. It has been demonstrated both theoretically and experimentally that the
statistical properties of the diagonal impedance Z11 of a high-loss (? ? 1) cavity cascade
system remain virtually unchanged regardless of the total number of connected cavities
in the chain [33]. The real and imaginary parts of measured Z11 values of multi-cavity
systems are shown in Fig. 4.2 (a) and (d). Direct observation of the frequency dependent
Z11 results cannot provide useful information to classify the number of connected cavities,
since the curves are essentially on top of each other and effectively indistinguishable vs.
frequency. The detailed Re(Z11) and Im(Z11) from 5.605-5.62GHz are shown in Fig.
45
Figure 4.2: Multi-cavity diagonal impedances Z11 from 2-port measurements on either
1-cavity (blue), 2-cavity (red) or 3-cavity (green) cascade. The inset in (c) is the schematic
diagram of the Z-parameter measurement of 1-, 2- and 3-cavity cascade structures. (a) and
(d): the real and imaginary parts of Z11 from 3.95-5.85GHz for a single realization of the
system, whose detailed views are shown in (b) and (e), and the corresponding statistical
analysis over a large ensemble of Re(Z11) and Im(Z11) are plotted in (c) and (f). The
single cavity loss parameter is ? = 9.7.
46
Figure 4.3: Generalized (recurrent) neural network architecture with two hidden layers.
The NN consists of the one input layer, hidden layers and one output layer. The RNN is
built upon the NN structure with the addition of ?memory effects? at each hidden layer
(shown as the dashed arrows).
4.2 (b) and (e). Even with the detailed view, it is still hard to differentiate the Z11 curves
from either the mode density, the level of fluctuation, or the averaged Z11 value levels.
Statistical analysis of the Re(Z11) and Im(Z11) measurements are shown in 4.2 (c) and
(f). The Z11 PDFs of 1-, 2- and 3-cavity cascade are very similar to each other and difficult
to systematically distinguish. The Fourier transforms of the multi-cavity Z11 data show
similar time-domain responses. We first wish to see if an ML algorithm can be trained to
correctly distinguish the number of cavities in the cascade simply from raw data such as
that shown in Fig. 4.2. The second objective is to see if an ML algorithm can predict the
evolution of the S- (or Z-) parameters as the cavity cascade is systematically perturbed.
47
4.3 ML Model I: Neural Network
The neural network (NN) and its modifications, inspired by the electric signal
propagation mechanism of brain neuron cells, are widely applied in various fields, such as
speech identification, pattern recognition and picture captioning [84]. A neural network
is one of the supervised learning algorithms whose goal is to utilize the given training
set information and establish a generic method to assign labels to the testing sets. As
a typical example of classifying cat and dog pictures, a supervised learning algorithm
refines its classification rule using thousands of images of correctly labeled cat and dog
images. More specifically, the algorithm identifies the label of a given picture through
the images pixel values. A typical structure of a NN is shown in Fig. 4.3. We will
next briefly introduce the operating principle of the NN with the example of an m?class
picture classification task (m = 2 for the cat and dog classification example).
In each iteration, the information of one picture from the training set is prepared into
a column vector and used as the input of the NN algorithm. For example, a picture of L?L
pixels is reshaped into a (L2, 1) vector x by simply concatenating the L? 1 lines L times.
The input vector is further modified as x? = x?<x> in order to improve the performance of
?x
the algorithm, where < x > and ?x are the mean and standard derivation of x. After the
input preparation, the ith input x?(i) is passed on to the first (represented in the subscripts)
hidden layer (i)h1 (with units) through
(i) (i) ? (i) (i) (i) (i)n1 h1 = ?(W1 x? +b1 ), where W1 and b1 are
a (n1, L2) coupling matrix and a (n1, 1) bias vector, respectively. The non-linear function
?(?) adds to the complexity of the NN and further expands the network?s functionality.
The first hidden layer vector (i)h1 is transferred to the second hidden layer through
(i)
h2 =
48
(i)
?(W2 ?
(i) (i)
h1 + b2 ), where
(i) and (i)W2 b2 are the coupling matrix and the bias of the second
hidden layer, and the same nonlinear function ?(?) is used. The algorithm labels the input
picture using (i) (i) ? (i) (i) . Similarly, (i) and (i)yNN = f [?(Wn+1 hn + bn+1)] Wn+1 bn+1 are the coupling
matrix and bias at the output coupling layer, and f is a normalization function which
transforms the value of (i) (i) (i)?(Wn+1 ? hn + bn+1) into a ?one-hot? label vector. The ?one-
hot? label is a one-column vector with m elements representing the m possible output
classes. For a given training data input, the normalization function f sets only one vector
element of (i)yNN to 1, which represents the label of the class for that particular data set,
and the other elements to zero. In the meantime, a cost function (i)f (i)cost = ||yNN ? y || is
calculated to describe the distance between the algorithm label (i)yNN and the known label
y(i) for that given input x(i). The goal of the training process is to minimize the value
of fcost through the so-called back propagation process, where the algorithm refines the
values of all coupling matrices W and biases b, utilizing the derivative values dfcost and
dW
dfcost . The back propagation process is repeated for all pictures in the training set. The
db
values of all network parameters are fixed after the training. In the testing phase, unseen
pictures are fed into the trained network with their machine predicted labels calculated,
while the back propagation process is not used. If the predicted labels agree well with the
correct labels, the classification NN algorithm is successfully trained.
The recurrent neural network (RNN) structure is based on the neural network with
the addition of a memory effect. As shown in Fig. 4.3, extra loops are employed in every
hidden layer. The RNN hidden layer h(i) has an additional influence from its previous
state h(i?1) through the h ? h coupling matrix (i)Wh?h. This feature grants a RNN the
ability to store information from previous iterations, and the potential to predict future
49
states. Differing from the NN algorithm, which is insensitive to the specific order of the
input information, a RNN algorithm requires that the inputs are fed into the algorithm
following a strict time-ordering or systematic evolution.
4.4 Task I: Multi-cavity System Classification Task
We utilize the neural network algorithm to classify different wave chaotic systems
based on raw scattering data. The goal is to train the machine to classify the total
number of cavities in the linear cascade chain by using the raw impedance (Z) matrix
measurements. As discussed in section I, the diagonal element Z11 of 1-, 2- and 3-cavity
systems is very difficult to distinguish by eye, and also shares nearly identical statistical
properties as shown in Fig. 4.2. In the 1-, 2- and 3-cavity cascade measurements, we
rotate the mode stirrer to 200 unique and distinct angles in order to generate a large
ensemble of measured data. The impedance measurement sweep is from 3.95-5.85 GHz
with 16001 frequency points. One-hot vectors of size (3, 1) are created as the correct
labels for all configurations. 80% of the total measurements (3? 200 = 600 realizations
from 1-, 2- and 3-cavity cascade cases) are used to train the algorithm, and 20% of the data
is reserved as a testing data set. The NN employs 4 hidden layers which have 25, 26, 33
and 18 units, respectively. The specific choice of the total number of layers and units per
layer can be varied. We use the Tensorflow package in the back propagation process. The
training and testing results are shown in Fig. 4.4. To speed up the algorithm, we select
subsets of size 2000 (Fig. 4.4 (a)) and 5000 (Fig. 4.4 (b)) uniformly chosen frequency
points out of the total 16001 Re(Z11) measurements as input vectors (the Im(Z11) input
50
algorithm also works just as well but the results are not shown). In both cases, the cost
function decreases sharply and the testing set accuracy reaches above 95% within 500
cost-function minimization iterations. We observe that training the algorithm with 5000
data points in the input vectors (Fig. 4.4 (b)) reaches a higher testing accuracy (98.3%) as
compared to the fewer data point case (95%). The testing accuracy reaches to 100% (Fig.
4.4 (c)) when all measured data points are fed into the NN algorithm. One can understand
the improvement of test performance from the fact that more information are delivered to
the algorithm with the increase of input vector size. The algorithm successfully identifies
the total number of connected cavities in the cascade utilizing only the knowledge of
the system diagonal impedance as a function of frequency. The result indicates that
the algorithm is able to detect and utilize ?hidden orders? embedded in the measured
data, completing a task that is hard to achieve with either visual inspection or by more
sophisticated analysis and statistical means [74].
We further explore the limitations and robustness of the developed multi-cavity
classification algorithm. We first test the algorithm?s ability to identify unseen measured
data. An algorithm is trained to successfully classify the Re(Z11) measured data from the
1- and 3-cavity cascade systems. The NN output has three classes: the 1- and 3-cavity
systems whose data are seen by the algorithm in the training phase, and an untrained
2-cavity system class. In the testing phase, we mix in the Re(Z11) data from 2-cavity
cascade measurements. We find that the algorithm classifies the 2-cavity data into the
all three categories with equal probabilities. This result indicates that the classification
algorithm applies strictly to the data drawn from the training ensemble. We next test the
algorithm?s tolerance level to noise. The algorithm is successfully trained to classify the
51
Figure 4.4: The NN classification algorithm performance for 1-, 2-, and 3-cavity Re(Z11)
data. (a-c): The cost function fcost (blue, left logarithmic axis) and test accuracy (red,
right linear axis) results are evaluated to 500 iterations of the algorithm. The dimension
of the input vector is 2000, 5000 and 16001 in (a), (b) and (c), respectively. (d): The test
accuracy results with and without additional Gaussian noise. The dimension of the input
vector is 5000. The signal to noise ratio values are -5, 5, 15 dB and the original data.
52
Re(Z11) data from 1-, 2-, and 3-cavity cascade measurements with minimal experimental
noise. In the testing phase, Gaussian noise are added to the Re(Z11) measurements to
achieve signal to noise ratios (SNR) of -5, 5 and 15 dB. The power level of the noise
Pn is calculated from SNR = 10 log10(Ps/Pn) based on the signal power level Ps. As
shown in Fig. 4.4 (d), the final test accuracy changed from 0.9 to 0.66 (saturating after
the 1000 iterations) with the decrease of SNR from +15 to -5 dB. The result shows that
the classification algorithm retains its ability to distinguish the cavity number when the
noise level is not too high.
The use of the frequency and realization-averaged Z11 values was investigated as a
possible alternative to the machine learning algorithm. The 1-, 2- and 3-cavity cascade
?Re(Z11)? are nearly identical (varying by 0.7% of the average value). The ML algorithm
fully retains the classification performance for the multi-cavity Re(Z11) data even after it
has been normalized by the average values. This demonstrates that the algorithm utilizes
details that are not easily summarized when making a high resolution distinction between
the three scattering scenarios based on statistics. With statistical methods, repetitive
independent measurements of the multi-cavity system can produce stable mean values
of ?Re(Z11)? to further distinguish the different cases. To obtain a good estimation
of ?Re(Z11)? and decrease the uncertainty, this approach involves a large number of
measurements of each different system in the data-gathering phase. In the prediction
phase, given an unknown system, a large amount of additional data is required to acquire
a good measurement of ?Re(Z11)? to perform classification. By contrast, a trained ML
algorithm can distinguish an unknown system with only one round of measurement.
Compared to the statistical method, the ML algorithm requires considerably less data
53
to produce a confident classification of the system.
4.5 Task II: Chaotic Cavity S-parameter Series Prediction
Inspired by the recent progress of predicting the future evolution of classically
chaotic systems with machine learning techniques [78, 79, 80], we utilize the RNN algorithm
to predict the evolution of a wave chaotic system under systematic perturbation. We
record the 2?2 S-parameter matrix of a single wave chaotic enclosure under a sequential
and systematic perturbation from a rotating mode stirrer. The dimension of the single
cavity is 0.038 ? 0.038 ? 0.089m3. EM waves from 75-110GHz are fed into the cavity
(? 104 modes at the operating frequency range) through single-mode WR10 waveguides
from Virginia Diodes VNA Extenders and the S-parameters are measured by a VNA.
The mode stirrer is rotated to 1000 unique angles with uniform step angle size. A full
rotation of the mode stirrer takes ? 120 steps. There exists some correlations between
the measurements. The task is to use a portion of the measured S11 and S21 data to
train the RNN algorithm, and then generate predictions of S21 by supplying additional
measured S11 data obtained later in the sequence. It is difficult to predict the future state
of a wave chaotic system based solely on its history information, since minute changes of
the system?s boundary condition can drastically affect the EM wave properties [58].
The RNN algorithm is implemented with the layer recurrent neural network function
from the Matlab Deep Learning Toolbox. Similar to the neural network training, we
first prepare the input vectors x? = x?<x> to improve the performance of the algorithm.
?x
In contrast with the input labeling process in the NN, one must prepare the input and
54
output vectors in the correct evolutionary order. We next define the desired network
parameters, such as the sizes of the hidden layer units and the method of back-propagation
optimization. In the chaotic cavity S-parameter prediction task, the input vectors are
assembled from measured |S11| at 50 frequency points with 175MHz spacing from 75-
110 GHz. The output vectors are from measured |S21| at 5 frequency points with 2.5GHz
spacing (? 1000 modes exist in this frequency interval). It is not shown here, but the
algorithm also works for the real/imaginary parts of the S-parameters as input and output.
The adopted RNN structure has 38 units in 1 hidden layer. The testing results from all
frequency points are shown in Fig. 4.5. The first 900 realization of the measured data are
used as the training set, and the testing begins with the 901st realization. Only measured
reflection information (|S11|) is fed to the algorithm, serving as the observer of the system
[78, 80], and we record the machine predicted transmission information (|S21|). As shown
in Fig. 4.5, we observe good agreement between the algorithm-generated and measured
transmission |S21| at least for the first 40 or so realizations for all 5 frequency points. The
prediction degrades beyond this point.
4.6 Discussion and Summary
The ML techniques are shown to successfully classify different wave chaotic systems
and to give predictions of future system states. Benefiting from its model-free property,
ML can be versatile for various tasks. However, the training process of ML methods
requires large amounts of measured or simulated data, and the tuning of algorithm hyper-
parameters, including the number of layers, the number of units per layer, etc. These
55
Figure 4.5: The |S21| prediction of a wave chaotic electromagnetic system future states
using measured |S11|. The measured |S21| data (colored solid lines) and algorithm
predictions (colored dashed lines) from 5 frequency points are compared for cavity
realizations 900 to 1000. The black lines show the absolute magnitude of the difference
between the measured data and the algorithm prediction.
56
choices must be guided by experience with the algorithms rather than knowledge of the
physical problem. The NN based wave chaotic classification task can be further developed
into a cavity loss parameter (?) detector, or an algorithm to distinguish different types of
perturbation of a wave chaotic system, for example. The application of RNN and other
techniques, such as the long short-term memory methods, the nonlinear autoregressive
neural network and reservoir computing, may give predictions for the future evolution
of the S-matrix when the boundary condition of the enclosure is subject to systematic
perturbation. The algorithm could be utilized to predict the future S-matrix of an evolving
system, thus allowing identification of coherent perfect absorption conditions of a wave
chaotic system, for example.
57
Chapter 5: Late-Time Time-Domain Random Coupling Model
5.1 Overview
The complex scattering of electromagnetic (EM) fields can be found in various
scenarios. In the short-wavelength limit, the wave dynamics of chaotic systems are
extremely sensitive to the geometrical details, which poses challenges for the modeling
of wave evolution. Benefiting from powerful numerical tools, accurate time-domain
modeling of the wave propagation can be simulated when the exact geometrical details
of a scattering system are known. However, the results obtained from deterministic
methods will not hold when perturbations are applied to the system, for example, a simple
change of system boundary condition. Under these circumstances, statistical studies are
advantageous for studying chaotic wave systems.
As we have introduced in the previous chapters, for frequency-domain quantities,
the Random Coupling Model (RCM) has found great success in characterizing the statistical
properties of various quantities, for example, the complete system scattering matrix and
the impedance matrix. The RCM method is based on the Random Plane Wave hypothesis,
which models the chaotic wave field as a superposition of many plane waves with random
phases, amplitudes, and directions. The coupling coefficients between the system ports
and the system modes are modeled as Gaussian random variables with unit variance.
58
Here, we aim to further expand the utilization of RCM from the frequency domain
to the time domain. By analogy with the frequency domain RCM, we construct the
time domain RCM (TD-RCM) by computing a set of system modes, each modeled by
a damped, driven, harmonic oscillator equation, using the finite-difference method. The
lossyness of the enclosure is modeled by the damping factor in each mode evolution
equation. System-specific quantities, such as the radiation resistance of the antennas, are
also included in the TD-RCM method. This chapter will first introduce the fundamental
equations of the TD-RCM method and then introduce the gigabox short-pulse experiments
used to verify the theory, followed by a comparison between experimental and TD-RCM
simulated results, and then finally I summarize our discussion. We would also note that,
as will be clear later, we are working towards a newer version of the TD-RCM method,
which includes the enclosure?s deterministic early-time (short orbit) behavior. Thus the
version of modeling discussed here should be taken as the late-time version of TD-RCM.
5.2 Details of TD-RCM Formulation
In this section, we will discuss the detailed formulation of the TD-RCM program.
The TD-RCM program models an enclosure with a volume V and supporting N eigenmodes
in the operating frequency bandwidth. The enclosure can have M ports, where each of
the ports can serve as a transmitting port or a receiving port. The type of load connected
at each of the ports can be different, and the degree of lossyness for each of the modes
can also be specified. The radiation resistance of the ports Rrad can be assigned with a
frequency dependence.
59
5.2.1 TD-RCM Equation of Motion
The time-domain Random Coupling Model (TD-RCM) computes the temporal evolution
of the cavity modes and port signals. The cavity modes (Un, n = 1, 2, ? ? ? , N ) are
represented by the voltages Un, where n is the mode index. The signal at port j is
represented by the port voltage Vj and port current Ij(j = 1, . . . ,M).
The cavity modes are described by the damped harmonic oscillator equation
M
d2 d ? d
U + v U + ?2n n n nUn = ?2 n cnjKnj Ij, (5.1)dt dt dt
j=1
where the port subscript j runs through all M system ports. The port voltages (Vj, j =
1, 2, ? ? ? ,M ) are modeled by a linear summation over all cavity modes,
?N
V = c ?j jnKjnUn ? Vjc. (5.2)
n=1
As such, the ports act to couple all the modes to each others.
2
In Eqn. 5.1, the damping term v = ??nn | | ? is derived from frequency domain RCM.?
The loss parameter is defined as ? = k2/(Q?k2). The loss parameter ? describes the
degree of lossyness of the system. With ? known, one is able to generate the statistical
ensemble of system impedance using RCM. The quantities k,Q,?k2 are the operating
wavenumber, enclosure quality factor, and the mean-mode-spacing, respectively; ?n and
??2n are the eigenmode frequency and mean-mode-spacing near ?n; cjn is the coupling
coefficient for mode n and port j (assumed to be Gaussian random variables); and Knj =
60
( ) [ ] 11
? ? 4 2Rrad,j(?n)??n 2 , where ? and ? are the (assumed uniform) permeability and
? ??n
permitivity of the enclosure volume material, and Rrad,j(?n) is the frequency domain
radiation resistance of port j at frequency ?n. In Eqn. 5.2, the random coupling coefficient
cjn(an)d m[ode voltages U]n are ?the same quantities as in Eqn. 5.1; K ?jn is defined as K ?jn =1 1
? ? 4 2Rrad,j(?n)??n 2 = ?Knj; Vjc represents the electrostatic mode contribution? ??n ?
of port j, defined as C dp Vjc = ?Ij; the term Cp is the port radiation reactance at lowdt
frequency (< 100MHz). We note that the order of subscript for the quantities Knj and
K ?jn is flipped. The reason for this setting is just putting the subscript that runs over the
summation in the second spot.
For an enclosure with N modes and M ports, we will have N sets of Eqns. 5.1 and M
sets of Eqns. 5.2. The unknown quantities that will be solved are the N Uns and M Vjs as a
function of time. Usually we will set the initial condition as Un(t = 0) = Vj(t = 0) = 0,
representing a quiet cavity. When the input pulses are turned on, the port voltage values
at the transmitting ports will have nonzero values as well as nonzero time derivatives.
As shown in the RHS of Eqn. 5.1, the change of port signals will drive cavity modes to
oscillate. The oscillating cavity modes will further contribute to the nonzero signals at
all ports (the RHS of Eqn. 5.2). With all equations and parameters given, it seems that
one is able to solve for all system mode voltages and port voltages using finite-difference
time-domain evolution methods. However, we still need to have an explicit relationship
between the port voltage and current signals, and this is related to the load attached at
each port.
61
Figure 5.1: Diagrams of the TD-RCM port load modeling. a. The linear load at port j
is represented with an impedance Zj . b. The schematic of a diode-loaded port. c. The
schematic of a multi-mode aperture. Note that the subscript ?in? corresponds to waves
following the arrow, going into the cavity.
5.2.2 System Port Treatments
In this subsection, we will discuss the treatment of two generic types of loads (Fig.
5.1): a linear impedance load and a nonlinear diode load. We first show the treatment of
ports with a linear load. The voltage and current at the port j can be expressed by
Vj = Vin + Vout,
(5.3)
Ij = (Vin ? Vout)/Zj,
where Vin (Vout) represents the incident (reflected) voltage wave coming into (out of) the
port from the cavity. The direction of wave is shown as arrows in Fig. 5.1. So that ?in?
corresponds to waves going into the cavity. The value of the load impedance Zj is usually
50?, which is represented by the load connected between the cavity and the ground in
Fig. 5.1 (a). Note that the TD-RCM allows the cavity to have multiple ports, shown
schematically with the lines connected to the left side of the cavity. For an input port, we
have I = Vin?Vout 2V= in?Vjj where the Vin is the waveform injected into the cavity. ForZj Zj
62
an output port, we have Ij = ?Vj/Zj since Vin = 0.
We may further combine the I-V relation of both input and output ports by the
following vector expression
I = Y V +D. (5.4)
Vector quantities are used here for simplicity: I and V are two M-by-1 column vectors for
all port currents and voltages, whose elements are just Ij and Vj . Here Y is a diagonal M-
by-M matrix where Yjj = 1/Zj , and D is a time-dependent vector with Dj = 2Vin/Zj .
Now we have obtained enough equations to solve for the port voltages. Note that the
values of port load impedance and the injected voltage signal values are known quantities,
and the above equation is only for one time step.
We next study the case with nonlinear output ports as illustrated in Fig. 5.1 (b). As
a true time-domain model, one of TD-RCM?s strengths is the ability to treat ports loaded
with nonlinear elements. As shown in Fig. 5.1 (b), the diode-loaded port is modeled with
a series resistance (R) and a shunt capacitance (C). The diode-port is modeled by the
following equations:
Vj = VD +RID, Ij = ?(IC + ID)
(5.5)
d eVD
IC = C Vj, ID = I0[e kT ? 1].
dt
The subscripts D and C represent the diode and capacitor. For convenience, we set ID = 0
for negative voltages. Eqn. 5.5 can be simplified to
d Vj VD
C Vj + + Ij = ? , (5.6)
dt R R
63
kT Vj ? VD
VD = ln(1 + ). (5.7)
e R I0
We will apply a finite difference expression for Eqn. 5.6. Port voltage and currents
are expressed by Vj(t) = [Vj(t + 1) + Vj(t)]/2, Ij(t) = [Ij(t + 1) + Ij(t)]/2, and
V ?j (t) = [Vj(t+1)?Vj(t)]/h, where h represents the time increment. The voltage across
the diode VD is computed by solving the transcendental equation 5.7. The Eqn. 5.6 can be
rewritten in a similar form as Eqn. 5.4. The matrix Y is still a diagonal M-by-M matrix
where Yjj = 1/R+2C/h, and the vector D is expressed by D = (2C ? 1j )V 2h R j?Ij? VR D.
5.2.3 Generating of System Eigenmodes with RMT
Note that one will need to provide a list of eigenmodes of the closed system to
run TD-RCM. It is possible to obtain all system eigenmodes with good accuracy using
numerical simulation tools. However, such a process would be time-consuming for an
electrically large system and not possible if the system?s geometrical details are not
accessible. Here we utilize the property of a chaotic system to simulate a list of system
modes by only knowing the macroscopic system information (i.e., volume, area) and the
operating frequency. The list of system eigenmode ?n is generated using Random Matrix
Theory (RMT). We start by computing the eigenvalues ?rmt of large random matrices.
The matrix size equals to the total number of system modes that we want to model, and
the matrix elements are drawn from a Gaussian distribution controlled by the symmetry
group of the system. The system eigenmode kn = ?n/c is computed from the relation
2 2
? = kc?knrmt 2 , where kc is the center frequency of the band that we want to model and?k
?k2 is the mean mode spacing. To sum-up, one first uses the total number of modes
64
Figure 5.2: The TD-RCM simulated RX port voltage signal from 0 to 1000ns. Five
realizations with distinct cavity eigenmode spectra and coupling are shown. The injected
waveform is a 5GHz center frequency and 5ns long sine pulse.
(N ) to generate N random matrix eigenvalues, and then uses the center frequency and
mean-mode-spacing to generate N system modes. Here we show several examples of the
TD-RCM simulated RX (receiver) port signals in Fig. 5.2. A diverse RX port response
can be clearly seen.
5.3 Gigabox Short-pulse Experiment
To test the applicability of the TD-RCM method, we conducted short pulse injection
experiments in a 3D enclosure (shown in Fig. 5.3). The gigabox is a metallic enclosure
with a volume of ? 1m3. A mode-stirrer, consisting of two thin metallic panels, is
65
Figure 5.3: Schematic of the short pulse time-domain experiment setup. The injected
waveform (left inset) is generated through IQ modulation using the arbitrary waveform
generator (AWG) and the programmable signal generator (PSG). The short pulse is
broadcast into the enclosure through the transmitting (TX) port (star). The data ensemble
is created by rotating a motorized mode-stirrer inside the enclosure. The receiver (RX)
port (star) induced voltage signal (right inset) is measured by the oscilloscope. A lab
computer is used for instrument control and data transmission/collection.
located inside the gigabox to randomize the ray trajectories. One may rotate the angle
of the mode-stirrer through many values to obtain an ensemble of data. To expedite the
measurement process, we installed a motor outside the gigabox to control the rotation of
the mode stirrer. We will conduct a 2-port time-domain experiment. One of the locations
serves as the input port, where a single center frequency modulated by a short pulse
envelope is injected into the gigabox. The injected signal will bounce around and excite
many gigabox eigenmodes within the pulse bandwidth. A non-zero induced voltage signal
will be measured at the receiver (RX) port. Typical input and output waveforms are shown
as insets in Fig. 5.3. Here we set the input wave as a single frequency sine wave with 20
ns duration and a center frequency fc. In the experiment we use a series of different fc
values ranging from 5 to 10 GHz. We rotate the mode-stirrer 200 times for each choice of
66
Figure 5.4: We plot the measured gigabox RX port voltage signals in (a), and the
corresponding TD-RCM simulation in (b). In each figure, the port signals from 20
realizations are overlaid to show the diversity of system response.
center frequency to create an ensemble of measurement systems and resulting transmitted
waveforms Vj(t), n = 1, 2, . . . , 200.
5.4 Statistical Analysis and Discussion
To test the performance of the proposed TD-RCM method, we have conducted
short-pulse injection experiments in the gigabox, and measured the induced voltage waveform
at a second part. Using realistic information for the experimental setup, namely the
gigabox mode-spacing at the input pulse center frequency fc, the loss parameter ?rcm,
and the measured I/O antenna radiation resistance, we can simulate our pulse injection
measurement with TD-RCM code at least as far as the statistical properties are concerned.
We show the simulated RX port signal and the data gathered from gigabox experiment in
Fig. 5.4 (b) and (a), respectively. Both data sets show a diverse response in terms of the
RX port pattern.
We studied the statistics of the maximum receiver (RX) port voltage in the time
67
Figure 5.5: The statistics of the maximum RX port voltage studied at different
frequencies merged in the gigabox. The input pulse is a single-frequency 20ns sine pulse.
The center frequency fc changes from 5GHz to 10GHz. Note that the input pulse has a
finite rise/fall time (1 ? 2ns). Data histograms (blue) are compared to the predictions of
the TD-RCM model (red).
domain. This quantity is defined as the maximum induced voltage signal at the port
after the input pulse is injected. The ensemble size of the TD-RCM is usually set as
1000 realizations, created by drawing a new list of system eigenmodes and the coupling
strength between the ports and the eigenmodes. As summarized in Figs. 5.5 and 5.6,
we have found relatively good agreement between the TD-RCM model prediction and
experimental data. Note that in Figs. 5.5, we find that the value of the RX port signal
decreases as the center frequency increases. This can be explained by the loss parameter
of the gigabox increases from ? ? 3 (5GHz) to ? ? 12 (10GHz).
68
Figure 5.6: Mean values of the receiver port maximum voltage PDFs shown in Fig. 5.5
as a function of the center frequency of the 20ns pulse in the Gigabox. The error bar is
one standard deviation.
69
5.5 Discussion and Summary
We studied the statistics of the maximum receiver (RX) port voltage in the time
domain. A relatively good agreement between the experimental data and the TD-RCM
simulation is observed. We have included the effect of SO by directly adding SO signal
to the TD-RCM simulated RX port voltages. To expand the applicability of TD-RCM
model, we have also worked out the formula for including nonlinear port loads. Note that
we have not tested these functionalities experimentally.
As mentioned earlier in the chapter introduction section, one of the limitations of
TD-RCM is an incomplete characterization of the early-time wave dynamics. In practice,
the proper setup of the reverberant field environment inside the chaotic enclosure usually
takes several Trt, where Trt represents the characteristic cavity time scale of the ray
completing one round-trip inside the enclosure. Before the setup of the reverberant fields,
the RX port-induced voltage signal will be mainly influenced by the short-orbit responses.
As shown in Fig. 5.7 (a), the experimental ensemble-averaged RX port signal shows a
sharp rise up at early time, and then starts to show an exponential decay trend at late-
time. The reason for the early-time RX signal rise-up comes from a highly repetitive and
stable prompt-response at the port location in many realizations of the cavity ensemble.
One can also find that the duration and amplitude of the rise-up increase is in proportion
to the injected pulse duration. However, as shown in Fig. 5.7 (b), our current TD-
RCM simulated port response does not have a similar short-time rise-up as seen in the
experimental case. Instead, one observes a smooth rise up for the ensemble-averaged port
signals.
70
Figure 5.7: The root-mean-square ensemble averaged RX port signal comparison. The
injected pulse is a 5GHz single frequency sine wave with a duration sweeping from 5 to
30ns. Experimental and TD-RCM simulation results are shown in a and b, respectively.
The insets are the blow-up of the early-time (before 100ns) curves.
There are several difficulties for including realistic early time response into TD-
RCM. Firstly, the TD-RCM has no spatial information. Early time responses are due to
short-orbits, and this requires knowing the interior geometry of the closed enclosure. This
in turn requires computing many short-orbit trajectories between the input and output
ports and determining how much of them survive the ensemble averaging process, and
considering the effect of the phase difference between the transmitted pulses travelling
along different trajectories [30]. This is similar to the Rician or Rayleigh fading scenario
in the wireless communication multi-path problems.
Another practical limitation is the trade-off between computing speed versus the
system mode-density. The TD-RCM models the N system eigenmodes. We can estimate
the required number of eigenmodes N = BW/?f(fc), where BW is the bandwidth of
the injected pulse, and ?f(fc) represents the system mean-mode-spacing at the pulse
center frequency fc. Here is an estimate of the scale of N in our gigabox short-pulse
experiment: there are ? 23.5k modes in the frequency range from 4.5 ? 5.5GHz, and
71
there are ? 93.8k modes in the frequency range from 9.5 ? 10.5GHz. Hence, the
computation speed of the gigabox short pulse experiments is limited by the total number
of modes. One may increase the pulse duration to decreases the BW . However, this trick
might be useful for proof-of-principle studies, and a method of systematically reducing
the total mode number is desired.
To sum up, the TD-RCM model computes the temporal evolution of the cavity
modes and port signals. A list of cavity modes is randomly generated using RMT by
knowing the realistic system mode-spacing. Each of the cavity modes is characterized
by a driven damped harmonic oscillator equation, where its equation parameters are
computed using the system loss parameter ?rcm and port radiation resistance Rrad. The
port signal is computed by summing up all cavity mode signals. A TD-RCM user will
need to input the basic enclosure information (mode-density and ?rcm), port information
(Rrad and voltage-current relationship for all nonlinear port loads), and the injected waveform.
The TD-RCM model will compute the temporal evolution of system modes and port
voltages. An ensemble can be created by changing the list of eigenmodes and the mode-
port random coupling coefficients. One can then study many statistical properties of the
time-dependent port voltages and eigenmode dynamics. As shown with the maximum RX
port voltage distribution study, the TD-RCM predictions show relatively good agreement
with the experimental results. Furthermore, we will be working towards improved versions
of TD-RCM to overcome the current model limitations. Topics to be addressed include
time-domain aspects of short orbits, as well as the inclusion of large multi-mode apertures.
72
Chapter 6: Wave-Based Reservoir Computing
6.1 Introduction
This work was a collaborative effort with the Maryland Wave Chaos group and the
preprint can be found in Ref. [85]. From recognizing images to becoming a Go world
champion, machine learning (ML) has changed our lives in myriad ways. These brain-
inspired techniques are able to imitate other processes without explicitly knowing their
underlying governing rules or equations of motion, and thus excel in starkly different
applications [84]. ML techniques and scientific exploration have become increasingly
intertwined [79, 86]. However, the computing speed and power efficiency of ML algorithms
largely depend on their computing substrates, which are usually computers based on von
Neumann architecture. Such machines only execute serial instructions and their speed is
limited by the ?von Neumann bottleneck?, which is induced by the serial transfer of data
between the memory and the processor [87, 88]. Secondly, CPU clock speed stagnation
has appeared as the Moore?s Law scaling continues to falter, and the Dennard Scaling
has broken down due to heat dissipation induced failures [89]. Thirdly, large amounts of
data are produced non-stop in both research and everyday activities that are hard to store,
both for physical and financial reasons. On the other hand, the ceaseless development of
ML has introduced massively parallel schemes which are clearly unfit for von Neumann
73
machines, whereas matching hardware to algorithms could lead to both fast and low-
power data analysis. Thus, researchers are prompted to explore alternative computing
paradigms, namely neuromorphic computing (NC), whose physical structures better fit
with the ML algorithms.
With a biologically-inspired physical structure, the NC approach aims to overcome
the challenges of current computers and to further realize the full potential of brain-like
ML algorithms [88, 90]. The construction of NC machines are deeply customized to
enable the parallel and low-cost execution of various ML schemes, ranging from the
spiking neural networks [91, 92], recurrent neural networks (RNN) [93], convolutional
neural networks [94] and reservoir computing (RC) systems [95, 96, 97, 98, 99]. The
application of NC can realize ?all-in-one? machines with fast data collection, storage
and analysis, and thus are suitable in various fields such as high-throughput experiments
[100], IoT and edge computing applications [101]. Compared to the electronic NC
approaches, photonics realizations promise an enhancement in computational speed and
power efficiency by utilizing speed-of-light information transmission and the intrinsic
parallelism of light [88, 102, 103, 104]. But current photonic approaches are still in
their infancy due to the challenges of creating and measuring high-dimensional physical
networks [93, 105]. For instance, the component number in a RNN device must grow
as N2 where N is the network size, just to maintain the parallel implementation of
matrices. Recently, the time-multiplexing concept has been implemented in delay systems
for creating high-dimensional RC-inspired hardware [106, 107]. However, this enhancement
of network scale directly results in a decrease of processing speed, as the virtual time-
delayed computing nodes must operate in a serial fashion [101, 108].
74
Here, we present an elegant solution to the ?speed vs. scale? dilemma with RC
inspired computing realized through a generic wave chaotic system. By exploiting the
chaotic dynamics of EM fields, our approach, in principle, could emulate arbitrarily
high-dimensional reservoir networks using only one physical structure. Promised by
the parallelism of light, the proposed architecture evolves autonomously and with high
speed. The computing power and versatility of the wave chaotic RC is demonstrated by
performing five different benchmark tests with only one device. We have also introduced
reservoir extension techniques (RETs) that allow the further enhancement of network
scale utilizing only modest measurement capabilities. More detailed content of the RC
setups are included in Appendix F.
6.2 Software Origin of RC
Reservoir computing is a general type of ML whose structure can be specified as
follows. Sequential (e.g., time-dependent) input variables, at time n denoted by the Nu-
dimensional vector column un (n = 1, 2, 3, ...), drive the evolution of a reservoir state
column vector r?n,
r?n+1 = f(r?n, un). (6.1)
In Eq. 6.1, it is assumed that f satisfies the ?echo-state? property: for an initial condition
r? = r?1 and a given sequence u1, u2, ..., un, ..., as n becomes large r?n becomes independent
of the initial condition r?1. The RC system output, denoted sn is taken to depend linearly
75
on a vector column rn (of dimension Nr) derived from the reservoir state r?n,
sn = W ? rn, rn = g(r?n) (6.2)
where W is a rectangular matrix. To promote effectiveness of the RC functionality,
Nr is typically taken to be much greater than Nu. The elements of the matrix W are
viewed as adjustable parameters whose values are chosen through a ?training? procedure,
whereby, based on training data consisting of many examples (un, s?n) of inputs un and the
corresponding desired resulting outputs s?n, W is determined by minimizing the deviations
overtime of s (the actual RC output) from s?. Heuristically, the basic idea of reservoir
computing is that, if the dimension Nr is large, and the Nr individual elements of the time-
dependent vector r evolve in diverse ways, then the best fit of s = Wr to s? will indeed be
a very good approximation to the time varying vector s?. The basic ML assumption is that,
following training, the outputs will continue to give a good approximation to the outputs
that would be desired for the post-training inputs. Operationally, this is promoted by the
application use of training regularization (typically, for RC, ridge regression) meant to
prevent over-fitting [78, 79]. In general, either the function f , or the function g, or both
should be nonlinear to allow the RC to perform a wide variety of nonlinear tasks [109].
We note that, by virtue of the relation (Eq. 6.2), the training of an RC reduced to a simple
linear regression [78, 79]. Thus the training of an RC can be very fast [110].
76
6.3 Hardware Realization of RC
As introduced above, the RC architecture is a type of RNN which consists of a
large number of nonlinear nodes connected through recurrent links. Compared with other
deep learning approaches, the key insight behind RC is its radically simplified training
process: only the output layer is optimised, by means of common regression methods
[78, 79]. In the meantime, RC demonstrates the same level of performance as compared
with other deep learning approaches [110]. The unique characteristics of RC make it an
ideal blueprint for neuromorphic computing applications [101, 108]. The training of a
hardware RC is completed in one shot by only optimising an out-coupling matrix, while
the entire physical structure remains un-disturbed and intact. In contrast, the training of
many inverse-designed [102, 103] or neural-network photonic devices [93, 111] relies
on the iterative tuning of all physical components, which requires extra time and energy
consumption. Moreover, the re-programming of the above devices is either impossible
[102, 103, 111] or involves a complete re-training process [93], while RC devices require
merely a simple switch of the output coupling matrix. RC-inspired hardware enjoys a
high training efficiency and task flexibility among different NC genres, and thus attracts
a great deal of research interest [101]. Earlier hardware RCs were built from distributed
physical computing nodes with high processing speed, but their network scales are limited
in practice [112]. Major efforts in the field were then directed to delay-based RCs which
circumvent such limitations with virtual temporal nodes. However, this approach is
palliative as its processing speed decreases with an increased network scale.
77
a Input b
Quater-bowtie Billiard 100
Read-out Random inputAWG Port signal
50
Oscillo-
0
scope
Amplifier -50
Analog -100
0 1000 2000 3000
negative positive Digital
Ez 10 cm Time (ps)
c Electric Dipole d e
120 1
Port signal FT
Top 0.860
Plate
0.6
Microwave 0 0.4
Diode
-60
0.2
Diode-port signal 4 GHz sine-wave
-120 0
Bottom Plate 0 500 1000 1500 2000 2500 0 2 4 6 8 10 12 14
Time (ps) Frequency (GHz)
Figure 6.1: Wave chaotic systems for reservoir computing. a, Schematic of the
experimental setup. The input information is transferred from a lab computer to the AWG
and broadcast into the reservoir. The amplified waveform is injected into the chaotic
enclosure through an electric dipole antenna. Several diode-loaded antennas are used to
probe the EM field, whose voltage signals are measured by the oscilloscope and further
transferred to the lab computer for training. b, Measured port signals (solid) under random
input stream (dashed). c, Schematics of the diode-loaded port. d, The dynamics of diode-
port voltages (red) under single frequency input at 4 GHz (black). e, The FT of the
diode-port signal shown in d.
6.4 Chaotic Wave-based RC
6.4.1 Operating Principles
The proposed wave chaotic RC solves the ?speed vs. scale? dilemma by exploiting
the parallelism of light and the richness of chaotic dynamics. Complex systems exist
in multiple physical and real-life scenarios. Examples include atomic nuclei, quantum
dots, and electromagnetic (EM) waves in electrically-large resonant systems. These
systems are chaotic in the short-wavelength limit where minute perturbations will lead
to a drastic change to the wave field configuration. Wave chaotic systems in the semi-
78
Port voltage (mV)
Amplitude (a.u.)
Port voltage (mV)
classical limit are successfully described with the Random Plane Wave hypothesis, where
the system wave properties can be simulated with a superposition of plane waves with
random amplitudes, phases and directions [37, 113]. Recently, researchers have exploited
the rich dynamics of chaotic systems to create novel phenomenon including wavefront
manipulations [114, 115, 116, 117] and ultra-fast random number generation [118]. Since
the software RCs are essentially constructed with complex networks, here we aim to
utilize such similarity and realize a new platform for RC-inspired hardware with nonlinear
wave chaotic systems.
Fig. 6.1a shows in detail the photonic signal processing and the reservoir dynamics
of the proposed RC. The digital input signals are transformed into analog waveforms by
the AWG (arbitrary waveform generator) for both the training and testing sets. Later, the
input pulse is amplified and injected into the chaotic quarter-bowtie billiard. The input
stream will excite hundreds of system eigenmodes, where each mode displays a distinct
eigenfunction pattern due to the chaotic nature of the enclosure in the short wavelength
limit. The EM field emerges as the real-time superposition of all system modes, sampled
at discrete locations in the wave chaotic system (as shown in Fig. 6.1b). Nonlinear
dynamics are introduced by the microwave diodes connected at the output antennas (Fig.
6.1c). The nonlinearity is demonstrated by the port signal distortion (Fig. 6.1d) and higher
harmonics responses (Fig. 6.1e). Thus, the complexity of the wave chaotic reservoir is
ensured by the ray-chaotic property and the nonlinear elements.
Here the induced voltage signals at the diode-loaded ports serve as the observable
reservoir nodes, and these are further transferred to a lab computer for off-line training
(Fig. 6.1a). In contrast to the existing hardware RCs which consist of physical or
79
virtual discrete-nodes, the dimension of the wave chaotic RC is arbitrarily high as its
reservoir layer is embodied by the continuum wave field inside the chaotic enclosure.
Therefore, the wave chaotic RC has brought the conventional RC structure to its optimal
stage: an arbitrarily high dimensional reservoir layer is made possible with a complex
scattering continuum. Constructing such a continuum reservoir is hardly practical in
software, but directly feasible in a photonic platform due to the true parallelism of light.
And high-speed evolution of the wave chaotic RC is always promised regardless of the
reservoir dimension. We note that similar continuum RC constructions have been studied
numerically and realized with a water-bucket or in the spatial light modulator systems
[119, 120, 121, 122]. However, the water-RC is nowhere suitable for application and the
evolution of the latter system requires significant assistance from von Neumann computers.
6.4.2 Reservoir Expansion Techniques
Challenges arise because the total number of measurement channels in the physical
RC is limited in practice. To overcome this limitation, we introduce novel Reservoir
Expansion Techniques (RETs) with the boundary condition perturbation method (Fig.
6.2a) and the frequency stirring method (Fig. 6.2b). Both methods aim to create a large
?combined? RC through the hybridization of many smaller ?unit? reservoirs. Fig. 6.2a
and b show snapshots of the electric field Ez landscape from time-domain simulations.
Small variations of system boundary condition, as well as the stimulating wavelengths, are
implemented by translating a metallic cylindrical perturber, and stretching or shrinking
the duration of the input pulse, respectively. Based on the intrinsic chaoticity, the application
80
a b u1
u u0 2
u3
c
Combined Reservoir
Observable node
Hidden node
Reservoir i-1 Reservoir i Reservoir i+1
Figure 6.2: Reservoir Expansion Techniques (RETs). a and b show the schematics
of the boundary condition perturbation method and the frequency-stirring method,
respectively. The first method may be realized by the translation of a metallic perturber
shown as the cylindrical rod. Under the same input waveform u0, uniquely different
evolutions of the EM fields inside the enclosure area are created by means of translating
the perturber to new locations. In b, the frequency stirring technique utilizes small
changes of center frequency of a given input waveform to create new wave field
configurations. In the experiment, the input wavelet is stretched/shrunken (shown as
u1 ? u3) in order to shift the center frequency. c, Schematics of the RET concept. The
combined reservoir consists of the observable nodes from various single realizations of
reservoirs.
81
Read-out
Read-out
of RETs results in a drastically different wave response, which is illustrated by the Ez
patterns (Fig. 6.2). New unit reservoirs each having distinct temporal dynamics are thus
created. Although a unit reservoir dimension might be limited in practice, its complex
dynamics are ensured by the high-dimensional continuum EM reservoir ?substrate?. As
illustrated in Fig. 6.2c, each unit continuum reservoir can be seen as a large reservoir
layer with many hidden nodes, and the port signals are just the observable nodes which
are sampled at discrete locations. Based on RETs, we have realized high-dimensional
combined RCs and demonstrated various tasks using a single chaotic enclosure with only
three measurement channels.
6.4.3 Experimental Set-up
The shape of the billiard is similar to a quarter of a bow-tie with an area of A =
0.115m2 [37, 123]. The characteristic length of the billiard is of A0.5 ? 35cm and
thus is electrically large compared to the operating wavelength (? 5cm). With a height
of d = 0.79cm, the billiard is considered to be quasi-2D because the electric field is
polarized in the z-direction for f < c/(2d) = 18.9GHz. The inputs are converted from
a digital time series to analog waveforms and then broadcasted by a high speed AWG
(Tektronix 70000A) with a sampling rate of 50 Gs/s. For a typical RC input waveform
(i.e., the Rossler observer task), we set the center frequency at fc ? 4GHz so that the
corresponding oscillation period is ? 0.25 ns. The BW of this Rossler observer task input
waveform is about 100 MHz where ? 3 resonator eigenmodes exist. We set the system
loss parameter at a range of ? ? 50, meaning that about 50 modes can be excited within
82
one mode?s 3dB BW. Thus we estimate that there are at most 150 modes excited in the
Rossler observer task. The input waveform is amplified by the RF-Lambda 2-18GHz
amplifier (RFLUPA0218G5). The output ports consist of SMA dipole antennas with high
switching speed diodes (Infineon BAS7004) connected between the center pin and the
billiard top plate (Fig. 6.1c) [123]. The voltage signals are measured with an oscilloscope
(Agilent DSO91304A) with a sampling rate of 40Gs/s. Ensembles of new uncorrelated
reservoirs are created by translating a metallic perturber (a 1.5cm radius, 0.75cm high
circular cylinder) and/or changing the signal carrier frequency from the AWG [27, 124].
For an input waveform centered at f0 = 4GHz, we translate the perturber by 1cm
(? 0.2?0) or shift the f0 by ?f = 100MHz (? 3 resonator eigenmodes) to create a new
reservoir. The scaling of f0 is equivalent to the stretching and shrinking of the enclosure
volume. Microwave absorbers are employed inside the billiard to alter the fading memory
of the system. For tasks running at ? GHz rates, the system decay time is in the ? ns
range which is common in most real-life settings [25].
6.4.4 Benchmark Test Results
Fig. 6.3 illustrates the computing ability and the effect of varying system parameters
of the wave chaotic RC. The observer task for the 3-component Ro?ssler chaotic attractor
(Fig. 6.3a-d) and the NARMA10 task (Fig. 6.3c and d) are executed with the wave chaotic
RC. See detailed information in Appendix F. For the observer task, the RC is expected
to produce continuous predictions for all of the variables of a complex dynamical system
(in our case, the Ro?ssler y and z component) based only on input from a subset of the
83
a b
2 Target Output 600 120
1 500
100
80
0 400
60
-1 300 40
-2 200 20
-3 0
0 5 10 15 20 25 500 750 1000 1250 1500
Time (T osc ) Tdecay (ps)
c 5 d
Target Output 10-1 Experiment 
4 Simulation
3
Ros
2 sler-
10-2 z
1
0 Rossler-y
-1
10-3
0 5 10 15 20 25 102 103
Time (T osc ) Reservoir size
e 0.7 f 160
Target Output
0.6 140
40
120 30
0.5
100
20
0.4 80
0.3 1060
0.2 40 0
0 20 40 60 80 100 500 750 1000 1250 1500
Bit index Tdecay (ps)
Figure 6.3: Benchmark tests and the RC performance map on system parameters.
a and c, Ro?ssler observer task testing set performance (blue) from a combined RC of 630
nodes. The true states are shown as red dashed lines in a, c and e. b, The deviation of
normalized mean square error for the Ro?ssler observer task: contour as a function of the
input duration Tosc and the system decay time Tdecay. The dashed boundary outlines the
optimal parameter island (error deviation < 20%) from where the combined 630-node
RC is built. e, NARMA10 task testing set performance at the optimal system parameters
from a 90-node RC, shown as the cross in f. f, The deviation of normalized mean square
error for the NARMA10 task.
84
Prediction Prediction Prediction
Tbit (ps) n.m.s error Tosc (ps)
error deviation (%) error deviation (%)
system (Ro?ssler x component). Like other ML methods, this prediction task is carried
out without knowledge of the equations of motion of the Ro?ssler system, but using only
a finite amount of data of all three components for training. The performance of the RC
pre?diction s
? = (x, y, z)?is evaluated with the normalized mean square defined as NMSE
= n[s
?(n) ? s(n)]2/ n s(n)2, where s is the target and n is the index of the testing
set data. With RET, a high-dimensional reservoir of 630 nodes is created and this RC
achieved an NMSE = (0.002, 0.017) for predictions of the Ro?ssler y and z components,
respectively. We visualize the high accuracy of the RC prediction by directly comparing
the testing set output to the target, shown as the blue and red curves in Fig. 6.3a and c,
respectively. As shown in Fig. 6.3a, the predicted Ro?ssler y component is practically
overlapped with the target. The predicted Ro?ssler z component covers the target curve
faithfully, where deviations only appear systematically at low-amplitude points due to
noise. For the NARMA10 task, the input stream is a random series drawn from the
interval [0, 0.5], and the goal is to infer a computed series based on the random input bits
run through a specific nonlinear, 10-state memory equation (Appendix F). The successful
execution of the NARMA10 posses a demanding requirement of both long-short term
memories in the ML system. With optimal system parameters (cross in Fig. 6.3f), a wave
chaotic RC with 90 nodes has achieved a record-setting performance of NMSE = 0.034
which is ? 5 times lower than reported hardware RCs [87]. Thus the complex memory
capacity of the wave chaotic RC is fully demonstrated.
Despite the phenomenal empirical successes of ML in many applications, its system
parameters and underlying mathematical mechanisms remain poorly understood. In contrast,
the governing parameters of the wave chaotic RC present a clear real-world physical
85
meaning: the input oscillation period Tosc (or bit duration Tbit) controls the RC input
speed, and the exponential decay time Tdecay of the reservoir dynamics determines the
fading-memory of the reservoir. The overall RC performance is thus a function of both
time-scales. As shown in Fig. 6.3b, the optimal reservoirs lie near the domain where
Tdecay ? 2.5Tosc, which indicates a fading-memory of ? 2.5 recent periods is preferred
for the Ro?ssler observer task. For the NARMA10 tests, the optimized performance island
arises when Tdecay ? 9Tbit (Fig. 6.3f). This empirical observation agrees nicely with the
setup of the NARMA10 task where each output bit is determined by its 10 previous inputs.
For the Ro?ssler observer task, a high-quality combined RC of 630 nodes is built from the
unit reservoir on the optimized island (the dashed area in Fig. 6.3b). As shown in Fig.
6.3d, the RC performance has enjoyed a substantial boost with the application of RET,
and will show further improvement through building a higher dimensional combined RC,
as demonstrated by the numerical tests.
Compared with the flexibility of electronic systems, neuromorphic photonics have
not yet developed into general-purpose applications [104]. Among different approaches,
the RC-inspired system enjoys a high task-flexibility benefited from its efficient training
procedure. Here we demonstrate the versatility of the wave chaotic RC by executing
benchmark tests with distinctly different goals utilizing just one physical device, namely
the nonlinear function simulator task, the He?non-map observer task, and the nonlinear
channel equalizer (NCE) task. As shown in Fig. 6.4a, the wave chaotic RC, serving as
an edge device, may switch between different tasks with a simple change of the output
coupling matrix. The wave chaotic RC is also able to incorporate with a wide range
of input speed. In the experiment, we have employed a Tosc ? 200ps for the function
86
simulator task and Tbit ? 80ps for the He?non-map observer task and NCE tasks. All
three tasks are successfully executed with high accuracy, shown by the overlapping RC
output and the target curves in Fig. 6.4b-d.
Besides experimental tests, we have also conducted numerical simulations of the
physical RC with the same shape billiard, identical degree of system loss, and realistic
model of the diodes (Appendix F). As shown in Figs. 6.3d and 6.4e-g, all of our experimental
cases are faithfully simulated with EM simulation numerical tools. The accurate simulation
capability of wave chaotic RC greatly benefits future RC optimization and follow-up
studies.
6.4.5 Discussions and Summary
As discussed earlier, the scale and the operation speed are the major measures of
the capability for both software and hardware ML. In practice, finding the ideal balance
between the two factors has been the center of the field of hardware RC, but current
realizations tend to focus only on one factor. Benefiting from the underlying ray chaos of
our wave system, we are able to construct an arbitrarily high-dimensional reservoir with
just one physical device. In the meantime, high-throughput computation is promised by
the autonomous evolution of the EM wave reservoir at the speed of light. Utilizing only
modest measurement resources (3 measuring channels), we introduce RETs and create
a large combined RC with over 600 nodes (Fig. 6.3d). The application of RETs would
introduce extra system operation time. However, the operation time can be linearly scaled
down by using higher frequency inputs, as well as enlarging the number of measurement
87
a
Cloud Analog
Task-dependent Digital
Input streams
Wave output couplers(microwave, Output 
optical, ...) chaotic RC streams
b 1 c 0.6 d 5
Target Output Target Output Target Output
3
0.5 0.3
1
0 0
-1
-0.5 -0.3
-3
-1 -0.6 -5
0 2 4 6 8 10 0 20 40 60 80 0 10 20 30 40 50
Time (T  ) Bit index Bit index osc
e f 0.25 g
0.3
Experiment Experiment Experiment 
10-3 Simulation 0.2 Simulation 0.25 Simulation
0.15 0.2
0.15
10-4 0.1
0.1
0.05 0.05
10-5 0 0
100 200 300 400 500 100 200 300 400 500 100 200 300 400 500
Reservoir size Reservoir size Reservoir size
Figure 6.4: Experimental demonstration of the task-versatility. a illustrates that a
single wave chaotic RC is able to execute different tasks with a simple switch of output
couplers, and thus serves as an ideal platform an edge device. The RC performances
(blue) for the function simulator, He?non-map observer, and the nonlinear channel
equalizer tasks are shown in b-d, respectively. The targets are shown as red dashed lines
in b-d. In e-g, we show the simulated (dashed) and experimental (solid) RC performance
as a function of the dimension of the reservoir.
88
n.m.s error Prediction
n.m.s error Prediction
Bit error rate Prediction
channels (see Appendix F). Under current experimental conditions, the combined wave
chaotic RC of 500 nodes performs on the order of 1014 OPS (operations per second)
which is close to the state-of-art supercomputers (1015 to 1017 OPS) [122]. The concept of
RC can be realized in photonic crystal cavity systems and operate at optical frequencies
[98, 120]. With an all-optical realization, the RC operation speed will enjoy a further
enhancement by changing the carrier frequency from from ? 1 GHz to ? 103 THz. Both
RETs can be executed in an all-electronic manner: the boundary condition perturbation
may be implemented with reconfigurable metasurfaces [116], and one can shift the input
waveform center frequency fc with AWG Sequencing which allows the fast switching
between different input waveforms.
Since our method applies for any reverberant environment for EM waves, the wave
chaotic RC may be realized in a broad range of wavelengths, spanning from WiFi waves
in real-life indoor environments [25] to on-chip nanophotonic chaotic cavities [114]. By
scaling up to the THz range, the wave chaotic RC will reach a much faster processing
speed with a highly compact device (Sup. Mat.). With state-of-art high speed cameras
[118], the number of observable nodes will also increase substantially by probing the
entire EM reservoir. Our RC paradigm can also be applied to other Random Plane Wave
systems such as multi-mode fibers.
Compared to existing hardware RCs, our chaotic RC shows a better robustness
against structural defects. The Random Plane Wave hypothesis further ensures that the RC
performance remains stable despite the specific location of the ports. The wave chaotic
RC also has an edge due to its high energy efficiency because the system is entirely
passive and the establishment of the reservoir layer is created by the EM radiation. All
89
of these characters are of great importance for practical utilities, such as edge computing
applications. Fig. 6.4 shows that a single wave chaotic RC can process different types
of tasks with a simple swap of output coupling matrix. As shown in Fig. 6.4a, the wave
chaotic RC, serving as an edge device, will carry out the in-situ data analysis and transfer
only the useful information back to the cloud. Because the device operates fully with
light, a wave chaotic RC combining with online output couplers allows high-throughput
analysis of raw inputs directly taken from various sensors. As an example, we have
demonstrated a 99.95% accuracy on the NCE task (Fig. 6.4d,g), which is a paradigm for
telecommunication applications. A stand-alone wave chaotic RC is also well-suited for
IoT applications, such as the on-site image classification and RF radar analysis for Smart
Cars. Its high-throughput computation is especially advantageous for applications which
demands fast and robust data analysis, for example the autonomous cars.
We note that the optimal ridge parameter region (Tosc ? [200, 500]ps) in the Fig.
6.3b contour does not extend to the lower and higher ends of Tosc. The values of Tosc
cannot be arbitrarily low considering the fundamental (lowest) frequency of the EM
enclosure (? 700 MHz). The higher end of Tosc is constrained by the practical limitation
of our equipment sampling rate. One is able to increase the mode density by letting the
RC operate at a higher speed. For example, the same bowtie billiard will host around 100
modes inside a 100 MHz bandwidth at 10 GHz. Utilizing an AWG with a higher sampling
rate, we should be able to increase the speed of the RC. The wave chaotic RC also faces
common practical limitations in nonlinear effects since the microwave diode requires a
sufficiently high field amplitude to display nonlinear properties [88, 123]. The existence
of noise, arising from the background noise or the limited analog-to-digital conversion
90
resolution, would also dampen the RC performance [78, 95]. We also note that the
analogue photonic systems are hard to compute with electronic devices in terms of high-
precision computation, but AI applications usually function nicely with low-precision
operation [104]. The efficacy of the output coupler may degrade as the RC scattering
properties loose fidelity over time. However, one can re-calibrate for the performance
drift with a simple update of the output coupling matrix because the training of the RC is
fast.
We have experimentally demonstrated a new physical platform for neuromorphic
computation utilizing the chaotic dynamics of waves, and found good agreement between
experiments and simulations. The computational power of the wave chaotic RC is clearly
presented by the successful execution of numerous benchmark tests, and the RC performance
is controlled by clearly defined and controllable physical properties. We present elegant
solutions to create a high-dimensional reservoir with minimal hardware requirements, a
substantial challenge for the field of NC. With all its unique qualities, we believe that the
wave chaotic RC will secure its place in future applications.
91
Chapter 7: Chaotic 1D Graphs with Photonic Crystal Waveguides
This work was a collaborative effort with the Maryland Wave Chaos group and
the preprint can be found in Ref. [125]. The statistical properties of wave chaotic
systems of varying dimensionalities and realizations have been studied extensively. These
systems are commonly characterized by the statistics of the short-range and long-range
eigenmode-spacing, and the one-point and two-point eigenfunction correlations. Here,
we propose photonic crystal (PC) defect waveguide graphs as a new physical setting for
chaotic graph studies. Photonic crystal graphs have two novel features, namely an unusual
dispersion relation for the propagating modes, and complex scattering properties of the
junctions and bends. Compared to traditional microwave cable graph systems, the PC-
graph bond lengths and node scattering properties are more easily modified. Moreover,
the electromagnetic fields in the entire graph, including both bonds and nodes, can be
more easily measured and analyzed in both experiments and numerical simulations. Here
we present numerically determined statistical properties of an ensemble of such PC-
graphs including both eigenfunction amplitude and eigenmode-spacing studies. Our proposed
system is compatible with silicon nanophotonic technology and opens chaotic graph
studies to a new community of researchers.
92
7.1 Introduction of PC Defect Waveguide
Wave-chaotic phenomena have been studied in various systems, ranging from 1D
graphs [117, 126, 127, 128, 129], 2D billiards [14, 37, 130, 131, 132, 133, 134] to 3D
enclosures [16, 58, 60, 66, 70, 116]. The statistical properties of many system quantities,
such as the closed system eigenvalues and the open system scattering/impedance matrices,
exhibit universal characteristics, which only depend on general symmetries (e.g., time-
reversal, symplectic) and the degree of system loss. Random Matrix Theory (RMT) has
enjoyed great success in describing the energy level spacing statistics of large nuclei [10,
135]. The statistics of complex systems that show time-reversal invariance (TRI) are
described by the Gaussian orthogonal ensemble (GOE) of random matrices, the statistics
of systems showing broken time-reversal invariance are described by the Gaussian unitary
ensemble (GUE), and the statistics of systems with TRI and antiunitary symmetry show
Gaussian symplectic ensemble (GSE) statistics [16, 18, 37]. It was later conjectured
by Bohigas, Giannoni, and Schmit that any system with chaotic dynamics in the classical
limit will also have wave properties whose statistics are governed by RMT associated with
one of these three ensembles [11]. In practice, identifying RMT statistical properties in
experimental data is difficult due to the presence of non-universal features that obscure the
universal fluctuations. The Random Coupling Model (RCM) has found great success in
characterizing the statistical properties of a variety of experimental systems by removing
the non-universal effects induced by port coupling and short-orbit effects [16, 27, 29, 30,
34, 37, 58, 60, 66, 123, 136].
Chaotic microwave graphs support complex scattering phenomena despite their
93
relatively simple structure [117, 126, 127, 128, 129]. A microwave graph structure can
be realized with coaxial cables as graph bonds connected at Tee-junctions. The graph
architecture allows for various useful circuit components (such as phase shifters and
attenuators) to be incorporated into the structure [117, 128, 129]. Recent studies show that
non-universal statistical features exist in chaotic graph systems, and these are hypothesized
to be caused by the non-zero reflection at the graph vertices. These reflections create
trapped modes that pollute the spectral statistics of the graph [129, 137, 138, 139]. However,
the experimental investigation of such phenomena is complicated by the fact that information
about the waves is accessible only at the nodes in cable graph systems, which are the
common choice for conducting graph experiments.
Here, we introduce an alternative type of chaotic graph built with photonic crystal
(PC) systems. Photonic crystals find extensive use in semiconductor-based nanophotonic
technology, where they act as low-loss waveguides for visible and infrared light [140].
The PC graph bonds are realized with defect waveguides and the nodes are formed by
the waveguide junctions. As TRI is preserved and spin-1/2 behavior is absent, chaotic
PC-graph systems are expected to fall into the GOE universality class of RMT. With
numerical simulation tools, we conduct a series of statistical tests of the chaotic photonic
crystal graph system including both eigenvalue and eigenfunction studies of closed graphs.
The advantages of using PC structures for quantum chaos studies are as follows. Compared
to other experimental graph systems (i.e., coaxial cables), one may probe the fields propagating
in PC waveguide graphs bonds easily. In addition, the bond lengths of PC-graphs can also
be altered by means of lithographic fabrication [140]. Secondly, the scattering properties
of PC-graph nodes can be engineered by changing the rod properties and geometry.
94
Figure 7.1: a. The open-plate view of the photonic crystal lattice with an L-shaped defect
waveguide region. b. The side view of the PC unit cell. A dielectric rod is sandwiched
between two PEC (perfect electric conductor) surfaces. The quantities h0 and d0 are
the height and diameter of the rod. c. The photonic band structure from a supercell
defect waveguide simulation. The waveguide modes appear in the bulk bandgap region
(delineated by the blue dashed lines). d. The E-field profile of one of the waveguide mode
solutions (the red circle in c). The black dashed line marks the center of the waveguide.
Moreover, PC waveguides and cavities are technologies that are widely used in the field of
integrated nanophotonics. To our knowledge, PC graphs have not been previously utilized
for quantum graph studies.
The paper is organized as follows. In Section II, we introduce the design details
of the proposed PC graph structure as well as its numerical implementation methods.
We present the chaotic closed-graph mode-spacing study in section III, and focus on
the discussion of closed-graph eigenfunction studies in section IV. Different methods of
conducting eigenfunction statistical studies are also discussed. We summarize the paper
and discuss the future applications of the proposed PC graphs in section V.
95
7.2 Photonic Crystal Graph
A photonic crystal system consists of a regular lattice of artificial atoms (or scatterers)
whose spacing is comparable to the operating wavelength [141, 142]. The material
properties and the geometrical details of the atoms are carefully designed in order to
achieve a specific functionality. A PC system is usually constructed as a 2D planar
structure, which makes it especially good for lithographic fabrication and photonics applications.
Importantly, a 2D PC can show a complete bulk bandgap in its electromagnetic excitation
spectrum [141, 142]. In PC-based devices, waveguides and cavities can be constructed for
example by making air defects (removing a certain number of the atoms) in the original
lattice. This creates guided propagating modes in the bulk bandgap, which ensures that
the modes are confined to the defect region. A variety of defect waveguides can be
realized by changing the atom properties [141]. Recent photonic topological insulator
(PTI) studies present an alternative form of PC waveguide using the interface between
two different topological domains [3, 143, 144, 145, 146]. Such PTI-based waveguides
can also be utilized as the bonds for realizing graphs. The interesting properties of the
topologically protected modes may further enhance the quantum graph studies. Here,
we will utilize the defect waveguide modes to build chaotic graph structures. Note that
these PC defect waveguide graphs are distinctly different from photonic crystal slabs
and billiards [134, 147, 148], and from coupled resonant dielectric cylinders [128, 149]
studied previously.
The construction of the chaotic PC graph starts by building a square lattice with
identical dielectric rods (the blue lattice in Fig. 7.1(a)). The lattice constant is a0 =
96
36.8mm. The dielectric rod lattice is located between two metallic surfaces and installed
in a vacuum background. The defect waveguide is created by simply removing one row
or column of the dielectric rods. The detailed shape of the dielectric rod is shown in
Fig. 7.1(b). The diameter of the dielectric rod is d0 = 13.2mm and the height is
h0 = 0.1a0 = 3.68mm. Because the PC lattice is thin in the vertical direction (z-
direction), the waveguide modes considered here are transverse magnetic (TM) polarized
(Ez ?= 0, Ex,y = 0, Hz = 0 and Hx,y ?= 0). The relative permitivity and permeability
of the dielectric rods are ?r = 11.56 and ?r = 1. We realize the proposed PC structure
numerically with COMSOL Multiphysics Software. The presence of the waveguide mode
is clearly demonstrated by the supercell photonic band structure (PBS) simulation (Fig.
7.1(c)). The supercell simulation model consists of a single column of PC lattice with
Floquet periodic boundaries on the two long sides. That is, fields on opposite long
sides differ by a phase factor ? = kxa0. The two short-end surfaces are assigned totally
absorbing boundary conditions. We remove the center rod to create the defect waveguide
region. The PBS simulation is conducted by computing the system eigenmodes while
varying the wavenumber kx in the range [??/a0, ?/a0]. As shown in Fig. 7.1(c), the
defect waveguide modes have emerged inside of the bulk bandgap region from 2.5 ?
3.6GHz. Here we see the first unique feature of PC graphs ? the PC waveguide modes
present a more complicated dispersion relationship between ? and kx (e.g., modes near
2.6GHz in Fig. 7.1(c)) as compared with coaxial cable (? = ckx). In Fig. 7.1(d), the |Ez|
profile of the entire simulation domain shows clearly that the mode solution is indeed a
guided wave because its amplitude is highly localized in the defect region.
Besides graph bonds, it is also important to characterize the scattering properties
97
Figure 7.2: a. The simulated S-parameters of the PC-graph right-angle junction (solid),
as well as the S-parameter from a 3.5mm uniform coaxial cable 2-port right-angle
connector (dashed), which serves as a reference. Inset: the schematics of the PC right-
angle junction and the 2-port coaxial cable junction. b. The simulated S-parameters of
the PC-graph T-junction (solid) and the S-parameters from a coaxial cable 3-port ?tee?
connector (dashed). Inset: the schematics of the PC T-junction and the 3-port coaxial
cable junction. All PC-junction S-parameters are frequency averaged with a 500MHz
window to remove the effect of spurious reflections.
of graph nodes [150]. The PC graph structure is realized by connecting multiple straight
defect waveguides with both right-angle and T-shaped junctions. We have characterized
the scattering matrix of both types of the junctions with COMSOL and CST Microwave
STUDIO. For the right-angle junctions, non-zero transmission is found only in the bulk
bandgap region. As shown in Fig. 7.2(a), The transmission values (S21) vary systematically
near but below 1 as a function of frequency, which qualifies the right-angle bends as
vertices with a non-trivial scattering matrix in a graph. The PC-graph T-junction presents
a complex scattering profile over the entire bandgap region. As shown in Fig. 7.2(b), the
magnitude of the reflection (transmission) coefficient deviates systematically from 1/3
(2/3) as a function of frequency. The scattering properties of coaxial cable right-angle
bends and tee-junctions are shown as red and blue dashed lines in Fig. 7.2, and they show
virtually no variation over the bandgap. These nontrivial junction scattering parameters
98
are a second unique feature of PC graphs, and can be expected to influence the properties
of the graph eigenmodes. We note that alternative methods of making waveguide bends,
or connecting waveguides, can be applied, for example by removing or adding dielectric
rods at the right-angle turn, in order to tune the transmission property of the waveguide
joints [151]. The engineering of the node reflection and transmission properties is beyond
the scope of this paper.
We simulate the closed PC graphs with the COMSOL eigenvalue solver. The
graph topology is that of a flattened tetrahedral graph having 14 straight segments and 13
junctions (including both 3-way junctions and right-angle junctions). The four exterior
faces of the parallel-plate PC structure (see Fig. 7.3(a)) are assigned totally absorbing
boundary conditions. The total length of the simulated graph is on the scale of ? 9m,
which hosts about 80 eigenmodes (within the bulk PC bandgap) in a typical realization.
We note that one is able to decrease the mode-spacing by enlarging the size of the PC
graph, and we choose the current system scale due to limited computational power. One
can create a statistical ensemble of chaotic PC graphs by changing the length of the bonds
for a given graph topology.
The E-field profile for a particular eigenmode solution of a graph (shown in Fig.
7.3(a)) is shown in Fig. 7.3(b). It is clear that the graph mode is localized to the defect
waveguide region and displays a longitudinal sinusoidal standing wave pattern. The mode
amplitude on a single bond is uniform but varies between different bonds (shown by the
color differences), where bonds are defined as straight waveguide regions between 3-way
junctions and right-angle bends.
99
Figure 7.3: a. Schematic diagram of the PC waveguide graph system. The graph bonds
are shown as the clear gray channels and the rectangular lattice is shown as a lattice
of dielectric rods (black circles). b. The simulated |Ez| profile of the graph system
eigenmode at 2.8GHz. c shows the statistics of the normalized mode spacing from graph
simulation (histogram). The theoretical predictions for the GOE, GUE and GSE systems
are shown in red, yellow and purple curves. d shows the spectral rigidity ?(L) of the
normalized spectrum from graph simulation (blue circles). The theoretical predictions for
the GOE and GUE systems are shown in red and yellow curves. e shows the statistics of
the consecutive mode spacing ratio r (histogram) and the theoretical predictions for the
GOE, GUE and GSE systems.
100
7.3 Eigenmode Spacing Statistical Analysis
We start by conducting the nearest neighbor mode spacing analysis of the proposed
PC graphs. An ensemble of 10 different graph realizations are studied numerically and
we obtained ? 800 eigenfrequency values from the ensemble. The graph topology ranges
from 13 ? 16 straight segments and 11 ? 14 junctions. We note that this ensemble of
graphs is conducted with the eigenvalue simulation in COMSOL. The simulation results
are used for mode-spacing studies here, as well as the eigenfunction statistics in Sec. IV.
The graph topology is kept fixed in the eigenfunction studies discussed below. Because
the topology and the total length Ltot of each graph realization is different, we normalized
the system eigenmodes solutions using the following method. We considered 10 different
graphs. For each graph we calculated eigenfrequencies in the range 2.8 GHz to 3.6 GHz
and placed the eigenfrequencies in ascending order. The eigenfrequencies were then
converted to eigen-wavenumbers (ki, i is the mode index) using the dispersion relation for
the PC-waveguide PBS shown in Fig. 7.1(c). This results in a tabulation of consecutive
wave-numbers. We then fit the mode-number index to the wave-number using a quadratic
fitting function n(k) = c 21k + c2k+ c3, where c1,2,3 are fitting parameters and n(k) is the
total number of modes below k [41]. For a typical graph we found c1 = 0.0054, c2 = 2.21
and c3 = ?55.28. The small value of c1 indicates the n(k) is nearly linear in k as expected
for chaotic graphs. We then evaluate the fitting functions at the calculated k-values for
each graph, creating a list of normalized eigenvalues ei = nfit(ki). The nearest neighbor
spacings are finally computed as si = ei+1 ? ei.
The distribution of the normalized mode spacing values of the entire ensemble is
101
shown as the histogram in Fig. 7.3(c). We have also included the approaximate theoretical
predictions of the mode-spacing statistics for the GOE, GUE, and GSE in the figure. The
three theoretical predictions are based on RMT [39, 152, 153], which is widely studied
as a seminal method of understanding the universal statistical properties of wave-chaotic
systems [11, 154]. As shown in Fig. 7.3(c), the distribution of the normalized mode-
spacing matches reasonably well with the theoretical prediction for the GOE universality
class. Good agreement between the graph nearest neighbor mode-spacing statistics and
the RMT theoretical prediction is also reported in various works on graphs, although long-
range statistical quantities tend to show non-universal behavior due to the trapped-mode
issue mentioned above [129, 137, 138, 139].
To test the long-range statistical properties of the eigenmodes, we compute the
spectral rigidity of the graph spectrum ?(L), where L is the length of a segment in
the normalized spectrum. Detailed discussion of the ?(L) computation and theoretical
predictions can be found in Refs. [135, 137, 138, 139]. Previous studies find that the long-
range statistics of quantum graphs usually shows some degree of agreement with the GOE
prediction in the range of roughly L ? [2, 5] [126, 137, 138]. The quantity ?(e, L) is the
least square deviation of n(e) from its best line fit inside the range [e, e + L]. Using
the spectrum from one graph realization, one can compute the ?(L) by averaging over
the ?(e, L) computed with non-overlapping segments of spectrum. Finally an ensemble
averaging process is conducted with the ?(L) data computed from all 10 graph realizations.
As shown in Fig. 7.3(d), the PC-graph spectrum actually shows a better agreement
with the GUE theory curve. Because the operating bandwidth of the proposed PC-graph
system is relatively narrow (Fig. 7.1 (c)), we believe the long-range statistics are corrupted
102
PC-Graph Random Order GOE GUE GSE
?r? 1.89 1.57 1.75 1.37 1.18
?r?? 0.51 0.56 0.54 0.60 0.68
Table 7.1: The summarized consecutive mode spacing ratios of the photonic crystal
graph (PC-Graph) system, the graph system with random ratio order, and the theoretical
predictions for the GOE, GUE and GSE systems [1].
by finite-size effects, unless the graph size is increased substantially.
In addition to the mode-spacing distribu(tion te)st, we note that the method of consecutive
mode spacing ratios r = sii and r?i = min r , 1i can also be adopted [1]. This type ofsi?1 ri
statistical study does not require unfolding the spectrum. For the PC graphs, the averaged
values of ?r? = 1.89 and ?r?? = 0.51, are closer to the GOE theoretical predictions than
the GUE and GSE predictions (Table. 7.1). We note that the high value of ?r? in the PC-
graph is related to the correlation between consecutive level-spacings. To eliminate the
correlation, we have randomized the order of original consecutive level spacings and then
computed the level spacing ratio values. The averaged values of the 1000 random draws
for the ratio are included in the Random Order column in Table. 7.1. We further present
the distribution of mode-spacing ratios r of the PC graphs and corresponding theoretical
predictions of the GOE, GUE and GSE systems in Fig. 7.3(e). The statistics of r match
reasonably well with the GOE prediction.
7.4 Eigenfunction Analysis
We next study the statistics of the closed PC graph eigenfunctions. The eigenfunction
statistics have been studied experimentally in 2D chaotic systems by probing the electromagnetic
(EM) standing wave field inside a microwave cavity [2, 14, 155, 156, 157, 158]. In those
103
studies, the experimental probability amplitude distribution and two-point correlation
function agree well with the random plane wave conjecture. Here, the wave properties
of the entire closed PC graph (nodes and bonds) can be faithfully simulated. We employ
the same eigenvalue simulation model used in the mode-spacing studies above. For a
graph system, the telegrapher?s equation is formally equivalent to the 1D Schro?dinger
equation, where the wavefunction ? is represented by the wave voltage. For thin parallel
plate waveguides, the wave voltage difference between the top and bottom metallic plates
is represented by the Ez value at the middle cutting plane (at the height of z = h0/2).
Because the PC graph bonds have a finite width, the bond eigenfunction will be evaluated
along a 1D line at the center of the waveguide (shown as the black dashed line in Fig.
7.1(d)). We next examine two PC graph eigenfunction characterization methods.
Method I: grid-wise representation. For each eigenmode of each graph realization,
we will use the entire set of the graph bond |Ez|2 vs. longitudinal position along the
center-line of the waveguide data points to represent the eigenfunction. The name ?grid-
wise? comes from the granular nature of this method, where the total number of eigenfunction
data points is inversely proportional to the computational grid size. We study the wavefunction
statistics of the PC graph by computing the distribution of the normalized probability
density v, which is defined as the square of the eigenfunction values vj = |?(rj)|2 =
?|Ez(r 2j)| ?Ltot
| |2? , where Ltot is the total length of the graph, rj and ?Lj are the locationj Ez(rj) ?Lj
and the grid size of the j?th grid point, |Ez(rj)| is the z-directed electric field magnitude
at the j?th grid point, and the summation in the denominator runs over every graph grid
point. In our simulation, the grid size ?Lj ? 0.05?op where ?op = 10cm is the operating
wavelength at 3GHz. The distribution of the probability density values is computed using
104
Figure 7.4: a. Probability amplitude distribution of photonic crystal defect waveguide
graph eigenmodes values (v) obtained from the grid-wise method (symbols connected by
blue line). The theoretical prediction for systems with GOE statistics are in red. b. shows
the statistics of the normalized bond amplitude value |a| (symbols connected by blue line)
and the Gaussian distribution (red dashed).
the data from all simulated eigenmode solutions (48 to be specific) for a single realization
of the graph, and the results are shown in Fig. 7.4(a), and discussed below.
Method II: bond-value representation. For each graph realization, the eigenfunction
of a mode is represented by a set of ?bond-values? Ez(bm) which is defined as the amplitude
of the standing-wave wavefunction on the graph bond bm. The quantity m is the index of
the bond and runs from 1 to 14. The standing wave on the bond is made up of two counter
propagation waves ?(x) = a eikx + a? e?ikxm m , where am is the wavefunction amplitude
at bond bm and x is the distance from a vertex along the bond. We first conduct a sine-fit
of the raw Ez(x) values on each bond, which yields the amplitude, and the value of |am|
is obtained as 1/2 of the amplitude value. T?he normalization process follows the same
method as in Ref. [159] which ensures that m L(bm)|a |2m = Ltot where L(bm) is the
length of bond bm. Here the distribution of the probability density values is computed
over 14 ? 78 data points, where 14 is the number of the bonds and 78 is the number
105
Figure 7.5: a. The inverse participation ratio (IPR) of graph eigenmodes obtained from
the bond-value method, computed based on the method IIa. The simulated |Ez| profile
of the red (black) circled mode in a is shown in c (d). b. The inverse participation ratio
(IPR) of graph eigenmodes computed based on the method IIb. The simulated |Ez| profile
of the red (black) circled mode in b is shown in e (f).
of eigenmode solutions from one realization of the graph, and the results for P (|a|) are
shown in Fig. 7.4(b).
We first discuss the eigenfunction statistics obtained using the grid-wise method in
Fig. 7.4(a). The theoretical prediction of the GOE systems, given by P (x) = ? 1 e?x/2,
2?x
is shown as a solid red line. We find that the simulated PC graph result matches reasonably
well with the GOE prediction. We notice that a mismatch between the graph data and the
GOE prediction exists at the small and large probability density values. One possibility is
that this method tends to over-count the appearance of medium-sized eigenfunction values
(similar to the systematic errors experienced in Ref. [14]). It may also indicate that the
data simply does not match the GOE prediction. In addition to the discrete eigenfunction
imaging method, we note that the eigenfunction statistics may also be tested through
resonance width distributions in transmission measurements (Porter-Thomas statistics)
106
[131, 160].
The random plane wave hypothesis underlying RMT eigenfunction treatments [39,
113] predicts that the complex coefficients am should have Gaussian random fluctuations
[159, 161]. Here we use Method II to study the distribution of the normalized bond-
values |a|, in Fig. 7.4(b), and find some degree of agreement with the random plane-wave
prediction. However, a clear deviation from Gaussian statistics is seen for low amplitudes,
similar to the deviation observed with Method I. Together, these results suggest that the
wavefunction statistics of this simple graph are not fully consistent with the random plane
wave hypothesis.
We next present the inverse participation ratio (IPR) computation based on the
bond-value method (Method II) in Fig. 7.5. IPR is a measure of the degree of localization
for a wave function [159, 161], and can be used to quantify the degree of deviation of
wavefunction statistics from the random plane wave hypothesis. Based on its IPR value,
the eigenfunction behavior varies between two limits, namely the maximum ergodic limit
where the wave function occupies each graph bond with equal chance, and the maximum
localization limit where the eigenmode is confined to only one bond. RMT predicts an
IPR value by assuming Gaussian random fluctuation of the eigenfunctions. Two different
IPR definitions are tested here. Method IIa follows the definition in Ref. [159], where
the IPR value for each graph mode is evaluated using the formula IPRIIa = ?|am|4?. As
shown in Fig. 7.5(a), the IPR values of the photonic crystal graph modes vary erratically
but lie between the maximum ergodic limit (IPRIIa = 1) and the RMT prediction limit
(IPRIIa = 2) [159]. In the maximum localization limit the IPRIIa = B, where B = 14
is the number of graph bonds, and the results are far from this limit. Method IIb follows
107
? [? ]2
the definition in Ref. [161] where IPRIIb = m |a?
4
m| / m |a?m|
2 . The quantity a?m
is the un-normalized bond wavefunction amplitude. Here, the graph IPR values also vary
erratically from mode to mode (Fig. 7.5(b)), but lie well below the maximum localization
limit (IPRIIb = 1) and closer to the RMT prediction limit (IPRIIb = 1/B = 0.07) [161].
The conclusions we draw from the above two methods are not exactly the same.
For Method IIa, two exemplary eigenfunction profiles are shown in Figs. 7.5(c) and d,
which correspond to the RMT and ergodic limits, respectively. One may directly spot
the different nature of these two modes based on their eigenfunction patterns, which is
a unique advantage of the photonic crystal graph system. For Method IIb, we present
the eigenfunction profile of two neighboring graph modes in Figs. 7.5(e) and (f). The
eigenmode in Fig. 7.5(e) has a larger value of IPRIIb and shows a strongly localized
distribution. That in Fig. 7.5(f) has a small value of IPRIIb and is more evenly distributed
over the bonds. The average value of IPR over all the modes are more meaningful in this
case, and we have ?IPRIIa? = 1.39 and ?IPRIIb? = 0.10 which are closer to the ergodic
and RMT limits, respectively, than to the localization limit. Previous work [126, 159] on
IPR on quantum graphs shows the surprising results that larger graphs, with the number
of vertices > 10, tend to show strong deviations from RMT predictions, while smaller
graphs show better agreement.
Taking into account the eigenfunction and IPR statistics, it would appear that PC
graphs are close to the random plane wave condition for tetrahedral graphs, but clear
systematic differences from RMT predictions remain. The complex vertex scattering
properties of right-angle bends and T-junctions, along with the unusual ?(k) dispersion
108
relation of PC defect waveguide modes, suggest that PC defect graphs may be very
effective for future wave chaotic studies.
To conclude, we have designed and simulated an alternative chaotic graph system
with photonic crystal defect waveguides that show an unusual dispersion relation. We
show that a series of statistical studies can be carried out on closed graphs, including
nearest-neighbor spacing statistics and eigenfunction statistics studies. Because both the
graph bonds and nodes can be probed directly in a convenient 2D geometry, one may
better analyze the non-universal features of chaotic graphs using a PC system. Because
both the defect region size and rods can be changed locally, one is able to change the
graph bond length easily by adding/removing rods, as well as change the graph node
scattering properties by engineering the rod property near the graph nodes. We note that
the PC 3-way and right-angle junctions have complicated scattering properties, and it
should be possible to adjust the degree of wave localization by engineering the scattering
properties of the PC waveguide junctions. These properties of the PC waveguides may
facilitate further studies of localization phenomena in graphs, for example the emergence
or suppression of trapped modes.
109
Chapter 8: Quantum Hall - Quantum Spin Hall Composite Photonic Topological
Insulator
This work was a collaborative effort with the Cornell group, and has been published
as Ref. [144]. I have also written a review article on microwave PTI (photonic topological
insulator) realizations, published as Ref. [145], but not included in this thesis. Photonic
topological systems, the electromagnetic analog of the topological materials in condensed
matter physics, create many opportunities to design optical devices with novel properties.
We present an experimental realization of the BMW (Bi-anisotropic Meta-Waveguide)
photonic system replicating both QH (quantum Hall) and QSH (quantum spin Hall) topological
insulating phases. With careful design, a composite QH-QSH photonic topological material
is created and experimentally shown to support reflection-free edgemodes, a heterogeneous
topological structure that is unprecedented in condensed matter physics. The effective
spin degree of freedom of such topologically protected modes determines their unique
pathways through these systems, free from backscattering and able to travel around sharp
corners. As an example of their novel properties, we experimentally demonstrate reflection-
less photonic devices including a 2-port isolator, a unique 3-port topological device, and
a full 4-port circulator based on composite QH and QSH structures.
110
8.1 Introduction to PTI Technology
During the past few years the realization of topological properties in photonic
systems has paved the way for many novel applications and the creation of new states
of light [162, 163, 164, 165, 166, 167, 168]. Translating concepts from condensed matter
physics into the optical domain has come from mapping electronic Hamiltonians, such
as the Kane-Mele Hamiltonian, into photonic platforms [165, 169, 170, 171, 172]. The
photonic topological insulators (PTIs) possess similar topological properties to their electronic
counterparts and allow light traveling through sharp corners without reflection [172]. The
observation of topologically protected surface waves (TPSWs) at microwave frequencies
has been reported with biased magneto-optical photonic crystals with broken time-reversal
invariance [173, 174, 175, 176, 177, 178, 179]. Due to the difficulty of achieving strong
magnetic effects at optical frequencies, a class of time-reversal invariant PTIs has been
proposed: these include coupled resonators [180, 181, 182, 183, 184] and bi-anisotropic
meta-waveguide (BMW) PTIs [143, 185, 186, 187, 188, 189, 190]. Recent studies revealed
new designs of time-reversal-invariant PTIs based on emulating the quantum valley Hall
(QVH) degree of freedom (DOF) in photonic systems [190, 191, 192, 193, 194, 195].
Elegant Floquet topological photonic designs are also reported [196, 197, 198, 199]. In
this manuscript, we focus on the experimental realization of quantum spin Hall (QSH)
and quantum Hall (QH) PTIs in the context of the BMW structure. We believe that this is
the first-ever realization of a time-reversal invariance breaking heterogeneous topological
system based on two different microscopic Hamiltonians, in either the photonic or electronic
domains.
111
The BMW approach to PTIs is unique in that it allows for integration of different
photonic topological paradigms based on the same unperturbed structure. Here we present
the first demonstration of a non-reciprocal BMW structure and integrate it with a reciprocal
PTI domain to create a completely new class of PTI materials. In such composite BMW
structures, a spin-momentum-locked TPSW can be spatially guided and filtered into desired
channels at the junction points between different topological domains with high efficiency
and free from backscattering [200]. We begin with the simulation and experimental
studies of the QH-PTIs followed by integration of QH-QSH PTIs. We are the first to
observe TPSWs at the interface of the QH-QSH composite system which effectively
serves as a 2-port PTI isolator. We then introduce experimental demonstrations of more
complex PTI structures, namely a unique topological Y-junction and a full 4-port circulator.
Our realized circulator structure has unique advantages over conventional photonic designs
where the guided modes are not topologically protected. The topologically trivial resonator-
based circulators generally suffer from narrow bandwidth and are susceptible to backscattering
due to impedance mismatches [201]. An alternative QH-PTI based circulator faces challenges
with respect to device size (tens of wavelengths) and the lack of a spin DOF [202]. On the
contrary, the topological protection and spin DOF ensure the propagation of guided waves
in heterogeneous BMW-PTI systems is backscatter-free, with higher directivity [200].
The utilization of topologically protected waves allows for a broad spectral operating
range, compact design, integration with other PTI paradigms and potential for high-power
operation.
112
Figure 8.1: (a) Quantum Hall (QH) unit cell model with h0 = 36.8mm, h1 = 1.45mm,
d0 = 12.7mm and d1 = 8.8mm. A QH rod consists of two ferrite disks (blue) straddling
a middle metallic rod, sandwiched by two metallic plates. (b) Photonic band structure for
a QH bulk region. The horizontal-axis is the kx vector mapping an excursion in the 1st
Brillouin zone and the vertical-axis is the eigenmode frequency. Under Hz = 1350Oe
the QH bulk bandgap (shaded area from 5.9 to 6.3 GHz) matches the quantum spin Hall
bandgap [3]. The inset offers an oblique view of a typical open-plate BMW hexagonal
lattice with lattice constant a0 = 36.8mm.
8.2 Design of QH-QSH BMW Structure
The BMW PTIs are based on an unperturbed structure of metallic rods arranged in
a hexagonal lattice sandwiched between two metal plates [3, 200, 203, 204, 205, 206].
The structure dimensions are carefully designed to ensure the TE and TM modes are
degenerate at the Dirac points with the same group velocity in the photonic band structure
(PBS). The two orbital DOF are emulated with left/right-hand circular polarizations and
the two synthetic spin DOF by the in-phase and out-of-phase linear combinations of the
TE and TM modes at the Dirac point [200]. This mode crossing at the Dirac point can
be destroyed by a broad range of spatial perturbations [187, 206] thus creating a band
gap, analogous to a bulk electronic insulator. In the present case, we create a QSH-
113
PTI by adding an air gap between the rods and one plate, resulting in a bi-anisotropic
response and an effective spin-orbit coupling (SOC) [3]. Ref. [200] Eqn. 6 shows
the Hamiltonian for the BMW system near the Dirac point, which has a same form
as the Kane-Mele Hamiltonian. More discussions on the BMW Hamiltonian can be
found in Eqn. 9.5. Alternatively one can create a QH-PTI by placing magneto-optical
materials in the gap. For QSH and QH domains, the spin-Chern numbers 2CSOCs,v =
?1 ? sgn(?SOC) and a global Chern number 2CTs,v = sgn(?T ), where s =?, ? is the
spin state label and v = K,K ?, are defined to quantify the topological properties of the
upper and lower bands [200]. The spin-orbit coupling ?SOC and time-reversal invariance
breaking ?T perturbations are defined as the overlap integrals of unperturbed modes?
field components within the perturbed volume in the QSH and QH regions, respectively.
An interface separating two PTIs with different Chern numbers will support edge states
that are spectrally localized in the bulk bandgap [162, 185, 206]. Both QSH and QVH
types of BMW have been experimentally realized [3, 190, 193], leaving the experimental
demonstration of the QH-BMW structure un-reported. Recent theoretical studies [200]
pave the way for realizing QH-PTIs in a BMW base by sandwiching biased gyromagnetic
materials on top and bottom of the metallic rod (Fig. 8.1(a)). In contrast with [200]
which characterized the ferrites using a dielectric permittivity ?r = 1 and a simplified
permeability tensor, here the QH region is re-designed in order to incorporate realistic
ferrite properties. Figure. 8.1 (a) and (b) show the bulk QH unit-cell model and its
PBS, respectively. This band diagram is obtained with first-principles electromagnetic
simulation using COMSOL Multiphysics. The ferrites are characterized by ?r = 14 and
the Polder permeability tensor [207] (see Appendix G). As shown in Fig. 8.1(b), under
114
an external bias magnetic field of Hz = 1350Oe, a bandgap opens throughout the first
Brillouin Zone, and matches closely with that of our existing QSH bulk structure [3, 206].
The degeneracy of TE and TM modes in the vicinity of the Dirac point is well-preserved
as shown in Fig. 8.1 (b).
8.3 Composite QH-QSH Systems
We conduct transmission measurements through the proposed QH bulk structure.
The gyromagnetic disks are made of commercial Yttrium Garnet ferrites which are magnetized
along the z axis using Nd-Fe-B permanent magnets located outside the structure on the top
and bottom plates. The QH region contains a 14 ? 5 ferrite-loaded lattice. Experimental
results are presented as transmission S-parameters in Fig. 8.2. The inset of Fig. 8.2
shows the measurement method: the two ports of the Vector Network Analyzer (VNA) are
connected to antennas at the center of a QH bulk structure and its perimeter. By switching
on the H-field to magnetize the ferrites, we observed an averaged 20 dB decrease of
transmission from port 1 to port 2 from 5.95 to 6.2 GHz which is slightly narrower
compared to the photonic bandgap predicted in Fig. 8.1(b).
We next construct a heterogeneous PTI structure with two topological phases: QH
and QSH PTIs. This composite PTI structure has an interface that will support topologically
protected edgemodes in the operating bandwidth (the shaded area in Fig. 8.1(b)) inside the
original bandgap frequencies for both QSH and QH bulk structures. The key requirement
is to maintain the spin degeneracy of the modes in the composite structure, which allows
edgemodes propagating without reflection [200].
115
Figure 8.2: |S21| transmission results through a bulk QH structure in the ? ? K(K ?)
direction for both magnetized and un-magnetized ferrites. With Hz on (magnetized
ferrites, red curve), we observe a decrease in transmission S-parameter at 5.9-6.3 GHz,
which is caused by the absence of bulk modes inside the bandgap region. Inset: schematic
of the experiment.
116
The composite QH-QSH PTI structure has 39 ? 5 BMW rods in both regions. We
use the VNA to conduct 2-port measurements for edgemode transmission along the QH-
QSH interface. In Fig. 8.3 (a)/(b), the left/right-going transmission SAB/SBA for the two
possible magnetic field biasing orientations is 35 dB higher than the transmission in the
other direction in the joint band gap frequency range. These experimental results show
that a QH-QSH interface can only support an edgemode with one synthetic spin direction.
The backscattering-free property of an edgemode is also examined in Appendix G.
The inversion of applied H-field in the z-direction will invert the sign of the normalized
bandgap ?T for QH-PTIs [200]. With the QSH-PTIs unchanged, the supported edgemode
at the QH-QSH interface will flip its spin, and more importantly invert its spin-locked
propagation direction. The topologically protected edgemodes ensure strong suppression
of the reflected flow of waves in the un-desired direction. The QH-QSH structure serves
as an effective 2-port PTI isolator which can be turned on/off by external H-field. The
isolator?s allowed propagating direction can be switched in real-time by the inversion of
the H-field.
8.4 Topological Y-junction
With the experimental realizations of both QSH-QSH [3] and QH-QSH waveguides,
we next propose a ?Y-junction? PTI structure which mixes three different topological
phases (Fig. 8.4(a)). The 3-port Y-junction possesses unique functions compared to
more traditional devices, such as a 3-port circulator. With its spin DOF and propagation
direction locked together, the edgemode transmission behavior inside the Y-junction can
117
Figure 8.3: (a) and (b) are the measured edgemode transmission along the QH-QSH
interface with opposite H-field directions applied to the QH region. The black-dot(red-
line) curves represent the right(left)-going waves. The insets show the locations of port A
and B along with the direction of the applied H-field. The white arrows mark the expected
edgemode propagating directions. We find that the edgemode propagation direction is
reversed when the biasing magnetic field direction is inverted.
be predicted theoretically. For instance, an edgemode with either spin up (??) or down
(??) can exist at interface ?ab?, with their propagating direction restricted to right and
left-going, respectively. Governed by the direction of the external H-field in the QH-
region, a ?Hz field will dictate that the Y-junction?s port 3 (2) becomes a forbidden source
(receiver). To illustrate a forbidden source (receiver) experimentally, note the decrease of
both S23 and S13 (S23 and S21) from 5.9-6.3 GHz shown in Fig. 8.4(b) (Fig. 8.4(b) and
(c)). Because the interface ?bc? can only support an up-spin mode ??, the signal injected at
port 3 can not travel left to port 1 or 2. In Fig. 8.4(c) the amplitude of S31 is 20 dB higher
than S21, indicating that a right-going wave transmitted from port 1 will travel to port 3
instead of port 2. This choice of path is dictated by the fact that such an up-spin TPSW can
only propagate along interface ?bc? but not ?bd?. For similar reasons, a down-spin wave
118
Figure 8.4: (a): The Y-junction design consist of three different topological domains:
two QSH domains (green and blue) and a QH region (red), creating three interfaces: ab,
bd and bc. The ??/? represents an edgemode with spin up/down. Purple/yellow arrows:
propagation directions of ??/? modes. Plot (b-d) are the transmission measurements of
the Y-junction with ?Hz in the QH region. The edgemode will follow the 2 ? 1 and the
1 ? 3 paths strictly without backscattering. Port 3 becomes a forbidden source as shown
in (b).
119
?? transmitted by port 2 will only go to port 1 as shown in?Fig. 8.4(d?). The ideal S-matrix
???0 1 0???
of the topological Y-junction can be written as S ? ?Topo Y = ???0 0 0????, which differs from
? 1 0 0
??
?
?0 1 0???
a traditional 3-port circulator S ?Circ = ???0 0 1???? which creates a 1 ? 3 ? 2 ? 1
1 0 0
circulation. The functionality of a topological Y-junction is similar to the active quasi-
circulator [208]. Inversion of the H-field direction will turn port 3 (2) into a forbidden
receiver (source), due to the change in propagating edgemode spin polarization at the two
QH-QSH interfaces (observed but not shown here) [200]. The topological Y-junction will
act as a building block for more sophisticated structures, such as the 4-port circulator.
The demonstrated Y-junction can serve as a versatile platform for many novel photonic
applications based on various combinations of topological domains. For example, one
can construct a backscattering-free photonic combiner with a QSH region and two QH
regions with opposite H-field directions that leads edgemodes from two paths to merge
onto the same third path with no backscattering. Furthermore, an inversion of biasing
H-fields in the QH domains will turn the photonic combiner into a spin filter that divides
photons with different synthetic spins into different paths, effectively creating a Stern-
Gerlach device for the spin DOF of light. Such junctions may find application in quantum
communication, Boolean networks models, and photonic devices based on manipulating
photons with different polarization states.
120
Figure 8.5: Experimental results of edge waves propagating in a topologically protected
4-port circulator. The QH (QSH with ?SOC > 0 and ?SOC < 0) domains are color-coded
as red (green and blue) for all insets. The applied H-field direction is inverted between
plots (a-b) and (c-d), creating a clockwise and counter-clockwise circulation around the
center QH island, respectively.
8.5 Topological 4-port Circulator
Novel non-reciprocal structure designs, like isolators and circulators, have been
proposed recently due to their important roles in photonic and microwave circuits [75,
209, 210, 211, 212, 213, 214, 215]. After the experimental demonstration of the Y-
junction structure, we are in position to integrate different PTIs into another practical
121
s?tructure, nam?ely a 4-port circulator. The ideal S-matrix of a 4-port circulator SCirc =???0 0 0 1???
??
? ?0 0 1 0?? ?
??
which creates a 1 ? 3 ? 2 ? 4 ? 1 circulation.
??
?
1 0 0 0????
0 1 0 0
In contrast to realizing optical circulators with only 2D gyromagnetic PTIs [201],
here we report the experimental realization of the combined PTI system with different
types of topological phases [172, 200]. Similar to the design in [200], the realized 4-
port BMW circulator consists of a center QH island and four surrounding QSH regions
with alternating ?SOC signs (Fig. 8.5). The directions of both QSH-QSH and QH-
QSH waveguides are chosen to be in the K(K ?) directions for optimized edgemode
propagation.
The operation of the BMW circulator is as follows. To start with, consider the
trajectory of the up-spin edgemode ?? which is launched from port 1 and propagating
right along the QSH(green)-QSH(blue) waveguide (inset of Fig. 8.5(a)). When the
edgemode arrives at the QSH-QH-QSH Y-junction, there exist two QH-QSH waveguides
(?red-blue? and ?red-green?) which support edgemodes with opposite spins. The topological
nature of these modes will lead the wave to only one of the two QH-QSH waveguides,
without backscattering [200]: at the first Y-junction, the spin up wave ?? will merge on to
the horizontal QH-QSH interface; at the second Y-junction, it will choose the QSH-QSH
interface over the vertical QH-QSH interface because that QH-QSH interface does not
support spin up edgemodes. As shown in Fig. 8.5(a), the measured transmission from
port 1 to 4 (S41) enjoys an averaged 15 dB boost compared to the competing vertical
122
path (port 1 to 3) in the bulk band gap from 5.9 to 6.3 GHz. Combined with the other
three guided edgemode paths, a clear circulation pattern is observed around the center
QH-island.
The inversion of the Hz direction will change the sign of ?T of the QH region
and further lead to the flip of spins of the propagating modes at all QH-QSH interfaces.
Considering again the same up-spin TPSW launched from port 1: instead of choosing
the horizontal QH-QSH interface, it will merge on to the vertical QH-QSH interface and
finally flow to port 3. As shown in Fig. 8.5(c), S31 (the transmission from port 1 to
port 3) is now higher than S41. To conclude: we observed the change of circulation
direction from counter-clockwise (Fig. 8.5 (a) and (b)) to clockwise (Fig. 8.5 (c) and
(d)) controlled by the external H-field biasing direction on the QH BMW region. The
isolation of the circulator design is 20 dB on average, the input power is +5 dBm limited
only by the VNA, and the frequency range is 5.9-6.3 GHz which can be tuned by scaling
the structural dimensions.
Note that the parallel QSH-QSH waveguides can have cross-coupling issues caused
by the finite-size limitation of the structure. The devices also face energy dissipation
due to the finite line-width of the ferrite, and enhanced insertion loss due to imperfect
coupling between the antenna and guided modes in the structure. Implementation of a
rotating dipole source will enhance the power delivery by directly and efficiently exciting
uni-directional edgemodes on the QSH-QSH interfaces [3].
123
Chapter 9: Chaotic Gaussian Symplecticc Ensemble Structures with Photonic
Topological Insulator Graphs
9.1 Overview
Unprecedented wave phenomena have been introduced with the creation of photonic
topological insulator (PTI) technologies [145]. One of them is the emulation of spin
in the photonic context, where a pair of topologically protected modes with opposite
effective spins have emerged. Mesoscopic condensed matter physics concepts, including
the quantum Hall, quantum spin Hall (QSH), and quantum valley Hall effects (QVH),
have been demonstrated in a series of landmark papers [162, 163, 164, 165, 166, 167,
193]. The bi-anisotropic metawaveguide (BMW) is a versatile platform for realizing
quantum Hall, quantum spin Hall, and quantum valley Hall effects for electromagnetic
waves. The remarkable transmission properties of PTI have excited a wide range of real-
life applications, for example, delay line devices, isolators, and circulators. However, the
spin-1/2 like degree of freedom of the topological modes is not as widely exploited.
Random Matrix Theory has enjoyed great success in describing the energy level
spacing statistics of large nuclei with various symmetries [10]. The statistics of complex
systems that show time-reversal invariance (TRI) are described by the Gaussian orthogonal
124
ensemble (GOE), the statistics of systems showing broken time-reversal invariance are
described by the Gaussian unitary ensemble (GUE), and the statistics of systems with
TRI and antiunitary symmetry show Gaussian symplectic ensemble (GSE) statistics [16,
18, 31, 37]. Both GOE and GUE statistics were observed in a wide variety of chaotic
microwave systems [18, 216]. Complex metric graphs made with standard coaxial cables
have been shown to have statistical properties described by the GOE and GUE random
matrix ensembles [126, 152]. However, due to the lack of a photonic analog of a spin-1/2
degree of freedom, the realization of a GSE system was not observed for a long time in
classical wave systems [41].
Here we aim to create a microwave simulator of GSE statistics by utilizing the
intrinsic physical properties of the QSH-BMW PTIs [3, 70, 206]. The bianisotropic effect
introduces an artificial gauge field and Berry connection for the two spin modes and gives
rise to a QSH-PTI. The system we consider here is a pyramidal graph with graph bonds
realized by PTI waveguides. Edge modes with both spin-up and spin-down indices can
travel through the entire graph so that there is an interference of waves at all of the nodes.
In addition to the clear observation of double (Kramers) degeneracy from photonic band
structure studies and edge mode transmission measurement, we will also prove that the
BMW Hamiltonian falls into the symplectic group and define the associated antiunitary
operator.
125
Figure 9.1: a. The schematic of the unit cell structure that define the QVH and QSH
regions. The quantities lv, gv, hv are the prism side length, the air gap distance, and the
prism height of one QVH rod. The quantities ds, hs are the diameter and height of the
QSH cylinder rod. b. The photonic band structure of the QSH-QVH interface waveguide.
Inset shows a schematic of the supercell simulation model. The red-dashed line marks
the interface between the QSH (left) and the QVH (right) region. c. The side-view of the
QSH-QVH structure (same structure as the b inset). The top and bottom metallic plates
are separated by h0. All QSH rods are attached to one of the metallic plates, and the QVH
rods are mounted on both plates, leaving an air gap.
126
9.2 Design of QVH-QSH Graphs
The BMW-PTI system can host domains with different topological orders inside
one composite structure. The first principle of the design of a composite QVH-QSH
structure is matching the bulk-bandgap of the single QVH and the QSH regions. There
have been several designs of QVH-QSH BMW structures [193, 217]. In principle, we
can adopt either QVH-QSH blueprint from Refs. [193, 217] to build the graph structure.
However, the design in Ref. [193] requires the QVH region to have metallic inclusions
that float between the top and bottom metal plates. Because we aim to build a graph
structure with many QVH (and QSH) rods, the proper fixation and alignment of these
?airborne? QVH rods is too tricky. The design in Ref. [217] uses QSH rods which
can change topological order through sliding the rod up and down mechanically. This
is a convenient design indeed. However, our focus is not on flipping the QSH domain
order but altering the shape of the QVH-QSH interface waveguide. This design goal
requires a uniform and robust mounting method for the QVH and QSH domains. Thus
we have designed a revised version of the QVH and QSH BMW region that operates
near 24.5GHz. As shown in Fig. 9.1 (a), the QVH region consists of two identical
metallic prisms with length lv = 4.30mm, height hv = 4.32mm, and the middle air
gap gv = 0.69mm. Each prism is firmly attached to one of the two metallic plates. The
QSH region is a metallic cylinder with diameter ds = 3.22mm and height hs = 8.3mm,
and it is attached to one of the two plates. Both regions are placed between two metallic
surfaces with a spacing of h0 = 9.33mm in the vertical direction. The lattice constant
of the triangular BMW lattice is a0 = h0 = 9.33mm. Fig. 9.1 (c) illustrates how
127
the QSH and QVH rods are mounted between the two metallic covering plates. We
examine the existence of edgemodes in the QVH-QSH waveguide through the so-called
supercell simulation, as shown in Fig. 9.1 (b). The supercell structure consists of a
slice of QSH and QVH rods (8 rods for each domain in our simulation), where the
Floquet boundary condition is applied to the two long sides to simulate an infinitely
long waveguide structure. It is clearly shown that the QVH-QSH waveguide can host
two counter-propagating edgemodes between 23.8 ? 24.8GHz. The two modes have
opposite spin and valley degrees of freedom (DOF).
The PTI graph can be thought of as the shape of a square pyramid realized in the
2D plane. Each bond of the graph is built with a QVH-QSH interface waveguide. The
graph will be electrically large and have an asymmetric shape so that the scattering will be
complex. The open-plate view of the QVH-QSH PTI-graph is shown in Fig. 9.2. All parts
are made with Aluminum 1070 alloy due to the high purity of Al (> 99.7%), a relatively
low percentage of iron (< 0.25%), and a reasonable cost. This choice of Aluminum alloy
is to enable a superconducting transition in the entire graph structure below a temperature
of 1.2K. Because of the relatively high operating frequency (? 24.5GHz), we can create
a graph structure with many lattice periods but moderate overall size. We note that
further scaling up of the operating frequency is limited by the bandwidth of the VNA,
and the further scale-up of the total graph size is limited by the room under the dilution
refrigerator mixing chamber plate (50cm diameter). To prevent the bending and flexure
of the covering plates, we have also installed several support rods (with diameter ds and
height h0) at the four corners and several spots in the horizontal plane. An ensemble of
PTI graphs can be realized by changing each of the graph bonds? shape and/or length.
128
Figure 9.2: The open-plate view of the PTI-graph structure. Shown are the top and
bottom metallic plates and the rods/prisms attached to each one of them. The top plate
has the top part of the QVH rods. On the bottom plate, we have installed the QSH rods as
well as the bottom part of the QVH rods. Several supporting rods are installed which will
spread out the weight of the top-plate evenly, and provide a fixed and uniform separation
between the plates, h0.
129
The conductive loss of the structure can be lowered by cooling down the graph using the
refrigerator.
9.3 Symplectic Group Analysis for BMW Systems
We are interested in using the QSH-QVH BMW system as a chaotic wave system
whose statistical properties fall in the universality class of the Gaussian symplectic ensemble
(GSE), which is appropriate for time-reversal invariant systems with a spin-1/2 degree of
freedom. An essential feature of a spin-1/2 system is an anti-unitary T-operator (T 2 =
?1). For a time-reversal invariant system, T 2 = ?1 leads to |?? and T |?? being two
different states with the same energy directly creates a Kramers degeneracy. An existing
microwave experiment on simulating GSE systems uses many elegant tricks to create
conditions that mimic a GSE graph [41]. As far as we know, no research group has
built an artificial physical system that replicates the full spin-1/2 analog. In Ref. [41],
the GSE statistics is emulated through a careful design graph which has two symmetry-
related sub-graphs. This approach is not general. Other than a limited number of works
[218, 219, 220], the derivation of the anti-unitary condition T 2 = ?1 is not widely
discussed for PTI systems and not discussed at all for the BMW system.
We started by proposing a guessed T-operator of the BMW system,
Tb = a ? ?x sy ?0K (9.1)
where a is a complex number coefficient with |a|2 = 1, K is the complex conjugation
operator (K2 = 1 and K i = ?iK). ?, s and ? are the Pauli matrices operating on
130
the space of the K and K ? valley, the spin states, and the Dirac degeneracy, respectively.
The reason for guessing Tb ? ?x and sy in Eqn. 9.1 comes from the valley and spin flip
associated with a time-reversal (TR) operation. The Kronecker product (?i?j ? ?i ? ?j)
shorthand is used.
A time reversal invariant system has [H(r), T ] = 0, or in the momentum representation
TH(k)T?1 = H(?k) where T represents the TR operator of the system. Here we will
examine TH(k)T?1 = H(?k). The effective Hamiltonians of the BMW lattice at the K
and K ? valleys are written as [206]:
H??K (?kx, ?ky) = vD s0(?kx?x + ?ky?y) + ?D?emsz ?z (9.2)
H??K?(?kx, ?ky) = vD s0(??kx?x + ?ky?y)? ?D?emsz ?z (9.3)
The quantities ?k = k?K refers to the distance to the Dirac point K, and ?k = (?kx, ?ky),
vD is the group velocity for both the degenerate TE and TM modes near the Dirac
point. The quantities s0 is identity matrix acting on the spin states, ?D is the Dirac
point frequency, and ?em is the bianisotropic terms which describes the coupling between
field components inside the air gap region. Note that the reason behind TE and TM
modes having the same group velocity is a process of carefully engineering the structure
geometry. The detailed introduction of the BMW system Hamiltonian theory is not the
goal here, and we kindly direct the readers to find more discussions in Refs. [200, 206].
With Tb = a ? ?x sy ?0K and K? K?1y = ??y, it is clear that
T H??(k)T?1 = T H??(?k , ?k )T?1 = H??b K b b K x y b K?(??kx,??k ) = H
??
y K?(?k) (9.4)
131
This result also applies to the K ?-Hamiltonian. One can also combine Eqns. 9.2 and 9.3
and write out the total Hamiltonian, which is [206]
H = vD s0(?z ?kx?x + ?0 ?ky?y) + ?D?em?z sz ?z, H? = ?? (9.5)
where ? is an 8-component spinor ? = [?K;? ]. Similarly we find that T H(k)T?1K? b b =
H(?k). We next study the effect of the TR operation on the spinor ?.
?? ???? ? ? ?0 1 ?K? ??K??
Tb? = a s ? ?? ?y?0K ? ?? ? = a sy?0K ?? ?? (9.6)
1 0 ?K? ?K
As shown in Eqn. 9.6, the two valley indices switch under Tb. The sy and ?0 will operate
on both the K and K ? valleys. For simplicity, we will only look at the K valley here.
The K-valley state ? ? ? ? R,? L,?K = [?K ;?K ] is a 4-component spinor with ?K = [?K ;?K ]
and ??K = [?
R,?
K ;?
L,?
K ] [206]. R(L) refers to the right (left) circular polarization. By
choosing a = i, we have Tb = i?xsy?0K and isy ?0K?K =
? ?? ? ? ?
? ? ? ? ? ? ? ? ??? 0 0 1 0 ?R,? ?R,?????? ??? ??? 0 1??? ???1 0??? ????? ??? ???? 0 0 0 1?
??????? K ? ? K ?L,?
K K ? ????????
?K ?
? ?
??? ???? ?
?
?L,?K ??
is ?y??0K = ? K = K ? = K ? ?
??K ?1 0 0 1 ?
?
K ????1 0 0 0???????R,????? ???K ???R,?K ????
0 ?1 0 0 ?L,? ??L,?K K
(9.7)
Thus the effective TRS operator Tb = i?xsy?0K meets its expectation: it results in
flipping both momentum (valley) and spin, as well as introduce a negative sign. However,
132
our job is not finished since we have not discussed the effect of the complex conjugate
operator K acting on the spinor ?. On the other hand, one might ask why the orbital
degree of freedom, represented by the circular polarization, is not discussed. Presumably
both spin and orbital angular momentum should flip sign under TR operation. It will be
clear later that the effect of K is to flip the orbital DOF.
The physical meaning of the K operator may be better understood in the original
TE and TM linear polarized basis AK [206], where
R/L,? R/L,?
?K = [1;?i; 1;?i], ?K =
[1;?i;?1;?i], and the superscripts R/L represent RCP/LCP. Note that the valley DOF is
linked to the circular polarity of the edgemode. Combining Eqns. 9.6 and 9.7, one may
sum-up the effect of Tb as:
R/L,? L/R,? R/L,? L/R,?
Tb?K,K? = ?K?,K , Tb?K,K? = ??K?,K . (9.8)
To sum-up, we first examined that the BMW Hamiltonian is invariant under Tb =
i?xsy ?0K. We then inspect the effect of Tb acting on the spinor ?. As shown in Eqn. 9.8,
we see that operator Tb flips the valley, spin and orbital degrees of freedom of the BMW
spin-mode. It is also clear that T 2b = (i?xsy ?0K)(i?xsy ?0K) = i?xs
2
y ?0 (?i)?x(?sy)?0K =
?1. Aside from deriving the effective Tb operator for the BMW Hamiltonian, we present
an alternative method of showing that the BMW belongs to the GSE group in Appendix
H.
133
9.4 Numerical Results
We first conduct studies of the eigenmodes of closed graphs using the numerical
simulation software CST (shown in Fig. 9.3). The eigenmode solver can obtain the
modal spectrum for closed systems with reasonable accuracy. It is clearly shown that,
inside the bulk bandgap frequencies (24 ? 24.8GHz as shown in Fig. 9.1 (b)), we
observe only doubly degenerated edgemodes. At the lower edge of the bulk bandgap
(? 23.5GHz), a co-existence of the degenerated edgemodes and the more continuous
bulkmodes is observed.
We then conduct studies of the eigenmodes of closed graphs using the numerical
simulation software HFSS. We used HFSS for this study because it is compatible with
the UMD High Performance Computing (HPC) cluster DeepThought2. For us, the ability
to run eigenmode simulations on an HPC is a must-have due to the PTI-graph models?
large amount of RAM consumption. We note that if RAM is not a limiting factor, using
COMSOL for eigenvalue simulations tends to be the best approach.
The simulated pyramidal graph has the area of ? 35 ? 60 photonic crystal lattice
constants, where ? 90 levels are sampled in the topological edgemode bandwidth. We
construct the graph structure with both QSH and QVH BMW-domains, and the two-fold
degeneracy (Kramers degeneracy) of system eigenmodes are observed from eigenvalue
simulations in HFSS. The normalized spacing statistics of the graph eigenfrequencies are
computed and compared to the theoretical predictions of GOE, GUE, and GSE systems.
The quantity si = ei+1 ? ei, where ei is the normalized (unfolded with the mean-mode-
spacing) graph eigenmode and the subscript i refers to the mode index. Details of the
134
Figure 9.3: The simulated eigenmode results of a PTI-graph realization. It is clearly
shown that the two-fold degeneracy is present at the bulk bandgap frequencies. At the
edge of the bandgap (? 23.5GHz), we observe a hybridization of the edgemode and bulk
modes. Inset shows the graph structure used in this numerical study, consisting of a QSH
region show by blue circles and a QVH region shown by grey triangles. The second inset
shows an enlarged region of the QVH-QSH boundary.
135
Figure 9.4: a The histogram approximation probability distribution of the PTI graph
nearest-neighbor mode-spacing s. The prediction of the GOE, GUE, and GSE groups are
plotted as solid lines for reference. b. The integrated probability function of s along with
the theory predictions of the three groups.
mode-spacing normalization process can be found in Chapter 7 Section 3. As shown in
Fig. 9.4 (a) and (b), the level spacing statistics of the pyramidal graph agree reasonably
well with the GSE statistics and show a distinctly different shape as compared to the
GOE and GUE systems. For small spacing values (s ? 0.5), we find that the graph
system deviates from the GSE class prediction.
9.5 Room-temperature Experimental Results
With the QVH-QSH region designed and verified by simulation tests, we fabricated
the QSH, QVH, supporting rods, and the top/bottom plates as shown in Fig. 9.2. Room
temperature S-parameter measurements are conducted with dipole antennas installed at
several of the graph nodes (shown as the black stars in Fig. 9.5). We use port 3 and
port 4 of the VNA to conduct S-parameter measurements. The reflection S-parameters
are shown in Fig. 9.6 (a), where several peaks can be spotted in the edgemode operating
frequency range (23.7 ? 24.8 GHz). From only one round of measurement, we have
136
Figure 9.5: Top view of the graph bottom plate, showing all QSH rods, and showing
the bottom half of the QVH rods. The white dashed lines are drawn to outline the
interface between the QVH and QSH domains, which represents the graph that hosts the
topologically-protected modes. The ports locations are marked by the black stars, along
with their numbering.
137
obtained a mean-mode-spacing of the graph of ?f = 0.09GHz, which is close to the
theoretical prediction ?f = c = 0.073GHz, where the quantity L
2L tot
= 2.05m is the
tot
total length of all graph bonds.
Our first goal is to identify the Kramers degeneracy expected in the graph modes.
To better examine the degenerate mode property of the data, we compute the Wigner
time-delay ? = ? i dW ln[det(S)] which helps us to identify scattering matrix poles andM df
zeros. The quantity M is the number of measurement channels, and S refers to the full
M ?M complex scattering matrix, measured as a function of frequency. The ?W around
groups of modes, as a function of frequency, is fit by Eqn. (7, 8) in Ref. [221]. For
convenience, we write down the ?W modeling equations here:
N [ ]
? 1
? ?n ? ? ? ?n + ?Re(?W ) = C (9.9)
M (f ? f )2n + (? ? ?)2n (f ? fn)2 + (? + ?)2n=1 n
?N [ ( )]1 1 1
Im(?W ) = ? (f ? fn) ?
M (f ? fn)2 + (?n ? ?)2 (f ? f )2n + (?n + ?)2n=1
(9.10)
where the quantities ?, fn and C are the uniform loss factor, the n?th system energy
level in frequency, and a constant term for the use of fitting to the data and accounts for
contributions to Re(?W ) from far-away resonances, respectively. The quantity N refers
to the total number of modes used during the fitting process. Since we are looking at
the Kramers degeneracy, we have N = 2 throughout the Chapter. Here we simplified
the fitting process by assuming the zeros and poles of the S-matrix, written as zn and
fn + i?n in Ref. [221], where n is the mode index, are just the complex conjugate
of each other. The quantities fn and ? ?n are the real and imaginary parts of the n th
138
pole. This assumption is good for the PTI-graph systems because there is no lumped loss
element. We note that the quality of the computed ?W curves is closely related to the
noise present in the S-parameter measurement. As advised by my colleague Lei Chen,
we adopted an intermediate frequency bandwidth of 10kHz on the VNA for a balance
between measurement speed and relatively low phase noise.
Based on Eqns. 9.9 and 9.10, we are able to fit the measured ?W (f) data. The
fitted curves are shown as the red solid lines in Fig. 9.6 (b) and (c). Note that we have
assumed that the two degenerate modes share the same uniform loss factor ?. It is clear
from the asymmetric ?W (f) curves in Fig. 9.6 (b) and (c) that the PTI-graph system
does host nearly degenerate modes. Our fitting results also back up this observation,
and high-quality fitting is achieved by assuming the existence of two nearly degenerate
modes. The best-fit parameters are (f1, f2) = (24.18, 24.20)GHz and (?1,?2, ?) =
(0.007, 0.0073, 0.0146)GHz for the data in Fig. 9.6 (b), and (f1, f2) = (24.45, 24.46)GHz
and (?1,?2, ?) = (0.0019, 0.005, 0.0098)GHz for Fig. 9.6 (c). One may argue that
the two nearby modes shown in Fig. 9.6 (b) and (c) might just be two closely located
un-related modes. However, the spacing between these two modes is on the order of
? 10MHz, which is much smaller than the ? 70MHz mean-mode spacing. Thus
we believe that we have identified some hints of 2-fold Kramers degeneracy in room-
temperature S-parameter experiments.
Note that we did not find 2-fold degeneracy characteristics for all modes evident
in the ?W (f) data. In fact, the two sets of curves shown in Fig. 9.6 (b) and (c) are the
only examples where the two degenerate modes poles/zeros have a large enough real-part
splitting. For other cases, the real part values of the degenerate mode pair might be almost
139
Figure 9.6: a. The measured reflection S-parameters as a function of frequency in the
bulk bandgap at Ports 3 and 4 of the composite QSH-QVH shown in Fig. 9.5 at room
temperature. The measured complex Wigner time delay curves near 24.20 and 24.46
GHz, shown in b and c, respectively. The fitting results are shown in solid red curves.
140
Figure 9.7: a and b show the real and imaginary parts of the ?W as a function of frequency
computed using different values of ?, respectively. The Re(?W ) and Im(?W ) show two
divergences as ? is varied (marked out by the blue arrows) and match the values of ?1 and
?2.
identical so that the corresponding ?W curves do not have an asymmetric shape. However,
we believe that one may utilize the difference in degenerate mode imaginary parts ?n
to reveal the 2-fold degeneracy. This can be accomplished by smoothly decreasing the
uniform loss ? value from high to low. For example, one may strategically set the initial
value of ? to be greater than max(?1,?2), and then decreasing ? until ? < min(?1,?2).
Note that this initial condition for ? is satisfied for the modes shown in Figs. 9.6 (b)
and (c). During this process, the value of ? will become equal to the ?n value of one of
the two degenerate modes. When this happens, we know from Eqns. 9.9 and 9.10 that
the time delay curve will diverge as the denominator ? 02 + 02. And since we usually
have ?1 ?= ?2, we can find that ?W curve ?diverges? twice during the process. One may
call this process a Lazarus time-delay behavior, where ? has made ?W ?die?verge twice to
demonstrate the power of degeneracy. We conducted a numerical exercise to illustrate this
process in Fig. 9.7. The ?W real and imaginary parts are computed using Eqns. 9.9 and
9.10, where we have set f1 = f2 = 24.2GHz (degenerate modes), ?1 = 0.007GHz,?2 =
141
Figure 9.8: a and b show the front and top views of the PTI-graph structure mounted
below the mixing chamber plate, respectively. c. The thermal connection (copper
connectors and heat straps) and cable connection on the PTI-graph structure. The inset
shows the heat-strap fixations on the mixing chamber plate.
0.0073GHz(?1 ?= ?2), and the uniform loss ? is swept from low to high through these
values of ?. It is shown that when the value of uniform loss ? equals either ?1 or ?2,
the time-delay curve diverges (the real-part flips sign, and the imaginary-part shows high
absolute values). This double time-delay divergence can not take place if only one mode
exists. We believe that a low-temperature experiment will allow us to gradually change
the degree of uniform loss by controlling the conductivity loss of the entire PTI structure
but keep the values of ?1 and ?2 approximately constant. In this way we can reveal the
existence of nearly degenerate modes from the evolution of complex time delay as ? is
systematically varied with temperature.
142
9.6 Low-temperature Experimental Results
Our first objective is to identify Kramers degenerate pairs of modes in the experiment.
We shall use the Lazarus method, created by cooling the graph to mK temperatures.
To vary the degree of loss ? in the PTI graph system, we have put the graph structure
inside a dilution refrigerator (dil-fridge). The dil-fridge cools the graph structure and
reaches a 30mK base temperature (measured by the mixing chamber plate thermometer).
During the cooldown from room temperature, the conductivity loss of the PTI graph will
decrease. Around the Tc of aluminum (1.2K), we should see a discontinuous change of
transmission measurement due to the superconducting transition of all the meral making
up the graph. The graph structure is mounted vertically to the mixing chamber plate (mxc
plate) using 4 copper blocks. As shown in Fig. 9.8 (a), the graph structure takes up nearly
the entire diameter of the mxc plate and can be considered as a large thermal load to the
dil-fridge. Two of the copper connection blocks can be seen in Fig. 9.8 (b). Fig. 9.8
(c) and its inset show the thermal braid connection and cable connection of the graph
structure. We used 5 braided heat straps to aid the thermal connection between the graph
and the mxc plate to ensure that the temperature of the graph is uniform and close to that
of the mxc plate.
In the first cooldown, we encountered several measurement issues that interrupted
the temperature sweep. Before discussing these practical issues, we first present the
analysis of measured ?W (f) curves during the cooldown. As discussed in the previous
section, we expect to find a double divergence of the ?W (f) as ? is swept through ?1
and ?2. We have found clear signs of the double divergence for the energy level near
143
Figure 9.9: a and b show the real and imaginary part of the measured PTI-graph ?W
versus temperature and frequency, respectively. c and d are the top view of a and b. To
guide the eye, white dashed lines are added to mark the temperatures at which Re(?W )
changes sign and Im(?W ) diverges.
144
24.45GHz. As shown in Fig. 9.9 (a) and (b), the ?W curves diverge twice when temperature
changes from 180K to 80K. To guide the eye, we use white dashed lines in 9.9 (c) and
(d) to mark out the locations of the ?W divergences (Re(?W ) changes sign and Im(?W )
shows high values). In the real part of ?W curves, we find a double sign flip, and in the
imaginary part of ?W curves, we find double high values. In Fig. 9.9 (c) and (d) we plot
the top view of the figures in Fig. 9.9 (a) and (b). It is clear that in Fig. 9.9 (c) and (d)
only one ridge region of ?W (f) exists during the measurement. This phenomenon means
that the frequencies of the two degenerate modes are almost identical. Fig. 9.9 (a) to (d)
shows exactly the case where the 2-fold degeneracy can not be identified in the different
real part of S-matrix poles/zeros (fi) but can be spotted from the difference in imaginary
parts (?i) during uniform loss (?) changes. We note that there exists an outlier event near
160K. It can be clearly spotted in 9.9 (a) and (c).
Here we discuss the issues discovered during the first cooldown. Firstly, it became
clear that we needed to install a thermometer onto the graph structure so that the PTI-
graph temperature could be better monitored. The rods and triangles in Fig. 9.5 are
attached to the aluminum plate using a mixture of non-magnetic brass and nylon screws.
We plan to change all of the screws from nylon to brass to achieve a better thermal
connection and less thermal contraction during the cooldown. Lastly, an odd calibration
issue was discovered. We have conducted a 4-port calibration at the end of the long
cables, which are connected at the top of the dil-fridge (Fig. 9.10 (a)). The resulting
graph measurements are shown in Fig. 9.10 (b) and (c). The S-parameter results (Fig.
9.10 (b)) are acceptable because we did see some differences between the measurements
taken at different temperatures. However, the computed ?W results (Fig. 9.10 (c)) show
145
Figure 9.10: a shows the 4 cable connections at the top of the fridge. These cables
connect to the VNA in room temperature, and go to the PTI-graphs in low temperature
environment. b and c show the measured S-parameter and the real part of ?W , respectively.
Different colors represent data measured at 130 different temperatures (from 80K to 4K).
146
rather noisy results. We hope to resolve these problems in the next cooldown.
9.7 Discussion
To sum up, we have proposed a PTI-graph system as a new platform for realizing
GSE statistics. An effective T-operator Tb for the BMW system is presented, and we
showed that the operator shows anti-unitary behavior (T 2b = ?1). A QVH-QSH composite
BMW structure is designed and fabricated. We first conducted a series of numerical
simulations for an ensemble of PTI-graph structures and conducted nearest neighbor level
spacing statistics. A relatively good agreement between the GSE theory prediction and
the simulated graph system data is found. We have identified the Kramers degeneracy
in experimental tests through the asymmetric shape of ?W (f) in room-temperature tests.
With the help of low-temperature measurements, we can systematically alter the system
uniform loss and find clear double divergence in the measured ?W . Thus, theoretical,
numerical, and experimental studies have presented clear evidence for Kramers degeneracy.
Our next objective will be conducting the nearest neighbor level spacing statistical test
with room-temperature experiments.
147
Chapter 10: Conclusions
In this chapter, I hope to summarize the thesis and discuss future directions for the
projects. Fig. 10.1 shows a Venn diagram which includes the projects in Chapters 2-9.
The general theme of the thesis is wave chaos. In Fig. 10.1, the wave chaos project
shows up as the biggest area circle. Under the theme of wave chaos, my projects focus
on studying the EM wave propagation in multiple connected chaotic cavities. From the
theoretical end, the Random Coupling Model (RCM) is successfully extended from single
enclosure systems to multiple coupled systems. I considered both few-mode coupling
as well as multi-mode coupling through different sizes of the aperture. Based on the
experimental structure designed by Bo Xiao, I conducted the experimental verification of
the multi-enclosure RCM theory. Building beyond this, we worked with our colleagues
at the University of Nottingham and applied the Power Balance Model to simplify the
multi-enclosure RCM model.
In addition to the frequency domain work, I have started to work on a time-domain
RCM (TD-RCM). Because the TD-RCM program is a finite difference method that operates
in the time-domain, we can include a variety of port-loads, including nonlinear elements.
We have also tested the agreement between TD-RCM and the experiment at late times
compared to the arrival of short orbits. For future work, the testing of early-time receiver
148
Figure 10.1: Summary of the projects. There are three themes in the thesis: wave chaos,
microwave topological photonics, and machine learning. The corresponding chapter
indices for the projects are also shown in this diagram. Wave chaos takes the center
stage. The topological photonics and the machine learning projects also intersects with
the wave chaos themes.
port signal and the experimental study of transmitting port reflected signal can be conducted.
We will also work on developing new types of statistical analysis that can also be simulated
with TD-RCM programs. An experimental study of the nonlinear port load cases is also
of interest. Another future direction can be packaging the TD-RCM codes into a graphical
software (through Matlab GUI).
As side projects, I have studied the photonic topological insulator (PTI) structures
and realized a composite PTI structure with both quantum Hall and quantum spin Hall
domains. I have also applied machine learning techniques to understand wave chaotic
system data details. These two fields of study are represented by two smaller circles in
Fig. 10.1. With the knowledge gained from these two projects centered around the theme
of wave chaos, I later started three projects: photonic crystal-graphs, PTI-graphs, and
wave-based reservoir computing. Because these projects are mixes of wave chaos and a
149
different field, they are pointed to the area where two circles intersected in Fig. 10.1. For
the wave-based reservoir computing project, I use the chaotic EM waves to operate as a
physical realization of reservoir computers. Future directions include the scaling down of
the resonator to 2D silicon photonics systems. One may also increase the mode-density
by changing the experimental structure from 2D to 3D. An increase of mode-density
will, in principle, enhance the degree of complexity of the wave-based reservoir. For the
PTI-graph project, I use a unique PTI waveguide to build a chaotic graph structure that
falls into the GSE (Gaussian Symplectic Ensemble) universality group. Future directions
include the low-temperature measurement of the PTI-graph system. By varying the bond
length of the PTI-graph, one can create an ensemble of graphs and study the system?s
eigenmode-spacing statistics, as well as scattering matrix statistics.
150
Appendix A: Multi-cavity RCM Modeling
Here we give a brief derivation of expressions for the input and trans-impedances of
a coupled multi-enclosure system through extension of the RCM. This material supplements
chapter 2.
A.1 Single Cavity Treatment
We write the radiation admittance matrix Yrad of a single cavity, as:
?
???Yrad,a 0 ?
?
Y = ?rad ? (A.1)
0 Yrad,b
Here Yrad is a block diagonal matrix, where the diagonal elements Yrad,a or Yrad,b
are either the radiation admittance of the port Yrad,port, or the radiation admittance matrix
of the aperture modes Yrad,aper, depending solely on whether a port or an aperture is
connected to that side of the cavity. The off-diagonal components are set to zero under
the assumption that there is no direct coupling between ports and apertures. We next
incorporate the RCM fluctuating quantities into the description of the fluctuating cavity
admittance. The cavity admittance matrix Ycav is:
151
Figure A.1: Schematic diagram of the ith cavity in a cascade with port a on the left and
port b on the right. Here U and I refer to the voltage and current at each port (a vector
quantity in general).
( ) ( )0.5 ( )0.5
Ycav = i ? Im Yrad + Re Yrad ? ? ?Re Yrad (A.2)
We represent the cavity admittance matrix Ycav (as shown in Fig. A.1) as
? ?
???Yaa YabY = ???cav (A.3)
Yba Ybb
We then have: ????Ia?
? ?
?? ?
? ? ? ? ?
= Y ?cav ?Ua??? ???Yaa Yab= ??? ???Ua??? (A.4)
Ib Ub Yba Ybb Ub
where the matrix ? is the normalized impedance which has the universal statistical properties
predicted by Random Matrix Theory (whose statistics depend only on the loss parameter
?). The dimension of the matrix ? is equal to the dimension of Yrad, which is Na + Nb,
where Na or Nb is the dimension of the radiation admittance matrices on the two sides of
152
this cavity.
A.2 Total Admittance of an N-Cavity Chain
After we construct the Ycav matrices for all cavities in the cascade, we are in position
to develop the cascade quantities for the enclosure chain. For the ith cavity in the cavity
cascade chain, with the information of how the ith cavity is connected to its left (previous)
and right (next) neighboring cavities, and knowledge of the load admittance presented by
the (i+1)(i+ 1)st cavity YL and everything beyond it, we will have an iterative approach to
calculating the total chain admittance matrix:
??????? ? ? ??????Ia??? ?
? ?Ua?
= Y (i) ? ? ?
cav ? ?
??? I U (A.5)? b b? ?I = Y (i+1)b L ? Ub
Eqn. A.5 (same as Eqn. 2.3) expresses the continuity of voltage and current in all
modes of the aperture. By solving Eqn. A.5 in matrix form, we have:
( )?1
Y (i) = Y (i) ? Y (i) Y (i) (i+1) (i)L aa ab bb + YL Yba (A.6)
( )?1
Ub = ? Y (i) (i+1) (i)bb + YL Yba Ua (A.7)
Eqn. A.6 and Eqn. A.7 are in the exact forms of Eqn. 2.5 and Eqn. 2.4 in the main
text. Eqn. A.6 gives the YL recursion relationship which calculates the load admittance
Y (i) thL of the i cavity given the Y (i+1)L of the (i + 1)st cavity and everything beyond
it. At the end of the cavity chain, the load impedance is the VNA measurement port
153
load impedance Z0 = 50?. We can now calculate the total load admittance of the entire
structure from the end load all the way back towards the input port at the first cavity.
A.3 Input and Trans-Impedances
We will investigate the trans-impedance Zt and the input impedance Zin of the
cascaded system from the voltage and YL recursion relationships (Eqn. A.6 and Eqn.
A.7). For the scalar Zin, we have
(1)
Ua 1
Zin = = (A.8)(1) (1)
Ia YL
where the (1)Ua and
(1)
Ia refers to the scalar voltage and current at the input side of the first
cavity.
From Eqn. A.6 and A.7, we have:
( )?1
(N) ? Y (N) + Y (N+1) (N) (N)U bb L Yba ? Ua
Z bt = =(1) ( (1)Ia Ia )?1
? Y (N) + Y (N+1) (N) ? (N?1)bb L Yba Ub
=
?[ ( (1)IaN ) ] (A.9)?1 (1)
? (i) (i+1) (i) ? Ua= ( Ybb )+ YL Yba (1)i=1 ? IaN?1
= ?(1) ?(i) ?(N) ? Zin
i=2
The quantities (N) (N) (1) (1)Ub , Ib , Ua , Ia are the scalar voltages and currents at the load
connected to the output side of the N th cavity (subscript b), and the input port at the first
cavity (subscript a). The multiplier ?(i) of the ith cavity is defined as
154
( )?1
?(i) = ? Y (i) + Y (i+1)bb L Y (i)ba . As discussed in the main text, the elements of ?(i)
are small because the matrix elements of Ybb are larger than Y (i)ba. The dimension of ? is
(Na, Nb) where Na,b refers to the number of propagating modes of the aperture at either
the a side or the b side of the cavity. Thus ?(i) can be either a matrix or a vector depending
on the specific location of the ith cavity in the entire cascade chain. For example, ?(i) is
in vector form at the first/last cavity (i = 1orN ) since the input/output of such cavity is a
single-mode port. ?(i) is a matrix at the intermediate cavities (i ? [2, N ? 1]) due to the
fact that these cavities are equipped with multi-mode apertures at both ends. With Eqn.
A.8 and A.9, we present the full theoretical formulas of Zin and Zt based on the RCM.
155
Appendix B: Multi-cavity RCM Experimental Details
B.1 Aperture and port radiation admittance studies
The radiation admittance Yrad and impedance Zrad refer to the cases where apertures
or ports radiate into free space [30, 37, 42]. Numerical simulations and experimental
methods are adopted in order to characterize the radiation information of the apertures
and ports employed in the experimental set-ups.
We use CST [222] to calculate the frequency dependent radiation admittance matrix
of both the scaled and full scale apertures. As shown in Fig. B.1, the aperture is carved
out of the surface of a Perfect Electrical Conductor (PEC) plate. The thickness of the
PEC plate equals to the thickness of the aperture of interest. The carved plate is attached
to a large vacuum box whose boundaries are assigned as radiating. The simulation is run
in the time domain solver mode. The waveguide port is assigned at the aperture surface
(blue circle in Fig. B.1). Given the operating frequency and the total number of port
(which is the aperture) modes, the simulation directly calculates the frequency dependent
complex radiation admittance matrix of the aperture. Similar techniques are employed in
the aperture cross section numerical study of Ref. [223].
The radiation impedance of the ports are obtained from experimental methods. For
the WR10 waveguide port used in the scaled set-up, the 1?1 radiation S-parameter (Srad)
156
Figure B.1: The model set-up for aperture admittance matrix calculations in the CST
simulations. A PEC plate is carved with the dimension of the aperture, and attached to a
vacuum box. The port is assigned at the 2D surface of the aperture. The wave enters the
aperture (blue arrow) and flows into the vacuum space before absorbed by the surrounding
radiation boundaries (orange arrows).
is measured by exposing the bare waveguide flange to a large radiating environment.
The Zrad is then transformed from the measured Srad. For the WR178 waveguide ports
used in the full scale measurements, time-gating techniques are adopted to calculate the
corresponding Zrad [224].
B.2 Convergence of the RCM formalism with inclusion of aperture modes
As introduced in the main text, the aperture radiation admittance Yrad is a (n, n)
matrix where n is the total number of considered aperture modes. Though the total
number of propagating modes can be calculated for given aperture dimensions and operating
frequency, the choice of the total number of non-propagating modes which are taken into
account is unclear. Here we demonstrate that the RCM cavity cascade formalism can find
convergence with an increasing number of considered non-propagating modes. As shown
157
Figure B.2: Statistics of the imaginary part of the trans-impedance for a scaled 3-cavity
cascade with circular aperture connections, illustrating the convergence of the theoretical
model with increasing mode number included in the aperture admittance matrix. The
circular aperture allows 100 propagating modes. The model PDFs curves are nearly
unchanged when the non-propagating modes are included.
in Fig. B.2, the convergence of the model is tested by adding circular aperture modes from
100 to 140 in the 3-cavity connection case. The statistics obtained from the theoretical
model with different aperture mode numbers remains unchanged beyond 100 modes, and
agrees well with the experimental results in that limit.
158
Figure B.3: (a) and (b): Statistics of the imaginary part of the trans-impedance of 2- and
3-cavity cascades when the loss of the system changes. (c): Comparison of the Re(Zt)
statistics experimental and model results for the scaled 2 cavity set-up, with circular
aperture connection. In (a), (b) and (c), the single cavity loss parameter ? is measured to
be ? = 5.7, 7.5 and 9.1, respectively.
B.3 Additional cavity cascade experimental studies
The 2- and 3-cavity cascade systems are studied experimentally with various single
cavity loss parameter values. As shown in Fig. B.3 (a) and (b), the statistics of 2- and
3-cavity model and experiment Im(Zt)s are shown with single cavity loss parameter ? =
5.7 and 7.5, respectively. The cascaded cavity structure studied in Fig. B.3 (a) and (b)
are full scale structures with circular shaped aperture connections. By placing 2 and 4
RF absorber cones in each individual cavity, single cavity loss parameters are measured
as ? = 5.7 and 7.5, respectively. We observed reasonably good agreement between the
model-generated Im(Zt) PDFs (solid lines) and the measured results (dashed lines). In
Fig. B.3 (c), the statistics of the Re(Zt) of the scaled-down 2-cavity cascade with circular
aperture connection are studied. We observed good agreement between the measured data
and model-generated results.
159
Appendix C: Multi-channel analysis for the cavity cascade system
Within the RCM description, the nth cavity in the cavity cascade chain can be
modeled by the following equation:
?? ?? ?? ? ? ???V i ?? ??Zii Zion n n ??? ???I in???= . (C.1)
V o oin Z Z
oo Io
n n n
The superscripts i and o refer to the input and output side of the nth cavity where we adopt
the convention that ?input? is on the side closest to the input TX port of the cavity chain
and ?output? is on the side with the RX port. The voltage vectors V and current vectors
I at the input and output sides of a cavity are connected through the cavity impedance
matrix as shown in Eq. ((C.1)). The Z-matrix of the nth cavity is written as a 2? 2 block
matrix. Assuming the nth cavity has m1 and m2 coupling channels at its input and output
Figure C.1: Schematic of a cavity cascade array with many-channel connecting apertures
between all cavities except that between cavity i and cavity i + 1. Note the definition of
input and output powers at each cavity aperture.
160
sides, the input vectors (V in, I
i
n) are m1 ? 1, and the output vectors (V on, Ion) are m2 ? 1.
Correspondingly, the dimension of the Z-matrix elements are m1?m1, m2?m2, m1?m2
and m ?m for Zii, Zoo2 1 , Zio and Zoi, respectively. For simplicity, we assume the RCMn n n n
description of the off-diagonal impedance matrices as
Zoi = (Ro )1/2? (Ri )1/2, (C.2)
n n n n
and the matrices Zoo,ii are diagonal. The Ri,o are the input and output radiation resistance
n n
matrices, and ? is the RCM normalized impedance matrix. Since the output side of the
n
nth cavity is connected to the input side of the (n+1)st cavity, the continuity relationships
of voltages and currents between the two cavities are:
V o = V i , Io = ?I in n+1 n n+1. (C.3)
In the high-loss limit, the output current at the (n+ 1)st cavity Ion+1 is much smaller than
the input current at the nth cavity I in. By solving Eqs. (C.1) and (C.3) and neglecting
terms Zio Ion+1 (small compared to the Z
ii terms), the relationship between the input
n+1
currents at the nth and (n+ 1)st cavity in the high-loss limit is written as
(Zii + Zoo)?1Zoi I i = I i .
n+1 n n n n+1
161
We then rewrite this equation as:
? = ? ? ? with
n+1 n?n+1 n n (C.4)
? = (Ri )1/2(Zii + Zoo)?1(Ro )1/2,
n?n+1 n+1 n+1 n n
a transition matrix which describes the coupling from the nth cavity to the (n+ 1)st with
Ri,o as in (C.2), and the ? = (Ri )1/2I in is a current-like vector that describes the inputn n n
signal of the cavity n. The power entering the nth cavity is given by,
1 [ ]i ? ii 1Pn = Re In Z I in = ??? (C.5)2 n 2 n n
We first assume the transition matrix ? to be diagonal. The input power at the (n + 1)st
cavity is written as:
1
P = ??n+1
2
1 (
?
n+1 n+1 )? ( )
= ?? ?( ? ? ? ? ) ? ?n n+1 n n n+1 n (C.6)2 n n1
= ?? ? 2n,i ?ji |?j| ?jk ?n,k2
i,j,k
The more general cases, where the con?straint o?f ? being diagonal is lifted, will be discussed
later in this section. With ?P ? = 1 ??? n ? ? ? being the average input power at the n
th
2 n n
cavity, and the fact that ?? ? ? ? = ??
?
ji j i ji?ji ?ii??jj? , we can further simplify Eq. (C.6) as:
? ? ? ? ??
?Pn+1? = ?Pn? ? |?|2 ? |?j|2 = ?P ? |?|2n Tj (C.7)
j j
162
where Tj = |?j|2 and the sum on j runs from all transmitting channels. Based on the
relationship in Eq. (C.7), we next analyze the fluctuation level of the power delivered
to the load (PL) by studying the ratio ? = ?P 2L? / ?PL?
2. To simplify the following
expression, X = ?? ?j n,j?n,j and Yj = ?n+1,j?n+?1,j are employed in the following calculations.
Equation (?C.6)?i?s now rewritten as Y = T ?
?
? j ? j ? i,i?? n,i?
?
n,i? ?ji?ji? , and further we have
? ? | |2 ? ? | |4 | |2 2
? ? ? ?
Yj = Tj ? i Xi . Utilizing ? = 2 ? and ?
? ? ? ? ?ji ji? ?ji?ji? = ?ji?ji ??ji? ?ji?? =
0 when i ?= i?, we have:
? ? ? ? [? ? ]
?YjYj?? = |?|2
2 ?X X ? T T + T 2i i?? j j? j (C.8)
jj? i,i?? jj? j
We then introduce a mean power transmission coefficient for the nth cavity T?n and an
effective total number of channels Mn, defined as:
? 2?? ? (? )2T?n = |?| T ?1 2j, Mn = Tj / Tj , (C.9)
j j j
respectively. Thus we have the power transmitted into the (n + 1)st cavity (Eq. (C.7))
written as:
?Pn+1? = T?n ?Pn? , (C.10)
and ? ? ( ) ? ?
P 2n+1 = T?
2
n 1 +M
?1
n P
2
n . (C.11)
163
The fluctuation level of the power after N cavities is:
?P 2? ?P 2? ?N ( ) ?N ( )
? = L 2 =
1
2 1 +M
?1 = 1 +M?1n n . (C.12)?PL? ?P1? n=1 n=1
Here we utilized the fact that the power injected into the first cavity (P1) is a fixed
scalar in the case of a steady input. Since (the last ca)vity is connected to the load with a
single-mode port in experiments, we have 1 +M?1N = 2 at the last cavity. The overall
fluctuating level of load power is
[ ]
? 2? NP ??1 ( )
? = L 2 = 2? 1 +M
?1
n . (C.13)?PL? n=1
For a system with a large number of connecting channels (Mn ? 1) between the neighbouring
cavities, such as the circular aperture connection which allows ? 100 propagating modes,
the factor (1 +M?1n ) ? 1. In this limit, the mean field methods are sufficient to describe
the power flow for cavities with strong coupling. On the contrary, the hybrid model is
not expected to work for system with few coupling channels inside a cavity cascade. The
study explains the experimental observations that hybrid models are successfully applied
to systems with multi-mode apertures connections, while fail to generate predictions for
systems with ?bottlenecks?.
We also note that the proposed hybrid model assumes each input current throughout
the chain is Gaussian distributed while the RCM predicts non-Gaussian statistics [33]. In
the presence of a ?bottleneck? connection, a deviation from Gaussian statistics of the
input current distribution will be present caused by the fluctuating nature of the narrow
164
aperture coupling.
We next expand the analysis to the more general cases where the transition matrix ?
is no longer a diagonal matrix. With Eqs. (C.6) and (C.8), the expression for Yj becomes:
? ?
Y = ? ? ? ?j jj??jj?? ?n,i??n,i?? ?j?i??j??i?? . (C.14)
j?j?? i?,i??
? ? ??
And correspon?ding?ly?the average of Yj is ?Yj? =
? 2
j? ?jj??jj? |?| i ?Xi?. With the
updated T? = |?|2 ?n jj? ?jj??jj? , we sum on j and have
? ? ?? ? ?
?Yj? = |?|2 ? ?jj??jj?? ?Xi? = T?n ?Xi? , (C.15)
j jj? i i
and the average of the product of Y ?s is
?
?YkYk?? = ? ? ? ?kj kj??k?j???k?j??? ?
j,j??,j??,j???? ? ? ? (C.16)
?? ? ? ? ?n,i n,i??n,i???n,i??? ?ji?j?i??j??i???j???i??? .
i,i?,i??,i???
With careful calculation, the double summation of Eq. (C.16) over k and k? gives:
? ? ? ? [? ]2
?YkYk??
2
/ |?|2 = ?XiX ?i?? ?kj?kj +
k,k? ?ii? ?jk[? ] (C.17)2
?XiXi?? ? ?kj?kj? .
ii? jj? k
With the updated T?n, we reach a similar expression of the power flow as in Eq. (C.11)
165
with the new M?1n written as
?[? ]2 [? ]2
M?1n = ?
?
kj?kj? / ?
?
jj??jj? . (C.18)
jj? k jj?
Thus the conclusions of the multi-channel analysis can be applied to the more general
cases with the refined ? , T?n and M?1n .
166
Appendix D: Applying the Hybrid Model to Generic Cavity Systems
Here we propose a scheme of applying the hybrid model to more general cases,
where there exists several limited-channel connections between the intermediate cavities.
Due to the large power flow fluctuations induced by these ?bottlenecks?, we apply RCM
to better characterize these connections. A general case is shown in Fig. C.1. We study
the treatment of limited channel connections between the ith and (i + 1)st cavities in a
N-cavity cascade chain. The treatment involves four steps:
1. Use the PWB N-cavity model to calculate the input, output power of the ith cavity
P iin| iPWB, Pout|PWB, and the input power at the last (N th) cavity PNin |PWB.
2. For the case of a bottleneck between cavity i and cavity i + 1 (Fig. C.1), use the
calculated P iin|PWB and RCM treatment to calculate an ensemble of the output
power of the ith cavity P iout|RCM .
3. Update the input power at the last cavity by
N
PN iin |RCM = Pout| ?
Pin
RCM |PWB.
P iout
Note that we therefore obtain an ensemble of the input power to the last cavity
PNin |RCM .
167
4. Use RCM to calculate an ensemble of the output power of the N th cavity PNout|RCM
with the ensemble of PNin |RCM obtained from (3).
For cases where there exist multiple ?bottlenecks? inside the chain (beyond cavity
i), one may repeat the procedure (2)-(4) by actively changing the cavity index N to the
next cavity connected through a bottleneck.
168
Appendix E: Simplified PWB-RCM Hybrid Method
Here we present a simplified treatment of the PWB-RCM hybrid model where
the full-wave simulated aperture radiation admittance matrix is no longer required. We
assume that all inter-cavity apertures are large on the scale of a wavelength, ie. the original
PWB-RCM treatment is applicable. We also will assume that all cavities are in the high-
loss (? > 1) limit. This Yaper-free treatment will further decrease the computational cost
of the hybrid model. We utilize Eqs. (C.4-C.6) from Appendix C to introduce the new
treatment.
Consider an N-cavity cascade chain, one may treat the load connected to the output
port of the N th cavity as the (N +1)st effective cavity. The induced power is equal to the
input power of the (N + 1)st cavity:
1
P ?load = PN+1 = ?(N+1 ?N+12 ) ( ) (E.1)1 ?
= ?N?N+1 ? ? ?N?N+1 ? ?
2 N n N n
where ? = (Ri )1/2(Zii + Zoo)?1(Ro )1/2N?N+1 N+1 N+1 N N is a scalar transition factor due to
the fact that a single-mode port acts as a single-mode ?aperture? connecting the N th and
(N + 1)st cavity (the load). As introduced in the main text, the radiation impedance
of the port, which is used to calculate ?N?N+1, is obtained from the scattering matrix
169
Figure E.1: The PDFs of load induced voltage |UL| of 2- and 3-cavity experiments
(solid), hybrid model (dashed), and the simplified hybrid model (dash-dotted). The single
cavity loss parameter is varied from 9.7, 7.5, 5.7 and 1.7 from (a-d), respectively.
170
measurements of the waveguide ports.
Equation (E.1) is further generalized as:
1 1
Pload = ?
?
N+1 ?N+1 = |? 2 ? ?N?N+1| ? ? ? ? . (E.2)
2 2 n N N n
The ensemble average of Pload is
1 ? ??P 2 ? ?load? = |?N?N+1| ? ?? ? ?2 n N N n? (E.3)1
= |? 2N?N+1| ??? ?? ? .
2 n n N N
In the above equation, the quantity ? is moved out of the ensemble averaging
n
because only ? is varying under ensemble average while ? and ?N?N+1 are not. AsN n
defined in Appendix C, we write 1/2? = Rii ? I in under the high-loss assumption. Heren N
Rii is set as the radiation impedance of the aperture which does not change value under
N
ensemble averaging. The input current vector I in is also considered to be not changing
since the cavity off-diagonal impedance element Zio (Eq. (C.1)) is small at high-loss limit,
n
and the input power of the N th cavity is fixed from the prior PWB calculation. Since the
input power of the N th cavity P = 1N ?? ? is given by th?e PWB?calculation, we may2 N N
f?urther ?simplify Eq. (E.3) as ?P ? = 1 |? |2 ? ? 2load N?N+1 ? ? ? ? = P2 n n N N N |?N?N+1| ?
?? ? . The vector element of ? , ?i, is generated using the same method as in SectionN N N
III. The vector product ?? ? is the summation of a large number of ??? s. Thus ?? ?
N N i i N N
and ??i ?i have the same expected value according to the law of large numbers. We finally
write ?P ? = P |? |2load N N?N+1 ???i ?i? and Eq. (E.2) is reduced to scalar components
Pload = P |? |2 ??N N?N+1 i ?i.
171
The comparison of the original hybrid model and the Yaper-free hybrid model is
presented in Fig. E.1. Good agreement between the simplified hybrid model predictions
and the experimental cases for the PDF of load induced voltages are found. The original
hybrid model slightly outperforms its simplified version for the 3-cavity cases (red curves
in Fig. E.1 (a-d)), and it is expected that the accuracy of predictions will decrease for
longer chains. Such an effect may be caused by the absence of an aperture admittance
matrix in the simplified hybrid model. The Yaper-free treatment would broaden the applicable
range of the PWB-RCM hybrid model to situations where the exact shape of the aperture
is unknown and the numerical simulation of the aperture is unavailable. However, this
shortcut shares the same high-loss (? > 1) assumption since it is built on the multi-
channel analysis (Appendix C). We also note that for situations discussed in Appendix
D, the admittance matrix of the narrow aperture (bottleneck) is required to compute the
transition matrix ? .
172
Appendix F: Wave-based Reservoir Computing Details
F.1 Simulation setup
We have conducted numerical simulations of the wave chaotic RC with the software
package CST Studio. The RC tasks are simulated with the time-domain solver, where the
bowtie-billiard is realized with the exact geometry and the same input waveforms used in
the experiment. The system parameters are matched to the experimental cases, including
the loss of the enclosure, the parameters of the microwave diodes, the amplitude and the
center frequency of the input signal. The high-dimensional combined RC is realized with
the RETs, where the boundary condition perturbation is realized using the cylindrical
metallic perturber with the same dimension as the experimental one. The simulated RCs
go through the same off-line training process using the simulated induced voltage signal
at the diode-loaded ports.
F.2 Observer tasks for dynamical systems
For the observer tasks, the input is a subset of a complex multi-variable dynamical
system and the output is the full set of variables of the system. We have applied the
observer task for two common chaotic attractors. The chaotic Ro?ssler system is governed
173
by
x? = ?y ? z,
y? = x+ ay, (F.1)
z? = b+ z(x? c),
where a = 0.5, b = 2.1, c = 3.5. The He?non-map is governed by
x 2n+1 = 1? a xn + yn,
(F.2)
yn+1 = b xn,
where a = 1.4, b = 0.3. The Ro?ssler x component is transformed into a continuously
waveform with a center frequency around 4GHz, and the RC is expected to produce
continuous predictions of both the y and z components without knowing the governing
equations. In the training period, the input is the Ro?ssler x component with ? 250
oscillation periods Tosc, and the testing period includes ? 50 periods. The input of the
discrete He?non-map is the x series where each bit is sampled for Tbit ? 60ps. Similarly,
the processing speed of the task can be further improved using equipment with higher
sampling rates. The training waveform of the He?non-map has a length of 4000 bits, and
the testing waveform is set to 1000 bits. The He?non-map observer task aims to infer the
y series solely with the knowledge of the x series. For both attractor observer tasks, the
input series are scaled to have a zero mean and unit standard deviation in order to balance
the positive and negative amplitudes. Depending on the nature of the chaotic attractors,
the optimal observer task performance is achieved under different RC parameters as
shown in Fig. 6.3. With a well-trained system, the wave chaotic RC is able to infer
174
the full set of chaotic attractors accurately as long as the observation of the system subset
is provided consistently.
F.3 Benchmark tests
For the NARMA10 task, the input is a random bit series u(n) drawn from a uniform
distribution over the interval [0, 0.5]. The target output is computed [fr?om the follow]ing
10th-order nonlinear relationship: y(n + 1) = 0.3y(n) + 0.05y(n) 9i=0 y(n? i) +
1.5u(n ? 9) ? u(n) + 0.1. Its complex behaviour and a 10-state memory requirement
make the NARMA10 task one of the most popular benchmark tests for both software and
hardware RC. In the experiment, the input is a time-domain waveform where each bit is
sampled for Tbit ? 60ps. We note that such choices of input waveform speeds are limited
by the sampling speed of both the AWG and the oscilloscope. The training waveform has
a length of 4000 random bits, and the testing waveform is set to 1000 bits.
For the function simulator task, the wave chaotic RC is expected to reproduce the
dynamics of the function s(t) = 3u(t)3, where u(t) is a continuous sine-wave input. In
4
the experiment, we employ a 4GHz sine wave input with a duration of ? 300 oscillation
periods. The length of the training and testing sets are set to a 80-20 ratio. We note that
the wave chaotic RC is also able to simulate other types of input functions, including
two-tone and three-tone signals.
For the nonlinear channel equalizer (NCE) task, the goal is to recover the random
4-level symbol sequence from a noisy series which simulates the received signal sent
through a nonlinear multi-path RF channel. The random input d(n) is chosen between
175
four levels {?3,?1, 1, 3}, and the channel output is defined as q(n) = 0.08d(n + 2) ?
0.12d(n+ 1) + d(n) + 0.18d(n? 1)? 0.1d(n? 2) + 0.091d(n? 3)? 0.05d(n? 4) +
0.04d(n ? 5) + 0.03d(n ? 6) + 0.01d(n ? 7). The input for NCE is the channel signal
q(n), and the task aims to retrieve the original four-level random d(n) series. In the
experiment, a similar input speed (Tbit ? 60ps) is adopted, and the training/testing set
includes 4000/1000 bits, respectively. The winner-takes-all method is applied where the
direct RC output is regularized to the nearest level.
F.4 RC operation speed analysis
Here we evaluate the effective operation speed of the wave chaotic reservoir computer.
The computation cost an RC task is estimated with a software RC problem [122]. Considering
a software RC with an N -dimensional reservoir layer, the execution of each input consumes
at least a ? N2 times digital operations (from a N ?N matrix multiplication). Thus the
total number of commands required for a M -point input signal is M ? N2. We note that
only the computation efficiency of the reservoir evolution is considered here, as the same
training procedure is adopted in both software and hardware RCs. We next consider the
same problem executed on a wave chaotic RC. For this M -point task, the RC operation
time is determined by the duration of the input series ?0 = M/fsam, where fsam is the
sampling rate of the AWG. A system with N0 output channels will reach a combined
N -dimensional reservoir with N/N0 times RET applications. Thus, the wave chaotic RC
completes the same problem with ?tot = ?0N/N0, and the effective operation speed is
M ? N2/?tot = N ? N0 ? fsam OPS. It is clear that the speed of the wave chaotic RC will
176
scale up linearly with larger N,N0, and fsam parameters.
177
Appendix G: Heterogeneous BMW Structure Details
G.1 Experimental Methods
The experimental set-up for the 4-port circulator is shown in an open-plate fashion
in Fig. G.1. The left (right) component are the top (bottom) plates with BMW rods
mounted. The QH island and QSH regions adopt the same color coding as the inset in
Fig. 5. The locations of the four ports are shown as ? signs in Fig. G.1. We used
Coilcraft 1812SMS-68NGLB loop inductor coils to excite the kink modes at the ports.
The magnetic coil antennas are adopted in all measurements shown from Fig. 2 to 5. We
can further use rotating magnetic dipole antennas to better excite the modes, which are
circularly polarized [29]. The S-parameter measurement is conducted with a Keysight
E8364C PNA Network Analyzer (2-port QH bulk region transmission measurement in
Fig. 2), and Keysight N5242A PNA-X Network Analyzer (2-port measurements in Fig. 3,
3-port measurements in Fig. 4 and 4-port measurement in Fig. 5). We mount the antenna
connected to port 2 (in Fig. 2) on the M-415.DG moving stage from Physik Instrumente
to measure the bulk transmission at various locations along the edge of the structure.
The QH region consists of Aluminum rods, YG-14 Yttrium-Garnet ferrites from HTD-
RF Technology and further magnetized by D150E Nd-Fe-B permanent magnets from
Amazing Magnets. The transmission results are first calibrated by an S-parameter measurement
178
of an empty BMW structure where the rods are removed (i.e., just parallel plates with
antennas) and then normalized by the maximum value. Such a calibration method (see SM
of [14]) is meant to remove the impedance mismatch between antenna and the experimental
structure.
Figure G.1: An open-plate view of the experimental set-up of the 4-port BMW circulator
structure. Each region is color coded to match the insets in Fig. 5.
G.2 Ferrite Material
The YG-14 Yttrium-Garnet ferrites have saturation magnetization Ms = 1400Oe
and resonance line-width ?H = 25Oe. The permeability tensor [47] adopted in the
simulation models is:
179
?? ??? ? (?) ?? 1 2
(?) 0???
?(?) = ?????2(?) ?1(?) 0???? (G.1)
0 0 1
where ? (?) = ? + ?m(?l+i??) and ? (?) = i??m1 ? 2? 2 2 . The Lamor frequency(?l+i??) ? (?l+i??)2??2
?l = ??0H0 where ? is the gyromagnetic ratio and H0 is the applied H-field. The
gyrotropic frequency is ?m = ?4?Ms where Ms is the saturation magnetization. The
damping factor ? is defined here as ? = ?0??H where ?H is the ferrite resonance line-
2?m
width as a function of magnetic field at a certain reference frequency ?m. The BMW
bandgap frequency range is near the ferromagnetic anti-resonance regime for the biasing
conditions used in this experiment. The Polder tensor values used in the COMSOL
simulations are assumed to be independent of frequency and evaluated at the Dirac point
(6.08 GHz).
G.3 Kink mode Back-scattering Studies
We first conduct band structure simulations in COMSOL for a QH-QSH composite
PTI structure. The simulation model shown in Fig. G.2 is a supercell structure: a strip
of BMW rods with 15 periods in both the QH and QSH regions, and an interface in the
middle. We will invert the applied H-field direction in the QH region and study the kink
mode properties at the QH-QSH interface. The band structures show that the propagating
direction of kink modes is the same in both the K and ?K direction. This kink mode
propagating direction is inverted when the H-field direction is inverted (compare Fig.
180
Figure G.2: The photonic band structure of a 15+15 period QH-QSH Supercell structure
from COMSOL simulations. The insets are the supercell structures near the QH-QSH
interface. The dashed lines represent the location of the QH-QSH interface. Under
applied magnetic field in opposite directions (?Hz/ +Hz in (a)/(b)), the supported kink
mode will have opposite synthetic spin (up/down-spin shown as ??/?) and group velocity.
Dirac points are at kxa0 = ?2?/3. Boxed tags: the Chern numbers of the bulk modes
below the bandgap for QSH and QH domains.
G.2(a) and (b)).
Transmission measurements at the QH-QSH interface are shown in Fig. G.4. We
excite the kink mode at one end of the interface and measure the kink mode transmission
at various pick-up locations along the interface. As shown in Fig. G.4(a), the kink mode
transmission will not be damped significantly when the measurement location moves
further from the source along the propagating direction. On the other end, the signal
amplitude on the non-propagating direction will remain 25 dB lower than the propagating
direction at various pick-up locations as shown in Fig. G.4(b).
181
Figure G.3: An open-plate view of the QH-QSH interface 2-port isolator experimental
set-up. Note that the biasing magnets on the QH region are not shown.
Figure G.4: Plots of kink mode transmission measurement in (a) propagating and (b)
non-propagating directions along a QH-QSH interface. In both plots, the S-parameters
represent the transmission between two ports: S2 is measured at a port further from the
source (port A) than S1. The plots show that the amplitude of the kink mode along the
propagating direction does not attenuate significantly, while the amplitude of kink mode
in the backward direction is always low.
182
G.4 Coupled waveguide issue
In transmission measurements of the 4-port circulator structure, we observed un-
expected increase of transmission between two parallel QSH-QSH waveguides near 6.1
GHz as shown in Fig. 5 (a) and (d). This phenomenon is caused by the finite-size
problem. As shown in the Fig. G.1, the middle QSH area (highlighted by blue boundaries)
is only 7 lattice periods in size. The close proximity between two parallel QSH-QSH
waveguides will induce coupling between them. In order to decrease the coupling effect
between two parallel QSH-QSH waveguides, we employed DD-13490 absorber from
ARC Technologies in the middle QSH (blue) area but away from the QH-QSH interfaces
(red boundaries in Fig. G.1). The absorber will decrease the transmission between two
coupled QSH-QSH waveguides, while not disturbing the QH-QSH waveguides around
the center QH island.
It is well known that two waveguides can couple to each other via their fringing
fields and lead to power transfer from one waveguide to another. This can be verified
from the PBS (photonic band structure) of a supercell consisting of three QSH region
with alternating ?SOC , where the central region contains 7 unit cells. Fig. G.5 shows
the PBS near the K? valley. There are 4 modes in the figure, originated from 2 TPSWs
(one forward, one backward propagating) at each waveguide (QSH-QSH interface). It
is evident from the figure that coupling between co-propagating modes and counter-
propagating modes both takes place. Most importantly, the coupling between counter-
propagating modes creates a secondary bandgap around 6.06-6.09 GHz. If one launches
a TPSW at one interface, it will exponentially decay as it propagates and gets reflected
183
Figure G.5: The 3-QSH supercell PBS near the Dirac point. The positive/negative
group velocity parts (blue/red) are mixture of the up/down-spin mode localized at the top
interface and down/up-spin mode at the bottom interface. Near the secondary bandgap at
about 6.06-6.09 GHz, there is a four-mode mixing.
184
Figure G.6: Supercell simulations results of the kink mode at the Dirac point. Color:
|E| of the eigenmodes near 6.1 GHz. The 5 and 7 periods of unit cell separation cases
cause E-field spreading in the entire middle QSH region, while the 9-period separation
case makes the kink modes better localized at the interfaces.
through the other interface. This is exactly what is happening in Fig. 5 near 6.08 GHz,
causing high transmission S31 and S42, in (a) and (d), respectively.
To illustrate the coupled waveguide effect, we conducted supercell parallel waveguide
simulations in COMSOL as shown in Fig. G.6. The supercell consists of three QSH bulk
regions with alternating ?SOC signs. We examine the E-field profile of the eigenmodes
while increasing the size of the middle QSH region. To overcome the finite-size problem,
we find that at least a 9-period lattice in the middle region is required to establish clear
isolation between two parallel waveguides.
G.5 QH-QSH structure with no external magnetic field
In contrast to the QH-QSH interface transmission experiments under external H-
field shown in Fig. 3, we conducted a zero H-field test using the same experimental
185
Figure G.7: Zero magnetic field transmission measurement between port A and B. Two
ports are in the same locations as in the measurement shown in Fig. 3, located on the QH-
QSH interface. We find that the transmission profile between two ports are very similar,
indicating the time-reversal symmetry in not broken in the un-magnetized ferrite case.
structure shown in Fig. G.7. The QH bulk region will no longer have a band gap in the
absence of an applied magnetic field (see QH bulk transmission in Fig. 2). Though the
QSH region is insensitive to the magnetic field, an interface between QSH and a unbiased
QH region is not able to host uni-directional TPSWs. Compared to the 35dB contrast
between left/right going waves in Fig. 3, the transmission between port A and B are
nearly the same in both directions.
186
G.6 Topological Y-junction
Figure G.8 shows an open-plate view of the topological Y-junction.
Figure G.8: The open-plate view of the topological Y-junction. The magnets for QH
region and the upper QSH rods and are not shown in this view. The three topological
regions are color-coded as in Fig. 4(a). The dashed lines mark the interface between
different topological phases. The stars are the locations of the ports.
187
Appendix H: Extended Symplectic Group Analysis of BMW Hamiltonian
Here we present a quick proof to show that the BMW Hamiltonian falls into the
symplectic group. The effective Hamiltonian of the BMW lattice can be written as
H(?k) = H0 + HSOC/P , where ?k = (?kx, ?ky) is the wavenumber measured from
the Dirac points K and K ?. H0 and HSOC/P are the unperturbed Hamiltonian and the
perturbation term, respectively. The subscripts represents the type of perturbations that
are introduced inside the BMW unit cell and which give rise to different topological
phases, namely the spin-orbit coupling (for QSH-PTI) and parity breaking (for QVH-
PTI). Here we note that both QSH and QVH-BMW systems are suitable as platforms to
realize GSE statistics. For simplicity, we focus on the QSH-BMW analysis in this paper,
while a similar analysis also applies to the QVH-BMW system. The effective Hamiltonian
of the QSH-BMW lattice at the K and K ? valleys are written as [206]:
H??K (?kx, ?ky) = vD s0(?kx?x + ?ky?y) + ?D?emsz ?z (H.1)
H??K?(?kx, ?ky) = vD s0(??kx?x + ?ky?y)? ?D?emsz ?z (H.2)
where s and ? are the Pauli matrices acting on the spin and the Dirac degeneracy, respectively.
The arrow ? (?) refers to the up (down) spin index. The quantities ?D and vD are the
188
frequency and the group velocity at the Dirac points. Note that the Kronecker product
(s ? ? s??) shorthand is used. The BMW Hamiltonian is Hermitian with real eigenvalues.
It is shown by Mehta and Dyson that a quaternion-real matrix H can be diagonalized by a
symplectic matrix [9, 10]. Our goal is to show that the BMW Hamiltonian is quaternion-
real (see the definition in Ref. [9]). For simplicity, we loo?k at the B?MW Hamiltonian at?q11 0 ?
the K-valley (Eqn. H.1). In matrix form, we have H??K = ?? ??, where q11, q22 are
0 q22
quaternion numbers written as q11 = 0 + ? ? q? and q22 = 0 + ? ? q?. The quantities
q? = (vD?kx, vD?ky, ?D?em?) and q? ?= (vD?kx, vD?ky,??D?em) are real. Thus it is?q11 0 ?
clear that the matrix H = ?? ??K is quaternion-real and fits into the symplectic
0 q22
group.
As discussed above, we have shown that the BMW Hamiltonian falls into the
GSE universality group, and further written out the anti-unitary T-operator for the BMW
system that leaves the Hamiltonian invariant.
189
Bibliography
[1] Y. Y. Atas, E. Bogomolny, O. Giraud, and G. Roux, Physical Review Letters 110,
084101 (2013).
[2] A. Gokirmak, D.-H. Wu, J. S. A. Bridgewater, and S. M. Anlage, Review of
Scientific Instruments 69, 3410 (1998).
[3] B. Xiao, T. M. Antonsen, E. Ott, and S. M. Anlage, Physical Review E 93, 052205
(2016).
[4] E. Ott, Chaos in Dynamical Systems (Cambridge University Press, Cambridge,
2002).
[5] E. N. Lorenz, (1972).
[6] O. E. Rossler, Physics Letters A 71, 155 (1979).
[7] E. N. Lorenz, Journal of the Atmospheric Sciences 20, 130 (1963).
[8] E. P. Wigner, The Annals of Mathematics 62, 548 (1955).
[9] F. J. Dyson, Journal of Mathematical Physics 3, 140 (1962).
[10] M. L. Mehta, Academic Press (2004) p. 706, arXiv:0509286 [hep-ph] .
[11] O. Bohigas, M. J. Giannoni, and C. Schmit, Physical Review Letters 52, 1 (1984).
[12] R. U. Haq, A. Pandey, and O. Bohigas, Physical Review Letters 48, 1086 (1982).
[13] Y. Alhassid, Reviews of Modern Physics 72, 895 (2000).
[14] D. H. Wu, J. S. A. Bridgewater, A. Gokirmak, and S. M. Anlage, Physical Review
Letters 81, 2890 (1998).
[15] B. Dietz and A. Richter, Chaos: An Interdisciplinary Journal of Nonlinear Science
25, 097601 (2015).
[16] S. Ma, B. Xiao, Z. Drikas, B. Addissie, R. Hong, T. M. Antonsen, E. Ott, and
S. M. Anlage, Physical Review E 101, 022201 (2020).
190
[17] E. Doron, U. Smilansky, and A. Frenkel, Physical Review Letters 65, 3072 (1990).
[18] P. So, S. M. Anlage, E. Ott, and R. N. Oerter, Physical Review Letters 74, 2662
(1995).
[19] U. Kuhl, M. Mart??nez-Mares, R. A. Me?ndez-Sa?nchez, and H.-J. Sto?ckmann,
Physical Review Letters 94, 144101 (2005).
[20] M. Hoijer and L. Kroon, IEEE Transactions on Electromagnetic Compatibility 55,
1328 (2013).
[21] A. Gifuni, I. D. Flintoft, S. J. Bale, G. C. R. Melia, and A. C. Marvin, IEEE
Transactions on Electromagnetic Compatibility 58, 678 (2016).
[22] B. Dietz, T. Klaus, M. Miski-Oglu, A. Richter, M. Wunderle, and C. Bouazza,
Physical Review Letters 116, 023901 (2016).
[23] M. Dupre?, P. Del Hougne, M. Fink, F. Lemoult, and G. Lerosey, Physical Review
Letters 115, 017701 (2015).
[24] N. Kaina, M. Dupre?, G. Lerosey, and M. Fink, Scientific Reports 4, 6693 (2015).
[25] P. del Hougne, M. F. Imani, T. Sleasman, J. N. Gollub, M. Fink, G. Lerosey, and
D. R. Smith, Scientific Reports 8, 6536 (2018).
[26] J. P. Parmantier, IEEE Transactions on Electromagnetic Compatibility 46, 359
(2004).
[27] X. Zheng, T. M. Antonsen, and E. Ott, Electromagnetics 26, 37 (2006).
[28] X. Zheng, T. M. Antonsen, and E. Ott, Electromagnetics 26, 3 (2006).
[29] J. A. Hart, T. M. Antonsen, and E. Ott, Physical Review E 80, 041109 (2009).
[30] J.-H. Yeh, J. A. Hart, E. Bradshaw, T. M. Antonsen, E. Ott, and S. M. Anlage,
Physical Review E 82, 041114 (2010).
[31] G. Gradoni, J.-H. Yeh, B. Xiao, T. M. Antonsen, S. M. Anlage, and E. Ott, Wave
Motion 51, 606 (2014).
[32] X. Zheng, S. Hemmady, T. M. Antonsen, S. M. Anlage, and E. Ott, Physical
Review E 73, 046208 (2006).
[33] G. Gradoni, T. M. Antonsen, and E. Ott, Physical Review E 86, 046204 (2012).
[34] X. Li, C. Meng, Y. Liu, E. Schamiloglu, and S. D. Hemmady, IEEE Transactions
on Electromagnetic Compatibility 57, 448 (2015).
[35] J. Q. Fan, Y. Pan, J. H. Hao, and H. Y. Zhang, Progress in Electromagnetics
Research Letters 65, 81 (2017).
191
[36] L. Luyao, M. Hongge, and W. Ming, in 2017 IEEE 2nd Advanced Information
Technology, Electronic and Automation Control Conference (IAEAC) (IEEE, 2017)
pp. 1526?1529.
[37] S. Hemmady, X. Zheng, E. Ott, T. M. Antonsen, and S. M. Anlage, Physical
Review Letters 94, 014102 (2005).
[38] S. Hemmady, X. Zheng, J. Hart, T. M. Antonsen, E. Ott, and S. M. Anlage,
Physical Review E 74, 036213 (2006).
[39] M. V. Berry, European Journal of Physics 2, 91 (1981).
[40] Y. V. Fyodorov and D. V. Savin, Journal of Experimental and Theoretical Physics
Letters 80, 725 (2004).
[41] A. Rehemanjiang, M. Allgaier, C. H. Joyner, S. Mu?ller, M. Sieber, U. Kuhl, and
H.-J. Sto?ckmann, Physical Review Letters 117, 064101 (2016).
[42] G. Gradoni, T. M. Antonsen, S. M. Anlage, and E. Ott, IEEE Transactions on
Electromagnetic Compatibility 57, 1049 (2015).
[43] M. Notomi, E. Kuramochi, and T. Tanabe, Nature Photonics 2, 741 (2008).
[44] A. Wallraff, D. I. Schuster, A. Blais, L. Frunzio, R.-S. Huang, J. Majer, S. Kumar,
S. M. Girvin, and R. J. Schoelkopf, Nature 431, 162 (2004).
[45] H. Alt, C. I. Barbosa, H.-D. Gra?f, T. Guhr, H. L. Harney, R. Hofferbert, H. Rehfeld,
and A. Richter, Physical Review Letters 81, 4847 (1998).
[46] A. O?. Kaya, L. J. Greenstein, and W. Trappe, IEEE Transactions on Wireless
Communications 8, 4165 (2009).
[47] W. Porod, Z.-A. Shao, and C. S. Lent, Applied Physics Letters 61, 1350 (1992).
[48] F. R. Waugh, M. J. Berry, D. J. Mar, R. M. Westervelt, K. L. Campman, and A. C.
Gossard, Physical Review Letters 75, 705 (1995).
[49] R. Brunner, R. Akis, D. K. Ferry, F. Kuchar, and R. Meisels, Physical Review
Letters 101, 024102 (2008).
[50] J. T. Chalker and P. D. Coddington, Journal of Physics C: Solid State Physics 21,
2665 (1988).
[51] J. D. Sau and S. D. Sarma, Nature Communications 3, 964 (2012).
[52] C. W. J. Beenakker, Reviews of Modern Physics 87, 1037 (2015).
[53] P. Zhang and F. Nori, New Journal of Physics 18, 043033 (2016).
[54] B. Dietz, T. Guhr, H. L. Harney, and A. Richter, Physical Review Letters 96,
254101 (2006).
192
[55] G. B. Tait, R. E. Richardson, M. B. Slocum, and M. O. Hatfield, IEEE Transactions
on Electromagnetic Compatibility 53, 846 (2011).
[56] G. B. Tait, R. E. Richardson, M. B. Slocum, M. O. Hatfield, and M. J. Rodriguez,
IEEE Transactions on Electromagnetic Compatibility 53, 229 (2011).
[57] I. Junqua, J.-P. Parmantier, and F. Issac, Electromagnetics 25, 603 (2005).
[58] S. Hemmady, T. M. Antonsen, E. Ott, and S. M. Anlage, IEEE Transactions on
Electromagnetic Compatibility 54, 758 (2012).
[59] L. Gagliardi, D. Micheli, G. Gradoni, F. Moglie, and V. M. Primiani, IEEE
Antennas and Wireless Propagation Letters 14, 1463 (2015).
[60] B. Xiao, T. M. Antonsen, E. Ott, Z. B. Drikas, J. G. Gil, and S. M. Anlage, Physical
Review E 97, 062220 (2018).
[61] D. A. Hill, Electromagnetic Fields in Cavities (IEEE, 2009) p. 280.
[62] J. Gil Gil, Z. B. Drikas, T. D. Andreadis, and S. M. Anlage, IEEE Transactions on
Electromagnetic Compatibility 58, 1535 (2016).
[63] M. Zhou, E. Ott, T. M. Antonsen, and S. M. Anlage, Chaos 27, 103114 (2017).
[64] J.-H. Yeh, Z. Drikas, J. Gil Gil, S. Hong, B. Taddese, E. Ott, T. Antonsen,
T. Andreadis, and S. Anlage, Acta Physica Polonica A 124, 1045 (2013).
[65] J.-H. Yeh and S. M. Anlage, Review of Scientific Instruments 84, 034706 (2013).
[66] S. Ma, S. Phang, Z. Drikas, B. Addissie, R. Hong, V. Blakaj, G. Gradoni, G. Tanner,
T. M. Antonsen, E. Ott, and S. M. Anlage, Physical Review Applied 14, 14022
(2020).
[67] D. Hill, M. Ma, A. Ondrejka, B. Riddle, M. Crawford, and R. Johnk, IEEE
Transactions on Electromagnetic Compatibility 36, 169 (1994).
[68] D. Hill, IEEE Transactions on Electromagnetic Compatibility 40, 209 (1998).
[69] D. Coppersmith and S. Winograd, Journal of Symbolic Computation 9, 251 (1990).
[70] S. Ma, B. Xiao, R. Hong, B. Addissie, Z. Drikas, T. Antonsen, E. Ott, and
S. Anlage, Acta Physica Polonica A 136, 757 (2019).
[71] P. Mehta, C.-H. Wang, A. G. R. Day, C. Richardson, M. Bukov, C. K. Fisher,
and D. J. Schwab, A high-bias, low-variance introduction to Machine Learning for
physicists, Tech. Rep. (2018) arXiv:1803.08823v1 .
[72] S. Das Sarma, D.-L. Deng, and L.-M. Duan, Physics Today 72, 48 (2019).
[73] J. Carrasquilla and R. G. Melko, Nature Physics 13, 431 (2017).
193
[74] J. Venderley, V. Khemani, and E.-A. Kim, Physical Review Letters 120, 257204
(2018).
[75] A. Seif, K. A. Landsman, N. M. Linke, C. Figgatt, C. Monroe, and M. Hafezi,
Journal of Physics B: Atomic, Molecular and Optical Physics 51, 174006 (2018).
[76] C. J. van Diepen, P. T. Eendebak, B. T. Buijtendorp, U. Mukhopadhyay, T. Fujita,
C. Reichl, W. Wegscheider, and L. M. K. Vandersypen, Applied Physics Letters
113, 033101 (2018).
[77] S. S. Kalantre, J. P. Zwolak, S. Ragole, X. Wu, N. M. Zimmerman, M. D. Stewart,
and J. M. Taylor, npj Quantum Information 5, 6 (2019).
[78] Z. Lu, J. Pathak, B. Hunt, M. Girvan, R. Brockett, and E. Ott, Chaos: An
Interdisciplinary Journal of Nonlinear Science 27, 041102 (2017).
[79] J. Pathak, B. Hunt, M. Girvan, Z. Lu, and E. Ott, Physical Review Letters 120,
024102 (2018).
[80] G. Neofotistos, M. Mattheakis, G. D. Barmparis, J. Hizanidis, G. P. Tsironis, and
E. Kaxiras, Frontiers in Physics 7, 1 (2019).
[81] M. Frazier, B. Taddese, T. Antonsen, and S. M. Anlage, Physical Review Letters
110, 063902 (2013).
[82] M. Frazier, B. Taddese, B. Xiao, T. Antonsen, E. Ott, and S. M. Anlage, Physical
Review E 88, 062910 (2013).
[83] Z. B. Drikas, J. Gil Gil, S. K. Hong, T. D. Andreadis, J.-H. Yeh, B. T. Taddese,
and S. M. Anlage, IEEE Transactions on Electromagnetic Compatibility 56, 1480
(2014).
[84] Y. Lecun, Y. Bengio, and G. Hinton, ?Deep learning,? (2015),
arXiv:arXiv:1312.6184v5 .
[85] S. Ma, T. Antonsen, S. Anlage, and E. Ott, (2021), 10.21203/rs.3.rs-783820/v1.
[86] Y. Zhang, A. Mesaros, K. Fujita, S. D. Edkins, M. H. Hamidian, K. Ch?ng,
H. Eisaki, S. Uchida, J. C. S. Davis, E. Khatami, and E.-A. Kim, Nature 570,
484 (2019).
[87] Y. K. Chembo, Chaos: An Interdisciplinary Journal of Nonlinear Science 30,
013111 (2020).
[88] B. J. Shastri, A. N. Tait, T. Ferreira de Lima, W. H. P. Pernice, H. Bhaskaran, C. D.
Wright, and P. R. Prucnal, Nature Photonics 15, 102 (2021).
[89] H. Esmaeilzadeh, E. Blem, R. St. Amant, K. Sankaralingam, and D. Burger, IEEE
Micro 32, 122 (2012).
194
[90] K. Roy, A. Jaiswal, and P. Panda, Nature 575, 607 (2019).
[91] P. A. Merolla, J. V. Arthur, R. Alvarez-Icaza, A. S. Cassidy, J. Sawada, F. Akopyan,
B. L. Jackson, N. Imam, C. Guo, Y. Nakamura, B. Brezzo, I. Vo, S. K. Esser,
R. Appuswamy, B. Taba, A. Amir, M. D. Flickner, W. P. Risk, R. Manohar, and
D. S. Modha, Science 345, 668 (2014).
[92] J. Feldmann, N. Youngblood, C. D. Wright, H. Bhaskaran, and W. H. P. Pernice,
Nature 569, 208 (2019).
[93] Y. Shen, N. C. Harris, S. Skirlo, M. Prabhu, T. Baehr-Jones, M. Hochberg, X. Sun,
S. Zhao, H. Larochelle, D. Englund, and M. Soljac?ic?, Nature Photonics 11, 441
(2017).
[94] P. Yao, H. Wu, B. Gao, J. Tang, Q. Zhang, W. Zhang, J. J. Yang, and H. Qian,
Nature 577, 641 (2020).
[95] J. Torrejon, M. Riou, F. A. Araujo, S. Tsunegi, G. Khalsa, D. Querlioz,
P. Bortolotti, V. Cros, K. Yakushiji, A. Fukushima, H. Kubota, S. Yuasa, M. D.
Stiles, and J. Grollier, Nature 547, 428 (2017).
[96] D. Canaday, A. Griffith, and D. J. Gauthier, Chaos: An Interdisciplinary Journal
of Nonlinear Science 28, 123119 (2018).
[97] K. Nakajima, K. Fujii, M. Negoro, K. Mitarai, and M. Kitagawa, Physical Review
Applied 11, 034021 (2019).
[98] F. Laporte, Novel architectures for brain-inspired photonic computers, Ph.D.
thesis, Universiteit Gent (2020).
[99] Y. Zhong, J. Tang, X. Li, B. Gao, H. Qian, and H. Wu, Nature Communications
12, 408 (2021).
[100] A. Sirunyan, A. Tumasyan, W. Adam, F. Ambrogi, and A. et al., Journal of
Instrumentation 15, P10017 (2020).
[101] G. Tanaka, T. Yamane, J. B. He?roux, R. Nakane, N. Kanazawa, S. Takeda,
H. Numata, D. Nakano, and A. Hirose, Neural Networks 115, 100 (2019).
[102] N. Mohammadi Estakhri, B. Edwards, and N. Engheta, Science 363, 1333 (2019).
[103] T. W. Hughes, I. A. D. Williamson, M. Minkov, and S. Fan, Science Advances 5,
eaay6946 (2019).
[104] G. Wetzstein, A. Ozcan, S. Gigan, S. Fan, D. Englund, M. Soljac?ic?, C. Denz,
D. A. B. Miller, and D. Psaltis, Nature 588, 39 (2020).
[105] J. Feldmann, N. Youngblood, M. Karpov, H. Gehring, X. Li, M. Stappers, M. Le
Gallo, X. Fu, A. Lukashchuk, A. S. Raja, J. Liu, C. D. Wright, A. Sebastian, T. J.
Kippenberg, W. H. P. Pernice, and H. Bhaskaran, Nature 589, 52 (2021).
195
[106] L. Appeltant, M. Soriano, G. Van der Sande, J. Danckaert, S. Massar, J. Dambre,
B. Schrauwen, C. Mirasso, and I. Fischer, Nature Communications 2, 468 (2011).
[107] Q. Vinckier, F. Duport, A. Smerieri, K. Vandoorne, P. Bienstman, M. Haelterman,
and S. Massar, Optica 2, 438 (2015).
[108] G. Van Der Sande, D. Brunner, and M. C. Soriano, Nanophotonics 6, 561 (2017).
[109] L. Grigoryeva and J.-P. Ortega, Neural Networks 108, 495 (2018).
[110] P. Vlachas, J. Pathak, B. Hunt, T. Sapsis, M. Girvan, E. Ott, and P. Koumoutsakos,
Neural Networks 126, 191 (2020).
[111] X. Lin, Y. Rivenson, N. T. Yardimci, M. Veli, Y. Luo, M. Jarrahi, and A. Ozcan,
Science 361, 1004 (2018).
[112] K. Vandoorne, P. Mechet, T. Van Vaerenbergh, M. Fiers, G. Morthier,
D. Verstraeten, B. Schrauwen, J. Dambre, and P. Bienstman, Nature
Communications 5, 3541 (2014).
[113] P. O?Connor, J. Gehlen, and E. J. Heller, Physical Review Letters 58, 1296 (1987).
[114] C. Liu, R. E. C. van der Wel, N. Rotenberg, L. Kuipers, T. F. Krauss, A. Di Falco,
and A. Fratalocchi, Nature Physics 11, 358 (2015).
[115] K. Pichler, M. Ku?hmayer, J. Bo?hm, A. Brandsto?tter, P. Ambichl, U. Kuhl, and
S. Rotter, Nature 567, 351 (2019).
[116] B. W. Frazier, T. M. Antonsen, S. M. Anlage, and E. Ott, Physical Review
Research 2, 043422 (2020).
[117] L. Chen, T. Kottos, and S. M. Anlage, Nature Communications 11, 5826 (2020).
[118] K. Kim, S. Bittner, Y. Zeng, S. Guazzotti, O. Hess, Q. J. Wang, and H. Cao,
Science 371, 948 (2021).
[119] C. Fernando and S. Sojakka, in Lecture Notes in Artificial Intelligence (Subseries
of Lecture Notes in Computer Science), Vol. 2801 (Springer Verlag, 2003) pp. 588?
597.
[120] F. Laporte, A. Katumba, J. Dambre, and P. Bienstman, Optics Express 26, 7955
(2018).
[121] G. Marcucci, D. Pierangeli, and C. Conti, Physical Review Letters 125, 093901
(2020).
[122] M. Rafayelyan, J. Dong, Y. Tan, F. Krzakala, and S. Gigan, Physical Review X 10,
041037 (2020).
[123] M. Zhou, E. Ott, T. M. Antonsen, and S. M. Anlage, Chaos: An Interdisciplinary
Journal of Nonlinear Science 29, 033113 (2019).
196
[124] S. Hemmady, X. Zheng, J. Hart, T. M. Antonsen, E. Ott, and S. M. Anlage,
Physical Review E 74, 036213 (2006).
[125] S. Ma, T. M. Antonsen, E. Ott, and S. M. Anlage, (2021),
10.48550/arxiv.2112.05306, arXiv:2112.05306 .
[126] O. Hul, S. Bauch, P. Pakon?ski, N. Savytskyy, K. Z?yczkowski, and L. Sirko,
Physical Review E 69, 056205 (2004).
[127] A. M. Mart??nez-Argu?ello, A. Rehemanjiang, M. Mart??nez-Mares, J. A. Me?ndez-
Bermu?dez, H.-J. Sto?ckmann, and U. Kuhl, Physical Review B 98, 075311 (2018).
[128] A. Rehemanjiang, M. Richter, U. Kuhl, and H.-J. Sto?ckmann, Physical Review
Letters 124, 116801 (2020).
[129] J. Lu, J. Che, X. Zhang, and B. Dietz, Physical Review E 102, 022309 (2020).
[130] C. Dembowski, B. Dietz, H.-D. Gra?f, A. Heine, F. Leyvraz, M. Miski-Oglu,
A. Richter, and T. H. Seligman, Physical Review Letters 90, 014102 (2003).
[131] C. Dembowski, B. Dietz, T. Friedrich, H.-D. Gra?f, H. L. Harney, A. Heine,
M. Miski-Oglu, and A. Richter, Physical Review E 71, 046202 (2005).
[132] U. Kuhl, R. Ho?hmann, J. Main, and H.-J. Sto?ckmann, Physical Review Letters
100, 254101 (2008).
[133] S. Shinohara, T. Harayama, T. Fukushima, M. Hentschel, T. Sasaki, and E. E.
Narimanov, Physical Review Letters 104, 163902 (2010).
[134] B. Dietz and A. Richter, Physica Scripta 94, 014002 (2019).
[135] B. Dietz, A. Heusler, K. H. Maier, A. Richter, and B. A. Brown, Physical Review
Letters 118, 012501 (2017).
[136] Y. Aure?gan and V. Pagneux, Acta Acustica united with Acustica 102, 869 (2016).
[137] M. Bia?ous, V. Yunko, S. Bauch, M. ?awniczak, B. Dietz, and L. Sirko, Physical
Review Letters 117, 144101 (2016).
[138] B. Dietz, V. Yunko, M. Bia?ous, S. Bauch, M. ?awniczak, and L. Sirko, Physical
Review E 95, 052202 (2017).
[139] J. Che, J. Lu, X. Zhang, B. Dietz, and G. Chai, Physical Review E 103, 042212
(2021).
[140] S. Johnson, A. Mekis, S. Fan, and J. Joannopoulos, Computing in Science
Engineering 3, 38 (2001).
[141] J. D. Joannopoulos, P. R. Villeneuve, and S. Fan, Nature 386, 143 (1997).
197
[142] J. D. Joannopoulos, S. G. Johnson, J. N. Winn, and R. D. Meade, Photonic Crystals
(Princeton University Press, 2008).
[143] K. Lai, T. Ma, X. Bo, S. Anlage, and G. Shvets, Scientific Reports 6, 28453 (2016),
1601.01311 .
[144] S. Ma, B. Xiao, Y. Yu, K. Lai, G. Shvets, and S. M. Anlage, Physical Review B
100, 085118 (2019).
[145] S. Ma and S. M. Anlage, Applied Physics Letters 116, 250502 (2020).
[146] Y. G. N. Liu, P. S. Jung, M. Parto, D. N. Christodoulides, and M. Khajavikhan,
Nature Physics 17, 704 (2021).
[147] W. Maimaiti, B. Dietz, and A. Andreanov, Physical Review B 102, 214301 (2020).
[148] W. Zhang and B. Dietz, Physical Review B 104, 064310 (2021).
[149] M. Reisner, M. Bellec, U. Kuhl, and F. Mortessagne, Optical Materials Express
11, 629 (2021).
[150] S. Bittner, B. Dietz, M. Miski-Oglu, A. Richter, C. Ripp, E. Sadurn??, and W. P.
Schleich, Physical Review E 87, 042912 (2013).
[151] A. Mekis, J. C. Chen, I. Kurland, S. Fan, P. R. Villeneuve, and J. D. Joannopoulos,
Physical Review Letters 77, 3787 (1996).
[152] T. Kottos and U. Smilansky, Physical Review Letters 79, 4794 (1997).
[153] T. Kottos and U. Smilansky, Annals of Physics 274, 76 (1999).
[154] G. Casati, F. Valz-Gris, and I. Guarnieri, Lettere al Nuovo Cimento 28, 279 (1980).
[155] A. Kudrolli, V. Kidambi, and S. Sridhar, Physical Review Letters 75, 822 (1995).
[156] S.-H. Chung, A. Gokirmak, D.-H. Wu, J. S. A. Bridgewater, E. Ott, T. M.
Antonsen, and S. M. Anlage, Physical Review Letters 85, 2482 (2000).
[157] U. Kuhl, E. Persson, M. Barth, and H.-J. Sto?ckmann, The European Physical
Journal B 17, 253 (2000).
[158] Y. Hlushchuk, L. Sirko, U. Kuhl, M. Barth, and H.-J. Sto?ckmann, Physical Review
E 63, 046208 (2001).
[159] L. Kaplan, Physical Review E 64, 036225 (2001).
[160] H. Alt, H. D. Gra?f, H. L. Harney, R. Hofferbert, H. Lengeler, A. Richter, P. Schardt,
and H. A. Weidenmu?ller, Physical Review Letters 74, 62 (1995).
[161] O. Hul, P. S?eba, and L. Sirko, Physical Review E 79, 066204 (2009).
198
[162] F. D. M. Haldane and S. Raghu, Physical Review Letters 100, 013904 (2008).
[163] M. Z. Hasan and C. L. Kane, Reviews of Modern Physics 82, 3045 (2010).
[164] Y. E. Kraus, Y. Lahini, Z. Ringel, M. Verbin, and O. Zilberberg, Physical Review
Letters 109, 106402 (2012).
[165] S. Raghu and F. D. M. Haldane, Physical Review A 78, 033834 (2008).
[166] B.-y. Xie, H.-F. Wang, X.-y. Zhu, M.-H. Lu, Z. D. Wang, and Y.-f. Chen, Optics
Express 26, 24531 (2018).
[167] F. Zhang, Science 358, 1075 (2017).
[168] J. Hou, Z. Li, X.-W. Luo, Q. Gu, and C. Zhang, (2018), arXiv:1808.06972 .
[169] C. L. Kane and E. J. Mele, Physical Review Letters 95, 226801 (2005).
[170] B. A. Bernevig and S.-C. Zhang, Physical Review Letters 96, 106802 (2006).
[171] X.-L. Qi and S.-C. Zhang, Reviews of Modern Physics 83, 1057 (2011).
[172] A. B. Khanikaev and G. Shvets, Nature Photonics 11, 763 (2017).
[173] Z. Wang, Y. D. Chong, J. D. Joannopoulos, and M. Soljac?ic?, Physical Review
Letters 100, 013905 (2008).
[174] Z. Wang, Y. Chong, J. D. Joannopoulos, and M. Soljac?ic?, Nature 461, 772 (2009).
[175] S. A. Skirlo, L. Lu, Y. Igarashi, Q. Yan, J. Joannopoulos, and M. Soljac?ic?, Physical
Review Letters 115, 253901 (2015).
[176] S. Longhi, D. Gatti, and G. D. Valle, Scientific Reports 5, 13376 (2015).
[177] L. Lu, J. D. Joannopoulos, and M. Soljac?ic?, Nature Physics 12, 626 (2016).
[178] Z.-G. Chen, J. Mei, X.-C. Sun, X. Zhang, J. Zhao, and Y. Wu, Physical Review A
95, 043827 (2017).
[179] B. Yang, H. Zhang, T. Wu, R. Dong, X. Yan, and X. Zhang, Physical Review B
99, 045307 (2019).
[180] M. Hafezi, E. A. Demler, M. D. Lukin, and J. M. Taylor, Nature Physics 7, 907
(2011).
[181] M. Hafezi, S. Mittal, J. Fan, A. Migdall, and J. M. Taylor, Nature Photonics 7,
1001 (2013).
[182] W. Gao, M. Lawrence, B. Yang, F. Liu, F. Fang, B. Be?ri, J. Li, and S. Zhang,
Physical Review Letters 114, 037402 (2015).
199
[183] S. Mittal, J. Fan, S. Faez, A. Migdall, J. M. Taylor, and M. Hafezi, Physical Review
Letters 113, 087403 (2014).
[184] D. Leykam, S. Mittal, M. Hafezi, and Y. D. Chong, Physical Review Letters 121,
023901 (2018).
[185] A. B. Khanikaev, S. Hossein Mousavi, W.-K. Tse, M. Kargarian, A. H. MacDonald,
and G. Shvets, Nature Materials 12, 233 (2013).
[186] W.-J. Chen, S.-J. Jiang, X.-D. Chen, B. Zhu, L. Zhou, J.-W. Dong, and C. T. Chan,
Nature Communications 5, 5782 (2014).
[187] X. Cheng, C. Jouvaud, X. Ni, S. Hossein Mousavi, A. Z. Genack, A. B. Khanikaev,
S. H. Mousavi, A. Z. Genack, and A. B. Khanikaev, Nature Materials 15 (2016).
[188] A. P. Slobozhanyuk, A. B. Khanikaev, D. S. Filonov, D. A. Smirnova, A. E.
Miroshnichenko, and Y. S. Kivshar, Scientific Reports 6, 22270 (2016).
[189] B. Xiao, K. Lai, Y. Yu, T. Ma, G. Shvets, and S. M. Anlage, Physical Review B
94, 195427 (2016).
[190] Y. Yang, Y. F. Xu, T. Xu, H.-X. Wang, J.-H. Jiang, X. Hu, and Z. H. Hang, Physical
Review Letters 120, 217401 (2018).
[191] X. D. Chen, W. M. Deng, J. C. Lu, and J. W. Dong, Physical Review B 97, 184201
(2018).
[192] J. Noh, S. Huang, K. P. Chen, and M. C. Rechtsman, Physical Review Letters 120,
063902 (2018).
[193] F. Gao, H. Xue, Z. Yang, K. Lai, Y. Yu, X. Lin, Y. Chong, G. Shvets, and B. Zhang,
Nature Physics 14, 140 (2017).
[194] R. G. Gladstone, M. Jung, and G. Shvets, (2018), arXiv:1809.02819 .
[195] X.-D. Chen, F.-L. Shi, H. Liu, J.-C. Lu, W.-M. Deng, J.-Y. Dai, Q. Cheng, and
J.-W. Dong, Physical Review Applied 10, 044002 (2018).
[196] K. Fang, Z. Yu, and S. Fan, Physical Review Letters 108, 153901 (2012).
[197] K. Fang, Z. Yu, and S. Fan, Nature Photonics 6, 782 (2012).
[198] M. C. Rechtsman, J. M. Zeuner, A. Tu?nnermann, S. Nolte, M. Segev, and
A. Szameit, Nature Photonics 7, 153 (2013).
[199] M. C. Rechtsman, J. M. Zeuner, Y. Plotnik, Y. Lumer, D. Podolsky, F. Dreisow,
S. Nolte, M. Segev, and A. Szameit, Nature 496, 196 (2013).
[200] T. Ma and G. Shvets, Physical Review B 95, 165102 (2017).
[201] Z. Wang and S. Fan, Optics Letters 30, 1989 (2005).
200
[202] W. Qiu, Z. Wang, and M. Soljac?ic?, Optics Express 19, 22248 (2011).
[203] S. R. Zandbergen and M. J. A. de Dood, Physical Review Letters 104, 043903
(2010).
[204] S. Bittner, B. Dietz, M. Miski-Oglu, P. Oria Iriarte, A. Richter, and F. Scha?fer,
Physical Review B 82, 014301 (2010).
[205] U. Kuhl, S. Barkhofen, T. Tudorovskiy, H.-J. Sto?ckmann, T. Hossain, L. de Forges
de Parny, and F. Mortessagne, Physical Review B 82, 094308 (2010).
[206] T. Ma, A. B. Khanikaev, S. H. Mousavi, and G. Shvets, Physical Review Letters
114, 127401 (2015).
[207] D. Polder, Phys. Rev. 73, 155 (1948).
[208] S. W. Mung and W. S. Chan, IEEE Microwave Magazine 20, 55 (2019).
[209] K. Fang, J. Luo, A. Metelmann, M. H. Matheny, F. Marquardt, A. A. Clerk, and
O. Painter, Nature Physics 13, 465 (2017).
[210] Z. Shen, Y.-L. Zhang, Y. Chen, F.-W. Sun, X.-B. Zou, G.-C. Guo, C.-L. Zou, and
C.-H. Dong, Nature Communications 9, 1797 (2018).
[211] X. Ni, M. Weiner, A. Alu?, and A. B. Khanikaev, Nature Materials 18 (2018),
10.1038/s41563-018-0252-9.
[212] H. Jia, R. Zhang, W. Gao, Q. Guo, B. Yang, J. Hu, Y. Bi, Y. Xiang, C. Liu, and
S. Zhang, Science 363, 148 (2019).
[213] R. Fleury, D. L. Sounas, C. F. Sieck, M. R. Haberman, and A. Alu?, Science 343,
516 (2014).
[214] S. Maayani, R. Dahan, Y. Kligerman, E. Moses, A. U. Hassan, H. Jing, F. Nori,
D. N. Christodoulides, and T. Carmon, Nature 558, 569 (2018).
[215] A. B. Khanikaev, R. Fleury, S. H. Mousavi, and A. Alu?, Nature Communications
6 (2015), 10.1038/ncomms9260.
[216] U. Stoffregen, J. Stein, H.-J. Sto?ckmann, M. Kus?, and F. Haake, Physical Review
Letters 74, 2666 (1995).
[217] Y. Kang, X. Ni, X. Cheng, A. B. Khanikaev, and A. Z. Genack, Nature
Communications 9, 3029 (2018), arXiv:1804.08707 .
[218] L.-H. Wu and X. Hu, Physical Review Letters 114, 223901 (2015).
[219] C. He, X.-C. Sun, X.-P. Liu, M.-H. Lu, Y. Chen, L. Feng, and Y.-F. Chen,
Proceedings of the National Academy of Sciences 113, 4924 (2016).
201
[220] X.-C. Sun, C. He, X.-P. Liu, Y. Zou, M.-H. Lu, X. Hu, and Y.-F. Chen, Crystals 9,
137 (2019).
[221] L. Chen, S. M. Anlage, and Y. V. Fyodorov, Physical Review E 103, L050203
(2021), 2101.08335 .
[222] Computer Simulation Technology, ?Weblink,? .
[223] R. Gunnarsson and M. Backstrom, in 2014 International Symposium on
Electromagnetic Compatibility (IEEE, 2014) pp. 169?174.
[224] B. D. Addissie, J. C. Rodgers, and T. M. Antonsen, in 2015 IEEE Metrology for
Aerospace (MetroAeroSpace) (IEEE, 2015) pp. 214?219.
202