Kernel polynomial method to Anderson transition in disordered $\beta$-graphyne

By means of variable moment kernel polynomial method, we analyze the localization properties of $\beta$-graphyne sheet subjected to the Anderson disorder. To detect the localization transition we focus on the scaling behavior of the normalized typical density of states. We find that there takes place a metal-insulator transition and the critical disorder strength is of the order of the bandwidth, which is contrary to the one-parameter scaling theory stating that for infinite two dimensional systems, all the electronic states are localized for an arbitrary strength of the Anderson disorder. As its particular localization properties, it is reasonable to predict there will exist dc conductivity for $\beta$-graphyne at zero temperature.


Introduction
Graphyne, a single layer of carbon sheet, was first predicted by Baughman et al. in 1987 as a layer consisting of sp and sp 2 hybridized carbon atoms [1,2]. Graphyne is similar in structure to graphene but it is constructed by inserting carbon triple bonds (−C≡C−) into single bonds (C−C) in graphene. Although graphyne is topologically equivalent to graphene, the presence of triple bonds introduces rich optical and electronic properties [3] that are quite different from those of graphene, such as higher carrier mobility [4], preferred chemical properties and strongly anisotropic optical adsorption [5,6]. Because of these unique properties, graphyne has inspired researchers to explore its potential applications in electronics, catalysis, energy storage and biomedical fields [7][8][9][10][11][12].
Due to different distributions of sp and sp 2 hybridized carbon atoms, there were proposed various kinds of graphyne. Specially, four different types of graphyne were identified, namely α-, β-, γand 6,6,12-graphyne [1,13]. In reference [13], it is reported that α-, βand 6,6,12-graphyne possess Dirac cone-like band structure around the Fermi level [13]. Unlike graphene, the Dirac points of β-graphyne not only appear at high-symmetry points but also at low-symmetry points. According to these results, it is reasonable to expect that these graphynes can be considered as competitors to graphene and the study of graphyne allotropes will open insights into designing new-generation electronics.
So far, bulk graphyne sample is experimentally unobtainable but can be considered chemically stable. Tremendous efforts have been devoted to synthesizing and assembling a precursor and subunit of graphyne, a class of molecules known as dehydrobenzoannulenes (DBAs) [14,15]. Graphdiyne [16], which is a sub-structure of graphyne with benzene rings linked by diacetylene (−C≡C−C≡C−), was successfully synthesized on a copper substrate via a cross-coupling reaction using hexaethynylbenzene [17]. Moreover, graphdiyne can be fabricated by adopting low-temperature chemical vapor deposition [18] or through a carbon-carbon coupling reaction among constitutive monomers at a liquid/liquid or gas/liquid interface [19]. These exciting experimental progresses indicate that the synthesis of graphyne and its allotropes can be expected in the near future.
Since in real materials disorder is unavoidably present, it is of particular importance to understand how disorder affects the physics of electrons in the fabricated samples. Disorder can arise for various reasons, such as impurities (charged impurities, chemical impurities, etc.), topological defects (vacancies, G.X. Wang edge disorder, etc.), ripples or ad-atoms [20][21][22][23][24]. For graphene, as the quality of samples is improving and synthesis techniques are developing, the dominant disorder varies from bulk disorder to edge defects. So far, researchers are devoted to fabrication of graphyne families and it is reasonable that bulk disorder will play a dominant role in these materials. When modelling a disorder, the Anderson model has been widely used in studying the basic localization features.
For graphene, due to the linear dispersion at the Dirac point, there was debated the issue whether the one-parameter scaling theory can hold [25][26][27][28]. Amini et al. study the effect of on-site uncorrelated disorder on the electronic properties of graphene [25], and find that weak disorder can decrease the velocities of Dirac fermions, with extended states remaining. However, when the disorder strength increases to a large enough amount, there will be localized states around the Fermi point, which is consistent with the results of Naumis [29]. Subsequently, Schileede pointed out for an arbitrary disorder strength that all the states are localized and there exists no mobility edge [30]. Amini et al. acknowledged that the mobility edge might be induced by the kernel [31], but they argued the comment that the localization can occur throughout the whole spectrum because an arbitrary weak disorder strength is unreasonable since there exists minimal conductivity in graphene samples, and Anderson transition was observed by Bostwick et al. [32]. Recently, further study has been undertaken [33][34][35][36][37].
Of late there exist plenty of well established numerical methods, such as exact diagonalization method, Monte Carlo simulations, Lanczos recursion method and kernel polynomial method (KPM). Every method is suited to particular problems and has severe limitations. For example, Lanczos algorithm is a powerful method to calculate extremal eigenstates of sparse matrices but it becomes unstable owing to the loss of orthogonality between the generated basis vectors. However, KPM, which expands the spectral function in terms of a complete set of orthogonal polynomials, can be stable and more convenient to deal with spectral densities and correlation functions. The spectral properties, such as LDOS, can be calculated to arbitrary precision without an exact diagonalization of the Hamiltonian [38].
In this Letter, we consider β-graphyne sheet with the tight-binding approximation and focus on the effect of the Anderson-type bulk disorder [39]. With the modified KPM-the variable moment kernel polynomial method [40] (VMKPM), we calculate the local density of states (LDOS) and analyze the scaling behavior of normalized typical density of states (DOS). To detect the Anderson transition, the normalized typical DOS can serve as an order parameter.

