Critical properties of three-dimensional many-flavour QEDs

We review several variants of three-dimensional quantum electrodynamics (QED$_3$) with $N_f$ fermion (or boson) flavors including fermionic (or spinorial) QED$_3$, bosonic (or scalar) QED$_3$, $\mathcal{N}=1$ supersymmetric QED and also models of reduced QED (supersymmetric or not). We begin with an introduction to these models and their flow to a stable infra-red fixed point in the large-$N_f$ limit. We then present detailed state-of-the-art computations of the critical exponents of these models within the dimensional regularization (and reduction) scheme(s), at the next-to-leading order in the $1/N_f$ expansion and in an arbitrary covariant gauge. We finally discuss dynamical (matter) mass generation and the current status of our understanding of the phase structure of these models.

A salient feature of both fQED 3 and bQED 3 is that they are super-renormalizable with a dimensionful coupling constant (e 2 has a dimension of mass). Remarkably, early studies [1,2,4] realized that, within a 1/N f expansion, an interacting fixed point emerges in the low-energy limit, and both models become effectively renormalizable deep in the infra-red (IR) with a dimensionless coupling constant, 1/N f . This led to the study of their critical properties with the help of large-N f techniques (see [34] for a review). Critical exponents such as the field and mass anomalous dimensions are particularly important. They encode the renormalization of the composite operator,ψψ [35,36], and play a crucial role in the study of fundamental quantum field theory mechanisms, such as dynamical symmetry breaking and electron mass generation. Precision calculations require the computation of higher-order corrections that often represent a major technical challenge. Beyond precision, these corrections sometimes reveal new physics that is missed by simple low-order estimates. They also allow the study of the stability of a non-trivial IR fixed point with respect to radiative corrections. Finally, while the above models are well behaved in the limit of a large number of flavors, the fate of IR singularities are ubiquitous in super-renormalizable models, and QED 3 often serves as a toy model for such studies [37][38][39][40][41] (see recent progress in [42][43][44]).

Supersymmetric QEDs
Another interesting variant of QED 3 corresponds to (minimal) N = 1 supersymmetric three-dimensional QED (SQED 3 ). This model can be obtained naively by combining the fermionic and bosonic QED 3 models described above, together with a superpartner for the photon, the photino. Mathematically, the degrees of freedom of N = 1 SQED 3 are the 2N f matter multiplets {ϕ j ,ψ j ,F j } and a gauge multiplet {A µ ,λ}. In the matter multiplet, the complex pseudo-scalars, ϕ j , are the superpartners of the two-component Dirac spinors, ψ j , and the F j are complex auxiliary scalar fields without any dynamics (they ensure the equality of the degrees of freedom in the matter and gauge multiplet). In the gauge multiplet, obtained after choosing the Wess-Zumino gauge [45], the photino, λ, is a two-component Majorana field. The action of N = 1 massless SQED 3 is then given by which, similarly to the two previous models, is super-renormalizable and has a non-trivial fixed point deep in the IR, at which it becomes effectively renormalizable [46], the more recent [47] and [33] for a review).
Supersymmetric variants of QED 3 have attracted continuous interest through the last decades. This has been partly motivated by the fact that the enhanced symmetry may simplify the resolution and, perhaps, even lead to an exact solution. As a matter of fact, the case of (non-minimal) N = 2 SQED 3 has been studied in an early seminal paper of Pisarski [3] by dimensional reduction from the case of (minimal) N = 1 four-dimensional supersymmetric QED (SQED 4 ), with focus on dynamical electron mass generation along the lines of the non-supersymmetric case. Actually, in N = 1 SQED 4 , a non-perturbative non-renormalization theorem forbids dynamical mass generation [48], and it was then argued in [46] that it, therefore, extends by dimensional reduction to N = 2 SQED 3 . Further evidence for the absence of dynamical mass generation in N = 2 SQED 3 came from numerical simulations [49] and a refined analytic treatment [50].
The situation in N = 1 SQED 3 is more subtle because of the absence of non-renormalization theorem in this case. The model was first considered by Koopmans and Steringa [46] along the lines set by Appelquist et al. for standard fQED 3 [5]. Their truncated (to leading order (LO) in 1/N f -expansion) Schwinger-Dyson equations approach resulted in a critical fermion flavor number, N c = 1. 62. This implies that a dynamical (parity-invariant) mass generation may occur for N f = 1, i.e., one four-component Dirac spinor. A decade later, additional evidence for the generation of dynamical electron mass in minimal SQED 3 was also found in [51]. There is, however, no rigorous statement for electron mass generation for minimally supersymmetric SQED 3 [46,51].
In the last two decades, N = 1 SQED 3 has attracted significant attention (together with other supersymmetric and non-supersymmetric gauge theories) in the context of the study of IR dualities and renormalization group flows (see, e.g., [52][53][54][55][56][57]). Interestingly, it was argued in [47] that N = 1 SQED 3 at N f = 1 is dual to a conformal field theory in the IR. This suggests that no dynamical mass for the electron should be generated in contrast to the previously mentioned early (leading order) calculations [46,51]. In this review, we will present a refined, next-to-leading order (NLO) analysis. We will show that, at NLO, N c = 0.39, which is strong evidence that no electron mass is radiatively generated in N = 1 SQED 3 , which is in agreement with the analysis based on dualities [47].
At the interface with condensed matter physics, there have also been proposals during the last years that SUSY may emerge in the low-energy limit of various lattice models (see, e.g., [58][59][60][61][62][63][64][65]). To this day, there is still no evidence that SUSY is realized in nature, and an emergent SUSY should certainly be difficult to detect in the lab [66]. Nevertheless, computing critical exponents in the corresponding models is certainly valuable in order to assess the potential impact of supersymmetry on experimentally measurable observables.

