{'fraction_non_alphanumeric': 0.036747630703571205, 'fraction_numerical': 0.014058871524508881, 'mean_word_length': 4.65719865080774, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Eigenfunction and eigenmode-spacing statistics in chaotic photonic crystal graphs

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 eigenmode-spacings and the statistics of the eigenfunctions. Here, we propose photonic crystal (PC) defect waveguide graphs as a new physical setting for chaotic graph studies. Photonic crystal waveguides possess a dispersion relation for the propagating modes which is engineerable. Graphs constructed by joining these waveguides possess junctions and bends with distinct scattering properties. 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.


I. INTRODUCTION
Wave-chaotic phenomena have been studied in various complex scattering systems, ranging from 1D graphs [1][2][3][4][5], 2D billiards [6][7][8][9][10][11][12] to 3D enclosures [13][14][15][16][17][18]. 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 [19,20]. The statistics of complex systems that show time-reversal invariance (TRI) are described by the Gaussian orthogonal ensemble (GOE) of random matrices, and 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 [9,16,21]. It was later conjectured by Bohigas, Giannoni, and Schmit that any system with chaotic dynamics in the classical limit would also have semi-classical wave properties whose statistics are governed by RMT associated with one of these three ensembles [22]. 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 [9,13,14,16,17,[23][24][25][26][27][28].
Chaotic microwave graphs support complex scattering phenomena despite their relatively simple structure [1][2][3][4][5]29]. We refer loosely to a chaotic graph as one in which there are multiple paths of differing lengths from one node to another, such that waves traveling along these paths interfere when arriving at nodes. A microwave graph structure can be realized with coaxial cables as graph bonds connected at T-junctions. The graph architecture allows for various useful circuit components (such as variable phase shifters and attenuators, as well as microwave circulators) to be incorporated into the structure [3][4][5]. 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 impact the spectral statistics of the graph [4,[30][31][32].
Here, we introduce an alternative type of chaotic graph built with photonic crystal (PC) waveguides. Photonic crystals find extensive use in semiconductor-based nanophotonic technology, where they act as low-loss waveguides for visible and infrared light [33]. 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 PCgraph systems are expected to fall into the GOE universality class of RMT. With numerical simulation tools, we conduct a series of statistical tests of chaotic photonic crystal graph systems including both eigenvalue and eigenfunction studies of closed graphs.
The unique properties of PC structures for quantum chaos studies are as follows. First, the bond lengths of PC-graphs can be altered by means of lithographic fabrication [33]. Second, the scattering properties of PCgraph nodes can be engineered by changing the rod properties and geometry. Moreover, PC defect waveguides and cavities are technologies that are widely used in the field of integrated nanophotonics. To our knowledge, PC waveguide 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 closedgraph 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.

II. 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 [34,35]. 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 [34,35]. 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 [34]. Recent photonic topological insulator (PTI) studies present an alternative form of PC waveguide using the interface between two different topological domains [36][37][38][39][40]. 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 [12,41,42], and from coupled resonant dielectric cylinders [3,43] studied previously.
The construction of the chaotic PC graph starts by building a square lattice with identical dielectric rods (the blue lattice in Fig. 1(a)). The lattice constant is a 0 = 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.  1(b). The diameter of the dielectric rod is d 0 = 0.36a 0 = 13.2mm and the height is h 0 = 0.1a 0 = 3.68mm. Because the PC lattice is thin in the vertical direction (zdirection), the waveguide modes considered here are polarized such that E z = 0, E x,y = 0, H z = 0 and H x,y = 0 (TM with respect to z), and these are the only modes that propagate for frequencies f < 12 GHz, well above the band gap of interest here. 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. 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 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 'T' connector (dashed). Inset: 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. long sides differ by a phase factor θ = k x a 0 . Fields on the opposite short sides 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 k x in the range [−π/a 0 , π/a 0 ]. As shown in Fig. 1(c), the defect waveguide modes have emerged inside of the bulk bandgap region from 2.5 ∼ 3.6GHz. The eigenfrequencies in the band-gap region are essentially real, as the absorbing boundaries are well separated from the defect creating the waveguide.
We see from Fig. 1(c) that the waveguide modes in the PC have a distinct dispersion relation (frequency versus k x ) that differs from that of a transmission line, ω = k x c. The dispersion relation is parabolic at the lower bandedge, as in a metallic waveguide. However, the group velocity vanishes at the upper band-edge, in contrast with that of a metallic waveguide for which there is no upper band-edge, rather just a boundary where new propagating modes are possible [29]. In Fig. 1(d), the |E z | 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.
In addition to the propagation on the graph bonds, it is also important to characterize the scattering properties of graph nodes [44]. 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 junctions with COMSOL and CST Microwave STUDIO. For the rightangle junctions, non-zero transmission is found only in the bulk bandgap region. As shown in Fig. 2(a), the transmission |S 21 | and the reflection |S 11 | for right angle bends vary systematically as functions of frequency. Thus, right-angle bends are vertices with a non-trivial scattering matrix in a graph. The PC-graph T-junction also presents a complex scattering profile over the entire bandgap region. As shown in Fig. 2(b), the magnitude of the reflection (transmission) coefficient deviates systematically from 1/3 (2/3) as a function of frequency.
It should be noted that for the PC graph junctions the sum of the squares of the displayed scattering coefficients is close to unity. This implies that at the frequencies of interest the photonic waveguides are operating in a single transverse mode. For an ideal transmission line there is no reflection at a bend. Also, at an ideal Tjunction the reflection coefficient |S 11 | = | − 1/3|, and the transmission coefficient is |S 12 | = |2/3|. The ideal scattering parameters apply to transmission lines whose transverse dimensions are much smaller than an axial wavelength (equivalently the frequency is well below the first cut-off frequency of non-TEM modes.). The calculated scattering properties of an air-filled coaxial (r o = 2 mm, r i = 0.86 mm) cable bend and T-junction are shown as well in Fig. 2, and, by contrast, they show minimal variation in the frequency range of interest. The cutoff frequency of the T E 11 mode in this cable is 33.34 GHz, well above the operating frequency. Consequently, the scattering parameters are close to the ideal values. Graphs constructed of metallic waveguides, as opposed to transmission lines, can be expected to have scattering parameters with characteristics similar to PC graphs for frequencies between the first and second lowest cutoff frequencies where single mode operation is possible [29]. These nontrivial junction scattering parameters 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 [45]. One may adjust the degree of wave localization of graph nodes by engineering the scattering properties of the PC waveguide junctions. 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 rightangle junctions). The four exterior faces of the parallelplate PC structure (see Fig. 3(a)) are assigned totally absorbing boundary conditions. The total length of the simulated graph is on the scale of ∼ 9 m, 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. We have created a statistical ensemble of chaotic PC graphs by changing the length of the bonds for a given graph topology.
A particular eigenmode solution of a graph (shown in Fig. 3(a)) is shown in Fig. 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.

III. EIGENMODE SPACING STATISTICAL ANALYSIS
In order to characterize the statistical properties of PC defect waveguide graphs, we start by conducting the nearest neighbor mode spacing analysis of the proposed PC graphs. An ensemble of 10 different graph realizations is 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 model is composed of lossless dielectric structures. However, the presence of radiating boundary conditions on the perimeter of the finite-size structure shown in Fig. 3(a) will give rise to a finite loss. The quality factors of the simulated graph modes are on the or-der of 10 5 to 10 7 , allowing for clear identification of all eigenfrequencies. The graph topology is kept fixed in the eigenfunction studies discussed below. Because the topology and the total length L tot of each graph realization are different, we normalized the system eigenmode 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 (k i , i is the mode index) using the dispersion relation for the PC-waveguide PBS shown in Fig. 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 2 k 2 + c 1 k + c 0 , where c 0,1,2 are fitting parameters and n(k) is the total number of modes below k [46]. For a typical graph we found c 2 = 0.0054m 2 , c 1 = 2.21m and c 0 = −55.28. The small value of c 2 indicates the n(k) is nearly linear in k as expected for quasi-one-dimensional structures. We then evaluate the fitting functions at the calculated k-values for each graph, creating a list of normalized eigenvalues e i = n f it (k i ). The nearest neighbor spacings are finally computed as s i = e i+1 − e i .
The distribution of the normalized nearest-neighbor mode spacing values of the entire ensemble is shown as the histogram in Fig. 3(c). We have included the approximate theoretical predictions of the mode-spacing statistics for the GOE and GUE systems in the figure. These theoretical predictions are based on RMT [47][48][49], which are methods of understanding the universal statistical properties of wave-chaotic systems [22,50]. As shown in Fig. 3(c), the distribution of the normalized modespacing 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 [4,[30][31][32].
To test the long-range statistical properties of the eigenmodes, we compute the spectral rigidity of the graph spectrum ∆ 3 (L), where L measures the length of a segment in the normalized spectrum. Our procedure is as follows. The quantity ∆ 3 (e, L) representing the deviation of the spectral staircase from a straight line, is evaluated for each value of L by where N (e) = n(k) and e is chosen at random. We use 100 values of e for each graph. Then we average the values of e over 10 graph realizations [51]. The resulting values of ∆ 3 (e, L) are then plotted in Fig. 3(d) along with theoretical predictions based on GOE, and GUE random a. Schematic diagram of the PC waveguide graph system. The graph bonds are shown as 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 and GUE systems are shown in red and yellow curves. d shows the spectral rigidity ∆3(L) of the normalized spectrum from graph simulation (blue circles). The theoretical predictions for the GOE, GUE, and the Poisson systems are shown in red, yellow, and black curves. e shows the statistics of the consecutive mode spacing ratio r (histogram) and the theoretical predictions for the GOE and GUE systems. matrices [32], and the predictions of Poisson statistics [52] (appropriate for classically integrable systems).
The ∆ 3 (L) curve in Fig. 3(d) is linear in L for small values of L. As shown, systems of all symmetries (Poisson, GOE, GUE) show a linear-in-L dependence for L < 1 [52], and the data is consistent with this trend. Chaotic systems are then expected to show a logarithmic increase in ∆ 3 (L) [52], and this is evident in the data up to approximately L = 5. It is also predicted that ∆ 3 (L) will saturate due to the presence of short orbits in any experimental realization of a wave chaotic system [52], and this is also seen in the data. Previous studies find that the long-range statistics of quantum graphs usually show some degree of agreement with the GOE prediction in the range of roughly L ∈ [2,5], along with deviations beyond this range [1,30,31]. We observe in Fig. 3(d) a deviation from the GOE RMT predictions, especially for larger L, similar to cable graph studies in Ref. [1]. Because the operating bandwidth of the proposed PC-graph system is relatively narrow (Fig. 1 (c)), we believe the longrange statistics may also be affected by finite-size effects because the number of available eigenmodes is bounded above and below by the edge of the PC bandgap. The solution is to increase the graph size substantially, but this is beyond our computational resources. In summary, the ∆ 3 (L) behavior is very similar to that observed in wave chaotic systems of higher dimensionality. In addition to the mode-spacing distribution test, we note that the method of consecutive mode spacing ratios can also be adopted [53]. Here the ratio statistic has the feature that it is a measure of the correlation between adjacent spacings. One advantage of this type of statistical study is that it does not require unfolding the spectrum [53]. However, here we use the unfolded spectrum s i to compute the spacing ratios. 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 prediction (Table. I). We note that the high value of r in the PC-graph is related to the correlation between consecutive level-spacings. We further present the distribution of mode-spacing ratios r of the PC graphs and corresponding theoretical predictions of the GOE and GUE systems in Fig. 3(e). The statistics of r match reasonably well with the GOE prediction.

IV. 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 [6,[54][55][56][57][58]. In those 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 Schrö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 E z value at the middle cutting plane (at the height of z = h 0 /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. 1(d)). We next examine two PC graph eigenfunction characterization methods.
Method I: grid-wise representation. For each eigen-mode of each graph realization, we will use the entire set of the graph bond |E z | 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 where L tot is the total length of the graph, r j and ∆L j are the location and the grid size of the j th grid point, |E z (r j )| 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 ∆L j ∼ 0.05λ op where λ op = 10cm is the operating wavelength at 3GHz. The distribution of the probability density values is computed using the data from all simulated eigenmode solutions (48 to be specific) from a single realization, and the results are shown in Fig. 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' E z (b m ) which is defined as the amplitude of the standing-wave wavefunction on the graph bond b m . 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 m e ikx + a * m e −ikx , where a m is the wavefunction amplitude at bond b m and x is the distance from a vertex along the bond. We first conduct a sine-fit of the raw E z (x) values on each bond, which yields the amplitude, and the value of |a m | is obtained as 1/2 of the amplitude value. The normalization process follows the same method as in Ref. [59] which ensures that 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 of eigenmode solutions from one graph configuration, and the results for P (|a|) are shown in Fig. 4(b).
The random plane wave hypothesis underlying chaotic eigenfunction treatments asserts that the eigenfunctions in a 2D or 3D cavity can be thought of as superpositions of plane waves with random directions and random phases. Consequently, the value of an eigenfunction at any point is a sum of many random waves and is thus a zero mean Gaussian random variable [47]. For systems with TRS the plane waves come in counter-propagating pairs, and the wave function is real. For systems with TRS broken, wave directions are completely random, and the wave function is complex with independent real and imaginary parts. For the PC graph the situation is changed in that at any point in the graph there are only two counter-propagating waves of equal amplitude. The amplitudes of these waves, |a m | will be different on the different bonds leading to statistical variations.
We first discuss the eigenfunction statistics obtained using the grid-wise method in Fig. 4(a). For waves in a 2D or 3D cavity the distribution of values of v based on the random plane wave hypothesis is P (v) = (2πv) −1/2 e −v/2 . This is shown as a solid red line in Fig. 4(a). We find that the simulated PC graph result matches to some degree the GOE prediction in that for v > 1 the curve approaches a straight line on a loglinear plot, and for v small there is an indication of square root singularity. We notice that the mismatch between the graph data and the GOE prediction exists at the small and large probability density values, indicating an absence of both very small and very large v 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. [6]). 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) [8,60].
Here we now use Method II to study the distribution of normalized bond values, |a| that are displayed in Fig.  4(b). Also plotted in Fig. 4(b) is a Gaussian distribution, which would be characteristic of an eigenfunction of a GOE random matrix. However, a clear deviation from Gaussian statistics is seen for low amplitudes, similar to the deviation observed with Method I. The inset to Fig.  4(b) shows a plot of P (|a| 2 ) compared to the Porter-Thomas distribution expected for GOE systems. There is an agreement between the two except for deviations between data and the model distribution at small values of |a| 2 , similar to that seen in Fig. 4(a). 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. 5. IPR is a measure of the degree of localization for a wave function [59,61], 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. [59], where the IPR value for each graph mode is evaluated using the formula IP R IIa = |a m | 4 . As shown in Fig. 5(a), the IPR values of the photonic crystal graph modes vary erratically but lie between the maximum ergodic limit (IP R IIa = 1) and the RMT prediction limit (IP R IIa = 2) [59]. In the maximum localization limit the IP R IIa = B, where B = 14 is the number of graph bonds, and the results are far from this limit. Method IIb follows the definition in Ref. [61] where . The quantityã m is the un-normalized bond wavefunction amplitude. Here, the graph IPR values also vary erratically from mode to mode ( Fig. 5(b)), but lie well below the maximum localization limit (IP R IIb = 1) and closer to the RMT prediction limit (IP R IIb = 1/B = 0.07) [61]. 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. 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 convenient feature of the photonic crystal graph system. For Method IIb, we present the eigenfunction profile of two neighboring graph modes in Figs. 5(e) and (f). The eigenmode in Fig. 5(e) has a larger value of IP R IIb and shows a strongly localized distribution. That in Fig.  5(f) has a small value of IP R IIb and is more evenly distributed over the bonds. The average value of IPR over all the modes is more meaningful in this case, and we have IP R IIa = 1.39 and IP R IIb = 0.10 which are closer to the ergodic and RMT limits, respectively, than to the localization limit. Previous work [1,59] 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. We note that by computing the IP R IIa on single loops rather than the entire graph, one may identify the modes trapped inside these loops from IP R IIa ∼ B sl , where B sl is the total bond number of the single loop. This study is discussed in the Appendix.
Taking into account the eigenfunction and IPR statistics, it would appear that PC graphs are close to the random plane wave condition for tetrahedral-like graphs, but clear systematic differences from RMT predictions remain. The complex vertex scattering properties of rightangle bends and T-junctions, along with the ω(k) dispersion relation of PC defect waveguide modes, suggest that PC defect graphs may be very effective for future wave chaotic studies.

V. CONCLUSION
To conclude, we have designed and simulated an alternative chaotic graph system with photonic crystal defect waveguides that show a dispersion relation that is different from that utilized in coaxial cable microwave graphs. We show that a series of statistical studies can be carried out on an ensemble of closed graphs, including nearestneighbor spacing statistics and eigenfunction statistics studies. Because both the graph bonds and nodes can be probed, one may better analyze the non-universal features of chaotic graphs using a PC system. It is possible to adjust the degree of wave localization by engineering the scattering properties of the PC waveguide junctions. This property of PC waveguides may facilitate further studies of localization phenomena in graphs, for example, the emergence or suppression of trapped modes.