The many facets of the Estrada indices of graphs and networks

The Estrada index of a graph/network is defined as the trace of the adjacency matrix exponential. It has been extended to other graph-theoretic matrices, such as the Laplacian, distance, Seidel adjacency, Harary, etc. Here, we describe many of these extensions, including new ones, such as Gaussian, Mittag–Leffler and Onsager ones. More importantly, we contextualize all of these indices in physico-mathematical frameworks which allow their interpretations and facilitate their extensions and further studies. We also describe several of the bounds and estimations of these indices reported in the literature and analyze many of them computationally for small graphs as well as large complex networks. This article is intended to formalize many of the Estrada indices proposed and studied in the mathematical literature serving as a guide for their further studies.


Introduction
At the dawn of the XXI century the current author proposed an index to quantify the "degree of folding" of a linear chain in a three-dimensional space [70]. The motivation of this work came from the fact that many scientific articles make claims like that the structure A "is more folded than" the structure B (see examples at: [44,67,128,237]), or that certain structure is "highly folded" (see for instance: [42,129,142,246]), etc. These expressions could be referring to protein or polymer structures, but also to brain regions or even geological structures (see previous refs.). However, in neither of these works there was an index that quantifies how folded a linear chain is. Thus, the author proposed the index I 3 = n j=1 exp λ j (W ) , where B Ernesto Estrada estrada@ifisc.uib-csic.es λ j (W ) are the eigenvalues of certain tridiagonal matrix W whose diagonal entries are related to the cosines of the dihedral angles between adjacent planes and W i,i+1 and W i+1,i are equal to one. This index characterizes very well the degree of folding of a geometric chain and it has been mainly applied to the study of the degree of folding of proteins (see for instance [71,73,211]), although it can be applied to the folding of any linear chain. Five years after the publication of the "folding degree" paper, the authors of [88] proposed the "subgraph centrality" as a way to characterize the importance of the nodes in a complex network. "Complex networks" are large graphs representing the skeleton of complex systems in social, ecological, cellular, molecular, infrastructural, semantic and other scenarios [78]. The subgraph centrality of a node v in a network is defined as SC v = n j=1 ψ 2 jv exp λ j (A) , where λ j = λ j (A) are the eigenvalues of the adjacency matrix of the graph and ψ jv is the vth entry of its jth normalized eigenvector. Then, the so-called subgraph centralization of the network is v SC v = n j=1 exp λ j (A) [88], which is similar to the folding degree I 3 . In June 2005 the current author presented the lecture "Topological characterization of complex networks" at the International Academy of Mathematical Chemistry in Dubrovnik, Croatia. As a consequence Ivan Gutman proposed to organize a small seminar at a park near the port of Dubrovnik to discuss some of the mathematical aspects of the index v SC v = n j=1 exp λ j for general graphs. As a result, a paper was published in 2006 in Croatica Chemica Acta introducing v SC v as a molecular structure descriptor [113]. A year later the paper "Estimating the Estrada index" was published, where the authors proposed to call E E (G) = n j=1 exp λ j the Estrada index [54]. The same year a statistical mechanics interpretation of E E (G) as the partition function of a graph [83] appeared. A year later, in 2008, there were more than 30 papers published in the mathematical literature containing "Estrada index" in the title.
It seems a priori that E E (G) has emerged in different, apparently unrelated, scenarios: folding of linear chains, subgraphs in networks, and partition function in statistical mechanics. This reminds us the story told by Eugene Wigner in the first paragraph of his paper "The unreasonable effectiveness of mathematics in the natural sciences" [233] where a fellow asked a former classmate, now a statistician, about a symbol in a paper dealing with population trends. The statistician replied that the symbol was "π" and to clarify the skepticism of the other he added that it is "the ratio of the circumference of the circle to its diameter." The fellow then replied more skeptical: "Well, now you are pushing your joke too far, surely the population has nothing to do with the circumference of the circle." The situation of the Estrada index seems murkier than the one in that story, particularly after the ad hoc definition of several other variations of the index based not on the eigenvalues of the adjacency matrix, but of the graph Laplacian, distance matrix, resolvent of the adjacency matrix, Hadamard pseudo-inverse of the distance matrix (a.k.a. Harary matrix), Mittag-Leffler matrix functions of A, etc.
The goal of this paper is to make an account of the different facets of the Estrada indices. In doing so we will provide contextualization of several of these indices, many of which have been proposed in an ad hoc way. Therefore, we will provide a physical and/or mathematical context and interpretation of these indices. They include a combinatorial interpretation based on counting subgraphs, a statistical mechanics approach, a probabilistic interpretation in the context of walk-regular graphs, an interpretation on the basis of oscillations in (quantum and classical) systems of ball-and-springs, a contextualization on the basis of epidemiological models (normal and fractional) on graphs, diffusive processes with negative diffusiveness, nonlocal processes on graphs, quantification of graph radius of gyration. Although this paper does not intend to describe all the results published in the literature on this topic we make an account of many of the different bounds and estimations of the Estrada, Seidel Estrada, Harary Estrada, Laplacian Estrada, resolvent Estrada, Mittag-Leffler Estrada, and distance Estrada indices. For this purpose we include some numerical analysis of these bounds in the set of 11,117 connected graphs with 8 nodes and in five real-world networks representing a variety of complex system scenarios. The paper is written in a way that intend to be selfcontained and make the necessary definitions for understanding the concepts used in it. The paper is then intended as a guide for further studies and developments in this area of spectral graph theory.

General definitions
Here we present some definitions which are used across the paper and settle down the notation. We consider here simple, connected graphs G = (V , E) with n nodes (vertices) and m edges.
Definition 1 A walk of length k in G is a set of nodes and edges v 1 , e 1,2 , v 2 · · · v k−1 , e k−1,k , v k such that for all 1 ≤ l ≤ k, (v l , v l+1 ) ∈ E. A closed walk is a walk for which v 1 = v k+1 .
Definition 2 A path of length k in G is a walk in which neither vertices nor edges are repeated. A cycle is a closed path. The length of the shortest path connecting two vertices v and w is the (topological) shortest path distance d vw between the two nodes. The diameter of G is the longest distance between two vertices of G.
An induced subgraph is a subgraph formed by a subset of the vertices of the graph and all of the edges connecting pairs of vertices in that subset.

Definition 4 A graph G = (V , E)
is connected if there is a path between every pair of nodes v, w ∈ V . If the graph is directed we said that it is strongly connected if there is a directed path between every pair of nodes v, w ∈ V . A (strongly) connected component in a (directed) graph is a subgraph in which any two vertices are connected to each other by (directed) paths, and which is connected to no additional vertices in the rest of the graph.

Definition 5
The degree of a node v is the number k v of edges incident with that node. A graph is regular if the degree of all its nodes is the same.
The following matrices will be considered ( Table 1): Other matrices such as the Seidel adjacency matrix and Harary matrix, are defined in situ in the corresponding sections of the paper. The following types of graphs are used in this work. Table 1 Definition of some matrices used in this paper Name Symbol Definition Spectrum Adjacency Distance D D i j = d i j i = j 0 i = j σ 1 ≥ · · · ≥ σ n -Complete graph of n vertices K n : the graph having an edge between every pair of vertices.
-Empty graph of n verticesK n : the graph having n vertices and no edges.
-Complete bipartite graph K n 1 ,n 2 : the graph with n = n 1 + n 2 vertices in which the vertex set is partitioned into two disjoint subsets of cardinalities n 1 and n 2 , respectively, such that every vertex in one set is connected to every vertex in the other set. -Star graph S n : the particular case of K n 1 ,n 2 in which n 1 = 1 and n 2 = n − 1.
-Path graph of n vertices P n : the connected graph in which every vertex has degree 2, except two vertices which have degree one. -Cycle C n : a connected graph in which every vertex has degree 2.
Finally we consider two kinds of random graphs.
-Erdős-Rényi (ER) G (n, p) [68] graph with n nodes: constructed by connecting nodes randomly in such a way that each edge is included in G (n, p) with probability p independent from every other edge. -Barabási and Albert (BA) one [21]: created on the basis of a preferential attachment process. The graph is constructed from an initial seed of m 0 vertices connected randomly like in an Erdős-Rényi G (n, p). Then, new nodes are added to the network in such a way that each new node is connected to c ≤ m 0 of the existing ones with a probability that is proportional to the degree of these existing nodes.
the following measures can be defined (3.2) where c k are coefficients which give more weight to the shorter than to the longer walks. Then, if c k = (k!) −1 : The previous result implies that we can express E E v as a weighted sum of subgraphs, which gives the index its name. However, as we are focused here on the Estrada index let us move to the fact that the Estrada index is the sum of the subgraph centralities of all nodes in the graph: (3.4) The sum of node centralities in a graph is known as the corresponding centralization of the graph, or simply as a graph-theoretic invariant. Therefore, the Estrada index of the graph can be seen as its subgraph centralization. Theorem 2 Let G be a (directed) graph and let F be the set of all (strongly) connected subgraphs of G, and let us designate the cardinality of the set F by η. Then, E E (G) = η l=1 c l F l , (3.5) where F l ∈ F and c l ∈ Q.
Proof Using Lemma 1 we can show that that M k = tr A k can be expressed as a weighted sum of (strongly) connected subgraphs. The weight of each subgraph is given by the number of closed walks of length k in the given subgraph. Then, grouping together all identical subgraphs and summing their weights we obtain the final result. For instance, let us consider the first seven powers of the adjacency matrix. Then, tr The expressions for calculating these subgraphs are given in the Appendix as adapted from [9]. The formula for F 20 is given here by the first time.