Reduced QEDs
Another model corresponds to the so-called reduced QED (QED dγ ,de ) that describes relativistic fermions in d e -dimensional space-time and interacting via the exchange of massless bosons in d γ -dimensions (d e ≤ d γ ). In a Minkowski space, the QED dγ ,de action [67][68][69] reads where ψ i are the 2N f flavors of two-component Dirac spinors in d e dimensions (µ e = 0,1,...,d e − 1). In (4), the coupling of the fermion to the gauge field is restricted to d e -dimensional space-time such that the gauge field is free in the d γ − d e co-dimensional space, i.e., D µe = ∂ µe + ieA µe , where A µe = A µe (z = 0) such that z is the collective coordinate in the (d γ − d e )-dimensional space. It was introduced in [67], motivated by the study of dynamical chiral symmetry breaking in brane-world theories (see also [70]). Soon after, a first application was devoted to the especially important case of conformal QED 4,3 (also known as pseudo-QED from [71] and mixed-dimensional QED from the recent [72]) in relation to graphene [73]. More precisely, QED 4,3 describes graphene [25] at its infra-red (IR) Lorentz-invariant fixed point [74]. Importantly, it has been shown in [75] that there is mapping between QED 4,3 and fermionic QED 3 in a large-N f limit.
Theoretically, there has been rather extensive studies on QED 4,3 during the last decade with primary applications to graphene-like systems, e.g., their transport and spectral properties [68,69,[76][77][78][79][80] and quantum Hall effect [72,81] (in [72], the model was invoked as an effective field theory describing half-filled fractional quantum Hall systems) and dynamical symmetry breaking [75,82], on which we will focus in the following. From a more field theoretic point of view, the model was shown to be unitary [83], and its properties were studied under Landau-Khalatnikov-Fradkin transformation [84,85] as well as under duality transformations [86]. Renewed interest in the model and its formal properties was triggered by a study [87] on interacting boundary conformal field theories (see, e.g., [88][89][90][91][92]).
Motivated by condensed matter applications, supersymmetric extensions of reduced QED were also constructed and analyzed in [93,94] via superconformal techniques on the boundary for both N = 1 and N = 2 cases. In [94], non-perturbative computations of transport properties in the N = 2 case were carried out with the help of localization techniques. Once again, the N = 1 case is more subtle, and no exact solution is, so far, known. It is this case that will be of interest to us in the following. In particular, the supersymmetric extension of QED 4,3 will be denoted as SQED 4,3 . In this ultra-relativistic super-graphene model, the matter fields (electrons and selectrons) are localized on a (2 + 1)-dimensional plane, while the gauge fields (photons and photinos) are (3 + 1)-dimensional. The corresponding action reads where d e = 3 and d γ = 4, Λ denotes a four-component Majorana field, Γ µγ are four 4 × 4 gamma matrices, γ µe are three 2 × 2 gamma matrices and D is a real auxiliary field without dynamics (see [93] for more details). As will be seen in the following, just as for the non-SUSY case, there exists mapping between SQED 3 in the large-N f expansion and SQED 4,3 in the small-coupling expansion. The mapping is very similar to the non-SUSY case up to a factor of two from SUSY.
Let us emphasize that, in this review, we always consider suspended (super)-graphene, as opposed to a model defined on the boundary, as considered in, e.g., [93]. In our case, the boundary is considered as a transparent interface, while the model of [93] considers a purely reflecting boundary (graphene on a substrate). Nevertheless, the two models can be related by doubling the interaction, α bdry = α/2.

Outline of the Review
In order to adopt a unifying approach, we will introduce in Section 2 a general model that encompasses all the above-described models. We will present the perturbative setup, consisting of the Feynman rules (including subtleties related to the presence of Majorana spinors) and our renormalization conventions. In Section 3, we will then perform the full LO and NLO computations in the large-N f expansion of all the polarizations and self-energies and derive all the corresponding anomalous dimensions. This will lead us to briefly discuss a generalized version of the Furry theorem for these models. In Section 4, we apply our general results to the models of interest, i.e., fQED 3 , bQED 3 , SQED 3 and, finally, QED 4,3 and SQED 4,3 , where we compute the optical conductivity of the (super-)graphene material. In Section 5, we will discuss the criteria for dynamical mass generation and the related phase structure of these models. The conclusion will be given in Section 6. Let us finally note that this review represents an extended version of our previous (short) paper [95].

The General gQED 3 Model
In the Introduction, we formally introduced a total of five models, namely, fQED 3 , bQED 3 , SQED 3 and also the reduced QED 4,3 and SQED 4,3 . In this section, we introduce a general model that encompasses the first three of them with the help of additional parameters. The reduced models will be analyzed from the results of the general model with mapping [75].
The general model is denoted as gQED 3 , and its corresponding action reads In (6), the first line corresponds to the fQED 3 action (1). The second line is the N = 1 SUSY content, where each superpartner field is associated with a tracking factor S ∈ {0,1} such that Φ → SΦ, ∀Φ ∈ {ϕ,λ,Ā µ } and S 2 = S. Hence, at any step of the calculation, we may turn on (respectively, off) SUSY by setting S = 1 (respectively, S = 0). This will highlight SUSY effects in our computations and allow us to check our expressions by recovering known results for the corresponding non-SUSY theories. The third line of (6) is due to the use of the dimensional reduction (DRED) scheme [96][97][98] (see also [99] for a review), which is the most convenient regularization scheme for practical calculations in supersymmetric theories. DRED allows for preserving SUSY at the perturbative level by the introduction of extra particles, the so-called ε-scalars, carried by the fieldĀ. These particles arise from the formal splitting of the gauge field as A µ =Â µ +Ā µ . Here, we use the notations of the review [100], where hatted (respectively, barred) quantities have d (respectively, 3 − d) components. To better appreciate the effects of DRED during the computations, the ε-scalar field will be associated with a tracking factor E ∈ {0,1} such thatĀ µ → EĀ µ and E 2 = E. Indeed, as we shall see in the following, although ε-scalars affect only few quantities at NLO, their effect is crucial in order to ensure the validity of supersymmetric identities.
In addition, we will work with 2N f arbitrary n-component spinors. In the SQED 3 case, n = 2 is necessary to ensure the equality of the matter and gauge degrees of freedom. Nevertheless, working with arbitrary n component spinors will allow us to take the case of n = 0-component spinors, i.e., no fermions, which corresponds to the case of bQED 3 . Indeed, by killing the spinorial degrees of freedom with n = 0, one exactly recovers the action of bQED 3 , (2). In order to keep track of both cases while limiting the complexity of our formulas, one can notice that the identity n 2 S = 2nS holds in both cases. We shall, therefore, use the constraint, n(n − 2)S = 0, to simplify our computations.
The action (6) is the general model we will work with in the rest of this article. It completely describes N = 1 supersymmetric QED in the DRED scheme with suitable parameters (S, E, n) that allow the recovery of the subcases of fQED 3 and bQED 3 (as well as QED 4,3 and SQED 4,3 via mapping) as well as the study of the effect of DRED by turning it on (or off) with E = 1 (or 0). These parametrizations are summarized in Table 1.
Let us remark that the action (6) is completely massless. In the following, in order to compute the mass anomalous dimensions of the model, we will introduce a mass term for the matter multiplet, i.e., for the electron and the selectron. Since we are interested in dynamical mass generation, we will focus on the parity-even mass terms, i.e., of the form, Moreover, we will work within the limit of small masses, i.e., with p E the Euclidean momentum. This limit will have, as a main advantage, to remove all the tadpoles in the theory. In bQED 3 and SQED 3 (and, therefore, SQED 4,3 ), the four-point bosonic coupling indeed gives rise to tadpoles, as opposed to the case of fQED 3 , where no tadpoles are present, since the theory only has three-point coupling. The masses will, therefore, enter the electron and selectron propagators as small IR regulative masses, allowing the computation of their corresponding mass anomalous dimensions. All our calculations will then be carried using massless techniques (see, e.g., the review in [101]). We postpone to Section 5 the study of a potential dynamical generation of small parity-even masses (7) in the matter multiplet.

Feynman Rules
The gQED 3 model (6) contains both Dirac and Majorana fermions. Therefore, one has to be careful to properly define the Feynman rules of the model in order to avoid sign mistakes. In the following, we will use a method based on the conventions of [102,103]. We first derive the bare gauge-multiplet propagators from the general action (6), reading It is important to remark that the photino (λ) Majorana line (9c) carries a fermion flow, but is not represented with a dedicated arrow. We also derive from (6) the bare matter-multiplet propagators, reading (10b) Note that the arrow on the Dirac fermion (ψ) and the pseudo-scalar (ϕ) propagators indicates the charge flow or, equivalently, the matter flow. As for the photino (9c), the fermion flow on the Dirac fermion line (10a) is not indicated. Together with these gauge and matter propagators comes additional rules for the loops • Each matter-field loop (ψ and ϕ field charge flow) gives a factor of 2N f , i.e., graphically • Each fermion loop (ψ and λ field fermion flow) gives a factor (−1) and a trace over the spinorial indices, i.e., graphically Lastly, we provide all the vertices of the action (6), yielding in graphical form Note that the first vertex (13a) is purely of fQED origin, the second line (13b) are the vertices of bQED origin, then the vertices (13c) come from the ε-scalar contributions and, finally, the vertices (13d) are of pure SUSY origin.
In addition to all these rules, one should be careful about fermion flows when both Dirac and Majorana fermions are present. This usually results in a multitude of additional Feynman rules to cope with all the possible flow cases in order to obtain the correct signs everywhere. In the following, we will use the compact Feynman rules of [102,103] that are based on assigning an additional fermion flow line on diagrams (when necessary) along fermionic lines to obtain the correct signs. The additional Feynman rules are then written down by specifying the fermion flow (arrow above) and, for the fermionic propagators (recalling that the middle arrow is the charge/matter flow and the bottom arrow is the momentum), they read which amounts to adding a minus sign on the flowing momentum for each opposite arrow. Similarly, for the fermionic Dirac vertices (fermion flow indicated with the arrow on the right), they read which amounts to a complex conjugation (charge conjugation) of the vertex if the fermion flow goes backward with respect to the charge/matter flow. Note that the other vertices mixing both Majorana and Dirac fermion (see (13d)) are real and are, therefore, unchanged under the inversion of the fermion flow.
Actually, in the vast majority of cases, the simple rules (9) to (13), i.e., without the additional fermion flow lines (14) and (15), are sufficient. This comes from the fact that most of the diagrams that we consider are such that the charge flow can follow naturally the fermion flow, both being continuous and unidirectional, i.e., graphically where the hidden fermion flow goes from left to right, i.e., through the Dirac fermion, then the Majorana fermion and then the Dirac fermion again so that all arrows are properly aligned (reversing any of the arrows in this diagram would generate non-trivial minus signs not accounted for in the simple Feynman rules above). In such a case, provided that the momentum arrows follow the (hidden) fermion flow and the charge flow, one can safely use the simple Feynman rules without the sign corrections shown above, i.e., (9) to (13).
Nevertheless, the advanced Feynman rules (14) and (15) will be required for a few diagrams, where one encounters a configuration of the type In such a case, we are forced to use the advanced Feynman rules (14) and (15). In the following, this will be the case for only one diagram, which is the seventh (labeled (g)) diagram of the two-loop contribution to the fermion self-energy at NLO, i.e., Σ ψ(g) 2 (see (104e)).
Note that, in principle, these advanced Feynman rules are also needed for the computation of the photino polarization because of the Majorana external legs. However, these diagrams are always appearing in pairs (with respect to opposite charge flows) that are exactly equal, such that we can consider only the case where all arrows are aligned and double the result. See the discussion below (56) for an example.
We conclude this section by a brief warning to the reader that would like to use the software Qgraf [104,105] (as we did) to generate the diagram expressions of any theory involving both Dirac and Majorana fermions. First, Qgraf does not seem to be able to provide the correct minus signs from the fermionic loops in (12). The simplest solution we found is to include additional trivial delta functions, δ αβ , in the propagators for ψ and λ, where α,β are the spinor indices, such that δ αα = −1.
Similarly, one can implement in an automated way the inclusion of the factors 2N f for (12) with similar delta functions on the fields ψ i and ϕ j , i.e., δ ij such that δ ii = 2N f . Moreover, Qgraf may have trouble in generating diagram expressions with continuous and unidirectional fermion flows in rare cases. More specifically, the software seems to always generate the flow properly (i.e., the indices generated by Qgraf that we use to orient the charge and fermion flows are aligned with the momenta arrows), except when there is an isolated fermion between two Majorana or the reverse, i.e., a chain of the form (17). In this particular case, we need additional routines to check the Qgraf output and possibly correct these fermionic flows by using the rules (14) and (15). As advertised before, our routine has corrected only one diagram in the NLO computations, which is Σ ψ(g) 2 (see (104e)).

Numerator Algebra
We work in a three-dimensional Minkowski space with the metric g µν = diag(+,−,−). The three n × n Dirac γ-matrices satisfy the usual Clifford algebra, {γ µ ,γ ν } = 2g µν I n , where Tr(I n ) = n. Since we work in the DRED scheme, the metric tensor and γ-matrices are decomposed as so that there are d = 3 − 2ε matricesγ µ and 2ε matricesγ µ , in order to keep a total integer number of three matrices γ µ . All of these matrices are of arbitrary size n × n to be able to take the limits n = 0 for bQED 3 , as well as n = 2 for SQED 3 and fQED 3 . In the DRED scheme, the following intuitive properties hold as well as the very important case of the mixed dimensional anticommutator As expected, the usual Dirac trace computations will be modified but in a somewhat trivial way thanks to the property (20). In the following, we will have to compute traces involving gamma matrices living in two different spacetimes, such as Tr(γ µγνγργσ ). This requires some care. In practice, one first sorts out the matrices, e.g., gathers hatted matrices to the left and barred ones to the right. This can be conducted using repetitively the anticommutation of the hatted and barred matrices (20), i.e., Tr(γ µ ···γ νγρ ···γ σ ) = −Tr(γ µ ···γ ργν ···γ σ ).
Once completely sorted, one splits the traces into two parts using the following crucial trace splitting formula where all matrices on the left are hatted and all matrices on the right are barred. Once sorted and split, both traces can be computed using the usual algorithm and the same algorithm for traces over only barred matrices,γ µ i . These recursive formulas allow us to reduce any trace methodically until reaching the fundamental ones Tr(I n ) = n, Tr(γ µ ) = 0, Tr(γ µ γ ν ) = ng µν , Tr(γ µ γ ν γ ρ ) = inT n ε µνρ .
At this point, we recall that in three-dimensional theories, the trace over three gammas may not be zero, depending on the choice of the representation for the γ matrices, but proportional to the fully antisymmetric tensor, ε µνρ . To this end, we introduce the additional parameter T n , such that T 2 = 1, T 4 = 0 and T 2 n = T n . Anyway, in large-N f massless three-dimensional QED 3 (fermionic, bosonic and supersymmetric), these odd traces do not contribute to any result, as expected from parity-even theory. We will explicitly check this fact by observing that T n will never appear in the rest of this article, even if we perform all the computations taking it into account. In the DRED scheme, the trace identities (24) split into two copies with the following intuitive properties Tr(γ µ ) = 0, Tr(γ µγν ) = nḡ µν , Tr(γ µγνγρ ) = 0.
Note that we takeε µνρ = 0, as it makes sense that the Levi-Civita tensor in 2ε dimensions vanishes as ε → 0. Using the (mixed dimensional) trace techniques described above allows for computing any fermionic trace in gQED 3 and its subcases.

Renormalization Setup
We now have sufficient background material to introduce the renormalization setup and conventions for the gQED 3 model. Upon turning on the interactions, the Feynman rules for the gauge multiplet (9) are dressed via their respective Dyson equations and read where the polarizations, Π x , for the photon (Π γ ), the ε-scalar (Π ε ) and the photino (Π λ ), are parameterized via the following projectionŝ respectively. Using this setup, all integrals can be carried out in the massless limit, i.e., m x → 0 for x = {ψ,ϕ}, as an IR rearrangement.
An important remark is that in (26a), the tensorial structure still yieldsd µν (p) =ĝ µν − (1 − ξ)(p µpν /p 2 ) because we are using a non-local gauge, i.e., we take where ξ will still be considered as the gauge-fixing parameter in the following. This trick is widely used in the QED 3 literature to keep computations light (see, e.g., [35]; see also, for the SUSY case, [106][107][108]). We recall that the use of a non-local gauge (28) is possible without affecting the physical results because the gauge-fixing parameter, ξ, is a mathematical artifact that does not appear in physical results.
As we will prove explicitly in the next sections, all the polarizations (27) are finite. Indeed, we recall that, in the large-N f limit, SQED 3 [46], similarly to bQED 3 [1, 2] and fQED 3 [3,4], is a nonrunning ("standing") gauge theory, i.e., the coupling is not renormalized, implying finite polarizations and, therefore, vanishing beta functions. This leads to the triviality of the renormalization constants for the coupling, gauge-multiplet fields and gauge-fixing parameter, formally Z x = 1, and γ x = 0 with x ∈ {e,γ,ε,λ,ξ}, which imply a trivial beta function for the running of the coupling e reading β = −2εᾱ, where α = e 2 /(4π) andᾱ = α/(4π). In this case, the coupling trivially renormalizes as α → µ 2ε α, where µ is the renormalization scale. In the following, we will work in the modified minimal subtraction scheme, where the renormalization scale is defined asμ 2 = 4πe −γ E µ 2 , and further, (MS scheme) subtracts 4π and γ E , the Euler-Mascheroni constant. We will refer to this modified version of the dimensional reduction scheme as DRED. Now considering the matter multiplet, turning on the interactions leads to the following dressed propagators where the matter-multiplet self-energies are parameterized as From these, the components p and m can be extracted with the following projectors where m x = {m ϕ ,m ψ }. As for the gauge polarizations, using this setup allows all integrals to be computed in the m x → 0 limit, i.e., completely massless, as an IR rearrangement.
The renormalization conventions for the non-trivial renormalization constants are defined as The renormalization constants can be extracted from the bare self-energies thanks to the expression of the renormalized self-energies leading to the following simple set of relations where "finite" means of the order of ε 0 , so that no additional counter diagrams needs to be computed. Finally, the associated anomalous dimensions are defined as and correspond to the critical exponents of the theory that we want to compute.

The Large-N f Expansion
In this section, we briefly introduce graphically the idea of the large-N f expansion (see, e.g., [34] for complete a review). Let us consider fQED 3 for simplicity. We first recall that, in the loop expansion, the Dyson equation (26a) for the photon can be written graphically as In the loop expansion (36), the perturbative series is well defined in the small-coupling e regime, just by vertex counting. When the coupling, e, is not suitable as the expansion parameter, such as in superrenormalizable theories like gQED 3 , one can use the so-called large-N f expansion technique. Naively, the series (36) is not perturbative in this regime since each fermion loop gives a factor N f , thereby increasing with the complexity of the diagram. The trick to perform the 1/N f expansion is then to resum the infinite chain of simple matter loops in force field propagators. Hence, considering the first term of each line in (36), i.e., the simple bubble chains, we can define the new propagator Going further, we recall that the bare photon propagator has momentum dependence, ∼ p −2 , and the fermionic simple bubble diagram reads ∼ e 2 N f p E (with p E the Euclidean momentum). Therefore, the new propagator (37), in the large-N f limit, reads This new photon propagator is then said to be softened [1,2] since its behavior in the infra-red is attenuated. Using this softened propagator, the first contribution to the electron self-energy is, therefore, where the (dimensionful) coupling, e 2 , drops in favor of 1/N f . Therefore, in the large-N f limit (that takes into account an infinite number of diagrams), fQED 3 becomes renormalizable with dimensionless coupling 1/N f . Moreover, the expansion for the dressed photon propagator (36) can be rewritten as Similarly, at the next-to-leading order (NLO), one can resum the two-loop contributions, yielding a new propagator, 2 ∼ 1/N 2 f , in the IR limit, which allows computing NLO corrections to the electron self-energy at 1/N 2 f , etc. So, the strategy goes as follows. At leading order (LO): (1) Compute the one-loop polarization using bare Feynman rules and compute the LO-softened photon by resumming the one-loop polarization.
(2) Compute the other diagrams of the theory at O(1/N f ) using the LO-softened photon only.
Then, at next-to-leading order (NLO): (3) Compute the two-loop polarizations using the LO-softened photon propagator and compute the new NLO-softened photon propagator by resumming the two-loop polarization.
(4) Compute the other diagrams of the theory at NLO, i.e., O(1/N 2 f ) using both the LO and NLOsoftened photon propagators.
and pursue similarly at NNLO if desired, which goes beyond the scope of this review.
This reasoning can be easily extended for the full gQED 3 model by resumming all polarizations of the gauge multiplet. In general, large-N f techniques are expected to be very powerful as they resum an infinite number of diagrams. Moreover, since the new coupling of the theory is 1/N f , the value of α = e 2 /(4π) can be arbitrarily large, which is extremely useful to study the critical properties of the corresponding field theories that originate from non-perturbative effects.

Gauge-Multiplet Polarizations at LO
In this first section, we compute in detail the first correction to the polarizations of the gauge multiplet, i.e., for the photon, the ε-scalar and the photino, at LO in the 1/N f expansion, i.e., at O(N f ).

Photon Polarization at LO
We first consider the photon propagator (26a) and compute the LO photon polarization operator, which consists of the following two contributionŝ Graphically, the corresponding two diagrams read Note that the first diagram (a) is of pure fermionic (QED 3 ) origin, while the second one (b) is of pure bosonic (bQED 3 ) origin. Therefore, at this order, the SQED 3 photon polarization directly appears as a simple sum of the fermionic (spinorial) and bosonic (scalar) results. Using the Feynman rules for the vertices (13) and the matter (electrons and selectrons) propagators (10) leads to the following expression Using the photonic polarization projector (27a) and performing the trace on the d = 3 − 2ε (hatted) space using the recursive Formula (23) gives the following expressions These integrals, once wick rotated to the Euclidean space, are then straightforward to compute using the results of Appendix A and yield in the DRED scheme where G(d,α,β) is known exactly and defined in Appendix A. Performing the ε-expansion yields where L p = log(−p 2 /µ 2 ). As expected, in the fQED 3 case (S = 0, n = 2), only the first diagram, which is purely fermionic, contributes. In contrast, in the bQED 3 case (S = 1, n = 0), only the second diagram, which is purely bosonic, contributes. The total photon polarization function is, therefore, given by and since it is finite in d = 3, its exact expression in this dimensionality reads Interestingly, in the cases of SQED 3 (S = 1, n = 2), fQED 3 (S = 0, n = 2) and bQED 3 (S = 1, n = 0), we deduce the following results In this very simple case, the SQED 3 photon polarization is simply the sum of the fermionic and bosonic parts since there is no one-loop diagram involving a mixture of both. Therefore, the SQED 3 photon polarization is twice the value found for fQED 3 , which was first obtained in [1,2]. Note that our result for SQED 3 coincides with the earlier one-loop result given in ref. [46] but now obtained in the dimensional reduction scheme.

ε-Scalar Polarization at LO
Next, we proceed similarly for the ε-scalar propagator (26b) and compute the LO ε-scalar polarization function, which consists of a single non-vanishing diagram, defined as Using the Feynman rules for the vertices (13) and the matter (electron and selectron) propagators (10) leads to the following expression Using the projector defined in (27b) and performing the trace in the 2ε-dimensional (barred) space with the help of the recursive formula (23) yields After wick rotation and using the results of Appendix A, we have Since this result is again finite in d = 3, it can be written as Let us note that in the case of bQED 3 (S = 1, n = 0, E = 0), as well as in the case of fQED 3 (S = 0, n = 2, E = 0), this polarization is obviously zero. Indeed, the ε-scalars are relevant only in the case of which is exactly equal to the polarization of the photon in the same case, Π SQED 3 1γ (p 2 ) calculated in (49). As we will comment later on, such an equality is expected from SUSY.

Photino Polarization at LO
Lastly, we proceed in the same way for the photino propagator (26c) and compute LO photino selfenergy, which consists of two non-vanishing diagrams with opposite charge flows. Since it is a photino (Majorana) polarization, in principle, we need to follow the advanced Feynman rules (14) and (15), where we assigned a continuous and unidirectional fermion flow that goes from left to right. In the case of the diagram (a), all flows are in the same direction so that there is no need for the advanced Feynman rules. In the case of diagram (b), the charge flow and the fermion flow are opposite to each others. However, the charge flow is also opposite to the momentum flow so that the momentum obtains an additional minus sign. All together, the Dirac propagator remains unchanged, and diagram (b) is equal to diagram (a). The resulting contribution is then defined as twice the configuration where all flows are aligned It turns out that this reasoning will apply to all the photino polarization diagrams at higher orders. Therefore, for automation purposes, one can always consider only the configuration where all flows are aligned and simply multiply it by two so that the advanced Feynman rules with fermion flow specification are almost never needed (see discussion below (17)). Using now the simple Feynman rules (11) and (13), its expression reads Then, using the projector (27c) and performing the fermionic trace, we have and, after wick rotation, using the results of Appendix A, it yields Since this result is again finite, we set it exactly in d = 3 and obtain Note that this result is relevant only in the SQED 3 case (S = 1, n = 2), reading which is exactly equal to both the one-loop photon (49) and ε-scalar (54) polarization functions. Summarizing, we find that for SQED 3 , the photon, ε-scalar and photino self-energies are all equal and finite at the LO of the 1/N f -expansion, reading which is first-order perturbative proof that the polarizations are all equal in the gauge multiplet, as expected from SUSY. Moreover, the finiteness of the polarizations is a first-order perturbative proof that the theory has no anomalous dimensions for the gauge fields in accordance with the fact that it is a standing gauge theory, as previously advertised.

IR-Softened Gauge Multiplet at LO
We are now in a position to compute the softened gauge propagators at the leading order of the 1/N f expansion (see (37)). By substituting the one-loop (LO) results obtained for the polarization of the photon (48), the ε-scalar (53) and the photino (61) into the definition of the dressed gauge propagators (26), the propagators soften in the large-N f limit, i.e., p E ≪ N f e 2 , and read where the tensorial structure of the photon is still given byd µν (p) =ĝ µν − (1 − ξ)(p µpν /p 2 ) thanks to the use of the non-local gauge (see (28)). These new softened propagators can then be used to compute the LO self-energies of both the electron and its superpartner.

Matter-Multiplet Self-Energies at LO
In this section, we compute in detail the first correction to the self-energies of the matter multiplet, i.e., for the electron and the selectron, at the LO in the 1/N f expansion, i.e., at O(1/N f ).

Electron Self-Energy at LO
We start with the electron propagator (30) and compute its LO correction, which consists of three contributions one for each gauge interaction, which are defined as where the photon, ε-scalar and photino propagators are indeed the IR-softened ones at first order (64). Using the simple Feynman rules here is enough, i.e., using Equations (9)-(13), and we obtain Note that in (67), the (dimensionful) electric constant, e, drops out in favor of the new coupling 1/N f thanks to the softening of the gauge-multiplet propagators. These diagrams are then split into the part proportional to the external momentum, p, (also called the vectorial part sinceit is proportional to / p) and the one proportional to the mass, m ψ , (also called the scalar part) using the projectors (31). First, focusing on the vectorial part, using the projector (31a) and computing these three diagrams with projection, trace calculation, wick rotation, integral evaluation and wick rotate back, we find the following exact results is finite due to the ε-scalar, while the two other contributions are singular in the limit d → 3. Secondly, focusing on the scalar part, using the projector (31b) and computing these three diagrams with the same approach yields the following exact results where the first contribution is singular in the limit d → 3, while the second diagram vanishes in d = 3, and the last graph (c) is exactly zero because of the gamma matrix trace. Summing all the contributions, the total vectorial and scalar electron self-energies are therefore, given, expanded in d = 3 − 2ε, by Note that some log(2) are resummed by adding a 4 next to the momentum p 2 . From this result, we extract straightforwardly, with (34), the LO electron wave function and mass renormalization As expected, the general mass renormalization factor is completely gauge-invariant, which is a strong check on our results. From the definition of the anomalous dimension (35), we derive the general anomalous dimensions In the relevant cases of SQED 3 (S = 1, n = 2) and fQED 3 (S = 0, n = 2), the general results simplify as Note that γ SQED 3 ψ vanishes in the Landau gauge (ξ = 0), which is then the so-called "good gauge" at LO. This is to be contrasted with the non-supersymmetric case, γ fQED 3 ψ , (first obtained in [35]) that vanishes in the so-called Nash gauge [6], ξ = 2/3. Note also that the bQED 3 case is obviously irrelevant here since we consider the anomalous dimension of the electron field and mass. We will further discuss the quantities (73) once we obtain their supersymmetric counterpart in the next section and their NLO correction after that.

Selectron Self-Energy at LO
We proceed similarly for the electron superpartner, the selectron, which is the scalar propagator and compute its LO scalar self-energy, which consists of the sum of the two diagrams that are defined as where the photon and photino propagators are indeed the IR-softened ones (64) and Σ ϕ 1b contains a hybrid (Dirac/Majorana) fermion loop. Note that, for the diagram (b), we can assign a counterclockwise fermion flow and momentum flows that follows the fermion loop consisting of the Dirac and Majorana fermions. Therefore, using the simple Feynman rules given by (9) to (13) is indeed enough and leads to −iΣ ϕ(a) 1 Performing the traces and using the projection (31b) yields the two LO contributions to the momentum part of the selectron self-energy as well as the two LO contributions to the mass part of the selectron self-energy The total selectron momentum and mass self-energy, in ε-expanded form, yields From these results, we extract the LO scalar wave function and mass renormalization As expected from such a physical quantity, the general mass renormalization factor is completely gauge-invariant. We can now derive the anomalous dimensions for the selectron field and mass using the definition (35), yielding Note that in the relevant cases of SQED 3 (S = 1, n = 2) and bQED 3 (S = 1, n = 0), the above general results simplify as A few remarks are necessary here. First, we observe that for SQED 3 , the mass anomalous dimension for the selectron (82a) is identical to the one of the electron (73a), i.e., as expected from supersymmetry. In striking contrast, the field anomalous dimensions for the selectron (82a) and the electron (73a) are different. This is due to the use of a gauge-fixing term that breaks supersymmetry (Wess-Zumino gauge). This is not an issue since the breaking of SUSY will occur only for gauge-dependent quantities that are, by definition, non-physical. Secondly, let us remark that in the SQED 3 case, the field anomalous dimension of the selectron vanishes for ξ = 2. Since for the fermionic part of SQED 3 , the good gauge was the Landau gauge ξ = 0, and it is, therefore, not possible to cancel both matter-field anomalous dimensions at the same time. This may cause trouble for computations of the critical properties of the model using the Schwinger-Dyson equations (see the devoted Section 5). As for the bQED case (82b), we see that the good gauge is then ξ = 8/3. We will further discuss these results again after improving them to NLO.

Vanishing Contributions and Generalized Furry Theorem
Before going to higher orders and computing any NLO diagrams, we first need to discuss some additional diagrams that may enter the incoming NLO computations as subdiagrams. These LO diagrams are made of matter bubbles and triangles and are of uttermost interest because a lot of them are vanishing, either exactly or in pairs. On the one hand, this will tremendously reduce the number of diagrams to be computed at NLO but will also ensure that matter bubbles are connected to each other in a way suitable for the large-N f expansion.
We first focus on three exactly-vanishing one-loop diagrams made of a matter bubble, as shown in Figure 1. The first vanishing bubble contribution is the mixed polarization B 1 (ongoing photon and outgoing ε-scalar), as displayed in Figure 1a. Since this diagram is proportional to Tr(γ µ ) = 0, it reads B 1 = 0.
More generally, we conjecture that every diagram with an odd number of ε-scalar external lines is exactly zero. The two other contributions, given by Figure 1b,c, are also exactly vanishing by parity on the internal momentum integral, i.e., B 2 = B 3 = 0.
Multiple other vanishing contributions come from matter triangles. These are built from triangles of electrons and selectrons together with external legs of any allowed kinds, i.e., taken in the gauge multiplet. In total, there are 8 triangles (disregarding the possible charge flows), as shown in Figure  2. In the following, we will briefly describe why all of these diagrams (that are proportional to N f ) are vanishing.
The first one is the pure fermionic QED matter triangle diagram, T 1 , and is vanishing because it always appears paired up with its mirror conjugate diagram with opposite charge/matter flow. An explicit computation is a check at the first order of the Furry theorem in QED [109], which is the all-order proof that in QED, and any diagram with an odd number of photon legs can be discarded since they will cancel with their opposite flow diagrams, as a direct consequence of the conservation of energy and charge conjugation symmetry. In the following, we will discuss how the Furry theorem, generalized for gQED 3 , holds at least at the leading order, i.e., for diagrams made of matter triangles and three external legs taken in the gauge multiplet. The second vanishing contribution is the pure bQED diagram, T 2 , and also vanishes with its opposite charge flow counterpart. It, therefore, generalizes Furry's theorem to the case of bQED 3 . The third and fourth vanishing diagrams are the supersymmetric triangles made of mixed electron and selectrons, thereby generalizing Furry's theorem to SQED 3 without ε-scalars. Finally, as for the diagrams containing ε-scalars, we have the three contributions T 5 , T 6 and T 7 that are exactly zero for any momentum or charge flow direction. This is because they contain an odd number of ε-scalar external legs, i.e., they are ultimately related to Tr(γ µ ) = 0. We are left with a last triangle, T 8 , made of an electron loop together with one photon plus two ε-scalars external legs. This diagram is different since it is not vanishing because of ε-scalars, as in this case it is proportional to Tr(γ µγν ) ̸ = 0. We then have to consider the diagrams with the two opposite charge flows, that also vanishes upon explicit computation. We have to then check explicitly that every matter triangle does indeed vanish, either exactly or with respect to their (reversed matter flow) twin diagrams. This completes the perturbative leading order proof that the generalized Furry theorem holds in SQED 3 withing the DRED scheme, and, therefore, as subcases, also in fQED 3 and bQED 3 . This means that every diagram containing a matter triangle can be set to zero in gQED 3 , i.e., in SQED 3 , fQED 3 and bQED 3 . Some prominent examples of diagrams that we can drop are the large number of Aslamazov-Larkin-type diagrams. Taking into account these various vanishing contributions tremendously reduces the number of diagrams that has to be computed at NLO. Indeed, as we shall see in the following, it ensures that up to NLO, not a single diagram of three-loop type needs to be computed. This is crucial because the three-loop master integrals with half integer indices are still unknown and are a big challenge to compute due to the inherent branch-cut structure of the integrals, which results in intricate hypergeometric functions and transcendental numbers (Catalan number, Clausen function, etc.; see, e.g., [110,111]). Moreover, the generalized Furry theorem at LO also guarantees that matter loops are connected by simple chains of force field propagators, like in the simpler fQED 3 case, in accordance with our starting assumption, ensuring that the large-N f expansion is reliable. We can now go forward and proceed with the NLO computations.

Gauge-Multiplet Polarizations at NLO
In this section, we compute the NLO polarizations of the gauge multiplet, i.e., for the photon, the ε-scalar and the photino at NLO in the 1/N f expansion, i.e., at O(1/N 0 f ). We will show that all of these polarizations are finite and gauge-invariant for gQED 3 . In the following, we shall use shorthand notation for the polarization results

Photon Polarization at NLO
We first consider the NLO correction to the photon polarization that consists of 20 Feynman diagrams labeled (a,b,...,t). Taking into account the fact that mirror conjugate graphs take the same value, we are left with 11 distinct graphs to evaluate. This can be conducted exactly for all the diagrams, following the same procedure as for the one-loop case. Their expressions read Summing all the contributions (85), all poles cancel, and the final result is finite, reading Several remarks are in order here. First, the ε-scalars do not contribute here (the corresponding tracking factor, E, is absent) because the two contributions,Π γ 2lm andΠ γ 2n , cancel each other. Second, the result is completely gauge-invariant, which provides a strong check on our result. Lastly, the finiteness of the results ensures that the theory is still standing at NLO, i.e., the coupling does not renormalize.
In the different cases of interest, i.e., SQED 3 (S = 1, n = 2), fQED 3 (S = 0, n = 2) and bQED 3 (S = 1, n = 0), it yields the correction coefficients As advertised in the previous Section 3. 3, it turns out that all the diagrams considered are twoloop. Indeed, since we are in the large-N f expansion, higher loop diagrams could have contributions at the same order in 1/N f . However, this is fortunately not the case. As proof, we have explicitly checked that, up to NLO, no three-loop diagram contributes to the photon polarization, either because they contain a vanishing contribution (see Section 3.3) or because they are of order 1/N 2 f or higher. This requires the check of 361 diagrams in an automated way. This is conducted by generating the expressions for each diagram and then computing only what is in the order of 1/N f , as well detecting subdiagram expressions that vanish because of the generalized Furry's theorem.

ε-Scalar Polarization at NLO
We now consider the NLO correction to the ε-scalar polarization that consists of nine Feynman diagrams labeled (a,b,...,i). Taking into account the fact that mirror conjugate graphs take the same value, we are left with six distinct graphs to evaluate. This can be conducted exactly for all the diagrams and reads which is, as expected, completely gauge-invariant and finite, providing a strong check on our result. Similar to the photon case, we can rewrite the LO + NLO ε-scalar polarization as and where the interaction correction coefficient to the ε-scalar polarization reads Note that in the only case of interest here, SQED 3 (S = 1, n = 2), this result trivially simplifies as which is exactly the same result as the photon correction coefficient in the SQED 3 case, as shown in (89), as expected from such a supersymmetric gauge-invariant quantity.
Again, we also have explicitly checked that none of the 147 three-loop diagrams contributes to the ε-scalar polarization thanks to the generalized Furry theorem and the already resummed one-loop contributions.

Photino Polarization at NLO
The last polarization to consider is the NLO correction to the photino polarization that consists of 14 Feynman diagrams labeled (a,b,...,n). Taking into account the fact that mirror conjugate graphs take the same value, we are left with seven distinct graphs to evaluate. This can be conducted exactly for all the diagrams and reads Summing all the contributions yields the gauge-invariant and finite result Note that this result depends non-trivially on the parameter E, which implies that the ε-scalars are crucial here to ensure that the result is correct, as we will see in the following. The LO + NLO result for the photino polarization can again be written in the form where the interaction coefficient to the photino polarization reads Note that this result is only of interest in the case of SQED 3 (S = 1, n = 2), where it reduces to provided that we allow for the ε-scalars (E = 1). This is again exactly the same result as for the photon and the ε-scalar correction coefficient in SQED 3 . Therefore, we have explicitly checked that, up to NLO, Π meaning that all polarizations of the gauge multiplet are equal, as expected from supersymmetry for such gauge-invariant and finite quantities.
Again, we also have explicitly checked that none of the 234 three-loop diagrams contribute to the photino polarization thanks to the generalized Furry theorem and the already resummed one-loop contributions.

IR-Softened Gauge Multiplet at NLO
We are now in a position to compute the NLO-softened propagators, i.e., of order 1/N 2 f . Their expressions readD µν where we take the infra-red limit p E ≪ e 2 N f , as advertised in (8). Interestingly, we observe the nice property that the LO (64) and NLO (101)-softened gauge-multiplet propagators are simply related via their polarization correction coefficients, i.e., where the tensorial structure of the photon is still given byd µν (p) =ĝ µν − (1 − ξ)(p µpν /p 2 ) thanks to the use of the non-local gauge (see (28)).

Matter-Multiplet Self-Energies at NLO
In this section, we compute the NLO self-energies of the matter multiplet, i.e., for the electron and the selectron at NLO in the 1/N f expansion, i.e., at O(1/N 2 f ) in gQED 3 . In the following, we shall use shorthand notation for the self-energies

Electron Self-Energy at NLO
We first consider the NLO correction to the electron self-energy that consists of 15 two-loop and 3 one-loop Feynman diagrams, altogether labeled (a,b,...,r). Indeed, contributions of the same order in N f with different loop orders are possible now that we have at our disposal both the LO (64) and NLO (101)-softened propagators. Taking into account the fact that mirror conjugate graphs take the same value, we are left with a total of 16 distinct graphs to evaluate. For each one of them, we extract both the momentum and mass parts using parametrization (30a). With all computations conducted, we obtain the results Note that the computation of the three last diagrams leads to the trivial result that they are simply their LO equivalents times their corresponding coefficient with a sign, −C x /N f , thanks to the equality (102). Moreover, C ε will not contribute to anomalous dimensions because the diagram (104o),Σ p2 , is finite. Similarly, C λ does not contribute to the mass anomalous dimension because the diagram, Σ ψ(r) m2 , is exactly zero. Again, we also have explicitly checked that none of the 390 three-loop diagrams contributes to the electron self-energy at NLO.
Summing all the NLO contributions (104), yields the following results where we introduced the useful notation We can now compute the renormalization functions up to NLO for the electron, reading From these, the anomalous dimensions read We will discuss these results once the anomalous dimensions of the superpartner are computed.

Selectron Self-Energy at NLO
We next consider the NLO correction to the selectron self-energy that consist of 15 two-loops and 2 one-loop Feynman diagrams labeled (a,b,...,p). Taking into account the fact that mirror conjugate graphs take the same value, we are left with a total of 14 distinct graphs to evaluate. For each of them, we extract both the momentum and mass parts using parametrization (30b). With all computations conducted, we obtain the results .

(109n)
Again, note that the computation of the last two graph leads to the trivial results that they are simply the one-loop diagram result times the corresponding interaction correction coefficient with a sign −C x thanks to the identities (101). Interestingly, C ε does not contribute at all to the selectron self-energy since there is no one-loop diagram containing an ε-scalar polarization at this order due to the absence of direct coupling between the selectron and the ε-scalar. Again, we also have explicitly checked that none of the 297 three-loop diagrams contributes to the electron self-energy at NLO.

Summing all the contributions (109) reads
where we again use the useful notation We note that ε-scalars contribute to the self-energies, in part, from the polarization correction C λ (this time for the selectron only) but not from C ε (see (93)). We can now compute the renormalization functions up to NLO for the selectron field and mass using defining Equation (34), reading The factors µ 2ε are discussed under Equation (107). Using the definition of the anomalous dimensions (35), we derive the anomalous dimensions for the selectron field and mass, reading We will discuss the results for the different cases in the next section.

Critical Exponents and Observables
In this section, we apply the general results obtained in the previous section to the various QEDs of interest. In particular, we will present the critical exponents and discuss the observables arising from the computation of the polarization operators of the gauge multiple. We will conclude with a study of the stability of the non-trivial IR fixed point at which all physical quantities are computed.

Results for Fermionic QED 3
As a check of our computations, we will first recover well-known results for large-N f fQED 3 . This can be achieved by considering our general results without supersymmetry (S = 0, n = 2), which yields: The field and mass anomalous dimensions correspond exactly to the expressions first found by Gracey with a different method. In particular, the field anomalous dimension was first derived in the Landau gauge in [35] and then in an arbitrary covariant gauge in [36]. The mass anomalous dimension was derived in [36]. The interaction correction coefficient, C fQED 3 γ , was first explicitly computed in [68,76,112]. Hence, our results are in complete agreement with those of the literature.

Results for N = 1 SQED 3
We now consider the novel case of N = 1 SQED 3 , i.e., taking S = 1, n = 2. First, it is interesting to consider the results with arbitrary E to study the effect of DRED. In this case, our general results yield Interestingly, the effect of the ε-scalar is stiff but crucial. Indeed, the quantities γ ψ and γ ϕ , as well as γ m ψ and C γ , are all E-independent up to NLO. Only C λ and γ m ϕ depend on E. Taking E = 1 so that DRED is allowed, these results simplify as Remarkably, the ε-scalars ensure the validity of identities for the polarization operators as well as for the mass anomalous dimensions both verified up to NLO in our calculations. This is a behavior expected from SUSY that physical (gauge invariant) quantities are identical in the same multiplet [113,114]. On the other hand, the field anomalous dimensions for the electron and the selectron are not equal, neither at LO nor at NLO. This is indeed due to the use of a gauge-fixing term that breaks supersymmetry (Wess-Zumino gauge). We recall here that this is expected and not an issue since the breaking of SUSY occurs only for gauge-dependent quantities that are, by definition, non-physical.

Results for Bosonic QED 3
We now consider the second subcase of interest in this article, which is bosonic QED 3 , i.e., taking S = 1, n = 0 and E = 0. In this case, our general results yield Note that the LO results are in accordance with those previously derived in [31][32][33], and to our knowledge, the NLO results (first published in our short paper [95]) are new.

Results for Reduced QED 4,3 (Graphene)
From the above results, we are now in a position to study QED 4,3 . Let us recall that this model is an effective description of graphene at its ultra-relativistic IR fixed point. We access its properties from those of fQED 3 with the help of mapping first proposed in [75]. Comparing the LO-softened photon propagator of fQED 3 (64a) with the corresponding bare photon propagator of QED 4,3 (see, e.g., [115]), yields the following naive map This map is enough to recover the results for the polarization at one and two-loop for QED 4,3 from the polarization of fQED 3 and, therefore, the corresponding correction coefficient C γ . This map is also sufficient to recover the one-loop anomalous dimensions of the QED 4,3 model from the LO result of the fQED 3 model. However, it breaks at two-loop. Indeed, these models, though very similar, have two major differences that manifest at NLO.
First, fQED 3 at NLO is expressed in a non-local gauge, while the QED 4,3 is not. To compensate this effect, it is enough to consider that if the two-loop polarization of fQED 3 is next to the gauge parameter, ξ, it should not be present in the QED 4,3 case. Since the two-loop polarization is proportional to C γ , one can use the additional rule [75] to recover the proper gauge dependence at two-loop in the anomalous dimensions of QED 4,3 .
Secondly, in fQED 3 , we have softened the photon propagator at NLO and computed additional (one loop but NLO) diagrams (see Equations (104n)-(104p), (109m) and (109n)). These diagrams are not present in QED 4,3 and are replaced by diagrams with a simple fermion loop. To take this into account, we should replace the NLO-softened propagator in fQED 3 with the LO one times the regular factor for a fermion loop in QED 4,3 , i.e., −N f . Since the factor between the two propagators is exactly 102a)), the additional needed rule reads [75] as to be applied on the anomalous dimensions only.
Performed carefully, this mapping yields the following results for QED 4,3 which perfectly recover the results of [68,69,76,80,115]. Note that the use of α instead ofᾱ for the polarization is on purpose.
Going further, following [68,76], we can use our results to compute the optical conductivity of graphene at its IR fixed point. Indeed, the polarization of the photon Π µν , finite and gauge invariant for QED 4,3 (hence, physical) can be related to the optical (AC) conductivity of graphene with the Kubo formula where p µ = (p 0 ,⃗ p). Since the parametrization for the photon polarization reads Π µν = (p 2 g µν − p µ p ν )Π and Π QED 4,3 ∼ 1/p E , Formula (125) simplifies as After momentarily restoring the constants ℏ, c and ε 0 for clarity, it can be written as where σ 0g is the well-known universal minimal AC conductivity of graphene. Moreover, following [116], the optical conductivity of graphene is related to its transmittance (T g ) and its absorbance (A g ) via the relation where α = e 2 /(4πε 0 ℏc) is the usual QED fine structure constant. At first order, since N f = 2 (i.e., 8 elementary spinors) and α = 1/137, we obtain an absorbance of Moreover, C QED 4,3 γ is the interaction correction coefficient to this quantity so that we can expand the leading order absorbance to compute corrections, reading From perturbation theory, we expect the next corrections to be even smaller so that the first one can be taken as an error bar for the next ones, i.e., multiplying the NLO factor by ±1, and since α = 1/137, we have numerically that pure standing relativistic graphene has an optical absorbance of Surprisingly, this result is very close to the one found in the experiments, A exp g = (2.3 ± 0.2)% [117,118], even though measurements are in the pseudo-relativistic limit (v F ≈c/300), while our value (131) is computed in the ultra-relativistic limit (with electrons traveling at the speed of light, v F = c) (see related discussions in [30]).

