# Edge density of bulk states due to relativity

Matthew D. Horner<EMAIL_ADDRESS>Jiannis K. Pachos  School of Physics and
Astronomy, University of Leeds, Leeds, LS2 9JT, United Kingdom

###### Abstract

The boundaries of quantum materials can host a variety of exotic effects such
as topologically robust edge states or anyonic quasiparticles. Here, we show
that fermionic systems such as graphene that admit a low energy Dirac
description can exhibit counterintuitive relativistic effects at their
boundaries. As an example, we consider carbon nanotubes and demonstrate that
relativistic bulk spinor states can have non zero charge density on the
boundaries, in contrast to the sinusoidal distribution of non-relativistic
wave functions that are necessarily zero at the boundaries. This unusual
property of relativistic spinors is complementary to the linear energy
dispersion relation exhibited by Dirac materials and can influence their
coupling to leads, transport properties or their response to external fields.

Introduction:– Several materials have low-energy quantum properties that are
faithfully described by the relativistic Dirac equation. The celebrated
example of graphene owes some of its unique properties, such as the half-
integer quantum Hall effect [1, 2, 3] and the Klein paradox effect [4, 5], to
the relativistic linear dispersion relation describing its low-energy sector.
This is by no means a singular case. A wide range of materials have been
recently identified that admit 1D, 2D or 3D relativistic Dirac description,
including many topological insulators and $d$-wave superconductors [6, 7, 8,
9, 10, 11]. The unusual dispersion relation of Dirac materials gives rise to
effective spinors, where the sublattice degree of freedom is encoded in the
pseudo-spin components. Nevertheless, the emerging excitations are spinor
quasiparticles that can exhibit novel transport properties or responses to
external fields akin only to relativistic physics [5].

Here we present another counter-intuitive aspect of relativistic physics in
Dirac materials manifested by the behaviour of bulk states at the boundaries.
In general, the choice of boundary conditions one imposes on single-particle
wavefunctions of a system must ensure its Hamiltonian remains Hermitian. For
the example of a non-relativistic particle in a box obeying the Schrödinger
equation, the boundary conditions are simply that the wavefunction vanishes on
the walls of the box. However, for spin-$1/2$ particles of mass $m$ obeying
the $(2+1)$D Dirac equation

$\left(-i\alpha^{i}\partial_{i}+\beta
m\right)\psi(\bm{r})=E\psi(\bm{r}),\,\,\psi(\bm{r})=\begin{pmatrix}\psi_{\uparrow}(\bm{r})\\\
\psi_{\downarrow}(\bm{r})\end{pmatrix},$ (1)

where $\alpha^{i}$ and $\beta$ are the $2\times 2$ Dirac alpha and beta
matrices, vanishing of the spinor $\psi(\bm{r})$ is not possible on all
boundaries without the solution being trivially zero everywhere. The
requirement that the Dirac Hamiltonian $h=-i\alpha^{i}\partial_{i}+\beta m$ is
Hermitian with respect to the inner product
$\int_{D}\mathrm{d}^{2}r\psi^{\dagger}(\bm{r})\phi(\bm{r})$ on a finite domain
$D$ is that the charge current
$J^{i}(\bm{r})=\psi^{\dagger}(\bm{r})\alpha^{i}\psi(\bm{r})$ normal to the
boundary $\partial D$ is zero for all spinors. In other words, if
$\hat{\bm{n}}$ is the outward pointing normal to the boundary, then

$\hat{\bm{n}}(\bm{r}_{0})\cdot\bm{J}(\bm{r}_{0})=0$ (2)

for all points $\bm{r}_{0}\in\partial D$ [12]. This condition ensures that
particles are trapped in $D$. In contrast to the non-relativistic case, the
zero flux condition of Eq. (2) allows for bulk solutions $\psi(\bm{r})$ whose
charge density $\rho(\bm{r}_{0})=\psi^{\dagger}(\bm{r}_{0})\psi(\bm{r}_{0})$
is non-zero on the boundaries [13, 14].

To exemplify our investigation, we consider how bulk spinor states behave at
the edges of a zig-zag carbon nanotube – a system which is described by the
Dirac equation of Eq. (1). We find that bulk states have support on the edges
of the nanotube depending on the size of the system. Importantly, these
relativistic effects become more dominant for gapless nanotubes, corresponding
to systems with a multiple of three unit cells in circumference, or when the
length of the nanotube is small. Such relativistic properties of spinor
eigenstates are expected to be present in all Dirac-like materials and are
complementary to the typically linear dispersion relation they exhibit. Bulk
states with non-zero density at the boundaries are expected to impact the
coupling of Dirac materials to external leads, their transport properties or
their response to external magnetic fields.

Relativistic description of zig-zag carbon nanotubes:– The honeycomb lattice
of graphene is formed from two triangular sublattices $A$ and $B$. We take the
two basis vectors $\bm{n}_{x}$,$\bm{n}_{y}$ and we take the vertical links as
our unit cells, as shown in Fig. 1. We label our lattice sites with the pair
$(\bm{r},\mu)$, where $\bm{r}=x\bm{n}_{x}+y\bm{n}_{y}$ labels the position of
the unit cell with non-Cartesian coordinates $x,y\in\mathbb{Z}$, while
$\mu\in\\{A,B\\}$ labels the site within the unit cell. The Hamiltonian of the
system is given by
$H=-t\sum_{\langle\bm{r},\bm{r}^{\prime}\rangle}a^{\dagger}_{\bm{r}}b_{\bm{r}^{\prime}}+\text{h.c.}$,
where $a^{\dagger}_{\bm{r}}$ ($b^{\dagger}_{\bm{r}}$) creates a fermion on
sublattice $A$ ($B$) of unit cell ${\bm{r}}$ [5]. Bloch momenta are given by
$\bm{k}=\frac{1}{2\pi}(k_{x}\bm{G}_{x}+k_{y}\bm{G}_{y})$, where $\bm{G}_{x}$,
$\bm{G}_{y}$ are the reciprocal basis vectors and $k_{x},k_{y}\in[-\pi,\pi]$
are the corresponding coordinates of the Brillouin zone (BZ) (see Appendix).

Figure 1: (Left) The honeycomb lattice comprising of two triangular
sublattices $A$ and $B$. The lattice basis vectors
$\bm{n}_{x}=\frac{1}{2}(1,-\sqrt{3})$ and $\bm{n}_{y}=(1,0)$ are depicted with
corresponding non-Cartesian coordinates $x$ and $y$. A nanotube with zig-zag
boundaries (red lines) is depicted with length $L$ corresponding to $L+1$ unit
cells (dashed ovals) in the $\bm{n}_{x}$ direction while the $\bm{n}_{y}$
direction is periodic with circumference $N$. The boundary condition is given
by the top $A$ sites and bottom $B$ sites having zero population. (Right) The
band structure of a zig-zag nanotube of circumference $N=10$ is displayed,
where $n=3$ (dashed lines) is one of the two minimum gap bands. All bands have
a single minimum except $n=N/2=5$ which is completely flat (dotted line).

