Strong-coupling phases of 3D Dirac and Weyl semimetals. A renormalization group approach

We investigate the strong-coupling phases that may arise in 3D Dirac and Weyl semimetals under the effect of the long-range Coulomb interaction, considering the many-body theory of these electron systems as a variant of the conventional fully relativistic Quantum Electrodynamics (QED). For this purpose, we apply two different nonperturba-tive approaches, consisting in the sum of ladder diagrams and taking the limit of a large number N of fermion flavors. We benefit from the renormalizability that the theory shows in both cases to compute the anomalous scaling dimensions of different operators exclusively in terms of the renormalized coupling constant, allowing us to determine the precise location of the singularities signaling the onset of the strong-coupling phases. We show then that the QED of 3D Dirac semimetals has two competing effects at strong coupling. One of them is the tendency to chiral symmetry breaking and dynamical mass generation, which are analogous to the same phenomena arising in the conventional QED at strong coupling. This trend is however outweighed by the strong suppression of electron quasiparticles that takes place at large N , leading to a different type of critical point at sufficiently large interaction strength, shared also by the 3D Weyl semimetals. Overall, the phase diagram of the 3D Dirac semimetals turns out to be richer than that of their 2D counterparts, displaying a transition to a phase with non-Fermi liquid behavior which may be observed in materials hosting a sufficiently large number of Dirac or Weyl points.


Introduction
In recent years, condensed matter physics has witnessed several remarkable discoveries that, starting from graphene [1], have shown the existence of electron systems where the quasiparticles have linear dispersion about degenerate points at the intersection between valence and conduction bands. This has led to introduce the class of so-called Dirac semimetals, in which the low-energy excitations can be described in terms of a number of Dirac spinors. Thus, there is already clear evidence that materials like Na 3 Bi or Cd 3 As 2 have Dirac points in their electronic dispersion [2][3][4], providing in that respect a 3D electronic analogue of graphene. More recently, there have been claims that materials like TaAs or the pyrochlore iridates should be examples of 3D Weyl semimetals [5,6], characterized by having Weyl points hosting each of them fermions with a given chirality.
Since the appearance of graphene, the possibility that the quasiparticles of a condensed matter system may behave as Dirac fermions has been very appealing, as the particular geometrical features of these mathematical objects have shown to be indeed at the origin of several unconventional effects. The Klein paradox [7] and the difficulty that such quasiparticles have to undergo backscattering [8] are a direct consequence of the additional pseudospin degree of freedom inherent to Dirac fermions. The topological insulators and their peculiar surface states can be taken also as an illustration of the fruitful interplay arising between the symmetry properties of Dirac fermions and the topology of the space [9,10].
Another important feature of known Dirac semimetals is that they need to be modeled as interacting electron systems which are naturally placed in the strong-coupling regime. This is so as the long-range Coulomb repulsion becomes the dominant interaction, while the Fermi velocity of the electrons takes values which are at least two orders of magnitude below the speed of light. Thus, the equivalent of the fine structure constant in the electron system can be more than one hundred times larger than its counterpart in the conventional fully relativistic Quantum Electrodynamics (QED). This leads to expect relevant manybody effects as, for instance, the imperfect screening of impurities with sufficiently large charge on top of graphene [11][12][13][14]. Another important property is the scale dependence of physical observables, like that predicted theoretically for the Fermi velocity in 2D Dirac semimetals [15,16] and measured recently in suspended graphene samples at very low doping levels [17]. We may think therefore of the Dirac semimetals as an ideal playground to observe scaling phenomena that would be otherwise confined to the investigation of field theories in particle physics.
In any event, a proper account of the scaling effects must rely on the renormalizability of the many-body theory. The formulation of the interacting theory of Dirac semimetals requires the introduction of a high-energy cutoff in the electronic spectrum, which can be set to a rather arbitrary level. In order to ensure the predictability of the theory, it becomes therefore crucial that all the dependences on the high-energy cutoff can be absorbed into the redefinition of a finite number of parameters. In the case of Dirac semimetals, such a renormalizability cannot be taken for granted, as the many-body theory does not enjoy the full covariance which enforces that property in typical relativistic field theories. Nevertheless, the investigation of the 2D Dirac semimetals has already provided evidence that the interacting field theory in these electron systems is renormalizable [18]. This condition has allowed for instance to compute very high order contributions to the renormalization of relevant order parameters, leading to an estimate of the critical coupling for chiral symmetry breaking that agrees with the value obtained with a quite different nonperturbative method like the resolution of the gap equation [19].
The aim of this paper is to investigate the strong-coupling phases that may arise in 3D Dirac and Weyl semimetals under the effect of the long-range Coulomb interaction. In this respect, the many-body theory of these electron systems can be viewed as a variant of conventional QED. As already pointed out, a main difference with this theory lies in that the Fermi velocity of the electrons does not coincide with the speed of light, which makes that quasiparticle parameter susceptible of being renormalized and therefore dependent on the energy scale. In three spatial dimensions, the electron charge needs also to be redefined to account for cutoff dependences of the bare theory, which renders the question of the renormalizability even more interesting in the present context.
In this paper we will apply two different nonperturbative approximations to the QED of 3D Dirac and Weyl semimetals, namely the ladder approach and the large-N approximation (N being the number of different fermion flavors), which can be viewed as complementary computational methods. In both cases, we will see that the field theory appears to be renormalizable, making possible to absorb all the dependences on the high-energy cutoff into a finite number of renormalization factors given only in terms of the renormalized coupling.
We will show that the QED of 3D Dirac semimetals has two competing effects at strong coupling. One of them is the tendency to chiral symmetry breaking and dynamical mass generation, which are analogous to the same phenomena arising in the conventional QED at strong coupling [20][21][22][23][24][25][26][27]. This trend is however outweighed by the strong suppression of electron quasiparticles that takes place at large N , leading then to a different type of critical point at sufficiently large interaction strength [28], shared also by the 3D Weyl semimetals. We will see that the nonperturbative approaches applied in the paper afford a precise characterization of the respective critical behaviors in terms of the anomalous dimensions of relevant operators. At the end, the phase diagram of the 3D Dirac semimetals turns out to be richer than that of their 2D counterparts, displaying a transition to a phase with non-Fermi liquid behavior (genuine also of the 3D Weyl semimetals) which may be observed in materials hosting a sufficiently large number of Dirac or Weyl points.

