Global Phase Diagram of a Dirty Weyl Liquid and Emergent Superuniversality

Pursuing complementary field-theoretic and numerical methods, we here paint the global phase diagram of a three-dimensional dirty Weyl system. The generalized Harris criterion, augmented by a perturbative renormalization-group (RG) analysis shows that weak disorder is an irrelevant perturbation at the Weyl semimetal(WSM)-insulator quantum critical point (QCP). But, a metallic phase sets in through a quantum phase transition (QPT) at strong disorder across a multicritical point (MCP). The field theoretic predictions for the correlation length exponent $\nu=2$ and dynamic scaling exponent $z=5/4$ at this MCP are in good agreement with the ones extracted numerically, yielding $\nu=1.98 \pm 0.10$ and $z=1.26 \pm 0.05$, from the scaling of the average density of states (DOS). Deep inside the WSM phase, generic disorder is also an irrelevant perturbation, while a metallic phase appears at strong disorder through a QPT. We here demonstrate that in the presence of generic, but strong disorder the WSM-metal QPT is ultimately always characterized by the exponents $\nu=1$ and $z=3/2$ (to one-loop order), originating from intra-node or chiral symmetric (e.g., regular and axial potential) disorder. We here anchor such emergent \emph{chiral superuniversality} through complementary RG calculations, controlled via $\epsilon$-expansions, and numerical analysis of average DOS across WSM-metal QPT. In addition, we also discuss a subsequent QPT (at even stronger disorder) of a Weyl metal into an Anderson insulator by numerically computing the typical DOS at zero energy. The scaling behavior of various physical observables, such as residue of quasiparticle pole, dynamic conductivity, specific heat, Gr$\ddot{\mbox{u}}$neisen ratio, inside various phases as well as across various QPTs in the global phase diagram of a dirty Weyl liquid are discussed.

