# Solitons induced by an in-plane magnetic field in rhombohedral multilayer
graphene

Max Tymczyszyn  Peter H. Cross  Edward McCann<EMAIL_ADDRESS>Department of Physics, Lancaster University, Lancaster, LA1 4YB, UK

###### Abstract

We model the influence of an in-plane magnetic field on the orbital motion of
electrons in rhombohedral graphene multilayers. For zero field, the low-energy
band structure includes a pair of flat bands near zero energy which are
localized on the surface layers of a finite thin film. For finite field, we
find that the zero-energy bands persist and that level bifurcations occur at
energies determined by the component of the in-plane wave vector $q$ that is
parallel to the external field. The occurrence of level bifurcations is
explained by invoking semiclassical quantization of the zero field Fermi
surface of rhombohedral graphite. We find parameter regions with a single
isoenergetic contour of Berry phase zero corresponding to a conventional
Landau level spectrum and regions with two isoenergetic contours, each of
Berry phase $\pi$, corresponding to a Dirac-like spectrum of levels. We write
down an analogous one-dimensional tight-binding model and relate the
persistence of the zero-energy bands in large magnetic fields to a soliton
texture supporting zero-energy states in the Su-Schrieffer-Heeger model. We
show that different states contributing to the zero-energy flat bands in
rhombohedral graphene multilayers in a large field, as determined by the wave
vector $q$, are localized on different bulk layers of the system, not just the
surfaces.

##  I Introduction

Advances in the production of thin films of rhombohedral multilayer graphene
(RMG) [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] recently culminated in the realization
of high-quality films with up to fifty layers [8]. Scanning tunneling
spectroscopy [1], magneto-Raman [2], and photoemission [1, 3] measurements
have confirmed the existence of flat bands near the Fermi surface, as
predicted theoretically [11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21].
Meanwhile, the observation of phases including ferromagnetism [22] and
superconductivity [23] in trilayers of rhombohedral graphene have been
attributed to flat bands with strong electronic interactions.

The presence of flat bands localized at the surfaces of RMG may be understood
by analogy with edge states within the bulk band gap of the one-dimensional
Su-Schrieffer-Heeger (SSH) model [24, 25, 26, 27], whereby the intra- and
interlayer couplings in RMG play the role of alternating hopping parameters in
the SSH model [18, 17], Fig. 1. It has been predicted [19, 2] that the bulk
band gap of RMG closes with layer number $N$ at wave vector
$q_{c}=\gamma_{1}/\hbar v$ near the K-points at the corner of the Brillouin
zone, where $\gamma_{1}$ is the interlayer hopping parameter and $v$ is the
intralayer velocity. In the limit of large $N$, i.e., bulk rhombohedral
graphite, the band gap is closed and zero energy states occupy a ‘Dirac
spiral’, which rotates, as a function of the perpendicular wave vector $k_{z}$
[28, 29, 30, 31], around a K-point at in-plane radius $q_{c}$.

The study of Landau level (LL) spectra and the integer quantum Hall effect
(QHE) in magnetic fields are key characterisation tools of graphene [32, 33]
and related nanomaterials. For out-of-plane (perpendicular) magnetic field of
magnitude as little as $B=1\,$T, Shubnikov–de Haas oscillations emerge in RMG
[8]. For slightly stronger fields, $B>3$T, discrete LL form [34, 35, 19],
leading to quantized Hall resistivity and the onset of the QHE observed
experimentally [8]. In-plane (parallel) magnetic fields have also been
explored in graphene-related systems [36, 37, 19, 38], specifically in
relation to the magnetic ratchet effect in bilayer graphene [39, 40],
superconductivity in trilayers [23, 41], and the energy spectrum of Bernal
stacked multilayers [42, 43].

Figure 1: Schematic side view of the unit cell for three layers of
rhombohedral multilayer graphene. Labels $A_{n}$ and $B_{n}$ denote the two
nonequivalent atomic sites on each layer, $\gamma_{0}$ and $\gamma_{1}$ are
the intralayer and interlayer hopping parameters, respectively. The in-plane
carbon-carbon bond length is $a/\sqrt{3}=1.42\AA$ where $a$ is the lattice
constant, while the interlayer distance $d=3.46\AA$. The out-of-plane $z$ axis
is shown on the left where $z_{n}$ is the coordinate of the
$n^{{\mathrm{th}}}$ layer.

Figure 2: Low-energy band structures for rhombohedral graphene with 20 layers,
with a visualisation of the corresponding semiclassical quantization showing
isoenergetic contours in bulk rhombohedral graphite in a plane perpendicular
to the magnetic field. (a) and (b) show band structures for no magnetic field
and $B=100$T respectively [44], they have been calculated through the
numerical diagonalization of the Hamiltonian (2). (c) shows an isoenergetic
tube (blue) centred around the Dirac nodal spiral, a plane (red) corresponding
to $q_{x}=q_{c}=\gamma_{1}/\hbar v$ intercepts the spiral (black), creating an
isoenergetic contour which is depicted in (d), indicating a singly degenerate
state (4-fold including spin and valley degeneracy). Figures (e) and (f) are
similar, except the plane now cuts through $q_{x}=0$, creating two contours
corresponding to a doubly degenerate state (8-fold including spin and valley
degeneracy) found at low energy for $q_{x}<q_{c}$ in (b).

In this paper, we model the influence of an in-plane magnetic field on the
orbital properties of electrons in RMG stacks with layer number $N\gg 1$,
using a tight binding model [11, 15, 19] and magnetic field incorporated via a
Peierls substitution [39]. We numerically study the effect of the magnetic
field on the band structure in the vicinity of the Dirac points. Fig. 2(a)
shows the band structure for zero field exhibiting two flat bands at zero
energy [11, 12, 13, 14, 15, 16, 17, 18, 19, 21], and Fig. 2(b) shows the
spectrum for large magnetic field where we plot energy levels as a function of
$q_{x}$ which is the component of the wave vector measured from the K point in
the direction of the applied field, $B_{x}$ [44]. We note two qualitative
features of the spectrum in an in-plane field: Firstly, there are a series of
bifurcations from doubly- to singly-degenerate energy levels and, secondly,
the zero-energy flat bands, observed at zero field, persist up to very high
field strengths.

For energies $\varepsilon$, we find that the bifurcations occur along a set of
points in the $\varepsilon$-$q_{x}$ plot, Fig. 2(b), given by

$\displaystyle|\varepsilon|\approx\hbar v(q_{c}-|q_{x}|),$ (1)

for $|q_{x}|<q_{c}$. The origin and position of the bifurcations (1) may be
understood by analysis of a two band model for bulk graphene in zero magnetic
field, in which zero energy states exist along the Dirac spiral in reciprocal
space near each K point [29, 30]. With an approximately linear dispersion near
zero energy, an isoenergetic surface forms a tube-like structure around the
spiral, Fig. 2(c). Semiclassical quantization of LL [45, 46, 47, 48] involves
quantization of the area of an isoenergetic surface in a plane perpendicular
to the applied field. For field with component $B_{x}$, we consider cuts of
the tube-like structure in the $y$-$z$ plane, where we find either one
contour, Fig. 2(d), or two contours, Fig. 2(f), depending on the $q_{x}$
value. We calculate the Berry phase [49] of these contours to be either zero
or $\pi$ (modulo $2\pi$) for one or two contours, respectively. Thus, the
region with one contour corresponds to a ‘conventional’ LL spectrum with
constant level spacing $\sim l$ (where $l$ is the level index) whereas the
region with two contours corresponds to a Dirac-like LL spectrum with doubly-
degenerate levels including zero energy and level spacing $\sim\sqrt{l}$.
Analysis of where the zero field spectrum transitions from one to two contours
leads to Eq. (1), as we discuss in further detail in Section II.2.

In the presence of an applied field, the spectrum retains two flat bands near
zero energy. Just as these bands at zero field may be related to edge states
of the SSH model [24, 18, 17, 25, 26, 27], then their presence in a finite
field may be understood by analogy to zero-energy states localized on solitons
in the SSH model as in the Jackiw-Rebbi mechanism [50, 51, 52, 26, 53], as
explained in Section IV. As the magnetic field strength increases, the
position of the soliton tends to move towards the center of the sample until,
for very large magnetic field strengths, the solitons are annihilated and the
energy levels separate away from zero energy. Similar bound states can also be
found in other systems, such as SSH models undergoing dynamic quenching [54]
or periodic driving [55]. Additionally, photonic SSH lattices have been shown
to host a family of stable solitons once nonlinearity is introduced to the
lattice, despite the breakdown of the bulk-boundary correspondence [56, 57,
58, 59].

