Introduction

The bilayer graphene with van der Waals interlayer interaction shows rich electronic properties that are depending on the stacking order and twist angle.1,2,3,4,5,6,7,8,9,10,11,12,13,14,15 Two graphene layers can be arranged in AA, AB or twisted configurations. AA and AB stacking are two well-known configurations with electronic properties drastically different from the monolayer.16,17,18,19 The twist angle between the layers offers an additional degree of freedom to tune the electronic properties. If the commensuration condition is satisfied,20 namely, the twist angle is \(\theta =co{s}^{-1}(\frac{3{q}^{2}\;-\;{p}^{2}}{3{q}^{2}\;+\;{p}^{2}})\), where \(p\) and \(q\) are integers, the twisted bilayer graphene (tBG) forms Moiré pattern, which is periodic in the space and has an elementary unit cell. The incommensurate tBG, on the other hand, has only quasi-periodicity without any translational symmetry.

By now, the tBG with 30\({}^{\circ }\) twist angle, i.e. dodecagonal graphene quasicrystal, which will be referred as graphene quasicrystal for simplicity in the rest of the paper, has been grown successfully on different substrates, including SiC(0001),21,22 Pt(111)23 and Cu-Ni(111)24 surfaces. This emergent quasicrystal consisting of two graphene sheets with perfect crystalline has attracted increasing attention because of the coexistence of the quasicrystalline nature and the relativistic properties.21,22,23,24,25,26,27,28,29,30 The 12-fold rotational symmetry has been determined by the Raman spectroscopy, low-energy electron microscopy/diffraction (LEEM/LEED), transmission electron microscopy (TEM) and scanning tunneling microscopy (STM) measurements.21,22,23,24,25 Angle-resolved photoemission spectroscopy (ARPES) measurements indicated that the interlayer interaction between the two layers leads to the emergence of the mirror-symmetric Dirac cones inside the Brillouin zone of each graphene layer21,22,23 and a gap opening at the zone boundary.23 The critical eigenstates27,29 and quantum oscillations28 were observed theoretically due to the quasi-periodicity. These experimental observations indicate that graphene quasicrystal is very different from graphene, although several studies proposed that the tBG with large twist angle (\(>\)15\({}^{\circ }\)) should behave like two decoupled graphene monolayers.31,32,33,34,35,36,37,38

These electronic properties behind the quasi-periodicity of graphene quasicrystal make it very urgent to propose a method to study its electronic properties further in theory because the band theory can not be applied directly. Recently, the k-space tight-binding method29 and continuum model39 have been proposed to understand its electronic structure. Both two methods have to introduce the k-wavevector cutoff to implement the calculations because k-wavector is no longer a good quantum number due to the interlayer interaction. However, a number of band theory based methods and programs can not be applied to post-process the data obtained from these two methods. Besides, the calculated quasi-band structure from k-space tight-binding method is supercell-like, and can not be compared directly with the ARPES measurements. In order to overcome these diffculties, we propose a series of approximants with transitional symmetry with the help of the large-scale tight-binding calculations in real space. The effective band structure, which can be compared with ARPES measurements directly, can be derived by unfolding the band structures of these approximants. For a quasicrystal, its approximant is a periodic structure that contains similar compositions and almost the same local atomic structure within a unit cell. For graphene quasicrystal, although the structure of an approximant can be different from the original quasicrystal, a good approximant should form the same rotational order within a unit cell and keep similar physical properties as the infinite quasicrystal. The approximant method has been widely used to understand the physics of a quasicrystal (for a review see ref. 40 and references therein). Importantly, theoretical studies on a periodic approximant are much more easier comparing to a real quasicrystal, since methods established for crystalline in condensed matter can be readily used. The main purpose in this paper is to construct the reliable and convincing approximants for graphene quasicrystal, and systematically study its electronic properties, including optical conductivity, effective band structure and Landau level spectrum in the presence of external magnetic fields.

To construct an accurate approximant, one needs first calculate properly the characteristics of a real quasicrystal without any further approximation. This is indeed a quite challenging problem from a computational point of view, as a quasicrystal is not periodic in the space and it contains, in principle, infinite number of sites. In this paper, graphene quasicrystal is modeled by a round disk of exactly 30\({}^{\circ }\) tBG. In order to reproduce the bulk properties of an infinite graphene quasicrystal, such as density of states and optical conductivities, a large enough round disk is needed to get rid of the influence of the edge states. The crucial issue is therefore how to efficiently calculate these properties of the sample with a large radius. It is for sure beyond the commonly used density-functional theory but may be accessible within the tight-binding model. In fact, for twisted bilayer graphene and transitional metal dichalcogenides, there are accurate tight-binding models developed from first principle calculations.41,42,43 The hopping parameters in these models have both distance and orientation dependence and can be used to describe properly interlayer interactions in quasi-periodic bilayer. We will use the one developed in ref. 43 which has been verified by comparing results with several experiments.43,44,45 Recent theorical study proposed that the electronic structure of tBG can be understood as Dirac fermions coupled with opposite pseudo-magnetic fields generated by the moiré pattern.46 The periodic intrinsic pseudo-magnetic fields in tBG have been observed in experiment and can be reproduced by this tight-binding model.45

