Asymptotic Spectra of Large (Grid) Graphs with a Uniform Local Structure (Part I): Theory

We are mainly concerned with sequences of graphs having a grid geometry, with a uniform local structure in a bounded domain Ω⊂Rd,d≥1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\Omega} {\subset} \mathbb{R}^{d}, d \geq 1$$\end{document}. When Ω=[0,1]\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega = [0, 1]$$\end{document}, such graphs include the standard Toeplitz graphs and, for Ω=[0,1]d\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega = [0, 1]^{d}$$\end{document}, the considered class includes d-level Toeplitz graphs. In the general case, the underlying sequence of adjacency matrices has a canonical eigenvalue distribution, in the Weyl sense, and we show that we can associate to it a symbol f\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathfrak{f}$$\end{document}. The knowledge of the symbol and of its basic analytical features provides many information on the eigenvalue structure, of localization, spectral gap, clustering, and distribution type. Few generalizations are also considered in connection with the notion of generalized locally Toeplitz sequences and applications are discussed, stemming e.g. from the approximation of differential operators via numerical schemes. Nevertheless, more applications can be taken into account, since the results presented here can be applied as well to study the spectral properties of adjacency matrices and Laplacian operators of general large graphs and networks.


Introduction
In this work we are interested in defining and studying a large class of graphs enjoying few structural properties: a) when we look at them from "far away", they should reconstruct approximately a given domain Ω ⊂ [0, 1] d , d ≥ 1, i.e., the larger is the number of the nodes the more accurate is the reconstruction of Ω;

Background and definitions
In this section we present some definitions, notations, and (spectral) properties associated with graphs (see [7] and references therein) and, in particular, with Toeplitz graphs [18].
between the nodes and a potential term κ : V → R. The non-zero values w(v i , v j ) of the weight function w are called weights associated to the edge (v i , v j ). Given an edge e = (v i , v j ) ∈ E, the nodes v i , v j are called end-nodes for the edge e. An edge e ∈ E is said to be incident to a node v i ∈ V if there exists a node v j = v i such that either e = (v i , v j ) or e = (v j , v i ). A walk of length k in G is a set of nodes v i1 , v i2 , . . . , v i k , v i k+1 such that for all 1 ≤ r ≤ k, (v ir , v ir+1 ) ∈ E. A closed walk is a walk for which v i1 = v i k+1 . A path is a walk with no repeated nodes. A graph is connected if there is a walk connecting every pair of nodes. A graph is said to be unweighted if w(v i , v j ) ∈ {0, 1} for every v i , v j ∈ V . In that case the weight function w is uniquely determined by the edges which belong to E.
A graph is said to be undirected if the weight function w is symmetric, i.e., for every edge In this case the edges (v i , v j ) and (v j , v i ) are considered equivalent and the edges are formed by unordered pairs of vertices. Two nodes v i , v j of an undirected graph are said to be neighbors if (v i , v j ) ∈ E and we will write v i ∼ v j . On the contrary, if (v i , v j ) / ∈ E then we will write v i ≁ v j .
An undirected graph with unweighted edges and no self-loops (edges from a node to itself ) is said to be simple.
Throughout this work, we will always consider undirected graphs without self-loops. Since the potential term κ will come into play only in Section 6, we will often omit it in the next subsections and in Section 4. Moreover, when dealing with simple graphs we will use the simplified notation G = (V, E). Definition 2.3 (Sub-graph, interior and boundary nodes) Given a graphḠ = (V ,Ē,w,κ) we say that a graph G = (V, E, w, κ) is a (proper) sub-graph ofḠ, and we write G ⊂Ḡ, if (iii) w =w |E ; (iv) κ(v i ) =κ(v i ) for every node v i such that there not existsv j ∈V \ V ,v j ∼ v i . Sometimes we will callḠ the mother graph. The set of nodes is called interior of V and its elements are called interior nodes. Vice-versa, the set of nodes ∂V := v i ∈ V | v i ∼v j for somev j ∈V \ V is called boundary of V and its elements are called boundary nodes. Therefore, condition ((iv)) can be restated saying that κ =κ |V . Observe that we do not request that κ =κ on the boundary of V . Definition 2.4 (Degree of a node) In an undirected graph, the degree of a node v i ∈ V , denoted by deg(v i ), is the sum of weights associated to the edges incident to v i , that is, Definition 2.5 (Adjacency matrix) Every graph G = (V, E, w, κ) with κ ≡ 0 can be represented as a matrix W = (w i,j ) n i,j=1 ∈ R n×n , called the adjacency matrix of the graph. In particular, there is a bijection between the set of weight functions w : V × V → R and the set of a adjacency matrices W ∈ R |V |×|V | .
The entries of the adjacency matrix W are In short, the adjacency matrix tells which nodes are connected and the 'weight'of the connection. If the graph does not admit self-loops, then the diagonal elements of the adjacency matrix are all equal to zero. In the particular case of an undirected graph, the associated adjacency matrix is symmetric, and thus its eigenvalues are real [3].
Definition 2.6 (Isomorphism between graphs) Given two graphs G = (V, E, w, κ), . . , v ′ m }, we say that G is isomorph to G ′ , and we write G ≃ G ′ , if i) n = m, i.e., |V | = |V ′ | where | · | is the cardinality of a set; ii) there exists a permutation P over the standard set [n] such that w(v i , v j ) = w ′ (v ′ P (i) , v ′ P (j) ), κ(v i ) = κ ′ (v ′ p(i) ). In short, two graphs are isomorphic if they contain the same number of vertices connected in the same way. Notice that an isomorphism between graphs is characterized by the permutation matrix P .
As an immediate consequence of the previous definition, it holds that G ≃ G ′ if and only if there exists a permutation matrix P such that W = P W ′ P −1 = P W ′ P T , where W, W ′ are the adjacency matrices of G and G ′ , respectively. Definition 2.7 (Linking-graph operator) Given ν ∈ N, we will call linking-graph operator for the reference node set [ν] any non-zero R ν×ν matrix, and we will indicate it with L. Namely, a linking-graph operator is the adjacency matrix for a (possibly not undirected) graph G = ([ν], E, l), with l a weight function. Trivially, note that L may have nonzero elements on the main diagonal, so it admits loops. When the entries of L are just in {0, 1} then we call it a simple linking-graph operator.
In Section 4, we will use the linking-graph operator to connect a (infinite) sequence of graphs G 1 ≃ G 2 ≃ . . . ≃ G n , and that will grant a uniform local structure on the graph G := ∞ n=1 G n . The set of real functions on V will be denoted as C(V ). Trivially, C(V ) is isomorph to R n . Of great importance for Section 6 will be the operator ∆ G : C(V ) → C(V ) defined below.
Definition 2.8 (Graph-Laplacian) Given an undirected graph with no loops G = (V, E, w, κ), the graph-Laplacian is the symmetric matrix ∆ G : C(V ) → C(V ) defined as where D is the degree matrix and K is the potential term matrix, that is, and W is the adjacency matrix of the graph G, that is, Namely,

