A Discrete symmetries, random matrix classification, and disorder

# Multifractal nature of the surface local density of states in three-dimensional topological insulators with magnetic and nonmagnetic disorder

## Abstract

We compute the multifractal spectra associated to local density of states (LDOS) fluctuations due to weak quenched disorder, for a single Dirac fermion in two spatial dimensions. Our results are relevant to the surfaces of topological insulators such as and , where LDOS modulations can be directly probed via scanning tunneling microscopy. We find a qualitative difference in spectra obtained for magnetic versus non-magnetic disorder. Randomly polarized magnetic impurities induce quadratic multifractality at first order in the impurity density; by contrast, no operator exhibits multifractal scaling at this order for a non-magnetic impurity profile. For the time-reversal invariant case, we compute the first non-trivial multifractal correction, which appears at two loops (impurity density squared). We discuss spectral enhancement approaching the Dirac point due to renormalization, and we survey known results for the opposite limit of strong disorder.

###### pacs:
73.20.-r, 73.20.Jc, 64.60.al, 72.15.Rn

## I Introduction

The defining attribute of a 3D topological insulator(1) (TI) is the presence of an odd number of 2D massless Dirac bands at the material surface.(2); (3) Unlike the Dirac electrons that can appear in a purely 2D system (notably in graphene), the surface states of a (strong) 3D TI are robustly protected from the opening of gap, so long as time-reversal symmetry is preserved. The protection can be viewed as a consequence of the parity anomaly,(4); (5); (6); (3) which “holographically” links surface states separated by a topologically non-trivial bulk, and gives rise to the signature properties of the TI state: the half-integer quantum Hall effect, quantized magnetoelectric coupling, “axion” electrodynamics, etc.(2); (3) As stressed by Schnyder et al. in Ref. (7), the robust character of the surface states in the presence of quenched disorder can also be taken as a principal characteristic of a topological insulator. In particular, these states are protected from Anderson localization,(8) even in the presence of a “strong” impurity potential, so long as time-reversal invariance is preserved.(10); (9)

With its 2D Dirac band pinned to an exposed surface, a 3D TI is ideally suited to local probes such as scanning tunneling microscopy (STM). In spectroscopic mode, an STM captures an areal map of the local density of states (LDOS). There are several ways of analyzing such data. One is to look for quasiparticle interference (QPI)(11); (14); (15); (12); (13) in the LDOS Fourier transform. This method is useful for determining short-distance details, and contains similar information as an analysis of LDOS Friedel oscillations in the presence of a single impurity.(16) It has been applied in TIs to experimental data and analyzed theoretically in Refs. (12); (13) and (14); (15), respectively. In QPI, the disorder is employed primarily as a facilitator to gleam information about the clean system.(11)

Multifractal analysis(18); (19); (17) provides a complementary method better suited to extracting large-distance, disorder-dominated features in the same LDOS data field. It is a standard tool for assaying quantum interference phenomena, and is employed in the analysis of wavefunctions near a metal-insulator transition(20); (21); (18); (19) as well as mesoscopic fluctuations in diffusive metallic systems.(22); (23) In this paper, we derive new results for LDOS multifractal spectra associated to disordered topological insulator surface states. In particular, we extend the pioneering results of Ref. (24) to the generic cases of time-reversal () preserving and breaking impurities. Our calculations are performed in the near-ballistic limit,(25) wherein weak disorder enters as a perturbation to the clean Dirac band structure. A key characteristic of 2D Dirac fermions is that this weak disorder regime is continuously connected to more conventional domains of multifractal analysis, i.e. the diffusive (symplectic) metal(20); (23) and the integer quantum Hall plateau transition.(18); (26); (27); (28) These appear at strong coupling (many impurities) for dirty Dirac fermions.(24); (10); (9); (29)

We consider the case of a single flavor Dirac surface band, relevant to (e.g.) the TIs BiSe and BiTe (Refs. (2); (3); (30)). The different kinds of -preserving and -breaking disorder are sketched in Fig. 1. We demonstrate that the LDOS multifractal spectra observed in the absence of time-reversal symmetry breaking (i.e., for non-magnetic disorder) is qualitatively weaker than that induced by magnetic impurities. In particular, the first multifractal correction obtains at first order in the impurity density for the case of broken , while the first non-trivial amplitude appears at second order in the -invariant case. We compute the leading terms via one- and two-loop calculations, respectively. We also compute unnormalized spectra for the spin LDOS(31) in the case of magnetic impurities. We show that renormalization effects can enhance multifractality near the Dirac point. Finally, we summarize prior results on various strong-coupling regimes. Our goal is to sketch the full portrait of quantum interference physics on the surface of a TI, valid when interparticle interactions can be neglected.

Our results indicate that the long-distance, disorder-dominated features captured by the multifractal analysis behave in many cases opposite to the short-distance characteristics that appear in quasiparticle interference.(14); (15); (12); (13) In Ref. (14), the authors observed that QPI is strongest for the spin LDOS response to magnetic impurities, while the unpolarized LDOS pattern vanishes for magnetic disorder (in the first Born approximation). The QPI response of the LDOS to non-magnetic disorder is weak but non-zero.(14) By contrast, in this work we find that the LDOS multifractality is strongest for magnetic impurities, while the spin LDOS spectrum comparatively exhibits the same or weaker strength fluctuations, depending upon the polarization direction.

The weak influence of non-magnetic disorder is tied to the intrinsic spin-orbit coupling that defines the massless Dirac kinetic term. Multifractality is suppressed at one loop due to interference mediated by the Dirac pseudospin, which is proportional to the physical spin on a insulator surface. The spin is also responsible for the suppression of backscattering from a single non-magnetic impurity.(32) On the TI surface, magnetic disorder Zeeman couples directly to the Dirac spin, enabling backscattering in near-ballistic transport, and inducing multifractal LDOS fluctuations at the lowest order in the impurity density.

