Criticality Across the Energy Spectrum from Random, Artificial Gravitational Lensing in Two-Dimensional Dirac Superconductors

TL;DR

This paper investigates how random gravitational-like disorder affects two-dimensional Dirac materials, revealing a universal quantum-critical state across energy levels, with implications for high-temperature superconductors and topological states.

Contribution

It introduces a novel model of quenched gravitational disorder in 2D Dirac systems, showing universal quantum-critical behavior and explaining experimental phenomena in cuprate superconductors.

Findings

Disorder induces a quantum-critical 'stack' of states with universal fluctuations.

Model with pseudospin rotations and nematic fluctuations best matches cuprate phenomenology.

Energy width of critical states increases with disorder strength.

Abstract

We numerically study weak, random, spatial velocity modulation ["quenched gravitational disorder" (QGD)] in two-dimensional massless Dirac materials. QGD couples to the spatial components of the stress tensor; the gauge-invariant disorder strength is encoded in the quenched curvature. Although expected to produce negligible effects, wave interference due to QGD transforms all but the lowest-energy states into a quantum-critical "stack" with universal, energy-independent spatial fluctuations. We study five variants of velocity disorder, incorporating three different local deformations of the Dirac cone: flattening or steepening of the cone, pseudospin rotations, and nematic deformation (squishing) of the cone. QGD should arise for nodal excitations in the $d$-wave cuprate superconductors (SCs), due to gap inhomogeneity. Our results may explain the division between low-energy "coherent" (plane-wave-like) and finite-energy "incoherent" (spatially inhomogeneous) excitations in the SC and pseudogap regimes. The model variant that best matches the cuprate phenomenology possesses quenched random pseudospin rotations and nematic fluctuations. This model variant and another with pure nematic randomness exhibit a robust energy swath of stacked critical states, the width of which increases with increasing disorder strength. By contrast, quenched fluctuations that isotropically flatten or steepen the Dirac cone tend to produce strong disorder effects, with more rarified wave functions at low- and high-energies. Our models also describe the surface states of class DIII topological SCs.

Figures

Tables

Table 1. Table 1: Key properties of the effective 2D dirty Dirac Bogoliubov-de Gennes Hamiltonian for quasiparticles in a d 𝑑 d -wave superconductor. The number of nodes coupled refers to elastic impurity scattering between and/or within nodes. The generic model is expected to Anderson localize at all energies. The models restricted to disorder that couples only 2 or 1 node(s) are all effectively topological, i.e. describe the surface states of 3D bulk TSCs. All of the restricted models have the low-energy thermal conductivity given by Eq. ( 13 ), independent of disorder Ludwig1994 ; Tsvelik95 ; Ostrovsky06 and interactions Xie15 . The exponent x 1 / z subscript 𝑥 1 𝑧 x_{1}/z governs the low-energy scaling of the density of states [Eq. ( 16 )]. The restricted models all exhibit “multifractal stacking” of quantum critical wave functions throughout the energy spectrum. The parameter θ ​ ( ε ) 𝜃 𝜀 \theta(\varepsilon) characterizes the multifractal spectrum of wave functions according to Eq. ( 18 ), near energy ε 𝜀 \varepsilon . The stacked critical states for classes CI and AIII apparently belong to the spin and integer quantum Hall plateau transitions (SQHPT and IQHPT, respectively). In the last entry, TQHPT refers to the thermal quantum Hall plateau transition; see Appendix C for a discussion.

# of				Effective
nodes		Class		winding			$x_1/z$		$\theta \left(\epsilon =0\right)$		$\theta \left(\epsilon \ne 0\right)$		Dirt type(s)
coupled				# $\left|\nu \right|$
4		CI		N/A			1 SenthilFisher99		N/A		N/A		Vector, potential,
									(localized) Senthil98 ; AltlandSimonsZirnbauer02		(localized) Evers2008		mass AltlandSimonsZirnbauer02
2		CI		2			1/7 Nersesyan1994		1/4 Mudry1996 ; Caux1996		$\simeq 1/8$ (SQHPT) Ghorashi18		SU(2) vector
1		AIII		1			$\frac{1-{\lambda }_A}{1+{\lambda }_A}$ Ludwig1994		${\lambda }_A$ Ludwig1994		$\simeq 1/4$ (IQHPT) Chou2014 ; Sbierski		U(1) vector
1		DIII		1			1 (clean)		0 (clean)		$\simeq 1/13$ (TQHPT?) Sec. III		Velocity/QGD

Equations

$$H=-\frac{1}{2}\sum _{a,b=1,2}\int d^2\mathrm{r}{v}_{ab}(\mathrm{r})\left(\bar{\psi }i{\hat{\sigma }}^a\overleftrightarrow{{\partial }_b}\psi \right),$$

$$\begin{array}{c}H=\underset{\mathrm{k}}{\int }{\Psi }^{†}(\mathrm{k})\hat{h}(\mathrm{k})\Psi (\mathrm{k}),{\Psi }^{\mathsf{T}}(\mathrm{k})\equiv \left[\begin{array}{c}{c}_{\uparrow }(\mathrm{k})c_{\downarrow }^{†}(-\mathrm{k})\end{array}\right],\\ \hat{h}(\mathrm{k})=\left[\begin{array}{cc}\tilde{\epsilon }(\mathrm{k})& \Delta (\mathrm{k})\\ {\Delta }^{\ast }(\mathrm{k})& -\tilde{\epsilon }(-\mathrm{k})\end{array}\right].\end{array}$$

$$T:$$

$$P:$$

$$T:$$

$$P:$$

$$H=\underset{\mathrm{k}}{\int }{\psi }^{†}(\mathrm{k}){\hat{h}}_0(\mathrm{k})\psi (\mathrm{k}).$$

$$T:$$

$$P:$$

$${\hat{h}}_0=$$

$$+{\hat{P}}_{-}\left[v_f{\hat{\alpha }}_1(-i{\partial }_y)+v_{\Delta }{\hat{\alpha }}_2(-i{\partial }_x)\right],$$

$${\hat{h}}_1=$$

$$+{\hat{P}}_{-}\left[{\hat{\alpha }}_1{A}_{y,i}^{(-)}(\mathrm{r})+{\hat{\alpha }}_2{A}_{x,i}^{(-)}(\mathrm{r})\right]{\hat{t}}^i,$$

$${\hat{h}}_2=$$

$$\hat{h}\equiv {\hat{h}}_0+{\hat{h}}_1+{\hat{h}}_2$$

$$\frac{\kappa }{T}=\frac{k_B^2}{3\hslash }\left(\frac{v_F}{v_{\Delta }}+\frac{v_{\Delta }}{v_F}\right).$$

$${\hat{h}}_{{\mathrm{CI}}_2}=\hat{\mathbf{\sigma }}\cdot \left[-i\nabla +{\mathrm{A}}_i(\mathrm{r}){\hat{\tau }}^i\right],$$

$$T:$$

$$P:$$

$$\underset{\epsilon \to 0}{\lim }\nu (\epsilon )\sim |\epsilon {|}^{x_1/z},$$

$$\sum _{i=1}^{N^2}({\mu }_i{)}^q\sim {\left(\frac{b}{L}\right)}^{\tau (q)}.$$

$$\begin{array}{rl}\tau (q)=& 2(q-1)+\Delta (q),\\ \Delta (q)=& \theta q(1-q),0\le |q|\le q_c,q_c\equiv \sqrt{\frac{2}{\theta }},\end{array}$$

$${\hat{h}}_{{\mathrm{AIII}}_1}=\hat{\mathbf{\sigma }}\cdot \left[-i\nabla +\mathrm{A}(\mathrm{r})\right].$$

$$x_1/z=(1-{\lambda }_A)/(1+{\lambda }_A).$$

$$\overline{A_{\bar{a}}(\mathrm{r})A_{{\bar{a}}^{\prime }}({\mathrm{r}}^{\prime })}=\pi {\lambda }_A{\delta }_{\bar{a},{\bar{a}}^{\prime }}{\delta }^{(2)}(\mathrm{r}-{\mathrm{r}}^{\prime }).$$

$$\theta ={\lambda }_A.$$

$$S=\int dt{d}^2\mathrm{r}\left[\bar{\psi }i{\partial }_t\psi +\sum _{a,b=1,2}\frac{v_{ab}(\mathrm{r})}{2}\left(\bar{\psi }i{\hat{\sigma }}^a\overleftrightarrow{{\partial }_b}\psi \right)\right].$$

$$S=$$

$$S=\int d^2\mathrm{r}\left\{-\left[{\eta }^{ij}+\delta G^{ij}(\mathrm{r})\right]T_{ij}+\epsilon \left(\bar{R}L+\bar{L}R\right)\right\},$$

$$T_{ij}=$$

$$=$$

$${\eta }^{ij}\to 2\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right].$$

$$\begin{array}{rl}\delta G^{+-}=& \left(\delta v_{11}+\delta v_{22}+i{v}_{12}-i{v}_{21}\right),\\ \delta G^{-+}=& \left(\delta v_{11}+\delta v_{22}-i{v}_{12}+i{v}_{21}\right),\\ \delta G^{++}=& \left(-\delta v_{11}+\delta v_{22}-i{v}_{12}-i{v}_{21}\right),\\ \delta G^{--}=& \left(-\delta v_{11}+\delta v_{22}+i{v}_{12}+i{v}_{21}\right).\end{array}$$

$$\begin{array}{c}{v}_{11}(\mathrm{r})=v_{22}(\mathrm{r})=v_0\left\{1+\vartheta \left[e{A}^0(\mathrm{r})/E_{\mathsf{bulk}}\right]\right\},\\ v_{12}=v_{21}=0,\end{array}$$

Full text

I Introduction

Understanding the interplay between strong correlations and quenched disorder in low-dimensional superconductors remains one of the key challenges in condensed matter physics TrivediReview . While pair-breaking due to elastic impurity scattering is detrimental to superconductivity, spatial inhomogeneity can locally enhance pairing. Near a superconductor-insulator transition, the latter can occur via the spatial accumulation of Cooper pairs, due to the multifractal rarification of single-particle eigenstates Feigelman07 ; Feigelman10 ; Foster12 ; Foster14 ; Burmistrov15 . Multifractality arises in quantum-critical states due to wave interference induced by multiple impurity scattering. Two-dimensional (2D) multifractal superconductivity has very recently come to the fore in studies of transition metal dichalcoginides Ji19 ; Ugeda18 ; Yeom19 .

