Spectrum-wide quantum criticality at the surface of class AIII topological phases: An"energy stack"of integer quantum Hall plateau transitions

In the absence of spin-orbit coupling, the conventional dogma of Anderson localization asserts that all states localize in two dimensions, with a glaring exception: the quantum Hall plateau transition (QHPT). In that case, the localization length diverges and interference-induced quantum-critical spatial fluctuations appear at all length scales. Normally QHPT states occur only at isolated energies; accessing them therefore requires fine-tuning of the electron density or magnetic field. In this paper we show that QHPT states can be realized throughout an energy continuum, i.e. as an"energy stack"of critical states wherein each state in the stack exhibits QHPT phenomenology. The stacking occurs without fine-tuning at the surface of a class AIII topological phase, where it is protected by U(1) and (anomalous) chiral or time-reversal symmetries. Spectrum-wide criticality is diagnosed by comparing numerics to universal results for the longitudinal Landauer conductance and wave function multifractality at the QHPT. Results are obtained from an effective 2D surface field theory and from a bulk 3D lattice model. We demonstrate that the stacking of quantum-critical QHPT states is a robust phenomenon that occurs for AIII topological phases with both odd and even winding numbers. The latter conclusion may have important implications for the still poorly-understood logarithmic conformal field theory believed to describe the QHPT.


