\section{1. Extended proofs}
\label{sec-periodic}

\noindent
In this section we prove the theorem and the proposition from the main text. For the convenience of the reader, we briefly recall our framework.

The Hamiltonian $H$ is a bounded hermitian operator on  $\ell^2(\Gamma;\mathcal{H})$, where  $\Gamma \subset \mathbb{R}^d$ is a countable set and $\mathcal{H}$ a separable Hilbert space. We assume that there is a maximum distance $r > 0$ such that any point in $\mathbb{R}^d$ is at most at distance~$r$ away from a point in $\Gamma$. Recall that distances in $\mathbb{R}^d$ are measured in the norm  $\|x\|_\infty := \max\{ |x_1|,\ldots, |x_d|\}$. The open cube with side-length $2r$ around $x\in\mathbb{R}^d$ is
\[
B_r(x) := \{ y\in \mathbb{R}^d\,|\, \|x-y\|_\infty <r\}
\]
and its closure $\overline{B_r}(x) := \{ y\in \mathbb{R}^d\,|\, \|x-y\|_\infty \leq r\}$.
For any set $A\subset \mathbb{R}^d$ and wave function $\psi\in \ell^2(\Gamma;\mathcal{H})$, we denote by  $\norm{\psi}_A^2 := \sum_{x\in \Gamma \cap A} \|\psi(x)\|_{\mathcal{H}}^2$ its $\ell^2$-mass  in the region~$A$.

The following definition makes precise the idea that a Hamiltonian is gapped in the bulk at a certain energy $\lambda$ and a certain length scale $L$,  if all $\epsilon$-quasimodes at that energy that are supported in a region of that scale \eqref{eq-prop} are disproportionately supported near the edge of the region \eqref{sup:eq-prop-conclusio}.

\begin{definition}
	\label{sup:def-eps-bulk-gapped}
	Let $\epsilon,L>0$ and $\lambda\in\mathbb{R}$. We say that a Hamiltonian $H$ is \textbf{locally $\epsilon$-bulk-gapped at energy $\lambda$ and scale $L$}, if there exist constants  $N \in \mathbb{N}$, $N\geq 2$ and $C<\frac{1}{N^d} $, such that for any $x\in\Gamma$ and for any $\psi \in \ell^2(\Gamma;\mathcal{H})$, we have the following implication: Whenever
	\begin{align}
	\norm{ (H - \lambda) \psi}_ {B_L(x)} \leq \epsilon \norm{\psi}_{B_{L}(x)}\,,
	\label{eq-prop}
	\end{align}
then for $l := \frac{L + r}{N} + r$ it holds that
\begin{align}
\label{sup:eq-prop-conclusio}
\norm{\psi}_{\overline{B_l}(x)}^2 \leq C \norm{\psi}_{B_L(x)}^2\,.
\end{align}
\end{definition}

In the following we call any   $\psi \in L^2(\mathbb{R}^d)$ that satisfies the inequality \eqref{eq-prop}   an \textbf{ $L$-local $\epsilon$-quasimode at $x$}. In these terms, Definition~\ref{sup:def-eps-bulk-gapped} says that every $L$-local $\epsilon$-quasimode at $x$ has most of its mass in $B_L(x)$ concentrated outside of $\overline{B_l}(x)$. This terminology is useful in the proof of the main theorem.

\begin{theorem}\label{main-theorem-recall}
If a Hamiltonian $H$ on $\ell^2(\Gamma; \mathcal{H})$ is locally $\epsilon$-bulk-gapped at energy $\lambda$ on some scale $L>0$, then the interval  $(\lambda - \epsilon, \lambda + \epsilon)$ is a gap in the spectrum of $H$, i.e.\
\begin{align*}
\sigma(H) \cap (\lambda - \epsilon, \lambda + \epsilon) = \varnothing.
\end{align*}
\end{theorem}

\begin{proof}
Assume that $H$ is $\epsilon$-bulk-gapped at energy $\lambda$ and scale $L>0$. Let $  l, N$ and $C$ be as in Definition~\ref{sup:def-eps-bulk-gapped}.
We generalise the proof strategy presented in the main text.
For every $q \in \{1, \dots, N\}^d$ let
\begin{align}
\widetilde Z_q := 2(L + r)(\mathbb{Z}^{d} + q/N)\label{eq-def-tilde-Zq}
\end{align}
and
\begin{align}
f_q: \widetilde  Z_q \to \Gamma\,,\; \widetilde x\mapsto f_q(\widetilde x) := \text{the point in $\Gamma$ closest to $\tilde x$}\,.
\end{align}
In the case where several   points  in $\Gamma$ minimize the distance to $\widetilde x$, we   arbitrarily  pick one of them in the definition of $f_q(\widetilde x)$.
By the assumption on $\Gamma$ we have that $\|\widetilde x- f_q(\widetilde x)\|_\infty \leq  r$ for all $\widetilde x\in\widetilde Z_q$. Since two distinct points in $\widetilde Z_q$ have at least distance $2(L + r)$, this implies that the map $f_q$ is one-to-one. Hence, as a map onto its image
\begin{align}
Z_q :=  f_q(\widetilde Z_q) \subset \Gamma\,,
\label{eq-def-Zq}
\end{align}
the map $f_q: \widetilde  Z_q \to Z_q$ is a bijection and we think of $Z_q$ as a deformation of~$\widetilde Z_q$.

Is is now straightforward to see that the large boxes $B_L(x)$ are still mutually disjoint when $x$ varies in one of the deformed sublattices $Z_q$ and that the small boxes  $\overline{B_l}(x)$ still cover all of $\mathbb{R}^d$  if $x$ varies in the union $\cup_q Z_q$:\\[2mm]
\noindent\textit{Disjointness:} For any $q \in \{1, \dots, N\}^d$ and two different points $x, y \in Z_q$
\begin{align}
B_L(x) \cap B_L(y) = \varnothing\,.
\label{eq-intersections-disjoint}
\end{align}
\begin{proof}Let $\widetilde x := f_q^{-1}(x)$ and $\widetilde y:= f_q^{-1}(y)$. Then $\widetilde x$ and $\widetilde y$ are two distinct points in the square lattice $\widetilde Z_q$ and thus their distance is at least $2(L + r)$.  Using the inverse triangle inequality we conclude that
\[
\|x-y\|_\infty = \|\tilde x -\tilde y + x-\tilde x +  \tilde y - y\|_\infty \geq  \|\tilde x -\tilde y\|_\infty - \|x-\tilde x\|_\infty -\|  \tilde y - y\|_\infty  \geq  2(L+r) - 2r=2L\,.
\]
Thus, the boxes $B_L(x)$ and  $B_L(y)$ do not overlap.
\end{proof}

\noindent \textit{Covering:}
\begin{align}
\bigcup_{q \in \{1, \dots, N\}^d} \bigcup_{x \in Z_q} \overline{B_l}(x) = \mathbb{R}^d.
\label{eq-double-union}
\end{align}
\begin{proof}
Note that the union
\begin{align*}
\widetilde Z  := \bigcup_{q \in \{1, \dots, N\}^d} \widetilde Z_q
\end{align*}
is a square lattice with side length $2(L + r) / N$. Thus, for every $p \in \mathbb{R}^d$, there exists a $q  \in \{1, \dots, N\}^d$ and  $\widetilde x \in \widetilde Z_q$ such that $\|p- \widetilde x\|_\infty \leq (L + r) / N$.
The triangle inequality now implies
\begin{align}
\| p- f_q(\widetilde x) \|_\infty \leq \|p-  \widetilde x  \|_\infty + \|\widetilde x- f_q(\widetilde x) \|_\infty \leq (L + r) / N + r=l\,
\end{align}
and therefore  $p \in \overline{B_l}(x)$ for $x := f_q(\widetilde x) \in Z_q$.
\end{proof}

 These two properties allow us to prove by contradiction that the spectrum $\sigma(H)$ does not contain the interval $(\lambda - \epsilon, \lambda + \epsilon)$.

  Suppose there exists   $\nu \in (\lambda - \epsilon, \lambda + \epsilon)\cap\sigma(H) $. Then, according to Weyl's criterion,  there exist arbitrarily precise quasimodes for the energy $\nu$. More precisely, for any $\delta>0$  there exists $\psi \in \ell^2(\Gamma, \mathcal{H})$   such that
\begin{align} \label{eq-quasimode}
\norm{\psi} = 1 \quad\mbox{and}\quad
\| (H - \nu) \psi \|  < \delta .
\end{align}
We now fix $\delta > 0$ such that
\begin{subequations}
\begin{gather}
\delta < (\epsilon - |\lambda - \nu|)^2 \label{eq-delta-bound-a}\\
\delta < N^{-d} - C \label{eq-delta-bound-b}
\end{gather}
\end{subequations}
and  choose a corresponding $\psi \in \ell^2(\Gamma, \mathcal{H})$ that satisfies \eqref{eq-quasimode}.

Notice that if $\psi$ is a $\delta$-quasimode on $\mathbb{R}^d$, this does not imply that $\psi$ is an $L$-local $\delta$-quasimode for all $x \in \mathbb{R}^d$. However, we have that the total mass of $\psi$ on those squares $B_L(x)$ for which $\psi$ is not an $L$-local $\delta$-quasimode is small (it is of order $\delta$). To see this, we split each $Z_q$ into two subsets
\begin{align*}
Z_q^+ := \big\{  x \in  Z_q  \, |\, \norm{(H - \nu) \psi}_{B_L(x)} \leq \sqrt{\delta} \norm{\psi}_{B_L(x)} \big\}\,,\\
Z_q^- := \big\{  x \in Z_q  \, |\, \norm{(H - \nu) \psi}_{B_L(x)} > \sqrt{\delta} \norm{\psi}_{B_L(x)}\big\}\,.
\end{align*}
Around the points $x \in Z_q^+$, $\psi$ is an $L$-local $\epsilon$-quasimode, since
\begin{align} \label{eq:LEpsQuasimode}
 \norm{(H - \lambda)\psi }_{B_L(x)}
 &\leq \|(H - \nu) \psi \|_{B_L(x)} \hspace{-3pt}+ |\lambda - \nu| \norm{\psi}_{B_L(x)} \nonumber
 \\ &\leq (\sqrt{\delta} + |\lambda - \nu|) \norm{\psi}_{B_L(x)}\nonumber
 \\ & < \epsilon \norm{\psi }_{B_L(x)}\,.
\end{align}
For the last step of the above, we rewrite \eqref{eq-delta-bound-a} as
\begin{align*}
\delta < (\epsilon - | \lambda - \nu|)^2 \quad \Rightarrow \quad \sqrt{\delta} < \epsilon - |\lambda - \nu |
\quad \Rightarrow \quad
\sqrt{\delta} + | \lambda - \nu | < \epsilon\,.
\end{align*}
Since  $H$ is $\epsilon$-bulk-gapped on the scale $L$, \eqref{eq:LEpsQuasimode} implies that for every $x \in Z^+_q$ the inequality \eqref{sup:eq-prop-conclusio} holds true. Summing over all $x \in Z_q^+$, we get
\begin{align}
\sum_{x \in Z_q^+} \norm{\psi}_{\overline{B_l}(x)}^2 \leq C \sum_{x \in Z_q^+} \norm{\psi}_{B_L(x)}^2\,.
\label{sum-Zplus}
\end{align}
By the disjointness condition shown before, we have that $B_L(x)$ and $B_L(y)$ are disjoint for different $x$ and $y$ in $Z^+_q$.
Since for disjoint sets $A$ and $B$, we have
$
\norm{\psi}_A^2 + \norm{\psi}_B^2 = \norm{\psi}_{A \cup B}^2
$, we obtain
\begin{equation}
\label{eq:aux1}
\sum_{x \in Z_q^+} \norm{\psi}_{B_L(x)}^2 = \norm{\psi}_U^2 \leq \norm{\psi}^2 = 1\quad\text{for } U := \bigcup_{x \in Z_q^+} B_L(x) .
\end{equation}
Combining \eqref{eq:aux1} with \eqref{sum-Zplus}, we get
\begin{align}
\sum_{x \in Z_q^+} \norm{\psi}_{\overline{B_l}(x)}^2 \leq C.
\label{eq-small-Zplus-sum}
\end{align}
Let now $x \in Z_q^-$, then  by definition of  $Z_q^-$ we have \begin{equation}
\label{eq:aux2}
\norm{\psi}^2_{B_L(x)} < \frac{1}{\delta}\norm{(H - \nu) \psi}^2_{B_L(x)}  \, .
\end{equation}
By taking the sum over all $x \in Z_q^-$ we get
\begin{align*}
\sum_{x \in Z_q^-} \norm{\psi}^2_{\overline{B_l}(x)} \leq \sum_{x \in Z_q^-}  \norm{\psi}^2_{B_L(x)} < \frac1\delta \sum_{x \in Z_q^-} \norm{(H - \nu) \psi}_{B_L(x)}^2 \leq \frac1\delta \norm{(H - \nu) \psi}^2 \,,
\end{align*}
  where we used the disjointness of the $B_L(z)$ for $z \in Z_q$ in the last step. Since $\psi$ is a $\delta$-quasimode, meaning $\norm{(H - \nu) \psi}^2 < \delta^2$, we get
\begin{align}
\sum_{x \in Z_q^-} \norm{\psi}^2_{\overline{B_l}(x)} < \delta \, .
\label{eq-small-Zminus-sum}
\end{align}
 Combining  \eqref{eq-small-Zplus-sum} and \eqref{eq-small-Zminus-sum} and using the covering property \eqref{eq-double-union}, we finally obtain
\begin{align*}
\norm{\psi}^2 &\leq \sum_{q \in \{ 1, \dots , N\}^d} \sum_{x \in Z_q} \norm{\psi}^2_{\overline{B_l}(x)} = \sum_{q \in \{ 1, \dots , N\}^d} \left( \sum_{x \in Z_q^+} \norm{\psi}^2_{\overline{B_l}(x)} + \sum_{x \in Z_q^-} \norm{\psi}^2_{\overline{B_l}(x)}  \right)\\
&\leq  \sum_{q \in \{ 1, \dots , N\}^d} (C + \delta) \;= \;N^d (C + \delta)  \; <\;1\,,
\end{align*}
where in the last inequality we used the hypothesis \eqref{eq-delta-bound-b} on $\delta$. Since $\norm{\psi}^2 < 1$ contradicts the  normalisation of $\psi$ assumed in \eqref{eq-quasimode}, such  a $\delta$-quasimode   cannot exist,  and $\nu$ is not in the spectrum of $H$.
\end{proof}

 \begin{remark} Notice that neither Definition \ref{sup:def-eps-bulk-gapped} nor Proposition \ref{sup:eq-prop-conclusio} require the boundedness of the Hamiltonian and can be easily generalized to unbounded operators taking into account only vectors that belong to the domain of the Hamiltonian.
