Quantum phases of disordered three-dimensional Majorana-Weyl fermions

The gapless Bogoliubov-de Gennes (BdG) quasiparticles of a clean three-dimensional spinless $p_x+i{p}_y$ superconductor provide an intriguing example of a thermal Hall semimetal (ThSM) phase of Majorana-Weyl fermions; such a phase can support a large anomalous thermal Hall conductivity and protected surface Majorana-Fermi arcs at zero energy. We study the effects of quenched disorder on such a gapless topological phase by carrying out extensive numerical and analytical calculations on a lattice model for a disordered, spinless $p_x+i{p}_y$ superconductor. Using the kernel polynomial method, we compute both average and typical density of states for the BdG quasiparticles, from which we construct the phase diagram of three-dimensional dirty $p_x+i{p}_y$ superconductors as a function of disorder strength and chemical potential of the underlying normal state. We establish that the power law quasilocalized states induced by rare statistical fluctuations of the disorder potential give rise to an exponentially small density of states at zero energy, and even infinitesimally weak disorder converts the ThSM into a thermal diffusive Hall metal (ThDM). Consequently, the phase diagram of the disordered model only consists of ThDM and thermal insulating phases. We show the existence of two types of thermal insulators: (i) a trivial thermal band insulator (ThBI) [or BEC phase] with a smeared gap that can occur for suitable band parameters and all strengths of disorder, supporting only exponentially localized Lifshitz states (at low energy) and (ii) a thermal Anderson insulator that only exists for large disorder strengths compared to all band parameters. We determine the nature of the two distinct localization-delocalization transitions between these two types of insulators and ThDM. Additionally, we establish the scaling properties of an avoided (or hidden) quantum critical point for moderate disorder strengths, which govern the crossover between ThSM and ThDM phases over a wide range of energy scales. We also discuss the experimental relevance of our findings for three-dimensional, time reversal symmetry breaking, triplet superconducting states.

Figure 1

(a) Phase diagram of the clean, lattice model for a three-dimensional $p_x+i{p}_y$ superconductor [see Eq. (1)], as a function of the ratio $\mu /t$, where $\mu$ and $t$ are, respectively, the chemical potential and the nearest neighbor hopping strength for the underlying normal state. Thermal Hall semimetals with two and four Weyl nodes are, respectively, denoted as ThSM2 and ThSM4, and ThBIs represent fully gapped, thermal insulators arising due to BEC/molecular pairing in the absence of an underlying Fermi surface. (b) Berry curvature induced low temperature thermal Hall conductivity ${\kappa }_{xy}$ for the clean model. The maximum value of ${\kappa }_{xy}$ is $k_B^{2}T/\left(6ha\right)$, where $a$ is the lattice spacing. This maximum value occurs when two Weyl nodes get shifted to the Brillouin zone boundary. Notice that one chiral, surface Majorana mode contributes $k_B^{2}T/\left(6h\right)$ to ${\kappa }_{xy}$, and the displayed bulk ${\kappa }_{xy}$ equals the net contribution from all chiral modes, reflecting the bulk-boundary correspondence. (c) Zero energy density of states $\rho \left(0\right)$ for a clean ($W=0$) system with linear dimension $L=60$ (in units of the lattice spacing). We can clearly discern each phase even though in the thermodynamic limit ($L\to \infty$), $\rho \left(0\right)=0$ in all of these phases. The $\delta$ functions in the DOS are broadened by a Gaussian of width $\sigma =D\pi /2^{11}$ for a bandwidth $D$ (Ref. [59]).

Figure 2

