Dodecagonal bilayer graphene quasicrystal and its approximants

TL;DR

This study investigates the electronic and optical properties of dodecagonal bilayer graphene quasicrystals using large-scale calculations and proposes periodic approximants to accurately model their unique non-periodic structure.

Contribution

The paper introduces a series of periodic approximants for graphene quasicrystals and demonstrates their effectiveness in reproducing quasicrystal properties through band-unfolding and tight-binding calculations.

Findings

Emergence of new Dirac points and band gap at M point

Agreement of effective band structure with recent experiments

Lattice mismatch identified as key factor for approximant accuracy

Abstract

Dodecagonal bilayer graphene quasicrystal has 12-fold rotational order but lacks translational symmetry which prevents the application of band theory. In this paper, we study the electronic and optical properties of graphene quasicrystal with large-scale tight-binding calculations involving more than ten million atoms. We propose a series of periodic approximants which reproduce accurately the properties of quasicrystal within a finite unit cell. By utilizing the band-unfolding method on the smallest approximant with only 2702 atoms, the effective band structure of graphene quasicrystal is derived. Novel features, such as the emergence of new Dirac points (especially the mirrored ones), the band gap at M point and the Fermi velocity are all in agreement with recent experiments. The properties of quasicrystal states are identified in the Landau level spectrum and optical excitations. Importantly, our results show that the lattice mismatch is the dominant factor determining the accuracy of layered approximants. The proposed approximants can be used directly for other layered materials in honeycomb lattice, and the design principles can be applied for any quasi-periodic incommensurate structures.

Equations

$$A(\mathbf{k},ϵ)=\sum _{I{\mathbf{k}}_{SC}}{P}_{I{\mathbf{k}}_{SC}}(\mathbf{k})\delta (ϵ-{ϵ}_{I{\mathbf{k}}_{SC}}),$$

$$P_{I{\mathbf{k}}_{SC}}(\mathbf{k})=\sum _{s=1}^2\sum _i{\left|⟨{\psi }_{i\mathbf{k}}^{P{C}_s}⟩{\Psi }_{I{\mathbf{k}}_{SC}}^{SC}\right|}^2=\sum _{s=1}^2{P}_{I{\mathbf{k}}_{SC}}^s(\mathbf{k}),$$

$$P_{I{\mathbf{k}}_{SC}}^s(\mathbf{k})=\frac{1}{n_s}\sum _{\alpha }\sum _{{\mathbf{l}}_{\mathbf{s}}{{\mathbf{l}}_{\mathbf{s}}}^{{}^{\prime }}}{e}^{i\mathbf{k}\cdot ({\mathbf{l}}_{\mathbf{s}}-{{\mathbf{l}}_{\mathbf{s}}}^{{}^{\prime }})}{U}_{I{\mathbf{k}}_{SC}}^{{\mathbf{l}}_{\mathbf{s}}{\alpha }^{\ast }}{U}_{I{\mathbf{k}}_{SC}}^{{{\mathbf{l}}_{\mathbf{s}}}^{{}^{\prime }}\alpha }.$$

$$|{\mathbf{k}}_{SC}{\mathbf{l}}_{\mathbf{s}}\alpha ⟩=\frac{1}{N}\sum _{\mathbf{L}}{e}^{i{\mathbf{k}}_{SC}\cdot \mathbf{L}}\varphi (\mathbf{r}-\mathbf{L}-{\mathbf{l}}_{\mathbf{s}}-{\mathbf{t}}_{\alpha }),$$

$$\delta N(\mathbf{k},ϵ)={\int }_{ϵ-\delta ϵ/2}^{ϵ+\delta ϵ/2}A(\mathbf{k},{ϵ}^{{}^{\prime }})d{ϵ}^{{}^{\prime }},$$

$$d(ϵ)=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }{e}^{iϵ\tau }⟨{\varphi }_0⟩e^{-iH\tau /\hslash }{\varphi }_{0}d\tau .$$

$$Re{\sigma }_{\alpha ,\beta }(\omega )=\underset{ϵ\to 0^{+}}{\lim }\frac{e^{-\hslash \omega /k_{B}T}-1}{\hslash \omega A}{\int }_0^{\infty }{e}^{-ϵ\tau }sin\omega \tau$$

