# Dirac Fermion Cloning, Moiré Flat Bands and Magic Lattice Constants in
Epitaxial Monolayer Graphene

Qiangsheng Lu  Department of Physics and Astronomy, University of Missouri,
Columbia, Missouri 65211, USA  Ching-Kai Chiu  Kavli Institute for Theoretical
Sciences, University of Chinese Academy of Sciences, Beijing 100190, China
Congcong Le  Max Planck Institute for Chemical Physics of Solids, 01187
Dresden, Germany  Jacob Cook  Department of Physics and Astronomy, University
of Missouri, Columbia, Missouri 65211, USA  Xiaoqian Zhang  Department of
Physics and Astronomy, University of Missouri, Columbia, Missouri 65211, USA
Xiaoqing He  Electron Microscopy Core Facility, University of Missouri,
Columbia, Missouri 65211, USA  Department of Mechanical and Aerospace
Engineering, University of Missouri, Columbia, MO 65211, USA  Mohammad Zarenia
Department of Physics and Astronomy, University of Missouri, Columbia,
Missouri 65211, USA  Mitchel Vaninger  Department of Physics and Astronomy,
University of Missouri, Columbia, Missouri 65211, USA  Paul F. Miceli
Department of Physics and Astronomy, University of Missouri, Columbia,
Missouri 65211, USA  Chang Liu  Department of Physics, Southern University of
Science and Technology, Shenzhen 518055, China  Tai-Chang Chiang  Department
of Physics, University of Illinois at Urbana-Champaign, 1110 West Green
Street, Urbana, IL 61801-3080, USA  Frederick Seitz Materials Research
Laboratory, University of Illinois at Urbana-Champaign, 104 South Goodwin
Avenue, Urbana, IL 61801-2902, USA  Giovanni Vignale  Department of Physics
and Astronomy, University of Missouri, Columbia, Missouri 65211, USA  Guang
Bian  Department of Physics and Astronomy, University of Missouri, Columbia,
Missouri 65211, USA

###### Abstract

Tuning interactions between Dirac states in graphene has attracted enormous
interest because it can modify the electronic spectrum of the two-dimensional
material, enhance electron correlations, and give rise to novel condensed-
matter phases such as superconductors, Mott insulators, Wigner crystals and
quantum anomalous Hall insulators. Previous works predominantly focus on the
flat band dispersion of coupled Dirac states from different graphene layers.
In this work, we propose a new route to realizing flat band physics in
monolayer graphene under a periodic modulation from substrates. We take
gaphene/SiC heterostructure as a role model and demonstrate experimentally the
substrate modulation leads to Dirac fermion cloning and consequently, the
proximity of the two Dirac cones of monolayer graphene in momentum space. Our
theoretical modeling captures the cloning mechanism of Dirac states and
indicates that flat bands can emerge at certain magic lattice constants of
substrate when the period of modulation becomes nearly commensurate with the
$\bm{(\sqrt{3}\times\sqrt{3})R30^{\circ}}$ supercell of graphene. The results
show that the epitaxial monolayer graphene is a promising platform for
exploring exotic many-body quantum phases arising from interactions between
Dirac electrons.

##  I Introduction

The discovery of graphene has revolutionized modern condensed matter physics
as it provides direct access to the physics of Dirac fermions in solid-state
systemsCastroNeto2009 ; Novoselov2005 ; Zhang2005 . It also sheds light on the
path towards a vast field of novel 2D materials including van der Waals (vdW)
materials and topological materials such as quantum spin Hall insulators
CastroNeto2009 ; Geim2013 ; Katsnelson2006 ; RevModPhys.86.959 ;
PhysRevLett.61.2015 ; RevModPhys.82.3045 ; RevModPhys.83.1057 ; Liu2011 . A
single layer graphene possesses two copies of Dirac cones residing at the
opposite corners of the Brillouin zone, leaving them essentially isolated from
each other. Stacking graphene layers in the default order duplicates Dirac
cones in the same valley and thus makes the two valleys remain decoupled.
Recently, new excitements in graphene-like systems have arisen as a
consequence of the creation of strongly coupled Dirac states in artificially
engineered structures such as twisted bilayer graphene (TBG). The interaction
between Dirac states from the two layers can be effectively tuned by the angle
mismatch and thus leads to emergent collective behaviors of electrons
including Mott insulating states, unconventional superconductivity, emergent
ferromagnetism, quantum anomalous Hall effects Bistritzer12233 ;
Ponomarenko2013 ; Dean2013 ; Cao20181 ; Cao20182 ; Wang2016 ; Sharpe605 ;
Serlin900 ; Kaxiras2020 . The essential ingredient for those new emergent
states is the nearly dispersionless bands at zero energy in the Moiré
Brillouin zone. Achieving the required large periodicity of Moiré pattern in
real space and the closeness of Dirac cones in momentum space generally
requires delicate controls over the twist between two mechanically exfoliated
graphene layers, which places stringent constraints on the techniques of
sample assembly. Therefore, there is a pressing need for accessing flat band
physics in systems without fine tuning on twist angles.

Here we report an alternative route to enable interactions between Dirac
electrons in a single layer of graphene and realize flat bands with aid from
the supporting substrate potential. We observed the cloning of Dirac bands in
monolayer graphene epitaxially grown on SiC substates by angle-resolved
photoemission spectroscopy (ARPES). The periodic substrate potential brings
closer the Dirac states from the two valleys and thus turns on intervalley
coupling, which is precisely captured by our tight-binding simulations. Our
theory further indicates that the perturbed graphene system with nearly
commensurate epitaxial relations hosts absolutely flat bands due to the fact
that the AB coupling between the two valleys is absent. The lack of AB
coupling yields a chiral symmetry of the low-energy effective Hamiltonian and
consequently makes the epitaxial graphene system a natural realization of the
chiral-symmetric continuum model proposed by Tarnopolsky et
al.PhysRevLett.122.106405 .

##  II Experimental Results

First, we study the substrate effects on the Dirac states of graphene. Our
graphene was grown epitaxially on the Si-face of a 6H-SiC substrate. The
lattice structure and epitaxial relation between the graphene overlayer and
the SiC substrate are plotted in Fig. 1a. The lattice constant of graphene and
SiC(0001) surface is 2.46 Å and 3.07 Å, respectively. The Brillouin zone of
graphene and SiC(0001) is depicted by red and blue lines, respectively, in
Fig. 1b. ${\rm K}_{\textrm{Gr}}$ and ${\rm K^{\prime}}_{\textrm{Gr}}$
represent the location of the two valleys of graphene Dirac states. The
atomic-resolution STM image of the graphene sample is shown in Fig. 1c. A
superhexagonal Moiré pattern with a period
$\lambda=6(1+\delta)a_{\textrm{SiC}}$ is observed due to the incommensurate
modulation of the SiC interface layer Conrad2017 ; PhysRevB.96.035411 . To
examine the structural quality, we performed high-resolution TEM measurements
on our graphene samples. A typical TEM image is shown in Fig. 1d. The sample
consists of four well-ordered graphene layers sitting on the carbonized
surface of the SiC substrate. A gap between the graphene and the SiC surface
is noticeable, indicating a sharp interface between the graphene layers and
the substrate. In our experiment, the thickness of graphene layers can be
precisely controlled down to a monolayer. We will focus on results obtained
from monolayer graphene samples in the following discussion, but the physics
discussed here also applies to thicker graphene films.

The fermi surface of the graphene sample mapped by ARPES is presented in Fig.
2a. The Brillouin zone of graphene is marked by the blue dashed lines. At the
corners of the Brillouin zone, the ${\rm K}_{\textrm{Gr}}$ points, we can see
bright fermi surface contours from the Dirac bands of graphene. Note that our
sample is $n$-typed doped with fermi level above the Dirac point. Therefore,
the fermi surface contours are small circles surrounding ${\rm
K}_{\textrm{Gr}}$, as marked by the black arrows. Beside the Dirac states of
graphene, there are extra circular contours located inside the Brillouin zone
of graphene, which are absent on the fermi surface of freestanding graphene.
These new contours marked by blue, green, yellow, and gray arrows are referred
to as blue, green, yellow, and gray contours (or cones) in the following
discussions. The geometric relations between the Brillouin zones of graphene
and SiC is depicted in Fig. 2b. The emergent contour with lesser photoemission
intensity can be considered as the clones of the graphene Dirac cone generated
by the periodic substrate perturbations. Shifting the ${\rm
K^{\prime}}_{\textrm{Gr}}$ by a reciprocal lattice vector of SiC as shown in
Fig. 2c, we can find the location of the blue and green cones. The yellow and
gray contours can be obtained through two successive shifts of the graphene
states, see Fig. 2d. Taking the substrate interaction as a perturbation, we
can attribute the blue (green) and yellow (gray) cones to the first-order and
second-order perturbation effects, which also explains that the intensity of
yellow (gray) cones is apparently lower than that of blue (green) cones. The
locations of all the first-order and second-order clones are summarized in
Fig. 2e. It is consistent with the experimental observations shown in Fig. 2a.