A WSM, the prime example of a gapless topological phase of matter, is constituted by so called Weyl nodes that in the reciprocal space (Brillouin zone) act as the source and sinks of Abelian Berry curvature, and thus always appear in pairs [26]. In a nutshell, the Abelian Berry flux enclosed by the system determines the integer topological invariant of a WSM and the degeneracy of topologically protected surface Fermi arc. A question of fundamental and practical importance in this context concerns the stability of such gapless topological phase against impurities or disorder, inevitably present in real materials. Combining complementary field theoretic renormalization group (RG) calculations and a numerical study the role of randomness in various regimes of the phase diagram of a Weyl system to arrive at the global phase diagram, schematically illustrated in Fig. 1. A WSM can be constructed by appropriately stacking two-dimensional layers of quantum anomalous Hall insulator (QAHI) in the momentum space along the k z direction, for example. Thus, by construction a WSM inherits the two dimensional integer topological invariant of constituting layers of QAHI, and the momentum space skyrmion number of QAHI jumps by an integer amount across two Weyl nodes. As a result, the Weyl nodes serve as the sources and sinks for Abelian Berry curvature, and in a clean system WSM is sandwiched between a topological Chern and a trivial insulating phase, as shown in Fig. 1. In an effective tight-binding model a WSMinsulator quantum phase transition (QPT) can be tuned by changing the effective hopping in the k z direction, as demonstrated in Sec. II. In this work we first assess the stability of such clean semimetal-insulator quantum critical point (QCP) in the presence of generic randomness in the system. In this regard we come to the following conclusions.
1. By generalizing the Harris criterion [27], we find that WSM-insulator QCP is stable against sufficiently weak, but otherwise generic disorder (see Sec. III). Such an outcome is further substantiated from the scaling analysis of disorder couplings, suggesting that any disorder is an irrelevant perturbation at such a clean QCP.
2. Subscribing to an appropriate controlledexpansion (see Sec. III), we demonstrate that a multicritical point (MCP) emerges at stronger disorder, where the WSM, a band insulator (either Chern or trivial) and a metallic phase meet. The critical semimetal residing at the phase boundary between a WSM and an insulator then becomes unstable toward the formation of a compressible metal through such a MCP. The exponents capturing the instability of critical excitations toward the onset of a metal are: (a) correlation length exponent (CLE) ν = 2, and (b) dynamic scaling exponent (DSE) z = 5/4 to the leading order in the -expansion. These two exponents also determine the scaling behavior of physical observables across the anisotropic critical semimetal-metal QPT.
3. By following the scaling of DOS along the phase boundary between the WSM and insulator with increasing randomness in the system, we numerically extract ν and z at the MCP across the critical semimetal-metal QPT [see Fig. 2]. Numerically extracted values of these two exponents are ν = 1.98 ± 0.05 and z = 1.23 ± 0.05 [see Sec. III B], which are in excellent agreement with our prediction from the leading order -expansion.
We now turn our focus on the WSM phase. The study of disorder effects in topological phases of matter has re- cently attracted a lot of attention, leading to a surge of analytical  and numerical [50][51][52][53][54][55][56][57][58][59][60][61][62][63] works. In particular, the focus has been concentrated on massless Dirac critical point separating two topologically distinct insulators (electrical or thermal), as well as inside a Dirac and Weyl semimetal phases. Even though the effects of generic disorders have been studied to some extent theoretically [30,36,[41][42][43], most of the numerical works solely focused on random charge impurities (for exception see Refs. [53,55]). By now there is both analytical and numerical evidence that chemical potential disorder when strong enough drives a QPT from the WSM to a diffusive metal, leaving its imprint on different observables, e.g., average DOS, specific heat and conductivity. Deep inside the WSM phase, the system possesses various emergent symmetries (see Table III), such as a continuous global chiral U (1) symmetry that is tied with the translational symmetry of a clean noninteracting WSM in the continuum limit [64]. In the absence of both inversion and timereversal symmetries, the simplest realization of a WSM with only two Weyl nodes is susceptible to sixteen possible sources of elastic scattering, displayed in Table III. They can be grouped in eight classes, among which only four preserve the emergent global chiral symmetry (intranode scattering), while the remaining ones directly mix two Weyl nodes with opposite (left and right) chiralities  (internode scattering) 1 . As we demonstrate in the paper, such characterization of disorders based on the chiral symmetry allows us to classify the WSM-metal QPTs in the presence of generic disorder.
To motivate our theoretical analysis, we now discuss the possible microscopic origin of disorders in the Weyl materials. Furthermore, knowing this in future may facilitate a control over randomness in experiments on these materials. For example, chemical potential disorder can be controlled by modifying the concentration of random charge impurities. Random asymmetric shifts of chemical potential between the left and right chiral Weyl cones correspond to the axial potential disorder. Therefore, in an inversion asymmetric WSM such disorder is always present. Magnetic disorder is yet another type of chiral symmetry preserving (CSP) disorder, and the strength of random magnetic scatterers can be efficiently tuned by systematically injecting magnetic ions in the system 2 . In contrast, all chiral symmetry breaking (CSB) disorders cause mixing of two Weyl nodes and in an effective model for WSMs they stem from various types of random bond disorder that also cause random fluctuation of band-width (see Appendix D). Therefore, strength of CSB disorder may be tuned by applying inhomogeneous pressure (hydrostatic or chemical) in the Weyl materials. Since the WSMs are found in strong spin-orbit coupled materials, a random spin-orbit coupling can be achieved when hopping (hybridization) between two orbitals with opposite parity acquires random spatial modulation. Yet another CSB but vector-like type of disorder is a random axial Zeeman coupling. Its source is the different  Table II: Numerically extracted critical strength of disorder for WSM-metal QPT (Wc), dynamic scaling exponent (z) and correlation length exponent (ν) in the presence of four individual disorder potentials that mix two Weyl nodes (nonchiral disorder), obtained from the scaling of average density of states. For field theoretic analysis of internode scatterers or non-chiral disorder see Sec. VII, and for data analysis of average density of states see Sec. VI and Figs. 14, 16. g-factor of two hybridizing bands that touch at the Weyl point [65][66][67]. Therefore, when magnetic impurities are injected in the system such disorder is naturally introduced, and depending on the relative strength of the g-factor in different bands, one can access regular (intranode) or axial (internode) random magnetic coupling. Finally, two different types of CSB mass disorders that tend to gap out the Weyl points are represented by random charge-or spin-density-wave order, depending on the microscopic details [68]. These disorders correspond to random scalar and pseudo-scalar mass in the field theory language. Due to their presence, Weyl nodes are gapped out in each disorder configuration, but the sign of the gap is random from realization to realization, and in the thermodynamic limit the nodes remain gapless. To the best of our knowledge, it is currently unknown how to tune the strength of all individual sources of elastic scattering in real Weyl materials. Nevertheless, we elucidate how all possible disorders can be obtained from a simple effective tight-binding model on a cubic lattice for a WSM with two nodes (see Appendix D), allowing us to numerically investigate the effects of generic disorder in this system.
Here we address the stability of a disordered WSM (i) in the field-theoretical framework by using two different renormalization-group (RG) schemes: (a) an mexpansion about a critical disorder distribution, where m = 1 − m, with the Gaussian white noise distribution realized as m → 0, and (b) d = 2 − d-expansion about d l = 2, the lower-critical spatial dimension for WSMmetal QPT, and (ii) lattice-based numerical evaluation of average DOS by using the kernel polynomial method (KPM) [69] in the presence of generic chiral symmetric disorder [see Fig. 3] as well as non-chiral disorder [see Fig. 4]. Comparisons between the field theoretic predictions and numerical findings for all chiral disorders are given in Table I Figure 3: Scaling of numerically evaluated [using the kernel polynomial method [69]] average density of states (ADOS) in dirty Weyl semimetals (WSMs) in the presence of (a) potential, (b) axial, (c) axial current and (d) current disorder for weak to strong disorder regime, in a cubic lattice of linear dimension L = 220. Notice that for weak enough disorder ADOS (E) ∼ |E| 2 for |E| 1. In the metallic phase, appearing for strong enough disorder, ADOS at zero energy (0) becomes finite. Around a (nonuniversal) critical strength of disorder W = Wc the ADOS scales as (E) ∼ |E| for |E| 1. Since (E) ∼ |E| d z −1 , the dynamic scaling exponent z ≈ 1.5 across the WSM-metal quantum phase transitions (QPTs), irrespective of the nature of the elastic scatterers. All four disorders preserve the emergent global chiral symmetry and represent intranode scattering, see Table III. Numerically extracted critical exponents across WSM-metal QPTs and their comparison with the field theoretic predictions are displayed in Table I, suggesting an excellent agreement between these two methods and emergence of a superuniversality across WSM-metal QPT. Lattice model is defined in Sec. II and for detailed analysis of ADOS see Sec. VI. In a suitably chosen lattice model the axial current disorder corresponds to magnetic impurities (Table III).

Figs. 3 and 4)
, depicting that DOS scales as (E) ∼ |E| 2 for small energy (E), when generic disorder is sufficiently weak.
2. We show in Sec. V that irrespectively of the details of two distinct -expansions, in the presence of a CSP disorder, the WSM-metal QPT takes place through either a QCP (when either potential or axial potential disorder is present) or a line of QCPs (when both types of scalar disorder are present simultaneously), characterized by critical exponents obtained from the leading order in -expansions, where = m or d , and = 1 corresponds to the physical situation. Therefore, irrespective of the nature of elastic scatterers, the universality class of the WSM-metal QPT in the presence of a CSP disorder is unique, and we name such universality class chiral superuniversality.
3. In Sec.VI we carry out a thorough numerical analysis of DOS in the presence of all four CSP disorders, obtained by using KPM from a lattice model (see Fig. 3). Within the numerical accuracy we find that z ≈ 1.5 and ν ≈ 1 across possible CSP disorder driven WSM-metal QPTs (see Fig. 15 and also Table I Table III for definition and field theoretic nomenclature]. All non-chiral disorders drive WSM-metal QPTs, around which (W = Wc) the ADOS scales roughly as (E) ∼ |E|, indicating z ≈ 1.5 even in the presence of strong internode scattering [see Table II  We point out that the notion of superuniversality is realized rather sparsely in condensed matter systems. Most prominent examples in this regard include the quantum Hall plateau transitions [70][71][72] and one-dimensional disordered superconducting wires [73]. Therefore, dirty Weyl semimetal represents, to the best of our knowledge, the only example of a three-dimensional system exhibiting superuniversality.
It is worth mentioning that for sufficiently strong disorder the metallic phase in a Weyl system undergoes a second continuous QPT into an Anderson insulating phase [28,53,74]. In Sec. IX, we address the metalinsulator Anderson transition (AT), but only in the presence of random charge impurities. Our central achievements regarding the fate of the AT in strongly disordered Weyl metal are the following: 1. We show that a Weyl metal undergoes a second transition at stronger disorder into an Anderson insulator (AI) phase. By numerically computing the typical density of states (TDOS) at zero energy [ t (0)] we show that t (0) vanishes smoothly across the Weyl metal-AI QPT, while displaying critical and single-paramter scaling. In particular, t (0) is pinned at zero in the WSM and AI phases, while it is finite inside the entire metallic phase. By contrast, the average DOS at zero energy [ (0)] remains finite in the metallic as well as AI phases, while being zero only in the weakly disordered WSM. Otherwise, (0) decreases smoothly and monotonically across the Weyl metal-AI QCP.
2. We demonstrate that TDOS at zero energy displays single-parameter scaling across both (a) WSMmetal and (b) metal-AI QPTs. Specifically the order-parameter exponent for t (0), β t , defined as t (0) ∼ |δ| βt , where δ defines the reduced distance from transition point, is β t = 1.80 ± 0.20 across the WSM-metal QPT (which is different from the one for the average DOS at zero energy for which β a = 1.50 ± 0.05).
3. We show that inside the metallic phase the mobility edge, separating the localized states from the extended ones reside at finite energy. With increasing strength of disorder the mobility edge slides down to smaller energy and across the AT the entire energy widow is occupied by localized states.
The rest of the paper is organized as follows. In Sec. II, we introduce a simple tight-binding model for a Weyl system and discuss possible phases and the phase transitions in the clean limit. In Sec. III, we demonstrate the effects of generic disorder near the clean WSM-insulator QCP, and perturbatively address the effects of strong disorder. In Sec. IV we set up the theoretical framework for addressing the role of randomness deep inside the WSM phase, and introduce the notion of m and d expansions for perturbative treatment of disorder. This section is rather technical and readers familiar with the formalism or interested in physical outcomes may wish to skip it. We devote Sec. V to the effects of CSP disorder and promote the notion of chiral superuniversality. Detailed numerical analysis of the scaling of DOS is presented in Sec. VI. Effects of CSB disorder are discussed in Sec VII and scaling of various physical observables, such as DOS, specific heat, conductivity, etc., across the WSM-metal QPT is discussed in Sec. VIII. We discuss the Anderson transition of the metallic phase at stronger disorder in Sec. IX. Concluding remarks and a summary of our main findings are presented in Sec. X. Some additional technical details have been relegated to the Appendices.

II. LATTICE MODEL FOR WEYL SYSTEM
Let us begin the discussion with a lattice realization of chiral Weyl fermions in a three-dimensional cubic lattice. Even though in most of the commonly known Weyl materials, such as the binary alloys TaAs and NbP, Weyl fermions emerge from complex band structures in noncentrosymmetric lattices, their salient features can be captured from a simple tight-binding model The two-component spinor is defined as ψ k = (c k,↑ , c k,↓ ), where c k,s is the fermionic annihilation operator with momentum k and spin/pseudospin projection s =↑, ↓, and σs are standard Pauli matrices. We here choose where a is the lattice spacing. The first term gives rise to two isolated Weyl nodes along the k z axis at k z = ±k 0 z , where with the following choice of pseudospin vectors The second term in Eq. (3), namely N M 3 (k) = t 0 [2 − cos(k x a) − cos(k y a)], plays the role of a momentum dependent Wilson mass [57,58]. We subscribe to this tight-binding model in Secs. III B, VI and IX to numerically study the effects of randomness in various regimes of a dirty Weyl system.  (5) and (3). Here, NI and CI respectively represents trivial (normal) and Chern insulators. Weyl nodes in the WSM phase are always located along the kz direction. Respectively WSM1,2,3 supports one, two and one pair of Weyl nodes. The projection of the Weyl nodes on the xy plane in these phases are at the (0, 0) point, (0, π) and (π, 0) points, and (π, π) point. This model therefore supports translationally active topological phases. The transitions between the WSM and insulating phases (solid lines) and the ones between two distinct WSM phases (dashed lines) are continuous. We emphasize that there is no symmetry distinction among these phases.
With the above chosen form of the Wilson mass only a single pair of Weyl nodes is realized at k 0 = 0, 0, ±k 0 z , whenm z ≤ 1, wherem z = m z /t 0 . For 1 ≤m z ≤ 2, two pairs of Weyl nodes are found at (0, π, ±k 0 z ) and (π, 0, ±k 0 z ). Whenm z ≥ 2, again a single pair of Weyl points is realized at (π, π, ±k 0 z ). Notice that whent z < 1 three distinct realizations of WSM phase are separated by intervening insulating phases, representing three dimensional Chern insulator, while the insulating phases realized form z < −1 andm z > 4 are topologically trivial. Hence, by tuningm z one can drive the system through WSM-insulator QPT as long ast z < 1, wherẽ t z = t z /t 0 . However, fort z > 1 the noninteracting model supports only trivial insulating phases whenm z /t z < −1 or (m z −4)/t z > 1. Direct and continuous transitions between two different WSMs take place without an intervening insulating phase whent z > 1. Therefore, our proposed model in Eqs. (5) and (3) supports translationally active topological phases with band-inversions at non-Γ points in the Brillouin zone [9,75] and the resulting phase diagram is shown in Fig. 5.
For the sake of simplicity, we hereafter only consider the parameter regime −1 <m z < 1 andt z ≤ 1, so that only a single pair of Weyl fermions is realized at k 0 = (0, 0, ± cos −1 |m z /t z |). In the vicinity of these two points the Weyl quasiparticles can be identified as left and right chiral fermions, respectively. A WSM can be found when |m z /t| ≤ 1 and the system becomes an insulator for |m z /t| > 1. Even though we here restrict our analysis within the aforementioned parameter regime, this analysis can be generalized to study the semimetal- insulator QPTs in various other regimes shown in Fig. 5.
Within this parameter regime, to capture the Weyl semimetal-insulator QPT which occurs along the linẽ t z /m z = 1, we expand the tight-binding model around the Γ = (0, 0, 0) point of the Brillouin zone to arrive at the effective low energy Hamiltonian where v = ta is the Fermi velocity in the xy plane and b = t z a 2 /2 bears the dimension of inverse mass. For ∆ = t z −m z < 0 the system becomes an insulator (Chern or trivial). On the other hand, when ∆ > 0, the lattice model describes a WSM. The QPT in this clean model between these two phases takes place at ∆ = 0. Hence, ∆ plays the role of a tuning parameter across the WSMinsulator QPT. The QCP separating these two phases is described by an anisotropic semimetal, captured by the Hamiltonian H Q (0) in Eq. (6), that in turn also determines the universality class of the transition. Notice that the expansion of the lattice Hamiltonian [see Eq. (5)] also yields terms ∼ k 2 x and ∼ k 2 y and higher order (from the Wilson mass), which are, however, irrelevant in the RG sense, and therefore do not affect the critical theory for the WSM-insulator QPT. Hence, we omit these higher gradient terms for now. We will discuss the paramount importance of such higher gradient terms close to the CSB disorder driven WSM-metal QPT in Sec. VII. Next we address the stability of this quantum critical semimetal against disorder in the system using scaling theory and RG analysis.

III. EFFECTS OF DISORDER ON SEMIMETAL-INSULATOR TRANSITION
The imaginary time (τ ) action associated with the low energy Hamiltonian [see Eq. (6)] reads as In the proximity to the Weyl semimetal-insulator QPT, the system can be susceptible to both random charge and random magnetic impurities, and their effect can be captured by the Euclidean action where V j (x) are random variables. The effect of random charge impurities is captured by V 0 (x), while V ⊥ (x) and V z (x) represents random magnetic impurities with the magnetic moment residing in the easy or xy plane and in the z direction (denoted here by x 3 for notational clarity), respectively, which we allow due to the anisotropy of the Hamiltonian [see Eq. (6)]. All types of disorder are assumed to be characterized by Gaussian white noise distributions.
The scale invariance of the noninteracting action [see Eq. (7)] mandates the following scaling ansatz: τ → e l τ , (x, y) → e l (x, y) and x 3 → e l/2 x 3 , followed by the rescaling of the field operator ψ → e −5l/4 ψ, where l is the logarithm of running RG scale. The scaling dimension of the tuning parameter ∆ is then given by [∆] = 1, implying that ∆ is a relevant perturbation at the WSMinsulator QCP, located at ∆ = 0. The scaling dimension of the tuning parameter ∆ plays the role of the correlation length exponent (ν) at this QCP, implying ν = 1. In the presence of disorder, as we show in Appendix A, the Harris stability criterion [27] can be generalized for the WSM-insulator QCP with the quantum-critical theory of the form given by Eq. (6), but in a system with the topological or monopole charge c [see Eq. (A1)]. The generalized Harris criterion then suggests that WSM-insulator QCP in clean system remains stable against sufficiently weak disorder only if and d * as the effective spatial dimensionality of the system under the coarse graining procedure. At the WSMinsulator QCP ν = 1, and the critical excitations residing at ∆ = 0 are therefore stable against weak disorder when c = 1 (regular WSM). We next analyze the effects of disorder on the WSM-insulator QCP using a RG approach. The same outcome can be arrived at from the computation of inverse scattering life-time (1/τ ) within the framework of self-consistent Born approximation [see Appendix H].

A. Perturbative RG analysis
After performing the disorder averaging in the action [see Eq. (8)] within the replica formalism, we arrive at the replicated Euclidean action where a, b are replica indices. Notice that here we have replaced k 2 3 → k n 3 , with n as an even integer so that such deformation of spectrum does not change the symmetry of the system. We we will show that such deformation of the quasiparticle spectrum allows us to control the perturbative RG calculation in terms of disorder coupling. The above imaginary-time action (S) remains invariant under the space-time scaling (x, y) → e l (x, y), x 3 → e l/n x 3 and τ → e zl τ . At the bare level the scale invariance of the free part of the action requires the field renormalization factor Z ψ = e −(2+1/n)l and ψ → Z −1/2 ψ ψ. From this scaling analysis we immediately find that the scaling dimension of disorder couplings is [∆ j ] = −1/n, for j = 0, ⊥, z. Therefore, at the WSMinsulator QCP, characterized by n = 2, disorder is an irrelevant perturbation, in accordance with the prediction from the generalized Harris criterion, implying the stability of this QCP against sufficiently weak randomness. Note that disorder couplings are marginal in a hypothetical limit n → ∞, for which the system effectively becomes a two-dimensional Weyl semimetal. Therefore, perturbative analysis in the presence of generic disorder is controlled via an n -expansion, where n = 1/n, about n → ∞, following the spirit of -expansions about upper or lower critical dimension [76] and infinite monopole charge [77,78].
Upon integrating out the fast Fourier modes within the momentum shell Λe −l < k ⊥ < Λ, where k ⊥ = k 2 x + k 2 y , 0 < k 2 3 < ∞ and accounting for pertubative corrections to one-loop order (see Fig. 6), we arrive at the following flow equations for X = v, b n , j = 0, ⊥, z, β Q ≡ dQ/dl is the β-function for the running parameter Q, and for brevity we omit the hat notation in Eq. (11). In the above flow equations, we have kept only the leading divergent contribution that survives as n → ∞. Inclusion of subleading divergences yields only nonuniversal corrections, as shown in Appendix B. The β−function for in-plane Fermi velocity (v) and b n leads to a scale dependent DSE Note that in this formalism the random charge-impurities do not generate any new disorder, allowing us to depict the RG flow in the (∆, ∆ 0 ) plane, as shown in Fig. 7(a). The coupled RG flow equations (11) support only two fixed points: • (∆, ∆ 0 , ∆ ⊥ , ∆ 3 ) = (0, 0, 0, 0), which has only one unstable direction along the ∆-direction that serves as the tuning parameter for WSM-insulator QPT. This fixed point stands as a QCP in the four dimensional coupling constant space. The correlation length exponent at this QCP is ν −1 = 1. All disorder couplings are irrelevant perturbations at this QCP [see the blue dot in Fig. 7(a)].
• (∆, ∆ 0 , ∆ ⊥ , ∆ 3 ) ≈ (0, n /2, 0, 0) stands as a multicritical point (MCP) with two unstable directions. At this MCP the WSM, an insulator and the metallic phase meet. Two correlation-length exponents  Table III: Various types of disorder represented by fermionic bilinears (j = 1, 2, 3), together with their symmetries under pseudo time-reversal (T ), parity (P), continuous chiral rotation (Uc) and charge-conjugation (C). The disorder couplings are represented by ∆N and Σµν = [γµ, γν ]/(2i). Note that true time-reversal symmetry in WSM in already broken. The pseudo time-reversal symmetry T is generated by an anti-unitary operator γ0γ2K, where K is complex conjugation, such that T 2 = −1 (The true time-reversal operator is γ1γ3K). The parity operator is P = γ0, while the chargeconjugation operator is C = γ2. The continuous chiral symmetry (Uc) is generated by γ5, the generator of translational symmetry in the continuum limit in a clean Weyl semimetal [64]. The Hermitain γ matrices satisfy standard anti-commutation relation {γµ, γν } = 2δµν for µ, ν = 0, 1, 2, 3, 5, and for explicit representation of γ-matrices see Sec. IV A. Here and × signify even and odd under a symmetry operation, respectively. With a slightly different tight-binding model, where (2)], the axial current corresponds to magnetization, temporal and spatial tensors to spin-orbit and axial magnetization, respectively. However, such microscopic details do not alter any physical outcome.
are ν −1 M = n determining the relevance of disorder coupling ∆ 0 , which drives the anisotropic critical semimetal [described byĤ Q (0)] into a diffusive metallic phase, and ν −1 = 1 that determines the relevance of the tuning parameter ∆, controlling the WSM-insulator transition. The DSE for critical semimmetal-metal QPT is z = 1 + n 2 . Therefore, for a three-dimensional anisotropic critical semimetal-metal QPT, setting n = 1/2, the critical exponents are ν M = 2 and z = 1.25.
The RG flow and the resulting phase diagrams are shown in Fig. 7(a) and 7(b), respectively. At the multicritical point the average DOS scales as (E) ∼ |E| d/z−1 ≈ |E| 1.4 to one-loop order. Beyond the critical strength of disorder system becomes a metal where the average DOS at zero energy [ (0)] is finite and the order parameter exponent β = (d − z)ν = 3.5 determines the scaling of (0) according to (0) ∼ δ β = δ 3.5 in the metallic phase, where δ = (∆ 0 − ∆ * 0 ) /∆ * 0 is the reduced disorder coupling from the critical one at ∆ 0 = ∆ * 0 . Next we numerically demonstrate (a) stability of WSM-insulator QCP at weak disorder, (b) emergence of a metallic phase through a MCP at finite disorder coupling that masks the direct transition between WSM and insulator by numerically computing the average DOS using the kernel polynomial method. As a natural outcome of this exercise, we will also show that numerically extracted values of the exponents, z and ν, at the MCP, associated with the critical excitations-metal QPT agree with the predictions from the leading order n -expansion 3 .
B. Scaling of density of states near WSM-insulator QCP: Numerical demonstration of the MCP Before we discuss the scaling behavior of the average DOS along the WSM-insulator phase boundary and inside the metallic phase, setting in through the instability of critical semimetallic phase, let us point out some crucial subtle issues associated with such analysis. Note that the average DOS of the critical semimetal [described bŷ H Q (0) in Eq. (6)] vanishes as (E) ∼ |E| 3/2 , while that in the WSM phase vanishes as (E) ∼ |E| 2 . But, in the insulating phase average DOS displays hard gap. Based on scaling analysis we expect WSM, insulator and the critical semimetal to be stable against sufficiently weak disorder. We exploit these characteristic features to pin the WSM-insulator phase boundary for weak disorder. On the other hand, for stronger disorder onset of a metallic phase can be identified from the existence of finite average DOS at zero energy. Following these diagnostic tools we arrive at the phase diagram of a Weyl materials re-siding in the close proximity to the WSM-insulator QPT; see Fig. 2 (left). We are ultimately interested in exposing the existence of a MCP in the (m z , t z ) plane [the red dot in Fig. 7(a)] which has two relevant directions. One of them controls critical semimetal-metal QPT, while the other one drives WSM-insulator QPT. Since we consider the former transition, our focus will be restricted on the black dashed line shown in Fig. 2.
More specifically, we here compute the average DOS by employing the kernel polynomial method [69] starting with the tight-binding model, introduced in Eqs. (2), (3) and (5), and staying in the close vicinity of m z /t 0 = 0.5 and t z /t 0 = 0.5 (see the phase diargam in Fig. 5). The tight-binding model is implemented on a cubic lattice with periodic boundary conditions in all three directions and the linear dimensionality of the system in each direction is L = 140. Even though average DOS is a selfaveraged quantity, we perform average over 20 random disorder realization to minimize the residual statistical fluctuations, compute 4096 Chebyshev moments and take trace over 12 random vector to obtain average DOS. For the sake of simplicity we here account for only random charge impurities. Potential disorder is distributed uniformly and randomly within the range [−W, W ]. The scaling of average DOS can be derived in the following way.
Since we are following only one relevant direction associated with the MCP, effectively it can be treated as a simple QCP across which various physical observables (such as average DOS) display single parameter scaling. Note that total number of states N (E, L) in a ddimensional system of linear dimension L, below the energy E is proportional to L d , and in general is a function of two dimensionless parameters L/ξ and E/E 0 . Here, ξ ∼ δ −ν is the correlation length that diverges at the QCP, located at δ = 0, where δ = W −Wc Wc is the reduced distance from the QCP, located at W = W c . Consequently, the correlation energy, defined as E 0 ∼ δ νz vanishes as the QCP is approached from either side of the transition [80]. Following the standard formalism of scaling theory we then can write where G is an universal but unknown scaling function. Therefore, from the definition of average DOS (E, L) = L −d dN (E, L)/dE we arrive at the following scaling form where F is yet another universal, but typically unknown scaling function. However, we can access the behavior of the scaling function in different regimes along the black dashed line shown in Figs. 2 (left), which we exploit to compute critical exponents characterizing the critical semimetal-metal QPT across the MCP. In the final step we have used the fact that average DOS remains particlehole symmetric, but on average. Note we will use exactly the same scaling function deep inside the WSM phase in the presence of generic disorder, discussed in Sec. VI. First of all, notice that average DOS (0) is pinned to zero along the phase boundary between the WSM and insulator for weak enough disorder, as shown in Fig. 8(a). Therefore, critical semimetal separating these two phases remains stable against weak disorder and the the nature of the WSM-insulator direct transition remains unchanged for weak enough randomness. However, beyond a critical strength of disorder, W c = 1.20 ± 0.05, (0) becomes finite and metallicity sets in through the MCP, see Figs. 2 (left) and 8(a). Beyond this point there exists no direct transition between the WSM and an insulator. Also note that for W W c , (E) ∼ |E| 1.5 as shown in Fig. 2 (right). Now we consider very close proximity to the MCP, located at W = W c along the disorder axis. At this MCP average DOS becomes independent of δ, yielding F (x) ∼ x d z −1 . By comparing (E) with E, obtain the DSE associated with critical semimetal-metal QPT to be z = 1.23 ± 0.05, see Fig. 8 Next we move into the metallic phase, but continue to follow the black dashed line from Fig. 2 (left). In the metallic phase (0) becomes finite [see Fig. 8(a)]. Thus to the leading order F (x) ∼ x 0 and consequently (0) ∼ δ (d−z)ν . With the prior notion of z = 1.23 ± 0.05, now by comparing (0) vs. δ we obtain the CLE at the MCP associated with the critical semimetal-metal QPT to be ν = 1.98 ± 0.05, as shown in Fig. 8(c).
Therefore, numerically extracted values of two critical exponents, namely ν and z, at the MCP associated with the critical semimetal-metal QPT match extremely well with the field theoretic prediction obtained from an n -expansion introduced in this work, which allows to  Figure 9: (a) Collapse of average DOS at finite energy (obtained in system with L = 140) across the multi-critical point (MCP) shown in Fig. 2 (left). All data collapse reasonably well onto two branches corresponding to anisotropic semimetal (upper branch) and metallic phase (lower branch), which tend to meet in the critical regime. (b) Data collapse of average DOS at zero energy for different system sizes inside the metallic phase, appearing across the MCP. These two data collapses are obtained with numerically extracted critical exponents z = 1.23 and ν = 1.98 [see Fig. 8].
control the RG calculation by tuning the flatness of the quasiparticle dispersion along the k z direction: a controlled ascent from two spatial dimension.
We now discuss two different types of data collapses across the disorder-driven MCP. The results are shown in Fig. 9. First we focus on the largest system with L = 140. From Eq. (14), upon neglecting the finite size effects, we compare (E)|δ| −(d−z)ν vs. |E||δ| −νz along the black line from Fig. 2(left). With numerically obtained values of ν and z we find that all data nicely collapse onto two branches (corresponding to the anisotropic semimetal and metallic sides of the QPT), which meet in the critical regime, as shown in Fig. 9(a). Next we compare the average DOS at zero energy in the metallic phase, namely (0)L d−z vs. L 1/ν δ, in systems of different sizes (L), as shown in Fig. 9(b). We also obtain excellent finite-size data collapse for a wide range of system sizes using already numerically extracted values of ν and z. Therefore, field-theoretic predictions and numerical findings across the disorder-driven MCP are in good agreement with each other. Next we address the effects of disorder inside the WSM phase by pursing complimentary field theoretic and numeric approaches.

