Experimental Realization of Anti-Unitary Wave-Chaotic Photonic Topological Insulator Graphs Showing Kramers Degeneracy and Symplectic Ensemble Statistics

Working in analogy with topological insulators in condensed matter, photonic topological insulators (PTI) have been experimentally realized, and protected electromagnetic edge-modes have been demonstrated in such systems. Moreover, PTI technology also emulates a synthetic spin-1/2 degree of freedom (DOF) in the reflectionless topological modes. The spin-1/2 DOF is carried by Quantum Valley Hall (QVH) / Quantum Spin Hall (QSH) interface modes created from the bianisotropic meta waveguide (BMW) platform, and realized both in simulation and experiment. We employ the PTI setting to build an ensemble of wave chaotic 1D metric graphs that display statistical properties consistent with Gaussian Symplectic Ensemble (GSE) statistics. The two critical ingredients required to create a physical system in the GSE universality class, the half-integer-spin DOF and preserved time-reversal invariance, are clearly realized in the QVH/QSH interface modes. We identify the anti-unitary T-operator for the PTI Hamiltonian underlying our experimental realization. An ensemble of PTI-edgemode metric graphs are proposed and experimentally demonstrated. We then demonstrate the Kramers degeneracy of eigenmodes of the PTI-graph systems with both numerical and experimental studies. We further conduct spectral statistical studies of the edgemode graphs and find good agreement with the GSE theoretical predictions. The PTI chaotic graph structures present an innovative and easily extendable platform for continued future investigation of GSE systems.