The cloning of Dirac states is deeply rooted in the periodic modulation
exerted by the substrate on the electrons in graphene. To understand the
mechanism of cloning, we performed a tight-binding simulation in which the
substrate effect is approximated by a periodic potential acting on the the
graphene electrons. The sample is n-type doped, so the fermi level is shifted
in the simulation to match the experimental results. There is a good agreement
between the experimental and theoretical fermi surfaces as shown in Figs. 3a
and 3b. All first-order and second-order clones in the simulation show up at
the locations observed in the experimental results. The spectrum along the
line of ‘cut1’ (marked in Fig. 3a) is plotted in Figs. 3c and 3d. The
brightness of bands indicates the photoemission intensity in the ARPES
spectrum and the spectral weight (the probability of finding the electron with
the corresponding energy and momentum) in the TB simulation. The clones show
the same dispersion as the primary Dirac cones but have lesser spectral weight
compared to the primary cone, consistent with the perturbative picture. This
is further corroborated by the ARPES spectrum taken along the lines of
‘cut2’-‘cut5’. Unlike the primary Dirac cones sitting far apart at the corners
of the Brillouin zone, the clones are much closer to each other in momentum
space. For example, the distance between the yellow and green clones is 0.21
Å-1 while the spacing between neighboring primary cones is 1.7 Å-1. The
smaller distance between clones enables them to overlap in momentum space, see
the iso-energy contours in the supplementary information. The yellow and green
cones cross at $E=-1.3$ eV while the blue and green cones cross at $E=-1.8$
eV. Such crossings of Dirac states are unavailable in freestanding graphene
films. It is worth noting that there is no gap opened at the crossings of the
green and blue bands. This is because both green and blue contours originate
from the same valley (see Fig. 2c) and thus they do not hybridize. Only Dirac
bands from different valleys can interact with other other and open
hybridization gaps, which will be shown in the following discussions. The
cloning of Dirac bands occurs not only in the monolayer graphene but also in
thicker films. Figure 3f shows the spectra of the primary Dirac cone from
samples with different thicknesses, namely, 1 ML, 3 ML, 5 ML. The number of
Dirac bands in the spectra indicates the thickness of the sample Zhou2006 ;
Liu2010 . Despite of the weaker intensity, the clone band (the blue cone
marked in Figs. 2a and 3c) are still observable in the 3-ML and 5-ML samples
as shown in Fig. 3f.

##  III Discussions

The perturbation theory of quantum mechanics can capture the essential physics
of Dirac fermion cloning in the graphene/SiC heterostructure. The SiC
substrate potential can be treated as a small perturbation and exerts a
periodic modulation to the Dirac states of graphene PhysRevB.96.035411 .
Taking $W({\bm{x}})$ as the SiC potential, the eigenfunction up to the second-
order perturbations is given by

$\displaystyle|\Psi_{\bm{k}}\rangle=$
$\displaystyle|\Psi_{\bm{k}}^{0}\rangle+\sum_{{\bm{p}}\neq{\bm{k}}}|\Psi_{\bm{p}}^{0}\rangle\frac{\langle\Psi_{\bm{p}}^{0}|W|\Psi_{\bm{k}}^{0}\rangle}{E_{\bm{k}}^{0}-E_{\bm{p}}^{0}}+\sum_{{\bm{p}}\neq{\bm{k}},{\bm{l}}\neq{\bm{k}}}|\Psi_{\bm{p}}^{0}\rangle\frac{\langle\Psi_{\bm{p}}^{0}|W|\Psi_{\bm{l}}^{0}\rangle\langle\Psi_{\bm{l}}^{0}|W|\Psi_{\bm{k}}^{0}\rangle}{(E_{\bm{k}}^{0}-E_{\bm{p}}^{0})(E_{\bm{k}}^{0}-E_{\bm{l}}^{0})}$
$\displaystyle-\sum_{{\bm{p}}\neq{\bm{k}}}|\Psi_{\bm{p}}^{0}\rangle\frac{\langle\Psi_{\bm{k}}^{0}|W|\Psi_{\bm{k}}^{0}\rangle\langle\Psi_{\bm{p}}^{0}|W|\Psi_{\bm{k}}^{0}\rangle}{(E_{\bm{k}}^{0}-E_{\bm{p}}^{0})^{2}}-\frac{1}{2}|\Psi_{\bm{k}}^{0}\rangle\sum_{{\bm{p}}\neq{\bm{k}}}\frac{|\langle\Psi_{\bm{p}}^{0}|W|\Psi_{\bm{k}}^{0}\rangle|^{2}}{(E_{\bm{k}}^{0}-E_{\bm{p}}^{0})^{2}},$
(1)

where ‘0’ indicates the original wavefunction of the graphene in the absence
of the SiC potential. The second term in Eq. (1) is the first-order
correction. To have non-vanishing coefficient
$\langle\Psi_{\bm{p}}^{0}|W|\Psi_{\bm{k}}^{0}\rangle$, the difference between
${\bm{p}}$ and ${\bm{k}}$ must align with the period of $W({\bm{x}})$. We note
that we simulate the SiC potential in terms of the simplest harmonic form,
that is, $W({\bm{x}})=w\big{(}\cos({\bm{b_{1}\cdot x}})+\cos({\bm{b_{2}\cdot
x}})+\cos(-({\bm{b_{1}}}+{\bm{b_{2}}){\bm{\cdot}}{\bm{x}}})\big{)}$. Since the
period of $W({\bm{x}})$ is described by the two reciprocal lattice vectors
${\bm{b}}_{1}$ and ${\bm{b}}_{2}$, the non-zero density of states appears only
at ${\bm{p}}={\bm{k}}\pm{\bm{b}}_{1}$, ${\bm{k}}\pm{\bm{b}}_{2}$, and
${\bm{k}}\pm({\bm{b}}_{1}+{\bm{b}}_{2})$. As ${\bm{k}}$ represents the
momentum of the Dirac states, there are 12 duplications of the first order
(the green and blue clones in Fig. 2e) within the Brillouin zone of graphene,
which is in agreement with the experimental observation.

The last three terms of $|\Psi_{\bm{k}}\rangle$ correspond to the second-order
perturbations. The second last term on the right side of Eq. (1) vanishes due
to the fact that $\langle\Psi_{\bm{k}}^{0}|W|\Psi_{\bm{k}}^{0}\rangle=0$. The
last term only induces a renormalization of the primary cone at ${\bm{k}}$.
Only the first term of the second-order survives in certain conditions and
give rise to clones in momentum space. To have non-vanishing
$\langle\Psi_{\bm{l}}^{0}|W|\Psi_{\bm{k}}^{0}\rangle$, the mediating momentum
obeys
${\bm{l}}={\bm{k}}\pm{\bm{b}}_{1},{\bm{k}}\pm{\bm{b}}_{2},{\bm{k}}\pm{\bm{b}}_{1}\pm{\bm{b}}_{2}$.
Likewise, the momentum of the final wavefunction satisfies
${\bm{p}}={\bm{k}}\pm 2{\bm{b}}_{1},{\bm{k}}\pm 2{\bm{b}}_{2},{\bm{k}}\pm
2({\bm{b}}_{1}+{\bm{b}}_{2}),{\bm{k}}\pm({\bm{b}}_{1}-{\bm{b}}_{2}),{\bm{k}}\pm(2{\bm{b}}_{1}+{\bm{b}}_{2}),{\bm{k}}\pm({\bm{b}}_{1}+2{\bm{b}}_{2})$.
In this regard, the second-order perturbations duplicate 24 Dirac cones (the
gray and yellow clones in Fig. 2e) at various ${\bm{p}}$, consistent with the
experimental observation.