\end{remark}

Next we prove the  numerically verifiable criterion for the $\epsilon$-bulk-gapped property. Recall that for any set $A\subset \mathbb{R}^d$ we denote by $\mathbf{1}_A$ the characteristic function of $A$ and we use the shorthand notation
$H_A :=\mathbf{1}_A H \mathbf{1}_A$.

\begin{proposition}
\label{proposition-recall}
In addition to the assumption of Theorem~\ref{main-theorem-recall},
 let $H$ have finite range, i.e.\ there exists a maximal hopping distance $m$ such that
\begin{align*}
H_{xy} = 0\quad\text{for all\; $x,y\in\Gamma$ \;with\; $\norm{x-y}_\infty > m$}.
\end{align*}
Let $L > 0$, $N \in \mathbb{N}$, $N\geq2$, $l := \frac{L + r}{N} + r$, $\lambda\in\mathbb{R}$.
Assume that for every $x \in \Gamma$   there exists   a set  $A(x)\subset \Gamma$ such that
$\overline{B_l}(x)\,\subseteq\, A(x)\subseteq B_{L-m}(x)$,  $\lambda \notin \sigma(H_{A(x)})$, and
\begin{align*}
D(x) =\norm{\mathbf{1}_{\overline{B_l}(x)} (H_{A(x)} - \lambda)^{-1}  \mathbf{1}_{A(x)} H \mathbf{1}_{B_L(x) \backslash A(x)}}_\mathrm{op}
\end{align*}
satisfies $D(x)<N^{-d/2}$.
Then $H$ is $\epsilon$-bulk-gapped at energy $\lambda$ and scale $L$ for any~$\epsilon>0$ with
\begin{align}
\epsilon < \inf_{x\in\Gamma} \frac{N^{-d/2} - D(x)}{\norm{\mathbf{1}_{\overline{B_l}(x)} (H_{A(x)}- \lambda)^{-1} \mathbf{1}_{A(x)}}_\mathrm{op}}\,.
\label{eq-conditions-epsilon-recall}
\end{align}
\end{proposition}

\begin{proof}
 Let $x \in \Gamma$ and suppose that for some $\psi  \in\ell^2(\Gamma;\mathcal{H}) $ property \eqref{eq-prop} holds, i.e.\
 \[
 u := \mathbf{1}_{B_L(x)} (H - \lambda) \psi
 \]
 satisfies
 \begin{equation}
 \norm{u} < \epsilon \norm{\psi}_{B_L(x)}\,,\label{eq-u-bound}
\end{equation}
with $\epsilon>0 $ satisfying \eqref{eq-conditions-epsilon-recall}. We need to show that this implies \eqref{sup:eq-prop-conclusio}.
Writing
\begin{align*}
(H_{A(x)} - \lambda) \mathbf{1}_{A(x)} \psi &= \mathbf{1}_{A(x)} (H - \lambda) \mathbf{1}_{A(x)} \psi\\
&= \mathbf{1}_{A(x)} (H - \lambda) \psi - \mathbf{1}_{A(x)} (H - \lambda) \mathbf{1}_{{A(x)}^c} \psi\\
&= \mathbf{1}_{A(x)} u -\mathbf{1}_{A(x)} H \mathbf{1}_{{A(x)}^c} \psi
\end{align*}
and multiplying this equality by  $(H_{A(x)} - \lambda)^{-1}$  gives
\begin{align*}
\mathbf{1}_{A(x)} \psi = (H_{A(x)} - \lambda)^{-1} (\mathbf{1}_{A(x)} u + \mathbf{1}_{A(x)} H \mathbf{1}_{{A(x)}^c} \psi)\,.
\end{align*}
Since we need to estimate $\norm{\psi}_{\overline{B_l}(x)}$, we can multiply by $\mathbf{1}_{B_l(x)}$ to obtain
\begin{align*}
\mathbf{1}_{B_l(x)}\psi &= \mathbf{1}_{\overline{B_l}(x)} (H_{A(x)} - \lambda)^{-1} (\mathbf{1}_{A(x)} u + \mathbf{1}_{A(x)} H \mathbf{1}_{{A(x)}^c} \psi)\,.
\end{align*}
Using the triangle inequality, we obtain
\begin{align}
\norm{\psi}_{\overline{B_l}(x)} &\leq \norm{\mathbf{1}_{\overline{B_l}(x)} (H_{A(x)} - \lambda)^{-1} \mathbf{1}_{A(x)} u} + \norm{\mathbf{1}_{\overline{B_l}(x)} (H_{A(x)} - \lambda)^{-1} \mathbf{1}_{A(x)} H \mathbf{1}_{{A(x)}^c} \psi} \,.
\label{eq-after-triangle}
\end{align}
The first term is easily bounded using  \eqref{eq-u-bound},
\begin{align*}
\norm{\mathbf{1}_{\overline{B_l}(x)} (H_{A(x)} - \lambda)^{-1} \mathbf{1}_{A(x)} u} < \norm{\mathbf{1}_{\overline{B_l}(x)} (H_{A(x)}- \lambda)^{-1} \mathbf{1}_{A(x)}}_\mathrm{op} \epsilon \norm{\psi}_{B_L(x)} =: M(x) \;\epsilon \norm{\psi}_{B_L(x)} \,.
\end{align*}
  Because of the finite range hypothesis on $H$ and the fact that ${A(x)} \subseteq B_{L-m}(x)$, we can rewrite the second term   as
$$
\mathbf{1}_{\overline{B_l}(x)} (H_{A(x)} - \lambda)^{-1} \mathbf{1}_{A(x)} H \mathbf{1}_{{A(x)}^c} \psi= \mathbf{1}_{\overline{B_l}(x)} (H_{A(x)} - \lambda)^{-1} \mathbf{1}_{A(x)} H \mathbf{1}_{B_L(x) \setminus {A(x)}} \psi.
$$
Thus, the second term in \eqref{eq-after-triangle} can be estimated using the assumption on $D(x)$ and we find that
\begin{align*}
\norm{\psi}_{\overline{B_l}(x)} < \epsilon M(x) \, \norm{\psi}_{B_L(x)} + D(x) \norm{\psi}_{B_L(x) \setminus {A(x)}} \leq (D(x) + \epsilon M(x)) \norm{\psi}_{B_L(x)}\,.
\end{align*}
Since   \eqref{eq-conditions-epsilon-recall} implies that
\[
\sup_{x\in\Gamma} (D(x) + \epsilon M(x))^2 < \sup_{x\in\Gamma} \left(D(x) + \frac{N^{-d/2} - D(x)}{M(x)} M(x)\right)^2   = N^{-d}\,,
\]
  \eqref{sup:eq-prop-conclusio} holds for any
 $C$ with $\sup_{x\in\Gamma} (D(x) + \epsilon M(x))^2 < C < N^{-d}$.
\end{proof}

\section{2. Description of the algorithms \label{sec-algorithms}}
\label{sm-sec:DescriptionAlg}
\noindent
In this section we explain in detail how to apply the general results of Theorem~\ref{main-theorem-recall}   to Hamiltonians modelled over the
Ammann-Beenker tiling by enumerating all local patches on a certain scale $L$ and then checking the conditions of Proposition~\ref{proposition-recall}. The method can be easily generalised to other cut-and-project tilings. Moreover, Proposition~\ref{proposition-recall} can be applied in principle to any Hamiltonian that has finite local   complexity. The problem of enumerating all possible local restrictions of $H_{B_L(x)}$ on a scale $L$ then needs to be solved in a way specific to the structure of $\Gamma$ and $H$.

We first describe an algorithm to determine all the different local patches that can occur the Ammann-Beenker quasicrystal $\Gamma_\mathrm{AB}$ defined with the cut-and-project method, that is to determine  the set   $\mathcal{C}_L:=
\{\mathcal{C}_L(x)\,|\, x\in\Gamma_\mathrm{AB}\}$ of all local patches
$
\mathcal{C}_L(x) := \{ y\in \mathbb{R}^d \,|\, x+ y \in \Gamma_\mathrm{AB}\cap B_L(x)\}
$.

Recall that the vertices of the Ammann-Beenker tiling are defined as the set
\begin{align*}
\Gamma_\mathrm{AB} = \big\{\,p(z) \;\big|\; z \in \mathbb{Z}^4,\,\kappa(z) \in R\,\big\}\,,
\end{align*}
where  $a=\frac{1}{\sqrt{2}}$, $R\subset \mathbb{R}^2$ is the regular axis-aligned octagon centered at $0$ with side length $1$, and
\begin{align*}
p = \begin{pmatrix}
1 & a & 0 & - a\\
0 & a & 1 & a
\end{pmatrix} \,, \quad
\kappa = \begin{pmatrix}
1 & -a & 0 & a\\
0 & a & -1 & a
\end{pmatrix}
\end{align*}
are the two ``projections'' as maps from $\mathbb{R}^4$ to $\mathbb{R}^2$.

The first algorithm we describe determines the set $V_L \subseteq \mathbb{Z}^4$ of ``candidate points'' defined as
\begin{align*}
V_L = \{\, v \in \mathbb{Z}^4 \;|\; p(v) \in  B_L(0) \;\text{and}\; \kappa(v) \in 2 R \,\}.
\end{align*}
As explained in the main text, this set is defined such that for every $z \in \mathbb{Z}^4$, the points of $\Gamma_\mathrm{AB} \cap B_L(p(z))$ are all of the form $p(z + v)$ for some $v \in V_L$. Thus to determine what points are part of the local patch around $p(z)$, we only need to check the points $z + v$ for $v \in V_L$.

The set  of candidate points $V_L$ only depends on $L$ and hence we only have to be compute it once at the beginning of the algorithm. To compute it, it would be possible in principle to simply check the two conditions $p(v) \in B_L(0)$ and $\kappa(v) \in 2R$ for all integer points in a four-dimensional cube around $0$ with radius $2L$, say, but this would be very inefficient as it requires checking $O(L^4)$ points. Because the condition $\kappa(v) \in R$ means that all points in $V_L$ lie close to the two-dimensional hyperplane defined by $\kappa(v) = 0$, it should only be necessary to check the conditions for $O(L^2)$~points.

To compute $V_L$ efficiently, consider the matrix
\begin{align*}
t = \begin{pmatrix}
1 & 0 \\
a & a \\
0 & 1 \\
-a & a
\end{pmatrix}\,.
\end{align*}
The matrix $t$ satisfies $\kappa t = 0$ and $p t = 2 \cdot \mathbf{1}$. Because the columns of $t$ together with the second and fourth canonical basis vectors $e_2$ and $e_4$ form a linear basis of $\mathbb{R}^4$, we can write any $v \in V_L$ as
\begin{align}
v = t\begin{pmatrix}w_1 \\ w_2\end{pmatrix} + q_1 e_2 + q_2 e_4 = \begin{pmatrix} w_1 \\ a(w_1+w_2) +q_1\\ w_2 \\ a(w_2-w_1) +q_2
\end{pmatrix}  \,. \label{eq-v-def}
\end{align}
Because $v \in \mathbb{Z}^4$, also   $w_1,w_2\in\mathbb{Z}$.
The condition $\kappa v \in 2 R$ that holds for all $v \in V_L$ can be used to bound $q_1$ and $q_2$. In fact, using $\kappa t = 0$  we conclude that
\begin{align}
\begin{pmatrix}a & a \\ -a & a\end{pmatrix} \begin{pmatrix}q_1 \\ q_2\end{pmatrix} = \kappa(q_1 e_2 + q_2 e_4) = \kappa(v) \in 2R\,.
\label{cond-q1-q2}
\end{align}
Since $\begin{pmatrix}a & a \\ -a & a\end{pmatrix}$ and hence also its inverse are symmetries of $R$,  \eqref{cond-q1-q2} just becomes
\begin{align*}
\begin{pmatrix}q_1 \\ q_2\end{pmatrix} \in 2R\,.
\end{align*}
In particular, this condition on $(q_1,q_2)$ is independent of $L$. This is what allows enumerating all points according to equation~\eqref{eq-v-def} with $q_1, q_2$ in this range to be $O(L^2)$.

We now compute the range of values that $w_1$ and $w_2$ can take such that  \eqref{eq-v-def}  defines  an element $v \in V_L$. Using $pt = 2\cdot\mathbf{1}$ we find
\begin{align*}
p(v) = p \left( t \begin{pmatrix}w_1\\w_2\end{pmatrix} + q_1 e_2 + q_2 e_4\right) = 2 \begin{pmatrix}w_1 \\ w_2\end{pmatrix} + \begin{pmatrix}a & -a \\ a & a\end{pmatrix} \begin{pmatrix}q_1 \\ q_2\end{pmatrix}\,.
\end{align*}
As also $\begin{pmatrix}a & -a \\ a & a\end{pmatrix}$  leaves $R$ invariant, $(q_1, q_2) \in 2R$ implies
\begin{align*}
p(v) - 2 \begin{pmatrix}w_1 \\ w_2\end{pmatrix} \in 2 R\,.
\end{align*}
Since   $v\in V_L$ requires  $pv \in B_L(0)$ and since the $\infty$-norm of $2R$ is bounded by $1 + \sqrt{2}$, we find that
\begin{align*}
2 \begin{pmatrix}w_1 \\ w_2\end{pmatrix} \in B_{L + 1 + \sqrt{2}}(0)\,.
\end{align*}
Thus $w_1$ and $w_2$ must both lie in the interval $(-L/2 - s, L/2 + s)$ with $s=\frac{1+\sqrt{2}}{2}$.

Thus, when determining whether $v\in\mathbb{Z}^4$ lies in $V_L$,  the values of $ v_1=w_1$ and $ v_3=w_2$ can be restricted to  integers in this range (line 1 in Algorithm 1 below). For $q_1$ and $q_2$, it suffices to consider all values in $2R\subset B_{2s}(0)$ such that the resulting values for $v_2$ and $v_4$ become integral (line 2 in Algorithm 1 below).