Field theory of 3D Dirac and Weyl semimetals
Our starting point is the field theory describing the 3D fermions with linear dependence of the energy ε on momentum p, ε(p) = v F |p|, and interacting through the scalar part φ of the electromagnetic potential. The fermion modes can be encoded into a number N of spinor fields ψ i , i = 1, . . . N , where the index i represents in general the way in which the modes are distributed into different Dirac or Weyl points in momentum space. For convenience, we will choose to pair the fermion chiralities into four-component spinors, extending the notation to think of each ψ i field as being made of two different chiralities from a given Dirac point, or from two different Weyl points in the case of a Weyl semimetal. The hamiltonian can be written then as In order to complete the field theory, we have to provide the propagator of the scalar field φ, for which we take the Lorentz gauge in the dynamics of the full electromagnetic potential. Thus, we get from the original relativistic theory We are interested in systems where the Fermi velocity v F is much smaller than the speed of light c, which makes possible to adopt a simpler description taking the limit c → ∞. This justifies that we can neglect magnetic interactions between our nonrelativistic fermions, and that we can just deal with the zero-frequency limit of (2.2). Thus, the free propagator for the scalar potential will be taken in what follows as A remarkable property is that the hamiltonian (2.1) together with the interaction mediated by (2.3) define an interacting field theory that is scale invariant at the classical level. That is, under a change in the scale of the space and time variables the fields can be taken to transform accordingly as This gives rise to a homogeneous scaling of the hamiltonian (2.1), consistent with the transformation of an energy variable. Such a scale invariance has important consequences, since it makes possible to give meaning to the divergences that the field theory develops at high energies. The computation of the quantum corrections requires indeed the introduction of a high-energy cutoff that spoils the classical scale invariance. However, provided the field theory is renormalizable, the singular dependences on the cutoff can be still absorbed into a finite number of renormalization factors redefining the parameters of the theory. Writing the action S corresponding to the hamiltonian (2.1), we may start with a formulation introducing suitable renormalization factors Z ψ , Z v , Z e and Z φ : The gauge invariance of the model implies that Z e Z φ = 1. The point is that, if the field theory is well-behaved, one has to be able to render all the physical observables cutoffindependent, when written in terms of the renormalized parameters v R and e.
The first example of the need to implement a high-energy cutoff can be taken from the computation of the electron self-energy. The first perturbative contribution is given by the first rainbow diagram at the right-hand-side in Fig. 1, which corresponds to the expression This contribution to the electron self-energy does not depend on frequency, but shows a divergent behavior at the upper end of the momentum integration, pointing at the renormalization of the Fermi velocity v F . In what follows, instead of dealing with a cutoff in momentum space, we will prefer to use a regularization method that is able to preserve the gauge invariance of the model. For this purpose, we will adopt the analytic continuation of the momentum integrals to dimension D = 3 − ǫ, which will convert all the high-energy divergences into poles in the ǫ parameter [29].
In the dimensional regularization method, the bare charge e 0 must get dimensions through an auxiliary momentum scale µ, being related to the physical dimensionless charge e by e 0 = µ ǫ/2 e (2.8) The computation of the first-order contribution to the self-energy (2.7) gives then The expression (2.9) shows the development of a 1/ǫ pole as ǫ → 0. This divergence can be indeed absorbed into a renormalization of the Fermi velocity v F , as the full fermion propagator G(k, ω) is related to the self-energy Σ(k, ω) by Thus, to first order in e 2 , the fermion propagator can be made finite in the limit ǫ → 0 by taking This renormalization of the fermion propagator already gives an idea of the effective dependence of the Fermi velocity on the energy scale. The bare Fermi velocity Z v v R cannot depend on the momentum scale µ, since the original theory does not know about that auxiliary scale. We have therefore For the same reason, the bare coupling e 0 cannot depend on µ, which leads to the equation Combining (2.12) and (2.13), we arrive at the scaling equation This is the expression of the growth of the renormalized Fermi velocity v R in the lowenergy limit, approached here as µ → 0. This scaling, which parallels the behavior of the Fermi velocity in 2D Dirac semimetals, is the main physical property deriving from the lowest-order renormalization of the 3D Dirac and Weyl semimetals [30].
To close this section, we comment that the result (2.9) for the electron self-energy has actually a wider range of validity, beyond the lowest-order approximation in which it has been obtained. This can be seen by exploiting a feature that is also shared with the 2D Dirac semimetals in the so-called ladder approximation, that is when corrections to the interaction vertex and the interaction propagator are neglected in the electron self-energy [18]. Then we remain with the sum of diagrams encoded in the self-consistent equation represented in Fig. 1. The electron self-energy in this approximation, Σ ladder (k), must have the form and therefore it is bound to satisfy the equation It is now easy to see that the second term at the right-hand-side of (2.16) identically vanishes. By performing a Wick rotation ω p = iω p , we get dω p 2π showing that Σ ladder (k) = Σ 1 (k) (2.18) The result (2.18) has important consequences for the ladder approximation, since it implies that the low-energy scaling of the Fermi velocity is again the main physical effect derived from the electron self-energy, limiting the possible electronic instabilities that can arise in that nonperturbative approach. = + Figure 1. Diagrammatic representation of the ladder approximation for the electron self-energy.