Results for Reduced
As a non-trivial application, we will map our results for SQED 3 to a model of super-graphene, i.e., for SQED 4,3 (see the action (5)). We recall that this model is an effective description of an eventual pure suspended super-graphene material at its ultra-relativistic fixed point. Following the non-SUSY case, the mapping arises from comparing the LO IR-softened gauge propagators of SQED 3 (64) with the propagators of SQED 4,3 derived from, e.g., [93]. These propagators read It is then straightforward to deduce the following naive mapping which is the same as the non-SUSY case up to a factor of two. This allows for accessing the polarization of SQED 4,3 up to two loops and also the anomalous dimensions up to one loop. In order to access the correct two-loop contribution to the anomalous dimensions for this model, as in the non-SUSY case, we first have to cancel the effect of the non-local gauge by using and then cancel the effect of the NLO softening of the gauge propagators by taking in the anomalous dimensions. Performed carefully, this mapping yields the following results Note that use of α instead ofᾱ for polarization is on purpose. These results are in accordance with [93] at one loop. To our knowledge, the two-loop contributions are a new result. Note that [93] considered a super-graphene model on the boundary (on a substrate) such that the coupling α bdry is twice as small than in our case, i.e., α bdry = α/2.
Similar to the non-supersymmetric case, we can derive the optical conductivity of the hypothetical super-graphene. Using the same Kubo formula, as in (125), yields the following minimal AC conductivity of super-graphene which is twice as big than the non-SUSY one (127). From here, a procedure similar to the non-SUSY case yields the following optical absorbance Therefore, the absorbance of ultra-relativistic freestanding super-graphene is twice the value of normal graphene in the same conditions. Amusingly, this result is exactly the same as for bilayer (non-SUSY) graphene, which is experimentally twice the absorbance of graphene, i.e., A g ≈ 4.6%, [116,119].

