# A continuum limit for the Laplace operator on metric graphs

Sidney Holden1
<EMAIL_ADDRESS>&Geoffrey Vasil1
<EMAIL_ADDRESS>
(1 School of Mathematics and Maxwell Institute for Mathematical Sciences
The University of Edinburgh, Edinburgh EH9 3JZ, United Kingdom
)

###### Abstract

Metric graphs — with continuous real-valued edges — have been studied
extensively as domains for solutions of PDEs. We study the action of the
Laplace operator on each edge as a metric graph fills in a $n-$dimensional
Riemannian manifold, e.g., a very dense spider web in the planar disc or ball
of wire within the $2-$sphere. We derive an $n$-dimensional limiting
differential operator whose low-order eigenmodes approximate the low-order
eigenmodes of the metric graph Laplace operator. High-density metric graphs
exhibit nontrivial, continuum-like behaviour at low frequencies and purely
graph-like behaviour at high frequencies. In the continuum limit, we find the
emergence of a new symmetric tensor field but with an unusual distance scaling
compared to the traditional Riemannian metric. The limiting operator from a
metric graph implies a possible new kind of continuous geometry.

_Keywords_ Quantum graph $\cdot$ Continuum model $\cdot$ Emergent geometry

##  1 Introduction

The concept of indivisible continuum media is a remarkably useful fiction. In
particular, we know that liquids, solids and gases comprise innumerable
interacting discrete elements, where atoms and molecules are traditionally the
relevant microscopic components. But this situation is not a requirement.
Instead, a good deal of modern modelling derives macroscopic “constitutive
relations" and “coarse-grained” dynamics for smooth functions representing the
approximate state of exotic, complex systems. These efforts include biological
“active matter”, where motile microorganisms take the role of molecules [1, 2,
3]. Another example is from finance, where the microscopic degrees of freedom
represent potential price histories of assets and commodities [4]. Coarse-
graining is often productive because the resulting macrophysics reduces to
local partial differential equations in many situations. Moreover, we usually
know much about the resulting PDEs, where the diffusion equation, wave
equation, or at least a Laplace-like operator dominates the story.

This paper considers a new kind of macroscopic media that emerges
systematically from a dense network of interconnected “wires” or “conduits”.
Our microphysics already occurs as a mixture of continuous and discrete
components. We begin with a collection of one-dimensional continuous edges,
each acting as the “spatial” domain for numerous possible local differential
equations. Multiple edges can meet at their mutual endpoints, with relevant
matching conditions at these junctures. Overall, finding the dynamics on each
edge demands information from the entire global system of vertices. In this
sense, we envision a kind of “mesoscopic” system.

The appropriate mathematical structure for our analysis is a “metric graph”: a
discrete graph with physical length along each edge. We embed such a graph
within some common space (e.g. 2D or 3D Euclidean space or possibly a smooth
manifold like the 2-sphere). We pose a typical dynamical model as a one-
dimensional differential equation on each edge. The overall question we focus
on is:

> What partial differential equations result from taking the “continuum limit”
> of the metric graph, where the vertices become dense in the embedded space,
> and the edge lengths shrink to zero?

In this paper, we solve this question for the eigenvalue equation of the one-
dimensional Laplace operator $\Delta=d^{2}/dx^{2}$ on each edge of a given
graph $G$.

##  2 Previous work

On the one hand, there has been significant work studying differential
equations on metric graphs. Much of the literature [5] focuses on quantum
Hamiltonians $\mathcal{H}=-\Delta+V(x)$, where the system
$\\{\mathcal{H},G,\text{boundary conditions}\\}$ is called a “quantum graph”
for its applications to the Schrödinger operator in quantum systems. On the
other hand, much more work has considered the continuum limit of discrete
lattices. As of yet, there has been minimal work combining the two.

###  2.1 Physical applications

Graphs, especially lattices, have garnered increased attention recently
regarding their applications in forefront scientific areas. In particular, for
nanomaterials due to their remarkable mechanical, electric and optical
properties [6, 7, 8, 9]. Quantum graphs model graphene and carbon nanotubes
[10, 11, 12], and also stacks of graphene using multi-layer quantum graphs
[13]. Recently, quantum graphs have shown promise in new approaches to
metamaterial design [14].

There has been significant research on the localization of charge carriers at
low energy electronic states in disordered semiconductors. Localization
landscape theory provides one promising model for this phenomenon [15, 16, 17,
18, 19], which has been also adapted in the quantum graph setting [20].
Anderson localization more generally has been studied extensively on quantum
graphs [21, 22, 23, 24, 25]. Recently, localization has been shown in random
laser networks both experimentally [26] and on quantum graphs [27].
Applications are for electrically pumped UV lasers [28], biosensing [29] and
optical information processing [30].

Neuroscience and brain dynamics are especially intriguing applications of
high-density graph behaviour [31, 32, 33]. Comparatively, little has included
quantum graph models; recent work typically builds discrete adjacency matrices
by solving inverse problems using brain-scan data. [34] models localised
neuronal activity in the brain using the spectrum of an estimate of the
Laplacian. Mapping the neuro-topology of animal brains is still a rapidly
developing area [35, 36].

Recently, an equal-edge 60-vertex quantum graph model provided the low-energy
eigenmodes of protonated methane [37]. The alternative approach involves
costly numerical simulations of the full Schrödinger operator decomposition
for the multi-particle system [38]. Quantum graphs result as the limit of
tubular networks [39]. They also describe the Casimir effect [40, 41, 42].
More relevant work is [43, 44, 45] on infinite metric graphs with applications
in the human circulatory system [46, 47] and in population dynamics [48]. [49]
model a “tissue membrane" with a dense equilateral network of strings.

###  2.2 Mathematical results

Hybrid monographs on both quantum graphs and PDEs on metric graphs [50, 51,
52] include studies of the non-linear Schrödinger and reaction-diffusion
equations. Metric graphs result from the limits of thin-branched structures
[53, 39]. Quantum graphs also apply to control problems [54] and semigroup
methods for evolution equations [55].

Other efforts have created a calculus on metric graphs, defining (up to)
second-order differential operators, integral theorems, a natural Laplacian,
its associated PDEs, and existence and uniqueness results [56, 57]. These
works also provide an alternative eigenvalue formula, linking the
combinatorial graph to equilateral quantum graphs [58, 59]. For a graph with
uniform edge lengths, $\ell$, then

$\displaystyle-k^{2}\ \in\ \sigma(\Delta)\qquad\iff\qquad-2\sin^{2}(k\ell/2)\
\in\ \sigma(D^{-1}L),$ (1)

where $\sigma$ is the operator spectrum, $\Delta$ is the metric-graph
Laplacian, $D$ is the diagonal degree matrix and $L=A-D$ is the discrete
Laplacian matrix, with adjacency matrix, $A$.

Notably, $-2\sin^{2}(\ell k/2)\sim k^{2}\ell^{2}/2$ as $\ell\to 0$. It seems
natural in that $k\sim\ell^{-1}$. However, the extra factor of $1/2$ is not
from a simple difference of scaling convention; it points to interesting
fundamental aspects of metric graphs in the continuum limit.

In the current work, we consider the limit as the graph becomes _dense_ in
$\mathbb{R}^{d}$. Edges are compact parametrised curves attached at vertices
$v\in\mathbb{R}^{d}$. The limiting approach is similar to that used for
finite-element meshes but with different motivations. Metric graphs rarely
appear as objects in Euclidean space. Little is known about arbitrarily dense
embedded graphs. _Infinite_ graphs generally correspond to the more general
class of graphs with infinite vertices, each of finite degree. Several works
on spectral problems for equilateral, symmetric and self-similar infinite
metric graphs. There have also been advances in spectral properties for
infinite metric graphs with edge lengths having no positive lower bound. For
infinite metric graphs, [60] finds self-adjointness for the Laplace operator
with variable vertex conditions.

Based on preliminary results in [61], a recent series of papers [62, 63] on
continuum limits of the quantum Hamiltonian on metric lattices proves the
continuum version of eq. (1). By linking the Euclidean Laplace operator to its
discrete form on a lattice and then to the metric graph Laplace operator, they
show that the dimension-scaled metric graph Laplace operator approximates the
Euclidean Laplace operator in the sense of norm resolvent convergence. We
introduce the relationship via our PDE and empirically, in 5.1. But [61]
approach still starts from a PDE and studies what emerges regarding discrete
graphs. Lastly, [64] shows Ollivier to Ricci curvature convergence in the
continuum limit of random geometric graphs by considering edges as having
simply real positive weights.

##  3 Mathematical setup

###  3.1 Metric graphs

A metric graph is just what one might expect intuitively, like a spider web,
electrical power grid, water-supply system or British Rail. Nevertheless, we
make a few precise definitions for conducting later mathematical derivations.
A metric graph, $G=(V,\,E)$ is a collection of vertices, $V$, and edges, $E$,
with the additional feature that edges are connected intervals of
$\mathbb{R}$, with the standard notion of Euclidean distance. A weighted and
oriented incidence matrix, $B_{v,e}$, encodes all necessary structural
information.

Two vertices, $v_{1},v_{2}$ are edge incident, $v_{1}\sim_{e}v_{2}$, if they
are connected by an edge $e$. The endpoints of any two edges are identically
equal if they meet at a vertex $v\in V$. Two edges, $e_{1},e_{2}$ are vertex
incident, $e_{1}\sim_{v}e_{2}$, if they share a common vertex, $v\in e_{1}$
and $v\in e_{2}$. Finally, we say that an edge, $e$, is incident to a vertex,
$e\sim v$, when $v$ is one of the endpoints of $e$. We omit subscripts and use
just $\sim$ when the implication is unambiguous.

Without loss of generality, we can assign local coordinates to each edge
$x_{1}\in[0,\ell_{1}]\simeq e_{1}$, $x_{2}\in[0,\ell_{2}]\simeq e_{2}$. An
oriented, weighted incidence matrix encodes which interval endpoints of
$e_{1}\sim_{v}e_{2}$ represent the common vertex, $v$. Assuming
$e\simeq[0_{e},\ell_{e}]$, where the subscript $e$ in $0_{e}=0$ distinguishes
edge starting points, then for the $v$th row of $B_{v,e}$ we define