IV. DIRTY WEYL SEMIMETAL: MODEL AND SCALING ANALYSIS
In this section, we set up the field theoretical framework to analyze the role of disorder when the system is deep inside the WSM phase. We will introduce the notion of two different -expansions: (a) an m −expansion about a critical disorder distribution, where m = 1 − m with Gaussian white noise distribution recovered as m → 0; (b) an d −expansion, with Gaussian white noise distribution from outset, about the lower critical dimension d c = 2 for WSM-metal QPT, where d = d − 2, and therefore for three spatial dimensions d = 1.

A. Hamiltonian and action
The effective low energy description of WSM can be obtained by expanding the lattice Hamiltonian [see Eq. (5)] around the Weyl nodes located at k 0 = (0, 0, ±k 0 z ), with k 0 z = cos −1 ( mz tz ). The resulting low energy Hamiltonian reads where v = ta, v z = a t 2 z − m 2 z , and the momentum is measured from the Weyl nodes. For simplicity we hereafter take the Fermi velocity to be isotropic, v = v z , so that the low energy Hamiltonian becomes rotationally symmetric. Upon performing a unitary rotation with U = σ 0 ⊕ σ 3 , the above Hamiltonian assumes a quasirel- Hermitian matrices, and summation over repeated spatial indices is assumed. To close the Clifford algebra of five mutually anticommuting matrices we define γ 5 = τ 3 ⊗σ 0 . Two sets of Pauli matrices σ µ and τ µ respectively operate on spin/pseudospin and valley or chiral (left and right) indices. The low energy effective Hamiltonian enjoys variety of emergent discrete and continuous symmetries. The above Hamiltonian is invariant under a pseudo-time-reversal symmetry, generated by antiunitary operator T = γ 0 γ 2 K, where K is the complex conjugation, a charge conjugation symmetry, generated by C = γ 2 , and parity or inversion symmetry generated by P = γ 0 . Furthermore, the Hamiltonian [see Eq. (15)] also possesses a global chiral U (1) symmetry, generated by γ 5 , which in the low energy limit corresponds to the generator of translational symmetry [64].
To incorporate the effects of disorder we consider the following minimal continuum action for a dirty WSM σ,τ is the fermionic creation operator near the Weyl point at τ k 0 for τ = ± (left/right) and with spin σ =↑, ↓, whileΨ = Ψ † γ 0 , as usual. Various disorder fields ϕ N , coupled to the fermion bilinears, are realized with different choices of 4 × 4 matrices, N , as shown in Table III. Notice that the matrices associated with four types of disorder anticommute with γ 5 and represent chiral symmetric disorder, while for the other four types of disorder [N, γ 5 ] = 0 and the corresponding disorder vertex breaks the U (1) chiral symmetry. As we demonstrate in this paper, this global chiral symmetry plays a fundamental role in classifying the disorder-driven WSM-metal QPTs.