An enduring mystery is the role of spatial inhomogeneity in the high-$T_{c}$ cuprate superconductors. These materials are characterized by two different energy scales in the underdoped regime: the pseudogap energy $\Delta_{1}$ and the smaller pairing energy $\Delta_{0}$; the former (latter) increases (decreases) with increased underdoping Sawatzky08 . Maps of $\Delta_{1}(\bm{\mathrm{r}})$ obtained via STM demonstrate strong nanoscale spatial inhomogeneity that increases with underdoping Davis01 ; Davis02 ; Davis05-a ; Davis05-b ; Davis08 ; DavisReview . A key observation in STM scans of the local density of states (LDOS) is that the lower energy scale $\Delta_{0}$ appears to divide fermionic excitations into two distinct classes. States with energies smaller than $\Delta_{0}$ behave like dispersive, “coherent” Bogoliubov-de Gennes quasiparticles, showing robust quasiparticle interference. States with energies above $\Delta_{0}$ are instead termed “incoherent,” as they exhibit strong spatial fluctuations that vary little with energy Davis08 ; DavisReview . Renewed urgency for understanding spatial inhomogeneity comes from Ref. Welp18 , which reported an increase in $T_{c}$ with increasing disorder.

Much recent work on the strange physics of the cuprates invokes gravitational descriptions of quantum criticality Sachdev19 . In this paper, we uncover a very different role possibly played by quantum criticality in these materials, due to a different type of “gravitational” physics. We show that the dichotomy between plane-wave-like, low-energy quasiparticle states and strongly inhomogeneous, finite-energy states can be reconciled in a simple model of noninteracting nodal Dirac quasiparticles, subject to random velocity disorder. Formally this can be cast as the effect of static spacetime curvature Ramond ; Volovik ; Juan12 [“quenched gravitational disorder” (QGD)]. For massless Dirac fermions, propagation is analogous to the lensing of starlight through a fixed, randomly gravitating spacetime Carroll . The gravitational language precisely defines the gauge-invariant content of the disorder, through the induced curvature.

QGD is realized whenever Dirac carriers arise from a correlation gap that is spatially inhomogeneous. In a $d$-wave superconductor, spatial gap fluctuations modulate the nodal quasiparticle velocity, see Fig. 1. In the context of conventional 2D Dirac materials such as graphene or the surface states of 3D topological insulators, QGD has so far received little attention Juan12 ; Nakai14 ; DiazFernandez17 ; Jafari19 . This is because disorder can usually couple in a more relevant fashion, through gauge and mass potentials Juan12 ; AltlandSimonsZirnbauer02 ; AleinerEfetov06 ; Evers2008 . These perturbations typically induce metallic or Anderson insulating behavior on the largest scales Bardarson07 ; Nomura07 ; Nomura08 ; Evers2008 .

In this paper, we use exact diagonalization to probe the effect of QGD on 2D Dirac carriers. While the low-energy states near the Dirac point may be only weakly affected, we find that most of the energy spectrum converges into a “stack” of critical wave states. These quantum-critical wave functions apparently exhibit universal multifractal LDOS fluctuations. A representative state is shown in Fig. 1(b), while LDOS maps of states with different energies appear in Figs. 2 and 3.

The Dirac model that we study is defined by the Hamiltonian

[TABLE]

where $\bm{\mathrm{r}}=\{x_{1},x_{2}\}$, $\psi=\psi_{\sigma}$ is a two-component spinor ($\sigma\in\{\uparrow,\downarrow\}$), $\{\hat{\sigma}^{1,2}\}$ are Pauli matrices in the usual basis, $A\!\!\stackrel{{\scriptstyle\leftrightarrow}}{{\vspace{-12pt}\partial}}\!\!B=A\partial B-(\partial A)B$, and the QGD is encoded in the four velocity components \big{\{}v_{11}(\bm{\mathrm{r}})\equiv 1+\delta v_{11}(\bm{\mathrm{r}}), $v_{22}(\bm{\mathrm{r}})\equiv 1+\delta v_{22}(\bm{\mathrm{r}}),$ $v_{12}(\bm{\mathrm{r}}),$ v_{21}(\bm{\mathrm{r}})\big{\}}. The four potentials couple to the spatial components of the energy-momentum tensor $T^{ab}$. We consider five different variants of the model, with different combinations of the random velocity perturbations. The variants are defined as follows,

1. (a)

Independent $\{\delta v_{11},\delta v_{22}\}$, $v_{12}=v_{21}=0$. Local isotropic flattening or steepening of the Dirac cone and nematic (squishing) of the cone. 2. (b)

Independent $\{v_{12},v_{21}\}$, $\delta v_{11}=\delta v_{22}=0$. Local pseudospin rotations and nematic (squishing) of the Dirac cone. 3. (c)

Independent $\{\delta v_{11}=\delta v_{22},\;v_{12}=-v_{21}\}$. Local isotropic flattening or steepening of the Dirac cone and pseudospin rotations. 4. (d)

Independent $\{\delta v_{11}=-\delta v_{22},\;v_{12}=v_{21}\}$. Local nematic (squishing) of the Dirac cone. 5. (e)

Independent $\{\delta v_{11},\delta v_{22},v_{12},v_{21}\}$. The generic model.

All nonzero random velocity components are taken to be short-range correlated, characterized by a dimensionless disorder strength $\lambda$. The physical disorder variance is expressed as $\lambda$ times the square of the inverse momentum cutoff (the disorder is formally irrelevant at the Dirac point—see below).

We find that critical-state stacking occurs for sufficiently weak disorder in all models (a)–(e). Stacking improves in models (b) and (d) with increasing disorder, i.e. more states become quantum-critical with identical, universal spectra throughout more of the spectrum. This is interesting because model (d) includes only pure nematic deformations of the Dirac cone. Nematicity has emerged as a key possible ingredient in the pseudogap phase of the cuprate superconductors NematicReview ; PDWReview , including quenched nematic fluctuations Davis2019 . By contrast, models (a) and (c) (which incorporate local flattening and steepening of the cone) are most susceptible to strong disorder effects, showing supercritical (more rarified) states at low- and high-energies for increasing randomness. The generic model (e) shows intermediate behavior, including a robust energy swath with critical stacking, but also supercritical behavior at low energies for stronger disorder. Critical stacking is exhibited for models (b) and (e) via the LDOS plots shown in Figs. 2 and 3.

Our results are very surprising, in that QGD is a strongly irrelevant perturbation at low energies. QGD was instead expected to preserve the quasi-ballistic nature of the clean system, even away from the Dirac point Nakai14 . Moreover, critical wave states typically arise only at isolated energies, such as a mobility edge or at the plateau transition of a quantum Hall effect Huckestein1995 ; EversSQHE03 ; MirlinSQHE03 ; Evers2008 . The population of critical states induced by QGD is compared to the total density of states in Figs. 5–9 for models (a)–(e). Our calculations are performed in momentum space; the largest size studied is a 193 $\times$ 193 grid. Since Dirac cones typically span $\sim 10\,\%$ of the 2D Brillouin zone, this corresponds roughly to a (200 nm)2 map. For models (b) and (d), more of the spectrum becomes more critical with increasing disorder (Figs. 6 and 8), with very little dependence on the system size (Figs. 14 and 15). The Dirac model with QGD in Eq. (1) resides in class DIII; despite this, we find no evidence for weak antilocalization at finite energy in any of the variants (a)–(e), contrary to expectations based upon standard symmetry arguments Evers2008 ; Ghorashi18 ; Sbierski .

Model (b) exhibits a phenomenology most similar to the LDOS maps observed in STM studies of BSCCO Davis08 ; DavisReview , with low-energy plane wave states and a robust quantum-critical stack at intermediate and higher energies, see Fig. 2. In addition to quenched nematic fluctuations, model (b) incorporates local rotations of the pseudospin $\{\hat{\sigma}^{1},\hat{\sigma}^{2}\}$ relative to the coordinate axes $\{x_{1},x_{2}\}$. As we review in Sec. II, other types of disorder afflicting $d$-wave quasiparticles are predicted to induce different physics: (I) One class of (effectively topological Foster14 ) perturbations should also produce critical scaling throughout the quasiparticle energy spectrum Ostrovsky2007 ; Chou2014 ; Ghorashi18 ; Sbierski . In this scenario, though, zero-energy states are also predicted to be multifractal, and the low-energy DOS $\nu(\varepsilon)$ would be enhanced via a sublinear power law, $\nu(\varepsilon\rightarrow 0)\sim|\varepsilon|^{\delta}$ with $\delta<1$ Ludwig1994 ; Nersesyan1994 ; Mudry1996 ; Caux1996 ; Bhaseen2001 ; AltlandSimonsZirnbauer02 . (II) Disorder that induces generic internode scattering is instead predicted to Anderson localize the entire quasiparticle spectrum Senthil98 ; AltlandSimonsZirnbauer02 . Neither scenario (I,II) is consistent with STM data DavisReview .

QGD should also arise at the surface of a topological superconductor (TSC). Due to “topological protection,” velocity modulation is the only allowed coupling of charged impurities to the 2D massless Majorana fluid expected to form at the surface of a class DIII TSC TSCRev1 ; TSCRev2 ; TSCRev3 ; TSCRev4 with winding number $|\nu|=1$ Nakai14 ; Roy19 , as we show in Appendix B.

The stack of critical finite-energy states found here is quite unusual, but not unprecedented. In Ref. Ghorashi18 , we observed stacks of critical, class C spin quantum Hall plateau transition (QHPT) states at finite energy in a surface model for a class CI TSC. Numerical studies Chou2014 ; Sbierski have found that the finite-energy surface states of class AIII TSCs sit at the class A integer QHPT Ludwig1994 ; Ostrovsky2007 . We expect that the finite-energy critical states identified here correspond to a version of the class D thermal QHPT D-1 ; D-2 ; D-3 ; D-4 ; D-5 ; D-6 ; D-7 .