$\displaystyle B_{v,e}\ =\ \begin{cases}+\ell_{e}&\text{if}\quad v\simeq
0_{e}\in[0_{e},\ell_{e}],\\\ -\ell_{e}&\text{if}\quad
v\simeq\ell_{e}\in[0_{e},\ell_{e}],\\\ \ 0&\text{otherwise.}\end{cases}$ (2)

The incidence matrix helps to define functions and their derivatives on $G$.
The signs of the incidence matrix entries,
$\hat{B}_{v,e}=\text{sign}(B_{v,e})$, are unit normal vectors pointing away
from $v$.

Global coordinates on $G$ are denoted $(e,x)$, where $x$ is the local
coordinate on $e$. When considering graphs embedded in regions
$\Omega\in\mathbb{R}^{d}$, we will drop the global/local distinction and
consider only $x\in G\subset\Omega$. The edges of an embedded metric graph are
parameterised curves, but with the caveat that they need not be straight. A
balled-up fishing net is the same metric graph as its laid-flat counterpart.

###  3.2 Functions and derivatives

We denote real-valued functions $f:G\to\mathbb{R}$ as $f(e,x)$ for $x\in e$,
e.g. figure 1. A function is continuous on $G$ if it is continuous on all
individual edges up to and including the endpoints. Recall that the endpoints
of adjacent edges are literally the same unique point $v\in G$. Therefore, a
continuous function on $G$ also has unique values at the vertices, denoted
$f(v)$.

Differentiability on metric graphs is straightforward in the middle of the
edges. For all $\varepsilon>0$ and $x$, $x+\varepsilon\,\hat{n}\in e$ where
$\hat{n}\in\left\\{\pm 1\right\\}$, a real-valued function $f:G\to\mathbb{R}$
is edgewise differentiable if the following exists,

$\displaystyle D_{\hat{n}}f(e,x)\ =\ \lim_{\varepsilon\to
0}\frac{f(e,x+\varepsilon\,\hat{n})-f(e,x)}{\varepsilon}\ =\
\hat{n}\,\frac{d}{dx}f(e,x).$ (3)

Functions need not have unique derivatives at vertices. The derivative at $v$
is an edgewise function

$\displaystyle D_{e}f(v)\ =\ \hat{B}_{v,e}\,\frac{d}{dx}f(e,x)\Big{|}_{x=v},$
(4)

which is now a coordinate-parameterisation invariant directional derivative
away from $v$.

Figure 1: Real, scalar-valued function (in red) on a star-shaped metric graph
(in black).

###  3.3 The Laplace operator

We define the $L_{2}$ inner product over the whole graph as the sum of $L_{2}$
inner products on each edge. For real-valued functions $f,g$ on $G$,

$\displaystyle\left<f,g\right>_{G}\ =\ \sum_{e\in
G}\int_{0_{e}}^{\ell_{e}}f(e,x)g(e,x)\,\text{d}x,$ (5)

with norm, $||f||_{G}^{2}=\left<f,f\right>_{G}$. For derivatives,

$\displaystyle||Df||_{G}^{2}\ =\ \sum_{e\in
G}\int_{0_{e}}^{\ell_{e}}\left|\frac{d}{dx}f(e,x)\right|^{2}\,\text{d}x,$ (6)

which is independent of the edge orientations because $\hat{B}_{v,e}^{2}=1$.
By analogy with standard manifolds, we define the eigenfunctions of the
Laplacian from the stationary points of the Rayleigh quotient, with eigenvalue

$\displaystyle k\ =\ \min_{||f||_{G}=1}||Df||_{G}.$ (7)

For any globally continuous and differentiable $g(e,x)$, the Euler-Lagrange
equations are,

$\displaystyle\sum_{e\in
G}\left\\{\left[g(e,x)\frac{d}{dx}f(e,x)\right]_{0_{e}}^{\ell_{e}}-\int_{0_{e}}^{\ell_{e}}g(e,x)\left(k^{2}+\frac{d^{2}}{dx^{2}}\right)f(e,x)\,\text{d}x\right\\}\
=\ 0.$ (8)

Taking $g(e,x)=0$ at all vertices implies the edge differential equations,

$\displaystyle-\frac{d^{2}}{dx^{2}}f(e,x)\ =\ k^{2}f(e,x).$ (9)

The remaining conditions are,

$\displaystyle\sum_{e\in
G}\left\\{g(e,\ell_{e})\frac{d}{dx}f(e,\ell_{e})-g(e,0_{e})\frac{d}{dx}f(e,0_{e})\right\\}\
=\ 0.$ (10)

At any vertex, $v\in V$ we can pick a preferred vertex, $v^{\prime}$, and
assume $g(e,v)=-\delta_{v,v^{\prime}}$, giving the celebrated Kirchhoff
condition,

$\displaystyle\sum_{e\sim v}D_{e}f(v)\ =\ 0.$ (11)

The edge orientations $\hat{B}_{v,e}$ account for the sign differences between
$\ell_{e}$ and $0_{e}$ in eq. (10).

###  3.4 A reduced vertex eigensystem

We transform the global eigensystem (posed on edges as differential equations
with vertex conditions) as a system of algebraic equations on the vertices
alone. We solve eq. (9) analytically and substitute the result into the
Kirchhoff conditions.

We focus on a typical vertex, $v$, with multiple edges $e\sim v$. Without loss
of generality, all edges are oriented outward from $v$, i.e..
$\hat{B}_{v,e}=+1$. Along a typical edge for $w\sim_{e}v$,

$\displaystyle f(e,x)\ =\
f(v)\,\frac{\sin(k(\ell_{e}-x))}{\sin(k\ell_{e})}+f(w)\,\frac{\sin(kx)}{\sin(k\ell_{e})},$
(12)

manifestly expressing continuity, $f(e,0)=f(v)$, $f(e,\ell_{e})=f(w)$. The
Kirchhoff condition is

$\displaystyle\sum_{e\sim v}D_{e}f(v)\ =\ k\sum_{w\sim
v}\left[\csc(k\ell_{vw})\,f(w)-\cot(k\ell_{vw})f(v)\right]\ =\ 0.$ (13)

Altogether, this is an algebraic system over each $v\in V$,

$\displaystyle\text{L}(k)\cdot\text{f}_{V}\ =\ 0,$ (14)

where $\text{f}_{V}$ is the column vector of values of $f$ at the vertices,
and the matrix

$\displaystyle\text{L}_{vw}(k)\ =\ \begin{cases}\sum_{v^{\prime}\sim
v}\cot\left(k\ell_{vv^{\prime}}\right)&\text{if}\quad w=v,\\\
-\csc\left(k\ell_{vw}\right)&\text{if}\quad w\sim v,\\\
0&\text{otherwise}.\end{cases}$ (15)

We solve the finite system numerically for eigenpairs $(k,\text{f}_{V})$ (see
Appendix D), fully specifying the eigenfunctions in eq. (12).

##  4 Analysis of the continuum limit

We now construct concrete graphs by embedding the vertices in a Riemannian
manifold, $G\subset\Omega$. The vertices are definite points in $v\in M$, with
more freedom for the edges. We could require that the edges be fully embedded
in $e\subset M$, e.g., finite geodesics. Or we could embed the manifold within
a higher Euclidean space and consider the edges as straight lines. In the
continuum limit, both options will give the same result: for any two vertices,
$w\sim v$, we take the limit as $\ell_{vw}\to 0$ but allow the number of
vertices to fill in densely in $M$.

We assume that $M$ has a Riemannian structure, e.g. $\nabla$, $\nabla\cdot$,
implying a natural volume measure,

$\displaystyle\int_{M}\,\text{d}{x}\ =\ \text{Vol}(M).$ (16)

An underlying Riemannian structure is not essential to our analysis. We find
our eventual results are independent of a particular underlying Riemannian
metric. However, the graph requires a vertex number density, $\mu(x)$ such
that integrating over the manifold gives the total number of vertices,

$\displaystyle\sum_{v\in G}\ =\ \int_{G\subset M}\mu(x)\,\text{d}{x}\ =\ |V|.$
(17)

The function, $\mu(x)$, is the Radon-Nikodym derivative with respect to the
geometric measure on the manifold.

Using straightforward trigonometric identities, eq. (13) is equivalent to

$\displaystyle\sum_{w\sim
v}\cot(k\ell_{vw})(f(v)-f(w))-\tan\\!\left(k\ell_{vw}/2\right)f(w)\ =\ 0,$
(18)

which is the stationary-point condition of the quadratic functional

$\displaystyle Q(\text{f}_{V})\ =\ \sum_{v}\sum_{w\sim
v}\left\\{\cot(k\ell_{vw})\left(f(v)-f(w)\right)^{2}-2\tan\left(k\ell_{vw}/2\right)f(v)f(w)\right\\}.$
(19)

As $z\to 0$, $\cot(z)\sim z^{-1}$, $2\tan(z/2)\sim z$. For $r_{vw}=w-v$ with
$|r_{vw}|=\ell_{vw}$,

$\displaystyle f(w)\ =\ f(v)+r_{vw}\cdot\nabla
f(v)+\mathcal{O}(\ell_{vw}^{2}).$ (20)

Put together, the discrete quadratic form becomes

$\displaystyle Q(\text{f}_{V})\ \sim\ \sum_{v}\sum_{w\sim
v}\frac{1}{|r_{vw}|}\left(r_{vw}\cdot\nabla
f(v)\right)^{2}-k^{2}|r_{vw}|f(v)^{2}+\mathcal{O}(\ell_{vw}^{2}).$ (21)

Now comes a major crux of our derivation. We define the following 2-tensor at
each vertex,

$\displaystyle R(v)\ \equiv\ \sum_{w\sim
v}|r_{vw}|\hat{r}_{vw}\otimes\hat{r}_{vw},\qquad\text{where}\qquad\text{tr}\\!\left({R(v)}\right)\
=\ \sum_{w\sim v}|r_{vw}|.$ (22)

Therefore,

$\displaystyle Q(\text{f}_{V})\ \sim\ \sum_{v\in G}\nabla f(v)\cdot
R(v)\cdot\nabla
f(v)-k^{2}\,\text{tr}\\!\left({R(v)}\right)f(v)^{2}+\mathcal{O}(\ell_{vw}^{2}).$
(23)

Finally, we replace the single sum with an integral over the vertex density,