3D Dirac semimetals in ladder approximation
The iterative sequence of interactions, illustrated in Fig. 1 for the electron self-energy, defines a partial sum of perturbation theory which can be also applied to analyze the renormalization of composite operators. This provides a suitable approach to investigate possible instabilities of the electron system, since some of those operators correspond to order parameters characterizing the breakdown of respective symmetries of the field theory. Once the renormalization program is accomplished, the singularities found in the correlators of the composite operators may signal the condensation of a given order parameter, pointing at the onset of a new phase of the electron system.
The most relevant composite operators are made of bilinears of the fermion fields, having associated vertex functions that can be used to measure different susceptibilities of the electron system. At this point, the discussion is confined to the 3D Dirac semimetals (excluding the case of Weyl semimetals), for which we can write different composite operators in terms of four-component spinors about the same Dirac point in momentum space. We pay attention in particular to those operators corresponding to the charge, the current, and the fermion mass, given by The respective one-particle-irreducible (1PI) vertices are The point is that, assuming the renormalizability of the field theory, there must exist a multiplicative renormalization of the vertices that renders Γ 0,ren , Γ c,ren and Γ m,ren independent of the high-energy cutoff. We check next that property of the field theory, dealing with the same iteration of the interaction applied before to the electron self-energy. In the present context, this defines the ladder approximation as given by the self-consistent diagrammatic equation represented for a generic vertex in Fig. 2.
We first observe that the renormalization of Γ 0 must be dictated by the gauge invariance of the model, since the composite operator ρ 0 (r) already appears in the action (2.6). This implies the result meaning that, according to the previous analysis of the self-energy, it must be Z 0 = 1 in the ladder approximation. Moreover, the composite operator ρ c may be introduced in the action multiplying it by the vector potential, showing that its renormalization can be related by gauge invariance to that of the kinetic term in (2.1). We conclude therefore that which, in the ladder approximation, leads to Z c = Z v . The results (3.10) and (3.11) can be also obtained more formally from the Ward idendities that reflect the gauge invariance of the model at the quantum level. They lead in particular to the equations These identities were used in Ref. [18] to check the suitability of the dimensional regularization method to preserve the gauge invariance of the field theory of 2D Dirac semimetals. It was shown there that, if the the electron self-energy is computed with the set of diagrams encoded in the equation of Fig. 1, then the Ward identities are satisfied when one takes for the vertex the series of ladder diagrams, but dressed with the lowest-order electron self-energy correction. With this scope of the nonperturbative approach, it can be also shown by direct renormalization of Γ c that Z c becomes as simple as the expression given by (2.11), which arises then from a remarkable cancellation of different contributions in the computation of the vertex. We conclude that the charge operator is not renormalized, while the current operator requires a multiplicative redefinition by the simple factor (2.11) in the ladder approximation. The absence of a singularity for any value of the coupling constant excludes the possibility of having an instability related to the condensation of the charge or the current in the electron system. We turn now to the case of the mass vertex (3.6), whose renormalization is not dictated by the previous analysis of the electron self-energy.