I.1 Outline

This paper is organized as follows. In Sec. II, we review results on impurity scattering in $d$-wave superconductors. We derive the low-energy 4-node Dirac theory from a microscopic model. We review how restrictions to 2- and 1-node models are effectively topological, i.e. correspond to models for surface states of bulk TSCs. We summarize results for the thermal conductivity, density of states, and multifractal spectra of low- and finite-energy wave functions. We discuss the quantum critical “stacking” of multifractal finite-energy states and the connection to plateau transitions in quantum Hall effects that has recently been found to occur in these models Ghorashi18 ; Sbierski .

In Sec. III we present the main numerical results of this paper for the five different variants (a)–(e) of the model defined by Eq. (1). Our results consist of the density of states and multifractal exponents throughout the energy spectrum, for a range of disorder strengths and system sizes. We contrast our results to conventional expectations, and discuss them in light of the multifractal stacking phenomenon.

The appendices detail more technical considerations. In Appendix A, we show how velocity randomness is embedded in a generally covariant framework for 2+1-D Dirac fermions propagating through curved spacetime. In Appendix B we show how electric potentials couple gravitationally to the single Majorana surface cone predicted to occur at the boundary of a class DIII TSC. Finally, in Appendix C we present a symmetry analysis of the Dirac fermion with QGD, and the possible connection to the thermal QHPT in class D.

II Disorder in 2D $d$-wave superconductors, and topological superconductor surface fluids

In this section we review a microscopic model for elastic impurity scattering in 2D $d$-wave superconductors. The generic model with elastic scattering between all four low-energy Dirac quasiparticle nodes Anderson localizes at all energies. We then show how restrictions on internode scattering produce independent, topologically protected sectors. The latter are equivalent to the surface states of TSCs in classes CI, AIII, and DIII. Finally, we review exact results for the thermal conductivity, density of states, and multifractal spectrum of local density of states fluctuations. With the exception of the recent spectrum-wide criticality numerical results obtained in Refs. Ghorashi18 ; Sbierski and in this paper, all of the following is well known AltlandSimonsZirnbauer02 ; SRFL08 ; Foster14 .

The discussion below is based on the linearized quasiparticle band structure in the vicinity of the Dirac nodes. For a recent numerical study of disorder effects in a lattice model for the cuprates, see e.g. Ref. Zaanen18 .

Readers primarily interested in new results for the Dirac model with QGD can skip this section, and proceed directly to Sec. III.

II.1 $d$-wave model

A 2D spin-singlet superconductor has the Bogoliubov-de Gennes Hamiltonian

[TABLE]

Here $\Psi(\bm{\mathrm{k}})$ is the Nambu spinor that annihilates spin-1/2 quasiparticles, $\mathsf{T}$ denotes the matrix transpose, $\tilde{\varepsilon}(\bm{\mathrm{k}})=\varepsilon(\bm{\mathrm{k}})-\mu$ is the bare energy dispersion relative to the chemical potential $\mu$, $\Delta(\bm{\mathrm{k}})$ is the mean-field BCS order parameter, and $\int\limits_{\bm{\mathrm{k}}}\equiv\int d^{2}\bm{\mathrm{k}}/(2\pi)^{2}$.

The two key symmetries we want to enforce are $T^{2}=-1$ time-reversal and spin SU(2) symmetry, the combination of which places the system in the Altland-Zirnbauer class CI AltlandSimonsZirnbauer02 ; Evers2008 . U(1) symmetry under $\hat{S}^{z}$-axis spin rotations is manifest in Eq. (4). A rotation around $\hat{S}^{x}$ or $\hat{S}^{y}$ generates a continuous particle-hole transformation in the Nambu language. To ensure SU(2) invariance, it is sufficient to impose invariance under a $\pi$ $\hat{S}^{x}$-rotation, which takes the form of a $P^{2}=-1$ effective particle-hole symmetry. These are encoded via the transformations

[TABLE]

where the Pauli matrices $\{\hat{\sigma}^{1,2,3}\}$ act on the particle-hole space grading the Nambu spinor. The $T$ ($P$) transformation is antiunitary (unitary) in second quantization. These translate into the following conditions on $\hat{h}(\bm{\mathrm{k}})$,

[TABLE]

Physical time-reversal invariance is thus represented by an effective unitary chiral symmetry in first quantization. Imposing both $T$ and spin SU(2) implies that $\tilde{\varepsilon}(-\bm{\mathrm{k}})=\tilde{\varepsilon}(\bm{\mathrm{k}})$, and that $\Delta(\bm{\mathrm{k}})=\Delta^{*}(\bm{\mathrm{k}})=\Delta(-\bm{\mathrm{k}})$.

For a $d_{x^{2}-y^{2}}$ superconductor, we take $\Delta(\bm{\mathrm{k}})=\Delta_{0}\,k_{x}k_{y}/k_{F}^{2}$, where $\Delta_{0}$ is a real amplitude, $k_{F}$ is the Fermi wave vector, and $(k_{x},k_{y})$ measure momentum along the nodal directions, see Fig. 4. Linearizing the Hamiltonian in the vicinity of each node gives an effective low-energy, anisotropic Dirac theory,

[TABLE]

Here $\psi(\bm{\mathrm{k}})$ is an 8-component spinor with long wavelength momentum $\bm{\mathrm{k}}$. The field $\psi$ is formed from the direct sum of $\Psi(\bm{\mathrm{k}}_{i}+\bm{\mathrm{k}})$ evaluated in the vicinity of the four nodes $i\in\{1,2,3,4\}$ depicted in Fig. 4, located at $\{\bm{\mathrm{k}}_{i}\}$. Let the Pauli matrices $\{\hat{\tau}^{1,2,3}\}$ act on the partners of a pair of nodes related by time-reversal, i.e. the matrix $\hat{\tau}^{1}$ exchanges states between nodes $(1\Leftrightarrow 2)$ and $(3\Leftrightarrow 4)$ in Fig. 4. Let the Pauli matrices $\{\hat{\kappa}^{1,2,3}\}$ act on pairs of nodes, i.e. the matrix $\hat{\kappa}^{1}$ exchanges $(1\Leftrightarrow 3)$ and $(2\Leftrightarrow 4)$. The symmetries in Eq. (6) become

[TABLE]

The Bogoliubov-de Gennes Hamiltonian takes the Dirac form

[TABLE]

where $\hat{P}_{\pm}\equiv(1/2)(\hat{1}\pm\hat{\kappa}^{3})$ projects onto nodal pair ($1,2$) ($+$) or ($3,4$) ($-$), and where the Clifford algebra matrices are $(\hat{\alpha}_{1},\hat{\alpha}_{2})=(\hat{\sigma}^{3}\hat{\tau}^{3},\hat{\sigma}^{1}\hat{\tau}^{3}).$ In Eq. (II.1), $v_{F}$ is the bare Fermi velocity, while the perpendicular dispersion $v_{\Delta}\simeq\Delta_{0}/k_{F}$ arises from the pairing.

Non-magnetic impurities can couple to the low-energy theory via local, Hermitian fermion bilinear operators. The most relevant of these in the RG sense have no derivatives, and take the form $\mathcal{O}^{{\scriptscriptstyle{({\alpha\beta\gamma})}}}(\bm{\mathrm{r}})\equiv\psi^{\dagger}(\bm{\mathrm{r}})\,\hat{\sigma}^{\alpha}\,\hat{\tau}^{\beta}\,\hat{\kappa}^{\gamma}\,\psi(\bm{\mathrm{r}}),$ where $\alpha,\beta,\gamma\in\{0,1,2,3\}$ and $\hat{\sigma}^{0}=\hat{\tau}^{0}=\hat{\kappa}^{0}=\hat{1}$ (the identity matrix). Imposing time-reversal and spin SU(2) symmetries from Eq. (8) reduces the number of allowed Hermitian bilinears to 20. This is consistent with the random matrix classification for class CI, where an $8\times 8$ matrix is formed from the generators of Sp$(8)/$U$(4)$. We can parameterize the 20 allowed perturbations via intrapair and interpair scattering Hamiltonians $\hat{h}_{1,2}$. Pure intrapair [$(1\Leftrightarrow 2)$, $(3\Leftrightarrow 4)$] scattering is encoded in

[TABLE]

where the repeated index $i$ is summed over $\{1,2,3\}$. The intrapair scattering terms take the form of SU(2) gauge potentials $\bm{\mathrm{A}}_{i}^{{\scriptscriptstyle{({\pm})}}}(\bm{\mathrm{r}})\,\hat{t}^{i}$, where the SU(2) generators are $\hat{t}^{i}=\{\hat{\sigma}^{2}\hat{\tau}^{1},\hat{\sigma}^{2}\hat{\tau}^{2},\hat{\tau}^{3}\}.$ The remaining 8 allowed perturbations scattering between pairs of nodes, and can be expressed through the Hamiltonian

[TABLE]

where repeated indices $i\in\{1,2,3\}$ and $\bar{a}\in\{1,2\}$ are summed. In contrast to the intrapair scattering encoded in $\hat{h}_{1}$, the interpair scattering mediated by $B_{\bar{a},i}(\bm{\mathrm{r}})$ and $C_{\bar{a}}(\bm{\mathrm{r}})$ take the form of scalar potentials and mass terms. The latter can be used to open up a gap without breaking $T$; this immediately implies that the generic $4$-node Hamiltonian

[TABLE]

cannot describe the isolated surface states of a strong TSC or insulator.

It is possible to derive the strengths and correlations of all 20 disorder potentials $\{A_{\bar{a},i}^{{\scriptscriptstyle{({\pm})}}},B_{\bar{a},i},C_{\bar{a}}\}$ in the effective low-energy field theory from microscopic perturbations of the original model in Eq. (4) AltlandSimonsZirnbauer02 , but we will not pursue this here. In fact, our main interest will be not in these potential perturbations, but in the weaker “quenched gravitational disorder” (QGD, nodal velocity randomness), which we argue below should be present in the cuprate superconductors DavisReview . In Appendix B, we provide a derivation of QGD from the coupling of the electric potential to Majorana surface states in a microscopic model of a class DIII TSC.