$\displaystyle Q(\text{f}_{V})\to\mathcal{Q}[f]\ =\ \int_{G\subset
M}\left[\nabla f(x)\cdot R(x)\cdot\nabla
f(x)-k^{2}\,\text{tr}\\!\left({R(x)}\right)f(x)^{2}\right]\mu(x)\,\text{d}{x}.$
(24)

Defining the vertex density can occur in several equivalent ways. Trivially,
we can use the empirical measure on the vertices

$\displaystyle\mu_{V}(x)\ =\ \sum_{v\in V}\delta(x-x_{v}).$ (25)

The empirical measure converges weakly to a continuum limit as the graph
becomes dense. However, we can also approximate the empirical measure with
several other measures converging weakly to the same eventual result. A
particularly useful alternative results from a dual graph. That is,

$\displaystyle\mu_{\widetilde{V}}(x)\ =\
\frac{1}{|\mathcal{A}_{v}|}\quad\text{for}\quad x\in\mathcal{A}_{v},$ (26)

where $\mathcal{A}_{v}$ is the dual cell area (volume in general) surrounding
a given vertex, $v$. At this stage, $\mu_{\widetilde{V}}(x)$ is not a
continuous function; however, smoother but more complicated interpolations
exist as well. The cell-area density approximates the empirical measure for
functions of bounded variation, which is also sufficient to define the weak
derivative [65].

Computing the variational derivative of the continuum functional produces the
eigenproblem in the dense limit. We construct $G$ to have a boundary $\partial
G$ such that $\partial G=V\cap\partial\Omega$. $\partial G$ fills in
$\partial\Omega$ in the same sense that $G$ fills in $\Omega$. For
$\varepsilon>0$ and $\varphi:G\to\mathbb{R}$ with $\varphi\big{|}_{\partial
G}=0$, then
$\frac{d}{d\varepsilon}\mathcal{Q}\left[f+\varepsilon\varphi\right]\big{|}_{\varepsilon=0}=0$
gives

$\displaystyle\frac{1}{\mu\,\text{tr}\\!\left({R}\right)}\nabla\cdot\left(\mu\,R\cdot\nabla
f\right)\ =\ -k^{2}f\quad\text{for}\quad x\in G\subset M.$ (27)

A “free” boundary condition results in starting with a graph with no
distinction between boundary and interior vertices. In that case,

$\displaystyle\hat{n}\cdot R\cdot\nabla f=0\quad\text{for}\quad x\in\partial
M,$ (28)

where $\hat{n}$ is the unit normal to the boundary. If the original graph has
pinned boundary vertices, then $f(x)=0$ for $x\in\partial M$.

The reduced system resembles the standard Laplace-Beltrami (LB) operator on a
manifold with an augmented local metric tensor, $g(x)$. However, in that case,

$\displaystyle\Delta_{LB}f\ =\ \nabla\cdot\nabla{f}\ =\
\frac{1}{\sqrt{\text{det}(g)}}\,\partial\cdot(\sqrt{\text{det}(g)}\,g^{-1}\partial
f).$ (29)

There is no clear way to put eq. (27) into one-to-one correspondence with an
LB operator, even if we impose conditions on the type of Riemannian metric.
The $R$ tensor is defined similarly to a discrete approximation of a
Riemannian metric, but with a different distance scaling, i.e.,

$\displaystyle g(v)^{-1}\ =\ \frac{d}{n_{v}}\sum_{w\sim
v}\hat{r}_{vw}\otimes\hat{r}_{vw},$ (30)

where $n_{v}$ is the number of neighbouring vertices (see appendix A).

As far as we know, it is impossible to relate $R$ in any straightforward way
to any other standard tensors (e.g. Ricci, Einstein) commonly defined on
Riemannian manifolds. Moreover, where the local volume measure on a Riemannian
manifold follows from the metric determinant, there appears no simple way to
produce the vertex density, $\mu(x)$, in terms of a function of $R(x)$ alone.
The underlying presence of a metric graph “substance” cannot be described
purely in geometrical terms.

An additional observation is that the Riemannian metrics included in the
$\nabla$ operators in eq. (27) and eq. (29) do not need to be the same. We
look at what PDEs on $\Omega$ arise from a metric graph $G$ filling in
$\Omega$. Though we don’t explicitly address it here, we have also looked at
the inverse problem of determining which metric graphs within $\Omega$ produce
a particular PDE on $\Omega$.

##  5 Specific graph examples

(a)

(b)

(c)

(d)

Figure 2: The examples of metric graphs $G$ and their embedding spaces
$\Omega$ that we study. From top to bottom, $|V|\ \approx\ 200\to 2500$. From
left to right, 2(a) and 2(b) are the square (homogeneous and isotropic) and
rectangular (homogeneous and anisotropic) lattices filling in the 2-cube, 2(c)
is the spider web filling in the 2-ball, and 2(d) is the Goldberg polyhedron
filling in the 2-sphere.

The following sections demonstrate empirically how the limiting PDE model
approximates the large-scale dynamics of graphs with various characteristics.
Starting with the motivation of previously discovered dispersion relations on
lattices, we show the power of our continuum model for increasingly complex
graphs. If the $R$ tensor of a graph is constant, then the graph is
homogeneous, and inhomogeneous otherwise. If $R$ is a multiple of the
identity, then the graph is isotropic, and otherwise anisotropic.

###  5.1 Fully solvable homogeneous and isotropic graphs

####  5.1.1 Two-dimensional regular toroidal lattice

A regular, toroidal lattice (e.g. figure 2(a) with periodic boundary
conditions) is one straightforward example that highlights much of the
essential nature of a metric graph and its continuum limit–without needing to
solve the limiting PDE. The Kirchhoff vertex condition 11 simplifies
considerably in this case with identical rectilinear edges,

$\displaystyle\cos(k\ell)f(v)-\frac{1}{4}\sum_{w\sim v}f(w)\ =\ 0.$ (31)

We consider vertices connected in the cardinal directions (N-S-E-W) at
distances $\ell$. Or alternatively, vertices connected in the ordinal
directions (NE-SW-SE-NW) at distances $\ell\to\sqrt{2}\,\ell$; see figure 3.

In each case, we assume a periodic lattice in the $x,y$

$\displaystyle f(x,y)\ =\ A\,e^{i(k_{x}x+k_{y}y)}+\text{c.c.}$ (32)

The total lattice periodicity quantises $k_{x},k_{y}$.

Cardinal connectivity — In this case, the total wavenumber dispersion
conditions reduces to

$\displaystyle\cos(k\ell)\ =\
\frac{1}{2}\left(\cos(k_{x}\ell)+\cos(k_{y}\ell)\right).$ (33)

To lowest order in $\ell\to 0$ we find

$\displaystyle 2k^{2}\ =\ k_{x}^{2}+k_{y}^{2}\ +\ \mathcal{O}(\ell^{2}).$ (34)

We find the usual Pythagorean relationship between wave numbers, just as in
the continuous periodic torus, albeit with an extra factor of $2$.

Ordinal connectivity — In this case, with neighbouring diagonal vertices,

$\displaystyle\cos(\sqrt{2}k\ell)\ =\ \cos(k_{x}\ell)\cos(k_{y}\ell).$ (35)

This result is particularly intriguing for the following reason:

> Eq. (35) is the Pythagorean Theorem for a right triangle on the surface of a
> sphere with spherical angles $a=k_{x}\ell$, $b=k_{y}\ell$ and
> $c=\sqrt{2}k\ell$.

We have no idea why the spherical geometry for wavenumbers arises in this
example. A spherical triangle automatically reduces to a Cartesian triangle as
$\ell\to 0$, and we again find eq. (34). Nontrivial geometric structure in
spectral space can lead to many interesting phenomena, including symmetry-
induced topological protection [66].

We can also see that the case with cardinal connectivity gives a similar
result only with the wavenumbers rotated by $\pi/4$, i.e.,

$\displaystyle\frac{\cos(a)+\cos(b)}{2}\ =\
\cos\\!\left(\frac{a+b}{2}\right)\cos\\!\left(\frac{a-b}{2}\right).$ (36)

Cardinal + ordinal connectivity — Finally (see figure 3(c)), we can combine
the above two examples into the total dispersion relation,

$\displaystyle\frac{\cos(k\ell)-\cos(\tfrac{(k_{x}+k_{y})\ell}{2})\cos(\tfrac{(k_{x}-k_{y})\ell}{2})}{\sin(k\ell)}+\frac{\cos(\sqrt{2}k\ell)-\cos(k_{x}\ell)\cos(k_{y}\ell)}{\sin(\sqrt{2}k\ell)}\
=\ 0.$ (37)

It is less clear how to interpret the combination in terms of the spherical
Pythagorean theorem. Nevertheless, in the limit $\ell\to 0$, we still recover
the Cartesian Pythagorean theorem in eq. (34), including the factor of $2$.

(a)

(b)

(c)

Figure 3: Three types of homogeneous and isotropic graph geometries: a)
cardinal, b) ordinal and c) both.

####  5.1.2 Homogeneity and isotropy in the continuum limit

The exactly solvable torus displays limiting behaviour that we expect from a
two-dimensional continuous Laplacian, albeit with a factor of 2. We better
understand the continuum limit by considering isotropic and homogeneous cases
in higher dimensions.

For an isotropic and homogeneous underlying graph, the only 2nd-order
isotropic and homogeneous tensor is proportional to the $d$-dimensional
identity matrix,

$\displaystyle R\ =\ \frac{\text{tr}\\!\left({R}\right)}{d}\,I_{d}.$ (38)

The only information about the underlying graph comes through
$\text{tr}\\!\left({R}\right)=c$ for some constant. For a vertex connected to
$N$ neighbours each a uniform distance $\ell$ away, $c=\ell\,N$. For a vertex
with cardinal and ordinal connectivity, we have
$\text{tr}\\!\left({R}\right)=4\ell+4\sqrt{2}\ell$. Moreover,
$\text{tr}\\!\left({R}\right)$ and the vertex density, $\mu$, are both
spatially constant and commute out of the derivatives. Therefore,

$\displaystyle\frac{1}{\mu\,\text{tr}\\!\left({R}\right)}\nabla\cdot\left(\mu
R\cdot\nabla f\right)\ =\ -k^{2}f\qquad\Longrightarrow\qquad\frac{1}{d}\Delta
f\ =\ -k^{2}f.$ (39)

The entire problem reduces to the standard Laplace-Beltrami operator, but
weighted by the spatial dimension of the graph embedding. It is clear why we
should expect the dispersion relations for the three toroidal lattices to have
the same limiting forms.

