Concrete method for recovering the Euler characteristic of quantum graphs

Trace formulas play a central role in the study of spectral geometry and in particular of quantum graphs. The basis of our work is the result by Kurasov which links the Euler characteristic $\chi$ of metric graphs to the spectrum of their standard Laplacian. These ideas were shown to be applicable even in an experimental context where only a finite number of eigenvalues from a physical realization of quantum graph can be measured. In the present work we analyse sufficient hypotheses which guarantee the successful recovery of $\chi$. We also study how to improve the efficiency of the method and in particular how to minimise the number of eigenvalues required. Finally, we compare our findings with numerical examples$-$surprisingly, just a few dozens of eigenvalues can be enough.


Physical model
Schrödinger operators on metric graphs have long been studied as simplified models of complex quantum mechanical systems, notably large molecules [1,2]. More recently, under the name of quantum graphs, they have been widely used to describe nano-wires, wave-guides and networks [3,Chapter 8]. They have also been investigated by mathematicians as a simplified setting in which to understand complex problems in spectral theory, such as quantum chaos [4], trace formulas [5,4,6,7,8] and nodal patterns [9]. From this point of view, they stand midway between one-dimensional and multi-dimensional Schrödinger operators. Indeed, metric graphs are one-dimensional objects with a non-trivial geometry. The central theme in studying quantum graphs is therefore the interaction of the spectrum with the geometry of the underlying metric graph.

Mathematical problem
We consider in this work a finite, compact and connected metric graph Γ. Following [10], we describe it by a set of N edges {e 1 , . . . , e N }, with e n = [x 2n−1 , x 2n ], and a set of M vertices {v 1 , . . . , v M }, which is a partition of the set of endpoints {x p : 1 ≤ p ≤ 2N }. The associated quantum graph is an operator in the Hilbert space More precisely, we define the standard Laplacian L Γ on Γ, which is the object of our study, in the following way. Its domain Dom(L Γ ) consists of functions f ∈ L 2 (Γ) satisfying 1. for every 1 ≤ n ≤ N , f ∈ W 2 2 (e n ); 2. for every 1 ≤ m ≤ M , (a) f (x p ) = f (x q ) whenever the endpoints x p and x q belong to the vertex v m (continuity condition), (b) xp∈vm ∂f (x p ) = 0 (balance condition).
In the balance condition, ∂f (x p ) is the derivative along the edge with endpoint x p , taken towards the inside of the edge. The conditions in (ii) are called the standard vertex conditions. The action of L Γ on f ∈ Dom(L Γ ) is then given by the differential expression L Γ f (x) = −f (x) for all n ∈ {1, . . . , N } and x in the interior of e n . It is known that the operator L Γ is self-adjoint, non-negative and has compact resolvent. Its spectrum therefore consists of a sequence of non-negative eigenvalues (λ j ) j≥1 of finite multiplicity tending to +∞ (we repeat them according to their multiplicities). Let us note that λ 1 is always 0 and that the hypothesis that Γ be connected implies that it is simple, with all the associated eigenfunctions constant. For the sake of convenience, we will work in the rest of the paper with the eigenfrequencies k j = λ j .
Our goal is to recover information about the geometry of the graph from the knowledge of some of its eigenfrequencies. This is a variation on Kac's famous problem "Can one hear the shape of a drum?" [11], with some specificities that we will address below. The analogue to Kac's problem in our setting would be showing that two graphs with the same spectrum are isometric. There are partial results in that direction. For instance, the eigenfrequencies of a quantum graph satisfy Weyl's law: as j → +∞, k j ∼ πj/L, where L is the total length of the graph, i.e. the sum of the lengths of all its edges (see for instance [10]). The total length of the graph can therefore be read off the complete spectrum. Furthermore, the answer to Kac's problem is positive provided all edge-lengths of the graph are rationally independent [6,7]. Nevertheless, if this last condition is not satisfied, we can find graphs with the same spectrum that are not isometric [6,12,13,14].
We will restrict ourselves to recovering the main topological invariant of a connected graph: its Euler characteristic χ, defined by χ = M − N . It is related to β 1 , the number of independent cycles in the graph, by β 1 = 1 − χ.