II.2 Topological restrictions and TSC surface states: Localization, quantum criticality, and physical properties

II.2.1 Four coupled nodes: Spectrum-wide Anderson localization

Quantum wave interference induced by generic elastic scattering between all four nodes in the model described by Eqs. (II.1)–(12) is expected to Anderson localize all quasiparticle states, at both zero and finite single-particle energies.

The localization near zero energy is understood as follows. The Dirac model described by Eq. (12) can be averaged over all 20 disorder potentials using the replica or supersymmetry (SUSY) trick. For Gaussian white-noise distributions, the model flows to strong coupling under the perturbative renormalization group AltlandSimonsZirnbauer02 . One then expects that the system can be described by a replicated or SUSY nonlinear sigma model in class CI. Although this sigma model admits a Wess-Zumino-Novikov-Witten (WZNW) term, it can be shown that no such term arises in the full four-node theory (e.g., it is forbidden by average parity invariance AltlandSimonsZirnbauer02 ; moreover, the WZNW model turns out to be topological SRFL08 , as reviewed below). Without the WZNW term, the 2D class CI sigma model also flows to strong coupling, interpreted as the tendency towards Anderson localization Senthil98 .

Anderson localization is anticipated at finite energy as well in the four-node model. The standard argument (which, however, needs to be revised for topological models as we emphasize below) goes as follows. At zero energy, 7 of the 10 Altland-Zirnbauer classes are characterized by a special particle-hole or chiral symmetry. This is true for all classes that describe spin-1/2, time-reversal-invariant superconductors in classes CI, AIII, and DIII SRFL08 ; Foster14 . This symmetry also relates states at positive and negative energies, but does not tell us anything about the character of the states at some particular fixed energy $\varepsilon\neq 0$. On the other hand, for quasiparticle states in a superconductor, one can always consider $|\varepsilon|\gg\Delta_{0}$, where $\Delta_{0}$ is the pairing energy. Then the quasiparticle states should reside in the standard Wigner-Dyson class of the normal metal that hosts the superconductivity. Since $\varepsilon=0$ is the only symmetry-distinguished energy, one concludes that all states with $\varepsilon\neq 0$ in the infinite-volume limit must reside in the orthogonal, unitary, or symplectic Wigner-Dyson class. For the four-node class CI model, the appropriate Wigner-Dyson class is the orthogonal metal class AI, which also preserves $T$ and spin SU(2) symmetry. This class is believed to always localize in 2D Evers2008 .

Localization of all quasiparticle states is not borne out by experiments in the high-$T_{c}$ cuprate superconductors, at or below optimal hole doping. First of all, localization is at variance with superconductivity itself, since $T_{c}$ is not protected by Anderson’s theorem for non-$s$-wave pairing AltlandSimonsZirnbauer02 . Second, STM data taken at energies less than a characteristic scale $\Delta_{0}$ (which is always below the pseudogap energy $\Delta_{1}$ on the underdoped side) show robust quasiparticle interference, a sign that there is just enough internode scattering to reveal the clean dispersion, but not enough to induce localization Davis08 ; DavisReview .

Finally, the experimentally measured low-temperature longitudinal thermal conductivity Taillefer00 ; Orenstein00 is nonzero, and close to the universal (disorder-independent) theoretical result

[TABLE]

This result was originally obtained via an approximate, self-consistent semiclassical calculation Lee93 ; Senthil98 ; LeeDurst00 . Eq. (13) is better understood as the exact result for the $T\rightarrow 0$ limit of the Landauer thermal conductivity in the clean four-node model, the analog (via Wiedemann-Franz) of the “ballistic” conductivity $\sigma^{xx}=4e^{2}/\pi h$ for pristine graphene doped exactly to charge neutrality. The absence of impurity scattering combined with the vanishing density of states produces a finite, universal result due to evanescent wave propagation Tworzydlo06 . A very special feature of the topological models reviewed below is that Eq. (13) remains exact in the presence of disorder Ludwig1994 ; Tsvelik95 ; Ostrovsky06 ; SRFL08 , and is even predicted to be independent of virtual interaction (Altshuler-Aronov) corrections in these special models Xie15 . By contrast, the Anderson localized model in Eq. (12) would have $\kappa/T\rightarrow 0$ as $T\rightarrow 0$.

II.2.2 Two coupled nodes: class CI TSC surface states and spectrum-wide spin quantum Hall criticality

If we restrict ourselves to intrapair scattering, such that the node pair (1,2) decouples from (3,4) (Fig. 4), then $\hat{h}=\hat{h}_{0}+\hat{h}_{1}$ [Eqs. (II.1) and (II.1)]. Since the pairs are independent, we focus on one pair (1,2). After a rescaling of $(x,y)$ and a basis change, the two-node Hamiltonian can be written as

[TABLE]

where $i\in\{1,2,3\}$ and $\hat{\bm{\sigma}}\equiv\hat{\sigma}^{1}\,\hat{x}+\hat{\sigma}^{2}\,\hat{y}$. Here $\hat{h}_{{\mathrm{CI}_{2}}}$ is a $4\times 4$ matrix differential operator acting in the composite (particle-hole)$\otimes$(valley) ($\sigma\otimes\tau$) space. Since nodes (1,2) are related by $T$, the model still resides in class CI. The symmetry operations in Eq. (8) become

[TABLE]

Given the (transformed) form of the Clifford algebra in Eq. (14), the effective chiral/physical time-reversal symmetry condition in Eq. (15a) is anomalous, and cannot be realized without fine-tuning in two spatial dimensions BernardLeClair02 ; SRFL08 . It is naturally realized on the surface of a class CI TSC SRFL08 ; Schnyder09 ; Foster12 ; Foster14 , with minimal winding number $|\nu|=2$. The sigma model takes the same form as the four-node model, except that it is now augmented with a WZNW term at level $k=1$ Nersesyan1994 ; Mudry1996 ; Caux1996 ; Bhaseen2001 ; Foster14 ; Schnyder09 .

The density of states $\nu(\varepsilon)$ exhibits critical scaling with energy $\varepsilon$ in the limit $|\varepsilon|\rightarrow 0$. In particular,

[TABLE]

where $x_{1}/z=1/7$ Nersesyan1994 . The zero-energy wave function $\psi_{0}(\bm{\mathrm{r}})$ exhibits quantum critical fluctuations on all length scales; these are characterized by the multifractal spectrum $\tau(q)$ (for a review, see e.g. Evers2008 ). If the system has size $L\times L$, one divides this up into $N^{2}$ boxes of size $b$, with $N\equiv L/b$. Then one introduces the box probability $\mu_{i}\equiv\int\limits_{b_{i}}d^{2}\bm{\mathrm{r}}\,|\psi_{0}(\bm{\mathrm{r}})|^{2}$, where the integral is performed over the $i^{\mathrm{th}}$ box. The multifractal spectrum is defined via the scaling behavior of moments of the box probability,

[TABLE]

Box probabilities can also be obtained by normalizing a spatial map of the local density of states in an STM experiment. In the limit $L\rightarrow\infty$, $\tau(q)$ is self-averaging. For the topological class CI model in Eq. (14), the spectrum is perfectly parabolic. The exact result is

[TABLE]

with $\theta=1/4$ Mudry1996 ; Caux1996 .

Very recently, the authors considered the question of the finite-energy states of the model in Eq. (14). On one hand, the same argument presented in Sec. II.2.1 leads to the conclusion that all finite-energy states should be Anderson localized in the orthogonal class AI; this was the “conventional wisdom” Evers2008 . On the other hand, Eq. (14) also describes a surface quasiparticle fluid that forms at the boundary of a bulk TSC SRFL08 ; Schnyder09 , protected by the anomalous form of time-reversal symmetry in Eq. (15a). From this perspective, the idea that only the zero-energy single-particle wave function $\psi_{0}(\bm{\mathrm{r}})$ escapes Anderson localization appears very strange. Indeed, it would correspond to a very weak form of “topological protection,” since in other topological phases such as quantum Hall liquids, 2D and 3D topological insulators, it is the entire band of edge or surface excitations that is protected from Anderson localization KaneHasan ; TSCRev1 .

The only alternative to localization was argued in Ref. Ghorashi18 to be a “stacking” of identical, quantum critical wave functions at all nonzero energies. It was argued that each such wave function should sit at the class C, spin quantum Hall plateau transition (SQHPT) SQHP1 ; SQHP2 ; SQHP3 . The latter shares a few critical exponents with classical 2D percolation SQHP2 ; EversSQHE03 ; MirlinSQHE03 ; Bhardwaj2015 , a logarithmic conformal field theory SQHP4 ; SQHP5 .

The numerical results of Ref. Ghorashi18 are consistent with the “critical stacking” scenario. Near zero energy, the DOS and multifractal spectrum confirm the predictions of the WZNW theory [Eqs. (16) and (18), with $x_{1}/z=1/7$ and $\theta=1/4$]. At intermediate energies, however, a wide swath of states is found to exhibit weaker universal multifractality, given approximately by Eq. (18), but now with $\theta\simeq 1/8$. The latter is consistent with the SQHPT EversSQHE03 ; MirlinSQHE03 . More states exhibit SQHPT phenomenology upon increasing the disorder strength or system size Ghorashi18 .

Finally, we note that the results described above are clearly incompatible with experiments in the cuprates. This is not surprising, because in the 2D $d$-wave model (Fig. 4) it is difficult to microscopically suppress scattering between node pairs $(1,2)\Leftrightarrow(3,4)$, while simultaneously retaining significant internode scattering between partners ($1\Leftrightarrow 2$) and ($3\Leftrightarrow 4$). The low-energy density of states vanishes as a strongly sublinear power-law $\nu(\varepsilon)\sim|\varepsilon|^{1/7}$. The two-node model exhibits the strongest multifractality at zero energy [Eq. (18) with $\theta=1/4$], and weaker multifractality at finite energies [Eq. (18) with $\theta\simeq 1/8$ Ghorashi18 ]. These features are opposite the observations in STM on BSCCO, which show minimal spatial inhomogeneity in low-energy LDOS maps, with stronger inhomogeneity above the energy scale $\Delta_{0}$; in addition, the low-energy DOS retains the linear character of the clean system DavisReview .