A notable problem in experiments probing topological insulator surface states has been the unintentional doping of carriers into the bulk bands, which then dominate transport measurements in large samples.(33) Even if the chemical potential is moved into the gap, it may reside far from the Dirac point, making it difficult to observe surface state carrier dynamics at low densities. In this respect, STM offers several advantages over transport experiments. First, the position of the chemical potential is no barrier to probing states at the Dirac point, since the latter can always be reached by tuning the bias voltage (although the Dirac point is not guaranteed to reside in the bulk gap).(3); (30) Assuming that the Dirac point or the low density regime can be accessed by tuning the tunneling bias, the advent of a finite, even large doping of the surface and/or bulk states may actually play a beneficial role in facilitating the observation of disorder-induced quantum interference effects. This is because a finite carrier density screens the long-range Coulomb potential introduced by charged defects. The potential landscape formed by screened impurities is short-range correlated on scales larger than the screening length. Good screening eliminates the problem of electron and hole puddle formation,(34); (35) which has until recently(36) occluded transport and other properties of Dirac carriers in graphene near the Dirac point. On the other hand, a low density of poorly-screened bulk dopants induces a long-range correlated potential and puddle formation, as in graphene.(13) LDOS fluctuations in the puddle regime are an important topic for future work.

Three-dimensional topological insulators provide us with an interesting paradigm flip for quantum interference phenomena. Isolating the surface state contribution in transport measurements is problematic. By comparison, direct LDOS imaging is easier than in conventional semiconductor systems, wherein the 2D electron gas is typically buried in a layered material stack. Moreover, the amount of surface disorder can to some extent be controlled; for example, magnetic impurities can be deposited across the surface of an otherwise high-quality bulk 3D sample. These can be charge-neutral adatoms or charged dopants; an example of the former (latter) is provided by iron (manganese)(37) in .

This paper is organized as follows. We begin in Sec. II with a lightning review of multifractal composite and spin LDOS measures. In Sec. III, we present new results for multifractal LDOS fluctuations in TI surface states, in the presence of weak disorder. We also show how renormalization can enhance multifractality close to the Dirac point. Finally, in Sec. IV, we review previous results on various strong disorder regimes relevant to the TI surface states and LDOS statistics. In particular, we discuss the symplectic metal, the integer quantum Hall plateau transition, and the Anderson insulator. Various technical details are relegated to appendices. In Appendix A, we review the symmetry classes of Anderson (de)localization that appear in the disordered Dirac surface theory. In Appendix B, we supply some details of our perturbative calculations.

## Ii Multifractal LDOS measures

### ii.1 Definitions

We suppose that the tunneling local density of states (LDOS) is imaged at a fixed energy over an field of view. The field is finely partitioned into a grid of boxes. The box edge length must be chosen larger than any “microscopic” scale , such as the correlation length of the random potential.(18) One introduces the box probability

 μn(ε)≡∫And2rν(ε,r)∑l[∫Ald2rν(ε,r)], (1)

where denotes the box. LDOS multifractality is defined through the inverse of the participation ratio (IPR),(20); (18)

 Pq(ε)≡∑nμqn(ε)∼(aL)τ(q,ε). (2)

The right-hand side (scaling limit) obtains when ; corrections are down by higher powers of . The exponent is the multifractal moment spectrum(38); (18) for LDOS fluctuations at energy .

The construction in Eqs. (1) and (2) is useful for characterizing a system with extended states, or for an Anderson localized system in which ; denotes the localization length. In what follows, we assume experiments are performed at sufficiently low temperatures so that inelastic cutoffs to quantum interference can be ignored.(8); (39) A clean system with plane wave states at energy has . Multifractality refers to the incorporation of corrections non-linear in . Physically, these arise due to quantum interference via multiple scattering of electron waves in a dirty environment, processes that serve as the precursor to Anderson localization.(20); (18); (22)

For weak disorder, the spectrum is typically dominated by the quadratic correction(20); (18)

 τ(q,ε)=2(q−1)−θ(ε)q(q−1), (3)

where gives a measure of the disorder strength. As an example, in a weakly disordered 2D metal (with for the orthogonal or unitary classes), one finds(20); (21); (23)

 θ=β−12π2N(ε)D, (4)

where denotes the average density of states, is the classical (Drude) diffusion constant, and , depending upon the presence or absence of time-reversal symmetry and spin-orbit scattering.(23); (40) At stronger disorder, higher order corrections in must be retained; for the diffusive metals, results are known to four loops.(20)

An alternative characterization of LDOS multifractality is provided by the singularity spectrum(38); (18) : Over a subset of the sample grid area that scales as , the box probability . The singularity spectrum is the Legendre transform of ,

 f(α)=qα−τ(q),dτ(q)dq=α.

For the quadratic spectrum in Eq. (3), one obtains

 f(α)=2−14θ(α−2−θ)2. (5)

In this “parabolic approximation,” the strength of the multifractality is encoded in the peak position [], and the width of the spectrum such that ,

 α0=2+θ,αW=4√2θ. (6)

Part of the power of multifractal analysis for disordered quantum systems derives from the fact that the spectra [ or ] typically depend only upon a few gross measures of the impurity potential. In the case of dirty metals, the entire spectrum can be computed as an expansion in one parameter, the inverse conductance (consistent with scaling theory).(20); (8); (22) At a non-interacting Anderson localization transition, and become universal functions, so that the critical point is characterized by an infinite set of critical exponents [e.g., the expansion coefficients for ].

The spectra above have been defined for data collected in a single fixed realization of the disorder. Strictly speaking, Eq. (3) then applies only for , where . Outside of this range, the associated to a fixed disorder realization is linear, a phenomenon known as spectral termination.(41); (42); (43) [This assumes that higher order corrections can be ignored for . Regardless, the spectrum is always linear for sufficiently large ]. Termination can be viewed as a consequence of the restriction to positive sample measures .(38); (19)

