Quantitative analytical theory for disordered nodal points [Article and Erratum]

TL;DR

This paper develops an analytical approach using the functional renormalization group to accurately evaluate the average density of states near disordered nodal points in materials like graphene and Weyl semimetals, matching numerical data.

Contribution

It introduces a novel analytical method for disordered nodal points that surpasses previous approaches in accuracy for both 2D and 3D systems.

Findings

Excellent agreement with numerical simulations in 2D

Significant improvement over existing methods in 3D

Provides a quantitative analytical framework for disorder effects

Abstract

Disorder effects are especially pronounced around nodal points in linearly dispersing bandstructures as present in graphene or Weyl semimetals. Despite the enormous experimental and numerical progress, even a simple quantity like the average density of states cannot be assessed quantitatively by analytical means. We demonstrate how this important problem can be solved employing the functional renormalization group method and, for the two dimensional case, demonstrate excellent agreement with reference data from numerical simulations based on tight-binding models. In three dimensions our analytic results also improve drastically on existing approaches.

Equations

$$H_d=H_{0,d}+U_d,$$

$${\mathcal{K}}_d\left(\mathbf{r}-{\mathbf{r}}^{\prime }\right)=⟨U_d(\mathbf{r})U_d({\mathbf{r}}^{\prime })⟩=K\frac{(\hslash v{)}^2}{(2\pi {)}^{d/2}{\xi }^2}{e}^{-|\mathbf{r}-{\mathbf{r}}^{\prime }{|}^2/2{\xi }^2},$$

$${\mathcal{K}}_d^{GWN}\left(\mathbf{r}\right)=K(\hslash v{)}^2{\xi }^{d-2}\delta \left(\mathbf{r}\right),$$

$$\nu \left(E\right)=-\frac{1}{\pi }\mathrm{Im}\mathrm{Tr}{\int }_{\mathbf{k}}{G}_{\mathbf{k}}^R\left(E\right),$$

$$H_{0,d}^{\mathrm{L}}$$

$$S_{d,\mathrm{dis}}$$

$$\times$$

$$d\tilde{K}/d\ln \lambda ={\tilde{K}}^2/\pi +{\tilde{K}}^3/(2{\pi }^2).$$

$$\Sigma \left(\mathbf{k}\right)$$

$$+K^2{\left(\hslash v\right)}^4{\int }_{\mathbf{q},\mathbf{p}}{e}^{-\frac{1}{2}{\xi }^2({|\mathbf{k}-\mathbf{p}|}^2+{|\mathbf{q}-\mathbf{p}|}^2)}$$

$$\times G(\mathbf{p})\cdot G(\mathbf{q})\cdot \left[G(\mathbf{k}+\mathbf{q}-\mathbf{p})+G(\mathbf{p})\right],$$

$${\partial }_{\Lambda }{\Sigma }_{\Lambda }\left(1^{\prime };1\right)=-{\int }_{2,2^{\prime }}{\left[{\dot{G}}_{\Lambda }\right]}_{2,2^{\prime }}{\Gamma }_{\Lambda }\left(1^{\prime },2^{\prime };1,2\right)$$

$${\partial }_{\Lambda }{\Gamma }_{\Lambda }\left(1^{\prime },2^{\prime };1,2\right)$$

$${\Gamma }_{\infty }\left(1^{\prime },2^{\prime };1,2\right)$$

$$A_{\infty }\left({\mathbf{k}}_1^{\prime }{\sigma }_1^{\prime },{\mathbf{k}}_2^{\prime }{\sigma }_2^{\prime };{\mathbf{k}}_1{\sigma }_1,{\mathbf{k}}_2{\sigma }_2\right)\equiv -2\pi {\delta }_{{\mathbf{k}}_1^{\prime }+{\mathbf{k}}_2^{\prime }-{\mathbf{k}}_1-{\mathbf{k}}_2}K{\left(\hslash v\right)}^2{e}^{-\frac{1}{2}{\xi }^2{|{\mathbf{k}}_1^{\prime }-{\mathbf{k}}_1|}^2}{\delta }_{{\sigma }_1^{\prime }{\sigma }_1}{\delta }_{{\sigma }_2^{\prime }{\sigma }_2},$$

$${\Gamma }_{\Lambda }\left(1_{{\alpha }_1^{\prime }i{\omega }_1^{\prime }}^{\prime },2_{{\alpha }_2^{\prime }i{\omega }_2^{\prime }}^{\prime };1_{{\alpha }_{1}i{\omega }_1},2_{{\alpha }_{2}i{\omega }_2}\right)$$

$${\partial }_{\Lambda }{\Gamma }_{\Lambda }\left(1_{{\alpha }_{1}i{\omega }_1}^{\prime },2_{{\alpha }_{2}i{\omega }_2}^{\prime };1_{{\alpha }_{1}i{\omega }_1},2_{{\alpha }_{2}i{\omega }_2}\right)={\int }_{3,3^{\prime },4,4^{\prime }}$$

$${\partial }_{\Lambda }{\Sigma }_{\Lambda }{\left({\alpha }_{1}i{\omega }_1{\mathbf{k}}_1\right)}_{{\sigma }_1^{\prime },{\sigma }_1}$$

$${\partial }_{\Lambda }{\Sigma }_{\Lambda }{\left(\alpha i\omega {\mathbf{k}}_1\right)}_{{\sigma }_1^{\prime },{\sigma }_1}=\sum _{{\sigma }_2,{\sigma }_2^{\prime }}{\int }_{{\mathbf{k}}_2}{\left[{\dot{G}}_{\Lambda }\left(\alpha i\omega {\mathbf{k}}_2\right)\right]}_{{\sigma }_2,{\sigma }_2^{\prime }}{A}_{\Lambda }^{\alpha i\omega ,\alpha i\omega }\left({\mathbf{k}}_1{\sigma }_1^{\prime },{\mathbf{k}}_2{\sigma }_2^{\prime };{\mathbf{k}}_2{\sigma }_2,{\mathbf{k}}_1{\sigma }_1\right).$$

$${\partial }_{\Lambda }{A}_{\Lambda }^{\alpha i\omega ,\alpha i\omega }\left({\mathbf{k}}_1^{\prime }{\sigma }_1^{\prime },{\mathbf{k}}_2^{\prime }{\sigma }_2^{\prime };{\mathbf{k}}_1{\sigma }_1,{\mathbf{k}}_2{\sigma }_2\right)={\int }_{{\mathbf{k}}_3,{\mathbf{k}}_4}\sum _{{\sigma }_3^{\prime }{\sigma }_3,{\sigma }_4^{\prime }{\sigma }_4}$$

$${\partial }_{\Lambda }{\Sigma }_{\Lambda }{\left(i\omega {\mathbf{k}}_1\right)}_{{\sigma }_1^{\prime },{\sigma }_1}=\delta \left(|i\omega |-\Lambda \right)\sum _{{\sigma }_2,{\sigma }_2^{\prime }}{\int }_{{\mathbf{k}}_2}{\tilde{G}}_{\Lambda }{\left(i\omega ,{\mathbf{k}}_2\right)}_{{\sigma }_2,{\sigma }_2^{\prime }}{A}_{\Lambda }^{i\omega ,i\omega }\left({\mathbf{k}}_1{\sigma }_1^{\prime },{\mathbf{k}}_2{\sigma }_2^{\prime };{\mathbf{k}}_2{\sigma }_2,{\mathbf{k}}_1{\sigma }_1\right),$$

$$\Sigma {\left(i\omega {\mathbf{k}}_1\right)}_{{\sigma }_1^{\prime },{\sigma }_1}=-\sum _{{\sigma }_2,{\sigma }_2^{\prime }}{\int }_{{\mathbf{k}}_2}G{\left(i\omega {\mathbf{k}}_2\right)}_{{\sigma }_2,{\sigma }_2^{\prime }}{A}^{i\omega ,i\omega }\left({\mathbf{k}}_1{\sigma }_1^{\prime },{\mathbf{k}}_2{\sigma }_2^{\prime };{\mathbf{k}}_2{\sigma }_2,{\mathbf{k}}_1{\sigma }_1\right),$$

$$A^{i\omega ,i\omega }\left({\mathbf{k}}_1^{\prime }{\sigma }_1^{\prime },{\mathbf{k}}_2^{\prime }{\sigma }_2^{\prime };{\mathbf{k}}_1{\sigma }_1,{\mathbf{k}}_2{\sigma }_2\right)=A_{\infty }\left({\mathbf{k}}_1^{\prime }{\sigma }_1^{\prime },{\mathbf{k}}_2^{\prime }{\sigma }_2^{\prime };{\mathbf{k}}_1{\sigma }_1,{\mathbf{k}}_2{\sigma }_2\right)-{\int }_{{\mathbf{k}}_3,{\mathbf{k}}_4}\sum _{{\sigma }_3^{\prime }{\sigma }_3,{\sigma }_4^{\prime }{\sigma }_4}$$

$$A^{i\omega ,i\omega }\left({\mathbf{k}}_1^{\prime }{\sigma }_1^{\prime },{\mathbf{k}}_2^{\prime }{\sigma }_2^{\prime };{\mathbf{k}}_3{\sigma }_3,{\mathbf{k}}_4{\sigma }_4\right)\left({\left[G\left(i\omega {\mathbf{k}}_3\right)\right]}_{{\sigma }_3,{\sigma }_3^{\prime }}{\left[G\left(i\omega {\mathbf{k}}_4\right)\right]}_{{\sigma }_4,{\sigma }_4^{\prime }}\right)A^{i\omega ,i\omega }\left({\mathbf{k}}_3{\sigma }_3^{\prime },{\mathbf{k}}_4{\sigma }_4^{\prime };{\mathbf{k}}_1{\sigma }_1,{\mathbf{k}}_2{\sigma }_2\right)$$

$$+A^{i\omega ,i\omega }\left({\mathbf{k}}_1^{\prime }{\sigma }_1^{\prime },{\mathbf{k}}_4{\sigma }_4^{\prime };{\mathbf{k}}_3{\sigma }_3,{\mathbf{k}}_2{\sigma }_2\right)\left({\left[G\left(i\omega {\mathbf{k}}_3\right)\right]}_{{\sigma }_3,{\sigma }_3^{\prime }}{\left[G\left(i\omega {\mathbf{k}}_4\right)\right]}_{{\sigma }_4,{\sigma }_4^{\prime }}\right)A^{i\omega ,i\omega }\left({\mathbf{k}}_3{\sigma }_3^{\prime },{\mathbf{k}}_2^{\prime }{\sigma }_2^{\prime };{\mathbf{k}}_1{\sigma }_1,{\mathbf{k}}_4{\sigma }_4\right).$$