Toeplitz matrices, d-level Toeplitz matrices, and symbol
Toeplitz matrices A n are characterized by the fact that all their diagonals parallel to the main diagonal have constant values a ij = a i−j , where i, j = 1, . . . , n, for given coefficients a k , k = 1 − n, . . . , n − 1: When every term a k is a matrix of fixed size ν, i.e., a k ∈ C ν×ν , the matrix A n is of block Toeplitz type. Owing to its intrinsic recursive nature, the definition of d-level (block) Toeplitz matrices is definitely more involved. More precisely, a d-level Toeplitz matrix is a Toeplitz matrix where each coefficient a k denotes a (d − 1)-level Toeplitz matrix and so on in a recursive manner. In a more formal detailed way, using a standard multi-index notation (see [33] and Remark 1 at the end of this section), a d-level Toeplitz matrix is of the form , with the multi-index n such that 0 < n = (n 1 , . . . , n d ) and a k ∈ C, −(n − 1) k n − 1. If the basic elements a k denote blocks of a fixed size ν ≥ 2, i.e. a k ∈ C ν×ν , then A n,ν is a d-level block Toeplitz matrix, For the sake of simplicity, we write down an example explicitly with d = 2 and ν = 3: Observe that each block A k1 has a (block) Toeplitz structure. When ν = 1 then we will just write A n,ν = A n .
Here we are interested in asymptotic results and thus it is important to a have a meaningful way for defining sequences of Toeplitz matrices, enjoying global common properties. A classical and successful possibility is given by the use of a fixed function, called the generating function, and by taking its Fourier coefficients as entries of all the matrices in the sequence.
More specifically, given a function f : [−π, π] d → C ν×ν belonging to L 1 ([−π, π] d ), we denote its Fourier coefficients bŷ (the integrals are done component-wise), and we associate to f the family of d-level block Toeplitz matrices We call {T n,ν (f )} n the family of multilevel block Toeplitz matrices associated with the function f , which is called the generating function of {T n,ν (f )} n . If f is Hermitian-valued, i.e. f (θ) is Hermitian for almost every θ, then it is plain to see that all the matrices T n,ν (f ) are Hermitian, simply because the Hermitian character of the generating function and relations (2) imply thatf −k =f H k for all k ∈ Z d . If, in addition, f (θ) = f (|θ|) for every θ, then all the matrices T n (f ) are real symmetric with real symmetric blockŝ Remark 1 (Multi-index notation) Given an integer d ≥ 1, a d-index k is an element of Z d , that is, k = (k 1 , . . . , k d ) with k r ∈ Z for every r = 1, . . . , d. Through this paper we will intend Z equipped with the lexicographic ordering, that is, given two d-indices i = (i 1 , . . . , i d ), j = (j 1 , . . . , j d ), then we write i ⊳ j if i r < j r for the first r = 1, 2, . . . , d such that i r = j r . The relations , ⊲, are defined accordingly.
Given two d-indices i, j, we write i < j if i r < j r for every r = 1, . . . , d. The relations ≤, >, ≥ are defined accordingly.