Stability of the IR Fixed Point
An important question is related to the stability of the non-trivial IR fixed point with respect to radiative corrections. As we have discussed in the Introduction, for all variants of QED 3 that we have studied, this fixed point arises in the large-N f limit and, more precisely, in the limit p E ≪ e 2 N f . As N f decreases, corrections in 1/N f increase, which calls for an examination of how the fixed point is affected.
Following the bQED 3 and fQED 3 cases (see [2,5], respectively), for all the QED 3 models studied in this article, one can define a dimensionless effective charge In the case of our general (gQED 3 ) model, the photon polarization operator takes the following form where C γ encodes the effects of interactions. From (139) and (140), one can define the beta function associated with the effective coupling, g r . Its expression is given by: and displays two fixed points. One of them is the (trivial) asymptotically free UV fixed point, g * r →0. The second one is the (non-trivial) interacting IR fixed point that we are interested in, g * r → − 1/X (see [1,2,4], as well as the more recent [43]). Summarizing g * r = 0 asymptotic UV fixed point, By combining the above results, the non-trivial IR fixed point of the various cases of interest read where the results are accurate up to the NLO in the 1/N f expansion. We see that fQED 3 is the least affected by radiative corrections and that the latter also weakly affects SQED 3 (though three times more than fQED 3 ). In the case of bQED 3 , however, the correction is of the order of 1. A resummation yields g * r bQED 3 so that, despite being shifted by an amount in the order of 1/N f , the fixed point still exists. It would be interesting to have an all-order proof of the existence of the fixed point, but this goes beyond the scope of the present paper.

