Two-loop mass anomalous dimension in reduced quantum electrodynamics and application to dynamical fermion mass generation

We consider reduced quantum electrodynamics (RQED$_{d_\gamma,d_e}$) a model describing fermions in a $d_e$-dimensional space-time and interacting via the exchange of massless bosons in $d_\gamma$-dimensions ($d_e \leq d_\gamma$). We compute the two-loop mass anomalous dimension, $\gamma_m$, in general RQED$_{4,d_e}$ with applications to RQED$_{4,3}$ and QED$_4$. We then proceed on studying dynamical (parity-even) fermion mass generation in RQED$_{4,d_e}$ by constructing a fully gauge-invariant gap equation for RQED$_{4,d_e}$ with $\gamma_m$ as the only input. This equation allows for a straightforward analytic computation of the gauge-invariant critical coupling constant, $\alpha_c$, which is such that a dynamical mass is generated for $\alpha_r>\alpha_c$, where $\alpha_r$ is the renormalized coupling constant, as well as the gauge-invariant critical number of fermion flavours, $N_c$, which is such that $\alpha_c \rightarrow \infty$ and a dynamical mass is generated for $N<N_c$. For RQED$_{4,3}$, our results are in perfect agreement with the more elaborate analysis based on the resolution of truncated Schwinger-Dyson equations at two-loop order. In the case of QED$_4$, our analytical results (that use state of the art five-loop expression for $\gamma_m$) are in good quantitative agreement with those obtained from numerical approaches.


Introduction
Reduced quantum electrodynamics (RQED dγ ,de ) is a model describing relativistic fermions in a d e -dimensional space-time and interacting via the exchange of massless bosons in d γdimensions (d e ≤ d γ ). It has been introduced in [1] motivated by the study of dynamical chiral symmetry breaking in brane-world theories (see also [2]). Soon after, a first application was devoted to the specially important case of conformal RQED 4,3 (also known as pseudo-QED from [3] and mixed-dimensional QED from the recent [4]) in relation with graphene [5]. More precisely, RQED 4,3 describes graphene [6] at its infra-red (IR) Lorentz invariant fixed point [7] (see also [8][9][10][11] for reviews on graphene). It has been shown in [12] that there is a mapping between RQED 4,3 and the well known three-dimensional QED (QED 3 ) in the large-N limit (where N is the number of 4-component Dirac spinors -see also [13] for a review on large-N quantum field theories) which is a celebrated effective field theory for planar condensed matter physics systems exhibiting Dirac-like low-energy excitations such as, e.g., high-T c superconductors [14][15][16][17], and quantum antiferromagnets [18]. The model of RQED 4,3 also captures some universal features of a broader range of so-called Dirac liquids that have been discovered experimentally during the last decade and are under active study such as, e.g., artificial graphene-like materials [19], surface states of topological insulators [20] and half-filled fractional quantum Hall systems [21].
Theoretically, there has been rather extensive studies of RQED 4,3 during the last decade with primary applications to Dirac liquids, e.g., their transport and spectral properties [22][23][24][25][26][27][28], quantum Hall effect [4,29] (in [4] the model was invoked as an effective field theory describing half-filled fractional quantum Hall systems) and dynamical symmetry breaking [12,30] which will be our main focus in this paper (see also [31] for a review on these topics). From a more field theoretic point of view, the model was shown to be unitary [32], it's properties were studied under the Landau-Khalatnikov-Fradkin transformation [33,34] as well as under duality transformations [35]. Renewed interest in the model and its formal properties was triggered by the study [36] of interacting boundary conformal field theories, see, e.g., [37][38][39][40][41] and supersymmetric extensions constructed and analyzed in [42,43].
In this paper we shall focus on (parity-conserving) dynamical fermion mass generation in (massless) RQED 4,de (d e ≤ 4) that include renormalizable models such as QED 4 and RQED 4,3 with dimensionless (renormalized) coupling constant α r (the fine structure constant). Let's recall that this non-perturbative mechanism relies on the existence of a critical coupling constant, α c , which is such that for α r < α c fermions are massless (assuming chiral invariance), while for α r > α c a dynamical fermion mass is generated at a given number of fermion flavours N (and chiral symmetry is dynamically broken). Alternatively, in the limit α c → ∞, a dynamical mass is generated for N < N c where N c is the critical number of fermion flavours. The knowledge of α c and N c therefore provides precious information on the phase structure of gauge theories. Moreover, in the condensed matter terminology, this mechanism corresponds to a (semi-)metal to insulator transition whereby a dynamical gap is generated at strong coupling. It's study is crucial, e.g., for the development of graphene-based transistors [44].
In the last four decades, the standard approaches to study dynamical mass generation in gauge theories were either based on lattice simulations [45][46][47][48][49][50] or solving Schwinger-Dyson (SD) equations on which we shall focus here, see [51] for an early detailed review as well as the manuscripts [52,53]. While initial interests were towards four-dimensional theories, the importance of QED 3 (and its large-N limit with no running coupling) was recognized very early [54][55][56] because of it's simpler UV structure 1 and similarity to quantum chromodynamics (QCD). Nevertheless, because of the formidable complexity of the task, calculations were generally carried out only at the leading order (LO) in the coupling constant [57] or by including non-perturbative ansätze for the vertex function in one-loop like SD equations [58][59][60]. Often, the resulting solution displayed residual gauge variance which is unsatisfactory for a physical quantity such as a critical coupling. Following early 1 Because of the absence of running of their coupling constants, large-N QED3 and RQED4,3 may be referred to as "standing" gauge theories. In contrast, QED4 has a running of the coupling constant. It is only in the quenched approximation (where fermion loops are neglected) that the coupling constant of QED4 does not renormalize. multi-loop works of Nash [61] and Kotikov [62,63], a complete gauge-invariant prescription up to next-to-leading (NLO) of the 1/N -expansion for QED 3 appeared only recently in [64] and [65,66] (see also [67] for a recent review). In [12] the results were then mapped to RQED 4,3 thereby extending the LO results of [1] to the NLO in α.
The gauge-invariant prescriptions found in [64] and [66] for large-N QED 3 alleviate doubts about the validity of the SD equation approach though a similar prescription still has to be implemented for QED 4 . Considering the complexity of the calculations, simpler approaches are worthwhile investigating. An argument often invoked in the recent literature on QED 3 , see, e.g., [68][69][70][71][72][73][74][75][76][77] (see also [78] for a review), is the fact that a fermion quadrilinear operator becomes relevant at criticality (this has actually been noted in the early literature on four-dimensional models, see, e.g., [79][80][81][82][83]). The computation of the anomalous dimension of the corresponding composite operator allows then to derive a marginality crossing equation -as referred to in [76][77][78] -in order to extract the value of the critical coupling. Actually, as noticed in [64], the SD gap equation incorporates such a criterion in a rather simple and efficient way, i.e., via the mass anomalous dimension (the anomalous dimension of the fermion bilinear mass operator) which is a gauge-invariant quantity governing the ultra-violet (UV) asymptotic behaviour of the fermion propagator [84][85][86][87] (see also the textbook [88]).
In the present paper we will compute the mass anomalous dimension, γ m , up to two loops, in a general theory of RQED 4,de generalizing recent results in RQED 4,3 [40] and recovering well known ones in QED 4 (see the lectures [89]). Following [64], we then construct a gap equation that allows us to study the critical properties of RQED 4,de . In the case of RQED 4,3 , we straightforwardly recover the results of [12] for the NLO critical coupling and flavour number. As for QED 4 , our approach is semi-phenomenological, but our results are in good quantitative agreement with those obtained from numerical solutions of SD equations, see, e.g., [52,90].
The paper is organized as follows. In Sec. 2 we present the model, the perturbative setup as well as our renormalization conventions. In Sec. 3, we present the one and two-loop calculations of γ m in RQED 4,de . The critical properties of the model are then analysed in Sec. 4 where the gap equation is derived and applications to RQED 4,3 and QED 4 (both in the quenched and unquenched cases) are provided. The conclusion is given in Sec. 5. Some conventions and results related to the master integrals entering our calculations are provided in App. A and exact results for all the computed diagrams are presented in App. B.