\begin{algorithm}[H]
\caption{Enumerate ``candidate set'' $V_L$
\label{algorithm-enumerate-candidates}}
\begin{algorithmic}[1]
\ForAll{integers $v_1, v_3$   in $(-L/2 - s,  L/2 + s)$}
\ForAll{integers $v_2 \in a(v_1 + v_3) + [-2s, 2s]$ and $v_4 \in a(v_1 - v_3) + [-2s, 2s]$}
\If{$p(v) \in B_L(0)$ and $\kappa(v) \in 2 R$}
\State add $v$ to $V_L$
\EndIf
\EndFor
\EndFor
\end{algorithmic}
\end{algorithm}
Next, we describe the algorithm to enumerate the local patches of the Ammann-Beenker tiling. As we have discussed in the main text, the local patch $\mathcal{C}_L(x)$ centered at $x=p(z)\in \Gamma_\mathrm{AB}$ can be characterized by which of the candidate points in $V_L$  ``become part of the tiling''. Namely, for every $v \in V_L$, we have
\begin{align}
p(z + v) \in \Gamma_\mathrm{AB} \quad\Leftrightarrow\quad \kappa(z + v) \in R \quad\Leftrightarrow\quad \kappa(z) \in R - \kappa(v)\,.
\label{eq-equivalences}
\end{align}
Defining for $v\in\mathbb{R}^4$ the following decomposition of the acceptance region,
\begin{gather*}
 P^1(v) := R \cap  (R - \kappa(v))\quad\text{and}\quad P^0(v) := R \setminus P^1(v)\,,
\end{gather*}
   \eqref{eq-equivalences} entails that  for every $v \in V_L$ it holds that
   \begin{equation}\label{condition1}
   p(v) \in \mathcal{C}_L(x) \quad\Leftrightarrow\quad
   \kappa(z)\in P^1(v)\,.
   \end{equation}
By labelling the points in the set $V_L$ with an index $i \in \{1,\dots, |V_L|\}$, we can uniquely associate to each bit string $b=(b_1,\dots,b_{|V_L|}) \in \{0,1\}^{|V_L|}$ a local patch
\begin{align*}
\mathcal{C}_{L,b} := \{ \,p(v_i) \, | \,    b_i = 1, i \in \{1,\dots, |V_L|\}\}\,.
\end{align*}
From \eqref{condition1} we conclude that
\begin{equation}
\mathcal{C}_L(x) = \mathcal{C}_{L,b} \quad\Leftrightarrow
\quad \kappa(z) \in P_{(b_i)} := \bigcap_{i = 1}^{|V_L|} P^{b_i}(v_i)\,.\label{eq-def-A}
\end{equation}

For many bit strings $b$ the set $P_{(b_i)}$ turns out to be empty and thus not all local patches $\mathcal{C}_{L,b}$ defined by bit strings $b$ of length $|V_L|$ correspond to actual local patches of the Ammann-Beenker tiling. In fact, the number of local patches (also referred to as the ``patch counting function'' or, in \cite{julien2010}, ``complexity'') of the Ammann-Beenker tiling is of order $O(L^2)$, cf.\  \cite{julien2010}, while the number of possible bit strings is $O(2^{|V_L|})$, with $|V_L|$ growing like $L^2$.

To enumerate those bit strings that correspond to actual local patches (i.e.,~for which the set $P_{(b_i)}$ is not empty), we can use a dynamic programming approach. To do so, we first extend the definition \eqref{eq-def-A} of $P_{(b_i)}$ to shorter bit strings $(b_i)$, namely  to $(b_i)\in \{0,1\}^m$, $1 \leq m \leq |V_L|$, by taking into account only the intersections up to the $m$-th place of the bit string.

Our algorithm then proceeds step by step by computing all nonempty $P_{(b_i)}$ for bit strings $(b_i)$ of length $1$, of length~$2$, and so on. By \eqref{eq-def-A}, we can go from a bit string of lenght $n$ to a bit string of lenght $n+1$ using the following recursion relation:
\begin{equation}
\begin{aligned}
P_{(b_i)\oplus 1} &= P_{(b_i)} \cap P^1(v_{i+1})\\
P_{(b_i)\oplus 0} &= P_{(b_i)} \cap P^0(v_{i+1})\,.
\end{aligned}
\label{P-split}
\end{equation}

Suppose we have computed a set $J_n$ of all bit strings $(b_i)$ of length $n$ and their associated sets $P_{(b_i)}$ for all $(b_i)$ where $P_{(b_i)}$ is not empty. We can then compute the set $J_{n+1}$ by computing, for every $(b_i) \in J_n$, the two intersections on the right hand sides of \eqref{P-split} and adding those of the two bit strings $(b_i) \oplus 0$ and $(b_i) \oplus 1$ for which the intersection is not zero. This procedure can be implemented algorithmically by simply maintaining a list $J$ of nonempty regions $P_{(b_i)}$ and the associated bit strings, as described in the pseudocode in the loop in line \ref{optimization-loop-to-avoid} of Algorithm \ref{algorithm-enumerate-tilings}. The following algorithm is formulated to decompose an area $R_0 \subseteq R$. While setting $R_0 = R$ will lead to an enumeration of all patches, it is sometimes advantageous to compute only the decompositions corresponding to $\kappa(z) \in R_0$ for some smaller $R_0$, as described below.