Model and method
β-graphyne consists of a honeycomb lattice of carbon atoms with eighteen sublattices as shown in figure 1. The tight-binding Hamiltonian of β-graphyne with Anderson disorder can be written as where c i (c † i ) is the annihilation (creation) operator destructing (creating) an electron at lattice site i; i, j represents the nearest-neighbor hopping with the hopping amplitude t i j , which is specifically denoted in figure 1; i is the random on-site potential induced by Anderson disorder. Here, we assume that the random on-site potential i has a uniform probability distribution as in which the parameter γ is a measure of disorder strength and θ is Heaviside step function. Many quantities can be used to describe localization properties, such as the averaged localization length, inverse participation number or the distribution of LDOS. Among these quantities, LDOS is the most favorable to study Anderson localization both theoretically and experimentally with the fact that it can be applied to more complex situations and can be directly measured by scanning tunneling microscopy. The LDOS at lattice site i can be defined as where |n is an eigenstate of H with energy E n and |i = c † i |0 . To distinguish a localized state from an extended one, it is convenient to compare two characteristic quantities: the typical DOS and the mean DOS, which are defined, respectively, as where K r denotes the number of disorder realizations and K s is the number of different lattice sites for a given realization. ρ ty and ρ me are both finite for extended states while for localized states ρ ty tends to zero but ρ me remains finite.
To detect the localization transition only for a single finite size system is insufficient [30]. It is reasonable to take into account the scaling behavior of typical DOS, ρ ty . Consider the ratio of the typical and mean DOS R(E), the normalized typical DOS, is positive for extended states whereas it is equal to zero for localized states. Therefore, R(E) can serve as an order parameter to reveal the position of the Anderson transition. By calculating R(E) for different disorder strength and lattice sizes, we can analyze the scaling behavior of the typical DOS. Theoretically, LDOS can be evaluated to arbitrary precision using the KPM [38], which is based on the expansion of the spectral function, such as ρ i (E), into a finite series of Chebyshev polynomials T n (x) = cos(n arccos x). Since Chebyshev polynomials are defined in the interval [-1,1], correspondingly, the Hamiltonian H needs to be rescaled before expansion with the rescaled HamiltonianH = (H − b)/a and rescaled energiesẼ = η is a small positive number and it is not essential to have accurate bounds on the spectrum.
In terms of Chebyshev polynomials, LDOS can be approximated as where the Chebyshev moments  Using the recursion relations between the Chebyshev polynomials, T m+1 (Ẽ) = 2ẼT m (Ẽ) −T m−1 (Ẽ), these moments can be efficiently obtained through a recursive procedure and the expression for µ m can be very simple: µ m = i| α m in which α m = T m (H) |i . The most expensive computation is the multiplication between a sparse Hamiltonian matrix and a vector. If the matrix is computed on the fly, without being stored, dimensions of the order can be carried out to 10 9 or more. Since the Hamiltonian is a sparse matrix, the calculation can be done even on a desktop computer. In contrast to Lanczos, which is unstable because of loss of orthogonality and spurious degeneracy, Chebyshev iteration is completely stable. Since the number of Chebyshev moments µ m in equation (2.7) is finite, it is necessary to convolute the approximated function with kernel polynomials in order to damp out the Gibbs oscillations and ensure properties such as positivity and normalization. In this paper, we use the strictly positive Jackson kernel [41] with φ = π/(M + 1). The broadened peak for an order M expansion of a δ function at position x 0 is a good approximation to a Gaussian of width σ [38] In order to distinguish between resolution and localization [40,42], it is essential to ensure a uniform energy resolution, namely constant σ, over the whole spectrum by restricting the number of moments in equation (2.10).

Calculation and discussion
In our calculation, the number of lattice sites is N = 4.5 × 10 4 − 1.8 × 10 5 with periodic boundary conditions. Following the reference [43], each carbon atom of β-graphyne is described by one 2p z orbital, and the corresponding hopping energies can be obtained, with t 1 = −2 eV, t 2 = −2.7 eV and t 3 = −4.3 eV. In figure 2, we directly calculate the band structure and its associated DOS of ordered β-graphyne without expanding in Chebyshev polynomials, which is in consistence with the results of the [13].
We calculate the typical DOS ρ ty and the mean DOS ρ me of a β-graphyne sheet with 1.8 × 10 5 lattice sites for different sites and disorder realizations K = 32 × 32. The results are shown in figure 3. As β-graphyne possesses the Dirac structure [13], the valence and conduction bands meet at the Dirac points and form the Dirac cone in the vicinity of the Dirac point. Consistently, the DOS is zero at the Fermi level, corresponding to E = 0 in figure 3 (a). For the ordered β-graphyne, ρ me (E) and ρ ty (E) coincide with each other, as shown in figure 3 (a). With disorder increasing, starting from two boundaries of the spectrum, ρ ty (E) is suppressed and when γ > 12t it vanishes, namely, at the critical disorder strength,  the entire spectrum is localized. At the same time, the value of ρ me (E) is always finite, no matter how disorder strength varies.

33706-4
As mentioned above, only a single finite size system is insufficient to study the localization properties. We need to consider the normalized typical DOS, R(E). Figure 4 presents R(E) of five different system sizes for the energy E = 0 as function of disorder strength. As Anderson transition can be influenced by the finite size of the system, and the resolution of KPM, for different system sizes R(E), does not coincide very well but we can find that the curves tend to converge. More importantly, with an increasing disorder strength, ρ ty → 0, which indicates Anderson localization. Then, by calculating R(E), it is possible to reliably determine the critical disorder strength γ c . For example, we obtain the W c 11t [E = 0 and R (E) = 0.05], which is of the order of bandwidth. The existence of critical disorder means that the localization behavior of β-graphyne is similar to 3D case, but contrary to the predicted 2D case of one-parameter theory. Figure 5 shows the contours of the normalized typical DOS R(E) in the energy-disorder plane. Since localized electrons do not contribute to the charge transport, the energy that separates localized states from the extended ones is called the mobility edge. From the contour map, we can find the mobility edge [ρ ty (E)/ρ me (E) = 0.05] and apparently all states become localized for γ > 12t. In figure 5, we can also find that the value of R(E) at two boundaries of the spectrum is finite. However, from figure 5 we know that ρ me (E) [which coincides with the standard density of states ρ(E)] is close to zero, which indicates that the probability of reaching these energy states is particularly small. Thus, these finite values of R(E) will make no sense and the corresponding curve is known as the Lifshitz boundaries.
For a given energy E, ρ i (E) can measure the local amplitude of the wave function at site i. In order to understand the internal structure of the extended and localized states, we calculate the LDOS in the band center (E = 0) for one realization of disorder in figure 6. For the ordered β-graphyne sheet, the 33706-5  distribution of LDOS is uniform [ figure 6 (a)]. When the disorder strength γ = 4t, from figure 6 (b) we can see that although the distribution of LDOS becomes nonuniform, the states still remain delocalized. As the disorder strength increases to γ = 8t, the states clump into some finite regions of clusters, as shown in figure 6 (c), which indicates that the localization takes place. However, between the clusters, there exist bridges where the LDOS is very small but nonvanishing. When the disorder strength reaches γ = 12t, the states are confined to many isolated islands, which implies that these states are totally localized and the metal turns into an insulator, that is to say, the Anderson transition occurs.

Conclusion
In summary, we use typical DOS ρ ty and the normalized typical DOS R(E), which can be efficiently and precisely calculated by means of variable moment kernel polynomial method, in order to detect the localization transition of disordered β-graphyne sheet. When modelling a disorder, the Anderson model is  widely used in studying the basic localization features. One parameter scaling theory predicts that disorder can induce localization transition in 3D systems while all states are localized for infinite 2D case. However, the metal-insulator transition in hydrogenated graphene [32], disordered GNRs [44], and Si-MOSFET inversion layers [45] have been reported, which are out of reach for Anderson model. For a single layer of disordered β-graphyne sheet, we have proved that it also exhibits metal-insulator transition and the critical disorder strength is of the order of the bandwidth, which might be explained by percolation-based approaches [44,46]. In our work, we only numerically analyze the localization properties of β-graphyne. More details, such as the localization length, conductivity [28,47], the contribution of Dirac nodes and the reduced back-scattering, and "resilience" properties of the Dirac nodes against weak to moderate disorder [25,48], need further investigation. Furthermore, we doubt that other allotropes, such as αand 6,6,12-graphyne, may possess similar localization properties, which will be investigated systematically later.

Acknowledgements
The work is supported by the Doctoral Program Foundation of Henan Institute of Technology.