To study the low-energy properties of a finite zig-zag nanotube, we first take
the continuum and thermodynamic limit in the the $\bm{n}_{x}$ direction only,
whilst keeping the the periodic $\bm{n}_{y}$ direction finite and discrete,
with $N$ unit cells in circumference. This gives rise to $N$ bands
parametrised by momenta $k_{y}=2n\pi/N$, where $n$ is an integer [15, 16]. The
$n$th band has the one-dimensional dispersion relation

$E_{n}(k_{x})=\pm t\sqrt{3+2g_{n}(k_{x})},$ (3)

where
$g_{n}(k_{x})=\cos\left(\frac{2n\pi}{N}\right)+\cos\left(\frac{2n\pi}{N}-k_{x}\right)+\cos(k_{x})$.
The zig-zag nanotube is typically gapped, unlike an infinite flat sheet of
graphene which is gapless. Each conduction band contains a single minima, as
seen in Fig. 1, which dictates the low-energy physics for that particular
band. Our model is a simplified version of a carbon nanotube as we ignore
effects due to curvature and spin orbit coupling that are not relevant to our
investigation [17, 18, 19, 20, 21, 5, 16].

Following the literature [19, 21], we expand the Hamiltonian about the minima
of Eq. (3) by letting $k_{x}=K_{\text{min}}+p$, for each band $n$, yielding
the $(1+1)$D massive Dirac Hamiltonian
$H_{n}=\int\mathrm{d}x\psi^{\dagger}_{n}h_{n}\psi_{n}$, with

$h_{n}=-ie_{n}\sigma^{y}\partial_{x}+\Delta_{n}\sigma^{x},\,\,\psi_{n}(x)=\begin{pmatrix}a_{n}(x)\\\
b_{n}(x)\end{pmatrix},$ (4)

where the sublattices $A$ and $B$ of the unit cell are encoded on the pseudo-
spin components, $e_{n}=2t\cos\left(\frac{n\pi}{N}\right)$ is the spatial
component of the zweibein and $\Delta_{n}$ is the energy gap of the $n$th
band, given by

$\Delta_{n}=t\left(2\cos\left(\frac{n\pi}{N}\right)-1\right).$ (5)

See Appendix A for a derivation. We now truncate the length of the nanotube to
a finite length $L$. We construct standing waves
$\psi(x)=(\psi_{A}(x),\psi_{B}(x))^{\text{T}}$ from forward and backward
propagating eigenstates of $h_{n}$ of Eq. (4). The zig-zag boundary conditions
are $\psi_{A}(0)=\psi_{B}(L)=0$, where $x=0$ and $x=L$ are the coordinates of
the unit cells of the top and bottom boundaries, as shown in Fig. 1 [19, 20].
These conditions obey the zero flux condition of Eq. (2). This gives the
solutions

$\displaystyle\psi_{n,p}(x)$
$\displaystyle=\mathcal{N}_{n,p}\begin{pmatrix}\sin(px)\\\
\sin(px+\theta_{n,p})\end{pmatrix},$ (6a) $\displaystyle\theta_{n,p}$
$\displaystyle=\arg(\Delta_{n}+ie_{n}p),$ (6b)

where $\mathcal{N}_{n,p}$ is a normalisation constant and $\theta_{n,p}$ is a
relative phase shift between the $A$ and $B$ sublattice wavefunctions. The
quantised momenta $p$ are solutions to the transcendental equation

$pL+\theta_{n,p}=m\pi,\ m\in\mathbb{N},$ (7)

which can be solved numerically (see Appendix A). Note that the wavefunctions
of Eq. (6a) correspond to bulk states, however graphene with zig-zag
boundaries also supports zero-energy states localised at the edges [11]. Edge
states correspond to complex solutions of Eq. (7) and are not considered here
[5, 19, 20].

Relativistic edge effects of bulk states:– The $U(1)$ electric charge density
of $(1+1)$D Dirac spinors $\psi=(\psi_{A},\psi_{B})^{T}$ is given by
$\rho=\psi^{\dagger}\psi=|\psi_{A}|^{2}+|\psi_{B}|^{2}$. With our
interpretation of the pseudo-spin components $\psi_{A}(x)$ and $\psi_{B}(x)$
as the sublattice wavefunctions, where $x$ labels the unit cell, $\rho$ is
therefore the charge density with respect to the unit cells. For the bulk
standing wave solutions of Eq. (6a), we have

$\rho_{n,p}(x)=|\mathcal{N}_{n,p}|^{2}\left(\sin^{2}(px)+\sin^{2}(px+\theta_{n,p})\right),$
(8)

which gives a charge density at the boundaries of

$\rho_{n,p}(0)=\rho_{n,p}(L)=|\mathcal{N}_{n,p}|^{2}\sin^{2}(\theta_{n,p}).$
(9)

We see that it is possible to have $\rho_{n,p}(0)\neq 0$ due to the phase
difference, $\theta_{n,p}$, which is purely a relativistic effect.

The edge charge density of bulk states is maximal when
$\theta_{n,p}=\pm\pi/2$. Referring to Eq. (6b), this is achieved when
$\Delta_{n}=0$, i.e., when the $n$th band is gapless. From Eq. (5) we see that
the gap closes if $n/N=\pm 1/3$ which is only possible if $N$ is a multiple of
three. Note that, for a gapless band, the charge density of Eq. (8) is also
completely uniform with

$\rho_{n,p}(x)=\frac{1}{L},$ (10)

which is independent of the momentum $p$, where we have chosen a $1$D
normalisation. On the other hand, when the system is gapped, then the density
oscillates along the length of the nanotube and becomes vanishingly small at
the edges. This shift in behaviour of the charge density reflects the expected
transition from the relativistic to non-relativistic regime witnessed in
confined Dirac particles as their mass increases [13].

Figure 2: (Left) The phase shifts $\theta\equiv\theta_{n,p}$ of Eq. (6b) and
the numerically measured edge densities $\rho(0)\equiv\rho_{n,p}(0)$ versus
circumference $N$ for the ground state of a system of fixed length $L=200$.
When $N$ is a multiple of three, i.e., when the system is gapless, the phase
shift is $\theta=\pi/2$ and the edge density is $\rho(0)=1/L=0.005$ confirming
Eq. (10). (Right) The phase shifts $\theta$ and edge densities $\rho(0)$
versus system length $L$ for the ground state of the gapless systems $N=30$
and its two neighbouring gapped systems $N=29$ and $N=31$. The solid line
represents the analytical formulas whilst the points represent numerics. The
edge density for the gapless system $N=30$ goes as $1/L$, while for gapped
systems $N=29$ and $N=31$ the edge density tends to zero quickly in agreement
with Eq. (9). (Inset) The integrated LDOS at the edge $x=0$ for a nanotube of
circumferences $N=29,30,31$ and lengths $L=25$ (dashed lines) and $L=100$
(solid lines). The LDOS displays the predicted behaviour of maximising for
gapless systems ($N=30$) and increasing with smaller system size $L$.