For a large round disk of graphene quasicrystal described by the tight-binding model, such as a system as large as ten million atoms, the calculations of the electronic properties is still challenging. A common approach one may consider is to diagonalize the Hamiltonian. However, in the diagonalization processes the costs of memory and CPU time scale as \(O({N}^{2})\) and \(O({N}^{3})\), respectively. In view of the potentially large number of arithmetic calculations, it is advisable to use 1315 digit floating-point arithmetic which corresponds to 8 bytes for a real number and 16 bytes for a complex number. Thus, for a tight-binding sample with \(N\) sites, we need in total \(16\;\times {N}^{2}\) bytes memory for the storage of all the eigenstates. Considering a common computational node with 256 GB memory, it has a maximum storage for eigenstates of about 126,491 sites, which is indeed too small to represent properly a quasicrystal. In fact, our numerical tests show that even a quasicrystal sample with two million sites is not large enough to neglect the influence of the edge states (please see section SI in supplementary information). One may use large supercomputers with more nodes and memory to increase the number of sites that can be reached, but it is computational too expensive as the method based on diagonalization does not scale linearly with the size of the sample. We will apply a different approach, tight-binding propagation method (TBPM),47 to overcome the difficulties raised by diagonalization. TBPM is based on the numerical solution of time-dependent Schrödinger equation without any diagonalization, and importantly, both memory and CPU costs scale linearly with the system size. The calculations of electronic, optical, transport and plasmonic properties can be easily implemented in TBPM without the requirement of any symmetry. In this paper, our main purpose is to build approximants and reproduce the electronic properties that have been observed in the experiments. More studies of other properties, which are not presented here, can be further explored by TBPM or applying band theory on these small approximants proposed.

Results and discussion

In this paper, graphene quasicrystal is simulated by the tight-binding model based on \({p}_{z}\) orbitals, where the hopping energy between site \(i\) and \(j\)48

$${t}_{ij}={n}^{2}{V}_{pp\sigma }(\left|{{\boldsymbol{r}}}_{{\boldsymbol{ij}}}\right|)+(1-{n}^{2}){V}_{pp\pi }(\left|{{\boldsymbol{r}}}_{{\boldsymbol{ij}}}\right|).$$
(1)

Here, \(n\) is the direction cosine of relative position vector \({{\boldsymbol{r}}}_{{\boldsymbol{ij}}}\) with respect to \({{\boldsymbol{e}}}_{{\boldsymbol{z}}}\). The Slater and Koster parameters \({V}_{pp\sigma }\) and \({V}_{pp\pi }\) have the following form:

$${V}_{pp\pi }(\left|{{\boldsymbol{r}}}_{{\boldsymbol{ij}}}\right|)=-{\gamma }_{0}{e}^{2.218(b-\left|{{\boldsymbol{r}}}_{{\boldsymbol{ij}}}\right|)}{F}_{c}(\left|{{\boldsymbol{r}}}_{{\boldsymbol{ij}}}\right|),$$
(2)
$${V}_{pp\sigma }(\left|{{\boldsymbol{r}}}_{{\boldsymbol{ij}}}\right|)={\gamma }_{1}{e}^{2.218(h-\left|{{\boldsymbol{r}}}_{{\boldsymbol{ij}}}\right|)}{F}_{c}(\left|{{\boldsymbol{r}}}_{{\boldsymbol{ij}}}\right|).$$
(3)

The interlayer distance \(h\) and nearest carbon-carbon distance \(b\) are chosen to be 3.349 and 1.418 Å respectively. \({\gamma }_{0}\) and \({\gamma }_{1}\) are turned to 3.12 and 0.48 eV to fit the experimental Fermi velocity (8.9 \(\sim\) 10.0\(\times 1{0}^{5}\) m/s), respectively.21,23 \({F}_{c}\) is a smooth function

$${F}_{c}(r)={(1+{e}^{(r-0.265)/5})}^{-1}.$$
(4)

Under a magnetic field, the hopping \({t}_{ij}\) will be replaced by a Peierls substitution.47

Approximant construction

The structure of graphene quasicrystal is shown in Fig. 1a. Along x direction, bottom layer (black) and top layer (red) have the period \(3b\) and \(a\) (\(=\!\!\sqrt{3}b\)), respectively. If a common period exists in such a structure, there should be two integers \(M\) and \(N\) that satisfy \(M\;\times 3b=N\times a\), i.e., \(N/M=\sqrt{3}\). Graphene quasicrystal posses the quasi-periodicity because \(\sqrt{3}\) is an irrational number, namely \(N/M=\sqrt{3}\) is never satisfied. The commensurate configurations of tBG with twist angle close to 30\({}^{\circ }\) can be used as the approximant, but the 12-fold rotational order will be destroyed. In this paper, we construct the approximant in the following way: the twist angle of 30\({}^{\circ }\) is fixed and the top graphene layer is compressed or stretched to satisfy the condition \(M\,\times\, 3b=N\,\times\, {a}_{t}\), where \({a}_{t}\) is the lattice constant of the top graphene layer with strain. This method was applied to construct the periodic structure to calculate the formation energy of graphene quasicrystal by using first principle calculations.23

Fig. 1: The atomic structures.
figure 1

a Graphene quasicrystal. b 4/7 approximant with four elementary unit cells. c Lattice constant \({a}_{t}\) of the top layer of the approximants with \(M\,<\,100\). The horizontal dashed line shows the lattice constant \(a=2.456\) Å of pristine graphene.

Full size image

The procedure to construct the approximant of graphene quasicrystal is given below: For a specific integer \(M\) (in this paper, all integers less than 100 are considered), \(N\) can be determined by