Mass vertex and dynamical mass generation
Our starting point for the analysis of the mass vertex in the ladder approximation is the self-consistent equation represented in Fig. 2. We are going to be mainly interested in the renormalization of Γ m and, for that purpose, it is enough to consider the limit of momentum transfer q → 0 and ω q → 0. The vertex must satisfy then the equation The resolution of (3.14) implies that Γ m has to be proportional to the identity matrix. Moreover, it turns out to be a function independent of the frequency ω k in this approximation. After a little algebra, we arrive then at the simplified equation Eq. (3.15) can be solved by means of an iterative procedure, expressing the vertex as a power series in the effective coupling λ 0 = e 2 0 /4πv F , Assuming that the term Γ (n) m (k) is proportional to 1/|k| nǫ , the next order can be computed consistently from (3.15), taking into account that Thus, the vertex can be written as an expansion where each order can be obtained from the previous one according to the recurrence relation We observe from (3.18) and (3.19) that Γ m develops higher order poles in the ǫ parameter as one progresses in the perturbative expansion. At this point, one has to check the renormalizability of the theory by ensuring that all the poles can be reabsorbed by means of a redefinition of the vertex like that in (3.9). In terms of the dimensionless coupling the renormalization factor Z m must have the structure with residues d i depending only on λ.
In the present case, it may be actually seen that the renormalized vertex Γ m,ren can be made finite in the limit ǫ → 0, with an appropriate choice of the functions d i (λ). We obtain for instance for the first perturbative orders An important point about the residues d i (λ) is that they can be chosen without having any dependence on the momentum k of the vertex. This is a signature of the renormalizability of the theory, by which we are able to render it independent of the high-energy cutoff, redefining just a finite number of local operators. On the other hand, the renormalization factor Z m contains important information about the behavior of the theory in the low-energy limit. This stems from the anomalous scale dependence that the vertex gets as a consequence of its renormalization [31]. The bare unrenormalized theory does not know about the momentum scale µ, but the factor Z m lends to Γ m,ren an anomalous scaling of the form The anomalous dimension γ m can be thus obtained as The renormalization factor Z m is given by an infinite series of poles in the ǫ parameter, and then it is highly nontrivial that the computation of γ m from (3.30) may provide a finite result in the limit ǫ → 0. At this stage, the dependence of Z m on µ arises from the scaling (2.13), since no self-energy corrections are taken into account yet. That leads to the equation The term with no poles in the ǫ parameter already gives the result for the anomalous dimension But it remains however to ensure that the rest of the poles identically cancel in Eq. (3.31), which requires the fulfillment of the hierarchy of equations Quite remarkably, we have checked that the set of equations (3.33) is satisfied, up to the order we have been able to compute the residues d i (λ) in the ladder approximation. In this task, we have managed to obtain numerically the first residues up to order λ 28 , with a precision of 36 digits. Besides verifying the hierarchy (3.33) at this level of approximation, we have also analyzed the trend of the power series for such large perturbative orders, finding that a singularity must exist at a certain critical coupling. Concentrating in particular on the function d 1 (λ), we have found that the terms in its perturbative expansion approach a geometric sequence at large n. The values we have obtained for the coefficients d (n) 1 are represented in Fig. 3 up to order n = 28. From the precise numerical computation of the coefficients, it is actually possible to verify that their ratio follows very accurately the asymptotic behavior From the estimate of the limit value d, we have obtained the finite radius of convergence for a value of the coupling constant The coupling λ c corresponds to a point where the anomalous dimension γ m diverges, which signals in turn the development of a singular susceptibility with respect to the mass parameter. That is, λ c has to be viewed as a critical coupling above which the electron system enters a new phase with dynamical generation of mass [32]. This characterization is also supported by the fact that, in the case of the 2D Dirac semimetals, a similar renormalization in the ladder approximation [18] has shown to lead to a value of the critical coupling that coincides very precisely with the point for dynamical mass generation determined from the resolution of the mass gap equation [19]. It is indeed remarkable that these quite different approaches give an identical result for the point of chiral symmetry breaking in the electron system. This reassures the predictability of our renormalization approach, that moreover allows to extend the analysis beyond the ladder approximation as we discuss in what follows.