$$\times 2Im⟨{\varphi }_2(\tau )⟩j_{\alpha }{{\varphi }_1(\tau )}_{\beta }d\tau .$$

$${|{\varphi }_1(\tau )⟩}_{\beta }=e^{-iH\tau /\hslash }[1-f(H)]j_{\beta }|{\varphi }_0⟩,$$

$$|{\varphi }_2(\tau )⟩=e^{-iH\tau /\hslash }f(H)|{\varphi }_0⟩,$$

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

$$V_{pp\pi }(\left|{\mathbf{r}}_{\mathbf{ij}}\right|)=-{\gamma }_0{e}^{2.218(b-|{\mathbf{r}}_{\mathbf{ij}}|)}{F}_c(\left|{\mathbf{r}}_{\mathbf{ij}}\right|),$$

$$V_{pp\sigma }(\left|{\mathbf{r}}_{\mathbf{ij}}\right|)={\gamma }_1{e}^{2.218(h-|{\mathbf{r}}_{\mathbf{ij}}|)}{F}_c(\left|{\mathbf{r}}_{\mathbf{ij}}\right|).$$

$$F_c(r)=(1+e^{(r-0.265)/5}{)}^{-1}.$$

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

$$a_t=\frac{M\times 3b}{N}.$$

$$\Delta \rho =\sqrt{\frac{1}{N}\sum _i{\left|{\rho }^{QC}({ϵ}_i)-\rho ({ϵ}_i)\right|}^2}$$

Full text

I INTRODUCTION

The bilayer graphene with van der Waals interlayer interaction shows rich electronic properties that are depending on the stacking order and twist angleTabert and Nicol (2013); Bistritzer and MacDonald (2011a); Lopes dos Santos et al. (2007); Sboychakov et al. (2015); Suárez Morell et al. (2010); Wang et al. (2012); Bistritzer and MacDonald (2011b); Po et al. (2018); Moon and Koshino (2013); Lee et al. (2011); de Gail et al. (2011); Li et al. (2009); Brihuega et al. (2012); Mele (2010); Cao et al. (2018). 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.Rakhmanov et al. (2012); Aoki and Amawashi (2007); Latil and Henrard (2006); Partoens and Peeters (2006). The twist angle between the layers offers an additional degree of freedom to tune the electronic properties. If the commensuration condition is satisfiedShallcross et al. (2010), namely, the twist angle is $\theta=cos^{-1}({{3q^{2}-p^{2}}\over{3q^{2}+p^{2}}})$ , where $p$ and $q$ are integers, the 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° 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 4H-SiC(0001)Ahn et al. (2018), Pt(111)Yao et al. (2018) and Cu-Ni(111)Takesaki et al. (2016) surfaces. The interplay between the quasi-periodicity and interlayer interaction results in the emergence of extra Dirac pointsAhn et al. (2018) and critical eigenstatesMoon et al. (2019). These novel properties are very different from those of graphene, although the twist angle is even more than 15°. Some of them have been explained by the k-space tight-binding method proposed in Ref. Moon et al. (2019). However, the calculated quasi-band structure is supercell-like, and can not be compared directly with the ARPES measurements. The overall understanding of the quasi-periodicity on the electronic and optical properties of graphene quasicrystal is still missing.

In this paper, we study the electronic properties of graphene quasicrystal by using large-scale tight-binding calculations in real space and propose a series of approximants with finite unit cell. 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. Goldman and Kelton (1993) 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. We are going to set up reliable and convincing approximants for recently fabricated graphene quasicrystal with 12-fold rotational order.

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° 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 calculationsFang et al. (2015); Cappelluti et al. (2013); Trambly de Laissardière et al. (2012). 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. Trambly de Laissardière et al. (2012) which has been verified by comparing results with several experiments Trambly de Laissardière et al. (2012); Huder et al. (2018); Shi et al. (2019). This model has been implemented in our home-made code of tight-binding propagation simulator (Tipsi) and reproduces successfully experimentally observed pseudo-magnetic fields in twisted and relaxed bilayer graphene lattice Shi et al. (2019).

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 13–15 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 126491 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 (data not shown). 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)Yuan et al. (2010), 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 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.