$$N=INT\left(\frac{M\,\times\, 3b}{a}\right)$$
(5)

where, \(INT(x)\) stands for the integer closest to \(x\). Then, \({a}_{t}\) can be obtained by

$${a}_{t}=\frac{M\,\times 3b}{N}.$$
(6)

We name such an approximant as \(M/N\), which has the lattice vectors \({{\boldsymbol{a}}}_{1}=(M\times 3b,0)\) and \({{\boldsymbol{a}}}_{2}=(\frac{1}{2}M\times 3b,\frac{\sqrt{3}}{2}M\times 3b)\). The structure of \(4/7\) approximant is shown in Fig. 1b. Note that if \(M\) and \(N\) share a common divisor, or \(N\) can be divisible by 3, a smaller elementary unit cell exists. Then, the approximants with such \(M\) and \(N\) will not be considered. The curve of the lattice constant \({a}_{t}\) with respective to the approximant size is given in Fig. 1c. It indicates that \({a}_{t}\) converges to \(a\) with the increase of the approximant size. Usually, the approximant with larger size is expected to reproduce more accurately the electronic properties of the quasicrystal. However, this is not always the case as we will show in the following, and it is important to figure out what are dominant factors that determine the accuracy of these approximants shown in Fig. 1c.

Density of states

The density of states (DOS) of graphene quasicrystal is shown in Fig. 2a. It is calculated from a round disk with radius about 204 nm and more than ten million atoms by using the TBPM.47 The previous calculation by using the k-space tight-binding method29 showed spiky peaks and dips in DOS, and all these Van Hove singularities they observed can be reproduced in our calculations. As the tight-binding parameter \({\gamma }_{0}\) we adopt (3.12 eV) is different from the value 2.7 eV used in ref., 29 the positions of these peaks in DOS are not exactly the same. Previous study23 and our calculations indicate that the experimental Fermi velocity \(9.1 \sim 10\times 1{0}^{5}\) m/s can be reproduced correctly when \({\gamma }_{0}\) = 3.1 \(\sim\) 3.12 eV. As shown in Fig. 2a, comparing the DOS of graphene quasicrystal with graphene, the significant difference is the emergence of some peaks in the spectrum, which is attributed to the interaction between layers.29 In the vicinity of Fermi level, the DOS of graphene quasicrystal is almost the same as pristine graphene, which indicates that their electronic and optical properties may behave similar at low energies.

Fig. 2: The comparison of DOS obtained from graphene quasicrystal (QC) and its approximants.
figure 2

For graphene quasicrystal, a round disk with more than ten million atoms is used. For approximants, the numbers of carbon atoms in a unit cell are given in brackets. a Accurate approximants with less than 60,000 carbon atoms in a unit cell, and the smallest one (15/26) contains only 2702 atoms. DOS of pristine graphene is also given for reference. b Examples of several inaccurate approximants. c The standard deviations of DOS of some approximants at different lattice mismatches, which are fitted by the function \(\Delta \rho =0.157\times {\delta }^{0.562}\).

Full size image

In order to figure out the influence of lattice constant (\({a}_{t}\)) and approximant size (\(M\)) on the accuracy of the approximants, we compare the DOS of graphene quasicrystal with these obtained from approximants in Fig. 2a, b. The results indicate that as long as \({a}_{t}\) is close enough to \(a\), even a small approximant, such as 15/26 which contains only 2702 carbon atoms, can reproduce the DOS of graphene quasicrystal with very high accuracy that one can not distinguish them by eye (see the results plotted in Fig. 2a). However, if \({a}_{t}\) is far from \(a\), even the sample contains much more atoms, and it can not be supposed to be a qualified approximant. For example, as the data shown in Fig. 2b, 80/139 contains 77042 carbon atoms in the unit cell, but the calculated DOS differs from the one obtained from quasicrystal, especially around Van Hove singularities.

The quantitative measurement on the accuracy of an approximant can be described by the standard deviation

$$\Delta \rho =\sqrt{\frac{1}{N}\sum _{i}{\left|{\rho }^{QC}({\epsilon }_{i})-\rho ({\epsilon }_{i})\right|}^{2}}$$
(7)

where \(N\) is the number of energy points, \({\rho }^{QC}({\epsilon }_{i})\) and \(\rho ({\epsilon }_{i})\) are the DOS of graphene quasicrystal and the approximant at energy \({\epsilon }_{i}\), respectively. We plot \(\Delta \rho\) of several approximants as a function of the relative lattice mismatch \(\delta\) (\(=\!\left|a-{a}_{t}\right|/a\)) in Fig. 2c, which show clearly that \(\Delta \rho\) increases monotonically with increasing lattice mismatch. It is obvious that the smaller lattice mismatch is, the more accurate DOS the approximant can reproduce. More precisely, the data can be fitted with a simple function \(\Delta \rho =0.157\times {\delta }^{0.562}\). Our numerical tests indicate that the lattice mismatch is a crucial factor to justify the accuracy of the approximant. In the following, we will only focus on these accurate approximants shown in Fig. 2a, such as 15/26, 41/71, 67/116 and so on. Actually, when a strain is applied to the top layer, both two factors, namely the change of the electronic structure of the top layer and the change of the interlayer interaction between the two layers, contribute to the deviation of DOS from the graphene quasicrystal. Our results (see section SII in Supplementary information) show that at small lattice mismatch the change of the electronic structure of the top layer is the dominating factor. But the change of the interlayer interaction becomes more important as the lattice mismatch increases and both two factors are comparable at large lattice mismatch.

