Efficient Linear Scaling Approach for Computing the Kubo Hall Conductivity
Abstract
We report an orderN approach to compute the Kubo Hall conductivity for disorderd twodimensional systems reaching tens of millions of orbitals, and realistic values of the applied external magnetic fields (as low as a few Tesla). A timeevolution scheme is employed to evaluate the Hall conductivity using a wavepacket propagation method and a continued fraction expansion for the computation of diagonal and offdiagonal matrix elements of the Green functions. The validity of the method is demonstrated by comparison of results with bruteforce diagonalization of the Kubo formula, using (disordered) graphene as system of study. This approach to mesoscopic system sizes is opening an unprecedented perspective for socalled reverse engineering in which the available experimental transport data are used to get a deeper understanding of the microscopic structure of the samples. Besides, this will not only allow addressing subtle issues in terms of resistance standardization of large scale materials (such as wafer scale polycrystalline graphene), but will also enable the discovery of new quantum transport phenomena in complex twodimensional materials, out of reach with classical methods.
pacs:
73.43.f,72.80.Vp, 73.22.PrI Introduction
The Hall effect is one of the central phenomenon in Condensed Matter physics with great importance for measuring carrier density and mobility in semiconductors. The long history started almost 150 years ago when E. H. Hall observed a transverse voltage drop over a conducting bar driven by a longitudinal current Hall (1879) and it experienced an important turn in 1980 when K. von Klitzing observed a quantized version of this effect, the socalled Quantum Hall Effect (QHE)v. Klitzing et al. (1980). This had an immediate impact and triggered huge amount of work which led to the definition of a resistance standard.Porier and Schopfer (2010)
The effect of disorder and related localization phenomenon is central to understand the integer QHE from a bulk perspective, where plateaus develop by varying the charge density, while the plateau region coincide with a region of vanishing . Aoki (1987) This is explained by Anderson localization of the wave functions in the Landau levels induced by disorder. The most standard method to treat bulk conductivities microscopically is the linearresponse theory, or the Kubo formula, which was first applied by Aoki and Ando Aoki and Ando (1981) to the QHE problem Aoki (1987). Such quantity can be connected to measurements in the Corbino geometry for which electrodes are attached to the inner and outer perimeters of an annular sample, which allow to investigate bulk transport physics in the high magnetic field regime.
Recently, the QHE has played a crucial role for the first unambigous proof of the existence of single atomic layers of graphene Novoselov et al. (2005); Zhang et al. (2005); Novoselov et al. (2007). Indeed, in graphene, the peculiar Berry phase related to the pseudospin quantum degree of freedom results in a unique spectrum with fourfold degenerate Landau levels (due to the spin and valley degeneracies), and halfinteger QHE with Hall plateaus emerging at , with integers . QHE related research in graphene is a vivid field, studying electronelectron interaction Nomura and MacDonald (2006); Young et al. (2012), superlattices Ponomarenko et al. (2013); Dean et al. (2013); Hunt et al. (2013); Yang et al. (2013), graphene functionalization and topological currents. Additional transverse transport phenomena are being studied such as (quantum) anomaleous Hall effect,Nagaosa et al. (2010); Chang et al. (2013), spin Hall effect Kato et al. (2004) and quantum spin Hall effect. König et al. (2007)
The role of disorder in graphene is also debated and the QHE features, with very robust behavior to large disorder Guillemette et al. (2013), also present additional peculiarities depending on the symmetry breaking aspects conveyed by defects.Ostrovsky et al. (2006); Tworzydło et al. (2006); Morpurgo and Guinea (2006); Peres et al. (2006); Gattenlöhner et al. (2014) Recently, we have studied the case of oxidized graphene, and have shown the formation of disorderinduced resonant critical states appearing in the zeroenergy Landau level with finite and suggesting a zeroenergy quantized plateauLeconte et al. (2014). The confirmation of such a plateau will demand to revise the description of topological invariant in disordered graphene, an exciting direction of work. However a real space calculation of is needed and demand for the development of new methodology.
Indeed, the importance of an efficient computational methodology for the Hall conductivity stands in sharp contrast to the limited possibilities to simulate this quantity with usual approaches. One of them starts with the Kubo formula and requires the knowledge of the full eigenstates of a given system. Aoki and Ando (1981). Such formulation allows for advanced analysisSheng et al. (2006); Goswami et al. (2008) but suffers from unfavourable cubic scaling with system size which restricts to small systems because diagonalization is necessary. All eigenstates are needed and have to be combined to calculate . This computational limitation makes the analysis of realistic models of disordered samples a priori impossible. Additionally, such approach is restricted to unreasonably strong magnetic fields that exceed by far what is achievable experimentally, and forces the magnetic length scale to dominate over all other length scales of the problem (mean free path, average distance between impurities,…) limiting the exploration of complex magnetotransport in intermediate or even lowmagnetic field regimes, thus calling for approaches with improved numerical performance.Ortmann and Roche (2013); J. H. Garcia (2014); Ishii et al. (2011)
In this paper, we present a highly efficient order algorithm for the computation of the Kubo Hall conductivity that allows us to circumvent aforementioned limitations, and offer fascinating perspectives for exploring transverse transport phenomena in disordered twodimensional materials with almost arbitrary complexity. A validation of the method on several models of clean or disordered graphene is made by comparing this new time evolution method with brute force diagonalization technique. Some preliminary illustration of the method have been published in a prior work Ortmann and Roche (2013). The method is inspired by the real space implementation of the dissipative Hall conductivity at zero frequency, which has proven undisputed computational efficiency and predictive power (see Ref.Roche and Mayou (1997); Foa Torres et al. (2014) for computational details, and Ref. Roche et al. (2012) for some illustration on realistic models of disordered graphene).
Ii Methodology
Kubo formalism – The general framework is provided by Kubo’s linear response theoryKubo (1957); Mahan (2000) in which the dc conductivity is given as
(1) 
Starting from Eq. (1), the evaluation of the conductivity proceeds by introducing an orthonormal manyparticle basis with which allows to rewrite Eq. (1) into
(2) 
where is a small parameter necessary to converge this expression.
In the case of noninteracting electrons one can proceed further towards the singleparticle picture and obtain
(3)  
where the FermiDirac distribution is described by and the chemical potential has been introduced. The singleparticle eigenvalues can in principle be obtained by matrix diagonalization of finite systems which has been widely employed in the literature for studying model systems with phenomenological description of disorder. Unfortunately such numerical operations quickly makes the methodology computationally prohibitive, in particular for simulating complex and large system sizes closer to the experimental ones and at moderate magnetic fields. Clearly, the computational problem then becomes untreatable with diagonalizationbased methods, i.e. it becomes impossible to calculate (and store) all eigenvalues and eigenvectors, so that timeevolution based approaches appear promising. By introducing we arrive at an alternative representation of the Hall conductivity
(4)  
which proves useful for the numerical implementation. In Eq. (4), the trace is replaced by an average over a randomphase state according to . Such states are normalized by using the number of sites . We further introduce the projection operator , as well as the timedependent fulfilling with
(5)  
where
(6) 
are the the matrix elements of the Green function (with ) and is the volume per site. Equation (5) is the basis for the numerical evaluation of the Hall conductivity for which one eventually takes the limit . One notes, however, that the result does not converge in all cases when taking this limit and an oscillatory behaviour of the integral is generally observed, particularly for the case of clean systems in high magnetic fields (and well separate Landau levels), which is the usual test reference for Kubo Hall transport simulations.
The reason is that the Kubo formula Eq. (1) assumes that at “infinite past” times, the system is in equilibrium, i.e. at zero electric field, the current should vanish. In other words, the currentcurrent correlation function should decay with time difference . This should be fulfilled also for the final result. In pristine systems this is not the case and must be imposed a posteriori. In order to take this into account an additional small parameter has to be introduced such that the integral
(7)  
converges.
We note that the same parameter has to be included in full analogy in the standard expression for the calculation of the Hall conductivity based upon the knowledge of the eigensystem. We here display such formula Aoki and Ando (1981); Aoki (1985)
(8)  
which is identical to Ref. Aoki (1985) for the reasonable choice that . These two equations (7) and (8) are comparable and will be evaluated for a couple of examples in Sect. III as numerical checks.
Numerical implementation – The calculation proceeds by splitting the timedependent and time independent parts () while the sum over in (7) controls the energy dependent convergence (see details below). Firstly, in order to evaluate Eq. (6), we make use of the Lanczos tridiagonalization of the Hamiltonian which allows us to use the continued fraction expansion for the diagonal Green function. From the diagonal elements one can obtain the necessary offdiagonal elements in a second iterative step. We find a simple recursion relation connecting both which reads for
(9) 
while the first step in the iteration is given as .
Secondly, the time evolutions of the quantities and can be performed efficiently through a Chebychev expansion of the timeevolution operator which usually converges quite rapidly. For the final calculation of the brackets in Eq. (7) we have to calculate off diagonal matrix elements with the advanced Green operator . For the calculation of this second set of offdiagonal Green functions in Eq. (7), we use the continued fraction expansion. All operations are linear in sample size and this holds also for the total expression since the scalar product in is restricted to a cutoff value . We use if not stated otherwise throughout the paper.
System description The electronic properties of graphene can be described by the standard orthogonal tightbinding model
(10) 
where is the nearest neighbor transfer integral and the onsite energy, which can be chosen to describe various disorder models. Spatially uncorrelated Anderson disorder is introduced through onsite energies taken at random from where gives the disorder strength. This is a commonly used disorder model for exploring the metalinsulator transition in lowdimensional systems Evers and Mirlin (2008); Sheng et al. (2006), which, in average, does not break electronhole symmetry or inversion symmetry. The inversion symmetry of the system related to the symmetry of AB sublattices can be artificially broken by a staggered potential of the form for A(B) sublattices (including all atoms). A less artificial variant of such potential can be introduced by diluting such sublattice terms for a weaker and random distribution, which mimicks some correlation in the potential beyond the Anderson model.
The magnetic field is implemented through a Peierls phase Luttinger (1951) added to which determines the magnetic flux to per plaquette. Ortmann et al. (2011) Spin degeneracy is assumed throughout the paper and an antisymmetrization procedure for has been performed consistent with electronhole symmetry for all considered disorder potentials.
Iii Results
Before we enter into the discussion of the physical results, we illustrate the convergence behaviour of the numerical algorithm which can be controled by the expression (6). A key quantity for the convergence is which can be studied prior to the time evolution of wavepackets and allows to estimate the required number of recursion steps which depend on the specific physical details (e.g. disorder, magnetic field), on the selected energy, and on the broadening . Fig. 1 presents the modulus of plotted versus for some selected energies in the case of pristine graphene (a) and disordered graphene (b). As a general trend, Fig. 1 shows that in case of pristine systems, the convergence of the Lanczos recursion can be reached with a low number of recursion steps (), while increasing disorder requires to be increased typically to a few thousands.
Turning to the simulations of the Hall conductivity, we first discuss the results for pristine graphene in Fig. 2 which shows the comparison of diagonalization method and timeevolution Kubo approach. Figure 1 (inset) shows the case of extremely high magnetic fields (fluxes). We recall that, in order to fulfill the boundary conditions of the periodic system, the limit of very high fields is usually considered by diagonalizationbased studies, because only small sample sizes can be treated by matrix diagonalization. The magnetic field is simultaneously given in the inset of Fig. 2 in terms of the total flux penetrating the sample. The corresponding results from the timeevolution Kubo method (TEK) is plotted as dashed lines for the same magnetic fields (values indicated in the inset legend). For the selected energies, the curves coincide visually. The square root scaling of their position with magnetic field is evident and Hallconductance steps occur at the expected positions and with the expected height corresponding to the halfinteger sequence .
Next, we consider the interesting case of moderately high magnetic fields such as frequently observed in current experimentsYoung et al. (2012); Ponomarenko et al. (2013); Dean et al. (2013); Hunt et al. (2013); Yang et al. (2013) in the main frame of Fig. 2. The Kubo approach allows to reduce magnetic fields down to only few Teslas where the first Landau levels appear at Fermi energies of few tens of meV. Still the plateau sequence is obtained with very good accuracy, which shows that such approach can be of great practical use. Most notably, the figure demonstrates conductance quantization for fields as low as 11 Tesla. As an effect of finite energy resolution, we observe that the relative sharpness of the steps is reduced with lowering the field in this regime. The conductance slope starts to become noticable as Landau levels get closer (as seen for 11.4 T at ). This consequence of the considered energy resolution (here ) can be systematically improved further by increasing the number of recursion steps and decreasing .
As a second example to illustrate the performance of the algorithm, we focus on the case of homogeneously disordered graphene. The case of Anderson disorder is chosen for Fig. 3, because this disorder is expected to induce a transition from the QHE regime to the conventional Anderson insulating state for sufficiently high value of . The mainframe indicates very good agreement (within few percents) obtained between the Kubo algorithm and exact diagonalization techniques at low energy for very high magnetic fields. At higher energy, allows for fast convergence, at the expense of small loss in information compared to the exact method when Landau levels are getting closer in energy. To illustrate the effect of increasing disorder, a zoom on the first two Hall steps is provided for T in the upper inset. The Landau levels are broadened with increasing around the critical states, as expected from perturbation theory. Simultaneously, the first and second plateaus remain at and indicating robust QHE. In the lower inset, a more realistic magnetic field is chosen ( T) to illustrate (i) the performance of the algorithm to probe low magnetic fields and (ii) the possibility to destroy conductance quantization at higher energy for . For such disorder, only the zeroenergy Landau level fully develops, while all states become localized for .
As another challenging example, we present the analysis of a weakly correlated potential. We chose a globally sublatticesymmetry breaking potential that is locally disordered. In contrast to the above Anderson disorder model, which shifts all onsite energies at random, the present potential is only locally present, i.e. only for a fraction of sites that are selected at random. The disorder is chosen such that it breaks globally the ABsublattice symmetry. For the results shown in Fig. 4 (main frame), we use a strength of with of sites affected. Besides ordinary conductance quantization at the values , an additional plateau is induced at zero energy with zero Hall conductivity.^{1}^{1}1The case of 45 T has been previously published in Ref. Ortmann and Roche (2013). The plateau is clearly visible for fields between 9 T and 45 T. We find the width of the zeroenergy plateau equal for all calculated magnetic fields. The corresponding energy value of the onset of the plateau at is determined by the strength of the potential which is given as (indicated by dashed vertical line). The transition to higher plateaus (i.e. from to ) is rather influenced by the magnetic fielddependence of Landau level position.
It is interesting to study the emergence of such plateau when gradually increasing the concentration of impurities. Results are plotted for the case of 9T in Fig. 4 (inset). The plateaulength scaling with is evident from the figure. For the considered energy resolution in the inset (), a trace of plateau onset is still visible for the lowest percentages of the impurity potential down to 0.5% (corresponding to or 2.7 meV).
Iv Conclusion
To conclude, an efficient real space algorithm to compute the Hall conductivity with the Kubo formula has been presented and validated on simple graphenebased systems. Such approach should become a useful computational tool to simulate QHE in very large size complex disordered materials such as polycrystalline grapheneCummings et al. (2014), graphene subjected to weak van der Waals interaction such as by a boronnitride substrateMartinezGordillo et al. (2014), or other types of disordered twodimensional materials. It should also allow to corroborate the formation of zeroenergy plateaus of driven by disorderinduced critical states, and disconnected from degeneracy lifting of Landau levelsLeconte et al. (2014), an issue of genuine fundamental interest in the context of topological interpretation of the quantized conductance.
Acknowledgements.
F.O. would like to acknowledge the Deutsche Forschungsgemeinschaft for financial support within the Emmy Noether scheme (grant OR 349/11). ICN2 acknowledges support from the Spanish Ministry of Economy and Competitiveness under contract MAT201233911 the Severo Ochoa Program (MINECO, Grant SEV 20130295). We acknowledge PRACE for awarding us access to the Curie supercomputing center based in France. The support of Mikolaj Szydlarski from MDLS HPC DevTeam, France to the technical work is gratefully acknowledged.References
 Hall (1879) E. H. Hall, Am. J. Math. 2, 287 (1879).
 v. Klitzing et al. (1980) K. v. Klitzing, G. Dorda, and M. Pepper, Phys. Rev. Lett. 45, 494 (1980).
 Porier and Schopfer (2010) W. Porier and F. Schopfer, Nature Nanotechnol. 5, 171 (2010).
 Aoki (1987) H. Aoki, Reports on Progress in Physics 50, 655 (1987).
 Aoki and Ando (1981) H. Aoki and T. Ando, Sol. State Commun. 38, 1079 (1981).
 Novoselov et al. (2005) K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos, and A. A. Firsov, Nature 438, 197 (2005).
 Zhang et al. (2005) Y. Zhang, Y. Tan, H. L. Störmer, and P. Kim, Nature 438, 201 (2005).
 Novoselov et al. (2007) K. S. Novoselov, Z. Jiang, Y. Zhang, S. V. Morozov, H. L. Stormer, U. Zeitler, J. C. Maan, G. S. Boebinger, P. Kim, and A. K. Geim, Science 315, 1379 (2007).
 Nomura and MacDonald (2006) K. Nomura and A. H. MacDonald, Phys. Rev. Lett. 96, 256602 (2006).
 Young et al. (2012) A. F. Young, C. R. Dean, L. Wang, H. Ren, P. CaddenZimansky, K. Watanabe, T. Taniguchi, J. Hone, K. L. Shepard, and P. Kim, Nature Phys. 8, 550 (2012).
 Ponomarenko et al. (2013) L. A. Ponomarenko, R. V. Gorbachev, G. L. Yu, D. C. Elias, R. Jalil, A. A. Patel, A. Mishchenko, A. S. Mayorov, C. R. Woods, J. R.Wallbank, et al., Nature 497, 594â597 (2013).
 Dean et al. (2013) C. R. Dean, L.Wang, P. Maher, C. Forsythe, F. Ghahari, Y. Gao, J. Katoch, M. Ishigami, P. Moon, M. Koshino, et al., Nature 497, 598 (2013).
 Hunt et al. (2013) B. Hunt, J. D. SanchezYamagishi, A. F. Young, M. Yankowitz, B. J. LeRoy, K. Watanabe, T. Taniguchi, P. Moon, M. Koshino, P. JarilloHerrero, et al., Science 340, 1427 (2013).
 Yang et al. (2013) W. Yang, G. Chen, Z. Shi, C.C. Liu, L. Zhang, G. Xie, M. Cheng, D. Wang, R. Yang, D. Shi, et al., Nature Mater. 12, 792 (2013).
 Nagaosa et al. (2010) N. Nagaosa, J. Sinova, S. Onoda, A. H. MacDonald, and N. P. Ong, Rev. Mod. Phys. 82, 1539 (2010).
 Chang et al. (2013) C.Z. Chang, J. Zhang, X. Feng, J. Shen, Z. Zhang, M. Guo, K. Li, Y. Ou, P. Wei, L.L. Wang, et al., Science 340, 167 (2013).
 Kato et al. (2004) Y. K. Kato, R. C. Myers, A. C. Gossard, and D. D. Awschalom, Science 306, 1910 (2004).
 König et al. (2007) M. König, S. Wiedmann, C. Brüne, A. Roth, H. Buhmann, L. W. Molenkamp, X.L. Qi, and S.C. Zhang, Science 318, 766â770 (2007).
 Guillemette et al. (2013) J. Guillemette, S. S. Sabri, B. Wu, K. Bennaceur, P. E. Gaskell, M. Savard, P. L. Levesque, F. Mahvash, A. Guermoune, M. Siaj, et al., Phys. Rev. Lett. 110, 176801 (2013).
 Ostrovsky et al. (2006) P. Ostrovsky, I. Gornyi, and A. Mirlin, Phys. Rev. B 74, 235443 (2006).
 Tworzydło et al. (2006) J. Tworzydło, B. Trauzettel, M. Titov, A. Rycerz, and C. Beenakker, Phys. Rev. Lett. 96, 246802 (2006).
 Morpurgo and Guinea (2006) A. Morpurgo and F. Guinea, Phys. Rev. Lett. 97, 196804 (2006).
 Peres et al. (2006) N. Peres, F. Guinea, and A. Castro Neto, Phys. Rev. B 73, 125411 (2006).
 Gattenlöhner et al. (2014) S. Gattenlöhner, W.R. Hannes, P. M. Ostrovsky, I. V. Gornyi, A. D. Mirlin, and M. Titov, Phys. Rev. Lett. 112, 026802 (2014).
 Leconte et al. (2014) N. Leconte, F. Ortmann, A. Cresti, J.C. Charlier, and S. Roche, 2D Mater. 1, 021001 (2014).
 Sheng et al. (2006) D. N. Sheng, L. Sheng, and Z. Y. Weng, Phys. Rev. B 73, 233406 (2006).
 Goswami et al. (2008) P. Goswami, X. Jia, and S. Chakravarty, Phys. Rev. B 78, 245406 (2008).
 Ortmann and Roche (2013) F. Ortmann and S. Roche, Phys. Rev. Lett. 110, 086602 (2013).
 J. H. Garcia (2014) T. G. R. J. H. Garcia, L. Covaci, arXiv:1410.8140 (2014).
 Ishii et al. (2011) H. Ishii, N. Kobayashi, and K. Hirose, Phys. Rev. B 83, 233403 (2011).
 Roche and Mayou (1997) S. Roche and D. Mayou, Phys. Rev. Lett. 79, 2518 (1997).

Foa Torres et al. (2014)
L. E. F. Foa Torres,
S. Roche, and
J. C. Charlier,
Introduction to GrapheneBased Nanomaterials: From
Electronic Structure to Quantum Transport (Cambridge University Press, Cambridge, 2014).  Roche et al. (2012) S. Roche, N. Leconte, F. Ortmann, A. Lherbier, D. Soriano, and J.C. Charlier, Solid State Communications 152, 1404 (2012).
 Kubo (1957) R. Kubo, J. Phys. Soc. Jap. 12, 570 (1957).
 Mahan (2000) G. D. Mahan, ManyParticle Physics (Kluwer Academic Publishers, New York, 2000).
 Aoki (1985) H. Aoki, Phys. Rev. Lett. 55, 1136 (1985).
 Evers and Mirlin (2008) F. Evers and A. D. Mirlin, Rev. Mod. Phys. 80, 1355 (2008).
 Luttinger (1951) J. Luttinger, Phys. Rev. 84, 814 (1951).
 Ortmann et al. (2011) F. Ortmann, A. Cresti, G. Montambaux, and S. Roche, EPL 94, 47006 (2011).
 Cummings et al. (2014) A. W. Cummings, A. Cresti, and S. Roche, Phys. Rev. B 90, 161401(R) (2014).
 MartinezGordillo et al. (2014) R. MartinezGordillo, S. Roche, F. Ortmann, and M. Pruneda, Phys. Rev. B 89, 161401 (2014).