Section II describes the methodology beginning with the minimal model in which
only nearest-neighbour hopping parameters are considered, and Section II.2
discusses the qualitative interpretation of level bifurcations in terms of the
semiclassical quantization of the zero field Fermi surface. Beyond the minimal
model, Section II.3, we find that the introduction of next-nearest neighbour
and next-nearest layer hopping parameters does not materially affect Eq. (1).
An alternative calculation is presented in Sec. III whereby we numerically
determine the energy spectrum for bulk rhombohedral graphite (i.e., an
infinite number of layers $N$) using a magnetic supercell. These calculations
recover the results presented in Fig. 2(b) for a stack with a finite number of
layers $N$. Section IV explores the topological nature of the band structure
through a comparison to the SSH model, relating the zero-energy flat bands in
a finite field to zero-energy states localized on solitons [50, 51, 52, 26,
53]. We also consider the influence of disorder on the spectrum by calculating
the disorder-averaged density of states (DOS) per unit energy, and we find
that chirality-preserving disorder does not affect the DOS in a magnetic field
to a significant degree, whereas the zero-energy flat band is shown to be
susceptible to chirality-destroying disorder [52, 21, 53].

##  II Methodology

###  II.1 Finite rhombohedral graphene layers with an in-plane magnetic field

For numerical determination of the energy spectrum, we employ a tight-binding
model [15, 19], initially using the minimal model including nearest-neighbour
intra- and inter-layer hopping only. We assume translational invariance within
each layer in order to Fourier transform to $\mathbf{k}_{||}=(k_{x},k_{y})$
space and create a $2N\times 2N$ Hamiltonian for $N$ layers and two atomic
sites (A and B) per layer. We modify the in-plane terms $H_{AB}$ of the
Hamiltonian with a Peierls substitution, which introduces a path integral of
the vector potential such that

$\displaystyle
H_{AB}(\mathbf{k}_{||})\\!=\\!-\gamma_{0}\sum^{3}_{j=1}\exp\\!\Bigg{(}\\!i\mathbf{k}_{||}\cdot\big{(}\mathbf{R}_{Bj}-\mathbf{R}_{A}\big{)}\\!-\\!\frac{ie}{\hbar}\int^{\mathbf{R_{A}}}_{\mathbf{R}_{Bj}}\mathbf{A}\cdot\mathbf{dl}\\!\Bigg{)},$

where $\gamma_{0}$ is the tight-binding parameter for intralayer hopping,
$\mathbf{R}_{A}$ is the position of an A atom, and $\mathbf{R}_{Bj}$ denotes
the positions of three adjacent B atoms. The magnetic vector potential
$\mathbf{A}=z(B_{y},-B_{x},0)$ is chosen so that translational invariance in
the in-plane direction is preserved [39], and $z$ is the coordinate
perpendicular to the graphene layers. Expanding and simplifying the above
yields

$\displaystyle
H_{AB}(\mathbf{k}_{||})=-\gamma_{0}\Bigg{[}\\!\exp\\!\Bigg{(}\\!\frac{ik_{y}a}{\sqrt{3}}-\frac{iezaB_{x}}{2\hbar}\\!\Bigg{)}\\!\hskip
68.28644pt$
$\displaystyle\quad\qquad+\,2\exp\\!\Bigg{(}\\!\frac{-ik_{y}a}{2\sqrt{3}}+\frac{iezaB_{x}}{4\sqrt{3}\hbar}\\!\Bigg{)}\\!\cos\\!\Bigg{(}\\!\frac{k_{x}a}{2}+\frac{ezaB_{y}}{4\hbar}\\!\Bigg{)}\\!\Bigg{]}\\!.$

Figure 3: Low-energy band structures for rhombohedral multilayer graphene with
an applied in-plane magnetic field, obtained through the numerical
diagonalization of Hamiltonian (2). (a) $N=300$ layers and $B_{x}=5$T, where,
for clearer presentation, only every $5$th pair of bands are depicted. (b)
$N=1000$ layers with $B_{x}=1$T, where only every $25$th pair of bands are
depicted. Bifurcations are still present at locations described by relation
(1), but are less visible when compared to Fig. 2(b) due to the larger layer
number.

We shift the in-plane wave vector $\mathbf{k}_{||}$ as
$\mathbf{q}=\hbar(\mathbf{k}_{||}-\mathbf{K_{\xi}})$ to be measured with
respect to the K-point at $\mathbf{K_{\xi}}=\xi(4\pi/3a,0)$ where $\xi=\pm 1$
is a valley index. Then, we expand the Hamiltonian for small $|\mathbf{q}|$
and small $|\mathbf{B}|$, which is valid for $|\mathbf{q}|a/\hbar\ll 1$ and
$e|z|a|\mathbf{B}|/\hbar\ll 1$, conditions which hold at low energy. The
resulting Hamiltonian reads

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!H^{N}(\mathbf{q})\\!=\\!\begin{pmatrix}0&\hbar
v\pi_{1}^{\dagger}&0&0&\cdots&0&0\\\ \hbar
v\pi_{1}&0&\gamma_{1}&0&\cdots&0&0\\\ 0&\gamma_{1}&0&\hbar
v\pi_{2}^{\dagger}&\cdots&0&0\\\ 0&0&\hbar v\pi_{2}&0&\cdots&0&0\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\ 0&0&0&0&\cdots&0&\hbar
v\pi_{N}^{\dagger}\\\ 0&0&0&0&\cdots&\hbar v\pi_{N}&0\end{pmatrix}\\!\\!\\!,$
(2)

where $v=\sqrt{3}\gamma_{0}a/2\hbar$ and

$\displaystyle\pi_{n}(\mathbf{q})=\xi(q_{x}+ez_{n}B_{y}/\hbar)+i(q_{y}-ez_{n}B_{x}/\hbar)\,,$

where $z_{n}$ is the position of the $n^{\mathrm{th}}$ layer on the $z$ axis,
Fig. 1. We choose $z=0$ to be in the middle of the layered system. In polar
coordinates we can write

$\displaystyle\pi_{n}(\mathbf{q})\\!=\\!\\!\sqrt{\\!(\\!q_{x}\\!+\\!ez_{n}B_{y}/\\!\hbar)^{2}\\!\\!+\\!\\!(q_{y}\\!\\!-\\!ez_{n}B_{x}/\\!\hbar)^{2}}\,e^{i\xi\varphi_{n}}\\!\\!=\\!r_{n}e^{i\xi\varphi_{n}},$

and it is clear the phase $\varphi_{n}$ can be gauged away (e.g., by a
redefinition of the atomic orbital wave functions). However, the magnetic
field retains its influence through the magnitude $r_{n}$, and can not be
gauged away. Without loss of generality, we set $B_{y}=0$ such that the
magnetic field $\mathbf{B}=\left(B_{x},0,0\right)$. By diagonalizing the
Hamiltonian (2) we obtain the band structure of Fig. 2(b), which, when
compared to the zero field band structure of Fig. 2(a), shows a spreading of
the energy levels along the energy axis and the appearance of bifurcations of
levels near the K-point, from doubly degenerate (8-fold including spin and
valley degeneracy) to singly degenerate states (4-fold including spin and
valley degeneracy). We see similar qualitative features for larger layer
numbers and smaller fields, Fig 3, although the bifurcations are less visible
for these parameter values. For this figure and other numerics in this paper
the following hopping parameter values have been used: $\gamma_{0}=3.16eV$,
$\gamma_{1}=0.381eV$, $\gamma_{2}=-0.017eV$, $\gamma_{3}=0.38eV$,
$\gamma_{4}=0.14eV$ [60, 61], $a=2.46$Å for the in-plane lattice constant, and
$d=3.46$Å for the interlayer spacing [62].

The qualitative features of Fig. 2(b) are still observable even at
comparatively low field strengths [44], although, at $q_{x}=0$, the energy
levels are only slightly spread from the zero-field degeneracy point of bulk
eigenvalues at $\varepsilon=\pm\gamma_{1}$. At low field, the bifurcations are
also localized in this region of the $\varepsilon$-$q_{x}$ plot, but still
obey relationship (1). This also holds true for small layer numbers $N\leq 10$
and larger fields, although the zero energy states do not extend to the band
gap closing point $q_{x}=q_{c}$. Additionally, bifurcations are not present
for $N\leq 4$ due to there being an insufficient number of bands to realise
doubly degenerate states for $\varepsilon<\gamma_{1}$. The bifurcations occur
at positions (1) that can be understood through a comparison to an infinite
layer graphene model, as described in the next subsection.

###  II.2 Semiclassical quantization of the zero-field Fermi surface

