Anomaly Non-renormalization in Interacting Weyl Semimetals

Weyl semimetals are 3D condensed matter systems characterized by a degenerate Fermi surface, consisting of a pair of ‘Weyl nodes’. Correspondingly, in the infrared limit, these systems behave effectively as Weyl fermions in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$3+1$$\end{document}3+1 dimensions. We consider a class of interacting 3D lattice models for Weyl semimetals and prove that the quadratic response of the quasi-particle flow between the Weyl nodes is universal, that is, independent of the interaction strength and form. Universality is the counterpart of the Adler–Bardeen non-renormalization property of the chiral anomaly for the infrared emergent description, which is proved here in the presence of a lattice and at a non-perturbative level. Our proof relies on constructive bounds for the Euclidean ground state correlations combined with lattice Ward Identities, and it is valid arbitrarily close to the critical point where the Weyl points merge and the relativistic description breaks down.


Introduction
The behaviour of relativistic fermionic particles is described by the Dirac equation and by its interacting extensions, such as Quantum ElectroDynamics (QED) or other standard Quantum Field Theory (QFT) models, whose spectacular and often counter-intuitive predictions have been tested for decades in high energy experiments. On the other hand, the electrons present in ordinary matter, like the conduction electrons in metals, which live at much lower energy scales, are described by the many-body Schrödinger equation, which breaks Lorentz invariance explicitly.
By scale separation, the Dirac high energy physics is not expected to play directly any role in the behaviour of conduction electrons. Nevertheless, there are cases in which the interaction of the conduction electrons with the underlying lattice can produce an effective description in terms of Dirac particles, with a velocity of propagation much smaller than the speed of light. Such emergent QFT description is valid in great generality for one dimensional metals, see [48,58,67] (see also [55] for a recent review), but is also possible in higher dimensional systems: notable examples were proposed in [37,68] and in [63]. The discovery of graphene [44,64] and of topological insulators [38] provides a condensed matter realization of systems of 2D Dirac fermions. In connection with their emergent QFT description, the transport coefficients of these systems are characterized by remarkable universality properties: examples include the optical conductivity of graphene [60] and the Hall conductance [11,66]. A mathematical proof of universality in the presence of interactions has been established in [6,12,13,26,[28][29][30][31]39,57]. More recently, 3D Dirac systems have been experimentally realized [16,46,47,62], following the theoretical predictions of [63] and [24,59,65], see [9,38] for reviews; these semimetals are dubbed 'Dirac' or 'Weyl' semimetals, depending on whether the Fermi points coincide or are at distinct locations in the Brillouin zone. Such experimental discoveries open the way to the observation of the counterpart of one of the most interesting phenomena characterizing 3D massless Dirac fermions: the Adler-Bell-Jackiw (ABJ) axial anomaly [2,14], which corresponds to the breaking of the classical chiral symmetry by quantum mechanical radiative corrections and results in the non-conservation of the axial current. Nielsen and Ninomiya [63] predicted that this quantum anomaly should have a realization, or 'simulation', in the quadratic response of the quasi-particle flow between the Weyl nodes to an external electromagnetic field, see also [18,69]; in particular, this quasi-particle flow is expected to be proportional to E · B, with coefficient given by the value of the axial anomaly in QED 4 . The experimental observation of the 'chiral anomaly' in Weyl semimetals has been reported in [41,43,45,71,72], see also [15] for an observation in the context of superfluid 3 He. Other applications of the chiral anomaly in condensed matter physics have been discussed in [23], in the context of a gauge theory of phases of matter, see also [22] for new proposed applications and [21] for a recent review.
The theoretical analysis in [63] neglects the effects of many-body interaction. What happens to the quadratic response of the quasi-particle flow in the presence of electron interactions, which are unavoidable in real systems? One of the distinctive features of the axial anomaly in QFT is the Adler-Bardeen (AB) non-renormalization property [3], which says that all the radiative corrections cancel out at all orders in the fine structure constant, so that the axial anomaly equals the value of the chiral triangle graph. However, the AB argument cannot be applied to Weyl semimetals: the analysis in [3] requires cancellations between classes of Feynman graphs based on Lorentz and chiral symmetries, which are broken explicitly in any lattice model of interacting electrons. Strangely enough, the issue of universality has not been addressed theoretically until now, despite the huge literature on the effects of interactions 1 [10,20,35,36,49,51,70]. In this paper, we consider a paradigmatic class of lattice models for 3D Weyl semimetals with short-range interactions, generalizing a widely used 'standard' model for these materials, see [9]. For this class of models, we prove that the quadratic response of the quasi-particle flow between the Weyl points, which is the counterpart of the chiral anomaly in the emergent QFT description of the system, is universal, that is, independent of the interaction strength and form.
Therefore, this quantity can be added to the limited list of universal transport coefficients in condensed matter, like the Hall conductivity and graphene's optical conductiv-ity, and is the only known example of universal quantum transport coefficient in three dimensions. This quadratic response has been measured [41,43,45,71,72] in different Weyl semimetal, such as Na 3 Bi, TaAs, GdPtBi, ZrTe 5 , and the proportionality between the quasi-particle flow between the Weyl points and E · B has been experimentally verified. The interaction does not appear to qualitatively modify such a response, in agreement with our result, even though a precise measurement of the proportionality coefficient remains to be performed. Remarkably, our result is valid arbitrarily close to the critical point where the Weyl nodes merge and the natural parameter measuring the strength of the interaction, that is, the coupling strength divided by the quasi-particle velocity, diverges: this is important for applications to real 3D Weyl semimetals, where typically the distance between the Weyl nodes is very small, as compared to the size of the Brillouin zone; see Sect. 2 for further discussion of this point.
In this paper, we focus our attention on the quadratic response coefficient of the 'chiral current' to an external electromagnetic field directly, without discussing the linear response coefficient, which is not expected to have any remarkable universality feature. We stress that in the experiments the 'chiral anomaly' is detected by measuring the peak of quasi-particle flow between the Weyl nodes as the angle between the external electric and magnetic fields E and B is varied, see e.g. [71]: such a procedure has the effect of automatically subtracting the linear response, which is uninteresting for the effect under investigation.
The paper is organized as follows: in Sect. 2, we introduce a class of lattice Weyl semimetals, we define the notions of lattice current and lattice 'chiral current', and state our main results, summarized in Theorem 2.1. In Sect. 3 we present the proof: first we describe the general strategy and the core argument (Sect. 3.1); next we discuss the technical details: the Grassmann representation of the generating function of correlations and its multiscale Renormalization Group (RG) computation are presented in Sects. 3.2 and 3.3; the RG flow of the effective couplings is studied in Sect. 3.4; the regularity properties of the Euclidean correlations is discussed in Sect. 3.5, and the relativistic contribution to the quadratic response of the chiral current to the electromagnetic field is computed in Sect. 3.6. A few additional technical aspects of the proof are deferred to the appendices: in Appendix A we prove a number of key symmetry properties of the kernels of the effective potentials produced by the RG iterations; in Appendix B we discuss a few technical aspects of the tree expansion for the effective potentials; in Appendix C we compute the quadratic response coefficient via time-dependent perturbation theory and prove a Wick rotation theorem that allows us to express it in terms of Euclidean correlations.

A basic non-interacting model for Weyl semimetals.
We start by describing a paradigmatic non-interacting lattice model for Weyl semimetal, see [9, Section II.B.1], with broken time-reversal symmetry, for which inversion symmetry is preserved. This example motivates the class of models which our analysis applies to. It describes an electron system on a simple cubic lattice subject to an external magnetic field and with the following features: the system is magnetically ordered, so that the bands have no spin degeneracy; moreover, every lattice site comes with two internal degrees of freedom, r = 1, 2, playing the role of pseudo-spin labels, corresponding to, e.g., orbital degrees of freedom; finally, these orbitals (say, s, p orbitals) have opposite parity. i.e., they behave differently under inversion. The model is defined in terms of the following Bloch Hamiltonian (see next section for details):Ĥ 0 (k) = α(k) β(k) β * (k) −α(k) , (2.1) where α(k) = 2 + ζ − cos k 1 − cos k 2 − cos k 3 , β(k) = t 1 sin k 1 − it 2 sin k 2 (2.2) and ζ, t 1 , t 2 are three real parameters. We assume that t 1 , t 2 > 0 and −1 < ζ < 1, so thatĤ 0 (k) is singular only at two Weyl nodes, or Fermi points, k = p ± F , with p ± F = (0, 0, ± arccos ζ ). We shall think of t 1 , t 2 as being fixed once and for all, while ζ is tunable: it controls the distance among the Weyl nodes, which tend to merge in the limit as |ζ | → 1 − . Note that the relative distance between the two nodes vanishes like √ 1 − |ζ | as |ζ | → 1 − .
Remark. Having a parameter tuning the distance between Weyl nodes is extremely natural: in real systems, magnetic fields in orthogonal directions are used to induce and control a 'topological' transitions from, e.g., a Dirac semimetal to a magnetic semiconductor, or to a Weyl semimetal; once a Weyl semimetal phase is entered, the intensity of the magnetic fields is used to tune the location of the Weyl nodes. See [9,Fig.2] and the corresponding discussion; see also [20]. In the limit as the Weyl nodes tend to merge, the valence and conduction bands generically tend to touch quadratically at the doubly-degenerate node. Experimentally, the Weyl nodes may be quite close to each other, compared to the size of the Brillouin zone. The fact that the Weyl points can be arbitrarily close is a peculiarity of Weyl semimetals, which distinguish them from other Dirac materials, like graphene.
In the vicinity of the Weyl nodes, k p ω F , the Bloch Hamiltonian can be linearized as:Ĥ 0 (k) = t 1 σ 1 k 1 + t 2 σ 2 k 2 + ω 1 − ζ 2 σ 3 (k 3 with σ i the Pauli matrices: The parameters t 1 , t 2 , ω 1 − ζ 2 in front of the three terms in the right side of (2.3) play the role of 'Fermi velocities' in the three coordinate directions. Note that, as |ζ | → 1, the third component of the Fermi velocity (the one along which the Weyl nodes are located), v 0 3,ω := ω 1 − ζ 2 , tends to zero at the same speed as p + F − p − F . This is no accident, it is a generic feature in the limit as two Weyl nodes merge.
In order to provide an asymptotic description of the system that holds unifromly in the limit as the Weyl nodes merge, it is convenient to explicitly extract the quadratic part of the Bloch Hamiltonian in the third coordinate direction. In our case: 3 ) 2 +O(k 2 1 , k 2 2 , |k 3 − p ω F,3 | 3 ). (2.5) Note that the prefactor of the quadratic term that we extracted, equal to ζ /2, tends to a non-zero constant in the limit as the Weyl nodes tend to merge, |ζ | → 1 − . This is again no accident, it is another generic feature in the limit as two Weyl nodes merge.