Weyl eigenvalue distribution and clustering
We say that a matrix-valued function g : D → C ν×ν , ν ≥ 1, defined on a measurable set D ⊆ R m , is measurable (resp. continuous, in L p (D)) if its components g ij : D → C, i, j = 1, . . . , ν, are measurable (resp. continuous, in L p (D)). Let µ m be the Lebesgue measure on R m and let C c (C) be the set of continuous functions with bounded support defined over C. Setting d n the dimension of the square matrix X n,ν , for F ∈ C c (C) we define Σ λ (F, X n,ν ) := 1 dn dn k=1 F (λ k (X n,ν )). Hereafter, with the symbol {X n,ν } we will indicate a sequence of square matrices with increasing dimensions, i.e., such that d n → ∞ as n → ∞, with ν fixed and independent of n. Definition 2.9 (Eigenvalue distribution of a sequence of matrices) Let {X n,ν } n be a sequence of matrices and let g : D → C ν×ν be a measurable Hermitian matrix-valued function defined on the measurable set D ⊂ R m , with 0 < µ m (D) < ∞.
We say that {X n,ν } n is distributed like g in the sense of the eigenvalues, in symbols {X n,ν } n ∼ λ g, if where λ 1 (g(y)), . . . , λ ν (g(y)) are the eigenvalues of g(y). Let us notice that in the case ν = 1, then identity (4) reduces to We will call g the (spectral) symbol of {X n,ν } n .
The following result on Toeplitz matrix sequences linking the definition of symbol function and generating function is due to P. Tilli.
that is the generating function of {T n,ν (f )} n coincides with its symbol according to Definition 2.9.
Since in this paper we work only with undirected (i.e., symmetric) graphs, and in virtue of the applications in Section 6, we will deal with symbol functions g such that λ k (g(y)) are real valued for every y ∈ D, and for every k = 1, . . . , ν. See for example Propositions 3.2,3.4 and Theorem 4.2.
The knowledge of the symbol function g can give valuable insights on the distribution of the eigenvalues of a sequence of matrices. We refer to the Appendix A where a collection of theoretical results are presented, and to Section 6 where numerical experiments are provided.

Diamond Toeplitz graphs
In this section we are going to present the main (local) graph-structure which will be used to build more general graphs as union of sequences of sub-graphs, i.e., diamond Toeplitz graphs. The resulting graphs will be then immersed in bounded regular domains of R d in Section 4. We will proceed step by step, gradually increasing the complexity of the graph structure.
As a matter of reference, we have the following scheme of inclusions, with the related variable coefficient versions:

Toeplitz graphs and d-level Toeplitz graphs
We first focus on a particular type of graphs, namely Toeplitz graphs. These are graphs whose adjacency matrices are Toeplitz matrices. Definition 3.1 (Toeplitz graph) Let n, m, t 1 , . . . , t m be positive integers such that 0 < t 1 < t 2 < . . . < t m < n, and fix m nonzero real numbers w t1 , . . . , w tm . A Toeplitz graph, denoted by T n (t 1 , w t1 ), . . . , (t m , w tm ) , is an undirected graph defined by a node set V n = {v 1 , . . . , v n } and a weight function w such that In the case of simple graphs, i.e., w t k = 1 for every k, then we will indicate the Toeplitz graph just as T n t 1 , . . . , t m . The number of edges in a Toeplitz graph is equal to m k=1 (n − t k ). By construction, the adjacency matrix W n = (w i,j ) n i,j=1 of a Toeplitz graph has a symmetric Toeplitz structure, i.e., w i,j = w |i−j| .
If we assume that m, t 1 , . . . , t m , are fixed (independent of n) and we let the size n grow, then the sequence of adjacency matrices W n can be related to a unique real integrable function f (the symbol) defined on [−π, π] and expanded periodically on R. In this case, according to (2), the entries w ij =f i−j of the matrix W n are defined via the Fourier coefficients of f , where the k-th Fourier coefficient of f is given byf We know that the Fourier coefficientsf k are all in {0, w t1 , . . . , w tm } and that the matrix is symmetric. Note that obviously any such graph is uniquely defined by the first row of its adjacency matrix. On the other hand, we know that w j−1 = w 1,j = w(v 1 , v j ) for j = 1, . . . , n, namely, w j−1 = 0 iff j − 1 ∈ {t 1 , . . . , t m }. From this conditions we can infer that the symbol has a special polynomial structure and in fact it is equal to In such a way, according to (5), our adjacency matrix W n is the matrix T n (f ) (real and symmetric) having the following structure and, as expected, the symbol f is real-valued and such that f (θ) = f (|θ|) for every θ. See Figure 1 for Figure 1: Example of a 1-level Toeplitz graph T 5 (1, w 1 ), (3, w 3 ) . On the left there is a visual representation of the graph while on the right it is explicated the associated adjacency matrix W n which presents the typical Toeplitz structure. In particular, W n has symbol function f (θ) = 2w 1 cos(θ) + 2w 3 cos(3θ).
Along the same lines, we can define d-level Toeplitz graphs as a generalization of the Toeplitz graphs, but beforehand we need to define the set of directions associated to a d-index.
Definition 3.2 Given a d-index t k = ((t k ) 1 , . . . , (t k ) d ) such that 0 t k and t = 0, then define Trivially, it holds that and fix m nonzero real vectors w 1 , . . . , w m , such that We then indicate the components of the vectors w k by the following index notation, is an undirected graph defined by a node set V n = {v k | 1 k n} and a weight function ω such that If there exist m nonzero real numbers such that w k = w k 1 for every k = 1, . . . , m, then the above relation translates into and we will indicate the d-level Toeplitz graph as T n {t 1 , w 1 }, . . . , {t m , w m } . In the case of simple graph, i.e., w k = 1 for every k, then we will indicate the d-level Toeplitz graph just as T n t 1 , . . . , t m . The number of nodes in a d-level Toeplitz graph is equal to D(n) with D(n) = d r=1 n r , while the number of edges is equal to m r=1 D(n − t r ).