Objectives
We focus on a situation that can be realized experimentally. A microwave network can be described to a good approximation by a quantum graph, whose eigenfrequencies correspond to resonances of the system. In [15] the authors describe the measurement of about 50 resonances of such networks and use this experimental result to recover χ. Contrary to expectations, they find this relatively small number of eigenfrequencies to be sufficient. Their computation is based on a trace formula as we will see below in more details. Extending the work in this paper, we want to check the correctness of the method and modify it in order to improve the efficiency, or namely to reduce the number of eigenfrequencies required.
The total length of the graph is a continuous parameter and an infinite number of eigenfrequencies are needed to compute it using Weyl's law, as described above. By contrast, χ is integer valued and we can hope to compute it with only a finite number of them. However, some assumptions on the graph are needed. Indeed, by attaching a very small loop to some edge in an existing graph, we decrease χ by 1 while altering the first 50 eigenfrequencies by an arbitrarily small amount. We come back to this point at the end of Section 5. Therefore, our goal will also be to determine a reasonable set of assumptions that guarantees the success of our methods.
2 Trace formula 2.1 Main tool All our methods for computing the Euler characteristic are based on the trace formula proved in [7,8], building on previous work in [5,4,6]. This formula connects the spectrum of L Γ , defined in the previous section, with the set of periodic orbits in Γ, to be described below. More precisely, the spectrum is described by the distribution which is a tempered distribution of order 0. We denote by P the set of closed continuous paths in Γ that are allowed to turn back on themselves only at vertices, which we call periodic orbits, or orbits. If p ∈ P, we denote by prim(p) any shortest orbit such that p can be obtained by repeating prim(p). Thus any two such orbits have the same length, hence (prim(p)) is well defined. The following quantity will play a crucial role in our methods: Proposition 2.1 (Theorem 1 from [7,8], also [10]). We have where S V (p) is the product of all the vertex scattering coefficients along the path p.
Let us point out that, in all the paper, we define the Fourier transform so thatf (k) = +∞ −∞ e ik f ( ) d for any f ∈ L 1 (R). Remark 2.2. Taking the Fourier transform on both sides of Equation (3), we obtain which is the form of the Trace Formula that we will be using in the rest of the paper.
By definition of the Fourier transform of a tempered distribution, for ϕ in the Schwartz class S(R),û[ϕ] = u[φ], which, using Equations (1) and (4), means Equation (5) holds for a class of functions larger than S(R). For our purpose, we will consider the set T of function ϕ satisfying the following assumptions: 1. ϕ is continuous on R and nonnegative; Let us make some comments on Condition (iv). The function F ϕ can be seen as the smallest non-increasing majorant of | ( ϕ)| on (0, +∞). Its integrability is used both to express χ as the sum of a convergent series and to prove error estimates.
We prove the above proposition in A. Let us note that whenever ϕ is in T , so is ϕ t (l) := tϕ(t ), for any t > 0.
We can now state the result by which we recover the Euler characteristic of Γ.
where min is defined in (2).
In the rest of the paper, we assume that t ≥ 1/ min . When we will analyse our method, it will accordingly be important to know a lower bound of min in order to guarantee a successful estimate of χ.
Remark 2.5. The above propositions hold for the standard Laplacian, specifically for standard vertex conditions imposed at all the vertices. In [10] is presented a general version of the trace formula for arbitrary vertex conditions where in particular the constant term inû may differ from χ. Hence our method might not be applicable in presence of vertex conditions other than standard. It is even possible that χ cannot be recovered knowing the whole spectrum of the graph: see [12,13,14] for an example of two graphs with mixed vertex conditions which have distinct values of χ despite being isospectral.