A class of non-interacting models for Weyl semimetals.
Motivated by the model introduced in the previous section, we now consider a class of non-interacting lattice models, including (2.1) as a special case. We will use adapted coordinates, such that, as in the model above, the Weyl nodes are located along the third coordinate direction. Let L be a cubic lattice of side L, with periodic boundary conditions; that is L = Z 3 L , with Z L = Z/LZ the integers modulo L. The coordinates of the particle on the lattice are specified by (x, r ), with x ∈ L the position and r = 1, 2 the internal degree of freedom, called 'color', or 'orbital'. Thus, the single-particle Hilbert space is h L = C L 3 ⊗ C 2 . We denote by H 0 the single-particle Hamiltonian on h L , which we assume to be Hermitian. We consider finite-ranged, translation invariant systems, such that H 0 (x, y) ≡ H 0 (y − x). We always assume that L is larger than the range of H 0 . We also assume that, after identifying y − x with its image in Z 3 , H 0 (y − x) is independent of L (i.e., the kernel of H 0 is independent of L, up to the periodicity condition). Due to periodicity and translation invariance, we can represent the single-particle Hamiltonian as: where B L is the finite-volume Brillouin zone, defined as: andĤ 0 (k) is the Bloch Hamiltonian, a two-by-two matrix parametrized by k ∈ B L . Under the above assumptions, the Bloch Hamiltonian is an L-independent function, Hermitian and analytic in k ∈ B ∞ ≡ T 3 . We also allow forĤ 0 to depend smoothly on an additional paremeter ζ ∈ [ζ 0 , ζ 1 ) ≡ I , which can be used to tune the location of the Weyl points, see below.
We assume the following symmetry properties, valid for all ζ ∈ I , which are all satisfied by the model in Eq. (2.1). The first two symmetries are the most physically significant: they characterize the class of Weyl semimetals under consideration. The other two are not crucial, but simplify the structure of the linearized Hamiltonian around the Fermi points.
In addition to these symmetries, we assume the following properties on the energy bands of the Bloch Hamiltonian. 2 We denote by z * , rather than byz, the complex conjugate of z ∈ C. Given a matrix A ∈ C n×n , we denote by A † its Hermitian conjugate, and by A * its complex conjugate, that is, the matrix whose elements satisfy (A * ) i j = (A i j ) * . 3 As mentioned in the previous section, we are assuming that the two orbitals r = 1 and r = 2 transform differently under inversion, that is, one is even and one is odd under x → −x, which explains the presence of the third Pauli matrix σ 3 .

Weyl points. We write the Bloch Hamiltonian in the form
(2.8) Using the symmetries (ii), (iii), (iv), we see that: a is odd in k 1 , even in k 2 and in k 3 ; b is odd in k 2 , even in k 1 and in k 3 ; c is even in k 1 , in k 2 and in k 3 ; d is odd in k 1 and in k 2 , even in k 3 . We let ε ± (k) be the k-dependent eigenvalues of the Bloch Hamiltonian, also known as the 'energy bands': For any ζ ∈ I , we assume that ε + (k) ≥ 0 ≥ ε − (k), with equalities only at two points, called the Weyl nodes, denoted by p ± F , at which we assume ε ± (k) to vanish linearly. Thanks to the symmetries (ii), (iii), (iv), the pair of Weyl nodes is invariant under reflections about any of the coordinate axes. Therefore, in order for the nodes to be exactly two, they need to be located along one of the three axes. Moreover, the requirement that ε ± (k) vanish linearly at p ± F , combined with the fact that a, b, c, d are all even functions of k 3 , implies that the Weyl nodes are located along the third axis: We allow for the possibility that the Weyl nodes merge in the limit ζ → ζ − 1 . For definiteness, we assume once and for all that If the limit were different from zero, the model could be treated in the same way as if ζ ∈ [ζ 0 , ζ 1 ] ⊂ [ζ 0 , ζ 1 ), which is a sub-case of the problem studied in this paper. 4 Let us now discuss the structure of the Bloch Hamiltonian in the vicinity of the Weyl nodes. We have a( p ± If we Taylor expand the Bloch Hamiltonian around p ω F , with ω ∈ {±}, at first order in k 1 , k 2 and at second order in k 3 − p ω F,3 , recalling the parity properties of a, b, c, d, we obtain the analogue of (2.5), namelyĤ where v 0 1 , v 0 2 , v 0 3 , b 0 are real variables, smoothly depending on ζ , and a R , b R , c R,ω , d R are the remainders of the Taylor expansion (note that d R = d). We assume that there exists a constant c 0 ≥ 1 such that uniformly in ζ , for ζ ∈ I . Thanks to their parity properties of a, b, c, d, the remainder terms satisfy the following bounds in the vicinity of p ω F : for a constant C that we assume to be independent of ζ . By the 'vicinity of p ω F ', we mean |k − p ω F | ≤ 2| p + F |. In the complementary region, min ω |k − p ω F | > 2| p + F |, we have: for a constant C that, without loss of generality, can be taken to be the same as in (2.14). Note that, in view of the bounds (2.13), (2.14), (2.15), it is possible to find c 1 > 0 such that (2.16) uniformly in ζ ∈ I . Moreover, by the smoothness ofĤ 0 (k) and the fact that detĤ 0 (k) vanishes only at π ± F , there exists c 2 > 0 such that uniformly in ζ ∈ I . Under these assumptions, we will be able to construct and analyze the ground state of a many-body interacting version of H 0 , uniformly in ζ , for ζ ∈ I . In particular, our analysis will be valid in the limit ζ → ζ − 1 as the Weyl nodes merge. For other examples of models in this class, in addition to the one discussed in the previous section, see [19,42].

The interacting model.
Let us now include the many-body interaction. We describe the system in the grand-canonical setting, in second quantization. The fermionic Fock space is defined as The many-body Hamiltonian is: where ρ x,r = a + x,r a − x,r is the density operator and N L = x∈ L (ρ x,1 − ρ x,2 ) is the staggered number operator. The first term in the right side of (2.19) is the hopping term, defined in terms of an H 0 satisfying the properties listed in the previous section. The second term is the many-body interaction, defined in terms of a w r,r (x, y), which we assume to be even, short-ranged, translational invariant, w r,r (x, y) = w r,r (y − x), and periodic over L . As for H 0 (x, y), we assume that w(x, y) is L-independent, up to the periodicity condition. Also, we suppose that the interaction is invariant under the reflections (iv) and (v) above: that is, we require that w r,r (x) is invariant under . The interaction strength λ will be assumed to be small, compared with the bandwidth max k H 0 (k) , uniformly in the system size and in the choice of the parameter ζ . The −1/2 appearing in the factors (ρ x,r − 1 2 ) and (ρ y,r − 1 2 ) correspond to a specific, convenient, choice of the chemical potential. Finally, the third term in the right side of (2.19) is a staggered chemical potential, with ν ≡ ν(λ) being a free parameter, such that ν(0) = 0, to be chosen in such a way that the Fermi points of the interacting theory do not move as the interaction varies: they will coincide with those of the noninteracting theory, p ± F , in a sense to be made precise later.
The Gibbs state of the model is: Given a local observable O x , even in the fermionic fields, we denote by We denote byÔ p its space-time Fourier transform, β Z is the set of bosonic Matsubara frequencies. If, instead, O x is an operator which is odd in the fermionic fields, such as a − x , similar definitions and formulas hold; however, in the definition of Fourier transform, the set of bosonic Matsubara frequencies is replaced by the set of fermionic Matsubara frequencies, M f β = 2π β (Z + 1 2 ). Given k ∈ M f β × B L , we denote the space-time Fourier transforms of the fermionic operators as: Later, we shall be interested in the Schwinger correlation functions of the model, defined, for x 0,i ∈ [0, β), as: with T the fermionic time-ordering, ordering the operators in decreasing imaginary-time order, see, e.g., [30,Eq. (3.3)]. We extend (2.23) anti-periodically to all imaginary times, x 0,i ∈ R. Of particular interest is the two-point Schwinger function, S β,L 2;r 1 ,r 2 (x, y) = Ta − x,r 1 a + y,r 2 β,L , together with its thermodynamic and zero temperature limit 6 Ta − x,r 1 a + y,r 2 β,L , (2.24) which is translationally invariant, whenever it exists. We also denote byŜ 2;r 1 ,r 2 (p) its Fourier transform.
In the absence of interaction, λ = ν = 0, the two-point function (thought of as a 2 × 2 matrix of elements S 2;r 1 ,r 2 ) can be written, for x 0 − y 0 / ∈ βZ, as: y). (2.25) We refer to g β,L as the free propagator of the model, we set g(·) = lim β,L→∞ g β,L (·), and we denote byĝ its Fourier transform. Let p ω F := (0, p ω F ). In the absence of interactions, as β, L → ∞, the Fourier transform of the free propagator reads, for k = k + p ω F and k small,ĝ We see that g ω (k ) := g(k +p ω F ) agrees at leading order with the propagator of chiral relativistic fermions, with chirality ω = ±, and anisotropic velocities. We will prove that, for λ small, by fixing ν = ν(λ) appropriately, the singularities of the Fourier transform of the interacting propagator are located at the same points, p ± F , as the non-interacting one; moreover, at those points, it has the same singularity structure: with v j = v 0 j (1 + O(λ)), j = 1, 2, 3, the interacting Fermi velocities, Z = 1 + O(λ) the wave function renormalization and, for any θ ∈ (0, 1), the remainder satisfies R ω (k ) = O(|k | θ ), non-uniformly in θ as θ → 1 − (and, possibly, non-uniformly in the distance between the Weyl nodes).

2.4.
Coupling to an external gauge field. Our analysis will focus on the transport properties of the model, after introducing an external electromagnetic field. The coupling of the model with an external vector potential is defined via the Peierls' substitution. This means that, both in the Hamiltonian and in the physical observables, any product of fermionic operators a + x,r a − y,r is replaced by: where A x ∈ R 3 , with x ∈ Q L := (R/LZ) 3 , and x→y A · d denotes the line integral: We use the following convention for the Fourier transform of A: if p ∈ 2π L Z 3 , Notice that the many-body interaction is not affected by the presence of the gauge field, due to the fact that the density operator is gauge invariant. Let us denote by H 0 L (A) the free Hamiltonian (λ = ν = 0) in the presence of the gauge field. The charge current operator is defined as: We also setĴ j, p ≡Ĵ j, p (0). A straightforward computation gives, letting η p (δ) := whereĴ j (k, p) should be understood as a 2 × 2 matrix, and a + k (resp. a − k ) as a row (resp. column) vector of components a + k,1 , a + k,2 (resp. a − k,1 , a − k,2 ). Note that, in the thermodynamic limit, the kernelĴ j (k, p) can be written more explicitly aŝ In particular, recalling (2.12), the kernel of the current at k = p ω F and p = 0 reads: (2.33) Equivalently, in real-space, where J j (y, z) = L −6 k∈B L p∈ 2π L Z 3Ĵ j (k, p)e ik·y e i p·z . In the thermodynamic limit, where δ(z − sy) is a Dirac delta. Note that the current satisfies a lattice continuity equation. In fact: which can be rewritten as We collect the density and the components of the lattice current to form a Euclidean 4-current (−iρ p ,Ĵ p ), whose components are denoted byĴ μ, p , μ ∈ {0, 1, 2, 3}:Ĵ whereĴ 0 = −i1 2 andĴ j , j = 1, 2, 3, are given by (2.31).