Lemma 3.1 A Toeplitz graph is a 1-level Toeplitz graph as in Definition 3.3.
Proof We simply note that, for d = 1, n, t 1 , . . . , t m and the associated w 1 , . . . , w m are scalars, so that the resulting graph has n points and weight function given by as in definition 3.1, completing the proof.
• If we assume that m, , w m }, are fixed (independent of n) and we let the sizes n j grow, j = 1, . . . , d, then the sequence of adjacency matrices can be related to a unique real integrable function f : [−π, π] d → R (the symbol) and expanded periodically on R d . In this case, the entries w i,j =f i−j of the adjacency matrix are defined via the Fourier coefficients of f , where the k-th Fourier coefficient of f is defined according to the equations in (2). Following the same considerations which led to Equation (5), we can summarize everything said till now in the following proposition.
, w m } are fixed and independent of n. Then the adjacency matrix W n of the graph is a symmetric matrix with a d-level Toeplitz structure (see Section 2.1), that is, is clear by definition 3.3, while, by direct computation of the Fourier coefficients of f , we see that As an example, the adjacency matrix of a 2-level Toeplitz graph has the form with j 1 = 0, . . . , n 1 − 1, j 2 = 0, . . . , n 2 − 1 and w 0 = w T 0 . See Figure 2 for an explicit graphic example.
On the left there is a visual representation of the graph while on the right there is the associated adjacency matrix W n , where we used the standard lexicographic ordering to sort the nodes {v k | (1, 1) (k 1 , k 2 ) (4, 3)}. Comparing the adjacency matrix of this example with W n of (9), the red blocks correspond to w 1 and w T 1 , and the blue blocks correspond to w 2 = w T 2 . All the blocks are Toeplitz matrices, and indeed W n is a matrix which possesses a block Toeplitz with Toeplitz blocks (BTTB) structure. In particular, W n has symbol function f (θ 1 , θ 2 ) = 2w 2,0 cos(2θ 1 ) + 2w 1,1 cos(θ 1 + θ 2 ) + 2w 1, 1 cos(θ 1 − θ 2 ): notice that the coefficient of the variable θ 1 refers to the diagonals block while the coefficient of θ 2 refers to the diagonals inside the block. Finally, observe that W n is not connected, since the graph can be decomposed into two disjoint subgraphs G 1 and G 2 having {v 1 , v 3 , v 5 , v 7 , v 9 , v 11 } and {v 2 , v 4 , v 6 , v 8 , v 10 , v 12 } as vertex sets, respectively.