B. m expansion in three dimensions
We assume that the disorder field obeys the distribution [38,79] or in the momentum space and the limit m → 0 corresponds to the Gaussian white noise distribution, which we are ultimately interested in. This form of the white noise distribution stems from the following representation of the d−dimensional δ−function [43] We now carry out the scaling analysis of the continuum action for a WSM given by Eq.
Due to linearly dispersing low-energy quasiparticles, a WSM corresponds to z = 1 fixed point, and in d = 3 the engineering dimension of the disorder strength is [∆ N ] = m − 1. A first implication of this result is that the white noise disorder, m = 0, is irrelevant close to the WSM ground state in d = 3. Second, for m = 1, the disorder is marginal and we use that to introduce the deviation from this value as an expansion parameter m = 1 − m.
The β−function (infrared) for the disorder coupling ∆ N in the m expansion is given in terms of its scaling dimension in Eq. (20), yielding in d = 3. Therefore, to obtain the explicit form of this β−function in terms of the disorder couplings, we have to compute the DSE and the anomalous dimension of the disorder field. The former is obtained from the fermion self-energy with the diagram shown in Fig. 10(a), while the latter is found from the vertex diagram in Fig. 10(b). Evaluation of these two diagrams have been carried out using field-theoretic method (see Appendix C). Alternatively, one may choose to integrate out the fast modes within the momentum shell Λe −l < k < Λ, with Λ as an ultraviolet cutoff in the momentum, to arrive at the RG flow equations for ∆ N .

Self-energy and dynamic scaling exponent
We first show the computation of the self-energy diagram, shown in Fig. 10(a), yielding the dynamical exponent and the anomalous dimension for the fermion field within the regularization scheme defined by the parameter m = 1 − m, the deviation from the critical disorder distribution. All the integrals are therefore performed in d = 3. The divergent part of the integral appears as a pole ∼ 1/ m , analogously to the case of the dimensional regularization where the deviation from the upper or lower critical space-time dimension plays the role of an expansion parameter. To find renormalization constants, we use minimal subtraction, i.e. we keep only divergent part appearing in the corresponding diagrams.
The action [see Eq. (16)] without the disorder yields the inverse free fermion propagator with v 0 as the bare Fermi velocity. Taking into account the self-energy correction, the inverse dressed fermion propagator is with Σ(iω, k) as the self-energy. After accounting for all possible disorders, we arrive at the following compact expression for the self-energy (see Appendix C for details) where as the dimensionless disorder strength, and for brevity we here drop the hat symbol in the final expression. From the above expression of the self-energy, together with the renormalization con- with v as the renormalized Fermi velocity, we arrive at the expression for the fermion-field renormalization (Z Ψ ) and velocity renormalization (Z v ) This equation then yields the anomalous dimension for the fermion field Furthermore, the renormalization factor Z v enters the renormalization condition for the Fermi velocity Finally, the β−function of the Fermi velocity is β v = (1 − z)v, which together with Eq. (28) determines the DSE

Vertex correction: Anomalous dimension of disorder field
We now turn to the vertex correction due to the disorder, shown in Fig. 10(b), which yields the anomalous dimension of the disorder field. As shown in Appendix C, the vertex represented by the matrix N receives the correction of the form The corresponding renormalization condition that determines the renormalization constant Z ϕ N for the disorder field reads with Z Ψ given by Eq. (26). The above condition in turn yields the anomalous dimension of the disorder field as which we then use to write the explicit form of the β−function Eq. (21) in terms of the disorder couplings.
Alternatively, one may take the Gaussian white noise distribution in Eq. (17) with m → 0 from the outset. In that case, the engineering dimension of the disorder coupling is equal to 2 − d, since z = 1 in a clean WSM. Therefore, d = 2 is the lower critical dimension in the problem and we can use d = d − 2 as an expansion parameter, following the spirit ofexpansion [30,36,37,40,41,44,55,76]. In this scheme, after performing the disorder averaging using the replica method, the imaginary time action assumes a similar form of Eq (10).
Within the framework of the d expansion only the temporal (frequency-dependent) component of self energy acquires a disorder-dependent correction to the leading order. The self-energy correction due to disorder reads as with the function f 1 (∆ j ) given by Eq. (24), and ∆ j Λ d /(2πv 2 ) → ∆ j . This result is obtained from Eq. (C1) with d = 2 + d and m = 0. As a result, the field renormalization factor Z Ψ = 1+f 1 (∆ j )/ d and the velocity renormalization factor is , we obtain the leading order RG flow equation for the Fermi velocity which yields a scale dependent dynamic exponent z = 1 + f 1 (∆ j ). The seemingly different expressions for the flow equation and DSE in these two schemes stems from underlying different methodology of capturing the ultraviolet divergences of various diagrams. However, such details do not alter any physical outcome.

V. CHIRAL SYMMETRIC OR INTRA-NODE DISORDER
We first focus on chiral-symmetric disorders. For a single pair of Weyl fermions there are four such disorders, namely chemical potential, axial potential, current and axial current disorders, as shown in Table III. With appropriate lattice model axial current disorder corresponds to magnetic impurities and from here onward we use this terminology. We will address the effect of weak and strong chiral symmetric disorder using both m and d expansions.

A. m expansion
Let us first analyze this problem pursuing the m expansion. Using Eqs. (21), (29), (31) and (32), we obtain the following RG flow equations for the coupling constants to the leading order in m The above set of flow equations supports a line of quantum critical points in the ∆ V − ∆ A plane, determined by where the quantities with subscript " * " represent their critical values for WSM-metal QPT. The RG flow in this plane is shown in Fig. 11(a). The line of QCPs also determines the WSM-metal phase boundary, and the corresponding phase diagram in the ∆ V − ∆ A plane is shown in Fig. 11(b). At each point of this line of QCPs the DSE and CLE are respectively given by Therefore, for the Gaussian white noise distribution, realized for m = 1, we obtain z = 3/2 and ν = 1 from the leading order m expansion. If the bare value of either the chemical potential or axial potential disorder strength is zero, the quantum-critical behavior is governed by the QCP corresponding to the disorder of a nonvanishing bare value [43]. This QCP features the critical exponents of the same value to the one-loop order as in the case of the quantum-critical line, given by Eq. (37). From the RG flow equations [see Eq. (35)], we find that both magnetic and current disorder are always irrelevant perturbations, at least to the leading order in the m -expansion. In the ∆ X − ∆ Y plane, where X = V, A and Y = M, C the RG flow diagram is shown in Fig. 12(a) and the corresponding phase diagram is shown in Fig. 12(b). Importantly, the QPT separating the metallic and the semimetallic phase in any ∆ X − ∆ Y plane is governed by the QCP located at ∆ X, * = 3 m /8. The phase boundary between these two phases is determined by the irrelevant direction at this QCP. Therefore, across the entire WSM-metal phase boundary in these planes the universality class of the QPT is identical and characterized by z = 1 + m /2 and ν −1 = m .