(2.4)
Upon turning on interactions, the dressed fermion propagator reads: where the fermion self-energy has the following parameterization appropriate to the massive case: and its vector and scalar parts can be extracted with the help of On the other hand, the dressed photon propagator can be expressed as: where the photon self-energy is parameterized as: .
Anticipating the RPA analysis of sec. 4.4.2, let's note that, in the case of RQED 4,3 , the transverse part of the photon propagator, (2.8b), reads: which can be usefully contrasted with the QED 4 case: With this perturbative setup, the renormalization constants are defined as: 12) and the following conventions will be used: where γ E is Euler's constant, µ is the renormalization scale appropriate to the modified minimal subtraction (MS) scheme in which all our calculations will be performed andᾱ is the bare reduced coupling constant. Notice that, from Ward identities, Z α = Z −1 A and Z ξ = Z A . Moreover, Z m (which is of special interest to us here) and Z ψ can be computed from the fermion self-energy parts with the help of 14) The corresponding anomalous dimensions and beta function can then be calculated in a standard way with the help of Notice that, in the MS-scheme, the renormalization constants are Laurent series in ε γ . They can be written as: 3 Perturbative calculations

Preliminaries
In this section we compute γ m in RQED 4,de . From (2.17) and (2.14), this requires the knowledge of the fermion self-energies Σ V and Σ S which can be extracted from (2.7). From [91], we recall that, in standard four-dimensional quantum field theories (that correspond to ε e = 0 and ε γ → 0), the anomalous dimensions cannot depend on external momenta and masses and, in the MS-scheme, on γ E and ζ 2 . Hence, as an IR rearrangement [91], we will compute the self-energies Σ V and Σ S in the massless limit, see [92] for a review on the corresponding techniques. Let's note that γ ψ was already computed in [25] for RQED 4,de and we shall follow the notations of that paper. In particular, the following parameters will be useful: where ψ(x) is the digamma function; see also App. A of [25] for the expansion of the relevant master integrals.
In the following, we will provide the ε γ -expansion of individual graphs contributing to Σ V and Σ S self-energies for arbitrary ε e (exact expressions for individual graphs contributing to Σ S and Σ V can be found in App. B of the present paper). 2 Combined with (2.14) and (2.17), this will allow us to compute the renormalization constants and anomalous dimensions in RQED 4,de . Actually, from the one-loop expansion of the polarization operator, we will see that the limits ε γ → 0 and ε e → 0 do not commute (a fact already noticed in [25]). In order to recover QED 4 from our general results, we will also present improved expressions for both γ ψ and γ m valid for RQED 4,de .