An interesting corollary is that the large-scale eigenfunctions of a dense
graph approximate well those from a continuous underlying geometry. However,
the eigenvalue will have an adjustment factor equal to the embedding-space
dimension.

As an illustrative example: Consider two drums, one made from a homogeneous
and isotropic wire mesh and the other made from whatever materials comprise
everyday drums. Both drums have the same total tension and mass density. The
shapes of the fundamental modes of both drums match nearly perfectly. However,
the frequencies of the fundamental modes will be different by a factor of
$\sqrt{d}=\sqrt{2}$. Overall, the continuum “graph material" can resemble many
other materials but is less stiff by the dimension of the embedding space.
Example 5.2 shows this behaviour in more detail.

###  5.2 Square lattice within the 2-cube

Consider now the locally homogeneous unit square
$\Omega=[0,1]\times[0,1]\subset\mathbb{R}^{2}$ (figure 2(a)) with clamped
boundaries ($\tilde{f}\big{|}_{\partial\Omega}=0$). Applying the same
reasoning as the above combinatorial torus, the continuum model in eq. (39)
with clamped boundaries has the recognisable solutions

$\displaystyle\tilde{k}_{m,n}\ =\
\frac{\pi}{\sqrt{2}}\sqrt{m^{2}+n^{2}}\;,\qquad\tilde{f}_{m,n}\ \propto\
\sin(m\pi x)\sin(n\pi y).$ (40)

Figure 10(a) shows the relative error
$\eta(k_{m,n},\tilde{k}_{m,n})=|(k_{m,n}-\tilde{k}_{m,n})/\tilde{k}_{m,n}|$
between corresponding eigenvalues. We use Newton’s method to solve
$\det(L(k))=0$, inputting a high-density uniform $k$-sampling of a given
interval. Each output $k$-value is matched to a particular $\tilde{k}_{m,n}$;
all of the first several eigenvalues match to within $10^{-2}$.

Figure 11(a) shows the metric graph $L^{2}$ error between functions from
corresponding eigenspaces. Comparing the eigenfunctions of the discrete and
continuous systems requires some care because of degeneracy owing to 4-fold
rotation symmetry.

We determine $\text{f}_{V}$ from the nullspace of $\text{L}(k_{m,n})$ using
the scipy.linalg package. However, because $k_{m,n}=k_{n,m}$ a given null
vector can correspond to an arbitrary linear combination,

$\displaystyle\text{f}_{V}\ \approx\
\alpha_{n,m}\tilde{f}_{n,m}\big{|}_{V}+\alpha_{m,n}\tilde{f}_{m,n}\big{|}_{V}.$
(41)

where $\alpha_{m,n}=\left<f_{m,n},\tilde{f}_{m,n}\right>_{\\!G}$ for
$||\tilde{f}_{m,n}||_{G}=1$ We measure the convergence between the ODE- and
PDE-eigenfunctions using the metric-graph norm. Figure 4(a) shows an $m\times
n=3\times 3$ grid of eigenfunctions for $|V|=\mathcal{O}(10^{2})$. Figure 4(b)
shows the PDE eigenfunctions, restricted to $G$.

(a) Graph ODE eigenfunctions $f_{m,n}$

(b) Graph PDE eigenfunctions $\tilde{f}_{m,n}$

Figure 4: Square lattice $G$ with $|V|\,=\,\mathcal{O}(10^{2})$ embedded
within the unit square $\Omega$. (rows, columns) $=(m,n)\in\\{1,2,3\\}^{2}$.
Figure 4(a) corresponds to the graph eigenmodes $f_{m,n}$, and 4(b) to the
continuum eigenmodes $\tilde{f}_{m,n}\propto\sin(m\pi x)\sin(n\pi y)$
restricted to $G\subset\Omega$.

###  5.3 Rectangular lattice within the 2-cube

For a homogeneous rectangular lattice (figure 2(b)) on the unit square
$\Omega$, $\mu$ is still constant, but now $R=\text{diag}(\text{r})$ for some
vector of constants $r>0$. Without loss of generality,
$\text{tr}\\!\left({R}\right)=1$. In local coordinates
$\sum_{i=1}^{d}r_{i}\partial_{i}^{2}f=-k^{2}f$. On the clamped unit square
with $r=(1/3,2/3)$,

$\displaystyle\tilde{k}_{m,n}\ =\
\pi\sqrt{\frac{m^{2}+2n^{2}}{3}}\;,\qquad\tilde{f}_{m,n}\ \propto\ \sin(m\pi
x)\sin(n\pi y).$ (42)

Again, we see how the local structure of the underlying graph leads to a
notable difference in the eigenvalues, even while eigenfunctions remain quite
innocent. Figures 10(b) and 11(b) show the convergence of $(k,f)$ for
increasing $|V|$. Figure 5(a) shows a $3\times 3$ grid of eigenfunctions for
$|V|=\mathcal{O}(10^{2})$. Figure 5(b) shows the PDE eigenfunctions,
restricted to $G$.

(a) Graph eigenfunctions $f_{m,n}$

(b) PDE eigenfunctions $\tilde{f}_{m,n}$

Figure 5: Rectangular lattice $G$ with $|V|\,=\,\mathcal{O}(10^{2})$ embedded
within the unit square $\Omega$. (rows, columns) $=(m,n)\in\\{1,2,3\\}^{2}$.
Figure 5(a) corresponds to the graph eigenmodes $f_{m,n}$, and 5(b) to the
continuum eigenmodes $\tilde{f}_{m,n}\propto\sin(m\pi x)\sin(n\pi y)$
restricted to $G\subset\Omega$.

###  5.4 Spider web within the 2-ball

We now study a subclass of inhomogeneous and anisotropic spider web-like
graphs.

The spider web (figure 2(c)) fills in the unit disc
$\Omega=\\{\sqrt{x^{2}+y^{2}}\leq 1\\}\in\mathbb{R}^{2}$. With vertex spacings
$\,\text{d}{r},\,\text{d}{\theta}$ for each coordinate $r,\theta$
respectively, the spider web has a central vertex at the origin, and for any
other non-boundary vertex; in Cartesian components,
$v(r,\theta)=r(\cos(\theta),\sin(\theta))$. Vertices on the unit circle lack
outer radial neighbours. In a polar basis,

$\displaystyle\begin{bmatrix}R_{rr}&R_{r\theta}\\\ R_{\theta
r}&R_{\theta\theta}\end{bmatrix}\ =\
\begin{bmatrix}2\,\text{d}{r}+4r\sin^{3}\\!\left(\\!\frac{\,\text{d}{\theta}}{2}\right)&0\\\
0&\frac{r\sin^{2}(\\!\,\text{d}{\theta})}{\sin\left(\frac{\,\text{d}{\theta}}{2}\right)}\end{bmatrix}\
\sim\ 2\begin{bmatrix}\,\text{d}{r}&0\\\ 0&r\,\text{d}\theta\end{bmatrix}.$
(43)

We determine the vertex number density from the inverse of the vertex cell
area, $A_{v}\sim r\,\text{d}{r}\,\text{d}{\theta}$. We require that
$\,\text{d}{\theta}=2\pi/M$ for some total number of angular vertices, $M$.
However, we allow $\,\text{d}{r}$ to depend smoothly on $r$. It is, therefore,
useful to define a function $\rho(r)$ where

$\displaystyle\,\text{d}{r}\ =\ \rho(r)\,\text{d}\theta\;,\qquad\mu(r)\ =\
\frac{1}{r\rho(r)\,\text{d}{\theta}^{2}}.$ (44)

Putting everything together,

$\displaystyle\frac{\nabla\cdot\left(\mu R\cdot\nabla
f\right)}{\mu\,\text{tr}\\!\left({R}\right)}\ =\
\frac{1}{r(r+\rho(r))}\left(r\rho(r)\frac{\partial^{2}}{\partial
r^{2}}+\frac{\partial^{2}}{\partial\theta^{2}}\right)f(r,\theta)\ =\
-k^{2}f(r,\theta).$ (45)

Notably, no form of $\rho(r)$ can render eq. (45) into a metric-based
Laplacian.

Angular symmetry allows Fourier solutions,
$f(r,\theta)=f_{m}(r)e^{im\theta}+\text{c.c.}$. Therefore,

$\displaystyle
r\rho(r)f_{m}^{\prime\prime}(r)+\left[r(r+\rho(r))k^{2}-m^{2}\right]f_{m}(r)\
=\ 0.$ (46)

The original discrete graph problem (13) also simplifies into decoupled radial
systems. In this case, $\,\text{d}{r}=\,\text{d}{r}(r)$ takes discrete values
on the edges. So at a vertex $v(r,\theta)$, edges to neighbours toward and
away from the central vertex have lengths $\,\text{d}{r}_{-}(r)$ and a
$\,\text{d}{r}_{+}(r)$ respectively,

$\displaystyle\frac{f_{m}(r+\,\text{d}{r}_{+})+f_{m}(r-\\!\,\text{d}{r}_{-})}{\sin(k\,\text{d}{r})}\hskip
227.62204pt$
$\displaystyle+\;\left[2\left(\frac{\cos(m\,\text{d}{\theta})}{\sin\\!\left(kr\,\text{d}{\sigma}\right)}-\cot(kr\,\text{d}{\sigma})\right)-\left(\cot(k\,\text{d}{r}_{+})+\cot(k\,\text{d}{r}_{-})\right)\right]f_{m}(r)\
=\ 0,$ (47)

where $\,\text{d}{\sigma}=2\sin(\,\text{d}{\theta}/2)$ and $|m|\leq
M=2\pi/\,\text{d}{\theta}$. As $\,\text{d}{r},\,\text{d}{\theta}\to 0$ with
$\rho=\,\text{d}{r}/\,\text{d}{\theta}\sim\mathcal{O}(1)$ eq. (47) reduces
exactly to eq. (46) at leading order. The vertex (with $M$ neighbours) at the
origin requires special care. For the discrete case,

$\displaystyle
f_{m}(0)\cos(k\,\text{d}{r})-f_{m}(\\!\,\text{d}{r})\,\delta_{m,0}\ =\ 0,$
(48)

with Kronecker-$\delta_{m,0}$. As $\,\text{d}{r}\to 0$, the centre conditions
are

