Field theoretic renormalization study of interaction corrections to the universal ac conductivity of graphene

The two-loop interaction correction coefficient to the universal ac conductivity of disorder-free intrinsic graphene is computed with the help of a field theoretic renormalization study using the BPHZ prescription. Non-standard Ward identities imply that divergent subgraphs (related to Fermi velocity renormalization) contribute to the renormalized optical conductivity. Proceeding either via density-density or via current-current correlation functions, a single well-defined value is obtained: $\mathcal{C}= (19-6\pi)/12) = 0.01$ in agreement with the result first obtained by Mishchenko and which is compatible with experimental uncertainties.


Introduction
Transport properties of graphene and similar planar Dirac liquids have been the subject of extensive studies for more than a decade now, see Ref. [1] for a review. Of central interest for charge transport is the conductivity, σ, which is in general a complicated function of, e.g., frequency (ω), momentum ( q ), temperature (T ), chemical potential (µ) and scattering rates (Γ). A remarkable feature of the ideal case of an intrinsic (µ = 0) disorder-free graphene monolayer is that, despite the vanishing density of states at the Fermi points, the relativistic-like nature of the charge carriers [2,3] yields a non-zero universal ac conductivity in the collisionless limit (ω T, Γ): 1 (1.1) This result, which was predicted long ago to hold for free Dirac fermions [4][5][6], agrees to within 1-2% with experiments in the optical regime, e.g., at ω ∼ 1eV (visible range of the spectrum) [7,8]. Still in the case of free fermions, adding other factors, e.g., T , µ, Γ, ..., is rather non-trivial but leads to results [9,10] which fit quite well the experimental 1 Throughout this paper, we work in units where = kB = 1. data [7,8]. This is rather surprising because the long-range Coulomb interaction among charge carriers is not only unscreened but also supposed to be strong as witnessed by the fine structure constant of suspended graphene: which is of the order of unity due to the fact that v ≈ c/300 where v and c are the Fermi and light velocities, respectively. Moreover, it is well known that Kohn's theorem [11] does not apply to pseudo-relativistic systems thereby allowing electron-electron interactions to affect Eq. (1.1), see also the recent Ref. [12] and references therein for more. There has therefore been extensive theoretical works during the past decade devoted to understanding the intriguing effect of interactions on the optical conductivity of graphene in the collisionless limit, see, e.g., Refs. [13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29], see also [30] for a short review. The latter can be defined via a density-density correlation function: where, in real time, Π 00 (t, q ) = T ρ(t, q )ρ(0, − q ) , ρ is the charge density and T the timeordering operator. Equivalently, from current conservation, it can also be defined via a current-current correlation function (Kubo formula): where, in real time, K ij (t, q ) = T j i (t, q )j j (0, − q ) and j is the charge current. In the case of weak short-range interactions, it was rigorously established that no interaction correction arises [13]. 2 However, no exact result is available in the more realistic case of long-range interactions. Analytically, though interactions are strong, the problem is generally considered with the help of perturbation theory with a focus on the lowest order corrections to Π 00 (q) and K ij (q): where the numbers C andC are the so-called first order interaction-correction coefficients.
On physical grounds, one expects that C =C, independent on the method used. It turns out that, despite the apparent simplicity of the task, different theoretical results have been found in the literature during the past ten years, see Tab. 1 for a summary of some results.
Following Refs. [23,25], the main three values read: As can be seen from Tab. 1, the most commonly accepted result is, up to date, the value C (2) of Mishchenko [15] since it has been recovered by a majority of groups. Incidentally, this is also the only result, among those of Eqs. (1.6), which is consistent with the experimental uncertainties [7,8] as C (2) α g ≈ 2%.
At this point, we note that there is a limit where the result for the interaction correction coefficient does not raise any doubt: the deep infra-red (IR) limit corresponding to the Lorentz-invariant fixed point [31]. The later arises from the long-range Coulomb interaction among the Dirac fermions which enforces the flow of the Fermi velocity, v ≈ c/300 at the present experimentally accessible scales, to the velocity of light, c, in the IR with a corresponding flow of the fine structure constant, α g ≈ 2.2, to the usual fine structure constant, α ≈ 1/137, in the IR. The IR fixed point therefore corresponds to an ultrarelativistic limit with fully-retarded interactions (v/c → 1). The corresponding conductivity reads [19,21,29,30]: Of course, at the fixed point α = 1/137 and the product C * α ≈ 10 −4 is very small leading to almost unobservable effects. However, it is interesting to note that C * = 0.056, a value which is of the same order of magnitude as C (2) . 3 The result of Eq. (1.7) was derived in Refs. [19,21] with the help of multi-loop techniques together with conventional renormalization. In Ref. [29], see also the review [30], the Bogoliubov-Parasiuk-Hepp-Zimmermann (BPHZ), see Refs. [32][33][34] as well as Ref. [35] for a textbook, prescription was used. 4 This is a powerful renormalization prescription, which allows to systematically construct counter-terms on a diagram by diagram basis and can be generalized to higher 3 The interest in comparing C (2) to C * (rather than C (2) αg to C * α) comes from the fact that there is a more general model (model I in the terminology of Ref. [29]) describing interactions in graphene which is valid for arbitrary v/c. So there is actually a non-trivial interaction correction function C(v/c) which encodes relativistic corrections (the dependence of the fine structure constant on v is trivial). Presently, only the limiting values C(v/c → 0) = C (2) (see the proof in the following pages) and C(v/c → 1) = C * are known. It is surprising that these values are of the same order as if C(v/c) was only weakly dependent on v/c. A study of C(v/c) for arbitrary v/c is beyond the scope of the present paper and we leave it for our future investigations. 4 Let's also note the more recent Hopf algebraic formulation [36] of renormalization (see Ref. [37] for a recent review). Its application to our model is beyond the scope of our present study. Table 1. Some values of C obtained over the years together with elements of the different methods used. In case of numerical simulations, we cite the numerical value obtained whenever available and when it slightly differs from the main three results found in the literature, C (i) (i = 1, 2, 3), Eqs. (1.6). For the sake of completeness, the value at the IR fixed point, C * , has been also added. DR is for dimensional regularization. CR is for conventional renormalization. The result derived in this paper appears in red.
orders. It was then found in Refs. [29,30] that, at the fixed point, divergent subgraphs (due to wave function renormalization) cancel each other due to standard, i.e., similar to those appearing in relativistic quantum electrodynamics (QED), Ward identities yielding the fixed point value of Eq. (1.7). Let's also remark that at the fixed point α 1 and the result (1.7) is reliable. But clearly, at higher energy scales such as those relevant to experiments, results based on the loop expansion (1.5) may receive strong corrections from higher orders and the perturbative approach is rightly questionable. A way to overcome this difficulty may be to reorder the perturbative series in the form of a 1/N -expansion or, in other words, to perform a random phase approximation-like resummation, see, e.g., Ref. [38] for an attempt to carry out next-to-leading order (NLO) computations. This is a very interesting suggestion which certainly requires further study. In the following, we shall pursue a more modest goal and restrict ourselves to the study of Eq. (1.5) without any additional resummation. The reason is that, in our opinion, the problem of an accurate evaluation of NLO corrections cannot reasonably be addressed before a full understanding of the first few orders of the loop expansion is achieved.
In the present paper, we reconsider the computation of the first order interaction correction to the minimal conductivity of graphene in the standard but more subtle nonrelativistic limit with instantaneous interactions (v/c → 0). Our analysis clarifies the fact that the origin of the discrepancy between the different results found in this limit, the C (i) (i = 1, 2, 3) of Eqs. (1.6), does not lie in the regularization method but, more simply, in the renormalization procedure itself. This was first revealed in Ref. [25] with a proof obtained with the help of conventional renormalization. In the following, we will give a stronger proof based on the BPHZ renormalization prescription. In particular, we will explicitly show that, in the present non-relativistic case, non-standard Ward identities imply that divergent subgraphs (related to Fermi velocity renormalization) do contribute to the renormalized optical conductivity. Proceeding either via density-density or currentcurrent correlation functions and properly taking into account of these counter-terms will provide a clear explanation for why radiative corrections to the optical conductivity of graphene are finite and perfectly well determined. 5 Anticipating the conclusion, we will show that our approach is in favour of the result first derived by Mishchenko [15]. Our final result reads: where r stands for renormalized and C originates from one-loop counterterms.
The paper is organized as follows. In Sec. 2, we introduce the basic model, the corresponding Feynman rules and the notations that will be used throughout the text. In Sec. 3, we briefly study the model at one-loop (self-energies and Ward identities) introducing the basic ingredients required for the two-loop computation. In Sec. 4, we focus on the optical conductivity at two-loop and compute the counter-terms with the help of the BPHZ prescription. In Sec. 5 we summarize our results and conclude. Some useful integrals are provided in App. A and App. B contains some additional illustration of the singular subgraphs in two-loop diagrams.