One-loop calculations
At one-loop, we have the following two simple contributions to the fermion and photon self-energies: that are displayed in Fig. 1. 2 The expansion of individual graphs contributing to ΣV in RQED 4,de was not explicit in [25]. For clarity and completeness, we provide such expressions here which can be usefully compared with the corresponding ones for ΣS.
We first consider the one-loop fermion self-energy (3.2a). Using the parameterization (2.7), we extract the scalar and vector parts that are then computed exactly (see (B.1) of App. B). The resulting ε γ -expansions at arbitrary ε e then read: whereL p was defined in (3.1a) and ψ (x) is the trigamma function. For completeness, we also provide these expansions in the RQED 4,3 case (ε γ → 0 and ε e = 1/2): With the help of (2.14), these results allow for a straightforward derivation of the one-loop renormalization constants where the parameterization (2.16) was used. From (2.17), the corresponding one-loop anomalous dimensions read: Similarly, we now consider the one-loop polarization operator (3.2b). Using the parameterization (2.9) and from the exact result (B.1c), it's ε γ -expansion at arbitrary ε e reads: which is valid only for d e = 3. In this case (ε e = 1/2), it is finite: which is related to the fact that the coupling constant does not run in RQED 4,3 . However, (3.7) does not reproduce the pole structure of QED 4 . The correct result is obtained from (3.2b) by first setting ε e = 0 and then performing the ε γ -expansion; this yields: This non-commutativity of the ε e and ε γ limits will appear in the two-loop fermion selfenergy graphs with a fermion loop insertion (see Fig. 2a) thereby affecting the anomalous dimensions. We will discuss improved expressions for them in Sec. 3.4.

Two-loop calculations
At two-loop, three diagrams contribute to the fermion self-energies: Fig. 2 where they are displayed.

Anomalous dimensions
Combining the above derived one-loop (3.6) and two-loop (3.16) contributions to the anomalous dimensions yields: where γ m is fully gauge invariant as it should be while the gauge-variance of γ ψ is at oneloop in accordance with the Landau-Khalatnikov-Fradkin transformation, see, e.g., [34]. Interestingly, (3.17a) and (3.17b) have similar structures. In particular, we see that in the limit of RQED 4,3 (ε e → 1/2) the factors of π 2 arise from K 1 , see (3.1a), and not from the ε γ -expansion of Γ-functions. It turns out however that, because of these K 1 terms, (3.17) (and similarly for the two-loop renormalization constants (3.15) and self-energies (3.11) and (3.13) above) do not apply to the case of QED 4 . As discussed below (3.7), this discrepancy (which is even more severe for γ m than γ ψ ) originates from the fact that the limits ε e → 0 and ε γ → 0 do not commute for the one-loop polarization operator and hence for (3.11c) and (3.13c) where it appears as a subgraph. 3 So (3.17) only apply to RQED 4,3 . The expressions (3.17) can be improved in order to cover both cases of QED 4 and RQED 4,3 at the expense of introducing an additional parameter z such that All calculations done, the two-loop contributions to the renormalization constants (3.15) now take the form: where, in QED 4 the running of the coupling constant at one-loop has been taken into account via: . Hence, we obtain one of the central results of this paper in the form of improved expressions for the anomalous dimensions in RQED 4,de : A remark is in order here in relation with (3.20b). Within the SD formalism, see, e.g., the review [51], a common procedure to simplify the equations and minimize the gauge-variance of the solutions is to consider the gauge for which the fermion anomalous dimension vanishes. Such a gauge is referred to as the "good gauge". At one-loop, from (3.20b), the "good gauge" is simply: which corresponds to the Landau gauge (ξ g = 0) in QED 4 and the Nash gauge (ξ g = 1/3, see [61]) in RQED 4,3 . Proceeding in a similar way at two-loop, with the help of (3.20b), we find that the "good gauge" becomes: (3.22) So, in the case of RQED 4,3 (z = 0 and ε e → 1/2), (3.22) yields: and in the case of QED 4 (z = 1 and ε e → 0), the (two-loop) "good gauge" reads: Note that, in the case of QED 4 , (3.24) is known since a long time (see, e.g., eq. (B.1) in [93] for an early derivation in the quenched case) and can presently be extended to 5 loops thanks to state of the art results for γ ψ [94,95].
We may now proceed in applying (3.19) and (3.20) to various cases of interest. Let's consider first the most important case of RQED 4,3 which amounts to set z = 0 and ε e → 1/2. From (3.5) and (3.19b), the renormalization constants read: And from (3.20), the anomalous dimensions read: where (3.26b) agrees with the result of [25]. Our result (3.26a) is new and corresponds to the parity-even mass anomalous dimension of RQED 4,3 for an arbitrary number N of 4-component spinors. 4 In order to get more confidence in our result for γ m , we proceed on showing that it can be mapped to the one of QED 3 at NLO in the large-N expansion with the help of [12]:ᾱ Let's note that, in [40], the parity-odd mass anomalous dimension of a single 2-component spinor was computed. The computation of such a quantity requires taking into account of the fact that the trace of three 2 × 2 Dirac gamma matrices is non-zero: Tr[γ µ γ ν γ ρ ] = 2i ε µνρ where ε µνρ is the rank 3 totally antisymmetric tensor. This affects the one-loop polarization operator and hence the coefficient of ζ2 in γm through the scalar part of diagram (a) of fig. 2 (and corresponding self-energy Σ2Sa(p 2 )). Recomputing this diagram with an arbitrary number, n, of 2-component spinors (which is such that n = 2N ) yields: where Tn arises from the epsilon tensors and is such that T 2k = 0 while T 2k+1 = 1 for any integer k. Eq. (3.27) is a general formula including both cases of an odd or an even number of 2-component spinors. Indeed, for n = 2N , we recover our eq. (3.26a) while for n = 1 we recover the result of [40] which, in our notations, reads: Note also that the model considered in [40] is that of a reflective boundary (say, with fine structure constant α bdry ) while our case corresponds to a transparent interface (with fine structure constant α); these models are related by: α bdry = α/2. 5 Notice that this result was also known to the authors of [62,63,97]. Though it did not appear explicitly in these papers, the knowledge ofΠ2 was required to perform the calculations carried out in these papers.
which corresponds exactly to the result obtained by Gracey in [98]. 6 Our results also allow us to consider the well-known case of QED 4 , see, e.g., the textbook [89], thereby strengthening their validity. From (3.5) and (3.19b), this amounts to set z = 1 and ε e → 0 yielding: Proceeding in a similar way from (3.20), the corresponding anomalous dimensions read: In (3.33a), the one-loop contribution was probably first computed in [93], the quenched two-loop one in [99] and the unquenched two-loop contribution can be extracted from the QCD calculations of [100,101] as we could learn from [102] where the three-loop QCD calculation was performed. Notice that, presently, γ m is known up to five loops in QCD [94,95] from which the corresponding QED result can be extracted (we shall use it in the next section). 7