II METHODS

II.1 Effective band structure

First of all, the spectral function at wavevector $\bm{k}$ and energy $\epsilon$ can be calculated byNishi et al. (2017)

[TABLE]

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

[TABLE]

where $\left|\psi^{PC_{s}}_{i\bm{k}}\right>$ and $\left|\Psi^{SC}_{I\bm{k}_{SC}}\right>$ 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

[TABLE]

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

[TABLE]

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

Then, the effective band structure can be obtained byMedeiros et al. (2014)

[TABLE]

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

II.2 Density of states

The density of states (DOS) are calculated by TBPMYuan et al. (2010) 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>$ with $\braket{\phi_{0}}{\phi_{0}}=1$. DOS is calculated as Fourier transform of the time-dependent correlation function

[TABLE]

During all calculations in TBPM, we always use systems with more than 10 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.

II.3 Optical conductivity

The optical conductivity is calculated by using the Kubo formula in TBPMYuan et al. (2010). The real part of the optical conductivity matrix $\sigma_{\alpha,\beta}$ at temperature $T$ reads

[TABLE]

Here, $A$ is the area of the unit cell per layer, and wave functions

[TABLE]

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.

III 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$ Slater and Koster (1954)

[TABLE]

Here, $n$ is the direction cosine of relative position vector $\bm{r_{ij}}$ with respect to $\bm{e_{z}}$. The Slater and Koster parameters $V_{pp\sigma}$ and $V_{pp\pi}$ have the following form:

[TABLE]

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 10^{5}$ $m/s$), respectivelyAhn et al. (2018); Yao et al. (2018). $F_{c}$ is a smooth function

[TABLE]

Under a magnetic field, the hopping $t_{ij}$ will be replaced by a Peierls substitutionYuan et al. (2010).

III.1 Approximant construction

The structure of graphene quasicrystal is shown in Fig. 1(a). 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$\degree$ 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$\degree$ 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 calculationsYao et al. (2018).

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

[TABLE]

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

[TABLE]

We name such an approximant as $M/N$, which has the lattice vectors $\bm{a}_{1}=(M\times 3b,0)$ and $\bm{a}_{2}=({{1}\over{2}}M\times 3b,{{\sqrt{3}}\over{2}}M\times 3b)$. The structure of $4/7$ approximant is shown in Fig. 1(b). 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. 1(c). 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. 1(c).

III.2 Density of states

The DOS of graphene quasicrystal is shown in Fig. 2(a). It is calculated from a round disk with radius about 204 nm by using the TBPMYuan et al. (2010), and is in good agreement with previous resultMoon et al. (2019). As shown in Fig. 2(a), 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 Moon et al. (2019). 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.

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. 2(a) and (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. 2(a)). 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. 2(b), 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 derivate

[TABLE]

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. 2(c), 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. 2(a), such as 15/26, 41/71, 67/116 and so on.

III.3 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 Method section for details) is based on tight-binding propagation methodYuan et al. (2010) 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. 2(a). 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. 2(a). 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.