In the localized regime, the states contributing to the LDOS at a given position in the sample have a discrete energy spectrum, quantized by the typical localization volume . As a result, all non-unity LDOS moments diverge in the absence of level smearing. In a tunneling experiment, smearing can appear due to inelastic scattering (temperature), open sample boundary conditions, or due to the finite energy resolution of the instrument. To characterize an Anderson insulating state over an field of view with , the full LDOS distribution should be examined;(44); (45) sensitive dependence of the distribution shape to smearing can serve as a telltale sign of the localized regime. LDOS fluctuations in the Anderson insulator are reviewed in more detail in Sec. IV.2.

### ii.2 Spin LDOS spectra

By restricting the character of the tunneling species, it may be possible to measure individual LDOS components separately. For example, in the case of a spin-polarized (ferromagnetic) STM tip,(31) the spin-projected components can be separately resolved. The use of an unpolarized tip recovers the composite LDOS . We define the spin LDOS along the spin space direction ,

 ν^ι(ε,r)≡ν^ι↑(ε,r)−ν^ι↓(ε,r). (7)

For a time-reversal invariant system (with or without spin-orbit scattering and/or disorder), one has . In a system with broken time-reversal (e.g., magnetic impurities), but zero average spin polarization, the integral of over a sufficiently large region becomes arbitrarily small; we cannot use the normalized construction in Eqs. (1) and (2) to characterize spin LDOS multifractals. Instead, we employ the un-normalized inverse spin participation ratio (ISPR)

 P^ιq(ε)≡∑n(μ^ιn)q,μ^ιn≡∫And2rν^ι(ε,r). (8)

In the scaling limit,

 P^ιq(ε)∼cq(aL)x^ιq−2, (9)

where the exponent is the scaling dimension for the corresponding moment operator in the disorder-averaged field theory description, and for even .

## Iii Weak disorder multifractality

### iii.1 Model. Short- and long-range correlated potential landscapes

The Dirac surface states of a topological insulator (TI) are guaranteed to appear in an odd number of flavors.(2); (3) In this paper, we consider the simplest case of a single flavor, relevant to (e.g.) and . The Hamiltonian is (in units such that )

 Unknown environment '% (10)

where , and repeated indices are summed. The coordinates chart the TI surface, while the topological bulk resides in the perpendicular direction. In Eq. (10), denotes the Fermi velocity, and the Dirac pseudospin Pauli matrices are related to the physical spin operators via . The vector, scalar, and mass potentials describe the effects of external electromagnetic fields and/or surface impurities. In the absence of time-reversal () symmetry breaking, . (See Appendix A for an enumeration of discrete symmetry operations.) Thus, a mass gap is explicitly forbidden so long as remains a good symmetry, a consequence of the protection afforded by the topologically nontrivial bulk. When is broken by an external magnetic field , the vector and mass potentials are

 Aμ =−eA(orb)μ−γ∥2vFϵμνB∥,ν, M =−γ⊥2Bz, (11)

where () denotes the Zeeman coupling to the in-plane field (out-of-plane field ), and the orbital effect is embedded in via .

Non-magnetic adatoms or charge traps are encoded in the scalar potential . In-plane (out-of-plane) polarized magnetic impurities additionally induce point exchange coupling to the vector [mass ] fields.(16) The different types of disorder leading to , , and are sketched in Fig. 1. Assuming a random surface distribution of impurities and spatial rotational invariance on average, the disorder potentials can be taken as Gaussian white noise distributed variables,

 Misplaced & (12)

The dimensionless variances quantify the disorder strength. In the first Born approximation, these are of the form

 Δv2F=nimp|~u(0)|2, (13)

where is the impurity density, and denotes the Fourier transform of the single impurity potential. We note that a net in-plane magnetization of the surface impurities can be removed by a gauge transformation, while the average scalar potential is absorbed into the chemical potential. We will assume that there is no net magnetization perpendicular to the surface, , or that we only probe LDOS fluctuations on energy scales much larger than the induced gap .

In 2D, the single impurity potential [Eq. (13)] must decay faster than (or oscillate rapidly enough) so that the limit exists; otherwise, the white noise assumption in Eq. (12) is invalidated by long range impurity potential correlations.(46) This causes a problem for charged impurities, which can become poorly screened for a small surface doping relative to the Dirac point. In graphene, the long-range correlated potential undulations induced by poorly-screened substrate impurities leads to a smearing of the Dirac point over an energy scale , and to the breakup of the sample into electron and hole puddles.(34); (35) The advent of electron-hole puddles has until recently prevented the observation of various “intrinsic” phenomena associated to the Dirac carriers in graphene experiments such as velocity renormalization(36) and hydrodynamic transport near the Dirac point. In this respect, a large surface or bulk doping actually improves the situation for STM measurement of disorder-induced quantum interference, since these carriers screen the potential of surface charges. The disorder potential can be considered short-range correlated for scales larger than the screening length.

If we consider only surface doping, with an insulating bulk, then the Thomas-Fermi wavelength due to a finite surface carrier density is given by

 λTF=1α√πn, (14)

where is the effective “fine structure constant.” The permittivity , the average of the bulk TI below and vacuum above the surface. For with a surface density of (corresponding to a doping level of eV relative to the Dirac point),(30) m/s (Ref. (30)), and permittivity(47) , one obtains nm. This is very large, and indicates that the surface state carrier density is inadequate to screen charged impurities. A smaller screening length is possible for bulk doping,(13) or by performing experiments on thin film samples exfoliated over a metallic gate. Alternatively, one can restrict the deposition of surface impurities to non-doping adatoms, e.g. iron in .(37) The disorder variance associated to Thomas-Fermi screened charged impurities is

 ΔV=πnimpn. (15)

Finally, we note that the appearance in isolation of any of the three disorder potentials in Eq. (10) realizes three different symmetry classes of Anderson (de)localization,(48); (49); (7) see Appendix A for a review. The -invariant case with belongs to the spin-orbit class AII, which is also the class of the topological bulk [Fig. 1(a)]. In the case of broken , realizes the random vector potential model in class AIII [Fig. 1(b)], while gives the random mass model in class D [Fig. 1(c)]. All three classes exhibit delocalized states in 2D, although this occurs only at the Dirac point for class AIII.(24) In the -invariant symplectic case, the unpaired single Dirac flavor avoids the usual spin-orbit metal-insulator transition,(9) remaining delocalized even for strong disorder due to a topological term.(10) The generic case of broken- with all three disorder potentials non-zero realizes the unitary class A, and is believed to flow under renormalization to the plateau transition in the integer quantum Hall effect.(24); (29) (See Sec. IV.1.2 for a review).