Critical properties
As explained in the Introduction, a standard approach to study the critical properties of gauge theories is by solving the SD equations. Truncating these equations at LO is generally unsatisfactory as the resulting critical coupling may be strongly gauge-variant. The fully gauge-invariant procedure advocated in the recent [64] and [66] for large-N QED 3 requires computing NLO corrections and then performing a so-called Nash resummation (following the seminal work of Nash [61]) in order to properly cancel the gauge dependence both at LO and NLO. The simplicity of the resulting gap equation (from which the gauge-invariant 6 In addition, the mapping of [12] allows to recover (3.26b) from the field anomalous dimension of QED3 at NLO in the large-N expansion which can also be extracted from Gracey's work [98] and that we provide here for completeness: whereΠ2 was defined in (3.28). Hence, all results are in complete agreement. Note that from (3.30) we may also compute the"good gauge" in large-N QED3 at NLO that reads: in agreement with eq. (4.8) in [64]. 7 Let's note that (3.20) seems to also apply to the case of QED4,1 (z = 0 and εe → 3/2) for which γm = 0 + O(ᾱ 3 r ) and γ ψ = 2ᾱr (ξr − 3) + O(ᾱ 3 r ). This limit corresponds to the quantum mechanics of a point particle coupled to a fully retarded four-dimensional electromagnetic field. critical coupling is extracted) drastically contrasts with the complexity of the calculations that need to be performed in order to derive it. Following [64], in this section we provide a semi-phenomenological construction of the gap equation for RQED 4,de and apply it to cases of interest. Before that, we would like nevertheless to apply the SD formalism at LO for RQED 4,de . This will allow us to illustrate some of the difficulties we have just mentioned but also to underline the importance of the "good gauge" introduced in the previous section.

Leading order solution of SD equations for RQED 4,de
The SD equation for the fermion self-energy in RQED 4,de reads: where S(k), Γ µ and D µν (k) are the full fermion propagator, vertex function and photon propagator, respectively. In general, (4.1) is coupled to the Dyson equation for the vertex function and the photon propagator. A dynamical mass arises as a non-trivial solution of this system of coupled equations for Σ S (p 2 ) that needs to be determined self-consistently.
Focusing on the LO approximation in the coupling constant, all equations decouple as both the vertex function and the photon propagators are taken as the free ones: Γ µ = Γ µ 0 and D µν (p) = D µν 0 (p). Moreover, the wave-function renormalization is neglected, i.e., Σ V = 0, which implies that S(p) = i( / p − Σ S ) −1 where Σ S is now the dynamically generated parity-conserving mass (for convenience, the mass parameter has been absorbed in Σ S ). With the help of the parametrization (2.7) (with m = 1), (4.1) significantly simplifies and can be written as: At the critical point, (4.2) can be linearized by taking the limit Σ 2 S (k 2 ) k 2 and a powerlaw ansatz can be taken for the mass function Σ S (p 2 ) such that where the mass is assumed to be dynamical in origin and the index a has to be selfconsistently determined. As first noticed by Kotikov [62,63], together with (4.3), the linearized eq. (4.2) is now a massless integral which is very easily solved in dimensional regularization with the help of (A.1). The solution being finite for all d e , one can set ε γ → 0 and straightforwardly derive the LO gap equation: Solving it yields two values for the index a: (4.5) Dynamical mass generation takes place for complex values of a, i.e., for α > α c with Notice that, in the case of RQED 4,3 (ε e = 1/2), (4.6) leads to: α c = π/(2(5 + ξ)) in agreement with [12]. On the other hand, in the case of QED 4 (ε e = 0), we obtain α c = π/(3 + ξ), which exactly corresponds to the result of Rembiesa [103] according to [104]. As anticipated, (4.6) has a strong gauge dependence which is not satisfactory. In order to minimize it, we consider the "good gauge" which is given by (3.21) at one-loop. This yields: (4.7) We shall come back to this result and discuss it in detail in the next subsections. At this point, let's note that the LO gap equation itself, (4.4), can also be written in the "good gauge" where it may be expressed in the form: Interestingly, the right hand side of (4.8) involves the one-loop mass anomalous dimension (3.6) -a gauge-invariant quantity -as the only input. The powerful technique of Kotikov [62,63], that we have used here to easily solve the LO SD equation in dimensional regularization, can possibly be extended to higher orders for QED 4 along the lines of the recent [66]. Of course, an NLO computation is more complicated and out of the scope of the present paper. Instead, in the following subsection, we will present general arguments allowing us to build a fully gauge-invariant gap equation that is valid at any order in the coupling constant thereby generalizing (4.8).