The stark contrast between gapped and gapless systems is confirmed numerically
(see Appendix B for numerical details). The left-hand column of Fig. 2 shows
the edge density $\rho$ of the ground state of a system of length $L=200$ for
varying circumferences $N$. When $N$ is a multiple of three, i.e., when the
system is gapless, the edge density spikes to the expected value of
$1/L=0.005$. On the other hand, when $N$ is not a multiple of three, i.e.,
when the system is gapped, the edge density is small. This behaviour is a
consequence of the highly oscillating phase shift $\theta$. When $N$ is a
multiple of three, the phase shift is $\pi/2$ exactly, maximising the edge
density according to Eq. (9). The smaller $N$ is, the stronger the effect as
the difference between gapped and gapless systems is much stronger due to the
gaps being larger. However, as $N$ increases, all zig-zag nanotubes tend
towards gapless systems even if $N$ is not a multiple of three, as there
exists a band $n$ such that $n/N\approx\pm 1/3$ when $N$ is large, so the gap
of Eq. (5) begins to close, so all systems begin to behave similarly.

The relativistic boundary effects also have a system length dependence [13].
The right-hand column of Fig. 2 shows the numerically measured edge density
$\rho$ of the ground state of a the gapless system $N=30$ and its two
neighbouring gapped systems $N=29$ and $N=31$ for varying system lengths $L$.
The edge density of the gapless system $N=30$ goes as $1/L$ whereas the edge
density gapped systems $N=29$ and $N=31$ tends to zero quickly, both in
accordance with Eq. (9). It is worth noting that, despite the fact that the
analytic results have been derived in the large $L$ limit where the continuum
approximation holds, the numerics and analytics are in surprisingly good
agreement even for very small $L$. This verifies the theoretically predicted
relativistic effects of nanotubes with small length $L$ where the violation of
the non-relativistic zero edge density is expected. Note that this behaviour
repeats itself for any $N$ that is a multiple of three and its two
neighbouring sizes above and below, which the left hand side of Fig. 2
demonstrates.

To explain the system size dependence of the charge density, note that for
very small $L$ the allowed momenta $p$ satisfying Eq. (7) become very large.
In this case, the imaginary contribution to the phase
$\theta_{n,p}=\arg(\Delta_{n}+ie_{n}p)$ dominates, giving
$\theta_{n,p}\approx\pi/2$ even if the gap is non-zero, as seen in the right-
hand column of Fig. 2. Hence, the edge density of Eq. (9) becomes significant
for small system sizes. For the gapless case, the phase is exactly equal to
$\pi/2$ regardless of the value of $p$ or system size $L$. This yields a
uniform charge density throughout the nanotube, resulting in the $1/L$ edge
density as observed.

Finally, Fig. 2 shows the integrated local density of states (LDOS) on the
edge at $x=0$ given by $N(E,\bm{r})=\sum_{m}\rho_{m}(\bm{r})\Theta(E-E_{m})$,
where $\rho_{m}(\bm{r})$ is the unit cell charge density of the $m$th
eigenstate of the $2$D model with eigenvalue $E_{m}$. We present this for
systems $N=29,30,31$ and $L=25,100$. The edge LDOS is maximised for a fixed
$L$ when the system is gapless, so for $N=30$ in this case. Moreover, the LDOS
increases as the system size decreases, which provides a clear signature for
the observation of the relativistic edge effect.

To summarise, the edge density is prominent if either the system is gapless,
so $N$ is a multiple of three, or the system length $L$ is small. The typical
lattice constant of a nanotube is given by $|\bm{n}_{x}|=|\bm{n}_{y}|\approx
2.46$Å [22, 16], so Fig. 2 applies to systems on the order of $1$nm in
diameter and $10$nm in length. However, the dependence on whether the system
is gapless or not is very strong, so this effect holds for much larger
circumferences $N$ and lengths $L$. Therefore, we expect these results to hold
for a wide range of experimentally accessible sizes.

Relativistic spinors from non-relativistic wavefunctions:– To explain the
emergence of relativistic boundary effects from a non-relativistic model, we
focus on the sublattice wavefunctions $\psi_{A}$ and $\psi_{B}$. For
concreteness, we examine a nanotube of dimension $(N,L)=(30,200)$ and
$(N,L)=(31,200)$ which have gapless and gapped spectra, respectively.

Figure 3: (Left) A comparison of the analytical wavefunctions $|\psi_{A}|^{2}$
and $|\psi_{B}|^{2}$ of Eq. (6a) and charge densities $\rho$ of Eq. (9) to the
numerical simulation (markers) for the gapless system $(N,L)=(30,200)$. We
present the first three excited states above the Fermi energy $E=0$. We
observe the wavefunctions act highly-relativistically, with a large boundary
support and uniform charge density of $1/L=0.005$. (Right) We present the same
information for the gapped system $(N,L)=(31,200)$. The behaviour contrasts
highly with the gapless system despite $N$ only being greater by one. The
wavefunctions and densities display a Schrödinger-like profile with a much
smaller edge density that tends to zero as $L$ increases as seen in Fig. 2.

In the left-hand column of Fig. 3 we compare the numerical sublattice
wavefunctions $\psi_{A}(x)$, $\psi_{B}(x)$ and the charge densities $\rho(x)$
to the analytical results of Eq. (6a) and Eq. (8) respectively, for the first
three excited states above the Fermi energy for the gapless system
$(N,L)=(30,200)$. We see that the sublattice wavefunctions $\psi_{A}$ and
$\psi_{B}$ are highly out of phase and maximise the edge support at $x=L$ and
$x=0$ respectively, yielding a charge density $\rho(x)$ with minor
oscillations about the predicted uniform value of $1/L=0.005$. These
oscillations are caused by finite-size effects.

In the right-hand column of Fig 3, we present the same information for the
gapped system $(N,L)=(31,200)$. Despite $N$ increasing only by $1$, the fact
the system now has a gap results in wavefunctions $\psi_{A}(x)$ and
$\psi_{B}(x)$ that contrast considerably to the gapless case, with a charge
density $\rho(x)$ that displays a more Schrödinger-like oscillatory profile.
As the system size $L$ increases, the relative phase shift $\theta$ modolo
$\pi$ between $\psi_{A}(x)$ and $\psi_{B}(x)$ decreases, as seen in Fig. 2,
and the wavefunctions begin to display the Schrödinger-like profile that tends
to zero on the boundaries. However, this is not the case for gapless systems
as the phase shift is always $\pi/2$ regardless of system size, as seen in
Fig. 2.