\begin{algorithm}[H]
\caption{ Enumerate the local patches $\mathcal{C}_L(x)$}
\label{algorithm-enumerate-tilings}
\begin{algorithmic}[1]
\State Initialize $J = \{(R_0, \text{` '})\}$.
\State Compute $V_L$ using Algorithm~\ref{algorithm-enumerate-candidates}
\ForAll{$v \in V_L$}
\State set $P^1(v) = R - \kappa(v)$.
\State Initialize $J_2 = \{\}$
\If{$R_0 \subseteq P^1(v)$}\label{optimization-start}
\State set $J \leftarrow \{(Q, s \oplus\text{`1'}) \,|\, (Q,s) \in J\}$
\State \textbf{continue} with next $v$
\EndIf
\If{$R_0 \cap P^1(v) = \varnothing$}
\State set $J \leftarrow \{(Q, s \oplus\text{`0'}) \,|\, (Q,s) \in J\}$
\State \textbf{continue} with next $v$
\EndIf\label{optimization-end}
\ForAll{$(Q, s) \in J$}\label{optimization-loop-to-avoid}
\If{$P^1(v) \cap Q \neq \varnothing$}
\State Add $(P \cap Q, s \oplus \text{`1'})$ to $J_{ 2}$
\EndIf
\If{$Q \backslash P^1(v) \neq \varnothing$}
\State Add $(Q \backslash P, s \oplus \text{`0'})$ to $J_{2}$
\EndIf
\EndFor
\State Update $J \leftarrow J_2$.
\EndFor
\end{algorithmic}
\end{algorithm}

Usually we decompose the entire region of acceptance $R$ according to the previous procedure. However, it can be useful to decompose only a smaller polygon $R_0 \subseteq R$. The reason for this is twofold. First, there are some symmetries of $R$ which can allow us to decompose a smaller region. For example, the tiling corresponding to $\kappa(z) = (k_1, k_2)$ is exactly the mirror image of the tiling corresponding to $(k_2, k_1)$. Therefore, considering
\begin{align*}
R_0 = R \cap \{ (x,y) \in \mathbb{R}^2 \,|\, x > 0, y > 0, y < x \}.
\end{align*}
is sufficient for enumerating all local patches up to all mirror symmetries. In that case one may replace   $R$   by $R_0$ in Algorithm~\ref{algorithm-enumerate-tilings}.

Apart from the case where one has to check only part of the acceptance region $R$ for symmetry reasons, splitting $R$ into sub-polygons can improve the performance of the algorithm. In the default case, for every candidate point $v_i \in V_L$, the intersection of the set $P^1(v_i)$ with all polygons $P_{(b_j)}$ distinguished up to that point has to be computed. The time taken for this appears to be quadratic in the final number $|\mathcal{C}_L|$ of local patches distinguished. If one only computes the decomposition of a smaller polygon $R' \subseteq R$, the intersection computation can be skipped in all cases where either $R' \subseteq P^1(v_i)$ or $R' \subseteq P^0(v_i)$, because the intersections are trivial (either empty or all of $P_{(b_j)}$). This is implemented in lines~\ref{optimization-start} to~\ref{optimization-end} of the above algorithm by replacing $R$ by $R'$.
 In practice, we have found that splitting the region $R$ into many convex polygonal pieces $R_j'$, such that $\bigcup_j R_j' = R_0$ (we used $80$ pieces $R_j'$ for $L = 100$) greatly improved the running time of the algorithm. The fact that some tilings could occur in several of the pieces was found to play a negligible role in terms of performance. Additionally, this method allows the algorithm to be parallelized across multiple cores or nodes, since the computation for every piece $R_i'$ is independent of all others.

Having enumerated the candidate set, we now  describe how to use   the enumeration Algorithm \ref{algorithm-enumerate-tilings}  to check whether the Hamiltonian on each patch satisfies the condition in   Proposition~\ref{proposition-recall}. Specifically, we have to compute the norms
\begin{equation}
\label{eq:Dx}
D(x) =  \norm{\mathbf{1}_{\overline{B_l}(x)} (H_{A(x)} - \lambda)^{-1}  \mathbf{1}_{A(x)} H \mathbf{1}_{B_L(x) \backslash {A(x)}}}
\end{equation}
and check that
\begin{align}
  D(x) < N^{-d/2}
\label{eq-resolvent-property}
\end{align}
for all $x \in \Gamma$. Let us first explain in detail how we choose $A(x)$. Since the norm of $(H_{A(x)} - \lambda)^{-1}$ is expected to fall off exponentially in the distance between $\mathbf{1}_{\overline{B_l}(x)}$ and $\mathbf{1}_{A(x)} H \mathbf{1}_{B_L(x) \backslash A(x)}$, $D(x)$ is expected to be minimal for $A(x) = B_{L - m}(x)$.
It turns out that this is true in most cases. However,
if $\lambda$ is very close to the spectrum of $H_{B_{L-m}(x)}$, which happened in our simulations only for a
 few local patches, then
 the norm of $(H_{B_{L-m}(x)} - \lambda)^{-1}$ and thus $D(x)$ may become too large.
 In those cases we  computed  \eqref{eq:Dx} again with a different choice of $A(x)$.
 We found that already removing one site chosen at random from the edge set $\Gamma_\mathrm{AB} \cap (B_{L - m}(x) \backslash B_{L - m -1}(x))$ from $B_{L-m}(x)$ was usually enough to perturb the spectrum of $H_{B_{L-m}(x)}$ sufficiently to remove outliers in the value of $D(x)$.

To actually compute the norm in \eqref{eq:Dx}, the most obvious method would be to invert $(H_{A(x)} - \lambda)$, compute the matrix products as written, and then compute the norm. However, computing matrix inverses is expensive and prone to numerical error, and it would not be feasible at all for matrices of the size we are considering in this paper. The method of computing $D(x)$ we   present now, by contrast, even works with sparse matrices, which makes it efficient enough that it can be employed for values of $L$ where an investigation of the spectrum of $H_{B_L(x)}$ by direct diagonalization would not be possible (although many algorithms are of course available for approximating the spectrum of sparse Hermitian operators \cite{saad2011}).

Instead of computing the matrix product as written in the right-hand side of \eqref{eq:Dx}, we compute the individual entries of the matrix. Let $M(x) := \mathbf{1}_{\overline{B_l}(x)} (H_{A(x)} - \lambda)^{-1}  \mathbf{1}_{A(x)} H \mathbf{1}_{B_L(x) \backslash A(x)}$.
Then, for $\tilde{x} \in \Gamma \cap\overline{B_l}(x)$ and $\tilde{y} \in \Gamma \cap (B_L(x) \backslash A(x))$, we compute the matrix entries
\begin{equation*}
M(x)_{\tilde{x}\tilde{y}} = \mathbf{1}_{\{\tilde{x}\}} (H_{A(x)} - \lambda)^{-1}  \mathbf{1}_{A(x)} H \mathbf{1}_{\{\tilde{y}\}} \, .
\end{equation*}
Depending on whether there are more $\tilde{x}$'s or $\tilde{y}$'s to consider,
we compute either the vector
\begin{align*}
\mathbf{1}_{\{\tilde{x}\}} (H_{A(x)} - \lambda)^{-1} \mathbf{1}_{A(x)} H
\quad\mbox{or}\quad
(H_{A(x)} - \lambda)^{-1} \mathbf{1}_{A(x)} H \mathbf{1}_{\tilde{y}}.
\end{align*}
This is easier because for a given vector $\mathbf{z}$ we can compute $(H_{A(x)} - \lambda)\mathbf{z}$ using the sparse matrix $H_{A(x)}$. The efficiency of the algorithm can be greatly increased by preparing a decomposition of $(H_{A(x)} - \lambda)$ that can be used to solve $(H_{A(x)} - \lambda)^{-1} \mathbf{z}$ for different vectors $\mathbf{z}$; we use a sparse LU decomposition \cite{superlu} for our computations. We then compute the norm $D(x)$ from these matrix entries. The following algorithm describes the computation we use to check \eqref{eq-resolvent-property}.

\begin{algorithm}[H]
\caption{Check equation \eqref{eq-resolvent-property}
\label{algorithm-check-resolvent}}
\begin{algorithmic}[1]
\State Compute the sparse LU factorization of $H_{A(x)} - \lambda$
\State Initialize an empty matrix $m \times n$ matrix $T$, where $m = |B_l(x)|_\Gamma$ and $n = |B_L(x) \backslash {A(x)}|_\Gamma$
\If{$n < m$} \label{if-resolvent-algorithm}
\ForAll{$\tilde{y} \in (B_L(x) \backslash A) \cap \Gamma$}
\State Set $b = (H\delta_{\tilde{y}})|_{A(x)}$ as a vector in $l^2({A(x)} \cap \Gamma)$
\State Solve $(H_{A(x)}- \lambda) y = b$ using the LU factorization
\State Put $y$ in the column of $T$ corresponding to $\tilde{y}$
\EndFor
\Else
\State Initialize an $m \times P$ matrix $\tilde T$, where $P = |{A(x)}|_\Gamma$.
\ForAll{$\tilde{x} \in B_l(x) \cap \Gamma$}
\State Set $b = \delta_{\tilde{x}}$ as a vector in $l^2(A \cap \Gamma)$
\State Solve $(H_{A(x)} - \lambda) y = b$ using the LU factorization
\State Put $y$ in the row of $\tilde T$ corresponding to $\tilde{x}$
\EndFor
\State Set $T = \tilde T H \mathbf{1}_{B_L(x) \backslash {A(x)}}$, an $m \times n$ matrix.
\State Set $D(x) = \norm{T}_{\mathrm{op}}$ and check $D(x) < N^{-d/2}$.
\EndIf
\end{algorithmic}
\end{algorithm}

The two branches of the ``if'' statement in line~\ref{if-resolvent-algorithm} of Algorithm \ref{algorithm-check-resolvent} always compute the same number, and the condition is only an optimization that transposes the matrix product in \eqref{eq-resolvent-property} in order to reduce the number of linear systems that have to be solved.

Finally, Algorithm \ref{algorithm-final} summarizes our general strategy, using the enumeration of the local patches from Algorithm~\ref{algorithm-enumerate-tilings} and checking the condition of Propositon~\ref{proposition-recall} on each of them using Algorithm~\ref{algorithm-check-resolvent}. The loop in line~\ref{outer-loop} can be performed on different computer nodes if necessary.

\begin{algorithm}[H]
\caption{Prove a gap at energy $\lambda$ \label{algorithm-final}}
\begin{algorithmic}[1]
\State Split the region $R$ into $n$ smaller polygons $R_1, \dots, R_j, \dots, R_n$ (for efficiency)
\ForAll{polygons $R_j$}\label{outer-loop}
\State Using Algorithm~\ref{algorithm-enumerate-tilings}, decompose $R_j$ into a number of polygons $P_{(b_i)}$ corresponding to local patches in $\mathcal{C}_L$.
\State Initialize $r_{\min, j} = \infty$, the minimum gap size
\ForAll{polygons $P_{(b_i)}$ corresponding to a bit string $(b_i)$}
\State Construct the Hamiltonian $H$ on $B_L(x)$, for an $x$ with $\kappa(p^{-1}(x)) \in P_{(b_i)}$
\State Generate a finite list $A_k$ of values of $A(x)$ to try. We always set $A_1= B_{L-m}(x)$, while $A_2, A_3, \dots, A_t$ are generated by removig random points from the edge of $A_1$.
\ForAll{$A = A_1, A_2, \dots, A_t$}
\State Set $A = B_{L-m}(x)$ \label{line-set-a}
\State Use Algorithm~\ref{algorithm-check-resolvent} to check equation \eqref{eq-resolvent-property}
\If{equation \eqref{eq-resolvent-property} is fulfilled}
\State Set $r_{\min,j} \leftarrow \min(r_{\min,j}, \epsilon],$ where $\epsilon$ is the maximum allowed by equation~\eqref{eq-conditions-epsilon-recall}.
\State \textbf{continue} with next polygon $P_{(b_i)}$
\EndIf
\EndFor
\If {Equation~\eqref{eq-resolvent-property} was not fulfilled for any $A$}
\State \textbf{end computation}, no gap could be proven
\EndIf
\EndFor
\State Let $r_{\min} = \min(r_{\min, 1}, \dots, r_{\min, n})$
\State \textbf{end computation}, the infinite Hamiltonian has a gap of size $r_{\min}$ at energy $\lambda$.
\EndFor
\end{algorithmic}
\end{algorithm}

\section{3. Fibonacci Quasicrystals}

\noindent
The Fibonacci quasicrystal is a simple one-dimensional quasicrystal that was studied even prior to the discovery of physical quasicrystals \cite{kohomoto1983, ostlund1983}.  In recent years, significant attention has been devoted to the mathematical rigorous study of the spectrum of the Hamiltonian associated to the Fibonacci quasicrystals \cite{jaga2021, even-dar_mandel_electronic_2008,damanik2014}. In particular, it has been proved that the spectrum of the  Fibonacci Hamiltonian is a Cantor set \cite{casdagli1986, suto1987}. In this section, we will describe how our method can be applied to systems in one dimension using the explicit example of the Fibonacci quasicrystal, which has the advantage that many of the constructions are easier to visualize in such case. In particular, we compute upper and lower bound for the distance to the spectrum for the Fibonacci Hamiltonian which clearly show the fractal structure of its spectrum, see Figure \ref{fig-fib-both-bounds}.

\begin{figure}[ht]
	\centering
	\includegraphics{sfig_gaps_wrapper.pdf}
	\caption{Upper and lower bounds for the distance to the spectrum for the Fibonacci Hamiltonian, defined by \eqref{FibonacciHamiltonian}, for $L = 500, N = 6, \alpha = 1$. The transparent blue intervals display the minimal and maximal sizes (defined as in Figure~4) of the gaps centered around the local maxima of the lower bound. For the Fibonacci quasicrystal, multiple gaps can be seen and proven to exist. As $L \to \infty$, the number of gaps will grow as the spectrum of the infinite Hamiltonian is a Cantor set.\label{fig-fib-both-bounds}}
\end{figure}

\subsection{Cut-and-project construction of the Fibonacci quasicrystal}

\noindent
As in the Ammann-Beenker case, we will define two projections, in this case from $\mathbb{R}^2 \rightarrow \mathbb{R}$, corresponding respectively to the real space and to the additional dimension:

\begin{align}
p = \begin{pmatrix}1 & \varphi\end{pmatrix} && \kappa = \begin{pmatrix}-\varphi & 1\end{pmatrix} ,
\end{align}
where $\varphi:=\frac{1+\sqrt5}{2}$ is the golden ration. Clearly the kernels of $p$ and $\kappa$ are again orthogonal.

The acceptance region in the case of the Fibonacci quasicrystal consists simply of the interval
\begin{align}
R = [0,1)\,.
\end{align}
This is the projection of the vertical interval $\{0\}\times[0,1)$ via $\kappa$.

We can then define the Fibonacci lattice as
\begin{align*}
\Gamma_{\text{Fib}} = \big\{\; p(z) \;\big|\; z \in \mathbb{Z}^2,\; \kappa(z) \in R \;\big\} \, .
\end{align*}
The condition $\kappa(z) \in R$ corresponds to the ``cutting'' step, the expression $p(z)$ is the ``projection'' step.

\begin{figure}[ht]
\includegraphics{sfig1_wrapper.pdf}
\caption{Cut-and-project construction of the Fibonacci quasicrystal. \label{fig-three-offsets}}
\end{figure}

Figure~\ref{fig-three-offsets} contains a pictorial representation of this cut-and-project construction. The yellow shaded area in Figure~\ref{fig-three-offsets} shows points in $\mathbb{R}^2$ for which $\kappa(z) \in R$.

This definition of the Fibonacci quasicrystal is equivalent to the more common one that uses the substitution rules
\begin{align*}
\mathrm{S} \rightarrow \mathrm{L} && \mathrm{L} \rightarrow \mathrm{LS}\,.
\end{align*}
Starting from the string ``L'', this substitution rule gives a sequence of strings ``LS'', ``LSL'', ``LSLLS'', ``LSLLSLSL''$\dots$ in which the $n$-th string is a prefix of the $(n+1)$-th and a suffix of the $(n+2)$-th. This allows one to define a Fibonacci string that is infinite in both directions.

The equivalence is based on the fact that neighbouring points in the Fibonacci lattice form either ``long'' or ``short'' distances, corresponding to the letters ``L'' and ``S'' in the substitution definition. As one can see in Figure~\ref{fig-three-offsets}, successive points $z$ with $\kappa(z) \in R$ always differ by either the vector $(1, 1)$ or the vector $(1,2)$. Under the projection $p$, these vectors map to offsets
\begin{align*}
(1 + \varphi)\quad\text{and}\quad   (1 + 2 \varphi) \, .
\end{align*}
The intervals of length $(1 + 2\varphi)$ are the ``long'' intervals L and the intervals of length $(1 + \varphi)$ are the ``short'' intervals S. The quotient between the lengths is again $\varphi$.

\subsection{Enumeration of local patches}

\noindent
Let us show how to enumerate all local patches $\mathcal{C}_L$ of the Fibonacci quasicrystal. As in the case of the Ammann-Beenker tiling, for any $x \in \Gamma_{\text{Fib}}$, the \textit{local patch around $x$} is the set
\begin{align*}
\mathcal{C}_L(x) = \left\{ \,\tilde x - x \,\middle|\, \tilde x \in \Gamma_{\text{Fib}} \cap B_L(x)\,\right\}.
\end{align*}
Since every point in $x \in \Gamma_{\text{Fib}}$ has exactly one preimage $z \in \mathbb{Z}^2$ with $p(z) = x$ and $\kappa(z) \in R$, by exploiting the linearity of $p$ and $\kappa$, the set of local patch can be rewritten in the usual form as
\begin{align*}
\mathcal{C}_L(x) = \left\{ \,p(v) \;\middle|\; v \in \mathbb{Z}^2, \kappa(v) \in R - \kappa(z), |p(v)| < L \,\right\}.
\end{align*}
This description of the local patch only depends on $\kappa(z)$, where $z$ is the integer preimage of $x$ under $p$ (which is unique by irrationality considerations). Our enumeration algorithm will decompose the interval $[0,1)$ in which $\kappa(z)$ lies into subintervals corresponding to different local patches.

It turns out that for any given $L$, there are only finitely many points $v \in \mathbb{Z}^2$ which can fulfill the two conditions $\kappa(v) \in R - \kappa(z)$ and $|p(v)| < L$, across all $z \in \mathbb{Z}^2$ with $\kappa(z) \in R$. Indeed, for a point $v \in \mathbb{Z}^2$ to be able to satisfy the condition $\kappa(v) \in R - \kappa(z)$ for any $z$ with $\kappa(z) \in R$, we must have
\begin{align*}
\kappa(v) \in R + (- R)
\end{align*}
where $R +(- R)$ denotes the Minkowski sum. This equals the interval $R +(- R) = (-1, 1)$.  As for the Ammann-Beenker tiling, we define the set of ``candidate points'' as
\begin{equation*}
V_L:=\{v \in \mathbb{Z}^2 \, | \, \kappa(v) \in (-1, 1) \, ; \, p(v) \in (- L, L)\, \}.
\end{equation*}

\begin{figure}[ht]
\centering \includegraphics{sfig2_wrapper.pdf}
\caption{The set of candidate points is defined by the two constraints $|p(v)| < L$ and $\kappa(v) \in (-1,1)$. In this picture, the constraint on $p(v)$ is indicated by the blue shaded area, and the restriction on $\kappa(v)$ by the orange shaded area. Because the projections have orthogonal kernel, it is clear that the intersection of both sets is compact and contains only finitely many integer points. \label{fig-candidates}}
\end{figure}

Because the two projections $p$ and $\kappa$ have orthogonal kernels, the two sets where the two conditions $\kappa(v) \in (-1, 1)$ and $p(v) \in (-L, L)$ respectively are fulfilled are two orthogonal strips, whose intersection is a rectangle, as shown in Figure~\ref{fig-candidates}. Any finite rectangle contains a finite number of points, therefore the set $V_L$ ist finite.

It is also easy to enumerate these candidate points algorithmically. To do so we describe here an alogorithm similar to Algorithm \ref{algorithm-enumerate-candidates} that we used for the Ammann-Beenker tiling.
In this case we can choose any vector in $\mathbb{R}^2$ to be a right inverse for $p$ (up to a constant factor) and we just need to impose that such vector is in the kernel of $\kappa$. Consider the vector $t:=(1 \;\; \varphi )$, we have $ pt=(1+\varphi^2)$ and we also have that $\kappa t =0$.  Therefore, every point $v \in V_L$ can be written as
\begin{equation*}
v=t y + q e_2
\end{equation*}
for some $y,q \in \mathbb{R}$. Moreover, $v \in \mathbb{Z}^2$ implies that $y \in \mathbb{Z}$, while $\kappa v \in (-1,1)$ implies that $q \in (-1,1)$. Thus, any candidate point is of the form $v=(v_1,v_2)$,  with $v_1 \in \mathbb{Z}$ and $v_2 \in \varphi v_1 + (-1,1)$. If we then consider also the condition $pv \in (-L,L)$, we get that $v_1 \in \frac{(-L,L) + (-1,1)}{1+\varphi^2}$. In particular, if we want to count all the candidate points $v$ such that $pv \in I$, with $I$ a given interval, we have to enumerate all the integers $v_1,v_2$ such that
\begin{equation*}
v_1 \in \frac{I+ (-1,1)}{1+\varphi^2} \qquad v_2 \in \varphi v_1 + (-1,1) \,
\end{equation*}
which results in two loops (over $v_1$ and $v_2$) similar two the loops in Algorithm~\ref{algorithm-enumerate-candidates}.

As we have explained in Section 2 \ref{sm-sec:DescriptionAlg}, once a set of candidate points has been computed, we can categorize the values $\kappa(z)$ by which local patch we get. The local patch is completely determined by which points $\tilde{z}:=z + v \in \mathbb{Z}^2$ become part of the tiling. All points in the candidate set $V_L$ fulfill
\begin{align*}
|p(v)| < L.
\end{align*}
Therefore $|p(z + v) - p(z)| < L$. Thus, the condition $|p(\tilde z) - p(z)| < L$ is always fulfilled for all points in our candidate set.

Instead, whether the other condition is fulfilled, namely $\kappa(\tilde z) \in R$, depends on the base point $z$, or more precisely on $\kappa(z)$. We have $\kappa(\tilde z) = \kappa(z + v) = \kappa(z)  + \kappa(v) \in R$, which implies
\begin{gather}
\kappa(z) \in R-\kappa(v) .\label{eq-whichisanint}
\end{gather}
For a given candidate point $v \in V_L$, the right hand side of \eqref{eq-whichisanint} is just an interval, thus for every candidate point $v \in V_L$, we get one such interval. The local patch of scale $L$ around $p(z)$ only depends on which of these intervals the point $\kappa(z)$ lies in.

Let us enumerate the candidate points as $V_L = \{v_i \,|\, i = 1, \dots, |V_L| \}$. This corresponds to an enumeration of the corresponding intervals $I_i =
R - \kappa(v_i)$. Now, a local patch can be described by a bit string $b \in \{ 0, 1 \}^{|V_L|}$, by setting $b_i = 1$ if $\kappa(z) \in I_i$ and $b_i = 0$ otherwise.

Of course, not every bit string will occur in this way. For every bit string $b$, the portion of the acceptance region for which this patch occurs is given by
\begin{gather*}
P_{(b_i)} = \bigcap_{i = 1}^{|V_L|} P^{b_i}(v_i)
\intertext{where}
P^{1}(v_i) = R \cap I_i, \quad
P^{0}(v_i) = R \setminus P^1(v_i).
\end{gather*}
For every $k = \kappa(z) \in P_{(b_i)}$, the bit string associated to the local patch around $x=p(z)$ is exactly $P_{(b_i)}$. If a bit string $b=(b_i)$ does not actually occur, the set $P_{(b_i)}$ will be empty.

To enumerate all actually possible bit strings $(b_i)$, we can start with only the first candidate point $\{v_1\}$, then with the first two $\{v_1, v_2\}$ and so on, proceeding iteratively. As in Section 2\ref{sm-sec:DescriptionAlg}, we extend the notation $P_{(b_i)}$ to bit strings of length $|b|\leq |V_L|$ simply by setting
\begin{align*}
P_{(b_i)} = \bigcap_{i = 1}^{|b|} P^{b_i}(v_i) \, .
\end{align*}
For every candidate we add, we can go from a bit string of lenght $n$ to a bit string of lenght $n+1$ using the following recursion relation:
\begin{equation*}
\begin{aligned}
P_{(b_i)\oplus 1} &= P_{(b_i)} \cap P^1(v_{i+1})\\
P_{(b_i)\oplus 0} &= P_{(b_i)} \cap P^0(v_{i+1})\,.
\end{aligned}
\end{equation*}

This suggests a simple algorithm for the enumeration of all nonempty $P_{(b_i)}$, with string lenght $|b|=|V_L|$, and their associated bit strings.

The algorithm consists of $|V_L|+1$ steps, $k = 0, \dots, |V_L|$. At every step, we maintain a list $\ell_k$ of all bit strings $b$ of length $k$ for which $P_{(b_i)}$ is nonempty, together with the associated $P_{(b_i)}$. Initially, the list only consists of the empty bit string $\varnothing$. The associated set is $S_\varnothing = [0,1)$.

In the $k$-th step, we are given the list $\ell_{k-1}$ and want to compute $\ell_k$. So every bit string $b \in \ell_{k-1}$, we compute the two sets $S_{b0}$ and $S_{b1}$ as described above, which simply amounts to an intersection of intervals. In many cases, one of these sets will be empty (but not both). We then add to the list $\ell_k$ all those strings $b'=(b'_i)$ for which $P_{(b'_i)}$ is not empty. In the end, the list $\ell_{|V_L|}$ will contain all the bit strings $b$ for which $P_{(b_i)} \neq \varnothing$. See Figure~\ref{fig-split} for a visualization of this interval splitting procedure.

\begin{figure}[ht]
\centering \includegraphics{sfig3_wrapper.pdf}
\caption{Successive splitting of the intervals as more candidate points are added, for $L = 10$. The intervals correspond to the sets $P_{(b_i)}$, the strings to the partial bit strings $b$. \label{fig-split}. The large blue interval on the left hand side is the interval $[0,1)$ associated to the empty string. Steps where no interval is split are not shown.}
\end{figure}

\subsection{Gap bounds}

\noindent
With the method of enumerating local patches of the tiling described above, we will now describe how our lower bound and the upper bound of \cite{colbrook_2019} can be computed in the Fibonacci case.

We consider the Fibonacci Hamiltonian given by the standard Laplacian on $\ell^2(\mathbb{Z})$ with an electric potential given by the Fibonacci sequence. That is, if we denote by $f(n)$ the sequence with $f(n) = 1$ if the $n$-th character of the infinite Fibonacci sequence is $L$ and $f(n) = 0$ otherwise, we can define an operator on $\ell^2(\mathbb{Z})$ by
\begin{equation}
\label{FibonacciHamiltonian}
\begin{aligned}
H_{xy} = \begin{cases}
\alpha f(x) & x = y\\
-1 & |x - y| = 1\\
0 & \text{otherwise}
\end{cases}\,,
\end{aligned}
\end{equation}
where $\alpha \in \mathbb{R}$ is a constant. For $\alpha = 0$, the operator is the standard Laplacian on $\ell^2(\mathbb{Z})$ with absolutely continuous spectrum. As $\alpha$ moves away from $0$, gaps will open in the spectrum of the Laplacian. An enumeration of the local patch of the Fibonacci string can be made similarly to what is described above.

The upper bound of \cite{colbrook_2019} becomes quite simple in a one-dimensional setting like ours. Indeed, following the procedure described in \cite{colbrook_2019}, in order to compute the upper bound to the distance to the spectrum of a given energy $\lambda$, one is reduced to analyze the rectangular matrix
\begin{align*}
B = \mathbf{1}_{[-(n+1), \dots, n+1]}(H - \lambda \mathbf{1}) \mathbf{1}_{[-n,n]}.
\end{align*}
Let $s$ be the lowest eigenvalue of $B^* B$, then the upper bound on the distance from $\lambda$ to the spectrum that is given in \cite{colbrook_2019} is simply $\sqrt{s}$. Although we could verify that the algorithm utilizing the Cholesky decomposition to compute $s$ given in \cite{colbrook_2019} is faster for large sparse matrices, a standard eigenvalue decomposition was sufficient to compute the values corresponding to the black line in Figure~\ref{fig-fib-both-bounds}.

To compute our lower bound, an enumeration of the possible tilings as outlined above was first made. Using $N = 6$ and $r = 1$, we split each interval $B_L(x)$ into an outer part $B_L(x) \backslash A(x)$, a middle part $A(x) \backslash B_l(x)$ and an inner part $B_l(x)$. It turned out that it was sufficient to always pick the maximal choice $A(x) = B_{L - 1}(x)$. This way, the outer part $B_L(x) \backslash A(x)$ has only two elements. Thus, to compute the norm in \eqref{eq-conditions-epsilon-recall}, we only need to compute
\begin{align*}
\mathbf{1}_{\overline{B_l}(x)} (H_{A(x)} - \lambda)^{-1} \mathbf{1}_{A(x)} H \mathbf{1}_{ B_L(x) \backslash A(x)}\mathbf{\tilde{x}}_i = \mathbf{1}_{\overline{B_l}(x)} (H_{A(x)} - \lambda)^{-1} \mathbf{1}_{A(x)} H \mathbf{\tilde{x}}_i
\end{align*}
for the two vectors $\tilde{\bf x}_1 = \delta_{x+L}$ and $\tilde{\bf x}_2 = \delta_{-L+x}$. The constant $M^{-1}(x)$ may be computed as the lowest eigenvalue of the matrix
\begin{align*}
\mathbf{1}_{\overline{B}_{l(x)}} (H - \lambda) \mathbf{1}_{B_{L-1}(x)}\,.
\end{align*}
The values of this lower bound are plotted as the red curve in Figure~\ref{fig-fib-both-bounds}. It can be seen that our method can resolve even the fine fractal gap structure of the Fibonacci Hamiltonian and gives exact bounds on the extent of gaps of different orders of magnitude.

\section{4. Computational complexity of spectral computations}

\noindent
The extistence of the algorithm presented in this paper seemingly contradicts the statement in \cite{colbrook_2019} that it is impossible to compute spectra with error control in a general setting. In fact, even the spectrum of a diagonal operator on an infinite Hilbert space cannot be computed by an algorithm accessing the matrix elements one-by-one. It is only by requiring the additional structure of finite local complexity that such an algorithm can be found. This can be further elucidated in the framework of the \textit{solvability complexity index} (SCI) \cite{hansen2011, benartzi2015, colbrook2019foundations}: We show that computing the spectrum of an operator of finite local complexity is a problem with solvability complexity index $1$, whereas it is known that the index is $2$ in the general case \cite{colbrook2019foundations}.

\begin{definition} A \textit{computational problem} consists of a tuple $(\Omega, \Lambda, (\mathcal{M}, d), \Xi)$. Here $\Omega$ is the domain, or the set of problems. (In our case, $\Omega$ will be a set of operators on a Hilbert space.) The metric space $(\mathcal M, d)$ is the set of possible solutions, which in our case will be the power set of $\mathbb{R}$ equipped with the Hausdorff metric. The \textit{problem function} $\Xi : \Omega \rightarrow \mathcal M$ describes the exact solution of the problem (for example, the function that maps every operator to its spectrum). Finally, $\Lambda$ is a set of functions, $f_i : \Omega \to \mathbb{R}$, the \textit{evaluation functions}, which the algorithm uses to access information on the given object in $\Omega$.
\end{definition}

\begin{definition}
\label{def-sci}
A computational problem $(\Omega, \Lambda, (\mathcal M, d), \Xi)$ is said to have \textit{solvability complexity index $1$} if and only if there exists a sequence $\Gamma_n$ of functions $\Gamma_n : \Omega \to \mathcal M$ such that:
\begin{enumerate}
\item Every function $\Gamma_n$ can be computed by a finite number of divisions, radicals, and comparisons from a finite number of evaluations $f_i : \Omega \to \mathbb R$. More precisely, we can ask for the computations to be performed by a machine in the Blum-Shub-Smale (BSS) model \cite{blum1998}. Which and how many evaluations are performed can be decided based on the previous evaluations.
\item For every $A \in \Omega$, the computations $\Gamma_n(\Omega)$ converge to $\Xi$ in a controlled way, that is
\begin{align*}
d(\Gamma_n(A), \Xi(A)) < 2^{-n}
\end{align*}
\end{enumerate}
\end{definition}

By passing to a subsequence, this definition of solvability complexity index $1$ is equivalent to giving an algorithm $\Gamma_n$ and an error control function $E_n(A)$ such that $d(\Gamma_n(A), \Xi(A)) < E_n(A)$ for all $A \in \Omega$ and $E_n(A) \to 0$ as $n \to \infty$~\cite{benartzi2015}.

\begin{definition}
A \textit{Hamiltonian on a set} $\Gamma$ is a bounded self-adjoint operator on the Hilbert space $\ell^2(\Gamma, \mathcal H)$, for some finite-dimensional Hilbert space $\mathcal H$.
\end{definition}

\begin{definition}
A set $\Gamma \subseteq \mathbb{R}^d$ is called \textit{uniformly discrete} if there exists an $\epsilon > 0$ such that every two points $x, y \in \Gamma$ have distance $\|x-y\|_{\infty} \geq \epsilon$. The \textit{local patch around} $x$ \textit{at scale} $L$ is the set $\mathcal{C}_L(x):=\left\{ y-x | y \in \Gamma \cap B_L(x) \right\}$.  A set $\Gamma$ is said to have \textit{finite local complexity} if for every $L$, the set
\begin{align*}
\mathcal{C}_L = \big\{ \mathcal{C}_L(x) \,|\, x \in \Gamma \big\}\,,
\end{align*}
called the set of \textit{local patches}, is finite.
\end{definition}

\begin{definition}
\label{definition-finite-local-complexity}
A Hamiltonian $H$ on a set $\Gamma$ is said to have finite local complexity if there is a function $W : \mathbb{N} \to \mathbb{N}$ such that, for every $L \in \mathbb{N}$ and every $x \in \Gamma$, there is an $y \in B_{W(L)}(0) \cap \Gamma$ such that $H_{B_L(x)}$ and $H_{B_L(y)}$ are equivalent in the following sense:
\begin{align*}
H_{B_L(x)} = T_{y-x}^* U^* H_{B_L(y)} U T_{y-x}\,
\end{align*}
where $T_{y-x}$ is the translation operator and $U$ is a diagonal unitary operator. (The operator $U$ corresponds to a change of gauge for the magnetic Laplacian.) We say that $H$ has \textit{maximum hopping length} $m$ if the matrix entries $H_{xy} = 0$ for all $x, y \in \Gamma$ with $d(x,y) > m$
\end{definition}

\noindent
Now, the following no-go theorem is proven in \cite{colbrook_2019}.

\begin{theorem}
For every $N_{\max} \in \mathbb{R}$ and $m \in \mathbb R$, define the set $\Omega$ as the set of all Hamiltonians $H$ on uniformly discrete sets $\Gamma$ with $\norm{H}_{\text{op}} \leq N_{\max}$ and maximum hopping length $m$. Then the computational problem $(\Omega, \Lambda, (\mathcal M, d), \Xi|_{\Omega})$ has solvability complexity index $> 1$.
\end{theorem}

The proof of this no-go theorem is surprisingly simple and generalizes across all computational architectures. In fact, the authors remark that it even generalizes to the class of \textit{diagonal operators}. It is based on the fact that any algorithm can only see a finite part of a given infinite matrix, and by changing the matrix in one place outside this range, we can change the spectrum in a way that cannot be detected by the given algorithm. (The counterexample matrix thus has to be adapted to the algorithm.)

It was suggested in \cite[Remark 9.1]{benartzi2015} that it might be a fruitful subject of research to find some additional structure under which one can reduce the SCI, and that this could lead to new algorithms. We show here, using the algorithm described in the previous sections, that adding the structure of \textit{finite local complexity} lowers the SCI of computing the spectrum of self-adjoint operators with bounded dispersion from two to one.

\newcommand{^{\mathrm{fc}}}{^{\mathrm{fc}}}
\begin{theorem}
\label{theorem-fc-complexity}
Let $\Omega^{\mathrm{fc}}$ be the set of all tuples $(\Gamma, H)$, where $\Gamma$ is a uniformly discrete subset of $\mathbb{R}^d$ and $H$ is a Hamiltonian on the set $\Gamma$ with finite local complexity and finite hopping length. We need to introduce several evaluation functions to allow the algorithm to access the geometry of $\Gamma$. We assume that the set $\Gamma$ comes with an enumeration $(\gamma_k)_{k \in \mathbb{N}}$ of its points such that $\norm{\gamma_k}$ is nondecreasing.
\begin{enumerate}
\item for every $k_1, k_2 \in \mathbb{N}$, let $f^{\mathrm{mat}}_{k_1, k_2}(\Gamma, H) := \langle \gamma_{k_1} | H | \gamma_{k_2} \rangle$ evaluate the matrix entries of $H$.
\item for every $k \in \mathbb{N}$ and $l \in \{1, \dots, d\}$, let $f^{\mathrm{pos}}_{k, l}(\Gamma, H) := (\gamma_k)_l$ be the $l$-th coordinate of the $k$-th point.
\item for any $n \in \mathbb{N}$, let $f^{\mathrm{flc}}_n(\Gamma, H) := R$ be smallest radius $R \in \mathbb{R}$ such that every local patch of size $n$ (that is, every set $(\Gamma \cap B_n(x)) - x$ for $x \in \Gamma$) occurs for some $x \in B_R(0)$. (This function quantifies the finite local complexity)
\item let $f^{\mathrm{hop}}(H, \Gamma)$ be the maximum hopping length of $H$.
\end{enumerate}
77Let the set $\Lambda^{\mathrm{fc}}$ comprise all the functions $f^{\mathrm{mat}}_{k_1, k_2}, f^{\mathrm{pos}}_{k, l}, f^{\mathrm{flc}}_n, f^{\mathrm{hop}}, $ and $f^{\mathrm{norm}}$. Finally, as before, define $(\mathcal{M}^{\mathrm{fc}}, d)$ as the power set of $\mathbb{R}$ equipped with the Hausdorff distance and let $\Xi^{\mathrm{fc}}(\Gamma, H) = \sigma(H)$, where $\sigma(H)$ denotes the spectrum of $H$.

The computational problem $(\Omega^{\mathrm{fc}}, \Lambda^{\mathrm{fc}}, (\mathcal M^{\mathrm{fc}}, d), \Xi^{\mathrm{fc}})$ so defined is $\in \Delta_1^A$. (That is, its SCI is $1$.)
\end{theorem}

A detailed proof of this theorem will be given in a future publication. Basically, one can combine the upper bound on the distance to the spectrum from \cite[Supplementary Information, Theorem 3]{colbrook_2019} with the computable lower bound we provide here. The pertinent theorem on the upper bound from \cite{colbrook_2019} can be forumulated as follows:

\newcommand{^{\mathrm{up}}}{^{\mathrm{up}}}
\newcommand{^{\mathrm{low}}}{^{\mathrm{low}}}
\begin{theorem}
\label{theorem-colbrook-bound-sci}
Let $A \in \Omega_2$ and let $\Gamma_n^{\mathrm{up}}(A), E^{\mathrm{up}}(n)$ be computed using the algorithm in \cite{colbrook_2019} with size parameter $n$. Then $\Gamma_n^{\mathrm{up}}(A) \to \sigma(A)$ and $E^{\mathrm{up}}(n) \to 0$ as $n \to \infty$ and $\Gamma_n^{\mathrm{up}}(n)$ is contained in the $E^{\mathrm{up}}(n)$ neighbourhood of $\sigma(A)$. Moreover, $\Gamma_n$ can be implemented using finitely many arithmetic operations and comparisons on the matrix elements of $A$.
\end{theorem}

Using the algorithm in this article, we can also give the following statement on a computable lower bound for operators of finite local complexity:

\begin{theorem}
\label{theorem-our-bound-sci}
For the computational problem $(\Omega^{\mathrm{fc}}, \Lambda^{\mathrm{fc}}, (\mathcal M^{\mathrm{fc}}, d), \Xi^{\mathrm{fc}})$, there exists an algorithm computing a spectral approximation $\Gamma^{\mathrm{low}}_n(H)$ and an error estimate $E^{\mathrm{low}}_n$ such that $\sigma(H)$ is contained in an $E^{\mathrm{low}}_n$-neighbourhood of $\Gamma^{\mathrm{low}}_n$.
\end{theorem}

From the combination of Theorems \ref{theorem-colbrook-bound-sci} and \ref{theorem-our-bound-sci}, we can clearly get an approximation of the spectrum that converges in Hausdorff distance with error control. (To get this set, simply compute the intersection of $\Gamma_n^{\mathrm{low}}$ with the $E_n^{\mathrm{up}}$ neighbourhood of $\Gamma_n^{\mathrm{up}}$ and set $E_n = \max(E_n^{\mathrm{up}}, E_n^{\mathrm{low}})$.) Hence the SCI of the problem $(\Omega^{\mathrm{fc}}, \Lambda^{\mathrm{fc}}, (\mathcal M^{\mathrm{fc}}, d), \Xi^{\mathrm{fc}})$ is 1.

The proof of Theorem \ref{theorem-our-bound-sci} is based on the following results.

\begin{theorem}
\label{theorem-always-bulk-gapped}
Let $H$ be an operator on a uniformly discrete set $\Gamma \subseteq \mathbb{R}^d$ with finite hopping length $m \in \mathbb{R}$. Then for every $\delta > 0$ there exists an $\epsilon > 0$ and an $L \in \mathbb{N}$, such that $H$ is $\epsilon$-bulk-gapped at scale $L$ for all energies $\lambda$ where the distance to the bulk spectrum $\rm{dist}(\lambda, \sigma(H)) > \delta$.
\end{theorem}

The tuple $(\epsilon, L)$ can be computed explicitly and  depends only on $\norm{H}$, on $m$ and on $\Gamma$. The proof of Theorem \ref{theorem-always-bulk-gapped}, which will be presented in a future publication, relies on a discrete version of the well-known Combes-Thomas estimates \cite{combes_thomas}.


\section{Introduction}
The connection between Graph Theory and Coding Theory has been considered by several researchers. The adjacency matrix of a simple graph
with $n$ vertices, labeled by $v_1, v_2,\ldots , v_n$, is a square matrix  of order $n$, whose $(i, j)$-entry is equal to $1$ if the
vertices $v_i$ and $v_j$ are adjacent and is equal to zero otherwise, i.e., the adjacency matrix of a simple graph is a symmetric binary
matrix with diagonal zero which has make it suitable for constructing binary codes, see \cite{rouayheb2012graph,jungnickel1996codes,mallik2021graph,tonchev2002error} and the references quoted therein. Different types of
codes are obtained from this construction depending of the structure of the involved graph. On the other hand, the
work of M. Tanner \cite{tanner1981recursive} introduced the construction of bipartite graphs from the parity-check matrix of codes. Special focus has been dedicated in the application of Tanner graphs to the encoding and decoding algorithms for low-density parity-check codes (LDPC). Intensive work has been done
to analyze the properties of these graphs and to study the characteristics of the related codes, see \cite{etzion1999codes,kschischang2003codes,banihashemi2001tanner,halford2006codes}.

Let $\mathbb{F}_2$ be the field with $2$ elements. A \emph{binary $[n,k]$-linear code} is a $k$-dimensional subspace of $\mathbb{F}_2^n$.
An element of a binary linear code is called a \emph{codeword}. The \emph{Hamming distance}, $d(\boldsymbol{x},\boldsymbol{y})$, between two codewords $\boldsymbol{x}=(x_1,\ldots,x_n), \boldsymbol{y}=(y_1,\ldots,y_n)\in \mathcal{C}\subseteq\mathbb{F}_2^n$ is the number of entries where they differ,
or equivalently,

$$
d(\boldsymbol{x},\boldsymbol{y})=d(\boldsymbol{x},\boldsymbol{y}) =|\{i:x_i\neq y_i, ~1\leq i\leq n\}|.
$$

For $\boldsymbol{x}\in\mathbb{F}_2^n$, the \emph{Hamming weight} of $\boldsymbol{x}$ is $\wt(\boldsymbol{x})=d(\boldsymbol{x},\boldsymbol{0})$, i.e.,
$\wt(\boldsymbol{x})$ is the number of non-zero coordinates in $\boldsymbol{x}$. The \emph{minimum distance} $d(\mathcal{C})=d$ of a linear
code $\mathcal{C}$ is defined as the minimum weight among all non-zero codewords, thus we called it a binary $[n,k,d]$-linear
code. A \emph{generator matrix} for an $[n, k]$-linear code $\mathcal{C}$ is any $k\times n$ matrix $G$ whose rows form a basis for
$\mathcal{C}$. So the code $\mathcal{C}$ can be seen as

\begin{equation}\label{generator_matrix}
\mathcal{C}=\{ \boldsymbol{x}G: \boldsymbol{x}\in \mathbb{F}_2^k\}.
\end{equation}

Also, as a binary linear code is a subspace of a vector space, it is the kernel of some linear transformation. In particular,
there is an $(n-k)\times n$ matrix $H$, called a \emph{parity-check matrix} for the $[n, k]$-linear code $\mathcal{C}$, such that

\begin{equation}\label{check_matrix}
\mathcal{C} =\{\boldsymbol{x} \in \mathbb{F}_2^n:  H\boldsymbol{x}^T = \boldsymbol{0}\}.
\end{equation}

By using elementary row and column operations, we can bring the generator matrix $G$ into the \emph{standard form} $[I_k|A]$ where
$I_k$ denotes the $k\times k$ identity matrix and $A$ is a $k\times (n-k)$ matrix and, then we can obtain a parity-check
matrix for the code $\mathcal{C}$ in the standard form, see \cite[Thm. 1.2.1]{Huffman2003}.

\begin{teor}
Let $\mathcal{C}$ be an $[n, k]$-linear code over $\mathbb{F}_2$ with generator matrix $G=[I_k|A]$. Then $H=[A^t|I_{n-k}]$
is a parity-check matrix for $\mathcal{C}$.
\end{teor}

In this paper, we present a graphical representation (\emph{Hasse diagram}) of a binary linear code using the ideas given in
Section 11.7 of \cite{Huffman2003}. Our main interest is to identify properties of the graph constructed:
we demonstrate that the graph defined here is a connected bipartite graph which no contain triangles. In the next section, we recall some basic concepts and results from Coding Theory and Graph Theory and we give our definition of the graph obtained from binary linear code. Finally, in Section \ref{sec3} we present our main results.

\section{Preliminaries}

In general, there are two basic process to do with a binary $[n,k]$-linear code $\mathcal{C}$: encoding and decoding.
The first one consists in convert a vector of $\mathbb{F}_2^k$ and get a codeword, this process can be done using
\eqref{generator_matrix}. In the other hand, to decode a received word it is necessary, in a first stage, to identify
whether it is a codeword or not. It can be solve by mean of \eqref{check_matrix}. But when we are sending a word it
could be appear some errors in the received word, so we need to detect them and (if it is possible) correct them. These
two properties can be estimated using the minimum distance of $\mathcal{C}$, see \cite[Thms. 2.5.6 and 2.5.10]{Ling2004}.

For an $[n,k,d]$-binary linear code $\mathcal{C}$, we can devise an algorithm using a table with $2^{n-k}$ rather than
$2^n$ elements where one can find the nearest codeword by looking up one of these $2^{n-k}$ entries. This general decoding
algorithm for linear codes is called \emph{syndrome decoding}. Because our code $\mathcal{C}$ is an elementary abelian subgroup
of the additive group of $\mathbb{F}_2^n$, its distinct cosets $\boldsymbol{x} + \mathcal{C}$ partition $\mathbb{F}_2^n$ into
$2^{n-k}$ sets of size $2^k$. Two vectors $\boldsymbol{x}$ and $\boldsymbol{y}$ belong to the same coset if and only if
$\boldsymbol{y}-\boldsymbol{x} \in \mathcal{C}$. We denote by
$\mathfrak{cl}(\mathcal{C})=\{\boldsymbol{x} + \mathcal{C}: \boldsymbol{x}\in \mathbb{F}_2^n\}$  the set of cosets of the code
$\mathcal{C}$.

The table with the $2^{n-k}$ elements can be constructed as follows: Let $H$ be a parity-check matrix for $\mathcal{C}$,
that is a generator matrix for the set $\mathcal{C}^{\perp}=\{\boldsymbol{x}\in \mathbb{F}_2^n: \boldsymbol{x}\cdot \boldsymbol{c}=0, \text{ for all $\boldsymbol{c}\in \mathcal{C}$}\}$. The \emph{syndrome} of a word $\boldsymbol{y} \in \mathbb{F}_2^n$ with respect to $H$ is the vector
$S(\boldsymbol{y}) = H\boldsymbol{y}^T \in \mathbb{F}_2^{n-k}$. Since the syndrome of a codeword is $\boldsymbol{0}$, if
$\boldsymbol{y_1}, \boldsymbol{y_2}\in \mathbb{F}_2^n$ are in the same coset of $\mathcal{C}$, then $\boldsymbol{y_1} - \boldsymbol{y_2}=\boldsymbol{c}\in \mathcal{C}$.
Therefore $S(\boldsymbol{y_1})=H(\boldsymbol{y_2}+\boldsymbol{c})^T=H\boldsymbol{y_2}^T+H\boldsymbol{c}^T=S(\boldsymbol{y_2})$. Conversely,
if $S(\boldsymbol{y_1})=S(\boldsymbol{y_2})$, then $H(\boldsymbol{y_2}-\boldsymbol{y_1})^T=\boldsymbol{0}$ and so $\boldsymbol{y_2}-\boldsymbol{y_1}\in \mathcal{C}$,
i.e., we have a way to test whether the word belongs to the code. Moreover, there is a one-to-one correspondence between cosets of
$\mathcal{C}$ and syndromes. In fact, let $\mathcal{C}$ be a binary linear code and  assume the codeword $\boldsymbol{v}$
is transmitted and the word $\boldsymbol{w}$ is received, resulting in the \emph{error pattern} (or error string)
$\boldsymbol{e} = \boldsymbol{w} - \boldsymbol{v} \in \boldsymbol{w} + \mathcal{C}$. Then $\boldsymbol{w} - \boldsymbol{e} = \boldsymbol{v} \in \mathcal{C}$,
so, the error pattern $\boldsymbol{e}$ and the received word $\boldsymbol{w}$ are in the same coset.

Since error patterns of small weight are the most likely to occur, nearest neighbor decoding works for a linear code $\mathcal{C}$ in the following manner. Upon receiving the word $\boldsymbol{w}$, we choose a word $\boldsymbol{e}$ of least weight in the coset $\boldsymbol{w} + \mathcal{C}$ and conclude that $\boldsymbol{v} = \boldsymbol{w}-\boldsymbol{e}$ was the codeword transmitted. To do this it is necessary to construct the (Slepian) \emph{standard array}, for which we need to find the weight of a coset which is the smallest weight of a vector in the coset. Any vector of this smallest weight in the coset is called a \emph{coset leader}. The zero vector is the unique coset leader of the code $\mathcal{C}$. More generally, every coset of weight at most $t = \lfloor(d-1)/2\rfloor$ has a unique coset leader.

There is a natural partial ordering $\preceq$ on the vectors in $\mathbb{F}_2^n$, which is defined as follows: for
$\boldsymbol{x},\boldsymbol{y}\in \mathbb{F}_2^n$, $\boldsymbol{x}\preceq \boldsymbol{y}$ provided that $\supp(\boldsymbol{x}) \subseteq \supp(\boldsymbol{y})$,
where $\supp(\boldsymbol{c})$ for $\boldsymbol{c}\in \mathbb{F}_2^n$ denote the {\it support} of the vector $\boldsymbol{c}$, i.e.,
the set of coordinates where $\boldsymbol{c}$ is non-zero. If $\boldsymbol{x}\preceq \boldsymbol{y}$, we will also say that
$\boldsymbol{y}$ covers $\boldsymbol{x}$. We now use this partial order on $\mathbb{F}_2^n$ to define a partial order, also
denoted $\preceq$, on the set of cosets of a binary linear code $\mathcal{C}$ of length $n$. If $\mathcal{C}_1$ and
$\mathcal{C}_2$ are two cosets of $\mathcal{C}$, then $\mathcal{C}_1 \preceq \mathcal{C}_2$ provided there are coset
leaders $\boldsymbol{x}_1$ of $\mathcal{C}_1$ and $\boldsymbol{x}_2$ of $\mathcal{C}_2$ such that $\boldsymbol{x}_1 \preceq \boldsymbol{x}_2$.
As usual, $\mathcal{C}_1 \prec\mathcal{C}_2$ means that $\mathcal{C}_1 \preceq \mathcal{C}_2$ but $\mathcal{C}_1\neq \mathcal{C}_2$.
Under this partial ordering the set of cosets of $\mathcal{C}$ has a unique minimal element, the code $\mathcal{C}$
itself. Since, the set $\mathfrak{cl}(\mathcal{C})$ is a  partially ordered set, \emph{poset} for short, with respect to
$\preceq$, we can use a Hasse diagram to represent the partial order on $\mathfrak{cl}(\mathcal{C})$.

Let  $\boldsymbol{x}\in \mathbb{F}_2^n$ and $r$ be a non-negative real number. The \emph{sphere of radius  $r$ and center in $\boldsymbol{x}$}
is the set $S(\boldsymbol{x},r)=\{\boldsymbol{y}\in \mathbb{F}_2^n: d(\boldsymbol{x},\boldsymbol{y})\leq r\}.$ The \emph{covering radius} of a code $\mathcal{C}$, denoted by $\rho = \rho(\mathcal{C})$, is the smallest integer
$s$ such that $\mathbb{F}_2^n$ is the union of the spheres of radius $s$ centered at the codewords of $\mathcal{C}$.
Actually, it can be proved that $$\rho(\mathcal{C})=\max_{\boldsymbol{x}\in \mathbb{F}_2^n}\min_{\boldsymbol{c}\in\mathcal{C}}d(\boldsymbol{x},\boldsymbol{c}).$$

Also, by  \cite[Thm. 11.1.2]{Huffman2003}, $\rho(\mathcal{C})$ is the largest value of the minimum weight of all cosets of $\mathcal{C}$, i.e. $\rho(\mathcal{C})=\max\{\wt(\boldsymbol{x}+\mathcal{C}): \boldsymbol{x}\in \mathbb{F}_2^n\}.$

A \emph{graph} $\Gamma$ with $n$ vertices is a pair $(V,E)$ where $V=\{v_1, v_2,\ldots ,v_n\}$ is the set of vertices and
$E \subseteq V\times V$ is the set of edges. Given two vertices $v_i$ and $v_j$, if $v_iv_j\in E$, then $v_i$ and
$v_j$ are said to be \emph{adjacent} or that $v_i$ and $v_j$ are \emph{neighbors}. In this case, $v_i$ and $v_j$ are said to be
the end vertices of the edge $v_iv_j$. If $v_iv_j\notin E$, then $v_i$ and $v_j$ are \emph{non-adjacent}. Furthermore, if
an edge $e$ has a vertex $v_i$ as an end vertex, we say that $v_i$ is incident with $e$. We denote the \emph{set of neighbors} or \emph{neighborhood} of the vertex $v_i$ by $\mathcal{N}(v_i)$,
i.e., $\mathcal{N}(v_i)=\{x\in V: v_ix\in E\}$, and the cardinality  of this set is the \emph{degree} of the vertex $v_i$,
which is denoted by $\mathop{\rm deg}\nolimits(v_i)$, in other words, the $\mathop{\rm deg}\nolimits(v_i)$ is the number of edges incident with $v_i$.

Given a binary linear code $\mathcal{C}$, we denote by  $\Gamma(\mathcal{C})=(V_{\mathcal{C}},E_{\mathcal{C}})$ the  graph
constructed on this wise: the set of vertices  $V_{\mathcal{C}}=\mathfrak{cl}(\mathcal{C})$  and $\mathcal{C}_1\mathcal{C}_2\in E_{\mathcal{C}}$
if $\mathcal{C}_1\prec \mathcal{C}_2$ and $\wt(\mathcal{C}_1)=\wt(\mathcal{C}_2)-1$. If $\mathcal{C}_1\mathcal{C}_2$ is an edge of $\Gamma(\mathcal{C})$, $\mathcal{C}_1$ is called a \emph{child} of $\mathcal{C}_2$, and $\mathcal{C}_2$ is a \emph{parent} of $\mathcal{C}_1$. Actually, the graph $\Gamma(\mathcal{C})$ is known as the \emph{Hasse diagram} of the poset $(\mathfrak{cl}(\mathcal{C}),\preceq)$.  It is clear that $|V_{\mathcal{C}}|=2^{n-k}$.
The main interest of this manuscript is to study the properties of the graph $\Gamma(\mathcal{C})$.

\begin{ejem}
\label{ejemplorelaciones}
	Let $\mathcal{C}$ be a $[5, 2, 3]$ binary linear code with generator matrix
	\begin{center}
		$G=\begin{pmatrix}
			1 & 1 & 1 & 1 & 0 \\
			0 & 0 & 1 & 1 & 1
		\end{pmatrix}.$
	\end{center}
	The cosets of $\mathcal{C}$ are
	\begin{align*}
		\mathcal{C}_0 &= 00000 + \mathcal{C},    &   \mathcal{C}_4 &= 00010 + \mathcal{C}, \\
		\mathcal{C}_1 &= 10000 + \mathcal{C},    &   \mathcal{C}_5 &= 00001 + \mathcal{C},\\
		\mathcal{C}_2 &= 01000 + \mathcal{C},    &   \mathcal{C}_6 &= 10100 + \mathcal{C} = 01010 + \mathcal{C},\\
		\mathcal{C}_3 &= 00100 + \mathcal{C},    &   \mathcal{C}_7 &= 10010 + \mathcal{C} = 01100 + \mathcal{C}.
    \end{align*}

The visual representation of the graph (Hasse diagram) $\Gamma({\mathcal{C}})$ is given in the Figure \ref{fig:1}.
	\begin{figure}[H]
		\centering
\hspace*{2.5cm}\includegraphics[scale=0.7]{Ejemplo_grafo1.pdf}
		\caption{Hasse diagram of the code $\mathcal{C}$.}
		\label{fig:1}
	\end{figure}
\end{ejem}

\textbf{Remark.} If $\mathcal{C}_1$ and $\mathcal{C}_2$ are cosets of $\mathcal{C}$ with $\mathcal{C}_1 \prec \mathcal{C}_2$, then $\mathcal{C}_1$ is called a \emph{descendant} of $\mathcal{C}_2$, and
$\mathcal{C}_2$ is an \emph{ancestor} of $\mathcal{C}_1$. Note that the coset $\mathcal{C}_0=\boldsymbol{0}+\mathcal{C}$ is always a descendant of any coset of the linear code $\mathcal{C}$, in other words $\mathcal{C}_0$ is minimal in the poset $(\mathfrak{cl}(\mathcal{C}),\preceq)$. Always, this vertex will be represented with a red dot in the figures. A coset of $\mathcal{C}$ is called an \emph{orphan} whenever it has no parents, that is, it is a maximal element in the poset $(\mathfrak{cl}(\mathcal{C}),\preceq)$. These vertices will be denoted, in the figures, with a blue dot, see Figure \ref{fig:1}.

\begin{ejem}
Let $\mathcal{C}$ be a binary linear code with parameters $[6,2,4]$ with generator matrix
$$G= \begin{pmatrix}
		1 & 1 & 1 & 1 & 0 & 0\\
		0 & 0 & 1 & 1 & 1 & 1\\
	\end{pmatrix}.$$
	It can be verified that $\mathcal{C}$ has 16 cosets, its graph $\Gamma(\mathcal{C})$ is given in the Figure \ref{ejem_grafo4} and the cosets $\mathcal{C}_7= +\mathcal{C}$, $\mathcal{C}_{14}=+ \mathcal{C}$ and $\mathcal{C}_{15}=+ \mathcal{C}$ are its orphans.

	\begin{figure}[H]
		\centering
		\includegraphics[scale=0.2]{Ejemplo_grafo4.pdf}
		\caption{Hasse diagram of the code $\mathcal{C}$.}
		\label{ejem_grafo4}
\end{figure}
\end{ejem}

Let us recall other definitions of graph theory, the reader can find them in \cite{chartrand2013first,B98}. Given a
graph $\Gamma$, a \emph{$u-v$ walk} $W$ in $\Gamma$ is an alternating sequence of vertices and edges
$v_0e_1v_1e_2\ldots v_{n-1}e_nv_n$ beginning with $v_0=u$ and ending at $v_n=v$, in which each
edge is incident with the two vertices immediately preceding and following it, i.e., $e_i=v_{i-1}v_i, ~1\leq i\leq n$.
This walk joins $v_0$ and $v_n$, and may also be denoted $W=v_0v_1v_2\cdots v_n$. The number of edges in
a walk (including multiple occurrences of an edge) is called the \emph{length} of the walk. We say that $W$ is \emph{closed} if $v_0=v_n$ and
is \emph{open} otherwise. It is a \emph{trail} if all the edges are distinct, and a \emph{path} if all the vertices (and thus necessarily
all the edges) do not repeat. Hence every path is a trail, but the converse is not true. If the walk is closed, then it
is a \emph{cycle} provided its $n$ vertices are distinct and $n\geq 3$. Hence a cycle begins and ends at the same vertex but
repeats no edges. A cycle of length $n$ is called a \emph{$n$-cycle}. As is usual, we denote by $C_n$ the graph consisting
of a cycle with $n$ vertices and by $P_n$ a path with $n$ vertices; $C_3$ is often called a \emph{triangle}. A cycle of odd length
is called an \emph{odd cycle}; while, otherwise it is called an \emph{even cycle}. A graph without triangles is
called a \emph{triangle-free graph}.

A graph $\Gamma$ is \emph{connected} if every two vertices of $\Gamma$ are joined by a path. A graph $\Gamma$ that is not connected
is called \emph{disconnected}. A maximal connected subgraph of $\Gamma$ is called a \emph{connected component} or simply a component of
$\Gamma$. Thus, a disconnected graph has at least two components. Let $\Gamma=(V,E)$ be a connected graph and $u, v\in V$ two vertices.
The \emph{distance} between $u$ and $v$ is the smallest length of any $u-v$ path in $G$ and is denoted by $d_\Gamma(u,v)$. The
greatest distance between any two vertices of a connected graph $\Gamma$ is called the \emph{diameter} of $\Gamma$ and is denoted by
$\mathop{ \rm diam}\nolimits(\Gamma)$.

 A \emph{coloring} of a graph $\Gamma$, is an assignment of colors (elements of some set) to the vertices of $\Gamma$, one
color to each vertex, in such a way that adjacent vertices have distinct colors. The smallest number of colors in any coloring
of a graph $\Gamma$ is called the \emph{chromatic number} of $\Gamma$ and is denoted by $\chi(\Gamma)$. A graph $\Gamma$ such
that $\chi(\Gamma)=2$ is a \emph{bipartite graph}, i.e., the set of vertices of $\Gamma$ can be decomposed into two disjoint sets
such that no two graph vertices within the same set are adjacent.

\section{Main results}\label{sec3}
In this section we prove the main results about some properties of the graph $\Gamma(\mathcal{C})$, for a given binary linear code $\mathcal{C}$.

\begin{teor}
Let $\mathcal{C}$ be a binary linear code. Then $\Gamma(\mathcal{C})$ is connected.
\end{teor}

\begin{proof}
	Let $\mathcal{C}_k$ and $\mathcal{C}_l$ be two  cosets of $\mathcal{C}$, such that $\wt(\mathcal{C}_k)\leq \wt(\mathcal{C}_l)$, with $\boldsymbol{x}_k$ and $\boldsymbol{x}_l$ their coset leaders, respectively. Then we study the following cases.
	\begin{enumerate}[(i)]
		\item Suppose that $\mathcal{C}_k$ is a descendant of $\mathcal{C}_l$. Take $A_1=\supp(\boldsymbol{x}_l)\setminus \supp(\boldsymbol{x}_k)$, $t=|A_1|$, $i_1=\min A_1$, $\boldsymbol{x}_{i_1}=\boldsymbol{x}_l+\boldsymbol{e}_{i_1}$ and for $2\leq s\leq t$
		\begin{align*}
			A_s&=A_{s-1}\setminus \{i_{s-1}\},\\
			i_s&=\min{A_s}, \\
			\boldsymbol{x}_{i_s}&=\boldsymbol{x}_{i_{s-1}}+\boldsymbol{e}_{i_s}.
		\end{align*}
		Additionally, consider the cosets $\mathcal{C}_{i_s}=\boldsymbol{x}_{i_s}+\mathcal{C}$, for $1\leq s\leq t$. By \cite[Corollary 11.7.7]{Huffman2003},  $\boldsymbol{x}_{i_1}$ is a coset leader of a child of $\mathcal{C}_l$ and also $\mathcal{C}_{i_s}$ is a child of $\mathcal{C}_{i_{s-1}}$. Thus  $\mathcal{C}_l\mathcal{C}_{i_1} \cdots \mathcal{C}_{i_t}$ is a $\mathcal{C}_l-\mathcal{C}_k$ path, where $\mathcal{C}_{i_t}=\mathcal{C}_k$. Observe that this path is the shortest between these two vertices, but it is not necessarily unique. As a consequence we get that $d_{\Gamma}(\mathcal{C}_{k},\mathcal{C}_{l})=t$.

		\item Assume that $\mathcal{C}_k$ is not a descendant of $\mathcal{C}_l$. Take $$\Omega=\{\boldsymbol{y}+\mathcal{C}: \boldsymbol{y}+\mathcal{C}\prec \mathcal{C}_k \text{ and } \boldsymbol{y}+\mathcal{C}\prec \mathcal{C}_l\},$$ that is $\Omega$ is the set of the common descendants of $\mathcal{C}_k$ and $\mathcal{C}_l$. Since that $\mathcal{C}_0=\boldsymbol{0}+\mathcal{C}$ always is a descendant of any coset of $\mathcal{C}$, then $\Omega\neq \emptyset$.  We take $\mathcal{C}''$ such that $\wt(\mathcal{C}'')=\max\{\wt(\mathcal{C}'): \mathcal{C}'\in \Omega\}$.

Since $\mathcal{C}''$ is a descendant of $\mathcal{C}_k$, by (i) there exists a  $\mathcal{C}_k-\mathcal{C}''$ path in $\Gamma(\mathcal{C})$, we call it $\mathcal{C}_k\cdots\mathcal{C}''$. Similarly, there is a $\mathcal{C}''-\mathcal{C}_l$ path: $\mathcal{C}''\mathcal{C}''_1\cdots \mathcal{C}_l$. Thus by the choose of $\mathcal{C}''$ we find a path from $\mathcal{C}_k$ to $\mathcal{C}_l$ in $\Gamma(\mathcal{C})$:  $\mathcal{C}_k\cdots \mathcal{C}''\mathcal{C}''_1\cdots \mathcal{C}_l$.
	\end{enumerate}
\end{proof}

For a coset $\mathcal{C}_k$ of a binary linear $\mathcal{C}$, with $CL(\mathcal{C}_k)$ we denote a complete set of the coset leaders of $\mathcal{C}_k$. Now, from the previous proof, we get the following result.

\begin{cor}\label{cor2.1}
	Let $\mathcal{C}$ be a binary linear code, $\mathcal{C}_k$ and $\mathcal{C}_l$ cosets of $\mathcal{C}$ such that $\boldsymbol{x}_k\in CL(\mathcal{C}_k)$, $\boldsymbol{x}_l\in CL(\mathcal{C}_l)$. Then
	$$
	d_{\Gamma}(\mathcal{C}_k,\mathcal{C}_l)=
	\begin{cases}
		d(\boldsymbol{x}_k,\boldsymbol{x}_l), & \text{if  $\mathcal{C}_k\prec \mathcal{C}_l$ and $\boldsymbol{x}_k\prec \boldsymbol{x}_l$},\\
		d(\boldsymbol{x}_k,\boldsymbol{x}'')+d(\boldsymbol{x}'',\boldsymbol{x}_l), & \text{if  $\mathcal{C}_k\not\prec \mathcal{C}_l$ and $\boldsymbol{x}''\in CL(\mathcal{C}'')$, $\mathcal{C}''\in \Omega$}\\
			& \text{such that $\wt(\mathcal{C}'')=\max\{\wt(\mathcal{C}'): \mathcal{C}'\in \Omega\}$.}
	\end{cases}
	$$
\end{cor}

Since $\mathcal{C}_0$ is a descendant of any coset of $\mathcal{C}$, if $\mathcal{C}_l$ is an orphan such that $\rho(\mathcal{C})=\wt(\mathcal{C}_l)$, that is $\mathcal{C}_l$ is a coset with the largest Hamming weight, see  \cite[Theorem 11.1.2]{Huffman2003} then we get the following result.

\begin{cor}
	Let  $\mathcal{C}$ be a binary linear code. Then  $\mathop{ \rm diam}\nolimits(\Gamma(\mathcal{C}))\geq \rho(\mathcal{C})$.
\end{cor}

The problem of deciding whether a given graph $\Gamma = (V, E)$ contains a
$k$-cycle is among the most natural and easily stated algorithmic graph
problems. In particular, the \emph{triangle finding problem} is the problem of determining whether a graph is triangle-free or not. It is possible to determine if a given graph with $m$ edges is triangle-free in time $O(m^{1.41})$, see \cite[Thm. 3.5]{AYZ97}. In the following result, we solve this problem to the family of  graphs $\Gamma(\mathcal{C})$.

\begin{teor}\label{triangle_free}
Let $\mathcal{C}$ be a binary linear code. Then $\Gamma(\mathcal{C})$ is a triangle-free graph.
\end{teor}
\begin{proof}
Assume that there is a triangle in $\Gamma(\mathcal{C})$ with vertices $\mathcal{C}_k$, $\mathcal{C}_q$ and $\mathcal{C}_m$. Suppose that $\mathcal{C}_k\mathcal{C}_q \in E_{\mathcal{C}}$. Without lost of generality, suppose that $\mathcal{C}_k$ is a child of $\mathcal{C}_q$. Since $\mathcal{C}_k\mathcal{C}_m\in E_{\mathcal{C}}$, then either $\mathcal{C}_k$ is a child of $\mathcal{C}_m$ or $\mathcal{C}_m$ is a child of $\mathcal{C}_k$.

If  $\mathcal{C}_k$ is a child of $\mathcal{C}_m$, then  $\wt(\mathcal{C}_m)=\wt(\mathcal{C}_q)$ and then $\mathcal{C}_q=\mathcal{C}_m$, which is a contradiction.	On the other hand, if  $\mathcal{C}_m$ is a child of $\mathcal{C}_k$, then $\wt(\mathcal{C}_m)=\wt(\mathcal{C}_k)-1$, furthermore $\wt(\mathcal{C}_k)=\wt(\mathcal{C}_q)-1$, thus $\wt(\mathcal{C}_m)=\wt(\mathcal{C}_q)-2$, thus there no exists an edge between $\mathcal{C}_m$ and $\mathcal{C}_q$, so again we obtain a contradiction. Consequently, $\mathcal{C}_k$ and  $\mathcal{C}_m$ are not adjacent vertices in  $\Gamma(\mathcal{C})$, and therefore it is a triangle-free graph.
\end{proof}

The study of the chromatic number of triangle‐free graphs is a classic topic, which has been studied extensively for several researchers and from many perspectives (for instance algebraic, probabilistic, quantitative Ramsey theory and algorithmic), see \cite{davies2020coloring} for more references about this problem.  In fact, triangle-free graphs can have large chromatic numbers, see \cite[Thm. 10.10]{chartrand2013first}. The next theorem demonstrates that if $\mathcal{C}$ is a proper vectorial subspace, then $\chi(\Gamma(\mathcal{C}))=2$.

\begin{teor}
Let $\mathcal{C}\subset \mathbb{F}_2^n$ be a binary linear code. Then $\Gamma(\mathcal{C})$ is a  bipartite graph.
\end{teor}

\begin{proof}
Since $\mathcal{C}\subset \mathbb{F}_2^n$, it is clear that  $|V_{\mathcal{C}}|\geq 2$. Consider $U$ and $V$, the set of cosets of $\mathcal{C}$ of odd weight and even weight, respectively. It is easy to prove that $V_{\mathcal{C}}=U\cup V$ and $U\cap V=\emptyset$. Now, by the definition of $\Gamma(\mathcal{C})$, the vertices of even weight are adjacent uniquely with vertices of odd weight and vice versa. Therefore, any edge in $E_{\mathcal{C}}$ connect nodes from $V_1$ with elements in $V_2$, but not connect vertices in the same set of vertices and the result follows.
\end{proof}
It is known that a nontrivial graph is bipartite  if and only if  it does not contains no odd cycles, see \cite[Theorem 1.12]{chartrand2013first}. Thus the last result generalize Theorem \ref{triangle_free}.

For $n\in \mathbb{Z}^+$, let $S_n$ denote the \emph{symmetric group on $n$ letters}. For   $\pi\in S_n$, a permutation on the coordinates of $\mathbb{F}_2^n$ is a function on $\mathbb{F}_2^n$ defined by $\mathcal{P}(u_1\ldots u_n)=u_{\pi(1)}\ldots u_{\pi(n)}$.  We remember that two binary linear codes $\mathcal{C}$ and $\mathcal{D}$ are \emph{permutation equivalent} if there is a permutation on the coordinates of $\mathbb{F}_2^n$ which sends $\mathcal{C}$ to $\mathcal{D}$, i.e. there exists a permutation $\pi \in S_n$ such that for any codeword $\boldsymbol{v}\in \mathcal{D}$ there exists a codeword $\boldsymbol{u}\in \mathcal{C}$ such that $\boldsymbol{v}=\mathcal{P}(\boldsymbol{u})$. Since that a permutation on the coordinates of $\mathbb{F}_2^n$ preserves the weight of any word, it can be proved that two permutation equivalent binary linear codes have the same parameters, thus they can be considered as ``the same code''. Deciding whether two codes are permutation equivalent codes is known as the \emph{Code equivalence problem}. This problem has been extensively studied in the last decades, see \cite{PR97,sendrier2000finding} and the references therein.

Recall that two graphs are \emph{isomorphic} if there is a correspondence between their vertex sets that preserves adjacency, i.e.,
$\Gamma=(V,E)$ is isomorphic to $\Gamma'=(V',E')$, denoted by $\Gamma\simeq \Gamma'$, if there is a a bijection $\varphi$ from
$V$ to $V'$ such that $xy\in E$ if and only if $\varphi(x)\varphi(y)\in E'$.  The \emph{graph isomorphism problem} requests to determine whether two given graphs are isomorphic. It is well known that deciding if two graphs are isomorphic is an algorithmic problem that has been widely studied, see \cite{Babai2018}.

\begin{teor}\label{grafos_isomorfos}
Let $\mathcal{C}$ and $\mathcal{D}$ two linear binary codes. If $\mathcal{C}$ and $\mathcal{D}$ are permutation equivalent codes, then $\Gamma(\mathcal{C})\simeq \Gamma(\mathcal{D})$.
\end{teor}
\begin{proof}
By the assumption there exists  $\mathcal{P}$ a permutation on the coordinates of $\mathbb{F}_2^n$ such that $\mathcal{D}=\{\mathcal{P}(\boldsymbol{u}):\boldsymbol{u}\in \mathcal{C}\}$. Take $\boldsymbol{u}$ a coset leader of $\boldsymbol{u}+\mathcal{C}$. Since $\mathcal{P}$ preserves Hamming distance, we have that for any $\boldsymbol{c}\in \mathcal{C}$, $\wt(\boldsymbol{u}+\boldsymbol{c})=\wt(\mathcal{P}(\boldsymbol{u})+\boldsymbol{v})$, where $\boldsymbol{v}=\mathcal{P}(\boldsymbol{c})$. Then $\mathcal{P}(\boldsymbol{u})$ is a leader of its coset in $\mathcal{D}$, that is $\mathcal{P}(\boldsymbol{u})+\mathcal{D}$ is a vertex of $\Gamma(\mathcal{D})$.

Let $\mathcal{C}_1=\boldsymbol{u}+\mathcal{C}$ and $\mathcal{C}_2=\boldsymbol{u}'+\mathcal{C}$ be two cosets of $\mathcal{C}$.  Now, if $\mathcal{C}_1\mathcal{C}_2$ is an edge of $\Gamma(\mathcal{C})$, then $\mathcal{C}_1\prec \mathcal{C}_2$. So  there is a coset leader $\boldsymbol{u}''$ of $\mathcal{C}_2$ such that $\boldsymbol{u}\prec \boldsymbol{u}''$; i.e. $\supp(\boldsymbol{u})\subset \supp(\boldsymbol{u}'')$. This last statement implies that $\supp(\mathcal{P}(\boldsymbol{u}))\subset \supp(\mathcal{P}(\boldsymbol{u}''))$ or equivalently we have that $\mathcal{P}(\boldsymbol{u})+\mathcal{D}$ and $\mathcal{P}(\boldsymbol{u}'')+\mathcal{D}=\mathcal{P}(\boldsymbol{u}')+\mathcal{D}$ are adjacent in $\Gamma(\mathcal{D})$. So the function $\varphi: V_{\mathcal{C}}\rightarrow V_{\mathcal{D}}$ given by $\varphi(\boldsymbol{u}+\mathcal{C})=\mathcal{P}(\boldsymbol{u})+\mathcal{D}$ is a graph isomorphism.
\end{proof}

\begin{ejem}
The reciprocal of Theorem \ref{grafos_isomorfos}  is false. In fact, consider the binary linear codes $\mathcal{C}$  and $\mathcal{D}$ with generator matrices $$M=\begin{pmatrix}
     	1 & 0 & 0 & 0 & 1 & 1\\
     	0 & 1 & 0 & 1 & 0 & 1\\
     	0 & 0 & 1 & 1 & 1 & 0\\
     \end{pmatrix} \text{ and }  N= \begin{pmatrix}
		1 & 0 & 0 & 0 & 0 & 1 & 1 & 0\\
		0 & 1 & 0 & 0 & 0 & 1 & 0 & 1\\
		0 & 0 & 1 & 0 & 0 & 0 & 1 & 0\\
		0 & 0 & 0 & 1 & 0 & 0 & 0 & 1\\
		0 & 0 & 0 & 0 & 1 & 1 & 1 & 1\\
	\end{pmatrix},$$
	respectively.  It can be verified that $d(\mathcal{C})=3$ and $d(\mathcal{D})=2$, fact that implies that these codes are non permutation equivalent codes. However, their Hasse diagrams are isomorphic, see Figure \ref{fig_grafos_isomorfos}.

\begin{figure}[H]
	\centering
	\hspace*{2.5cm}\includegraphics[scale=0.25]{Ejemplo_grafo3.pdf}
	\caption{Hasse diagram of $\mathcal{C}$ and $\mathcal{D}$
.}
	\label{fig_grafos_isomorfos}
\end{figure}
\end{ejem}

Theorem \ref{grafos_isomorfos} and the above example motivate us to ask when two non permutation equivalent binary codes $\mathcal{C}$ and $\mathcal{D}$ satisfy that $\Gamma(\mathcal{C})\simeq \Gamma(\mathcal{D})$. More  precisely:

\begin{ques}
What conditions must satisfy two non permutation equivalent binary codes to guarantee that their Hasse diagrams are structurally identical?
\end{ques}

Now that we know that $\Gamma(\mathcal{C})$ is a connected bipartite graph, we can give a particular answer to the last question. In graph theory, for a positive integer $k$, the \emph{star $S_k$} is the complete bipartite graph $K_{1,k}$. Note that if $\Gamma(\mathcal{C})\simeq S_{2^{n-k}}$, this implies that all the cosets of $\mathcal{C}$ have weight at most 1; that is for all $\boldsymbol{u}\notin \mathcal{C}$ there is $i\in \{1,\ldots,n\}$ such that $\boldsymbol{u}\in \boldsymbol{e}_i+\mathcal{C}$, or equivalently  $\boldsymbol{u}=\boldsymbol{e}_i+\boldsymbol{c}$ for some $\boldsymbol{c}\in \mathcal{C}$, thus $d(\boldsymbol{u},\boldsymbol{c})=1$. It means that  for all $\boldsymbol{u}\in \mathbb{F}_2^n$,  there exists $\boldsymbol{c}\in \mathcal{C}$ such that  $d(\boldsymbol{u},\boldsymbol{c})\leq 1$. Therefore, $\mathbb{F}_2^n=\bigcup_{\boldsymbol{c}\in \mathcal{C}}S(\boldsymbol{c},1)$; namely $\rho(\mathcal{C})=1$. It is easy to see that if $\rho(\mathcal{C})=1$, then $\mathfrak{cl}(\mathcal{C})=\{\boldsymbol{e}_i+\mathcal{C}: \text{ for some } 1\leq i\leq n\}$. Consequently, the Hasse diagram  $\Gamma(\mathcal{C})$ is isomorphic to a graph as that one given in Figure \ref{fig:estrellainfinita}.

\begin{figure}[H]
		\centering
		\hspace*{2.5cm}\includegraphics[scale=0.4]{Estrella_infinita.pdf}
		\caption{$\Gamma(\mathcal{C})$ with $\rho(\mathcal{C})=1$.}
		\label{fig:estrellainfinita}
	\end{figure}

Actually, in the last paragraph we have proven the following result.
\begin{teor} Let $\mathcal{C}$ a binary linear code. Then
$\Gamma(\mathcal{C})\simeq S_{2^{n-k}}$ the star graph if and only if $\rho(\mathcal{C})=1$.

\end{teor}
Binary codes with covering radius 1 has been studied in the literature, see \cite{habsieger1997} and the references quoted therein.