Dynamical (Matter) Mass Generation
As an application of our results, we now turn to an estimate of N c , the critical number of (s)electron flavors, which is such that for N f > N c , the (s)electron is massless, while for N f < N c , a dynamical mass, with a Miransky scaling [5], is generated, reading As discussed in the model section, at the level of the action, the potentially generated parity-even mass terms (parity-odd masses cannot be dynamically generated [120]) are of the form (7). Let us remark that only the electron mass term breaks the global flavor symmetry. From SUSY, we also expect that m dyn ψ = m dyn ϕ , which we will simply call m dyn .

The (Semi-Phenomenological) Gap Equation
In principle, the critical number of fermion flavors should be derived via the self-consistent resolution of properly truncated (coupled) Schwinger-Dyson (SD) equations. Due to the complexity of the calculations, for decades, this task has been carried out only at the LO, which has resulted in several inconsistencies, such as severe gauge dependence [121,122] and/or broken Ward identities (see also the thesis in [123]). In the case of fQED 3 , following early multiloop works of Nash [6] and Kotikov [9,10], a complete gauge-invariant prescription up to the NLO of the 1/N f -expansion appeared only rather recently in [15] and [16,17] (see also [19] for a recent review). In [75], the results were then mapped to QED 4,3 , thereby extending the LO results of [67] to the NLO in α.
The systematic approach reviewed in [19] alleviates doubts about the validity of the SD equation approach to access the non-perturbative regime of dynamical mass generation. Nevertheless, it is very technical and difficult to apply to, e.g., SQED 3 where SUSY leads to a dramatic increase in the number of graphs with respect to fQED 3 where γ m = γ m1 + γ m2 + ··· expanded either in loop or large-N f . In (146), the electron dynamical mass scaling is m dyn ∼ p −b E , and dynamical generation occurs when b becomes complex. Actually, as was already noted in the early literature on four-dimensional models (see, e.g., [124][125][126][127][128]) and reconsidered recently [15,115], the form of this gap equation can be deduced from the UV asymptotic behavior of the fermion propagator. In [115], it was argued by the present authors that the gap equation may be quadratic in b at all loop orders and, therefore, semi-phenomenologically written non-pertubatively in with the gauge-invariant γ m as the only input. In this case, b becomes complex for (b − (d e − 2)/2) 2 < 0, yielding the criterion from which the corresponding critical number of fermion N c (and possibly the critical coupling α c ) can be computed. Note that, if γ m would be known exactly, the gap equation would then simply yield γ m (N c ) = (d e − 2)/2. However, when the mass anomalous dimension is known only perturbatively up to a certain order, the gap Equation (148) accordingly needs to be properly truncated, i.e., where γ m = γ 1m + γ 2m + ···. Hence, we should use and then solve K(N c ) = 0 (or possibly K(α c ) = 0). Since γ m is gauge-invariant by construction, the resulting N c will automatically be gauge-invariant too. Moreover, as it is built from the SD formalism, it can be truncated to the accuracy at which γ m is known (Equation (147) reduces to (146) at the LO in 1/N f ). From this polynomial equation, we will obtain multiple solutions for N c . The physical N c will be taken as the largest real solution that is found, which is in accordance with perturbation theory. Though semi-phenomenological, such an approach is straightforward and completely gauge invariant.
For completeness, we provide numerically the mass anomalous dimensions that we obtained in Table 2.