$$A\left({\mathbf{k}}_1^{\prime }{\sigma }_1^{\prime },{\mathbf{k}}_2^{\prime }{\sigma }_2^{\prime };{\mathbf{k}}_1{\sigma }_1,{\mathbf{k}}_2{\sigma }_2\right)=A_{\infty }\left({\mathbf{k}}_1^{\prime }{\sigma }_1^{\prime },{\mathbf{k}}_2^{\prime }{\sigma }_2^{\prime };{\mathbf{k}}_1{\sigma }_1,{\mathbf{k}}_2{\sigma }_2\right)-K^2{\left(\hslash v\right)}^4{\int }_{{\mathbf{k}}_3}$$

$$A\left({\mathbf{k}}_1{\sigma }_1^{\prime },{\mathbf{k}}_2{\sigma }_2^{\prime };{\mathbf{k}}_2{\sigma }_2,{\mathbf{k}}_1{\sigma }_1\right)=-K{\left(\hslash v\right)}^2{e}^{-\frac{1}{2}{\xi }^2{|{\mathbf{k}}_2-{\mathbf{k}}_1|}^2}{\delta }_{{\sigma }_1^{\prime }{\sigma }_1}{\delta }_{{\sigma }_2^{\prime }{\sigma }_2}-K^2{\left(\hslash v\right)}^4{\int }_{{\mathbf{k}}_3}$$

$$\Sigma \left(\mathbf{k}\right)$$

$$+K^2{\left(\hslash v\right)}^4{\int }_{\mathbf{q},\mathbf{p}}{e}^{-\frac{1}{2}{\xi }^2({|\mathbf{k}-\mathbf{p}|}^2+{|\mathbf{q}-\mathbf{p}|}^2)}G(\mathbf{p})\cdot G(\mathbf{q})\cdot \left[G(\mathbf{k}+\mathbf{q}-\mathbf{p})+G(\mathbf{p})\right].$$

$$\frac{{\Sigma }_{d=2}\left(\mathbf{x}=\mathbf{k}\xi \right)}{\hslash v/\xi }=m_2\left(x\right)\left\{{\sigma }_x\cos \left[\varphi \right]+{\sigma }_y\sin \left[\varphi \right]\right\}+i{M}_2\left(x\right),$$

$$\frac{{\Sigma }_{d=3}\left(\mathbf{x}=\mathbf{k}\xi \right)}{\hslash v/\xi }=m_3\left(x\right)\left(\sin \left[\theta \right]\left\{{\sigma }_x\cos \left[\varphi \right]+{\sigma }_y\sin \left[\varphi \right]\right\}+{\sigma }_z\cos \left[\theta \right]\right)+i{M}_3\left(x\right).$$

Full text

I Introduction

Two dimensional graphene Novoselov et al. (2005) and three dimensional Weyl materials Bernevig (2015) are important examples of Dirac type semimetals. Their electronic structure features a nodal degeneracy point where two linearly dispersing Bloch bands meet. Due to the vanishing density of states (DOS), disorder effects can be expected to be particularly pronounced in these materials and have been actively studied, for reviews see Refs. Das Sarma et al. (2011); Syzranov and Radzihovsky (2016). Despite all this effort on the disorder problem for nodal points, analytical results, even for a quantity as simple as the DOS, are at best qualitatively correct but fail widely in their quantitative predictions, even for weak disorder. This is surprising insofar as exact answers can be obtained with ease from numerical simulations of non-interacting lattice Hamiltonians. The scope of this work is to show how tremendous progress on this long-standing problem can be achieved by employing a variant of the functional renormalization group (fRG).

We consider the minimal continuum model of a single disordered node in $d\!=\!2,3$ dimensions,

[TABLE]

where $H_{0,2}=\hbar v(\sigma_{x}k_{x}+\sigma_{y}k_{y})$ is a $d\!=\!2$ Dirac Hamiltonian and $H_{0,3}=\hbar v(\sigma_{x}k_{x}+\sigma_{y}k_{y}+\sigma_{z}k_{z})$ a $d\!=\!3$ Weyl Hamiltonian written with the standard Pauli matrices $\sigma_{i=x,y,z}$. The disorder potential $U_{d}(\mathbf{r})$, taken to be proportional to the unit matrix, is commonly assumed to have Gaussian correlations and zero mean. Explicitly, we assume a smooth form of the correlator

[TABLE]

where $\left\langle...\right\rangle$ denotes the disorder average. As $H_{0,d}$ is lacking any scale, the disorder correlation length $\xi$ serves as the fundamental scale in the problem. The dimensionless parameter $K$ measures the disorder strength. In the Brillouin zone of real materials, nodal points usually come in pairs. This is enforced by symmetry (graphene) or topology (Weyl). However, these pairs can have a sizable k-space separation $\Delta k$. If $\xi\Delta k\gg 1$ the intra-node scattering dominates over inter-node scattering and the model (1) is a reasonable low-energy approximation for realistic materials.

While Eq. (1) with the correlator (2) has the advantage that it can be easily approximated in tight-binding models if $\xi\gg a$ ($a$ being the lattice scale) another common choice for $\mathcal{K}_{d}$ more convenient for analytical calculations is the white noise limit $\xi\rightarrow 0$,

[TABLE]

along with the prescription that $1/\xi$ serves as an ultraviolet cutoff for the clean dispersion $H_{0,d}$. We will use the white noise approximation to make contact with known results.

The bulk DOS can be calculated as

[TABLE]