Graphs with uniform local structure: introducing the "diamond"
The idea here is that each node in Definition 3.3 is replaced by a subgraph of fixed dimension ν. For instance, fix a reference simple graph with adjacency matrix W and where [ν] is the standard set of cardinality ν ∈ N. Consider 0 < n ∈ N copies of such a graph, i.e., G(k) = (V (k), E(k)) such that G(k) ≃ G for every k = 1, . . . , n. Indicating the distinct elements of each V (k), k = 1, . . . , n, with the notation v (k,r) , for r = 1, . . . , ν, we can define a new node set V n as the disjoint union of the sets V (k), i.e., Fix now m integers 0 < t 1 < . . . < t m with 1 ≤ m ≤ n − 1, and moreover fix L t1 , . . . , L tm simple linking-graph operators for the reference node set [ν], as in Definition 2.7, along with their uniquely Namely, E n is the disjoint union of all the edge sets E(k) plus all the edges which possibly connect nodes in a graph G(i) with nodes in a graph G(j): two graphs G(i), G(j) are connected iff |i − j| ∈ {t 1 , . . . , t m } and in that case the connection between the nodes of the two graphs is determined by the linking-graph operator L t k (and by its transpose L T t k ). We can define then a kind of symmetric 'weight-graph function ' w : It is not difficult then to prove that the adjacency matrix W G n,ν of the graph (V n , E n ) is of the form Trivially, W G n,ν is a symmetric matrix with a block-Toeplitz structure and symbol function f given by Let us observe that f (θ) is an Hermitian matrix in C ν×ν for every θ ∈ [0, π], and therefore λ j (f (θ)) are real for every j = 1, . . . , ν, as we requested at the end of Subsection 2.2. We will call T G n (t 1 , L t1 ) , . . . , (t m , L tm ) := (V n , E n ) a (simple) diamond Toeplitz graph associated to the graph G. A copy G(k) of the graph G will be called k-th diamond.
See Figure 3 for an example. We can now generalize everything said till now.
be a fixed undirected graph which we call mold graph.
We then indicate the elements of the sets L k by the following index notation, Finally, consider n copies G(k) ≃ G of the mold graph, which we will call diamonds.
and characterized by the weight function w n : V n × V n → R such that otherwise.
The number of nodes in a d-level diamond Toeplitz graph is equal to νD(n) with D((n) = d r=1 n r , while the number of edges is equal to ν m r=1 D(n − t r ).

Corollary 3.3 A d-level
Toeplitz graph is a special case of a d-level diamond Toeplitz graph.

Grid graphs with uniform local structure and main spectral results
The section is divided into two parts. In the first we give the definition of grid graphs with uniform local structure. In the second part we show the links of the above notions with Toeplitz and GLT sequences and we use the latter for proving the main spectral results.

Sequence of grid graphs with uniform local structure
The main idea in this section is to immerse the graphs presented in Section 2.1 inside a bounded regular domain Ω ⊂ R d . We start with a series of definitions in order to give a mathematical rigor to our derivations.
and w is the weight function defined in (6). With abuse of notation we will identify V n =Ṽ n and we will write Observe that now w p k , for k = 1, . . . , m, are not anymore constant vectors as w k but vector-valued functions w p k : otherwise, for α = 1, . . . , c k . It is then not difficult to see that we can express the weight function w p as In other words, taking in mind the role of the reference domain [0, 1] d , x i can be connected to x j only if , for all r = 1, . . . , d. From this property we derive the name of 'grid graphs with local structure'. Naturally, the above notion can be generalized to any domain Ω ⊂ [0, 1] d : as we will see in the next subsection, the only restriction in order to have meaningful spectral properties of the related sequences, is that Ω is regular.
Clearly, V Ω n ′ = n ′ ≤ d r=1 n r = |V n |. Nevertheless, n ′ = n ′ (n) → ∞ as n → ∞. Therefore, with abuse of notation we will keep writing n instead of n ′ . We will indicate such a graph with the notation In the application, as we will see in Section 6, once it is chosen the domain Ω and the kind of discretization technique to solve a differential equation, then the weight function w is fixed accordingly, and consequently the coefficients w 1 , . . . , w m . In particular, it is important to remark that the weight function of T n {[t 1 ], w 1 }, . . . , {[t m ], w m } , will depend on the differential equation and on the discretization technique.  With abuse of notation we will write Remark 3 While in the case of a d-level Toeplitz graph the immersion map ι was introduced naturally as the Hadamard product between the indices of the graph nodes and the natural Cartesian representation of points in R d , the diamond Toeplitz graphs grant another degree of freedom for the immersion map. In Definition 4.3 we decided for the simplest choice, namely lining-up all the nodes of the diamonds along the first axis. Clearly, other choices of the immersion map ι would be able to describe more complex grid geometries.

Definition 4.4 (d-level diamond Toeplitz grid graphs in Ω)
The same definition as in Definition 4.2 where the d-level Toeplitz grid graph is replaced by a d-level diamond Toeplitz grid graph. We will indicate such a graph with the notation

Asymptotic spectral results
We start this section containing the spectral results, by giving the distribution theorem in the Weyl sense in its maximal generality, i.e. for a sequence of weighted (diamond) local grid graphs in Ω, according to the case depicted in Definition 4.4.
where f (θ) is the symbol function defined in (8).
Proof We note that, in the case where ν = 1, a d-level diamond Toeplitz grid graph reduces to a d-level Toeplitz grid graph according to definition 4.2. The conclusion of the theorem is then obvious once we prove our next result, Theorem 4.2.
Proof First of all we observe that our assumption of Ω regular is equivalent to require Ω to be measurable according to Peano-Jordan, which is the fundamental assumption to apply the GLT theory (see [12]).
Assume Ω = [0, 1] d and p(x) ≡ 1 over Ω. Then our sequence of graphs reduces to a sequence of d-level diamond Toeplitz graphs and the proof is over using Proposition 3.4. Assuming now that Ω = [0, 1] d and p is just a Riemann-integrable function over Ω, we decompose the adjacency matrix W G,p n,ν as W G,p n,ν = diag n (p)T n (f ) + E n . The only observation needed here is that {diag n (p)} is a multilevel block GLT with symbol function p, while {T n (f )} is a multilevel block GLT with symbol function f (see item GLT2 in [11]). Moreover, by direct calculation, we see that E n , for n large, can be written as a term of small spectral norm, plus a term of relatively small rank. Therefore, E n is a zero-distributed sequence of matrices and hence a multilevel block GLT with symbol function 0. Summing up we have Now, by the structure of algebra of multilevel block GLT sequences and using the symmetry of the sequence (see item GLT3 in [11]), we conclude that W G,p n,ν n ∼ p(x)f (θ) over is constructed according to Definition 4.3 and the function p is substituted by p |Ω := p(i(x)), where i : Ω ֒→ [0, 1] d is the inclusion map. Since Ω is regular we conclude that p |Ω is Riemann-integrable, p |Ω ≡ p over Ω, and that where c ∈ R is a fixed constant, I is the identity matrix and p : Ω ⊂ [0, 1] d → R is a continuous a.e. function as in Definition 4.4. Then it holds that where ∆ Gn,ν is the graph-Laplacian as in Definition 2.8 and f (θ) is defined in (11).
Proof From Definition 2.8, ∆ Gn,ν = D + K − W G,Ω,p n,ν . By assumption, note that D + K is a GLT sequence with symbol function cp, so that the conclusion follows once again by the algebra structure of the GLT sequences. •

Spectral gaps
In this section we report general results on the gaps between the extreme eigenvalues of adjacency and graph-Laplacian matrices of the type considered so far. The spectral properties of a Toeplitz matrix T n (f ) are well understood by considering f ; in fact, it is known (see e.g. [19,20]) that the spectrum of Recall the following result due to Kac, Murdoch, Szegö, Parter, Widom and later third author, Böttcher, Grudsky, Garoni (see also [4] and references therein for more details and for the history of such results).
At the light of some of the new results presented in [1] and reported in the Appendix A, we can make the following statements within the scope of this paper.
dn ≤ max R g definitely for n → ∞; (ii)g is piecewise Lipschitz and differentiable at x = 1; then it holds that Proof It is immediate from Corollary A.4.

Examples
As a first example we consider the Toeplitz graph T n (1, 1), (2, −6), (3, 1), (4, 1) with corresponding symbol function f (θ) = 2 cos(θ) − 12 cos(2θ) + 2 cos(3θ) + 2 cos(4θ), according to Proposition 3.2. In this case, by symmetry of the symbol f over [−π, π] we can restrict it to θ ∈ [0, π] without affecting the validity of the identity (4). It is easy to verify that the graph is connected, and in Figure 4 and Table 1 it is possible to check the numerical validity of Theorem 5.2, Corollary 5.3 and Theorem A.2.    Figure 4b     k(n) is chosen such that k(n)/n is constant for every fixed n = 10 2 , 5 · 10 3 , 10 3 , 2 · 10 3 . As it can be seen, the relative errors decrease as n increases, in accordance with Theorem A.2. The convergence to zero is slow and not uniform: one of the reasons is that we are using a linear approximation ofg instead ofg itself. In the last row we show the computation of the gap between the biggest and the second-biggest eigenvalues, confirming the prediction of Corollary 5.3.
Due to the symmetry of the symbol f over [−π, π] 2 we can restrict it to (θ 1 , θ 2 ) ∈ [0, π] 2 without affecting the validity of the identity (4). In the following Table 3 Table 3: The first four rows present the relative errors between the eigenvalue λ (n) k(n) and the sampling of the monotone rearrangementf k(n) dn . The index k(n) is chosen such that k(n)/d n is constant for every fixed n = (n, n) with n = 10, 50, 100. The relative errors decrease as n increases, in accordance with Theorem A.2. The convergence is slow and not uniform: one of the reasons is that we are using a simple linear approximation off . In the last row we show relative errors between the gap of the biggest and the second-biggest eigenvalues and the gap between the maximum and second-maximum values of a uniform sampling off , confirming the prediction of Theorem 5.2. In particular, d n λ
In this case, by symmetry of the symbol f over [−π, π] we can restrict it to θ ∈ [0, π] without affecting the validity of the identity (4). Moreover, due to Proposition A.1 and taking into account that we restricted θ to [0, π], with abuse of notation we write where The map θ → 4θ − (k − 1)π is a diffeomorphism between I k and [0, π], and the eigenvalues functions are ordered by magnitude, namely, λ k (f k (θ)) is the k-th eigenvalue function of the matrix (17) over θ ∈ [0, π]. In particular, in this case it holds that In the following Table 4 and Figure 6 it is possible to check numerically the validity of Theorem A.  Table 4: The first four rows present the relative errors between the eigenvalue λ (n) k(n) and the sampling of the monotone rearrangementf k(n) 4n , i.e.: is chosen such that k(n)/(4n) is constant for every fixed 4n = 10 2 , 10 3 , 5 · 10 3 . As it can be seen, the relative errors decrease as n increases, in accordance with Theorem A.2. The convergence is slow and not uniform: one of the reasons is that we are using a linear approximation ofg instead ofg itself. In the last row we show the computation of the gap between the biggest and the second-biggest eigenvalues, confirming the prediction of Corollary 5.3.

Applications to PDEs approximation
The section is divided into three parts where we show that the approximation of a model differential problem by three celebrated approximation techniques leads to sequences of matrices that fall in the theory developed in the previous sections.   (18), for θ ∈ [0, π], and of the eigenvalues λ (n) k of the adjacency matrix W G n , for k = 1, . . . , 10 3 . We have that θ 1 = π/4, θ 2 = π/2, θ 3 = 3π/4 are points of discontinuity since the inequalities in (19) are strict. With the notations of Definition A.1, notice that in this example R f is the union of four disjoint intervals and that λ [0,π] λ 1 (f (θ)), max [0,π] λ 4 (f (θ))]. In Figure 4b, instead, it is possible to observe the validity of Theorem A.2, comparing the distribution of an approximation of the monotone rearrangementf , for x ∈ [0, 1], and the distribution of λ mation with (4m + 1)-points. That is, our model operator with Dirichlet BCs is given by Fixing the diffusion term p(x, y) = 1 + (x − 1/2) 2 + (y − 1/2) 2 and the potential term q(x, y) = e xy , then the operator L is self-adjoint and has purely discrete spectrum, see [8].
Now, if we fix m = 1, n ∈ N and i, j ∈ Z, then the uniform second-order 5-point FD approximation of the (negative) Laplacian operator (i.e., ∆[·] = −(∂ 2 ) is given by Extend now continuously the diffusion term p(x, y) outside B 1/2 , that is, and define the graph G n = T as a sub-graph of G n = T n [1,0], wp , [0, 1], wp ,κ whereV n = (x i , y j ) ∈ [0, 1] 2 (24) p as in (22), like in Definition 2.3. Namely, the mother graphḠ n is the 2-level grid graph on [0, 1] 2 obtained by extending continuously the diffusivity function p to [0, 1] 2 and adding a nontrivial potential termκ which naturally depends on q. On the other hand, the potential term k of the sub-graph G n describes the edge deficiency of nodes in G n compared to the same nodes inḠ n , κ(x i , y j ) = h 2 q(x i , y j ) + (xr,ys)∼(xi,yj) (xr,ys)∈Vn\Vn wp((x r , y s ), (x i , x j )), |κ(x i , y j ) −κ(x i , y j )| = (xr,ys)∼(xi,yj) (xr,ys)∈Vn\Vn wp((x r , y s ), (x i , x j )), see [24, pg.197]. It is easy to check that the boundary points (x i , x j ) ∈ ∂V n are connected at most with two points ofV n \ V n , therefore the potential term κ can be split into three terms Figure 7. The white nodes are the nodes of V n while the gray nodes belong toV n \ V n . The potential termκ of the mother graphḠ n is determined only by the potential term q from (20) while the potential term κ of the sub-graph G n is influenced by the nodes inV n \ V n on the boundary set ∂V n . This influence is due to the presence of Dirichlet BCs in (21). The green connections represent the weighted edges whose end-nodes are both interior nodes of V n , while the red connections represent the weighted edges which have at least one end-node that belongs tō V n \ V n . The potential term κ sums the weight of a red edge to every of its end-nodes which belong to V n .
The given graph-Laplacian ∆ Gn approximates the weighted operator h 2 L. Moreover, by Corollary 4.3 it holds that By the symmetry of g over [−π, π] 2 we can restrict it to B 1/2 × [0, π] 2 without affecting the validity of the identity (4). If we consider now the monotone rearrangementg : [0, 1] → [0, 10] of the symbol g as in Definition A.1, then we can see from Table 5 and Figure 8 that for any index sequence {k(n)}, 1 ≤ k(n) ≤ d n , such that k(n) dn → x ∈ [0, 1] as n → ∞, where d n is the dimension of the graph-Laplacian ∆ Gn . We want to stress out that d n < n 2 since V n V n , but clearly it holds that d n → ∞ as n → ∞, and the Hausdorff distance between the node set V n and the disk B 1/2 is going to zero. Sinceg does not posses an easy analytical expression to calculate, it has been approximated by an equispaced sampling of g over B 1/2 × [0, π] 2 by d n -points and then rearranging it in non-decreasing order. The approximation converges tog as n → ∞, see for example [30]. Finally, see Remark 8. For other applications of the GLT techniques for approximation of partial differential operators see [2] relative errors n = 10 n = 50 n = 80 k(n) dn 0.1 0.0788 0.0094 0.0020 0.5 0.0055 3.3995e-04 1.2565e-04 0.8 0.0100 7.2173e-04 3.9353e-05 1 0.0443 0.0440 0.0360  k(n) is chosen such that k(n)/d n is constant for every fixed n = 10, 50, 80. As it can be seen, the relative errors decrease as n increases, in accordance with Theorem A.2. The convergence to zero is not uniform and slower in some subintervals. k } (orange-dotted line), for n = 80, of the graph-Laplacian ∆ Gn associated to the graph G n defined in (23). In this case we have that d n = 5140 < 80 2 .