II.2.3 One node, vector potential disorder: class AIII TSC surface states and spectrum-wide quantum Hall criticality

We can further restrict to pure intranode scattering. The effective Hamiltonian is “half” of Eq. (14), i.e. a single 2-component Dirac fermion subject to U(1) vector potential randomness,

[TABLE]

Formally this single-node Hamiltonian is the same as the surface fluid of a class AIII TSC with minimal winding number $|\nu|=1$. The anomalous (topological) time-reversal symmetry is still encoded by Eq. (15a).

The density of states scales as in Eq. (16), with exponent Ludwig1994

[TABLE]

Here $\lambda_{A}$ denotes the variance of the assumed white-noise vector disorder potential,

[TABLE]

At zero energy, the multifractal spectrum is given by Eq. (18), with Ludwig1994

[TABLE]

For a review of this model in the context of TSCs and its higher-winding number generalizations, see e.g. Ref. Foster14 .

At finite energy, the eigenstates of Eq. (19) form a stack of quantum-critical wave functions as in the CI case. These were expected to reside at the ordinary class A integer quantum Hall plateau transition (IQHPT) Ludwig1994 ; Ostrovsky2007 , and this result has been confirmed numerically in Refs. Chou2014 ; Sbierski . The IQHPT has an approximately parabolic multifractal spectrum as in Eq. (18), with $\theta\simeq 1/4$ Huckestein1995 ; Evers2008 .

The spectrum-wide “stacked” IQHPT multifractality is quite strong. Unlike the $|\nu|=2$ class CI model, the low-energy class AIII model predictions for the DOS and multifractal spectrum depend on the nonuniversal parameter $\lambda_{A}$. When the latter is strong enough to render the finite-energy states critical over a length scale that is not too large, one would expect to see multifractality extending all the way down to zero energy [Eq. (21)], as well as a nonlinear enhancement of the DOS [Eqs. (16) and (20)]. Neither of these features are seen in STM data on BSCCO DavisReview .

II.2.4 One node, gravitational disorder due to spatial gap inhomogeneity: class DIII TSC surface states and spectrum-wide criticality

The simplest possible model neglects all forms of potential scattering. However, even in this case it is possible for disorder to produce a nonzero effect. In particular, a slow spatial modulation of $d$-wave the gap amplitude $\Delta_{0}=\Delta_{0}(\bm{\mathrm{r}})$ should induce spatial modulation of $v_{\Delta}\sim\Delta_{0}/k_{F}$ in Eq. (II.1), referred to in the sequel as QGD. Because the disorder couples to an operator with a spatial derivative, it is formally irrelevant in an RG sense at zero energy, in contrast with the potential perturbations in Eqs. (II.1) and (11). The latter are marginal (at tree level); these conclusions assume short-range correlated disorder. As a result, for sufficiently weak QGD in Eq. (1), the low-energy DOS $\nu(\varepsilon)\sim|\varepsilon|$ as in the clean limit, and the states near zero energy are not multifractal. Although velocity disorder modifies the definition of the (spin) current operator, the low-temperature thermal conductivity should be unchanged from the clean Landauer result in Eq. (13), due to the irrelevance of disorder near $\varepsilon=0$.

The effects of such velocity disorder at finite $\varepsilon\neq 0$ are the main focus of this paper; the setup and results are discussed in the next section. Formally, a single node with QGD resides in class DIII, and could be realized as a dirty Majorana cone on the surface a bulk TSC (such as the candidate material CuxBi2Se3 TSCRev1 ). This is discussed in detail in Appendix B. A summary of the results discussed in this section appears in Table 1.

III 2D Dirac fermions with quenched gravitational disorder: Results

III.1 Model, formulations, and applications

The single 2D Dirac fermion with quenched velocity disorder described by the Hamiltonian in Eq. (1) can be associated to a (2+1)-D action

[TABLE]

The action can be interpreted in terms of fermions propagating through curved spacetime, with a static metric $g_{\mu\nu}(\bm{\mathrm{r}})$ defined explicitly in terms of the velocity fields $\{v_{ab}\}$ in Appendix A. Eq. (22) obtains from $g_{\mu\nu}(\bm{\mathrm{r}})$ and from the covariant action

[TABLE]

where $\mu\in\{t,x_{1},x_{2}\}$ and $A,B,C\in\{0,1,2\}$; repeated indices are summed. Here $\sqrt{|g|}\,d^{3}x$ is the volume measure, $\{\hat{\gamma}^{A}\}$ are the gamma matrices, $E^{\mu}_{A}$ is the “dreibein,” $\omega_{\mu}{}^{BC}$ is the spin connection, and $\hat{S}_{BC}$ generates local Lorentz transformations Carroll . Since the velocity modulation enters through the effects of spacetime curvature in Eq. (23), we alternatively refer to this as “quenched gravitational disorder” (QGD). The gauge-invariant content of the disorder can be characterized via the induced curvature, as shown in Appendix A.

The Dirac model in Eqs. (1) and (22) is noninteracting. We can therefore alternatively cast the problem in terms of (2+0)-D action designed to compute Green’s functions at a given, fixed single particle energy $\varepsilon$. Using the chiral decomposition $\psi\equiv\begin{bmatrix}L&R\end{bmatrix}$, $\bar{\psi}\equiv\begin{bmatrix}\bar{R}&\bar{L}\end{bmatrix}$, we get a perturbed two-dimensional free-fermion conformal field theory

[TABLE]

where repeated indices are summed over $i,j\in\{+,-\}$, and where the free-fermion stress tensor is

[TABLE]

with $z=x+iy$ and $\partial=(1/2)(\partial_{x}-i\partial_{y})$. The Euclidean metric in conformal coordinates is

[TABLE]

Eq. (24) shows that the effects of disorder couple gravitationally via the metric perturbation $\delta G^{ij}$, while nonzero energy $\varepsilon$ couples to an (imaginary or tachyonic) mass operator for the fermions. The mass must be imaginary in order to ensure plane-wave correlations in the clean limit. The components of the perturbing 2D metric are

[TABLE]

The five variants of the model (a)–(e) were defined and described below Eq. (1). In the conformal language [Eq. (24)], models (c) and (d) are special. Model (c) has $\delta G^{++}=\delta G^{--}=0$; in this case, disorder couples only to the “diagonal” components of the stress tensor $\{T_{+-},T_{-+}\}$. These are operators with conformal spin $s=0$; the disorder induces isotropic flattening or steepening of the Dirac cone, as well as local rotations of the Dirac pseudospin $\{\hat{\sigma}^{1,2}\}$ relative to the coordinate axes $\{x_{1,2}\}$. By contrast, model (d) has $\delta G^{+-}=\delta G^{-+}=0$. In this case disorder couples only to the purely chiral $T(z)$ and $\bar{T}(\bar{z})$ fields, which carry conformal spins $s=\pm 2$. These “$d$-wave” operators describe nematic deformations of the Dirac cone. We will see that the nematic-only model (d) shows very robust “critical stacking” behavior, while stacking in model (c) is disrupted by increasing the disorder strength.

QGD can affect massless 2+1-D Dirac carriers whenever the latter arise from a correlation gap. Strong spatial inhomogeneity has been observed in gap maps of the superconducting and pseudogap regimes in the $d$-wave cuprate superconductor BSCCO Davis01 ; Davis02 ; Davis05-a ; Davis05-b ; Davis08 ; DavisReview . This should in turn imply the modulation of the nodal quasiparticle velocities along the Fermi surface [see Fig. 1(c)].

Additional sources of internode scattering in the cuprates can arise due to short-ranged impurities such as interstitial oxygen dopants DavisReview ; AltlandSimonsZirnbauer02 ; Hirschfeld02 ; LeeMott06 ; these can dominate over the effects of velocity modulation. Depending on whether one, two, or four nodes are coupled in the $d$-wave problem, one can get critical states all the way down to zero energy Ostrovsky2007 ; Chou2014 ; Ghorashi18 (with a concomitant sublinear scaling of the low-energy global density of states) Ludwig1994 ; Nersesyan1994 ; Mudry1996 ; Caux1996 ; Bhaseen2001 . The other possibility is Anderson localization across the energy spectrum for the most generic model of impurity scattering Senthil98 ; AltlandSimonsZirnbauer02 ; see Sec. II.2 and Table 1 for a review of these models and scenarios. We only emphasize two points here. First, the linear low-energy density of states and the coherent, plane-wave-like nature of the lowest-energy quasiparticles observed in STM studies Davis05-a ; Davis05-b ; Davis08 ; DavisReview appear inconsistent with strong internode scattering. Second, the experimentally-observed thermal conductivity Taillefer00 ; Orenstein00 is not consistent with localization; in fact, the “universal result” in Eq. (13) is predicted to hold for any topologically restricted model of scattering, including Eq. (22) (see Sec. II.2).

Eq. (22) could be also realized on the surface of a class DIII topological superconductor (TSC) TSCRev1 ; TSCRev2 ; TSCRev3 ; TSCRev4 with winding number $|\nu|=1$ Nakai14 ; Roy19 . In this case $\psi$ is a real Majorana field with $\bar{\psi}_{\sigma}=i\psi_{\sigma^{\prime}}\left(\hat{\sigma}^{1}\right)_{\sigma^{\prime}\sigma}$, and one can show that (see Appendix B)

[TABLE]

where $v_{0}$ is the bare surface Majorana velocity, $A^{0}(\bm{\mathrm{r}})$ is the screened electric potential due to (e.g.) static charged impurities, $\vartheta$ is a constant, and $E_{\mathsf{bulk}}$ is the bulk gap energy of the TSC. The disorder can be characterized by a variance $\tilde{\lambda}\propto n_{\mathsf{imp}}\left[e^{2}/\left(k_{F}k_{\mathsf{TF}}\right)\right]^{2},$ where $n_{\mathsf{imp}}$ is the surface areal impurity density, and $k_{F}$ ($k_{\mathsf{TF}}$) is the bulk Fermi (Thomas-Fermi screening) wave vector. In units such that $v_{0}=1$, $\tilde{\lambda}$ is a squared-length; weak, short-range correlated QGD is therefore strongly irrelevant KPZ–Note at zero energy on the surface. Low-energy states are thus expected to be only weakly affected by QGD Nakai14 . The fate of the finite-energy states is not obvious, however, since energy is itself strongly relevant Ludwig1994 ; Ghorashi18 .