Choice of the test function
Let us apply Formula (7) to a triangular function: Everywhere in the paper, 1 A stands for the characteristic function of the set A. We find Proposition (2.4) then gives us the following expression for χ, assuming t ≥ 1/ min : We can also apply the same formula to a truncated trigonometric function: We deduce another expression for χ, valid whenever t ≥ 1/ min : This formula is the one employed in [15]. In order to recover χ from a finite set of eigenfrequencies, we have to replace the series on the right-hand side of Formulas (9) and (10) by partial sums. Since the terms in the second series go to zero more quickly than in the first, we can expect that fewer eigenfrequencies are necessary to find a close approximation of χ, which we will then be rounded (up or down) to the nearest integer. Both numerical trials and the analysis of the experimental results in [15] confirm this. The more regular ϕ is, the faster ϕ decays at infinity. We can then wonder how much we can decrease the number of eigenfrequencies used to compute χ. To study this question, we will consider the one-parameter family of test functions where d is a positive integer and C d is chosen so that +∞ −∞ ϕ d ( ) d = 1. We easily check that the function ϕ d is of class C 2d−1 on R. Our goal is to give rigorous error estimates and to perform numerical trials for this family of functions.
3 General estimate

Calculation with exact spectrum
For a fixed graph Γ, we assume that only the smallest J eigenvalues of the standard Laplacian are known. We consider the distribution (1) evaluated on the Fourier transform of a generic function ϕ ∈ T and split the expression in the following way We call S J (t) the truncated sum and R J (t) the remainder or tail. We have Let ϕ ∈ T such that T = 1 in Condition (ii) be fixed. Because of Proposition 2.
where F ϕ , is defined in Condition (iv) of T . In particular, F ϕ is positive decreasing to zero 1 . At this point we employ the classical eigenvalue estimate k j ≥ (j − M )π/L (see for instance [16,10]) and we bound the summation over the integers of a monotone function by its integral.
The function E ϕ : R 2 + → R + is continuous, positive, and monotone with respect to both arguments. In particular E ϕ (J − M, ·) is strictly monotone increasing, while E ϕ (·, Lt) is strictly monotone decreasing to zero. Since χ is an integer number (generally negative) it is enough to ensure that |R J (t)| < 1/2 in order to compute χ. Therefore we define J min (M, Lt) := min{J ∈ N : E ϕ (J − M, Lt) < 1/2}, this is the minimum number of eigenvalues needed for computing χ. This can be further extended to the case where only partial information on Γ is known. The result is summarised in • has total length at most L, i.e. L ≤ L; • has its shortest orbits at least min long, i.e. min ≥ min .
Let ϕ ∈ T with T = 1 and let t = 1/ min . Then J ≥ J min = J min (M , Lt) eigenvalues of the standard Laplacian on Γ suffice for recovering the Euler characteristic of Γ via the formula where nint(x) denotes the nearest integer function. Proof . Hence for J ≥ J min it follows that |χ − S J (t)| < 1/2, and thus we obtain (18).
Notice that the cases S J (t) = χ ± 1/2 is excluded, hence we do not need to worry of the definition of the function nint[·] over the half integers.
We notice that the hypothesis of the above theorem cannot be relaxed too much. In particular, the condition min ≥ min is crucial, because adding small loops to any graph perturbs its smallest eigenvalues only a little, but at the same time it does dramatically affect its topological structure, decreasing χ by 1.
We also notice that the quantity L can be obtained from other information on the graph when available. For instance if the number of edges of the graph is at most N then L ≤ min 2≤j≤J (J + N )π/k j .
Remark 3.2. The above Theorem 3.1 can be reproduced after replacing F ϕ by any other majorant function monotone decreasing. The J min computed will be different. We will use this remark in Section 4.