Because in-plane (out-of-plane) Zeeman coupling appears in the vector (mass) potential [Eq. (III.1)], one is tempted to identify class AIII (class D) with the limit of an otherwise clean surface, dusted with charge neutral magnetic impurities randomly polarized in-plane (perpendicular to the TI surface). However, a magnetic adatom is expected to also induce a local scalar potential deformation . For example, it can dope the surface or bulk, as occurs for a manganese impurity in (Ref. (37))]. As discussed in Appendix A, the advent of any two flavors of disorder destroys the additional discrete symmetries enjoyed by the special class D and AIII Hamiltonians. The asymptotic long-distance LDOS scaling is then governed by the unitary class A, discussed above. Nevertheless, depending upon the relative microscopic strength of the magnetic versus potential perturbations induced by polarized magnetic impurities, the class AIII or D model may provide an adequate approximation for broken- LDOS fluctuations on intermediate scales.

### iii.2 Results

To compute the scaling of LDOS moments in a quantum theory with quenched disorder, one employs a path integral to express products of fermion Green’s functions as functionals of the disorder configuration. Using a trick (replicas,(20); (8); (22) supersymmetry,(50) or Keldysh(51)) to normalize , the Green’s functions are formally averaged over disorder configurations (typically with a Gaussian weight). The result is a translationally-invariant, but “interacting” field theory, where the disorder strength appears as a coupling constant.(8); (50) Perturbative calculations are controlled via loop expansion for small .

To determine the scaling, one decomposes the LDOS moment into projections upon the renormalization group (RG) eigenoperators of the disorder-averaged theory.(20); (21); (22) The multifractal spectrum is determined by the most relevant (negative)(52) scaling dimension exhibited by an eigenoperator in this decomposition, and is given by(19); (43)

 τ(q)=2(q−1)+xq−qx1. (16)

#### Broken T: random vector potential disorder (Class AIII)

The properties of the model in Eq. (10) with short-range correlated disorder [Eq. (12)] were originally studied in Ref. (24). In this work, the exact multifractal spectrum was calculated for the broken-, random vector potential ( in-plane polarized magnetic impurity)(53) class AIII model, to all orders in . Technically, this result obtains because the disorder-averaged AIII model is conformally invariant at the Dirac point, and the exact LDOS moment spectra can be extracted using an Abelian bosonization treatment. The exact spectrum(24) is quadratic in , and takes the form of Eq. (3), with

 θA=ΔAπ. (17)

Subsequent work(41); (42) on the random vector potential model elucidated the mechanisms of termination and freezing, transitions that occur in the spectral statistics for large moments or strong disorder .

For this broken- class, we can also examine the spin LDOS fluctuations, utilizing the same nonperturbative bosonization treatment employed in Ref. (24). The spin LDOS taken along an axis in spin space was defined by Eq. (7). Moment fluctuations are characterized by the inverse spin participation ratio (ISPR) in Eq. (8), the scaling limit of which is controlled by the dimension that appears in Eq. (9). The out-of-plane ISPR is associated to the “mass” fermion bilinear . For the random vector potential model, the most relevant contribution to carries the same scaling dimension that gives the composite LDOS scaling in Eqs. (3) and (17),

 x^3q=q−ΔAπq2. (18)

The chiral components of the in-plane spin LDOS are the energy-resolved Dirac current operators

 ν±≡ν^1±iν^2=ψ†^σ±ψ. (19)

Moments of these are RG eigenoperators that receive no corrections. The scaling of the associated ISPR is governed by the disorder-independent (tree level) exponent

 x±q=q. (20)

Eqs. (17), (18), and (20) are exact results that hold to all orders in .

#### Broken T: random mass disorder (Class D)

In the rest of this section, we provide new results for the broken-, random mass ( out-of-plane polarized magnetic impurity)(53) class D model, the -invariant class AII model, and the generic broken- unitary class A model. For weak disorder, none of these are conformally invariant, and we resort to perturbation theory. In this section we summarize results; some technical aspects are sketched in Appendix B. The results obtained below hold only for small , wherein the disorder appears as a weak marginal perturbation (at tree level) to the clean Dirac surface band structure.

For the broken- case of random mass disorder (with ), one obtains quadratic multifractality at one loop, again governed by Eq. (3), with

 θM=ΔM2π+O(Δ2M). (21)

Moments of the out-of-plane spin LDOS operator , as well as of the chiral in-plane [ current] operators constitute RG eigenoperators at one loop, with scaling dimensions

 x^3q= q+ΔM2πq+O(Δ2M), (22) x±q= q+O(Δ2M). (23)

Note that the first correction in Eq. (22) is positive (and linear in ); this should be contrasted with the AIII case, Eq. (18) above. On general grounds, the anomalous scaling dimension associated to the moment of the composite LDOS, or any projected component thereof, must appear with a negative sign. The reason is that this quantity is associated to a moment of a normalized probability distribution(18); (52) through Eqs. (1) and (2). For a quadratic spectrum, this leads in particular to in Eq. (3) [consistent with a positive, real disorder variance—c.f. Eqs. (17), (21), and (24)]. By contrast, the spin LDOS is defined as the difference between two orthogonal projections [Eq. (7)]; for this reason, the first disorder correction to the scaling dimension in Eq. (22) is not required to appear with a particular sign.

#### Non-magnetic disorder (Class AII)