Figure 4: Dirac spirals of Eq. (4) in ${\bf k}$-space located at the bulk
rhombohedral graphite valley corresponding to $\xi=-1$. (a) shows the zero
energy Dirac points forming a nodal line, and (b) the points of an arbitrary
small energy $\varepsilon<\gamma_{1}$ centerd around zero energy, forming a
‘squashed’ tube-like structure.

To develop a qualitative understanding of the numerically-obtained spectrum,
we consider bulk rhombohedral graphite for zero magnetic field. Assuming
translational invariance in both the intra- and inter-layer directions, we
introduce perpendicular wave vector $k_{z}$ in addition to the small in-plane
wave vector $\mathbf{q}$ measured from the K-point. By doing so the
Hamiltonian [63, 29, 30] is reduced to a $2\times 2$ matrix,

$\displaystyle H^{\mathrm{inf}}(\mathbf{q},k_{z})=\begin{pmatrix}0&\hbar
v\pi^{\dagger}-\gamma_{1}e^{ik_{z}d}\\\ \hbar
v\pi-\gamma_{1}e^{-ik_{z}d}&0\end{pmatrix},$ (3)

where $d$ is the interlayer spacing and $\pi=\xi q_{x}+iq_{y}$. In-plane wave
vectors at which the energy is zero are represented in polar coordinates as
$\mathbf{q}_{D}=(q_{D}\cos\varphi_{D},q_{D}\sin\varphi_{D})$ where
$q_{D}=|\mathbf{q}_{D}|$ with

$\displaystyle q_{D}=q_{c}=\frac{\gamma_{1}}{v\hbar}\,,$
$\displaystyle\varphi_{D}=-\xi\left(k_{z}d-\frac{\pi}{2}\right)-\frac{\pi}{2}\,.$
(4)

The dependence of $\varphi_{D}$ on $k_{z}$ shows that these solutions
represent a spiral through ${\bf k}$-space, shown in Fig. 4(a). The energy
dispersion very close to the Dirac spiral is linear as $\varepsilon=\hbar
v|\mathbf{\kappa}|$, where $\mathbf{\kappa}=(\kappa_{x},\kappa_{y})$ is the
wave vector measured from the spiral for a given $k_{z}$ value. Isoenergetic
contours are circles of constant radius $\mathbf{\kappa}$ centerd around the
spiral for each $k_{z}$ value, producing the tubular structure shown in Fig.
4(b).

Onsager semiclassical quantization of LL [45, 46, 47, 48] involves
quantization of the area $S(\varepsilon)$ enclosed by isoenergetic contours in
${\bf k}$-space in a plane perpendicular to the magnetic field, and is
dependent on the Berry phase $\gamma$ of these contours,

$\displaystyle
S(\varepsilon)=\frac{2\pi}{\lambda_{B}^{2}}\left(l+\frac{1}{2}-\frac{\gamma}{2\pi}\right),$
(5)

where $\lambda_{B}=\sqrt{\hbar/eB_{x}}$ is the magnetic length, and integer
$l\geq 0$. For example, a ‘conventional’ 2D semiconductor has quadratic
dispersion $\varepsilon=\hbar^{2}k^{2}/(2m_{*})$, effective mass $m_{*}$, and
zero Berry phase. With area $S(\varepsilon)=\pi k^{2}=2\pi
m_{*}\varepsilon/\hbar$, Eq. (5) gives the harmonic oscillator relation
$\varepsilon_{l}=\hbar\omega_{c}(l+1/2)$ with cyclotron frequency
$\omega_{c}=eB/m_{*}$. By way of contrast, monolayer graphene has a linear
dispersion relation $\varepsilon=\pm\hbar vk$, Berry phase $\pi$, and area
$S(\varepsilon)=\pi k^{2}=\pi\varepsilon^{2}/(\hbar v)^{2}$. The Berry phase
cancels the factor of $1/2$ in Eq. (5), producing
$\varepsilon_{l,\pm}=\pm(\hbar v/\lambda_{B})\sqrt{2l}$, which does not have
zero point energy $\hbar\omega_{c}/2$, but a level fixed at
$\varepsilon_{0}=0$. We assume that $l=0$ is admissible in Eq. (5), although
it may relate to contours enclosing zero area, and is therefore somewhat ill-
defined [46, 47].

In this paper, we apply the semiclassical quantization Eq. (5) to infinite
layer graphite, Eq. (3), in an in-plane field in the $x$ direction.
Perpendicular fields have already been studied, where the relevant contours
are circles of radius $|\kappa|$ in the $x$-$y$ plane, resulting in a LL
spectrum that it essentially the same as monolayer graphene [63]. For an in-
plane field in the $x$ direction, the contours instead lie in the $y$-$z$
plane, Fig. 2. The number of contours varies for different values of $q_{x}$,
and we find

$\displaystyle\begin{cases}\text{0 contours}&\text{if
}|q_{x}|>q_{c}+|\varepsilon|/(\hbar v)\\\ \text{1 contour}&\text{if
}q_{c}-|\varepsilon|/(\hbar v)<|q_{x}|<q_{c}+|\varepsilon|/(\hbar v)\\\
\text{2 contours}&\text{if }|q_{x}|<q_{c}-|\varepsilon|/(\hbar v)\\\
\end{cases}.$

Moreover, when there are two contours, they have the same area. Thus, one and
two contours correspond to singly and doubly degenerate energy levels,
respectively, i.e., the same degeneracies observed in the energy spectrum
under an applied in-plane magnetic field. The location at which the level
bifurcations occur in the finite field $\varepsilon$-$q_{x}$ plot are the same
as the points of transition between one and two contours in the zero field
quantization.

Figure 5: Low-energy band structures for rhombohedral graphene with an applied
in-plane magnetic field, with an overlay (red) of the linear relationship for
the location of the bifurcation points in the LL spectrum, Eq. (1). (a) shows
the band structure obtained from the semiclassical quantization, Eq. (5), of
the numerically obtained areas enclosed by isoenergetic contours in the
infinite bulk graphite model, Fig. 2. (b) depicts a similar band structure for
a finite 40 layer rhombohedral graphene stack with an applied in-plane
magnetic field of $B_{x}=50$T, obtained from the numerical diagonalization of
Hamiltonian (2).

We find that a path around a single contour acquires a Berry phase of zero
whereas the double contours each correspond to Berry phase $\pi$ (modulo
$2\pi$). The latter case corresponds to a Dirac-like LL spectrum as observed
in the finite field band structure at $q_{x}=0$, Fig. 2(b). At zero field and
$q_{x}=0$, each of the double contours has a dispersion at low energy given by

$\displaystyle\varepsilon\approx\pm\hbar
v\sqrt{\kappa_{y}^{2}+(q_{c}d\kappa_{z})^{2}}$ (6)

where $\kappa_{z}$ is the wave vector in the $z$ direction measured from the
center of the Dirac spiral. Contours described by Eq. (6) are approximately
circular with radius $\varepsilon$ and the corresponding enclosed area is
$S(\varepsilon)=\pi\varepsilon^{2}$. Substituting this into Eq. (5) returns
the Dirac-like LL spacing $\sim\sqrt{l}$, and, as there are two identical
contours, these levels are doubly degenerate. We can extend this to finite
values of $q_{x}$ by numerically calculating the contour areas, resulting in
LL corresponding to an infinite system, Fig. 5(a), while maintaining the same
qualitative features as the band structure obtained numerically from the model
of a system with a finite number of layers, Fig. 5(b). The location of
bifurcations of contours in the zero field spectrum can be written as Eq. (1),
and superimposing this linear $\varepsilon$-$q_{x}$ relationship over the
numerically obtained band structures, Fig. 5, visually confirms it.

###  II.3 Inclusion of long-range hopping parameters

We consider whether the qualitative features discussed in the context of the
minimal model still hold when considering up to next-nearest layer (long-
range) hopping parameters in the Hamiltonian. The full Hamiltonian including
these extra parameters $\gamma_{2}$, $\gamma_{3}$, $\gamma_{4}$ [11, 14, 15,
19] is

$\displaystyle
H^{N}(\mathbf{q})=\begin{pmatrix}\mathcal{A}_{1}&\mathcal{B}&\mathcal{C}&\cdots&0&0&0\\\
\mathcal{B}^{\dagger}&\mathcal{A}_{2}&\mathcal{B}&\cdots&0&0&0\\\
\mathcal{C}^{\dagger}&\mathcal{B}^{\dagger}&\mathcal{A}_{3}&\cdots&0&0&0\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
0&0&0&\cdots&\mathcal{A}_{N-2}&\mathcal{B}&\mathcal{C}\\\
0&0&0&\cdots&\mathcal{B}^{\dagger}&\mathcal{A}_{N-1}&\mathcal{B}\\\
0&0&0&\cdots&\mathcal{C}^{\dagger}&\mathcal{B}^{\dagger}&\mathcal{A}_{N}\end{pmatrix},$
(7)

where each element represents a 2x2 matrix,

$\displaystyle\mathcal{A}_{n}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}0&\hbar v\pi^{\dagger}_{n}\\\ \hbar
v\pi_{n}&0\end{pmatrix};\qquad\mathcal{B}=\begin{pmatrix}-\hbar
v_{4}\pi^{\dagger}&\hbar v_{3}\pi\\\ \gamma_{1}&-\hbar
v_{4}\pi^{\dagger}\end{pmatrix},$ $\displaystyle\mathcal{C}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}0&\frac{\gamma_{2}}{2}\\\ 0&0\end{pmatrix};\qquad
v_{i}=\frac{\sqrt{3}a\gamma_{i}}{2\hbar}.$

The band structure for this Hamiltonian, Fig. 6, retains significant
similarities to the minimal model, but $\gamma_{4}$ introduces finite
dispersion with electron-hole asymmetry and $\gamma_{3}$ creates trigonal
warping (anisotropy of the dispersion).

Figure 6: A comparison of band structures for 20 layer rhombohedral graphene
with up to next-nearest layer hopping parameters $\gamma_{0,1,2,3,4}$
included, according to the numerical diagonalization of Hamiltonian (7),
located at the valley corresponding to $\xi=1$. (a) Shows the band structure
for a 20 layer system for $B=0$ to illustrate the difference compared to Fig.
1(a), (b) shows the same system with the overlay of numerically obtained
bifurcation points for $B=100$T. Negative $q_{x}$ values have been included
here to highlight anisotropy of the dispersion (trigonal warping) due to the
presence of parameter $\gamma_{3}$.

For the model (3) of infinite layer graphite at zero field, the additional
parameters modify matrix elements [63] as

$\displaystyle H^{\mathrm{inf}}_{11}$ $\displaystyle=$ $\displaystyle
H^{\mathrm{inf}}_{22}=2\hbar v_{4}|\mathbf{q}|\cos(\varphi+k_{z}d),$
$\displaystyle H^{\mathrm{inf}}_{12}$ $\displaystyle=$
$\displaystyle\left(H^{\mathrm{inf}}_{21}\right)^{*}=-\hbar
v|\mathbf{q}|e^{-i\varphi}+\gamma_{1}e^{ik_{z}d}$
$\displaystyle\qquad\qquad\qquad+\,\gamma_{2}e^{-2ik_{z}d}+\hbar
v_{3}|\mathbf{q}|e^{i(\varphi-k_{z}d)},$

where $\varphi$ is the polar angle. Since zero-energy states coincide with a
band degeneracy, their location may be determined by setting
$H_{12}^{\mathrm{inf}}=0$, and $\gamma_{4}$ has no effect on the location of
the nodal line. The solution may be written as a perturbation to Eq. (4) and,
up to $\mathcal{O}(v_{3}/v)$, it can be written as

$\displaystyle q_{D}$ $\displaystyle=$
$\displaystyle\frac{\gamma_{1}}{v\hbar}\left(1+\frac{v_{3}}{v}\cos(3k_{z}d)\right),$
(8) $\displaystyle\Phi_{D}$ $\displaystyle=$
$\displaystyle-\xi\left(k_{z}d+\frac{v_{3}}{v}\sin(3k_{z}d)\right),$ (9)

such that the original terms still dominate [63], and parameter $\gamma_{2}$
may be neglected. The presence of $\gamma_{3}$ (trigonal warping) slightly
modifies the semiclassical estimate (1) of the position of the bifurcations
which we determine numerically and show as the dashed overlaid line in Fig.
5(b). It has recently been suggested that the sign of $v_{3}$ is negative [61,
64], in which case one should invert the sign of $q_{x}$ in Fig. 6(b).

##  III Numerics for infinite bulk rhombohedral graphite

In Section II, we described methodology for numerically determining the
spectrum of a system with a finite number of layers in an external in-plane
field. Here, instead, we describe how to numerically determine the spectrum of
bulk rhombohedral graphite (with an infinite number of layers) using a finite-
sized magnetic supercell. Specifically, we consider a magnetic supercell
described by Hamiltonian

$\displaystyle\mbox{\normalsize$H_{\mathrm{supercell}}^{M}(\mathbf{q})$
}=\begin{pmatrix}m_{1}&0&0&\cdots&\Gamma_{1}\\\ 0&m_{2}&0&\cdots&0\\\
0&0&m_{3}&\cdots&0\\\ \vdots&\vdots&\vdots&\vdots&\vdots\\\
\Gamma_{1}^{\dagger}&0&0&\cdots&m_{M}\\\ \end{pmatrix}.$ (10)

The magnetic supercell is composed of an integer number $M$ of three-layer
hexagonal unit cells, each described by a $6\times 6$ matrix
$m_{n}(\mathbf{q})$ which is similar to the previously considered Hamiltonian
(2),

$\mbox{\normalsize$m_{n}(\mathbf{q})$ }=\\\
\\!\\!\\!\\!\\!\\!\begin{pmatrix}0&\hbar v\pi_{3n-2}^{\dagger}&0&0&0&0\\\
\hbar v\pi_{3n-2}&0&\gamma_{1}&0&0&0\\\ 0&\gamma_{1}&0&\hbar
v\pi^{\dagger}_{3n-1}&0&0\\\ 0&0&\hbar v\pi_{3n-1}&0&\gamma_{1}&0\\\
0&0&0&\gamma_{1}&0&\hbar v\pi_{3n}^{\dagger}\\\ 0&0&0&0&\hbar v\pi_{3n}&0\\\
\end{pmatrix}\\!.$

In the Hamiltonian (10), $\Gamma_{1}$ represents a $6\times 6$ matrix with
$\gamma_{1}$ as its top right entry and zeros everywhere else, i.e.
$(\Gamma_{1})_{ij}=0$ for all $i$, $j$, except $(\Gamma_{1})_{16}=\gamma_{1}$.
This accounts for periodicity of the supercell.

We can choose the in-plane component of the magnetic field to lie in the $x$
direction without loss of generality (for the minimal model). It can only take
certain discrete values related to the number $M$ of hexagonal unit cells, as
determined by periodicity of the Hamiltonian. In particular, for vector
potential ${\bf A}=z(0,-B_{x},0)$, a matrix element connecting positions
$\mathbf{R}$ and $\mathbf{R^{\prime}}$ has a field-dependent phase factor such
as $\Delta=(B_{x}e/\hbar)\int^{\mathbf{R^{\prime}}}_{\mathbf{R}}z\;dy$. This
should be equal to $2\pi$ for a vertical translation of $3Md$ and a horizontal
translation of $\sqrt{3}a/2$, so that allowed values of the magnetic field are

$\displaystyle B_{x}=\frac{2h}{3\sqrt{3}adeM},$ (11)

for integer $M$. For example, when $M=100$, then $B_{x}=96.6$T, and the
numerical calculation requires the diagonalization of a $600\times 600$ matrix
(10). Despite differences, this method produces similar qualitative results to
those for a finite system, section II, including the location of bifurcations
from doubly to singly degenerate states in the energy spectrum as shown in Fig
7.

Figure 7: Low-energy band structures for infinite bulk rhombohedral graphene
with an applied in-plane magnetic field, obtained through the numerical
diagonalization of the supercell Hamiltonian (10). (a) Shows the band
structure for a supercell of size $M=200$, and corresponding field
$B_{x}=96.6$T. (b) Is the band structure for a supercell of size $M=400$, and
corresponding field $B_{x}=48.3$T.

##  IV Analogy to the SSH model

For zero magnetic field, the spectrum of RMG in the minimal model may be
related to the energy levels of the SSH model [24, 18, 17, 25, 26, 27] by
dimensional reduction [65], i.e., by treating the in-plane wave vector
$\mathbf{q}$ as a parameter. In particular, alternating hoppings $t$ and $w$
of the SSH model are equivalent to $\hbar v|{\bf q}|$ and $\gamma_{1}$,
respectively, of RMG. Here we generalize this analogy in the presence of an
in-plane magnetic field, resulting in a one-dimensional tight-binding model
with spatially-dependent hopping parameters that bears some similarity to the
commensurate off-diagonal Aubry-André-Harper model [66].

Starting from the Hamiltonian (2) for $N$ layers of RMG, we gauge away the
phase of the intralayer hopping term for layer $n$ leaving the magnitude,

$\displaystyle\hbar v|\pi_{n}(q_{x},z_{n})|=\hbar vr_{n}=\\!\sqrt{(\hbar
vq_{x})^{2}+(evz_{n}B_{x})^{2}}\,,$ (12)