Results for (S)QED 3
In the following, we shall only focus on the electron mass generation and not its superpartner. Indeed, in the case of bQED 3 with N f scalars, we did not find any evidence for dynamical scalar mass generation in bQED 3 , suggesting that for that model, either via the SD method or via the effective gap equation method. Note that the picture seems different if one allows a non-zero quartic coupling λ(|ϕ| 2 ) 2 in three dimensions (see, e.g., [129]), where they obtained N bQED 3 c (λ ̸ = 0) = 6.1 ± 1.95 from fixed-point collision in a four-loop expansion combined with advanced resummation techniques. The situation seems to be also different in four dimensions (see [130]).
On the other hand, for SQED 3 (similar to the four-dimensional case (see [131,132])), we find a possibility that a selectron mass can be induced by the electron condensate, if the latter exists. As we will see in the following, our results suggest that electrons do not condense in SQED 3 .
coinciding with the Landau gauge result of [46], which is indeed the good gauge for SQED 3 . The LO result suggests that an electron mass is generated for N f = 1, thus seemingly breaking both flavor and SUSY symmetries. We find that higher-order corrections dramatically change this picture. Indeed, truncating the gap equation at the NLO of the 1/N f expansion, we find that Such a complex value arises because of the negative NLO contribution (due to the selectron) to the mass anomalous dimension (116c), as shown in Table 2, which prevents the gap equation from having any real valued solution. This calls for a 1/N 3 f computation that is clearly outside the scope of this article. So, in order to overcome this difficulty, we shall proceed with a resummation of the seemingly alternating asymptotic series. A simple Padé approximant [1/1] of (116c) leads to Using this new improved value to solve the gap equation non-perturbatively, i.e., γ m ψ (N c ) = 1/2, yields This result is strong evidence that beyond the LO of the 1/N f expansion, no dynamical (parity-even) mass is generated for the electron in N = 1 SQED 3 . Though a dynamical breaking of SUSY may take place in SQED 3 (the Witten index is not well-defined with massless matter fields; see, e.g., ref. [133] and references therein), the absence of any electron condensate suggests that SUSY is preserved, in accordance with our perturbative result γ m ψ = γ m ϕ up to NLO.
We then focus on the case of fQED 3 (S = 0, n = 2), for which the gap equation is known exactly up to NLO [15,17,19]. The same procedure, this time using (114b) for the mass anomalous dimension, leads at LO to and at NLO to which is in accordance with [15,17,19]. Although the problem of a complex N c is not encountered in this case (because the NLO term in (114b) is positive, as shown in Table 2), we still provide for completeness the improved N c value obtained with resummation, i.e., As expected from radiative correction effects, this value is smaller than the exact NLO one but still quite close, in accordance with the stability of the critical point. In striking contrast with both SQED 3 and bQED 3 , this suggests that a dynamical (flavor-breaking and parity-even) mass is radiatively generated for the electron in fQED 3 for small values of N f , i.e., for N f = 1 and 2. This new improved value (157) is to compare with the extensive fQED 3 DχSB literature, see Table 3, where seemingly all values between 0 and 4 (even infinite in some early studies) have been obtained over four decades.  On the other hand, results from lattice simulations are inconsistent. This may partly be due to the fact that, as N f = 2 is close to N c , the dynamically generated mass is so small (see estimate and discussion in [15]) that it is difficult to extract from lattice simulations.