The perturbative corrections to the wavefunction in Eq. (1) give rise to the
clones of Dirac cones. The clones represent a redistribution of spectral
weight of the primary Dirac cone in the momentum space. That is why the
observed clones share the same band dispersion as the primary Dirac cones. The
clones derived from the same primary Dirac cone do not hybridize with each
other. On the other hand, the hybridization are allowed for the clone or
primary contours from different valleys. The hybridization between two valleys
is mediated by the substrate potential. Here the two valleys of monolayer
graphene behave like the two sets of Dirac cones from the two layers of TBG.
Flat bands can be created under certain substrate conditions, which can be
seen in the following discussions.

Our ARPES and tight-binding results indicate that the substrate potential
places a periodic modulation to the graphene band structure, produces clones
of Dirac states, and effectively shortens the distance between the two valleys
by the reciporical vectors of the substrate. This machanism can enable a
direct coupling between the Dirac states from the two valleys when the
substrate lattice is nearly commensurate with graphene. Up to date, various
graphene-based heterostructures such as graphene/metals Rotenberg2015 ;
HERNANDEZRODRIGUEZ201558 ; YU2019633 , graphene/boron nitride Yankowitz2019 ;
Dean2010 ; Haigh2012 , and graphene/chalcogenide compounds grapheneVDW ;
Aeschlimanneaay0761 ; Mao2020 ; Geim2013 ; Ponomarenko2011 ; Georgiou2013
have been experimentally realized. To investigate the substrate effects in the
nearly commensurate condition, we performed tight-binding simulation for a
generic graphene heterostructure with a hexagonal substrate rotated by 30∘
relative to the graphene unit cell. The substrate lattice constant is chosen
to be 3.8 Å, which is about 10% smaller than the commensurate value
$\sqrt{3}a_{\textrm{Gr}}=4.26\ \textrm{\AA}$. The cacluated band structure is
shown in Fig. 4a. At the Fermi level, there are two primary Dirac points (DP)
denoted by D and D’ and six duplicated DPs denoted by C1, C1’, C2, C2’, C3,
and C3’. The clones of C1-C3 are from the valley of the ‘D’ Dirac cone while
those of C1’-C3’ are from the other valley. When two Dirac bands from
different valleys (for example, C1’ and D, or C2 and C2’) intersect, an energy
gap is opened at the crossing point. The gapped band structure give rise to
Van Hove singularities (VHS) in the density of states (DOS) as marked by the
yellow and red arrows. The iso-energy contours at $E=-0.08$ eV is plotted in
Fig. 4b. Close to the zero energy, all the primary and cloned contours are
isolated in momentum space and thus contribute to the DOS as independent Dirac
cones. Therefore, the DOS vanishes at zero energy as shown in Fig. 4a. The
effective distance between the two primary Dirac points in the presence of
substrate perturbations is

$q=\big{|}|{\textrm{K}^{\prime}}_{\textrm{Gr}}{\textrm{K}}_{\textrm{Gr}}|-b\big{|}=\big{|}|\Gamma{\textrm{K}}_{\textrm{Gr}}|-b\big{|},$
(2)

where $b$ is the length of the substrate reciprocal lattice vector, as
schematically shown in Fig. 4c. That is also the separation between DPs C1’
and D in Fig. 4b. As the substrate constant approaches the commensurate value
$\sqrt{3}a_{\textrm{Gr}}$, the clone contours move closer to the the primary
cone and enhance the coupling between the two valleys, since the effective
coupling is described by a dimensionless parameter $\alpha=\frac{w}{\hbar
v_{\rm F}q}$, where $w$ is the amplitude of the substrate potential and
$v_{\rm F}$ is the Fermi velocity of electrons in graphene Bistritzer12233 .
$w$ is set to be $0.05t=140$ meV (comparable to that of TBG, $w\approx 110$
meV)Bistritzer12233 , where $t=2.8$ eV is the nearest-neighbor hopping
parameter of grapheneCastroNeto2009 . When the substrate lattice constant is
equal to 4.166 Å and 4.370 Å, a sharp peak shows up at zero energy in DOS and
an energy gap of size $\sim 2w$ emerges between the conduction and valence
bands, as shown in Figs. 4e,f. The zero-energy peaks in DOS cannot be
described by isolated Dirac cones, therefore there must be dispersionless
bands emerging at low energy as a consequence of hybridization of Dirac states
from the two valleys.

The effective separation between the two DPs of monolayer gaphene in the
presence of substrate modulation is described by three vectors, ${\bm{q}_{1}}$
= $q$(1, 0), ${\bm{q}_{2}}$ = $q$($\frac{-1}{2}$, $\frac{\sqrt{3}}{2}$), and
${\bm{q}_{3}}$ = $q$($\frac{-1}{2}$, $\frac{-\sqrt{3}}{2}$). Repeated hopping
between the two valleys generates a $k$-space honeycomb lattice shown in Fig.
5a. The unit vectors of this lattice are same as reciprocal vectors of the
Moire pattern (MP) of this hybrid structure, namely, ${\bm{b}}^{\rm
MP}_{1}={\bm{q}_{1}}-{\bm{q}_{2}}$ and ${\bm{b}}^{\rm
MP}_{2}={\bm{q}_{1}}-{\bm{q}_{3}}$. The low-energy electron dynamics can be
described by an effective Hamiltonian,

$H_{\rm eff}({\bm{k}})=\begin{bmatrix}D({\bm{k}})&T_{1}&T_{2}&T_{3}\\\
T^{\dagger}_{1}&\overline{D}({\bm{k}-{\bm{q}_{1}}})&0&0\\\
T^{\dagger}_{2}&0&\overline{D}({\bm{k}-{\bm{q}_{2}}})&0\\\
T^{\dagger}_{3}&0&0&\overline{D}({\bm{k}-{\bm{q}_{3}}})\end{bmatrix},$ (3)

where $D({\bm{k}})=\hbar v_{F}{\bm{k}}{\bm{\cdot}}{\bm{\sigma}}$,
$\overline{D}({\bm{k}})=-\hbar v_{F}{\bm{k}}{\bm{\cdot}}{\bm{\sigma}^{*}}$,
and $T_{m}=w\exp({i\frac{2(1-m)\pi}{3}}\sigma_{z})$, $m=$1, 2, and 3. Here
${\bm{\sigma}}=(\sigma_{x},\sigma_{y}$) and $\sigma_{z}$ are Pauli matrices
associated with the A and B sublattices of graphene. We note that
$D({\bm{k}})$ and $\overline{D}({\bm{k}})$ describes the Dirac cone at $\rm K$
and $\rm K^{\prime}$, respectively, as the two cones are time-reversal
partners. The derivation of the effective Hamiltonian can be found in the
supplementary information. Using this effective Hamiltonian, we calculated the
band structure along the path A-B-C-D-A and the DOS, see Fig. 5b. For $\alpha$
= 0.1, the DPs at B and C remained isolated meanwhile an energy gap is opened
at the crossing point of two Dirac bands. The DOS shows a nearly linear
dependence of energy and several peaks from VHSs, which is consistent with the
tight-binding result in Fig. 4a. When $\alpha$ increases to 0.586, a pair of
absolutely flat bands exist at zero energy inside the band gap ($\Delta E\sim
2w$). This leads to a sharp zero-energy peak in DOS, which agrees remarkably
well with the tight-binding DOS with $a_{\rm sub}$ = 4.370 Å and $w=0.05t$ as
shown in the inset of Fig. 5b. The perfect flatness of zero-energy bands is
due to the chiral symmetry of the Hamiltonian, since $H_{\rm eff}$ is
equivalent to the chiral-symmetric continuum model proposed by Tarnopolsky et
al.PhysRevLett.122.106405  (see the proof in the supplementary information).
The intervalley transition matrices $T_{m}$ ($m$ = 1, 2, and 3) in $H_{\rm
eff}$ contain only the diagonal AA and BB couplings due to the on-site
substrate potential and the fact that the wavefunctions of Dirac states at K
and K’ are defined with respect to the same A and B sublattices of monolayer
graphene. The absence of AB couplings yields the chiral symmetry of the
effective Hamiltonian and creates the absolutely flat bands. In other words,
this system is a natural realization of the chiral-symmetric model of flat
bandsPhysRevLett.122.106405 . The unique coupling $\alpha^{*}=0.586$
corresponds to two ‘magic’ lattice constants according to Eq. 2,