The lattice chiral current.
We now introduce a lattice current for the quasi-particle flow between the Weyl nodes, playing the role of the 'chiral current' for our lattice model: where p ∈ B L and, letting σ 0 := −i1 2 , (2.40) Note that the kernel Z 5 μ,bareĴ 5 μ (k, p) of the chiral current at k = p ω F and p = 0 reads: Comparing with (2.33), we see that at the Weyl nodes the different components of the chiral current behave like those of the total current, up to an additional 'chirality sign' ω = ± and different multiplicative factors. This implies that, in the 'infrared regime' of p small and k close to the Weyl nodes, the 'chiral density'Ĵ 5 0, p defined by (2.39) with μ = 0 reduces to the difference between the quasi-particle densities around the Weyl nodes, up to a multiplicative pre-factor; similarly, the spatial components of the chiral current reduce to the difference between the quasi-particle currents around the Weyl nodes. Making a parallel with QED 4 , this infrared chiral current coincides with the standard QED 4 axial current with different velocities.
The multiplicative normalizations Z 5 μ,bare are fixed as follows. Let us define the vertex functionsˆ where we recall thatŜ 2 (k) is the interacting propagator, see (2.24) and (2.27). Existence of the limits in (2.42)-(2.43), in the sense of footnote 6 above, is part of the main results of this paper, stated in the following section. We will impose that the amputated vertex functions agree in the infrared limit, up to a sign: that is, if p ω F is the 4-vector where, for technical convenience, we take the limit in such a way that k − p ω F , p, k + p − p ω F are all of the same order as they tend to zero. The normalization condition (2.46) can be interpreted as the requirement that, at the level of the amputated vertex correlation functions, for momenta close to a Weyl node, the insertion of the operator J 5 μ,p is equivalent to the insertion of ωĴ μ,p , where ω is the chirality of the Weyl node. Finally, we denote by J 5 μ (A) the chiral current coupled to the external gauge field, defined again via the Peierls substitution (2.28). We have: is the inverse Fourier transform ofĴ 5 μ (k, p).
Remark. 1) The limiting values ofγ μ (k, p) and ofγ 5 μ (k, p) as k → p ω F and p → 0 divided by Z v μ , with Z the wave function renormalization and v μ the dressed velocity, 7 see (2.27), have the physical meaning of (dressed) charges of the quasi-particles around the Weyl nodes. Gauge invariance ensures that the charge of quasi-particles associated with the electromagnetic current is not renormalized by the interaction, that is, its dressed and bare values coincide, see (3.4) below. In contrast, the charge associated with the chiral current is renormalized non-trivially: this should come as no surprise, because the model is not invariant under chiral gauge transformations and, therefore, there is no symmetry protecting such a charge from being dressed by the interaction. The choice of Z 5 μ,bare is used to impose, via (2.46), the correct physical normalization of the chiral current, namely the one guaranteeing that the charges of the quasi-particles around each Weyl node computed either via the electromagnetic or the chiral current are the same. An analogous phenomenon takes place in QED 4 , where it is well known that the axial vertex renormalization is different from the vectorial one, see e.g. [1, eq.(2)]: the discrepancy between the values of the two vertex renormalizations is a manifestation of the chiral anomaly, since a formal use of chiral gauge invariance would naively suggest their identification.
2) The definition of the components of the chiral current is largely arbitrary, as long as their kernels have the right asymptotic form at the Weyl nodes, see (2.41), and the right discrete symmetries, discussed in Appendix A.1. Changing the specific definition of the chiral current would only affect the specific value of the parameters Z 5 μ,bare . In view of the universality result we prove, this arbitrariness does not affect the quadratic response at dominant order.
3) The introduction of a lattice current for Weyl semimetals, generalizing the proposal of [63] (see also [40]) to lattice interacting models, is an original contribution of this paper. Despite the formal similarity between the electromagnetic and the chiral lattice current, it is important to highlight some basic differences. In the infrared limit ( p small and k close to the Weyl nodes), they reduce to the electromagnetic and axial currents of QED 4 , respectively. In QED 4 , both these currents are conserved, and their conservation is associated with basic gauge symmetries of the model (total and chiral, respectively). On the contrary, in our lattice realization, only the electromagnetic current is associated with a gauge symmetry, from which exact Ward Identities between correlations follow: these imply, in particular, that the dressed charge is not renormalized, see (3.4) below, and that the 4-divergence of the current vanishes even in the presence of an external electromagnetic field, see, e.g., (3.11) below. Neither of these properties holds for the lattice chiral current.
2.6. Main result: condensed matter simulation of the chiral anomaly. We are interested in the response of the expectation of the chiral 4-currentĴ 5 μ, p to an adiabatic external gauge field of the form A x (t) = e ηt A x , where η > 0 plays the role of adiabatic parameter. We denote by H L (A(t)) the Hamiltonian of the interacting system, coupled to the external gauge field via the Peierls substitution (2.28). We will consider the timeevolution of the many-body system from t = −∞ to t = 0.
Let · β,L;t ≡ Tr · ρ(t) be the time-dependent state of the system, where ρ(t) is the solution of the von Neumann equation i∂ t ρ(t) = [H L (A(t)), ρ(t)], with boundary condition ρ(−∞) = ρ β,L . Let Ĵ 5 μ, p (A(t)) β,L;t be the time-dependent average of the chiral lattice current. As discussed in the introduction, we focus our attention on the quadratic variation of this quantity with respect to the external field, denoted by Ĵ 5 μ, p (A(t)) (2) β,L;t : with p i = (η, p i ). The quadratic responseˆ 5;β,L μ,i, j (p 1 , p 2 ) is the analogue, in our condensed matter context, of the chiral anomaly in QED 4 , see (2.51) and following discussion for more details. It can be expressed in terms of equilibrium correlation functions, via second-order time-dependent perturbation theory; see Appendix C, Eq. (C.12). A rigorous application of the Wick rotation, proven in Appendix C, allows us to express the β, L → ∞ limit of this quantity in terms of Euclidean correlation functions. Letˆ 5 μ,i, j (p 1 , p 2 ) = lim β,L→∞ˆ 5;β,L μ,i, j (p 1 , p 2 ). For p 1 = (η, p 1 ), p 2 = (η, p 2 ), where · ∞ = lim β,L→∞ (β L 3 ) −1 · β,L . We refer the reader to Eq. (C. 15) for the precise form of the Schwinger terms. More generally, we denote byˆ 5 μ,ν,σ (p 1 , p 2 ) the extension of (2.49) to space-time current insertions, the labels ν, σ = 0 corresponding to insertions of a density operatorĴ 0,p i . The next theorem gives the explicit expression of the interacting quadratic response of the chiral current,ˆ 5 μ,ν,σ (p 1 , p 2 ), in the limit of low momenta.
The proof of this theorem, presented in the rest of this paper, provides a constructive algorithm for computing the functions ν, Z 5 μ , v j , as well as the β, L → ∞ limit of the Euclidean correlation functions of the system. The construction and proof of analyticity of the staggered chemical potential ν and of the Fermi velocities v j , j = 1, 2, 3, uniformly in ζ , is not new, see [51,52]. An explicit computation at the level of first order perturbation theory shows that the interacting Fermi velocities v j , j = 1, 2, 3, are non-trivial functions of λ and of all the other parameters entering the definition of the Hamiltonian, generically in the choice of H 0 and w. In this sense, the Fermi velocities are non-universal quantities: their value depends on the details of the microscopic Hamiltonian. The same is true for the longitudinal Kubo conductivity, see [51], and for several other physical observables.
The important new piece of information contained in Theorem 2.1 is Eq. (2.50), which shows that the quadratic response of the chiral current is universal in the low momentum limit. It plays the same role as the quantized transverse conductivity in quantum Hall fluids. Remarkably, (2.50) is valid for all the choices of the parameter ζ controlling the distance among the Weyl points; it is valid, in particular, arbitrarily close to the critical point ζ = ζ 1 where the Weyl points merge and the valence and conduction bands touch quadratically, rather than linearly. The situation where the Weyl nodes are very close and the quasi-particle velocities are very small is very common in the actual experimental realization of Weyl semimetals, see e.g. [9,Figs. 2 and 3] and [20, Fig. 1]. In such a situation, the natural parameter measuring the strength of the interaction is the coupling strength divided by the quasi-particle velocity; surprisingly, even if this parameter becomes huge, the interaction still appears not to affect the anomaly coefficient.
Note that the universality of the 'chiral anomaly' in (2.50) holds, provided that the lattice chiral current verifies the condition (2.46). The subtle cancellation underlying our universality result can also be stated in an alternative, yet equivalent, way: let the chiral current be defined as in (2.39), but without apriori fixing Z 5 μ,bare in any special way; then the divergence of the quadratic response coefficientˆ 5 μ,ν,σ divided by the ratio of the axial and vector vertex renormalizations (i.e., by the left side of (2.46) times the chirality index ω) is universal at leading order. The fact that the quadratic response of the axial current to the electromagnetic field divided by the axial renormalization is universal (and not the quadratic response itself, as often stated incorrectly) is the very content of the Adler-Bardeen theorem in massless QED 4 , see the discussion in [1] after Eqs. (6) and (7). Indeed, a naive perturbative computation of the chiral anomaly in QED 4 would apparently give rise to higher order corrections beyond the chiral triangle graph, see [5]: however, as discussed in [1], these corrections are cancelled exactly by the axial vertex renormalization, which is different from the vectorial one.
Note also that, of course, the quadratic response of the electromagnetic, rather than axial, current vanishes, by current conservation.
In order to make the connection between the quadratic responseˆ 5 μ,ν,σ (p 1 , p 2 ) and the transport experiments in Weyl semimetals, we choose the vector potential to be of compact support in space, and set N 5 β,L;t . We now want to compute ∂ t N 5 t , which represents the quasi-particle flow from the first Weyl node to the second, at quadratic order in A. From the construction of the interacting Gibbs state of the system in the β, L → ∞ limit, which Theorem 2.1 is based on, one finds that the time derivative commutes with the β, L → ∞ limit. Moreover, note that, by differentiating (2.48) with respect to time, the right side simply gets multiplied by 2η, because its time dependence is all inÂ i, p 1 (t) = e ηtÂ i, p 1 (0) andÂ j, p 2 (t) = e ηtÂ j, p 2 (0). Note also that 2η is minus the zero-th component of p = −(2η, 0) (recall that in our case p = −p 1 − p 2 = 0). Therefore, by using (2.50) in the expression obtained from (2.48) by differentiating with respect to time and by taking β, L → ∞, we find: where, in the first line, summation over repeated indices is understood (α, β run from 0 to 3, while i, j from 1 to 3), and ∂ 0 denotes the time derivative. The error term e A (t), due to the error term R 5 ν,σ in (2.50), collects contributions involving a higher number of derivatives on the vector potential. It is subdominant for a vector potential slowly varying in space. For the second identity, we used the definitions of electric and magnetic fields: [41,43,45,71,72] in different Weyl semimetals reported evidence for an anomalous flow of quasi-particles between the Weyl nodes, which causes a negative, highly anisotropic, magnetoresistance. For instance, in [71], the resistivity tensor was measured for different values of the angle between E and B, and the response proved to be strongest for E and B parallel, see e.g. [71, Fig. 6A and 6B], in agreement with the E · B dependence in (2.51). Comparing measures at different angles is a way to access directly the quadratic response, which is the one related to the anomaly. We stress that the E · B dependence found in real (interacting) Weyl semimetals is the same predicted originally for noninteracting systems [63]. This is a first important instance of universality, even though a precise, quantitative, experimental verification of the interaction-independence of the chiral anomaly coefficient has still to come.
The dominant term in the right side of (2.51) is the usual chiral anomaly of QED 4 , with the same universal prefactor, irrespective of the presence of the interaction. The interaction independence of the coefficient is the analogue, in our context, of the Adler-Bardeen anomaly non-renormalization theorem [3], which our result is a rigorous version of. In our context, the lattice plays the role of a fixed ultraviolet regularization of an effective QFT model: chiral and Lorentz symmetry are emerging and approximate but nevertheless the chiral anomaly satisfies the AB non-renormalization property. We expect that our methods can be used to prove a generalization of the Adler-Bardeen's theorem for interacting lattice gauge theories at finite cutoff, such as infrared lattice QED 4 , see [50,53].

Proof of Theorem 2.1
3.1. General strategy and core argument of the proof. The proof of Theorem 2.1 is based on the following strategy: 1. We first compute the generating function of correlations via RG methods; the output of the RG construction is that the correlation functions are expressed in terms of a series expansion in λ and in a sequence of effective parameters, which is convergent, provided that |λ| is sufficiently small and the effective parameters are uniformly bounded under the iterations of the RG map. 2. Next, we show that these effective parameters are uniformly bounded, as desired, provided that the staggered chemical potential ν is chosen appropriately. We also show how to fix the bare parameters Z 5 μ,bare in such a way that (2.46) is verified. 3. The RG expansion also allow us to identify an explicit dominant contribution to the correlations (for momenta close to the Weyl nodes), which can be written explicitly in terms of the 'dressed parameters' of the model. The subdominant contributions to the correlations admit improved dimensional bounds, which imply better regularity properties in momentum space, as compared to their dominant counterparts. 4. Once that the correlations have been written as a dominant contribution, plus a betterbehaved remainder, we are in business for proving (2.50): in fact, the left side of (2.50) can be written as the contribution from the dominant part, which can be computed explicitly, plus a remainder, whose value at small momenta is fixed by Ward Identities.
In the next three subsections, Sects. 3.1.1-3.1.3, containing the core argument of the proof, we will describe more technically: (1) the outcome of items 1 to 3 for the quadratic responseˆ 5 μ,ν,σ , see in particular (3.1) below; (2) the proof of our main result, Eq. (2.50), starting from the representation formula (3.1) forˆ 5 μ,ν,σ . The main ingredients involved in these steps are the following: (i) The decomposition of the quadratic response into an explicit, non-differentiable, term, proportional to the usual chiral triangle graph of QED 4 with momentum cutoff (the function I μ,ν,σ in (3.1) below), plus a correction (the function H 5 μ,ν,σ in (3.1) below). Note: the correction is not subdominant; rather, it gives a contribution to the chiral anomaly comparable with the one from the first term. However, since it comes from the integration of the infrared-irrelevant terms, it is more regular than the first one: more precisely, it is continuously differentiable in the momenta, and its derivative is Hölder continuous (contrary to the first term, whose derivative is discontinuous at p 1 = p 2 = 0); this allows us to expand it in Taylor series up to first order included. Let us stress that the decomposition (3.1) is, obviously, nonunique: it is natural in a Wilsonian RG computation of the correlation functions, the first term being obtained by neglecting irrelevant terms in the RG sense. Due to the presence of a momentum cutoff, this term breaks gauge invariance, as apparent from (3.7) below. Such a term has the correct tensorial structure (i.e., its divergence with respect to the first index is proportional to α,β p 1,α p 2,β ε α,β,ν,σ , see (2.50)). However, if we neglected the remainder, it would produce an incorrect prefactor in the right side of (2.50) (even in the absence of interactions). (ii) The explicit computation of the first term in the right side of (3.1), i.e., of the chiral triangle graph of QED 4 with momentum cutoff, leading to (3.6)-(3.7) below. (iii) The use of Ward Identities, guaranteeing two crucial facts: (1) the vectorial vertex renormalization is proportional to the Fermi velocity, see (3.4); (2) the first order Taylor coefficients of the correction term to the quadratic response, which are well defined thanks to its differentiability, are uniquely fixed by the conservation of the vectorial current, namely by the condition (3.2). (iv) The normalization of the chiral current, i.e., condition (2.46), which guarantees that the chiral vertex renormalization equals the vectorial one, see (3.5).
Let us now discuss the core argument of the proof in more detail.
3.1.1. The splitting ofˆ 5 μ,ν,σ into chiral triangle graph + differentiable correction As an outcome of the RG analysis mentioned in item 1, of the choice of ν mentioned in item 2 and of the splitting of the correlation functions into dominant plus subdominant parts mentioned in item 3, we will prove that the quadratic response of the chiral current to the external electromagnetic field can be written aŝ where: the first term in the right side is defined in terms of the effective parameters Z 5 μ , Z μ (the chiral/non-chiral vertex renormalizations 9 ), of Z (the wave function renormalization, see (2.27)), of v l (the dressed velocities, see (2.27)), and of the explicit function I μ,ν,σ (the 'chiral triangle graph' of QED 4 with momentum cutoff, see (3.141) below), computed at the second term in the right side is a correction term, which is differentiable in the external momenta in a sufficiently small neighbourhood around the origin, with derivatives that are Hölder continuous of order θ , for all θ ∈ (0, 1).

Ward identities and the chiral triangle graph
Before we proceed further, let us specify a few additional properties of the quadratic responseˆ 5 μ,ν,σ and of the first term in the right side of (3.1). 9  The first key property we wish to emphasize is a remarkable consequence of the conservation (2.37), namely the Ward Identity We stress that there is no similar identity for the contraction of the μ index with the external momenta, due to the lack of lattice chiral gauge invariance: while the formal infrared limit of both the electromagnetic and chiral currents are conserved, only the first one is associated with a lattice gauge invariance principle. Of course, gauge invariance implies the validity of Ward Identities for (infinitely many) other correlations. A crucial one is the 'vertex Ward Identity' relating the vertex which implies the following relation between the renormalized parameters: with v 0 := 1. The condition (3.4) can be interpreted by saying that the electric charge transported by the electromagnetic current is not renormalized, thanks to lattice Ward Identities. The analogous identity does not hold for the lattice chiral current, because there is no lattice chiral symmetry protecting the 'chiral charge' from renormalizing. This is the reason why we need to impose the identity of the charges transported by the electromagnetic and chiral currents via the constants Z 5 μ,bare , as discussed in Sect. 2.5; in fact, by imposing (2.46), we obtain that, for μ = 0, 1, 2, 3, Finally, an explicit computation of the chiral triangle graph I μ,ν,σ (p 1 , p 2 ) shows that where ε αβνσ is the Levi-Civita symbol, see footnote 8, and where both R and R are smaller than C P 3 /| p + F | for P = max{|p 1 |, |p 2 |} sufficiently small, as compared with | p + F |. Note that the fact that the 4-divergences in the left side of (3.7) is non-zero, contrary to (3.2), is ultimately due to the fact that the triangle graph I μ,ν,σ (p 1 , p 2 ) is computed in the presence of an ultraviolet momentum cutoff, see (3.141) below, which breaks global and chiral gauge symmetries explicitly.

Roadmap
In view of the previous subsection, in order to complete the proof of (2.50), we 'just' need to prove the validity of (3.1), with Z μ , Z 5 μ , Z , v l satisfying (3.5)-(3.4), and with I μ,ν,σ satisfying (3.6)-(3.7); we also need to establish (3.2). The proof of all these claims will be given below, together with the proof of the other claims in Theorem 2.1. We will follow the strategy outlined in items 1 to 3 at the beginning of this section. More in detail, the proof in the next sections is organized as follows: • In Sect. 3.2, we represent the generating function of Euclidean correlations in terms of a Grassmann functional integral. In particular, in Sect. 3.2.1 we state the gaugeinvariance property of the Grassmann generating function, which implies a hierarchy of Ward Identities for the correlation functions, including (3.2). • In Sect. 3.3 we describe the iterative RG computation of the Grassmann generating function, whose output is a convergent expansion for the generating function, in terms of a sequence of effective coupling constants Z μ,h , Z 5 μ,h , Z h , v l,h , labelled by the step h of the RG iteration (these are nothing but the finite-h analogues of the parameters Z μ , Z 5 μ , Z , v l in (3.1)). If these parameters are suitably bounded, uniformly in h, and if they admit a limit as the number of RG iterations tends to infinity, the expansion for the correlation functions is convergent uniformly in the β, L → ∞ limit, for |λ| small enough, and the limit of the correlations as β, L → ∞ exists. In order for the bounds on the radius of convergence to be uniform in ζ ∈ I , the RG iteration must be defined differently, depending on whether the momentum scale under consideration in the RG step is larger or smaller than the separation between the Weyl nodes; the two regimes are described in Sects. 3.3.1 and 3.3.2, respectively. • In Sect. 3.4 we explain how to fix the bare staggered chemical potential ν in such a way that the sequence of effective coupling constants Z μ,h , Z 5 μ,h , Z h , v l,h satisfies the desired bounds. We also show how to fix the bare parameters Z 5 μ,bare in such a way that (3.5) holds. • In Sect. 3.5 we explain how to use the convergent expansion provided by the iterative RG computation of the generating function, in order to define a splitting of the main correlation functions of interest into an explicit dominant part plus a remainder, whose Fourier transforms admit improved dimensional bounds at low external momenta. In particular, we prove (3.1) and (2.27), we show that (3.5) implies (2.46) and, by using the vertex Ward Identity, we prove (3.4).
Further technical details are deferred to the appendices: in Appendix A we discuss the symmetries of the Grassmann action and their consequences for the effective action obtained at the h-th step of the RG iteration; in Appendix B we provide further technical details on the 'tree expansion' for the kernels of the effective action and their dimensional bounds; in Appendix C we prove that the quadratic response coefficients we are interested in can be expressed in terms of Euclidean correlation functions.