Calculation with approximate spectrum
Since we have in mind to apply our method in an experimental context we need to consider the case where the eigenvalues are known only up to some approximation. We denote by { k j } J j=1 the approximate spectrum where each eigenfrequency is measured with precision δ, meaning that k j − k j ≤ δ for all 1 ≤ j ≤ J. Hence we define S J (t) = 2 j≤J ϕ t (k j ) the truncated sum evaluated at the approximated eigenvalues, hence our new aim is to approximate χ from S J (t). Therefore to the error contribution E ϕ (J − M, Lt) we need to add a new term This second term can be bounded linearly in δ in the following way Therefore we define For this new function does not necessarily always exists J such that E ϕ (J − M, L, t) < 1/2. For simplicity we consider separately the following two inequalities, which together are sufficient: • E ϕ (J − M, Lt) < 1/4; • 2δJ/t < 1/4.
We derive from them the following two conditions: If the first J eigenvalues are known such that both the above conditions are satisfied then we are able to compute the Euler characteristic via the approximated spectrum. The above analysis can be summarised in the following Theorem 3.3. Let Γ, ϕ and t be as in the hypotheses of Theorem 3.1, and moreover let the first J eigenvalues of the standard Laplacian on Γ be known up to error δ. Assume that J and δ satisfy both conditions (26) and (27) with M, L replaced by M , L. Then the Euler characteristic of Γ can be recovered via the following formula 4 Example of a family of test functions

Fourier transform
Let us recall that we are considering the family of test functions with d a positive integer, called the order. To apply Proposition 2.4, we need to compute the Fourier transform of ϕ d , and also to determine the normalisation constant C d . Let us set, for k ∈ R,

Standard trigonometric identities give us
We can evaluate quickly the integral on the right using complex integration, and assuming k > 2πd. Let us denote by γ the closed path in the complex plane C obtained by adding the path γ 0 , consisting of the segment [−1, 1] traversed in the positive direction, and γ 1 , the semi-circle centered at 0 traversed in the anti-clockwise direction from 1 to −1.
Since γ f (z) dz = 0, we obtain We find the two terms on the right-hand side being the integral on the positive and negative part of γ 0 , respectively. Using the change of variable t = √ s, we transform the expression into where B denotes the Euler beta function. We have where the last equality results from using, 2d + 1 times, the identity Γ(x + 1) = xΓ(x). We obtain Both sides of the equation define a function of k which is holomorphic in C. The equality, which we have proved for k real and positive, therefore holds for k ∈ C by the principle of isolated zeros. In particular, we find and finallŷ e ik/2 sin(k/2).

Error estimate and optimal order
We now apply the analysis in Section 3 to the present case. In the rest of the section, we assume that t = 1/ min . We have , so that, when k > 2πd, Let us write R d,J (t) for the remainder denoted in Section 3 by R J (t), in the case whereφ =φ d . According to Remark 3.2, we find, for J > 2tLd + M , Using the inequality d! ≤ √ 2πdd d e −d+1/6 (see for instance Formula (3.9) in [17]), we find that with ρ := 2tL and Let us note that we have replaced the discrete (integer) values d and J by the continuous parameters α and β.   Using the substitution s := 1 2α ln e 1/6 ρ ε − 1, we see that this derivative vanishes if, and only if, se s = 1. This equation has a unique positive solution which is, by definition, s = W (1). Therefore, the derivative vanishes only at α = α * (ε, ρ), with α * (ε, ρ) = 1 2(1 + W (1)) ln e 1/6 ρ ε .
We can easily check that the derivative is negative below that value and positive above, so that α → β(ε, M, ρ, α) is first decreasing and then increasing. The function d → β(ε, M, ρ, d) takes its smallest value among all positive integers for the largest integer below α * (ε, ρ) or the smallest integer above. Equations (31) and (32) then follow easily.

Numerical examples
We have numerically tested our method on some graphs for which we have computed the exact spectrum up to machine error precision and then added an error term uniformly distributed with prescribed variance. In this section we give an overview of the implemented algorithm and show a few examples of its application.

Code overview
The code consists of three distinct parts each of them developed to run on the software Mathematica [19].
1. The first part is a GUI for building a graph with a drawing, drag-and-drop method and inserting the lengths of the edge.
2. The second part contains an algorithm which computes the eigenvalues of the standard Laplacian of the graph in input. This is based on the von Below formula [20] which holds for equilateral metric graphs, more details can be found in [10]. This code shares parts with an other study in [21] on how this method can be extended to approximate eigenvalues of non equilateral metric graphs. Several samples of approximate spectra { k j } j∈N are realized by summing a pseudo-random number uniformly distributed over a prescribed interval to the exact spectrum k j = k j + e j , e j ∼ U([−δ, +δ]).
3. The last part computes the parameters d * , J * , and Formula (18), plots of the maps t, J → S J (t), S J (t), and compares the error terms (30) with the numerical error | S J (t) − S J (t)|.