$a^{*}_{\rm sub}=\sqrt{3}a_{\rm Gr}\pm\frac{3wa_{\rm Gr}}{2\pi t\alpha^{*}}.$
(4)

Plug in $\alpha^{*}=0.586$, $a_{\rm Gr}$ = 2.46 Å and $w=0.05t$, we find that
$a^{*}_{\rm sub}$ = 4.161 and 4.361 Å. The magic lattice constants echo the
values (4.166 and 4.370 Å) we found in the tight-binding simulations. The
small discrepancy between two results can be attributed to the finite size of
the supercell we used in the tight-binding simulations. The effective
Hamiltonian gives rise to a series of magic coupling $a^{*}$ with a
periodicity of $\Delta\alpha\simeq 1.5$ Bistritzer12233 ;
PhysRevLett.126.016404 . The second magic coupling is $a^{*}$ = 2.221. The
band structure with this magic coupling (the bottom panel of Fig. 5b) exhibits
absolute flat bands at zero energy and a smaller band gap ($\Delta E\sim
0.2w$). We notice that other bands such as the ones at $\pm$0.25$w$ are also
flattened, leading sharp spikes in the DOS. The second and higher magic
couplings correspond to a very small deviation from the commensurate lattice
constant, $\Delta a\lesssim 0.02$ Å. It is technically challenging to detect
such small lattice deviations in experiments, just like the smaller twists
corresponding the higher magic couplings in the TBG systems.

The low-energy dynamics in the graphene heterostructures is essentially
governed by the coupling $\alpha$, and $\alpha$ is determined by the
separation $q$ between two adjacent Dirac cones in the Moiré lattice. For
$w=0.05t$ ($\sim$140 meV) and $\alpha^{*}=0.586$, $q$ = 0.04
${\rm\AA}^{-1}\ll|\Gamma{\rm K}_{\rm Gr}|=1.7\ {\rm\AA}^{-1}$, which means the
reciprocal lattice vectors of substrate must almost connect the two DPs at K
and K’ in momentum space, see Fig. 4c. This places a constraint on the
possible substrate lattice constants and orientations. To have flat bands, the
largest possible substrate lattice constant corresponds to a
$(\sqrt{3}\times\sqrt{3})R30^{\circ}$ supercell of graphene, $i.e.$, $a_{\rm
sub}\approx\sqrt{3}a_{\rm Gr}=4.26$ Å. This commensurate relation is also
known as the Kekulé superlatticePhysRevB.88.155415 . For other commensurate
relations between graphene and substrates, the substrate lattice constant has
to be $\lesssim\frac{\sqrt{3}}{2}a_{\rm Gr}=2.13$ Å, which is very rare in
real materials. Therefore, the substrate materials for the flat-band
heterostructure must have a surface with $C_{3}$ rotation symmetry and lattice
constant close to 4.26 Å. In Table 1, we list several materials which can be
potentially employed in the proposed heterostructures. In addition to the
conventional assembly method by growing or transferring graphene onto the
substrate surface, the proposed flat-band heterostructures can be readily
synthesized via a “top-down” approach, as schematically plotted in Fig. 5c.
High-quality graphene layers can be epitxially grown on the SiC(0001) surface,
and the Dirac states of graphene remain isolated due to the large lattice
mismatch. Therefore, the graphene/SiC structure can serve as a supporting
substrate for the growth of various materials with nearly commensurate
relations as suggested in Table 1. For example, Bi2Se3 ($a$ = 4.136 Å) thin
layers have been grown on the graphene/SiC(0001) surface Jin2013 .

In conclusion, our ARPES experiments demonstrated the cloning of Dirac
fermions in the graphene/SiC(0001) heterostructure due to the periodic
modulation of the substrate potential. Our theoretical calculations showed
this modulation effect from the substrate can effectively couple the two
valleys of Dirac states in monolayer graphene in the nearly commensurate
condition. The graphene heterostructures can be a promising alternate system
for exploring the intriguing flat-band physics that was found in TBG. The
criterion for realizing flat bands is a matchup of the surface potential
strength and the periodicity of the substrate to reach the magic effective
coupling $\alpha^{*}$. There are a vast selections of substrate materials that
can be potentially used in this hybrid structure. In addition, the charge and
spin orderings in the substrates such as antiferromagnetism in
MnTePhysRevB.96.214418 , superconductivity in PdTe2PhysRevB.97.014523 , and
topological surface states in Bi2Se3Jin2013  can further enrich flat-band
physics in graphene via proximity effects, systematic investigations of which
is left to future works.

Materials |  $a_{\rm sub}$ (Å) |  Space group |  Materials |  $a_{\rm sub}$ (Å) |  Space group
---|---|---|---|---|---
CuSe | 3.980 | $P6_{3}/mmc$ |  Bi2Se3 | 4.136 |  $R$-$3mH$
InSe | 4.000 |  $R$-$3mH$ |  CdS | 4.137 | $P6_{3}/mc$
CrTe | 4.005 | $P6_{3}/mmc$ |  MnTe | 4.148 | $P6_{3}/mmc$
PdTe2 | 4.024 |  $P$-$3m1$ |  GeTe | 4.156 |  $R$-$3mH$
In2Se3 | 4.026 |  $R$-$3mH$ |  PdTe | 4.200 | $P6_{3}/mmc$
PtTe2 | 4.026 |  $P$-$3m1$ |  CdSe | 4.232 | $P6_{3}/mc$
InSe | 4.050 | $P6_{3}/mmc$ |  Cu2Te | 4.237 | $P6/mmm$
As2Te3 | 4.058 | $R3mH$ |  Sb2Te3 | 4.264 |  $R$-$3mH$
GaTe | 4.060 | $P6_{3}/mmc$ |  SiTe2 | 4.289 |  $P$-$3m1$
ZnTe | 4.092 | $P3_{1}21$ |  MgSe | 4.319 | $P6_{3}/mmc$
ScTe | 4.097 | $P6_{3}/mmc$ |  HgSe | 4.320 | $P3_{2}21$
AuTe2 | 4.107 |  $P$-$3m1$ |  Sb(111) | 4.332 |  $R$-$3mH$
PtTe | 4.111 | $P6_{3}/mmc$ |  HgTe | 4.392 | $P3_{1}21$
AuSe | 4.120 | $P6_{3}/mmc$ |  Bi2Te3 | 4.403 |  $R$-$3mH$
MnSe | 4.120 | $P6_{3}/mc$ |  MgTe | 4.531 | $P6_{3}/mmc$
Cu2Se | 4.132 |  $R$-$3mH$ |  Bi(111) | 4.546 |  $R$-$3mH$
Table 1: Material candidates for flat-band heterostructures. The materials are
ordered according to their in-plane lattice constant. Bi(111) and Sb(111) are
Bi and Sb films grown in the rhombohedral (111) directionBIAN2019109 .

##  IV Methods

###  IV.1 Tight-bind modeling

To simulate the observed spectrum of the graphene on the top of the SiC film,
we use the simplest graphene model with the inclusion of only the nearest
neighbor hopping. Since SiC is an insulator with large gaps, the main low
energy physics can be described by the graphene monolayer embedded in the
electric potential on the SiC surface. Although the details of the SiC surface
potential is unknown, it suffices to employ an approximate potential that
preserves the crystalline symmetry. The SiC surface belongs to the wallpaper
group $P3m1$, thus the SiC potential can be written in this simplest form

$\displaystyle W(x,y)=$ $\displaystyle
w\big{(}\cos(\pi\frac{4y}{3b})+\cos(\pi\frac{-2\sqrt{3}x-2y}{3b})+\cos(\pi\frac{2\sqrt{3}x-2y}{3b})\big{)}$
(5)

where $b$ (1.77 Å) is the in-plane silicon-carbon distance on the SiC(0001)
surface and ${\bm{b}}_{1,2}=\frac{2\pi}{3b}(\pm\sqrt{3},-1)$. That is, the SiC
potential includes only the first order of the Fourier series and breaks the
$C_{6}$ rotation symmetry due to the different potentials stemming from Si and
C atoms. We can write the effective low-energy Hamiltonian in the second
quantization form

$\displaystyle\hat{H}=$
$\displaystyle\sum_{n,m}\Big{[}t(a^{\dagger}_{n,m}b_{n,m}+a^{\dagger}_{n+1,m}b_{n,m}+a^{\dagger}_{n,m+1}b_{n,m}+{\rm{h.c.}})+W(x_{nm}^{+},y_{nm}^{+})a^{\dagger}_{n,m}a_{n,m}$
$\displaystyle+W(x_{nm}^{-},y_{nm}^{-})b^{\dagger}_{n,m}b_{n,m}\Big{]},$ (6)

where $x_{nm}^{\pm}=(\mp 1+3n-3m)a/2,\ y_{nm}^{\pm}=(\sqrt{3}n+\sqrt{3}m)a/2$
and $a\approx 1.42$ Å is the carbon-carbon distance in graphene. A
400$\times$400 supercell is employed in the calculations of band structure and
DOS. By diagonalizing the Hamiltonian, we find the eigenstates within
$E\pm\delta$. Then we transform the eigenstates to momentum space so that at
energy $E$, the density of the wavefunctions can be plotted in momentum space.

###  IV.2 Sample synthesis and characterizations

The graphene films were prepared by annealing a 6H-SiC(0001) substrate at 1150
°C in an integrated MBE-STM-ARPES ultrahigh vacuum (UHV) system with base
pressure below 2$\times$10-10 mbar. After the growth, the graphene samples
were in-situ transferred the ARPES stage. ARPES measurements were performed at
100 K using a SPECS PHOIBOS-150 hemisphere analyzer with a SPECS UVS-300
helium discharge lamp (photon energy = 21.2 eV). The size of the beam spot on
the sample was $\sim$1.5 mm. The topography of the sample surface was mapped
in-situ by an Aarhus STM equipped in the growth chamber.

The TEM samples were prepared by a lift-out method in a ThermoFisher
Scientific Scios focused ion beam (FIB) instrument at room temperature, and
imaged in the ThermoFisher Scientific G2 Tecnai F30 FEG high resolution TEM
operated at 300 kV. The SiC substrate was tilted to the [100] zone axis and
the lattice fringes from both the graphene and the SiC can be clearly
resolved. Great care has been taken to reduce the beam damage on the thin film
samples both during the FIB lift out and during the sample tilting and high-
resolution image acquisition process.

Figure 1: Lattice structure and Brillouin zone of graphene/SiC
heterostructure. a, Lattice structure of graphene and SiC(0001) surface. The
unit cell of graphene and SiC surface is depicted by the green and black
parallelograms, respectively. b, Brillouin zone of graphene (blue lines) and
SiC(0001) surface (red lines). c, STM image of graphene with atomic
resolution. A Moiré pattern is observed with a period approximately equal to 6
times the lattice constant of SiC(0001). d, Cross-sectional TEM image of
graphene/SiC heterostructure.

Figure 2: ARPES spectrum and momentum-space analysis of Dirac fermion cloning.
a, Fermi surface of graphene/SiC(0001) heterostructure measured by ARPES. b,
Geometrical relation between the Brillouin zones of graphene and SiC(0001)
surface. c, The location of green and blue clones from the first-order
perturbations. d, The location of yellow and gray clones from the second order
perturbations. e, Distribution of first-order and second-order clones in the
Brillouin zone of graphene.

Figure 3: ARPES spectra and tight-binding simulation of Dirac bands and their
clones. a, ARPES and b, tight-binding fermi surface of graphene/SiC
heterostructure. c, ARPES and d, tight-binding spectrum taken along ‘cut1’
marked in a. e, ARPES spectra taken along ‘cut2-5’. f, ARPES spectrum of the
primary Dirac cone taken along ‘cut2’ from 3-, 5-monolayer graphene films
grown on SiC(0001) surface. g, ARPES spectrum of cloned Dirac bands taken
along ‘cut4’ from 3-, 5-monolayer graphene films.

Figure 4: Tight-binding simulation of epitaxial graphene with various
substrate lattice constants. a, The tight-binding band structure and density
of states with $a_{\textrm{sub}}=3.8\ \textrm{\AA}$. b, Calculated iso-energy
contours at $E=-0.08\ eV$. c, Coupling mechanism between the two Dirac cones
at ${\textrm{K}_{\textrm{Gr}}}$ and ${\textrm{K}^{\prime}_{\textrm{Gr}}}$. d,
Schematic of the movement of clone contours C1’ as the substrate approaches to
the commensurate value $\sqrt{3}a_{\textrm{Gr}}=4.26\ \textrm{\AA}$. e,
Calculated density of states for $a_{\textrm{sub}}=4.166\ \textrm{\AA}$ and
$w=0.05t$. f, Same as e but for $a_{\textrm{sub}}=4.370\ \textrm{\AA}$.

Figure 5: Effective model of Moiré mini lattice and flat bands. a, Moiré
Brillouin zone of a nearly commensurate graphene heterostructure. b,
Calculated band structure and density of states with $\alpha=0.1,0.586$, and
$2.221$. The DOS plot in the middle panel includes an inset showing the tight-
binding DOS (the blue curve) from Fig. 4f. c, Schematic of the two ways of
assembling heterostructures with flat bands.