Grassmann representation.
In this section we introduce a Grassmann integral formulation of the model, which will be the starting point for our RG analysis. In the following, C, C , . . ., denote universal constants (in particular, independent of the distance between the Weyl nodes), whose specific values may change from line to line.
Let B L be the finite-volume Brillouin zone, Eq. (2.7), and let D β,L ,N := M f β,N × B L . We consider the finite Grassmann algebra generated by the Grassmann variables {ψ ± k,r } with k = (k 0 , k) ∈ D β,L ,N and r = 1, 2.
The Grassmann Gaussian integration P N (dψ)(·) is a linear functional acting on the Grassmann algebra as follows. Its action on a given Grassmann monomial n j=1ψ ε j k j ,q j is zero unless |{ j : ε j = +}| = |{ j : ε j = −}|, in which case: We represent the free propagator as: , these functions satisfy the bounds in (2.14) and (2.15).
The Grassmann Gaussian integration might also be rewritten as: with dψ(·) is the usual Berezin integration over the Grassmann algebra and N β,L ,N a normalization factor: Let us introduce the configuration space Grassmann field as follows, for x ∈ [0, β)× L : we then have: Denoting by g β,L (x) the two-point function of the Fock-space model, it is well-known (and easy to check) that lim N →∞ g β,L , Next, we introduce the interaction of the Grassmann field theory as: From now on, we assume that ν is real-analytic in λ, for λ small, and |ν| ≤ C|λ|. A posteriori, we will see that this assumption is compatible with the other requirements on ν we will need. More generally, the interaction of the Grassmann field theory in the presence of complex external fields A μ,x , A 5 ν,x (depending now on coordinates x = (x 0 , x), with x 0 an imaginary time in [0, β)) and of a Grassmann external field φ ± x,r is: where (ψ + , φ − ) := 2 r =1 dxψ + x,r ψ − x,r and (3.27) with the understanding that (A 0 , j 0 ) = dx A 0,x j 0,x , with j 0,x = −i 2 r =1 ψ + x,r ψ − x,r , and similarly for (A 5 μ , j 5 μ ); here j 5 μ denotes the Grassmann counterpart of the chiral lattice current, recall Eq. (2.47): where x→y A z 0 · d denotes the line integral 1 0 A (z 0 ,x+s(y−x)) · (y − x) ds. Finally, we introduce the generating functional of the correlations as: (3.28) The existence of the ultraviolet limit N → ∞ is essentially model independent, and has already been proven in a number of places, see e.g. [27,29] for the Hubbard model on the honeycomb lattice. It is well-known, see e.g. Section 5.1 of [30], that the Fockspace Euclidean correlation functions can be expressed in terms of the derivatives of the generating functional with respect to the external fields. At non-coinciding points: 5;β,L μ;r,r (z, x, y) = By translation invariance, all the correlation functions in (3.29) only depend on the relative differences of the configuration-space arguments.  0, x, y), (3.33) and their β, L → ∞ limits, denoted byŜ 2;r,r (k),ˆ μ;r,r (k, p),ˆ 5 μ;r,r (k, p),ˆ 5 μ,ν,σ (p 1 , p 2 ), respectively. For the purpose of proving Theorem 2.1, we can limit ourselves to computing these functions at sufficiently low momenta (more precisely, at k sufficiently close to the Weyl nodes and p, p 1 , p 2 sufficiently close to 0). In particular, in the proof below, we can freely assume thatÂ μ,p andÂ 5 μ.p are supported in a sufficiently small neighbourhood of the origin, and we shall do so in the following.
for any smooth function α x on R 3 , periodic in x 0 of period β and in x of period L (here φe iα is a shorthand for {e +iα . Differentiating Eq. (3.34) with respect to α and with respect to the external fields, we obtain a hierarchy of Ward identities. As already emphasized in Sect. 3.1.2, two Ward Identities that are of particular importance for us are (3.2) and (3.3): the first is obtained by deriving once with respect to α, once with respect to A, once with respect to A 5 , then setting the external fields A, A 5 , φ to zero, and taking the β, L → ∞ limit; the second is obtained by deriving once with respect to α, twice with respect to φ, then setting the external fields A, A 5 , φ to zero, and taking the β, L → ∞ limit.

Renormalization group analysis.
In this section we sketch the RG analysis of the model in the presence of external fields A and A 5 . Our analysis follows and extends the one in [51,52], where RG methods have been used to prove the analyticity of the free energy of a specific lattice model of Weyl semimetals within the class of models considered in this paper, and to compute the corresponding two-point functionŜ β,L 2 . In the discussion below, we limit ourselves to describe the general scheme and to highlight the main differences compared to [51,52], referring to those papers, and in particular to [52, Sections 2 and 3], for additional technical details. Moreover, for simplicity, we focus on generating functional W β,L (A, A 5 , 0). The generalization to the presence of an external fermionic field φ is straightforward and will not be explicitly worked out, see, e.g., [25,Section 12] or [29, Section 2.2] for a discussion of the necessary modifications.
Remark. If h * = 0, then the infrared analysis simplifies, in that the first regime, discussed in the next subsection, disappears. This is the case when ζ is far from the value ζ 1 at which the Weyl points merge (or, similarly, the case, mentioned after (2.11), in which the Weyl nodes are well separated uniformly in ζ ).

3.3.1.
Regime h > h * Assume h * < 0 (otherwise the reader can pass directly to next subsection). We inductively assume that, for all h * ≤ h ≤ 0, the generating function of correlations can be written as: (3.42) The Gaussian Grassmann field ψ (≤h) has Fourier-space covarianceĝ (≤h) (k) given by: where: with v l,h real-analytic in λ, such that |v l,h − v 0 l | ≤ C|λ|, for l = 1, 2, and |ζ h | ≤ C|λ|. The effective potential V (h) has the form: which is similar to (3.39), with the difference that now there are the operators∂ q i , with i = 1, . . . , 2n, acting on the Grassmann variables; here the labels q i are multi-indices of the form , 1, 2, 3, 4}, and∂ q i is a pseudo-differential operator, equal to the identity if q i = (0, 0, 0, 0), and dimensionally equivalent to the composition of a derivative of order q 0 i in direction 0, of a derivative of order q 1 i in direction 1, etc, otherwise 10 (if h = 0, the only non-vanishing contribution to the right side of (3.47) is the one with q i = (0, 0, 0, 0), for all i = 1, . . . , 2n, and (3.47) reduces to (3.39)). The kernels W (h) n,m 1 ,m 2 ,q;r ,μ,ν belong to a suitable weighted L 1 space, see (3.74) below, and are analytic in λ for |λ| small enough, uniformly in h (recall that ν is assumed to be analytic in λ, as well, and of order λ). We stress that the representation in (3.47) is not unique: the claim is that there exists such a representation, with the kernels satisfying natural dimensional estimates, see (3.74) below. The generating functional of correlations on scale h, W (h) , admits a similar representation, except that only kernels with n = 0 contribute. These assumptions are true at scale zero, with v l,0 = v 0 l , ζ 0 = 0 and Z 0 = 1.
The inductive assumption (3.42) is verified at scale h = 0, as an outcome of the integration of the ultraviolet degrees of freedom. We now assume that (3.42) is valid at scale h * < h ≤ 0 and prove it for h − 1. For this purpose, we intend to integrate out the degrees of freedom on scale h, after having properly rewritten the effective potential in the right side of (3.42). As a first step, we split where L is the localization operator, acting on V (h) as follows 11 : where dk (2π) 4 and dp (2π) 4  whereŴ (h),∞ 1,0,0 (k) indicates the β, L → ∞ limit ofŴ (h) 1,0,0 (k), and similarly for W (h),∞ 1,1,0;μ (k, p) andŴ (h) 1,0,1;μ (k, p), whose existence follows from the inductive construction of the kernels and from the corresponding bounds, uniform in β, L described below. The renormalization operator is defined as R := 1 − L. Notice that RL = LR = 0. 10 In Fourier space, the Fourier symbol of∂ q is a function of k behaving like k 3 for k small. When writing the analogue of (3.47) in Fourier space, we can re-absorb the Fourier symbols of the operatorŝ ∂ q i , with i = 1, . . . , 2n, into the Fourier symbol of the kernel W (h) n,m 1 ,m 2 ,q;r ,μ,ν , and we shall do so; after summation over q, we will denote the resulting modified Fourier symbol of the kernels byŴ (h) n,m 1 ,m 2 ;r ,μ,ν . 11 Whenever it will be convenient, from now on, the dependence upon the indices r 1 , r 2 , . . . , will be left implicit.ψ + k (resp.ψ − k ) will be thought of as a row (resp. columnn) vector with componentsψ + k,1 ,ψ + k,2 (resp. ψ − k,1 ,ψ − k,2 ). When writingψ + k M(k)ψ − k with M(k) a 2×2 matrix, we will mean 2 r,r =1ψ + k,r M r,r (k)ψ − k,r .

Remarks. 1) There are potentially other terms in the Taylor expansion ofŴ
(h),∞ 1,0,0 (k) that we could include 'for free' in the definition in the first line, namely k 3 1, {k μ } μ=0,1,2,3 , {k μ k 3 } μ=0,1,2,3 , and k 3 3 , where some of the terms have not been written explicitly, simply because they are zero. Consequently, RŴ (h) 1,0,0 (k) is equal to the corresponding Taylor remainder ofŴ (h),∞ 1,0,0 (k), which is quadratic in {k μ } μ=0,1,2 and, at k 0 = k 1 = k 2 = 0, is quartic in k 3 ; plus a finite size correction, proportional toŴ which is exponentially small in β, L as β, L → ∞, and can be bounded as discussed in [34,Appendix B]. We anticipate the fact that the scaling dimension ofŴ (h),∞ n,m 1 ,m 2 ;μ,ν is 7 2 − 5 2 n − m 1 − m 2 , see (3.74) (the convention here is that kernels with positive/zero/negative scaling dimensions correspond to relevant/marginal/irrelevant operators in the RG sense); any additional derivative with respect to k μ decreases the scaling dimension by 1, if μ = 0, 1, 2, and by 1/2, if μ = 3. Therefore, from the comments above, it follows that RŴ 1,0,1;μ (k, p) are equal to Taylor remainders that are linear in k μ , p μ with μ = 0, 1, 2, and quadratic in k 3 , p 3 , up to a (subdominant) finite size correction. From the formula of the scaling dimension anticipated in the previous item, it follows that they are irrelevant with scaling dimension −1. Note that the linear terms (k 3 1,0,1;μ (0, 0) in the second and third lines of (3.49) are irrelevant with dimension −1/2; the reason why we prefer to include them into the local part is to guarantee that the irrelevant part has largest scaling dimension equal to −1. This guarantees that the improved dimensional bound discussed in the paragraph after (3.74) has an additional factor γ θ h , with θ any positive constant smaller than 1 (rather than 1/2), see also Appendix B; this fact is useful in the study of the flow of the effective couplings discussed in Sect. 3.4, see in particular case μ = 3 in (3.107) and case μ = 0, 1, 2 in (3.114), where having θ > 1/2 has important consequences for the resulting flow of Z μ,h and Z 5 μ,h .
The next lemma shows that the lattice symmetries of the model impose non-trivial constraints on the kernels of LV (h) .
Let us now go back to (3.54). By using the addition principle for Gaussian Grassmann variables, we can rewrite it as: where P h (dψ (h) ) has covariance given by g (h) , P ≤h−1 (dψ (≤h−1) ) has covariance given by g (≤h− 1) , and, in the exponent in the second line, ψ (≤h) = ψ (≤h−1) + ψ (h) . We now integrate the field ψ (h) and denote the logarithm of the result of the integration in the sec- This iterative integration procedure goes on until we reach scale h * , at which point the procedure is changed, as described in the next subsection. As discussed, e.g., in [25,27,29,52], the effective potential and generating function at scale h, obtained via such an iterative procedure, can be represented as convergent series over Gallavotti-Nicolò (GN) trees, see in particular [52, Section 2], where the tree expansion for a model of Weyl semimetal in the same regime as the one considered here is discussed in detail.
The kernels of the single-scale contribution to the generating function, W (h) (A, A 5 ) satisfy an estimate analogous to (3.74) with n = 0. The combination 7 2 − 5 2 n − m 1 − m 2 − d(q) is the scaling dimension of the kernels with 2n Grassmann fields, derivatives of order q acting on them, and m 1 + m 2 external fields of type A or A 5 in the first regime h ≥ h * . Note, in particular, that the effective quartic interaction, i.e., the kernel with n = 2, no derivatives, and m 1 = m 2 = 0, is irrelevant, with scaling dimension −3/2. The irrelevance of the quartic interaction allows us to derived improved bounds on all the contributions to the kernels associated with GN trees containing at least one interaction endpoint: this is the same as in models of graphene with short range interactions, see [27] and, in particular, [27,Theorem 2]. More precisely, we can split W n,m 1 ,m 2 ,q , where 'd' stands for dominant and 'r' for remainder, and where the second term collects the contribution of all GN trees with at least one endpoint of type ν k or on scale 1 (in particular, it includes all the contributions from GN trees with at least one endpoint associated with a quartic interaction, which is necessarily on scale 1); the term W (h);r n,m 1 ,m 2 ,q admits an improved dimensional bound, analogous to (3.74), but with an extra factor γ θ h in the right side, for any fixed θ smaller than 1 (and the constant C replaced by a θ -dependent constant C θ ). The proof of (3.74) and of its improved analogue for W (h);r n,m 1 ,m 2 ,q are completely analogous to those of Eqs. (60) and (61) in Lemma 1 of [52], respectively, modulo a few minor differences discussed in Appendix B.
By using these bounds, we can also prove the inductive assumptions (3.62). In fact, assume that (3.62) and (3.63) are valid for k ≥ h. Then the bound (3.74) is valid at scale h − 1, and similarly for its improved analogue for W (h−1);r n,m 1 ,m 2 ,q . By using the definitions of z μ,h−1 with μ = 0, 1, 2, 3, and the bound on W (h−1);r 1,0,0,q , we find that (3.62) is valid with k = h − 1, as well.
Concerning the proof of (3.63), let β ν h := ν h−1 − 2ν h be the beta function for ν h . From the bound on W (h−1);int 1,0,0 , we find that for any θ ∈ (0, 1). This does not imply that any solution of the beta function equation (3.63). However, it is enough that we find one special choice of ν 0 (or, equivalently, of the staggered chemical potential ν in (2.19)) for which the solution satisfies such a bound. We temporarily proceeding by assuming the existence of such a good initial datum, and we will prove this assumption in Sect. 3.4. In that section, we will also derive bounds on the effective couplings Z μ,h , Z 5 μ,h .

3.3.2.
Regime h ≤ h * In order to integrate the remaining scales, we proceed as follows.
We choose α 0 in the definition of χ (see beginning of Sect. 3.2) small enough that the support of χ h * (k) consists of two 'well-separated' sets, 12 centered at p + F and p − F , respectively. We decompose the fermionic field in terms of two independent quasiparticle fields, associated with the two Weyl nodes (as in the previous section, we will often leave the r indices implicit, and think ψ + x,ω resp. ψ − x,ω as two-component row resp. column vectors): where in the second line the sum over k runs over the four-dimensional vectors such As in the previous regime, we integrate the scales h ≤ h * in a multiscale fashion: at each step, we decompose the quasi-particle Grassmann field ψ (≤h) ω as the sum of two independent fields, one describing the fluctuations at scale h, and the other at smaller scales: ψ ω ; we decompose the effective potential as the sum of a localized part plus an 'irrelevant' remainder; we combine part of the localized part of the effective potential with the Grassmann measure; we integrate out the field at scale h; we iterate. More precisely, we inductively assume that, for any h β ≤ h ≤ h * , with h β = log γ (π/β) , the generating functional of the correlations can be written as: (3.78) Let us explain the meaning of the various objects involved: P ≤h (dψ (≤h) ) denotes a Gaussian Grassmannn integration with covariance P ≤h (dψ (≤h) )ψ −(≤h) and: , Z h is a real-analytic function of λ, to be inductively defined below, (3.83) and, recalling that c R,ω was defined by (2.12), we letc R,ω (k ) : Moreover, the function V (h) in the right side of (3.78) is the effective potential, which can be represented in a way analogous to (3.47), namely n,m 1 ,m 2 ,q;ω,r ,μ,ν (X, Y, Z), (3.84) and W (h) (A, A 5 ) is the finite-scale contribution to the generating function, which admits a representation analogous to (3.84), with the difference that only terms with n = 0 contribute to the sums. The inductive assumption (3.78) and following equations is verified at scale h = h * , as an outcome of the integration of the first regime, provided we let (note that the representation (3.84) at scale h * follows from (3.47) at the same scale, in combination with the definition of quasi-particle fields (3.77)). We now assume that (3.78) is valid at scale h and prove it for h − 1. We decompose with L the localization operator, defined as follows: 4 Â μ,pψ   (k , p).
Remark. We anticipate the fact that the scaling dimension ofŴ  (k , p) are all irrelevant, with scaling dimension −1. The reader may recognize that the scaling dimensions are the same as those of QED 4 : our theory can be seen as a lattice realization of lattice infrared QED 4 , with the important difference that the speed of light is replaced by the anisotropic running velocities v μ,h . Since v 3,h (or, equivalently, v 0 3 ) vanishes in the limit as the Weyl nodes merge, the dimensional bounds on the effective potentials, see (3.102) below, are potentially affected by dangerous 1/v 3,h factor, which may a priori have an impact on the convergence properties of the infrared expansion of the observables of interest. However, by carefully tracking the loss and gain of factors v 3,h , due to the non-uniform dependence of the propagator upon v 3,h and to the difference of scaling dimensions between the first and second regimes, one finds that the bad and good dependences upon these factors compensates, and lead to the overall factor |v 0 3 | n−1+ q 3 1 in (3.102) below; see [52,Section 3] and Appendix B.
As in the regime h > h * , we shall use the localized term LV (h) to redefine the Grassmann integration and the coupling of the fermions with the external fields. The next lemma establishes important symmetry properties of the kernels of LV (h) (see Appendix A.2 for the proof). with Z μ,h , Z 5 μ,h ∈ R.

Lemma 3.2. The following identities hold:
We now: rewrite V (h) = LV (h) + RV (h) in the right side of (3.78); then write LV (h) explicitly, in view of Lemma 3.2, in terms of n h , z μ,h , Z μ,h and Z 5 μ,h ; then re-absorb the quadratic part of LV (h) associated with μ k μ ∂ μŴ (h),∞ 1,0,0;ω (0) ≡ μ k μ σ μ,ω z μ,h in the Gaussian Grassmann integration, thus getting (3.92) and the new Grassmann Gaussian integration P ≤h (dψ (≤h) ) has covariance given by: Similarly to the previous regime, we inductively assume that, for all scales h ≤ k ≤ h * and any θ ∈ (0, 1): where v 0 0 := 1. We will check the validity of these bounds on scale k = h − 1. We also assume that ν k satisfies (3.63), for all h ≤ k ≤ h * . Notice that (3.96), in combination with the bounds on Z h , v l,h derived in Sect. 3.3.1, implies: We also let and, of course, these effective parameters satisfy the same bounds as (3.97). To check the inductive assumption, we decompose the Grassmann field as: and we used the fact that, on the support of On the support of f h,ω , |k μ | ≤ Cγ h for μ = 0, 1, 2, and |k 3 | ≤ Cγ h /|v 0 3 |. Therefore, for these values of k , recalling the definition of h * , which implies ĝ (h) (k) ∞ ≤ Cγ −h and ĝ(k) 1 ≤ Cγ 3h /|v 0 3 |. By a similar discussion, we can also dimensionally bound the derivatives ofã where α = (α 0 , α 1 , α 2 , α 3 ) and |α| = μ |α μ |. These in turn imply the following bound in configuration space: At this point, we integrate the field ψ (h) out, as done in Eqs. (3.72) and (3.73), thus obtaining the analogue of (3.73), where the effective action V (h−1) can be represented as in (3.84), with h replaced by h − 1. Also in this case, the kernels of the effective potential can be represented in terms of a convergent GN tree expansion, and are analytic in λ, for |λ| small, uniformly in ζ , i.e., in the distance between the Weyl nodes. They satisfy the following weighted L 1 estimates: for all h β ≤ h < h * , n ≥ 1, and |λ| small enough, uniformly in ζ , if (3.96) and (3.63) are valid for h < k ≤ h * , then, see [52, Lemma 2, Eq. (93)], where δ * (Q) is the Steiner diameter measured by using the norm d(x) defined after (3.101), k were defined after (3.74)). The kernels of the single-scale contribution to the generating function, W (h) (A, A 5 ) satisfy an estimate analogous to (3.102) with n = 0. The combination 4 − 3n − m 1 − m 2 − q 1 is the scaling dimension of the kernels with 2n Grassmann fields, derivatives of order q acting on them, and m 1 + m 2 external fields of type A or A 5 in the second regime h < h * . Also in this case, the quartic interaction (n = 2, m 1 = m 2 = 0, q 1 = 0) is irrelevant, and its scaling dimension is −2. Also in the second regime, we can obtain improved bounds on all the contributions to the kernels associated with GN trees containing at least one interaction endpoint: more precisely, if we split W n,m 1 ,m 2 ,q , where the second term collects the contribution of all GN trees with at least one endpoint of type ν k or one endpoint on scale h * + 1 (including, in particular, all the contributions from GN trees with endpoints associated with quartic interactions), the term W (h);r n,m 1 ,m 2 ,q admits an improved dimensional bound, analogous to (3.74), but with an extra factor in the right side equal to: γ θ h , if m 1 = m 2 = 0 or n ≥ 2; γ θ(h−h * ) , if m 1 + m 2 ≥ 1 and n = 1 (again, θ can be chosen to be any positive exponent smaller than 1, at the cost of replacing the constant C in the right side of (3.102) by C θ , which possibly diverges as θ → 1 − ). The proof is completely analogous to that of Eq. (94) of Lemma 2 in [52], modulo a few minor differences discussed in Appendix B.
The definition of z μ,h−1 and the bound on W (h−1);r n,m 1 ,m 2 ,q;ω,r ,μ,ν , readily imply the validity of (3.96) at scale k = h − 1, as desired. Concerning (3.63), also in this second regime we find that ν h satisfies the beta function equation (3.76), with β ν h satisfying (3.75) even for h < h * . The choice of the initial datum ν 0 (or, equivalently, of the staggered chemical potential ν) for which ν k satisfies (3.63) for all h ≤ 0 will be discussed in the next section, where we will also derive bounds on the flow of the effective couplings Z μ,h , Z 5 μ,h . The iterative RG integration of the second regime goes on until we reach scale h β . At that point, we integrate 'in one shot' all the remaining Grassmann degrees of freedom, thus obtaining the desired generating function of correlations, in the form of a sum of single-scale contributions from h = h β to h = 0 (the covariance of ψ (≤h β ) ω admits the same dimensional bounds as the one of ψ (h β +1) ω ; therefore, the result of the integration of ψ (≤h β ) ω admits the same qualitative bounds as the one at scale h β + 1, so that the last iteration step does not give any additional difficulty).
The uniformity in β, L and ζ of the convergence of the GN tree expansion for the single-scale contributions to the generating function, as well as the dimensional bounds (3.74) and (3.102), imply that the correlation functions at fixed space-time positions are analytic functions of λ and {ν h } h≤0 , in a sufficiently small neighbourhood of the origin, uniformly in β, L and ζ . Elaborating on this and on the fact that all the Taylor coefficients of such convergent expansions admit a limit as β, L → ∞ implies the existence of the β, L → ∞ limit of the Euclidean correlations, as well as the analyticity of the limits, stated in Theorem 2.1, provided that the sequence {ν h } h≤0 satisfies the promised bounds, see (3.63), and is itself analytic in λ. This will be proved in the next section, together with the bounds on the flow of the effective vertex functions Z h , Z 5 h .