(a) The zero energy (or temperature) schematic phase diagram of Hamiltonian Eq. (1) with $\mu$ as the chemical potential for the electrons, $W$ is the strength of disorder, and $t$ is the hopping strength. ThDM is a diffusive thermal Hall metal with finite density of states at zero energy, ThBI is the thermal band insulator supporting exponentially localized Lifshitz states, and AI denotes an Anderson insulator phase of class D, BdG quasiparticles. The thin solid orange line indicates the crossover from thermal Hall semimetals (ThSMs) and ThDM which occur at finite energies and are governed by an avoided QCP [as shown in (b)]. The dashed line at large $W/t$ and $\mu /t$ represents the possible separation between the two Anderson insulating phases. We expect that a ThDM regime (of a width on the order of the hopping strength) could exist along this line due to the different nature of the two localized phases and their respective transitions. (b) A sketch of what each regime means [taken as a cut in (a) for either fixed $\mu$ or $W$ such that one passes through both the ThBI to ThDM localization transition and the avoided transition], which we verify in each relevant section below. Due to the exponentially small DOS at $E=0$ there is only a ThDM at low energies; despite this the DOS still resembles that of a ThSM above a crossover energy scale. The perturbatively predicted critical point between ThSM and ThDM becomes avoided, but its quantum critical fan can be probed numerically over a wide range of energies. The ThBI phase is insulating and an actual critical point separates it from the ThDM phase with its own critical fan. Nonetheless, $\rho \left(0\right)>0$ in the ThBI phase purely from localized rare states (represented by the grayed region in the ThBI phase). $\mathrm{\Lambda }$ is some characteristic energy scale, above which short distance or lattice effects are important. This figure does not depict the Anderson localization transition occurring at a very large disorder strength (or energy).

Figure 3

The stochastic trace evolution of the density of states produces an artificial background as seen by the use of “unnormalized” vectors. However, if one just uses a small number of “normalized vectors” much more accurate results at low values of $\rho \left(E\right)$ can be found. The calculation in this plot was done with 100 disorder realizations and Gaussian disorder at $\mu =0$. This effect persists at lower system sizes where it is confirmed by computing the exact trace.

Figure 4

Deep within the insulating phase, we fit the Lifshitz Tail due to rare region effects: $\rho \left(E\right)-\rho \left(0\right)\approx A{e}^{b|E-E_0{|}^{-3/2}}$ for box disorder (with $A\approx 5.31\times {10}^{-12}$, $b\approx 0.209$, and $E_0\approx 0.233t$). The thin individual lines are averages over 1000 realizations while the thick (red) line is an average over 28,000 realizations—all calculations are done with $N_c=2^{13}$. The thin lines clearly show rare contributions to the DOS deep in the average band gap. $E_0$ is a good approximation for the average band edge, drawn as a vertical line in both plot and inset. (Inset) Same data on a linear-linear scale.

Figure 5

The finite size DOS at weak disorder. (a) Depending on the physics we are looking at, it is useful to sometimes keep one Weyl peak at zero energy as shown here. (b) The DOS for various $L$ with Gaussian disorder as a function of $LE$. At weak disorder the peaks are spaced like $1/L$ and therefore are composed of perturbatively dressed Weyl states that are smoothly connected to their $W=0$ counterparts. There are no states at zero energy; the finite value of the DOS at $E=0$ is solely due to the overlap of the two Weyl peaks at $EL/t\approx 5$ broadened by the KPM. (c) The DOS versus $E$ computed with a gaussian disorder distribution displaying the crossover regimes at low $E$. At very low energies the density of states is nonzero and essentially $E$ independent displaying the ThDM regime. At finite $E$, the ThSM regime is intact with characteristic $\rho \left(E\right)\sim E^2$ dispersion, and at large $E$, we start to see the critical fan $\rho \left(E\right)\sim \left|E\right|$.

Figure 6

(a) Projected probability density ${\sum }_z{\left|\psi (x,y,z)\right|}^2$ in the $xy$ -plane of a rare wave function that corresponds to a state in the low energy tail of the DOS with $W/t=0.8$, $L=18$, and periodic boundary conditions. (b) A scatter plot of the wave function decay away from its maximum value on a log-log plot, showing a clear power-law trend. (c) The scatter plot data is binned and is fit to $a/r^x$ showing a clear power law decay at small $r$ with $x=1.6$ and the rare wave function is indeed quasilocalized.

Figure 7

