Unconventional Features in the Quantum Hall Regime of Disordered Graphene: Percolating Impurity States and Hall Conductance Quantization

We report on the formation of critical states in disordered graphene, at the origin of variable and unconventional transport properties in the quantum Hall regime, such as a zero-energy Hall conductance plateau in the absence of an energy bandgap and Landau level degeneracy breaking. By using efficient real-space transport methodologies, we compute both the dissipative and Hall conductivities of large size graphene sheets with random distribution of model single and double vacancies. By analyzing the scaling of transport coefficients with defect density, system size and magnetic length, we elucidate the origin of anomalous quantum Hall features as magnetic-field dependent impurity states, which percolate at some critical energies. These findings shed light on unidentified states and quantum transport anomalies reported experimentally.


I. INTRODUCTION
The role of disorder in the Quantum Hall Effect (QHE) 1 has been essentially related to the existence of a localization/delocalization transition between electronic states, with the formation of critical (extended) states at the center of Landau Levels (LLs) 2 . In very clean samples, the presence of this transition is ensured by the system edges, which force the formation of extended states, while bulk states are localized by the magnetic field. The robustness of the QHE in the bulk limit is guaranteed by the contribution of either weak impurity potentials satisfying the so-called weakness condition 3,4 , strong scattering centers sufficiently far away from each other [5][6][7] , or smooth potentials with long-range spatial variation [7][8][9][10][11][12][13][14][15] . Whenever disorder becomes too strong, all QHE features eventually vanish away.
Under high enough magnetic fields, the electronic properties of graphene are characterized by the presence of a four-fold degenerate zero energy LL (where electrons and holes coexist) together with non-equidistant LLs at energies E n = sgn(n) 2hv F 2 eB|n| (v F is the Fermi velocity, B is the magnetic field and n is the integer LL index) [16][17][18][19][20] . This electronic spectrum results in a Hall conductance quantization σ xy = 4e 2 h (n + 1 2 ) 20 , which is weakly affected by electron-hole puddles or weak surface disorder, but can exhibit further fragmentation of the plateau structure whenever additional symmetrybreaking mechanisms lift the four-fold degeneracy 21,22 . The presence of an additional quantized Hall plateau σ xy = 0 at low energy in high-mobility samples has been for instance assigned to Zeeman splitting or to the formation of quantum Hall ferromagnetism [23][24][25] , with B-field dependent transport scaling behavior conveyed by the dominant symmetry breaking mechanism at play 25,26 .
Recently, several experiments have reported unexplained QHE features in disordered graphene, including sets of extended states in Hall measurements and the formation of a zero energy Hall plateau [27][28][29][30][31] . These features do not fit the usual energy quantization scheme of massless and massive Dirac charge carriers, and are, as such, often generically attributed to disorder. Additionally, the observation of a quantized Hall conductance in highly resistive (millimeter-scale) hydrogenated graphene, with mobility less than 10cm 2 /V.s and estimated mean free path far beyond the Ioffe-Regel limit 32 , suggests some unprecedented robustness of the QHE in damaged graphene 33 .
These findings are considered unconventional in the sense that common belief often states that high concentrations of strong disorder are detrimental to the Hall quantization in 2DEGs. In graphene, however, a variety of literature suggests that things are different and more subtle. The topological contribution to the Berry phase 34 , which is basically a winding number of the pseudo-spin 1/2 35,36 , is predicted to be more robust under disorder, as it should persist even in the presence of sub-lattice symmetry breaking and associated gap-opening. Such behavior has been demonstrated experimentally in hydrogenated graphene 37 . Also, a robust QHE is expected in graphene, even when strong impurities are at a distance smaller than the magnetic length from each other. For instance, for dense impurity concentrations, depending on the symmetry class and impurity strength, a splitting of the critical energy within a single Landau level is predicted when disorder introduces valley mixing [38][39][40][41][42] , similar to splittings al-ready discussed for 2DEG under the influence of certain types of smooth potential 43 . In (quasi-)periodic systems, impurity-engineered Landau levels have been proposed to exist as well 44 .
Up to now, a proper quantitative description of these phenomena for a completely random distribution of disorder and different types of disorder has been lacking, mostly due to computational limitations.
In this Article, by using efficient computational methods, we provide tangible numerical insight into the rich physics of the QHE in disordered graphene, backing up the single-particle scenarios [38][39][40][41][42]44 in which dense distributions of defects can explain the formation of a zero energy plateau in disordered graphene 29,30 and the presence of sets of extended states in Hall measurements 28,31 .
In the presence of single vacancies (SV) and double vacancies (DV), critical states are found to preclude the formation of the usual graphene LLs. Rather, two sets of extended states form at energies different from E n , within each LL. Consequently, at low energy, by tuning the magnetic field and the impurity concentration, a zero energy Hall plateau can be engineered. We extend the present knowledge by characterizing these states following their real-space behavior. We find that the critical states are predominantly located in the impurity dense regions, while the localized states are trapped inside the pristine-like regions. This suggests that the disorder is triggering a percolation of states mechanism, which is also supported by our estimation of the critical exponent. Furthermore, by calculating the transverse conductivity numerically, we prove that σ xy (E) retains quantized values between Landau levels, even in a highly disordered environment. Our numerical approach circumvents the problems commonly associated with Chern number calculation from (i) a computational point of view, namely we don't have to diagonalize matrices containing millions of elements, and (ii) a physical point of view, as the Chern number approach might become ill-defined when the gaps between extended states close due to increasing disorder contributions.