Electron self-energy corrections to the ladder approximation
The most relevant way to improve the ladder approximation corresponds to including electron self-energy corrections in the diagrams encoded in the equation of Fig. 2. With these additional contributions we may account for an important feature of the electron system, which is the growth of the Fermi velocity in the low-energy limit. This leads to a reduction of the effective interaction strength, from which we can expect the need to push the nominal coupling to larger values in order to reach the phase with dynamical mass generation.
A fully consistent approach can be devised by considering the self-energy contributions in the same ladder approximation discussed in Sec. II. In this case we know that the electron self-energy coincides with the lowest-order result given by Eq. (2.9). This can be translated into a redefinition of v F , leading in the fermion propagator to an effective Fermi velocitỹ The electron self-energy corrections are then accounted for automatically after replacing the parameter v F byṽ F (k) in the self-consistent equation (3.15) for the mass vertex. It can be easily checked that the perturbative expansion of 1/ṽ F (k) introduced in that equation corresponds to the iteration of self-energy rainbow diagrams correcting the fermion propagators in the original ladder approximation. As in the previous subsection, we can assume that the mass vertex admits now the expansion We can apply a recurrent procedure as before to compute the different orders in (3.39), with a main difference in that now each t n is going to depend on all the lower orders in the expansion. This is so as the n-th order, when inserted at the improved right-hand-side of (3.15), can give rise to contributions to any higher order as the factor 1/ṽ F (k) is also expanded. At the end, we arrive at the result that We have to check again that the vertex can be made finite in the limit ǫ → 0 by a suitable multiplicative renormalization. In this case, we have first to express all quantities in terms of the renormalized Fermi velocity v R arising after subtraction of the pole at the right-hand-side of (3.37).
where Z v has already appeared in (2.11). The renormalized coupling is given now by When expressed in terms of λ, the new renormalization factor Z m must have the structure The first terms in the expansions of the residuesd i (λ) can be obtained analytically, with the result that The important point is again that the functionsd i (λ) can be chosen with no nonlocal dependence (in fact, with no dependence) on the external momentum of the vertex, which means that the renormalization can be accomplished with the redefinition of a finite number of purely local operators. The expansion of the residued 1 (λ) allows us to estimate the effect of the self-energy corrections on the anomalous dimension of the vertex. This is given as before by Eq. (3.30), while now we have to account for the change in the scaling of the renormalized coupling λ. This can be obtained from Eqs. (2.13) and (2.14), with the result that The computation of the anomalous dimension proceeds then as Eq. (3.51) provides an expression for γ m that contains poles of all orders in the ǫ parameter. As in the previous subsection, it is therefore highly nontrivial that a finite result may be obtained for the anomalous dimension in the limit ǫ → 0. From (3.51), the equation that we have to inspect in this case is The zeroth order in ǫ leads to the equation for the anomalous dimension  However, this derivation makes sense only when the cancellation of all the poles in (3.52) can be guaranteed, which implies the set of equations It is reassuring to see that the analytic expressions given in (3.44)-(3.49) satisfy the hierarchy of equations (3.54). A more comprehensive check can be performed however with the numerical computation of the expansions of the first residues, that we have carried out up to order λ 30 and with a precision of 40 digits. In this way, we have been able to certify that the conditions (3.54) are also fulfilled at that level of approximation.
The residued 1 (λ) has in particular an expansioñ with coefficients that have been represented in Fig. 4 up to order n = 30. We observe that thed (n) 1 series approaches a geometric sequence at large n. The ratio between consecutive orders can be fitted indeed with great accuracy by a behavior like (3.35), allowing to compute the asymptotic valued with an error estimated to be (as in (3.36)) in the last digit. As remarked before, the critical coupling λ c corresponds to the point at which the anomalous dimension γ m diverges. This is in turn the signature of the development of a nonvanishing expectation value of the mass operator (3.3). Thus we see that, even after taking into account the effect of the Fermi velocity renormalization, there still remains a strong-coupling phase in the electron system characterized by the dynamical generation of mass for the Dirac fermions. The value of the critical coupling (3.57) is sensibly larger than that obtained in the previous subsection, which is consistent with the fact that the self-energy corrections effectively reduce the interaction strength. The situation is in this respect rather similar to the case of the 2D Dirac semimetals, where the effective growth of the Fermi velocity at low energies has been invoked as the reason why no gap is observed in the electronic spectrum of graphene [33,34], even in the free-standing samples prepared at very low doping levels about the Dirac point [17].

