Implications of Spectral Interlacing for Quantum Graphs

Quantum graphs are ideally suited to studying the spectral statistics of chaotic systems. Depending on the boundary conditions at the vertices, there are Neumann and Dirichlet graphs. The latter ones correspond to totally disassembled graphs with a spectrum being the superposition of the spectra of the individual bonds. According to the interlacing theorem, Neumann and Dirichlet eigenvalues on average alternate as a function of the wave number, with the consequence that the Neumann spectral statistics deviate from random matrix predictions. There is, e.g., a strict upper bound for the spacing of neighboring Neumann eigenvalues given by the number of bonds (in units of the mean level spacing). Here, we present analytic expressions for level spacing distribution and number variance for ensemble averaged spectra of Dirichlet graphs in dependence of the bond number, and compare them with numerical results. For a number of small Neumann graphs, numerical results for the same quantities are shown, and their deviations from random matrix predictions are discussed.


Motivation
Quantum graphs are composed of bonds which are connected with each other at vertices. Along the bonds wave propagation is governed by the Schrödinger equation without potential and boundary conditions depending on the details of the coupling at the vertices. Quantum graphs were first introduced by Pauling [1] in the context of free electron models of organic molecules. Later, they were studied intensely in physics [2] and mathematics [3], and experimentally implemented in correspondingly-shaped microwave networks [4]. They are conceptually simple, but still complex, and there is a straightforward symbolic alphabet to classify the periodic orbits. Casati and coworkers [5] suggested that the universal features of the spectra of chaotic systems might be described by random matrix theory (RMT), which later was expressed by Bohigas, Giannoni, and Schmit [6] in the form of a conjecture. Using supersymmetry techniques, Gnutzmann and Altland [7] proved the conjecture for the two-point correlation function for fully connected graphs with incommensurate bond lengths. Their result was generalized to all correlation functions by Pluhař and Weidenmüller [8]. Just as for billiard systems [9], there is a one-to-one correspondence between a quantum graph and the corresponding microwave networks, which has been used, in particular, by Sirko and coworkers in numerous experiments to study spectral and scattering properties of microwave graphs (see e.g., Ref. [4]).
In a recent microwave experiment in tetrahedral graphs [10], however, we noticed that one important aspect is missing in the above scenario. It is hidden in the structure of the equation system determining the graph spectrum. Using energy and current conservation (the Kirchhoff rules in experimental networks), one arrives at a secular equation [2] V ∑ m=1 h nm ϕ m = 0 , where the sum runs over all vertices V, and ϕ m is the potential at vertex m. In the experiment, the bonds are connected by ordinary T junctions corresponding to Neumann boundary conditions. For this situation, the elements of the secular matrix h are given by where the l nm are the lengths of the bonds connecting vertices n and m, and k is the wave number. For the homogeneous equation system (1) to have non-trivial solutions, the determinant of h(k) has to vanish, The roots k n of the equation generate the spectrum of the graph. It will be called "Neumann" spectrum in the following since the T junctions at all vertices obey Neumann boundary conditions. On the other hand, h nm becomes singular, whenever kl nm is an integer multiple of π. This situation corresponds to a totally disassembled graph with a spectrum being the sum of the spectra of all individual bonds with Dirichlet boundary conditions at both ends, thus the vertices have no influence any longer. This "Dirichlet" spectrum hence appears via the poles of |h(k)|, whereas the Neumann spectrum is given by the zeros of |h(k)|. In the following, all lengths will be assumed to be incommensurable to avoid degeneracies of the Dirichlet spectrum. This "spectral duality", as we termed it in our previous publication [10], has important consequences for the spectral statistics. The cause is the interlacing theorem (see, e.g., Chapter 3.11 of Ref. [3]): If the boundary conditions at one vertex of a graph are changed from Neumann to Dirichlet, or somewhere in between, the eigenvalues of the original and the new graph appear strictly alternating.
To move from the Neumann to the Dirichlet spectrum for a complete graph, the boundary conditions have to be changed one after the other at all vertices, not just at one of them. Now, there is no longer a strict alternation in the sequence of the respective eigenvalues, but still a strong correlation remains-the maximum number of Neumann eigenvalues confined between two successive Dirichlet ones is given by the number of vertices V, and vice versa.
The mean density of states for a graph of a total length of l tot is given bȳ For a graph with B bonds and a given total length, the maximum level spacing is found for the limiting case where all bonds are equal. For this case, the Dirichlet spectrum is B-fold degenerate, and the maximum distance between neighboring eigenvalues, in units of the mean level spacing, is just s max = B. Due to the interlacing theorem, the same must be true for the Neumann resonances. There is hence a cut-off in the level spacing distribution p(s) at s max = B, at the latest, both for the Dirichlet and the Neumann spectrum. Consequences of spectral interlacing for the number variance have been discussed already in our previous paper [10]. Thus, there are clear deviations from the RMT expectation for small graphs. This is not in contradiction with the proofs mentioned in the beginning that the spectra of graphs with irrational length ratio do exhibit RMT behavior, since these proofs work in the limit of infinitely large graphs only. From the practical point of view this is of little help since numerical, as well as experimental, studies are necessarily restricted to comparatively small graphs.
Therefore, an understanding of the impact of Dirichlet-Neumann interlacing is mandatory for the correct interpretation of the spectral statistics in small graphs. Since, in the moment, a good idea to approach Neumann spectral statistics is still missing, we start with a more modest task-the interpretation of Dirichlet spectral statistics. Analytic results are given for level spacing distribution and number variance for a random superposition of lattice fence spectra and compared with numerical results. For the Neumann spectra, we restrict ourselves to an illustration of the fingerprints of spectral interlacing in level spacing distribution and number variance, but have to leave the theoretical interpretation to future papers. We do not discuss experimental results from microwave graphs in the present paper. This remark may be necessary since probably this is exactly what readers do expect from our group.

