Nature of single-particle states in disordered graphene
We analyze the nature of the single particle states, away from the Dirac point, in the presence of long-range charge impurities in a tight-binding model for electrons on a two-dimensional honeycomb lattice which is of direct relevance for graphene. For a disorder potential , we demonstrate that not only the Dirac state but all the single particle states remain extended for weak enough disorder. Based on our numerical calculations of inverse participation ratio, dc conductivity, diffusion coefficient and the localization length from time evolution dynamics of the wave packet, we show that the threshold required to localize a single particle state of energy is minimum for the states near the band edge and is maximum for states near the band center, implying a mobility edge starting from the band edge for weak disorder and moving towards the band center as the disorder strength increases. This can be explained in terms of the low energy Hamiltonian at any point which has the same nature as that at the Dirac point. From the nature of the eigenfunctions it follows that a weak long range impurity will cause weak anti localization effects, which can be suppressed, giving localization if the strength of impurities is sufficiently large to cause inter-valley scattering. The inter valley spacing increases as one moves in from the band edge towards the band center, which is reflected in the behavior of and the mobility edge.
pacs:72.80.VP,73.20.Fz, 72.15.Rn, 73.22.Pr
Graphene provides a two-dimensional (2D) electron system which is different from conventional 2D systems in many ways. Graphene can be described by a tight- binding model on a 2D honeycomb lattice (1) and the energy dispersion is linear in momentum near the Fermi points at half-filling resulting in Dirac like quasi-particles (2); (3). It is well known that the Dirac states in graphene evade Anderson localization in the presence of weak long range charge impurities (4); (6); (5). Only short range impurities or very strong long range impurities, either of which can cause scattering from one Dirac valley to another, show weak localization correction to conductivity rather than weak-anti localization (7). But the nature of single particle states away from the Dirac point has not been explored in detail so far. Does a higher energy state, away from the Dirac point, get localized in the presence of an infinitesimal strength of disorder, as expected for a conventional 2D system from the scaling theory of localization (8) or does it evade localization like a Dirac state? In this paper we focus on these questions in detail.
An impurity which is long range in real space can not induce large momentum transfer in a scattering event. Therefore, in the presence of weak long range impurities, only intra valley scattering around a Dirac point is possible. A backscattering event around a Dirac point is further suppressed due to Berry phase effect of Dirac states (7). A near backscattering event, then gives a positive correction to the conductivity, that is, resulting in a postitive function and weak-anti localization effect. Thus, although in a conventional 2-dimensional system, all single particle states are localized in the presence of an infinitesimal disorder (8), Dirac states in graphene stay extended for weak long range impurities. Though “trigonal warping” of the spectrum, obtained by keeping terms up to second order in momentum, is necessary to understand the physics of WAL in this model fully (9) and can suppress WAL effects under certain conditions, for graphene WAL has been observed experimentally (10) and numerically (4); (11); (5); (12); (13). When the long range impurity potential, increases in strength, such that the scattering amplitude for intervalley scattering becomes comparable to the intra valley scattering amplitude (as shown in (14)), transition from WAL to weak localization takes place (4); (13).
In this work, we study the nature of all the single particle states of the tight binding model on the honeycomb lattice in the presence of long range charge impurities. For the disorder potential (, which is a widely used model of impurities for graphene (4); (5); (14); (12); (11); (15), we demonstrate that all single particle states remain extended for weak enough strength of disorder. We show clear evidence in support of this from analysis of inverse participation ratio (IPR) and the localization length obtained from the time evolution dynamics of the wave packet. We show that the conductivity increases with the disorder strength and the system size in the weak disorder regime, not only for the charge neutral graphene but also for highly doped graphene. There is a threshold disorder beyond which an energy state shows localized trend. is small for states near the band edge and is maximum for states near the band center which implies the appearance of a mobility edge (ME) in this system. We also provide an intuitive explanation for this surprising result from an effective Hamiltonian calculation. We derive the low energy Hamiltonian for any point on the Brillouin zone and show that it has the same nature as the low energy Hamiltonian at the Dirac point. From the nature of the Hamiltonian and the eigenfunctions, it follows that (as shown e.g. by Ando (7) even before the discovery of graphene) a random potential will cause weak anti localization effects; these can be suppressed, giving localization (7) if the impurities can cause sufficient inter-valley () scattering. Since the inter valley spacing increases as one moves in from the band edge towards the band center, required to localize a single particle state also increases implying the formation of the mobility edge starting from the band edge and moving in, with increasing overall strength of the long range potential.
The rest of the paper is organized as follows. In section I, we describe the model and the method used to study followed by the details of numerical results on IPR, conductivity and diffusion coefficient in section II. In section III, we describe the effective Hamiltonian picture away from the Dirac point in support of our numerical results and end this paper with conclusions and discussions.
Ii Model and Method
The disorder model we study is
Here is the nearest neighbor hopping amplitude, is the chemical potential which can be tuned to fix the average particle density in the system. At half filling . is the random on site potential due to smooth scatterers considered. is the sum of contributions from impurities located randomly at among lattice sites, with being the length of one sub lattice, so that , where is randomly distributed within . We exactly diagonalize the Hamiltonian matrix for a given disorder configuration. Once we have all the eigenvalues, we can calculate physical quantities of interest exactly for that disorder configuration. Data is obtained for various independent disorder configurations and all physical quantities are obtained by doing averaging over many independent disorder configurations. We mainly focus on inverse participation ratio and dc conductivity which is calculated using the Kubo formula. We also study time evolution dynamics of the wave packet to analyze the nature of single particle states.
The results described in the following sections are for and the range of the impurity potential where is the inter atomic distance. We first present results for inverse participation ratio which is crucial to analyze whether a single particle state is extended or localized. In the following subsections we present results for the dc conductivity and time evolution dynamics study respectively.
iii.1 Inverse participation ratio (IPR)
To decide whether a single particle state is localized or not, we calculate the IPR which is defined as
where is an eigenfunction (normalized) with eigenvalue of the Hamiltonian in Eq. 1 and indicates the configuration averaging. In our numerical calculation of IPR(E), we replace the delta function by a box distribution of finite width . IPR of a state is inversely proportional to the portion of space which “participates” in that eigenstate. For plane waves for a d dimensional system of length . But for a general extended state where (16). Therefore, in the thermodynamic limit, IPR vanishes for an extended state while it is non-zero, having the form , for a localized state. Fig. 1 shows for various energy states as a function of . At , for all the states goes as implying that all these states are extended (for details, see Appendix A). As increases, first states near the band edge get localized and then the states in the middle of the band get localized. Eventually at , IPR in the thermodynamic limit becomes non zero for all the states. Our scaling analysis, which shows that the Dirac state () is extended below and is localized for larger values, is consistent with earlier work on magneto conductance calculation where WAL to WL transition was observed for charge neutral graphene at (13). This provides a check for our scaling analysis.
The IPR analysis provides a clear evidence that not only the Dirac states but all the single particle states remain extended for very weak long range disorder. A threshold disorder is required to localize a state at energy . is larger for a state close to the band center as compared to its value for a state near the band edge. Fig. 2 shows how the mobility edge , which separates localized states from extended states, moves as increases. The results shown are for quite big system sizes (maximum ) and disorder averaging is done over 50-100 configurations. Density plots of the wave functions also support the IPR analysis.
Right panel of Fig. 3 shows the spatial extent of the wave functions for an eigenstate with energy using the density plots of for a specific disorder configuration. One can see that at weak disorder, e.g., for all the states, including those near the band edge (e.g., ), is extended over the entire lattice. Further, states close to the band edge show fractal nature even at while states near the band center are very uniformly distributed over the entire lattice, which is consistent with the values of obtained from IPR analysis. For , for has weight only at a few points on the lattice as shown in Fig. 3. Also notice that fractal nature of the wave functions becomes dominant for all the energy states as the disorder strength increases. Energy state at gets localized at comparatively larger disorder value as shown in the bottom panel. This is consistent with our IPR analysis. Left panel of Fig. 2 shows the schematic plot of the mobility edge obtained from the IPR analysis.
In this section we present results for the regular part of the dc conductivity calculated using the Kubo formula details of which are given in Appendix B. Conventionally, in a metallic system scattering with impurities reduces the conductivity. On the contrary our calculation shows an increase in with disorder in the limit of weak disorder.
As shown in Fig. 4, first increases with , reaches a maximum followed up by its decrease on further increase in . Increase in with has been observed for graphene at the charge neutral point , where the Fermi surface consists of Dirac points (4); (11); (12); (5), and is explained in terms of impurity assisted tunneling (17). This relies upon the symplectic symmetry of the Dirac Hamiltonian (17) which is preserved in the presence of weak long range impurities due to the absence of inter-valley scattering. Our results show that even when the system is doped far away from the charge neutral point, shows a trend similar to that of the charge neutral graphene. As we will show in section III, this happens because the effective Hamiltonian away from the Dirac point still belongs to the symplectic symmetry class in the presence of weak long range impurities.
Panel [a] in Fig. 5 shows vs the system size for . In this extreme weak disorder limit, increases with the system size, not only for , which has been seen in earlier work (4); (12); (5); (7), but also for highly doped systems which demonstrates that higher energy states are similar in nature to the Dirac states. Like Dirac states, higher energy states also remain extended for weak disorder and show weak anti localization effects.
Upon increasing further, for small values, continues to show an increase with up to intermediate disorder strength (), implying a positive function, (see Fig. 5[b]) which is consistent with conclusions from earlier work (7); (4); (15); (5). We see signatures of positive function for as well for . But for higher values of , for intermediate disorder remains constant w.r.t . This is because for higher states is smaller than that for the lower states, and is not weak compared to for these states. For higher disorder strength, for all values, is more or less constant for the system sizes studied as shown in Fig. 5. But eventually for large disorder should decrease with as the system enters in the localized regime which could not be captured by the Kubo formula though it is very clearly visible in the IPR analysis shown earlier and also in the diffusion coefficient study which we will discuss in the next section.
iii.3 Time Evolution Dynamics
In this section we study the time evolution dynamics of a wave packet to analyze the effect of disorder on it (18); (19). We start with an initial wave packet , which is an eigenstate of the Hamiltonian in Eq. 1 for a given value of and a particular disorder configuration on a small lattice for which the sub-lattice length is . Here is the site index and is the time (in units of ). We place this wave packet in the middle of a bigger lattice of sub lattice length , which corresponds to certain expectation value of the Hamiltonian, and evolve it in time with respect to the Hamiltonian of Eq. 1 for bigger system size. Then the wave function at any subsequent time is given by where is an eigenvalue of the Hamiltonian for system and is the corresponding eigenvector. Note that the energy of the evolving wave packet is determined as the expectation value of the Hamiltonian operator .
We calculate the expectation value of where is the position operator;
where denotes the position of the th lattice point measured from the center of the lattice. For the clean case, where all the states are extended, . For , the state is diffusive with the constant of proportionality being the diffusion coefficient. For a localized state, the wave packet spreads up to the localization length and then remains constant in time in the limit of large time.
Diffusion Coefficient: We define a generalized diffusion coefficient following (20) as
which depends upon time . We calculated for various values of and the disorder strength . Data shown in Fig. 6 has been averaged over 15 independent disorder configurations.
As shown in Fig. 6, for , for all values keeps increasing with indicating their extended nature in consistency with our earlier results. Note that for , first increases with followed by a saturation of at a maximum value which indicates the diffusive behavior. For intermediate values of , as shown in panel [b] for , , in the large time regime, shows diffusive trend for lower energy states while it decreases with for (close to the band edge) due to localization effects. The time after which starts showing localization trend decreases with increase in disorder strength. Eventually for very high disorder values, as shown in panels [c] and [d], decreases with in most of the time regime for all the energy states.
Elastic Mean Free path: For the diffusive regime, where saturates to a maximum value , one can estimate the elastic mean free path using the relation , where is the Fermi velocity. We extend this definition following (13) and calculate using the maximum value for all disorder regimes. Fig. 7 shows vs for various values of . The data shown is for and . As expected, decreases with increase in the disorder strength. For weaker disorder values and , shows a dip for while the energy dependence gets smoothened out for stronger disorder values. This energy dependence of for weak disorder is consistent with what has been observed earlier (13) for long range disorder.
Most important point to notice is that even for weak disorder, for states near the band edge. This ensures that our analysis of IPR, which suggests that for weak disorder not only the Dirac state but also the states near the band edge remain extended, is reliable. For , for all values of , is smaller than all the system sizes studied for IPR analysis.
Localization Length: We also calculated the localization length from time evolution dynamics of the wave packet. Since gives the spread of the wave packet of energy at time , maximum value of provides a good estimate of . Fig. 8 shows vs for various disorder values. This data is obtained for . For , is of the order of the system size for all states. As the disorder strength increases, decreases eventually becoming much smaller than the system sizes studied indicating the localization behavior.
To ensure whether an energy state is extended for weak enough disorder or not, we did system size analysis of . In Fig. 9, we have plotted vs , in analogy with IPR plots. Since gives the area enclosed by a wave packet, is proportional to IPR. As seen in Fig. 9, for , in the thermodynamic limit for all states except for and . Let is the extrapolated value of in the thermodynamic limit. Hence, for and , is finite indicating the localized nature of these states while for all other energy states is infinite indicating their extended nature. As increases, more energy states in the band show localized behavior. For example for , even shows finite . Eventually for very large disorder , is finite for all the energy states.
Fig. 10 shows vs for various values of disorder. Here for the sake of presentation, wherever is infinity, we represent it by 1000a. As one can see, for , diverges after . For , is finite for a longer range of energy inside the band and shows divergence only for . This analysis supports our proposal and IPR analysis that all energy states remain extended below certain threshold disorder strength resulting in the appearance of a mobility edge in this system (22).
Iv Effective Hamiltonian
In this section we provide an explanation for our numerical findings by looking at the effective low energy Hamiltonian. The tight binding Hamiltonian on the honeycomb lattice can be written as
with . Since the spin symmetry is preserved in the presence of impurities, we will ignore real spin index from now onwards. Here are the nearest neighbor lattice vectors. In the vicinity of a point on an equal energy contour of energy , one can do the Taylor series expansion to get for small where and are functions of . Thus the low energy Hamiltonian around a point on an equal energy contour corresponding to is given by
where are the Pauli matrices and . Note that like the Dirac Hamiltonian, this Hamiltonian is invariant under pseudo-spin reversal symmetry (6); (5). Around the diagonally opposite point on this energy contour , the coefficients change sign while remain same. Eigenfunctions for the Hamiltonian are given by
where correspond to positive and negative eigenvalues and with and . Just like for the Dirac states, due to the non trivial Berry phase effect of the wavefunction, a back scattering process around the point and its time reversed process are out of phase with each other causing suppression of the back scattering process. A near back scattering event, again due to the Berry phase, then gives rise to the opposite sign of the conductivity correction (7) resulting in WAL effect. Note that since we have kept only first order terms in the Taylor expansion, this picture will hold only in close vicinity of the point . Thus only for very weak disorder, which can induce small momentum transfer scatterings, this picture holds true. But it conveys the main message that for extremely weak disorder, all states will show WAL effect.
A strong long range impurity which can result in momentum transfer can scatter an electron from point on the energy contour to the diagonally opposite point , and give rise to weak localization effects. The strength of disorder required to cause scattering from to depends on the energy contour. For energies close to the band edge (), where the energy contours have very small radius (Fig. 11), small disorder is sufficient to cause this scattering and this value of disorder strength increases as one moves away from the band edge towards the band center (). This intuitive picture supports the results obtained from our numerical calculations.
V Conclusions and Discussions
In conclusion, we have analyzed the nature of single particle states in the tight binding honeycomb lattice in the presence of long range charge impurities and demonstrated that not only the Dirac state but all the single particle states remain extended in the presence of weak long range impurities due to generic anti localization effects and prohibited backscattering events. This is because the nature of low energy Hamiltonian and the wave functions are basically of the same nature as for the Dirac states. A threshold strength of disorder is required to localize a state with energy since the random potential needs to scatter an electron from point on an equal energy contour of energy to point on it. Therefore, though states near the band center remain extended for weak to intermediate values of disorder strength, states near the band edge show extended nature only for weak disorder. We do see clear indications of this in the analysis of inverse participation ratio (IPR), charge conductivity, the generalized diffusion coefficient and the localization length obtained from time evolution dynamics of the wave packet.
In this work we presented all the results for a specific choice of the range of impurity potential and the concentration of impurities . Location of the mobility edge for a given value of also depends upon and . For larger values of , we should expect ME to be moving slowly with the disorder strength . In principle, since screening of the potential depends upon the density of states, effective range of the potential itself is a function of the doping. It will be interesting to carry out a full self-consistent calculation for a model of long range impurity potential where range of the potential is decided by the doping. These are questions for future work.
At this end, we would like to mention that for simplicity we have worked with the nearest neighbor tight binding (NNTB) model of graphene which is a very good model at low energy (near Dirac states). Also the 3rd NNTB model does not make much of a difference in the dispersion with respect to what we have for the NNTB model (21). Only the exact band structure calculation shows that there are new somewhat non dispersive high and low energy bands near the gamma point (21). There is no specific reason why these new sub bands should be particularly effective for localization. There cannot be any decay a la Fermi Golden rule, for states at the band edges, which are our main points.
We hope our study will motivate experiments on highly doped graphene to look for extended nature of single particle states away from the Dirac point.
S. N would like to acknowledge Kausik Das (SINP) for computational help.
Vii Appendix A
In this Appendix, we give details of IPR analysis. For plane waves for a d dimensional system of length . But for a general extended state where . Therefore, in the thermodynamic limit, IPR vanishes for an extended state while it is non-zero, having the form , for a localized state. Below we tabulate the values of for various values of energy states and the disorder strength .
As shown in this table, for , for states away from the band edge because is weak disorder for the single particle states away from the band edge, but states near the band edge have values quite off from . As increases, for almost all the energy states, decreases and becomes significantly smaller than .
Viii Appendix B
In this Appendix, we will provide details of calculation of conductivity within the Kubo formalism (23) which is given by
Here is the current-current correlation function, given by
and being the x component of the nearest neighbor vector between sites and . Here is the index for electron spin and since there is a spin degeneracy in the system, each spin makes equal contribution to the charge conductivity. Thus from now onwards, we will skip the spin index.
The tight binding Hamiltonian in Eq. in the presence of disorder can be diagonalized using the transformation
and the diagonalized form is
Here ’s are positive eigenvalues and ’s are negative eigenvalues.
In the transformed basis, the current operator can be written as
Here and . Using this expression for the current operator and with some simple algebra, we get the following Lehmann representation for the current-current correlation function
where is given by
Similarly is given by
can be obtained from by replacing in the indices for wave functions by which also implies replacing by and by respectively. Similarly can be obtained from by replacing by and by . Correspondingly and are replaced by and and vice versa. This is the expression after taking the limit. Here is the Fermi function and . Finally we substitute for and calculate the conductivity using Eq..
- P. R. Wallace, Phys. Rev. 71, 622 (1947).
- A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov and A. K. Geim, Rev. Mod. Phys., 81, 109 (2009).
- N. M. R. Peres, Rev. Mod. Phys. 82, 2673 (2010).
- J. H. Bardarson, J. Tworzydlo, P.W. Brouwer, and C. W. J. Beenakker, Phys. Rev. Lett. 99. 106801 (2007).
- S. Das Sarma, S. Adam, E. H. Hwang, and E. Rossi, Rev. Mod. Phys. 83, 407 (2011) and references there in.
- P. M. Ostrovsky, I. V. Gornyi, and A. D. Mirlin, Phys. Rev. B 74, 235443 (2006), F. Evers, and A. D. Mirlin, Rev. Mod. Phys. 80, 1355 (2008).
- H. Suzuura and T. Ando, Phys. Rev. lett. 89, 266603 (2002).
- E. Abrahams, P. W. Anderson, D. C. Licciardello, and T. V. Ramakrishnan, Phys. Rev. Lett. 42, 673 (1979).
- E. McCann, K. Kechedzhi, V. I. Falko, H. Surat, T. Ando, and B. L. Altshuler, Phys. Rev. Lett. 97, 146805 (2006).
- F. V. Tikhonenko, A. A. Kozikov, A. K. Savchenko, and R. V. Gorbachev, Phys. Rev. lett. 103, 226801 (2009).
- A. Rycerz, J. Tworzydlo, and C. W. J. Beenakker, Euro Phys. Lett. 79, 57003 (2007).
- E. R. Mucciolo, and C. H. Lewenkopf, J. Phys: Condens Matter 22, 273201 (2010).
- F. Ortmann, A. Cresti, G. Montambaux, and S. Roche, Euro Phys. Lett. 94, 47006 (2011).
- Y. Zhang, H. Jiangping, B. A. Bernevig, X. R. Wang, X. C. Xie, and W. M. Liu, Phys. Rev. Lett. 102, 106401 (2009).
- E. R. Mucciolo, and C. H. Lewenkopf, J. Phys: Condens Matter 22, 273201 (2010).
- B. Kramer and A. Mackinnon, Rep. Prog. Phys. 56, 1469 (1993).
- M. Titov, Euro Phys. Lett. 79, 17004 (2007).
- P. Markos, Acta Physica Univ. Comenianae, L-LI, 45 (2009).
- A. Lherbier, B. Biel, Y. Niquet, and S. Roche, Phys. Rev. Lett. 100, 036803 (2008).
- T. M. Radchenko, A. A. Shylau and I. V. Zozoulenko, Phys. Rev. B 86, 035418 (2012).
- A. Cresti, N. Nemec, B. Biel, G. Niebler, F. Trizon, G. Cuniberti, and S. Roche, Nano Res. 1, 361 (2008).
- Though from the time evolution dynamics provides qualitatively the same picture as we obtained from the IPR analysis, there are quantitative differences. The location of the mobility edge for a given disorder are slighly different in two calculations. Main reason for this might be less number of configuration averaging for the time evolution dynamics (as it is very time consuming) compared to the IPR analysis and different bin sizes used in two calculations.
- Many Particle Physics, by G. D. Mahan, 2nd edition (1990).