Do the surface Fermi arcs in Weyl semimetals survive disorder?

We theoretically study the topological robustness of the surface physics induced by Weyl Fermi-arc surface states in the presence of short-ranged quenched disorder and surface-bulk hybridization. This is investigated with numerically exact calculations on a lattice model exhibiting Weyl Fermi arcs. We find that the Fermi-arc surface states, in addition to having a finite lifetime from disorder broadening, hybridize with nonperturbative bulk rare states making them no longer bound to the surface (i.e., they lose their purely surface spectral character). Thus, we provide strong numerical evidence that the Weyl Fermi arcs are not topologically protected from disorder. Nonetheless, the surface chiral velocity is robust and survives in the presence of strong disorder, persisting all the way to the Anderson-localized phase by forming localized current loops that live within the localization length of the surface. Thus, the Weyl semimetal is not topologically robust to the presence of disorder, but the surface chiral velocity is.

Figure 1

(a) Schematic of a Weyl semimetal with two cones in the bulk and therefore a single Fermi arc. The chiral velocity is defined perpendicular to the arc. (b) Schematically, how the density of states and (log of) the surface chiral velocity change in the phases and regimes this model exhibits (horizontal axis). Setting $m=(3/2)t$, both quantities are evaluated as disorder averages at $E=0$. The diffusive metal phase lies within $0<W<W_l$ and the Anderson insulator for $W>W_l$. (c) and (d) For $m=(3/2)t$, we plot cuts of the dispersion for Eq. (1) in the clean limit with open boundary conditions displaying the bulk bands and the topological surface Fermi-arc states (red and blue) dispersing like $E(k_y,k_z)=\pm tsin\left(k_y\right)$ in the pseudogap. (c) $E(k_y,k_z)$ versus $k_y$ with $k_z=0$, and (d) is $E(k_y,k_z)$ versus $k_z$ with $k_y=0$; Weyl points at ${\mathbf{K}}_W=(0,0,\pm 2\pi /3)$ can be seen.

Figure 2

(a) Surface electronic dispersion curves (EDCs) where each value of $A_S(k_y,k_z,\omega )$ for $k_z=0$ with $k=k_y$ is shifted by $(L/2\pi )k$; blue is the top surface and the bottom is red. (b) The spectral function versus $\omega$ on the surface at ${\mathbf{k}}_{\perp }=0$ for various disorder strengths. We see a smooth broadening of the Fermi-arc peak with disorder, as captured by the width of the spectral function $\mathrm{\Gamma }({\mathbf{k}}_{\perp }=0)$ shown in the inset. (c) The typical DOS on the surface for strong disorder and weak disorder (inset). As the system size $L$ increases, the typical DOS on the surface converges to the average surface DOS implying the arcs do not localize for weak disorder; for strong disorder we find the bulk and surface localization transitions agree. (d) Average spectral weight for states of definite momentum on the surface to tunnel into the bulk for three representative surface momenta on the arc: ${\mathbf{k}}_{\perp }=0$, at the Weyl node projection ${\mathbf{k}}_{W,S}=(0,2\pi /3)$, and off the arc ${\mathbf{k}}_{\perp }=(0,\pi )$, all computed at $W=0.5t$ and $L=30$. The finite value of the spectral weight in the middle of the sample indicates surface-bulk hybridization.

Figure 3

(a) The low-energy eigenstates as a function of a twist in the $z$ direction for a disordered sample without any rare states; total weight of the eigenstate on the $x=1$ surface is indicated by the color scale. The indicated state shown in (c) represents hybridization between bulk Weyl states and surface states. Green represents the bulk states found with periodic BCs. (b) The low energy eigenstates as a function of a twist in the $y$ direction for a sample with a rare bulk state; The weight of the wave function on the rare state is indicated by the color scale; green again represents bulk states found with periodic BCs. Opposite chiral velocities represent states on opposing surfaces. The rare state hybridizes with both surfaces (d)–(f), strongly renormalizing the dispersion (b). The density plots are partially summed $\rho (x,y)={\sum }_z|\langle x,y,z|\psi \rangle {|}^2$, and all plots are at weak disorder $W/t=0.5$ and have $L=18$.