B. d expansion
The RG flow equations for the chiral symmetric disorder coupling constants within the framework of dexpansion are . These coupled flow equations also support only a line of QCPs in the ∆ V −∆ A plane, as we previously found from Eq. (35) using m -expansion, now determined by In the presence of only potential disorder we find z = 3/2 and ν = 1 [30,36,37,40,41,44]. Notice that if we start with only magnetic or current disorder, the axial disorder gets generated from Feynman diagrams (c) and (d) in Fig. 6. Thus, to close the RG flow equations, we need to account for ∆ A coupling from the outset, and the resulting RG flow equations read  Table I and II]. The numerical methodology is described in details in Sec. VI.

C. Chiral superuniversality
From the discussion in previous two subsections, we can conclude that in the presence of chiral-symmetric disorder in a WSM, the semimetal-metal QPT takes place either through a QCP or a line of QCPs. The location of the line of QCPs and the resulting phase boundaries are nonuniversal and thus dependent on the RG scheme. However, the universal quantum critical behavior with chiral symmetric disorder couplings is insensitive of these details, at least to the leading order in the expansion parameter, and all QPTs in the four-dimensional hyperplane of disorder coupling constants, are characterized by an identical set of critical exponents, namely z = 1 + /2 and ν −1 = , with = 1. Therefore, emergent quantum critical behavior for strong chiral-symmetric disorder stands as a rare example of superuniversality, and we name it chiral superuniversality. Next we demonstrate emergence of such superuniversality across WSM-metal QPT by numerically analyzing the scaling of average DOS in the presence of generic chiral symmetric disorder.

VI. NUMERICAL DEMONSTRATION OF CHIRAL SUPERUNIVERSALITY
Motivated by the field-theoretic prediction of emergent chiral superuniversality across the WSM-metal QPTs driven by CSP disorder, next we numerically investigate the scaling of average DOS across such QPTs. Since (0) vanishes and is finite in the WSM and metallic phases, respectively, it can be promoted as a bonafide orderparameter across the WSM-metal QPT [51,53,55,56,58,62]. In addition, such analysis endows an opportunity to extract the critical exponents for the transition nonperturbatively and, at the same time, test the validity of the proposed scenario for chiral superuniversality. The WSM phase is realized from the tight-binding model, defined through Eqs. (3) and (5), which we implement on a cubic lattice of linear dimension L. For numerical analysis we always set m z = 0, and for current disorder take t = t z = 1 = t 0 , while t = 1 = t 0 , t z = 1 2 for remaining seven types of elastic scatterers [see Table III], in the clean model, given by Eqs. (2)-(5). We use lattice realizations of disorder introduced in Appendix D. We impose periodic boundary condition in all three directions. The average DOS is computed by using the kernel polynomial method [69]. The average is taken over 20 random real-  Table III]. First column shows the scaling of ADOS (E) vs. E around the critical strength of disorder (W = Wc). The second column depicts the scaling of ADOS at zero energy (0) vs. δ, the reduced distance from the critical disorder defined as δ = W −Wc Wc . In the third column we display (E)δ −(d−z)ν vs. |E||δ| −νz for weak (W < Wc) and strong (W > Wc) disorder and |E| t(= 1). All data collapse onto two branches. The top branch represents the metallic phase, while the lower one WSM. Note that these two branches meet at large values of |E||δ| −νz , corresponding to the quantum critical regime. All data in first three columns are obtained from a system of linear dimension L = 220. The finite size data collapse inside the metallic phase is shown in the forth column, where we compare (0)L d−z vs. δL 1/ν for 100 ≤ L ≤ 220. Notice that all data collapse onto one branch for small to moderate values of δL 1/ν , with the numerically extracted values of the critical exponents z and ν, quoted in the figure and summarized in Table I. The quality of the data collapse progressively worsens for larger values of δL 1/ν due to the existence of a second QPT of a three-dimensional dirty Weyl metal into the Anderson insulating phase, discussed in Sec. IX. Scaling of ADOS and data analysis are discussed in details in Sec. VI.
ization of disorder that minimizes the residual statistical error in average DOS, which is a self-averaged quantity. We typically compute 4096 Chebyshev moments and take trace over ∼ 12 random vectors to compute average DOS. All types of disorder are distributed uniformly and randomly within the range [−W, W ]. The scaling theory for average DOS has already been discussed in Sec. III B. Thus, we can readily start from the final expression of the general scaling form of the average DOS, presented in Eq. (14) and continue with our numerical analysis.  Table III for definitions]. In a specific lattice model the former two sources of elastic scattering respectively correspond to random spin-orbit coupling and random axial-magnetization. Methods of analysis and system size are identical to the ones in Fig. 15. Final results of our analysis are quoted in Table II.

A. Numerical analysis with random intra-node scatterers or chiral-symmetric disorder
We begin the discussion on the effects of randomness on WSM by focusing on the intra-node or chiral symmetric disorder. Let us first focus on the quantum critical regime and for now we assume that the system size is sufficiently large so that we can neglect the L-dependence in Eq. (14). In this regime the scaling function must be independent of δ, dictating F (x) ∼ x  Table I. Within the numerical accuracy, we always find z ≈ 1.5 in excellent agreement with the field-theoretic result, obtained from the leading order expansions.
Next we proceed to the metallic side of the transition, where average DOS at zero energy becomes finite. From the scaling function in Eq. (14), we obtain (0) ∼ δ (d−z)ν . Thus by comparing (0) vs. δ, we extract the CLE ν, using already obtained value of the DSE z, as shown in the second column of Fig. 15. The numerically found CLE is also quoted in Table I, and within numerical accuracy ν ≈ 1 always, irrespective of the nature of CSP disorder. Once again we find an excellent agreement of numerically extracted values of ν with the one obtained from the leading order -expansions. These two results strongly support the picture of chiral superuniversality.
To test the quality of our numerical analysis we search for two types of data collapse. First, we compare (E)|δ| −ν(d−z) vs. |δ| −νz |E|, motivated by the scaling form of average DOS, displayed in Eq. (14). Using numerically obtained values of ν and z, we find that for energies much smaller than the bandwidth (|E| 1), all data collapse onto two separate branches for all four disorders, as shown in the third column of Fig. 15. While the top branch corresponds to the metallic phase, the lower one stems from the WSM phase and eventually these two branches meet in the quantum critical regime.
Finally, we demonstrate a finite size data collapse for (0) for different system sizes (L) by focusing on the metallic side of the transition. Setting E = 0 in Eq. (14), we obtain (0) = L z−d F (0, δL 1/ν ). Hence, we compare (0)L d−z vs. δL 1/ν and find an excellent data collapse for 100 < L < 220, using numerically obtained values of ν and z for all four disorders, as shown in the fourth column of Fig. 15. The data collapse becomes systematically worse for large values of δ or stronger disorder due to the existence of a second transition that takes the system from a metallic phase to an Anderson insulator. Therefore, our thorough numerical analysis provides a valuable and unprecedented insight into the nature of the WSM-metal QPTs driven by generic chiral symmetric disorder, and staunchly supports the proposal of an emergent chiral superuniversality across such QPTs.

B. Numerical analysis with random inter-node scatterers or non-chiral disorder
Motivated by the intriguing possibility of realizing an emergent superuniversality we further seek to examine its robustness in the presence of inter-node scattering (also referred as non-chiral disorder). In the simplest version of a Weyl semimetal comprised of only two Weyl nodes there are four sources of internode scattering, highlighted in Table III, and their lattice realization is shown in Appendix D. We rely on the scaling of average DOS in the presence of non-chiral disorder as well, and all the parameters and numerical strategies are identical to the ones pursued for chiral symmetric (intranode) disorder. The analyses of average DOS in various regimes of the phase diagram of disordered WSM are performed in the same fashion. The locations of WSM-metal QPT are shown in Fig. 14 (lower row), and numerically extracted values of two critical exponents ν and z are reported in Table II. The details of the data analysis are displayed in Fig. 16.
Within the numerical accuracy we find that the WSMmetal QPT driven by CSB disorder is also characterized by ν ≈ 1 and z ≈ 1, 5. Therefore, the chiral superuniversality appears to be generic in a dirty WSM, and the WSM-metal QPTs belong to the same universality class, irrespective of the nature of impurities. Hence, there is a substantial evidence supporting emergence of superuniversality across WSM-metal QPT. Such an intriguing outcome further motivates us to understand the effect of internode scattering in a WSM from a field theoretic point of view, which we present in the following section by carrying out two different -expansions, described in Secs. IV B and IV C.

VII. CHIRAL SYMMETRY BREAKING OR INTER-NODE DISORDER
In a WSM constituted by a single pair of Weyl nodes, there are four CSB disorders, namely temporal and spatial components of a tensor disorder, which in a suitable lattice model respectively represents spin-orbit and axial magnetic disorder, as well as scalar and pseudoscalar mass disorder, see Table III. We will address the effects of weak and strong CSB disorder by using both m and d expansions.

A. m expansion
Within the framework of an m expansion the RG flow equations to one-loop order read as Therefore, individually each CSB disorder is always an irrelevant perturbation, at least to the leading order in the m expansion, and as such does not lead to any QPTs. However, in the absence of chiral symmetry all four disorder couplings are present and to address the critical properties in this situation we recast the above flow equations in terms of newly defined coupling constants as where ∆ ± V = ∆ SO ± ∆ AM , ∆ ± S = ∆ S ± ∆ P S . The above set of RG flow equations supports a line of QCPs determined by the equation Notice that if we tune the CSB disorders, so that ∆ − V = ∆ − S = 0, these two coupling constants do not get generated through quantum corrections, and the plane with ∆ − V = ∆ − S = 0, shown in Fig. 17, remains invariant under the RG. The RG flow in this plane is shown in Fig. 17(a), and the corresponding phase diagram is presented in Fig. 17(b). The WSM-metal phase boundary in the ∆ + V − ∆ + S plane is determined by the line of QCPs, given by Eq. (43), qualitatively similar to the situation in the presence of potential and axial disorders, as shown in Fig. 11. However, these two scenarios are fundamentally different in the sense that while the DSE z = 1 + /2, with = m or d , is fixed along the entire line of QCPs in the ∆ V − ∆ A plane, it varies continuously along the