with $q_{y}=B_{y}=0$. For an even number of layers $N$, the vertical
coordinate of the $n$th layer is $z_{n}=d(n-(N+1)/2)$ with the origin ($z=0$)
at the center of the system.

Treating $q_{x}$ as a parameter, the Hamiltonian (2) may be written as a one-
dimensional tight-binding model which is a generalisation of the SSH model
with $2N$ orbitals,

$\displaystyle H_{\mathrm{SSH}}=\begin{pmatrix}0&t_{1}&0&0&0&\cdots&0&0\\\
t_{1}&0&w&0&0&\cdots&0&0\\\ 0&w&0&t_{2}&0&\cdots&0&0\\\
0&0&t_{2}&0&w&\cdots&0&0\\\ 0&0&0&w&0&\cdots&0&0\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
0&0&0&0&0&\cdots&0&t_{N}\\\ 0&0&0&0&0&\cdots&t_{N}&0\\\ \end{pmatrix}\,,$ (13)

where

$\displaystyle
t_{n}=\sqrt{t^{2}+\left(\left[n-\tfrac{1}{2}(N+1)\right]\chi\right)^{2}}\,,$
(14)

and $\chi$ is a parameter describing the strength of the modulated hopping.
This has the following equivalence with RMG,

$\displaystyle t$ $\displaystyle\equiv$ $\displaystyle\hbar v|q_{x}|,$ (15)
$\displaystyle w$ $\displaystyle\equiv$ $\displaystyle\gamma_{1},$ (16)
$\displaystyle\chi$ $\displaystyle\equiv$ $\displaystyle evd|B_{x}|.$ (17)

The Hamiltonian (13) has constant intercell hopping $w$ with modulated hopping
$t_{n}$ on every intracell bond. The modulation is inversion symmetric about
the center of the system, increasing in magnitude from
$t_{N/2}=\sqrt{t^{2}+(\chi/2)^{2}}\sim t$ at the center to
$t_{N}\sim\sqrt{t^{2}+(N\chi/2)^{2}}$ at the edge.

Figure 8: Probability density of the wavefunctions close to the K-point,
$q_{x}=q_{y}=0$, for rhombohedral graphene multilayers on carbon atoms at
positions $j=1,2,3,\ldots,2N$. (a) shows edge state localization at zero
energy while (b) is a generic bulk state lying within the valence band. Both
(a) and (b) are for a $N=10$ layer system with $20$ atoms in total and zero
magnetic field. (c) shows the distribution of a soliton zero energy state for
a $N=100$ layer system and an in-plane field $B_{x}=50$T. We use a larger
layer number to better depict the two soliton states, $n^{\ast}$ denotes the
center of the soliton.

We can interpret the effect of $\chi$ as producing a texture of $t$ in the
system, i.e. a soliton and anti-soliton pair placed symmetrically about the
center, each of which supports a state near zero energy. For $t<w$ and
$\chi=0$, the Hamiltonian (2) describes the SSH model with a bulk band gap of
$2(w-t)$ and states near zero energy localized at the edges, Fig. 8(a). By
contrast, Fig. 8(b) shows the spatial distribution of a state lying within the
bulk valence band. For nonzero $\chi$, and large enough $N$, $t_{N}$ at the
edge becomes larger than $w$ so that there are no edge states. However, at
some position $n=n^{\ast}$ there will be $t_{n^{\ast}}\approx w$ separating
regions with $t_{n}>w$ towards the edge and with $t_{n}<w$ towards the center.
Location $n^{\ast}$ is at the center of a soliton and supports a zero energy
state, and there is a corresponding antisoliton with the opposite texture of
$t_{n}$ at the other side of the system which also supports a zero energy
state, Fig. 8(c). Thus, as $\chi$ increases, the soliton-antisoliton pair move
from the edges towards the center of the system, always supporting zero energy
states. Finally, for large enough $\chi$, they will reach the center and
annihilate, when $t_{N/2}=\sqrt{t^{2}+(\chi/2)^{2}}\approx w$, and there will
be no zero energy states.

For fixed $\chi$, the location of the soliton $n^{\ast}$ is dependent on the
value of $t$. From $t_{n^{\ast}}\approx w$ it may be estimated at the right
side of the system ($n^{\ast}>(N+1)/2$), say, to be
$n^{\ast}\approx\mathrm{min}\big{(}N,(N+1)/2+\chi^{-1}\sqrt{w^{2}-t^{2}}\big{)}$.
Using the equivalence to RMG, Eqs. (15-17), we can express the layer position
$z^{\ast}\equiv\big{(}n^{\ast}-(N+1)/2\big{)}d$ where the flat band state is
localized as

$\displaystyle\frac{z^{\ast}}{d}\approx\mathrm{min}\bigg{(}\frac{N-1}{2},\frac{a\sqrt{q_{c}^{2}-|q_{x}|^{2}}}{2\pi(\phi_{d}/\phi_{0})}\bigg{)},$
(18)

for $|q_{x}|<q_{c}$ where $q_{c}=\gamma_{1}/(\hbar v)$, $\phi_{d}=ad|B_{x}|$
is the magnetic flux per unit cell, and $\phi_{0}=h/e$ is the single particle
flux quantum. This expression is for the top half of the thin film
($z^{\ast}>0$), and the position for the lower half is $-z^{\ast}$. Fig 9
shows a plot of $z^{\ast}$ obtained by numerically determining the peak
positions of the wavefunctions, with the estimate Eq. (18) shown for
comparison as the dashed line. It shows that the zero-energy flat bands are
localized on different bulk layers of the system, not just the surfaces, as
determined by the wave vector $q_{x}$. This accounts for the behavior observed
for the flat bands at zero energy in RMG in the presence of an in-plane
magnetic field. As the field strength is increased, the zero energy states
persist, but the extent of the flat bands in $q_{x}$ is reduced because, for
large $q_{x}$ values, the solitons have disappeared. Finally, at very large
field strengths, they separate and move into the bulk (when the solitons for
all $q_{x}$ values have vanished).

Figure 9: Vertical position $z^{\ast}$, as a function of the in-plane wave
vector $q_{x}$, of the localized states with zero-energy flat bands in
rhombohedral multilayer graphene with $N=50$ layers and an applied in-plane
magnetic field $B_{x}=50$T. Solid lines are obtained by numerically
determining the peak positions of the wavefunctions, the dashed line is the
estimate Eq. (18). States at $q_{x}=0$ are localized at high $|z|$ values
towards the surface of the thin film. As $q_{x}$ increases, they move inwards
to the center at which point they begin to annihilate one another,
asymptotically approaching the minimum of $|z|=d/2$.

By equating the two terms on the right side of Eq. (18), we can estimate the
extent of the flat band, in terms of $q_{x}$, that remains on the surface in
the presence of the in-plane field as

$\displaystyle q_{x}a=\sqrt{q_{c}^{2}a^{2}-(\pi\phi_{N}/\phi_{0})^{2}},$ (19)

where the magnetic flux summed over all layers is
$\phi_{N}=(N-1)\phi_{d}=ad(N-1)B_{x}$.

The limit $t=0$ for the SSH model corresponds to $q_{x}=0$ in RMG, in which
case $t_{N}\sim N\chi/2$ at the edge. Thus, the zero energy states for $t=0$
remain at the edge until $t_{N}\approx w$, i.e., $\chi\sim 2w/N$. In terms of
RMG, this means that states at $q_{x}=0$ remain on the surfaces until the
magnetic flux $\phi_{N}$ reaches values

$\displaystyle\frac{\phi_{N}}{\phi_{0}}\sim\frac{2\gamma_{1}}{\sqrt{3}\pi\gamma_{0}},$
(20)

where $\phi_{0}=h/e$ is the single particle flux quantum; numerical values
give $\phi_{N}\sim\phi_{0}/20$. We can alternatively write this estimate as
$\phi_{N}/\phi_{0}\sim q_{c}a/\pi$. Up to this point, we can use the zero
field estimate for the energy of the surface states,
$\varepsilon\approx\pm(\hbar v|{\bf q}|)^{N}/\gamma_{1}^{N-1}$, with the
substitution $|{\bf q}|\rightarrow B_{x}$ to estimate that the surface states
at $q_{x}=0$ behave as $\varepsilon\sim\pm B_{x}^{N}$. This is confirmed
numerically in Fig. 10, which also depicts larger magnetic fields breaking
this relationship as the localized states move from the edges towards the
center. The localized states will eventually be destroyed when
$t_{N/2}=\sqrt{t^{2}+(\chi/2)^{2}}\approx w$ and, for $t=0$, this gives an
estimate $\chi\sim 2w$. In terms of RMG, this means that the zero energy
states remain until