where $\int_{\mathbf{k}}=\left(2\pi\right)^{-d}\int d\mathbf{k}$ and $G_{\mathbf{k}}^{R}\left(E\right)$ is the retarded (matrix-valued) Green function. For the clean Hamiltonian $H_{0,d}$, one has $\nu_{0,d}\left(E\right)=|E|^{d-1}/(2\pi)^{d-1}\left(\hbar v\right)^{d}$, vanishing at the degeneracy point. If disorder is thought of as a local chemical potential creating carriers from conduction or valence bands, a finite $\nu_{d}\left(E\!=\!0\right)$ can be expected (since disorder is a self-averaging quantity, we omit $\left\langle...\right\rangle$). In the following, we distinguish between ’numerical’ approaches based on explicit generation of a large number of random disorder realizations $U_{d}(\mathbf{r})$ in Eq. (1) and ’analytical’ methods starting from Eq. (2). While the former are well established, up to now there is no known analytical method that could reproduce numerical results with reasonable accuracy, not even for small $K$.

The scope of this work is to show how this long-standing problem can be solved by a variant of the functional renormalization group (fRG) which allows to rewrite the disorder problem as an — in principle infinite — hierarchy of coupled self-consistency equations for vertex functions. We apply this technique to calculate the DOS $\nu_{d}\left(E=0\right)$ and find that even a simple truncation of the above hierarchy yields results in very good quantitative agreement with numerically exact data obtained from the kernel polynomial method at much higher computational costs. We acknowledge an earlier study by Katanin Katanin (2013) with similar objectives but a different variant of the fRG. However, our results go significantly beyond those of Ref. Katanin (2013), where only $d\!=\!2$ was investigated without comparison to numerically exact results.

II Exact numerical DOS

To gauge the quality of analytical approaches discussed in the remaining sections, let us start by obtaining numerically exact DOS data for the Dirac and Weyl systems with smooth disorder, described by Eqns. (1) and (2). We apply the kernel polynomial method (KPM) Weisse et al. (2006), a numerically efficient tool to approximate the DOS of large lattice Hamiltonians $H$ represented as sparse matrices. The DOS $\nu\left(E\right)$ as a function of energy $E$ is expanded in Chebyshev polynomials and the expansion coefficients $\mu^{(n)}$ are expressed as a trace over a polynomial in $H$. Using recursion properties of Chebyshev polynomials, the $\mu^{(n)}$ can be efficiently computed (up to order $N$) involving only sparse matrix-vector products and a statistical evaluation of the trace.

The clean nodal Hamiltonian $H_{0,d}$ is approximated as the low energy theory of the following tight-binding models on a square/cubic lattice (with constant $a$, size $L^{d}$)

[TABLE]

which feature four/eight nodal points for $d\!=\!2$ and $d\!=\!3$, respectively, with minimal mutual distance $\Delta k=\frac{\pi}{a}$. We apply periodic boundary conditions and add a correlated disorder potential as in Eq. (2). If our disordered lattice model would faithfully emulate the continuum Hamiltonian (1), the DOS at zero energy must be of the scaling form $\nu_{d}\left(E=0\right)=\left(\hbar v\right)^{-1}\xi^{1-d}f\left(K\right)$ with $f(K)$ a dimensionless function. We have checked that the KPM data based on the lattice Hamiltonian Eq. (5) fulfills this scaling condition once $\xi\gg a$ so that (i) the smooth disorder correlations are well represented on the discrete lattice, (ii) the disorder induced energy scale is well below the scale of order $\hbar v/a$ where $H_{0,d}^{\mathrm{L}}$ deviates from $H_{0,d}$ and (iii) the inter-node scattering rate is sufficiently suppressed compared to the intra-node rate (the factor is $\exp[-(\Delta k)^{2}\xi^{2}/2]$). Moreover, we require $L\gg\xi$ to suppress finite-size effects. Thus, the KPM data (normalized to a single node) shown as dots in Fig. 1 ($d=2$) and Fig. 2 ($d=3$) can be regarded as the exact zero energy DOS of the continuum model Eq. (1). Simulation parameters are given in the figure captions. In spite of the abundant literature on similar numerical studies for the DOS of disordered 2d Dirac (see Refs. Peres et al. (2006); Pereira et al. (2008); Wu et al. (2008); Ziegler et al. (2009)) and 3d Weyl systems (see Refs. Kobayashi et al. (2014); Pixley et al. (2015); Trescher et al. (2017)), we are not aware of existing high-precision data obtained for a smooth disorder correlator and with the required scaling properties fulfilled.

III Disordered $d\!=\!2$ Dirac node

We proceed by discussing existing analytical approaches to the disorder problem in the $d\!=\!2$ Dirac case. The self-consistent Born approximation (SCBA) determines the disorder induced self-energy $\Sigma\equiv G^{-1}-G_{0}^{-1}$ (where $G_{0}$ is the Green function of the clean system) according to the diagram in Fig. 3(i) Shon and Ando (1998); Ostrovsky et al. (2006); Noro et al. (2010). The corresponding self-consistent equation can be solved in closed form for the white noise correlator Eq. (3) and yields a disorder induced scale $\Gamma=\frac{\hbar v}{\xi}e^{-2\pi/K}$ (for $K\lesssim 1$) exponentially small in $K$ appearing in the imaginary self-energy $\Sigma=\pm i\Gamma$ and a DOS $\nu_{2}\left(E=0\right)\hbar v\xi\propto e^{-2\pi/K}/K$ Ostrovsky et al. (2006). In Fig. 1 (bottom panel), this result (dashed line) compares well to the DOS obtained from the SCBA with smooth disorder correlator (2) (blue line). However, comparing to the exact KPM-DOS in Fig. 1 (dots), we find that albeit the exponential form is correctly predicted by the SCBA, the slope (prefactor in the exponent) is roughly a factor 2 off.

The failure of the SCBA can be attributed to interference corrections from multiple disorder scattering events Aleiner and Efetov (2006), see diagrams (ii.a) and (ii.b) for the lowest order corrections. While unimportant in ordinary metals (where $1/k_{F}l\ll 1$ with $k_{F}$ Fermi wavevector and $l$ the mean free path serves as a small parameter), for Dirac materials these diagrams provide corrections of order $\ln\left[\hbar v/\xi\Gamma\right]$. Accordingly, their contribution vanishes for strong disorder where the SCBA becomes reliable, c.f. Fig. 1.

To go beyond the SCBA, Refs. Ostrovsky et al. (2006); Aleiner and Efetov (2006) used the super-symmetry method. Alternatively, the replica trick Altland and Simons (2006) can be employed: It takes a disorder average over $R$ copies (replicas) of the original problem seeing the same disorder potential. The resulting action $S=S_{d,0}+S_{d,dis}$ is translational invariant but contains, besides the free part $S_{d,0}=\sum_{\alpha=1}^{R}\int_{\omega}\int_{\mathbf{k}}\sum_{\sigma,\sigma^{\prime}}\bar{\psi}_{\omega\mathbf{k}\sigma^{\prime}}^{\alpha}(i\omega-H_{0,d})_{\sigma^{\prime}\sigma}\psi_{\omega\mathbf{k}\sigma}^{\alpha}$ an attractive inter-replica interaction which is elastic (i.e. without frequency transfer)

[TABLE]

Assuming the white noise correlator (3) that comes with the UV cutoff $1/\xi$ in k-space, this action is susceptible to a Wilsonian momentum-shell RG analysis Ostrovsky et al. (2006); Aleiner and Efetov (2006); Schuessler et al. (2009). Successively integrating out high energy modes down to $\lambda^{-1}/\xi$ ($\lambda\geq 1$) perturbatively, the action can be approximately mapped to itself with rescaled momenta, fields and coupling constants. If the velocity is kept constant, the two-loop RG equation for the flowing disorder strength $\tilde{K}\left(\lambda\right)$ reads Schuessler et al. (2009)

[TABLE]

Starting with the initial condition $\tilde{K}\left(1\right)\!=\!K$ the flow is to strong coupling where the perturbation theory leading to Eq. (7) breaks down. To find the energy scale $\Gamma$ where this happens (and below which the DOS is presumably constant), let us assert $\tilde{K}\left(\hbar v/\Gamma\xi\right)\sim 1$ which, in the limit of $K\ll 1$, leads to $\Gamma\propto\frac{\hbar v}{\xi}\sqrt{\frac{1}{K}}e^{-\pi/K}$ Schuessler et al. (2009) correcting for the factor 2 in the exponent as found from the SCBA. The DOS at the nodal point is expected to be governed by this emergent energy scale $\nu_{2}\left(E=0\right)\hbar v\xi\propto\Gamma$, in agreement with the KPM results in Fig. 1.

The Wilsonian RG calculation gave the correct exponential scale governing the disorder problem. However, it is not quantitative in the sense that numerical estimates for, say, the DOS could be obtained in the strong coupling limit. We will now show how the fRG method overcomes the difficulties mentioned above and use it to obtain quantitative results for the disorder induced DOS at the nodal point without any fitting parameters.

IV fRG approach

The fRG Metzner et al. (2012) introduces a flow parameter $\Lambda$ in the bare propagator and rewrites the many-body problem in a hierarchy of coupled flow equations for vertex functions with respect to $\Lambda$. The flow parameter is chosen such that for $\Lambda=\infty$, the vertex functions are known exactly and for $\Lambda=0$ the original problem is retained. We relegate a detailed discussion of technicalities to the appendix and only highlight the most important points and modifications related to use of the fRG with the replicated action.

To actually calculate expectation values and vertex functions from the replicated action, the replica limit $\left\langle O\right\rangle=\underset{R\rightarrow 0}{\mathrm{lim}}\frac{1}{R}\sum_{\alpha=1}^{R}\left\langle\mathcal{O}\left(\bar{\psi}^{\alpha},\psi^{\alpha}\right)\right\rangle_{\psi}$ is required, where $\left\langle\mathcal{O}\left(\bar{\psi},\psi\right)\right\rangle_{\psi}=\int D(\bar{\psi},\psi)\mathcal{O}\left(\bar{\psi},\psi\right)e^{-S\left[\bar{\psi},\psi\right]}$ stands for the standard functional average over a polynomial of fields $\mathcal{O}\left(\bar{\psi},\psi\right)$ Altland and Simons (2006). In a peturbative expansion (which is also at the heart of the fRG flow equations), thus only diagrams without closed fermion loops have a finite contribution in the replica limit. This also means that mixing of replica indices in the relevant diagrams is avoided. One can also show that the elastic nature of the interaction vertex derived from (6) is maintained along the flow. As a consequence, on the right hand side (rhs) of the flow equations the frequency integral as required for inelastic (true) interactions, is absent. Thus introducing $\Lambda$ via a Matsubara frequency cutoff scheme results in a Dirac delta function on the rhs which allows for a direct integration of the corresponding flow equations and results in a self-consistent hierarchy of equations for the vertices. So far no approximations have been made. To proceed, we truncate the hierarchy to order $K^{2}$. This is a pragmatic choice, that still goes beyond all diagrammatic schemes previously applied to disordered Dirac materials explicitly. Subsequently, we eliminate the interaction vertex in favor of the self-energy. The remaining self-consistency equation reads

[TABLE]

and is displayed in Fig. 3 diagrammatically: The term of order $K$ represents the SCBA approximation, c.f. diagram (i), the two second order terms are shown in diagrams (ii.a) and (ii.b) respectively. Although these diagrams would also appear in perturbation theory, the fRG approach (i) rigorously justifies the use of the self-energy dressed propagators and (ii) indicates how we could consistently go beyond order $K^{2}$ by allowing feedback for the vertex self-consistency equation.

To solve Eq. (8), we parameterize the self-energy using polar ($d\!=\!2$) or spherical ($d\!=\!3$) coordinates and proceed by iteration. We compute the DOS from Eq. (4). Further details are given in the appendix. In the $d=2$ Dirac case, the resulting DOS (red line) shows excellent agreement with the numerically exact KPM data and justifies the used order $K^{2}$ truncation a posteriori, well capable of capturing the exponential scale derived from Eq. (7).

On the pragmatic side, let us note that our fRG method also has advantages over the KPM method besides being analytic. For example, in Fig. 1, the KPM data for $\nu_{2}\left(E\right)$ shows a dip around $E=0$ that can only be resolved for small $K$ if the system size $L$ and expansion order $N$ is taken large. In comparison, the solution of Eq. (8) requires only a small fraction of computational effort.

V Disordered $d\!=\!3$ Weyl node

We now turn to the disorder induced DOS for a $d=3$ Weyl node. Here, weak disorder is irrelevant so that the DOS is maintained at zero. Only for $K>K_{c}$, disorder induces a finite DOS, see Fig. 2 for the KPM data (dots). These qualitative features were correctly predicted by the SCBA (blue line, see Refs. Fradkin (1986a, b); Biswas and Ryu (2014); Ominato and Koshino (2015)) and by the momentum shell RG treatment, see Refs. Syzranov et al. (2015, 2016). From the KPM, we find $K_{c}^{KPM}=4\pm 0.5$ (the precision is limited by finite size effects) while $K_{c}^{SCBA}\simeq 11$ (blue line) is off by more than a factor two. The one-loop RG result $K_{c}^{RG_{1}}=\pi^{2}\simeq 10$ can be improved with respect to the KPM value by adding two-loop corrections $K_{c}^{RG_{2}}=\pi^{2}/2\simeq 5$. However, quantitative predictions for the DOS in the strong-disorder phase cannot be obtained with the RG approach.

When compared to the $d\!=\!2$ case, the additional challenge for the fRG approach in the Weyl case is that the interesting disorder strengths $K\apprge K_{c}$ are not numerically small. Thus we assume that our $\mathcal{O}(K^{2})$ truncation of the fRG equations might cause a sizable error. Surprisingly, the fRG results (red line) yield $K_{c}^{fRG}\simeq 6$ and predict the available exact DOS for $K>7$ within an error of a few percent. On the one hand, we expect that the remaining numerical error of the fRG method could be systematically reduced by considering the fRG flow of the interaction vertex, which we leave for future research. On the other hand, this might not improve the accuracy for $K\simeq K_{c}$ where rare region effects which lie beyond any order of perturbation theory, are expected to dominate the DOS Nandkishore et al. (2014); Pixley et al. (2016a, 2017); Holder et al. (2017). However, it is known that their influence can be suppressed by choosing a different disorder model Pixley et al. (2016b).

VI Conclusion

We applied the fRG to treat the disorder problem at nodal points in two and three dimensions. From the resulting hierarchy of self-consistency equations, we calculate the bulk DOS and show that it is superior in accuracy to any other existing analytical approach. Suprisingly, for two dimensions, a truncation of the self-consistency equation at second order of $K$ is sufficient, while in three dimensions the accuracy could probably be increased with increasing order. We leave this suggestion for future work, along with the calculation of other experimentally relevant transport properties from fRG. More complicated disorder models, in particular vector disorder in two dimensions and its characteristic $\nu(E)$ behavior or scattering between multiple nodal points, as present in realistic materials, could also be studied in the future.

Acknowledgments

We thank Pavel Ostrovsky and Voker Meden for useful discussions. Numerical computations were done on the HPC cluster of Fachbereich Physik at FU Berlin. Financial support was granted by the Deutsche Forschungsgemeinschaft through the Emmy Noether program (KA 3360/2-1) and the CRC/Transregio 183 (Project A02).* *

Appendix: fRG with replicated action and solution of self-consistency

equation

fRG flow equations and vertex structure for the replica interaction.—

The fRG flow equations for the self energy $\Sigma$ and the interaction vertex $\Gamma$ have the form Metzner et al. (2012)

[TABLE]

and, in three-particle vertex truncation,

[TABLE]

where $\dot{G}_{\Lambda}=G_{\Lambda}(\partial_{\Lambda}[G_{0,\Lambda}^{-1}])G_{\Lambda}$ is the single-scale propagator and the multi-index $\left\{\alpha_{1}i\omega_{1}\mathbf{k}_{1}\sigma_{1}\right\}\equiv 1$ includes the relevant single-particle indices: replica index, Matsubara frequency, momentum and spin, respectively. We also use the notation $1_{\alpha_{j}}\equiv\left\{\alpha_{j}i\omega_{1}\mathbf{k}_{1}\sigma_{1}\right\}$ and $1_{\alpha_{j}i\omega_{k}}\equiv\left\{\alpha_{j}i\omega_{k}\mathbf{k}_{1}\sigma_{1}\right\}$ at our convenience and also abbreviate integrals and sums on the rhs as $\int_{1}\equiv\sum_{\alpha_{1}}\,\frac{1}{2\pi}\int d\omega_{1}\,\frac{1}{(2\pi)^{d}}\int\mathrm{d}\mathbf{k}_{1}\,\sum_{\sigma_{1}}$.

Starting from the inter-replica interaction $S_{\mathrm{dis}}$, Eq. (6) in the main text, we find the bare vertex by anti-symmetrization

[TABLE]

where we defined

[TABLE]

symmetric under the simultaneous exchange $\mathbf{k}_{1}^{\prime}\sigma_{1}^{\prime}\leftrightarrow\mathbf{k}_{2}^{\prime}\sigma_{2}^{\prime}$ and $\mathbf{k}_{1}\sigma_{1}\leftrightarrow\mathbf{k}_{2}\sigma_{2}$.

It is easy to see from the vertex flow equations (fRG flow equations and vertex structure for the replica interaction.—) that the locking of the replica and frequency indices, as present in the bare vertex (fRG flow equations and vertex structure for the replica interaction.—), is preserved in the flow (since the Green functions are frequency- and replica-diagonal). This means the flowing vertex is always of the form $\Gamma_{\Lambda}(1_{\alpha_{1}i\omega_{1}}^{\prime},2_{\alpha_{2}i\omega_{2}}^{\prime};1_{\alpha_{1}i\omega_{1}},2_{\alpha_{2}i\omega_{2}})$ or $\Gamma_{\Lambda}(1_{\alpha_{1}i\omega_{1}}^{\prime},2_{\alpha_{2}i\omega_{2}}^{\prime};1_{\alpha_{2}i\omega_{2}},2_{\alpha_{1}i\omega_{1}})$. Hence, in analogy to the bare vertex, we can write

[TABLE]

with $A_{\Lambda}^{\alpha_{1}^{\prime}i\omega_{1}^{\prime},\alpha_{2}^{\prime}i\omega_{2}^{\prime}}\left(\mathbf{k}_{1}^{\prime}\sigma_{1}^{\prime},\mathbf{k}_{2}^{\prime}\sigma_{2}^{\prime};\mathbf{k}_{1}\sigma_{1},\mathbf{k}_{2}\sigma_{2}\right)$ symmetric under the simultaneous exchange $\alpha_{1}^{\prime}i\omega_{1}^{\prime}\leftrightarrow\alpha_{2}^{\prime}i\omega_{2}^{\prime}$ as well as $\mathbf{k}_{1}^{\prime}\sigma_{1}^{\prime}\leftrightarrow\mathbf{k}_{2}^{\prime}\sigma_{2}^{\prime}$ and $\mathbf{k}_{1}\sigma_{1}\leftrightarrow\mathbf{k}_{2}\sigma_{2}$ like $A_{\infty}$.

Next, we need to leave out all terms on the rhs of Eqs. (9) and (fRG flow equations and vertex structure for the replica interaction.—) where the sums $\sum_{\alpha_{3},\alpha_{4}}$ on the rhs provide an extra factor of $R$ as these vanish in the replica limit, $\underset{R\rightarrow 0}{\mathrm{lim}}\frac{1}{R}\sum_{\alpha=1}^{R}\left\langle\mathcal{O}\left(\bar{\psi}^{\alpha},\psi^{\alpha}\right)\right\rangle_{\psi}\propto\underset{R\rightarrow 0}{\mathrm{lim}}\frac{1}{R}\sum_{\alpha=1}^{R}R\propto\underset{R\rightarrow 0}{\mathrm{lim}}R=0$. Note that fixing $1_{\alpha_{2}i\omega_{2}},2_{\alpha_{1}i\omega_{1}}$ in the first line of Eq. (fRG flow equations and vertex structure for the replica interaction.—) we can associate $\alpha_{1}$ with the replica index from multi-index $3$ or $4$, in the second line we have no such choice and the third line is always $\propto R$ and vanishes. If we would draw diagrams to represent Eq. (fRG flow equations and vertex structure for the replica interaction.—), the replica limit condition is equivalent of leaving out diagrams with internal fermion loops. We find

[TABLE]

Self-energy and vertex flow.—

Eventually, for the DOS we are interested in the Green function which involves the self-energy. Employing the replica-frequency locking of the vertex for the self-energy flow Eq. (9), we find

[TABLE]

Applying Eq. (fRG flow equations and vertex structure for the replica interaction.—), we find that only the second part avoids the replica sum leading to $\propto R$. The Green function locks all frequencies and replica indices appearing on the rhs of the self-energy flow equation,

[TABLE]

In Eq. (15), the function $A$ only appears with equal replica and frequency indices. We insert this structure in Eq. (fRG flow equations and vertex structure for the replica interaction.—) and obtain

[TABLE]

We can now drop the replica index $\alpha$ from our intermediate flow equations (15) and (Self-energy and vertex flow.—) and proceed to specify the flow parameter $\Lambda$ which was general so far.

Matsubara frequency cutoff.—

In its standard application to systems with inelastic (true) interactions, the fRG flow equations contain frequency integrals on the rhs Metzner et al. (2012). This integral is absent in Eq. (Self-energy and vertex flow.—) due to the elastic structure of the disorder induced interaction vertex. We can take this to our advantage and choose a Matsubara cutoff scheme which will allow exact integration of the flow equations. In the Matsubara cutoff scheme a multiplicative cutoff to the bare Green function is employed $G_{0,\Lambda}(1_{i\omega_{1}})=\theta\left(|i\omega_{1}|-\Lambda\right)G_{0}(1_{i\omega_{1}})$, the corresponding single scale propagator reads $\dot{G}_{\Lambda}(1_{i\omega_{1}})=\delta\left(|i\omega_{1}|-\Lambda\right)\tilde{G}_{\Lambda}(1_{i\omega_{1}})$ and $\dot{G}_{\Lambda}(1_{i\omega_{1}})G_{\Lambda}(1_{i\omega_{2}})=\delta\left(|i\omega_{1}|-\Lambda\right)\Theta\left(|i\omega_{2}|-\Lambda\right)\tilde{G}_{\Lambda}(1_{i\omega_{1}})\tilde{G}_{\Lambda}(1_{i\omega_{2}})$ where $\tilde{G}_{\Lambda}(1_{i\omega_{1}})=\left[G_{0}^{-1}(1_{i\omega_{1}})-\Sigma_{\Lambda}(1_{i\omega_{1}})\right]^{-1}$ Metzner et al. (2012) and $\theta(0)=1/2$ is understood by Morris Lemma Morris (1994).

We find

[TABLE]

Assuming $|\omega|>0$, we can now integrate both flow equations exactly over $\Lambda$ from $\Lambda=\infty$ to $\Lambda=0$ to find the physical self-energy $\Sigma=\Sigma_{\Lambda=0}$ and vertex function $A=A_{\Lambda=0}$. The initial condition for the interaction vertex is the bare interaction. Writing simply $G$ instead of $\tilde{G}_{\Lambda=0}$, we find

[TABLE]

Instead of the usual coupled fRG flow equations that have to be integrated, we thus have rephrased the disorder problem in terms of the coupled self-consistent Eqns. (17) and (18). Note that the above derivation did not depend on the three (or $N$-) particle vertex truncation in Eq. (fRG flow equations and vertex structure for the replica interaction.—) and thus, an extended set of coupled self-consistency equations would still be exact.

We turn back to our initial goal to find the DOS at the nodal point $E=0$. For this, we need the retarded real frequency self-energy, see Eq. (4) in the main text, that is connected to $\Sigma\left(i\omega\right)$ by an analytical continuation $i\omega=0+i0^{+}$ where $0^{+}$ is a positive real infinitesimal. After this step, we drop the frequency variable from now on. Let us emphasize that the appearance of a single frequency in the hierarchy of self-consistent equations is a remnant of the elastic nature of disorder scattering.

Solution correct to order $K^{2}$.—

Even the set of self-consistency equations (17) and (18) (with the three-particle vertex dropped) is difficult to solve without further approximations. To obtain the self-energy correct to at least $\mathcal{O}(K^{2})$, on the rhs of Eq. (18), is is sufficient to use the bare vertex Eq. (12). This is a pragmatic approach, which, however still goes beyond existing studies in the literature. We obtain from Eqns. (12) and (18)

[TABLE]

and further specialize to the spin-momentum structure needed for the self-energy flow Eq. (15)

[TABLE]

We combine Eq. (Solution correct to order $K^{2}$.—) with (17) and find the final self-consistency equation. Relabeling $\mathbf{k}_{1}\rightarrow\mathbf{k}$, $\mathbf{k}_{2}\rightarrow\mathbf{q}$ and $\mathbf{k}_{3}\rightarrow\mathbf{p}$ and using $``\cdot$” to indicate matrix products for the 2x2 matrix-valued Green functions, we arrive at Eq. (8) from the main text:

[TABLE]

If the feedback of the flowing vertex $A$ to the rhs of its own flow equation would be considered, this would yield two equations for $\Sigma$ and $A$ to be solved self-consistently.

Numerical solution of self-consistency equations.—

The self-consistency equation (20) can be solved numerically by iteration. We use dimensionless units (measuring momenta in $1/\xi$ and energies in $\hbar v/\xi$) and the dimensionless self-energy in $d=2$ (at the nodal point) can be parametrized as

[TABLE]

with $x,\phi$ polar coordinates. The term $M_{2}\left(x\right)$ has to be purely real (to avoid a spontaneous creation of chemical potential) and $>0$ for the retarded self energy. As a result, on the rhs of Eq. (20), we can chose $\mathbf{k}$ in say, the x-direction and also take only the $\sigma_{x}$ component of the product of Green functions (it can be checked that all other components vanish). The final self-consistency loop is then only for the functions $m_{2}\left(x\right)$ and $M_{2}(x)$, which turn out to be rather smooth. They can be discretized on a geometric grid for the variable $x$, the angular integrations can be done using a linearly spaced integration grid for the angles. We made sure that our results are converged with respect to the resolution of the discretization grids. Once $m_{2},\,M_{2}$ do not change any more under insertion on the rhs of Eq. (20), the DOS is computed from Eq. (4) using interpolation of the integrand and quadrature integration. Likewise, in $d=3$, the same strategy is applied using a parametrization in spherical coordinates $x,\phi,\theta$:

[TABLE]