III.2 Density of critical states (DOCS) and disorder-strength scaling

In order to work directly in the continuum, we diagonalize Eq. (1) exactly in momentum space,

[TABLE]

We set the parameter $\nu=1$ in Eq. (III.2), appropriate to the relativistic clean system. A higher odd-integral value of $\nu\in\{3,5,\ldots\}$ can be used to represent a 2D class DIII Majorana surface fluid that arises from a bulk topological superconductor with corresponding winding number $\nu$ Roy19 . We used Eq. (III.2) to analyze the multifractal spectra of low-energy surface states with $\nu\in\{3,5,7\}$ in Ref. Roy19 , where we tested predictions of the class DIII SO($2n$)ν ($n\rightarrow 0$) conformal field theory expected to describe the zero-energy surface states in these cases Foster14 .

The impurity velocity potentials in Eq. (III.2) are each taken to be a composition of random phases in momentum space, e.g.

[TABLE]

Here $\theta_{11}(-\bm{\mathrm{k}})=-\theta_{11}(\bm{\mathrm{k}})$; these are otherwise independent, uniformly distributed phase angles. We choose a short correlation length so as to approximate white noise disorder, appropriate (e.g.) to model the nanoscale gap inhomogeneity observed in the cuprates DavisReview , or efficient screening of Coulomb impurities on the surface of a TSC: $\zeta\equiv(0.25)(2\pi/\Lambda)$, where $\Lambda$ is the momentum cutoff. Disorder becomes “strong” when the local variance of the velocity components in position space ($\equiv\Delta v_{ij}$) becomes of order one. Disorder beyond this threshold regularly tips the velocity components through zero, which can create curvature horizons Volovik ; Jafari19 [see Eq. (44)]. Since $\Delta v_{ij}=\sqrt{\tilde{\lambda}_{ij}/(2\pi\zeta^{2})}$, this corresponds to the condition $\lambda_{ij}\equiv\tilde{\lambda}_{ij}(\Lambda/2\pi)^{2}=\pi/8\simeq 0.393$. Here $\lambda_{ij}$ denotes the dimensionless disorder strength.

We study the five model variants (a)–(e) defined below Eq. (1). For simplicity, we assign the same dimensionless disorder variance $\lambda$ to all independent, nonzero velocity potentials in each variant.

Representative plots of the LDOS for states at different energies in models (b) and (e) appear in Figs. 2 and 3. Wave function quantum criticality is characterized by the spectrum of exponents $\tau(q)$, reviewed in Sec. II.2.2 [Eqs. (17) and (18)]. We calculate the multifractal spectrum $\tau(q)$ by the usual box-counting method; the reader is referred to Ref. Ghorashi18 for technical details. In order to quantify the degree of criticality throughout the energy spectrum, we employ the following criterion. We compare the computed $\tau(q)$ spectrum for every state in regularly spaced energy bins to a quadratic ansatz Ghorashi18 , $\tau(q)=2(q-1)(1-q/q_{c}^{2})$ for $|q|\leq q_{c}$. We employ the “fitness” criteria, defined as follows Ghorashi18 . For each eigenstate $\psi(\bm{\mathrm{r}})$, we compute the error between the numerical spectrum [$\equiv\tau_{N}(q)$] and the appropriate analytical prediction [$\equiv\tau_{A}(q)$], error$(q)\equiv|\tau_{N}(q)-\tau_{A}(q)|/\tau_{A}(q)$. If the error is less than or equal to 4% for 85% of the evaluated $q$-points in the interval $0\leq q\leq q_{c}$, we keep the state. We consider 200 total bins across the energy spectrum. In each bin we analyze 15 states with equally spaced eigenenergies; the DOCS is computed from the fraction of these that satisfy the fitness criterion.

We empirically choose $q_{c}=5.1$ for the parabolic ansatz; this corresponds to $\theta\simeq 1/13$ in Eq. (18). Representative anomalous multifractal spectra $\Delta(q)\equiv\tau(q)-2(q-1)$ are shown in Figs. 10, 11, and 12 for the models (b), (d), and (e) that exhibit robust critical stacking. When we scan through states at different energies to evaluate their “fitness” relative to the universal critical spectrum conjectured above, we exclude negative moments $q<0$, since evaluating these accurately requires significant coarse-graining; for this reason negative moments are typically not reported. We are unable to determine if the deviation seen between the ansatz and the data for $q<0$ in Figs. 10–12 is intrinsic, or simply a finite-size limitation.

We define the density of critical states (DOCS) as the number of states within an energy bin satisfying the above criterion. The ratio of the DOCS to the total density of states (DOS) is the effective energy-resolved distribution function for critical states Ghorashi18 . For models (a)–(e), the DOS and the DOCS are shown in Figs. 5–9, respectively. In each case, we show results for the maximum system size (193$\times$193 momentum grid) and six different disorder strengths $\lambda\in\{0.05,0.1,0.2,0.3,0.4,0.5\}$; $\lambda\gtrsim 0.393$ corresponds to strong disorder (see above). In these plots, we also exhibit the bare second inverse participation ratio (IPR)

[TABLE]

This gives a raw measure of the spatial extent of each wave function. For a plane wave $\mathcal{P}_{2}\sim 1/L^{2}$, while for a localized state with $\zeta_{\mathsf{loc}}\ll L$, $\mathcal{P}_{2}\sim 1/\zeta_{\mathsf{loc}}^{2}$, where $L$ and $\zeta_{\mathsf{loc}}$ denote the system size and localization length, respectively. More generally, $(L^{2})\mathcal{P}_{2}$ is enhanced for any rarification of the wave function (relative to the uniform plane wave case).

All five models exhibit a swath of critical states at finite energy for sufficiently weak disorder. As we discuss below, we find very weak system-size dependence for these results. A crucial difference is how the critical swath varies as the strength of the disorder is increased. Figs. 5 and 7 show that the swath shrinks for models (a) and (c) with increasing disorder; this is especially pronounced for model (c). The second IPR $\mathcal{P}_{2}$ in these plots indicates that states outside of the critical swath are more rarified (supercritical), except near zero energy for the weakest disorder strength. As described below Eq. (1), models (a) and (c) both feature disorder that locally flattens or steepens the Dirac cone in an isotropic fashion; model (a) also incorporates nematic fluctuations.

By contrast, Figs. 6 and 8 show that the width of the critical swath for models (b) and (d) increases with increasing disorder strength. In model (b), in particular, the main effect of increasing the disorder is to lower the crossover energy ($\equiv\Delta_{0}$) between low-energy plane-wave-like states and finite-energy critical states. The coexistence of low-energy plane-waves and finite-energy critical states (with fixed, universal multifractal fluctuations in the latter) over a wide range of disorder strengths means that model (b) is closest to the phenomenology of the cuprates, i.e. to the energy-dependence of the LDOS maps observed in STM studies of BSCCO Davis08 ; DavisReview . Model (d) behaves in a similar fashion, except that the states near zero energy begin to show more rarified behavior (larger $\mathcal{P}_{2}$) at intermediate disorder strengths, even while the transition to “stacked” criticality decreases in energy. In model (d), this is concomitant with a filling in of the total DOS at zero energy. Models (b) and (d) feature nematic disorder; model (d) corresponds to the pure “$T\bar{T}$” deformation in the CFT language of Eq. (24) (after averaging over disorder, using e.g. replicas). Model (b) additionally incorporates local rotations.

The DOS and DOCS for the generic model (e) are depicted in Figs. 9. In this case the behavior is intermediate: a wide critical swath appears for all disorder strengths, but is pushed to higher and higher energies with increasing $\lambda$. Except for the weakest disorder, states near zero energy become more rarified than those in the critical swath, as indicated by the enhancement of $\mathcal{P}_{2}$.

The anomalous multifractal spectra selected from an energy bin where the ratio of the DOCS to DOS is maximized for models (b), (d), and (e) are plotted for the same six disorder strengths in Figs. 10–12. Note that despite the variation of $\lambda$ by an order of magnitude, there is little change in the spectra. This is completely different from the conventional symmetry-based prediction reviewed in Appendix C, which argues that finite-energy states of the class DIII model in Eq. (1) should reside in the symplectic class Evers2008 . In that case, although multifractality would be expected for weakly antilocalizing states in finite volume, the curvature of the spectrum is predicted to scale linearly in the disorder strength or inverse bare (semiclassical) conductance [see Eqs. (35) and (36), below].

III.3 Plane-wave to critical-state crossover energy “$\Delta_{0}$” versus disorder strength

A key observation in STM studies of BSCCO in the optimal to underdoped regime is the qualitative separation between local density of states (LDOS) spectra obtained below and above a weakly doping-dependent scale $\Delta_{0}$. At energies below (above) $\Delta_{0}$, the LDOS maps exhibit robust energy-dispersing quasiparticle interference (strong energy-independent, nanometer-scale spatial inhomogeneity), interpreted as separating “coherent” low-energy quasiparticle states from “incoherent” intermediate-energy excitations Davis08 ; DavisReview .

For the QGD model (b) [defined below Eq. (1)] that best fits the cuprate phenomenology (see Figs. 2 and 6), we can try to compute the scale $\Delta_{0}$ in terms of the energy cutoff $\Lambda$, as a function of the dimensionless disorder strength $\lambda$. We define $\Delta_{0}$ as the threshold where the DOCS becomes larger than a certain percentage of the total DOS. We choose the (arbitrary) criterion

[TABLE]