$\displaystyle\frac{\phi_{N}}{\phi_{0}}\sim\frac{2N\gamma_{1}}{\sqrt{3}\pi\gamma_{0}},$
(21)

We can alternatively write this estimate as $\phi_{d}/\phi_{0}\sim q_{c}a/\pi$
where $\phi_{d}=ad|B_{x}|$ is the magnetic flux per unit cell. Clearly, this
is a huge magnetic field, a large factor (of $N\gg 1$) greater than the field
(20) for which the zero energy states move away from the edge.

Figure 10: Energy of the states near zero energy at $q_{x}=0$ in $N=20$ layer
rhombohedral graphene as a function of magnetic field $B_{x}$, represented as
a log-log plot. States localized on surfaces of the RMG film at small magnetic
field are characterized by a linear relation between energy and magnetic field
in the log-log plot, reflecting the analytical solution $\varepsilon\propto
B_{x}^{N}$.

We determine the density of states per unit energy per unit area (DOS) to show
how the flat bands at zero energy contribute to sharp features in the DOS.
Since the Hamiltonian (2) satisfies sublattice chiral symmetry in the same way
as the SSH model, we check that the zero energy states are robust to chiral-
preserving disorder by calculating the disorder-averaged DOS in the presence
of disorder. We calculate the DOS per unit energy numerically by approximating
the Dirac delta function as a Lorentzian with finite width
$\zeta=0.005\gamma_{1}$,

$\displaystyle g(\varepsilon)=\frac{1}{\pi
L^{2}}\sum_{n}\frac{\zeta}{(\varepsilon-\varepsilon_{n})^{2}+\zeta^{2}},$ (22)

where $L^{2}$ is the area of the system.

Figure 11: The DOS as a function of energy for $N=20$ layer rhombohedral
graphene systems calculated numerically using the minimal model Hamiltonian
(2) and DOS definition (22). (a) depicts DOS for pristine RMG at zero field.
(b) pristine RMG with an applied in-plane magnetic field of $B_{x}=50\,$T. (c)
is the disorder-averaged DOS at $B_{x}=50\,$T created by averaging over 20
randomized disorder realizations for chirality-destroying disorder. (d)
disorder-averaged DOS at $B_{x}=50\,$T for 20 realizations of chirality-
preserving disorder.

For pristine RMG at zero magnetic field, the DOS is characterized by a series
of Van Hove singularities at energies corresponding to turning points of the
bulk bands, Fig. 11a, with a prominent central peak at zero energy caused by
the flat bands. The DOS in a magnetic field of $B_{x}=50$T is shown in Fig.
11b, and is characterized by the softening of Van Hove singularities, but the
peak corresponding to flat bands remains as prominent as before.

We introduce chirality-breaking disorder through the addition of onsite
energies randomly chosen from a distribution $[-\delta,\delta]$, where
$\delta=50\mbox{meV}\ll\gamma_{1}$ [21]. The resulting disorder-averaged DOS
in the presence of a magnetic field is shown in Fig. 11c with averaging over
20 disorder realisations. We also consider chirality-preserving disorder in
the form of random additions $\delta$ to the interlayer coupling $\gamma_{1}$
as $\gamma_{1}+\delta$, chosen from a distribution $[-\delta,\delta]$, where
$\delta=50\mbox{meV}\ll\gamma_{1}$. The result of averaging over 20 disorder
realisations in a finite magnetic field is shown in Fig. 11d.

The DOS for both types of disorder are characterized by further softening of
the Van Hove singularities at finite energies. However, chirality-destroying
disorder destroys the zero energy flat band peak, whereas chirality-preserving
disorder has no visible effect on the zero-energy peak. This is expected
because the zero energy states rely on the protection of chiral symmetry [67,
25, 52, 68]. In the plot of Fig. 11c, for chirality-destroying disorder, there
is also electron-hole asymmetry in the DOS, but this is only an artefact of
having a finite number of disorder realisations.

##  V Conclusion

We numerically determine the effect of an in-plane magnetic field on the
electronic spectrum of rhombohedral multilayer graphene (RMG). The spectrum
shows two qualitative features: Firstly, there are a series of bifurcations
from doubly- to singly-degenerate energy levels and, secondly, the zero-energy
flat bands, observed at zero field, persist up to very high field strengths.
Using semiclassical quantization [45, 46, 47, 48] of the zero field Fermi
surface of rhombohedral graphite, we are able to relate the presence of level
bifurcations to bifurcations of zero field isoenergetic contours from one
contour, with zero Berry phase, to two contours each with Berry phase $\pi$.
By writing down a one-dimensional tight-binding model that is analogous to RMG
in an in-plane field, we are able to relate the persistence of the zero-energy
bands in large magnetic fields to a soliton texture supporting zero-energy
states in the Su-Schrieffer-Heeger model. Whereas the zero-energy flat bands
at zero field are localized on the surfaces of the RMG thin film, in a finite
magnetic field they are generally localized on different layers, ranging from
the surface to the bulk as a function of the wave vector.

The electronic properties of RMG films may be experimentally accessed by a
variety of methods including scanning tunneling spectroscopy [1], magneto-
Raman [69, 70, 71, 72, 73, 74, 2], photoemission [1, 3], and magneto-transport
[75, 76, 8, 77]. For example, scanning tunneling microscopy was recently used
to probe the surface states of RMG with up to $N=17$ layers [10], and magneto-
transport properties, probing the Landau level spectra, have been measured in
RMG with up to $N=50$ layers [8] and in Bernal-stacked graphitic films with up
to a few hundred layers [76, 77].

We neglected spin splitting and assumed a fourfold degeneracy due to spins and
valleys throughout. In an external in-plane magnetic field, splin splitting of
the energy levels will occur. With $g$-factor $g=2$, the spin splitting should
be $\Delta\varepsilon=2\mu_{B}B$ for Bohr magneton $\mu_{B}$ and field
strength $B$. Thus, for $B=100\,$T, we expect $\Delta\varepsilon=0.012\,$eV
and $\Delta\varepsilon/\gamma_{1}=0.030$. Experiments [78] in monolayer and
bilayer graphene reported an effective $g$-factor, $g^{\ast}=2.7\pm 0.2$, and
attributed its enhancement to electron-electron interaction effects.
Nevertheless, even with such an enhancement, we expect the spin splitting to
be a small effect on the scale of our typical plots, e.g. Fig. 2(b).

All relevant data presented in this paper can be accessed [79].

###### Acknowledgements.

The authors thank A. Mishchenko and S. Ozdemir for helpful discussions.