In the 1950s, it was proposed that Random Matrix Theory (RMT) describes the energy level spacing statistics of large nuclei with various symmetries [1][2][3].Based on the BGS (Bohigas-Giannoni-Schmit) conjecture, RMT was later applied to describe the statistical properties of a broader range of the quantum/wave properties of systems that show chaos in the classical limit [4][5][6].The Dyson three-fold way [2,7] classifies the statistics of chaotic systems in three different universality classes.Systems that show time-reversal invariance (TRI) are described by the Gaussian orthogonal ensemble (GOE), the statistics of systems showing broken timereversal invariance are described by the Gaussian unitary ensemble (GUE), and the statistics of systems with TRI and anti-unitary symmetry show Gaussian symplectic ensemble (GSE) statistics [2,3,[8][9][10][11].Both GOE and GUE statistics were observed in a wide variety of chaotic systems, and in particular, microwave billiard systems [12][13][14][15][16][17][18][19].Complex metric graphs made with standard microwave coaxial cables have been shown to have statistical properties described by the GOE and GUE random matrix ensembles to a good approximation, although deviations are seen in some properties [20][21][22][23][24][25].
Due to the lack of a photonic analog of a spin-1/2 degree of freedom in polarization-preserving wave propagation, the realization of a GSE system was not observed for a long time in classical wave systems.Recently, a microwave experiment utilizing a carefully designed cable * skma@umd.edunetwork shows clear evidence of GSE statistics [26,27].This work relies on a specially engineered construction involving two geometrically identical sub-graphs, along with the creation of phase shifts inside and between the sub-graphs that enables an anti-unitary symmetry for the eigenmodes.The secular matrices of the network, written based on Kirchhoff's laws, directly resemble the form of the Hamiltonian of a spin-1/2 system, and a demonstration of the GSE energy level spacing statistics is demonstrated.However, this method lacks generality since delicate tricks are required: the creation of two sub-graphs with exactly equal cable lengths and carefully tuned 0 and π phase shifts for waves traveling between the two sub-graphs [26][27][28][29][30][31][32].As far as we know, there still awaits a physical realization that displays a clear and recognizable analog of full spin-1/2 physics in an intrinsic manner for the excitations of the system.In particular, we seek a physical realization in which eigenfunction and other properties can be studied in a manner similar to those employed for other wave/quantum chaotic systems, such as billiards.
Unprecedented wave phenomena have been introduced with the creation of photonic topological insulator (PTI) technologies [33,34].One of these phenomena is the synthetic spin degree of freedom (DOF) emulation in the photonic context.In particular, a pair of topologically protected modes with opposite effective spins has emerged.Mesoscopic condensed matter physics concepts, including the quantum Hall (QH), quantum spin Hall (QSH), and quantum valley Hall effects (QVH), have been demonstrated as classical analogs in the photonic context in a series of landmark papers [35][36][37][38][39][40].Researchers have emulated nontrivial topological edgemodes in many photonic structures [41][42][43][44].The bi-anisotropic metawaveguide (BMW) is one type of PTI system which can host all QH, QSH, and QVH analog phases for electromagnetic waves as perturbations of one basic structure [33,40,[45][46][47].The remarkable topologically-protected transmission properties of PTI have excited a wide range of real-life applications, for example, delay line devices, isolators, and circulators [33,46,48].However, the spin-1/2-like degree of freedom of the topological modes is not as widely exploited.Here we aim to create a microwave realization of a symplectic Hamiltonian system by using BMW structures with both QSH and QVH domains [40,49].Because the QSH-QVH interface waveguide has TRI and can host edgemodes with opposite synthetic spin DOFs, such a structure is an ideal building block for realizing complex structures that fall within the GSE universality class.We have recently demonstrated a chaotic 1D graph with topologically trivial photonic crystal defect waveguides and found good agreement with the corresponding RMT theoretical predictions [50].Our next step is to consider a chaotic graph structure where the BMW QSH-QVH interface waveguides realize the graph bonds, and represent the first use of topological waveguides for wave chaos studies.
In this paper, we first present the anti-unitary Toperator for the BMW Hamiltonian, directly linked to the Kramers degeneracy.We then introduce the design of the QSH-QVH eigenmode topological graph system.We present the clear observation of double (Kramers) degeneracy using numerical tools from photonic band structure studies.Electrically large metric graphs based on QSH-QVH interface modes are then studied for their mode-spacing statistics.We also conducted transmission measurements of the experimentally realized graph structure and find evidence of the Kramers degeneracy in the eigenmodes of the PTI graph.

II. II. ANTI-UNITARY T-OPERATOR
We are interested in using the QSH-QVH interface BMW waveguides as the building block of chaotic graphs that fall in the Gaussian symplectic ensemble (GSE) universality class, 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 (satisfying T 2 = −1).For a time-reversal invariant system, T 2 = −1 leads to |ψ⟩ and T |ψ⟩ being two different states with the same energy, creating what is known as a Kramers degeneracy.Other than a limited number of works [51][52][53], the derivation of the anti-unitary condition T 2 = −1 is not widely discussed for PTI systems, and is not discussed at all in the literature for the BMW system, to our knowledge.
As introduced before, researchers have demonstrated experimentally that the BMW system can host different topological domains in one composite structure [40,47,54].The BMW structure starts with a 2D hexagonal photonic crystal (PC) sandwiched between two con-ducting plates.The physical dimension of the lattice is carefully engineered so that the TE and TM modes are degenerate at the Dirac points (K and K ′ points) of the PC band structure.At the Dirac points, a synthetic spin-1/2 DOF is emulated by the orthogonal linear combinations of the TE and TM modes.For QSH, QH, and QVH BMW systems, the effective Hamiltonian near the Dirac points is written as [45] where τ , s, and σ are the Pauli matrices acting on the space of K and K ′ valley, the spin, and the Dirac degeneracy, respectively.The quantity ω D is the frequency at the Dirac points, and v D is the group velocity for both the degenerate TE and TM modes near the Dirac point.Note that the Kronecker product (s σ ≡ s ⊗ σ) shorthand is used.The quantity δk = (δk x , δk y ) refers to the distance to the Dirac point in k-space.The second term in Eqn. 1 represents one of several types of perturbations that can be made to the hexagonal PC lattice, in this case creating an effective spin-orbit interaction, discussed in more detail below.The BMW Hamiltonian has a similar form to the Kane-Mele Hamiltonian, which is used to model graphene-like topological insulators [45].The eigenvectors are written as Ψ , where the superscript ↑ (↓) refers to the spin-up (-down) index and the subscript represents the valley index [45].The second term introduces bianisotropy and contains ∆ em which is the overlap integral of unperturbed TE and TM modes (subscripts e, m) inside the perturbed unitcell volume.Different topological orders are then achieved by introducing different perturbations to the basic unit cell.For the two topological domains used in this paper, the QSH and QVH phases are realized by breaking the out-of-plane mirror symmetry, and the in-plane rotation symmetry, respectively [49].These symmetry-breaking perturbations will lift the mode degeneracy at the Dirac point, and create a full band gap in the PC band structure.The detailed introduction of the BMW system Hamiltonian theory is not the goal here, and we kindly direct the readers to find more detailed discussions in Refs.[45,49].
Here we propose an effective T-operator T b = iτ x s y σ 0 C for the BMW system, where C is the complex conjugation operator (C i = −i C), and τ is the Pauli matrix that acts on the valley degree of freedom (K and K ′ ).We want to first point out that the proposed T-operator, T b = iτ x s y σ 0 C, clearly satisfies the anti-unitary criterion T 2 b = −1, leading to the Kramers degeneracy.The physical meaning of the operator T b lies in that it flips momentum (valley) and the spin degrees of freedom of the edgemode, expressed as satisfies this criterion.

III. III. PTI-GRAPH DESIGN
Recent works have reported that the BMW-PTI system can host both QSH and QVH domains inside one composite structure [40,54].Here we present a new version of the QVH and QSH BMW-PTI structure and experimentally construct chaotic graph systems which consume ∼ 1000 unit cells.We designed the structure to operate around 24.5GHz to shrink down the physical size of the graph structure based on the electromagnetic scaling law.As shown in Fig. 1 (a), the QVH region unit cell consists of two identical metallic prisms with length l v = 4.30 mm, height h v = 4.32 mm, and the middle air gap g v = 0.69 mm.Each prism is firmly attached to one of the two parallel metallic plates.The QSH region unit cell is a metallic cylinder with diameter d s = 3.22 mm and height h s = 8.3 mm.Both regions are placed between two metallic surfaces with a spacing of h 0 = 9.33 mm in the vertical direction.The cylinder is firmly attached to one metallic plate and a gap of thickness g s = h 0 − h s = 1.01 mm is made with the second metallic plate.The lattice constant of the triangular BMW lattice is a 0 = h 0 = 9.33 mm.These physical dimensions are carefully engineered to ensure the QSH and QVH bulk bandgaps are matched in both center frequency and bandwidth.This bandgap occurs when the guided wavelength of waves in the interface region is comparable to a 0 .Fig. 1 (b) illustrates how the QSH and QVH rods are mounted between the two metallic covering plates.
We examine the existence of edgemodes in the QSH-QVH interface waveguide through the so-called supercell simulation in COMSOL, as shown in Fig. 1 (c).The supercell structure consists of a slice of QSH and QVH rods (8 rods for each domain in our simulation).The Floquet boundary condition is applied to the two long sides to simulate an infinitely long waveguide structure.It is shown that the QSH-QVH waveguide can host two counter-propagating edgemodes between 23.8 ∼ 24.8 GHz, and the two modes have opposite spin and valley DOFs.Near the Dirac point, the synthetic spin-1/2 DOF is simulated with the plus and minus linear combinations of the TE and TM degenerate modes [45].These propagating modes show up inside the bulk bandgap region in Fig. 1 (c), with opposite propagation directions.In real-space, the edgemodes and bulkmodes also show distinct eigenfunction distributions.The edgemodes reside strictly at the interface of different topological domains (marked by the red dashed line in Fig. 1 (b) and inset of (c)), and the mode amplitude decays exponentially into the surrounding domain regions.The eigenfunction dis-tribution of the bulkmodes reside in the interior of one or more topological domains.The clear distinction in the mode eigenfunction is another indicator for distinguishing between edgemodes and bulkmodes.
The field of quantum chaos originated from nuclear physics studies, where the spectral statistics of the Hamiltonian of large nuclei are found to be well described by the statistics of the eigenvalues of random matrices [2,3,[55][56][57].It was later conjectured that all complicated quantum/wave systems that have chaotic classical/ray counterparts possess universal statistical properties described by RMT [5].Such chaotic phenomena have been manifested in wave systems with different dimensionalities [17-19, 58, 59].The quantum graph system is an inter-connected network of transmission lines that meet at nodes/junctions [60][61][62].Waves propagate on the network according to the telegrapher equations which resemble the Schrodinger equation, and the conservation of probability current imposes a Neumann boundary condition at each node.Previous work has shown that photonic crystal defect waveguides with trivial topology can be used as quantum graphs, and display many of the same statistical properties of more conventional cable graphs [50].The proposed PTI symplectic quantum graph structure can be considered a square pyramid realized in the 2D plane (see the inset in Fig. 2(c)).Each bond of the graph is built with a QSH-QVH interface waveguide.The interface waveguides meet at nodes with valance 2 or 4. The graph will be electrically large (i.e. each bond will be many wavelengths long) and create multiple paths between any two nodes, so that the scattering will be complex.
The open-plate view of the QSH-QVH PTI-graph experimental realization is shown in Fig. 1 (d).Because of the relatively high operating frequency (∼ 24.5 GHz), we can create a graph structure with 60 lattice periods on a side but with moderate overall size, on the scale of 50 cm × 50 cm.To prevent the bending and flexure of the covering plates, we have also installed several support rods (with diameter d s and height h 0 ) at several locations in the horizontal plane in the bulk regions, away from the graph bonds.
The closed graph has numerous eigenfrequencies and corresponding eigenfunctions.When the graph is connected to one or more outside leads at discrete ports, one can define a frequency-dependent scattering matrix S. The ports are installed at the graph nodes where multiple bonds meet, and each port consists of a coaxial cable mounted perpendicular to one plate whose center conductor extends ∼ 2 mm into the graph below the top metal plate.The scattering matrix is measured as a function of frequency using a Keysight N5242A vector network analyzer.
An ensemble of PTI graphs can be realized by changing each graph bond's length, similar to studies of topologically-trivial PC [50] and waveguide-based [63] graphs.The losses for edge modes propagating in the structure are dominated by conductor loss in the plates, prisms, and cylinders, and their points of contact with the plates.A small amount of radiation loss also exists out of the 4 open edges of the finite-size structure.The conductive loss of the structure can be approximated by a uniform attenuation η, and its value can be reduced by cooling down the graph from room temperature using a dilution refrigerator.

IV. IV. NUMERICAL SIMULATION STUDIES
We first present the eigenmode study of closed PTI graphs using the numerical simulation software (CST).The simulated graphs have zero uniform loss, although loss is present in the simulation due to leakage through the boundaries of the finite-size simulation domain, where radiating boundary conditions are assumed.Inside the bulk bandgap frequency range (24 ∼ 24.8 GHz as shown in Fig. 1 (b)), we observe only doubly degenerated edgemodes, as shown in Fig. 2 (a).At the lower edge of the bulk bandgap (∼ 23.5 GHz), a co-existence of the degenerate edgemodes and the higher spectral-density bulkmodes is observed, where the spectrum becomes nearly continuous and the doubly degenerate modes mix with the bulk modes [64].Thus the double (Kramers) degeneracy of the QVH-QSH edge modes is clearly demonstrated with numerical tools.
We next study the nearest neighbor level spacing statistics of the proposed pyramidal graphs.The nearest neighbor level spacing distribution is an excellent diagnostic of the universality class that governs the fluctuations present in the properties of a wave chaotic system [2,3,10,14].The eigenvalue simulations are conducted using the software HFSS running on the UMD DeepThought2 HPC cluster.The simulated graph has the area of ∼ 35 × 60 photonic crystal lattice constants, where ∼ 12 levels (2 degenerate modes per level) are found in the topological edgemode bandwidth.For the nearest neighbor level spacing study, the quantity s i = e i+1 − e i , where e i is the graph eigenmode level and the subscript i refers to the mode index.The quantity s i is then normalized as s i /⟨s i ⟩ where ⟨s i ⟩ is the averaged level spacing.We simulate 30 different configurations of the PTI graph and obtain 342 level-spacing values.In Fig. 2 (c) inset we show one configuration of the pyramidal graph.New graph configurations are created by varying the lengths of graph bonds, which is implemented by adding/removing zig-zag paths on graph bonds.Due to the topological protection, the waveguide modes can travel through these zig-zag paths without back-scattering.By varying the length of the bonds, we are able to change the spectrum of resonant modes of the entire graph structure.
In the quantum chaos community, it is of interest to study the statistics of the system level-spacing s through its distribution function, P (s), and its integrated value I(s) = s 0 P (s ′ )ds ′ .By comparing the functions P (s) and I(s) with the theoretical predictions of GOE, GUE, and   The ensemble-averaged consecutive mode spacing ratios of the photonic topological insulator graph (PTI-Graph) system, the theoretical predictions for the GOE, GUE, and GSE systems [66], and the mode spacing ratios of the topologically-trivial photonic crystal graph (PC-Graph) [50].

PTI-Graph GSE GUE GOE PC-
GSE systems, one is able to classify the symmetry group of the system of interest [2,3,14].As shown in Fig. 2 (b) and (c), the level spacing statistics of the pyramidal graph agree reasonably well with the GSE statistics at larger s values (s > 0.5), and show a distinctly different shape as compared to the GOE and GUE systems, especially in terms of the peak and tail of P (s).For small spacing values (s < 0.5), we find that the graph system deviates from the GSE class prediction.We attribute this to the small statistical ensemble (∼ 12 levels per realization), and the proximity of many bulk modes to the bandgap.We have examined the statistics of nearest neighbor spacings in finite-size GOE matrices in which alternating eigenvalues are ignored to replicate GSE statistics, as demonstrated before [3,65].We find that taking ten 30 × 30 GOE matrices, yielding 14 nearest-neighbor pairs, show statistical properties similar to our data (see details in Supp.Mat.).
To augment the mode spacing distribution test, we also investigated the method of consecutive mode spacing ratios, where one computes two parameters r i = si si−1 and ri = min r i , 1  ri [66].The mode spacing ratio ensembleaveraged data are reported in Table.I.For the PTI graph system, the averaged values of ⟨r⟩ = 1.05 and ⟨r⟩ = 0.89 are closer to the GSE theoretical predictions than either the GUE or GOE predictions.Note that we also include the mode spacing ratio values of the GOE graph system consisting of topologically-trivial photonic crystal waveguides in the last column [50].It is clear that a GOE graph system mode spacing ratios are at the other extreme of values, and closer to the GOE theoretical prediction values.From this mode spacing study of numerical realizations of the PTI graph system we conclude that it is more consistent with GSE statistics than any other Gaussian ensemble.

V. V. EXPERIMENTAL STUDIES
With the QSH-QVH waveguide carefully designed by means of numerical simulations, we fabricated a physical realization of the chaotic PTI-graph structure, as shown in Fig. 1 (d).The PTI rods and the two covering flat plates are made of aluminum 1070 alloy (aluminum > 99.7% and iron < 0.25%).The structure is designed to be compatible with low-temperature measurements in a dilution refrigerator (DR) so that the entire graph can undergo uniform changes in temperature to systematically vary the uniform attenuation η of the edgemodes.We have installed numerous support rods to prevent the drooping/bending of the covering plates to ensure uniform air gaps throughout the structure.The structure is measured in a vertical geometry to minimize the effects of gravity on the graph properties.Note that the center frequency and the width of the bandgap are sensitive to the set of parameters h 0 , g v , and g s , and the triangular prism QVH rods on the two cover plates need to be carefully aligned when assembling the graph structure.It should be noted that small variations in h 0 , g v , and g s , and the small variations in the orientation of the QVH rods, have no effect on the topological properties of the PTI-guided waves as long as they do not close the bandgap [67].
We conduct 2-port S-parameter measurements using a Keysight N5242A PNA-X microwave network analyzer.Electric dipole antennas are installed at two of the graph nodes (shown as the stars in Fig. 1 (d)).The electric field inside the graph structure is coupled to the dipoles which allows one to probe the eigenmodes of the closed graphs.The lengths of the dipoles are chosen to be short so that the coupling is weak, aiming to minimize the perturbation to the modes.For room-temperature S-parameter measurements, the microwave network analyzer calibration planes are at the end of the dipole antenna input terminals.For low-temperature tests, we conduct uncalibrated S-parameter measurements because the cables inside the DR experience varying temperatures during the measurement and the calibrations cannot be performed at each temperature.
Our main objective in this experiment is to identify the Kramers degeneracy expected in the QSH-QVH interface graph modes.However, the degenerate modes are difficult to identify by simply studying the resonances in the S-parameter spectrum.Several methods to identify degenerate modes in microwave resonant systems have been developed in the past.For example, the perturbations caused by the ports used to measure a billiard can be used to identify two degenerate modes [68].By exciting a microwave cavity at two locations with a phase-shifted signal, one can also identify degenerate modes [69].In the graph structure with anti-unitary symmetry [26,27], the system can be perturbed away from the interference condition that underlies anti-unitarity, and observe the splitting of the Kramers degenerate modes.However, the topologically-protected nature of the QVH-QSH edge modes prevents such studies here.
To better examine the question of degenerate modes in the data, we compute the complex Wigner-Smith timedelay [70][71][72][73] which can be used to identify the location of the scattering matrix poles and zeros in the complex frequency plane [73,74].The quantity M is the number of measurement channels, and S(f ) refers to the full M ×M complex scattering matrix, measured as a function of frequency.We adopt the Heidelberg picture of scattering theory [75][76][77].In this picture, we regard the closed graph to be described by a Hamiltonian H, having real eigenvalues that describe the resonant frequencies of the isolated graph.Coupling of the graph to the outside world through M ports introduces a coupling matrix W whose matrix elements W i,j denote the coupling of mode i to port j.This creates an effective (non-Hermitian) Hamiltonian H ef f = H − iΓ W , where Γ W = πW W + .The complex eigenvalues of H ef f are the poles of the scattering matrix, written as f n −iΓ n , with the convention that Γ n > 0 for passive lossy systems.Because the losses are uniform in our system, we shall assume that the zeros of the S-matrix are simply the complex conjugates of the poles.The complex time delay τ W around an isolated group of N modes suffering uniform attenuation η > 0, as a function of frequency, can be fit by the following relationships [73]: (3) where the quantities η, f n − iΓ n and C are the uniform loss rate, the n th system complex resonant frequency (Smatrix pole), and a constant term to account for contributions to Re(τ W ) from far-away resonances, respectively.For a single isolated resonance mode in the PTIgraph τ W spectrum, we expect N = 2 because there are 2 degenerate modes for each energy level.Note that we have simplified the fitting process by assuming the zeros and poles of the S-matrix, written as z n and f n − iΓ n in Ref. [73], where n is the mode index, are complex conjugates of each other (i.e.z n = f n + iΓ n ).This assumption is suitable for the PTI-graph systems because it has only uniform loss and there is no lumped loss element present in the graphs [73,74].The value of η is dictated by resistive losses in the Aluminum plates and rods, and can be considered independent of frequency over the range of fits to Eqs. 3 and 4, in the bandgap of the PTI graph structures.
Based on Eqns. 3 and 4, we are able to fit the complex Wigner-Smith time delay τ W (f ) obtained from the scattering matrix data for the PTI eigenmodes.Two examples of the measured (blue dots) and the fitted (solid red line) complex τ W curves are shown in Fig. 3 for an isolated degenerate mode of the PTI graph.Note that these τ W (f ) curves have an asymmetric shape as a function of frequency.From the forms of Eqns. 3 and 4, it is clear that the asymmetric τ W (f ) curves must be the result of two extremely nearby modes.To further prove that the nearby modes are degenerate, we conduct high-quality numerical fitting of the measured data.The  d)) might just be two closely located but unrelated modes.However, the spacing between these two modes is on the order of ∼ 10 MHz, which is much smaller than the ∼ 70 MHz mean mode-spacing for a graph of this size.The mean-mode-spacing is computed by c 2Ltot = 0.073GHz, where c is the speed of light and L tot = 2.05m is the total length of all graph bonds.Note that the imaginary parts of the poles, as well as the uniform attenuation, are all much less than the mean spacing between edgemodes, allowing for a definitive determination of all the scattering singularities.Thus, in the case of some of the modes, we believe that we have identified the expected 2fold Kramers degeneracy through the measured complex time-delay spectrum.
For most resonant modes of the graph, the f n values of the degenerate mode pair appear to be almost identical so that the corresponding τ W (f ) curves do not have an asymmetric shape.However, in these situations, one can utilize the difference in the imaginary parts of the Smatrix poles (Γ n ) to reveal the 2-fold degeneracy as follows.Two nearly-degenerate S-matrix poles at The difference in imaginary parts is due to the different overlap of the modes with the measurement ports, and this difference is a generic feature of chaotic graphs (note that this assumption is the basis of the Random Coupling model for the statistical properties of complex scattering systems [16,18,78,79]).To find nearly degenerate modes, we shall exploit the properties of complex time delay (Eqns.3 and 4) by smoothly sweeping the uniform loss η value so that it sequentially passes through the values of both Γ 1 and Γ 2 .For example, one may strategically set the initial value (room temperature) of η to be greater than max[Γ 1 , Γ 2 ], and then systematically decrease η until η < min[Γ 1 , Γ 2 ].During this process, the value of η will first equal the Γ n value of one of the two degenerate modes (max[Γ 1 , Γ 2 ]).When this happens, the τ W (f ) spectrum will diverge as one of the denominators of Eqs. 3 and 4 go to zero when f = f n .When τ W (f ) diverges, the Re(τ W ) will change sign and the Im(τ W ) will show high absolute values just above and below resonance.Moreover, because one generically has Γ 1 ̸ = Γ 2 , one expects to find that the τ W (f ) curve 'diverges' twice during the process of systematically varying η.We call this process a Lazarus time-delay behavior, where varying η has made τ W 'die'verge twice to demonstrate the presence of an eigenmode degeneracy.This double-timedelay divergence will not take place if only one mode exists in the frequency range of interest.Hence, cooling the graph in a dilution refrigerator (DR) from room temperature to low-temperatures will allow us to gradually decrease the degree of uniform loss η by controlling the PTI structure global conductivity loss, while keeping the values of Γ 1 and Γ 2 approximately constant since they are dictated by the structure of the port.In this way, we can reveal the existence of nearly degenerate modes from the evolution of complex time delay as η is systematically varied by means of the temperature variation of the PTI graph.
To identify Kramers degenerate pairs of modes, we adopted the Lazarus method experimentally by cooling the graph from room temperature to ∼ 30 mK.The conductivity loss of the system, represented by η in Eqns. 3 and 4, will gradually and monotonically decrease during the cooldown process [80].The coupling at M = 2 ports is made weak enough such that max[Γ 1 , Γ 2 ] < η at room temperature (this is the case for both of the modes shown in Fig. 3).The graph structure is suspended vertically below the mixing chamber plate (mxc plate) of a BlueFors XLD-400 DR, as shown in Fig. 4 (a).The graph structure takes up nearly the entire diameter of the mxc plate and is a large thermal mass attached to the DR.We used braided copper wire straps (not shown) to aid the thermal connection between the graph and the mxc plate, and to make the temperature of the entire structure uniform.The residual resistivity ratio of Aluminum typically exceeds 10 3 [80] so a reduction of η by at least a factor of 30 is expected between room temperature and the base temperature of the DR.Note that the S-parameter measurement is not calibrated in this case, but this does not affect our ability to identify degenerate modes from complex time delays.
Here we present the analysis of measured τ W (f ) curves during the cooldown.Figures. 4 (b) and (c) show the measured τ W (f ) spectrum for a single pair of modes when the temperature changes from 180 K to 80 K.As discussed above, we expect a double divergence of the τ W (f ) as the temperature-dependent value of η is swept through Γ 1 and Γ 2 .To better illustrate the evolution of τ W under uniform loss variation, we also present a numerical example of similar character in Fig. 4 (d) and (e).The numerically-generated τ W (f ) data in Figs. 4 (d) and (e) is obtained using Eqs. 3 and 4, with parameters (f 1 , f 2 ) = (24.382,24.382) GHz, (Γ 1 , Γ 2 ) = (0.0054, 0.0057) GHz, and η varying from 0.0062 to 0.0050 GHz.As shown in both the measurement and the numerical example, when τ W diverges, we observe that Re(τ W ) changes sign and Im(τ W ) shows high magnitudes near resonance.A clear resemblance between the measurement and the numerical example is found.In Figs. 4 (f) and (g) we show the top view of the real and imaginary parts of the measured τ W (same data as in (b) and (c)).It is clear that only one dispersing feature of enhanced τ W (f ) exists throughout the measurement, as opposed to two distinct dispersing features, as one would expect for accidentally degenerate modes.This observation means that the frequencies of the two degenerate modes remain essentially identical as the PTI graph contracts and becomes less lossy through cooling.No other singular features in τ W (f ) are observed from this set of degenerate modes outside of this temperature range.Figs. 4 (b) to (e) show exactly the case where the 2-fold degeneracy can not be identified through a difference in the real part of S-matrix poles/zeros (f n ), but can be identified from the difference in imaginary parts (Γ n ) during uniform loss (η) variation.Effectively, we identify the presence of degenerate modes by sweeping 'imaginary frequency' through the S-matrix zeros.We note in passing that the divergence of Re[τ W ] is an indication that an S-matrix zero has been brought to the real frequency axis [73] and that coherent perfect absorption can be achieved by exciting the system with the corresponding eigenvector [81,82].

VI. DISCUSSION
Through the use of complex time delay and temperature variation of the uniform attenuation of the graph, we have identified the Kramers degeneracy in the proposed PTI graph structure.Our work presents a clear experimental method for identifying Kramers degeneracy within PTI systems.Additionally, our innovative use of complex time delay and variable uniform attenuation, leading to what we term 'Lazarus time-delay behavior', can potentially be applied to other fields of study.We also note that in addition to the demonstrated method of uniform attenuation tuning, the manipulation of other parameters, such as coupling strength, can lead to similar 'Lazarus' time delay divergent events.Future exploration of Kramers degeneracy in other contexts is thus enabled by our experiment.
The introduction of Kramers degeneracy to microwave and optical resonant systems plays a pivotal role in extending the advantages of PTI systems.The presence of Kramers degeneracy within PTI systems is intricately related to the system's band structure, endowing the modes with topological protection as long as the bandgap remains intact.This remarkable property enables PTI modes to travel seamlessly through sharply bent waveguides and propagate without reflection, making PTI systems invaluable for a variety of applications [33].
The QVH/QSH PTI graph incorporates the spin-1/2 degree of freedom directly into the wave excitation that propagates throughout the system.This is in contrast to previous GSE experimental studies which relied on connecting two GUE graphs in a particular way to create a composite system that had anti-unitary symmetry.As a result of the simpler and distributed nature of our GSE realization, our approach also enables studies of many other wave chaotic properties of GSE systems, such as eigenfunction statistics, for the first time.This will enable a complete investigation of the statistical properties of systems making up Dyson's three-fold way [2,7].

VII. CONCLUSION
To conclude, we have proposed a PTI-graph structure as a new platform for realizing wave chaotic systems that display symplectic Hamiltonian behavior in a natural and general manner.We present an effective time-reversal operator T b for the BMW system and show that it is antiunitary.A new composite QSH-QVH BMW interface structure is designed, fabricated, and employed as the building block for creating chaotic graphs in the Gaussian symplectic ensemble universality class.We applied numerical tools to study the proposed graphs and find clear evidence of Kramers degeneracy with eigenvalue simulations of the closed graph structure.We studied nearest-neighbor level spacing statistics by simulating an ensemble of such graphs.A relatively good agreement between the GSE prediction and the graph system data is found.We then realized an equivalent QSH-QVH BMW interface graph structure in Aluminum.In the experimental test, we identified the Kramers degeneracy through the asymmetric shape of the measured complex Wigner-Smith time delay τ W (f ) spectrum.With the help of variable-temperature measurements, we systematically varied the system uniform loss η and found clear double divergence in the measured τ W (f ), demonstrating the presence of two degenerate modes.Thus, theoretical, numerical, and experimental studies have accumulated clear evidence for Kramers degeneracy in this unique setting.Compared to the reflection-less waveguide transmission in PTI, its delicate emulation of a spin-1/2-like DOF of light is highly under-utilized in scientific explorations [33].This work serves as an example of how the synthetic spin DOF can be employed to address a long-standing puzzle in another field of physics, namely quantum chaos.
FIG. 1. a. Schematic of the unit cell structure that defines 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 and hs are the diameter and height of the QSH cylinder rod.b.The side-view of the QSH-QVH structure with an interface (red dashed line).The top and bottom metallic plates are separated by h0.All QSH rods for a given region are attached to one of the metallic plates, and the QVH rods are mounted on both plates, leaving an air gap.c.The photonic band structure of the QSH-QVH interface waveguide.Inset shows a schematic of the supercell simulation model (same structure as in b).The red dashed line marks the interface between the QSH (left) and the QVH (right) region.d. shows the top and bottom metallic plates and the rods/prisms attached to each one of them in the experimental realization of a GSE graph.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 spread out the weight of the top plate evenly, and provide a fixed and uniform separation between the plates, h0.The inset introduces a node structure in the pyramidal graph, where the QSH and QVH rods are colored in grey and blue, respectively.

FIG. 2 .
FIG. 2. a.The simulated eigenmode results of a closed 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.5 GHz), we observe a hybridization of the edgemode and bulk modes.b.The histogram approximation probability distribution of the PTI graph nearest-neighbor mode-spacing s.The prediction of the GOE, GUE, and GSE ensembles are plotted as solid lines for reference.c.The integrated probability function I(s) along with the theory predictions of the three groups.The inset shows the schematic of the square pyramidal graph.The black and blue lines represent the graph bonds and cables connected to the ports, respectively.

FIG. 3 .
FIG. 3. The real and imaginary parts of the measured complex Wigner time delay as a function of frequency for nearly degenerate pairs of modes around 24.20 and 24.46 GHz at room temperature, shown in (a, b) and (c, d), respectively.The fitting results are shown in solid red curves, and the fit parameters are given in the text.The deviations between data and the fits at the extreme edges of (a), (c), and (d) are due to the presence of neighboring levels, which are not included in the fit.

FIG. 4 .
FIG. 4. a shows the front view of the closed PTI-graph structure mounted to the DR mixing chamber plate.b and c show the real and imaginary part of the measured PTI-graph complex Wigner time delay τW versus temperature and frequency, respectively.d and e show the real and imaginary part of a numeric example of complex Wigner time delay τW versus varying uniform attenuation η and frequency, respectively.The yellow arrows mark the location of τW divergence in b through g. f and g show top views of the data in b and c, plotted on the same color scale.