in order to define $\Delta_{0}$, measured relative to the Dirac point. Results are shown in Fig. 13. We find a weak decrease of $\Delta_{0}$ with increasing disorder strength. This is qualitatively similar to cuprate data Davis08 ; DavisReview with increased underdoping, if the magnitude of the fluctuations in $\Delta_{1}(\bm{\mathrm{r}})$ is taken as a proxy for the effective disorder strength (instead of the doping level).

For the class DIII model (b), we are unable to extract a clear trend of $\Delta_{0}$ with respect to the available system sizes $N$ (see Fig. 14). In the other studied instances of quantum-critical wave-function stacking, the disorder is of vector-potential type Foster14 (see Table 1 in Sec. II). For short-range correlated vector-potential randomness, the disorder strength is intrinsically dimensionless. This leads to the expectation that the crossover energy scale analogous to $\Delta_{0}$ must go to zero as the system size $N\rightarrow\infty$, consistent with the numerical results obtained for class CI and AIII surface states in Refs. Ghorashi18 ; Sbierski . For QGD studied here, however, the intrinsic disorder strength $\tilde{\lambda}$ has units of length squared. Then we can write

[TABLE]

where $L$ is the dimensionful linear system size, and we have assumed that the scaling function $F$ does not depend explicitly on the cutoff $\Lambda$ (as usual GoldenfeldBook ). Clearly we must have $\Delta_{0}\rightarrow\infty$ as $\tilde{\lambda}\rightarrow 0$. We cannot rule out that the limit $F(x\rightarrow 0)>0$, which would imply a finite $\Delta_{0}\sim\tilde{\lambda}^{-1/2}$ in the $L\rightarrow\infty$ limit. As evidenced by the roughness of the trend indicated in Fig. 13, much more work is needed to establish the behavior of $\Delta_{0}$ with respect to disorder and system size.

In the QGD stacking scenario, a finite $\Delta_{0}$ would seemingly imply a finite-energy semimetal to quantum-critical stacked transition. A crucial goal for future work is to establish whether this transition exists at $\varepsilon>0$, and if so what independent detectable signatures could be associated to it. In experiments on BSCCO, $\Delta_{0}$ can also be extracted from a “kink” in the energy-dependence of the spatially averaged tunneling spectrum DavisReview . The precise variation of the DOS at energies away from $\varepsilon=0$ is however likely nonuniversal and model-dependent Evers2008 , in contrast to the universal spatial fluctuations of the critical states in the stacking scenario Ghorashi18 ; Sbierski . Finally, we note that a kink in the DOS observed in experiments at the threshold of critical stacking might arise due to the interplay with quasiparticle-quasiparticle interactions (i.e., an Altshuler-Aronov-type correction) LeeRamakrishnan .

III.4 System-size scaling of the DOS and DOCS

As emphasized above, weak QGD is strongly irrelevant near the Dirac point. For an electrically charged Dirac field $\psi$, kinetic theory predicts a divergent dc conductivity for all nonzero temperatures DavisFoster . The finite-energy critical-state stack found here instead suggests a finite dc conductivity.

On the other hand, conventional symmetry arguments imply that the finite-energy states of a class DIII model should reside in the symplectic class AII Evers2008 , see Appendix C. This class can exhibit weak antilocalization, leading to a “supermetallic” phase in 2D Bardarson07 ; Nomura07 . In a finite-size system, it is therefore not surprising to find multifractal wave functions in 2D, since the scaling is only logarithmic in the system size. However, the curvature of the multifractal spectrum in that case should be related to the bare disorder strength. Moreover, more of the spectrum should become more weakly multifractal with increasing system size.

The leading order anomalous multifractal spectrum for a 2D symplectic metal is given by HoefWegner86

[TABLE]

where the dimensionless conductance $G(L)$ scales according to weak antilocalization LeeRamakrishnan

[TABLE]

Here $L$ is the system size and $l_{\mathsf{el}}$ denotes the elastic impurity scattering length. The bare conductance $G_{0}=\nu(\varepsilon)D(\varepsilon)$, where $\nu(\varepsilon)$ is the density of states and $D(\varepsilon)$ is the diffusion constant, at the single-particle energy $\varepsilon$. Fixing the disorder strength fixes $G_{0}$ at each energy. Increasing the system size $L$ should then reduce multifractality according to Eq. (35), as the system flows slowly towards a “supermetallic” phase at infinite $L$.

We plot the evolution of the density of states (DOS) and density of critical states (DOCS) as a function of system size for models (b), (d), and (e) in Figs. 14–16. These plots exhibit no discernible systematic variation with system size. This is not consistent with weak antilocalization behavior expected from the symplectic metal class AII in Eqs. (35) and (36).

Finally, in Figs. 17 and 18, for the robust stacking models (b) and (d) we exhibit the finite-size scaling of two particular multifractal dimensions with system size. Results are shown for the same 6 disorder strengths $\lambda$ employed in Figs. 6, 8, 10, and 11. At each system size $N$, the dimensions are computed from an energy bin wherein the DOCS is maximized relative to the DOS. These figures indicate a possible very weak dependence of the dimensions on the disorder strength. This is very different from the symplectic metal prediction in Eqs. (35) and (36), which would instead imply a linear dependence of the dimensions on the disorder strength (proportional to $1/G_{0}$).

III.5 Discussion and open questions

We have uncovered a third instance of quantum-critical wave function stacking Ghorashi18 , due here to QGD. In both Ref. Ghorashi18 (class CI) and the present paper (class DIII), the stacking occurs in disordered Dirac models that can describe the 2D surface states of 3D topological superconductors. Another instance (for topological superconductor surface states in class AIII) was touched upon in Ref. Chou2014 , and extensively explored in Ref. Sbierski .

We argued that QGD naturally obtains whenever 2D Dirac quasiparticles arise from a spatially inhomogeneous gap. The critical-state stack found in this paper might provide a simple explanation for the division between plane-wave-like and spatially inhomogeneous LDOS fluctuations measured at low and finite energies, respectively, in the high-$T_{c}$ cuprates.

We studied the five model variants (a)–(e), defined below Eq. (1). We found the best qualitative agreement with cuprate STM phenomenology DavisReview for model (b), which includes random pseudospin rotations and quenched nematic fluctuations of the Dirac cone. Robust stacking that increases with increasing disorder strength is observed in both model (b) and model (d), the latter which includes only pure nematic disorder. By contrast, isotropic flattening or steepening of the Dirac cone appears to suppress critical stacking for intermediate and strong disorder [models (a) and (c)]. The generic model (e) shows intermediate behavior, with a robust critical stack at all disorders, but low-energy states that become more rarified (supercritical) with increasing disorder strength.

Avenues and questions for future work include the following. (1) The anomalous multifractal spectrum for spatial LDOS fluctuations computed for the critical stack, exhibited in Figs. 10–12, could be directly compared to corresponding spectra computed from experimental STM data on BSCCO DavisReview . (2) A crucial question is whether the crossover scale $\Delta_{0}$ [Fig. 13 and Eq. (34)] remains finite or vanishes in the infinite system-size $L\rightarrow\infty$ limit. A nonzero $\Delta_{0}(L\rightarrow\infty)$ would indicate a finite-energy transition between ballistic and quantum-critical stacked behavior. (3) Another question regards transport versus temperature in the Dirac model with QGD. In particular, could the transition from “ballistic” low-energy states to critical, finite-energy ones explain the ubiquitous linear-in-$T$ resistivity observed in the strange metal phase above $T_{c}$ LeeMott06 ? (4) What is the role of dephasing on transport AAK ; IQHP-Deph due to inelastic quasiparticle scattering, which is presumably important in the cuprates due to the strong correlations ($U\gg t$), and finally (5) Does the multifractal-stacking phenomenon enhance superconductivity in a self-consistent calculation of the gap?

Acknowledgements.

We thank Alex Altland, Mustafa Amin, Seamus Davis, Ilya Gruzberg, Victor Gurarie, Christopher Hooley, Enrico Rossi, and Björn Sbierski for useful discussions. S.A.A.G. was supported by the U.S. Army Research Office Grant No. W911NF-18-1-0290, and acknowledges partial support from NSF CAREER Grant No. DMR-1455233 and ONR Grant No. ONR-N00014-16-1-3158. J.F.K. acknowledges the KHYS research travel grant and Graduate Funding from the German States. S.M.D. and M.S.F. acknowledge support by the Welch Foundation Grant No. C-1809 and by NSF CAREER Grant No. DMR-1552327. M.S.F. also acknowledges support by the U.S. Army Research Office Grant No. W911NF-17-1-0259. M.S.F. thanks the Aspen Center for Physics, which is supported by the NSF Grant No. PHY-1607611, for its hospitality while part of this work was performed.

Appendix A Random velocity modulation from Dirac fermions in curved spacetime

The generally covariant action for (2+1)-D Dirac fermions propagating through curved spacetime is given by Eq. (23). In the context of disordered Dirac superconductors with QGD, however, there is a naturally preferred coordinate system, where $\bm{\mathrm{r}}=\{x_{1},x_{2}\}$ measure physical Euclidean distances across the sample. In this “physical” coordinate system, the Berry phase term $\bar{\psi}\hat{\gamma}^{0}i\partial_{t}\psi$ arises from Trotterization in the path integral, and should therefore have unit coefficient. We call this the “temporal flatness condition.” Then, we can take the dreibein $E_{A}^{\mu}$ in Eq. (23) to be off-diagonal only in the spatial-spatial sector. Let $i,j\in\{1,2\}$. We define

[TABLE]

Temporal flatness requires

[TABLE]

Shifting $\bar{\psi}\rightarrow\bar{\psi}\hat{\gamma}^{0}$ and enforcing Eqs. (37) and (38), Eq. (23) reduces to Eq. (22). The spin connection is eliminated via integration-by-parts in the spatially modulated kinetic terms; this leads to the manifestly Hermitian form of the Hamiltonian in Eqs. (1).

The inverse metric can be expressed as

[TABLE]

where $\eta^{AB}\rightarrow\operatorname{diag}(-1,1,1)$ is the Minkowski metric, and

[TABLE]

Imposing Eq. (38), the metric is

[TABLE]

This can be specialized for the five model variants (a)–(e), defined below Eq. (1). For example, model (a) has $v_{12}=v_{21}=0$, leading to