Model, Feynman rules and renormalization
The effective low-energy action (model II in the terminology of Ref. [29]) that we wish to consider reads: where v and e are the bare Fermi velocity and charge, respectively, ψ σ is a four-component spinor field describing a fermion of specie σ (σ = 1, · · · , N F and for graphene N F = 2) and A 0 is the gauge field mediating the instantaneous Coulomb interaction. The Dirac matrices, γ µ = (γ 0 , γ ) satisfy the usual algebra {γ µ , γ ν } = 2g µν , with metric tensor g µν = diag(+, −, −). From Eq. (2.1), the bare momentum space fermion propagator reads: where we use pseudo-relativistic notations with p µ = (p 0 , v p ). Because of the absence of any retardation effect (non-relativistic limit corresponding to v/c → 0), the effective photon propagator reduces to the instantaneous Coulomb interaction and reads: Moreover, vector photons decouple and the bare vertex reduces to the temporal part: For our future purposes, it will be nevertheless convenient to define Γ µ = (Γ 0 , Γ ) such that Γ µ 0 = γ µ = (γ 0 , γ ) and employ the following graphical notations: In the following, we will use conventional dimensional regularization [40]. The above Feynman rules stay the same but momenta, Dirac matrices and metric tensor are extended to span a D e -dimensional space (keeping Tr [1] All bare parameters and fields are then related to renormalized ones via renormalization constants in a standard way: where: and µ is the so-called renormalization scale in the minimal subtraction (MS) scheme. The latter is related to the corresponding parameter µ in the modified minimal subtraction (MS) scheme with the help of: where γ E is Euler's constant. The attractive feature of the MS scheme is that the renormalization constants take the simple form: is the renormalized coupling constant and l runs over the number of loops at which ultra-violet (UV) singularities are subtracted. The Z x do not depend on momentum; furthermore, the dependence on µ is only through α gr . So the Z x depend only on α gr (µ) and ε. The basic correlation functions of the model are renormalized as follows: where the notation implicitly assumes that, besides the Z x , only the renormalized quantities depend on µ.
For the model of Eq. (2.1), once integrated over the third spacial direction the effective action of the free gauge field becomes non-local. This implies that the gauge field is not renormalized: Z A 0 = 1 [41][42][43]. As a consequence, the charge is not renormalized: Z e = 1. The renormalization of the coupling constant is therefore entirely due to the renormalization of the velocity which is the only running parameter of the model: The velocity β-function is then defined as: where the subscript B indicates that bare parameters, which do not depend on µ, are fixed.

One-loop analysis
We start by analyzing model of Eq. (2.1) at one-loop.

One-loop fermion and photon self-energies
Let's first recall some elementary results for clarity, see Ref. [25] for more details. The one-loop fermion self-energy, Fig. 1a, is defined as: where d e = 1 + D e is the space-time dimension. The following parametrization is useful: which is such that: Using this parametrization together with the standard rules for integrating massless Feynman diagrams straightforwardly yields: where G(α, β) is defined in Eq.(A.1) and Σ 1ω (k 2 ) = 0. The factor G(1/2, 1/2) contains the pole 1/ε γ at ε γ → 0 and is such that: Performing the ε γ expansion [25] then yields: where L k = log(| k | 2 /µ 2 ). From the above results, we see that there is a momentum but no frequency-dependence of the one-loop fermion self-energy due the instantaneous nature of the interaction. At one-loop, there is therefore a Fermi velocity renormalization [31] but no wave-function renormalization: From Eq. (2.13b) we see that the velocity beta-function is negative: implying that Fermi velocity grows in the infrared [31]. Graphically, these results can be summarized as follows: where the pseudo Lorentz structure of the diagram in the brackets has been projected out, i.e., the displayed graph corresponds to Σ 1k (k 2 ) and not Σ 1 (k), and the K operator is defined as: (3.10) Similar notations and conventions will be used in the following. The one-loop photon self-energy, Fig. 1b, is defined in the usual way as: Focusing on Π 00 , performing the integrations and the ε γ expansion [25] yields: which is finite as expected and where L q 0 = log(−q 2 0 /(4v 2 µ 2 )). Combining Eqs. (1.3) and (3.12), we arrive at the (renormalized) one-loop conductivity: which, for N F = 2, corresponds to the well known universal ac conductivity σ 0 , Eq. (1.1). We may proceed in a similar way with the help of the Kubo formula Eq. (1.4). The transversality of Π µν allows to parametrize it as follows: . (3.14) Then:σ Performing the integrations and the ε γ expansion [25] yields:

One-loop fermion-photon vertex and Ward identities
At this point it is also instructive to look at the vertex part: 1 is the one-loop correction, see Fig. 2a. The latter is defined as: We first consider the temporal part, Fig. 2b. In order to single out it's UV divergent part we will evaluate it for k = q = 0. Performing the trace and going to euclidean space (q 0 = iq E0 ) yields: where the frequency integral vanishes identically, see App. A for some useful integrals. The temporal part of the vertex is therefore not renormalized at one-loop. Graphically, this result can be summarized as follows: Together with Eq. (3.7b), this implies that: Z 1ψ = Z 1Γ 0 = 1 and the Ward identity, Z ψ = Z Γ 0 , is (trivially) satisfied at one-loop. Let's now turn to the vector part of the vertex, see Fig. 2c, focusing again of the case where k = q = 0. Performing the trace and going to euclidean space yields: where now the frequency integral is non-zero. Performing the remaining integrations, we arrive at: which shows that the vector part of the vertex is UV singular (m is just an arbitrary IR regulator). Extracting the pole part, the corresponding renormalization constant together Figure 3. Two-loop vacuum polarization, Π µν 2 , diagrams. with its graphical representation are given by: At one-loop, this result is consistent with the Ward identity: Z Γ = Z v Z ψ which may be graphically represented as: Notice that the peculiar Ward identities Eqs. (3.20) and (3.24) are rather unusual with respect to those which can be found in the case of usual (Lorentz-invariant) QEDs. This will play a crucial role in deriving the correct interaction correction to the optical conductivity in the non-relativistic limit.

Optical conductivity at two loops
We now proceed on computing the 2-loop corrections displayed on Fig. 3: Π µν 2 (q) = 2Π µν 2a (q) + Π µν 2b (q) where Π 2a is the so-called self-energy correction and Π 2b is the so-called vertex correction. The latter are defined in the usual way as: where the fermion self-energy was defined in Eq. (3.1) and the fermion-photon vertex in Eq. (3.18). For completeness, the diagrams associated to Π 00 (a) are displayed on Fig. 4. Our convention for the conductivity will be that: and similarly forσ 2 .

Density-density correlation function approach
The computation of the self-energy diagrams displayed on Fig. 4 has been performed in Ref. [25] (see also App. B). We briefly recall the final results for clarity: where, with two-loop accuracy, the coupling constant is the renormalized one. Combining these results with Eqs. (4.2) yields: from which the total (bare) optical conductivity reads: with a (bare) interaction correction coefficient corresponding to C (3) . We now proceed on computing the renormalized conductivity using the BPHZ prescription. According to the latter, the renormalized diagrams contributing to the density-density correlation function, Π 00 2 α r (q) (α = a, b), are related to the bare ones as follows: where R is the so-called incomplete R-operation because it subtracts only the subdivergences and Π 00 2α may be graphically represented as: The peculiarity of the present non-relativistic theory is that the one-loop fermion selfenergy and fermion-photon vertex subgraphs appearing in Eq. (4.7) are not related by a Ward identity and therefore do not cancel each other (contrarily to what happens in usual QED). The case of Π 00 2b is trivial: this diagram is finite overall so KΠ 00 2b = 0 and, from Eq. (3.20), its subgraphs is also finite so: Π 00 2b = 0. This diagram is therefore absolutely convergent in Weinberg's sense [44] so that: The case of Π 00 2a is more interesting: this diagram is also finite overall, so: KΠ 00 2a = 0. However, its subgraph is divergent, see Eq. (3.9), and needs to be subtracted. In order to compute Π 00 2a we go to the integral representation of Eq. (4.7a) as the operation does not reduce to a simple multiplication (the diagram is not logarithmic and care must be taken in projecting out its pseudo-Lorentz structure). Performing the trace, going to euclidean space, integrating over frequencies and taking the q → 0 limit leads to: where m = q E0 /2v. Substituting Eq. (3.9) and computing the integral using the formulas of App. A yields: D e e 2γ E εγ B(1, 1/2) , (4.10) which is finite in the limit ε γ → 0 and reads: The fact that both Π 00 2a and Π 00 2a are finite implies that K R Π 00 2 a = 0 which was expected since there is no global divergence to subtract. However, the subtraction of the subdivergence brings a finite contribution to the renormalized function: The corresponding renormalized conductivity reads: and is decreased with respect to its bare value Eq. (4.4a) by a quantity: C = −1/4. Hence, the total two-loop renormalized conductivity reads: 14) in accordance with the advertised result for the interaction correction coefficient, Eq. (1.8).

Current-current correlation function approach
We now proceed in a similar way with the self-energy diagrams displayed on Fig. 3. The latter were computed in Ref. [25] (see also App. B). We briefly recall the final results for clarity: where again the coupling constant can be replaced by the renormalized one with two-loop accuracy and L q 0 = log(−q 2 0 /(4v 2 µ 2 )). Notice that the diagrams are now individually divergent contrarily to those of Eqs. (4.3). Singular terms cancel from the sum of Eqs. (4.15a) and (4.15b). Combining these equations with Eq. (3.15) yields the total (bare) optical conductivity at two-loops:σ with a (bare) interaction correction coefficient corresponding once again to C (3) . We now proceed on computing the renormalized conductivity using the BPHZ prescription. Similarly to the density-density case, the renormalized diagrams contributing to the current-current correlation function, K 2 α r (q 0 ) (α = a, b), are related to the bare ones as follows: where K 2α may be graphically represented as: Contrarily to the density-density case, both K 2a and K 2b need to be properly renormalized. Indeed, in the present case, both one-loop fermion self-energy and fermion-photon vertex subgraphs appearing in Eqs. (4.18) are singular. One may wonder at this point if these contributions will cancel each other due to the Ward identity Eq. (3.24). As will be shown in the following, this is not the case. The proof requires some care as the operation does not (necessarily) reduce to a simple multiplication. Notice for example that only the vector part of the fermion-photon vertex subgraph is singular; but the full vertex appears in Eq. (4.18b). So, one needs to be careful in projecting out only the vector component.
In order to do this we will use the integral representation of Eqs. (4.18) We first consider the integral representation of K 2a . Performing the trace, going to euclidean space, integrating over frequencies and taking the q → 0 limit leads to: Substituting Eq. (3.9) and computing the integral yields: The ε γ -expansion reads: Adding the contribution of 2K 2a the pole terms cancel each other and we find that: The fact that KR K 2 a = 0 means that the overall counter-term is zero in accordance with the fact that the only singularity which needs to be subtracted is the one associated with the fermion self-energy subgraph. From Eq. (4.22), we may now derive the expression for the renormalized conductivity associated with Fig. 3a: We now consider the case of K 2b . In order to compute the latter, we find it convenient to go back to the initial definition Eq. (4.1b) replacing Λ µ 1 by γ K Λ 1 / γ with the appropriate γ-matrix contractions. Performing the trace, going to euclidean space, integrating over frequencies and taking the q → 0 limit leads to: The integral is easily computed with the rules for evaluation semi-massive tadpoles, see App. A, and the result reads: Substituting Eq. (3.23) and performing the ε γ -expansion then yields: Adding the contribution of K 2b , the pole terms cancel each other and we find that: (4.27) The fact that KR K 2 b = 0 means that the overall counter-term is zero in accordance with the fact that the only singularity which needs to be subtracted is the one associated with the vector photon-fermion vertex subgraph. From Eq. (4.27), we may now derive the expression for the renormalized conductivity associated with Fig. 3b: Hence, adding Eqs. (4.23) and (4.28), the total two-loop renormalized conductivity reads:σ 29) in accordance with the advertised result for the interaction correction coefficient, Eq. (1.8), and therefore with the result obtained from the density-density correlation function, Eq. (4.14).