Approximations of PDEs vs sequences of weighted diamond graphs: FEM
Consider the model boundary value problem where Ω = (0, 1) and f ∈ L 2 (Ω). We approximate (26) by using the quadratic C 0 B-spline discretization on the uniform mesh with stepsize 1 n , where the basis functions are chosen as the C 0 B-spline of degree 2 defined on the knot sequence { 1 n+1 , 1 n+1 , 2 n+1 , 2 n+1 , ..., n n+1 , n n+1 }. Proceeding as in [15], we trace the problem back to solving a linear system whose stiffness matrix reads as follows: We note that 1 n A n can be seen as the graph-Laplacian of a 1-level diamond Toeplitz graph with nonzero killing term. Namely, according to Definition 2.8, we have that where K n is the diagonal matrix given by if i is even 8 3 otherwise and W G n,2 is the adjacency matrix of the 1-level diamond Toeplitz graph T G 3 (1, L 1 ) , with Combining now Proposition 3.4 and Corollary 4.3 we get that 1 n A n has asymptotic spectral distribution with symbol function f : [−π, π] → C 2×2 given by Remark 5 It is possible to study the multi-dimensional case of the problem in the example above using the fact that for every m, s ∈ N d there exists a permutation matrix Γ m,s of dimension d j=1 m j s j such that for any choice of trigonometric polynomials p : [−π, π] → C sj ×sj , j = 1, ..., d, as stated in Lemma 4 of [15], where ⊗ denotes the usual tensor product. In the d-dimensional case (see, again, [15] and references therein), the discretizing matrix is given by where A n is defined as in the example above and M n is a 2n × 2n matrix with the same block Toeplitz structure as A n and, hence, an analogous symbol function, which we denote by h. Assuming now that the multi-index n = νn = (ν 1 n, ν 2 n, ..., ν d n) for a fixed ν ∈ Q d >0 := {(ν 1 , ..., ν d ) ∈ Q : ν 1 , ..., ν d > 0}, it is immediate to see by the considerations above that n d−2 A n is the linear combination of graph-Laplacians of d-level diamond Toeplitz graphs with spectral distribution given by the following symbol function g : with c k (ν) = ν k ν 1 · · · ν k−1 ν k+1 · · · ν d , k = 1, ..., d.