Figure 4: The full single-particle ground state $\Psi_{0,j}$ and the first
state above the Fermi energy $\Psi_{j}$ for a system of size $(N,L)=(30,9)$.
The emergent relativistic physics near the Fermi energy can be seen clearly as
a result of aliasing of a high frequency Schrödinger wavefunction. Comparing
this with the left-hand column of Fig. 3, we see how the sublattice
wavefunctions $\psi_{A}$ and $\psi_{B}$ described by the spinor Eq. (6a)
emerges.

We now analyse the total wavefunctions $\Psi_{j}$ of the lattice fermions,
where $j\in\mathbb{N}$ is the real space coordinate of the bipartite lattice,
alternating between sublattices $A$ and $B$. This coordinate should be
contrasted to the unit cell coordinate $x$ of the spinor $\psi(x)$. Fig. 4
shows the wavefunctions of the single-particle eigenstate with the most
negative energy below the Fermi energy, $\Psi_{0,j}$, and the first single-
particle eigenstate above the Fermi energy, $\Psi_{j}$, for a system of
dimension $(N,L)=(30,9)$.

The wavefunctions $\Psi_{0,j}$ and $\Psi_{j}$ are both non-relativistic
wavefunctions which vanish on the boundaries. This is to be expected as the
microscopic model is non-relativistic. However, due to high frequency
oscillations, the support of $\Psi_{j}$ on each sublattice is highly out of
phase. Comparing with the left-hand column of Fig. 3, we see that these
oscillations give the impression of two separate wavefunctions faithfully
described by the components of a Dirac spinor. Non-relativistically, we expect
the system to behave like a particle in a box, so we take the ansatz
wavefunction $\Psi_{j}\propto\sin(pj)$. From inspection, we see that this
matches the numerics for momenta $p=(l+1)\pi/2l$, where $l=2L+1$ is the total
length of the bipartite chain and $L+1$ is the number of unit cells as defined
in Fig. 1. This gives a wavelength comparable to the lattice spacing.
Therefore, the emergent relativistic physics described by the spinor of Eq.
(6a) is a consequence of aliasing from sampling a high frequency non-
relativistic wavefunction at discrete intervals. This effect is independent of
length $L$. Such high frequency wavefunctions correspond to the middle of the
spectrum where the relativistic linear dispersion is present.

On the other hand, gapped systems display a Schrödinger-like wavefunction for
both $\Psi_{0,j}$ and $\Psi_{j}$ if the system length $L$ is large. This can
be seen clearly in the right-hand column of Fig. 3 where the sublattice
wavefunctions are almost in phase. The total wavefunctions that describe these
can also be described by the ansatz wavefunction of a particle in a box, but
for a small $p$ instead, so $A$ and $B$ sublattices are now more in phase,
similar to the left hand side of Fig. 4. We also see this in the right hand
side of Fig. 2 where the edge densities drop to zero on the walls, suggesting
a non-relativistic behaviour.

Conclusion:– Our analysis demonstrates that relativistic effects can dominate
certain geometries of Dirac materials, resulting in large edge support. We
studied this effect analytically and numerically for zig-zag carbon nanotubes
and demonstrated that it holds strongly for a wide range of experimentally
accessible sizes. We found that the effect is dominant when the system is
either gapless or has a small length on the order of $10$nm. Nevertheless,
this relativistic effect is general and it is expected to be present in 1D, 2D
and 3D materials with the same qualitative properties presented here. While
high edge densities of bulk states should be measurable with STM [23, 24, 25,
26], it is expected to have a significant effect on the conductivity of the
material when attaching leads to its boundaries or its response to a magnetic
field [27, 19, 20, 28]. In addition, determining if such effects will be
present in 2D materials containing a finite density of defects which
effectively imposes boundary conditions on the wavefunctions within the
material will be intriguing [29, 30, 5, 31]. We leave these questions for a
future work.

###### Acknowledgements.

Acknowledgements:– We would like to thank Oscar Cespedes, Jamie Lake, Alex
Little and Satoshi Sasaki for inspiring conversations. This work was supported
by the EPSRC grant EP/R020612/1. Statement of compliance with EPSRC policy
framework on research data: This publication is theoretical work that does not
require supporting research data.