(a) The DOS at zero energy $\rho \left(0\right)$ as a function of disorder strength $W$. A Weyl state at zero energy produces a very large finite size effect in $\rho \left(0\right)$, whereas the DOS without a Weyl state becomes $L$ independent at weak disorder ($W\gtrsim 0.75t$). With a Weyl state at zero energy, we can pinpoint the crossover due to the avoided transition to occur roughly around $W/t\approx 1.0$ (for Gaussian disorder by finite size scaling). (b) Without a Weyl state at $E=0$, rare states begin populating the low-$\left|E\right|$ DOS, and they appear as an $L$ independent DOS. The data is well fit to the rare region form $\rho \left(0\right)\sim exp(-a/W^2)$, and our results are consistent with only the ThDM phase persisting at zero energy for weak disorder. Even without an $E=0$ state, we hit an artificial background DOS due to the Gaussian broadening of the Dirac-delta function in the DOS of nearby states around $W/t\approx 0.7$.

Figure 8

(a) The second derivative of the DOS at zero energy at $\mu =0$ and as a function of disorder strength $W/t$ for box disorder. The peak describes the avoided QCP and remains finite as we saturate the value of $L$ at each particular $N_C$. (b) Focusing at $W/t=3.4$, and increasing with $N_C$, we eventually see the value of the peak saturating suggesting the critical point is indeed avoided. This data has been averaged over the twisted boundary conditions.

Figure 9

(a) The second derivative of the DOS across the transition from ThSM to ThBI. At $W=0$, this peak is nonanalytic, but for finite disorder, we find that it saturates as we increase $N_C$. We see this saturating clearly in (b) where the peak's height is independent of the size for $N_C=2^{10}$. (c) The second derivative of the DOS as a function of $\mu /t$ across the transition from ThDM to ThBI. We saturate the second derivative of the DOS across the transition. The results are suggestive that the DOS remains analytic across the transition but are inconclusive as to whether an intermediate ThSM regime remains between the ThBI phase and ThDM regime. Each of these was computed with 1000 disorder realizations.

Figure 10

Phase diagram for axial chemical potential disorder in this system with $\mathrm{\Delta }=t$. The respective phases are “metallic” or “insulating” in terms of the Bogoliubov quasiparticles (leading to thermal insulators or thermal conductors), but thermally metallic phases host different regimes that can appear semimetallic at nonzero energy or diffusive at low energies with avoided quantum critical lines that “separate” the two as we plot here in dashed orange. $\mathrm{ThSM}n$ represents a thermal semimetallic regime with $n$ nodes. ThBI is the thermal band insulating phase while ThDM is the thermal diffusive metal regime. The dashed line in the ThDM phase is the onset of a zero-energy antilocalization peak in the DOS associated with the class D system [52]. The color plot is the value of the order parameter $\rho \left(0\right)$ and is interpolated from all of the data computed in this paper. The noisy nature of $\rho \left(0\right)$ in the semimetallic phases is due to finite size effects and sparsity of data. The color scale uses a lower bound of $\rho \left(0\right)\ge 3\times {10}^{-4}$.

Figure 11

The crossover from ThSM4 to ThDM at $\mu =0$ for box disorder. (a) The rise of $\rho \left(0\right)$ as a function of disorder on a linear-linear scale. In the inset, we see the saturation of $\rho \left(0\right)$ in terms in the ThDM regime for large $L$. (b) The DOS changes across this transition, going through a quantum critical region where the DOS is linear near $E=0$. Notice that a peak begins to develop for the higher disorder. We investigate this peak in Appendix pp2.

Figure 12

(a) The crossover from ThSM4 (left of peak) to ThSM2 (right of peak) by observing the finite-size induced peak for $\rho \left(0\right)$ as a function of $\mu$. Notice that away from the peak the data becomes noisy. This is the expected behavior for $\rho \left(0\right)$ in the semimetallic regime at finite size [see Fig. 1]. This data is taken at $L=60$. An estimate for the full-width at half-maximum is used for the error. (b) The avoided multicritical point can be visually captured here by considering $\rho \left(0\right)$ as a function of $\mu /t$. The peak that began at $\mu =2t$ for $W=0$ smoothly transitions to the peak seen here around $\mu \sim 2.4t$. On the right side, we have either ThDM or ThSM4 depending on the value of $W/t$. On the right side, it quickly becomes ThSM2.

Figure 13