III.4 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 $60\degree n+30\degree$ followed by the mirror reflection with the mirror plane in the middle of the two graphene layers plus just rotation of $60\degree n$ ($n$ is any integer), because after these operations the structure remains unchanged. Two 12-folded rotational symmetric eigenstates of 15/26 approximant are shown in Fig. 4(a) and (b), and the corresponding states in 41/71 approximant are shown in Fig. 4(c) and (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. 4(e) and (f). Such a result is reasonable. For a real graphene quasicrystal, some critical eigenstates expand more than 20 nm in spaceMoon et al. (2019), which can not be simulated by a small approximant, for instance, the approximant with 29.84° twist angle in RefMoon et al. (2019) and the 15/26 approximant in our model. But as the approximant size increased to be large enough, these critical states appear again.

III.5 Effective band structure

According to the generalized Umklapp scattering theoryKoshino and Moon (2015); Koshino (2015), the interlayer interaction will couple the wavevector $\bm{k}$ in the bottom layer and the wavevector $\widetilde{\bm{k}}$ in the top layer if $\bm{k}+\bm{G}=\widetilde{\bm{k}}+\widetilde{\bm{K}}$ is satisfied, where $\bm{G}$ and $\widetilde{\bm{G}}$ are the reciprocal lattice vectors of bottom and top layers respectively. The coupling matrix element $\left|\braket{\widetilde{k},\widetilde{X}}{U}{k,X}\right|$ can also be understood as the scattering from the $k$ state of bottom layer to $\widetilde{k}$ state of top layer, where $\left|k,X\right>$ ($\left|\widetilde{k},\widetilde{X}\right>$) 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|\bm{q}\right|=\left|\bm{k}+\bm{G}\right|=\left|\widetilde{\bm{k}}+\widetilde{\bm{G}}\right|$, because the hopping energy between two layers is isotropic along the in-plane direction. And smaller $|\bm{q}|$ corresponds to the stronger scattering process. The Brillouin zones (BZs) of two graphene layers are shown in Fig. 5(a). 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 $K_{1}$ and $\widetilde{K}_{1}$ after the strongest scatterings and $K_{2}$ and $\widetilde{K}_{2}$ after the second strongest scatterings, respectively. As an example, the first two strongest scattering paths from $K^{{}^{\prime}}$ are shown in Fig. 5(b) and (c) respectively.

Previous ARPES measurementsAhn et al. (2018); Yao et al. (2018) and theoretical studiesAhn et al. (2018) 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 10^{5}$ m/s, respectivelyAhn et al. (2018). 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) substrateYao et al. (2018).

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 measures. The effective band structure (EBS) derived by applying the band-unfold 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. 5(a) is given in Fig. 5(d), 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 theroyAhn et al. (2018) and the ARPES measurementsYao et al. (2018). 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 10^{6}$ respectively in Fig. 5. The Fermi surfaces and the detailed EBS around the six kpoints are shown in Fig. 5(e) and (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 10^{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. 5(d), the band gap at $M$ point below Fermi level, which has been observed in experimentYao et al. (2018), 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 detectedYao et al. (2018). 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 differentAhn et al. (2018). The signal is even stronger for some mirrored Dirac cones than for the original onesAhn et al. (2018). Besides, graphene quasicrystal on 4H-SiC(0001) should be n-type doped according to the ARPES measurementsAhn et al. (2018). 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 theoryAhn et al. (2018) and main experimental resultsAhn et al. (2018).

III.6 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 levelsGuinea et al. (2006); Jiang et al. (2007); Zhang et al. (2006, 2005). The two dimensional Dirac fermion behavior of graphene can be expressed by the Landau levels which satisfies $E_{n}=v_{F}\sqrt{2e\hbar 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. 6(a)). 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 10^{5}$ m/s for hole and electron, respectively, which are obtained by fitting Landau levels under 50 T magnetic field (see Fig. 6(b)). It agrees well with the values from the EBS (9.1$\times 10^{5}$ m/s) and experiments (8.9 $\sim$ 9.3$\times 10^{5}$ m/s)Ahn et al. (2018). Moreover, it is almost the same as but slightly smaller than that in graphene (9.25 and 9.05 $\times 10^{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 10^{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.

IV CONCLUSION

In this paper, 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 10^{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\hbar nB}$ but with a reduced Fermi velocity of $\sim 5.2\times 10^{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., $MoS_{2}$, $WS_{2}$, $MoSe_{2}$ and $WSe_{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 modelCappelluti et al. (2013); Fang et al. (2015), therefore, an approximation with finite unit cell will dramatically simplify the complexities of the modelling. 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.

ACKNOWLEDGEMENTS

This work is supported by the National Key R$\&$D Program of China (Grant No. 2018FYA0305800) and China Postdoctoral Science Foundation (Grant No. 2018M632902). MIK acknowledges a support by the JTC-FLAGERA Project GRANSPORT. Numerical calculations presented in this paper have been performed on the supercomputing system in the Supercomputing Center of Wuhan University.

AUTHOR CONTRIBUTIONS

G. Yu and Z. Wu performed the calculations. G. Yu wrote the manuscript. All authors analyzed the results and revised the manuscript. S. Yuan conducted the project.

Competing interests: The authors declare no competing interests.