We also consider the influence of the structure relaxation on DOS (see section SIII in Supplementary information). For one limiting case that the atoms in the bottom layer are fixed, which corresponds to the strongest control from the substrate, the top layer almost keep flat, and the DOS is almost the same as the non-relaxed structure. For the other limiting case, namely a free-standing sample, which corresponds to the weakest control from the substrate, the sample undergoes a strong lattice deformation around the edge and there are also some small ripples around the center. With such lattice deformation the change of the DOS of the graphene quasicrystal becomes more obvious, but all the Van Hove singularities appeared in non-relaxed structure still exist and their energies don’t change. As the deformation of the lattice structure may lead to the emergence of both scale and vector potentials, and subsequently peudo-magnetic fields in some parts of the sample,45,46,49,50 it is worth to study the relaxation of graphene quasicrystal in specified conditions (such as substrate and boundary), and discuss them in details as a separated paper.

Optical conductivity

Density of states is the counting of states at certain energy, and it does not contain the details of the wave function, such as the amplitude and phase distribution in the space. To make the verification of our approximants more complete, the comparison of optical conductivities among graphene quasicrystal, its approximants and graphene is shown in Fig. 3. The calculation of optical conductivity (see the “Methods” section for details) is based on tight-binding propagation method47 without diagonalization of the Hamiltonian matrix. Although the method itself does not using directly eigenstates of the Hamiltonian, it is equivalent to the standard Kubo formulation which calculates excitations between occupied and unoccupied eigenstates. Optical conductivity is indeed a bulk property determined by wave functions, and can be used to check the accuracy of our approximants. Indeed, as we see from the results plotted in Fig. 3, there is a perfect agreement between graphene quasicrystal and its approximants with small mismatch as these shown in Fig. 2a. Importantly, by comparing with the results of pristine graphene, there are emerged peaks around 4.0 \(\sim\) 4.6 eV in the optical spectrum of graphene quasicrystal, which are reproduced exactly with the same energies and amplitudes by proposed approximants. These peaks in the optical spectrum are attributed to Van Hove singularities appeared in density of states shown in Fig. 2a. It is possible to identify the one-to-one correspondence of the peaks in optical spectrum and density of states by using the band structure of proposed approximant. We leave these detailed studies in future, together with other optical and plasmonic properties of graphene quasicrystals. The peak at 5.63 eV corresponds to the transition between the singularities of graphene monolayer, which shifts towards higher frequency by about 0.13 eV than graphene.

Fig. 3: The comparison of optical conductivities of graphene quasicrystal (QC), its approximants and monolayer graphene.
figure 3

The inset shows a zoom of the peaks associated with the quasi-periodic states appeared in graphene quasicrystal, which is all accurately reproduced by proposed approximants.

Full size image

Eigenstates with 12-fold rotational symmetry

In spite of the periodicity of approximants, the quasi-periodicity still remains inside each unit cell. The electronic properties related to the quasi-periodicity can be verified by the extistance of the 12-fold rotational symmetric eigenstates within the unit cell. In this paper, 12-fold rotational symmetry means a rotation of \(6{0}^{\circ }n+3{0}^{\circ }\) followed by the mirror reflection with the mirror plane in the middle of the two graphene layers plus just rotation of \(6{0}^{\circ }n\) (\(n\) is any integer), because after these operations the structure remains unchanged. Two 12-fold rotational symmetric eigenstates of 15/26 approximant are shown in Fig. 4a, b, and the corresponding states in 41/71 approximant are shown in Fig. 4c, d, respectively. Comparing with 15/26 approximant, some new 12-fold rotational symmetric eigenstates appear for the 41/71 approximant, two of which are given in Fig. 4e, f. Such a result is reasonable. For a real graphene quasicrystal, some critical eigenstates expand more than 20 nm in space,29 which can not be simulated by a small approximant, for instance, the approximant with 29.84\({}^{\circ }\) twist angle in ref. 29 and the 15/26 approximant in our model. But as the approximant size increased to be large enough, these critical states appear again. It is worth noting that the top layer still has the tiny strain in 15/26 and 41/71 approximants, so the eigenstates we discussed here have approximately but not exactly the 12-fold rotational symmetry.

Fig. 4: Eigenstates with 12-fold rotational symmetry.
figure 4

a, b The eigenstates of 15/26 approximant at −4.2 and −2.76 eV, respectively. c–f The eigenstates of 41/71 approximant at −4.2, −2.76, −2.4, and −2.23 eV, respectively. Red and blue circles represent the projection on the top and bottom layers, respectively.

Full size image

Effective band structure