In the -invariant case of scalar potential disorder, it turns out that no local operator (without derivatives) exhibits multifractal scaling to first order in . For Dirac fermions, this applies to both LDOS and energy-resolved current moments. Physically, the weak influence of non-magnetic disorder is due to interference mediated by the Dirac pseudospin (equivalent to physical spin on the TI surface). The Dirac pseudospin is also responsible for the suppression of backscattering from a single non-magnetic impurity.(32) Technically, this result is derived by mapping the one-loop renormalization process of local operators to the action of a certain spin- Hamiltonian , and identifying renormalization group eigenoperators with states that diagonalize (see Appendix B). As a result, to lowest order one observes plane wave scaling in the LDOS IPR [Eq. (2)]. The spin LDOS vanishes exactly, due to .

The first non-trivial correction to the LDOS appears at two loops. To this order, the spectrum is again quadratic as in Eq. (3). A straight-forward but laborious calculation gives the coefficient in this equation,

 θV=3Δ2V8π2+O(Δ3V). (24)

Since [Eqs. (13) or (15)], we find that the non-trivial multifractal scaling begins at second order in the impurity density. This is qualitatively weaker than any of the broken- regimes, where the quadratic multifractality appears already at first order, Eqs. (17) and (21). This distinction between -invariant and -broken surfaces is our primary result, and can be tested directly in STM experiments by varying the concentration of deposited surface disorder. Although the -invariant case is not conformally invariant (for a discussion of renormalization effects, see Sec. III.3, below), the multifractal and spectra depend only upon a single parameter, the variance . Eq. (24) can be extended to higher loops, allowing ever more precise tests against numerics or experimental data within the perturbatively accessible regime. The multifractal spectrum therefore provides a unique fingerprint for the time-reversal invariant Dirac surface state of the topological insulator, in the presence of weak but otherwise generic non-magnetic disorder. The opposite limit of strong disorder for the -invariant case is discussed below in Sec. IV.1.1.

#### Broken T: generic disorder (Class A)

When is broken and any two disorder flavors appear, the system resides in the unitary class A. The third disorder flavor is always generated under renormalization—see Sec. III.3, below. The results of Eqs. (17) and (21) for the LDOS spectrum in the random vector and mass potential models suggest that the unitary case also exhibits multifractality to first order in the impurity density , since .

With multiple flavors of the disorder, solving the operator mixing problem for the LDOS moment requires the diagonalization of an effective spin Hamiltonian , transcribed in Eq. (53) of Appendix B. In Figs. 3 and 4, we present the results obtained by numerically diagonalizing this matrix for various combinations of . In these figures we plot the Renyi dimension(17) , defined for via

 dq≡τ(q)q−1. (25)

Figs. 3 and 4 show that the generic broken- case is multifractal at one loop, and easily distinguished from the two-loop -invariant result, in the limit of weak disorder. [Note that Fig. 4 indicates that the spectrum is not purely quadratic in this general case.] It should therefore be possible to precisely distinguish the broken- and -invariant spectra experimentally, by observing the dependence of the deviation on . The single-disorder flavor results for comparable strengths are plotted in Fig. 2 for reference.

For the multidisorder unitary model, the same RG eigenoperators dominate the scaling of composite and out-of-plane spin LDOS moments. The dimension that determines the spin LDOS scaling via Eq. (9) also enters into the LDOS spectrum in Eq. (16), leading to Figs. 3 and 4. By contrast, moments of the chiral spin LDOS components [Eq. (19)] remain eigenoperators that acquire no corrections at one loop,

 x±q=q+O(ΔαΔβ), (26)

.

As reviewed in the subsequent Sec. IV.1.2, for , the generic broken- model is believed to flow to the critical state at the integer quantum Hall plateau transition.(24); (29) This state exhibits strong multifractality that has been extensively studied in numerics.(18); (19); (26); (27); (28)

### iii.3 Renormalization effects

As discussed at the beginning of the previous section, the disorder-averaged Dirac surface state theory used to compute LDOS multifractal spectra is an “interacting” field theory, wherein the disorder strengths appear as coupling constants (c.f. Appendix B). Because these parameters are dimensionless, at weak coupling the disorder constitutes a marginal perturbation of the clean Dirac band structure. The one-loop RG equations for these parameters are given by(24); (54)

 dΔAdl =1πΔMΔV, (27a) dΔMdl =1π(2ΔA−ΔM)(ΔM+ΔV), (27b) dΔVdl =1π(2ΔA+ΔV)(ΔM+ΔV), (27c)

where denotes the log of the RG length scale (e.g., the system size). Energy scales as

 dlnεdl=z(l), (28)

where the (scale-dependent) dynamic critical exponent is

 z=1+12π(2ΔA+ΔM+ΔV)+O(ΔαΔβ), (29)

.

In this section, we use Eqs. (27)–(29) to derive the dynamical scaling of the disorder parameters ; here energy is measured relative to the Dirac point, not the Fermi energy. (From the point-of-view of the disordered Dirac theory, a finite energy above the Dirac point constitutes a relevant perturbation.(24)) Using the results obtained in the previous section, we thereby determine the enhancement or suppression of LDOS multifractality approaching the Dirac point, due to renormalization.

#### Broken T: random vector potential disorder (Class AIII)

For the random vector potential model with , Eq. (27) implies

so that (constant), where is the “microscopic” value derived from the randomly polarized in-plane magnetic impurity distribution.(53) This result in fact holds to all orders in ;(24) in this case, the theory describing LDOS fluctuations at the Dirac point is conformally invariant. Multifractality is neither enhanced nor suppressed as one moves away from the Dirac point, defined as . However, for non-zero energies , in an infinite size sample all states are in fact localized.(24) The localization length diverges upon approaching the band center as , with [Eq. (29)]. Eqs. (3) and (17) for hold on scales smaller than .

#### Broken T: random mass disorder (Class D)

For the random mass model with ,

 dΔMdl=−Δ2Mπ+O(Δ3M),

so that the disorder is marginally irrelevant at weak coupling.(24) Integrating this equation and using Eqs. (28) and (29), we can compute the scaling of with energy. At energy scale , we define ; then for the smaller energy scale (relative to the Dirac point), we obtain the logarithmic suppression

 ΔM(ε≲Υ)∼ Δ(∘)M−(Δ(∘)M)2πlog(Υε) +O{[Δ(∘)M(1−εΥ)]2,(Δ(∘)M)3}. (30)

