Criticality of two-dimensional disordered Dirac fermions in the unitary class and universality of the integer quantum Hall transition
Björn Sbierski, Elizabeth J. Dresselhaus, Joel E. Moore, Ilya A. Gruzberg
CCriticality of two-dimensional disordered Dirac fermions in the unitary class anduniversality of the integer quantum Hall transition
Bj¨orn Sbierski, Elizabeth J. Dresselhaus, Joel E. Moore, and Ilya A. Gruzberg Department of Physics, University of California, Berkeley, California 94720, USA Ohio State University, Department of Physics, 191 West Woodruff Ave, Columbus OH, 43210 (Dated: August 12, 2020)Two-dimensional (2D) Dirac fermions are a central paradigm of modern condensed matter physics,describing low-energy excitations in graphene, in certain classes of superconductors, and on surfacesof 3D topological insulators. At zero energy E = 0, Dirac fermions with mass m are band insulators,with the Chern number jumping by unity at m = 0. This observation motivated Ludwig et al. [Phys.Rev. B , 7526 (1994)] to suggest a relation between 2D disordered Dirac fermions (DDF) andthe integer quantum Hall transition (IQHT). They conjectured that the transitions in both systemsare controlled by the same fixed point and possess the same universal critical properties. Given thefar reaching implications for modern condensed matter physics and our understanding of disorderedcritical points in general, it is surprising that the above conjecture has never been tested numerically.Here, we report the results of extensive numerics on criticality and energy-mass phase diagram of2D-DDF in the unitary symmetry class. We find a critical line at m = 0, with energy dependentlocalization length exponent. At large energies, our results for the DDF are consistent with state-of-the-art numerical results ν IQH = 2 . .
62 from models of the IQHT. At E = 0 however, weobtain ν = 2 . . incompatible with ν IQH . Our result has practical importance for a variety ofexperimental systems but also challenges conjectured relations between models of the IQHT.
Introduction.
The integer quantum Hall effect appearswhen a two-dimensional (2D) electron gas is placed ina strong perpendicular magnetic field. Without disor-der, the electron eigenstates form Landau levels, and eachfilled level contributes unity to the total Chern number C . Disorder is essential for experimental observation ofthe (dimensionless) quantized Hall conductivity σ xy = C .Disorder broadens the Landau levels into bands and lo-calizes eigenstates on a scale ξ ( E ) that diverges in apower law fashion at a critical energy E c [1]: ξ ( E ) ∼ | E − E c | − ν IQH . (1)For Fermi energies E (cid:54) = E c and system sizes L (cid:29) ξ ( E )the Hall conductivity is quantized. The integer quan-tum Hall transition (IQHT) at E = E c is the moststudied Anderson transition [2] because of its conceptualsimplicity, low dimensionality, and relevance to exper-iments. However, critical properties at the IQHT arenotoriously difficult to compute analytically; they aremostly known from numerical studies which employedthe Chalker-Coddington (CC) network model [3–13], mi-croscopic continuous [14, 15], lattice [10, 14–17], andFloquet Hamiltonians [18]. In recent works, the criti-cal properties agree between models, indicating univer-sality of the IQHT. They include the localization lengthexponent ν IQH = 2 . .
62 and the leading irrelevant ex-ponent y (cid:39) . y describes the approach of the dimensionless quasi-1DLyapunov exponent Γ to its limiting value at infinite sys-tem size Γ IQH0 = 0 . .
82 [5–7, 9, 11–13, 16]. A similarexponent y was found for the average conductance g ofa square sample with limiting value g IQH = 0 . . (a) cleanbandinsulator (b) disorderedbandinsulator (Anderson) ins.(Anderson) ins.metal metal QHP critical (delocalized) FIG. 1. Schematic phase diagram for massive 2D Diracfermions. In the clean case (a), a metal intervenes be-tween two band insulators with different Chern numbers C at | m | > | E | . With disorder in the unitary class (b), thecritical line m = 0 separates Anderson insulators. The local-ization length exponent ν at E = 0 is significantly differentfrom the one at E > E = 0 , . , . [21–23] and the discussion below. The IQHT has alsobeen discussed recently in the context of exotic topolog-ical superconductor surface states [24].A longstanding conjecture by Ludwig et al. [25] statesthat the IQHT fixed point also controls the criticalityof 2D disordered Dirac fermions (DDF) in the unitarysymmetry class [26, 27]. The clean Dirac Hamiltonian is H = (cid:126) v ( − iσ x ∂ x − iσ y ∂ y ) + mσ z , (2)with 2 × σ µ , mass m , and velocity v .The spectrum of H has a gap of size 2 | m | symmetricaround E = 0. For Fermi energies E within the gap, thesystem is a band insulator with a half-integer quantized σ xy = C ( m ) = − sgn( m ) [25], see Fig. 1(a). If the a r X i v : . [ c ond - m a t . m e s - h a ll ] A ug Dirac fermion is regularized on a lattice as in the Haldanemodel [28] or Eq. (7) below, the Hamiltonian (2) onlydescribes the low-energy excitations near a certain pointin the Brillouin zone. Bloch states elsewhere contributeanother 1/2 to C , such that | σ xy | jumps between zeroand one as m changes sign.With m taking the role of energy, the superficial sim-ilarity of this transition to the IQHT motivated Ludwiget al. [25] to study effects of disorder on the Dirac Hamil-tonian, placing the model in the unitary symmetry class, H = H + (cid:88) µ =0 ,x,y,z U µ ( x, y ) σ µ . (3)The random scalar ( U ) and vector ( U x,y ) potentials,and the random part of the mass ( U z ) are taken tobe independent Gaussian fields with the correlators U µ ( r ) U ν ( r (cid:48) ) = δ µν K µ ( | r (cid:48) − r | ) and zero mean. We ex-pect that for m (cid:54) = 0, the eigenstates or H are localizedwith the localization length ξ ( m ) ∼ | m | − ν E , c.f. Eq. (1),with a possibly E -dependent critical exponent.Although model (3) is not solvable analytically, theconjecture ν E = ν IQH is consistent with several plau-sible arguments which we discuss below. However, de-spite the importance of 2D Dirac Hamiltonians in moderncondensed matter physics, the conjectured emergence ofIQHT criticality in DDF was never checked numerically.In what follows, we address this issue with extensive nu-merical simulations employing different microscopic mod-els and scaling observables. We start with the continuummodel (3) and use the transfer matrix (TM) approach inquasi-1D (q1D) geometry to find the critical behaviornear the line m = 0 in the m - E plane, see Fig. 1(b). Atlarge E our results are consistent with ν E = ν IQH , butthe critical exponent at zero energy ν E =0 = 2 . .
36 isclose to, but strikingly incompatible with ν IQH . We cor-roborate our results in a lattice model of DDF, employingan alternative 2D scaling observable [29].
Continuum model and disorder induced length scale.
We start with Hamiltonian (3) at E = 0 and smoothdisorder with correlation length a , K µ ( r ) = W e − r / a / π. (4)We use a and (cid:126) v/a as units of length and energy. Thenthe disorder strength W , taken to be the same for allfour disorder fields, is dimensionless and is the bare en-ergy scale in the model. The mean free path l W equalsthe quasiparticle decay time, l W ≡ − / Im Σ ↑↑ (0 ,
0) de-fined in terms of the disorder-averaged Green function G ( k , ω ) = [ ω − H ( k ) − Σ( k , ω )] − . For weak dis-order W (cid:28)
1, a perturbative renormalization group(RG) [25, 30] gives, for m = 0, l W ∝ e c/W , with c = O (1). To ensure that our system sizes L (cid:29) l W ,we work with strong disorder W ≥ . l W =1 . = 1 .
54. Wealso observe that for kl W >
1, the peaks in the spectral function A ( k , ω ) = − π tr Im G ( k , ω ) occur at frequencies ω (cid:39) ± (cid:126) vk , i.e. the velocity v is almost un-renormalized.We conclude that for W = 1 .
5, system sizes L (cid:38) O (10)are large enough to exhibit disorder-dominated physics. Lyapunov exponent (LE).
A common method to ana-lyze critical behavior in disordered systems employs theself-averaging of the LEs γ i of q1D samples of width L y and length L x → ∞ [32]. The smallest γ i > γL y , which increases (decreases) with L y in a local-ized (extended) phase and is scale-invariant at a criticalpoint. Following Ref. [33], we use finite L x = O (10 ) andfind Γ as the average over (cid:63)
200 disorder realizations, seethe supplemental material (SM) for details.The eigenvalue problem for the DDF (3) can be rewrit-ten as ∂ x ψ ( x, k y ) = f ( ψ ( x, k (cid:48) y )). The right hand sidecontains scattering between transversal wavevectors k y but is local in x , which allows us to express the TM inexponential form. We impose periodic boundary condi-tions (BC) in the y direction. We discretize the x di-rection and stabilize the TM multiplication by repeatedQR-decompositions [1] (to obtain Γ) or in a scatteringmatrix formalism [34] (for the conductance of moderatelysized systems). Both methods are numerically exact andfaithfully treat a single Dirac node without band bend-ing or node doubling, and maintain the statistical sym-metry m → − m . The only approximations are relatedto the cutoff | k y | ≤ k max and the x -discretization. Theassociated length scales (taken equal) were chosen muchsmaller than a , and the results are converged with respectto these parameters.Results for the dimensionless LE Γ at E = 0, W = 1 . m and system widths L y are presented inFig. 2. The solid lines are fits to the scaling functionΓ( m, L y ) = Γ + α L − yy + α m L /νy , (5)which is the lowest-order polynomial ansatz allowed bysymmetry, including an irrelevant contribution. The fitgives the following critical properties: ν E =0 = 2 . , y = 0 . , Γ = 0 . , (6)the number in parentheses denotes one standard devia-tion. In the SM, we give a detailed account for the fittingprocedure and show its stability with respect to higherorder terms in Eq. (5) and a removal of data pointsof large m and small L y . There, we also present datafor an increased disorder strength W = 2 .
0, which yields ν E =0 = 2 . y = 0 . = 0 . Lattice model and alternative scaling observable.
Wenow confirm the value of ν E =0 using a square-lattice reg-ularization of the DDF allowing access to an alternativescaling observable introduced by Fulga et al. [35]. Inmomentum space, the clean model reads [36] H L = σ x sin k x + σ y sin k y + σ z ( m − k x + cos k y ) , (7) m = L y = 2.32 ± 0.01 L y ( m = ) y = 0.51 ± 0.03 = 0.84 ± 0.01 FIG. 2. Top: LEs Γ for E = 0 and W = 1 . m . System widths are L y = 40 , , , , , , , ≤ . m = 0) with extrapolation toinfinite system size determining Γ (cross). where the lattice spacing and the energy scale have beenset to unity. For | k | (cid:28)
1, this model reduces to theDirac Hamiltonian (2), with a topological transition at m = 0 where the Chern number C changes by 1. Weadd on-site disorder potentials, V = (cid:80) r i ,µ U µ ( r i ) σ µ with U µ ( r i ) uniformly drawn from the interval [ − w/ , w/ r i and µ = 0 , x, y, z .Transport calculations use the kwant software package[37] and employ two extended leads in the x -direction,described by the Hamiltonian H Lead ( k x , k y ) = σ x sin k x + σ z (1 + cos k x ) . (8)The lattice model (7) has no symmetry that ensures m c = 0 in the presence of disorder. However, the Diracnode energy is not renormalized away from E = 0 due tothe sum of non-trivial Pauli matrices in Eq. (7). Bandbending effects are important for momenta and energiesabove unity. Transport simulations of a 200 ×
200 squaresystem with periodic BC in the y direction and disorderstrength w = 2 . (cid:39) . m (cid:39) . y direction shows a smooth transition from 0to 1 within the same mass range, revealing the integrityof chiral edge states. This shows that the localizationlengths for w = 2 . m . Thus,we use w = 2 . ν E =0 , we consider the re- m = l o g | z | = 2.33 ± 0.03 m c = 0.0911 ± 0.0002 FIG. 3. Scaling plot of the variable Λ for the lattice model atenergy E = 0 and disorder strength w = 2 .
5. Dots representaverages over at least 10 disorder realizations for system sizes L = 60 , , , ,
200 and the solid curves are results of thefitting procedure described in the SM. flection matrix r ( φ ) of the left lead as a function of thephase φ of twisted BC in the y direction. For a given dis-order realization, the critical value m = m c occurs whenthere exist a φ such that r ( φ ) has a zero eigenvalue anddet r ( φ ) = 0. Fulga et al. [35] showed that a scaling ob-servable Λ can be obtained by working with generalizedtwisted BC ψ x,y = L − = z ψ x,y =0 for all x = 0 , , ..., L − z ∈ C . Now, det r ( z ) has zeros z even for m (cid:54) = m c but with | z | (cid:54) = 1. Then for the z closest to the unitcircle we define Λ = log | z | as a measure of distance tocriticality Λ = 0. For the CC model, Fulga et al. demon-strated scaling of Λ with system size L , giving ν = 2 . H L + V for m around 0 . L = 60 and 200, seeFig. 3 for the results and the SM for details of the fit.We find ν E =0 = 2 . w = 2 .
25 and 2 .
75 (not shown)yields compatible ν within the given error bars. Results for finite energy (
E > ). We now consider thecontinuum model (3) with smooth disorder (4) at finiteenergy
E >
E < E = 0 . E = 0 . W = 2. As for the E = 0 case we find localizing behaviorfor any m (cid:54) = 0. The exponent ν E =0 . = 2 . ν E =0 . However, ν E =0 . = 2 . ν E =0 and is much closer to ν IQH . Othercritical properties (Γ and y ) for E > E = 0 results.To further probe the critical line m = 0, we computedthe critical distributions of the Landauer conductance g of L × L systems with periodic BC in the y direction,and metallic leads modeled as highly doped Dirac nodes[38], see SM for details at E = 0. Such distributionsand their moments are expected to be scale-invariant anduniversal [2, 4]. In Fig. 4 we present results for the E g L = 100, W = 1.5 L = 200, W = 1.5 L = 100, W = 2.5 L = 200, W = 2.5 FIG. 4. Landauer conductance g of square samples at van-ishing mass m = 0, size L = 100 ,
200 and periodic BC intransversal direction averaged over at least 10000 disorder re-alizations. The disorder strengths are W = 1 . . average g . We observe that for E (cid:46) . g (cid:39) . E , which weinterpret as evidence of proximity to the underlying fixedpoint. With increasing L , g slightly increases, consistentwith decreasing Γ( m = 0) in Fig. 2 (bottom).For 0 . (cid:46) E , g begins to depend on W , and varies with E by ∼
50% for W = 1 . ∼
10% for W = 2 . W = 1 . E > . g slightly decreases when L grows form 100 to 200. We interpret this as a remnant ofthe crossover from the diffusive to the critical behavior.It is consistent that LEs obtained in this regime (notshown) cease to obey critical scaling. Discussion.
In summary, our numerical results forDDF are consistent with localized behavior anywhere inthe mass-energy plane except on a critical line m = 0,see Fig. 1(b). At m = 0, both the dimensionless LE ex-trapolated to infinite system size Γ = 0 . .
85 and theirrelevant exponent y do not vary significantly with en-ergy or disorder strength below E (cid:39)
1, while the averageconductance g of fixed-size square samples at strongerdisorder varies at most by ∼ ν E sig-nificantly depends on energy. While ν E =0 . = 2 . ν IQH = 2 . .
62, the value ν E =0 is significantlysmaller. Taking a union over error bars for the two mod-els and two scaling methods that we used, we obtain thefollowing conservative estimate: ν E =0 = 2 . . . (9)The value ν E =0 . = 2 . E and low W characterized by a large Drude conductivity σ D xx (cid:29) σ D xx controls thederivation of an effective field theory for the DDF withshort-range disorder [39] as it justifies the required sad-dle point approximation. The resulting non-linear sigmamodel with θ -term can also be derived for other models of the IQHT: the Schr¨odinger equation with short-rangedisorder and strong magnetic field [40, 41], and the CCnetwork model [42, 43]. These relations rationalize ourfinding of IQHT-like criticality in the DDF at E = 0 . σ D xx , and the derivation of the sigmamodel for it is uncontrolled, as well as for the DDF at E (cid:39)
0, where σ xx < E = 0 and discuss possible reasonsfor the inequality ν E =0 (cid:54) = ν IQH in two scenarios. (a) Insufficient system size.
Given the long and tortu-ous history of numerical work on the IQHT, where refinedfitting functions and the ability to study larger systemsshifted the value of ν considerably over time, we cannotexclude the scenario that ν E =0 in Eq. (9) is not the trueasymptotic value, and further increase in system sizeswould lead to a crossover to ν IQH . Indeed, taking thedifference Γ( m = 0 , L y, max ) − Γ as a proxy for the dis-tance to the critical point, its value at E = 0 is roughlythree times larger than at E = 0 . ν E =0 if the minimal L y involved in the fitis increased from 40 to 68, see SM. Finally, we corrobo-rated our result (9) at a different disorder strength andusing an alternative scaling observable in a lattice modelregularization of the DDF. Our finding for ν E =0 is alsosupported by numerical results on a massless DDF in amagnetic field [44]. At strong enough potential disorder,only the critical state deriving from the Landau level at E = 0 persists, separating localized states at E ≶
0. Thescaling of dσ xy /dE | E =0 and the width of the conductancepeak around E = 0 with system size L gave ν ≈ .
3, butno error bars were provided. (b) Different fixed points.
In a more intriguing scenarioour results can be explained if we conjecture the exis-tence of two different fixed points. One of them is theconventional IQHT fixed point that controls the criticalbehavior of DDF at
E >
0, while the other fixed pointcontrols the system at m = 0 , E = 0. We conjecturethat this fixed point is multicritical, where both m and E are relevant, with the RG eigenvalues y m = 1 /ν E =0 and y E . The RG flow near this point would resemblethat near the tricritical point in the Ising model with va-cancies [45] or near the point σ xx = 0, σ xy = 1 / E > E = 0 . E is marginally relevant) value ofthe crossover exponent φ = y E /y m at the multicriticalpoint, resulting in the cusp-like shape of the crossoverline in Fig. 1(b). However, this scenario cannot explainthe apparent energy independence of Γ .Finally, we review the standard (non-rigorous) argu-ments for ν E =0 = ν IQH , and identify their possible flaws.(i) The original argument of Ludwig et al. [25] assumessufficiently smooth disorder in the DDF model. Then thelow energy states are chiral Jackiw-Rebbi fermions mov-ing along lines of zero mass m + U z ( x, y ) = 0. Otherdisorders lead to random phases accumulated betweensaddle points in U z , where scattering controlled by thevalue of m occurs. This argument parallels the one thatleads to the CC network model [49]. A possible flaw isthat unlike in the large B -field case of the IQHT, non-chiral higher energy states might be important in DDF.Indeed, solving the Dirac equation for a linear domainwall U z ( x ) = cx arising from our smooth disorder poten-tial with c ∼ W/a , the first pair of counterpropagatingedge modes appears at E ∼ ±√ W . Naively, it takes adisorder strength W (cid:63) U fluctuations, but their relevance for the network modelis presently unclear. In addition, the effects of geometricdisorder might further modify the above argument, seethe discussion below.(ii) Another argument in favor for the equivalence ofthe CC network and DDF was given by Ho and Chalker[50]. The authors showed that the spectrum of quasiener-gies of the clean CC model (viewed as a Floquet system)has a Dirac point with mass m proportional to the de-viation from the critical point. Then they argued thatdisorder in the CC model leads to all possible types ofdisorder in the Dirac model. However, the mapping as-sumes smooth disorder in the network model, and maybe invalid for strong disorder.In either case, our results raise the issue of universal-ity of the IQHT, which was also challenged by two recentresults. Refs. [51, 52] numerically studied the role ofstructural (geometric) disorder in the CC model and re-ported ν ≈ . ν seems to indicate therelevance of geometric disorder, and is consistenly repro-duced by the alternative scaling method of Ref. [35] inongoing work [53]. In a different development, Zirnbauer[23] proposed a solvable conformal field theory for theIQHT, predicting ν = ∞ and y = 0. If true, this wouldimply apparent critical exponents that show logarithmicflow with L , invalidate all existing numerical studies ofthe exponent ν at the IQHT, and exclude any possibilityto find universal critical behavior numerically. Outlook.
We hope that our findings will prompt a care-ful re-examination of criticality at the IQHT and otherAnderson transitions. For future work on critical DDF, itis interesting to numerically study the multifractal prop-erties of wavefunctions and compare them to establishedresults for the IQHT [2]. Further, extension of our meth-ods to DDF in the symmetry classes of the spin and ther-mal quantum Hall effects is worthwhile.
Acknowledgements.
We acknowledge useful discussionswith Jens Bardarson, Matt Foster, Cosma Fulga, IgorGornyi, Alexander Mirlin and Pavel Ostrovsky. Com-putations were performed at the Ohio Supercomputer Center and the Lawrencium cluster at Lawrence Berke-ley National Lab. B.S. acknowledges financial supportby the German National Academy of Sciences Leopold-ina through grant LPDS 2018-12. E.J.D was supportedby NSF Graduate Research Fellowship Program, NSFDGE 1752814. J.E.M. acknowledges support by TIMESat Lawrence Berkeley National Laboratory supported bythe U.S. Department of Energy, Office of Basic EnergySciences, Division of Materials Sciences and Engineering,under Contract No. DE-AC02-76SF00515 and a SimonsInvestigatorship. [1] B. Huckestein.
Scaling theory of the integer quantum Halleffect . Rev. Mod. Phys., , 357 (1995).[2] F. Evers and A. D. Mirlin. Anderson transitions . Rev.Mod. Phys., , 1355 (2008).[3] J. T. Chalker and P. D. Coddington. Percolation, quan-tum tunnelling and the integer Hall effect . J. Phys. C, , 2665–2679 (1988).[4] B. Kramer, T. Ohtsuki, and S. Kettemann. Random net-work models and quantum phase transitions in two di-mensions . Phys. Rep., , 211–342 (2005).[5] H. Obuse, A. R. Subramaniam, A. Furusaki, I. A.Gruzberg, and A. W. W. Ludwig.
Boundary Multifrac-tality at the Integer Quantum Hall Plateau Transition:Implications for the Critical Theory . Phys. Rev. Lett., , 116802 (2008).[6] F. Evers, A. Mildenberger, and A. D. Mirlin.
Multi-fractality at the Quantum Hall Transition: Beyond theParabolic Paradigm . Phys. Rev. Lett., , 116803(2008).[7] K. Slevin and T. Ohtsuki.
Critical exponent for the quan-tum Hall transition . Phys. Rev. B, , 041304 (2009).[8] H. Obuse, A. R. Subramaniam, A. Furusaki, I. A.Gruzberg, and A. W. W. Ludwig. Conformal invari-ance, multifractality, and finite-size scaling at Andersonlocalization transitions in two dimensions . Phys. Rev. B, , 035309 (2010).[9] M. Amado, A. V. Malyshev, A. Sedrakyan, andF. Dom´ınguez-Adame. Numerical study of the localiza-tion length critical index in a network model of plateau-plateau transitions in the quantum Hall effect . Phys. Rev.Lett., , 066402 (2011).[10] I. C. Fulga, F. Hassler, A. R. Akhmerov, and C. W. J.Beenakker.
Topological quantum number and critical ex-ponent from conductance fluctuations at the quantumHall plateau transition . Phys. Rev. B, , 245447 (2011).[11] H. Obuse, I. A. Gruzberg, and F. Evers. Finite-size ef-fects and irrelevant corrections to scaling near the integerquantum Hall transition . Phys. Rev. Lett., , 206804(2012).[12] K. Slevin and T. Ohtsuki.
Finite Size Scaling of theChalker-Coddington Model . Int. J. Mod. Phys. Conf. Ser., , 60–69 (2012).[13] W. Nuding, A. Kl¨umper, and A. Sedrakyan. Localiza-tion length index and subleading corrections in a Chalker-Coddington model: A numerical study . Phys. Rev. B, ,115107 (2015).[14] M. Ippoliti, S. D. Geraedts, and R. N. Bhatt. Integer quantum Hall transition in a fraction of a Landau level .Phys. Rev. B, , 014205 (2018).[15] Q. Zhu, P. Wu, R. N. Bhatt, and X. Wan. Localization-length exponent in two models of quantum Hall plateautransitions . Phys. Rev. B, , 024205 (2019).[16] M. Puschmann, P. Cain, M. Schreiber, and T. Vojta. Integer quantum Hall transition on a tight-binding lattice .Phys. Rev. B, , 121301 (2019).[17] M. Puschmann, P. Cain, M. Schreiber, and T. Vojta. Edge state critical behavior of the integer quantum Halltransition . arXiv e-prints, arXiv:2004.01611 (2020).[18] J. P. Dahlhaus, J. M. Edge, J. Tworzyd(cid:32)lo, and C. W. J.Beenakker.
Quantum Hall effect in a one-dimensionaldynamical system . Phys. Rev. B, , 115133 (2011).[19] Z. Wang, B. Jovanovi´c, and D.-H. Lee. Critical Conduc-tance and Its Fluctuations at Integer Hall Plateau Tran-sitions . Phys. Rev. Lett., , 4426 (1996).[20] L. Schweitzer and P. Markoˇs. Universal conductance andconductivity at critical points in integer quantum hall sys-tems . Phys. Rev. Lett., , 256805 (2005).[21] R. Bondesan, D. Wieczorek, and M. R. Zirnbauer. Purescaling operators at the integer quantum Hall plateautransition . Phys. Rev. Lett., , 186803 (2014).[22] R. Bondesan, D. Wieczorek, and M. R. Zirnbauer.
Gaus-sian free fields at the integer quantum Hall plateau tran-sition . Nuclear Physics B, , 52–90 (2017).[23] M. R. Zirnbauer.
The integer quantum Hall plateau tran-sition is a current algebra after all . Nucl. Phys. B, ,458 (2019).[24] B. Sbierski, J. F. Karcher, and M. S. Foster.
Spectrum-Wide Quantum Criticality at the Surface of Class AIIITopological Phases: An “Energy Stack” of Integer Quan-tum Hall Plateau Transitions . Phys. Rev. X, , 021025(2020).[25] A. W. W. Ludwig, M. P. A. Fisher, R. Shankar, andG. Grinstein. Integer quantum Hall transition: An al-ternative approach and exact results . Phys. Rev. B, ,7526 (1994).[26] M. R. Zirnbauer. Riemannian symmetric superspaces andtheir origin in random-matrix theory . J. Math. Phys., ,4986 (1996).[27] A. Altland and M. R. Zirnbauer. Nonstandard symme-try classes in mesoscopic normal-superconducting hybridstructures . Phys. Rev. B, , 1142 (1997).[28] F. Haldane. Model for a quantum Hall effect withoutLandau levels: Condensed-matter realization of the ”par-ity anomaly” . Phys. Rev. Lett., , 2015 (1988).[29] I. Fulga, F. Hassler, A. Akhmerov, and C. Beenakker. Topological quantum number and critical exponent fromconductance fluctuations at the quantum Hall plateautransition . Phys. Rev. B, , 245447 (2011).[30] A. Schuessler, P. Ostrovsky, I. Gornyi, and A. Mirlin. An-alytic theory of ballistic transport in disordered graphene .Phys. Rev. B, , 075405 (2009).[31] B. Sbierski and C. Fr¨aßdorf. Strong disorder in nodalsemimetals: Schwinger-Dyson-Ward approach . Phys.Rev. B, , 020201 (2019).[32] B. Kramer and A. McKinnon. Localization: theory andexperiment . Rep. Prog. Phys., , 1469 (1993).[33] M. Amado, A. V. Malyshev, A. Sedrakyan, andF. Dom´ınguez-Adame. Numerical study of the localiza-tion length critical index in a network model of plateau-plateau transitions in the quantum hall effect . Phys. Rev.Lett., , 066402 (2011). [34] J. Bardarson, J. Tworzydlo, P. Brouwer, andC. Beenakker.
One-Parameter Scaling at the DiracPoint in Graphene . Phys. Rev. Lett., , 106801 (2007).[35] I. C. Fulga, F. Hassler, a. R. Akhmerov, and C. W. J.Beenakker. Scattering formula for the topological quan-tum number of a disordered multimode wire . Phys. Rev.B, , 155429 (2011).[36] X.-l. Qi, Y.-s. Wu, and S.-c. Zhang. Topological quantiza-tion of the spin Hall effect in two-dimensional paramag-netic semiconductors . Phys. Rev. B, , 085308 (2006).[37] C. Groth, M. Wimmer, A. Akhmerov, and X. Waintal. Kwant: a software package for quantum transport . NewJ. Phys., , 063065 (2014).[38] J. Tworzydlo, B. Trauzettel, M. Titov, A. Rycerz, andC. W. J. Beenakker. Quantum-limited shot noise ingraphene . Phys. Rev. Lett., , 246802 (2006).[39] P. M. Ostrovsky, I. V. Gornyi, and A. D. Mirlin. Quan-tum criticality and minimal conductivity in graphene withlong-range disorder . Phys. Rev. Lett., , 256801 (2007).[40] A. M. M. Pruisken. On localization in the theory of thequantized hall effect: A two-dimensional realization of the θ -vacuum . Nucl. Phys. B., , 277 (1984).[41] H. Weidenm¨uller. Single electron in a random potentialand a strong magnetic field . Nucl. Phys. B., , 87(1987).[42] N. Read (1991). Unpublished.[43] M. R. Zirnbauer.
Towards a theory of the integer quan-tum Hall transition: From the nonlinear sigma model tosuperspin chains . Annalen Phys., , 513–577 (1994). [Er-ratum: Annalen Phys. , 89 (1995)].[44] K. Nomura, S. Ryu, M. Koshino, C. Mudry, and A. Fu-rusaki. Quantum Hall Effect of Massless Dirac Fermionsin a Vanishing Magnetic Field . Phys. Rev. Lett., ,246806 (2008).[45] J. Cardy.
Scaling and Renormalization in StatisticalPhysics (1996).[46] D. E. Khmel’Nitskiˇi.
Quantization of Hall conductivity .Soviet Journal of Experimental and Theoretical PhysicsLetters, , 552–556 (1983).[47] D. E. Khmelnitskii. Quantum hall effect and addi-tional oscillations of conductivity in weak magnetic fields .Physics Letters A, , 182–183 (1984).[48] A. M. M. Pruisken.
Dilute instanton gas as the precursorto the integral quantum Hall effect . Phys. Rev. B, ,2636–2639 (1985).[49] J. T. Chalker and P. D. Coddington. Percolation, quan-tum tunnelling and the integer Hall effect . J. Phys. CSolid State Phys., , 2665 (1988).[50] C. Ho and J. Chalker. Models for the integer quantumHall effect: The network model, the Dirac equation, anda tight-binding Hamiltonian . Phys. Rev. B., , 8708(1996).[51] I. A. Gruzberg, W. Nuding, and A. Sedrakyan. Geometri-cally disordered network models, quenched quantum grav-ity, and critical behavior at quantum Hall plateau transi-tions . Phys. Rev. B, , 125414 (2017).[52] A. Kl¨umper, W. Nuding, and A. Sedrakyan. Randomnetwork models with variable disorder of geometry . Phys.Rev. B, , 140201 (2019).[53] E. J. Dresselhaus, B. Sbierski, I. A. Gruzberg, and J. E.Moore.
Manuscript in preparation. (2020).[54] M. Janssen, M. Metzler, and M. Zirnbauer.
Point-contactconductances at the quantum Hall transition . Phys. Rev.B., , 15836 (1999). [55] R. Klesse and M. Zirnbauer. Point-Contact Conduc-tances from Density Correlations . Phys. Rev. Lett., ,2094 (2000).[56] I. A. Gruzberg, A. D. Mirlin, and M. R. Zirnbauer. Classification and symmetry properties of scaling dimen-sions at Anderson transitions . Phys. Rev. B, , 125144(2013).[57] B. Jovanovic and Z. Wang. Conductance Correlationsnear Integer Quantum Hall Transition . Phys. Rev. Lett., , 2767 (1998). Supplemental material for “Criticality of two-dimensional disordered Dirac fermions in the uni-tary class and universality of the integer quantumHall transition”
Quasi-1D Lyapunov exponents: Additional data,fitting procedure
Here, we give a detailed description of the ensem-ble of data sets, fitting procedure and results for thequasi-1D Lyapunov exponents Γ of the disordered Diracfermion. Table I summarizes all data sets used in thisstudy. The first two columns denote energy E and disor-der strength W , respectively. The third column givesthe length L x of the quasi-1D slabs, the widths are L y = 40 , , , , , , ,
160 for all data sets. Weuse a linearly spaced grid of squared Dirac masses m as given in the respective column. The minimal num-ber of disorder realizations per tuple ( m, L y ) is shownin the next column. In Fig. S1 we show a typical his-togram of the ensemble of finite-length Γ approximants,for m = 0 . L y = 102 and E = 0, W = 1 .
5. Fol-lowing Ref. [33], we confirm their gaussian distribution(red line) and approximate Γ as the mean, the error σ Γ is the standard deviation of this ensemble divided by thesquare root of the number of disorder realizations. Themaximal error so obtained in given in the respective col-umn of Table I. For a few representative data points wehave increased L x by a factor of five and checked the Γonly vary within error bars, moreover note the factor ofthree between the two values for L x used for the two datasets at E = 0. The so obtained dimensionless Lyapunovexponents Γ are shown as dots in Fig. 2 of the main text( E = 0, W = 1 .
5) and in Figs. S2-S4, for the remaining(
E, W ) parameter pairs of table I.We now consider the scaling ansatz for fitting the P () FIG. S1. Histogram of finite-length Lyapunov exponents for E = 0, W = 1 . L y = 102 and m = 0 .
05, involving 1000disorder realizations. The red line is a gaussian fit. Γ( m, L y ) data. Following Ref. [12], we take into accountthe relevant and the leading irrelevant scaling variables, x R ( m ) and x I ( m ), respectively:Γ( m, L y ) = F (cid:16) X R = x R ( m ) L /νy , X I = x I ( m ) L − yy (cid:17) . (S1)The relevant variable vanishes at the critical point, x R (0) = 0. To satisfy the (statistical) symmetry m →− m , we follow the customary choice [12, 33] an expand x R ( m ) as an odd-power polynomial of order R , and x I ( m ) as an even-power polynomial of order I . Thenthe series expansion of Eq. (S1) must contain only evenpowers of X R :Γ( m, L y ) = Γ + α X I + α X R + α X I + α X R X I + ... (S2)We fix the expansion orders R, I and the truncation ofthe series in Eq. (S2) and find the parameters ν, y, Γ , α ij using non-linear least squares fitting. The quality of a fitis reported by the value of ˜ χ = χ / ( N − N c ) where N is the number of data points Γ j (j stands for the tuple( m, L y )), N c is the number of fitting parameters in thechosen scaling function, and χ = (cid:80) Nj =1 (Γ j − Γ) /σ j isthe sum of squared residues weighted by the variances σ j .The errors of the fit parameters (one standard deviation)are obtained from repeated fitting of ∼
100 syntheticdata sets generated from a normal distribution of meanΓ j and variance σ j . The fit results, i.e. the optimizeduniversal parameters of the ansatz with the overall low-est ˜ χ and all relative errors of fit parameters below 30%,are given in the last three columns of Table I and alsoin Fig. 2 and S2-S4. It turns out that for all data setsconsidered, the minimal ansatz which allows for an irrel-evant contribution (first three terms on the rhs of Eq. S2with R = 1, I = 0 and N c = 5, given explicitly in Eq. (5)of the main text) is chosen by our fitting algorithm. Asshown in the stability plots of Figs. S5 to S8, while some-times a higher order fitting function can give a slightlysmaller value of ˜ χ , this always comes at the cost of arelative error exceeding our stabilitiy bound above. Asa further stability criterion, we request that with one ortwo small system widths L y or a few largest masses m removed from the fit, critical properties still overlap inerror bars. This is also confirmed in Table I, see columntitled “restriction”. E W L x ( · ) m min realizations max σ Γ / Γ restriction ν y Γ ˜ χ N m ≤ .
005 2.32(1) 0.54(2) 0.85(1) 0.996 72 L y ≥
68 2.32(1) 0.39(9) 0.82(4) 0.978 660.0 2.0 3 0,0.0005,...,0.005 200 0.2% none 2.31(2) 0.51(3) 0.84(1) 0.960 88 m ≤ .
004 2.32(3) 0.51(5) 0.84(1) 1.163 72 L y ≥
68 2.30(3) 0.47(10) 0.84(2) 0.919 660.3 2.0 1 0,0.000625,...,0.005 200 0.25% none 2.34(2) 0.50(5) 0.84(1) 0.803 72 m ≤ . L y ≥
50 2.34(3) 0.42(6) 0.82(2) 0.786 630.7 2.0 1 0,0.00125,...,0.01 550 0.2% none 2.53(2) 0.60(9) 0.83(1) 0.993 72 m ≤ . L y = 40 , , , , , , ,
160 for all data sets unless restricted for a stabilitycheck. m = L y = 2.31 ± 0.02 N = 88, N c = 5 = 0.96 L y ( m = ) y = 0.51 ± 0.03 = 0.84 ± 0.01 FIG. S2. Quasi-1D Lyapunov exponents as in Fig. 2 of themain text, but for E = 0 and W = 2 .
0. Solid lines denote thebest fit with the ansatz in Eq. (2) of the main text.
Details for calculation and fitting of the alternativescaling observable
We provide details here of calculating the alternativescaling observable Λ proposed by Fulga et al. [35]. Fortight-binding models, the concept of virtual leads extend-ing in y -direction is used to apply the generalized twistedtransversal boundary conditions mentioned in the maintext [35]. When extending the method to the DDF lat-tice model, and more generally any model with multipleorbitals per site, it is crucial to match the propagatingmodes in the virtual leads on top and bottom. With the m = L y = 2.34 ± 0.02 N = 72, N c = 5 = 0.803 L y ( m = ) y = 0.5 ± 0.05 = 0.84 ± 0.01 FIG. S3. Quasi-1D Lyapunov exponents as in Fig. 2 of themain text, but for E = 0 . W = 2 .
0. Solid lines denotethe best fit with the ansatz in Eq. (2) of the main text. kwant software package [37] this is not automaticallyensured, but can be corrected through the following pro-cedure.With the
PropagatingModes method of kwant , thepropagating modes in the top and bottom leads can becompared. They are beyond the user’s control and arefound to be different generically. We construct unitarymatrices that rotate the, say, ingoing top lead modesto match the outgoing modes of the bottom lead, ψ inn = (cid:80) k U innk ˜ ψ k and likewise for the outgoing modes of the toplead. Applying these transformations in the Schr¨odingerequation defining the scattering matrix S returned by0 m = L y = 2.53 ± 0.02 N = 72, N c = 5 = 0.993 L y ( m = ) y = 0.6 ± 0.09 = 0.83 ± 0.01 FIG. S4. Quasi-1D Lyapunov exponents as in Fig. 2 of themain text, but for E = 0 . W = 2 .
0. Solid lines denotethe best fit with the ansatz in Eq. (2) of the main text. kwant , we can transform to a scattering matrix ˜ S =( U out ) T S ( U in ) ∗ defined in a matching mode basis for topand bottom leads which can now be safely short-circuitedincluding the z -factor as detailed in Ref. [35].As explained in the main text, we find the z = z wherethe determinant of the reflection matrix r vanishes. InFig. S9 we show the histogram of log | z | for subcriti-cal m = 0 .
05 and three different system sizes. Like inthe case of Lyapunov exponents, they are gaussian dis-tributed and their mean and standard deviation deter-mine the scaling observable Λ reported in Fig. 3 of themain text.The fitting procedure and scaling function Λ( m, L )parallels the discussion for the Lyapunov exponentsΓ( m, L y ) discussed above. However, as there is no sym-metry fixing m c = 0, we have to introduce a fitting pa-rameter m c so that δ = m − m c appears in the expansionof the relevant scaling field and general polynomials areallowed. We find the ansatz Λ( m, L ) = β X R + β X R with X R = L /ν [ δ + αδ ]) to be optimal in the sensedefined above. The fit in Fig. 3 of the main text ischaracterized by N c = 5, N = 50 and ˜ χ = 0 . E = 0.0, W = 1.5, L y = 40 160, m m a x . r e l . e rr . max. rel. err. = 30% y _ _ R _ I _ _ R _ I _ _ R _ I _ _ R _ I _ _ _ R _ I _ _ _ R _ I _ _ _ R _ I _ _ _ R _ I _ _ _ R _ I _ _ _ R _ I _ _ _ R _ I _ _ _ R _ I _ _ _ R _ I _ _ _ R _ I _ _ _ R _ I _ _ _ R _ I FIG. S5. Stability plot showing the performance of differentfitting functions for the quasi-1D Lyapunov exponents ob-tained for E = 0 and W = 1 . m, L y ) datapoints. The labels on the x -axis encode the chosen fittingfunction, c.f. Eq. (S2). The green color indicates the bestfit for which the maximum relative error of all fit parametersis below 30%, see dahsed line in the second panel. In thethird panel, horizontal lines denote ν = 2 .
57 as obtained re-cently for the CC network model in Ref. [52] and our estimate ν E =0 = 2 . We finally remark that Fulga et al. [35] claim that theirscaling observable requires samples of large aspect ratio.However, we find that calculations with square samplesyield converging numerics and save considerable compu-tational resources. We discern no underlying argumentin favor of particular aspect ratios.1 E = 0.0, W = 2.0, L y = 40 160, m m a x . r e l . e rr . max. rel. err. = 30% y _ _ R _ I _ _ R _ I _ _ R _ I _ _ R _ I _ _ _ R _ I _ _ _ R _ I _ _ _ R _ I _ _ _ R _ I _ _ _ R _ I _ _ _ R _ I _ _ _ R _ I _ _ _ R _ I _ _ _ R _ I _ _ _ R _ I _ _ _ R _ I _ _ _ R _ I FIG. S6. Stability plot as in Fig. S5, but for E = 0, W = 2. Landauer conductance for square system at m = 0 , E = 0 Critical conductance distributions depend on the sam-ple shape and boundary conditions, and are typicallyvery broad [2, 4]. In addition, unlike the self-averagingLyapunov exponents, critical transport observables ex-hibit multifractality: their different moments are de-scribed by different scaling dimensions [54–56].In the main text, we have presented disorder averagedsquare system Landauer conductance g for m = 0 anda range of energies. Here, focusing on the case E = 0,we present the histogram of log g , see Fig. S10 (top) fora range of system sizes L = 28–280. It is qualitativelysimilar to the IQHT case [20, 57].In Fig. S10 (middle), we attempt a power law fit of the E = 0.3, W = 2.0, L y = 40 160, 0.0 m m a x . r e l . e rr . max. rel. err. = 30% y _ _ R _ I _ _ R _ I _ _ R _ I _ _ R _ I _ _ _ R _ I _ _ _ R _ I _ _ _ R _ I _ _ _ R _ I _ _ _ R _ I _ _ _ R _ I _ _ _ R _ I _ _ _ R _ I _ _ _ R _ I _ _ _ R _ I _ _ _ R _ I _ _ _ R _ I FIG. S7. Stability plot as in Fig. S5, but for E = 0 . W = 2. average conductance, using the ansatz g ( L ) = g − αL − ¯ y .It yields ¯ y = 0 . g = 0 . g . The bottom panelpresents a logarithmic fit g ( L ) = g ( ∞ ) − b log( λL ) as hasbeen suggested in literature as well [16].Notice that due to the above-mentioned multifractalityof conductances at criticality, we do not expect that theirrelevant exponent y describing the approach of g ( L ) toits limiting value is the same as the exponent y for theself-averaging Lyapunov exponent Γ, which is related tothe conductance of a strongly localized quasi-1D system.2 E = 0.7, W = 2.0, L y = 40 160, m m a x . r e l . e rr . max. rel. err. = 30% y _ _ R _ I _ _ R _ I _ _ R _ I _ _ R _ I _ _ _ R _ I _ _ _ R _ I _ _ _ R _ I _ _ _ R _ I _ _ _ R _ I _ _ _ R _ I _ _ _ R _ I _ _ _ R _ I _ _ _ R _ I _ _ _ R _ I _ _ _ R _ I _ _ _ R _ I FIG. S8. Stability plot as in Fig. S5, but for E = 0 . W = 2. Drude conductivity
In Ref. [39], Ostrovsky et al. derived Pruisken’s non-linear sigma model as an effective long-range theory forthe DDF. As emphasized in the main text, this derviationis only controlled if the longitudinal Drude conductivity σ D xx is much larger than unity. In Fig. S11, we plot thenumerically obtained Drude conductivity for the contin-uum DDF with smooth disorder at m = 0 for a rangeof energies E and disorder strengths W , see caption fordetails of the procedure. In the metallic limit of large E and small W , the data obeys σ Drude xx ∼ /W [39]. Forthe parameter combinations ( E, W ) used for the scalinganalysis above, the Drude conductivity in Fig. S11 is oforder unity such that the derivation of Ref. [39] is not P ( l o g | z | ) L = 60 P ( l o g | z | ) L = 110 z |0.00.2 P ( l o g | z | ) L = 200 FIG. S9. Histograms for the quantity log | z | underlying Fig.3 of the main text for m = 0 .
05. The solid lines are Gaussiansfor comparison. controlled. Note that this does not rule out the suggestedphase diagram in Fig. 1(b) of the main text.3 g P ( l o g g ) L = 28 L = 50 L = 140 L = 280 L g ( L ) g ( ) = 0.6 ± 0.04, y = 0.18 ± 0.06 L g ( L ) g ( ) = 0.73 ± 0.08, = 776 ± 2614 FIG. S10. Landauer conductance of square continuum Diracsystems, Eq. (3) of the main text, with m = 0, E = 0 andsmooth disorder, Eq. (4), of strength W = 2 .