## References

  * Novoselov  _et al._ [2005] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos, and A. A. Firsov, Two-dimensional gas of massless dirac fermions in graphene, Nature 438, 197–200 (2005).
  * Novoselov  _et al._ [2007] K. S. Novoselov, Z. Jiang, Y. Zhang, S. V. Morozov, H. L. Stormer, U. Zeitler, J. C. Maan, G. S. Boebinger, P. Kim, and A. K. Geim, Room-temperature quantum hall effect in graphene, Science 315, 1379–1379 (2007).
  * Fujita and Suzuki [2016] S. Fujita and A. Suzuki, Theory of the half-integer quantum hall effect in graphene, International Journal of Theoretical Physics 55, 4830 (2016).
  * Katsnelson  _et al._ [2006] M. I. Katsnelson, K. S. Novoselov, and A. K. Geim, Chiral tunnelling and the klein paradox in graphene, Nature Physics 2, 620–625 (2006).
  * Castro Neto  _et al._ [2009] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, The electronic properties of graphene, Rev. Mod. Phys. 81, 109 (2009).
  * Wehling  _et al._ [2014] T. Wehling, A. Black-Schaffer, and A. Balatsky, Dirac materials, Advances in Physics 63, 1–76 (2014).
  * Moore [2010] J. E. Moore, The birth of topological insulators, Nature 464, 194 (2010).
  * Jia  _et al._ [2016] S. Jia, S.-Y. Xu, and M. Z. Hasan, Weyl semimetals, fermi arcs and chiral anomalies, Nature Materials 15, 1140 (2016).
  * Hasan and Kane [2010] M. Z. Hasan and C. L. Kane, Colloquium: Topological insulators, Reviews of Modern Physics 82, 3045–3067 (2010).
  * Hasan and Moore [2011] M. Z. Hasan and J. E. Moore, Three-dimensional topological insulators, Annual Review of Condensed Matter Physics 2, 55 (2011).
  * Bernevig and Hughes [2013] B. A. Bernevig and T. L. Hughes,  _Topological Insulators and Topological Superconductors_ (Princeton University Press, 2013) pp. 80–90.
  * Berry and Mondragon [1987] M. V. Berry and R. Mondragon, Neutrino billiards: time-reversal symmetry-breaking without magnetic fields, Proc. R. Soc. Lond. A 412, 53 (1987).
  * Alberto  _et al._ [1996] P. Alberto, C. Fiolhais, and V. M. S. Gil, Relativistic particle in a box, European Journal of Physics 17, 19 (1996).
  * Alonso  _et al._ [1997] V. Alonso, S. D. Vincenzo, and L. Mondino, On the boundary conditions for the dirac equation, European Journal of Physics 18, 315 (1997).
  * Charlier  _et al._ [2007] J.-C. Charlier, X. Blase, and S. Roche, Electronic and transport properties of nanotubes, Rev. Mod. Phys. 79, 677 (2007).
  * Saito  _et al._ [1998] R. Saito, G. Dresselhaus, and M. S. Dresselhaus,  _Physical Properties of Carbon Nanotubes_ (Imperial College Press, 1998).
  * Kane and Mele [1997] C. L. Kane and E. J. Mele, Size, shape, and low energy electronic structure of carbon nanotubes, Phys. Rev. Lett. 78, 1932 (1997).
  * Ando [2000] T. Ando, Spin-orbit interaction in carbon nanotubes, Journal of the Physical Society of Japan 69, 1757 (2000).
  * Margańska  _et al._ [2019] M. Margańska, D. R. Schmid, A. Dirnaichner, P. L. Stiller, C. Strunk, M. Grifoni, and A. K. Hüttel, Shaping electron wave functions in a carbon nanotube with a parallel magnetic field, Phys. Rev. Lett. 122, 086802 (2019).
  * Margańska  _et al._ [2011] M. Margańska, M. del Valle, S. H. Jhang, C. Strunk, and M. Grifoni, Localization induced by magnetic fields in carbon nanotubes, Phys. Rev. B 83, 193407 (2011).
  * Akhmerov and Beenakker [2008] A. R. Akhmerov and C. W. J. Beenakker, Boundary conditions for dirac fermions on a terminated honeycomb lattice, Phys. Rev. B 77, 085423 (2008).
  * Altland and Simons [2010] A. Altland and B. D. Simons,  _Condensed Matter Field Theory_ (Cambridge University Press, 2010) pp. 55–58.
  * Andrei  _et al._ [2012] E. Y. Andrei, G. Li, and X. Du, Electronic properties of graphene: a perspective from scanning tunneling microscopy and magnetotransport, Reports on Progress in Physics 75, 056501 (2012).
  * Kim  _et al._ [2000] P. Kim, T. W. Odom, J. Huang, and C. M. Lieber, Stm study of single-walled carbon nanotubes, Carbon 38, 1741 (2000), fullerenes ’99.
  * Venema  _et al._ [2000] L. C. Venema, V. Meunier, P. Lambin, and C. Dekker, Atomic structure of carbon nanotubes from scanning tunneling microscopy, Phys. Rev. B 61, 2991 (2000).
  * Hassanien  _et al._ [1998] A. Hassanien, M. Tokumoto, Y. Kumazawa, H. Kataura, Y. Maniwa, S. Suzuki, and Y. Achiba, Atomic structure and electronic properties of single-wall carbon nanotubes probed by scanning tunneling microscope at room temperature, Applied Physics Letters 73, 3839 (1998).
  * Laird  _et al._ [2015] E. A. Laird, F. Kuemmeth, G. A. Steele, K. Grove-Rasmussen, J. Nygård, K. Flensberg, and L. P. Kouwenhoven, Quantum transport in carbon nanotubes, Rev. Mod. Phys. 87, 703 (2015).
  * Ajiki and Ando [1993] H. Ajiki and T. Ando, Electronic states of carbon nanotubes, Journal of the Physical Society of Japan 62, 1255 (1993).
  * Algharagholy [2019] L. A. Algharagholy, Defects in carbon nanotubes and their impact on the electronic transport properties, Journal of Electronic Materials 48, 2301 (2019).
  * Araujo  _et al._ [2012] P. T. Araujo, M. Terrones, and M. S. Dresselhaus, Defects and impurities in graphene-like materials, Materials Today 15, 98 (2012).
  * Dutreix  _et al._ [2019] C. Dutreix, H. González-Herrero, I. Brihuega, M. I. Katsnelson, C. Chapelier, and V. T. Renard, Measuring the berry phase of graphene from wavefront dislocations in friedel oscillations, Nature 574, 219 (2019).

##  I Appendix A: Continuum limit of zig-zag carbon nanotubes

The honeycomb lattice of graphene is formed from a triangular Bravais lattice
with a unit cell containing two sites, one on sublattice $A$ and the other on
sublattice $B$, as shown in Fig. 1 of the main text. The Bravais lattice is
generated by the two basis vectors

$\bm{n}_{x}=\frac{a}{2}(1,-\sqrt{3}),\ \bm{n}_{y}=a(1,0),$ (11)

where $|\bm{n}_{x}|=|\bm{n}_{y}|=a$ is the lattice spacing. In this Letter, we
take $a=1$, but we leave it in this supplementary material for completeness.
We label our lattice sites with the pair $(\bm{r},\mu)$, where
$\bm{r}=x\bm{n}_{x}+y\bm{n}_{y}$ labels the position of the unit cell which we
take to coincide with sublattice $A$, where $x,y\in\mathbb{Z}$ are the non-
Cartesian coordinates, and $\mu\in\\{A,B\\}$ labels the site within the unit
cell. The corresponding reciprocal basis is given by

$\bm{G}_{x}=\frac{2\pi}{a\sqrt{3}}(0,-2),\
\bm{G}_{y}=\frac{2\pi}{a\sqrt{3}}(\sqrt{3},1).$ (12)

where $\bm{n}_{i}\cdot\bm{G}_{j}=2\pi\delta_{ij}$. With this reciprocal basis,
the Bloch momenta are given by
$\bm{k}=\frac{a}{2\pi}(k_{x}\bm{G}_{x}+k_{y}\bm{G}_{y})$, where
$k_{x},k_{y}\in[-\pi/a,\pi/a]$ defines the Brillouin zone (BZ) which is square
in this coordinate system. The components of momenta in the $\bm{n}_{i}$
direction are given by $k_{i}=\bm{k}\cdot\hat{\bm{n}}_{i}$, where
$\hat{\bm{n}}_{i}=\bm{n}_{i}/a$ are unit vectors in these directions.

Figure 5: The nearest-neighbour hoppings of the Hamiltonian can be tiled
across the honeycomb with this particular unit cell of hoppings along $X$, $Y$
and $Z$ links represented by the thick black links. The unfilled circles
represent sites that are not contained in this choice of unit cell.

We take the tight-binding Hamiltonian of graphene with nearest-neighbour
hoppings only. With our choice of unit cell, basis vectors and labelling
convention, the tight-binding Hamiltonian can be written as

$\displaystyle H$
$\displaystyle=-t\sum_{\langle\bm{r},\bm{r}^{\prime}\rangle}a^{\dagger}_{\bm{r}}b_{\bm{r}^{\prime}}+\text{h.c.}$
(13) $\displaystyle=-t\sum_{{\bm{r}}\in
A}a^{\dagger}_{\bm{r}}\left(b_{\bm{r}}+b_{\bm{r}-{\bm{n}}_{x}}+b_{\bm{r}-{\bm{n}}_{x}+{\bm{n}}_{y}}\right)+\text{h.c.}$

