Disorder-driven itinerant quantum criticality of three-dimensional massless Dirac fermions

Progress in the understanding of quantum critical properties of itinerant electrons has been hindered by the lack of effective models which are amenable to controlled analytical and numerically exact calculations. Here we establish that the disorder-driven semimetal to metal quantum phase transition of three-dimensional massless Dirac fermions could serve as a paradigmatic toy model for studying itinerant quantum criticality, which is solved in this work by exact numerical and approximate field-theoretic calculations. As a result, we establish the robust existence of a non-Gaussian universality class, and also construct the relevant low-energy effective field theory that could guide the understanding of quantum critical scaling for many strange metals. Using the kernel polynomial method (KPM), we provide numerical results for the calculated dynamical exponent ($z$) and correlation length exponent ($\nu$) for the disorder-driven semimetal (SM) to diffusive metal (DM) quantum phase transition at the Dirac point for several types of disorder, establishing its universal nature and obtaining the numerical scaling functions in agreement with our field-theoretical analysis.

Figure 1

Numerically calculated finite temperature $\left(T\right)$ phase diagram as a function of disorder $\left(W\right)$ computed for a linear system size $L=110$ (with $t$, the band hopping energy, setting the energy scale). The Dirac semimetal (SM), the diffusive metal (DM), and the QCP separating them only exist at zero temperature and there is no finite temperature phase transition. The squares mark the crossover out of the SM phase (green) into the quantum critical (QC) region (orange), which is anchored by the QCP at the critical disorder value $W_c/t=2.55\pm 0.05$ (star). This crossover is determined by where the specific heat fails to be described by $C_V\sim T^3$. The circles mark the crossover out of the DM phase (blue) into the QC region, which is marked by where the specific heat fails to be described by $C_V\sim T$. The dashed lines connecting the squares/circles are a guide to the eye. The dashed line at $\mathrm{\Lambda }/t$ marks the high-energy cutoff above which the lowenergy continuum description in terms of massless Dirac fermions is inapplicable. In a narrow region inside the quantum critical regime the specific heat manifests non-Fermi-liquid behavior $C_V\sim T^2$ which gradually crosses over to the power-law behaviors of SM and DM. The dashed horizontal line $E_F/t$ is a schematic for the (generally unknown) Fermi energy associated with any residual doping that drives the system slightly away from the precise Dirac point, indicating the low-energy cutoff above which the zero-doping Dirac point theory is applicable. Note that as long as $E_F$ is not too large, quantum criticality manifests itself in the orange region of the phase diagram, even though the QCP itself (the star at $W_c$ and $T=0$) is “hidden” by the residual doping. At temperature and energy scales smaller than $E_F$, the system behaves as a conventional metal or diffusive Fermi liquid. In a nominally undoped system, the inequality $\mathrm{\Lambda }\gg E_F$ applies, and quantum criticality is observable for $k_{B}T>E_F$.

Figure 2

Numerically (KPM) calculated DOS as a function of disorder strength for a linear system size $L=110$. (a) DOS at zero energy $\rho \left(0\right)$ as a function of disorder strength for each type of disorder matrix $A_W$, which shows that axial, mass, and potential disorder yield identical results. Inset: Corresponding calculation for axial disorder in two dimensions [i.e., setting the hopping in the $z$ direction to zero in Eq. (1)] clearly establishing that there is no stable SM phase and $\rho \left(0\right)$ is finite even for an infinitesimal amount of disorder. (b) DOS for axial disorder as a function of energy passing through the QCP with the arrow indicating increasing disorder strength. Specific heat (solid lines) with axial disorder for $W<W_c$ (c) and $W>W_c$ (d), the fits (dashed lines) are to the leading $C_V\sim T^3$ in the SM and $C_V\sim T$ in the DM. The temperatures at which the fits no longer apply determine the crossover boundary shown in Fig. 1.

Figure 3

Determination of the critical exponent $z$ for potential disorder with the box distribution and $L=130$. (a) DOS as a function of energy, and (b) the specific heat as a function of temperature, both obtained at the QCP, $W_c=2.55\pm 0.05$.

Figure 4

Determination of the critical exponent $\nu$ for potential disorder with the box distribution [(a) and (c)] and the Gaussian distribution [(b) and (d)]. Panels (a) and (b) show the numerical DOS as a function of the distance to the QCP. The dashed lines are power-law fits to extract $\nu$, the thick (thin) dashed lines are fits closer to (further from) the QCP. Panels (c) and (d) show the numerical finite-size scaling and data collapse of $\rho (0,L)$; the labels for $L$ are shared between (c) and (d). Insets: The local-linearity function [37], which yields the best scaling collapse at its minimum value. Here $\delta =|W-W_c|/W_c$ is the usual dimensionless tuning parameter defining the QPT.

Figure 5

Scaling collapse of the DOS and the specific heat in the quantum critical fan ${\delta }^{\nu z}\ll E/t,T/t\ll \mathrm{\Lambda }/t$ (with $\delta =|W-W_c|/W_c$) establishing that both exhibit single parameter scaling as a function of the distance to the QCP given by Eqs. (7) and (8). The dashed line is the approximate crossover functions analytically determined from the one-loop $\beta$ functions (see Appendix pp1); (a) and (c) are for $W<W_c$, and (b) and (d) are for $W>W_c$. We stress that the dashed lines are not a fit, and only the nonuniversal amplitudes are adjusted. Despite the significant difference between the numerically obtained value for $\nu$ and the one-loop result ($\nu =1$), we find considerable agreement between the analytical and numerical crossover functions. We are using $z=1.46$ and the value of $\nu$ from the finite-size scaling for the collapse of the numerical data ${\nu }_L\approx 1.46$. These results are for the axial disorder using the box distribution, but the results for other distributions and disorder are very similar.

Figure 6

(a) Finite-size scaling of $\rho (0,L,N_c)$ on a log-log scale, for $N_c=1028$ in the clean limit ($W=0$) yields a power-law decay in system size $1/L^{2.85}$ consistent with $1/L^3$. (b) Scaling of $\rho (0,L,N_c)$ on a log-log scale with $N_c$ for $L=60$ yielding a linear in $N_c$ scaling. Here $N_c$ is the number of terms kept in the Chebyshev expansion in the KPM numerics.

Figure 7

DOS for $L=60$ and $W=2.0t<W_c$ (a) $N_c$ dependence of $\rho (0,L,N_c)$ converges in $N_c$ for $N_c\ge 4096$. (b) Low-energy dependence of $\rho (E,L,N_c)$ for various $N_c$ showing that the deviations in the DOS for larger $N_c$ in the SM phase is a result of the KPM beginning to resolve individual eigenstates. The low-energy peak is a result of the eight Dirac states that are there in the clean limit that have now been broadened by disorder.

Figure 8

DOS for $L=60$ and $W=3.0t>W_c$; (a) $N_c$ dependence of $\rho (0,L,N_c)$ converges in $N_c$ for $N_c\ge 4096$. (b) Low-energy dependence of $\rho (E,L,N_c)$ for various $N_c$ showing that there is a very weak dependence in the DM phase.

Figure 9

DOS at zero energy $\rho \left(0\right)$ as a function of of the inverse system size $1/L$ on a log-log scale in both the SM and DM phases for potential disorder up to $L=110$. Results for the box distribution (a) and the Gaussian distribution (b). In the SM phase $\rho \left(0\right)$ goes to zero as the system size is increased, while it saturates in the DM phase, and at the QCP $\rho \left(0\right)$ is still decreasing for increasing $L$. In the DM phase the results are saturating in $1/L$.