Approximations of PDEs vs sequences of weighted d-level graphs: IgA approach
In this section we consider the approximation of the same differential equation (26) by using the Galerkin B-splines approximation of degree ν in every direction: while the approximation of standard Q ν Lagrangian FEM considered in Section 6.2 leads to a symbol which is a linear trigonometric polynomial in every variable, but C ν d ×ν d Hermitian matrix valued, here the symbol is scalar valued, but the degree of the trigonometric polynomial is much higher. For example, by means of cubic C 2 B-spline discretization the (normalized) stiffness matrix is given by 0 · · · 0 3 60 9 360 Observe that the principal R (n−4)×(n−4) -submatrix (highlighted in blue) is an exact symmetric Toeplitz matrix while globally A n is not Toeplitz due to the presence of perturbations near the boundary points (highlighted in yellow). This behavior is influenced by the presence of BCs and the specific choice for the test-functions (B-spline with C 3 local regularity and C 2 global regularity). For those reasons, A n can not be representative of the graph-Laplacian of a graph in the form G n = (T n (t 1 , w 1 ), . . . , (t m , w m ) , κ), even if it is clearly the graph-Laplacian for another kind of graph G n which does not own globally the symmetries of the graphs studied in sections 3 and 4. Nevertheless, since the perturbations are local, it happens that {A n } n ∼ λ f (θ) where f is the same symbol function of the graph-Laplacian of G n = T n 1, 30 240 , 1, 48 240 , 1, 2 240 , κ , κ ≡ 0, namely, f (θ) = 160 240 − 60 240 cos(θ) − 96 240 cos(2θ) − 4 240 cos(3θ). For a complete treatment of the study of the eigenvalue distribution for Galerkin B-splines approximations we refer to [17], where many examples are provided along the exposure.