The flow of effective couplings.
In this section we discuss the RG flow of the effective parameters ν h , Z μ,h , Z 5 μ,h (as well as the one ofZ 3,h andZ 5 μ,h with μ = 0, 1, 2, which appear in the first regime). In particular, we explain how to fix the bare staggered chemical potential ν in (2.19) in such a way that the running staggered chemical potential ν h satisfies (3.63) at all scales h ≤ 0. We also explain how to fix Z 5 μ,bare , see (2.39), in such a way that their dressed counterpart, Z 5 μ,h , flow to any prescribed 4-tuple of values as h → −∞. The bounds discussed in the following use a few properties of the GN trees contributing to the effective potential, discussed in [52, Sections 2 and 3], see also Appendix B.

Flow of ν Starting from the beta function equation for ν h , (3.76), we find
for any k < 0. Here ν 0 is an analytic function of λ and ν that, assuming ν to be of order λ, is of the form ν 0 = ν + O(λ 2 ). If we now require that ν k tends to zero as k → −∞, we get and, more in general, Recall that β ν h is a function of λ and of the effective parameters ν h on scales h ≥ h. Therefore, we can regard the right side of (3.105) as a function of the whole sequence ν := {ν k } k≤0 , which we denote by T k (ν, λ), so that (3.105) can be read as a fixed point equation ν k = T k (ν, λ) for the sequence ν. By proceeding as done in many other papers before, see e.g. [34, Section 6.4.2], we look for a fixed point in the Banach space of sequences ν such that ν θ := sup h≤0 |ν h |γ −θ h ≤ K |λ|, for θ = 3/4 (say) and K a suitable (sufficiently large) constant. Following the same strategy as [34,Section 6.4.2] (in a much simpler setting), the reader can check that the map T : ν → {T k (ν, λ)} k≤0 is a contraction on such a Banach space, which implies the existence of a unique fixed point in that space. The value of ν 0 of such a fixed point sequence corresponds to the 'right' initial datum to assign in order for (3.63) to be satisfied at all scales. Finally, note that fixing ν 0 is equivalent to fixing ν (recall that ν 0 = ν 0 (ν, λ) is analytic in ν, λ and that ν 0 (ν, λ) = ν 0 + O(λ 2 ); appealing to the implicit function theorem, we can analytically invert ν 0 with respect to ν). The value of the bare staggered chemical potential ν fixed via this strategy is the one stated in Theorem 2.1.