where $t$ is the hopping parameter and $a^{\dagger}_{\bm{r}}$
($b^{\dagger}_{\bm{r}}$) are fermionic operators which create an electron on
sublattice $A$ ($B$) of the unit cell located at $\bm{r}$ [5]. The second
equality corresponds to the fact one can construct the honeycomb lattice by
tiling sublattice $A$ with the “Y” shape of links that originate from a single
$A$ site. Careful consideration must be taken with the top and bottom row when
we have boundaries as the external vertical links are missing, as shown in
Fig. 5.

We impose periodic boundary conditions upon $H$ in the $\bm{n}_{x}$ and
$\bm{n}_{y}$ directions, with $L+1$ and $N$ unit cells in these directions
respectively, which gives us a carbon nanotube of length $L$. Therefore, the
Hamiltonian can be diagonalised by first taking the discrete Fourier transform

$a_{\bm{r}}=\frac{1}{\sqrt{N_{c}}}\sum_{\bm{k}\in\text{BZ}}e^{i\bm{k}\cdot\bm{r}}a_{\bm{k}}$
(14)

and similarly for $b_{\bm{r}}$, where $N_{c}=N(L+1)$ is the number of unit
cells in the honeycomb lattice. With this, the Hamiltonian takes the form

$H=\sum_{\bm{k}\in\text{BZ}}\psi^{\dagger}_{\bm{k}}h(\bm{k})\psi_{\bm{k}},\
h(\bm{k})=\begin{pmatrix}0&f(\bm{k})\\\ f^{*}(\bm{k})&0\end{pmatrix},$ (15)

where $\psi_{\bm{k}}=(a_{\bm{k}},b_{\bm{k}})^{T}$ is a two-dimensional spinor,
where the sublattice degrees of freedom appear as the “spin” degrees of
freedom of the spinor, and

$f(\bm{k})=-t\left(1+e^{-i(k_{x}-k_{y})a}+e^{-ik_{x}a}\right).$ (16)

The single-particle dispersion relation of $H$ is given by
$E(\bm{k})=\pm|f(\bm{k})|$ which is shown in the 2D colour plot of Fig. 6.
This dispersion is gapless and contains two zero energy Dirac points about
which the system acts relativistically, as shown by the crosses in Fig. 6.

To construct a zig-zag nanotube, we let the length $L$ in the $\bm{n}_{x}$
direction tend to infinity whilst letting the periodic periodic length in the
$\bm{n}_{y}$ direction remain finite, with $N$ unit cells in circumference,
which creates an infinitely long nanotube. In this case, the Bloch momenta are
semi-quantised within the BZ, with $k_{x}$ unconstrained and

$k_{y}=\frac{2n\pi}{Na},$ (17)

where $n$ is an integer. The quantisation of $k_{y}$ means that the system
only has access to one-dimensional bands of momentum states within the BZ
labelled by the integer $n$ [15]. The $n$th band has the dispersion relation

$E_{n}(k_{x})=\pm t\sqrt{3+2g_{n}(k_{x})},$ (18)

where

$g_{n}(k_{x})=\cos\left(\frac{2n\pi}{N}\right)+\cos\left(\frac{2n\pi}{N}-k_{x}a\right)+\cos(k_{x}a),$
(19)

which is obtained by simply substituting the quantised values of $k_{y}$ into
the dispersion $E(\bm{k})$ of graphene. For each value of $n$ we have an
energy band which is gapped in general, unlike an infinite flat sheet of
graphene which is always gapless. This is due to the finite circumference of
$N$ unit cells. Each conduction band contains a single minima which will
describe the low-energy physics for that particular band. We stress that these
minima are not the two zero-energy Dirac points of an infinite sheet of
graphene, as these two points are inaccessible to the system in general. Only
for special values of $N$ will these points be accessible, yielding a gapless
nanotube.

Figure 6: The Brillouin zone for a system zig-zag nanotube of $N=10$ with
lattice spacing $a=1$. The horizontal lines are the paths each band
$E_{n}(k_{x})$ takes through the BZ. The two red crosses mark the positions of
the zero energy Dirac points for an infinite sheet of graphene. We see that,
for this particular nanotube, the bands do not pass through the Dirac points,
resulting in a gapped system. The $n$th band will pass through a Dirac point,
resulting in a gapless system, if $n/N=\pm 1/3$.

We study the low-energy properties for the $n$th band by Taylor expanding the
Hamiltonian about the band minima, following a similar route to that of
references [19, 21]. These minima are located at the same position as the
minima of $F_{n}(k_{x})=E_{n}^{2}(k_{x})$, so we use this as it is easier to
work with. The turning points $K_{\text{min}}$ obey
$F_{n}^{\prime}(K_{\text{min}})=0$, which gives us the equation

$\sin\left(\frac{2n\pi}{N}-K_{\text{min}}a\right)=\sin\left(K_{\text{min}}a\right).$
(20)

Using the result that $\sin(x)=\sin(y)\Rightarrow x=m\pi+(-1)^{m}y$ for
$m\in\mathbb{Z}$, this implies

$K_{\text{min}}=\frac{\pi}{(1+(-1)^{m})a}\left(\frac{2n}{N}-m\right).$ (21)

Due to the denominator we have finite solutions only when $m$ is an even
number, so we take $m=-2l$ which gives us the turning points

$K_{\text{min}}=\frac{\pi}{a}\left(\frac{n}{N}+l\right),\quad l\in\mathbb{Z}.$
(22)

Note that as $K_{\text{min}}\in[-\pi,\pi]$ and $n/N\in[-1/2,1/2]$, the only
possibilities are that $l=-1,0,1$. We now need to identify which of these
turning points are minima, where $F^{\prime\prime}_{n}(K_{\text{min}})>0$.
Substituting in $K_{\text{min}}$ from above, we require

$(-1)^{l}\cos\left(\frac{n\pi}{N}\right)<0.$ (23)

As $n/N\in[-1/2,1/2]$, then $\cos(\frac{n\pi}{N})\in[0,1]$, therefore in order
to satisfy this constraint whilst ensuring $K_{\text{min}}\in[-\pi,\pi]$, we
require

$l=\begin{cases}-1&\text{if}\ n>0\\\ 1&\text{if}\ n\leq 0\end{cases}\
\Rightarrow\
K_{\text{min}}=\begin{cases}\left(\frac{n}{N}-1\right)\frac{\pi}{a}&n>0\\\
\left(\frac{n}{N}+1\right)\frac{\pi}{a}&n\leq 0\end{cases}.$ (24)