Dirichlet Graphs
For Dirichlet graphs, there are Dirichlet boundary conditions at each end of all bonds, thus the bonds are not coupled at the vertices and the spectrum corresponds to a superposition of B separated bonds. Here, we present analytical and numerical results of the spectral statistics for ensemble-averaged Dirichlet graphs. Following the usual practice, the mean density of statesρ(k) = l tot /π was kept constant and normalized to one, meaning a total length of l tot = π for all graphs entering the average. For the numerics, the lengths had been created by generating B − 1 random numbers r n between 0 and π, and by taking the appearing B segments as lengths l n . The procedure yields p B (l 1 , . . . , l B ) = 1 π B−1 δ(∑ l n − π) for the joint length probability. The l n are hence uniformly distributed on the interval 0 to π with the constraint ∑ l n = π. Integrating out all l n but one obtains the distribution for the remaining l, being constant only for B = 2. The derivation and a plot of the length distributions can be found in Appendix A. Alternatively, one could think of taking B lengths l n from an interval between 0 and 1, and afterward, normalizing each length via l n = πl n / ∑ B n=1 l n to a mean density of one, i.e., ∑ l n = π. The resulting joint length probability is non-uniform, in contrast to the one above. For the sake of conciseness, we shall refer in the following to the two respective ensembles as the uniform and the non-uniform one. The non-uniform approach would be more in the spirit of the usual unfolding technique used in quantum chaos to make spectra taken from different systems comparable. For the non-uniform ensemble again, numerical length distributions are presented in Appendix A. Since it would be hard to obtain analytical results for the non-uniform ensemble, all analytics and numerics, if not explicitly stated differently, are for the uniform one.
In the next two subsections, theoretical expressions for nearest neighbor spacing distribution p(s) and number variance Σ 2 are given and compared with numerical data.