Results for (S)QED 4,3 (Graphene and Super-Graphene)
From the results, we have obtained for QED 4,3 , in particular, the mass anomalous dimension at two loops (124b), we may apply our semi-phenomenological gap equation formalism and derive the critical coupling constant and critical fermion flavor number. This was conducted in [115] with results that are in complete agreement with those derived from the SD equation formalism [75]. We review them in the following and then carry on with the supersymmetric case.
The computations require the use of an RPA-like procedure, which consist of resumming the two-loop N f dependency (see [75,115] for more details), yielding Similarly, at two loops, the following can be obtained which, for the range of allowed non-zero values of N f , (160) yield We recall that, in the case of graphene, we are interested in N f = 2 because graphene has a total of eight spinors (two cones/sub-lattices × two valley/chirality × two spins). Moreover, in the ultra-relativistic limit we are interested in, graphene has a weak coupling constant, α ∼ 1/137. The result N QED 4,3 c = 3 is unreachable because it is valid in the limit α→∞. Moreover, the critical coupling α c (N f = 2) ≈ 1.2 is much larger than α ∼ 1/137. Hence, graphene remains (semi-)metallic in the ultra-relativistic limit. This agrees with the results originally derived in [75,115]. It is also compatible with experiments on graphene that do not find any evidence for a metal-to-insulator transition.
Similarly, we derive the results for the SQED 4,3 critical coupling and fermion flavor number, corresponding to a phase transition in the ultra-relativistic limit of freestanding super-graphene. Following the RPA-like procedure introduced in [75,115], we obtain For the range of allowed values of N f , it leads numerically to Unfortunately, at two loops, the result is non-physical, despite trying RPA-like or Padé resummations. This is probably a parity effect (like the four-loop approach in QED 4 , see [115]), we then settle for the one-loop approach. Since the relevant value for super-graphene is also N f = 2, and that N SQED 4,3 c = 1.6211 at one loop, and that we expect higher-order corrections to lower N c , we expect that N f = 2 will always be above N c in SQED 4,3 . Hence, super-graphene is even further away from the insulating phase than graphene. For completeness, we provide in Figure 3 the phase diagram of (super-)graphene.  Here, insulator refers to an excitonic insulating phase, while metal refers to a semimetallic phase.