II. MODEL AND METHODS
A single-orbital first-neighbor tight-binding (TB) model restricted to p z orbitals is used to describe graphene, with hopping terms equal to γ 0 and zero on-site energies. Model SV and DV are described by removing the corresponding carbon orbitals, and are randomly distributed on graphene samples containing up to 12 million atoms. SV and DV are short-range scatterers that entail inter-valley scattering [38][39][40] , but with genuine differences. Indeed, SV locally break the sublattice symmetry and induce stronger localization effects, whereas DV locally preserve the sublattice symmetry. Both of these defect models retain electron-hole symmetry, which simplifies calculations and analysis of the physics at play, unhindered by the full complexity of DFT-fitted models, such as for oxygenated graphene 45 .
The order-N method to obtain the dissipative bulk conductivity by wavepacket evolution is already well established 46,47 . By following the time evolution of the wavepackets, length-dependent conductivities σ xx (L) can be extracted, probing diffusive and localization regimes. The effect of a perpendicular magnetic field is modeled through a Peierls phase substitution 48 . As the spin polarization is neglected in present simulations (a factor two is included to take into account the spin degree of freedom), Zeeman splitting is not considered. We don't expect this to alter our conclusions, as is discussed at the end of the paper.
The non dissipative Hall conductivity is calculated using a newly developed efficient real-space algorithm 49,50 [see also Ref. 51 ] as follows: with V the volume of the system, the current operator j x , f the Fermi function and η → 0 a small parameter required for numerical convergence. |φ RP is a random phase state that allows to drastically limit the computation time.
To simulate two-terminal transport, we consider a standard configuration composed of a graphene ribbon with a central region of length L, where DV are randomly distributed. To mimic source and drain contacts, graphene is highly doped outside this region. An on-site energy shift of γ 0 in the Hamiltonian describing the contacts accounts for the doping. We obtain the differential conductance of the system and the spatial distribution of the spectral current by means of the Green's function approach 52 .

A. Density of States
As a first step, we calculate the density of states (DOS) for different vacancy densities in the presence of a magnetic field of 80 T, see Fig. 1. For the chosen densities, the average distance d between defects is in the order of the magnetic length (l B ≈ 25 nm / B/T ), thus favoring strong coupling between impurity states. For the weakest impurity concentrations considered in Fig. 1 (0.05% for DV and 0.125% for SV), LLs still exist at conventional quantization energies E n . Yet, they are broadened due to lingering coupling between impurities (the transition to the non-coupled case is discussed in Sect. III E). The LLs at conventional energies gradually disappear for larger density, when impurity states emerge at energies below and above each pristine LL energy E n for the DV case, and at energies larger than E n for the SV case, in agreement with the predictions for a (quasi-) periodic model 44 . In contrast, the robustness of the zero energy states in the SV case in Fig. 1 is explained by the ranknullity theorem 53 . This theorem predicts the existence of zero-energy modes (ZEM) for dilute concentrations of SV, while they should not form for DV 54 .