According to the generalized Umklapp scattering theory,51,52 the interlayer interaction will couple the wavevector \({\boldsymbol{k}}\) in the bottom layer and the wavevector \(\widetilde{{\boldsymbol{k}}}\) in the top layer if \({\boldsymbol{k}}+{\boldsymbol{G}}=\widetilde{{\boldsymbol{k}}}+\widetilde{{\boldsymbol{G}}}\) is satisfied, where \({\boldsymbol{G}}\) and \(\widetilde{{\boldsymbol{G}}}\) are the reciprocal lattice vectors of bottom and top layers respectively. The coupling matrix element \(\left|\left\langle \widetilde{{\boldsymbol{k}}},\widetilde{X}| U| {\boldsymbol{k}},X\right\rangle \right|\) can also be understood as the scattering from the \({\boldsymbol{k}}\) state of bottom layer to \(\widetilde{{\boldsymbol{k}}}\) state of top layer, where \(\left|{\boldsymbol{k}},X\right.\rangle\) (\(\left|\widetilde{{\boldsymbol{k}}},\widetilde{X}\right.\rangle\)) is the Bloch basis function in sublattice \(X\) (\(\widetilde{X}\)) for bottom (top) layer, \(U\) is the interlayer interaction. The scattering strength just depends on the length of \(\left|{\boldsymbol{q}}\right|=\left|{\boldsymbol{k}}+{\boldsymbol{G}}\right|=\left|\widetilde{{\boldsymbol{k}}}+\widetilde{{\boldsymbol{G}}}\right|\), because the hopping energy between two layers is isotropic along the in-plane direction. The smaller \(| {\boldsymbol{q}}|\) corresponds to the stronger scattering process. The Brillouin zones (BZs) of two graphene layers are shown in Fig. 5a. If we take \(K\), \({K}^{^{\prime} }\), \(\widetilde{K}\), and \({\widetilde{K}}^{^{\prime} }\) as the original points for scattering, where the Dirac cones exist for two layers, they will be scattered to their mirrored points \({\widetilde{K}}_{1}\) and \({K}_{1}\) after the strongest scatterings and \({\widetilde{K}}_{2}\) and \({K}_{2}\) after the second strongest scatterings, respectively. As an example, the first two strongest scattering paths from \({K}^{^{\prime} }\) are shown in Fig. 5b, c, respectively.

Fig. 5: The effective band structure of graphene quasicrystal.
figure 5

a Brillouin zones of the two layers. b, c The first two strongest scattering paths for \({K}^{^{\prime} }\) (\({K}^{^{\prime} }\to {\widetilde{K}}_{1}\) and \({K}^{^{\prime} }\to {\widetilde{K}}_{2}\)). d The effective band structure of graphene quasicrystal obtained by unfolding the band structure of 15/26 approximant. The insets at the bottom left and right corners show the gaps at \(M\) and \(\widetilde{M}\) clearly. e Fermi surfaces around \(K\), \({K}_{1}\), \({K}_{2}\), \(\widetilde{K}\), \({\widetilde{K}}_{1}\), and \({\widetilde{K}}_{2}\) above Fermi energy by 35 meV. f Effective band structures around \(K\), \({K}_{1}\), \({K}_{2}\), \(\widetilde{K}\), \({\widetilde{K}}_{1}\), and \({\widetilde{K}}_{2}\) along four directions. The larger spectral function is denoted by larger black dot in (d, f) and lighter color in (e), respectively. The spectral functions around \({K}_{1}\)/\({\widetilde{K}}_{1}\) and \({K}_{2}\)/\({\widetilde{K}}_{2}\) are multiplied by 300 and \(5\times 1{0}^{6}\), respectively.

Full size image

Previous ARPES measurements21,23 and theoretical studies21 show that the band structures around these six kpoints (\(K\), \(\widetilde{K}\), \({K}_{1}\), \({\widetilde{K}}_{1}\), \({K}_{2}\), and \({\widetilde{K}}_{2}\)) are Dirac cones. The Fermi velocities for \(\widetilde{K}\), \(K\), \({K}_{1}\), and \({\widetilde{K}}_{1}\) Dirac cones are 9.3, 9.2, 8.9, and 9.1 \(\times 1{0}^{5}\) m/s, respectively.21 Due to the hybridization between the Dirac cones at \(K\) and \({\widetilde{K}}_{1}\), a \(\sim\)200 meV gap can be observed at \(M\) point below the Fermi level for the graphene quasicrystal on Pt(111) substrate.23

Now, we check whether the approximants we proposed can reproduce the experimental results. Because the approximant contains lots of unit cells of the two layers, the band structure calculated directly from the approximant can not be used to compare with the ARPES measurements. The effective band structure (EBS) derived by applying the band-unfolding method can overcome this problem. Here we focus on only the smallest approximant (15/27). The EBS along the path plotted by blue dashed line in Fig. 5a is given in Fig. 5d, which contains about four bands in the whole energy region. That is because only two \({p}_{z}\) orbitals exist in the unit cell of each monolayer. The strength of the spectral function becomes weaker at the Fermi level from \(K\) (\(\widetilde{K}\)) to \({\widetilde{K}}_{2}\)(\({K}_{2}\)) via \({\widetilde{K}}_{1}\) (\({K}_{1}\)). It just follows the order of scattering strength, namely, the weaker scattering will lead to weaker strength of the spectral function of the Dirac cones at the ending points. Such results are in accordance with the generalized Umklapp scattering theroy21 and the ARPES measurements.23 In order to show the results clear, the spectral functions around \({K}_{1}\)(\({\widetilde{K}}_{1}\)) and \({K}_{2}\) (\({\widetilde{K}}_{2}\)) are always timed by 300 and 5\(\times 1{0}^{6}\), respectively in Fig. 5. The Fermi surfaces and the detailed EBS around the six kpoints are shown in Fig. 5e, f, respectively, which show the band structures are all Dirac cones clearly. Although our results indicate that the four Dirac cones at \(\widetilde{K}\), \(K\), \({\widetilde{K}}_{1}\) and \({K}_{1}\) have the same Fermi velocity of \(9.08\times 1{0}^{5}\) m/s, the value is still very close to that observed in experiment. Moreover, the spectral function shows strong anisotropic intensity, which leads to the strong contrast of two bands along y direction for \({\widetilde{K}}_{1}\) and \({\widetilde{K}}_{2}\) and along x direction for \({K}_{1}\) and \({K}_{2}\), and almost half circles of the Fermi surfaces at \({K}_{1}\), \({\widetilde{K}}_{1}\), \({K}_{2}\) and \({\widetilde{K}}_{2}\). As shown in the insets in Fig. 5d, the band gap at \(M\) point below Fermi level, which has been observed in experiment,23 can also be obtained by our model, although the gap value (\(\sim\)130 meV) is a little smaller than that in experiment (\(\sim\)200 meV).