Conclusions
In this article, we have reviewed the critical properties of several variants of three-dimensional QED with fermionic (fQED 3 ), bosonic (bQED 3 ) and minimally supersymmetric (SQED 3 ) charged matter. All these cases were considered in a unified way with the help of a general gQED 3 model. Reduced QED models and their supersymmetric extension were also studied in relation to graphene (QED 4,3 ) and super-graphene (SQED 4,3 ).
In the general framework provided by the gQED 3 model, we performed a complete analytical perturbative computation of matter and gauge field anomalous dimensions at the LO and the NLO in the large-N f expansion, in the DRED scheme and for arbitrary covariant gauge fixing. All these quantities correspond to the critical exponents of the considered models at the non-trivial IR fixed point that arises in the large-N f limit. Expanding on our previous (short) paper [95], all calculation details were provided. Along the way, we added perturbative proof of the Furry theorem, generalized for these models. We also studied thoroughly the effect of DRED and showed its crucial importance in ensuring that the theory remains SUSY invariant.
All of our results have a transcendental structure that is similar to that known in the case of fQED 3 . There are, however, noticeable quantitative differences, with radiative corrections having a tendency to increase vacuum polarization in bQED 3 with respect to fQED 3 while acting oppositely on the mass anomalous dimension. The case of SQED 3 is, somehow, intermediate between fQED 3 and bQED 3 with, in particular, a tendency of the bosonic contribution from the selectron to, on the one hand, increase the overall photon polarization and, on the other hand, decrease the overall electron mass anomalous dimension.
As a first application of our results, we computed the optical conductivity of super-graphene (SQED 4,3 ) and showed that it leads to an optical absorbance of ∼ 4.6%. This result is exactly twice the absorbance of usual (non-SUSY) graphene (QED 4,3 ) and is a direct consequence of the enhanced effect of interactions that, if ever realized, SUSY would bring on the optical absorbance. Another application was devoted to the study of a potential dynamical (matter) mass generation. This nonperturbative phenomenon again revealed a marked difference between fQED 3 and bQED 3 . In fQED 3 , a flavor-breaking parity-invariant mass is generated for N f ≤ 2 (in terms of four-component spinors), while in bQED 3 , we did not find any evidence for a dynamically generated scalar mass. In all cases, our results indicate that the addition of SUSY to an Abelian gauge theory seems to suppress (rather than enhance) dynamical mass generation. In the case of SQED 3 , the value found for the critical electron flavor number, N c , that is such that for N f < N c a dynamical mass for the electron would be generated, is given by N c = 0.39 (in terms of four-component spinors). Contrary to fQED 3 , this strongly suggests that N = 1 SQED 3 remains in an interacting conformal phase for all values of N f . In the case of SQED 4,3 , we found N c = 1.62, which is indeed lower than the corresponding value for QED 4,3 , N c = 3.04. Note that these results hold only at strong coupling. Graphene at its IR fixed point is weakly coupled (α ∼ 1/137) and, thus, is deep in a semimetallic phase, which is in qualitative agreement with the experiments in actual samples.
This diamond integral is the only one we have to consider at two loops for two-point massless diagrams as it encompasses all the other two-loop topologies. Indeed, since lines with zero index contract like we have that all possible two-loop sub-topologies can be written using the diamond diagram with well-chosen zero indices, reading where these three sub-topologies have been reduced to the exactly known one-loop master topology, G(d,α,β). Indeed, the first one (double bubble) is obviously the multiplications of two one-loop bubbles, and upon closer inspection, the two others (eye and sunset) are convolutions of one-loop bubbles.
We are then left with the case where all α i are non-zero, like, for example, J(d,⃗ p,1,1,1,1,1). In this case, IBP techniques can be used [152][153][154] and allow the derivation of identities of the form where α ± i = α i ± 1. Several similar IBP identities can be derived, and altogether, they form a powerful reduction algorithm. One can show that, ultimately, every two-loop integral J(d,⃗ p,α 1 ,α 2 ,α 3 ,α 4 ,α 5 ), with integer indices α i , can be reduced as a linear combination of a small set of master integrals. The implementation of the IBP identities and the reduction process can be conveniently automated with the Mathematica versatile package LiteRed by Roman Lee [155,156].