B. Longitudinal and transverse conductivity
To study the nature of these impurity states, we perform conductivity calculations for 1% of DV, see Fig. 2 (a), and for 0.25% of SV in Fig. 3(a). By computing both the σ xx (E) (dark color) and σ xy (E) (light color), the occurrence of localized and extended states for different energies is clarified.
Focusing on the low energy region, Figs. 2 and 3 show that, for both SV and DV, the impurity states of Fig. 1 are more extended at the energies E + c and E − c , where σ xx (E) peaks appear. For DV, the DOS does not exhibit any double peak structure (dashed line in Fig. 2), while the conductivity clearly resolves it. The different position of the peaks for SV and DV is reminiscent of their behavior when arranged in periodic arrays 44 .
For both SV and DV, the height of the two peaks (σ xx (E) 1.2e 2 /h) resembles the value for the case where inter-valley mixing leads to two sets of extended states within the same Landau level 40 . This explains the modified transition between quantized Hall plateaus at −2e 2 /h and 2e 2 /h, with an additional plateau at E = 0, observed in Fig. 2 and Fig. 3 for DV and SV, respectively. The difference in strength of DV compared to SV leads to σ xy (E) profiles with varying slope between −2e 2 /h and 2e 2 /h. Such gradient in the localization strength has also been observed in Ref. 42 , where the considered magnetic fields are several orders of magnitude larger. For the selected DV case in Fig. 2(a), the slope at E 0 is not completely equal to zero, thus still providing a limited contribution to the longitudinal conductivity (around 0.7e 2 /h for the calculated length-scale).
The impurity states at E − c and E + c below and above E 0 (position of the conventional zero energy LL) only contribute 2e 2 /h each to the transverse conductance, indicating that they originate from the original LL0, which would contribute 4e 2 /h (spin degeneracy included). Because of the neglected spin polarization in our simulations, these results support the valley-mixing scenario predicted theoretically [38][39][40] , where the existence of the two extended states at E − c and E + c results from the mixing between K and K valleys, within the same LL.
To illustrate the spatial distribution of extended impurity states at low energy, we calculate the projected DOS (PDOS) at E ± c and E 0 in insets (c) and (d) respectively of Fig. 2, and compare it to the local impurity density defined as W i = In addition, to extend these results to the higher energies and to LLs different from LL0, length-dependent conductivities σ xx (L) are considered, whose decay is related to the strength of localization effects, see Fig. 3

(b).
The SV case is chosen as it depicts better energy resolution between extended and delocalized states at high energy. New sets of extended impurity states clearly develop also away from the Dirac point up to −0.2γ 0 , witnessed by the longitudinal conductivity. Extended state energies in σ xx (L) are confirmed by the plateau transitions in σ xy (E) [ Fig. 3(a)]. Contrary to the LL0 case, extended impurity states at higher energy contribute to σ xy (E) with integer multiples of 4e 2 /h; no clear step is observed at σ xy = ±4e 2 /h for instance. Traces of quantization are found at ±6e 2 /h and ±10e 2 /h for SV in Fig. 3(a) . Similarly to the low energy case, these states form at energies different than the ones predicted from conventional Hall quantization. The full 4e 2 /h steps at higher energies (in contrast with the two 2e 2 /h steps in LL0) are rationalized by the fact that, even for the fully periodic case where random disorder broadening is absent, the two new states within each higher-energy-LL are very close in energy 44 . Thus, unrealistically small broadening and a complete absence of disorder would be required to resolve the energy lifting of extended states for LLs different than LL0. Finally, from this σ xx (L) plot, we also note that, although SV significantly contribute to the DOS at E 0 , these states (ZEM) turn out to be strongly localized, in sharp contrast with the DV case where the conductivity at E = 0 remains finite for much larger length scale. In the next section, this different behavior is further scrutinized by simulating different densities.