The transition from ThDM to ThBI. (a) $\rho \left(0\right)$ as it tends to zero from the ThDM to the ThBI (left) and the average gap that forms in the insulator (right). The dashed line represents a power-law fit of the gap for $L=120$. (b) The density of states across the transition. Notice that roughly when the gap closes, $\rho \left(0\right)>0$ indicating the ThDM phase. (c) The transition from the ThSM2 regime (left of peak) to the ThBI phase (right of peak) up to disorder $W/t=4$. As with the ThSM4 to ThSM2 crossover in Fig. 12, the peak at the anisotropic Weyl node is expected and the subsequent fall to $\rho \left(0\right)=0$ even at this finite size of $L=60$ is expected in the ThBI phase. The transition and error bars are read off the plot by the position of the peak and resolution of the data.

Figure 14

When fitting $\rho \left(E\right)\sim E^{d/z^{*}-1}$ across the phase transition from ThSM2 to ThBI, we can observe a sudden increase in the power $d/z^{*}-1$ due to entering the ThBI phase. This can be used to numerically give an upper bound to the transition ${\mu }_I\left(W\right)$.

Figure 15

In this figure, we determine the critical exponents for the quantum critical region between the ThSM4 to ThDM regimes. The exponent $\nu$ is determined by two methods: (a) by $\rho \left(0\right)\propto {\delta }^{\nu (d-z)}$ where $\delta =(W-W_c)/W_c$ and (b) by data collapse minimizing the collapse function $S\left(\nu \right)$. This is illustrated here for $W_c=3.525$. The exponent $z$ is found simply in (c) by fitting $\rho \left(E\right)\propto {(E-b)}^{\frac{3}{z}-1}$. The dashed lines are numerical fits. The computed values for $\nu$ and $z$ are quoted in Table 2.

Figure 16

The critical exponents $z$ (a) and $\nu$ (b) along the critical line connecting the ThSM phases to the ThDM. It becomes more difficult to pin down a value of $z$ the closer we get to the ThBI phase. In fact, the uncertainty in $z$ is related to our inability to pin down the phase transition to a good accuracy in those regions. Similarly, it becomes difficult to pin down $\nu$ (when $\mu >8t$, $\nu$ cannot even be determined reliably due to proximity to the transition to ThBI). The $\nu$ data was found by finding the inflection point of $log\rho \left(0\right)$ vs $log\delta$ by cubic interpolation and determining the tangent at that point. Most of this data was found with $\rho \left(E\right)$ computed with $L=120$ (for $\mu =0$, 2.32, 2.42, 2.52, and 5, $L=140$ data is used).

Figure 17

The critical exponents on the critical line ${\mu }_I\left(W\right)$. (a) The critical exponent $z^{*}$ extracted from fitting $\rho \left(E\right)\sim E^{d/z^{*}-1}$. In the ThBI phase $z^{*}\to 0^{+}$. For $W=0$ (the blue data point), this value is known exactly from analytics. For $1\le W\le 4$ (the orange data points), $\rho \left(E\right)$ is noisy, so we find this exponent for larger bounds that could see higher band effects. For $W\ge 5$, we can fit $\rho \left(E\right)$ with $E/t\ll 1$ and these are the exponents we find. (b) The critical exponent $\gamma =\nu z$, that describes how the band insulating gap increases in the ThBI phase: ${\mathrm{\Delta }}_g\sim {(\mu -{\mu }_c)}^{\gamma }$. Systematics that could lead to appreciable error in this quantity cannot be reliably estimated.

Figure 18

The fitting of the DOS at the transition from the ThSM2 regime to the ThBI phase. In (a) we have a smaller uncertainty for the transition, and in (b) we can fit a power law within the entire region to get a smooth change from $z^{*}\approx 1.6$ to $z^{*}\approx 0.4$ within this window.

Figure 19

The QCP associated with the ThBI to ThDM transition has a critical fan that can be observed here in the DOS. (a) Moving from the ThBI to the critical point (represented by red), we see the fan sharply as the gap closes. (b) On the ThSM side of the transition, the fan is a bit more subtle since the critical fan is a slight change in exponent $\rho \left(E\right)\sim E^2$ (or $z=1$) for the ThSM and $\rho \left(E\right)\sim E^{1.12}$ (or $z=1.4$) for the other (see the inset). The critical fan appears at higher energies the farther from the transition we go. (c) At higher disorder, the critical fan from the avoided QCP swamps all effects and we begin seeing that quantum critical fan appear. For all plots the red line represents the approximate critical point and we have taken $L=60$ and box potential disorder for these plots.