Gap equation and criterion for dynamical mass generation
We start by recalling that, from the operator product expansion, the scalar part of the fermion self-energy has two asymptotes in the UV [84][85][86][87] (see also the textbook [88] sec. 12.3): where p 2 E = −p 2 is the Euclidean momentum, m is the bare mass and m dyn the dynamical one. As noticed in [87], dynamical mass generation arises from the coalescence of these two asymptotes. In particular, deep in the UV, p 2 E → ∞, the dynamical mass is favoured over the bare mass provided the mass anomalous dimension is large enough: γ m > (d e − 2)/2. The non-perturbative criterion for dynamical mass generation is therefore given by γ m (α c ) = (d e − 2)/2 which requires the knowledge of the exact γ m .
In our case, we only have access to a perturbative expansion for γ m . Hence, we would like to find a criterion for dynamical mass generation (which is intrinsically a nonperturbative mechanism) that could be truncated at a given order of the perturbative expansion of γ m . This can be achieved with the help of the gap equation found from the SD formalism. Actually, as we saw in the previous subsection, we know that the all order ansatz (valid at the critical point) for the fermion self-mass is given by Σ S (p) ∼ p −a E (see (4.3)) where the index a is determined self-consistently. Comparing this ansatz to (4.9), and in particular to the second asymptote of this equation, we see that a = d e − 2 − γ m and is therefore related to the mass anomalous dimension. From [64], we learn that the gap equation (which is quadratic in a) is built up from the two asymptotes of (4.9) (at least for large-N QED 3 at NLO in the 1/N -expansion). Extending the result of Gusynin and Pyatkovskiy [64] to arbitrary d e , we therefore assume that, for RQED 4,de , the gap equation is quadratic in a at all loop orders and takes the form: with γ m as the only input. As we saw in the previous subsection, in the SD formalism, the dynamical mass is generated when a becomes complex, i.e., for (a − (d e − 2)/2) 2 < 0. In terms of the mass anomalous dimension and at the critical point, this translates into: from which the critical coupling α c can be computed. Since γ m is gauge invariant, the resulting critical coupling will automatically be gauge invariant too. And as it is built from the SD formalism, it can be truncated to the accuracy at which γ m is known (Eq. (4.10) reduces to (4.8) at the LO in α). From this polynomial equation, we will obtain multiple solutions for α c . In the following, the physical α c will be taken as the the smallest positive real solution that is found. At this point, we would like to comment on the fact that the large mass anomalous dimension required for dynamical mass generation also affects the dimension of the quartic fermion operator ∆[(ψψ) 2 ] = 2d e − 2 − γ (ψψ) 2 where γ (ψψ) 2 is the associated anomalous dimension. Indeed, assuming that γ (ψψ) 2 = 2γ m , we have ∆[(ψψ) 2 ] = 2d e − 2 − 2γ m ≤ d e for γ m ≥ (d e − 2)/2, i.e., the quartic fermion operator is marginal at the critical point and becomes relevant once the dynamical mass is generated. However, from these arguments, the marginality of the quartic fermion operator at criticality appears as a consequence of (4.11) and holds only in an approximate way. According to the literature on QED 4 , see, e.g., the review [51] as well as [79][80][81][82], the assumption γ (ψψ) 2 = 2γ m is supposed to be valid in the quenched and rainbow approximation. 8 In the case of QED 4 , there is evidence that it still holds beyond the rainbow approximation [104]. Note also that the quenched approximation significantly affects QED 4 because in this approximation the coupling does not run. But its effect is weaker for a "standing" theory such as RQED 4,3 . These arguments suggest that the marginality of the quartic fermion operator at criticality may hold even beyond the quenched (especially for RQED 4,3 ) and rainbow approximations. A more careful study of the validity of this approximation goes beyond the scope of this paper.
In the next sections, we will apply (4.11) to the study of the critical properties of RQED 4,de with applications to RQED 4,3 and QED 4 . For the case of RQED 4,3 , we will find perfect agreement with the SD formalism [12] which is natural since (4.11) is built from the SD formalism for large-N QED 3 [64] which can be mapped to RQED 4,3 [12]. As for QED 4 , our approach is more phenomenological since it amounts to extrapolate an equation valid for RQED 4,3 to the more subtle case of a "running" theory.
In the case of QED 4 (ε e = 0), we recover the celebrated result: α c = π/3 = 1.0472, which was obtained via the rainbow approximation in the Landau gauge (which is the "good gauge" for QED 4 at one loop) in the early papers [105][106][107]. In [104], a numerical solution of SD equations with Curtis-Pennington vertex (and hard cut-off regularization) led to α c (ξ) ≈ 0.93 with variations of the critical coupling of only a couple of percent over a wide range of the gauge fixing parameter; such a result deviates by 11% from π/3. In [108], calculations using dimensional regularization, which is a gauge-invariant regularization scheme, were found to agree with the hard cut-off regularization ones to within numerical precision. Moreover, in [109] the result α c = π/3 was extracted analytically from dimensionally regularized SD equations in the Landau gauge in perfect agreement with our own calculation of sec. 4.1. Based on (4.11), we obtained in a very simple way that α c = π/3 is actually the fully gauge-invariant result at LO. This motivates us to include higher order corrections which is quite straightforward in our formalism provided γ m is known.
We therefore consider RQED 4,de at two-loop. Substituting the two-loop expression of γ m , (3.20a), in (4.11), and selecting as a physical critical coupling the smallest value obtained, yields the following result: (4.14) It turns out that α c in (4.13) is actually complex and hence unphysical for all integer values of N except N = 0. This comes from the fact that ∆ N in (4.14) is positive only for values  Table 1: Critical couplings of quenched QED 4 computed from (4.11) up to 5 loops (the symbol "-" indicates that no physical solution is found).