Figure 4

(a) Average surface chiral velocity $v_c$ as a function of disorder on a linear and (inset) log scale for various system sizes. $v_S$ can be computed from the pole of the surface Green function only at weak disorder $W\le 0.6t$, we show $L=60$. (b) Broadness of the distribution of the surface velocity given by its standard deviation divided by the mean. The dashed lines are power-law fits to the two distinct power-law regimes, for $W\lesssim 1t$ we find $\sigma \left[v_{c,S}\right]/v_{c,s}\sim W^{1.08}$ and for $W>1t$ it crosses over to $\sigma \left[v_{c,S}\right]/v_{c,s}\sim W^{3.9}$. (c) Average velocity $v_c$ in each $y-z$ sheet located at $x$ for various disorder strengths. Disorder leads to a completely random velocity in the bulk but the surface chiral velocity persists to large disorder (inset) depicting disorder strengths $W=2.0t$ up to $7.0t$ in steps of $1.0t$.

Figure 5

Plots of characteristic near-surface eigenfunctions: Their bulk density profile (left) and surface current (right). The density profile is partially summed $\rho (x,y)={\sum }_z|\langle x,y,z|\psi \rangle {|}^2$, and the surface is at $x=L=18$. The current density profiles are normalized by the largest current value on any bond, with colored squares indicating current flowing onto or out of the surface (i.e., the surface divergence of the current). In (a) a rare state is located and moved close to the surface (indicated by a blue X in the current density profile). Notice that each wave function has more chiral velocity (red, left moving) on the surface, even when the state is delocalized throughout the bulk (c). The fully localized state in (d) has a surface chiral velocity but only as small current loop near the surface (and only near the plane $z=15$).

Figure 6

Tracking the chiral velocity at $E=0$, we see that a bend characterized by $t^{\prime \prime }$ does not affect the surface chiral velocity. These are results on a system of size $L=10$.

Figure 7

(a) The density of states at $E=0$ vs disorder strength $W$ of this system as found from via the KPM method and (b) the saturation of the second derivative of the density of states ${\rho }^{\prime \prime }\left(0\right)$ vs W. The peak of the former characterizes the location of the avoided quantum critical point.

Figure 8

Plot of the typical density of states at $E=0, {\rho }_t\left(0\right)$. (a) Shows the decrease in the typical density of states which mirrors that seen in Fig. 2. (b) Furthermore, notice that dependence on $N_C$ begins to set in around $W_l/t\approx 5.6–6.0$.

Figure 9

(a) The level statistics as they change for higher disorder with periodic boundary conditions. Notice the crossing roughly around $W\approx 6.0t$. (b) The level statistics as they change for higher disorder with open boundary conditions. Notice the crossing roughly around $W\approx 6.0t$. The bins in energy space to determine both plots are roughly $4\%$ of the total bandwidth and symmetric about $E=0$.

Figure 10

Saturation of the width of the peak of the surface spectral function in $\omega$ defined as $\mathrm{\Gamma }\left({\mathbf{k}}_{\perp }\right)=1/\mathrm{Im}{G}_S({\mathbf{k}}_{\perp },\omega =0)$ with the surface wave-vector ${\mathbf{k}}_{\perp }=0$. We are able to converge our results for $W\ge 0.1t$.

Figure 11

(a) Twisting the boundary conditions moves the bulk Weyl states away from zero energy, revealing low-lying rare states; here we focus on the lowest one pictured, $r=309$. It has energy $E\approx 0.067t$ and remains stable when boundary conditions are twisted (b) and (c). This data is for disorder strength $W=0.5t$ and $L=18$.

Figure 12

(a) The absolute size of the wave function on sites a distance $r=|\mathbf{r}-{\mathbf{r}}_{\max }|$ from the maximum. (b) The result of binning the wave function and the red line is a power-law fit to the resulting data: $1/r^{1.83}$.