## References

  * (1) Castro Neto, A. H., Guinea, F., Peres, N. M. R., Novoselov, K. S. & Geim, A. K.  The electronic properties of graphene.  _Reviews of Modern Physics_ 81, 109–162 (2009).
  * (2) Novoselov, K. S. _et al._ Two-dimensional gas of massless Dirac fermions in graphene.  _Nature_ 438, 197–200 (2005).
  * (3) Zhang, Y., Tan, Y.-W., Stormer, H. L. & Kim, P.  Experimental observation of the quantum Hall effect and Berry’s phase in graphene.  _Nature_ 438, 201–204 (2005).
  * (4) Geim, A. K. & Grigorieva, I. V.  Van der Waals heterostructures.  _Nature_ 499, 419–425 (2013).
  * (5) Katsnelson, M. I., Novoselov, K. S. & Geim, A. K.  Chiral tunnelling and the Klein paradox in graphene.  _Nature Physics_ 2, 620–625 (2006).
  * (6) Basov, D. N., Fogler, M. M., Lanzara, A., Wang, F. & Zhang, Y.  Colloquium: Graphene spectroscopy.  _Rev. Mod. Phys._ 86, 959–994 (2014).
  * (7) Haldane, F. D. M.  Model for a quantum hall effect without landau levels: Condensed-matter realization of the ”parity anomaly”.  _Phys. Rev. Lett._ 61, 2015–2018 (1988).
  * (8) Hasan, M. Z. & Kane, C. L.  Colloquium: Topological insulators.  _Rev. Mod. Phys._ 82, 3045–3067 (2010).
  * (9) Qi, X.-L. & Zhang, S.-C.  Topological insulators and superconductors.  _Rev. Mod. Phys._ 83, 1057–1110 (2011).
  * (10) Liu, Y., Bian, G., Miller, T. & Chiang, T.-C.  Visualizing Electronic Chirality and Berry Phases in Graphene Systems Using Photoemission with Circularly Polarized Light.  _Physical Review Letters_ 107, 166803 (2011).
  * (11) Bistritzer, R. & MacDonald, A. H.  Moiré bands in twisted double-layer graphene.  _Proceedings of the National Academy of Sciences_ 108, 12233–12237 (2011).
  * (12) Ponomarenko, L. A. _et al._ Cloning of Dirac fermions in graphene superlattices.  _Nature_ 497, 594–597 (2013).
  * (13) Dean, C. R. _et al._ Hofstadter’s butterfly and the fractal quantum Hall effect in moiré superlattices.  _Nature_ 497, 598–602 (2013).
  * (14) Cao, Y. _et al._ Correlated insulator behaviour at half-filling in magic-angle graphene superlattices.  _Nature_ 556, 80–84 (2018).
  * (15) Cao, Y. _et al._ Unconventional superconductivity in magic-angle graphene superlattices.  _Nature_ 556, 43–50 (2018).
  * (16) Wang, E. _et al._ Gaps induced by inversion symmetry breaking and second-generation Dirac cones in graphene/hexagonal boron nitride.  _Nature Physics_ 12, 1111–1115 (2016).
  * (17) Sharpe, A. L. _et al._ Emergent ferromagnetism near three-quarters filling in twisted bilayer graphene.  _Science_ 365, 605–608 (2019).
  * (18) Serlin, M. _et al._ Intrinsic quantized anomalous hall effect in a moiré heterostructure.  _Science_ 367, 900–903 (2020).
  * (19) Tritsaris, G. A. _et al._ Electronic structure calculations of twisted multi-layer graphene superlattices.  _2D Materials_ 7, 035028 (2020).
  * (20) Tarnopolsky, G., Kruchkov, A. J. & Vishwanath, A.  Origin of magic angles in twisted bilayer graphene.  _Phys. Rev. Lett._ 122, 106405 (2019).
  * (21) Conrad, M. _et al._ Wide band gap semiconductor from a hidden 2D incommensurate graphene phase.  _Nano Letters_ 17, 341–347 (2017).
  * (22) Huang, L. _et al._ Effects of moiré lattice structure on electronic properties of graphene.  _Phys. Rev. B_ 96, 035411 (2017).
  * (23) Zhou, S. Y. _et al._ First direct observation of Dirac fermions in graphite.  _Nature Physics_ 2, 595–599 (2006).
  * (24) Liu, Y. _et al._ Phonon-induced gaps in graphene and graphite observed by angle-resolved photoemission.  _Physical Review Letters_ 105, 1–4 (2010).
  * (25) Rotenberg, E. & Bostwick, A.  Superlattice effects in graphene on SiC(0001) and Ir(111) probed by ARPES.  _Synthetic Metals_ 210, 85–94 (2015).
  * (26) Hernández-Rodríguez, I., García, J. M., Martín-Gago, J. A., de Andrés, P. L. & Méndez, J.  Graphene growth on pt(111) and au(111) using a mbe carbon solid-source.  _Diamond and Related Materials_ 57, 58 – 62 (2015).
  * (27) Yu, J. _et al._ Study on aln buffer layer for gan on graphene/copper sheet grown by mbe at low growth temperature.  _Journal of Alloys and Compounds_ 783, 633 – 642 (2019).
  * (28) Yankowitz, M., Ma, Q., Jarillo-Herrero, P. & LeRoy, B. J.  van der Waals heterostructures combining graphene and hexagonal boron nitride.  _Nature Reviews Physics_ 1, 112–125 (2019).
  * (29) Dean, C. R. _et al._ Boron nitride substrates for high-quality graphene electronics.  _Nature Nanotechnology_ 5, 722–726 (2010).
  * (30) Haigh, S. J. _et al._ Cross-sectional imaging of individual layers and buried interfaces of graphene-based heterostructures and superlattices.  _Nature Materials_ 11, 764–767 (2012).
  * (31) Zhang, Z. _et al._ Graphene-based mixed-dimensional van der waals heterostructures for advanced optoelectronics.  _Advanced Materials_ 31, 1806411 (2019).
  * (32) Aeschlimann, S. _et al._ Direct evidence for efficient ultrafast charge separation in epitaxial ws2/graphene heterostructures.  _Science Advances_ 6 (2020).
  * (33) Mao, J. _et al._ Evidence of flat bands and correlated states in buckled graphene superlattices.  _Nature_ 584, 215–220 (2020).
  * (34) Ponomarenko, L. A. _et al._ Tunable metal–insulator transition in double-layer graphene heterostructures.  _Nature Physics_ 7, 958–961 (2011).
  * (35) Georgiou, T. _et al._ Vertical field-effect transistor based on graphene–WS2 heterostructures for flexible and transparent electronics.  _Nature Nanotechnology_ 8, 100–103 (2013).
  * (36) Ren, Y., Gao, Q., MacDonald, A. H. & Niu, Q.  Wkb estimate of bilayer graphene’s magic twist angles.  _Phys. Rev. Lett._ 126, 016404 (2021).
  * (37) Wallbank, J. R., Mucha-Kruczyński, M. & Fal’ko, V. I.  Moiré minibands in graphene heterostructures with almost commensurate $\sqrt{3}\times{}\sqrt{3}$ hexagonal crystals.  _Phys. Rev. B_ 88, 155415 (2013).
  * (38) Jin, K.-H. & Jhi, S.-H.  Proximity-induced giant spin-orbit interaction in epitaxial graphene on a topological insulator.  _Physical Review B_ 87, 075442 (2013).
  * (39) Kriegner, D. _et al._ Magnetic anisotropy in antiferromagnetic hexagonal mnte.  _Phys. Rev. B_ 96, 214418 (2017).
  * (40) Das, S. _et al._ Conventional superconductivity in the type-ii dirac semimetal ${\mathrm{pdte}}_{2}$.  _Phys. Rev. B_ 97, 014523 (2018).
  * (41) Bian, G. _et al._ Survey of electronic structure of bi and sb thin films by first-principles calculations and photoemission measurements.  _Journal of Physics and Chemistry of Solids_ 128, 109–117 (2019).


# A Gallery of Gaussian Periods

Ellen Eischen1 and Stephan Ramon Garcia2
1Department of Mathematics, University of Oregon<EMAIL_ADDRESS>
2Department of Mathematics, Pomona College<EMAIL_ADDRESS>
###### Abstract

Gaussian periods are certain sums of roots of unity whose study dates back to
Gauss’s seminal work in algebra and number theory. Recently, large scale plots
of Gaussian periods have been revealed to exhibit striking visual patterns,
some of which have been explored in the second named author’s prior work. In
2020, the first named author produced a new app, Gaussian periods, which
allows anyone to create these plots much more efficiently and at a larger
scale than before. In this paper, we introduce Gaussian periods, present
illustrations created with the new app, and summarize how mathematics controls
some visual features, including colorings left unexplained in earlier work.

## Introduction

Gaussian periods, certain sums of roots of unity introduced by Gauss, have
played a key role in several mathematical developments. For example, Gauss
employed them in his work on constructibility of regular polygons by unmarked
straightedge and compass, as well as in number theory. In the past few years,
realizations of their role in the supercharacter theory of Diaconis and Isaacs
have led in new directions [3].

The fast computations afforded by modern technology provide new insights into
Gaussian periods by enabling us to study them at a scale that was until
recently unfathomable. In particular, large-scale plots of Gaussian periods
display striking visual properties that were previously unknown (see, e.g.,
Figure 1). The goals of this paper lie in the intersection of the mathematical
and artistic aspects of Gaussian periods:

  * •

Introduce Gaussian periods, together with new illustrations representing some
of their visual features.

  * •

Introduce a new app, created in 2020, that quickly produces illustrations of
Gaussian periods and is appropriate as a tool for art, exploratory
mathematical research, illustration at scale, and pedagogy.

  * •

Summarize how mathematics controls some features, including colorings
unexplained in earlier work.

(a) $n=1481151$, $\omega=54184$, $c=21$

(b) $n=255255$, $\omega=254$, $c=7$

Figure 1: Examples of $\mathsf{G}(n,\omega)$ for different choices of input.

Before proceeding, we clarify what we mean by the term Gaussian period. Given
a positive integer $n$, an integer $\omega$ coprime to $n$, and an integer
$k$, we set

$\displaystyle\eta_{n,\omega,k}:=\sum_{j=0}^{d-1}e^{\frac{2\pi
i\omega^{j}k}{n}},$

where $d$ is the multiplicative order of $\omega$ mod $n$; that is, $d$ is the
smallest positive integer such that $\omega^{d}\equiv
1\,(\operatorname{mod}n)$. When $k$ is relatively prime to $n$, we say that
$\eta_{n,\omega,k}$ is a Gaussian period of modulus $n$ and generator
$\omega$. In this paper, by a slight abuse of notation that coincides with the
conventions employed in [7, 4, 8], we still call $\eta_{n,\omega,k}$ a
Gaussian period of modulus $n$ and generator $\omega$ even when $k$ is not
coprime to $n$. As we explain later in this paper, $\eta_{n,\omega,k}$ is a
positive integer multiple of the Gaussian period
$\eta_{\frac{n}{\gcd(n,k)},\omega,k}$, which allows us to preserve certain key
structures of interest.

## Creating and viewing plots of Gaussian periods

Each of the images in this paper is a plot in the complex plane of

$\displaystyle\mathsf{G}(n,\omega):=\left\\{\eta_{n,\omega,k}:k=1,2,\ldots,n\right\\}\subset\mathbb{C}.$