Flow of Z μ,h
In the first regime, recalling that Z μ,h is defined via (3.51), we write where β μ,h includes the contributions from GN trees that have at least one endpoint of type ν h or one endpoint on scale 1 of order λ. Therefore, β μ,h is bounded in the same way as W (h);r 1,1,0,q , with q 1 = q 3 1 equal to 0 or 1, depending on whether μ = 0, 1, 2 or μ = 3, respectively; see the paragraph after (3.74). We thus get, for any 0 < θ < 1, where Z μ,0 for μ = 0, 1, 2, 3 andZ 3,0 are analytic functions of λ, bounded as |Z μ,0 − v 0 μ | ≤ C|λ| |v 0 μ | and |Z 3,0 | ≤ C|λ| |v 0 3 |. Remark. The importance of having a gain factor γ θ h with θ any positive constant smaller than 1 (rather than 1/2, as other simpler choices of the localization operator would have implied), is apparent from (3.107). In fact, choosing θ larger than 1/2 makes all the components of the beta function summable in h, uniformly in v 0 3 . A posteriori, this motivates the definitions (3.49), see also the remarks following it. Similar considerations are valid for the beta function for Z 5 μ,h , see in particular (3.114) below.
Similarly, the flow equation forZ 3,h is with From these bounds on the beta function, we readily find for any h ≥ h * , uniformly in h.
In the second regime, the flow of Z μ,h is controlled by a flow equation of the same form as (3.106), where β μ,h is bounded in the same way as W (h);r 1,1,0,0 , see the paragraph after (3.102): (3.111) Using this bound on the beta function and (3.110), we readily find that for any h ≤ 0, that Z μ := lim h→−∞ Z μ,h exists and is analytically close to Z μ,0 , and that the limit is reached at an exponential rate:
In the second regime, the flow of Z 5 μ,h is controlled by a flow equation of the same form as (3.113), where β 5 μ,h is bounded in the same way as W (h);r 1,0,1,0 , see the paragraph after (3.102): Using this bound on the beta function and (3.117), we readily find that for any h ≤ 0, that Z 5 μ := lim h→−∞ Z 5 μ,h exists and is analytically close to Z 5 μ,bare , and that the limit is reached at an exponential rate: . Note that these facts imply that the relation between Z 5 μ and Z 5 μ,bare can be inverted into Z 5 μ,bare = (1 + O(λ))Z 5 μ , thus allowing us to fix the dressed chiral renormalizations as desired. In particular, in the following we will need Z 5 μ ≡ Z μ ; therefore, we will fix the bare chiral parameters in such a way that this condition is satisfied.

Asymptotic infrared behaviour of the correlation functions.
In this section we use the RG construction described in the previous sections to study the β, L → ∞ limit of the correlation functions in Equations (3.30) to (3.33). Our goal is to isolate the dominant part, at large distances and/or momenta close to the Weyl nodes, which has an explicit structure in terms of the dressed parameters Thanks to the bound (3.102) (or, better to its analogue with n = 0) and to the boundedness of Z := sup k≤0 Z k and of Z 5 := sup k≤0 Z 5 k , proved in the previous section, we see that the h-th term in the sum in the right side of (3.122) is bounded by CZ 5 (Z) 2 γ h /|v 0 3 |. Therefore, the sum in the right side of (3.122) is absolutely convergent. Moreover, each term in the sum is continuous with respect to p 1 , p 2 , uniformly in h. Hence, (3.122) is continuous as a function of the external momenta.
Note, however, that the bound (3.102) does not imply differentiability of (3.122) in the external momenta. In fact, it implies that the derivatives in p 1 or p 2 ofŴ The summands and their derivatives are bounded via (3.74), (3.102), and its improved version for the remainder term, discussed after (3.102). These bounds imply that (3.123) is not only continuous in the external momenta, but also differentiable: in fact, they imply that the derivatives in p 1 or p 2 ofŴ , and that those ofŴ . Therefore, the derivatives in p 1 or p 2 of (3.123) are finite, and bounded by CZ 5 (Z) 2 |v 0 3 | −1 . Moreover, they are continuous in a sufficiently small neighbourhood of the origin; more precisely, the bound proportional to γ θ(h−h * ) on the derivatives ofŴ (h),∞;r 0,2,1;μ,ν,σ (p 1 , p 2 ) with h < h * implies that the derivatives in p 1 , p 2 of (3.123) are Hölder continuous of exponent θ > 0, for any θ < 1.
The dressed propagator and the vertex functions. We now consider the dressed prop-agatorŜ 2 (k) and the vertex functionsˆ μ (k, p) andˆ 5 μ (k, p), which are obtained by an RG construction analogous to the one described in Sect. 3.3, in the presence of the Grassmann source field φ, see (3.26) and (3.28). We decomposeŜ 2 (k) (resp.ˆ μ (k, p), resp.ˆ 5 μ (k, p)) in a way analogous to (3.121): that is, we distinguish a dominant part, analogous to (3.122), which consists of the sum of the single-scale contributions from the GN trees with one endpoint of type φ + ψ − and one of type ψ + φ − (resp. one endpoint of type Aψ + ψ − , one of type φ + ψ − , one of type ψ + φ − , resp. one endpoint of type A 5 ψ + ψ − , one of type φ + ψ − , one of type ψ + φ − ), all on scales smaller than h * , from a remainder term, analogous to (3.123), which includes all the other contributions (those from GN trees with root on scale smaller than h * and at least one endpoint on scale h * + 1, and those from GN with root on scale h * or larger). The dominant term is further decomposed in analogy with (3.124) into a 'relativistic' contribution, obtained via the substitutions spelled in items (i) and (ii) after (3.124), plus an additional remainder. The relativistic contribution is explicit: for the dressed propagator with k close to p ω F , letting where χ h,ω = k≤h f h,ω , andĝ rel ω (k) is defined in the same way as (3.125), with f h,ω (k) replaced by 1, that is,ĝ rel where v 0 = 1 and σ μ,ω was defined in (3.89) (note, a posteriori, that the cutoff function χ h * −1,ω can be dropped from the right side of (3.129) for k small enough, because χ h * −1,ω ≡ 1 in a sufficiently small neighbourhood of the origin). Similarly, for the vertex functions with k close to p ω F , letting again k = k − p F , the relativistic contribution readsˆ On the other hand, the remainder terms (the analogue of (3.123) and the analogue of the non-relativistic term in (3.124)) admit an improved dimensional bound: if k − p ω F , p and k + p − p ω F are all small and of the order γ h , the remainder terms have an additional factor γ θ h , as compared to the corresponding dominant, relativistic, term.
Summarizing, letting k = k − p ω F and assuming that |k |, |k + p|, |p| are all in [cγ h , Cγ h ], for some h < h * and c, C > 0, we find thatŜ 2 (3.133) and (2.27) into the vertex Ward Identity (3.3), we obtain Eq. (3.4), that is, Z μ = v μ Z . As already mentioned, this identity says that the charge carried by the electromagnetic current is not renormalized by the interaction.
There is no analogue of the vertex Ward Identity for the chiral vertexˆ 5 μ (k, p). Therefore, the condition (2.46) must be enforced by properly fixing Z 5 μ,bare . By using (3.133) and (3.134), we find that the left side of (2.46) equals (Z 5 μ /Z μ )ω1 2 , so that (2.46) reduces to (3.5), Z 5 μ = Z μ . This condition is enforced by fixing the bare parameters Z 5 μ,bare in the way discussed at the end of Sect. 3.4.3.
Now, A νσ (p 1 , p 2 ) = B νσ (p 1 , p 2 ) = D νσ (p 1 , p 2 ) = 0 by simple parity reasons: in fact, after the computation of the trace, letting χ be the derivative of χ with respect to |k|, k μ = k μ /|k|, and ανβσ the Levi-Civita symbol (see footnote 8), which are all zero by the anti-symmetry in α ← →β. The only non-trivial term we are left with is where, using that the angular integration in the 4D integral over k gives 2π 2 , we can rewrite where in the second equality we used that ∞ 0 χ 2 (ρ)χ (ρ)dρ = −1/3, where, with some abuse of notation, we denoted χ(|k|) ≡ χ(k). Plugging this back into (3.154) and using again the anti-symmetry in α ← →β, we find (recalling that q = p 1 + p 2 ) C νσ (p 1 , p 2 ) = − 1 12π 2 p 1,α p 2,β ε ανβσ , (3.155) so that, putting things together, which proves (3.6) (note that the order of the indices in ε αβνσ in the right side of (3.6) is different from the one in the right side of (3.156), which explains the different sign).