Each band has a single minima. If $n=N/2$, there does not exist a minima as
the band is completely flat so we do not consider this.

The continuum limit Hamiltonian of the $n$th band is defined as
$h_{n}(p)=h(K_{\text{min}}+p,2n\pi/Na)$ to first order in $p$, where $h(p)$ is
given in Eq. (15). Therefore, substituting in $k_{y}=2n\pi/Na$ into
$f(\bm{p})$ gives us the set of $N$ functions

$f_{n}(p)\equiv
f(p,2n\pi/Na)=-t\left(1+e^{-ipa}\left(e^{i\frac{2n\pi}{N}}+1\right)\right),$
(25)

which are enumerated by the band index $n$. We Taylor expand $f_{n}(p)$ about
$K_{\text{min}}$:

$f_{n}(K_{\text{min}}+p)=f_{n}(K_{\text{min}})+pf^{\prime}_{n}(K_{\text{min}})+O(p^{2}).$
(26)

We have

$\displaystyle f_{n}(K_{\text{min}})$
$\displaystyle=-t\left(1+e^{-iK_{\text{min}}a}\left(e^{i\frac{2n\pi}{N}}+1\right)\right)$
(27) $\displaystyle=-t\left(1+e^{-i\frac{n\pi}{N}}e^{\mp
i\pi}\left(e^{i\frac{2n\pi}{N}}+1\right)\right)$
$\displaystyle=-t\left(1-e^{-i\frac{n\pi}{N}}\left(e^{i\frac{2n\pi}{N}}+1\right)\right)$
$\displaystyle=-t\left(1-\left(e^{i\frac{n\pi}{N}}+e^{-i\frac{n\pi}{N}}\right)\right)$
$\displaystyle=-t\left(1-2\cos\left(\frac{n\pi}{N}\right)\right)$
$\displaystyle\equiv\Delta_{n}.$

We also have

$\displaystyle f^{\prime}_{n}(p)$
$\displaystyle=iate^{-ipa}\left(e^{i\frac{2n\pi}{N}}+1\right)$ (28)
$\displaystyle\Rightarrow f^{\prime}_{n}(K_{\text{min}})$
$\displaystyle=iate^{-iK_{\text{min}}a}\left(e^{i\frac{2n\pi}{N}}+1\right)$
$\displaystyle=iate^{-i\frac{n\pi}{N}}e^{\mp
i\pi}\left(e^{i\frac{2n\pi}{N}}+1\right)$
$\displaystyle=-iat\left(e^{i\frac{n\pi}{N}}+e^{-i\frac{n\pi}{N}}\right)$
$\displaystyle=-2iat\cos\left(\frac{n\pi}{N}\right)$ $\displaystyle\equiv-
ie_{n}.$

Pulling everything together, we get

$h_{n}(p)=\begin{pmatrix}0&f_{n}(K_{\text{min}}+p)\\\
f_{n}^{*}(K_{\text{min}}+p)&0\end{pmatrix}=\begin{pmatrix}0&\Delta_{n}-ie_{n}p\\\
\Delta_{n}+ie_{n}p&0\end{pmatrix}=e_{n}\sigma^{y}p+\Delta_{n}\sigma^{x}+O(a^{2}p^{2}).$
(29)

We choose to interpret $e_{n}$ as a zweibein as it allows us to generalise to
curved spacetimes where $t$ is space-dependent. We see that the failure of the
continuum limit to describe the flat band $n=N/2$ is encoded in the zweibein
as it vanishes here. Taking the limit that $a\rightarrow 0$ whilst keeping
$at$ fixed, we can safely ignore the $O(a^{2}p^{2})$ terms and we now have our
continuum limit/low-energy Hamiltonian.

The above Hamiltonian is a $(1+1)$D Dirac Hamiltonian with mass $\Delta_{n}$
using the representation $\alpha^{x}=\sigma^{y}$ and $\beta=\sigma^{x}$. The
relativistic dispersion relation is given by
$E_{n}(p)=\pm\sqrt{\Delta^{2}_{n}+e_{n}^{2}p^{2}}$. Note the gap closes if
$n/N=\pm 1/3$, which is only possible if $N$ is a multiple of three. The
corresponding Dirac equation reads

$\begin{pmatrix}0&\Delta_{n}-ie_{n}p\\\
\Delta_{n}+ie_{n}p&0\end{pmatrix}\begin{pmatrix}\phi_{A}\\\
\phi_{B}\end{pmatrix}=E_{n}\begin{pmatrix}\phi_{A}\\\ \phi_{B}\end{pmatrix}.$
(30)

This yields two equivalent equations, both implying

$\phi_{B}=\left(\frac{\Delta_{n}+ie_{n}p}{E_{n}}\right)\phi_{A}\equiv
se^{i\theta_{n,p}}\phi_{A},$ (31)

where $\theta_{n,p}=\arg(\Delta_{n}+ie_{n}p)$ and $s=\mathrm{sgn}(E_{n})$. The
corresponding un-normalised eigenvectors in one-dimensional position space are
given by

$\phi_{n,p}(x)=\begin{pmatrix}1\\\ se^{i\theta_{n,p}}\end{pmatrix}e^{ipx},$
(32)

where we now we rename $ax\mapsto x$ which is our continuum coordinate system
when $a\rightarrow 0$. We interpret the top and bottom components of our
spinors as the wavefunction on sublattices $A$ and $B$ respectively.

With the continuum limit approximation, we now study a nanotube of finite
length $L$ in the $\bm{n}_{x}$ direction by imposing suitable boundary
conditions. Note that for the purposes of numerically encoding this, this
requires $L+1$ unit cells in the $\bm{n}_{x}$ direction. First, we build
positive energy ($s=1$) standing waves by superimposing forward and backward
propagating waves as $\psi_{n,p}=\phi_{n,p}+R\phi_{n,-p}$, where
$R\in\mathbb{C}$ as

$\psi_{n,p}=\begin{pmatrix}1\\\
e^{i\theta_{n,p}}\end{pmatrix}e^{ipx}+R\begin{pmatrix}1\\\
e^{-i\theta_{n,p}}\end{pmatrix}e^{-ipx},$ (33)

where $\theta_{n,-p}=-\theta_{n,p}$. The zig-zag boundary conditions are given
by $\psi_{A}(0)=\psi_{B}(L)=0$. These boundary conditions can be seen clearly
in Fig. 7 as the unit cells of the top and bottom row, where $x=0$ and $x=L$,
each contain a “missing” site that is outside of the system (recall that our
coordinates $x$ label the unit cell and not the individual sites). The non-
relativistic wavefunction must vanish on these sites so the corresponding
components of the spinor must vanish. Note that, in our representation of the
Dirac alpha and beta matrices, the zero-flux condition of Eq. (2) reads
$\mathrm{Im}(\psi^{*}_{A}\psi_{B})=0$ on the boundaries which the zig-zag
boundary conditions satisfy. The first boundary condition gives $R=-1$, so our
solutions take the form