Mathematical principles guarantee that $\mathsf{G}(n,\omega)$ exhibits certain
basic symmetries. For large values of $n$ (which we refer to as large scale),
however, plots of $\mathsf{G}(n,\omega)$ also exhibit striking patterns and
intricacies whose aesthetic properties can be appreciated even by those
without mathematical training.

That large-scale plots of Gaussian periods exhibit such variety and intricate
patterns was a surprise. This was discovered by B. Lutz, an undergraduate, in
the course of his senior thesis. In his explorations, he also introduced a
coloring scheme, which is employed in the plots in this paper and discussed in
the next section.

A new app, Gaussian Periods (written in Swift for Apple computers and freely
available [6]) was produced in 2020 by the first named author, with assistance
from R. Lipshitz. This app:

  * •

plots Gaussian periods faster than previous code (including the aforementioned
Mathematica code), taking seconds to produce images that used to take hours;

  * •

allows larger scale plots than previously possible (e.g., Figure 4, which
contains over 9 million points), which can be useful for exploring or
illustrating asymptotic behavior;

  * •

does not require programming experience or mathematical expertise;

  * •

allows one to quickly modify values of $n$ and $\omega$, as well as a coloring
parameter $c$; and

  * •

includes an option to save layers suitable for further steps, e.g.,
manipulation in Adobe Photoshop to customize color choices (as was done to
produce the images in this paper).

As a result, the app is suitable for projects ranging from art to exploratory
mathematical investigations to illustration. We also used it to produce all
the images in this paper. To improve the image quality in this paper, we
layered the different color components (plotting all the points corresponding
to a color at the same time), which is assisted by the layers option in the
app.

While notions of “beauty” and “aesthetic appeal” are subjective, symmetry has
often appeared in discussions of beauty in both art and math, dating back to
the ancient Babylonians. In a historical context, given the algebraic origins
of Gaussian periods, it is fitting to note that symmetry originally became a
prominent concept in mathematics not through geometry, but through algebra [9,
Preface]. Gaussian periods, together with the various symmetries they exhibit,
can also be considered of independent artistic merit and be appreciated in
their own right (even by those with no mathematical training). More broadly,
patterns in plots of certain other families of algebraic numbers have also
been recognized for their beauty [1]. For those who wish to delve further into
philosophical considerations of the artistic merits of Gaussian periods, a
brief survey of visual aesthetics in similar mathematical contexts can be
found in [2, pp. 121-4].

## Using color to reveal structures

Whether focusing on art or math, a coloring scheme can be used to highlight
some structures in $\mathsf{G}(n,\omega)$ for a given $n$ and $\omega$.
Following [7, §3], fix a positive integer $c\mid n$ and, for $j=1,2,\ldots,c$,
assign the same color to all points in the set $\\{\eta_{n,\omega,k}:k\equiv
j\,(\operatorname{mod}c)\\}$. Two points $\eta_{n,\omega,k}$ and
$\eta_{n,\omega,\ell}$ might have the same color even if
$k\nequiv\ell\,(\operatorname{mod}c)$, since
$\eta_{n,\omega,k}=\eta_{n,\omega,k\omega^{j}}$ for all integers $j$.

Since this is the coloring scheme employed in each of the images here (and is
implemented in Gaussian Periods), we begin with small-scale examples to help
readers grasp the subtleties of this coloring scheme.

######  Example 1.

Suppose $n=27$, $\omega=2$, and $c=9$. There are three orbits of
$\langle\omega\rangle$ acting on $\mathbb{Z}/27\mathbb{Z}$:
$\mathcal{O}_{1}:=\mathbb{Z}/27\mathbb{Z}^{\times}$,
$\mathcal{O}_{0}:=\left\\{0\right\\}$, and
$\mathcal{O}_{18}:=\\{[a]\in\mathbb{Z}/27\mathbb{Z}:27>\gcd(a,27)>1\\}$, and
$G(27,2)$ consists of the three points $\eta_{27,2,1},$ $\eta_{27,2,0}$, and
$\eta_{27,2,18}$. So by our rules for coloring, since $c=9$ and $18\equiv
0\,(\operatorname{mod}9)$, we must assign the same colors to $\eta_{27,2,0}$
and $\eta_{27,2,9}$. Since elements in $(\mathbb{Z}/27\mathbb{Z})^{\times}$
are not congruent modulo $9$ to elements divisible by $3$, there are no
restrictions on the color of $\eta_{27,2,1}$. As seen in Figure 2, for this
input, our recipe produces a plot with $3$ points and just $2$ colors, even
though $c=9$.

(a) $n=27$, $\omega=2$, $c=9$

(b) $n=12$, $\omega=5$, $c=3$

(c) $n=12$, $\omega=5$, $c=4$

Figure 2: Small scale illustrations for Examples 1 and 2.

This approach to coloring was chosen not as a way to illustrate a particular
mathematical principle but rather as a recipe Lutz discovered that produced
visually appealing pictures. In fact, in the final paragraph of [7], he and
his coauthors note that this approach to coloring is _ad hoc_ and that “a
general theory is necessary to formalize our intuition.”

## Algebraic structures

Through the lens of Galois theory (which was never considered in prior work on
illustrating $\mathsf{G}(n,\omega)$), the above coloring scheme can be linked
to mathematical structure. The main theorem of Galois theory gives a bijection
between the subfields of $\mathbb{Q}\left(\zeta_{n}\right)$, with $\zeta_{n}$
a primitive $n$th of unity, and the subgroups of the Galois group
$G:=\mathrm{Gal}\left(\mathbb{Q}\left(\zeta_{n}\right)/\mathbb{Q}\right)=\operatorname{Aut}\left(\mathbb{Q}\left(\zeta_{n}\right)\right)$,
which is identified with $\left(\mathbb{Z}/n\mathbb{Z}\right)^{\times}$ via

$\displaystyle\left(\mathbb{Z}/n\mathbb{Z}\right)^{\times}$
$\displaystyle\overset{\sim}{\rightarrow}G$ $\displaystyle[m]$
$\displaystyle\mapsto\psi_{m}$
$\displaystyle\psi_{m}\Big{(}\sum_{j}a_{j}\zeta_{n}^{j}\Big{)}$
$\displaystyle:=\sum_{j}a_{j}\left(\zeta_{n}^{m}\right)^{j},\hskip
36.135pta_{j}\in\mathbb{Q}.$

Each subgroup $H\subseteq G$ corresponds to the fixed field
$\mathbb{Q}\left(\zeta_{n}\right)^{H}:=\left\\{z\in\mathbb{Q}\left(\zeta_{n}\right):\psi(z)=z\mbox{
for all }\psi\in H\right\\}$. Then $H$ is the Galois group of
$\mathbb{Q}\left(\zeta_{n}\right)$ over
$\mathbb{Q}\left(\zeta_{n}\right)^{H}$, i.e., the group of field automorphisms
of $\mathbb{Q}\left(\zeta_{n}\right)$ fixing
$\mathbb{Q}\left(\zeta_{n}\right)^{H}$, and $G/H$ is the Galois group of
$\mathbb{Q}\left(\zeta_{n}\right)^{H}$ over $\mathbb{Q}$, i.e., the group of
field automorphisms of $\mathbb{Q}\left(\zeta_{n}\right)^{H}$.

To help identify elements of the subfields of
$\mathbb{Q}\left(\zeta_{n}\right)$, students in Galois theory courses are
sometimes assigned to determine Gaussian periods (e.g., see [5, §14.5]).
Indeed, if $H=\langle\omega\rangle$, then

$\eta_{n,\omega,k}=\sum_{\psi\in
H}{\left(\psi\left(\zeta_{n}\right)^{k}\right)}=\sum_{\psi\in
H}\psi\big{(}\zeta_{n}^{k}\big{)}\in\mathbb{Q}\left(\zeta_{n}\right)^{H},$

So the plot of $\mathsf{G}(n,\omega)$ contains a (rescaled) plot of
$\mathsf{G}(\frac{n}{\gcd(n,k)},\omega)$. If $c$ is the coloring number, then
the plot of $\mathsf{G}(n,\omega)$ contains a single-colored plot of
$\mathsf{G}(\frac{n}{c},\omega)$, rescaled by a factor of
$\frac{\mathrm{ord}_{n}(\omega)}{\mathrm{ord}_{n/c}(\omega)}$, with
$\mathrm{ord}_{m}(\omega)$ denoting the multiplicative order of
$\omega\,(\operatorname{mod}m)$. For clarification, we illustrate this in a
simple example.