## References

  * [1] D. Pierucci, H. Sediri, M. Hajlaoui, J.-C. Girard, T. Brumme, M. Calandra, E. Velez-Fort, G. Patriarche, M. G. Silly, G. Ferro, V Soulière, M. Marangolo, F. Sirotti, F. Mauri, and A. Ouerghi, Evidence for Flat Bands near the Fermi Level in Epitaxial Rhombohedral Multilayer Graphene, ACS Nano 9, 5432 (2015).
  * [2] Y. Henni, H. P. O. Collado, K. Nogajewski, M. R. Molas, G. Usaj, C. A. Balseiro, M. Orlita, M. Potemski, and C. Faugeras, Rhombohedral Multilayer Graphene: A Magneto-Raman Scattering Study, Nano Lett. 16, 3710 (2016).
  * [3] H. Henck, J. Avila, Z. Ben Aziza, D. Pierucci, J. Baima, B. Pamuk, J. Chaste, D. Utt, M. Bartos, K. Nogajewski, B. A. Piot, M. Orlita, M. Potemski, M. Calandra, M. C. Asensio, F. Mauri, C. Faugeras, and A. Ouerghi, Flat electronic bands in long sequences of rhombohedral-stacked graphene, Phys. Rev. B 97, 245421 (2018).
  * [4] T. Latychevskaia, S.-K. Son, Y. Yang, D. Chancellor, M. Brown, S. Ozdemir, I. Madan, G. Berruto, F. Carbone, A. Mishchenko and K. S. Novoselov, Stacking transition in rhombohedral graphite, Front. Phys. 14, 13608 (2019).
  * [5] Y. Yang, Y.-C. Zou, C. R. Woods, Y. Shi, J. Yin, S. Xu, S. Ozdemir, T. Taniguchi, K. Watanabe, A. K. Geim, K. S. Novoselov, S. J. Haigh, and A. Mishchenko, Stacking Order in Graphite Films Controlled by van der Waals Technology, Nano Lett. 19, 8526 (2019).
  * [6] F. R. Geisenhof, F. Winterer, S. Wakolbinger, T. D. Gokus, Y. C. Durmaz, D. Priesack, J. Lenz, F. Keilmann, K. Watanabe, T. Taniguchi, R. Guerrero-Avilés, M. Pelc, A. Ayuela and R. Thomas Weitz, Anisotropic Strain-Induced Soliton Movement Changes Stacking Order and Band Structure of Graphene Multilayers: Implications for Charge Transport, ACS Appl. Nano Mater. 2, 6067 (2019).
  * [7] C. Bouhafs, S. Pezzini, F. R. Geisenhof, N. Mishra, V. Mis̆eikis, Y. Niu, C. Struzzi, R. T. Weitz, A. A. Zakharov, S. Forti, and C. Coletti, Synthesis of large-area rhombohedral few-layer graphene by chemical vapor deposition on copper, Carbon 177, 282 (2021).
  * [8] Y. Shi, S. Xu, Y. Yang, S. Slizovskiy, S. V. Morozov, S.-K. Son, S. Ozdemir, C. Mullan, J. Barrier, J. Yin, A. I. Berdyugin, B. A. Piot, T. Taniguchi, K. Watanabe, V. I. Fal’ko, K. S. Novoselov, A. K. Geim, and A. Mishchenko, Electronic phase separation in multilayer rhombohedral graphite, Nature 584, 210 (2020).
  * [9] A. Kerelsky, C. Rubio-Verdú, L. Xian, D. M. Kennes, D. Halbertal, N. Finney, L. Song, S. Turkel, L. Wanga, K. Watanabe, T. Taniguchi, J. Hone, C. Dean, D. N. Basov, A. Rubio, and A. N. Pasupathy, Moiréless correlations in ABCA graphene, Proc. Natl. Acad. Sci. U.S.A. 118, 2017366118 (2021).
  * [10] I. Hagymási, M. S. Mohd Isa, Z. Tajkov, K. Márity, L. Oroszlány, J. Koltai, A. Alassaf, P. Kun, K. Kandrai, A. Pálinkás, P. Vancsó, L. Tapasztó, and P. Nemes-Incze, Observation of competing, correlated ground states in the flat band of rhombohedral graphite, Sci. Adv. 8, eabo6879 (2022).
  * [11] J. W. McClure, Electron energy band structure and electronic properties of rhombohedral graphite, Carbon 7, 425 (1969).
  * [12] S. Latil and L. Henrard, Charge Carriers in Few-Layer Graphene Films, Phys. Rev. Lett. 97, 036803 (2006).
  * [13] M. Aoki and H. Amawashi, Dependence of band structures on stacking and field in layered graphene, Solid State Commun. 142, 123 (2007).
  * [14] D. P. Arovas and F. Guinea, Stacking faults, bound states, and quantum Hall plateaus in crystalline graphite, Phys. Rev. B 78, 245416 (2008).
  * [15] M. Koshino and E. McCann, Trigonal warping and Berry’s phase $N\pi$ in ABC-stacked multilayer graphene Phys. Rev. B 80, 165409 (2009).
  * [16] F. Zhang, J. Jung, G. A. Fiete, Q. Niu, and A. H. MacDonald, Spontaneous Quantum Hall States in Chirally Stacked Few-Layer Graphene Systems, Phys. Rev. Lett. 106, 156801 (2011).
  * [17] T. T. Heikkilä and G. E. Volovik, Dimensional crossover in topological matter: Evolution of the multiple Dirac point in the layered system to the flat band on the surface, JETP Lett. 93, 59 (2011).
  * [18] R. Xiao, F. Tasnádi, K. Koepernik, J. W. F. Venderbos, M. Richter, and M. Taut, Density functional investigation of rhombohedral stacks of graphene: Topological surface states, nonlinear dielectric response, and bulk limit, Phys. Rev. B 84, 165404 (2011).
  * [19] S. Slizovskiy, E. McCann, M. Koshino, and V. I. Fal’ko, Films of rhombohedral graphite as two-dimensional topological semimetals, Commun. Phys. 2, 164 (2019).
  * [20] A. Garcia-Ruiz, S. Slizovskiy, M. Mucha-Kruczyński, and V. I. Fal’ko, Spectroscopic Signatures of Electronic Excitations in Raman Scattering in Thin Films of Rhombohedral Graphite, Nano Lett. 19, 6152 (2019).
  * [21] J. H. Muten, A. J. Copeland, and E. McCann, Exchange interaction, disorder, and stacking faults in rhombohedral graphene multilayers, Phys. Rev. B 104, 035404 (2021).
  * [22] H. Zhou, T. Xie, A. Ghazaryan, T. Holder, J. R. Ehrets, E. M. Spanton, T. Taniguchi, K. Watanabe, E. Berg, M. Serbyn, and A. F. Young, Half- and quarter-metals in rhombohedral trilayer graphene, Nature 598, 429 (2021).
  * [23] H. Zhou, T. Xie, T. Taniguchi, K. Watanabe, and A. F. Young, Superconductivity in rhombohedral trilayer graphene, Nature 598, 434 (2021).
  * [24] W. P. Su, J. R. Schrieffer, and A. J. Heeger, Solitons in polyacetylene, Phys. Rev. Lett. 42, 1698 (1979).
  * [25] J. K. Asbóth, L. Oroszlány, and A. Pályi, A Short Course on Topological Insulators (Springer, Switzerland, 2016).
  * [26] J. Cayssol and J.-N. Fuchs, Topological and geometrical aspects of band theory, J. Phys. Mater. 4, 034007 (2021).
  * [27] E. McCann, Catalog of noninteracting tight-binding models with two energy bands in one dimension, Phys. Rev. B 107, 245401 (2023).
  * [28] T. T. Heikkilä, N. B. Kopnin, and G. E. Volovik, Flat bands in topological media, JETP Lett. 94 233 (2011).
  * [29] C.-H. Ho, C.-P. Chang, and M.-F. Lin, Landau subband wave functions and chirality manifestation in rhombohedral graphite, Solid State Commun. 197, 11 (2014).
  * [30] C.-H. Ho, C.-P. Chang, and M.-F. Lin, Evolution and dimensional crossover from the bulk subbands in ABC-stacked graphene to a three-dimensional Dirac cone structure in rhombohedral graphite, Phys. Rev. B 93, 075437 (2016).
  * [31] W. Chen, H.-Z. Lu, and O. Zilberberg, Weak Localization and Antilocalization in Nodal-Line Semimetals: Dimensionality and Topological Effects, Phys. Rev. Lett. 122, 196603 (2019).
  * [32] 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 (2005).
  * [33] Y. Zhang, Y.-W. Tan, H. L. Stormer, and P. Kim, Experimental observation of the quantum Hall effect and Berry’s phase in graphene, Nature 438, 201 (2005).
  * [34] C. H. Ho, Y. H. Ho, Y. H. Chiu, Y. N. Chen, and M. F. Lin, Magneto-electronic properties of rhombohedral trilayer graphene: Peierls tight-binding model, Ann. Phys. (N. Y.) 326, 721 (2011).
  * [35] C.-H Ho, Y.-H. Ho, Y.-Y. Liao, Y.-H. Chiu, C.-P. Chang, and M.-F. Lin, Diagonalization of Landau Level Spectra in Rhombohedral Graphite, J. Phys. Soc. Jpn 81, 024701 (2012).
  * [36] M. B. Lundeberg and J. A. Folk, Rippled Graphene in an In-Plane Magnetic Field: Effects of a Random Vector Potential, Phys. Rev. Lett. 105, 146804 (2010).
  * [37] M. Van der Donck, F. M. Peeters, and B. Van Duppen, Transport properties of bilayer graphene in a strong in-plane magnetic field, Phys. Rev. B 93, 115423 (2016).
  * [38] W. Qin and A. H. MacDonald, In-Plane Critical Magnetic Fields in Magic-Angle Twisted Trilayer Graphene, Phys. Rev. Lett. 127, 097001 (2021).
  * [39] N. Kheirabadi, E. McCann, and V. I. Fal’ko, Magnetic ratchet effect in bilayer graphene, Phys. Rev. B 94, 165404 (2016).
  * [40] N. Kheirabadi, E. McCann, and V. I. Fal’ko, Cyclotron resonance of the magnetic ratchet effect and second harmonic generation in bilayer graphene, Phys. Rev. B 97, 075415 (2018).
  * [41] A. L. Szabó and B. Roy, Metals, fractional metals, and superconductivity in rhombohedral trilayer graphene, Phys. Rev. B 105, L081407 (2022).
  * [42] S. S. Pershoguba and V. M. Yakovenko, Energy spectrum of graphene multilayers in a parallel magnetic field, Phys. Rev. B 82, 205408 (2010).
  * [43] N. A Goncharuk and L. Smrčka, Tight-binding description of Landau levels of graphite in tilted magnetic fields, J. Phys.: Condens. Matter 24, 185503 (2012).
  * [44] For clarity of presentation, we show numerical data for systems with a relatively small number of layers and very large magnetic fields, see Fig. 1(b). However, similar qualitative features appear for lower magnetic field strengths in samples with more layers, Fig 3, although bifurcations in the energy spectrum are less visible for these parameter values.
  * [45] L. Onsager, Interpretation of the de Haas-van Alphen effect, Philos. Mag. 43, 1006 (1952).
  * [46] A. Kormányos, P. Rakyta, L. Oroszlány, and J. Cserti, Bound states in inhomogeneous magnetic field in graphene: Semiclassical approach, Phys. Rev. B 78, 045430 (2008).
  * [47] J. N. Fuchs, F. Piéchon, M. O. Goerbig, and G. Montambaux, Topological Berry phase and semiclassical quantization of cyclotron orbits for two dimensional electrons in coupled band models, Eur. Phys. J. B 77, 351 (2010).
  * [48] D. Xiao, M.-C. Chang, and Q. Niu, Berry phase effects on electronic properties, Rev. Mod. Phys. 82, 1959 (2010).
  * [49] M. V. Berry, Quantal phase factors accompanying adiabatic changes, Proc. R. Soc. Lond. A 392, 45 (1984).
  * [50] R. Jackiw and C. Rebbi, Solitons with fermion number $1/2$, Phys. Rev. D 13, 3398 (1976).
  * [51] R. Jackiw, Fractional and Majorana fermions: the physics of zero-energy modes, Phys. Scr. 2012, 014005 (2012).
  * [52] M. Scollon and M. P. Kennett, Persistence of chirality in the Su-Schrieffer-Heeger model in the presence of on-site disorder, Phys. Rev. B 101, 144204 (2020).
  * [53] R. E. J. Allen, H. V. Gibbons, A. M. Sherlock, H. R. M. Stanfield, and E. McCann, Nonsymmorphic chiral symmetry and solitons in the Rice-Mele model, Phys. Rev. B 106, 165409 (2022).
  * [54] L. Rossi, F. Rossi, F. Dolcini, Real-space effects of a quench in the Su–Schrieffer–Heeger model and elusive dynamical appearance of the topological edge states, New J. Phys. 24, 013011 (2022).
  * [55] Z. Fedorova (Cherpakova), C. Jörg, C. Dauer, F. Letscher, M. Fleischhauer, S. Eggert, S. Linden, and G. von Freymann, Limits of topological protection under local periodic driving. Light Sci. Appl. 8, 63 (2019).
  * [56] D. Leykam and Y. D. Chong, Edge Solitons in Nonlinear-Photonic Topological Insulators, Phys. Rev. Lett. 117, 143901 (2016).
  * [57] M. A. Gorlach, A. P. Slobozhanyuk, Nonlinear topological states in the Su– Schrieffer–Heeger model, Nanosyst.: Phys. Chem. Math. 8, 695 (2017).
  * [58] M. Guo, S. Xia, N. Wang, D. Song, Z. Chen, and J. Yang, Weakly nonlinear topological gap solitons in Su–Schrieffer–Heeger photonic lattices, Opt. Lett. 45, 6466 (2020).
  * [59] Y.-P. Ma and H. Susanto, Topological edge solitons and their stability in a nonlinear Su-Schrieffer-Heeger model, Phys. Rev. E 104, 054206 (2021).
  * [60] A. B. Kuzmenko, I. Crassee, D. van der Marel, P. Blake, and K. S. Novoselov, Determination of the gate-tunable band gap and tight-binding parameters in bilayer graphene using infrared spectroscopy, Phys. Rev. B 80, 165406 (2009).
  * [61] C. Bao, W. Yao, E. Wang, C. Chen, J. Avila, M. C. Asensio, and S. Zhou, Stacking-Dependent Electronic Structure of Trilayer Graphene Resolved by Nanospot Angle-Resolved Photoemission Spectroscopy, Nano Lett. 17, 1564 (2017).
  * [62] J. D. Bernal, The Structure of Graphite, Proc. R. Soc. Lond. A 106, 749 (1924).
  * [63] C. H. Ho, C. P. Chang, W. P. Su, and M. F. Lin, Precessing anisotropic Dirac cone and Landau subbands along a nodal spiral, New J. Phys. 15, 053032 (2013).
  * [64] A. Garcia-Ruiz, S. Slizovskiy, and V. I. Fal’ko, Flat Bands for Electrons in Rhombohedral Graphene Multilayers with a Twin Boundary, Adv. Mater. Interfaces 10, 2202221 (2023).
  * [65] S. Ryu, A. P. Schnyder, A. Furusaki, and A. W. W. Ludwig, Topological insulators and superconductors: tenfold way and dimensional hierarchy, New J. Phys. 12, 065010 (2010).
  * [66] S. Ganeshan, K. Sun, and S. Das Sarma, Topological Zero-Energy Modes in Gapless Commensurate Aubry-André-Harper Models, Phys. Rev. Lett. 110, 180403 (2013).
  * [67] M. Inui, S. A. Trugman, and E. Abrahams, Unusual properties of midband states in systems with off-diagonal disorder, Phys. Rev. B 49, 3190 (1994).
  * [68] M. Bellec, C. Poli, U. Kuhl, F. Mortessagne, and H. Schomerus, Observation of supersymmetric pseudo-Landau levels in strained microwave graphene, Light Sci. Appl. 9, 146 (2020).
  * [69] C. Faugeras, M. Amado, P. Kossacki, M. Orlita, M. Sprinkle, C. Berger, W. A. de Heer, and M. Potemski, Tuning the Electron-Phonon Coupling in Multilayer Graphene with Magnetic Fields, Phys. Rev. Lett. 103, 186803 (2009).
  * [70] J. Yan, S. Goler, T. D. Rhone, M. Han, R. He, P. Kim, V. Pellegrini, and A. Pinczuk, Observation of Magnetophonon Resonance of Dirac Fermions in Graphite, Phys. Rev. Lett. 105, 227401 (2010).
  * [71] C. Faugeras, M. Amado, P. Kossacki, M. Orlita, M. Kühne, A. A. L. Nicolet, Yu. I. Latyshev, and M. Potemski, Magneto-Raman Scattering of Graphene on Graphite: Electronic and Phonon Excitations, Phys. Rev. Lett. 107, 036807 (2011).
  * [72] M. Kühne, C. Faugeras, P. Kossacki, A. A. L. Nicolet, M. Orlita, Yu. I. Latyshev, and M. Potemski, Polarization-resolved magneto-Raman scattering of graphenelike domains on natural graphite, Phys. Rev. B 85, 195406 (2012).
  * [73] S. Goler, J. Yan, V. Pellegrini, A. Pinczuk, Raman spectroscopy of magneto-phonon resonances in graphene and graphite, Solid State Commun. 152, 1289 (2012).
  * [74] C. Qiu, X. Shen, B. Cao, C. Cong, R. Saito, J. Yu, M. S. Dresselhaus, and T. Yu, Strong magnetophonon resonance induced triple G-mode splitting in graphene on graphite probed by micromagneto Raman spectroscopy, Phys. Rev. B 88, 165407 (2013).
  * [75] F. Chiappini, S. Wiedmann, M. Titov, A. K. Geim, R. V. Gorbachev, E. Khestanova, A. Mishchenko, K. S. Novoselov, J. C. Maan, and U. Zeitler, Magnetotransport in single-layer graphene in a large parallel magnetic field, Phys. Rev. B 94, 085302 (2016).
  * [76] J. Yin, S. Slizovskiy, Y. Cao, S. Hu, Y. Yang, I. Lobanova, B. A. Piot, S.-K. Son, S. Ozdemir, T. Taniguchi, K. Watanabe, K. S. Novoselov, F. Guinea, A. K. Geim, V. Fal’ko, and A. Mishchenko, Dimensional reduction, quantum Hall effect and layer parity in graphite films, Nat. Phys. 15, 437 (2019).
  * [77] C. Mullan, S. Slizovskiy, J. Yin, Z. Wang, Q. Yang, S. Xu, Y. Yang, B. A. Piot, S. Hu, T. Taniguchi, K. Watanabe, K. S. Novoselov, A. K. Geim, V. I. Fal’ko, and A. Mishchenko, Mixing of moiré-surface and bulk states in graphite, Nature 620, 756 (2023).
  * [78] E. V. Kurganova, H. J. van Elferen, A. McCollam, L. A. Ponomarenko, K. S. Novoselov, A. Veligura, B. J. van Wees, J. C. Maan, and U. Zeitler, Spin splitting in graphene studied by means of tilted magnetic-field experiments, Phys. Rev. B 84, 121407(R) (2011).
  * [79] https://doi.org/10.17635/lancaster/researchdata/614