$\displaystyle f_{m}(0)\ =\ 0\quad\text{if}\quad m\neq 0\qquad\text{or}\qquad
f_{0}^{\prime}(0)\ =\ 0\quad\text{if}\quad m=0$ (49)

Analysing the PDE within and $\varepsilon$-ball around the origin produces
identical centre conditions.

Bessel-like solutions — Some specific forms of $\rho(r)$ allow analytical
solutions to eq. (46). In particular, if $\rho(r)=r/\gamma$, then for $m\neq
0$,

$\displaystyle
f_{m}(r)\propto\sqrt{r}J_{\nu}(\sqrt{1+\gamma}\,kr)\qquad\text{where}\qquad\nu\
=\ \frac{\sqrt{4\gamma\,m^{2}+1}}{2},$ (50)

where $J_{\nu}$ is the non-integer order regular Bessel function, the spectral
parameter, $k=\zeta_{n,\nu}/\sqrt{1+\gamma}$, and $J_{\nu}(\zeta_{n,\nu})=0$
results from the clamped outer boundary condition; i.e., $f_{m}(r=1)=0$. The
$m=0$ solution uses the Neumann function branch, which reduces to

$\displaystyle
f_{0}(r)\propto\cos(\sqrt{1+\gamma}\,kr)\qquad\text{where}\qquad k\ =\
\frac{(2n+1)\pi}{2\sqrt{1+\gamma}}.$ (51)

The spider web with $\gamma=1$ is in figure 2(c). Figures 10(c) and 11(c) show
the convergence of $(k,f)$ for increasing $|V|$. Figure 6(a) shows a $3\times
3$ grid of eigenfunctions for $|V|=\mathcal{O}(10^{2})$. Figure 6(b) shows the
PDE eigenfunctions restricted to $G$.

(a) Graph eigenfunctions $f_{m,n}$

(b) PDE eigenfunctions $\tilde{f}_{m,n}$

Figure 6: Spider web $G$ with $|V|\,=\,\mathcal{O}(10^{2})$ embedded within
the unit disc $\Omega$. (rows, columns)
$=(m,n)\in(\\{0,1,2\\}\times\\{1,2,3\\})$. Figure 6(a) corresponds to the
graph eigenmodes $f_{m,n}$, and 6(b) to the continuum eigenmodes
$\tilde{f}_{m,n}\propto\cos(\sqrt{2}\tilde{k}_{m,n})\delta_{m,0}+\sqrt{r}J_{\nu_{m}}(\sqrt{2}\,\tilde{k}_{m,n}r)\left(1-\delta_{m,0}\right)$
restricted to $G\subset\Omega$.

###  5.5 Goldberg polyhedron within the 2-sphere

The 2-sphere is the first embedding domain $\Omega$ with Euler characteristic
$\chi=2>1$. There are no trivial examples to introduce as with the periodic
lattices within the unit square. The only regular tilings of the 2-sphere are
the _low-density_ Platonic solids: from the tetrahedron ($|V|=4$) to the
icosahedron ($|V|=12$). As a result, industry standards for high-density
tilings focus on preserving the shapes of faces. For example, video games and
film animations generally use either geodesic or Goldberg polyhedra, which are
built strictly from triangles or pentagons+hexagons, respectively. We
construct polyhedra with icosahedral symmetry by starting with the icosahedron
and successively perform the Conway truncate $``t"$ and dual $``d"$
operations, projecting the vertices to the 2-sphere at each step. This results
in a sequence of Goldberg polyhedra and their dual geodesic polyhedra (see
figure 7).

Figure 7: Construction of geodesic polyhedra and their dual Goldberg polyhedra
by successively applying the Conway operations dual $``d"$ and truncate $``t"$
using the icosahedron as the seed.

The geodesic polyhedra contain the 12 vertices of the icosahedron as a subset,
and similarly, each of the 12 pentagons of the Goldberg polyhedra can be
identified with a vertex. We study the resulting sequence of Goldberg
polyhedra, starting with the soccer ball. For increasing vertex density, $R$
and $\mu$ are close to but never homogeneous over $G$. If they were indeed
homogeneous, we would obtain yet again the rescaled version of the eigenvalue
equation for the Laplace operator eq. (39), this time the solutions being
spherical harmonics

$\displaystyle\tilde{k}_{j}\ =\
\frac{1}{\sqrt{2}}\sqrt{j(j+1)}\;,\qquad\tilde{f}_{j,m,}\ \propto\
Y_{j,m}(\theta,\phi).$ (52)

where $j=0,1,2,\ldots$ and $-j\leq m\leq j$. Figure 10(d) shows the relative
error between the graph eigenvalues and spherical harmonic eigenvalues as
solid lines. The graph deformations are such that while solutions $(k,f)$
accord with $(\tilde{k},\tilde{f})$, there is an offset that does not vanish
in the continuum limit. To study this offset, we consider solutions $(k,f)$ as
perturbations of $(\tilde{k},\tilde{f})$. We explicitly treat the soccer ball
and then show results for increasing vertex density. Eigenfunctions for
$j=0,1,2$ are shown in figure 8.

(a) Graph eigenfunctions $f_{m,n}$

(b) PDE eigenfunctions $\tilde{f}_{m,n}$

Figure 8: Soccer ball $G$ with $|V|\,=\,60$ embedded within the unit disc
$\Omega$. (rows, columns) $=(j,m)=(j,2j+1)$ for $j=0,1,2$. Figure 8(a)
corresponds to the graph eigenmodes $f_{m,n}$, and 8(b) to the continuum
eigenmodes $\tilde{f}_{m,n}\propto Y_{j}^{m}(\theta,\phi)$ restricted to
$G\subset\Omega$.

####  5.5.1 First order perturbation theory for the soccer ball

Each vertex $v$ of the soccer ball has three neighbours at a set distance
$\ell=\sqrt{\frac{2}{109}\left(29-9\sqrt{5}\right)}$. As a result, the vertex
condition eq. (13) can be expressed in terms of the unweighted adjacency
matrix $(A)_{vw}=\delta_{vw}$, where now if $\gamma$ is an eigenvalue of $A$,
we have

$\displaystyle k\ =\ \frac{1}{\ell}\arccos\left(\frac{1}{3}\gamma\right).$
(53)

We can write down the characteristic polynomial of $A$ explicitly

$\displaystyle p(\gamma)\ =\
(\gamma-3)(\gamma-1)^{9}(\gamma+2)^{4}\left(\gamma^{2}-\gamma-3\right)^{5}\left(\gamma^{2}+\gamma-4\right)^{4}\hskip
56.9055pt$
$\displaystyle\left(\gamma^{2}+\gamma-1\right)^{5}\left(\gamma^{2}+3\gamma+1\right)^{3}\left(\gamma^{4}-3\gamma^{3}-2\gamma^{2}+7\gamma+1\right)^{3},$
(54)

and so we end up with explicit formulae and degeneracies for each $k$. We show
$k_{jm}^{2}$ for $0\leq j\leq 3$ in figure 9(a).

(a)

(b)

Figure 9: Visual comparison between $k_{jm}^{2}$ and $\lambda_{j}^{(0)}$
(left) and $\lambda_{j}^{(0)}+\lambda_{jm}^{(1)}$ (right). The
$\lambda_{j}^{(0)}=\frac{j(j+1)}{d}$ are the eigenvalues of $\Delta_{S^{2}}$
being perturbed, inversely scaled by the dimension $d=2$ (as in eq. (39)).
Each $\lambda_{j}^{(0)}$ has degeneracy $2j+1=1,3,5,7,9$. We see in 9(b) that
the splittings in the graph eigenvalues $k_{jm}^{2}$ are approximated well by
the $\lambda_{jm}^{(1)}$.

To calculate the splittings in 9(a) we use the first order perturbation
condition eq. (78) derived in Appendix E. For the soccer ball (and the higher
density Goldberg polyhedra more generally), we first define
$S^{(0)}=\mu^{(0)}R^{(0)}$, where $R^{(0)},\mu^{(0)}$ are the graph tensor and
vertex density of the theoretical homogeneous and isotropic graph. Here,
$R^{(0)}=r_{0}I$, where $r_{0}$ is the mean diagonal entry of $R$ over $G$

$\displaystyle r_{0}\ =\ \frac{1}{|V|}\sum_{v}\sum_{w\sim
v}\frac{|r_{vw}|}{\text{deg}(v)}.$ (55)

By definition, $\mu^{(0)}$ is the mean volume per vertex

$\displaystyle\mu^{(0)}\ =\ \frac{|V|}{|\Omega|}.$ (56)

$\lambda^{(0)}$ and $f^{(0)}$ are given by eq. (52), where we know that the
$f^{(0)}$ form a complete orthonormal basis of $L^{2}(S^{2})$. Substituting
everything into the continuum PDE, we first obtain the leading order equation
eq. (39). The first order equation eq. (78) we solve, finding for each $j$,
$2j+1$ first order corrections $\lambda_{jm}^{(1)}$ for $-j\leq m\leq j$. We
see in 9(b) that $\lambda_{jm}^{(1)}$ accurately reflects the splittings in
$k_{jm}^{2}$. At such a low vertex density, perturbations for increasing $j$
become less reliable. For $j=5,6$ for example,
$\lambda_{jm}^{(1)}=\mathcal{O}(\lambda_{j}^{(0)})$, and so these are not
plotted. However, for increasing vertex densities, we see significant
improvements in convergence.

The dashed lines in figure 10(d) correspond to the mean relative errors

$\displaystyle\eta_{j}\ =\
\frac{1}{2j+1}\sum_{m=-j}^{j}\eta\left(k_{jm}\;,\;\sqrt{-\tilde{k}_{j}^{2}-\lambda_{jm}^{(1)}}\right).$
(57)

For the icosphere, the mean splitting $\lambda_{jm}^{(1)}$ for a given $j$ is
0, because of the following identity

$\displaystyle\sum_{m=-j}^{j}\sum_{v\in V}\nabla
Y_{jm}\nabla\bar{Y}_{jm}:R(v)+\lambda_{j}^{(0)}\text{tr}\\!\left({R}\right)Y_{jm}\bar{Y}_{jm}\
=\ 0,$ (58)

which follows because of the following two properties of spherical harmonics,