C. Coupled impurity states
In the very dilute limit (around 0.05%), localized impurity states can be sufficiently separated to have no significantly overlap. In this limit, the conventional QHE is essentially preserved 55 . This limiting case will be considered separately in Sect. III E. By increasing the impurity density, noticeably extended states start to appear when l B is of the order of d as in the case for 1% of DV in Fig. 2(a) and for 0.25% of SV in Fig. 3(a), when localized impurity states couple. This condition can be satisfied by increasing the impurity concentration or the magnetic length (by decreasing the magnetic field), which we demonstrate for several cases in Fig. 4(a) and Fig. 5(a). Once in the coupled regime, the energy of the extended states increases with B 44,56 and concentration, see for instance Fig. 4(a) for 0.5% DV concentration at 320 T.
To analyze the percolation of states driven by impurities, the extracted localization lengths in the strong localization regime, using 57 are plotted in Fig. 4(b) for an impurity concentration of 0.5% (symbols) (similar results are obtained for other concentrations, not shown here). Percolation theory 38,58 predicts a critical exponent ν = 2.34 for ξ ∼ |E − E c | −ν (solid line). Visual agreement is obtained between numerics and theory for the right tail of E + c . However, for the left tail, towards E 0 , agreement cannot be claimed. This behavior results from the set of remnant states in the low impurity density regions of the samples, which are not fully localized for DV at the considered length scale. To further characterize the puzzling behavior of states at E = 0, we plot σ xx (L) in Fig. 5(b) for SV and DV. On one hand, for SV the states are strongly localized. This is in agreement with the localization behavior of ZEM predicted at zero magnetic field 59,60 , following a power-law behavior σ xx (L) ∼ L −2 . On the other hand for the DV case, σ xx (L) exhibits a linear decay. This explains the finite conductivity contributions observed in Fig. 2(a) and Fig. 4(a). Actually, the highest density concentration (2%), with an increased energy split between extended impurity states, even allows resolving the three sets of states at E ± c and E 0 in the σ xx (E) curve. This is counterintuitive in the sense that increasing the disorder decreases the amount of clean patches in the sample, the habitat for delocalized states at E 0 . σ xx (E) at E 0 remains nevertheless surprisingly robust up to long length scales, which is not observed for SV. This could either be explained by the weaker DV disorder strength allowing states to propagate more easily from one pristine patch to the other, or by referring to Ostrovski et al. 55,61 providing an argument against localization at E = 0 for point-like chiral disorder. Ref. 62 recently demonstrated through simulations that it is very difficult and computationally much more demanding to accurately capture the singularity associated to the ZEM, even more so for stronger SV compared to DV. We thus remain cautious about making conclusions on the exact localization behavior at the E = 0 point.

D. Two-terminal Calculations
To gain complementary insight in the length scaling and percolation behavior of impurity states, we also use a different methodology, namely by calculating the twoterminal conductance of a 100 nm wide ribbon for B = 80 T. In the pristine case, the conductance shows a 2e 2 /h plateau in the energy region of Fig. 6(a) (see the inset for a larger energy range). The corresponding spectral current injected from the left contact and transmitted along the chiral top edge channel of the ribbon is indicated by the arrows in Fig. 6(b). For Figs 6(c)-(e), 0.5% of DV are distributed over a ribbon of length L from 25 nm to 100 nm. In quantitative agreement with above calculations in 2D geometry, Fig. 6(a) shows that DV induce bulk states within a certain energy window (from −0.015γ 0 to 0.015γ 0 ), which modify the pristine conductance. For L=25 nm, most of these states are extended enough to connect source and drain contacts and induce a conductance increase well above the pristine value. In fact, the observed conductance for this case is reminiscent of the DOS curve in Fig. 1(a) and the bulk spectral current distribution in Fig. 6(c) illustrates in real space that electrodes are connected by the states and explains the high conductance also seen in 2D simulations. When increasing L to 50 nm, less bulk states are sufficiently extended to allow the electrons reaching the drain contact. As a consequence, the conductance decreases and narrow peaks appear corresponding to the energy regions where more extended states are concentrated, contributing to the transport. The current distribution of Fig. 6(d) indicates that, for this energy, the electron penetration decreases, with a less efficient bridging between electrodes. Analogous behavior is observed for L = 100 nm in Fig. 6(e), with a more pronounced fragmentation into peaks and a weaker conductance decrease at the center and the sides of the DV energy window [ Fig. 6(a)], in agreement with the 2D simulations. For L > 100 nm, the peaks reduce to isolated resonances with conductance below or almost 2e 2 /h, and a transport gap progressively opens around E = 0 (not shown here).
Note that, in the 2D results, the narrow resonances are masked by self-averaging effects. On the other hand, twoterminal simulations should be performed over a large ensemble of disorder realizations to recover the statistical information provided by 2D bulk conductivity simulations, such as the exact position of the more extended states corresponding to the critical energies. The broadness of energy distribution of the DV-induced states and their localization is inversely proportional to the defect density. In the limit of a periodic DV distribution 44 , the bulk states are completely delocalized and concentrated around very specific energies.