Large-N approximation for 3D Dirac and Weyl semimetals
We deal next with an approach complementary to the ladder approximation, paying attention to the effect of the photon self-energy corrections on the electron quasiparticles. This will take into account the renormalization of the electron charge, which is a relevant feature in the field theory of 3D semimetals as well as in the fully relativistic 3D QED [35]. In order to include a consistent set of quantum corrections, we will dress the interaction with the sum of bubble diagrams obtained by iteration of the electron-hole polarization, as represented in Fig. 5. This approximation can be then considered as the result of taking the leading order in a 1/N expansion, providing a resolution of the theory in a well-defined limit as N → ∞. The corrections to the bare interaction are represented by the polarization Π(q, ω q ), from which the full interaction propagator D(q, ω q ) is obtained according to In the present approximation, the polarization is given by the dominant contribution in the large-N limit with the free Dirac propagator Computing the integrals in analytic continuation to D = 3 − ǫ, we get the result The polarization (4.4) diverges in the limit ǫ → 0, which points at the need to renormalize the scalar field φ mediating the Coulomb interaction. As already mentioned, the gauge invariance of the action (2.6) implies that Z e Z φ = 1, which means that the electron charge is consequently renormalized. In the present approximation, these effects can be taken into account at once in the large-N expression of the electron self-energy Thus, we can determine the electron charge renormalization by devising a finite limit of the effective e-e interaction in (4.6) as ǫ → 0. This can be achieved by absorbing the pole in (4.4) into the bare charge e 0 according to the redefinition In terms of the renormalized charge e, the electron self-energy becomes The renormalization program has to be completed anyhow by accounting for the divergences that arise when performing the integrals in (4.8). In order to make sense of the theory in the large-N limit, we can assume formally that N e 2 /2π 2 v F ∼ O(N 0 ). The difference between v F and its renormalized counterpart v R is then a quantity of order ∼ 1/N , and the electron self-energy gets a natural dependence on the effective coupling We may actually resort to a perturbative expansion It may be seen from (4.10) that, to order g n , the self-energy can develop divergences as large as ∼ 1/ǫ n . If the theory is renormalizable, it must be possible to absorb these poles into the definition of suitable renormalization factors in the full fermion propagator (2.10). We look then for the large-N limit of Z ψ and Z v , which must have the general structure with coefficients b i and c i depending only on the renormalized coupling g. It is not difficult to compute the integrals that are needed to get in general the n-th order of Σ(k, iω k ). We are interested in terms that are linear in ω k and k, and the only technical point is that these contributions must be also regularized in the infrared, using for instance the own external frequency or momentum, or some other suitable parameter. It can be checked that the choice of a particular infrared regulator is not important when extracting the high-energy divergences as ǫ → 0. For this reason, we have resorted to introduce a fictitious Dirac mass ν, which has the ability of keeping well-behaved both types of contributions linear in ω k and k. We rely then on the results for the generic integrals and With the help of (4.13) and (4.14), one can obtain analytically for instance the first terms in the expansions of the residues in (4.11) and (4.12), by imposing that the renor-malized fermion propagator becomes finite in the limit ǫ → 0. We find in this way and for the Fermi velocity renormalization The important point is that these expansions of the residues do not show any dependence on the auxiliary scales µ and ν (or on the external frequency and momentum when these are used alternatively to regularize the self-energy in the infrared). This is a nice check of the renormalizability of the theory which guarantees that, at least in the large-N approximation, the high-energy divergences can be absorbed into a finite number of renormalization factors depending only on the renormalized coupling g.
The renormalization factors Z ψ and Z v provide important information about the behavior of the electron system in the low-energy limit. In the renormalized theory, the Dirac fermion field gets an anomalous scaling dimension γ ψ (g), which can be computed as [31] This anomalous dimension governs the scaling of the correlators involving the Dirac field. When sitting at a fixed-point of the renormalized parameters, we would get for instance for the Dirac propagator G(sk, sω) ≈ s −1+γ ψ G(k, ω) (4.27) On the other hand, the renormalized Fermi velocity v R gets also a scaling dependence which can be assessed in terms of a function γ v (g) such that Thus, the knowledge of γ ψ (g) and γ v (g) allows us to inspect the theory in the limit of long wavelengths and low energies as µ → 0. The anomalous dimension γ ψ (g) can be computed from the dependence of the renormalized coupling g on the auxiliary scale µ, taking into account that The scaling of g can be obtained in turn by differentiating the expression (4.7) with respect to µ, bearing in mind the independence of e 0 with respect to that auxiliary parameter. This leads to µ ∂ ∂µ As already mentioned, the difference between v F and v R can be taken formally as a quantity of order ∼ 1/N , so that we end up in the large-N limit with the equation This expression can be plugged into (4.29) to get Working to leading order in the 1/N expansion, we can set Z ψ = 1 at the right-handside of (4.32). The term with no poles in that equation leads to the result This expression of the anomalous dimension makes only sense, however, provided that one can certify the cancellation of the rest of terms carrying poles of all orders in the ǫ parameter, which implies the hierarchy of equations It is remarkable that the power series of the residues given in (4.15)-(4.19) satisfy indeed the hierarchy (4.34). Beyond the analytic approach, we have computed numerically the expansion of the first residues c i (g) up to order g 32 , with a precision of 40 digits. In this way, we have been able to check that the conditions (4.34) are fulfilled at that level of approximation, reassuring the consistency of the large-N approach for the present field theory.
Another important detail of the numerical calculation of the residues is the evidence that the coefficients of each perturbative expansion approach a geometric sequence for large orders of the coupling. In the case of the residue c 1 (g), we have for instance the power series with coefficients that we have represented in Fig. 6 up to the order we have carried out the numerical computation. The behavior observed for the series c implies that the perturbative expansion must have a finite radius of convergence. This can be obtained by approaching the ratio between consecutive coefficients according to the dependence which provides indeed an excellent fit at large n. We get in this way the estimate r ≈ 0.3333333 (4.37) where the error lies in the last digit. We find then the singular point for the coupling g (the radius of convergence) at The critical value g c corresponds to a point where the anomalous scaling dimension γ ψ (g) diverges. This can be appreciated in Fig. 7, where we have represented that function from the results of our numerical calculation. The singularity found at g c has a precise physical meaning, since γ ψ governs the scaling of all the correlators involving the Dirac fermion field. Away from a fixed-point in the renormalized parameters, the scaling is not so simple as in (4.27), but from that expression we already get the idea that the divergence of γ ψ implies the decay of the Dirac propagator [36]. In the limit g → g c , the singularity in the scaling dimension leads indeed to the suppression of the quasiparticle weight. This characterizes a particular form of correlated behavior, which has been identified in several other instances of interacting electrons and has led to constitute the class of socalled non-Fermi liquids [37][38][39][40]. The distinctive feature of this class is the absence of a quasiparticle pole in the electron propagator, which leads to a very appealing paradigm to explain unconventional properties like those found in the normal state of copper-oxide superconductors.
In order to ensure the relevance of the divergence at g c , we have anyhow to verify that such a singular behavior is not prevented by some other instability in the low-energy scaling of the electron system. We pay attention in particular to the renormalization of the Fermi velocity, which has a natural tendency to grow in the low-energy limit. The function γ v (g) that dictates the scaling in Eq. (4.28) can be obtained from the independence of the bare Fermi velocity on the auxiliary scale µ: The dependence of Z v on µ comes only from the coupling g so that, using Eq. (4.31), we get Working to leading order in the 1/N expansion, the term free of poles in Eq. (4.40) leads to the result As in the case of the anomalous scaling dimension, one has to make sure however that the rest of terms carrying poles of all orders in Eq. (4.40) cancel out identically, which is enforced by the conditions The expressions of the residues in (4.20)-(4.25) satisfy indeed the set of equations (4.42), and a more extensive numerical computation confirms that this is also the case for the power series of the b i (g) evaluated up to order g 32 . This analysis also shows that these expansions do not lead to any singularity in the range of couplings up to the critical point g c . The corresponding function γ v (g) obtained from b 1 (g) is represented in Fig. 8. The regular behavior observed in the plot implies that the scaling of the Fermi velocity is not an obstacle for the suppression of the fermion quasiparticles as the interaction strength becomes sufficiently large to hit the critical point at g c .
As a last check, we look also for the possible tendency towards chiral symmetry breaking and dynamical mass generation in the large-N approximation. For this purpose, we may analyze the scaling of the vertex for the mass operator represented in Fig. 9. Computing in the limit where both frequency and momentum transfer vanish, we get We have to account as before for the renormalization of the charge, which leads to an expression in terms of the renormalized coupling g ψψ Figure 9. Corrections to the vertex for the mass operator in the large-N approximation. The cross represents the operator ψψ and the thick wiggly line stands for the dressed interaction propagator as defined in Fig. 5.
As already done in the case of the electron self-energy, we can compute the high-energy divergences of the vertex in the limit of vanishing k and ω k , regularizing the integrals in the infrared with an auxiliary mass ν in the Dirac propagators. In this way, the different terms in the expansion (4.44) can be obtained using the general result We see that the vertex can develop in general divergences of order ∼ 1/ǫ n at the level g n in the perturbative expansion. These have to be reabsorbed into a multiplicative renormalization of the type shown in (3.9). The point to bear in mind is that now Z m is made of two different factors, coming independently from the renormalization of the composite mass operator (3.3) and the Dirac fermion fields in the vertex [31]: The renormalization factor Z ψ 2 for the mass operator can have the general structure Under the assumption that the theory is renormalizable, it must be possible to choose appropriately the residuesd i (g) to end up with a renormalized vertex, finite in the limit ǫ → 0, given by Γ m,ren = Z ψ Z ψ 2 Γ m (4.48) We have checked that the vertex (4.48) can be made free of poles, at least up to order g 32 we have been able to compute it numerically, with a set of residuesd i (g) that depend only on the renormalized coupling g. We have also seen that these functions have a regular behavior in the range up to the critical coupling g c where the anomalous dimension γ ψ diverges. This makes clear that the dominant instability in the large-N approximation is indeed characterized by the suppression of the electron quasiparticles.
Of all the residuesd i (g), the first of them conveys the most relevant piece of information, since it is related to the anomalous dimension of the composite mass operator, defined by Paralleling the above derivation of similar scaling dimensions, one arrives at the result The plot of this function obtained from our numerical computation ofd 1 (g) is represented in Fig. 10. The regular behavior observed in the figure is the evidence that no singularity can be expected in the correlators of the mass operator ρ m , whose magnitude has got to be bound to the scaling dictated by the anomalous dimension γ ψ 2 in the low-energy limit.