Comparison of the theoretical bounds with the experimental error
Here we present some computations of χ performed for the graphs Lasso and K 5 (see Figure 1). Let us explain how the numerical results are displayed in Figures 2 and 3.
Top left : the function t → S J (t) (in black) is compared with a sample of t → S J (t) (in purple) with δ indicated in caption. The curves are not visibly distinguishable despite the different error terms.
Top right : the function t → S J (t) is compared for two different test functions, ϕ d (in black) and the triangular function ψ (in blue). 1. In the two plots at the top right we notice that the black curve approximates χ with higher precision than the blue curve, i.e. it stays within distance 1/2 from χ for a certain interval in t. This shows that the method with ϕ d is significantly more efficient compared with ψ in terms of number of eigenfrequencies sufficient in order to compute χ.

2.
The two plots at the bottom show that the theoretical bound exceeds the numerical error by at least a factor 10 2 . This phenomenon can be explained by the several cancellations occurring inside the tail R d,J (t) which we did not take into account in (29).
3. In the bottom left plots it can be observed that both the theoretical bound and the numerical error grow with respect to the parameter t, indicating that it is important for the effectiveness of the method to have a good lower bound of min , corresponding to a smaller t.

4.
In the bottom right plot we observe that for J growing the numerical error stabilizes at small positive values different for different samples. This suggests that the precision does not improve by taking J J * and that the error mainly depends on the contribution from the first terms in the summation S J (t).
5. The theoretical bound shows to be even worse with respect to the numerical error for J J * , this can be explained by the decreasing magnitude of the Fourier coefficients of the test function ϕ inside | S J (t) − S J (t)|, which we did not take into account in (20-24).

Additional examples
We compare the numerical computation of χ for the following three equilateral graphs K 5 , K 5 with one edge made pendant, i.e. disconnected at one end-point, and K 3,3 . These graphs have Euler characteristic respectively χ = −5, −4 and −3, they have all three been constructed equilateral with each edge of length one. The first two graphs are chosen to have very close spectra. We can numerically observe that only 1 in 5 of their eigenfrequencies differ significantly.   Figure 3: Plots for the graph K 5 : at the top left is t → S J (t), S J (t) for ϕ d ; at the top right t → S J (t) for ϕ d (black) and ψ (blue); at the bottom left t → |χ − S J (t)| (purple) compared with the theoretical bounds (red); at the bottom right the same quantity, but in function of J. In the first three plots J = J * = 41 is fixed and the value t = t = 0.5 is shown with the dashed vertical line. In the bottom right plot t = t = 0.5 is fixed. The approximate spectrum is computed with δ ≈ 1.5 · 10 −3 .  Figure 4: From left to right: K 5 , K 5 with one endpoint disconnected and K 3,3 . On the right, we plot the functions t → S J (t) for each graph, with J = J * = 41, 42, 37 respectively.

Influence of the degree of the test function
• The distribution of the edge length has large variance, i.e. the length of the shortest edge is very small compared to the total length. For instance, a Lasso graph with the loop at least 8 times shorter than the pendant edge.
• The distribution of the edge length has null, or very small variance, but the number of edges is relatively large, at least 17 if the graph does not have loops. For instance, the equilateral K 7 .

A Proof of Proposition 2.4
A priori, we only know that Formula (5) holds for test functions in the Schwartz class. In order to extend it to functions in the set T , we use a standard argument of regularization by convolution. Let us give more details. We choose a function ρ : R → R of class C ∞ , non-negative, whose support is in [−1, 1]. We also require that ρ is even and that +∞ −∞ ρ( ) d = 1. It follows thatρ is real-valued and even, withρ(0) = 1 and |ρ(k)| ≤ 1 for all k ∈ R. For any ε > 0, we define ρ ε ( ) := (1/ε)ρ( /ε). We have ρ ε (k) =ρ(εk) for all k ∈ R.