Nearest Neighbor Spacing Distribution for Dirichlet Graphs
To calculate the distribution of nearest neighbors spacings p(s) for the Dirichlet spectrum of a graph, we apply a strategy that had been used already by Berry and Robnik [11] to calculate p(s) for an uncorrelated superposition of two spectra, one associated with the chaotic part, the other with the regular part of a mixed phase-space system. A key element in the calculation is the gap probability e(s) describing the probability for a spectral range of length s to be empty of eigenvalues. The gap probability is related to the level spacing distribution via p(s) = e (s) , where a mean level spacing of one has been assumed. Expression (6) is well-known to those working in the field, but for readers not familiar with the subject, a didactic derivation is given in Appendix B. For a picket fence spectrum with a mean level spacing of ∆s = 1, the gap probability is given by e(s) is, in contrast to p(s), multiplicative for superimposed uncorrelated spectra, whence follows for the Dirichlet spectrum of a graph with B bonds of lengths l n , n = 1, . . . , B From Equation (9) for e B (k), the Dirichlet level spacing distribution can now be obtained by taking the second derivative, see Equation (6). The Dirichlet level spacing distribution has already been calculated by Barra and Gaspard [12], who did not follow, however, the approach of Berry and Robnik [11]. Their derivation therefore is much less concise and considerably longer then the present one, see the appendix of [12]. Expression (9) has to be averaged over all different realizations of l n with the constraint that the total length l tot is constant with the substitution s n = l n π k, the constraint becomes Thus, s is the wave number in units of the mean level spacing π/l tot . In the following we shall use the letter s exclusively for spectra with a mean level spacing of ∆s = 1. Now the average can written as where w B (s) is given by The factorial in Equation (12) reflects the number of possible l sequences to do the average. Equation (13) can be used to calculate w B (s) iteratively. For B = 2, e.g., one obtains where the limits of integration in the different s windows take care of the cut-off of e(s).
With help of Equations (12) and (6), we now obtain for the level spacing distribution In this way, the p B (s) may be obtained iteratively, resulting in formulas with a complexity increasing step by step.
A more direct approach takes advantage of the fact that the integral in Equation (13) is nothing but a convolution. In such a situation, Laplace transform techniques are the method of choice. Applying a Laplace transform to Equation (13), the convolution theorem yieldŝ whereŵ are the Laplace transforms of w(s), and e(s), respectively. Iterating Equation (16), one getŝ whence w B (s) is obtained via an inverse Laplace transform The inverse Laplace transform can be done with the result where s denotes the largest integer ≤ s and Further details can be found in Appendix C.
To verify our results, we compare the analytical results with numerical simulations. In Figure 1, the histograms for B = 2 to 6 and 100 are shown together with the corresponding theoretical predictions in a linear, and in Figure 2 in a logarithmic scale. All distributions show the expected cut-off at s max = B. There is a perfect agreement between numerics and theory. Note the discontinuity for B = 2 at s = 1, which for larger B is smoothed out and vanishes for B → ∞, where an an exponential decay is expected, corresponding to a Poisson distribution. This can be seen in Figures 1f and 2f, showing the results for B = 100. There are still deviations from the exponential behavior as can be seen in the inset of Figure 2f (6) and (12). In (f), the orange dashed line corresponds to a Poisson distribution, i.e., an exponential.   Figure 1 but in logarithm scale. In the inset in (f), the abscissa ranges from 0 to 16, and the ordinate from 10 −9 to 5.

Number Variance for Dirichlet Graphs
The number variance, defined as where n is the number of eigenvalues in an interval of length s, yields for the lattice fence spectrum of a single bond of length l It is convenient to express Σ 2 (s) in terms of its Fourier transform, For B bonds with independent bond lengths l n , the spectrum is just the superposition of the B spectra with bond lengths l 1 , l 2 , . . . , l B . The number variation Σ 2 (s) is additive for uncorrelated spectra leading to for the ensemble averaged number variance, where s k = kl k π , and . . . l means the average over all l n with the constraint ∑ l k = l tot , i.e., with α = 2mk and p B (l) given by Equation (5). In Figure 3, the ensemble averaged number variance for Dirichlet graphs are shown for a number of different bonds. For a single bond (B = 1), there is just one lattice fence spectrum with a spacing of one. Hence, one observes a periodic modulation with an average value of 1/6, as described by Equation (25). With increasing B, these oscillations are damped out more and more, until Σ 2 (s) eventually approaches the linear increase expected for a Poissonian ensemble. A good agreement between the simulations and the analytical predictions given by Equation (27) is found.