######  Example 2.

Let $n=12$ and $\omega=5$. Note that $5$ has multiplicative order $2$ mod $12$
and $1$ mod $4$. So $\eta_{12,5,3k}=2\eta_{4,5,k}=2e^{\frac{2\pi ik}{4}}$ for
all $k$. So $G(4,5)=\\{1,-1,i,-i\\}$, and $G(12,5)\supseteq 2\cdot
G(4,5)=\\{2,-2,2i,-2i\\}$. If we choose $c=3$, then the points in $2\cdot
G(4,5)$ all must be the same color as each other (the red, outer diamond in
Figure 2), while the other points need not be colored that color (the blue,
inner diamond). On the other hand, selecting $c=4$ forces the pair of points
in $G(3,5)\subset G(12,5)$ to be the same color (shown in Figure 2 in red). As
illustrated in Figure 2, the pair gets rotated by $2\pi/4=2\pi\omega/4$, with
a new color allowed at each rotation.

Regarding the current coloring scheme, if two elements of
$\mathbb{Z}/n\mathbb{Z}^{\times}\cong\mathrm{Gal}(\mathbb{Q}(\zeta_{n})/\mathbb{Q})$
are congruent $\mod c$, then the corresponding elements of the Galois group
restrict to the same element of
$\mathrm{Gal}(\mathbb{Q}(\zeta_{c})/\mathbb{Q})$. So, for example, given an
element $\psi\in\mathrm{Gal}(\mathbb{Q}(\zeta_{c})/\mathbb{Q})$ and an integer
$a$, all points $\Psi(\eta_{n,\omega,a})$, as $\Psi$ ranges over the
extensions of $\psi$ to $\mathbb{Q}(\zeta_{n}),$ are colored the same. There
are also other natural coloring schemes. For example, an option (called
“period squared”) in Gaussian Periods is to color $\eta_{n,\omega,k}$ and
$\eta_{n,\omega,-k}$ the same, thus creating symmetry across the real axis.
This corresponds to coloring the points in the Galois orbit of the complex
conjugation automorphism the same. More generally, one might color all points
in some other given Galois orbits the same. Furthermore, our reformulation in
terms of elements of fixed fields should naturally generalize beyond Gaussian
periods to the illustration of symmetries in other settings beyond the scope
of this short paper (such as finite nonabelian extensions of $\mathbb{Q}$).

## Asymptotic behavior

Symmetry is the most obvious feature in typical Gaussian-period plots. We say
that $\mathsf{G}(n,\omega)$ has _$k$ -fold dihedral symmetry_ if
$\mathsf{G}(n,\omega)$ is invariant under the action of the dihedral group of
order $2k$. That is, $\mathsf{G}(n,\omega)$ is invariant under complex
conjugation and rotation by $2\pi/k$ about the origin. It turns out that
$\mathsf{G}(n,\omega)$ has at least $\gcd(\omega-1,n)$-fold dihedral symmetry
[4, Prop. 3.1]. This symmetry refers to the uncolored graph; the colors
highlight additional features beyond the initial symmetry. This is illustrated
in Figure 3.

(a) $n=29070$, $\omega=1189$, $\gcd(1188,29070)=18$, $c=3$

(b) $n=70091$, $\omega=21792$, $\gcd(21791,70091)=7$, $c=7$

(c) $n=255255$, $\omega=254$, $\gcd(255255,253)=11$, $c=7$

Figure 3: Dihedral symmetry of $\mathsf{G}(n,\omega)$.

Gaussian period plots often demonstrate great structural coherence if the
parameters $n$ and $\omega$ vary in the appropriate manner. The following
“filling out” of various shapes was discovered in [4, Thm. 6.3]. See the
thorough exposition in [7, Thm. 1]. Let $q=p^{a}$ be a nonzero power of an odd
prime and let $\omega=\omega(q)$ be such that $d=|\langle\omega\rangle|$
divides $p-1$. Then $\mathsf{G}(q,\omega)$ is contained in the image of the
Laurent polynomial function $g:\mathbb{T}^{\varphi(d)}\to\mathbb{C}$ defined
by

$g(z_{1},z_{2},\ldots,z_{\varphi(d)})=\sum_{k=0}^{d-1}\prod_{j=0}^{\varphi(d)-1}z_{j+1}^{b_{k,j}},$

where the integers $b_{k,j}$ are determined by
$t^{k}\equiv\sum_{j=0}^{\varphi(d)-1}b_{k,j}t^{j}\,(\operatorname{mod}\Phi_{d}(t)).$
Here $\mathbb{T}$ denotes the unit circle in $\mathbb{C}$, $\varphi$ is the
Euler totient function, and $\Phi_{d}$ denotes the $d$th cyclotomic
polynomial. For a fixed $d$, as $q$ becomes large, $\mathsf{G}(q,\omega)$
“fills out” the image of $g$; see Figure 4.

(a) $n=3019$, $\omega=239$

(b) $n=13063$, $\omega=1347$

(c) $n=9114361$, $\omega=3082638$

Figure 4: Sometimes $\mathsf{G}(n,\omega)$ appears to “fill out” the image of
a Laurent polynomial $g:\mathbb{T}^{\varphi(d)}\to\mathbb{C}$. Here
$g(z_{1},z_{2})=z_{1}+z_{1}+1/(z_{1}z_{2})$.

With ellipses, hypocycloids, and so forth as primitive graphical elements, one
can use the Chinese remainder theorem to produce new images of startling
complexity (as explained in [8]). We content ourselves here with a few
aesthetically pleasing images produced in such a manner; see Figure 5. We
encourage the reader to enjoy more examples by experimenting with our app
Gaussian Periods (freely available at
http://www.elleneischen.com/gaussianperiods.html).

(a) $n=37367$, $\omega=608$, $c=11$

(b) $n=185925$, $\omega=766$, $c=25$

(c) $n=82677$, $\omega=8147$, $c=21$

Figure 5: Primitive graphical elements (here, a “filled” deltoid) can form
more elaborate plots.

## Acknowledgments

Ellen Eischen was partially funded by NSF grant DMS-1751281. Stephan Ramon
Garcia was partially funded by NSF grant DMS-1800123. This material is based
partly upon work supported by the NSF grant DMS-1439786 and the Alfred P.
Sloan Foundation award G-2019-11406 while the authors were in residence at the
Institute for Computational and Experimental Research in Mathematics (ICERM)
in Providence, RI, during the Illustrating Mathematics program (Fall 2019). We
thank R. Lipshitz for help with the code.

## References

  * [1] J. Baez. The beauty of roots. 2011\. http://math.ucr.edu/home/baez/roots/.
  * [2] T. M. Brown, “Kaleidoscopes, Chessboards, and Symmetry,” Journal of Humanistic Mathematics, vol. 6, no. 1, 2016, pp. 110–126. https://scholarship.claremont.edu/jhm/vol6/iss1/8
  * [3] J. L. Brumbaugh, M. Bulkow, P. S. Fleming, L. A. G. German, S. R. Garcia, G. Karaali, M. Michal, A. P. Turner, and H. Suh. “Supercharacters, exponential sums, and the uncertainty principle.” Journal of Number Theory, vol. 144, 2014, pp. 151–175.
  * [4] W. Duke, S. R. Garcia, and B. Lutz. “The graphic nature of Gaussian periods.” Proceedings of the American Mathematical Society, vol. 143, no. 5, 2015, pp. 1849–1863.
  * [5] D. S. Dummit and R. M. Foote. Abstract Algebra, 3rd ed., John Wiley & Sons, 2004.
  * [6] E. E. Eischen and R. Lipshitz. Gaussian Periods. http://www.elleneischen.com/gaussianperiods.html, 2020.
  * [7] S. R. Garcia, T. Hyde, and B. Lutz. “Gauss’s hidden menagerie: from cyclotomy to supercharacters.” Notices of the American Mathematical Society, vol. 62, no. 8, 2015, pp. 878–888.
  * [8] B. Lutz. “Graphical cyclic supercharacters for composite moduli.” Proceedings of the American Mathematical Society, vol. 147, no. 9, 2019, pp. 3649–3663.
  * [9] I. Stewart. Why Beauty is Truth: The History of Symmetry. Basic Books, 2007.