Conclusion
In this paper, we have elaborated a field theoretic renormalization approach to the computation of electron-electron interactions in graphene and related planar Dirac liquids in the non-relativistic limit (v/c → 0) with instantaneous Coulomb interaction. As we have seen, the broken Lorentz invariance brings some peculiarities with respect to the relativistic case (non-standard Ward identities) as well as some complications (semi-massive Feynman diagrams instead of the massless ones). Focusing on the optical conductivity of these materials, we have provided a clear proof that radiative corrections affecting this experimentally relevant observable are finite and well determined. Our proof makes use of the powerful BPHZ prescription. It confirms the validity of our previous approach based on conventional renormalization [25] but is considerably more robust and allowed us to perform a refined diagram by diagram analysis. Both the density (where individual diagrams are overall finite) and Kubo formula (where individual diagrams are explicitly singular) approaches were shown to yield a single well-defined result for the two-loop interaction correction to the conductivity: C (2) = (19 − 6π)/12 ≈ 0.013, a result first found by Mishchenko [15] and which is compatible with experimental uncertainties [7,8]. At this point, let's recall that Mishchenko's analysis [15] warns against the use of the Kubo formula with a hard cut-off regulator. Reassuringly, in dimensional regularization no such problem is encountered and both density-density and current-current approaches can be safely used. Our formalism, can be extended to higher orders, 6 other quantities and/or systems of other dimensionality, e.g., the optical conductivity of 3D Dirac materials. It also constitutes a solid base on which the perturbation theory may be optimized in order to deal with the strong-coupling problem. 7 We leave these issues for our future investigations. 6 As we have discussed in the Introduction, for coupling constant αg ≈ 1, perturbation theory is highly questionable. For the polarization operator, the appearance of a small numerical constant C (2) in factor of αg brings a small parameter C (2) αg. This may restore the validity of perturbation theory for the optical conductivity provided that higher order terms may be neglected as well. However, beyond first order, the value of Cn at order n is unknown and it is therefore an open question as to whether or not Cnα n g is small. 7 The computation of NLO corrections in Ref. [38] was approximate as an NLO diagram with two-loop polarization insertion was neglected (see lines below Eq. (9) in that paper). It turns out that this important where, as in the main text, m = q E0 /2v. Notice that this result can be represented in the following form: 3) The one-loop self-energy Σ 1k (| k 2 |) is singular at ε γ → 0 (see (3.6)), but the integral J + (â) ∼ O(ε γ ) in this limit. Indeed: 2)) and in Π 00 2a (q 0 , q → 0) (see Eq. (4.9)) are related with the tadpoles J + (â = 1) and J + (â = 0), respectively. Since at ε γ → 0 J + (â) ∼ (â + 1)ε γ (see Eq. (B.4)), then the counter-term Π 00 2a (q 0 , q → 0) equals half the value of Π 00 2a (q 0 , q → 0) and thus the final result of Π 00 2ar (q 0 , q → 0) in (4.12) is two times less than the one of Π 00 2a (q 0 , q → 0).
Performing the trace, going to euclidean space, integrating over frequencies and taking the q → 0 limit leads to [25]: As it was for Π 00 2a (q 0 , q → 0), this result can be represented in the following form: The one-loop self-energy Σ 1k (| k 2 |) is singular at ε γ → 0 (see (3.6)) and the integral J − (â) is finite in this limit: