Anomalous Ballistic Transport in Disordered Bilayer Graphene: Dimer Vacancies induced Dirac Semimetal

We report anomalous quantum transport features in bilayer graphene in presence of a random distribution of structural vacancies. By using an efficient real-space Kubo-Greenwood transport methodology, the impact of a varying density of dimer versus non-dimer vacancies is investigated in very large scale disordered models. While non-dimer vacancies are shown to induce localization regimes, dimer vacancies result in an unexpected ballistic regime whose energy window surprisingly enlarges with increasing impurity density. Such counterintuitive phenomenon is explained by the formation of an effective linear dispersion in the bilayer bandstructure, which roots in the symmetry breaking effects driven by dimer vacancies, and provides a novel realization of Dirac semimetals in high dimension.

We report anomalous quantum transport features in bilayer graphene in presence of a random distribution of structural vacancies. By using an efficient real-space Kubo-Greenwood transport methodology, the impact of a varying density of dimer versus non-dimer vacancies is investigated in very large scale disordered models. While non-dimer vacancies are shown to induce localization regimes, dimer vacancies result in an unexpected ballistic regime whose energy window surprisingly enlarges with increasing impurity density. Such counterintuitive phenomenon is explained by the formation of an effective linear dispersion in the bilayer bandstructure, which roots in the symmetry breaking effects driven by dimer vacancies, and provides a novel realization of Dirac semimetals in high dimension. Introduction.-Single layer graphene (SLG) has attracted a great attention owing to its remarkable electrical, chemical and mechanical properties providing an endless list of novel opportunities for practical applications [1]. SLG possesses a linear electronic spectrum with chiral (A-B sublattice) symmetry which lead to exotic low-energy transport such as Klein tunneling [2,3], weak antilocalization [4,5], half-integer quantum Hall effect [6,7].
Bilayer graphene (BLG) differs from SLG by the parabolic band dispersion which however retains the chiral nature of low-energy electronic excitations. One of the salient and unique property of BLG is the possibility of creating an electronic bandgap by applying external gate voltage [8,9]. However surprisingly, the understanding of quantum transport in disordered BLG remains far less understood than the SLG case, because of the enhanced structural complexity. Experimental studies evidence critical differences in transport behaviors of SLG and BLG [10]. Some generalization of the localization theory in BLG has been derived [11], while a minimum conductivity σ min ≃ 4e 2 /h is predicted at the charge neutrality point (CNP) [12][13][14]. A recent scanning tunneling microscopy (STM) study shows that vacancies in graphite induce peculiar impurity states known as zero-energy modes (ZEMs) which are maximally localized at the defect position and then decay as the inverse of the distance from the vacancy [15]. The impact of ZEMs on BLG is so far poorly understood, especially the role played by dimer and non-dimer vacancies which strongly differ in terms of symmetry breaking characteristics.
In this Letter, we start by analysing the localization features of ZEMs in BLG for all types of vacancies, and found a highly inhomogenous sublattice state population (pseudospin polarization) as reported in STM experiments [15]. The depletion of low energy states in one sublattice is additionally further classified into two different classes depending on the vacancy position. Then by using efficient computational methods, we explore quantum transport in disor-dered BLG and analyze how the local nature of the vacancy (dimer versus non-dimer) impact on scattering and localization phenomena. The type and concentration of vacancies are found to dictate the nature of the transport regime which ranges from weak localization to anomalous ballistic conduction.
Electronic features induced by vacancies in graphene.-BLG ( Fig.1(a) and (c)) can be considered as two coupled SLGs with the top layer (in red) shifted a carbon bond from the bottom layer (in black). Consequently, BLG consists of four carbon atoms in its unit cell, two carbons A 1 , B 1 in the bottom SLG unit cell and A 2 , B 2 in the top layer where B 2 lies on the top of A 1 , namely dimer or α sites whereas B 1 , A 2 are called non-dimer or β sites. The tight-binding Hamiltonian model for BLG reads [18,22] where l = 1, 2 labels the bottom and top layer respectively. The annihilation (creation) operators acting on . The first term in H describes the intralayer hopping between nearestneighbor π-orbitals. The second term denotes the interlayer hopping (γ 1 = 340 meV) between dimer sites {A 1 , B 2 }, while the third term gives the interlayer coupling between B 1 and its closet A 2 with γ 3 = 280 meV. The fourth term corresponds to the hopping from A 1 to its nearest A 2 site and from B 1 to its nearest B 2 neighbors (γ 4 = 145 meV). The energy asymmetry between dimer and non-dimer sites is taken into account by introducing in the final term of H an on-site energy difference ∆ = 9.6 meV between dimer sites. All these parameters are derived from the ab-initio calculations [16,18].
The Local Density of States (LDOS) of pristine BLG are first scrutinized (dash curves in the inset of Fig.1(b) or (d)).
Since the two layers are identical one can restrict the discussion to α (dimer) and β (non-dimer) sites. The LDOSs on both sites show a sudden change in the slope at E = ±γ 1 , especially for α sites (black dashed line) owing to the contribution of higher energy bands [22]. The LDOSs also clearly exhibit fingerprints of pseudospin polarization on each layer, in the sense that the state mainly populate β sites (LDOS β = 0 and LDOS α = 0).
Bottom layer A monovacancy on BLG is simulated by removing a carbon atom from the pristine BLG which can be either chosen on the dimer or the non-dimer site. Fig.1(a) and (c) show the sketch of BLG with dimer and non-dimer vacancies (dashed circle), while the decay of the impurity states are given in Fig. 1(b) and (d), respectively. We first consider the effect of a single vacancy on its layer (here the bottom layer). Our results ( Fig.1(b) and (d) inset) show that both dimer and non-dimer vacancies have a strong impact on the layer where the vacancy resides with a strong depletion of ZEM on the vacancy sublattice and a state abundance in the other sublattice (similarly to SLG, see Supp. Mat. [20]). This is seen in Fig.1 (b) and (d) (insets) by comparing the LDOSs of the two nearest sites in the bottom layer (black and red solid lines) with the pristine BLG case (dashed lines). Dimer vacancies also more strongly impact at high energy as seen in the decay of the LDOS at E = ±γ 1 . The spatial distribution of the sublattice abundant state induced by the vacancy in its layer is further investigated in Fig.1(b) and (d) where the data (closed circle) for both dimer ( Fig.1 (b)) and non-dimer vacancy ( Fig.1 (d)) are fitted by power law r −2 (red lines). The depletion of the charge density in the vacancy sublattice can be actually classified into two different classes. The first class involves six second-nearest neighbors of the vacancy and forms a hexagonal lattice (green dashed line in Fig.1(c)) with an enlarged graphene lattice constant by a factor of √ 3, whereas the second class together with the vacancy populate the centers of the previous hexagons. Such long range nature of the impurity state distribution will be key in understanding anomalous ballistic transport. We next investigate the effect of a single BLG vacancy on its adjacent layer. While the dimer vacancy introduces ZEMs on both sublattices of the adjacent layer (the peaks at E = 0 in green and blue curves in Fig.1(b) inset), the non-dimer vacancy leaves the second layer almost unaffected. Indeed, the LDOSs for α 2 (blue solid line) and β 2 (green solid line) site in Fig.1(d) inset are almost identical to the ones for pristine BLG (dashed lines).
Transport properties in the dilute vacancy limit.-Charge transport properties of BLG are investigated for a finite density of vacancies, differentiating the cases with only dimer or non-dimer vacancies (uncompensated cases) from the case with equally distributed mixture of both types (compensated case). Here we also assume an equal distribution among top and bottom layers. We use a real-space order-N wave packet evolution approach [23]. The Kubo-Greenwood conductivity is written as σ(E, t) = e 2 ρ(E)∆X 2 (E, t)/t, where ρ(E) is the DOS and ∆X 2 (E, t) is the mean quadratic displacement of the wave packet and gives the diffusion coefficient D(E, t) = ∆X 2 (E, t)/t. With disorder, D(t) changes from a ballistic motion to a saturation regime, from which the mean free path ℓ e is deduced through ℓ e (E) = D max (E)/2v(E) (v(E) is the velocity and D max the maximum value). At long times, D(E, t) eventually decay owing to quantum interferences which drive the system either to weak or strong (Anderson) localization regime [21]. Fig.2(a) shows the total DOS for BLG with 0.05% of dimer vacancies α (red solid curve), non-dimer β (blue solid curve) and mixed α − β case (black solid curve), together with the total DOS for pristine case (dashed curve). As already suggested in Fig.1 (b) and (d) (insets), especially at the CNP, dimer vacancies more strongly affect the DOS than non-dimer ones, whereas the depletion of state at E = ±γ 1 is unobservable. However, a marked difference is observed on the mean free path l e (Fig.2(c)) and semiclassical conductivities σ sc = e 2 ρ(E)D max (Fig.2(d)) with two significant peaks at E = ±γ 1 . The values of l e (resp. σ sc ) at such energy are higher for the dimer than for the non-dimer vacancies by a factor of 3 (resp. a factor of 2) indicating a stronger scattering efficiency for non-dimer vacancies. This likely originates from the imbalance of state at E = ±γ 1 on the two sublattices around the non-dimer vacancies ( Fig.1(b),(d) insets). The stronger impurity scattering strength at high energy for non-dimer vacancies is also confirmed by the fact that l e and σ sc are almost the same for the compensated case (black solid curves in Fig.2(c),(d)) and the non-dimer vacancy case (blue solid curves). For three cases l e are almost the same but σ sc for the dimer vacancy is larger due to enhanced DOS. D(E, t) further evidence much stronger localization effects for the dimer vacancies ( Fig.2(b)), consistently with the DOS results ( Fig.1 (b) and (d)). Importantly in the low impurity density limit (≤ 0.05%), regardless the nature of vacancies, a localization regime is obtained for the whole energy spectrum. However, surprising features are obtained for high enough vacancy density when the affected spatial areas around vacancies start to overlap. One first scrutinizes the ratios of semiclassical conductivity σ sc for the compensated case in the dilute limit with 0.02%, 0.05% of vacancy, and dense limit with 0.5% of vacancies (Fig.2(d) (inset)). In the dilute limit σ sc perfectly obeys the Fermi's golden rule ∼ 1/n i in the whole energy band (black solid curve in Fig.2(d) (inset). The Fermi's golden rule, however, underestimates (overestimates) the values of σ sc at CNP (E = ±γ 1 ) in the dense impurity limit (red solid curve in Fig.2(d) (inset), i.e. when increasing the vacancy coverage σ sc decreases faster at E = ±γ 1 and slower at the CNP. Transport properties in dense vacancy limit.- Fig.3 and Fig.4 show results for 0.5% of non-dimer vacancies and compensated case (Fig.3 ) and dimer vacancies (Fig.4). Charge transport in presence of non-dimer vacancies and for the compensated case similarly exhibit a localization behavior (as seen in Fig.3(a)). The low-energy semiclassical con- ductivities also clearly saturates at σ min sc ≃ 4e 2 /h (dashed curve) [14].
The calculation of the quantum correction of the semiclassical conductivity ∆σ(L) = σ(L) − σ sc at CNP (Fig.3(c)) further confirms that the non-dimer vacancies induce localization, since a very good fit ∆σ(L) ≃ a 2e 2 h ln(L/l * e ) is obtained with a = −0.295 ∼ −1/π in agreement with 2D weak localization theory (∆σ(L) ≃ − 2e 2 πh ln(L/l * e )). The case of high density of dimer vacancy is remarkably different from the non-dimer one. Fig.4(a) shows D(E, t) for 0.5% dimer vacancies at CNP (red solid line) and at high energy E = 0.05γ 0 (red dashed line). While high energy charge transport quickly enters the localization regime (see D(E = 0.05γ 0 , t)), in sharp contrast D(E = 0, t) increases linearly without any sign of saturation, as expected in a truly ballistic motion. This remarkable behavior is captured in the quantum conductivity scaling σ(t) in the vicinity of CNP. Fig.4(b) shows that σ(t) at t = 5 ps (dashed lines) is smaller than its value at t = 40 ps (solid lines) around CNP. A crossover from the ballistic regime to the localization regime actually occurs at a critical energy E c = 0.008γ 0 (red dotted lines). This anomalous behavior is also observed for a higher dimer vacancy density of 5% ( Fig.4(a, b) (blue lines))., for which a significant enlargement of the ballistic-transport energy range is obtained (E c = 0.017γ 0 (blue dotted lines)). This is a remarkably counterintuitive phenomenon in which the quantum conductivity increases with disorder density. As explained below, the interlayer hopping between nondimer sites γ 3 governs such anomalous behavior. First, by setting γ 3 = 0 for the case of 5% dimer vacancies, D(t) is seen to exhibit a strong localization with a very fast decay at long times and a value reduced by two orders of magnitude when compared to the case with γ 3 = 0 ( Fig.4(a) (green solid line)).
To rationalize such effect, let us then consider a simple tight-binding model by suppressing the terms γ 4 and ∆ in Eq.(2), which is equivalent to the continuum model [9,22] It describes two possible interlayer hopping events between non-dimer sites, the ones with dominant charge density ( Fig.1(b) and (d) inset (dashed lines)). The first term activates hopping via dimer sites which generate a mass term m = γ 1 /2v 2 and a parabolic band [9]. The second term describes a direct hopping γ 3 between non-dimer sites in the bottom and top layers with velocity v 3 = ( √ 3/2)aγ 3 /h. This term is dominant at low energy and generates a linear energy dispersion [9,22] where ξ = +1 (−1) corresponds to the K (K ′ ) valley and φ labels the momentum direction. This linear dispersion occurs for ǫ < 1 2 γ 1 (v 3 /v) 2 . Hereafter we will show that this value is enhanced in presence of dimer vacancies and yield anomalous ballistic motion. The impurity state created by a dimer vacancy is mainly localized on non-dimer sites, and has almost no weight on dimer sites around the defect. Such sublattice polarization induced around the vacancy is shown in Fig.1(b). Fig.4 (c) further evidences the depletion of LDOS on the dimer sites of the layer where the vacancy is lying. Similarly to the SLG case (see [20]), the LDOSs on dimer sites can also be separated into two classes (circle and square symbols) which are both well fitted with the scaling law ρ = ρ 0 − br −0.9 (black solid lines) where ρ 0 is the LDOS far away from the vacancy. From Eq.(2), one observes that the abundance of charge density on the non-dimer sites leads to the enhancement of the second term in Eq.(2) which corresponds to the direct hopping between non-dimer sites. The reduction of electronic state on the dimer sites, on the other hand, leads to a reduction of the first term in Eq. (2). Indeed, this term involves a hopping process in which electron has to hop into a dimer site. This process is limited by dimer vacancies because it creates around it, in dimer sites, a wide area of depletion of electron density, i.e. prohibiting electron residence on the dimer sites in the vacancy vicinity. We can thus consider BLG with a high enough density of dimer vacancies as an effective pristine BLG with renormalized parameters v * < v and v * 3 > v 3 . The more dimer vacancies in the BLG, the smaller the renormalized value of v * and the larger the renormalized value of v * 3 . This thus leads to the expansion of the energy range ( 1 2 γ 1 (v * 3 /v * ) 2 ) where the linear dispersion dominates. This also explains the increase of critical energy E c with the dimer vacancy density in Fig.4(b). Such chiral electrons in the renormalized linear BLG energy dispersion share many similarities with low-energy excitations propagating in SLG (the sublattice index is replaced by the layer index) yielding backscattering suppression. This situation thus provides a novel realization of a Dirac semimetal in high dimensionality, which is an unconventional transport regime provoked by symmetry effects of impurities. This could motivate a systematic exploration of quantum transport in irradiated BLG, following pioneering experimental studies [15].