I. INTRODUCTION
Non-interacting topological quantum phases of matter feature robust gapless edge or surface states [1][2][3]. Remarkably, these states are protected from Anderson localization [4][5][6], defying the central dogma of conventional localization physics according to which all states localize in one dimension, as do most states in two dimensions [7]. An exception occurs for 2D systems with strong spin-orbit coupling and sufficiently weak disorder; in this case Anderson localization occurs near the band edges, but states near the band center can remain delocalized due to weak antilocalization [8]. By contrast, for the chiral [helical] edge states of the quantum Hall effect [2D topological insulators (TIs)], as well as the isolated 2D Dirac cone at the boundary of a 3D TI, all states spanning the bulk gap are prevented from localizing.
While the protection of chiral and helical 1D edge states is easily understood as the absence of elastic backscattering [1,2,9,10], the effects of disorder on 2D surface states of bulk topological phases is more subtle. Indeed, although pure backscattering is suppressed for the single Dirac fermion cone that forms at the boundary of the simplest 3D TI, elastic impurity scattering does occur at all other angles. The protection of the 2D surface states throughout the bulk energy gap is understood from the combination of weak antilocalization, as well as a Z 2 topological term in the effective field theory for the surface with quenched disorder (a nonlinear sigma model in the symplectic class) [7,[11][12][13][14][15].
In this paper, we consider the effects of disorder on Dirac node in the presence of quenched vector potential disorder A(r) of correlation length ξ and strength λ, modeling the surface state of a topological phase in symmetry class AIII. It is well known that the zero energy field theory is critical while finite energies E feature a power-law density of states (b) and an emergent length scale ζ ∼ E −1/z . (c) By studying the real-space characteristics of finite energy eigenstates (conductivity σ and the curvature of the multifractality spectrum θ), we identify ζ as a crossover scale beyond which criticality of the quantum Hall plateau transition (QHPT) type emerges. Our results imply that the finite-energy states form a "stack" of critical QHPT wave functions; this is in sharp contrast to the conventional expectation of Anderson localization [16]. We also show that this unusual energy-continuum of critical states is stable under the addition of a second Dirac node.
2D Dirac surface states of class AIII topological phases, at energies throughout the bulk gap. (See Table I for a review of the 10-fold classification). In the absence of strong interactions [17], class AIII is characterized by a Z-valued winding number ν in three dimensions, reflecting the number of 2D Dirac nodes on the surface. In this case, due to the absence of full spin-orbit coupling, protection from Anderson localization cannot be attributed to weak antilocalization. Class AIII could be realized either as a topological insulator on a bipartite lattice, with complex hopping that preserves sublattice symmetry but breaks time-reversal, or as a time-reversal invariant, spintriplet topological superconductor (TSC) that preserves a U(1) remnant of spin SU(2) symmetry [4,[18][19][20]. The effects of quenched disorder for class AIII surface states at energies close to zero are well known; here energy is measured relative to the surface Dirac point and the disorder comes in the form of a random vector potential. Topologically-protected continuum Dirac models were originally studied in the contexts of the quantum Hall effect [16] and d-wave superconductor quasiparticles [21][22][23][24][25][26]. The zero-energy states always escape Anderson localization [16,21]. Instead, these states are quantum critical, exhibiting strongly inhomogeneous spatial fluctuations on all length scales. These fluctuations can be characterized by the multifractal spectrum of wave function intensities [7], which is known exactly in these models [16,22,23,27]. Two other key attributes include a critical low-energy density of states ρ(E) that vanishes or diverges as E → 0, and a precisely quantized longitudinal conductivity at zero energy, independent of disorder [7,16,28,29] and (at zero temperature) interactions [30].
Despite this plethora of exact results at zero energy, very little was known about the character of finite-energy states or of the finite-temperature response. Ludwig et al. [16] argued that all finite-energy states of a 2D Dirac model with one node (corresponding to the bulk winding number ν = 1) should Anderson localize. Later, Ostrovsky et al. [14] instead conjectured that finite-energy states in this model could escape localization via a remarkable route: each of these states would sit exactly at the integer quantum Hall plateau transition (QHPT). The QHPT is a quantum phase transition that divides Anderson topological-insulating quantum Hall plateaux; it is known to also be quantum critical, with apparently universal multifractal spectra and average longitudinal conductivity [7,[31][32][33][34][35][36]. Ostrovsky et al. suggested however that this would only hold very close to the Dirac point [14]. Moreover, the derivation would seem to imply an "even-odd" effect, similar to the Haldane conjecture for half-integer versus integer Heisenberg spin chains. Specifically, finite-energy class AIII topological surface states with odd winding numbers are predicted to exhibit QHPT criticality, while those with even winding numbers would simply localize [5,14,15].
In this paper, we present strong numerical evidence that finite-energy surface states throughout the bulk gap of a class AIII topological phase exhibit QHPT criticality beyond the crossover length scale ζ ∼ E −1/z of the zero energy critical theory, see Fig. 1 for a graphical summary. Moreover, we find no evidence of an even-odd effect. Calculations are performed for winding numbers ν = {1, 2} using two different model types: (1) effective low-energy 2D surface field theories and (2) surface states of 3D topological lattice models in the slab geometry. For the continuum theories, we compute the energy-resolved Landauer conductivity, Landauer conductance distribution, multifractal spectrum, and density of states. We observe the crossover of the conductivity between the exact zeroenergy class AIII result and the known average conductivity of the QHPT [32][33][34]36]; concomitantly we observe QHPT multifractality. We demonstrate the crossover between two copies of the ν = 1 model to ν = 2 as a function of internode scattering. For square systems, we obtain a Landauer conductance distribution that scales towards the known QHPT result [35] with increasing system size for both ν = 1, 2. Our results confirm and significantly extend a previous glimpse of QHPT multifractality at finite energy reported in Ref. [37]. For the lattice models, we provide evidence that states throughout the gap avoid localization, and that finite-energy states exhibit multifractal spectra consistent with the QHPT. We summarize our main findings in Fig. 1.
We conclude that class AIII appears to solve the problem of topologically protecting its surface state spectrum throughout the bulk gap by forming an energy continuum of QHPT states. This is surprising, given that in the 2D quantum Hall effect, critical states appear only at isolated energies; accessing them in experiment requires fine-tuning of the electron density or magnetic field to a quantum critical point. The phenomenon of spectrumwide criticality observed here is very unusual, but has been previously reported in two other cases: models for surface states of class CI and DIII topological superconductors (TSCs) [38,39]. The finite-energy states in class CI [38] were found to exhibit critical statistics consistent with the class C spin quantum Hall plateau transition [7,[40][41][42][43][44][45][46]. A single surface Majorana fermion cone is predicted to occur at the boundary of a class DIII TSC, such as the candidate material Cu x Bi 2 Se 3 [2,47]. In the presence of disorder, finite-energy class DIII Majorana surface states also appear to exhibit robust, universal criticality [39], conjectured to be related to the thermal quantum Hall plateau transition in class D [48][49][50][51][52][53][54]. The 10-fold way symmetry classification is reviewed in Table I, with an emphasis on the implied connections between surface states of the 3D topological classes CI, AIII, DIII, and 2D topological phase transitions in classes C, A, D.
Together, the results obtained in this paper and in Refs. [38,39] imply that the ordinary 3D class AII, Z 2 topological insulator is actually the exceptional case. The metallic surface states of a class AII TI do not exhibit universal criticality; instead, the classical surface conductivity at each energy in the gap is determined by the impurity density and the density of states, and this  [4]. The 10 classes are defined by different combinations of the three effective discrete symmetries T (time-reversal), P (particle-hole), and S (chiral or sublattice). For a d-dimensional bulk, any deformation of the clean band structure that preserves T , P , and S and does not close a gap preserves the topological winding number. For a (d − 1)-dimensional edge or surface theory, the equivalent statement is that any static deformation of the surface (quenched disorder) that preserves T , P , and S also preserves the "topological protection" against Anderson localization. Of particular interest here are classes C, A, D on one hand, and classes CI, AIII, DIII on the other. Classes C, A, and D are topological in d = 2, and describe the spin (SQHE), integer (IQHE), and thermal (TQHE) quantum Hall effects; all three can be realized as TSCs with broken T . Classes CI, AIII, and DIII are topological in d = 3, and can describe 3D time-reversal-invariant TSCs. (In this case, the physical time-reversal symmetry appears as the effective chiral symmetry S in the table [4,18].) The column "spin sym." denotes the amount of spin SU(2) symmetry preserved for TSC realizations of these 6 classes. In this work, we show that the finite-energy surface states of the 3D class AIII topological phase appear to form a "stack" of quantum-critical, Anderson delocalized wave functions. We provide evidence that each state in the stack is statistically identical on large length scales, and corresponds precisely to the topological quantum phase transition in class A (the integer quantum Hall plateau transition). Previous numerical results also uncovered a gap-spanning "stack" of critical states at the surface of the class CI TSC, which match the statistics of the class C spin quantum Hall plateau transition [38]. A third study [39] revealed spectrum-wide criticality at the surface of class DIII, conjectured to be related to the thermal quantum Hall plateau transition in class D. Thus the numerical results presented in this paper and in Refs. [38,39] appear to connect the quantum Hall plateau transitions in classes C, A, and D to gap-spanning stacks of critical states at the surfaces of class CI, AIII, and DIII 3D topological phases. The conclusion of gap-spanning surface criticality, locked to plateau transitions in classes C, A, D, is in sharp contrast to the conventional expectation for finite-energy 2D surface states in classes CI, AIII, DIII. The conventional expectation is that finite energy always breaks the defining chiral S symmetry, producing a standard Wigner-Dyson class, so that [7] CI → AI (Anderson localized), AIII → A (Anderson localized), and DIII → AII (Anderson localized or weak antilocalization) (see however [55]).
nonuniversal, energy-dependent value is continuously enhanced at larger distances due to weak antilocalization [7,[11][12][13][14]. On the other hand, three of the five [4,20] topological classes in 3D exhibit spectrum-wide quantum criticality. There appears to be an unexpected, integral connection between Z-graded topological phases in two and three spatial dimensions: the quantum phase transitions in the former (classes C, A, D) are bundled into gap-spanning stacks of surface eigenstates in the latter (classes CI, AIII, DIII). Our results should have a deeper topological underpinning, which however remains yet to be uncovered.
Our results may also have implications for the theory of the class A integer QHPT itself. The transition has been extensively studied numerically [7,31], with universal signatures consistent with a conformally invariant critical point. Despite 30 years of effort, very little is known about this critical point analytically, which is believed to be governed by a logarithmic conformal field theory (LCFT) [56]. A recent proposal by Zirnbauer [57] in fact argues that the QHPT is itself governed by a zeroenergy, class AIII theory with winding number ν = 4. An additional special feature is that the purely marginal abelian disorder strength (a generic feature for class AIII theories, reviewed below) is fine-tuned to a special value so as to render the density of states non-critical. A very surprising consequence of this proposal is that the correlation length exponent ν QHPT is actually predicted to be infinite; the claim being that existing numerical studies showing ν QHPT in the range between 2.3-2.6 [7,31,[58][59][60][61][62][63][64][65][66] are beset by finite-size effects. We note that for the finite-energy class AIII surface states that exhibit QHPT critical statistics, the exponent ν QHPT does not play a role because all states are delocalized.
If Zirnbauer's proposal is correct, then our numerical results suggest that both the zero-and finite-energy states of the class AIII bulk topological phase with the special winding number ν = 4 are governed by the same conformal field theory. In this scenario, the energy perturbation serves only to fine-tune the abelian disorder strength to the correct value so as to achieve QHPT criticality. At least this proposal should be falsifiable using the class AIII surface state theory, which is a relatively well-understood LCFT (a Wess-Zumino-Novikov-Witten model [21][22][23][24]67]).

A. Outline
This paper is organized as follows. In Sec. II, we transcribe a disordered Dirac Hamiltonian for two nodes. We enumerate topological (anomalous) and regular symmetry operations and identify ν = 1, 2 class AIII surface state models. We also define two non-topological models that will be studied for comparison: a Dirac version of the Gade "sublattice hopping" class AIII model [68][69][70], and a topologically trivial version of spinless graphene in the unitary class A. We summarize key known results for these models, as well as for the quantum Hall plateau transition.
In Sec. III, we present numerical results for the ν = 1 topological AIII model at zero and finite energies. These include the density of states, the multifractal spectrum for wave function and local density of states fluctuations, the Landauer conductivity, and the distribution function for the Landauer conductance. In Sec. IV, we repeat this analysis for the ν = 2 topological surface state model. In Sec. V, we present analogous results for the non-topological Dirac models defined in Sec. II.
Results for multifractal spectra of surface states obtained from a 3D topological lattice model in the slab geometry, with bulk winding numbers ν = {1, 2}, are presented in Sec. VI. We discuss our results in the context of similar findings in classes CI and DIII [38,39] in Sec. VII, and in particular highlight possible implications for the analytical understanding of the logarithmic conformal field theory [56] governing the QHPT. The latter is explored in light of a very recent proposal by Zirnbauer [57] linking the QHPT to an effective AIII surface state model with winding number ν = 4.

II. 2D DIRAC MODELS, SYMMETRY CLASSES, AND TOPOLOGICAL SURFACE THEORIES
The low-energy physics of the class AIII topological surface states that we study in this paper in Secs. III and IV, along with two non-topological models employed for comparison (Sec. V), can be encoded in a 2D, two-valley Dirac-fermion Hamiltonian. The most generic model can be considered the low-energy effective field theory for spinless graphene, with the two cones corresponding to the K and K valleys. The Hamiltonian is [7] For spinless fermions hopping on the honeycomb lattice, three key symmetries are T 2 = +1 time-reversal, P 2 = +1 particle-hole, and sublattice S symmetry. The latter is the product of T and P . With respect to the single-particle Hamiltonianĥ, these symmetries are encoded in the conditions [72] T :σ 2τ2ĥ * σ 2τ2 =ĥ, (3a) In the graphene context, {m 1,2 }, m 3 , and m 0 respectively denote Kekulé, CDW, and Haldane masses [3,73]. The particle-hole symmetric vector potentials {Aā ,3 } can encode elastic strain (via pseudomagnetic fields) [74]. The scalar potential is V 0 , while real magnetic fields can be encoded in {Aā ,0 }.

A. Non-topological class A and AIII models
In the absence of T , P , and S, the Dirac Hamiltonian in Eq. (2) resides in the unitary class A. If we take all 16 bilinear perturbations to be Gaussian random variables with zero average (including a vanishing average Haldane/Chern insulator mass m 0 = 0), then this theory realizes a topologically trivial Anderson insulator, with localized states at all energies [7,29].
We can get a purely two-dimensional, nontopologicalsurface-state version of class AIII by imposing only sublattice symmetry S. In the honeycomb lattice model, sublattice symmetry is implemented via This is a symmetry for real or complex nearest-neighbor hopping on any bipartite lattice at half-filling [68,69]. The allowed bilinear perturbations implied by Eq. (3c) . The vector potentials act independently on the two nodes, while the scalar and mass potentials scatter between them. This class AIII Hamiltonian can alternatively describe Bogoliubov-de Gennes quasiparticles in a 2D p x -wave, spin-triplet ("polar" [75,76]) superconductor [18]. At the band center (≡ zero energy), sublattice symmetry introduces additional Goldstone modes in the nonlinear sigma model target manifold relative to ordinary metals. As a result, states near zero energy evade Anderson localization, as originally shown by Gade and Wegner [68,69]. Localization is still possible for strong random hopping [77]. On the other hand, at energies away from zero we expect eigenstates to reside in an ordinary Wigner-Dyson class. For a system possessing only sublattice symmetry, the finite-energy states are expected to reside in the unitary class A, which localizes without fine-tuning to the quantum Hall plateau transition. (See Sec. VII A for a precise nonlinear sigma model formulation in the topological case.) Thus all finite-energy states of a non-topological class AIII model in two dimensions are expected to localize.
There are in fact two unitarily inequivalent types of chiral symmetry available for 2D Dirac fermions [78]. The sublattice symmetry in Eq. (3c) defines a 2D class AIII Dirac theory with mass, scalar potential, and vector potential perturbations. Since one can gap out the system with a nonzero average mass without breaking the defining chiral symmetry, this cannot be a topological surface state theory. On the other hand, for the same 2-node Dirac Hamiltonian in Eq. (2), we can introduce an "anomalous" version of the chiral symmetry, This differs from Eq. (3c) due to the absence of a grading in valley space. In the context of a class AIII topological superconductor (TSC), Eq. (5) represents physical timereversal symmetry realized at the boundary of the sample [4,20,67]. In a chiral bulk topological insulator, it is the surface projection of the bulk sublattice symmetry [19]. Eq. (5) is anomalous because it cannot be realized via a local symmetry transformation in a strictly 2D lattice model [unlike Eq. (3c), which is the continuum version of the lattice operation in Eq. (4)]. Imposing Eq. (5) on the Dirac model, only vector potential perturbations are allowed. I.e., the Hamiltonian in Eq. (2) reduces tô where the repeated indicesā ∈ {1, 2} and i ∈ {1, 2, 3} are summed. Eq. (6) describes 2D Dirac fermions perturbed by quenched random U(1) (Aā ,0 ) and SU(2) (Aā ,i ) vector potentials. If we further suppress scattering between the two nodes, we get two copies of the simpler single-node theoryĥ In this case only abelian vector potential disorder appears. Eqs. (6) and (7) can be realized as surface state theories for class AIII TSCs with bulk winding numbers ν = 2 and ν = 1, respectively. Alternatively, two independent copies of Eq. (7) obtain from the 2D class AIII Gade-Dirac model defined by Eq. (3c), when internode scattering (mediated by {m 1 , m 2 , V 1 , V 2 }) is suppressed by hand. Similar fine-tuning of impurity amplitudes in a model for 2D d-wave superconductor quasiparticle scattering can realize 2 or 4 copies ofĥ (2) AIII orĥ (1) AIII , respectively [25,39].
The effective field theory describing the wave function and transport statistics near zero energy for the surface state models in Eqs. (6) and (7) is a class AIII Wess-Zumino-Novikov-Witten (WZNW) nonlinear sigma model [21][22][23][24]67] (see Sec. VII A). It is identical in form to the principal chiral sigma model relevant to Gade-Wegner physics in class AIII [68][69][70], except that the latter lacks the WZNW term. The presence or absence of this term completely changes the low-energy physics [7,70,78]. The key low-energy properties of Eqs. (6) and (7) are as follows.
The non-topological 2D Gade-Wegner models also feature an abelian parameter λ A ; this is related to the variance of the vector potentials {Aā ,0 , Aā ,3 } in the model defined by Eq. (3c). Unlike the topological surfacestate theory, however, the λ A parameter always flows to infinity under the RG [68][69][70]. For this reason, the low-energy physics of 2D Gade-Wegner models is always frozen [81,82].
2. Finite energy: correlation length and density of states-A finite energy E = 0 introduces a scale to the conformally invariant zero-energy theory. Formally, energy is a relevant coupling strength that flows to large values under the renormalization group [16,38,39]. We can associate to this relevant coupling a correlation length ζ(E) ∼ |E| −1/z , with the dynamical critical exponent [21,22,67] We will see that ζ(E) plays an unconventional role in the characterization of finite-energy eigenstates for topological class AIII Hamiltonians [Eqs. (6) and (7)], in that it divides two different quantum critical scaling regimes, see Fig. 1 and Secs. III, IV and VI. The length ζ(E) thus plays the role of a crossover scale between different critical fixed points, instead of dividing a critical point (the delocalized zero-energy theory) from a massive fixed point (which would be an Anderson insulator at finite energy). The dynamical critical exponent z also determines the scaling form of the low-energy density of states (DoS) ρ(E) for the Hamiltonians in Eqs. (6) and (7), which vanishes or diverges with a power-law, Eq. (12) assumes weak λ A , such that the system is not "frozen" at zero energy (in the sense discussed above).
3. Zero-temperature, zero-energy conductivity-At zero temperature and zero doping relative to the Dirac point, the Landauer conductivity ofĥ (ν) AIII is independent of the disorder [16,28,29], and equivalent to that of the clean limit [7,85] in units such that the conductance quantum e 2 /h = 1. This result also describes charge conduction through a wide, perfectly clean nanoscopic flake of graphene, precisely doped to the Dirac point [85]. For the topological WZNW class AIII theory, Eq. (14) appears to hold even in the presence of both disorder and interactions [30]. By contrast, in the Gade-Wegner class AIII Dirac model, the conductivity at zero energy is non-universal, depending on the microscopic strengths of the mass and potential disorders [7,70]. Interactions further modify the conductivity via Altshuler-Aronov corrections [18,86]; the latter are excised by the WZNW term in the topological case [30].

C. Integer quantum Hall plateau transition
In Secs. III and IV, we will look for signatures of the class A quantum Hall plateau transition (QHPT) in the finite-energy eigenstates of Eqs. (6) and (7). The key attributes we detect are the multifractal spectrum and Landauer transport properties. From numerical studies of the QHPT, the τ (q) spectrum is known to be approximately parabolic, given by Eq. (8) with θ QHPT 1/4 [7,31]. At the plateau transition, the disorder-averaged conductivity has been obtained from a Kubo-formula calculation of a disordered tight-binding model in a magnetic field [36]. Extrapolating results on very long samples with finite width to infinite width, the result is again in units such that e 2 /h = 1. The distribution of the critical Landauer conductance g has been computed for square samples of the Chalker-Coddington network model [35] and the tight-binding model [36]. In the limit of large system sizes, g (QHPT) = 0.60 ± 0.02 was obtained [36]. For earlier studies see Refs. [31][32][33][34]. In this section, we present numerical evidence for a continuous band of critical QHPT states in a single Dirac node with a random U(1) vector potential. This corresponds to the ν = 1 class AIII surface state Hamiltonian in Eq. (7), which is a 2 × 2 matrix differential operator, protected by the "topological" version of chiral symmetry in Eq. (5). For analytical approaches, the white noise disorder correlator Eq. (10) along with an UV cutoff on the Dirac dispersion is a convenient choice. In contrast, for our numerics, we find it useful to work with finite-range disorder correlations and take the UV-cutoff for the dispersion to infinity. We assume the disorder fields Aā ,0 (r) (ā ∈ {1, 2}) to vary smoothly on a scale ξ, and these fields are taken to have zero average value over the sample area. The disorder statistics are taken to be Gaussian with the correlator with disorder correlation length ξ; in this and the following two sections we restore . The dimensionless dis- In the limit of ξ being the smallest scale, the exact analytical results [16] for the white noise disorder case should hold at low energies. The dynamic critical exponent is predicted to be given by [Eq. (11) with ν = 1] This corresponds to the purely marginal parameter λ A = W 2 in Eq. (10), which parameterizes a line of fixed points [16]. These results are valid in the low-energy |E| → 0 limit, below the freezing transition that occurs at W c = Landauer conductance g (in units of e 2 /h) of square samples of length L for the topological class AIII surface model with a single Dirac node (ν = 1), as in Fig. 2. The probability distribution p(ln g) of the logarithm of the conductance for L/ξ = 60, 100, 180, 300 is depicted in panels (a)-(d) for energies Eξ/ v = 0.01, 0.1, 0.8, 1.0, respectively. The data is based on between 1000 and 10000 disorder realizations, depending on system size. The grey filled area depicts the result obtained for a square-shaped Chalker-Coddington network model of size 128 reproduced from Ref. [35]. Panel (e) shows the scaling of the mean square conductance g(L) for the above lengths and energies and the asymptotic value for the QHPT obtained from Ref. [36]. For the horizontal axis, the exponent on L has been chosen arbitrarily. √ 2π 2.5 [27,37,79,80]. In the following, we will work with W = 2.3 in order to keep disorder-induced length scales short while staying below W c . We then expect z = 2.68.
We start with an exact diagonalization (ED) study of the Hamiltonian in Eq. (7), regularized on a lattice in momentum space, assuming periodic boundary conditions in real space. Compared with tight-binding models in real space, the momentum space approach avoids fermion doubling and band bending effects on the one hand, but results in a dense Hamiltonian matrix which limits the available system sizes on the other hand. We choose a linear dimension of L = 60ξ and a momentum cutoff ξ|kā| ≤ R (ā ∈ {1, 2}) with R = 5 (corresponding to a matrix Hamiltonian of size 18050), which gives sufficient real space resolution to resolve the smooth variation of the disorder field. We checked that the results are converged with respect to R and the number of disorder realizations (typically a few hundred).
In Fig. 2(a) we plot the ED-DoS ρ versus energy E > 0. The scaling form ρ(E) ∼ E 2/z−1 implies a divergence at zero energy when z > 2 as in our case. In our finite size system, the divergence is replaced by a peak. At larger energies above ∼ 2 v/ξ, the DoS is unaffected by disorder scattering and asymptotically approaches the clean value ρ 0 (E) = E 2π v (dashed line). In Fig. 2(b) we quantify the DoS power law by plotting N (E) = E 0 d ρ( ) (solid) and confirm that it goes like N (E) ∼ E 2/z for small E (dashed line). This agreement with non-trival analytical predictions validates our numerical model and demonstrates precise control over the disorder strength.
We now turn to the results of a transport calculation for which we assume x to be the transport direction. The time independent Schrödinger equationĥ (1) AIII ψ = E ψ can be written as whereÛ ≡ A 1,0σ1 +A 2,0σ2 and k F ≡ E/ v. This can be solved in terms of the transfer matrix, using the method of Ref. [11]. We assume periodic boundary conditions in the transverse direction with |k y | = 2π|n|/L y ≤ R/ξ (n ∈ Z) with L y = 400ξ, and again use the mode cutoff R large enough so that the results below do not depend on it. Assuming clean, highly doped leads attached at x = 0 and x = L x , the propagating lead modes areσ 1 eigenstates [11]. This allows the definition of a scattering matrix S. From the transmission block t, the longitudinal conductance can be found as G = e 2 h tr t † t . We take the disorder average G and plot the resistance, normalized to the sample width L y /G in Fig. 2(c) as a function of L x . For a diffusive sample, we expect L y /G = L y R 0 + 1 σ L x where R 0 is some contact resistance and σ the bulk conductivity. The data for E = 0 agrees with σ = e 2 /hπ [Eq. (14)] and R 0 = 0 for each disorder realization, as required by gauge invariance and chiral symmetry [29,87].
For any finite energy, there is a crossover of the resistance trace to diffusive behavior with a larger conductivity. We define the crossover scale ζ as the length L x where L y /G deviates by 5% from the E = 0 result. In Fig. 2(d) we show that the crossover scale for E v/ξ indeed follows the scaling ζ ∼ E −1/z expected for the correlation length.
In Fig. 2(e), we plot the conductivities as fitted from the resistance data well above the crossover scale for L x > ζ + 10ξ. As the resistance trace is not perfectly linear in L x , there is a few percent ambiguity in the definition of σ that we further address in the next paragraph. For 0 < E v/ξ, the conductivities show a plateau at σ 0.55 e 2 h in fair agreement with the value σ (QHPT) x,x = 0.58 ± 0.02 e 2 h obtained by Schweitzer and Markoš [36] via the Kubo formula for a lattice model at the QHPT transition. At larger energies E, the conductivity at the accessible length scales increases with energy. This is as expected for the crossover to the semiclassical Drude conductivity, which goes as σ ∼ (e 2 /h)(1/W 2 ) [72].
To get further insight into the behavior of quantum transport as a function of system size, we next study the probability distribution of the Landauer conductance g(L) for square samples with increasing length L, keeping the disorder strength the same as above. The results for energies Eξ/ v = 0.01, 0.1, 0.8, 1.0 are shown in Fig. 3, panels (a)-(d). In all cases, the distributions p(ln g) agree reasonably well with the established distribution obtained from a square shaped Chalker-Coddington network model of size 128, reproduced from Ref. [35] (grey shaded area). For the energies Eξ/ v = 0.01, 0.1, an evolution of the distribution towards the Chalker-Coddington result can be observed with increased system size L, whereas for Eξ/ v = 0.8, 1.0, the distribution does not change visibly on the available length scales. The evolution of the average square conductance g(L) with L is depicted in Fig. 3(e). For energies Eξ/ v = 0.01, 0.1, we find an increase with system size, whereas there is a weak decrease for Eξ/ v = 0.8, 1.0 (not discernible in the conductance distributions). The data for all energies is consistent with a limiting value of g (QHPT) = 0.60 ± 0.02, as obtained in Ref. [36] by extrapolating results for a disordered lattice model to large system sizes (grey line). For the horizontal axis in panel (e), the exponent −0.5 on L has been chosen arbitrarily. Our range of length values is not sufficient to determine if there is indeed a power law approach of g(L) to its limiting value with a characteristic irrelevant exponent, as established for the fine-tuned QHPT in Ref. [36]. Our data suggest that the relevant length scale for the approach to the QHPT fixed point strongly increases with energy beyond E v/ξ. Interestingly, this parallels the behavior of the Boltzmann transport mean free path for our correlated disorder model where the scattering range in momentum space is restricted to 1/ξ.
Further evidence for the presence of a stack of QHPT critical states at finite energies comes from the anomalous part of the multifractal spectrum ∆(q) = τ (q) − 2(q − 1), as shown in Fig. 2(f). Here, τ (q) is defined via the scaling of the disorder averaged generalized inverse participation ratios [7], where, for a single disorder realization, and the outer sum is over square regions B with linear size b covering the real space lattice. To find P q , we use the ED eigenstate ψ closest to energy E. It turns out that for small energies E 0.1 v/ξ, the power law in Eq. (19) is governed by two different multifractal spectrums τ (q), with the crossover at b ζ as found in Fig. 2(d). To capture only the long distance physics beyond the correlation length, we restrict to b > ζ, which is only possible for E ≥ 0.01 v/ξ due to the overall system size restriction. For a wide range of energies E v/ξ, we find ∆(q) in good agreement with the established form for QHPT states [7,31] with θ QHPT 0.25. The class AIII multifractality from Eq. (9) expected at low energies and for box sizes below the length scale ζ is not clearly visible in our numerics. This is due to insufficient system size for the high degree of rarification associated to our large disorder strength (chosen just below the freezing transition [27,79,80]). For a similar study at weaker disorder strength confirming Eq. (9) at low energies, see Ref. [37].
As above for ν = 1, we start with an ED study of Hamiltonianĥ (2) AIII and use L = 60ξ and R = 5 as a UVcutoff parameter (corresponding to a matrix Hamiltonian of size 36100). Fig. 4(a)  expected peak at energy E = 0 and a crossover to the clean DoS of two Dirac cones (dashed line) at large energies. In Fig. 4(b) we confirm the expected power law of the DoS ρ(E) ∼ E 2/z−1 by plotting N (E) ∼ E 2/z . We now turn to quantum transport results, obtained as in the ν = 1 case. We use samples of width L y = 400ξ and increase the length L x up to 350ξ. Fig. 4(c) shows the resistance normalized to sample width as a function of length L x for a range of energies between E = 0 and E = 1.2 v/ξ. We define a crossover length ζ from the resistance data as above for ν = 1 and confirm the expected scaling ζ ∼ E −1/z in Fig. 4(d) for energies E 0.3 v/ξ. Beyond this crossover length, the slope of the low energy resistance curves in Fig. 4(c) increases and then roughly saturates for L x 200ξ. The conductivity fitted from this large L x data is plotted versus energy in Fig. 4(e); we find values in close vicinity to σ (QHPT) x,x [Eq. (15)] for small E and conductivities in-creasing with energy for E 0.3 v/ξ. To further investigate the scaling of transport properties with system size, we study the probability distribution of the Landauer conductance g(L) of square samples of length L, keeping the disorder strength the same as above. The results for energies Eξ/ v = 0.01, 0.1, 0.8 are shown in Fig. 5, panels (a)-(c), and the evolution of the average conductance g(L) with L is depicted in Fig. 5(d). The situation is similar to the ν = 1 case. While small energies Eξ/ v = 0.01, 0.1 result in distributions and average square conductances visibly approaching the QHPT case, at higher energy Eξ/ v = 0.8 the relevant length scale for this process exceeds our available system size and only a slight trend of the average conductance is visible.
Finally, Fig. 4(f) shows the anomalous part of the multifractal spectrum, as calculated from the ED eigenstates in the vicinity of the energies used for the transport study. In analogy with the conductivities, we find good agree- (d)  Fig. 4. The probability distribution p(ln g) of the logarithm of the conductance for L/ξ = 60, 100, 180, 300 is depicted in panels (a)-(c) for energies Eξ/ v = 0.01, 0.1, 0.8, respectively. The data is based on between 3000 and 600 disorder realizations, depending on system size. The grey filled area depicts the result obtained for a square shaped Chalker-Coddington network model of size 128, reproduced from Ref. [35]. Panel (d) shows the scaling of the mean square conductance g(L) for the above lengths and energies and the asymptotic value for the QHPT obtained from Ref. [36]. For the horizontal axis, the exponent on L has been chosen arbitrarily. x,x in the absence of internode scattering, to a plateau with value 1 × σ (QHPT) x,x in its presence. The data is obtained for systems of width Ly = 400ξ and is averaged over 600 disorder realizations.
fixed energy E = 0.03 v/ξ and present transport data in Fig. 6. For W N = 0, the two nodes are decoupled and we find a conductivity close to 2 × σ (QHPT) x,x in agreement with the single node case presented in Sec. III. For W N = 1, 2, 3 the nodes are coupled and the conductivity is close to the value 1 × σ (QHPT) x,x .

V. NON-TOPOLOGICAL TWO-DIRAC NODE MODELS: GADE-WEGNER AIII, LOCALIZED GRAPHENE
We now contrast our findings for the topologicalsurface-state class AIII models presented in Secs. III and IV with two non-topological Dirac models, both of which  [7,29]. All disorder potentials are taken to carry the same strength W = 1.5. In Fig. 7, we show numerical results for this model in the top row, panels (a,b,c). The numerical results are obtained analogously to those presented in the previous sections. The DoS in panel (a) is finite and featureless at and around zero energy. In panel (b), we show the IPR P q=2 , where we use a box size corresponding to the disorder correlation length b ξ. While the IPR is of the same order as for class AIII-WZWN (crosses), the localizing behavior is clearly observed in the resistance data in panel (c) probing larger length scales. For small energies E ≤ 0.3 v/ξ the resistance is growing with increasing slope to values much above the clean zero energy resistance (dashed line).
We also consider the non-topological Gade-Wegner class AIII Dirac model, defined by imposing the sublattice symmetry in Eq. (3c). Just as in the topological case studied in Sec. IV, there are eight allowed disorder perturbations. The crucial difference from the topological model in Eq. (6) is that Eq. (3c) permits mass and scalar potentials, in addition to vector potential disorder. The allowed perturbations from Eq. (2) are 70]. Numerical results are shown in Fig. 7(d,e,f). The DoS in panel (d) shows the characteristic Gade-type singularity at E = 0 [Eq. (13)], and the IPR in panel (e) is increased as compared to the topological version from Fig. 4. Concurrently, the tendency to localization in the resistance plot in (f) is clearly observable. However, note that the conductivity at zero energy in the Gade class AIII Dirac model is finite and slightly lower than the clean conductivity σ = 2e 2 /(hπ) (dashed line), in agreement with theoretical predictions [29,70].

VI. DIRTY SURFACE STATES OF A CLASS AIII 3D TOPOLOGICAL LATTICE MODEL
We complement the results from the former sections dealing with pure surface theories by exactdiagonalization studies of a slab topological superconductor system. We investigate the multifractal spectra of surface states in the superconducting gap. For this purpose, we take the model from Ref. [30]. Although not microscopically realistic, this model encodes a class AIII topological superconductor on the diamond lattice, as explicated in Fig. 8 and defined below. By tuning of parameters one can realize winding numbers of ν = −2, −1, 0, 1, 2. See Fig. 9 for a sketch of the topological phase diagram. This means we can study both even and odd winding numbers by tuning parameters in this model.
The lattice Hamiltonian can be decomposed into 3 parts, H = H 1 +H 2 +H 3 . We transcribe the model in the absence of disorder as follows (see also Fig. 8). The first kinetic term consists of nearest-neighbor hopping from one sublattice (A) to the other (B): (21) where C iν (R) annihilates an electron at site R on sublattice i ∈ {A, B}, with spin polarization ν ∈ {↑, ↓}. Here the displacement vectors {δ} connect nearest-neighbor A and B sites, and "H.c." is the Hermitian conjugate. A second kinetic term involves a staggered chemical potential and next-nearest-neighbor hopping in the xzand yz-planes, where the vectors {δ} connect next-nearest-neighbor (same-sublattice) sites on the diamond lattice, and where The Hamiltonian in Eqs. (21) and (22) [7,31] within 4% for 75% of the q ∈ [0, qc] below freezing qc = 2.83 [38,39]. The total DoS is shown in blue, the critical fraction is colored red. (c) Fits for extracting the multifractal exponents are shown for an exemplary state at E = 0.47. Commensurate box sizes are shown as red dots, incommensurate ones (padded with zeros) are blue. When fitting only commensurate boxings are taken into account, and these satisfy the expected linear relation very well even for large q. One can observe that log P b q depends linearly on log b for a single state.
In the last equation, the Pauli matrices {μ 1,2,3 } act on spin-1/2 space. The key takeaway is that we observe no even-odd effect: both ν = 1, 2 exhibit critical surface states (with approximate QHPT multifractality) spanning the bulk energy gap. An even-odd effect would mean critical (Anderson localized) finite-energy states for odd (even) winding numbers. Such an effect is predicted from a naive deformation of the WZNW term in the theory governing zero-energy critical states, see Sec. VII A. The ED studies on the lattice model agree with S-matrix and momentum space (Secs. III, IV) studies in that both ν = 1, 2 match the QHPT critical behavior at finite energies.
We implement the Hamiltonian for a slab geometry. In the x-direction we terminate the diamond lattice in two surfaces, with W pairs of (L × L) A, B layers in between. We choose periodic boundary conditions in y-and z-directions to obtain a quasi-infinite slab of thickness W . Disorder is implemented by choosing each coupling λ ∈ {t , µ s , t, ∆} [Eqs. (21)- (26)] from a box distribution λ ∈ [λ−δ xλ ,λ+δ xλ ], whereλ is the average value of the coupling strength. In order to avoid localization of the bulk states, we choose δ x that quickly decays away from the two surfaces, so that the disorder effectively resides only near sample boundaries. Suppressing bulk disorder further prevents shifting of topological phase boundaries as function of the disorder strength.
To obtain all the in gap states of the sparse Hamiltonian, we employ the FEAST algorithm [88]. In numbers, this enables us to compute O(10 3 ) eigenstates in the gap for matrix sizes up to 0.5×10 6 to accuracy of 10 −10 times the matrix norm. As one can see by comparing to the data in Refs. 59, 63, and 89 on QH critical wave functions on the network model, we can access the IPRs with reasonable precision using the system sizes under consideration. To investigate parabolicity and critical exponents, i.e. accessing more than a decade of systems sizes with averaging, our computational power is insufficient. We therefore leave precision finite size scaling studies of the topological surface wave functions to future work.
The lattice wave function carries crystal coordinates {x, y, z}, sublattice index i ∈ {A, B}, and spin index ν ∈ {↑, ↓}. The surface-resolved probability density of an eigenstate wave function ψ x,y,z;i,ν is defined via where we trace over the slab depth x, as well as sublattice and spin spaces. The surface states decay quickly (exponentially) into the bulk, and our results are not dependent on how many slices orthogonal to the x-direction are taken for states deeply in the gap. We define the box probability µ n and the inverse participation ratio P b q as usual, We subdivide the (L × L) surface into square boxes {A b n } of linear size b. For the linear regression to extract the multifractal exponents τ q , only commensurate box-sizes are taken into account: For 0 < q < 1, the restriction to commensurate box-sizes is crucial. For a large system with linear surface size L = 72 and even winding number ν = 2, we show the density of states and eigenstate multifractality in Fig. 10. Energy E is measured in units of the nearest-neighbor hopping strength t . For the analysis presented in Fig. 10, we only extract the in-gap states that comprise less than 1% of the entire energy spectrum. Throughout the gap, a power-law integrated density of states N (E) ∝ E α , with α close to unity is obeyed [ Fig. 10(a)]. The dimensionless disorder strength δ implemented on all surface couplings of the microscopic lattice model is related to the abelian (λ A ) and nonabelian disorder strengths of the surface theory [Eqs. (6) and (10)] in a non-trivial way. For this reason, a comparison of bare and effective disorder via Eq. (12) as in the preceding sections is therefore not revealing in this case.
In Fig. 10(b), a histogram exhibiting the density of critical versus the density of all states is shown. A state is critical, once its anomalous multifractal exponents ∆ q match the class A QHPT prediction ∆ (QHPT) q (0.25)q(1 − q) [7,31] within 4% for 75% of the q ∈ [0, q c ] below freezing q c = 2.83 [38,39]. Almost a third of the in-gap states satisfies this. We verify that the inverse participation ratios log P b q depend linear on log b and show this for an exemplary state in Fig. 10(c). Further in Fig. 10(d), we check that the wave functions in the gap consist primarily of surface states by analyzing their behavior integrated over slices in the open direction x. Here, the total surface probabilities are defined via Finally, panels (e) and (f) in Fig. 10 exhibit anomalous multifractal spectra for states spanning the bulk gap. Although states near zero energy exhibit stronger multifractality, most of the in-gap states show weaker fractality, consistent with the QHPT prediction.
In Fig. 11, we show disorder-averaged multifractal spectra for smaller (L = 48) systems. We compare the two winding numbers ν = 1 and ν = 2. In agreement with the preceding conductivity and multifractality study of the surface theory (Secs. III and IV), we find the same behavior at finite energies for both winding numbers. In accordance with the expectation [Eq. (9)] for the topological AIII surface theory, there is strong multifractality for strong disorder at zero energy. For E > 0.5 , there is very good agreement between the class A QHPT parabola and the observed average multifractal spectra. At higher energies (not shown), the surface states start extending into the clean bulk.
The first three lines of Eq. (32) describe the topologically trivial "Gade" AIII Dirac model [18,70,77,86], defined by the restriction of Eq. (2) via the sublattice (chiral) symmetry in Eq. (3c). As we have seen in Sec. V and Fig. 7, the finite-energy states of this model are strongly Anderson localized.
A symmetry-based argument for the effect of ω = 0 in Eq. (32) is the following. With ω = 0, the model exhibits U(2n)⊗U(2n) symmetry, defined as invariance under independent left-and right-group transformations, After absorbing the matrixΛ by left-group translation, ω = 0 in Eq. (32) constrainsÛ L =Û R , reducing the symmetry to the diagonal U(2n). Although ω constitutes an imaginary mass, if we assume that the oscillations induced by fluctuations in the larger U(2n)⊗U(2n) space FIG. 12. Schematic "physical" renormalization group flow picture connecting zero-energy and finite-energy fixed-point descriptions of class AIII topological surface states. In the zero-energy horizontal plane, the system state is determined by the winding number ν [equal to π times the longitudinal conductivity, Eq. (14)] and by the abelian disorder strength λA, Eqs. (10) and (32). The latter parameterizes the parallel zero-energy fixed lines, indicated in black. The numerical results presented in Secs. III, IV, and VI strongly suggest that nonzero energy induces a flow (green arrows) to the logarithmic conformal field theory that governs the plateau transition of the ordinary quantum Hall effect (QHPT). For λA = 0, finite energy produces special behavior for winding numbers ν = 1 (clean Dirac) and ν = 2. The latter flows instead to the plateau transition of the spin quantum Hall plateau transition (SQHPT) [38]. Also indicated is the schematic flow to the QHPT critical point from the ultraviolet limit of the quantum Hall plateau transition, as captured by the Pruisken sigma model in Eq. (35) with a large bare conductivity σx,x 1 and half-integer σx,y.
The λ A and frequency terms have now disappeared due to Eq. (34). [We perform the group translation that replaceŝ Λ →1 2n in Eq. (32), before implementing the constraint.] The topological theta term on the second line of Eq. (35) has the coefficient σ x,y , given by Eq. (36). In the context of a single massless Dirac fermion, it has been argued that the parameter σ x,y is not the physical Hall conductivity [15]. The latter should vanish for the on-average parity-invariant class AIII surface [25].
The main point is that Eq. (36) implies a nontrivial even-odd effect. For odd values of the winding number ν, the Pruisken term in Eq. (35) has ϑ ≡ 2πσ x,y = π (modulo 2π), corresponding to the plateau transition of the integer quantum Hall effect [7,14,90]. On the other hand, the parameter ϑ = 2π for even winding numbers. This corresponds to the Anderson insulating quantum Hall plateau. We conclude that the "derivation" of Eq. (35) from Eq. (32) predicts that class AIII surface states at finite energy are critically delocalized at the QHPT (Anderson localized in the plateau) for odd (even) winding numbers. Although we have obtained Eqs. (35) and (36) using nonabelian bosonization, the same conclusion is seemingly implied by a straight-forward gradient expansion [5,14].
A key result of the numerical results presented in Secs. III, IV, and VI is that no such even-odd effect is observed. The Landauer conductivity and multifractal spectra results for surface theories with ν = 1, 2 shown in Figs. 2 and 4, as well as the results for the surface states of the class AIII bulk lattice model exhibited in Fig. 11, show no difference between even and odd winding numbers. In both cases, finite-energy states are observed to form a "stack" of quantum critical QHPT states.
We emphasize that Eqs. (35) and (36) do not obtain by following a physical renormalization group (RG) flow, starting from the ω = 0 perturbed WZNW model in Eq. (32). Instead, the imposition of the constraint in Eq. (34) is essentially a "mean field" guess for how the zero-energy field theory responds to the finite-energy perturbation. It is certainly possible that a physical RG flow connects Eq. (32) with ω = 0 in the ultraviolet to Eq. (35), but with σ x,y = 1/2 for all nonzero ν. At the same time, the Pruisken theory is also an ultraviolet description of quantum Hall physics, i.e. Eq. (35) can only be derived in a controlled fashion for large bare σ x,x 1 [7,90]. The Pruisken model flows to strong coupling σ x,x → O (1); for σ x,y = 1/2, this flow is presumed to terminate at the logarithmic conformal fixed point governing the QHPT.
A schematic RG flow diagram is shown in Fig. 12, connecting the zero-energy fixed lines of Eq. (32) (with ω = 0, parameterized by the purely marginal abelian disorder variance λ A ), to the putative QHPT critical point. A conjectured flow from the ultraviolet Pruisken model (with ϑ = π) is also indicated in this figure. In this case, Eq. (35) is the effective Euclidean spacetime field theory for the 1+1-D antiferromagnetic SU (2) Heisenberg spin-s chain, with σ x,y = s [91]. According to Haldane's conjecture, for half-integer s this model is equivalent to the SU(2) sector of the U(2) WZNW model in Eq. (32), at level ν = 1. Although the "ultraviolet" Pruisken theory has only diagonal SU(2) symmetry, the WZNW model is invariant under the larger chiral symmetry group SU(2)⊗SU (2). Enlarged chiral symmetry is a generic property of rational 2D conformal field theories. By contrast, spin chains with integer s have ϑ = 2π, and flow to a gapped paramagnetic phase: this is the famous even/odd effect predicted by Haldane's conjecture [91].
Although the replica limit n → 0 produces nonunitary, c = 0 logarithmic conformal field theories [24,56], one might hope that a version of the Haldane conjecture survives in the replica limit. The QHPT features a noncritical density of states, as well as an approximately parabolic multifractal spectrum τ (q), Eq. (8) with θ 0.25 [7,31]. We can apparently satisfy both of these conditions in the class AIII WZNW model [Eqs. (9), (11) and (12)] by choosing Very recently, Zirnbauer has argued that consideration of the Chalker-Coddington network model independently implies the existence of a class AIII current algebra with level ν = 4 [57]. Owing to the special way that the WZNW model emerges in the quantum Hall context, Zirnbauer argues that the mean conductance at the plateau transition is exactly half the predicted value given by Eq. (14), σ x,x = 2/π 0.64. This is not too different from previous numerics [Eq. (15)] [7,[31][32][33][34]36], and the value we find at finite energy for ν = 1, 2 class AIII topological surface states in Figs. 2, 4, and 6.
A very surprising aspect of Zirnbauer's proposal [57] is that the correlation length exponent ν QHPT (not to be confused with the class AIII winding number ν) for tuning away from the transition into the quantum Hall plateau is actually predicted to be infinite. The claim is that existing numerical studies showing ν QHPT in the range between 2.3-2.6 [7,31,[58][59][60][61][62][63][64][65][66] are beset by finite-size effects. We emphasize that for the finite-energy class AIII surface states that exhibit QHPT critical statistics studied here, the exponent ν QHPT does not play a role. This is because all states are delocalized, and the only relevant length scale is the crossover scale ζ(E) ∼ E −1/z between the low-energy class AIII states (which exhibit winding number ν-dependent statistics) and the finite-energy QHPT states. Here z is the dynamic critical exponent associated to the class AIII, low-energy states [Eq. (11)].
The numerical results obtained in this paper open up a potential avenue to test Zirnbauer's ν = 4 theory. In particular, if we assume that the ν = 4 class AIII WZNW model [Eq. (32)] also exhibits a "stack" of QHPT critical eigenstates at finite energy, then the only difference between zero-and finite-energies for ν = 4 is the finetuning of the abelian parameter λ A . The latter must satisfy Eq. (37) in order to realize QHPT multifractality at finite energy. Thus, both zero-and finite-energy states would be governed by the same quantum field theory, and the only missing piece is a controlled renormalization group mechanism that fine tunes λ A → π/16 for ω = 0 in Eq. (32). This is an interesting avenue for future analytical and/or numerical work.

C. Conclusion
In this paper we have presented substantial numerical evidence that most finite-energy surface states of a class AIII topological superconductor (or chiral topological insulator) form a "stack" of quantum critical, QHPT states, in the presence of surface disorder that preserves the defining class AIII symmetries. Similar results based on the multifractal analysis of surface states with disorder were obtained in Refs. [38,39], for class CI and DIII topological superconductors. For class AIII treated here, we have provided two additional sources of evidence for the stacking of quantum Hall plateau transition states: (1) Landauer conductance of the finite-energy states in the continuum surface-only theory, and (2) multifractal analysis of all in-gap states of a class AIII topological diamond lattice model in the slab geometry.
In total, the results here and in Refs. [38,39] strongly argue for an unexpected, new connection between Zgraded topological phases in two and three spatial dimensions. The 2D topological phases in classes C, A, and D describe the spin, charge (ordinary integer), and thermal quantum Hall effects, with broken time-reversal symmetry. The 3D topological phases in classes CI, AIII, and DIII can describe time-reversal invariant topological superconductors. The numerics suggest that the quantum critical plateau transitions in classes C, A, and D, which occur only at isolated energies in two dimensions, reappear in gap-spanning "stacks" at the surfaces of class CI, AIII, and DIII topological phases, respectively, in the presence of symmetry-preserving disorder.
Much of the physics of these newly uncovered quantum-critical, multifractal surface metals remains to be explored. The role of interparticle interactionmediated instabilities in the presence of gap-spanning multifractality is an obvious example. Another avenue is the investigation of Zirnbauer's proposal for the field theory of the plateau transition as a particular zero-energy class AIII topological surface state [57], see Sec. VII B for a discussion.
Classes CI, AIII, and DIII are also notable for exhibiting a precisely quantized longitudinal surface spin or thermal conductivity, which holds even in the presence of disorder and interactions [7,16,[28][29][30]. The connection between 2D class C, A, and D quantum Hall effects and 3D class CI, AIII, and DIII topological superconductors likely reflects a deep, topological origin. The latter remains to be uncovered in future studies.