$\displaystyle\sum_{m=-j}^{j}\nabla Y_{jm}\otimes\nabla\bar{Y}_{jm}\ =\
\frac{(2j+1)}{4\pi}\frac{j(j+1)}{2}I,\qquad\sum_{m=-j}^{j}Y_{jm}\bar{Y}_{jm}\
=\ \frac{(2j+1)}{4\pi}.$ (59)

For $j=1,2$, there are no splittings–the $k_{jm}$ are the same for each $m$.
This explains why we see that $\lambda_{jm}^{(1)}=0$. For $j>2$, we see an
order of magnitude improvement in the convergence. For the eigenfunctions, we
see a similar offset in convergence to leading order as for the eigenvalues.
We expect that a first order perturbation calculation would track this offset
just as for the eigenvalues. We plan to address this in future work.

##  6 Discussion and conclusions

These limiting structures present a new kind of continuum geometry encoding
information about underlying semi-discrete graph structures. Our work derives
a PDE whose solutions approximate the eigenmodes of the Laplace operator on
metric graphs in the continuum limit. A simple extension to Schrödinger
operators, with nontrivial potential functions, would transition to the
continuum with virtually no modifications to the analysis. The potential part
of a Schrödinger operator contains no derivatives and therefore limits to a
locally averaged version of the underlying graph potential function.

There are many possible directions for future study. Applications of the
continuum model likely apply to materials science and condensed matter
physics. We also expect to study applications to systems as diverse as force
networks in sand piles and, most interestingly, to dense neuron networks
within the brain. It is exciting to think about which dynamics may be modelled
on our graph material. At the same time, it is clear that our mathematical
setup is highly idealised. However, the Laplace operator with its natural
boundary conditions is a ubiquitous object in mathematical modelling. We hope
that our work motivates research in new directions.

##  Appendix A Metric tensor derivation

Suppose we are in a manifold with distance and consider two nearby points,
$x\approx y$. Then consider the ratio,

$\displaystyle S(x,y)\ =\ \frac{(f(x)-f(y))^{2}}{\text{dist}(x,y)^{2}}.$ (60)

If $\gamma(t)$ is a parameterised curve connecting $x=\gamma(0)$ to
$y=\gamma(t)$. Then as $t\to 0$,

$\displaystyle S(x,y)\ \sim\ \frac{(v\cdot\,\text{d}{f(x)})^{2}}{v\cdot
g(x)\cdot v},\quad\text{as}\quad y\to x,$ (61)

where the tangent vector $v=\gamma^{\prime}(0)$. The local metric tensor at
the point $x$ is $g(x)$ and $\,\text{d}{f(x)}$ is the local gradient 1-form.

The tangent vector, $v$, does its best to “point” from $x$ to $y$; it is the
closest analogue to an edge in a metric graph. Because of the definition of
the ratio, the magnitude of $v$ cancels out of the total expression.

In a local orthogonal basis,

$\displaystyle g\ =\
\sum_{i=1}^{d}h_{i}\,\,\text{d}{x}_{i}\otimes\,\text{d}{x}_{i},$ (62)

with scale factors $h_{i}>0$. A general unit vector is

$\displaystyle
v_{\hat{n}}=\sum_{i=1}^{d}\frac{\hat{n}_{i}}{\sqrt{h_{i}}}\,\partial_{i},\qquad\text{where}\qquad\sum_{i=1}^{d}\hat{n}_{i}^{2}\
=\ 1,$ (63)

where $\,\text{d}{x}_{i}\cdot\partial_{j}=\delta_{i,j}$. That is,
$\hat{v}_{\hat{n}}\cdot g\cdot\hat{v}_{\hat{n}}=1$.

For a typical point, $y$, there is some direction, $\hat{n}$ “pointing” at it.
Therefore,

$\displaystyle S(x,y)\ \sim\ (\hat{v}_{\hat{n}}\cdot\,\text{d}{f(x)})^{2}$
(64)

If we average this expression over the sphere $S^{d-1}$, then

$\displaystyle\lim_{y\to x}\left<S(x,y)\right>_{S^{d-1}}$ $\displaystyle=$
$\displaystyle\frac{1}{|S^{d-1}|}\int_{S^{d-1}}(\hat{v}_{\hat{n}}\cdot\,\text{d}{f(x)})^{2}\,\text{d}{\hat{n}}$
(65) $\displaystyle=$
$\displaystyle\sum_{i,j}\frac{1}{\sqrt{h_{i}h_{j}}}\,\left<\hat{n}_{i}\hat{n}_{j}\right>_{S^{d-1}}\,\partial_{i}f\,\partial_{j}f$

However, for a unit vector in the sphere,
$\left<\hat{n}_{i}\hat{n}_{j}\right>_{S^{d-1}}=d^{-1}\,\delta_{i,j}$.
Therefore,

$\displaystyle\lim_{y\to x}\left<S(x,y)\right>_{S^{d-1}}$ $\displaystyle=$
$\displaystyle\frac{1}{d}\,\text{d}{f(x)}\cdot
g(x)^{-1}\cdot\,\text{d}{f(x)},$ (66)

where

$\displaystyle g(x)^{-1}\ =\
\sum_{i,j}\frac{1}{h_{i}}\partial_{i}\otimes\partial_{i}.$ (67)

To make an analogy with the metric graph, with normal vectors along the edges,
$\hat{r}_{vw}$,

$\displaystyle\frac{d}{n_{v}}\sum_{w\sim v}\hat{r}_{vw}\otimes\hat{r}_{vw}\
\approx\ g(x_{v})^{-1},$ (68)

where $n_{v}$ is the number of neighbours for vertex, $v$, and $d$ is the
dimension of space.

##  Appendix B Eigenvalue comparison

(a) Square lattice

(b) Rectangular lattice

(c) Spider web

(d) Goldberg polyhedron

Figure 10: We plot the relative error
$\eta(k_{\alpha},\tilde{k}_{\alpha})=|(k_{\alpha}-\tilde{k}_{\alpha})/\tilde{k}_{\alpha}|$
vs. $|V|$. For the lattices and spider web, $\alpha=(m,n)$, and for the
Goldberg polyhedra, $\alpha=j$. Dashed lines in 10(d) correspond to the
relative error between $k$ and the _first_ order approximation.

##  Appendix C Eigenfunction comparison

In the following, if the $k_{m,n}-$eigenspace is 2-fold degenerate, then the
graph $L^{2}$ error is calculated for only one $f_{m,n}$.

(a) Square lattice

(b) Rectangular lattice

(c) Spider web

Figure 11: We plot the graph $L^{2}$ error
$\chi(f_{\alpha},\tilde{f}_{\alpha})=||f_{\alpha}-P_{\tilde{U}}(\tilde{f}_{\alpha})||_{G}$
vs. $|V|$. For the lattices and spider web, $\alpha=(m,n)$.

##  Appendix D Numerical algorithm

While the original differential eq. (9) contained $-k^{2}$ as a traditional
eigenvalue, the matrix $L(k)$ in eq. (14) depends trigonometrically on $k$,
and we cannot use standard numerical linear algebra libraries to solve a
finite eigenvalue problem. Instead, $L(k)$ has a nontrivial null vector of
vertex values if its determinant vanishes,

$\displaystyle D(k)\ :=\ \det(L(k))\ =\ 0.$ (69)

We find $k$ solutions by iterating Newton’s method to convergence,
$k\leftarrow k-D(k)/D^{\prime}(k)$. Unfortunately, computing large
determinants is numerically prohibitive and unstable. However, the presence of
the logarithmic derivative in Newton’s method is fortuitous because of
Jacobi’s trace formula,

$\displaystyle\frac{d}{dk}\log D(k)\ =\ \text{tr}\\!\left({L(k)^{-1}\cdot
L^{\prime}(k)}\right).$ (70)

Therefore, we iterate the following

$\displaystyle k\ \leftarrow\
k-\frac{1}{\text{tr}\\!\left({X}\right)},\qquad\text{where}\qquad L(k)\cdot X\
=\ \frac{d}{dk}L(k).$ (71)

For moderately-sized problems (up to $|V|=\mathcal{O}(10^{3})$ for our
purposes), we can use standard techniques to solve for $X$. For larger
problems, we use stochastic trace estimation on
$\text{tr}\\!\left({L^{-1}L^{\prime}}\right)=\text{tr}\\!\left({X}\right)$,
which drastically reduces the computation time. The essential idea is as
follows: if $u_{i}$ is taken from any $V$-dimensional probability distribution
$U$ with expected value $\mu=0$ and covariance matrix $\Sigma=I$, then the
quadratic form $u_{i}^{T}Xu_{i}$ has expected value given by,

$\displaystyle\mathbb{E}\left[u_{i}^{T}Xu_{i}\right]\ =\
\text{tr}\\!\left({X\Sigma}\right)+\mu^{T}X\mu\ =\
\text{tr}\\!\left({X}\right).$ (72)

Then sampling $m$ $u_{i}$’s and taking a simple mean
$\text{tr}_{m}\\!\left[X\right]:=\frac{1}{m}\sum_{i=1}^{m}u_{i}^{T}Xu_{i}$
reduces the variance. This idea dates back to [67] and [68], and there are
many modern extensions in cite and cite. For our purposes, calculating
$u_{i}^{T}L^{-1}L^{\prime}u_{i}$ amounts to solving $L\cdot x=u_{i}$, for some
_vector_ $x$.

##  Appendix E First order perturbation condition

Here, we derive a first order perturbation condition of eq. (27) for when we
do not have explicit formulae for $\mu,R$ but can consider $S:=\mu R$ as
$S^{(0)}$ with some first order correction $S^{(1)}$ simply defined as
$S^{(1)}=S-S^{(0)}$. Analogously, if each $k_{m}$ is $N-$fold degenerate, then
to first order we express $(k_{m},f_{mn})$ as

$\displaystyle-k_{m}^{2}=\lambda_{mn}\ =\
\lambda_{m}^{(0)}+\lambda_{mn}^{(1)}\quad,\quad f_{mn}\ =\
f_{mn}^{(0)}+f_{mn}^{(1)}\quad,\quad n\ =\ 1\ldots N.$ (73)

Substituting everything into the continuum PDE, we obtain the leading order
equation by design

$\displaystyle\nabla\cdot\left(S^{(0)}\cdot\nabla f_{mn}^{(0)}\right)\ =\
\lambda^{(0)}\text{tr}\\!\left({S^{(0)}}\right)f_{mn}^{(0)}.$ (74)

If we then define