line of QCPs in the ∆
where the quantity with subscript " * " denote the critical value for WSM-metal transition. Such continuously varying DSE leaves its signature in critical scaling of various physical observables, as we discuss below, and qualitatively mimics the picture of Kosterlitz-Thouless transition. Notice that the end point of such line of QCPs on the ∆ + V axis reside in the ∆ SO − ∆ AM plane at ∆ SO = ∆ AM = 3 m /4, and the RG flow in this plane is shown in Fig. 18(a). The phase diagram of a dirty WSM containing only spin-orbit and axial magnetic disorder in this plane is shown in Fig. 18(b), with z = 1+5 m , which is directly obtained from Eq. (44) by setting ∆ + S = 0. It is worth pointing out that in the ∆ SO − ∆ AM plane the phase boundary between the WSM and metallic phase is set by the irrelevant parameter associated with the QCP, while when such QCP percolates through ∆ + V −∆ + S plane in the form of a line of QCPs, it is determined by the relevant direction at each point on the line of QCPs.

B. d expansion
Next let us address the effects of CSB disorder within the framework of an d expansion. In this method the RG flow equations become very complicated due to the ultraviolet divergent contribution arising from the class of the Feynman diagrams shown in Fig. 6 (c) and (d), and it is challenging to decode the emergent quantumcritical phenomena. Thus we attempt to unearth critical properties by focusing on various coupling constant subspaces that remain closed under the RG, at least to the leading order. Let us first focus on spin-orbit or axial magnetic disorder. The RG flow equations read where X = SO, AM . Notice that even though the bare theory contains only spin-orbit or axial magnetic disorders, the CSP axial disorder gets generated and in order to keep the RG flow equations closed we need to include the latter from the outset. The coupled flow equations support one QCP, located at ∆ X, * = 9 d /10, ∆ A, * = 6 d /5 [30,41]. The RG flow diagram is shown in Fig. 19(a), and the resulting phase diagram is displayed in Fig. 19(b). Note that QCP obtained in the absence of the CSB disorders, located at ∆ A, * = d /2 now becomes unstable in the presence of either spin-orbit or axial magnetic disorder, and a new QCP results from the competition between these two disorders, as mentioned above. This outcome although bears contrast with our previously reported results obtained from m expansion, still shows some qualitative similarities, as we argue below. Notice that the DSE and CLE at the new QCP, shown in Fig. 19(a), are respectively given by As a result the mean DOS at the QCP diverges as (E) ∼ |E| −5/11 for d = 1 or d = 3, since z > d. Hence, both -expansions give rise to diverging DOS at the QCP controlled via spin orbit and axial magnetic disorder. Although the calculated values of DSE depend on RG scheme, to the leading order they do not differ significantly, z = 6 for m = 1, and z = 11/2 for d = 1, while ν = 1, is independent of the RG scheme.

C. Mass disorder
We now discuss the role of mass disorder in WSMs. It should be noted that a WSM can be susceptible to two different types of mass disorder (a) scalar mass disorder and (b) pseudo-scalar mass disorder. Both of them break the chiral symmetry, but can be rotated into each other by the generator of the chiral symmetry γ 5 . The flow equation for mass disorder within the framework of an expansion reads as for X = S, P S, α m = 8/3 and α d = 2, j = m, d corresponds to m and d expansions, respectively. Hence, by itself scalar or pseudoscalar mass disorder does not drive any WSM-metal QCP, at least within the leading order in -expansions. In this regard both m and d expansions yield an identical result. Finally, we discuss yet another interesting aspect of mass disorder, when it coexists with the axial one. The flow equations in the presence of these two disorders are for X = S, P S, where ∆ − = ∆ A − ∆ X ,α m = 8/3, α d = 2, and respectively j = m, d corresponds to m and d expansions. These two flow equations support a line of QCPs, determined by The location of such line of QCPs is regularization dependent (throughα j ), along which the DSE and CLE, given by are identical in both -expansion schemes. Therefore, in a WSM with these two disorders the DSE continuously increases from z = 3/2 in an unbounded fashion, while the CLE remains fixed. The numerical investigation of such interesting possibility is left for a future work.

D. Why chiral superuniversality is so robust?
Leaving aside the interesting possibilities of realizing such as line of QCPs with continuously verying critical exponents, perhaps the most urgent issue to be addressed is the following: Why does the disorder-driven WSMmetal QPT always display same universality class, characterized by ν ≈ 1 and z ≈ 1.5?
Answer to question in presence of intra-node or chiralsymmetric disorders has already been provided in Sec. V. Note that scaling dimension of any disorder coupling in a d-dimensional WSM is [∆ a ] = 2z − d. But at all CSB disorder driven QCPs controlling the WSM-metal QPT z > d irrespective of the RG methodology. Therefore, even though the bare values of CSP disorders in latticebased simulations are set to be zero, as done in the ones discussed in Sec VI B, they do get generated as we approach the Weyl points through the coarse graining procedure. Ultimately the CSP disorder becomes relevant at CSB disorder driven WSM-metal QCPs. As a result, the dirty system even though initially tends to flow toward the QCPs with z > d, described in this section, it flows back toward the chiral symmetric QCP or line of QCPs shown in Fig. 11(a). This is the reason why the WSM-metal QPTs are always characterized by CLE ν ≈ 1 and DSE z ≈ 1.5 (within numerical accuracy), the characteristics of the proposed chiral superuniversality. The above argument is very generic and does not depend on the number of Weyl nodes. Therefore, in any lattice system, we expect WSM-metal QPT to always belong to the chiral superuniversality. This outcome can be anchored from the RG calculation in the presence of all eight possible disorder couplings (since in strong disorder regime all disorders get generated even if the bare coupling for some specific channel is set to be zero), as shown in Appendix E within the framework of both m and d expansions. Such analysis confirms that only the line of QCPs, defined through Eq. (36) or Eq. (39), and shown in Fig. 11(a), ultimately controls the quantumcritical behavior. This strongly supports the above argument in favor of chiral superuniversality under generic circumstances 4 .
The specific tight-binding model we subscribe in this work [see Sec. II] also contains Wilson mass that bears higher gradient terms, such that τ 3 b ⊥ (k 2 x + k 2 y ), with b ⊥ = t 0 a 2 /2. The scaling dimension of such operator is [b ⊥ ] = z − 2. Hence, the higher gradient terms are irrelevant at clean WSM fixed point ([b ⊥ ] = −1) as well as at the chiral symmetric line of QCPs ([b ⊥ ] = −1/2), but becomes relevant at pure CSB disorder-driven QCPs (since z > d > 2). This is also the reason why chiral superuniversality is such a generic and utmost stable situation.
Nevertheless, we believe pure CSB disorder driven QCPs (with z > d) can in principle be realized in a numerical simulation performed in momentum space, where forward/intranode/CSP scattering processes can be suppressed deliberately and higher gradient terms can be avoided completely. Such an analysis is an interesting exercise of a pure academic interest, and we leave it for a future investigation.