Some elementary properties of the Estrada index
Before proceeding to more complex properties of the Estrada index let us state a few elementary ones that could be helpful in understanding the structural nature of this index. The reader is referred to the following references [54,57,112,116] for details and references. Lemma 3 Let G be a simple graph and let G − e the same graph from which edge e has been removed. Then (3.13)

Corollary 1 Let G be a simple graph and let T be a tree with the same number of nodes as G.
Then (3.14) Theorem 3 [53,56] Let G be a simple connected graph with n nodes. Then Theorem 4 Let G be a simple graph with n nodes. Then The Estrada indices of some elementary graphs are given below.

Numerical analysis
We consider here two datasets which will be used in the rest of the paper for the numerical evaluation of the different indices and bounds. The first one consists of the 11,117 connected graphs with 8 nodes. The second one is formed by five real-world networks, which correspond to a food web at Stony stream, a network of the neurons in the worm C. elegans, the proteinprotein interaction network of yeast, a representation of the Internet at the autonomous system (AS) level, and a network of the USA western power grid system. The number of nodes n, of edges m, the maximum degree of the nodes k max , and the diameter d max of each network are given in Table 2.
The main goal of these numerical experiments is to show how close the bounds reported in the literature are to the actual values of the Estrada index. This is done because in most of the papers where these bounds are proposed there are no numerical experiments to illustrate this relation. When possible we will find some connection between structural characteristics of the networks studied and the corresponding bounds analyzed to understand why are they close or far away the actual values of the Estrada index. First, we consider the deviation of the bound from the actual value as |E E exact − E E bound | /E E exact expressed as percentage. We do this calculation considering the bound given in Lemma 2 for all the connected graphs with 8 nodes. The histogram illustrating the number of graphs having a given relative deviation (frequency) among the 11,117 connected graphs with 8 nodes is illustrated in Fig. 2. We should remark that we use here the terms "good bound" or refer to a bound as "better than" another just on the basis of the deviation of this bound relative to the actual value of the index. This is used only as a guide as for many cases there is large room for improvement as some of the bounds reported are orders of magnitude further from the real values of the indices.
The mean deviation is 5.768 ± 4.169, which indicates that this bound is a good estimation of the Estrada index for these small graphs. The largest deviation is 40.352 obtained for the complete graph K 8 . In general, the most densely connected graphs are richer in small subgraphs than the poorly dense ones, which increases the relative deviation of this bound for these graphs.
In Table 3 we illustrate the results for the five real-world networks. The largest deviation occurs for the Internet at AS indicating that in this network there are many larger subgraphs with important contribution to the Estrada index. On the other hand, the bound is very close to the actual value for the power grid of western USA, which points out that the Estrada index of this network is well approximated by counting the number of the 21 subgraphs

Estrada index and matrix functions
Soon after the definition of the Estrada index and the subgraph centrality several authors started to be interested in these indices due to their clear relation to functions of the adjacency matrix. The study of matrix functions is an active area of research in (numerical) linear algebra [25,97,127,222]. The topic of matrix functions in network theory has been recently reviewed by the authors of [28]. Therefore, we will not give too many details here and the interested reader is directed to the excellent review [28]. The goal of this section is then to establish the connection between the Estrada indices and functions of the corresponding matrices which pave the way for further sections of the article. Here we will follow the book [127]. Let M be any graph-theoretic matrix, e.g., adjacency, Laplacian, distance, etc. Then, its Jordan canonical form is given by where where Z is nonsingular and m 1 + m 2 + · · · + m p = n.

Definition 6
Let λ 1 , . . . , λ s be the distinct eigenvalues of M and let and let n i be the order of the largest Jordan block in which λ i appears, which is called the index of λ i . The function f is defined on the spectrum of M if the values Another, equivalent, definition is given via the Cauchy integral.
here f is analytic on and inside a closed contour Γ that encloses the spectrum of M.

Estrada index and spectral graph theory
An obvious connection exists between the Estrada index and the area of algebraic graph theory. Algebraic graph theory [24,30,105] deals with the use of algebraic methods to solve problems about graphs. Of particular interest is the use of the spectra of graph theoretic matrices to understand the structure of graphs, which is known as spectral graph theory [46,[50][51][52]213,214]. This area of research started in an applied context when Collatz and Sinogowitz published their paper entitled: "Spektren endlicher grafen" motivated by application problems such as the vibrations of a membrane [223]. Let us consider a simple example of the connections between structural properties of graphs and their spectra: counting triangles in a graph. The number of triangles, which is a combinatorial property of the graph, can be obtained from the spectrum of the adjacency matrix as: 1 6 n j=1 λ 3 j , where λ j are the eigenvalues of the adjacency matrix. The field of spectral graph theory had a tremendous impulse in the 1970s due to its connection with electronic properties of conjugated molecules [59,95,124,215,216,219].
The relation between the trace of a matrix and its eigenvalues immediately implies that the Estrada index of a graph can be expressed in terms of the eigenvalues of A as follows: In general, the exponentiation of A enlarges the spectral gap λ 1 − λ 2 and contracts the negative part of the spectrum. On the contrary, exp (−A) largely contracts the positive part of the spectrum and enlarges its negative part. These simple dilation/contraction effects of the main parts of the spectrum of A have important consequences on the Estrada index of a graph as we will see in the next parts of this review.
The analysis of the relation between the spectrum of a graph, i.e., the eigenvalues of its adjacency matrix, and the structure of the graph is the main goal of spectral graph theory.
One of the first results on spectral graph theory related to the Estrada index was the following bounds obtained by the authors of [54].

Theorem 5 Let G be a simple graph with n nodes and m edges. Then, the Estrada index of G is bounded as
with equality attained if and only if G ∼ =Kn.
These bounds were further improved in [166] where the following was proved.

Theorem 6
Let G be a simple graph with n nodes and m ≥ 1 edges. Then, the Estrada index of G is bounded as Based on Gauss-Radau quadrature rule the authors of [27] obtained the following bounds.

Theorem 7
Let G be a simple graph and let a, b ∈ R be such that the spectrum of A is contained in [a, b]. Then, the Estrada index of G is bounded as where k i is the degree of the node i.  Other bounds have been proposed, specially lower bounds, for the Estrada index. Some examples are given below.
Theorem 9 [247] Let G be a simple graph with n nodes and let Z = n i=1 k 2 i . Then, the Estrada index of G is bounded as with equality attained if and only if G ∼ = K n or G ∼ =Kn.
Theorem 10 [110] Let G be a simple graph with n nodes and m edges either without isolated vertices or having the property 2m/n > 1, then, the Estrada index of G is bounded as with equality if and only if G is a regular graph of degree 1.
Theorem 11 [110] Let G be a simple graph with n nodes and m edges, such that 2m/n < 1.

Then, the Estrada index of G is bounded as
where equality holds if and only if G consists of n − 2m isolated vertices and m copies of K 2 .
Theorem 12 [110,119] Let G be a simple graph with n nodes, m edges and graph nullity η 0 . Then, the Estrada index of G is bounded as where equality is attained if and only if n − η 0 is even, and if G consists of copies of complete bipartite graphs K r i ,s i , i = 1, . . . , (n − η 0 ) /2, such that all products r i · s i are mutually equal.
Theorem 13 [190] Let G be a simple graph with n nodes, m edges and minimum degree k min . Then, the Estrada index of G is bounded as with equality if and only if G ∼ = K p, p ∪ K 1 with n = 2 p + 1.
Theorem 14 [190] Let G be a simple graph with n nodes, m edges and minimum degree k min . Then, the Estrada index of G is bounded as E E (G) ≥ 2 cosh 2 cos π n + 1 + n − 2, (5.11) with equality if and only if G ∼ = P 2 or G ∼ = P 4 .
Theorem 15 [19] Let G be a simple graph with n nodes, m edges and t triangles. Then, the Estrada index of G is bounded as E E (G) ≥ n 2 + mn + 2nt, (5.12) with equality if and only if G ∼ =Kn.
Other bounds reported in the literature are based on different graph-theoretic indices and properties or for specific classes of graphs. A non-exhaustive resume is provided in Table 4.