This equation holds for . In the limit as , the disorder strength vanishes as

 ΔM(ε→0)∼ π⎡⎢⎣ln⎛⎜⎝√πΔ(∘)MΥε⎞⎟⎠⎤⎥⎦−1 +O⎡⎢⎣1Δ(∘)Mln−2⎛⎜⎝√πΔ(∘)MΥε⎞⎟⎠⎤⎥⎦. (31)

For small , Eq. (III.3.2) applies only at very small energies .

#### Non-magnetic disorder (Class AII)

Now we consider the -invariant model. The flow equation for is

 dΔVdl=Δ2Vπ+O(Δ3V). (32)

In contrast to the random mass, the random scalar potential is a marginally relevant perturbation to the clean band structure.(24) Examining lower and lower energy scales approaching the Dirac point, one observes stronger effects of the disorder. In the asymptotic scaling limit wherein the impurity potential strength becomes “large” (), numerical(9) results and analytical(10) arguments imply that the disordered -invariant Dirac theory renormalizes into the “conventional” symplectic metal. The metal is distinguished from the Dirac theory by its non-zero (and non-critical) density of states at zero energy,(8) and by its spectrum.(20); (23) We discuss the strong coupling LDOS multifractality below in Sec. IV.1.1. If at energy , , then for a somewhat smaller energy we obtain the logarithmic enhancement

 ΔV(ε≲Υ)∼ Δ(∘)V+(Δ(∘)V)2πlog(Υε) +O{[Δ(∘)V(1−εΥ)]2,(Δ(∘)V)3}. (33)

Eq. (III.3.3) implies that renormalization strengthens multifractality approaching the Dirac point , for the -invariant case. We emphasize that this has nothing to do with weak (anti-)localization. The latter occurs in the diffusive metallic regime with , where denotes the elastic mean free path. The diffusive regime obtains at strong coupling(9) near the Dirac point [Sec. IV.1.1, below]. The impurity strength renormalization in Eq. (III.3.3) is a quantum effect deriving from the clean band structure, in the “near-ballistic” regime.(25)

#### Broken T: generic disorder (Class A)

In the generic case of broken , with multiple disorder coupling strengths non-zero, the system flows toward strong coupling . As a result, multifractality is enhanced approaching the Dirac point. The RG flow ultimately terminates at a strong coupling critical point, or in the Anderson insulator, discussed in the next section.

## Iv Strong coupling regimes

In this section, we review prior results on strong coupling regimes relevant to the disordered Dirac topological insulator surface states, LDOS fluctuations and associated multifractal spectra. These are not new, but provide complimentary information to the new results derived in the previous section.

In both generic cases of -invariant, and -breaking impurities, the disordered Dirac description used in Sec. III fails on the largest length and lowest energy scales (approaching charge neutrality). For a sufficiently dilute concentration impurities, the results obtained in the previous section characterize the start of the scaling regime, over energy and length scales such that the disorder strengths remain weak . When these parameters become order one (due to renormalization down to lower energies and longer lengths), the system crosses over to one of the strong coupling regimes discussed below.

### iv.1 Delocalized states at strong disorder

#### T-invariant case: diffusive metal via strong disorder

In a random scalar potential field, the Dirac point vacillates in energy with spatial location; as a result, the density of states near charge neutrality is enhanced by the disorder. Due to the suppression of pure backscattering for Dirac fermions,(32) the state density enhancement more than compensates for the increased scattering introduced by the additional impurities. As a result, scalar potential disorder actually increases the (zero temperature, Landauer) conductance at charge neutrality beyond the clean ballistic result, (Refs. (9); (25)). As in Sec. III, here we assume short-range correlated disorder, due either to charge neutral impurities or efficient screening by bulk and/or surface carriers. We do not discuss the puddle regime(34); (35) in the present paper.

The effective disorder strength is enhanced by renormalization, as indicated by the runaway flow implied by Eq. (32). The concomitant density of states and conductance growth suggests that the disordered Dirac theory ultimately crosses over to the ordinary diffusive symplectic metal, a result born out by numerics.(9) In the absence of time-reversal symmetry breaking, Anderson localization is prohibited on the surface of a topological insulator.(7); (10)

The symplectic metal possesses a finite (non-critical) average density of states at charge neutrality, and a distinct spectrum. For a large effective diffusion constant (induced for a Dirac fermion subject to sufficiently strong disorder,(9) or for a weakly disordered system examined on large length scales), the lowest order result for the multifractal spectrum appears in Eqs. (3) and Eq. (4), above. In the latter equation, for the symplectic class.(23); (40)

For the -invariant case, the strongest multifractality is expected at intermediate coupling. Weak disorder induces weak multifractality in the Dirac language [ in Eqs. (3) and (24)], while strong disorder ultimately pushes the system into the symplectic metal, where a large diffusion constant suppresses the first correction in Eq. (4).

#### Broken T: IQHP transition

For generic -breaking disorder, i.e. all three non-zero, the disordered Dirac theory is also unstable under renormalization. When the average mass is zero (see below), the flow in Eq. (27) is believed to terminate at the critical point of the integer quantum Hall plateau transition.(24); (29) This is the delocalized state separating adjacent Hall plateaux; it exhibits strong multifractality that has been extensively studied in numerics.(18); (26); (27); (28) The spectrum is believed to be universal,(18) and is approximately(28) parabolic as in Eqs. (3) and (5), with (Refs. (27); (28)).

### iv.2 Anderson insulator

At zero chemical potential relative to the Dirac point, an average out-of-plane spin magnetization at the surface of a TI corresponds to the presence of a non-zero Dirac mass for the surface carriers. This insulating state resides in a quantum Hall plateau [with ].(6); (24); (2); (3) In the presence of surface disorder, the plateau state will assume the character of a localized Anderson insulator. In this section we review LDOS fluctuations in the Anderson insulator. The discussion is relevant not only to the magnetized surface of a 3D TI, but also to localized states populating the bulk gap of a disordered TI. Proposals exploiting localization to realize so-called “topological Anderson insulators” by adding impurities to clean hosts include those in Ref. (55).

To understand local density of states fluctuations in an Anderson insulator, it is useful to first study a toy problem. Consider a tight-binding model on a dimensional lattice, subject to nearest-neighbor hopping , and a random on-site potential , distributed uniformly over the region . We assume the absence of spatial correlations in the disorder potential. The inverse relative strength of the disorder is measured by the ratio . We consider first the extreme limit of zero hopping, . In that case, the LDOS is the on-site operator

 νi(ε,Vi)=η/π(ε−Vi)2+η2,

where denotes an energy-smearing parameter. Physically, smearing is determined by inelastic scattering, open sample boundary conditions, or due to the finite energy resolution of the probing instrument.

At the “band center” , the distribution function for disorder-averaged LDOS moments evaluates to

 p(ν)≡ ∫W/2−W/2dVWδ[ν−ν(ε,V)] = 1πν2W√ννmax−ν. (34)

In this equation, the LDOS is constrained to the interval , where

 νmin=4ηπ(W2+4η2),νmax=1πη. (35)

Using Eq. (IV.2), one can compute the disorder-averaged moments of the LDOS,

 ¯¯¯¯¯νq=Γ(q−12)Wπq−1/2Γ(q)η1−q. (36)

The average LDOS is ; all higher moments are proportional to , and thus diverge in the limit of zero energy smearing . This not surprising, because the energy spectrum in our trivial toy model is discrete, so that the LDOS operator becomes a delta function with ill-defined moments as . The moments are dominated by the power-law (Pareto) tail of the distribution, accumulating at the upper limit . By contrast, the typical LDOS, defined as is dominated by the infrared

 νtyp=4ηe2πW2.

This vanishes in the limit .

We see that observables exhibit broad statistics in the single site model, governed by the power-law distribution in Eq. (IV.2). The moments are rendered finite only by the non-zero energy smearing parameter . This should be compared to the LDOS statistics in a system with extended states and weak multifractality, e.g. that characterized by the quadratic spectrum in Eq. (3), with . It is known(22) that the corresponding LDOS distribution has a Gaussian bulk, with a small amplitude log-normal tail responsible for the weak multifractality. For the metallic system, the result is independent of energy smearing, provided that the thermodynamic limit is taken before the smearing is set to zero. Returning to the toy insulator model, we observe that the global density of states (GDOS) is self-averaging in the same limit. The GDOS is defined via

 νG≡1NN∑i=1νi(Vi),

where denotes the number of sites. In the large -limit, the cumulant expansion can be evaluated via the saddle-point. The cumulants of the GDOS then take the form

 ¯¯¯¯¯¯¯¯¯¯¯[νG]qc=N1−q(¯¯¯¯¯νq+…),

where denotes the cumulant, and the omitted terms are smaller by positive powers of . Taking the infinite system size limit before sending the energy smearing to zero leads to the vanishing of all GDOS cumulants.

The calculations above can be extended to non-zero hopping via a locator expansion in small , as performed by Anderson in his original 1958 paper.(56) This expansion can be formally summed to all orders in 1D and on the Bethe lattice,(57) but an explicit solution for the LDOS statistics is difficult to obtain this way; see Ref. (50) for an alternative approach.

Altshuler and Prigodin(44) succeeded in deriving the distribution generating disorder-averaged LDOS moments in a 1D system, which is exponentially localized for arbitrarily weak disorder.(8) In the thermodynamic limit for a closed sample, they obtain the “inverse Gaussian” distribution

 p(~ν)=√4ηπϵ1~ν3/2exp[−4ηϵ(~ν−1)2~ν], (37)

where , and is the typical energy level spacing in a localization volume; is also the elastic scattering lifetime.(44) In the limit of small smearing , this distribution has moments

 ¯¯¯¯¯~νq =41−qΓ(q−12)√π(ηϵ)1−q. (38)

The exact result in Eq. (38) for the 1D Anderson insulator exhibits the same singular dependence on the energy smearing as the single site model moment in Eq. (36). In fact, the distributions in Eqs. (IV.2) and (37) are very similar: both feature a power law at intermediate , while the exponential factor in Eq. (37) plays the role of the hard cutoffs in Eqs. (IV.2) and (35). The close resemblance of the exact 1D and single site model results can be attributed to the discrete spectrum of energy levels contributing to the LDOS in an Anderson insulator, with an energy level spacing determined by the localization volume in spatial dimensions.

The take away is that the LDOS distribution in an Anderson insulator becomes very broad, with a power-law tail yielding divergent moments, in the limit of vanishingly small energy smearing. In an STM experiment performed at ultra-low temperature, on a large, isolated Anderson localized surface, the collected LDOS statistics should be very sensitive to the smearing induced by the energy resolution of the measurement itself.

The locally discrete energy spectrum of the LDOS in the Anderson insulator invalidates the use of Eq. (2) as a tool to compute the multifractal spectrum. As advocated above, the shape of the LDOS distribution function and its sensitivity to smearing can best reveal the insulating phase. If one insists upon computing moments, one must employ(20)

 τ(loc)(q) ≡−ddlnL¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ln[1¯¯¯ν∫Ldddr∑i|ψi(r)|2qδ(ε−εi)]. (39)

Since the levels are discrete,

 limη→0(πη)q−1νq(ε,r)=∑i|ψi(r)|2qδ(ε−εi). (40)

Thus,

 τ(loc)(q)= −ddlnLln{∫Ldddr[limη→0(πη)q−1¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯νq(ε,r)]} +ddlnLln{∫Ldddr¯¯¯¯¯¯¯¯¯¯¯¯¯¯ν(ε,r)}. (41)

In this equation, we replace averages-of-the-logs with logs-of-the-average, a procedure that is legitimate here because spatial and disorder-averaging are expected to yield the same results on the insulating side. Noting that the LDOS moments are -independent in the insulator for , we obtain the expected result(18) for localized states

 τ(loc)(q)=0, (42)

computed in a well-defined limit.

## V Acknowledgments