General remarks
In order to be able to consider the case of unquenched RQED 4,de , we shall follow [1,5,12,112,113] and include a dynamical screening of the interaction. In the case of RQED 4,3 , this can be done with the help of the random phase approximation (RPA) that amounts to a simple redefinition of the coupling constant [1,22] because, for this "standing" model, the polarization operator is finite (see (3.8) combined with (2.10)). This procedure resums exactly all the N -dependence in γ m allowing us to go beyond the quenched approximation [12].
Following the RQED 4,3 case [12], we shall proceed with a generalization to RQED 4,de by enforcing the resummation of the N -dependence of γ m at each loop order in an effective coupling constant. As we shall see below, such a procedure does not correspond to RPA in the case of QED 4 . Indeed, applying RPA to QED 4 within our formalism does not allow to fully resum the N -dependence in γ m . This originates from the running of the coupling constant which arises from the UV divergent polarization operator (see (3.9)) and renders the QED 4 case more subtle to treat than the RQED 4,3 one. Enforcing the resummation of the N -dependence of γ m in QED 4 is a semi-phenomenological prescription. But, as we shall see, it will provide results that are in quantitative agreement with those obtained from numerical solutions of SD equations. Moreover, this effective coupling method will allow us to compute a critical fermion flavour number, N c , such that α c → ∞, i.e., a clearly non-perturbative quantity that requires the resummation of the N -dependence of γ m .
Within our formalism, an interesting fact about studying QED 4 is that it allows to consider higher orders in the loop expansion which are presently inaccessible in RQED 4,3 . We already exploited this in the last paragraph related to the quenched case. As we saw there, in both cases of RQED 4,3 and QED 4 , there is no solution to the 2-loop gap equation for N > 0. It turns out that, for QED 4 , solutions of (4.11) appear at higher loops for non-zero N values without any need to introduce an effective coupling constant. We do not find these results satisfactory but, for completeness, we summarize them in Tab. 2 and discuss them briefly.
First, as can be seen from Tab. 2, a parity effect is observed as no solution is obtained at 2-and 4-loop for N > 0. At 3-and 5-loops, we obtain solutions over a range of N -values smaller that some N max which is now larger than 1. However, this N max is not related to N c and it's physical interpretation is not clear; it is close to the N c found from the effective coupling approach at 3-loop but deviates substantially from it at 5-loop, see Tab. 3 for the results of the effective coupling approach. Moreover, at 5-loop the critical coupling obtained this way decreases with increasing N while the opposite behaviour is observed at 3-loop. Actually, on physical grounds we expect that α c should increase with increasing N . This comes from the fact that screening increases with N thus effectively weakening interaction effects which in turn requires a larger value of the coupling constant in order to dynamically generate a mass.
The effective coupling approach seems to overcome these difficulties (as far as our semiphenomenological approach can tell) and we shall therefore focus on it in the following.