VIII. QUANTUM CRITICAL SCALING OF PHYSICAL OBSERVABLES
As demonstrated in the previous two sections that QPT from a WSM to a diffusive metal can be driven by different types of elastic scatters, and the critical exponents are remarkably independent of the actual nature of randomness. We here highlight how these exponents can affect the scaling behavior of measurable quantities as the Weyl material undergoes this QPT 5 . 4 We note that the quality of data collapses for CSB disorders, shown in Fig. 16, is slightly less pronouced than those for CSP disorder, displayed in Fig. 15, which can qualitatively be understood in the following way. In the presence of only inter-node scatterers system first tends to flow toward the line of QCPs set by purely CSB disorder, discussed early in this section. Only when disorder gets sufficiently strong the intra-node disorder becomes relevant and the system starts flowing toward the line of QCPs discussed in Sec. V. The system then gets stuck in the crossover regime dominated by CSB disorder, and consequently the data collapse (involving finite energy states) becomes slightly less prominent. To achieve equally good quality data collapse even in the presence of CSB disorder we therefore need to subscribe to larger systems, which can be numerically challenging. 5 In spite of the emergent superuniversality, the putative line of QCPs driven by CSB disorders with continuously varying DSE z > d may leave its imprint on the physical observables in the crossover regime before the CSP disorders take over and ultimately the system flows toward the chiral symmetric quantum-• Residue of quasiparticle pole: As the WSM-metal QCP is approached from the semimetallic phase, the residue of quasiparticle pole vanishes and beyond the critical strength of disorder Weyl fermions cease to exist as sharp quasiparticle excitations, similar to the situation for two-dimensional Dirac fermion-Mott insulator QPT in the presence of a strong Hubbard interaction [81,82]. The residue of quasiparticle pole (Z) vanishes according to where η Ψ is the fermionic anomalous dimension at the critical point located at the disorder strength ∆ = ∆ * . Within the framework of an d expansion η Ψ = 0 to the leading order in d , and one needs to account for twoloop diagrams to obtain finite η Ψ . In contrast, in the m expansion we obtain nontrivial fermionic anomalous dimension even to the one-loop order, and η ψ ∼ , as shown in Eq. (27). Therefore, at the WSM-metal QCP, the quasiparticle spectrum displays a branch-cut and the critical point represents a strongly coupled non-Fermi liquid. Alternatively, the residue of quasiparticle pole plays the role of a bonafide order parameter on the semimetallic side. It is worth mentioning that the disappearance of residue of quasiparticle pole has recently been tracked in quantum Monte Carlo simulations for Hubbard model in two-dimensional honeycomb lattice [82], and we can expect that future numerical work can verify our proposed scaling form in Eq. (51) across the disorder driven WSMmetal QPTs. The Fermi velocity scales as v ∼ |δ| ν(z−1) , and since z > 1 at the QCP or the quantum-critical line, the Fermi velocity vanishes at the transition to the metallic phase.
• Average density of states: The most widely studied physical quantity in numerical simulations across the WSM-metal QPT is the average DOS [51,53,55,56,58,62]. Since throughout the paper we have already extensively used the average DOS to characterize phases, for the sake of completeness we here review only its salient features. We can infer the scaling form of the average DOS in the thermodynamic limit L → ∞ in different phases by using its scaling function [see Eq. (14)]. In the quantum critical regime (E) should be independent of δ, yielding Q (E) ∼ E d/z−1 . Inside the WSM phase, the average DOS scales as W (E) ∼ δ (1−z)dν |E| 2 . In the metallic phase average DOS at zero energy is finite and scales as (0) ∼ δ (d−z)ν . From the quoted values of DSE and CLE, it is straightforward to find the scaling of average DOS in these three regimes of the phase diagram in a dirty WSM, which we have used in the numerical analysis of this observable in the previous sections.
critical line with z = 3/2 and ν = 1. In that sense the physical observables we discuss in this section can also distinguish between different types of disorder (inter-node vs intra-node).
• Conductivity: The optical conductivity (σ) at T = 0 can as well serve as an order parameter across the WSMmetal QPT, and assumes the following scaling ansatz for frequency (Ω) much smaller than the bandwidth [42] where G is yet another unknown universal scaling function. This scaling form remains operative even at finite temperature as long as Ω T , i.e., in the collisionless regime. In the collision dominated regime at T Ω, the dc conductivity also assumes a similar scaling form as in Eq. (52), upon replacing the frequency (Ω) by temperature (T ) [37,52,83]. In the WSM side of the transition, the OC vanishes linearly with Ω and scales as Inside the critical regime the OC scales as σ Q (Ω) ∼ Ω (d−2)/z . In the presence of strong CSP disorder z = 3/2, and the optical conductivity inside the quantum critical regime thus vanishes as σ Q (Ω) ∼ Ω 2/3 . Since for non-chiral disorder the DSE is typically much bigger than in the presence of chiral symmetric one, the optical conductivity vanishes with a weaker power as Ω → 0 when the system is still dominated by CSB disorder before CSP disorder takes over. Hence, in this regime the system becomes more metallic in the presence of CSB disorder than with only CSP disorder. Inside the metallic phase, the optical conductivity becomes finite and scales as σ M (0) ∼ δ ν(d−2) as Ω → 0. Within the leading order m or d expansions, the conductivity of the metal is therefore always independent of the actual nature of elastic scatterers, since ν −1 = m or d , and m = 1, d = 1. Otherwise, weak disorder (such as potential) causes enhancement of optical conductivity without altering σ ∼ Ω scaling [42] [see also Appendix F for a simple derivation].
• Specific heat: The specific heat (C v ) also displays distinct scaling behavior in three regimes of the phase diagram of a dirty WSM. The scaling of specific heat at temperature much smaller than bandwidth follows the ansatz [36] where H is also an unknown universal scaling function. In the WSM phase, H(x) ∼ x d(z−1)/z and the specific heat scales as C V ∼ δ dν(1−z) T d , so that we recover T 3 dependence for three dimensional Weyl fermion. Inside the metallic phase, H(x) ∼ x 1−d/z , yielding C V ∼ δ ν(d−z) T and we obtain T -linear specific heat, similar to the situation in Fermi liquids. By contrast, inside the critical regime H(x) ∼ x 0 , yielding C V ∼ T 3/z . Therefore, the specific heat, analogous to the conductivity, displays distinct power-law dependence on temperature inside the quantum critical regime depending on the dominant source of disorder, while its scaling inside the WSM and metallic phases is insensitive to the nature of random impurities. Hence, the scaling of specific heat can be used to extract the extent of the critical regime and crossover boundaries among different phases of a dirty Weyl system at finite temperature [55].
• Mean-free path: The quasiparticle mean-free path (L ) also follows the critical scaling where J is a universal, but unknown scaling function, with energy much smaller than bandwidth. At the QCP (δ = 0) the mean-free path should be independent of δ, implying J (x) ∼ x −1/z . Therefore, inside the quantum critical fan, the mean-free path at zero energy diverges as L (E) ∼ E −1/z . In the metallic phase, J (x) ∼ x 0 as x → 0, leading to finite mean-free path at zero energy, and L (0) ∼ δ −ν . On the other hand, in the WSM phase, the mean-free path Since at all disorder driven QCPs z > 1, L W (E) decreases with increasing disorder, indicating propensity toward the onset of a metallicity in the system. • Grüneisen parameter : Yet another directly measurable quantity is the Grüneisen parameter, defined as γ = α/C P , where α is the thermal expansion parameter, and C P is the specific heat measured at constant pressure. The Grüneisen ratio in the WSM phase γ W ∼ T −4 , while inside the critical regime γ Q ∼ T −(1+d/z) . Inside the metallic phase γ M ∼ T −2 . Therefore, the Grüneisen parameter displays distinct power law behavior in three different phases of a dirty WSM.
Fascinating scaling behavior can also be observed for the magnetic Grüneisen ratio, defined as Γ H = (∂M/∂T ) H /C H , where M ∝ H is magnetization, C H is the molar specific heat, and H is the magnetic field strength. In the presence of sufficiently weak randomness when Landau quantization is sharp (ω c τ 1, where ω c is cyclotron frequency and τ is scattering lifetime) and dominates over the Zeeman coupling, leading to Γ H ∼ T −4/z . On the other hand, in the presence of strong elastic scattering when ω c τ 1 the Landau levels are sufficiently broadened and the dominant energy scale is set by Zeeman coupling, yielding Γ H ∼ T −2 , which is independent of dimensionality (d) or DSE (z). Therefore, for a fixed weak magnetic field, as the strength of impurities is gradually increased, the magnetic Grüneisen ratio should display a smooth crossover from T −4 to T −2 dependence. Note that such a crossover will take place even before the system enters the quantum critical regime and will persist in the metallic regime as well, since elastic scattering is strong in these two phases. i,α Here L 3 is the system size, |i, α is the eigenstate with site index i and orbital index α(= 1, 2) at energy E i,α . As previously discussed, average DOS is a self-averaging quantity so to minimize statistical fluctuations we only extract the disorder-averaged smoothened data, which we carry out by computing N m = 1024 Chebyshev moments and performing disorder average over 20 random disorder realizations. On the other hand, LDOS and TDOS are not self-averaging quantities. Therefore, numerical extraction of TDOS is extremely demanding for which we compute N m = 8192 moments and perform disorder average over 100 random disorder realization to construct the TDOS. To further suppress statistical fluctuations in TDOS we average over a small cube of size N s = L 3 s L 3 , and we here take L s = 4. Such averaging is justified since translational symmetry gets restored after disorder averaging has been performed.
The scaling of average DOS and TDOS over a wide range of disorder strength is shown in Fig. 20(a). Note that in the WSM phase both average DOS and TDOS at zero energy are pinned to zero, which then become finite across the WSM-metal QPT at W c,1 = 1.65 ± 0.05. Therefore, either average DOS or TDOS can be identified as a bonafide order-parameter to pin the WSM-metal QCP. Respectively these two quantities scale as near the WSM-metal QCP, with β a = 1.50 ± 0.05, β t = 1.80 ± 0.20, as shown in Fig. 20(b). Even though the numerical errorbar for β t is quite large, in general, we expect it to be different from β a , as their difference, ∆β = β t −β, is intimately tied with the multifractal dimension of the wavefunction across a disorder-driven QPT [84][85][86][87]. However, more precise determination of β t requires additional extensive numerical simulation. Therefore, we leave this issue as a subject for a future investigation. Inside the compressible diffusive metallic phase these two quantities increase monotonically and follow each each other up to a moderate strength of disorder W * ≈ 3.5. Upon further increasing strength of disorder the TDOS smoothly vanishes around W c,2 = 9.30 ± 0.25. Therefore, a metal-insulator transition (MIT) takes place at W = W c,2 , commonly known as AT. Note that the average DOS decreases monotonically across the AT, but remains non-critical, as shown in Fig. 20(a). In Fig. 21(a) we present the scaling of TDOS with the number of Chebyshev moments (N m ). We explicitly compute TDOS from moderate to strong disorder regime (6 ≤ W ≤ 10), in the close vicinity of the AT, for N m = 2048, 4096 and 8192. From the scaling of t (0) vs. N m we conclude that AT (identified with t (0) → 0) takes place around W c,2 = 9.30 in the N m → ∞ limit. Therefore, we can conclude that a three-dimensional diffusive Weyl metal is a stable phase of matter for moderately strong disorder, which ultimately undergoes a QPT into the Anderson insulator phase for sufficiently strong disorder. Across the AT the TDOS at zero energy display single-parameter scaling with β = 1.5 ± 0.1.
Recall that for weak disorder average DOS a (E) ∼ |E| 2 and around the WSM-metal QCP it scales as (E) ∼ |E|. Inside the metallic phase a (0) is finite. In Fig. 22(a), we show that within the range of disorder strength 0.50(weak) ≤ W ≤ 3.5(moderate) the TDOS also display the same scaling behavior as average DOS. This observation confirms that TDOS can also be subscribed as a bonafide order-parameter across the WSMmetal QPT. On the other hand, for strong enough disorder the TDOS t (E) decreases monotonically for any energy E, and ultimately t (0) becomes zero across the AT. Therefore, TDOS can serve as the order-parameter across all possible disorder-driven QPTs considered here.
Finally, we focus on the evolution of the location of the mobility edge in a dirty Weyl metal as a function of disorder strength by numerically computing the mobility edge, defined as In particular, the mobility edge defines the boundary between the extended and localized states, and we here focus on this quantity in the strong disorder regime W ≥ 2 > W c,1 . The results are shown in Fig. 21(b). For weak disorder the mobility edge resides at high-energy, indicating the metallic nature of a moderately dirty Weyl system. However, the mobility edge progressively slides down toward smaller energy with increasing randomness in the system. Finally, across the AT the mobility edge comes down to zero energy, indicating that all states inside the Anderson insulator are localized. Notice that the shape of the mobility edge is quite distinct in a Weyl metal in comparison to its counterpart in conventional metal [88], which however can solely be attributed to the linear dispersion of Weyl quasiparticles in the clean system.