Remark 6
There is another quite important difference between this case and the case of Section 6.1. In the FD case, the nodes of the graph were representative of the physical domain while in this case, even if the node set can be immersed in [0, 1], the nodes represent the base functions of the test-functions set. That said, given A n , it will be of interest to calculate the corresponding weight function w on the node set of the physical domain: this approach could lead some insight about the problem of the presence of a fixed number of outliers in the spectrum of A n , see [6, Chapter 5.1.2 p. 153].

Conclusions, open problems, and future work
We have defined general classes of graph sequences having a grid geometry with a uniform local structure in a domain Ω ⊂ [0, 1] d , d ≥ 1. With the only weak requirement that Ω is Lebesgue measurable with boundary of zero Lebesgue measure, we have shown that the underlying sequences of adjacency matrices have a canonical eigenvalue distribution, in the Weyl sense, with a symbol f being a trigonometric polynomial in the d Fourier variables: as specific cases, we mention standard Toeplitz graphs, when Ω = [0, 1], and d-level Toeplitz graphs when Ω = [0, 1] d , but also matrices coming from the approximation of differential operators by local techniques, including Finite Differences, Finite Elements, Isogeometric Analysis etc. In such a case we considered block structures and weighted graphs, from the perspective of GLT sequences, where the tools taken from the latter field have resulted crucial for deducing all the asymptotic spectral results. In particular, the knowledge of the symbol and of basic analytical features have been employed for deducing a lot of information on the eigenvalue structure, including precise asymptotics on the gaps between the largest eigenvalues.
Many open problems remain, ranging from a deeper analysis of the matrix sequences arising from different families of Finite Element approximations of multidimensional differential problems to the study of the convergence features of the ordered asymptotic spectra to the rearrangement of the corresponding symbol (see also the study and discussions in [17,1]).

A Appendix
Fix a square matrix sequence {X n,ν } n of dimension d n , with symbol function g : D ⊂ R m → C ν×ν as Definition 2.9. Observe that g is not unique and in general not univariate. To avoid this, we will soon introduce in Definition A.1 the notion of monotone rearrangement of the symbol. In order to simplify the notations and since all the cases we investigate in this paper can be lead back to this situation, we make the following assumptions: Assumptions (AS1) D is compact and of the form Ω × [−π, π] d with Ω ⊆ [0, 1] d , and therefore m = 2d; (AS2) g(y) = g(x, θ) = p(x)f (θ), with (x, θ) ∈ Ω × (−π, π) d and p : Ω → R, f : [−π, π] d → C ν×ν ; (AS3) p : Ω → R piecewise continuous; (AS4) every component f i,j : [−π, π] d → C of f is continuous; (AS5) f (θ) has real eigenvalues for every θ ∈ [−π, π] d .
Since we are assuming the eigenvalues to be real, then for convenience notation we will order the eigenvalue functions λ k (p(x)f (θ)) by magnitude, namely λ 1 (p(x)f (θ)) ≤ . . . ≤ λ ν (p(x)f (θ)). This kind of ordering could affect the global regularity of the eigenvalue functions, but it does not affect the global regularity of the monotone rearrangement of the symbol. Nevertheless, by well-known results (see [21]), items (AS3) and (AS4) imply that λ k (p(x)f (θ)) is at least piecewise continuous for every k = 1, . . . , ν. We have the following result.
• Remark 7 Trivially, the maps I k ∋ θ → νθ − (2k − 1 − ν)π are diffeomorphism between I k and [−π, π] d . Therefore, the image set of ν k=1 λ k (f k (θ)) over I k is exactly the image set of λ k (f (θ)) over [−π, π] d . The next definition of monotone rearrangementà la Hardy-Littlewood is crucial for the understanding of the asymptotic distribution of the eigenvalues of {X n,ν } n .
(29b) We callg the monotone rearrangement of g.
Clearly,g is well-defined, univariate, monotone strictly increasing and right-continuous. Within our assumptions on g, it is then easily possible to extend [9,Theorem 3.4] for this multi-variate matrixvalued case, and it holds that (i)g(0) = min R g ,g(1) = max R g ; (ii) lim n→∞ Σ λ (F, X n,ν ) = 1 0 F (g(x))dµ 1 (x) . We have the following results, see [1,Section 3.4]. The statements and proofs are exactly the same, we just adjusted them to fit with the multi-index notation we are adopting here.