$\displaystyle\mathcal{L}^{(\alpha)}(f)\ =\
\nabla\cdot\left(S^{(\alpha)}\cdot\nabla
f\right)-\lambda_{m}^{(0)}\text{tr}\\!\left({S^{(\alpha)}}\right)f,$ (75)

then to first order

$\displaystyle\mathcal{L}^{(0)}\left(f_{mn}^{(1)}\right)+\mathcal{L}^{(1)}\left(f_{mn}^{(0)}\right)\
=\ \lambda_{mn}^{(1)}\text{tr}\\!\left({S^{(0)}}\right)f_{mn}^{(0)}.$ (76)

Re-expressing $S^{(1)}=S-S^{(0)}$, we obtain the terms in eq. (74) again which
cancel. Next, let
$\overline{f}_{m}^{(0)}=\sum_{n^{\prime}=1}^{N}u_{mn^{\prime}}f_{mn^{\prime}}^{(0)}$
be some linear combination of the $\lambda_{m}^{(0)}-$eigenbasis with some
arbitrary vector of coefficients $(u_{m})_{n^{\prime}}$. After projecting onto
$\overline{f}_{m}^{(0)}$ and integrating by parts, we see that

$\displaystyle\int_{\Omega}\overline{f}_{m}^{(0)}\mathcal{L}^{(0)}\left(f_{mn}^{(1)}\right)dx\
=\ 0.$ (77)

The terms on the LHS can be written as sums by using
$\int_{\Omega}\mu(x)dx=\sum_{v\in V}$. Finally, to calculate each of the $N$
splittings, we form the system

$\displaystyle Au_{m}\ =\ \lambda_{mn}^{(1)}Bu_{m}$ (78)

where

$\displaystyle\left(A\right)_{ij}\ =\ \frac{1}{|V|}\sum_{v\in V}\nabla
f_{mi}^{(0)}\cdot R\cdot\nabla
f_{mj}^{(0)}+\lambda_{m}^{(0)}\text{tr}\\!\left({R}\right)f_{mi}^{(0)}f_{mj}^{(0)}$
(79)

and

$\displaystyle\left(B\right)_{ij}\ =\ -\sum_{n^{\prime}=1}^{N}\frac{d\cdot
r_{0}}{|\Omega|}\int_{\Omega}f_{mi}^{(0)}f_{mj}^{(0)}dx.$ (80)