X. SUMMARY AND DISCUSSION
In this paper we have studied the role of generic disorder in a Weyl semimetal, by considering its simplest realization, comprised of only two Weyl nodes. When the system resides in the proximity of semimetal-insulator quantum phase transition, the generalized Harris criterion suggests that such critical point is stable in the presence of weak but generic disorder. By contrast, a multicritical point appears in the phase diagram for strong disorder, where the Weyl semimetal, an insulator and a metallic phase meet. Within the framework of an appropriate -expansion we show that the critical exponents at such multicritical point are (i) dynamic scaling exponent z = 1 + n /2, and (ii) correlation length exponent ν = 1/ n that controls the relevance of disorder coupling, where n = 1/2 for physical system. These findings are in good agreement with the ones obtained numerically, yielding ν = 1.98 ± 0.05 and z = 1.23 ± 0.05.
On the other hand, when the system is deep inside the Weyl semimetal phase, we have shown that the continuous global chiral U (1) symmetry plays a fundamental rule in classifying the disorder-driven Weyl semimetal-metal quantum phase transitions. The simplest realization of a Weyl semimetal is susceptible to eight types of disorder, among which only four preserve such chiral symmetry. Using two different -expansions, we have shown that the chiral symmetric disorder driven semimetal-metal transition takes place through either a quantum critical point or a line of quantum critical points. Irrespective of details, the critical exponents are always z = 1 + /2 and ν = −1 , and = 1 corresponds to the physical situation. Such unique set of exponents in the presence of generic chiral symmetric disorder gives birth to an emergent chiral superuniversality across the Weyl semimetalmetal quantum phase transition. Furthermore, we have performed a thorough numerical analysis of average density of states in Weyl semimetals with chiral symmetric disorder. The emergence of chiral superuniversality has been demonstrated through numerical analysis of average density of states near zero energy. We show that for any such disorder Weyl semimetal undergoes a continuous quantum phase transition into a diffusive metallic phase. Within the numerical accuracy, we find that across this transition z ≈ 1.5 and ν ≈ 1, in excellent agreement with our field theoretic predictions obtained from leading order -expansions (see Table I for comparison). The quality as well as reliability of our numerical analysis has been anchored through two completely different types of high-quality data collapses, shown in Fig. 15, in the entire phase diagram of a dirty Weyl semimetal for all possible chiral disorder.
For chiral symmetry breaking disorder, the Weyl semimetal-metal quantum phase transition also takes place through a critical point or a line of critical points, but the critical exponents are significantly different from the ones reported in the presence of chiral disorder. Even though the critical exponents across such semimetalmetal transition turn out to be slightly dependent on the renormalization group scheme, we always find z > d and ν = 1/ from leading order -expansions. Consequently, all chiral symmetric or intra-node disorder (as well as higher gradient terms that are inevitably present in a lattice) become relevant at such putative line of critical points. As a result, inter-node disorder driven semimetal-metal phase transition is ultimately always governed by the chiral symmetric disorder, yielding ν ≈ 1 and z ≈ 3/2, characteristic of chiral superuniversality. We anchor these outcomes by numerically extracting the scaling of average density of states in the presence of inter-node disorder, and the results are shown in Table II and Figs. 4 and 16.
Even though we promoted such classification scheme in a Weyl semimetal with only two nodes, our prescription can easily be generalized to Weyl systems with multiple flavors, as well as topological Dirac semimetals with bonafide time-reversal symmetry that has recently been found in Cd 2 As 3 [89] and Na 3 Bi [90] and the ones at the quantum critical point residing between two topologically distinct insulating vacua.
We here mention that d expansion can be problematic beyond the leading order in d , since the contribution from diagrams (c) and (d) in Fig. 6 and their higherloop cousins are typically ultraviolet divergent and one looses the order by order control over the perturbative calculation [40,44]. Such a class of diagrams is, however, ultraviolet finite and thus does not contribute to renormalization group flow equations in the m expansion scheme. We, therefore, strongly believe that higher order perturbation theory can only be reliable in the m expansion scheme, and our results can serve as an ideal testbed to investigate the reliability of these two schemes.
In addition to the Weyl semimetal-metal quantum phase transition, we also establish that a compressible Weyl metal undergoes a a subsequent transition st stronger disorder into a Anderson insulator. We track the typical density of states to pin the onset of such insulating phase that only accommodates localized states. In particular, we show that across the Weyl metal-insulator transition the typical density of states at zero energy ( t (0)) smoothly vanishes, and thus serving as bonafide orderparameter, while the average density of states remains non-critical across this transition. In addition, we also find that t (0) remains pinned in the Weyl semimetal phase and becomes finite in the metallic phase. Therefore, typical density of states at zero energy serves as a unified order-parameter across all possible disorderdriven quantum phase transition in a Weyl semimetal.
Finally we comment on some non-perturbative effects of disorder in Weyl semimetals, such as puddles [91], Lishiftz tail [92], and rare-region states and Griffiths physics [34,60]. Puddles are inevitable in real materials as there are always density fluctuations that locally shift the chemical potential away from the Weyl nodes, while maintaining the overall charge neutrality of the system. In addition, presence of disorder can also support quasi-localized rare states at zero-energy even for subcritical strength of disorder. Although such effects are important and interesting, they do not affect the quantum critical behavior. This is so because rare and critical excitations appear to be decoupled from each other. Also a recent numerical work has demonstrated that such non-perturbative effects can be systematically suppressed with a suitable choice of the distribution of disorder [62]. Therefore, these effects do not alter any physical outcome we reported in this paper.
insulator for ∆ < 0 and (ii) WSM for ∆ > 0, with c representing the monopole charge of the Weyl nodes. Respectively for c = 1, 2 and 3, single, double and triple WSMs are realized in a crystalline environment [93][94][95]. The effective dimensionality (d * ) of such critical semimetallic phase can be found from the corresponding imaginary time Euclidean action where ψ is a two component spinor, describing the critical excitations residing at the WSM-insulator QCP. All parameters, such as α c and b, remain invariant under the rescaling of space-time(imaginary) co-ordinates according to τ → e l τ , (x, y) → e l/c (x, y), x 3 → e l/2 x 3 , when accompanied by the field normalization ψ → Z Therefore, only the single (c = 1) WSM-insulator QCP is stable against sufficiently weak mass/bond disorder. Furthermore, the stability of the WSM-insulator QCP in the presence of generic disorder, which appears similar to V z (x) in Eq. (8), can be established from the generalized Harris criterion (A3). Hence, a single WSM-insulator QCP is guaranteed to be stable against generic disorder. In this regard a comment is due. Our derivation of generalized Harris criterion differs from the original one in Ref. [27], where d * is replaced by the physical dimensionality of the system (d) and the CLE ν varies depending on the nature of the phase transition. On the other hand, within the framework of anisotropic scaling of spatial co-ordinates we always find ν = 1, but actual spatial dimension gets replaced by an effective dimensionality of the system (d * ) under the process of coarse graining. We believe that these two methods are complimentary to each other.
(B2) Therefore, as n → ∞ contribution only from h 1 (n) survives and for any finite n, h 2 (n) and h 3 (n) give rise to subleading divergences. The RG flow equations obtained by keeping only the leading divergence are shown in Eq. (11) of the main text. In order to reliably extract the critical exponents, in particular the DSE z at the multicritical point shown as a red dot in Fig.7, we should only keep the leading order (n → ∞) contribution. At least to the leading order in n expansion, inclusion of subleading divergences [arising from h 2 (n) and h 3 (n)] does not alter the value of CLE, and we always find ν −1 = n = 1/2.

2.
d expansion about lower critical dimension d l = 5 2 In this section we demonstrate the role of disorder in the vicinity of WSM-insulator QPT perturbatively using an d expansion near the lower critical dimension d l = 5/2 in the theory, where d = d − 5/2. As we will see the outcomes are qualitatively the same as in the n regularization scheme. The exact values of the critical exponents are, however, different from the ones announced in Sec. III, although only slightly so, at least to the one-loop order. Upon integrating the fast modes within the shell E c e −l < v 2 k 2 ⊥ + b 2 k 4 z < E c , where E c is the ultraviolet energy cutoff for critical excitations residing the WSM-insulator QCP, we arrive at the following flow equations to the leading order in d expansion for X = v, b, after defining the dimensionless disorder coupling constant as ∆ j α → ∆ j for j = 0, ⊥, z, where α = E d c / 20π 2 v 2 b 1/2 and ∆/E c → ∆. Then, β−function for v and b in the presence of disorder yields a scale dependent dynamic scaling exponent z(l) = 1 + 5 [∆ 0 + 2∆ ⊥ + ∆ 3 ] (l). (B4) The coupled RG flow equations from Eq. (B3) also support only two fixed points: (i) (∆, ∆ 0 , ∆ ⊥ , ∆ 3 ) = (0, 0, 0, 0), representing the WSM-insulator QCP in the clean limit (the blue dot in which is extremely close to the ones reported in Sec. III for d = 1/2. Therefore, both methods produce qualitatively similar results near WSM-insulator QPT, and the obtained critical exponents for anisotropic semimetalmetal transition are extremely close to each other, at least to the leading order. The resulting phase diagram is shown in Fig. 23. In this appendix we display the detailed analysis of various one-loop diagrams, shown in Fig. 10, within the framework of an m expansion.

Self-energy
Let us first consider the self energy diagram in Fig. 10(a). The expression for the self-energy reads with d = 3, the summation is taken over all eight types of disorder (see Table III) and q ≡ |q|. The contribution from one-loop self-energy diagram from the disorder represented by the matrix N reads We will evaluate the temporal and spatial components of the self-energy diagram separately. Let us first set k = 0, for which , where x is the Feynman parameter. Upon completing the integrals over q and x, and setting m = 1 − (for brevity, we use here shorthand notation m → ) we obtain Next we set ω = 0 and the spatial component of selfenergy correction is then given by Hence, the total self energy correction reads where we have redefined ∆ N k /(2π 2 v 2 ) → ∆ N , which is Eq. (23) in the main text.

Vertex
The vertex correction for the disorder vertex shown in Fig. 10(b) with the matrix N reads where we kept only one external momentum as an infrared regulator. The last expression can be compactly written as where We now present the evaluation of the above integral After shifting the momentum variable as q − xk → q, we obtain after taking m = 1 − , since only the q−dependent part in the numerator of the integrand yields a divergent contribution. We use the last expression to obtain Eq. (30) in the main text.
The above set of coupled flow equations supports only a line of QCPs, given by Eq. (39), in the ∆ V − ∆ A plane, shown in Fig. 11. Along the entire line of QPCs, the exponents are ν −1 = d and z = 1 + d /2 (to the leading order in d ). Therefore, in a three-dimensional WSM ( d = 1) the semimetal-metal QPT driven by arbitrary disorder potential is always characterized by ν = 1 and z = 3/2, thus strongly supporting the proposed emergent chiral superuniversality.
Appendix F: Alternative derivation of correction to optical conductivity Direct computation of the correction to the OC due to arbitrary disorder by using the Kubo formula has already been presented in Ref. [42]. Specifically, we compute the disorder driven correction to the current-current correlation function (involving computation of two-loop diagrams) and then via analytic continuation we found the OC at frequency Ω in a weakly disordered WSM to be σ(Ω) = N e 2 0 Ω 12hv 1 + ∆ V Λ π 2 v 2 ≡ σ 0 (Ω) 1 + where N is the number of Weyl nodes, e 0 is the electron charge in vacuum [see Eq.
(3) of Ref. [42]]. For con-creteness, we here restrict ourselves to potential disorder or random charge impurities (∆ V ), possessing Gaussian white noise distribution in three dimensions. In the absence of disorder (∆ V = 0), we recover the OC in a clean WSM, σ 0 (Ω), [30,32,96]. We now present an alternative derivation of the same expression. The OC is given by where Z Ψ = 1 + ∆ V Λ/(2π 2 v 2 ) is the field renormaliza-tion factor, as presented in Sec. IV C, for d = 1. The same expression for the field-renormalization factor can directly be obtained by integrating over the entire Weylband with 0 ≤ |k| ≤ Λ, which is legitimate since we are interested in the OC of a weakly disordered WSM for which sharp quasiparticle excitations persists all the way down to zero energy or momentum. Upon substituting Z ψ in the above expression we immediately recover Eq. (F1).
where δ = W − W c /(W W c ) measures the reduced disorder strength from the critical one (W = W c ). After regularizing the ultraviolet divergence we can take the limit E Λ → ∞ without encountering any divergence. The self-consistent solution of the scattering life-time is then obtained from the following universal scaling form which immediately implies that τ −1 s is finite only when δ > 0 or W > W c , and for W < W c we get τ −1 s = 0. Therefore, critical fermions separating a WSM and an insulator retain its ballistic nature upto a critical strength of disorder W c ∼ E 1/2 Λ . Only for strong disorder W > W c a metallic phase emerges where τ −1 s is finite. Therefore, conclusion from self-consistent Born approximation is in qualitative agreement with our results found by field theoretic RG analysis and numerical calculation, presented in Sec. III.