Neumann Graphs
Here, we present the results of the Neumann graphs shown in Table 1. We restricted ourselves to graphs where bonds are connected at least to two other bonds and where there are no disconnected parts. In addition the verticity V B , the number of bonds that connect at a vertex, has been assumed to be the same for all vertices. The smallest graph of interest is the tetrahedron which has been used repeatedly for RMT studies in numerics [13], as well as in experiments [4,10]. Table 1. The investigated graphs for Neumann boundary conditions at the vertices. The lower part of the table shows the results of the numerics: (i) s c is the s value, where p(s) drops below 10 −5 , i. e., p(s c ) = 10 −5 , (ii) s m is the s value, where Σ 2 (s) takes its maximal value, and (iii) Σ 2 sat is the limit of Σ 2 (s) for s → ∞, obtained by taking the average of Σ 2 (s) in the range of s between 10 and 20.  Figure 4 shows level spacing distributions for the graphs presented in Table 1. In addition, the exact RMT nearest neighbor distribution is shown [14,15]. Whereas the linear plot suggests a reasonable agreement with the RMT prediction, in the logarithmic plot for all four graphs, a suppression for large values becomes obvious, with the strongest suppression for the tetrahedron having the smallest number of bonds. This is in accordance with the expectation; due to the interlacing theorem, the largest possible distance is given by the number of bonds, B = 6 for the tetrahedron, and 10 or 12 for the other graphs. The decay, however, does not increase monotonously with B-it is faster for the hexahedron than for the octahedron, though the number of bonds is the same, and the decay for the fully connected five-vertices graph is as fast as for the octahedron, though the number of bonds is not the same.

Nearest Neighbor Spacing Distribution for Neumann Graphs
Similar deviations of p(s) from RMT have also been observed by Barra and Gaspard [12] but not discussed in detail.
To quantify these findings, we determined s c , the value where p(s) drops below 10 −5 , i.e., p(s c ) = 10 −5 , close to the limit of our statistical precision. In Figure 4, this s value reflects the point where p(s) crosses the abscissa. The extracted values are presented in Table 1. Regrettably, it is impossible to fix the cut-off point by the numerics, which should be, e.g., s max = 6 for the tetrahedron. This would need more than 10 11 spacings, by far beyond our computer resources.  Table 1 in a logarithmic scale, with the tetrahedron (blue), the fully connected graph with 5 vertices (orange), the octahedron (green), and the hexahedron (red). The plots were generated by superimposing the results from 52.2 · 10 6 , 14.6 · 10 6 , 6.7 · 10 6 , 35.0 · 10 6 spacings. The inset shows the same data in a linear scale. In addition, the exact RMT distribution is plotted (dashed black).

Number Variance for Neumann Graphs
In Figure 5, the ensemble-averaged number variances for the graphs shown in Table 1 are plotted, exhibiting a saturation at about s = 2 in contrast to the behavior predicted by RMT. Similar deviations from RMT predictions for the number variance in graphs have been reported in Refs. [16,17], and for the spectral rigidity in Ref. [4]. For non-experts, we mention that the spectral rigidity may be looked upon as a smoothed version of the number variance (the exact definition is technical and not of relevance here [14]). Already in 1985, Casati and coworkers [18] discovered a saturation of the spectral rigidity in the spectra of rectangular billiards, which could be traced back by Berry to the influence of the shortest periodic orbit [19]. In the present case this can not be the explanation. From the shortest periodic orbit, there should be a saturation of Σ 2 (s) at about s sat = π/l min = l tot /l min . Since necessarily l min ≤ l tot /B, periodic orbit theory predicts a saturation of the number variance not until s sat ≥ B, whereas actually for all graphs the saturation is observed much earlier at about s = 2. In the lower part of Table 1, the saturation values are given, as well as s m corresponding to the s value where Σ 2 (s) is maximal. s m is a convenient tool to quantify the point of cross-over from a linear increase of Σ 2 (s) for small s values to a saturation for s → ∞. There is a clear correlation between s c , s m , and Σ 2 sat , collected in the lower part of Table 1. From the interlacing theorem one would expect a correlation of these quantities with the bond number, which is observed for the first three graphs presented in the table, the tetrahedron, the totally connected graphs with five vertices, and the octahedron, but the hexahedron does not fit into the sequence. In fact, there is a stronger correlation with the vertex valency, the number of bonds meeting at a vertex. Obviously, the interlacing theorem alone is not sufficient to describe all these features. A quantitative explanation, in particular, of the saturation values, has to be postponed to further studies, but there is a qualitative explanation. Semiclassical theory relates RMT to periodic orbits [20]. Essential ingredients are correlations between orbits and their timereversed partners, in case there is time-reversal symmetry [19], and between various types of self-intersecting orbits with their non-intersecting partners [21,22]. Apart from this, all orbit lengths are assumed as uncorrelated. This assumption is severely violated in graphs, where all orbits are composed from a finite number of elements.