B. Dimensional Bounds on the Kernels of the Effective Potential
In this appendix we provide some additional details on the proofs of (3.74), (3.102) and of their improved analogues, discussed after (3.74) and (3.102), respectively. As anticipated above, the proofs of these bounds are completely analogous to those of Lemma 1 and Lemma 2 in [52], where a specific model of Weyl semimetal within the class of models considered in this paper was analyzed. Since there a few minor differences in the bounds (3.74), (3.102) and in their improved analogues, as compared to those stated in [52, Lemma 1 and 2], here we highlight where these differences rely and explain their origin. A first macroscopic difference is that [52] explicitly treats only the free energy (i.e., A = A 5 = 0), which obviously has an impact on the definition of the iterative step. In addition to this, there are other more technical differences, which we discuss separately for the first and second regimes. First regime. The RG construction in the first regime is virtually the same as the one described in [52, Section 2], modulo a slightly different definition of localization (compare [52,Eq. (43)] with the definition in (3.48)-(3.49)) that, in particular, takes into account the presence of the external fields A and A 5 . At each scale h * ≤ h ≤ 0, the effective potential can be expressed in terms of a GN tree expansion virtually equivalent to the one described in [52, Section 2.1]. As already mentioned, the proof of (3.74) and of its improved analogue (see discussion after (3.74)) goes along the same lines as the proof of Eqs. (60) and (61) of [52], respectively. However, there are a few technical differences that deserve comments: • The L 1 norm in the left side of (3.74) has a stretched exponential weight, contrary to the one used in [52,Eq. (60)]. However, the inclusion of the stretched exponential weight involving the tree distance can be accomodate without extra difficulties, thanks to the stretched exponential decay of the propagator (3.101). See [8,32] for two recent works on fermionic RG where similar exponentially weighted L 1 norms are studied via the same kind of methods of this paper. operator is −1, rather than −1/2; therefore, by the same line of reasoning, the GN trees with at least one endpoint on scale 1 admit a bound with an additional 'short memory factor' γ θ h with θ a positive constant strictly smaller than 1 (rather than 'just' 1/2).
Second regime. The RG construction in the first regime is virtually the same as the one described in [52,Section 3], modulo a slightly different definition of localization (compare [52,Eq. (88)] with the definition in (3.86)-(3.87)) that, in particular, takes into account the presence of the external fields A and A 5 . At each scale h < h * , the effective potential can be expressed in terms of a GN tree expansion virtually equivalent to the one described in [52, Section 3.1]. As already mentioned, the proof of (3.102) and of its improved analogue (see discussion in the paragraph after (3.102)) goes along the same lines as the proof of Eqs. (93) and (94) of [52], respectively. Also in this case, there are a few technical differences that deserve comments: • As in the first regime, in (3.102) we use a weighted L 1 norm, rather than the standard one, and we take into account the possibile presence of the derivative labels q; as commented above, the proof of the bound remains essentially unaltered by these two modifications. • Concerning (3.102) and its improved version, an important difference with respect to [52,Eqs. (93)-(94)] is the presence of the pre-factor |v 0 3 | n−1+ q 3 1 . Even though such factor does not appear in [52,Eqs. (93)-(94)], its presence is actually proved in [52] (in the case q 3 1 = 0, but the proof in the general case is essentially unchanged), see l.1 after [52,Eq. (107)].
• Concerning the improved dimensional bound on W (h);r n,m 1 ,m 2 ,q , in the paragraph after (3.102) we claim that it has an additional factor γ θ h if m 1 = m 2 = 0 or n ≥ 2, or γ θ(h−h * ) if m 1 + m 2 ≥ 1 and n = 1, with θ any positive constant strictly smaller than 1; this has to be compared with the factor γ 1 8 (h−h * ) stated in [52,Eq. (94)]. First of all, for the same reasons explained above for the first regime (see in particular the third item of the dotted list), the factor γ [52,Eq. (94)] can be straightforwardly improved to γ θ(h−h * ) , for any 0 < θ < 1, because the largest scaling dimension of an irrelevant operator is −1. We still need to discuss the origin of an additional gain factor γ θ h * for the terms with m 1 = m 2 = 0 or n ≥ 2. The point is that these terms are defined in terms of GN trees with at least one endpoint on scale h * + 1 with m 1 = m 2 = 0 or n ≥ 2; in turn, each such endpoint has a kernel generated by the tree expansion of the first regime, and comes from GN trees with at least one endpoint on scale 1: therefore, it is bounded via the improved dimensional bound discussed in the paragraph after (3.74), which is characterized by an additional factor γ θ h * as compared to the 'basic' bound (3.74), as desired.

C. The Quadratic Response
In this appendix we compute the quadratic response coefficient in the right side of (2.48) and prove that its β, L → ∞ limit can be expressed in terms of Euclidean correlation functions. We compute I and II separately, at second order in the external field. To this end, it is convenient to study the evolution of the state in the interaction picture. Let ρ int (t) := e iH L t ρ(t)e −iH L t . Then: We denote by I (2) the second order contribution in A. If we write O (k) t (A) for the k-th order of O t (A) in A, and similarly for P (k) t (A), we get: k,k =1,2,3 dp 1 (2π) 3 dp 2 (2π) 3Â k, p 1Âk , p 2ˆ k,k ( p 1 , p 2 ).

C.2. Wick rotation for correlations of three observables.
In this section we prove (C.14). Let be the function of interest, thought of as a function of η. We add and subtract T β,L (η β ), where η β = 2π β βη 2π ∈ 2π β N. We will prove below that T β,L (η) − T β,L (η β ) is bounded by (const.)β −1 , uniformly in L. As for T β,L (η β ), it is equal to a double integral over imaginary times of the appropriate Euclidean correlation function, as implied by the following proposition. with O X even in the fermionic operators, commuting with the total number operator. We assume that the sum in (C.16) runs over subsets X such that |X | is bounded uniformly in L. Let A(z) = e iH L z Ae −iH L z , for z ∈ C. Let η i ∈ (2π/β)N (with the convention that N is the set of positive integers) and consider: Before giving the proof of Proposition A.1, let us explain how to adapt it to the case at hand. We let: C =Ĵ 5 μ, p − Ĵ 5 μ, p β,L , A(s 1 ) =Ĵ k, p 1 ,s 1 − Ĵ k, p 1 β,L , and B(s 2 ) = J k , p 2 ,s 2 − Ĵ k, p 2 β,L , so that, for η 1 = η 2 = η β , T β,L (η β ) = I β,L ABC . Using the proposition, we get withp i = (η β , p i ) and p = (−2η β , − p 1 − p 2 ). If we now take β, L → ∞ with η β → η, this expression tends to TĴ k,p 1 ;Ĵ k ,p 2 ;Ĵ 5 μ,p ∞ = TĴ 5 μ,p ;Ĵ k,p 1 ;Ĵ k ,p 2 ∞ , as desired (existence of the limit follows from the construction of the Euclidean correlation functions of Sect. 3.3).
Proof. By using the definition of T, we rewrite (C.18) as: Consider I. We apply Cauchy theorem to rewrite the integral over s 2 as follows: where we used the fact that B η 2 (z) → 0 as Rez → −∞, thanks to the factor e η 2 Rez and the fact that η 2 > 0; see Fig. 1.
We now rewrite the integral over s 1 in a similar way, thus getting − A η 1 (t 1 )B η 2 (t 1 + t 2 )C β,L + A η 1 (t 1 − iβ)B η 2 (t 1 + t 2 − iβ)C β,L . (C.24) Recalling that e iη 1 β = 1, we have A η 1 (t 1 − iβ)B β,L = B A η 1 (t 1 ) β,L , and similarly for B η 2 (t 1 + t 2 − iβ), so that (C.23) can be rewritten as: By proceeding in the same way, we can rewrite II analogously: In the first term in the second line we split the integral over t 2 as 0 −∞ dt 2 = t 1 −∞ dt 2 + 0 t 1 dt 2 . If we now combine together the two terms with the integral over t 2 from −∞ to t 1 and use the Jacobi identity again, we get: which is the same as the right side of (C.17).
In view of the discussion at the beginning of this subsection, as well as of (C.19) and following lines, in order to conclude the proof of (C.14) we are left with proving that T β,L (η) − T β,L (η β ) is bounded from above by (const.)β −1 , uniformly in L. Using the fact that 0 ≤ η β − η ≤ (2π)/β, we find that with A, B, C as defined before (C. 19). By the Lieb-Robinson bounds for multi-commutators [17], see in particular [17, item (ii) of Corollary 4.12], one finds that there exist C ABC > 0 independent of L such that the norm [C, A(s 1 )], B(s 2 ) is bounded from above by C ABC L 3 (1 + |s 1 | + |s 2 |) 6 , where 6 should be understood as twice the spatial dimension. Therefore, which vanishes in the limit as β → ∞, uniformly in L, for any η > 0. This concludes the proof of (C.14).