E. Level-condensation
In the previous sections, we have always considered a distance between defects short enough to allow their coupling and the formation of impurity states at new critical energies. In Ref. 55 , the transition from interacting impurities to non-interacting impurities in the very low concentration regime has been considered. A distance criterion is provided at which so-called level-condensation should occur, namely For the case of 80T and the lowest concentration considered so far (0.05%), r = 0.63. By reducing the concentration to 0.01%, one gets r = 0.28, for which level condensation should occur. Similarly, by increasing the magnetic field, and keeping 0.05% of DV, one can achieve values of r = 0.32 (for 320T) and r = 0.16 (for 1280T), respectively. The latter approach (high magnetic fieldlow concentration) is simpler from a computational point of view, the former (lower magnetic field -higher concentration) is more realistic from an experimental point of view.
When calculating the DOS, very small impurity concentrations can be considered, as reported in Fig. 7. focuses on the energy region around the LL0. A logarithmic scale was used for the zoom, for the sake of clarity. With the broadening energy set to 0.0005γ 0 , the concentration-dependent width of the peak (∆E) is reported in panel (c) for two arbitrary peak heights. For the dotted line (0.0017 arb. unit), the value of the splitting does not vary for concentrations below 0.01%. This confirms qualitatively and quantitatively (with an uncertainty related to height and numerical broadening) the level condensation from a DOS point of view. A conductivity analysis for the condensated regime is out of the scope of the present study, because all states become localized due to the magnetic field. Either edges or additional long-range disorder would have to be included to allow for extended states to develop at the energies E n of LLs in the conventional pristine quantization.

IV. CONCLUSION
We have investigated numerically the possible origin of anomalous features reported in the quantum Hall regime of low mobility graphene samples [28][29][30] , such as resonances in the dissipative conductivity and a zero-energy Hall plateau. They result from the formation of disorderinduced percolating bulk states, whose density and extension is maximal around two critical energies that depend on the magnetic field and on the impurity density. The presence of defect-induced critical states on novel QHE properties does not seem to depend on the local symmetry breaking (as induced by SV), as they also form for DV impurities. Rather, valley-mixing is required. Finally, a vanishing contribution of these defect-induced critical states is observed in the very low impurity density limit, characterized by so-called level condensation 41 . Results on highly disordered graphene using more realistic impurities, with TB parameters extracted from ab initio simulations, suggest a certain universality of these impurityinduced extended states 45,63,64 , even in the presence of electron-hole asymmetry. This asymmetry can provide additional signature on the nature of the contaminating impurities 65 . The energy dependence of impurity states with magnetic field can be inferred from existing DOS literature for different types of impurities 44,45,66 . The study of electron transport, as we performed in this paper, is nevertheless required to precisely assess the extendedness of states, as is apparent from the energetically unresolved static DOS features for the DV case.
The inclusion of a Zeeman term in the present form of our Hamiltonian, while desirable for high magnetic fields, is not expected to alter our conclusions, as we do not consider spin-spin or spin-orbit interactions at this point. Such term would then simply induce two spindependent copies of the same physics (and an additional trivial splitting with the Zeeman energy). We also note that the present conclusions at high magnetic fields (80 T in this work), which are computationally less demanding for numerical convergence in the transverse conductivity than at low magnetic fields, are expected to be robust for much smaller magnetic fields (where Zeeman interaction is weakened), as is demonstrated in a separate work on oxygenated graphene 45 . In our work, we demonstrate that the unconventional transport features can be explained even with this simplification in neglecting the Zeeman term. Finally, we propose two future lines of research. First, the way the Chern number is modified or not in highly disordered graphene should be investigated (as well as simply explore if Chern number classification is still appropriate). Our method has the advantage to predict the Hall conductivity even for large disorder, but we presently do not have any tool to calculate the Chern number without going through exact diagonalization of a system containing millions of elements. Second, we comment on how interaction effects might play a role. Indeed, high magnetic fields may induce strongly localized states. The interaction between particles may thus become more relevant. Possible influences could be an interaction-induced localization of the critical states or an additional splitting as in the case of quantum Hall ferromagnetism 25 ; or a competition between both mechanisms.