The author thanks Kostya Kechedzhi for a collaboration that lead to this work, and thanks Andreas Ludwig, Igor Aleiner, Pedram Roushan, Haim Beidenkopf, Emil Yuzbashyan, and Deepak Iyer for useful discussions. The author acknowledges support by the National Science Foundation under Grant No. DMR-0547769, and by the David and Lucile Packard Foundation.

## Appendix A Discrete symmetries, random matrix classification, and disorder

The 10 symmetry classes of disordered Hamiltonians (Hermitian random matrices) can be efficiently distinguished by the presence or absence of time-reversal , particle-hole , and chiral/“sublattice” symmetry .(48); (49); (7) For the two-component Dirac Hamiltonian in Eq. (10), the definitions of these symmetries are essentially unique. In terms of the two-component Dirac spinor , these appear as

 T: ψ →−i^σ2ψ,i→−i (43a) P: ψ →^σ1[ψ†]T, (43b) C: ψ →^σ3[ψ†]T,i→−i. (43c)

In the second quantized language, and are antiunitary transformations; the unitary can be taken as the product of these.

The imposition of any one of the discrete symmetries upon the Hamiltonian in Eq. (10) in every disorder realization restricts its form, and selects a particular random matrix symmetry class.(24); (48); (49); (7) (1) -invariance: , only potential disorder is allowed. Since , this is the symplectic (spin-orbit) class AII, which is also the symmetry class of the (presumed -invariant) topological bulk. (2) -invariance: , only random mass disorder is allowed. Because , this is the broken time-reversal class D. (3) -invariance: , only random vector potential disorder is allowed. This is the broken time-reversal class AIII. (Technically, it is the “topological”/WZW class in the language of Refs. (7); (49).)

Class AII is generically realized whenever time-reversal is unbroken. Magnetic impurities randomly polarized parallel (perpendicular) to the TI surface manifest as point exchange sources in the vector (mass) potentials of Eq. (10); we are thus tempted to identify symmetry classes D and AIII with these two limits. However, a magnetic impurity will typically induce a local potential fluctuation as well. As a consequence, the generic case of broken-time reversal symmetry corresponds to the absence of , and , which gives the unitary class A.(48); (49); (7) In fact, for a vanishing average mass , the surface of a topological insulator with generic time-reversal breaking disorder is expected to flow under renormalization to the critical point of the integer quantum Hall plateau transition.(24); (29) On a different note, the class and class D versions of in Eq. (10) can be realized on the surface of a bulk -invariant 3D topological superconductor, where time-reversal is respectively preserved or broken at the surface.(7)

## Appendix B Perturbation theory

### b.1 Chiral Decomposition and one-loop renormalization

We write a 2+0-D fermion path integral to represent correlation functions in the disordered Dirac Hamiltonian transcribed in Eq. (10). The fermion operators are replaced with the Grassmann fields ; here denotes a replica index, and we are to send at the end of the calculation.(20); (8) We employ a “chiral decomposition” of the two-component spinors,

 ψi=[LiRi],¯ψi=[¯Ri¯Li] (44)

Then the action of the replicated theory is

 S=∫d2r[ε(¯RiLi+¯LiRi)+¯RiLiϕ+¯LiRi¯ϕ+¯Ri(−i∂+A)Ri+¯Li(−i¯∂+¯A)Li], (45)

where we have introduced complex coordinates , , and disorder potentials , . The energy is a fixed parameter. In Eq. (45), repeated replica indices are summed. Assuming the Gaussian white-noise variances for the disorder potentials enumerated in Eq. (12), the replicated theory can be averaged over disorder configurations. The post-ensemble averaged action is

 ¯¯¯¯S=S0+¯¯¯¯SA+¯¯¯¯S1+¯¯¯¯S2, (46)

where is the clean Dirac action, and

 ¯¯¯¯SA= −2ΔA∫d2r¯RiRi¯LjLj, (47a) ¯¯¯¯S1= −ΔV+ΔM2∫d2r(¯RiLi¯RjLj+¯LiRi¯LjRj), (47b) ¯¯¯¯S2= −(ΔV−ΔM)∫d2r¯RiLi¯LjRj. (47c)

Different replicas become coupled through the disorder.(20); (8)

The disorder-averaged composite LDOS corresponds to the fermion bilinear expectation

 ¯¯¯ν=⟨¯ψψ⟩=⟨¯RL+¯LR⟩. (48)

The spin LDOS was defined by Eq. (7). For the out-of-plane and in-plane (chiral) components, one has

 ¯¯¯¯¯ν^3= ⟨¯ψ^σ3ψ⟩=⟨¯RL−¯LR⟩, (49a) ¯¯¯¯¯¯ν±≡ ⟨¯ψ^σ±ψ⟩=2{⟨¯RR⟩,⟨¯LL⟩}, (49b)

The overlines appearing in the left-hand sides of Eqs. (48) and (49) denote disorder-averaging, whereas the angle brackets on the right-hand sides represent integration in the fermion path integral, using the action in Eq. (46).

A generic local operator corresponding to the moment of some fermion bilinear can be viewed as sum of “strings,” where each string consists of total right (R) and left (L) mover labels, arranged in some order. For example, the disorder-averaged moment of the LDOS is represented by the the composite operator expectation value

 ¯¯¯¯¯¯¯¯¯¯¯¯νq(r)=⟨q∏i=1[¯RiLi(r)+¯LiRi(r)]⟩. (50)

In this equation, a product is taken over operators carrying indices in the first replicas. The LDOS moment is computed by placing one copy of the LDOS operator into each of different replicas; before averaging, this gives in a fixed realization of the disorder. (Placing instead the copies into the same replica would give the disorder-averaged first moment of a -point Green’s function.) The operator in Eq. (50) is an even weight sum of “strings”, all of length . I.e.,

 ¯¯¯¯¯¯¯¯¯¯¯¯νq(r)= {¯RL;¯RL;…;¯RL} +{¯LR;¯RL;¯RL;…;¯RL} +{¯RL;¯LR;¯RL;…;¯RL} +… +{