Dirty higher-order Dirac semimetal: Quantum criticality and bulk-boundary correspondence

We analyze the stability of three-dimensional higher-order topological (HOT) Dirac semimetals (DSMs) and the associated one-dimensional hinge modes in the presence of random pointlike charge impurities.~Complementary numerical and renormalization group (RG) analyses suggest that a HOTDSM, while being a stable phase of matter for sufficiently weak disorder, undergoes a continuous quantum phase transition into a metallic phase at finite disorder. However, the corresponding critical exponents (obtained from the scaling of the density of states) are extremely close to the ones found in a dirty, but first-order DSM (supporting two Fermi arcs). This observation suggests an emergent superuniversality in the entire family of dirty DSMs, as also predicted by the RG analysis. As a direct consequence of the bulk-boundary correspondence, the hinge modes gradually fade away with increasing randomness in the system, and completely dissolve in the trivial metallic phase.

Introduction. Traditionally, the bulk-boundary correspondence in a d-dimensional topological system refers to the boundary modes residing on a (d − 1)-dimensional surface [1][2][3][4][5], also characterized by the codimension d c = d − (d − 1) = 1. This notion has recently been generalized to encompass boundary modes with d c ∈ Z (integers)> 1, such as zero-dimensional corner (d c = d) and one-dimensional hinge (d c = d − 1) modes, which led to the construction of higher order topological (HOT) phases in insulating (electric and thermal) as well as nodal systems . A question regarding the stability of such exotic phase of matter in the presence of interactions and/or disorder arises naturally. Due to a finite gap in the quasiparticle spectra, while the HOT insulators enjoy robustness against sufficiently weak interactions and disorder, their influences on gapless HOT phases demand more careful analyses.
Here we focus on a three-dimensional HOT Dirac semimetal (HOTDSM), supporting one-dimensional hinge modes, and scrutinize its stability when littered with pointlike quenched random charge impurities. Using numerical and field-theoretic renormalization group (RG) analyses, we find that HOTDSM is a stable phase of matter in the presence of sufficiently weak disorder. But, it undergoes a continuous quantum phase transition (QPT) into a diffusive metallic phase at finite disorder, see Fig. 1. We arrive at these conclusions from the scaling of average density of states (DOS), computed using the kernel polynomial method (KPM) [28]. And as a direct consequence of the bulk-boundary correspondence, the topological hinge modes gradually fade away with increasing randomness and completely melt into a trivial metallic bulk for sufficiently strong disorder, see Fig. 2.
Moreover, numerically extracted values for the dynamical scaling exponent z (DSE) and correlation length exponent ν (CLE) across a broad range of parameters (that also includes first-order DSMs) appear to be almost the same (within the numerical accuracy), see Table I, and yield satisfactory single-parameter data collapses across as the ultraviolet cutoff (Fermi velocity). When disorder (W or g) is weak HOT and first-order DSMs are stable. But, at stronger disorder they undergo a QPT into a metallic phase, where DOS at zero energy is finite. The blue dots in (a) are numerically obtained transition points between DSM and metal, across which we find single-parameter scaling of DOS, yielding the data collapses shown in Fig. 3.
the (HOT)DSM-metal continuous QPT, see Fig. 3. This observation strongly promotes the notion of an emergent superuniversality near a diffusive quantum critical point (QCP) in the entire family of dirty DSMs, irrespective of their topological order. We also substantiate these findings from a leading-order RG analysis. Lattice model. We set out by considering a tightbinding model for a three-dimensional HOTDSM on a cubic latticeĥ k =ĥ k 0 +ĥ k 1 , where [16,17,24,26] C j ≡ cos(k j a), S j ≡ sin(k j a), k j being the compo-  (1)], which is not responsible for the hinge states when the crystal is cleaved such that its corners are at (± L 2 , ± L 2 ) for any z. The critical disorder Wc = 3.7 ± 0.1 is obtained from the scaling of average DOS in a periodic system of L = 200. The dissolution of the hinges into the bulk sets in for W < Wc due to the finite size effects.
On the other hand,ĥ k 1 acts as a discrete-symmetry breaking, momentum-dependent or Wilson mass in each layer of QSHI. It changes sign under the four-fold rotation and thus acts as a domain wall mass; resulting in four corner-localized zero-energy states (with d c = 2) in each two-dimensional insulating layer, according to a generalized Jackiw-Rebbi index theorem [29]. Stacking such layers of two-dimensional HOT insulators in the momentum space along the k z axis gives rise to the one-dimensional hinge modes along the z direction [26], see Fig. 2. We then realize a second-order DSM, asĥ k 1 vanishes at the Dirac nodes (±K).
As the Wilson massĥ k 1 vanishes quadratically with momentum around the Dirac points, it only reduces the DOS without altering its overall ρ(E) ∼ |E| 2 scaling at sufficiently low energies, see Fig. 4(a). The subdomi-nant influence of the Wilson mass on the DOS suggests that the HOTDSM is stable in the presence of sufficiently weak disorder, and enters into a metallic phase via a QPT at finite disorder, qualitatively similar to the situation in first-order Dirac and Weyl semimetals , up to exponentially small, but debated rare region effect [53][54][55]. Next we anchor this anticipation by numerically computing the DOS using the KPM, and investigate the critical properties of such a QPT in HOTDSM. For concreteness, here we focus on pointlike charge impurities, which is the dominant source of elastic scattering in any real material, distributed uniformly and independently within the range [−W/2, W/2] at each site of the cubic lattice. Numerically obtained critical exponents are summarized in Table I, which we use to perform a finite energy (size) data collapse, shown in Fig. 3 top (bottom) [56].
Scaling. To formulate the scaling theory for DOS, we concentrate around the dirty QCP at W = W c , and parametrize the distance from it by δ = (W − W c )/W c . The number of states N (E, L) in a d-dimensional system of linear size L below some energy E is in general a function of two dimensionless parameters, L/ξ and E/E 0 . Here ξ is the correlation length, which diverges at the QCP as ξ ∼ δ −ν , and E 0 ∼ ξ −z is the corresponding energy scale. Since the number of states is proportional to L d , the functional form of N (E, L) ought to be [33,40,42] N (E, where F (x, y) is an unknown, but universal function of its arguments. The DOS is then given by To investigate the scaling behavior of the universal function G, we consider the scaling of ρ(E) at low energies  inside the HOTDSM and metallic phases, as well as inside the critical regime around the QCP at W = W c . For now we assume L to be sufficiently large, such that the L-dependence of G can be neglected. For linearly dispersing HOT Dirac fermions the DOS at low energies scales as ρ(E) ∼ |E| d−1 [see Figs. 4(a) and 4(b)] when δ < 0 or W < W c , thus On the other hand, inside the metallic phase the DOS at zero energy is finite, leading to for δ > 0 or W > W c . Lastly, at the critical point (δ = 0) ξ diverges, and therefore the ξ independence of G implies The fact that the DOS at zero energy vanishes for W ≤ W c and becomes finite in the metallic phase, allows one to treat ρ(0) as an order-parameter in dirty HOTDSM. Numerically we reconstruct the average DOS by using the KPM in a cubic system of linear dimension L = 200 for various choices of t and t 1 . First, we identify the critical disorder strength W c , where ρ(0) deviates from zero. Subsequently, we compute (a) the DSE z from the scaling of ρ(E) around W = W c [see Eq.  Table I and the details of the numerical analysis are shown in the Supplemental Materials [56]. Across a wide range of hopping parameters (t and t 1 ) we find that z ≈ 1.5 and ν ≈ 1.0, which are fairly close (within the numerical accuracy) to the ones found for first-order Dirac and Weyl semimetals.
With the numerically extracted values of the critical exponents, we obtain convincing data collapses by comparing ρ(E)|E| −ν(d−z) vs Eδ −νz , see Fig. 3 (top row). All data points collapse onto three curves, corresponding to the DSM (lower ones), a metal (upper left ones) and the quantum critical regime (upper right ones). Finally, from the same set of exponents we obtain excellent finite size data collapses by comparing ρ(0)L d−z with L 1/ν δ inside the metallic phase for 100 ≤ L ≤ 200, as shown in Fig. 3 (bottom row). All together our extensive numerical analyses strongly suggest an emergent superuniversality across the DSM-metal QPT irrespective of their topological order (first or second), which now we substantiate from a leading-order RG analysis.
RG analysis. To perform the RG analysis in a dirty HOTDSM, we consider the following low-energy model obtained by expandingĥ k around one of the Dirac points in general their RG flows are different when |b| > 0, as it breaks the rotational symmetry between p ⊥ = (p 1 , p 2 ) and p 3 . To incorporate the effects of disorder we consider the following imaginary time (τ ) Eucledian action in d dimensions [45,48] Here Φ is the disorder field that minimally couples to the four-component fermionic fields (Ψ) like a gauge field. up to the leading-order, the RG flow equations read [56] where g = ∆Λ /(2π 2 v 2 ) is the dimensionless disorder coupling, x = bΛ/v is also dimensionless, and with S θ ≡ sin θ, C θ ≡ cos θ. Since we are interested in the leading-order RG analysis, (1) the quantum corrections are computed by setting v ⊥ = v 3 = v, and (2) only the engineering dimension of x has been taken into account. The flow equation of v 3 is fixed by its scaling dimension, yielding dynamic scaling exponent ]. The flow equation for g then supports two fixed points: (1) an infrared stable one at g = 0, describing a clean HOTDSM, and (2) an infrared unstable QCP at ]. The latter one controls the QPT into a metallic phase at finite disorder, where z = 1 + /2 = 3/2 for Gaussian white noise distribution ( = 1) for any x. The CLE, defined as ν −1 = d(dg/d )/dg| g=g * = at the dirty QCP. Note that (g * , x) determines the phase boundary between the DSMs and a metal [the blue line in Fig. 1 (b)], which is symmetric about x = 0, as all the functions in Eq. (11) are symmetric under x → −x. In principle, one can also obtain DSE z from the flow equation of v ⊥ , which however does not alter any outcome qualitatively [56]. Finally we note that the four-fold symmetry breaking Wilson mass (x) is an irrelevant parameter. Thus in the deep infrared regime ( → ∞) x → 0, where f ⊥ (0) = f 3 (0) = 2/3, and the velocity anisotropy becomes marginal. Consequently, the ratio v ⊥ /v 3 ultimately flows to its bare value, set by the hopping parameters t and t z . Therefore, the Wilson mass (x) only changes the location of the dirty QCP (g * ) without altering the universality class of the semimetal-metal QPT (determined by ν and z, see Table I), giving rise to a superuniversality in the entire family of dirty DSMs, that includes its first-and second-order cousins.
Discussion. Here we investigate the stability of a HOTDSM in the presence of quenched charge impurities, and identify a semimetal-metal QPT at finite disorder. While the topological hinge modes gradually melt across this transition (similar to the Fermi arcs in conventional dirty Weyl semimetals [57]), see Fig. 2, we come to the conclusion that the topological order (first or second) does not affect its universality class, see Table I. Since real materials are inherently dirty, the stability of HOTDSM is critical for its experimental realization and observation of the hinge modes (via scanning tunneling microscopy, for example) in sufficiently clean systems. We also note that HOTDSM in general is more stable than its first-order counterpart (due to the suppression of DOS by the Wilson mass), see Fig. 1. In the future, it will be worthwhile to investigate the nature of the critical wavefunctions in dirty HOTDSM, and search for measurable signatures of the discrete symmetry breaking.
An interesting limit of the lattice model in Eq. (1) is when t = 0. The system then describes a double firstorder DSM, supporting four Fermi arcs in the k x k z or k y k z plane. However, due to the |E|-linear DOS this system is unstable in the presence of infinitesimal disorder. A finite t converts the double first-order DSM into a second-order DSM, where DOS vanishes as ρ ∼ |E| 2 , and stabilizes the system against sufficiently weak disorder.