It is worth noting that our calculations show almost equivalent results for the two graphene layers because of their almost the same lattice constants. However, for the graphene quasicrystal grown on Pt(111) surface, mirrored Dirac points of only top layer can be detected.23 But for the one grown on 4H-SiC(0001) surface, the mirrored Dirac points of both layers can be detected, although the ARPES signals are obviously different.21 The signal is even stronger for some mirrored Dirac cones than for the original ones.21 Besides, graphene quasicrystal on 4H-SiC(0001) should be n-type doped according to the ARPES measurements.21 Therefore, both Pt(111) and 4H-SiC(0001) substrates should impact the electronic properties of graphene quasicrystal. The existence of substrates may be the reason why our results are not exactly the same as the experimental results. However, our theoretical results from approximant 15/27 are in accordance with generalized Umklapp scattering theory21 and main experimental results.21

Landau levels

For graphene under a strong magnetic field perpendicular to the graphene plane, there is a quantization of the energy states, the so-called Landau levels.53,54,55,56 The two-dimensional Dirac fermion behavior of graphene can be expressed by the Landau levels which satisfies \({E}_{n}={v}_{F}\sqrt{2e\hslash B\left|N\right|}\). In order to study the electronic properties of graphene quasicrystal under magnetic field, we calculate the Hofstadter’s butterfly using 41/71 approximant under the magnetic field less than 50 T (see Fig. 6a). Similar to graphene, the DOS of the graphene quasicrystal under magnetic field also shows the Landau levels. In the vicinity of the Fermi energy, they follow the two-dimensional Dirac fermion behavior but with the Fermi velocity \({v}_{F}\) of 9.20 and 8.84\(\times 1{0}^{5}\) m/s for hole and electron, respectively, which are obtained by fitting Landau levels under 50 T magnetic field (see Fig. 6b). It agrees well with the values from the EBS (9.1\(\times 1{0}^{5}\) m/s) and experiments (8.9 \(\sim\) 9.3\(\times 1{0}^{5}\) m/s).21 Moreover, it is almost the same as but slightly smaller than that in graphene (9.25 and 9.05 \(\times 1{0}^{5}\) m/s for hole and electron, respectively, which are calculated by using the same tight-binding parameters as one layer in graphene quasicrystal). It means the electronic properties of graphene quasicrystal should be similar to graphene at low energies. More interestingly, some new Landau levels appear below Fermi level by about 1.6 eV when magnetic field is more than 10 T. By fitting the new Landau levels, it can be found that they also follow the two-dimensional Dirac fermion behavior but with a reduced Fermi velocity 5.21 \(\times 1{0}^{5}\) m/s. Besides, the Landau level of n = 0 doesn’t exist, but its position is predicted to be around 1.49 eV below Fermi level by interpolation. It is also around the position where the band gap at \(M\) appears, and the \(K\)(\(\widetilde{K}\)) and \({\widetilde{K}}_{1}\)(\({K}_{1}\)) valleys hybridize strongest.

Fig. 6: The Landau levels.
figure 6

a Hofstadter’s butterflies of 41/71 approximant with magnetic field less than 50 T. Upper and lower panels correspond to the energy regions of −1.0 \(\sim\) 1.0 eV and −2.0 \(\sim\) −1.5 eV, respectively. Color bar stands for the value of DOS. The blue numbers in the lower panel indicate the indexes of the corresponding Landau levels. b The DOS of 41/71 approximant (left) and Landau levels fitting (right) under a 50 T magnetic field. The inset is the zoom in of the new emerging Landau levels in the energy range about −2.0 \(\sim\) −1.6 eV. For the Landau levels fitting, the top and middle panels correspond to the Landau levels of the holes and electrons, respectively, in the vicinity of the Fermi energy, and the bottom panel corresponds to the Landau levels in the inset.

Full size image

In summary, we performed a systematic study of the electronic properties of graphene quasicrystal in the framework of tight-binding approximation. Large-scale calculations of round disks of graphene quasicyrstal with more than ten million atoms have been implemented to model the real graphene quasicrystal. Furthermore, we proposed a series of approximants with translational symmetry to represent graphene quasicrystal, and the accuracies of these approximants have been verified by comparing their density of states and optical conductivities to the large round disk of graphene quasicrystal. The number of atoms in these approximants are only few thousands or tens of thousands. The lattice mismatch between two layers in the approximant is found to be the dominant factor, which determines the accuracy of this approximation. An approximant with smaller mismatch can approximate graphene quasicrystal better than the one with larger mismatch, independent on the size of the approximant. This is indeed a quite surprising result as one would expect that larger approximant would leads to better accuracy. In fact, decreasing lattice mismatch between layers should be a designing principle when building approximants for any incommensurate layered structure.

Furthermore, by applying band-unfolding procedure to the smallest approximant with 2702 atoms, the effective band structure of graphene quasicrystal can be derived and compared directly with recent ARPES measurements. Such a comparison indicates that its properties agree well with the main experimental results, such as: (1) the emergence of new Dirac points, especially the mirrored ones, (2) the appearance of a band gap \(\sim\)130 meV (\(\sim\)200 meV in experiment) at the \(M\) points below the Fermi energy, and (3) the Fermi velocity of \(\sim \!9.1\times 1{0}^{5}\) m/s. Besides, our results show a strong anisotropic intensity in spectral function around these new Dirac points. By calculating the density of states in the presence of strong perpendicular magnetic fields (B \(>\) 10 T), we found that besides the usual Landau level spectrum in the vicinity of Fermi energy, which is almost the same as monolayer graphene, a group of new Landau levels appear with energies about 1.6 eV below the Fermi energy. They also follow the property of two-dimensional Dirac fermion \({E}_{n}={v}_{f}\sqrt{2e\hslash nB}\) but with a reduced Fermi velocity of \(\sim \!5.1\times 1{0}^{5}\) m/s. Interestingly, the zero order Landau level does not appear in the spectrum, no matter how strong the magnetic field is. The optical conductivities of graphene quasicrystal and its approximants have been studied numerically by using the tight-binding propagation method within the linear response theory. The optical excitations associated with quasicrystal states have been observed between 4.0 \(\sim\) 4.6 eV, which should be measurable in optical experiments.

Importantly, our proposal of the approximants for graphene quasicrystal can be applied directly for any bilayer or multilayer quasicrystals formed by two-dimensional honeycomb lattices. It works not only for graphene, but also for other materials such as hexagonal transition metal dichalcogenides, e.g., \(Mo{S}_{2}\), \(W{S}_{2}\), \(MoS{e}_{2}\) and \(WS{e}_{2}\), etc. In fact, to model accurately these materials, there are much more orbitals and hoppings that need to be considered in the tight-binding model,41,42 therefore, an approximant with finite unit cell will dramatically simplify the complexities of the modeling. The numerical costs such as memory and CPU time in the numerical calculations will be reduced significantly. The studies of layered quasicrystals with other rotational order and/or other materials by using the principles proposed in this paper will be continued in future works.

Methods

Effective band structure

First of all, the spectral function at wavevector \({\boldsymbol{k}}\) and energy \(\epsilon\) can be calculated by57

$$A({\boldsymbol{k}},\epsilon)=\sum _{I{{\boldsymbol{k}}}_{SC}}{P}_{I{{\boldsymbol{k}}}_{SC}}({\boldsymbol{k}})\delta (\epsilon -{\epsilon }_{I{{\boldsymbol{k}}}_{SC}}),$$
(8)

where \({\epsilon }_{I{{\boldsymbol{k}}}_{SC}}\) is the energy for \({I}^{\mathrm{th}}\) band at wavevector \({{\boldsymbol{k}}}_{SC}\) for the approximant. Actually, only one \({{\boldsymbol{k}}}_{SC}\), namely \({{\boldsymbol{k}}}_{SC}={\boldsymbol{k}}+{\boldsymbol{G}}\) being \({\boldsymbol{G}}\) the reciprocal lattice vector of the approximant, contributes to the spectral function. The spectral weight is defined by

$${P}_{I{{\boldsymbol{k}}}_{SC}}({\boldsymbol{k}})=\sum _{s=1}^{2}\sum _{i}{\left|\left\langle {\psi }_{i{\boldsymbol{k}}}^{P{C}_{s}}| {\Psi }_{I{{\boldsymbol{k}}}_{SC}}^{SC}\right\rangle \right|}^{2}=\sum _{s=1}^{2}{P}_{I{{\boldsymbol{k}}}_{SC}}^{s}({\boldsymbol{k}}),$$
(9)

where \(\left|{\psi }_{i{\boldsymbol{k}}}^{P{C}_{s}}\right.\rangle\) and \(\left|{\Psi }_{I{{\boldsymbol{k}}}_{SC}}^{SC}\right.\rangle\) are the eigenstates of layer \(s\) and the approximant, respectively. Under the tight-binding method, the spectral weight contributed from layer \(s\) can be described by

$${P}_{I{{\boldsymbol{k}}}_{SC}}^{s}({\boldsymbol{k}})=\frac{1}{{n}_{s}}\sum _{\alpha }\sum _{{{\boldsymbol{l}}}_{{\boldsymbol{s}}}{{\boldsymbol{l}}}_{{\boldsymbol{s}}}^{^{\prime} }}{e}^{i{\boldsymbol{k}}\cdot ({{\boldsymbol{l}}}_{{\boldsymbol{s}}}-{{\boldsymbol{l}}}_{{\boldsymbol{s}}}^{^{\prime} })}{U}_{I{{\boldsymbol{k}}}_{SC}}^{{{\boldsymbol{l}}}_{{\boldsymbol{s}}}{\alpha }^{* }}{U}_{I{{\boldsymbol{k}}}_{SC}}^{{{\boldsymbol{l}}}_{{\boldsymbol{s}}}^{^{\prime} }\alpha }.$$
(10)