$\psi_{n,p}(x)=\mathcal{N}_{n,p}\begin{pmatrix}\sin(px)\\\
\sin(px+\theta_{n,p})\end{pmatrix},$ (34)

where $\mathcal{N}_{n,p}$ is a normalisation constant. The second boundary
condition gives $\sin(pL+\theta_{n,p})=0$, giving the transcendental equation
for the allowed momenta

$pL+\theta_{n,p}=m\pi,\ m\in\mathbb{N},$ (35)

which can be solved numerically by minimising the function
$f_{nm}(p)=|pL+\theta_{n,p}-m\pi|$ with respect to $p$ for a fixed $n$,$m$.

##  II Appendix B: Numerical simulation of a zig-zag carbon nanotube

In order to numerically simulate the zig-zag carbon nanotube, we must modify
the Hamiltonian Eq. (13) slightly to take into account the open boundaries of
the system. We take the Hamiltonian

$H=-t\sum_{{\bf r}\in A}a^{\dagger}_{\bf r}\left(x_{\bf
r}b_{\bm{r}-{\bm{n}}_{x}}+y_{\bm{r}}b_{\bm{r}-{\bm{n}}_{x}+{\bm{n}}_{y}}+z_{\bm{r}}b_{\bm{r}}\right)+\text{h.c.},$
(36)

where $x_{\bm{r}},y_{\bm{r}},z_{\bm{r}}\in\\{0,1\\}$ are numerical factors
that take into account the top and bottom boundaries of the system. These
terms “switch off” the external $X$, $Y$ and $Z$ links of the Hamiltonian
respectively, as seen in Fig. 7, to ensure the nanotube has zig-zag boundaries
represented by the red links.

Figure 7: Due to our choice of unit cell in the Hamiltonian as seen in Fig. 5,
the numerical factors $x_{\bm{r}}$, $y_{\bm{r}}$ and $z_{\bm{r}}$ are required
to “switch off” the links that are outside of the system to create the red
zig-zag boundaries. At the top boundary where $x=0$, all grey $X$, $Y$ and $Z$
links are outside of the system so the numerical factors are $0$ here. Along
the bottom row at $x=L$, only the $Z$ link is switched off.

In order to fix which band we are in numerically, we derive the corresponding
tight-binding model. We Fourier transform with respect to the $\bm{n}_{y}$
direction only, with the definition

$a_{\mathbf{r}}=\frac{1}{\sqrt{N}}\sum_{p_{y}}e^{ip_{y}y}a_{x}(p_{y}).$ (37)

Substituting this into the tight-binding Hamiltonian Eq. (36) gives us

$H=\sum_{p_{y}}H(p_{y}),$ (38)

where our one-dimensional chains have the Hamiltonian

$H(p_{y})=-t\sum_{i=0}^{L}a^{\dagger}_{i}\left(z_{i}b_{i}+(y_{i}+x_{i}e^{ip_{y}})b_{i-1}\right)+\text{h.c.},$
(39)

where we have swtiched to the index $i\in\mathbb{N}$ to label the sites of our
one-dimensional chain.

In order to describe a nanotube, we now set $k_{y}=2n\pi/N$ to give us $N$
Hamiltonians $H_{n}\equiv H(2n\pi/N)$ which describe each band $n$ of the
nanotube. Note that numerically we take the lattice spacing $a=1$. In order to
encode this numerically, we need the single-particle Hamiltonian. We can write
this Hamiltonian as

$H_{n}=\sum_{i,j}a^{\dagger}_{i}(h_{n})_{ij}b_{j}+b_{j}^{\dagger}(h_{n}^{\dagger})_{ij}a_{i},$
(40)

where

$(h_{n})_{ij}=-t\left(z_{i}\delta_{ij}+\left(1+e^{i\frac{2n\pi}{N}}\right)\delta_{i-1,j}\right).$
(41)

Note that $x_{i}$ and $y_{i}$ are no longer needed here, as the
$\delta_{i-1,j}$ removes the $X$ and $Y$ links on the top of the cylinder as
required.

If we define the $2(L+1)$-dimensional spinor $\psi=(a_{0},a_{1},\ldots
a_{L};b_{0},b_{1},\ldots b_{L})^{T}\equiv(\bm{a},\bm{b})^{T}$, then the many-
body Hamiltonian can be written as

$H_{n}=\begin{pmatrix}\bm{a}^{\dagger}&\bm{b}^{\dagger}\end{pmatrix}\begin{pmatrix}0&h_{n}\\\
h^{\dagger}_{n}&0\end{pmatrix}\begin{pmatrix}\bm{a}\\\
\bm{b}\end{pmatrix}\equiv\psi^{\dagger}\mathcal{H}_{n}\psi$ (42)

Numerically, we diagonalise the matrix $2(L+1)\times 2(L+1)$-dimensional
matrix $\mathcal{H}_{n}$ which his our single-particle Hamiltonian. In our
basis, the first (last) $L+1$ components will be the wavefunctions $\psi_{A}$
($\psi_{B}$) of sublattice $A$ ($B$).

##  III Appendix C: Band-dependence of densities

The contrasting behaviour between gapped and gapless systems can also be
confirmed numerically by studying each band $n$ of the system, as shown in
Fig. 8. In this figure, the edge density of the ground state and first four
excited states of a particular band is plotted for all bands $n$ of a system
of size $(N,L)=(30,50)$, whereby ground state we mean the lowest energy state
of that particular band. We observe that if the nanotube is in an eigenstate
of the $n$th band where $n=\pm N/3=\pm 10$, which are the two gapless bands,
then the edge density takes the predicted value of $1/L=0.02$ for all
eigenstates. Away from these special bands a gap opens up, the edge density
falls off sharply and a dependence of the density on eigenstate emerges. The
edge density behaviour is a consequence of the phase shift,
$\theta_{n,p}=\arg(\Delta_{n}+ie_{n}p)$, shown in Fig. 8. If $n=\pm N/3\pm
10$, then $\Delta_{n}=0$ exactly, so $\theta_{n,p}=\pi/2$ for all $p$ and
consequently all eigenstates behave identically with an edge density of
$1/L=0.02$ according to Eq. (10) of the main text. Away from these points the
phase rapidly jumps to $0$ or $\pi$, which explains the sharpness of the
density peaks.

Figure 8:  The analytical phase shifts $\theta_{n,p}$ and numerical edge
densities $\rho_{n,p}(0)=\rho_{n,p}(L)$ for the ground state and first four
excited states of all $N$ bands $n$ of the system of dimension
$(N,L)=(30,50)$. We see that the edge density is maximised for bands $n=\pm
N/3=\pm 10$ coinciding with where the phase shift is $\theta_{n,p}=\pi/2$.