Conclusions
The implications of spectral interlacing in quantum graphs have been discussed. For Dirichlet graphs, explicit analytic expressions have been obtained for level spacing distribution and number variance, and compared with numerical results. For Neumann graphs, numerical results for the same quantities have been presented, showing clear deviations from RMT predictions due to spectral interlacing. For Neumann graphs, an analytic description of these features is still missing and has to be left to future work.

Appendix A. Length Distributions
In most of our numerical studies, the lengths were created by generating B − 1 random numbers r n between 0 and π, and taking the appearing B segments as lengths l i , resulting in the uniform ensemble introduced in Section 2. This results in a joint probability for the lengths l 1 , . . . , l B . The brackets denote an average over the r i , with the constraint 0 ≤ r 1 ≤ · · · ≤ r B−1 ≤ π. Due to the delta functions, the integrations are done easily step-by-step with the result The l n are thus equally distributed over the interval 0 to π with the constraint ∑ l n = π. Integrating over all l n but one, the distribution for a single length l is obtained, where the factorial takes account of the number of possible sequences to do the integrations. The integrations can be performed iteratively, the result is Equation (5). Figure A1 shows numerically obtained length distributions together with the theoretical curves from Equation (5). A perfect agreement is found. Figure A2 shows the same for the non-uniform ensemble but without the analytic. There are clear differences to the uniform ensemble. In particular, for B = 2, there is a uniform length distribution for the first ensemble, whereas for the second one, there is a cusp in the distribution for l = π/2. We did not try to calculate analytic expressions for the length distribution for the latter case, they are not needed in the present context.   Figure A2. The same as in Figure A1, but for the non-uniform ensemble, and without analytical expression.

Appendix B. Derivation of Equation (6)
The relation between gap probability e(k) and level spacing distribution p(k) is obtained as follows (see, e.g., Section 3.2.2 of Ref. [23]): The difference (e(k) − e(k + ∆k)) may be interpreted as Here, h(k) is the one-sided gap probability, namely the probability for an interval of length s to be empty of eigenvalue, and an eigenvalue at one end of the interval, let us assume the lower one. ρ∆k is the probability to have an eigenvalue in the interval ∆s, where ρ is the mean density of states. It follows Furthermore, [h(k) − h(k + ∆k)] may be interpreted as the probability to have one eigenvalue at the lower end and the next nearest one between k and k + ∆k, i.e., whence follows where p(k) is the probability density for an interval to be empty of eigenvalues, and an eigenvalue on each side, i.e., the nearest neighbor spacing distribution. Combining Equations (A6) and (A8), one obtains the wanted relation between e(k) and p(k), reducing to Equation (6) for a spectrum with a mean density ρ = 1. For a graphical illustration see Figure A3.

Appendix C. Calculation of the Inverse Laplace Transform (20)
Using the standard formula for the inverse Laplace transform Equation (20) yields The path of integration has to be performed to the right of all singularities of the integrand. The only possible singularity is at λ = 0, hence a has to be positive real. The path of integration can be closed at infinity by a large semi-circle either to the left, for x > k, or to the right, for x < k.
For x < k, the integration loop encloses no singularities, i.e., dλ · · · = 0, for x > k the integration can be performed by means of the residuum method resulting in where k max = Min(B, x ). For x > B, Equation (A11) yields since (e −λ − 1 + λ) = O(λ 2 ). For x < B, the sum in Equation (A11) runs from 0 to x , and after a number of elementary steps one arrives at Equation (21).

Appendix D. Dirichlet Level Spacing Distributions for the Non-Uniform Ensemble
To illustrate the difference between the two ensembles introduced in Section 2, we present here numerical results for the level spacing distribution for the non-uniform ensembles. Figures A4 and A5 show the results in linear and logarithmic scale. In addition, the analytical curves for the uniform ensemble are plotted. Not surprisingly, there are deviations in detail, but the qualitative behavior, in particular, the discontinuities for integer values of s, and the cut-off at s = B is the same in both cases. Note that the non-uniform ensemble approaches the Poisson limit faster than the uniform one, see inset of Figure Figure A5. The same as Figure A4, but in a logarithmic scale. In the inset in (f), the abscissa ranges from 0 to 16, and the ordinate from 10 −9 to 5.