Here, \({n}_{s}\) is the number of primitive unit cell of layer \(s\) in one elementary unit cell of the approximant. \({U}_{I{{\boldsymbol{k}}}_{SC}}^{{{\boldsymbol{l}}}_{{\boldsymbol{s}}}\alpha }\) is the projection of \(\left|{\Psi }_{I{{\boldsymbol{k}}}_{SC}}^{SC}\right.\rangle\) (the eigenstate of the approximant) on \(\left|{{\boldsymbol{k}}}_{SC}{{\boldsymbol{l}}}_{{\boldsymbol{s}}}\alpha \right.\rangle\) (the Bloch basis function of approximant), which can be constructed by all atomic orbitals

$$\left|{{\boldsymbol{k}}}_{SC}{{\boldsymbol{l}}}_{{\boldsymbol{s}}}\alpha \right.\rangle=\frac{1}{N}\sum _{{\boldsymbol{L}}}{e}^{i{{\boldsymbol{k}}}_{SC}\cdot {\boldsymbol{L}}}\phi ({\boldsymbol{r}}-{\boldsymbol{L}}-{{\boldsymbol{l}}}_{{\boldsymbol{s}}}-{{\boldsymbol{t}}}_{\alpha }),$$
(11)

where \(N\) is the number of the SC. \(\phi ({\boldsymbol{r}}-{\boldsymbol{L}}-{{\boldsymbol{l}}}_{{\boldsymbol{s}}}-{{\boldsymbol{t}}}_{\alpha })\) is the \({p}_{z}\) orbital located at \({\boldsymbol{L}}+{{\boldsymbol{l}}}_{{\boldsymbol{s}}}+{{\boldsymbol{t}}}_{\alpha }\). Equation 10 indicates that only the eigenstates of approximant are necessary to obtain the spectral function.

Then, the effective band structure can be obtained by58

$$\delta N({\boldsymbol{k}},\epsilon)=\int_{\epsilon -\delta \epsilon /2} ^{\epsilon +\delta \epsilon /2} A({\boldsymbol{k}},{\epsilon }{^\prime})d{\epsilon }{^\prime},$$
(12)

where \(\delta \epsilon\) is the bin width in energy sampling.

Tight-binding propagation method

The density of states (DOS) are calculated by TBPM47 based on the numerical solution of the time-dependent Schrödinger equation. In order to calculate the DOS under a magnetic field and then Hofstader’s butterfly, the hopping is replaced by the Peierls substitution. In TBPM, a random superposition of the \({p}_{z}\) orbitals at all sites is used as the initial state \(\left|{\phi }_{0}\right.\rangle\) with \(\left\langle {\phi }_{0}| {\phi }_{0}\right\rangle =1\). DOS is calculated as Fourier transform of the time-dependent correlation function

$$d(\epsilon)=\frac{1}{2\pi }{\int_{-\infty}^{\infty}}{e}^{i\epsilon \tau }\left\langle {\phi }_{0}| {e}^{-iH\tau /\hslash }| {\phi }_{0}\right\rangle d\tau .$$
(13)

During all calculations in TBPM, we always use systems with more than ten million atoms, for graphene quasicrystal, it is a disk with radius about 204 nm. The open and periodic boundary conditions are used for graphene quasicrystal and the approximants, respectively.

The optical conductivity is calculated by using the Kubo formula in TBPM.47 The real part of the optical conductivity matrix \({\sigma }_{\alpha ,\beta }\) at temperature \(T\) reads

$$\begin{array}{ll}Re{\sigma }_{\alpha ,\beta }(\omega)=\mathop{\mathrm{lim}}\limits_{\epsilon \to {0}^{+}}&\frac{{e}^{-\hslash \omega /{k}_{B}T}-1}{\hslash \omega A} \int_{0} ^{\infty }{e}^{-\epsilon \tau }{\mathrm{sin}}\omega \tau \\ &\times 2Im{\left\langle {\phi }_{2}(\tau)| {j}_{\alpha }| {\phi }_{1}(\tau)\right\rangle }_{\beta }d\tau .\end{array}$$
(14)

Here, \(A\) is the area of the unit cell per layer, and wave functions

$$\begin{array}{lll}{\left|{\phi }_{1}(\tau)\right.\rangle}_{\beta}&=&{e}^{-iH\tau /\hslash }[1-f(H)]{j}_{\beta }\left|{\phi }_{0}\right.\rangle,\\ &&\left|{\phi }_{2}(\tau)\right.\rangle={e}^{-iH\tau /\hslash }f(H)\left|{\phi }_{0}\right.\rangle,\end{array}$$
(15)

where \(f(H)=1/({e}^{\beta (H-\mu)}+1)\) is the Fermi-Dirac distribution operator.

For both density of states and optical conductivities, exponential decayed time windows are multiplied to the time-dependent correlations before performing the Fourier transform, in order to improve the approximation of the integrals. All the results are averaged with a number of different realizations of random sequences in the initial states to reduce the fluctuations appeared in the spectrum due to limit size of the sample.

More detailed description of TBPM is collected in section SIV of Supplementary information.

Structure relaxation

The atomistic model based on the classical REBO59 (intra-layer interaction) and Kolmogorov–Crespi60 (inter-layer interaction) potentials is implemented in LAMMPS software.61 For Kolmogerov-Crespi potential, the kolmogorov/crespi/z version is used in LAMMPS. This method has been used to study the effect of the atomic relaxation on the structure of moiré patterns in twisted graphene on graphite and bilayer graphene62 and it can also accurately reproduce the potential energy surface of the graphite substrate.63