Conclusions
We have studied the development of the strong-coupling phases in the QED of 3D Dirac and Weyl semimetals by means of two different nonperturbative approaches, consisting in the sum of ladder diagrams on one hand, and taking the limit of large number of fermion flavors on the other hand. We have benefited from the renormalizability that the theory shows in both cases, which makes possible to render all the renormalized quantities independent of the high-energy cutoff. Thus we have been able to compute the anomalous scaling dimensions of a number of operators exclusively in terms of the renormalized coupling constant, which has allowed us to determine the precise location of the singularities signaling the onset of the strong-coupling phases.
We have seen that, in the ladder approximation, the 3D Dirac semimetals have a critical point marking the boundary of a phase with chiral symmetry breaking and dynamical generation of mass, in analogy with the same strong-coupling phenomenon in conventional QED [20][21][22][23][24][25][26][27]. We have found however that such a breakdown of symmetry does not persist in the large-N limit of the theory, which is instead characterized by the growth of the anomalous dimension of the electron quasiparticles at large interaction strength. The picture that emerges by combining the results from the two nonperturbative approaches is that chiral symmetry breaking must govern the strong-coupling physics of the 3D Dirac semimetals for sufficiently small N , while there has to be a transition to a different phase with strong suppression of electron quasiparticles prevailing above a certain value of N .
With the results obtained for the different critical points we can draw an approximate phase diagram in terms of the number N of Dirac fermions and the renormalized coupling g defined in (4.9). In principle, the critical coupling (4.38) can provide a reliable estimate for the onset of non-Fermi liquid behavior at sufficiently large N . On the other hand, the critical value of g deriving from (3.57) may lead to a sensible map of the phase with chiral symmetry breaking as long as its magnitude does not become larger than that of (4.38). The resulting phase diagram of the 3D Dirac semimetals, represented in Fig. 11, shows the regions where the behavior of the phase boundaries may be captured by our alternative approximations, away from the intermediate regime about N = 4 where the competition between the two strong-coupling phases cannot be reliably described within our analytic framework. It is interesting to compare at this point the phase diagram in Fig. 11 with that obtained with the resolution of the Schwinger-Dyson equations, which can be trusted for all values of N . We can see from the results of Ref. [34] that such a numerical approach sets indeed the interplay between the phases with chiral symmetry breaking and non-Fermi liquid behavior at N = 4. For N > 4, we observe that the critical line for that latter phase is not straight, although it approaches a constant asymptotic limit at large N . The selfconsistent resolution leads also to a critical line for chiral symmetry breaking in the regime N ≤ 4 that has an approximate linear behavior as a function of N . It is worth to mention that the values of the critical couplings found in the numerical resolution of Ref. [34] are at first sight much larger than their counterparts in the diagram of Fig. 11. We have to bear in mind, however, that the critical values in that paper were given for the bare couplings, in the theory with a finite high-energy cutoff, while critical points like (3.57) and (4.38) are referred in the present context to the renormalized couplings. Overall, we may conclude that there is good qualitative agreement between the location of the phases in the diagram of Fig. 11 and in the more complete map obtained from the resolution of the Schwinger-Dyson equations.
We end up remarking that our analysis can be applied to characterize not only the strong-coupling regime of the 3D Dirac semimetals, but also of the Weyl semimetals. In these materials, each Weyl point hosts a fermion with a definite chirality, which means that the electron system cannot undergo chiral symmetry breaking through condensation of some order parameter at zero momentum [41]. This obstruction does not hold however for the strong-coupling phase that we have identified in terms of the suppression of fermion quasiparticles. It is clear that the large-N approach of Sec. IV applies equally well for Weyl and Dirac fermions, so that Weyl semimetals are susceptible of developing the phase with non-Fermi liquid behavior that we have mapped at large N in the diagram of Fig. 11. In this respect, there should be good prospects to observe such an unconventional behavior in present candidates for Weyl semimetals, like TaAs or the pyrochlore iridates, which have up to 12 pairs of Weyl points [5,6]. These considerations show that the strong-coupling phases studied here are not beyond reach, and that they may be actually realized in 3D Dirac or Weyl semimetals with suitably small values of the Fermi velocity or with the large number of fermion flavors already exhibited by known materials.