Figure 20

Plots of the mobility edge defined by $M_{\mathrm{edge}}\left(E\right)={\rho }_t\left(E\right)/\rho \left(E\right)$ tuning disorder $W$ while keeping $\mu =0$ constant. We see a large antilocalization peak in the DOS develop in addition to the standard behavior that states far from $E=0$ localize first. The function $\rho \left(E\right)$ is calculated at $L=60$ and ${\rho }_t\left(E\right)$ is calculated at $L=30$ and $N_C=2^{13}$. In (a) we have the whole range given by a color plot, and in (b) we see some cuts illustrating how drastic the peak is.

Figure 21

Typical DOS displaying the Anderson localization transition at large disorder for $\mu =0$ and $L=30$. (a) The $N_C$ dependence of the typical DOS. In the localized phase, ${\rho }_t\left(0\right)$ decreases as $N_C$ is increased whereas in the ThDM phase ${\rho }_t\left(0\right)$ is relatively insensitive to $N_C$. While this is plotted for $\mu =0$, the plot is nearly identical for any $\mu$ [see Fig. 22]. (b) The typical DOS goes to zero while the average DOS remains finite. The average DOS is computed with $L=60$. (c) The peak in the typical DOS before the transition to the AI. One sees little or even positive $N_C$ dependence at $E=0$. (d) The peak splits and now ${\rho }_t\left(E\right)$ has a strong $N_C$ dependence. With the theory for the ThDM breaking down and other measures, we conclude this is the localized phase.

Figure 22

(a) The Anderson insulator transition line at large disorder extracted from the $N_C$ dependence of the typical DOS. (b) The fit of the power law for the typical DOS across the Anderson localization transition in the regime that is $N_C$ independent. We find $\beta =2.9\pm 0.1$ and the fit is the red line. All data here is taken at $L=30$.

Figure 23

Renormalized disorder couplings vs RG flow time $l=log({\mathrm{\Lambda }}_0/\mathrm{\Lambda })$. The disorder couplings ${\mathrm{\Delta }}_1$, ${\mathrm{\Delta }}_2$, ${\mathrm{\Delta }}_3$, and ${\mathrm{\Delta }}_4$ are, respectively, shown as the green, orange, blue, and purple lines. (a) for the initial values ${\mathrm{\Delta }}_1\left(0\right)={\mathrm{\Delta }}_2\left(0\right)=0.8701$, ${\mathrm{\Delta }}_3\left(0\right)={\mathrm{\Delta }}_4\left(0\right)=0$ the renormalized disorder couplings flow to zero at long wavelength limit, signifying a perturbatively stable thermal semimetal phase. (b) For the initial values ${\mathrm{\Delta }}_1\left(0\right)={\mathrm{\Delta }}_2\left(0\right)=0.8702$, ${\mathrm{\Delta }}_3\left(0\right)={\mathrm{\Delta }}_4\left(0\right)=0$, the renormalized axial chemical potential coupling ${\mathrm{\Delta }}_4\left(l\right)$ diverges at long wavelength limit, indicating a diffusive metal phase. Therefore the universality class of quantum phase transition between semimetal and diffusive metal phases is described by the axial chemical potential disorder controlled quantum critical point, even though we are explicitly tuning the strength of random chemical potential for normal quasiparticles.

Figure 24

The peak in both the normal DOS (a) and typical DOS (b) around $E=0$ as predicted by perturbation theory for a ThDM. The peak becomes well defined at finite disorder and persists up to all disorders considered here. The DOS is calculated at $L=60$ and the typical DOS is calculated at $L=30$. (c) We show here two example fits to the expected analytical form for the antilocalization peak in $\rho \left(E\right)$ to the numerical data defined with the integral in Eq. (B1) for small $E$. (d) The fit (the dashed red line) works surprisingly well to even lower energies for ${\rho }_t\left(E\right)$.