Numerical analysis
We now do some calculations to show how close to the actual values of the Estrada index are some of the bounds studied in the previous sections. In particular, we consider the following five bounds: Bound 1 (Theorem 5); Bound 2 (Theorem 6; Bound 3 (Theorem 7 using a = −λ 1 and b = −λ n ); Bound 4 (Theorem 7 using a = −k max and b = k max ); Bound 5 (Theorem 8). First, we study these bounds for the 11,117 connected graphs with 8 nodes. The histograms of the relative deviations of these bounds are illustrated in Fig. 3, where the lower bound is always drawn in blue and the upper one in red. The means and standard deviations of the General [10,38,63,101,121,189,198,201,238] Weighted general [197,200] Trees [55,62,159,188,244] Molecular trees [115,134] Unicyclic [64] Bicyclic [228] Tricyclic [252] Tetracyclic [186] Pentacyclic [185] Bipartite [91,120,245,250] Line graphs [4,208] Strongly quotients [33] Folded hypercubes [165] Cacti [157] Cayley [103] Specific graphs [104] Ramanujan [199] Benzenoids [118] Phenylenes [187] Fullerenes [14] Möbius [96] lower, upper bounds are as follow: In Fig. 4 we illustrate the results for the five real-world networks considered in this work. In general, with the exception of Bound 3, which is based on eigenvalues, and Bound 5, which uses tr A 4 , the rest of the bounds are very far from the actual values for these four networks. With these two exceptions, the upper bounds exaggerate dramatically the estimation, in particular the Bound 1. Bound 4, performs very badly when the maximum degree of the network is very high and not close to the spectral radius, which is the case for instance of Internet, but also of many real-world networks. All in all, these results point out to the necessity of improving the bounds for the Estrada index of large graphs.
We then consider simple bounds based on the spectral radius of the adjacency matrix λ 1 . That is, The results are given in Table 5. As can be seen the bounds are very close to the actual values of the Estrada index. This is a consequence of the relatively large values of the spectral radius and of the spectral gap observed in most of the real-world networks, which when exponentiated are significantly enlarged. Notice that the largest deviation is obtained for powergrid, where the spectral radius is significantly smaller than in the rest of the networks and the spectral gap is very small.

Random graphs
In the study of real-world networks it is desired to investigate how unique are their structural and dynamical properties in relation to some null model. For instance, suppose that we have found that certain network displays relatively large Estrada index in relation to other networks of the same size. Is this a characteristic feature of the topological organization of this network or just an artifact emerging from a random interconnection of their nodes? A way to investigate this is by comparing the Estrada index of these networks with those of random realizations of such networks with the same number of nodes and edges. Then, the use of random graphs is frequent in the analysis of real-world networks [220]. Two classical models, although not the only ones, to do such studies are the Erdős-Rényi random graphs [68] and the Barabási-Albert preferential attachment model [21]. For instance, the Estrada index of the network "neurons" studied here is E E (G real ) ≈ 1.3062 × 10 10 and that of an Erdős-Rényi random graph with the same number of nodes and edges is E E (G ER ) ≈ 3.4688×10 6 , which indicates that the large Estrada index of this network is not due to a random interconnection of the neurons of C. elegans. However, the consideration of a Barabási-Albert network with the same number of nodes and edges than those in the network "neurons" gives E E (G BA ) ≈ 1.2131 × 10 10 , which clearly points out that the relatively large Estrada index of this network may be explained by its skewed degree distribution. For the Estrada index of random graphs, only the Erdős-Rényi model has been considered so far, indicating the necessity of extending these studied to other classes of random graphs such as the Barabási-Albert one. The following estimates were found for Erdős-Rényi random graphs based on the number of nodes and the probability of connection.
Lemma 4 [196] Let G n, p be an Erdős-Rényi random graph with n nodes and probability ln n n p < 1 − ln n n . (5.14) Then, the Estrada index is almost surely as n → ∞.

Theorem 16 [43] Let G n, p be e an Erdős-Rényi random graph with n nodes and probability p. Then, the Estrada index is
almost surely (a.s.) if and only if lim n→∞ n 2 /n 1 = 1.
In the case of Erdős-Rényi random bipartite graphs the author of [206] proved the following bounds for the Estrada index.

Estrada index and statistical mechanics
The analogy of the Estrada index E E (G) = tr e A with the partition function of a quantum system Z = tr e −τ H (see further for definitions) is remarkable, and was noticed soon after the definition of this index [83]. The importance of establishing this connection is twofold. On the one hand, the index can be interpreted in a physical context which at the same time facilitates its interpretation in other contexts where it is applied. On the other hand, new tools and techniques from statistical mechanics can be used to enrich the theory behind this index.
Here, we will describe the statistical mechanics interpretation of the Estrada index. Let us consider a physical system S that can be represented by a graph G, such that the total energy E of S can be obtained by the time-independent Schrödinger equation:ĤΨ = EΨ , where Ψ is the wavefunction andĤ is the Hamiltonian describing the interactions between the elements of S. In certain approaches in physics and chemistry, it is customary to use an effective Hamiltonian which describes the interaction between nearest-neighbors (NN) in the systemĤ where α is a self-energy parameter for the elements of S and t N N is the energy of the interaction between pairs of adjacent elements. In Chemistry this model is known as the Hückel Molecular Orbital (HMO) method [154,239], while in Physics it is better known as the tight-binding approach [184]. The parameter t N N is negative as it is supposed to be an attractive interaction. Therefore, it is common to set α = 0 and t N N = −1, such that H = −A. Therefore, the energy levels of the system are E j = −λ j and the wavefunctions correspond to the eigenvectors associated to the eigenvalues of A.
In the statistical mechanics framework [23,69], the Boltzmann probability p j (τ ) of finding the system in a state with energy E j when the inverse temperature of the system is τ = (k B T ) −1 > 0 with k B being a constant and T being the temperature 1 is where Z = tr e −τĤ N N . Therefore, the Boltzmann probability of the system is given by where the Estrada index plays the role of the partition function of the graph. We now can define the entropy of the graph as [83] S which in general is bounded as follows.

Lemma 5 Let G be a simple graph. Then, the free energy of G is bounded as
where the upper bound is attained for the null graph K n and the lower bound is reached for the complete graph K n .
From the general expression of the entropy one can obtain the graph "enthalpy" H (G, τ ) = − j p j λ j and the graph free energy, which is sometimes named the natural connectivity of the network [83]: We can write the logarithm of the Estrada index as follows, which implies that where = λ 1 − λ 2 is the spectral gap. Therefore, we have proved the following.

Lemma 6 Let G be a simple graph. Then, the free energy of G is bounded as
More generally, the free energy of a graph can be bounded by using the many bounds obtained for the Estrada index which have been previously reported in the literature. One important example is the following [83].

Lemma 7 Let G be a simple graph. Then, the free energy of G is bounded as
where the lower bound is obtained for the complete graph K n and the upper bound for the null graph K n .

Numerical analysis
We consider here numerical experiments to illustrate some general characteristics of the indices described in the previous section. We report the change of the entropy, enthalpy and free energy of all connected graphs with the increase of the number of edges in the connected graphs with 8 nodes, i.e., its edge density. It can be seen in Fig. 5, as expected, that the three parameters decay with the increase in the edge density. However, it should be noticed that for graphs with exactly the same number of edges there is a wide variability in these parameters, particularly for the entropy. The readers interested in more details about the implications of these parameters on the structure of graphs are referred to [83]. We then computed the three statistical mechanics parameters for the five networks studied here. The results are in Table 6 where we also give the values of the edge density of these graphs: δ (G) = 2m/ (n (n − 1)) where n and m are the number of nodes and edges of the graph. The most densely connected network, Stony, displays the lowest entropy and the least dense, powergrid, displays the largest one. However, as can be seen for the intermediate values of δ (G) this trend is not always followed as there are other structural factors influencing these statistical mechanics parameters. For instance, the network of Internet at AS displays the second smaller entropy of all the networks and the lowest free energy of all, although it is not very dense. Table 6 Values of the entropy, enthalpy and free energy of the five real-world networks analyzed here

Marginal probability, walk entropy and walk regularity
Having in mind the importance that the probability p j (τ ) has in the definition of statistical mechanics properties of networks we propose to explore it further in this section. That is, we consider here the role of the Estrada index in defining some probability-based measures for graphs. Let us start with two definitions from basic statistics (see for instance [49,Ch. 2]).

Definition 9
The conditional probability P (A |B ) is the probability that the event A occurs given that the event B occurs.

Definition 10
The marginal probability is the unconditional probability of one event A. That is, the probability that A occurs regardless of whether B occurs or not.
To obtain the marginal probability of an event A one should sum all possible configurations of the other event to obtain a weighted average probability Let us then return to the time-independent Schrödinger equation: where E j are the energy levels of the system and ψ j are the corresponding eigenfunctions. As usual, ψ j,v 2 represents the probability of finding a quantum particle at a given vertex v and time conditional to the system to be at the energy level described by the wave function ψ j . That is, ψ j,v 2 = P (v | j ) using the notation defined before.
On the other hand, p j (τ ) which was defined in the previous section accounts for the probability that the system is at the jth energy level for a given τ . Then, fixing τ , p j (τ ) = P ( j) . Therefore, the marginal probability that the node v is occupied by the quantum particle independently of the energy level in which the system is, is given by: which can be expressed as [86]: The corresponding entropy, known as the walk-entropy of the graph [86], is defined using Shannon formula: We now consider a graph property known as walk-regularity and the role that the walk entropy play in its characterization. Let us introduce the concept of walk regularity first (see for instance [106]).

Definition 11
A graph is walk-regular if ∀i, j ∈ V and for every nonnegative integer r , The following conjecture was formulated in [86] as an extension of the conjecture related to the subgraph centrality which had been previously stated in [88].

Conjecture 1 A graph is walk-regular if and only if S w (τ ) = ln n for all τ > 0.
Let us then introduce some necessary concepts for the further developments in the proof of this conjecture.

Definition 12 Two vertices
Definition 13 A graph is τ -subgraph regular if all pairs of vertices are τ -subgraph equivalent.
The following result was a step forwards the proof of Conjecture 1.
Theorem 18 [26] A graph G is walk-regular if and only if G is τ -subgraph regular for all τ ∈ I ⊆ R, where I is any set of real numbers containing an accumulation point.
In the saga, in [151] the authors found some counterexamples to a new conjecture proposed in [26] and stated a new conjecture. The final proof of Conjecture 1 came from an elegant Theorem in 2021 [18] where the authors used results from the Lindemann-Weierstrass Theorem.
Theorem 19 [18] Let τ > 0 be an algebraic number and let G be a connected undirected graph with adjacency matrix A.

G is τ -subgraph regular if and only if G is walk-regular.
2. If two vertices i, j are τ -subgraph equivalent, then the degree and eigenvector centralities of i and j are equal.

If G is τ -subgraph regular, then the degree and eigenvector centralities are also identical for all nodes.
Walk regular graphs can be constructed by using Kronecker product of the adjacency matrices of two walk-regular graphs [106]. That is, if G 1 and G 2 are walk regular graphs, then G 1 ⊗G 2 is also walk regular. Therefore, we have the following result.
Proposition 1 [86] Let G 1 and G 2 be two simple graphs with n 1 and n 2 vertices, respectively. Then, for all τ > 0 if G 1 and G 2 are walk-regular.

Bipartivity, signed graphs and Seidel Estrada index
is bipartite if its set of nodes V can be split into two subsets V 1 and V 2 such that there are edges only between the two sets but no edge connects vertices in neither V 1 nor V 2 . Therefore, a graph is or is not bipartite. However, in certain real-world situations a graph can be "close to bipartite", meaning that by removing very few edges the graph become bipartite. This is the case, for instance, of human sexual contact networks and human romance or partnership networks as remarked in [130]. In 2003 the authors of [130] proposed to quantify the "bipartivity" of a graph. The first of their measures is defined by where m f is the number of edges that if removed the network becomes bipartite. 2 The calculation of this index is computationally intractable as it is NP complete. The authors [130] then proposed another index in which m f is assessed computationally. Here we will show how the use of the Estrada index of graphs allows the calculation of an index of bipartivity which depends only on the eigenvalues of the graph. The first of these approaches was published in [87] and will not be discussed here. Instead we will consider the index studied in [82]. Another measure of bipartivity was also proposed in [180]. We will start with some basic definitions for the sake of completeness of this section. A bipartite graph is characterized by the following result proved by Konig in 1916 [153] (see also [15]).

Corollary 2 A graph G is bipartite if and only if it contains no closed walks of odd length.
The Estrada index of a graph can be expressed in terms of the hyperbolic matrix functions as The tr (sinh (A)) counts the odd-length closed walks in the graph: Similarly, tr (cosh (A)) counts the even-length closed walks. An odd closed walk of any length in the graph exists if and only if the graph contains at least one odd-length cycle. Therefore, we can reformulate the previous Corollary as.

Corollary 3 A graph G is bipartite if and only if tr (sinh (A)) = 0.
Based on this result the authors of [82] proposed the following.

Definition 14
The bipartivity of a graph is defined as the relative difference between the number of closed walks of even and odd length, where G − is the graph in which all the edges are weighted by −1.
It is easy to see that tr (exp (−A)) reaches its minimum for the complete graph, which is also the graph for which E E (G) is maximum (see an example in Fig. 6). In this figure the reader can also visualize how the bipartivity index changes monotonically with the increase of the number of edges "frustrating" the bipartition. Then, we have the following result.

Lemma 8 Let G be a simple graph. Then, its bipartivity is bounded as
where the upper bound is attained for any bipartite graph and the lower bound is reached for G ∼ = K n .
Therefore, we have that

Signed graphs
In order to understand why the index b (G) quantifies the bipartivity of a graph we should start by recognizing that the numerator of b (G) is the trace of the adjacency matrix of a fully-negative signed graph. For an exhaustive compilation of mathematical results about signed graphs the reader is referred to [241]. Let us introduce here the necessary concepts for understanding the connections between bipartivity and signed graphs. We will start with the following. Therefore, the numerator of b (G) counts the number of negative cycles in G, where a negative cycle is any cycle in which the product of the sign of its edges is negative. In a fully-negative graph, the negative cycles are the odd-length cycles, which are indeed those that break the bipartivity of the graph. In the theory of signed graphs we have the following important concept (for a list of references and some critical account see [79]).

Definition 16
A signed graph G +− is balanced if all its cycles are positive. Then, it is obvious that a fully-negative graph is balanced if and only if it is bipartite. In the general case of any signed graph the following result is well-known.

Theorem 21 A signed graph G +− is balanced if and only if its nodes can be separated into two mutually disjoint sets, such that positive edges joint nodes only inside the subsets and negative edges joint nodes from different subsets.
The adjacency matrix of a signed graph can be expressed as: A = A + − A − , where A + represents the adjacency between pairs of nodes connected by positive edges, and A − represents the adjacency between pairs of nodes connected by negative edges. Definition 17 [79,81] The balance of a signed network with adjacency matrix A = A + − A − can be quantified by where |·| represents the entrywise absolute of the corresponding matrix.
The following result was proved in 1980 [2].

Theorem 22
For any signed graph, the matrices A + − A − and A + − A − are isospectral (cospectral) if and only if the signed graph is balanced.
Then, we have the following.
where the upper bound is attained for any balanced graph and the lower bound is reached for a fully-negative complete graph.
Then, we also have that which is a maximally unbalanced graph.

Seidel Estrada index
Let us focus now on a particular kind of signed graph. Let J and I be the all-ones and identity matrices, respectively. The following matrix was introduced in [221] and it is nowadays known as the Seidel matrix.

Definition 18
The Seidel matrix of a simple graph G with adjacency matrix A is defined as Therefore, we have the following result.

Theorem 24 Let G +− be a signed graph with adjacency matrix S (G). Then, G is balanced if and only if S (G) is isospectral to A (K n ).
Proof The balance index of a signed graph with adjacency matrix S (G) is (8.11) which immediately implies the result.

Remark 3
The term tr (exp (S (G))) =: S E E (G) was denoted in [122] as the Seidel Estrada index of the graph. The name Seidel honors mathematician Johan Jacob Seidel (1919-2001). 4 We can prove here the following result.
Theorem 25 Let K n 1 ,n 2 be a complete bipartite graph. Let S K n 1 ,n 2 be the Seidel matrix of K n 1 ,n 2 . Then, S K n 1 ,n 2 and A K n 1 +n 2 are cospectral.
Proof Using the structural balance theorem we can show that the signed graph whose adjacency matrix is S K n 1 ,n 2 is balanced. That is, we can split the set of nodes into two disjoint sets containing n 1 and n 2 nodes, respectively, in which the inter-set edges are negative and all intra-set edges are positive. Then, using Theorem 24 we prove the result.

Remark 4
The previous result implies that any signed graph with adjacency matrix In [122] it was proved the following results for the Seidel Estrada index.

Theorem 26 Let G be a simple graph with n ≥ 2 nodes, m edges, t triangles and Z
Theorem 27 Let G be a simple k-regular graph. Then, Theorem 28 Let G be a simple k-regular bipartite graph. Then, In this subsection we have shown that although the so-called Seidel Estrada index was proposed and studied in a completely ad hoc way, it can be connected with the theory of signed graphs. This may facilitate further studies of this index, its extension to consider statistical mechanics parameters and its applications to the study of real-world signed graphs.

Negative absolute temperatures and the Onsager Estrada index
In the definition of the bipartivity index we have considered in the numerator of Eq.
In the context of statistical mechanics which we have analyzed in Sect. 6 this corresponds to consider the inverse temperature τ = −1. So far, we have considered the inverse temperature τ to be positive. So, what a negative inverse temperature could mean? Let us first analyze what is the physical definition of τ . Let S be the statistical entropy, which is a function of the possible microstates of the system, and let E be the system's energy. Then, the absolute temperature is defined as: Graphically, it corresponds to the slope of the curve of entropy versus energy at a given point. Therefore, as can be seen in Fig. 7 the inverse temperature can be negative. In a system at negative temperature the high-energy states are more occupied than low-energy states. Such systems have been created by physicists in the real-world [35].
From a graph-theory perspective what it means that "the high-energy states are more occupied than low-energy states"? In the Sect. 6 we have considered that the Hamiltonian describing the graph as a quantum system is given by the negative of the adjacency matrix H N N = −A, such that the energy levels of the system are E j = −λ j and the wavefunctions are the eigenvectors associated to the eigenvalues of A. In this case the partition function of the graph is given by Z = n j=1 e τ λ j with τ > 0. Therefore, for τ → ∞, we have that Z = e τ λ 1 . In the current case, where τ < 0, we have that when τ → −∞, the partition function is: Z = e τ λ n . This means that we have changed the "importance" given to the different eigenvalues in the Estrada index, giving now more weight to the contributions of the smallest ones. Because Lars Onsager (1903Onsager ( -1976 was the scientist who first study the negative absolute temperatures in [178] we propose to name the following index in his honor. 5  -O E E K n 1 ,n 2 = E E K n 1 ,n 2 = n 1 + n 2 − 2 + 2 cosh √ n 1 n 2 ;

Definition 19 The Onsager Estrada index of G is defined as
-lim n→∞ E E (C n ) = n I 0 , n even, where I 0 = 1 π π 0 e 2 cos x dx; -lim n→∞ E E (P n ) = (n − 1) − 2 cosh (2). As we can see only bipartite subgraphs make a positive contribution to the Onsager Estrada index.

Numerical analysis
Here we compute the bipartivity index for all connected graphs with 8 nodes. We select two other network parameters to compare with the bipartivity. The first is the edge density δ (G) = 2m/ (n (n − 1)) where m is the number of edges. The reason for selecting this parameter is that as the density of the graph increases the number of cycles of any length will also increase. For instance, in Erdős-Rényi random graphs we can find that the number of triangles F 4 (see Fig. 1) is bounded as The second parameter is the clustering coefficient C (G), which is defined as C (G) = 3F 4 /F 3 , where F 3 is the number of paths of length 2 in the graph (see [78]). Here again we would expect that the bipartivity and the clustering coefficient are negatively correlated due to the fact that the increase in clustering means the relative increase in the number of triangles. However, bipartivity is also related to other odd-cycles in the graphs and we want to investigate their influence of this network parameter.
In Fig. 8 we plot the results of the bipartivity vs. the clustering coefficient where the points are colored according to the number of edges that the graph has. As can be seen the most dense graphs also have the highest clustering and lowest bipartivity, as expected. Although there is a decaying trend between the bipartivity and the clustering coefficient, it is clear that even for these small graphs, the contribution of longer cycles to the bipartivity is very important.
In Table 7 we give the values of the bipartivity for the five networks studied in this work. The networks of Stony and powergrid have significant bipartivity, while neurons and yeast are highly non-bipartite. As can be seen in the Table there is not a clear trend between bipartivity and edge density nor to the clustering coefficient of these graphs.  Table 7 Values of the bipartivity, clustering coefficient and edge density of the five real-world networks studied in this paper In the case of Stony we have obtained a bipartition of the network using a technique also based on matrix exponentials. The result is illustrated in Fig. 9 where the edges colored in red or in blue are those that frustrate the bipartition of the network, i.e., those that, if removed, make the graph bipartite.

Gaussian Estrada indices
As we have seen in the previous analysis there are situations in which the Estrada index of a graph is mainly determined by the spectral radius of the adjacency matrix. That is, when λ 1 λ 2 1 the sum j exp λ j is approximated very well by exp (λ 1 ) . From the structural point of view, this means that most of the information contained in the eigenvalues λ j for j > 1 is making almost no contribution to the Estrada index. It is well-known that structural information encoded by some other eigenvalues other than λ 1 is very important for several kinds of problems [46,[50][51][52]213,214]. For instance, the nullity of the graph (see [111] for a review), i.e., the multiplicity of the zero eigenvalue of the adjacency matrix, plays a fundamental role in explaining magnetic properties of materials [230]. In general, many real-world networks have large multiplicity of λ j = 0 (nullity) and of λ j = −1 which points to the fact that some important structural information on these networks is encoded in eigenvalues different from λ 1 .
In this section we investigate Estrada indices that give higher weights to the contribution of eigenvalues other than the spectral radius. In particular we use here a technique known as spectral folding [36,229] to produce Gaussian Estrada indices [5,80]. In the following let λ be a given reference eigenvalue, I (z) be the modified Bessel function of the first kinds, Fig. 9 Illustration of a bipartition of the network of Stony stream using the method developed by [82]. The dotted lines joints the two partitions and continuous lines connect vertices inside the same partition, i.e., they frustrate the bipartition of the network erf (z) be the error function and erfc (z) = 1 − erf (z) be the complimentary error function [5,80].

Definition 20 The Gaussian Estrada index of G is defined as
(9.1) The idea behind this Gaussian Estrada index is explained graphically in Fig. 10. The name Gaussian honors Carl Friedrich Gauss (1777-1855). 6 First we give a few general results for the Gaussian Estrada index (see [5,80]).

Lemma 12 Let G be any graph. Then,
Theorem 29 Let G be a graph with n nodes and m edges. Then,

3)
where k i is the degree of the node i in the graph G .
with equality if and only if G ∼ =Kn.

Random graphs
In this subsection we consider the estimation of the Gaussian Estrada indices of random graphs. The reasons for studying random graphs have been explained in Sect. 5.2. Here we will consider both Erdős-Rényi and Barabási-Albert random graphs.

Double Gaussian Estrada index
Another important situation appearing in many molecular systems is the existence of two reference eigenvalues, typically located around the mid part of the spectrum, which are of great relevance for understanding the behavior of these systems. In 1952, Fukui et al. [99] calculated the chemical reactivity of molecules by using molecular orbital theory, but their method neglects all molecular orbitals except two, the occupied one of higher energy (HOMO) and the vacant one of lowest energy (LUMO). According to Fukui the HOMO gives a molecule a character of electron donor, whereas the LUMO acts as an electron acceptor. The theory was further applied by Woodward and Hoffmann [234] in the interpretation of the stereochemistry of electrocyclic organic reactions. Both, the Frontiers Molecular Orbital (FMO) theory of Fukui and the Woodward-Hoffmann rules are paradigmatic examples of success of theoretical approaches in Chemistry. Both Fukui and Hoffmann won the Nobel Prize in Chemistry for such works. Since then [98], FMO is widely applied for studying chemical reactivity [176]. Let us consider here, for instance, molecular systems S where the energy E is obtained by the time-independent Schrödinger equation: (α I + t N N A) Ψ = EΨ , as described before. Then, when α = 0 and t N N = −1, the energy levels of the system are E j = −λ j . Typically, the states with energy levels E j < 0 are occupied by electrons, while those with energy E j ≥ 0 are empty. Then, the energy level just below E j = 0 is known as the highest occupied molecular orbital (HOMO) and the one just over E j = 0 is the lowest unoccuppied molecular orbital (LUMO). These two molecular orbitals are fundamental in understanding the chemical reactivity of these molecular systems [182]. They can be described in the current approach by the negative of two references eigenvaluesλ 1 andλ 2 of the adjacency matrix. Then, we have the following [6] (Fig. 11). Let K n , K n 1 ,n 2 and K 1,n−1 be the complete, bicomplete and star graphs of n nodes, respectively. Then

Numerical analysis
We consider here the bounds given in Theorem 30 and in Theorem 32 for all connected graphs with 8 nodes. The bound given in Theorem 31 is not applicable in all the cases and we do not considered it for this general case. We show in Fig. 12 the histogram of the relative deviations for these two bounds in these small graphs. The mean relative deviations (in %) of the two bounds are, respectively 89.84 ± 2.51 and 66.51 ± 6.55, which points to the fact that the second bound is a better approximation than the first one to the Gaussian Estrada index. In Table 8 we give the values of the three bounds for the five networks studied here as well as the values of the actual Gaussian Estrada index forλ = 0. The bound given in Theorem 30 is extremely far away from the actual values and practically says the same as the trivial bound G E E 0 (G) > 0. The same happens for Theorem 32 in the cases of Stony and neurons, but it gives more decent estimations for the cases of the bigger networks of Internet and powergrid.

Mittag-Leffler Estrada indices
As we have seen in previous sections of this paper, the Estrada indices of a graph may arise as the solution of the linear dynamical system where M is a given graph matrix, with initial condition u (0) = u 0 . The solution of this system is given by u (t) = exp (t M) u 0 . The case in which M = A is the adjacency matrix of the graph has been analyzed in the paper [172]. Let us consider that, instead of using the first derivative of u (t) respect to time, we use a fractional derivative. Then we have a system of the form: for 0 < α < 1, where the Caputo fractional derivative [37] is given by where f (k) represents the kth derivative of f and Γ (z) is the Euler gamma function: The solution of this system is given by , (10.6) which are the Mittag-Leffler matrix functions (for properties of Mittag-Leffler matrix function the reader is referred to [102,183]).
To catch the analogy with the standard Estrada index of a graph we can write is as due to the fact that Γ (k + 1) = k!, Therefore we can generalize the Estrada index to account for a wider class of penalization functions, such that we write . (10.8) At the same time we keep in mind that E E α,β (G) is the solution of a system of equations of the form D α t u (t) = Au (t) as we will explore later. We propose the name Mittag-Leffler Estrada index for E E α,β (G) in honor to the mathematician Gösta Mittag-Leffler (1846-1927). 7 Some examples of closed formulas are illustrated in Table 9. One important aspect of these functions in general is that by controlling the parameters α and β we can penalize the walks of k length in different ways. For instance, if ((α − 1) k + β) < 0 for all k, then the walks of any length are penalized less than in E E 1,1 (G). This is for instance, the case of E E 1/2,1 (G) (see Table 9). In those cases where ((α − 1) k + β) > 0 for all k, the penalization of all walks is heavier than in the exponential, which are for instance the cases of E E α>1,β (G). There is a third case which occurs when (α − 1) k + β is negative for 0 ≤ k ≤ k c and positive for k > k c , where k c is a given integer. This is the case, for instance, of the matrix functions where k c < − β α − 1 .
Let us first consider the Estrada index E E 1/2,1 (G) = tr exp A 2 (I + erf (A)) . A similar index was defined and studied in the paper [89] in the following form: where k!! is the double factorial of k. The goal of defining such index was to account for less penalization of longer walks, which may play an important role in several applications (for some examples see [1,89] Fig. 13 Examples of three graphs with 8 nodes and chordless cycles of different lengths: G I has a triangle and an heptagon; G I I has a square and a hexagon; G I I I has two pentagons G I , G I I and G I I I , respectively, which represent 19.6% of difference between the first pair and 6.2% between the second one. Finally, we consider the Mittag-Leffler Estrada indices defined as follow: These indices were developed and studied in 2010 in the paper [75], as a way to penalize more heavily the longer walks than in the matrix exponential. The indices E E 1,β+1 (G) are also the trace of the so-called matrix Ψ functions: (10.12) where (10.13) which obey the following recurrence formula: When the adjacency matrix is not singular we can represent these Estrada indices as follow

Resolvent Estrada index
The context of Mittag-Leffler Estrada indices also allow the consideration of other indices that were previously proposed in the literature. This is the case of an index proposed in 2010 in [85]. The goal of introducing this index was to change the penalization of the different powers of the adjacency matrix from k! to (n − 1) k to increase the contribution of walks of longer lengths. The index proposed in [85] corresponds to the trace of the resolvent of the adjacency matrix: which eventually was proposed in [27] to be named as the resolvent Estrada index of the graph. It can be seen that the resolvent Estrada index is a particular case of Mittag-Leffler Estrada index:

Remark 7
The use of the normalization c k = 1/ (n − 1) in (c k A) k is just one of the many possibilities that exist. In reality this normalization is not a good one, because the corresponding Estrada index is very close to the number of nodes of the graph, as can be inferred from the bounds presented before. Then, other general choices of the type ( A) k where < (λ 1 ) −1 are more appropriate here.
A nice result relating the resolvent Estrada index and the characteristic polynomial of the adjacency matrix was proved by in [41] and is given below.

Theorem 35 Let G any graph with n nodes and let P (G, x) be the characteristic polynomial of the adjacency matrix of G, aka its characteristic polynomial. Then,
where P (G, n) is the first derivative of P (G, x) evaluated at x = n.
To illustrate the previous result let us consider the three graphs in Fig. 13. Their characteristic polynomials are, respectively: which give E E 0,1 (G I ) = 4023/484,E E 0,1 (G I I ) = 2980/359, and E E 0,1 (G I I I ) = 2191/264. That is, the difference between the first pair of graphs is only 0.13% and between the second pair is only 0.02%. This is a direct consequence of penalizing more heavily the longer walks than in the exponential matrix function. Some other inequalities have been reported for the resolvent Estrada index in terms of the number of nodes, edges, maximum degree, etc.
Lemma 20 [41] Let G be a simple graph with n nodes and m edges. Then, with equality if and only if G ∼ =Kn.
Lemma 21 [114] Let G be a simple noncomplete graph with n > 3 nodes and m edges. Then, , (10.24) with equality if and only if G ∼ =Kn. Lemma 22 [175] Let G be a simple graph with n nodes, m edges and maximum degree k max = n − 1. Then, . (10.26)

Estrada indices and network of oscillators
The study of vibrations on regular graphs, known as lattices, is standard in solid-state physics (see for instance Chapter 4 in [144]). The techniques of classical as well as of quantum mechanics are used in the analysis of such vibrational problems. In 2003 this analysis was extended to consider non-regular networks [146] where the vibrations where analyzed in the context of a quantum system. Here we investigate the connections existing between some of the Estrada indices and the network vibrations, used in a metaphorical sense. That is, although some physical systems represented by non-regular graphs can be analyzed using the techniques developed here we consider the current approach as an appropriate tool for giving a physical meaning to the indices involved. Let us consider a system S consisting of ball of mass M which are connected by springs with the spring constant Mω 2 . Let us consider that the ball-spring system is submerged into a thermal bath at the temperature τ . Then the balls in the complex network oscillate under thermal disturbances. We will consider that every ball is tied to the ground by a spring which has spring constant satisfying K max v k v (see Fig. 14). This guarantees that the system can oscillates but do not translate from a fixed position. In this way we can analyze how a given ball can transmit small oscillations to the rest of the balls of the system.
The general Hamiltonian of this system is written as where the first term represents the kinetic energy of the corresponding balls and the second term represents the potential energy of the system, with p v being the momentum and x v the coordinate of the ball v.

Quantum oscillators
In this setup we consider that the system obeys the laws of quantum mechanics. Then, the momenta p w and the coordinates x v are not independent variables but they are operators that satisfy the commutation relation: [x v , p w ] = i δ vw ,where i = √ −1, is the Dirac constant and δ vw is Dirac delta. Additionally we will use second quantization to express the creation and annihilation of oscillations at the given balls of the system. That is, we use the boson creation and annihilation operators defined by [184] With the use of these operators, the Hamiltonian of a network of quantum harmonic oscillators is given by [84] (11.4) where Ω = √ K /Mω and K is a constant such that K max v k v . Since A is symmetric, we can diagonalize it by means of an orthogonal matrix O as in (11.5) where Λ is the diagonal matrix with the eigenvalues λ μ of (K I − A) on the diagonal. This generates a new set of boson creation and annihilation operators as We can then decouple the Hamiltonian aŝ (11.9) We now introduce an approximation in which each mode of oscillation does not get excited beyond the first excited state. In other words, we restrict ourselves to the space spanned by the ground state (the vacuum) |vac and the first excited states b † μ |vac . Then the second term in the last line of the Hamiltonian (11.9) equals zero and we have (11.10)

Remark 8
This approximation is justified when the energy level spacing Ω is much greater than the energy scale of external disturbances, (specifically the temperature fluctuation k B T = 1/τ ), as well as than the energy of the network springs ω, i.e. τ Ω >> 1 and Ω >> ω. This happens when the mass of each oscillator is small, when the springs to the ground, MΩ 2 , are strong, and when the network springs Mω 2 are weak. Then an oscillation of tiny amplitude propagates over the network.
We are going to work in this limit hereafter. The thermal bath represents here an 'external situation' which affects all the links in the network at the same time, e.g., economic crisis, social agitation, extreme physiological conditions, etc. After equilibration, all links in the network are weighted by the parameter τ = (k B T ) −1 .
Let us now compute how much an excitation at the node p propagates throughout the network before coming back to the same node and being annihilated. This information is obtained through the diagonal thermal Green's function, which is given in the framework of quantum mechanics by (11.11) where the partition function is given by The diagonal thermal Green's function can then be obtained as [84] (11.13) where we have used the spectral decomposition of A in the last line.

Remark 9
In [84] it is remarked that (S, τ ) , (11.15) which indicates that the Estrada index represents the sum of the excitations started at every node of a graph, which propagate throughout the network before coming back to the same node and being annihilated in a network of quantum harmonic oscillators.

Classical oscillators
Here we consider a system S like the one described before but obeying the laws of classical mechanics. In this case we can write the Hamiltonian of the system by considering only the potential energy (for justification see [84]): (11.16) where x = (x 1 , x 2 , . . . , x n ) T and I is the n × n identity matrix.
We can now diagonalize A as before and by taking a sufficiently large value of the constant K , we can make all eigenvalues λ μ positive. By defining a new set of variables y μ as y = Ox and x = O T y, we can transform the Hamiltonian in the form Here again we focus of the quantification of those oscillations that start at a given ball of the system, navigates the whole system and return to the corresponding ball. Namely, (11.18) where the partition function is given by (11.19) where the integral is n-fold. Now, because the Jacobian of the orthogonal matrix O is unity we have i dx i = μ dy μ . Therefore, the multi-fold integration in the partition function is decoupled to give which can be written in matrix form as Since we have made all the eigenvalues of (K I − A) positive, its determinant is positive. Similarly, we have μ dy μ . (11.22) In the integrand, odd functions with respect to y μ vanish. Therefore, only the terms of y 2 σ survive after integration in the expansion of the square parentheses in the integrand. This gives (11.23) Comparing this expression with that of the partition function we have [84] (11.24)

Remark 10
In [84] it is remarked that if K = n − 1 then we have (11.25) which indicates that the resolvent Estrada index represents the sum of the excitations started at every node of a graph, which propagate throughout the network before coming back to the same node in a network of classical harmonic oscillators.

Estrada indices and epidemics on networks
In continuation with the previous line of research in which the Estrada index is derived from a given dynamical systems context we analyze here its connection with epidemiological models on networks. The field of mathematical epidemiology has a long tradition in applied mathematics (see for instance [7,34,169]. In 2001, the authors of the seminal work [179] discovered the tremendous influence of network topology on epidemic spreadings. Since then, the use of network-theoretic approaches together with epidemiological models have become a necessary combination [143,149]. Here we will show that such networked epidemiological models have a clear connection with the Estrada index of a graph. For that we will briefly introduce the Susceptible/Infected (SI) model on networks. The reader should be aware that this is a generalist model that can be used in many different scenarios, not only on the analysis of diseases propagating on a network. Let G be a graphs whose nodes can be in either of two states: either susceptible or infected one. An infected node can transmit the infection to any other node in the graph to which it interacts with. Then, if ζ is the rate at which such infection is transmitted between nodes, and if s v (t) and x v (t) are the probabilities that the node v is susceptible or get perturbed at time t, respectively, we can write the dynamics [170]: where A vw are the entries of the adjacency matrix of the graph for the pair of nodes v and w, and N is the set of nearest neighbors of v. In matrix-vector form becomes [170]: with initial condition x (0) = x 0 . The SI model can be rewritten as which is equivalent to Lee et al. [155] have considered the following linearized version of the previous nonlinear equation . They have found that the solution to this linearized model is [155]: When t 0 = 0, x i (0) = c/N , i = 1, 2, . . . , N for some c, the previous equation is transformed toŷ where γ = 1 − c/N and 1 is the all-ones vector. In [155] the authors proved that this solution is an upper bound to the exact solution of the SI model. Therefore, if we take the sum of the entries ofŷ (t) at a given t we have where C 1 = (1/γ − 1) and C 2 = (1/γ − 1 + log (γ )) . Obviously the first term in the square bracket is the Estrada index of the graph in which edges are weighted by γ ζ t. This term represents the circulability of the infection around the nodes of the graph, while the term tr (J − I ) e γ ζ t A , where J is the all-ones matrix, accounts for the transmissibility of the disease between the nodes.

Fractional SI model on networks
In recent years there have been an explosion of works in which the classical derivatives used in the epidemiological models have been substituted by fractional ones [8,11,12,133,192]. There have been several reasons for adopting fractional epidemiological models. They include for instance 1. the fact that the fractional parameter can be tuned to adjust the output of the model to real data [8] and so they can be more accurate that models using standard derivatives; 2. the fact that a fractional differential operators may be derived from epidemiological models whenever the recovery time from the disease is power-law distributed [11]; 3. the fact that fractional derivatives capture the history of the variable, that is, they have memory, and the effect of recent memory is more important than the effect of older memory [12,192].
In general, fractional derivatives are nowadays widely used to model biological processes [136] to incorporate different aspects of the dynamics in such systems. Here, we will describe a model which naturally gives rise to the Mittag-Leffler Estrada index in the context of epidemiological models. We proceed by considering fractional time-derivatives in the modified SI model proposed in [155]. That is, in [1] the authors considered the following linearized fractional SI equation wherex (t) = f ŷ (t) in whichx (t) is an approximate solution to the fractional SI model,ŷ is the solution of (12.9) with initial conditionŷ (0) = g (x (0)) and b (x) Here D α t f (t) is the fractional time derivative in the Caputo formulation [37], which was previously given in Eq. (10.3).
For convenience, we write Ω := diag (1 − x 0 ) , andÂ = AΩ. It was then proved that this solution is an upper bound to the exact fractional SI model.
Let us fix the following notation. Let x and y be two vectors of the same length n. Then, we say that x y is x i ≤ y i for all i = 1, . . . n. Letx(t) be the solution of the linearized fractional models of the form: D α tx (t) = ζ Ax(t), which is exponential unstable.
Therefore, here again if we take the sum of the entries ofŷ (t) at a given t we have Thus, again the Mittag-Leffler Estrada index, which is the first term in the squared bracket, represents the circulability of the infection around the nodes of the graph in the fractional SI model, while the second term represents the transmissibility of the disease between the nodes.

Estrada indices from piecewise walk penalization
In the same work [75] in which the author proposed the use of the matrix Ψ functions as a way to increase the penalization of longer walks in graphs, a different strategy was proposed to drop such penalization relative to the exponential matrix function. This strategy can be formulated as a piecewise penalization as follows. Suppose that we do not want to penalize the walks of lengths smaller than certain value t ∈ Z. Then, we define the following stepwise function: such that the piecewise Estrada index of the graph G is defined as: In the case that the adjacency matrix has no unity eigenvalue we can write this Estrada index as [75]: For computational purposes this expression can be adapted for any network as follows. Let r be a constant sufficiently close to one, such that r = 1/λ for all λ, which are the eigenvalues of A. Then,

Nonlocal adjacency, Harary Estrada index and beyond
There are physical situations in which the entities of a system not only interact if they are nearest neighbors, but also through nonlocal interactions. These long-range interactions have been documented in physical, chemical and biological systems [3,39,148,164,177,194,195,235,240]. In a physical context, like the tight-binding kind of models described before, these nonlocal interactions corresponds to the case where the Hamiltonian of the system describes not only NN interactions but also next-nearest-neighbor (NNN) and other interactions beyond them [31,168,174,193]:Ĥ In this framework we have that the system can be described by the weighted sum of higher-order adjacency matrices: where A 2 is a matrix with entries (A 2 ) i j equal to one if i and j are not adjacent and are separated by two edges or zero otherwise. We can extend this concept to any other separation, such that [77] where d i j is the length of the shortest path between the two nodes. The parameters t N N , t N N N , etc. are expected to decay with the length of the separation between the corresponding entities. That is, the strength of the interaction decays with a given law of their separation d, i.e., f (d). In this way we can write [72,218] where diam is the diameter of the graph. Let us see how we can constructĤ.
Here we will use a min-plus algebra to define what otherwise is the shortest path distance matrix of the graph. We do that because it is a mathematically elegant approach, which may also open some possibilities for studying other kinds of functions for graphs.
Let (R ∪ {+∞} , ⊕, ⊗) be the min tropical semiring with the operations [32,131,139]: x ⊕ y := min {x, y} , x ⊗ y := x + y. (14.5) The identity element for ⊕ is +∞ and that for ⊗ is 0. Then, we can define the tropical adjacency matrix power as where A ⊗0 =Î , which is the tropical identity matrix, i.e., a matrix with zeros in the main diagonal and ∞ outside it. Here, the function f could be an exponential, a trigonometric function or simply the power function. Let us hereafter focus only on the negative power function, such that (−s) represents the entrywise power. We can now write: . (14.8) The tropical sum is carried out up to infinity as it converges in all cases where there are no negative cycles in the graph. A negative cycle is a cycle where the product of the weights of all its edges is negative. Typically, except for signed graphs, we consider positive edge weights, which always avoid such negative cycles. The infinite sum ∞ k=0 A ⊗k is known as the Kleene star operator of A [32,131,139]. Obviously, A s = A i j (s) , where [72,218] (14.9) which are the entry-wise powers of all nondiagonal entries of the shortest path distance matrix of the graph. The parameter s accounts for the strength of the nonlocal interaction. Notice that Here again, in the statistical physics context, the partition function of the system containing nonlocal interactions is: where τ is the inverse temperature as before. Because the parameters t N N , t N N N , etc., are negative we have that Z = tr e τ A s =: E E (A s , β) , (14.11) Then, using the same definitions as the ones given before we can define the entropy, enthalpy and free energy of the system having local and nonlocal interactions. It is important to notice that When s = 1, the corresponding matrix A s=1 is known in mathematical chemistry as the Harary matrix [137,167,181] in honor to mathematician Frank Harary . 8 Nowadays there are not many results about the HEE index. Hereafter we collect some of the existing ones for simple graphs [109], H E E (G) = E E (A s=1 , τ = 1) .

Theorem 37 [109] Let G be a simple graph with n vertices and m edges. Then, the Harary Estrada index of G is bounded as
with equalities attained if and only if G ∼ =K1.
Theorem 38 [138] Let G be a simple graph with n ≥ 2 vertices and let κ = 1 2 tr A 2 s=1 . Then, the Harary Estrada index is bounded as .

Numerical analysis
In Fig. 15  indices for these small graphs, which is an extremely poor performance of this bound.
In Table 10 we give the values of the bounds previously considered for the five real-world networks analyzed here as well as the actual values of H E E (G). As can be seen both lower and upper bounds are extremely far from the actual values of the Harary Estrada indices of these real-world networks. In particular, the upper bound is extremely higher than the actual values.

Laplacian Estrada index and backward diffusion
In the study of graph properties, the function K (G) = e −t L , where L is the graph Laplacian, has found many applications [22,45,150,236]. The name "Laplacian" honors  mathematician Pierre-Simon Laplace (1749-1827). 9 The function K (G) is known as the heat kernel of the graph [125,191] and appears naturally in the solution u (t) = exp (−tD L) u 0 of the diffusion equation on graphs: where D is the diffusivity (see Sect. 16). Therefore, the trace of the heat kernel would correspond to a sort of diffusion Estrada index. However, in [92] the following index was proposed and named "the Laplacian Estrada index" of the graph where μ j are the corresponding eigenvalues of L. Therefore, what the authors of [92] have proposed can be though as an index related to the solution of the backward diffusion equation, i.e., negative time, or as a diffusion equation with negative diffusivity D < 0 [65]. There are physical situations in which such negative diffusivity appears [48,140,173,225,226]. For instance, in the simultaneous diffusion of boron and point defect in silicon, the diffusivities of interstitial could be negative [225]. That is, the diffusion process of interstitial or vacancy could be a backward diffusion in silicon. In other scenarios, a backward diffusive model is used to detect the potential location of sources in spreading processes.
In [92] the authors proved the following result.
Proposition 2 Let G be a simple graph with n nodes and m edges. Let Z = i k 2 i be the first Zagreb index of G. Then, with equality if and only if G ∼ =Kn.
Further, in [249] the authors proved the following results.

Proposition 3
Let G be a simple graph with n nodes and m edges. Let Z = i k 2 i be the first Zagreb index of G. Then,

4)
with equality if and only if G ∼ = K 2 ∪K n−2 or G ∼ =Kn.

Proposition 4 Let G be a simple graph with n nodes and m edges. Then,
L E E (G) ≥ 2 + n (n − 1) e 4m + 4 − 3n − 4m, (15.5) with equality if and only if G ∼ =Kn.
Other bounds were obtained in [248] on the basis of the degree sequence of a graph.
with equality if and only if G is a vertex disjoint union of complete subgraphs.
Several bounds have been proposed on the basis of the maximum and minimum degrees of a graph. We resume some of them here.
Theorem 39 [160] Let G be a simple graph with n nodes and m edges. Let k max and k min be the maximum and minimum node degrees of G. Then, L E E (G) ≥ e k max +1−2m/n + (n − 2) e 4m/n−k max −1 1/(n−2) + e −2m/n , (15.8) with equality if and only if G ∼ = K n or G ∼ = S n .
Theorem 40 [40] Let G be a simple graph with n nodes and m edges. Let k max and k min be the maximum and minimum node degrees of G. Then, with equality if and only if G ∼ = 2K 1 ∨ K n−2 or G ∼ = K 1,n−1 or G ∼ = K (n−1)/2 ∪ K (n−1)/2 (n is odd).
Theorem 41 [161] Let G be a simple graph with n nodes and m edges. Let k max and k min be the maximum and minimum node degrees of G. Then, with equality if and only if G ∼ =Kn.
The following is an upper bound found in [163].

Theorem 42
Let G be a simple graph with n nodes and m edges. Then, (15.11) with equality if and only if G ∼ = K n or G ∼ = K n − e.
Finally, we present the estimation made in [100] for the Laplacian Estrada index of Erdős-Rényi random graphs.
Theorem 43 Let G n, p be an Erdős-Rényi random graph with n nodes and probability p. Then, the Laplacian Estrada index is given by L E E G n, p = e np (n − 1) e o(1)n + o (1) , a.s. (15.12) In [132] the authors find estimations for the Laplacian Estrada index of random multipartite graphs.

Remark 11
Other bounds and estimations have been reported for the Laplacian Estrada index of specific graphs, or based on other graph parameters not considered here. Some nonexhaustive examples are: [20,58,60,61,117,135,145,162,207,242,243,251].

Remark 12
The normalized Laplacian Estrada index defined as where K is the diagonal matrix of node degree has been studied in [47,123,158,203,207].

Remark 13
The signless Laplacian Estrada index defined as has been also studied in [16,117,227].

Numerical analysis
In Fig. 16 we illustrate the histograms of the relative deviations (in %) of the lower bounds (Proposition 2), (Proposition 4) (Proposition 5), (Proposition 6), (Theorem 40) and (Theorem 41) for all connected graphs with 8 nodes. The best performance is obtained from the bound (Proposition 6) followed by (Proposition 5). We also analyzed the upper bounds given in Proposition 2, Proposition 3, Theorem 41 and Theorem 42 for the same set of graphs. In these cases the best performances were obtained for Theorem 41 and Theorem 42, while 2 give very high upper bounds (Fig. 17).
In Table 11 we give the lower bounds for the Laplacian Estrada index of five real-world networks. In general, the bounds (Proposition 5), (Proposition 6) and (Theorem 40) perform very well, while (Proposition 2), (Proposition 4) and (Proposition 41) are several orders of magnitude below the actual values of the Laplacian Estrada indices of these networks.
The case of the upper bound is much more contrasting with values several orders of magnitude over the actual values of the Laplacian Estrada indices of these five networks. We have used variable-precision floating-point arithmetic" (VPA) to evaluate each element of the symbolic input in Matlab for these calculations. It is used to evaluate symbolic inputs with variable-precision floating-point arithmetic, calculating values to 32 significant digits.   The results are given in Table 12.
The following bounds based on the largest Laplacian eigenvalue μ 1 perform very well for the four real-world networks analyzed as can be seen in Table 13. The reason is that the largest eigenvalue of the Laplacian matrix dominates the spectrum of this matrix, i.e., it is very large and separated from the second largest eigenvalue.

Radius of gyration and distance Estrada index
When presenting the diffusion equation on graphs, Eq. (15.1), we mentioned in passing the diffusion coefficient D, which appears in the equation and in its solution. The diffusion coefficient is related to the radius r of the spherical particle diffusing on a medium of viscosity where τ is the inverse temperature. In the case of small molecules like drugs, or macromolecular systems like proteins, the particles cannot longer be considered spherical. In these cases it is customary to replace the radius of the spherical particle by the radius of gyration of the corresponding molecule [94,107,126,152,171]. The radius of gyration is defined as follows. Let S = ( p 1 , . . . , p n ) be a system formed by n particles or points p i , which are located in a given region of the three-dimensional Euclidean space. Let r i j be the Euclidean distance between the particles p i and p j . Then, the radius of gyration of S is defined as [94]. However, it has been shown that even when the radius of gyration based on Euclidean distances is used, there are cases of undesired degeneration of the index for pairs of clusters [74]. That is, there are pairs of nonisomorphic clusters which have the same radius of gyration. Some examples in 2-and in 3-dimensions are given in Fig. 18. The radius of gyration is widely used in organic chemistry, polymer sciences, proteins, and RNA, in general, for the study of their compactness. Most of these molecular systems can be represented as graphs. For instance, molecules are typically represented by molecular graphs [217], proteins can be represented by "protein residue networks" [76,212], and the secondary structure of RNA is also represented by graphs [147]. Then, it is important to extend the concept of "radius of gyration" to graphs.

Definition 22
Let G be a simple graph. Let d i j be the shortest path distance between the nodes i and j. The graph radius of gyration is defined as Let D be the shortest path distance matrix of G. Then, it is straightforward to realize that Then, if σ 1 ≥ σ 2 ≥ · · · ≥ σ n are the eigenvalues of D, the graph radius of gyration is the second spectral moment of D, i.e., Therefore, we can say that the second moment of the shortest path distance matrix is a measure of the packing of the graph. In the case of the graph we can have a similar degeneracy of the index R 2 G for nonisomorphic graphs. For instance in Fig. 19 we give an example of four nonisomorphic graphs with the same value of R 2 G = 9 50 . In order to ameliorate this degeneracy problem we can think on extending the packing measure to higher moments of D as P (G) = c 2 tr D 2 + c 3 tr D 3 + · · · . (16.5) Let us include the term c 0 tr D 0 + c 1 tr D 1 (the first is the weighted number of vertices and the second is zero in simple graphs) and let us consider c k = (k!) −1 . Then, we get the following index, first defined in [108].

Definition 23
Let G be a simple graph. Let D be the shortest path distance matrix of G with eigenvalues σ 1 ≥ σ 2 ≥ · · · ≥ σ n . Then the distance Estrada index of the graph is which is an index of the packing of the graph.
For instance, for the graphs in Fig. 19 Fig. 19, (d) is the most "packed" one in terms of its shortest path distance and the one in (a) is the least packed. Several bound have been obtained for the distance Estrada index, some of which are resumed below.
Theorem 44 [108] Let G be a simple graph with n nodes and m edges. Then, if the diameter is d max , the distance Estrada index is bounded as Theorem 47 [204] Let G be a simple graph with n nodes and m edges, the distance Estrada index is bounded as DE E (G) ≥ e 2(n−1)−2m/n + e −(2(n−1)−2m/n) + n − 2 (16.11) where equality is attained if and only if G ∼ = K 2 .

Random graphs
The distance Estrada index has been studied for random graphs where some bounds have been reported in [205,209].

Theorem 48
Let G n, p be an Erdős-Rényi random graph with n nodes and probability p. Then, the distance Estrada index is bounded as (16.14)

Numerical analysis
We analyze here the lower bounds in Theorems 44, 45, 46 and 47. The relative deviations (in %) are illustrated in Fig. 20. As can be seen the closest values are obtained by the bound given in Theorem 45. The distributions of the relative deviations for Theorems 46 and 47 appears to show some dependencies with the structure of the graphs, which produce the multi-peak structures observed in the histograms. We also considered the upper bounds given in Theorems 44, 45 and 46 were we observe that these bounds are several order of magnitude over the actual values of DE E (G) even for small graphs like the ones studied here (Fig. 21).
In Table 14 we give the values of the lower and upper bounds as well as the actual values calculated with very precise arithmetic (vpa) in Matlab for the five real-world networks studied. As can be seen in the Table 14 the bound given in Theorem 45 gives the best lower and upper estimates of the distance Estrada index. It is also interesting to remark that the network of the western USA power grid displays an extremely large value of DE E (G), indicating that it is a very poorly packed network. Indeed, this network is planar as the power stations are embedded in the landscape of the western USA The results obtained for Erdős-Rényi random graphs G n, p with 1000 ≤ n ≤ 4000 and p = 0.5 are illustrated in Table 15, showing good agreement between the actual values and those predicted by Theorem 48. The values were computed in [209] using variable-precision floating-point arithmetic (VPA) in Matlab.

Conclusions
We presented an account of the many different facets of the Estrada indices of graphs. Starting from the "classical" Estrada index we give several interpretations of the index based   on (i) combinatorics of subgraphs, (ii) statistical mechanics, (iii) marginal probability in a quantum system, (iv) oscillations models on networks, and (v) epidemiological models on networks. Then we move forward to the analysis of other kinds of Estrada indices. First we contextualize these indices originally introduced in an ad hoc way in the mathematical literature. For instance, the Seidel Estrada index is placed in the context of signed graphs, the theory of balance and the concept of network bipartivity. The resolvent Estrada index is analyzed as a case of Mittag-Leffler Estrada indices which appear in the context of fractional epidemiological models on graphs. The Harary Estrada index is understood as a particular case of nonlocal operator on graphs. The Laplacian Estrada index is now pondered on the basis of the diffusion equation with negative diffusivity or a backward diffusive process. Finally, the distance Estrada index is considered in the context of the radius of gyration of a graph, which can be connected to the diffusion coefficient of graphs via the Stokes-Einstein equation. In all cases we have provided numerical analysis of several of the bounds and estimations made for these indices. Such results have revealed the necessity of investigating more robust bounds, particularly upper bounds, for most of the indices studied. In many cases the bounds, although correct, are very far away from the actual values of the indices, which leaves large rooms for improvements. We encourage authors searching for new bounds to compare them with the existing ones with the challenge of improving them for general classes of graphs. Finally, we have not considered many of the results obtained in the literature for specific classes of graphs, which would make this paper too long to be digested. We advice the reader that such bounds exist for several of the indices described in this paper and for several classes of graphs of importance in specific areas of applications.
Acknowledgements The author thanks financial support from Ministerio de Ciencia, Innovacion y Universidades, Spain for the grant PID2019-107603GB-I00 "Hubs-repelling/attracting Laplacian operators and related dynamics on graphs/networks". The author thanks the referees for exhaustive and construtive revision of the manuscript.
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. number of subgraphs illustrated in Fig. 1 are obtained as follow: 2) 3)