Effective coupling constant method
In order to include dynamical screening in RQED 4,de , we define the following effective (reduced) coupling constant:  Table 2: Critical couplings of unquenched QED 4 computed from (4.11) up to 5 loops (the symbol "-" indicates that no physical solution is found).
whereḡ r = g r /(4π) and all the N -dependence, i.e., the effect of fermion loops (or screening), is by definition inΠ (4,de) 1 eff . Substituting (4.17) in the mass anomalous dimension, (3.20a), and expanding the resulting expression up to second order inḡ r , we require thatΠ (4,de) 1 eff cancels the two-loop N -dependent terms in (3.20a). This yields: together with was defined in (3.8), i.e., the procedure exactly corresponds to RPA in agreement with (2.10). In the case of QED 4 (z = 1 and ε e → 0), we find from (4.18) that whereΠ (4) 1 was defined in (3.9). As we anticipated, in the case of QED 4 , the procedure is not RPA (see (2.11)) as it corresponds to a resummation of half of an effective one-loop renormalized polarization operator (at L p = 0).
We may now proceed on solving (4.11) with (4.19) truncated at one-loop order. The effective critical coupling g c that we obtain is equal to the one-loop critical coupling given by (4.12). With the help of (4.17), this result can be generalized to the unquenched case and yields: Note that we can deduce from (4.22) a critical fermion flavour number, N c , such that α c → ∞. Its expression reads: We now apply our general one-loop results to specific cases starting from RQED 4,3 (z = 0 and ε e → 1/2). In this case, (4.22) and (4.23) yield: Note that the one-loop SD equation in the rainbow approximation and including the polarization operator in the LAK-approximation (following the work of Landau, Abrikosov and Khalatnikov [114]) was approximately solved in the Landau gauge in [115] with the results α c (N = 1) = 1.95 (see also [52] for a review of other results as well as discussions below). It deviates by about 69% from our result, α c (N = 1) = 1.1541. Next, in order to access the two-loop case, we solve (4.11) with the full (4.19) as the input. This yields the critical coupling g c = α c (N = 0) where α c (N = 0) is the quenched critical coupling given by (4.16). With the help of (4.17), this result can be generalized to the unquenched case and yields: Proceeding as in the one-loop case, we deduce from (4.28) a critical fermion flavour number, N c , such that α c → ∞. Its expression reads: We now apply our general two-loop results to specific cases starting from RQED 4,3 (z = 0 and ε e → 1/2). In this case, (4.28) and (4.29) yield: thereby recovering in a simple and straightforward way all the results of [12]. As already noticed in [12], we see from (4.30) that the 2-loop N c in RQED 4,3 is slightly higher than the corresponding one in NLO large-N QED 3 (for which N c = 2.8470, see discussion below Thanks to the knowledge of γ m up to 5 loops in QED 4 , we can extend the above method to higher orders. At each order, the effective dynamical coupling is determined in such a way that it resums all the N -dependence of γ m . This allows us to deduce a critical flavour fermion number, N c , and the critical values of the coupling, α c (N ) for N < N c . Our results for unquenched QED 4 (as well as the results of RQED 4,3 for clarity) are summarized in Tab. 3. As in previous approaches, no physical solution is found at 4 loops (all solutions are complex) and therefore we effectively have N c = 0 in this case. At 1-, 2-, 3-and 5-loops, a finite non-zero value of N c is obtained which decreases with increasing loop order. In agreement with our physical intuition, at 1-, 2-, 3-and 5-loops the critical coupling is seen to increase with N (because of the increase of screening) and diverges at N c . Moreover, provided that the apparent decrease of N c with increasing loop order is smooth, we are able to extrapolate its value to infinite loop order using a simple decreasing exponential fit (discarding the 4-loop order). This yields: We therefore conjecture that, in the non-perturbative regime (infinite loop order), the critical coupling is defined only up to N = 2. Moreover, we can also infer that, α c (N = 2) may be extremely large since N = 2 is very close to N c = 2.0663. Unfortunately, the values of the critical coupling oscillate too much thus preventing us from attempting a similar fit. In order to numerically compare our results with the ones of the literature, we should first discuss some "good gauge" issues. Indeed, at a given loop order greater than 1, a  Table 3: Critical couplings of unquenched QED 4,3 and QED 4 , computed from (4.11) together with the effective dynamical coupling (4.17) up to 5 loops (the symbol "-" indicates that no physical solution is found and the symbol "??" indicates that we could not compute the corresponding value).
numerical evaluation of the "good gauge" parameter, ξ g (α c , N ), shows substantial deviations from the Landau gauge as N increases and even a divergence as N approaches N c . At N = 1, a comparison with the Landau gauge results is justified since the "good gauge" is reasonably small, 0. Taking as a reference the review part of [90], table III, a wide range of results have been found for α c (N = 1) in the Landau gauge over the last 30 years. Most of the first order approaches (that neglect one or more equations of the SD system) presented in [90], lead to results that are far from the results of our paper, e.g., a relative difference of 60% up to 101% by comparing with our 5 loop approach. However, better agreements are found when comparing our results with more sophisticated numerical approaches, especially the results of Bloch [52], that solve the full system of one-loop like SD equations with various vertex ansätze (bare, Ball-Chiu and modified CP as referred to in [52]). These various approaches lead to results that are very close to each other in the N = 1 case. Of course, the non-perturbative nature of the vertex ansätze used in these approaches does not have (to the best of our knowledge) a clear correspondence in terms of loops (such as in our case). Nevertheless, from these results, we find deviations of 2%-5% with respect to our two-loop results and deviations of 31%-42% with respect to our 4 and 5-loop results. We also note a smaller deviation of 22% upon comparing our 5-loop N = 1 result with the most recent results of [116].

Conclusion
In this paper we have computed the two-loop mass anomalous dimension of RQED 4,de with N flavours of four-component Dirac fermions. Our main formulas, (3.20), for γ m and also the field anomalous dimension γ ψ , take into account the non-commutativity of the ε γ → 0 and ε e → 0 limits and are valid in both cases of RQED 4,3 (which applies to graphene at its IR Lorentz invariant fixed point) and QED 4 . For RQED 4,3 , our formula for γ m , (3.27), generalizes the one obtained in [40] to an arbitrary number of (2-component) fermion flavours but our interest was mostly focused on the parity-even case (3.26a). When the later is mapped to large-N QED 3 [12], it corresponds exactly to the result obtained by Gracey in [98] thereby strengthening our result.
We then proceeded on studying the critical properties (dynamical generation of a parity-conserving fermion mass) of RQED 4,de (both in the quenched and the unquenched cases) with the help of the gap equation (4.10) (and its solution (4.11)). The latter was derived in a semi-phenomenological way on the basis of the modern approach of Gusynin and Pyatkovskiy [64]. It matches exactly the NLO gap equation of large-N QED 3 [64,66] as well of that of RQED 4,3 [12] and its range of application was extended to cover the case of QED 4 . Its only input being the mass anomalous dimension, it is fully gauge invariant by construction and allows to derive gauge-invariant critical coupling constants and critical fermion flavour numbers of RQED 4,de models with the help of (3.20a). We also underlined the fact that quantities derived with the help of this gauge-invariant method do match (at one-loop from our SD analysis) with gauge-dependant computations evaluated in the "good gauge" (a gauge where the fermion anomalous dimension vanishes order by order in perturbation theory). In the case of RQED 4,3 , we have explicitly checked that (4.11) does allow us to recover in a simple and straightforward way all the NLO results of [12]. In the case of QED 4 , we have first recovered the celebrated α c = π/3 (now exactly gauge-invariant, see (4.12) with ε e = 0) at LO in agreement with the result obtained from solving the SD equations in the "good gauge" (see (4.7) with ε e = 0). We have then used state of the art expressions for γ m up to 5-loops in QED 4 , to access the critical properties of the model. The unquenched case was considered first and our results are summarized in Tab. 1. The latter shows that, at 5-loops, our value for the critical coupling, α c (N = 0) = 1.0941, deviates by only 4% from π/3 and by 15% from the results of [52,90] where numerical solution of SD equations in the Landau gauge with various vertex ansätze lead to α c (N = 0) ≈ 0.93. The unquenched case is more subtle due to the running of the QED 4 coupling constant (in contrast to RQED 4,3 which is a "standing" gauge theory). Nevertheless, a resummation of the N -dependence of γ m into an effective coupling constant gave us access to N c which is such that α c → ∞ together with all values of α c (N ) for N < N c . Our results are summarized in Tab. 3. At 5 loops, we find that N c = 3.3954 with α c (N = 1) = 1.2221. Comparing our results to the ones of the literature, we find deviations of 2%-5% of our two-loop results with respect to the results of [52] and deviations of 31%-42% of our 4 and 5-loop results with respect to the results of [52]. We also note an interestingly smaller deviation of 22% upon comparing our 5-loop N = 1 result with the most recent results of [116]. From our results, we could extrapolate N c to infinite loop order finding N c = 2.0663, see (4.34). This suggests that dynamical mass generation in QED 4 may take place only for N ≤ 2 (and hence for a possibly very large value of α c (N = 2)).
Our work has shown that (4.10) (and its solution (4.11)) allow for a simple and straightforward study of the critical properties of RQED 4,de . In the case of RQED 4,3 , we were limited to two-loop order as, to the best of our knowledge, γ m is still unknown for this model at 3-loop and beyond. It would be very interesting to compute γ m at 3-loop for RQED 4,3 and use it as an input to (4.11). It would also be very instructive to solve the NNLO SD equations for RQED 4,3 along the lines of [66] for the NLO case, in order to have a solid proof that (4.10) does hold at this order. But such analytic computations may be quite tedious at this order and it remains to be seen if they can even be carried out (let's recall that it took approximately 30 years since the seminal work of Nash [61] for the gauge-invariant NLO calculation in large-N QED 3 to be completed in [64,66]). In the case of QED 4 , our straightforward solution of the LO SD equations with the method of Kotikov [62,63] is a good indication that this method might be extended to NLO along the lines of [66] (but now for a running theory in the unquenched case). We leave all these projects for future investigations.

Acknowledgments
We would like to warmly thank Valery Gusynin and Anatoly Kotikov for stimulating discussions and comments on our work as well as Lorenzo Di Pietro and Jingxiang Wu for fruitful discussions on [40]. We would also like to warmly thank Roman Lee for kindly providing us access to the most recent version of LiteRed.

A Master integrals
Following standard notations and conventions, see, e.g., the review [92], the one-loop Euclidean massless propagator-type integral is defined as: Similarly, the two-loop Euclidean massless propagator-type integral is defined as: J(d e , p, α 1 , α 2 , α 3 , α 4 , The dimensionless function G(d e , α 1 , α 2 , α 3 , α 4 , α 5 ) is presently unknown for arbitrary values of the indices α i (i = 1 − 5) but has a long history (see more in [117] and [92]). As can be seen from our exact results below, the only complicated two-loop master integral is G(d e , 1 − ε e , 1, 1 − ε e , 1, 1) with arbitrary indices α 1 = α 3 = 1 − ε e . In the case of QED 4 (ε e = 0), it's expression is known for a long time [118][119][120][121]. For arbitrary ε e , it has been represented as a combination of two 3 F 2 -hypergeometric functions of argument 1 in [25] following the general method of [122] (but the corresponding ε γ -expansion was not provided). The case ε e = 1/2, relevant to RQED 4,3 , has been recently encountered in [123] and it's leading-order epsilon expansion was computed in that paper. From the results of [123], it is straightforward to find that: where C ≈ 0.91597 is the Catalan constant and Cl s (θ) is the Clausen function with Cl 4 (π/2) ≈ 0.98895 . As noticed in [25], the finite diagram (A.4) enters Σ V at two-loop with a factor proportional to ε γ (see the last term in (B.3c) below). Hence, it does not contribute to Σ 2V in the limit ε γ → 0 relevant to RQED 4,de . On the other hand, as can be seen from the last term in (B.2c) below, (A.4) does contribute to Σ S at two-loop order in the limit ε γ → 0. But it only affects it's finite part. Hence, it does not contribute to the two-loop Z m and corresponding anomalous dimension γ m that we focus on in the main text.

B Exact results
In this appendix, we provide exact results for all the diagrams that where considered in the main text. All calculations were done by hand and crossed checked with the help of LiteRed [124,125]. These results will be expressed in terms of the one and two-loop massless propagator-type master integrals presented in A. Except for Σ 1S , all one-loop diagrams were previously computed in [25]. The expression of Σ 1S and, for completeness, the ones of Σ 1V and Π 1 are given by: