Efficient Linear Scaling Approach for Computing the Kubo Hall Conductivity

We report an order-N approach to compute the Kubo Hall conductivity for disorderd two-dimensional systems reaching tens of millions of orbitals, and realistic values of the applied external magnetic fields (as low as a few Tesla). A time-evolution scheme is employed to evaluate the Hall conductivity $\sigma_{xy}$ using a wavepacket propagation method and a continued fraction expansion for the computation of diagonal and off-diagonal matrix elements of the Green functions. The validity of the method is demonstrated by comparison of results with brute-force diagonalization of the Kubo formula, using (disordered) graphene as system of study. This approach to mesoscopic system sizes is opening an unprecedented perspective for so-called 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 two-dimensional materials, out of reach with classical methods.


I. 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 1 and it experienced an important turn in 1980 when K. von Klitzing observed a quantized version of this effect, the so-called Quantum Hall Effect (QHE) 2 . This had an immediate impact and triggered huge amount of work which led to the definition of a resistance standard. 3 The effect of disorder and related localization phenomenon is central to understand the integer QHE from a bulk perspective, where σ xy plateaus develop by varying the charge density, while the plateau region coincide with a region of vanishing σ xx . 4 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 linear-response theory, or the Kubo formula, which was first applied by Aoki and Ando 5 to the QHE problem 4 . 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 [6][7][8] . 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 half-integer QHE with Hall plateaus emerging at σ xy = 4e 2 h (n + 1/2), with integers n. QHE related research in graphene is a vivid field, studying electron-electron interaction 9,10 , superlattices [11][12][13][14] , graphene functionalization and topological currents. Additional transverse transport phenomena are being studied such as (quantum-) anomaleous Hall effect, 15,16 , spin Hall effect 17 and quantum spin Hall effect. 18 The role of disorder in graphene is also debated and the QHE features, with very robust behavior to large disorder 19 , also present additional peculiarities depending on the symmetry breaking aspects conveyed by defects. [20][21][22][23][24] Recently, we have studied the case of oxidized graphene, and have shown the formation of disorder-induced resonant critical states appearing in the zero-energy Landau level with finite σ xx and suggesting a zero-energy σ xy quantized plateau 25 . 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 σ xy is needed and demand for the development of new methodology.
Indeed, the importance of an efficient computational methodology for the Hall conductivity σ xy 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. 5 . Such formulation allows for advanced analysis 26,27 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 σ xy . This computational limitation makes the analysis of realistic models of disordered samples a priori im-possible. 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 low-magnetic field regimes, thus calling for approaches with improved numerical performance. [28][29][30] In this paper, we present a highly efficient order-N 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 28 . 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. 31,32 for computational details, and Ref. 33 for some illustration on realistic models of disordered graphene).

II. METHODOLOGY
Kubo formalism -The general framework is provided by Kubo's linear response theory 34,35 in which the dc conductivity is given as Starting from Eq. (1), the evaluation of the conductivity σ xy proceeds by introducing an orthonormal manyparticle basis with n|H|n = E n which allows to rewrite Eq. (1) into (2) where η is a small parameter necessary to converge this expression.
In the case of non-interacting electrons one can proceed further towards the single-particle picture and obtain where the Fermi-Dirac distribution is described by f (ε k − ζ) and the chemical potential ζ has been introduced. The single-particle eigenvalues ε k 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 diagonalization-based methods, i.e. it becomes impossible to calculate (and store) all eigenvalues and eigenvectors, so that time-evolution based approaches appear promising. By introducing ∞ −∞ dEδ(E − H) 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 random-phase state |Ψ 1 according to k k|...|k → N s Ψ 1 |...|Ψ 1 . Such states are normalized by using the number of sites N s . We further introduce the projection operator NR j=1 |Ψ j Ψ j |, as well as the time-dependent σ(t) fulfilling σ xy = lim t→∞ σ xy (t) with where are the the matrix elements of the Green function (with z = E + iη) and V s = V /N s 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 t ′ → ∞. 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 current-current correlation function should decay with time difference t. 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 δ > 0 has to be introduced such that the integral 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 5,36 which is identical to Ref. 36 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 time-dependent and time independent parts (κ j (E)) while the sum over j 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 H which allows us to use the continued fraction expansion for the diagonal Green function. From the diagonal elements one can obtain the necessary off-diagonal elements in a second iterative step. We find a simple recursion relation connecting both which reads for n > 1 while the first step in the iteration is given as Secondly, the time evolutions of the quantities |Ψ 1 and |j y Ψ j can be performed efficiently through a Chebychev expansion of the time-evolution operator U (t) 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 (z − H) −1 . 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 j is restricted to a cutoff value N R . We use δ = 0.002 if not stated otherwise throughout the paper.
System description The electronic properties of graphene can be described by the standard orthogonal

tight-binding model
where γ 0 = 2.7eV is the nearest neighbor transfer integral and V α the on-site energy, which can be chosen to describe various disorder models. Spatially uncorrelated Anderson disorder is introduced through on-site energies taken at random from [−W γ 0 /2, W γ 0 /2] where W gives the disorder strength. This is a commonly used disorder model for exploring the metal-insulator transition in low-dimensional systems 26,37 , which, in average, does not break electron-hole 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 V α = +(−)V AB 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 38 added to γ 0 which determines the magnetic flux to φ = A · dl = h/e hexagon ϕ αβ per plaquette. 39 Spin degeneracy is assumed throughout the paper and an anti-symmetrization procedure for σ xy (E) has been performed consistent with electron-hole 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 κ j (E) 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 j 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 (N R ), while increasing disorder requires N R 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 time-evolution 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 diagonalization-based 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 time-evolution 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 Hall-conductance steps occur at the expected positions and with the expected height corresponding to the half-integer sequence σ xy = ±(1/2 + n)4e 2 /h.
Next, we consider the interesting case of moderately high magnetic fields such as frequently observed in current experiments [10][11][12][13][14] 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 σ xy = ±10e 2 /h). This consequence of the considered energy resolution (here η = 0.00025) 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 W . The main-frame 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 963 T in the upper inset. The Landau levels are broadened with increasing W around the critical states, as expected from perturbation theory. Simultaneously, the first and second plateaus remain at ±2e 2 /h and ±6e 2 /h indicating robust QHE. In the lower inset, a more realistic magnetic field is chosen (45 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 W = 1. For such disorder, only the zero-energy Landau level fully develops, while all states become localized for E > E 1 .
As another challenging example, we present the analysis of a weakly correlated potential. We chose a globally sublattice-symmetry breaking potential that is locally disordered. In contrast to the above Anderson disorder model, which shifts all on-site 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 V AB = 0.2γ 0 with p = 2.5% of sites affected. Besides ordinary conductance quantization at the values ±2e 2 /h, ±6e 2 /h, ..., an additional plateau is induced at zero energy with zero Hall conductivity. 42 The plateau is clearly visible for fields between 9 T and 45 T. We find the width of the zero-energy plateau equal for all calculated magnetic fields. The corresponding energy value of the onset of the plateau at σ = 0 is determined by the strength of the potential which is given as pV AB = 0.005γ 0 (indicated by dashed vertical line). The transition to higher plateaus (i.e. from ±2e 2 /h to ±6e 2 /h) is rather influenced by the magnetic field-dependence 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 plateau-length scaling with p is evident from the figure. For the considered energy resolution in the inset (η = 0.0004), a trace of plateau onset is still visible for the lowest percentages of the impurity potential down to 0.5% (corresponding to pV AB = 0.001γ 0 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 graphene-based systems. Such approach should become a useful computational tool to simulate QHE in very large size complex disordered materials such as polycrystalline graphene 40 , graphene subjected to weak van der Waals interaction such as by a boron-nitride substrate 41 , or other types of disordered two-dimensional materials. It should also allow to corroborate the formation of zero-energy plateaus of σ xy driven by disorder-induced critical states, and disconnected from degeneracy lifting of Landau levels 25 , an issue of genuine fundamental interest in the context of topological interpretation of the quantized conductance.