[TABLE]

while model (b) has $v_{11}=v_{22}=1$, and

[TABLE]

In the generic case, the scalar curvature can be written has

[TABLE]

where $N$ is a homogeneous quadratic in spatial derivatives of $\{v_{ij}\}$. As an example, for model (a) with $v_{11}=v_{22}\equiv v(\bm{\mathrm{r}})$,

[TABLE]

Appendix B “Gravitational” coupling of electric potentials to the surface Majorana fluid of a class DIII topological superconductor

In this section we derive the form of the velocity modulation in Eq. (29), corresponding to the effect of an electric potential $A^{0}(\bm{\mathrm{r}})$ on the 2D Majorana fluid expected to form at the surface of a class DIII topological superconductor.

B.1 Bulk and surface states for solid-state 3He-$B$

As a simple model for a class DIII bulk topological superconductor with winding number $|\nu|=1$, we consider a solid-state analog of 3He-$B$ Volovik ; TSCRev2 . The bulk Bogoliubov-de Gennes Hamiltonian features isotropic $\hat{\bm{\sigma}}\cdot\bm{\mathrm{k}}$ pairing, where $\hat{\bm{\sigma}}=\hat{\sigma}^{a}\hat{x_{a}}$ is the vector of Pauli matrices acting on the physical spin-1/2 components $\sigma\in\{\uparrow,\downarrow\}$, and $a$ is summed over $\{1,2,3\}$. The mean-field Hamiltonian is

[TABLE]

where $\mu>0$ is the chemical potential and $\Delta$ is the $p$-wave pairing amplitude. Here we have introduced the Balian-Werthammer (“Majorana”) spinor

[TABLE]

The Pauli matrices $\{\hat{\tau}^{1,2,3}\}$ act on particle-hole space. Particle-hole $P$, time-reversal $T$, and chiral $S$ ($\equiv T\times P$) symmetries are defined via

[TABLE]

consistent with class DIII SRFL08 ; Foster14 .

As in Ghorashi17 , we implement hardwall boundary conditions at $z=0$ in order to get surface states. Eq. (46) separates into a $\bm{\mathrm{k}}=0$ piece, and a nonzero $\bm{\mathrm{k}}$ piece, where $\bm{\mathrm{k}}\equiv\{k_{x},k_{y}\}$ now accounts only for conserved transverse momenta.

[TABLE]

We have included a scalar electric potential $A^{0}(\bm{\mathrm{r}})$ as a perturbation in $\hat{h}_{1}$.

The Hamiltonian $\hat{h}_{0}$ has a pair of zero-energy Majorana bound states,

[TABLE]

where $|{\tau^{1}=m_{s}}\rangle$ is the particle-hole ($\tau$) space spinor, which is “locked” to the $\hat{\sigma}^{z}$-spin projection $m_{s}$ (in the plus and minus $\hat{\tau}^{1}$-direction for $m_{s}=\uparrow$ and $\downarrow$, respectively). The spatial profile of the bound state is

[TABLE]

where ${\mathcal{N}_{{0}}}$ is a normalization factor.

B.2 $k\cdot p$ perturbation theory

An effective Hamiltonian for the surface theory at $\bm{\mathrm{k}}\neq 0$ obtains by taking matrix elements of the following operator between the $m_{s}=\pm 1$ zero modes in Eq. (50),

[TABLE]

where we have expanded the second term via the $\bm{\mathrm{k}}=0$ resolution of the identity. The operator $\hat{P}_{1}$ on the first line is the projection out of the degenerate eigenspace of the zero modes, $\hat{P}_{1}=\hat{1}-\hat{P}_{0}$, where

[TABLE]

The state $|{\psi_{q,m_{s}}}\rangle$ is a positive-energy (gapped) bulk eigenstate of $\hat{h}_{0}$, parameterized by the standing wave momentum $q$, while $\hat{\tau}^{1}|{\psi_{q,m_{s}}}\rangle$ is its negative-energy chiral conjugate [Eq. (48)]. $E_{q}=\sqrt{(q^{2}/2m-\mu)^{2}+q^{2}\Delta^{2}}$ denotes the positive eigenenergy.

The first term in Eq. (B.2) gives the relativistic dispersion for the Majorana surface fluid,

[TABLE]

where $\bm{\mathrm{A}}\wedge\bm{\mathrm{B}}=A_{x}B_{y}-A_{y}B_{x}$. This is consistent with the surface projection of the symmetry conditions in Eq. (48),

[TABLE]

Coupling to the vector potential $A^{0}$ obtains from the second term in Eq. (B.2). Working to linear order in the potential, the relevant $+-$ matrix elements take the form

[TABLE]

and

[TABLE]

where $\mathsf{k}\equiv k_{x}-ik_{y}$. Evaluating these leads to perturbation of the form

[TABLE]

where $E_{\mathsf{bulk}}\simeq k_{F}\Delta$ is the bulk excitation gap, and $\vartheta$ is a pure order-one number. Since the bare Majorana fluid velocity is $\Delta$, we recover Eq. (29).

Appendix C Symmetry class for finite-energy states; Connection to the class D thermal quantum Hall plateau transition (?)

The 2D velocity-randomized Dirac Hamiltonian given by the sum $\hat{h}_{\scriptscriptstyle{\mathsf{S}}}^{{\scriptscriptstyle{({1})}}}+\hat{h}_{\scriptscriptstyle{\mathsf{S}}}^{{\scriptscriptstyle{({2})}}}$ in Eqs. (53) and (57) resides in class DIII, due to the $T^{2}=-1$ and $P^{2}=+1$ time-reversal and particle-hole symmetries encoded in Eq. (54). These symmetries hold irrespective of whether the fermion field $\psi$ in Eq. (22) is a complex-valued Dirac or real-valued Majorana spinor.

Typically, we can associate an effective field theory, the nonlinear sigma model, to describe the wave functions of any single-particle Hamiltonian at some particular fixed energy. The nonlinear sigma model employs local operators to encode the probability statistics of the extended, critical, or localized states. Using fermionic replicas to perform disorder-averaging, class DIII is associated to a sigma model with the target manifold O$(2n)$, where $n$ is proportional to the number of replicas Evers2008 . In two spatial dimensions, this can be seen via the nonabelian bosonization of the clean, zero-energy Majorana fermion field theory Foster14 .

For the DIII theory, nonzero energy is a relevant perturbation that couples to the principal chiral field Foster14 ; Ghorashi18 in the nonlinear sigma model. In this case, the O($2n$)$\times$O($2n$) symmetry of the zero-energy theory is broken down to the diagonal subgroup O($2n$). The energy perturbation will induce an RG flow to a new fixed point, which should be associated to a different sigma model with target manifold O($2n$)/$H$. There are only two possibilities:

[TABLE]

Class AII, the “symplectic” class associated to disordered metals with time-reversal symmetry and strong spin-orbit coupling. This class typically exhibits weak antilocalization in two dimensions. Class AII also describes the 2D surface states of 3D topological insulators. The target manifold is Evers2008

[TABLE] 2. 2.

Class D, typically associated to superconductors with broken time-reversal symmetry and strong spin-orbit coupling. Class D should be realized in dirty $p+ip$ superconductors D-1 ; D-2 ; D-5 ; D-6 ; D-7 . The target manifold is

[TABLE]

We argued in Sec. III that class AII appears incompatible with the numerical results obtained here for finite-energy states of the velocity-modulated Dirac Hamiltonian in Eq. (1). We therefore expect that our finding of critical states throughout the energy spectrum with universal multifractal spectra should instead be associated to class D.

There are three classes of time-reversal invariant topological superconductors (TSCs) in 3D, differing by the amount of spin rotational symmetry: CI [SU(2)], AIII [U(1)], and DIII (no spin symmetry, strong spin-orbit coupling) ClassRevNote . There are also three classes of time-reversal broken quantum Hall topological insulators or TSCs in 2D: class C [spin quantum Hall (SQH) effect], class A [ordinary integer quantum Hall (IQH) effect], and class D [thermal quantum Hall (TQH) effect] Evers2008 . All of these classes are characterized by integer (AIII, DIII, A, D) or twice-integer (CI, C) bulk topological winding numbers ClassRevNote .

In our previous work Ghorashi18 , we performed a similar population analysis as presented in this paper for the finite-energy 2D surface states of a 3D class CI topological superconductor. Based on the numerical results obtained there, we concluded that the finite-energy states of the 2D class CI model take the form of a “stack” of critical wave functions with universal multifractal spectra. The spectra are energy independent, and consistent with the plateau transition of the spin quantum Hall effect in class C SQHP3 ; EversSQHE03 ; MirlinSQHE03 . A connection between classes CI and C based on symmetry considerations similar to Eq. (58) was also presented in Ref. Ghorashi18 .

Finite-energy states in class AIII must reside in class A Ludwig1994 . For the 2D surface states of a 3D class AIII topological superconductor, whether these are critically delocalized or Anderson localized depends upon the presence or absence of a topological theta term (at $\theta=\pi$) in the effective nonlinear sigma model Pruisken . In Ref. Sbierski , strong numerical evidence was presented that 2D finite-energy surface states of 3D class AIII TSCs form a “stack” of critical wave functions [see also Ref. Chou2014 ]. The spectra are energy independent, and consistent with the plateau transition of the ordinary class A integer quantum Hall effect Huckestein1995 . For classes CI and AIII, the only alternative to the plateau-transition-stacking scenario is Anderson localization, but this is not observed in the numerical studies Ghorashi18 ; Chou2014 ; Sbierski .

We therefore expect that the “stack” of critical states found in the present paper can be associated to the thermal quantum Hall plateau transition in class D. In comparison to classes C and A, relatively little is known about the thermal quantum Hall plateau transition. The global phase diagram for a 2D system in class D is complicated by the advent of a thermal metal phase D-1 ; D-2 ; D-5 ; D-6 ; D-7 , in addition to Anderson localized thermal Hall plateaux Evers2008 . Further studies of the possible multicritical point in the phase diagram of class D D-1 ; D-5 ; D-6 could shed light on the nature of the finite-energy states found here.