## References

  * Saintillan and Shelley [2012] David Saintillan and Michael J Shelley.  Emergence of coherent structures and large-scale flows in motile suspensions.  _Journal of the Royal Society Interface_ , 9(68):571–585, 2012.
  * Bililign et al. [2022] Ephraim S Bililign, Florencio Balboa Usabiaga, Yehuda A Ganan, Alexis Poncet, Vishal Soni, Sofia Magkiriadou, Michael J Shelley, Denis Bartolo, and William Irvine.  Motile dislocations knead odd crystals into whorls.  _Nature Physics_ , 18(2):212–218, 2022.
  * Shahzad et al. [2022] Faisal Shahzad, Wasim Jamshed, Rabha W Ibrahim, Farheen Aslam, Tag El Din, M El Sayed, Hamiden Abd El-Wahed Khalifa, Fayza Abdel Aziz ElSeabee, et al.  Galerkin finite element analysis for magnetized radiative-reactive walters-b nanofluid with motile microorganisms on a riga plate.  _Scientific Reports_ , 12(1):1–20, 2022.
  * MacBeth and Merville [1979] James D MacBeth and Larry J Merville.  An empirical examination of the black-scholes call option pricing model.  _The journal of finance_ , 34(5):1173–1186, 1979\.
  * Berkolaiko and Kuchment [2013] Gregory Berkolaiko and Peter Kuchment.  _Introduction to quantum graphs_ , volume 186 of _Mathematical Surveys and Monographs_.  American Mathematical Soc., 2013.
  * Chiu et al. [2022] Christie S Chiu, Annette N Carroll, Nicolas Regnault, and Andrew A Houck.  Line-graph-lattice crystal structures of stoichiometric materials.  _Physical Review Research_ , 4(2):023063, 2022\.
  * Hagel et al. [2021] Joakim Hagel, Samuel Brem, Christopher Linderälv, Paul Erhart, and Ermin Malic.  Exciton landscape in van der waals heterostructures.  _Physical Review Research_ , 3(4):043217, 2021\.
  * Cao et al. [2018a] Yuan Cao, Valla Fatemi, Shiang Fang, Kenji Watanabe, Takashi Taniguchi, Efthimios Kaxiras, and Pablo Jarillo-Herrero.  Unconventional superconductivity in magic-angle graphene superlattices.  _Nature_ , 556(7699):43–50, 2018a.
  * Cao et al. [2018b] Yuan Cao, Valla Fatemi, Ahmet Demir, Shiang Fang, Spencer L Tomarken, Jason Y Luo, Javier D Sanchez-Yamagishi, Kenji Watanabe, Takashi Taniguchi, Efthimios Kaxiras, et al.  Correlated insulator behaviour at half-filling in magic-angle graphene superlattices.  _Nature_ , 556(7699):80–84, 2018b.
  * Amovilli et al. [2004] Claudio Amovilli, Frederik E Leys, and Norman H March.  Electronic energy spectrum of two-dimensional solids and a chain of c atoms from a quantum network model.  _Journal of mathematical chemistry_ , 36(2):93–112, 2004.
  * Leys et al. [2004] Frederik E Leys, Claudio Amovilli, and Norman H March.  Topology, connectivity, and electronic structure of c and b cages and the corresponding nanotubes.  _Journal of chemical information and computer sciences_ , 44(1):122–135, 2004.
  * de Oliveira and Rocha [2022] César R de Oliveira and Vinícius L Rocha.  Bidimensional honeycomb materials: A graph model through dirac operator.  _Reports on Mathematical Physics_ , 89(2):231–252, 2022.
  * Fisher et al. [2021] Lee Fisher, Wei Li, and Stephen P Shipman.  Reducible fermi surface for multi-layer quantum graphs including stacked graphene.  _Communications in Mathematical Physics_ , 385(3):1499–1534, 2021.
  * Lawrie et al. [2022] Tristan Lawrie, Gregor Tanner, and Dimitrios Chronopoulos.  A quantum graph approach to metamaterial design.  _Scientific Reports_ , 12(1):1–14, 2022.
  * Arnold et al. [2016] Douglas N. Arnold, Guy David, David Jerison, Svitlana Mayboroda, and Marcel Filoche.  Effective Confining Potential of Quantum States in Disordered Media.  _Physical Review Letters_ , 116(5):056602, February 2016.  ISSN 0031-9007, 1079-7114.  doi:10.1103/PhysRevLett.116.056602.  URL https://link.aps.org/doi/10.1103/PhysRevLett.116.056602.
  * Filoche et al. [2017] Marcel Filoche, Marco Piccardo, Yuh-Renn Wu, Chi-Kang Li, Claude Weisbuch, and Svitlana Mayboroda.  Localization landscape theory of disorder in semiconductors. I. Theory and modeling.  _Physical Review B_ , 95(14):144204, April 2017\.  ISSN 2469-9950, 2469-9969.  doi:10.1103/PhysRevB.95.144204.  URL http://link.aps.org/doi/10.1103/PhysRevB.95.144204.
  * Piccardo et al. [2017] Marco Piccardo, Chi-Kang Li, Yuh-Renn Wu, James S. Speck, Bastien Bonef, Robert M. Farrell, Marcel Filoche, Lucio Martinelli, Jacques Peretti, and Claude Weisbuch.  Localization landscape theory of disorder in semiconductors. II. Urbach tails of disordered quantum well layers.  _Physical Review B_ , 95(14):144205, April 2017\.  ISSN 2469-9950, 2469-9969.  doi:10.1103/PhysRevB.95.144205.  URL http://link.aps.org/doi/10.1103/PhysRevB.95.144205.
  * Li et al. [2017] Chi-Kang Li, Marco Piccardo, Li-Shuo Lu, Svitlana Mayboroda, Lucio Martinelli, Jacques Peretti, James S. Speck, Claude Weisbuch, Marcel Filoche, and Yuh-Renn Wu.  Localization landscape theory of disorder in semiconductors. III. Application to carrier transport and recombination in light emitting diodes.  _Physical Review B_ , 95(14):144206, April 2017\.  ISSN 2469-9950, 2469-9969.  doi:10.1103/PhysRevB.95.144206.  URL http://link.aps.org/doi/10.1103/PhysRevB.95.144206.
  * Gebhard et al. [2022] F. Gebhard, A. V. Nenashev, K. Meerholz, and S. D. Baranovskii.  Quantum states in disordered media. i. low-pass filter approach, 2022\.
  * Harrell II and Maltsev [2019] Evans M. Harrell II and Anna V. Maltsev.  Localization and landscape functions on quantum graphs.  _Transactions of the American Mathematical Society_ , 373(3):1701–1729, November 2019.  ISSN 0002-9947, 1088-6850.  doi:10.1090/tran/7908.  URL https://www.ams.org/tran/2020-373-03/S0002-9947-2019-07908-4/.
  * Schanz and Smilansky [2000] Holger Schanz and Uzy Smilansky.  Periodic-orbit theory of anderson localization on graphs.  _Physical Review Letters_ , 84(7):1427, 2000.
  * Hislop and Post [2009] Peter D Hislop and Olaf Post.  Anderson localization for radial tree-like random quantum graphs.  _Waves in Random and Complex Media_ , 19(2):216–261, 2009.
  * Klopp and Pankrashkin [2009] Frédéric Klopp and Konstantin Pankrashkin.  Localization on quantum graphs with random edge lengths.  _Letters in Mathematical Physics_ , 87(1):99–114, 2009.
  * Sabri [2014] Mostafa Sabri.  Anderson localization for a multi-particle quantum graph.  _Reviews in Mathematical Physics_ , 26(01):1350020, 2014.
  * Damanik et al. [2020] David Damanik, Jake Fillman, and Selim Sukhtaiev.  Localization for anderson models on metric and discrete tree graphs.  _Mathematische Annalen_ , 376(3):1337–1393, 2020\.
  * Gaio et al. [2019] Michele Gaio, Dhruv Saxena, Jacopo Bertolotti, Dario Pisignano, Andrea Camposeo, and Riccardo Sapienza.  A nanophotonic laser on a graph.  _Nature Communications_ , 10(1), January 2019.  doi:10.1038/s41467-018-08132-7.  URL https://doi.org/10.1038/s41467-018-08132-7.
  * Brio et al. [2022] M. Brio, J. G. Caputo, and H. Kravitz.  Localized eigenvectors on metric graphs, 2022.
  * Suja et al. [2017] Mohammad Suja, Sunayna Binte Bashar, Bishwajit Debnath, Longxing Su, Wenhao Shi, Roger Lake, and Jianlin Liu.  Electrically driven deep ultraviolet MgZnO lasers at room temperature.  _Scientific Reports_ , 7(1), June 2017.  doi:10.1038/s41598-017-02791-0.  URL https://doi.org/10.1038/s41598-017-02791-0.
  * Jeong et al. [2022] Bongjin Jeong, Rashida Akter, Jeonghyun Oh, Dong-Gi Lee, Chang-Geun Ahn, Jong-Soon Choi, Md Rahman, et al.  Novel electrochemical pmi marker biosensor based on quantum dot dissolution using a double-label strategy.  _Scientific Reports_ , 12(1):1–8, 2022.
  * Kulce et al. [2021] Onur Kulce, Deniz Mengu, Yair Rivenson, and Aydogan Ozcan.  All-optical information-processing capacity of diffractive surfaces.  _Light: Science & Applications_, 10(1), January 2021.  doi:10.1038/s41377-020-00439-9.  URL https://doi.org/10.1038/s41377-020-00439-9.
  * Montbrió et al. [2015] Ernest Montbrió, Diego Pazó, and Alex Roxin.  Macroscopic description for networks of spiking neurons.  _Physical Review X_ , 5(2):021028, 2015.
  * Martinello et al. [2017] Matteo Martinello, Jorge Hidalgo, Amos Maritan, Serena Di Santo, Dietmar Plenz, and Miguel A Muñoz.  Neutral theory and scale-free neural dynamics.  _Physical Review X_ , 7(4):041071, 2017.
  * Kadmon and Sompolinsky [2015] Jonathan Kadmon and Haim Sompolinsky.  Transition to chaos in random neuronal networks.  _Physical Review X_ , 5(4):041030, 2015.
  * Esfandiary et al. [2020] Samaneh Esfandiary, Ali Safari, Jakob Renner, Paolo Moretti, and Miguel A Muñoz.  Anomalous lifshitz dimension in hierarchical networks of brain connectivity.  _Physical Review Research_ , 2(4):043291, 2020\.
  * Menzel et al. [2020] Miriam Menzel, Markus Axer, Hans De Raedt, Irene Costantini, Ludovico Silvestri, Francesco S Pavone, Katrin Amunts, and Kristel Michielsen.  Toward a high-resolution reconstruction of 3d nerve fiber architectures and crossings in the brain using light scattering measurements and finite-difference time-domain simulations.  _Physical Review X_ , 10(2):021002, 2020.
  * Elam et al. [2021] Jennifer Stine Elam, Matthew F Glasser, Michael P Harms, Stamatios N Sotiropoulos, Jesper LR Andersson, Gregory C Burgess, Sandra W Curtiss, Robert Oostenveld, Linda J Larson-Prior, Jan-Mathijs Schoffelen, et al.  The human connectome project: A retrospective.  _NeuroImage_ , 244:118543, 2021.
  * Rawlinson et al. [2021] Jonathan I Rawlinson, Csaba Fábri, and Attila G Császár.  Exactly solvable 1d model explains the low-energy vibrational level structure of protonated methane.  _Chemical Communications_ , 57(39):4827–4830, 2021.
  * Fábri et al. [2017] Csaba Fábri, Martin Quack, and Attila G Császár.  On the use of nonrigid-molecular symmetry in nuclear motion computations employing a discrete variable representation: A case study of the bending energy levels of ch 5+.  _The Journal of Chemical Physics_ , 147(13):134101, 2017.
  * Post [2012] Olaf Post.  _Spectral analysis on graph-like spaces_ , volume 2039.  Springer Science & Business Media, 2012.
  * Berkolaiko et al. [2008] G Berkolaiko, JM Harrison, and JH Wilson.  Mathematical aspects of vacuum energy on quantum graphs.  _Journal of Physics A: Mathematical and Theoretical_ , 42(2):025204, 2008.
  * Fulling et al. [2007] SA Fulling, L Kaplan, and JH Wilson.  Vacuum energy and repulsive casimir forces in quantum star graphs.  _Physical review A_ , 76(1):012118, 2007.
  * Harrison and Kirsten [2011] JM Harrison and K Kirsten.  Zeta functions of quantum graphs.  _Journal of Physics A: Mathematical and Theoretical_ , 44(23):235301, 2011.
  * Carlson [2008] Robert Carlson.  Boundary value problems for infinite metric graphs.  _Analysis on graphs and its applications_ , 77:355–368, 2008.
  * Carlson [2011] Robert Carlson.  Dirichlet to neumann maps for infinite quantum graphs.  _arXiv preprint arXiv:1109.3132_ , 2011.
  * Carlson [2020] Robert Carlson.  Metric graphs with totally disconnected boundary.  _arXiv preprint arXiv:2011.00145_ , 2020.
  * Carlson [2006] Robert Carlson.  Linear network models related to blood flow.  _Contemporary Mathematics_ , 415:65–80, 2006.
  * Maury et al. [2009] Bertrand Maury, Delphine Salort, and Christine Vannier.  Trace theorems for trees and application to the human lungs.  _Networks & Heterogeneous Media_, 4(3):469, 2009\.
  * Carlson [2013] Robert Carlson.  Myopic models of population dynamics on infinite networks.  _arXiv preprint arXiv:1309.5378_ , 2013.
  * Pokornyi and Borovskikh [2004] Yu V Pokornyi and AV Borovskikh.  Differential equations on networks (geometric graphs).  _Journal of Mathematical Sciences_ , 119(6):691–718, 2004.
  * Exner et al. [2008] Pavel Exner, Jonathan P Keating, Peter Kuchment, Alexander Teplyaev, and Toshikazu Sunada.  _Analysis on Graphs and Its Applications: Isaac Newton Institute for Mathematical Sciences, Cambridge, UK, January 8-June 29, 2007_ , volume 77.  American Mathematical Soc., 2008.
  * Berkolaiko et al. [2006] Gregory Berkolaiko et al.  _Quantum Graphs and Their Applications: Proceedings of an AMS-IMS-SIAM Joint Summer Research Conference on Quantum Graphs and Their Applications, June 19-23, 2005, Snowbird, Utah_ , volume 415.  American Mathematical Soc., 2006.
  * Mugnolo [2015] Delio Mugnolo.  _Mathematical Technology of Networks: Bielefeld, December 2013_ , volume 128.  Springer, 2015.
  * Exner and Kovařík [2015] Pavel Exner and Hynek Kovařík.  _Quantum waveguides_.  Springer, 2015.
  * Dáger and Zuazua [2006] René Dáger and Enrique Zuazua.  _Wave propagation, observation and control in 1-d flexible multi-structures_ , volume 50.  Springer Science & Business Media, 2006.
  * Mugnolo [2014] Delio Mugnolo.  _Semigroup methods for evolution equations on networks_ , volume 20.  Springer, 2014.
  * Friedman and Tillich [2004a] Joel Friedman and Jean-Pierre Tillich.  Calculus on graphs.  _arXiv preprint cs/0408028_ , 2004a.
  * Friedman and Tillich [2004b] Joel Friedman and Jean-Pierre Tillich.  Wave equations for graphs and the edge-based laplacian.  _Pacific Journal of Mathematics_ , 216(2):229–266, 2004b.
  * Cattaneo [1997] Carla Cattaneo.  The spectrum of the continuous laplacian on a graph.  _Monatshefte für Mathematik_ , 124(3):215–235, 1997.
  * Pankrashkin [2006] Konstantin Pankrashkin.  Spectra of schrödinger operators on equilateral quantum graphs.  _Letters in Mathematical Physics_ , 77(2):139–154, 2006.
  * Exner et al. [2018] Pavel Exner, Aleksey Kostenko, Mark Malamud, and Hagen Neidhardt.  Spectral theory of infinite quantum graphs.  In _Annales Henri Poincaré_ , volume 19, pages 3457–3510. Springer, 2018.
  * Exner et al. [2006] Pavel Exner, Pavel Hejčík, and Petr Šeba.  Approximations by graphs and emergence of global structures.  _Reports on Mathematical Physics_ , 57(3):445–455, June 2006.  doi:10.1016/s0034-4877(06)80031-5.  URL https://doi.org/10.1016/s0034-4877(06)80031-5.
  * Nakamura and Tadano [2021] Shu Nakamura and Yukihide Tadano.  On a continuum limit of discrete schrödinger operators on square lattice.  _Journal of Spectral Theory_ , 11(1):355–367, March 2021.  doi:10.4171/jst/343.  URL https://doi.org/10.4171/jst/343.
  * Exner et al. [2022] Pavel Exner, Shu Nakamura, and Yukihide Tadano.  Continuum limit of the lattice quantum graph hamiltonian.  _Letters in Mathematical Physics_ , 112(4), August 2022\.  doi:10.1007/s11005-022-01576-5.  URL https://doi.org/10.1007/s11005-022-01576-5.
  * van der Hoorn et al. [2021] Pim van der Hoorn, William J Cunningham, Gabor Lippner, Carlo Trugenberger, and Dmitri Krioukov.  Ollivier-ricci curvature convergence in random geometric graphs.  _Physical Review Research_ , 3(1):013211, 2021\.
  * Evans and Garzepy [2018] Lawrence C Evans and Ronald F Garzepy.  _Measure theory and fine properties of functions_.  Routledge, 2018.
  * Berry [1984] Michael Victor Berry.  Quantal phase factors accompanying adiabatic changes.  _Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences_ , 392(1802):45–57, 1984.
  * Girard [1987] D Girard.  Un algorithme simple et rapide pour la validation croisée généralisée sur des problèmes de grande taille.  _Informatique et Mathématiques Appliquées de Grenoble: Grenoble, France_ , 1987.
  * Hutchinson [1989] Michael F Hutchinson.  A stochastic estimator of the trace of the influence matrix for laplacian smoothing splines.  _Communications in Statistics-Simulation and Computation_ , 18(3):1059–1076, 1989.