# Graphene conductivity: Kubo model versus QFT-based model

Pablo Rodriguez-Lopez, Mauro Antezza

► To cite this version:

Pablo Rodriguez-Lopez, Mauro Antezza. Graphene conductivity: Kubo model versus QFT-based model. 2024. hal-04523215

HAL Id: hal-04523215

https://hal.science/hal-04523215

Preprint submitted on 19 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Graphene conductivity: Kubo model versus QFT-based model

Pablo Rodriguez-Lopez, Mauro Antezza

<figure><img src="image_1.png" /></figure>

To cite this version:

Pablo Rodriguez-Lopez, Mauro Antezza. Graphene conductivity: Kubo model versus QFT-based model. 2024. hal-04523215

HAL Id: hal-04523215

https://hal.science/hal-04523215v1

Preprint submitted on 19 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Graphene conductivity: Kubo model versus QFT-based model

Pablo Rodriguez-Lopez¹,²,\* and Mauro Antezza²,³,†

¹Área de Electromagnetismo and Grupo Interdisciplinar de Sistemas Complejos (GISC), Universidad Rey Juan Carlos, 28933, Móstoles, Madrid, Spain

²Laboratoire Charles Coulomb (L2C), UMR 5221 CNRS-University of Montpellier, F-34095 Montpellier, France

³Institut Universitaire de France, Ministère de l'Enseignement Supérieur et de la Recherche, 1 rue Descartes, F-75231, Paris, France

(Dated: May 7, 2024)

We compare three available models of graphene conductivity: a non-local Kubo model, a local model derived by Falkovsky, and finally a non-local Quantum Field Theory based (QFT-b) model. The first two models are extensively used in the nanophotonic community. All these models are not ab-initio since they contain phenomenological parameters (like Fermi velocity, chemical potential and/or mass gap parameters that depend on the chosen material and possibly on external perturbations), and are supposed to provide coherent results since they are derived from the same starting Hamiltonian. While we confirm that the local model is a proper limit of the non-local Kubo model, we find some inconsistencies in the QFT-b model as derived and used in the literature. In particular, differently from the Kubo model, the QFT-b model shows a plasma-like behavior for the interband transversal conductivity at low frequencies instead of the expected behavior (an almost constant conductivity as a function of frequency ω with a gap for frequencies $\hbar\omega < \sqrt{(\hbar v_F q)^2 + 4m^2}$). We show how to correctly regularize the QFT-b model in order to satisfy the gauge invariance and, once also losses are correctly included, we show that the Kubo and QFT-b model exactly coincide. Our finding can be of relevant interest for both theory, predictions and experimental tests in both the nanophotonic and Casimir effect communities.

# I. INTRODUCTION

Since it was isolated in 2004 [1], the conductivity of graphene has been of great interest due to its potential applications [2][3][4]. There are several different models for the conductivity of graphene that can be classified into at least three different kinds: based on the Kubo formula [5][6][7], Quantum Field Theory based (QFT-b) models [8][9] and based on a hydrodynamic description [10][11][12]. Here we will focus on the comparison between the Kubo formula and QFT-b models. These two families have a deep connection, but are not necessarily equivalent. In Nanophotonics, the rule of thumb is to use the Kubo formula [13], on the other hand, in Casimir physics, some groups use the Kubo formula [14][15][16][17] while another use the QFT-b models [8, 9, 18] that has been used in several dozen papers in last few years, and claimed as much fundamental than the Kubo one since derived "from first principles" [18][19]. In this article we compare Kubo-based and QFT-based models for the conductivity of graphene. We find that, by construction, the Kubo formula provides regularized results that guarantees the fulfillment of the condition $\lim_{\omega \to 0} \Pi^{\mu\nu}(\omega) A_\nu(\omega) = 0$ for constant static $A_\nu(\omega)$ [20] and includes the effect of dissipation of electronic quasiparticles in the conductivity. On the other hand, the QFT-b model as it was developed and used in the literature [8][9][21] not only does not consider unavoidable effects

of losses on the conductivity, but also predicts an additional divergent dissipation-less Plasma behavior that cannot be cured by adding losses. We show that the QFT-b model is derived using a not correctly regularized theory and that, if correctly regularized not only does not provide any non-physical results, but it becomes *exactly identical* to the Kubo model. The Kubo and the QFT-b model are not different, none of them has some particular "first-principle" advantage, they start from exactly the same Hamiltonian and they provide exactly the same final conductivity/polarization result.

This dissipation-less Plasma behavior that appears in the non-local transverse interband conductivity [8][9][21] even when the chemical potential $\mu$ is inside the mass gap of the band spectrum clearly is not an acceptable result in normal materials and it would lead to unobserved dissipation-less currents in graphene, irrespectively of the mass gap, chemical potential, temperature and dissipation, in close analogy to superconductivity [22][23], but without a proper microscopic theory [24][25]. This is the main difference between the two models. However, we find a range of parameters where the 2 models give equivalent results: 1) The longitudinal conductivity for all frequencies; 2) in the local limit and 3) the transverse conductivity, for a sufficiently large (real and/or complex) frequency $\omega$.

In [26] it was shown that the thermal Casimir energy between graphene monolayers is corrected by $\xi_T = \beta\hbar v_F$ instead of $\lambda_T = \beta\hbar c (\beta = (k_B T)^{-1}, v_F \approx c/300$ is the Fermi velocity of electronic excitations in graphene, $c$ is the speed of light, $k_B$ is the Boltzmann constant and $\hbar$ is the reduced Planck constant); the study of the Casimir effect between graphene layers has been studied by the

<footnote>\* pablo.ropez@urjc.es</footnote>

<footnote>† mauro.antezza@umontpellier.fr</footnote>
Graphene conductivity: Kubo model versus QFT-based model

Pablo Rodriguez-Lopez¹,²,\* and Mauro Antezza²,³,†

¹Área de Electromagnetismo and Grupo Interdisciplinar de Sistemas Complejos (GISC),
Universidad Rey Juan Carlos, 28933, Móstoles, Madrid, Spain

²Laboratoire Charles Coulomb (L2C), UMR 5221 CNRS-University of Montpellier, F-34095 Montpellier, France

³Institut Universitaire de France, Ministère de l'Enseignement
Supérieur et de la Recherche, 1 rue Descartes, F-75231, Paris, France

(Dated: May 7, 2024)

We compare three available models of graphene conductivity: a non-local Kubo model, a local model derived by Falkovsky, and finally a non-local Quantum Field Theory based (QFT-b) model. The first two models are extensively used in the nanophotonic community. All these models are not ab-initio since they contain phenomenological parameters (like Fermi velocity, chemical potential and/or mass gap parameters that depend on the chosen material and possibly on external perturbations), and are supposed to provide coherent results since they are derived from the same starting Hamiltonian. While we confirm that the local model is a proper limit of the non-local Kubo model, we find some inconsistencies in the QFT-b model as derived and used in the literature. In particular, differently from the Kubo model, the QFT-b model shows a plasma-like behavior for the interband transversal conductivity at low frequencies instead of the expected behavior (an almost constant conductivity as a function of frequency ω with a gap for frequencies $\hbar\omega < \sqrt{(\hbar v_F q)^2 + 4m^2}$). We show how to correctly regularize the QFT-b model in order to satisfy the gauge invariance and, once also losses are correctly included, we show that the Kubo and QFT-b model exactly coincide. Our finding can be of relevant interest for both theory, predictions and experimental tests in both the nanophotonic and Casimir effect communities.

# I. INTRODUCTION

Since it was isolated in 2004 [1], the conductivity of graphene has been of great interest due to its potential applications [2][3][4]. There are several different models for the conductivity of graphene that can be classified into at least three different kinds: based on the Kubo formula [5][6][7], Quantum Field Theory based (QFT-b) models [8][9] and based on a hydrodynamic description [10][11][12]. Here we will focus on the comparison between the Kubo formula and QFT-b models. These two families have a deep connection, but are not necessarily equivalent. In Nanophotonics, the rule of thumb is to use the Kubo formula [13], on the other hand, in Casimir physics, some groups use the Kubo formula [14][15][16][17] while another use the QFT-b models [8, 9, 18] that has been used in several dozen papers in last few years, and claimed as much fundamental than the Kubo one since derived "from first principles" [18][19]. In this article we compare Kubo-based and QFT-based models for the conductivity of graphene. We find that, by construction, the Kubo formula provides regularized results that guarantees the fulfillment of the condition $\lim_{\omega \to 0} \Pi^{\mu\nu}(\omega) A_\nu(\omega) = 0$ for constant static $A_\nu(\omega)$ [20] and includes the effect of dissipation of electronic quasiparticles in the conductivity. On the other hand, the QFT-b model as it was developed and used in the literature [8][9][21] not only does not consider unavoidable effects

of losses on the conductivity, but also predicts an additional divergent dissipation-less Plasma behavior that cannot be cured by adding losses. We show that the QFT-b model is derived using a not correctly regularized theory and that, if correctly regularized not only does not provide any non-physical results, but it becomes *exactly identical* to the Kubo model. The Kubo and the QFT-b model are not different, none of them has some particular "first-principle" advantage, they start from exactly the same Hamiltonian and they provide exactly the same final conductivity/polarization result.

This dissipation-less Plasma behavior that appears in the non-local transverse interband conductivity [8][9][21] even when the chemical potential $\mu$ is inside the mass gap of the band spectrum clearly is not an acceptable result in normal materials and it would lead to unobserved dissipation-less currents in graphene, irrespectively of the mass gap, chemical potential, temperature and dissipation, in close analogy to superconductivity [22][23], but without a proper microscopic theory [24][25]. This is the main difference between the two models. However, we find a range of parameters where the 2 models give equivalent results: 1) The longitudinal conductivity for all frequencies; 2) in the local limit and 3) the transverse conductivity, for a sufficiently large (real and/or complex) frequency $\omega$.

In [26] it was shown that the thermal Casimir energy between graphene monolayers is corrected by $\xi_T = \beta\hbar v_F$ instead of $\lambda_T = \beta\hbar c (\beta = (k_B T)^{-1}, v_F \approx c/300$ is the Fermi velocity of electronic excitations in graphene, c is the speed of light, $k_B$ is the Boltzmann constant and $\hbar$ is the reduced Planck constant); the study of the Casimir effect between graphene layers has been studied by the

<footnote>\* pablo.ropez@urjc.es</footnote>

<footnote>† mauro.antezza@umontpellier.fr</footnote>
community as a platform to study the effect of temperature on Casimir effect [3].

Several different models of the em response have been used to study the Casimir effect of graphene [2][3][27]. For safety small frequencies, as the Dirac point is close to the chemical potential $\mu$, the tight binding model of graphene can be approached to two $(2+1)D$ massless four-spinor or to a sum of four two-spinors. Due to its simplicity and adequacy to experimental results, the local limit of the Kubo formula, derived by Falkovsky et. al. [28] (see also [29] and [30]) has been widely used [31][32][33]. This model takes into account the (real or imaginary) frequency $\omega$, the chemical potential $\mu$, the temperature $T$, and the dissipation rate $\Gamma = \tau^{-1}$ of the electronic quasiparticles for the Drude conductivity. However, the dissipation for interband transitions and the non-zero mass gap cases are not considered in the model.

In [6], by using the Kubo formula [34] and the two-spinor representation, the generalization to non-local conductivities of [28] for finite mass gaps and non-zero dissipation rate of the interband conductivity was performed. The authors presented closed analytical results for all complex frequencies of the imaginary positive complex plane for the zero temperature limit. From these results, the conductivity for finite temperature is easily obtained.

Another different approach based on Quantum Field Theory (QFT-b) of the four spinor in $(2+1)D$ and on the RPA, like in [35][36][37][38][39][40], gives the conductivity from the polarization operator [8][9][21][33][41][42][43][44][45][46][47][19] to cite a very small set of the literature, in this case, the dissipation rate $\Gamma = \tau^{-1}$ of the electronic quasiparticles is not considered (being equivalent to be set equal to zero), but the results are valid for finite chemical potential $\mu$, temperature $T$ and non-topological Dirac masses $m$ (note that, in [46], the effect of topological Dirac masses was described). The results of those models are "obtained on the solid foundation of quantum field theory and do not use any phenomenology" [19], as a consequence, they lead to the presence of double poles at zero frequency for the dielectric susceptibility "of doubtless physical significance." [19]. Here, as this double pole is translated into a single pole of the conductivity, we argue that the inclusion of electronic quasiparticle's dissipation and of the proper definition of the conductivity tensor in terms of the polarization tensor eliminate the spurious pole.

In this article, we compare the three different derivations of the conductivity of graphene in the small $(k \cdot p)$ limit, and show that the local result of Falkovsky et. al. can be derived from the non-local Kubo result. We show how the QFT-b results are related with the non-local Kubo results, that those results do not coincide and why it is the case.

We hope that this study will clarify the kind of approximations used in each different model, their similarities and differences.

The article is organized as follows: In Sect. II, we introduce the tight-binding model of graphene and the approximations used in the article, we also introduce the notation we will be using throughout the article. In Sect. III we derive and present the formulas used to obtain the polarization and conductivity of graphene in the three different models. In Sect. IV we show how to relate the different quantities obtained in the non-local Kubo model (the longitudinal and transversal conductivities) with the quantities obtained in the QFT-b model (the pure temporal term and trace of the polarization operator). In Sect. V, the non-local model of conductivity derived from the Kubo formula is shown. In Sect. VI the Falkovsky local model of conductivity is presented, and its convergence of the non-local Kubo model is shown. In Sect. VII, the QFT-b model for the polarization (and therefore the conductivity) of graphene is presented. We re-derive the results shown in other articles and explicitly show what is the relation of this model to the non-local Kubo model, when the results coincide and when and why they do not. In Sect. VIII we compare numerically the three different models, highlighting their similarities and differences. We finish in Sect. IX with the conclusions.

## II. TIGHT-BINDING MODEL OF GRAPHENE

In this section we are going to derive the tight-binding model of graphene. The goal is to show what approximations are needed to obtain the $(2+1)D$ Dirac Hamiltonian and the sum of four $(2+1)D$ 2-spinor Hamiltonian. The relation between the two formulas for the conductivity we are discussing pivots around those 2 different representations and their Green functions.

Graphene is a 2D material with a honeycomb lattice, whose unit cell consists on 2 nonequivalent carbon atoms in $sp^2$ electronic configuration. We model the electronic excitations of graphene in the macrocanonical ensemble with a tight-binding model of a bidimensional honeycomb lattice [2][48][49][50][51]. This lattice consists on 2 nonequivalent triangular lattices (denoted as A and B here). The position of the atoms in each sublattice (or of the unit cell) can be specified by a vector $R_{n_1,n_2} = n_1a_1 + n_2a_2$ ($n_i \in \mathbb{Z}$), with lattice vectors

$$a_1 = \frac{\sqrt{3}a}{2} \begin{pmatrix} 1 \\ \sqrt{3} \end{pmatrix}, \quad a_2 = \frac{\sqrt{3}a}{2} \begin{pmatrix} 1 \\ -\sqrt{3} \end{pmatrix}, \quad (1)$$

being $a = 1.42$ Å the carbon-carbon interatomic distance in graphene. The nearest neighbors of an atom of the sublattice A are given by the vectors

$$\delta_1 = a \begin{pmatrix} 0 \\ 1 \end{pmatrix}, \quad \delta_2 = \frac{a}{2} \begin{pmatrix} \sqrt{3} \\ -1 \end{pmatrix}, \quad \delta_3 = \frac{a}{2} \begin{pmatrix} -\sqrt{3} \\ -1 \end{pmatrix}, \quad (2)$$

The reciprocal lattice is also a honeycomb lattice, whose fundamental translation vectors $b_j$ are defined by the re-
community as a platform to study the effect of temperature on Casimir effect [3].

Several different models of the em response have been used to study the Casimir effect of graphene [2][3][27]. For safety small frequencies, as the Dirac point is close to the chemical potential $\mu$, the tight binding model of graphene can be approached to two $(2+1)D$ massless four-spinor or to a sum of four two-spinors. Due to its simplicity and adequacy to experimental results, the local limit of the Kubo formula, derived by Falkovsky et. al. [28] (see also [29] and [30]) has been widely used [31][32][33]. This model takes into account the (real or imaginary) frequency $\omega$, the chemical potential $\mu$, the temperature $T$, and the dissipation rate $\Gamma = \tau^{-1}$ of the electronic quasiparticles for the Drude conductivity. However, the dissipation for interband transitions and the non-zero mass gap cases are not considered in the model.

In [6], by using the Kubo formula [34] and the two-spinor representation, the generalization to non-local conductivities of [28] for finite mass gaps and non-zero dissipation rate of the interband conductivity was performed. The authors presented closed analytical results for all complex frequencies of the imaginary positive complex plane for the zero temperature limit. From these results, the conductivity for finite temperature is easily obtained.

Another different approach based on Quantum Field Theory (QFT-b) of the four spinor in $(2+1)D$ and on the RPA, like in [35][36][37][38][39][40], gives the conductivity from the polarization operator [8][9][21][33][41][42][43][44][45][46][47][19] to cite a very small set of the literature, in this case, the dissipation rate $\Gamma = \tau^{-1}$ of the electronic quasiparticles is not considered (being equivalent to be set equal to zero), but the results are valid for finite chemical potential $\mu$, temperature $T$ and non-topological Dirac masses $m$ (note that, in [46], the effect of topological Dirac masses was described). The results of those models are "obtained on the solid foundation of quantum field theory and do not use any phenomenology" [19], as a consequence, they lead to the presence of double poles at zero frequency for the dielectric susceptibility "of doubtless physical significance." [19]. Here, as this double pole is translated into a single pole of the conductivity, we argue that the inclusion of electronic quasiparticle's dissipation and of the proper definition of the conductivity tensor in terms of the polarization tensor eliminate the spurious pole.

In this article, we compare the three different derivations of the conductivity of graphene in the small $(k \cdot p)$ limit, and show that the local result of Falkovsky et. al. can be derived from the non-local Kubo result. We show how the QFT-b results are related with the non-local Kubo results, that those results do not coincide and why it is the case.

We hope that this study will clarify the kind of approximations used in each different model, their similarities and differences.

The article is organized as follows: In Sect. II, we introduce the tight-binding model of graphene and the approximations used in the article, we also introduce the notation we will be using throughout the article. In Sect. III we derive and present the formulas used to obtain the polarization and conductivity of graphene in the three different models. In Sect. IV we show how to relate the different quantities obtained in the non-local Kubo model (the longitudinal and transversal conductivities) with the quantities obtained in the QFT-b model (the pure temporal term and trace of the polarization operator). In Sect. V, the non-local model of conductivity derived from the Kubo formula is shown. In Sect. VI the Falkovsky local model of conductivity is presented, and its convergence of the non-local Kubo model is shown. In Sect. VII, the QFT-b model for the polarization (and therefore the conductivity) of graphene is presented. We re-derive the results shown in other articles and explicitly show what is the relation of this model to the non-local Kubo model, when the results coincide and when and why they do not. In Sect. VIII we compare numerically the three different models, highlighting their similarities and differences. We finish in Sect. IX with the conclusions.

## II. TIGHT-BINDING MODEL OF GRAPHENE

In this section we are going to derive the tight-binding model of graphene. The goal is to show what approximations are needed to obtain the $(2+1)D$ Dirac Hamiltonian and the sum of four $(2+1)D$ 2-spinor Hamiltonian. The relation between the two formulas for the conductivity we are discussing pivots around those 2 different representations and their Green functions.

Graphene is a 2D material with a honeycomb lattice, whose unit cell consists on 2 nonequivalent carbon atoms in $sp^2$ electronic configuration. We model the electronic excitations of graphene in the macrocanonical ensemble with a tight-binding model of a bidimensional honeycomb lattice [2][48][49][50][51]. This lattice consists on 2 nonequivalent triangular lattices (denoted as A and B here). The position of the atoms in each sublattice (or of the unit cell) can be specified by a vector $R_{n_1,n_2} = n_1a_1 + n_2a_2$ ($n_i \in \mathbb{Z}$), with lattice vectors

$$a_1 = \frac{\sqrt{3}a}{2} \begin{pmatrix} 1 \\ \sqrt{3} \end{pmatrix}, \quad a_2 = \frac{\sqrt{3}a}{2} \begin{pmatrix} 1 \\ -\sqrt{3} \end{pmatrix}, \quad (1)$$

being $a = 1.42$ Å the carbon-carbon interatomic distance in graphene. The nearest neighbors of an atom of the sublattice A are given by the vectors

$$\delta_1 = a \begin{pmatrix} 0 \\ 1 \end{pmatrix}, \quad \delta_2 = \frac{a}{2} \begin{pmatrix} \sqrt{3} \\ -1 \end{pmatrix}, \quad \delta_3 = \frac{a}{2} \begin{pmatrix} -\sqrt{3} \\ -1 \end{pmatrix}, \quad (2)$$

The reciprocal lattice is also a honeycomb lattice, whose fundamental translation vectors $b_j$ are defined by the re-
<table><caption>Table I. Table with the notation used along this article.</caption>
  <tbody>
    <tr>
      <td>k<sub>0</sub> = &omega;</td>
      <td>q<sub>0</sub> = &Omega; = &omega; + i&Gamma;</td>
    </tr>
    <tr>
      <td>k&#x0303;<sub>i</sub> = &hbar;v<sub>F</sub>k<sub>i</sub></td>
      <td>q&#x0303;<sub>0</sub> = &hbar;&Omega; = &hbar;&omega; + i&hbar;&Gamma;</td>
    </tr>
    <tr>
      <td>k&#x0303;<sub>&mu;</sub> = k&#x0303;<sub>0</sub> + &mu;</td>
      <td>k&#x0303;<sub>0</sub> = &hbar;&omega; = i&hbar;&xi; = i&Xi;</td>
    </tr>
    <tr>
      <td>k&#x0303;<sub>&parallel;</sub> = &hbar;v<sub>F</sub>k<sub>&parallel;</sub></td>
      <td>k<sub>&parallel;</sub> = (k<sub>1</sub>, k<sub>2</sub>)</td>
    </tr>
    <tr>
      <td>k&#x0303;<sub>&parallel;</sub> = &hbar;v<sub>F</sub>k<sub>&parallel;</sub> = &radic;(k&#x0303;<sub>1</sub><sup>2</sup> + k&#x0303;<sub>2</sub><sup>2</sup>)</td>
      <td>k<sub>&parallel;</sub> = &radic;(k&#x0303;<sub>1</sub><sup>2</sup> + k&#x0303;<sub>2</sub><sup>2</sup>)</td>
    </tr>
    <tr>
      <td>q&#x0303;<sub>&parallel;</sub> = &hbar;v<sub>F</sub>q<sub>&parallel;</sub></td>
      <td>q<sub>&parallel;</sub> = (q<sub>1</sub>, q<sub>2</sub>)</td>
    </tr>
    <tr>
      <td>q&#x0303;<sub>&parallel;</sub> = &hbar;v<sub>F</sub>q<sub>&parallel;</sub> = &radic;(q&#x0303;<sub>1</sub><sup>2</sup> + q&#x0303;<sub>2</sub><sup>2</sup>)</td>
      <td>q<sub>&parallel;</sub> = &radic;(q&#x0303;<sub>1</sub><sup>2</sup> + q&#x0303;<sub>2</sub><sup>2</sup>)</td>
    </tr>
    <tr>
      <td>k&#x0303;<sub>z</sub> = &hbar;v<sub>F</sub>k<sub>z</sub> = &radic;(k<sub>0</sub><sup>2</sup> - k&#x0303;<sub>&parallel;</sub><sup>2</sup>)</td>
      <td>k<sub>z</sub> = &radic;((<sup>&omega;</sup>&frasl;<sub>v<sub>F</sub></sub>)<sup>2</sup> - k<sub>&parallel;</sub><sup>2</sup>)</td>
    </tr>
    <tr>
      <td>q&#x0303;<sub>z</sub> = &hbar;v<sub>F</sub>q<sub>z</sub> = &radic;(q<sub>0</sub><sup>2</sup> - q&#x0303;<sub>&parallel;</sub><sup>2</sup>)</td>
      <td>q<sub>z</sub> = &radic;((<sup>&Omega;</sup>&frasl;<sub>v<sub>F</sub></sub>)<sup>2</sup> - q<sub>&parallel;</sub><sup>2</sup>)</td>
    </tr>
    <tr>
      <td>&kappa;&#x0303;<sub>z</sub> = i&#x0303;k<sub>z</sub> = &radic;(&Xi;<sup>2</sup> + k&#x0303;<sub>&parallel;</sub><sup>2</sup>)</td>
      <td>k = k<sub>&mu;</sub> = (k<sub>0</sub>, k&#x0303;<sub>&parallel;</sub>)</td>
    </tr>
    <tr>
      <td>&theta;&#x0303;<sub>z</sub> = i&#x0303;q&#x0303;<sub>z</sub> = &radic;(&Xi;<sup>2</sup> + q&#x0303;<sub>&parallel;</sub><sup>2</sup>)</td>
      <td>q = q<sub>&mu;</sub> = (q<sub>0</sub>, q&#x0303;<sub>&parallel;</sub>)</td>
    </tr>
    <tr>
      <td>&gamma;&#x0303;<sub>&mu;</sub> = (&gamma;&#x0303;<sub>0</sub>, &gamma;&#x0303;) = (&hbar;&gamma;<sub>0</sub>, &hbar;v<sub>F</sub>&gamma;)</td>
      <td>&gamma;&#x0303;<sub>i</sub> = v<sub>F</sub>&gamma;<sub>i</sub></td>
    </tr>
    <tr>
      <td>&int;<sub>k</sub> = &int;<sub>-&infin;</sub><sup>&infin;</sup> (dk<sub>0</sub> / 2&pi;) &int;<sub>k</sub></td>
      <td>&int;<sub>k</sub> = &int;<sub>BZ</sub> (d<sup>2</sup>k<sub>&parallel;</sub> / (2&pi;)<sup>2</sup>)</td>
    </tr>
    <tr>
      <td>&delta; = (2m / &theta;&#x0303;<sub>z</sub>)</td>
      <td>&gamma; = (&Xi; / &theta;&#x0303;<sub>z</sub>)</td>
    </tr>
    <tr>
      <td colspan="2">K<sub>&nu;</sub> = k + (&mu;, 0) = (k<sub>&mu;</sub>, k<sub>&parallel;</sub>) = (k<sub>0</sub> + &mu;, &hbar;v<sub>F</sub>k<sub>&parallel;</sub>)</td>
    </tr>
    <tr>
      <td colspan="2">S<sub>&mu;</sub> = K<sub>&mu;</sub> + q<sub>&mu;</sub> = (s&#x0303;<sub>0</sub>, s&#x0303;<sub>&parallel;</sub>) = (k<sub>0</sub> + q<sub>0</sub> + &mu;, k<sub>&parallel;</sub> + q&#x0303;<sub>&parallel;</sub>)</td>
    </tr>
  </tbody>
</table>

lation $\mathbf{a}_i \cdot \mathbf{b}_j = 2\pi\delta_{ij}$, resulting in

$$ \mathbf{b}_1 = \frac{2\pi}{3a} \begin{pmatrix} \sqrt{3} \\ 1 \end{pmatrix}, \quad \mathbf{b}_2 = \frac{4\pi}{3a} \begin{pmatrix} -\sqrt{3} \\ 1 \end{pmatrix}, \qquad (3) $$

The tight-binding Hamiltonian of graphene in real space is

$$ H = \sum_{\mathbf{n} \in \mathbb{Z}^d} \sum_{\lambda, \lambda'} \sum_{j=1}^{|\delta_j|<\delta} \hat{c}_{\lambda, \mathbf{R}_n}^{\dagger} t_{\lambda, \lambda'} (\delta_j) \hat{c}_{\lambda', \mathbf{R}_n + \delta_j}, \quad (4) $$

where $\hat{c}_{\lambda, \mathbf{R}_n}^{\dagger}$ is the creator operator of an electron placed at $\mathbf{R}_n = \mathbf{a}_1 n_1 + \mathbf{a}_2 n_2$, with spin $s = \{\uparrow, \downarrow\}$, triangular sub-lattice (A or B) and orbital ($2p_z$ only in our case) labelled by $\lambda$ in the unit cell. $\hat{c}_{\lambda, \mathbf{R}_n}$ is the annihilation operator of the electron. $t_{\lambda, \lambda'}(\delta_j)$ is the tight-binding coupling between an electron placed at $\mathbf{R}_n$, with spin and orbital $\lambda$ and another electron placed at $\mathbf{R}_n + \delta_j$, with spin and orbital $\lambda'$. Those coefficients can be calculated for each particular case. $\delta$ is the maximum hopping distance between atoms we consider in the model. As the chemical potential $\mu$ is close to the Dirac points, and we are interested in relatively small frequencies, we only take into account the $\pi$ and $\pi^*$ bands in our model and first neighbors coupling only, therefore, the tight-binding Hamiltonian operator of graphene in real space is reduced to

$$ H = t \sum_{s=\pm 1} \sum_{\langle i | j \rangle} \hat{c}_{i, s}^{\dagger} \hat{c}_{j, s}. \qquad (5) $$

For clarity, we will add the contribution of the chemical potential to the Hamiltonian later. $t = V_{pp\pi} \approx 2.8$ eV is the nearest-neighbor hopping energy [2][50], $\hat{c}_{i,s}^{\dagger}$ and $\hat{c}_{j,s}$

are the creator and annihilation operators of electronic excitations with spin $s = \{\uparrow, \downarrow\}$ at site $i$, and in $\langle i | j \rangle$ $i$ run to all the atoms of the lattice while $j$ run all over the nearest neighbors hopping sites of $i$.

In the momentum space, the Hamiltonian becomes

$$ H = \sum_{s=\pm 1} \int_{\mathbf{p}} \hat{\chi}_{s}^{\dagger}(\mathbf{p}) \hat{H}_{s}(\mathbf{p}) \hat{\chi}_{s}(\mathbf{p}), \qquad (6) $$

where we have defined $\int_{\mathbf{p}} = \int_{BZ} \frac{d^2\mathbf{p}}{(2\pi)^2}$ as the momentum integral defined over the Brillouin Zone BZ (see tab. I),

$$ \hat{H}_s(\mathbf{p}) = t \begin{pmatrix} 0 & f(\mathbf{p}) \\ f^*(\mathbf{p}) & 0 \end{pmatrix}, \qquad (7) $$

and the bi-spinor in the valley sub-space as

$$ \hat{\chi}_s(\mathbf{p}) = \begin{pmatrix} \hat{c}_{A,s}(\mathbf{p}) \\ \hat{c}_{B,s}(\mathbf{p}) \end{pmatrix}, \qquad (8) $$

where $\hat{c}_{A,s}(\mathbf{p})$ is the annihilation operator of electronic excitations in the sublattice A with spin $s$ and momentum $\mathbf{p}$, and

$$ \begin{aligned} f(\mathbf{p}) &= -\sum_{j=1}^{3} e^{i\mathbf{p}\cdot\delta_j} \\ &= -e^{-i\mathbf{ap}_y} - 2e^{i\frac{\mathbf{ap}_y}{2}} \cos\left(\frac{\sqrt{3}}{2}ap_x\right). \end{aligned} \qquad (9) $$

Diagonalizing $\hat{H}$ in momentum space gives the energy spectrum as

$$ E_{\lambda}(\mathbf{p}) = \lambda t |f(\mathbf{p})|, \qquad (10) $$

with $\lambda = \pm 1$ representing the conduction ($\lambda = +1$) and valence ($\lambda = -1$) bands respectively, and,

$$ |f(\mathbf{p})|^2 = 1 + 4 \cos^2 \left( \frac{\sqrt{3}ap_x}{2} \right) \\ +4 \cos \left( \frac{3ap_x}{2} \right) \cos \left( \frac{\sqrt{3}ap_y}{2} \right). \qquad (11) $$

Inside the Brillouin Zone defined by the parallelogram $\mathbf{b}_1 \otimes \mathbf{b}_2$, this function is zero at the $K_{\pm}$ points defined as

$$ K_{\eta} = \frac{2\pi}{\sqrt{3}a} \begin{pmatrix} 1 - \frac{\eta}{3} \\ 0 \end{pmatrix}, \qquad (12) $$

with $\eta = \pm 1$ the valley index. As the chemical potential $\mu$ is close to the crossing points between $\pi$ and $\pi^*$ bands at the $K_{\eta}$ points, the dispersion of the bands can be approached as [52]

$$ \hat{H}_s(\mathbf{K}_\eta + \mathbf{k}) = \frac{3at}{2} (\eta k_1 \tau_1 + k_2 \tau_2), \qquad (13) $$
<table><caption>Table I. Table with the notation used along this article.</caption>
  <tbody>
    <tr>
      <td>k<sub>0</sub> = &omega;</td>
      <td>q<sub>0</sub> = &Omega; = &omega; + i&Gamma;</td>
    </tr>
    <tr>
      <td>k&#773;<sub>i</sub> = &hbar;v<sub>F</sub>k<sub>i</sub></td>
      <td>q&#773;<sub>0</sub> = &hbar;&Omega; = &hbar;&omega; + i&hbar;&Gamma;</td>
    </tr>
    <tr>
      <td>k&#773;<sub>&mu;</sub> = k&#773;<sub>0</sub> + &mu;</td>
      <td>k&#773;<sub>0</sub> = &hbar;&omega; = i&hbar;&xi; = i&Xi;</td>
    </tr>
    <tr>
      <td>k&#773;<sub>&parallel;</sub> = &hbar;v<sub>F</sub>k<sub>&parallel;</sub></td>
      <td>k<sub>&parallel;</sub> = (k<sub>1</sub>, k<sub>2</sub>)</td>
    </tr>
    <tr>
      <td>k&#773;<sub>&parallel;</sub> = &hbar;v<sub>F</sub>k<sub>&parallel;</sub> = &radic;(k&#773;<sub>1</sub><sup>2</sup> + k&#773;<sub>2</sub><sup>2</sup>)</td>
      <td>k<sub>&parallel;</sub> = &radic;(k&#773;<sub>1</sub><sup>2</sup> + k&#773;<sub>2</sub><sup>2</sup>)</td>
    </tr>
    <tr>
      <td>q&#773;<sub>&parallel;</sub> = &hbar;v<sub>F</sub>q<sub>&parallel;</sub></td>
      <td>q<sub>&parallel;</sub> = (q<sub>1</sub>, q<sub>2</sub>)</td>
    </tr>
    <tr>
      <td>q&#773;<sub>&parallel;</sub> = &hbar;v<sub>F</sub>q<sub>&parallel;</sub> = &radic;(q&#773;<sub>1</sub><sup>2</sup> + q&#773;<sub>2</sub><sup>2</sup>)</td>
      <td>q<sub>&parallel;</sub> = &radic;(q&#773;<sub>1</sub><sup>2</sup> + q&#773;<sub>2</sub><sup>2</sup>)</td>
    </tr>
    <tr>
      <td>k&#773;<sub>z</sub> = &hbar;v<sub>F</sub>k<sub>z</sub> = &radic;(k<sub>0</sub><sup>2</sup> - k&#773;<sub>&parallel;</sub><sup>2</sup>)</td>
      <td>k<sub>z</sub> = &radic;((<sup>&omega;</sup>&frasl;<sub>v<sub>F</sub></sub>)<sup>2</sup> - k<sub>&parallel;</sub><sup>2</sup>)</td>
    </tr>
    <tr>
      <td>q&#773;<sub>z</sub> = &hbar;v<sub>F</sub>q<sub>z</sub> = &radic;(q<sub>0</sub><sup>2</sup> - q&#773;<sub>&parallel;</sub><sup>2</sup>)</td>
      <td>q<sub>z</sub> = &radic;((<sup>&Omega;</sup>&frasl;<sub>v<sub>F</sub></sub>)<sup>2</sup> - q<sub>&parallel;</sub><sup>2</sup>)</td>
    </tr>
    <tr>
      <td>&kappa;&#773;<sub>z</sub> = i&#773;k<sub>z</sub> = &radic;(&Xi;<sup>2</sup> + k&#773;<sub>&parallel;</sub><sup>2</sup>)</td>
      <td>k = k<sub>&mu;</sub> = (k<sub>0</sub>, k&#773;<sub>&parallel;</sub>)</td>
    </tr>
    <tr>
      <td>&theta;&#773;<sub>z</sub> = i&#773;q&#773;<sub>z</sub> = &radic;(&Xi;<sup>2</sup> + q&#773;<sub>&parallel;</sub><sup>2</sup>)</td>
      <td>q = q<sub>&mu;</sub> = (q<sub>0</sub>, q&#773;<sub>&parallel;</sub>)</td>
    </tr>
    <tr>
      <td>&gamma;&#773;<sub>&mu;</sub> = (&gamma;&#773;<sub>0</sub>, &gamma;&#773;) = (&hbar;&gamma;<sub>0</sub>, &hbar;v<sub>F</sub>&gamma;)</td>
      <td>&gamma;&#773;<sub>i</sub> = v<sub>F</sub>&gamma;<sub>i</sub></td>
    </tr>
    <tr>
      <td>&int;<sub>k</sub> = &int;<sub>-&infin;</sub><sup>&infin;</sup> (dk<sub>0</sub> / 2&pi;) &int;<sub>k</sub></td>
      <td>&int;<sub>k</sub> = &int;<sub>BZ</sub> (d<sup>2</sup>k<sub>&parallel;</sub> / (2&pi;)<sup>2</sup>)</td>
    </tr>
    <tr>
      <td>&delta; = 2m / &theta;&#773;<sub>z</sub></td>
      <td>&gamma; = &theta;&#773;<sub>z</sub></td>
    </tr>
    <tr>
      <td colspan="2">
        K<sub>&nu;</sub> = k + (&mu;, 0) = (k<sub>&mu;</sub>, k<sub>&parallel;</sub>) = (k<sub>0</sub> + &mu;, &hbar;v<sub>F</sub>k<sub>&parallel;</sub>)<br/>
        S<sub>&mu;</sub> = K<sub>&mu;</sub> + q<sub>&mu;</sub> = (s&#773;<sub>0</sub>, s&#773;<sub>&parallel;</sub>) = (k<sub>0</sub> + q<sub>0</sub> + &mu;, k<sub>&parallel;</sub> + q&#773;<sub>&parallel;</sub>)
      </td>
    </tr>
  </tbody>
</table>

lation $\mathbf{a}_i \cdot \mathbf{b}_j = 2\pi\delta_{ij}$, resulting in

$$
\mathbf{b}_1 = \frac{2\pi}{3a} \begin{pmatrix} \sqrt{3} \\ 1 \end{pmatrix}, \quad \mathbf{b}_2 = \frac{4\pi}{3a} \begin{pmatrix} -\sqrt{3} \\ 1 \end{pmatrix}, \qquad (3)
$$

The tight-binding Hamiltonian of graphene in real space
is

$$
H = \sum_{\mathbf{n} \in \mathbb{Z}^d} \sum_{\lambda, \lambda'} \sum_{j=1}^{|\delta_j|<\delta} \hat{c}_{\lambda, \mathbf{R}_n}^{\dagger} t_{\lambda, \lambda'} (\delta_j) \hat{c}_{\lambda', \mathbf{R}_n + \delta_j}, \quad (4)
$$

where $\hat{c}_{\lambda, \mathbf{R}_n}^{\dagger}$ is the creator operator of an electron placed at $\mathbf{R}_n = \mathbf{a}_1 n_1 + \mathbf{a}_2 n_2$, with spin $s = \{\uparrow, \downarrow\}$, triangular sub-lattice (A or B) and orbital ($2p_z$ only in our case) labelled by $\lambda$ in the unit cell. $\hat{c}_{\lambda, \mathbf{R}_n}$ is the annihilation operator of the electron. $t_{\lambda, \lambda'}(\delta_j)$ is the tight-binding coupling between an electron placed at $\mathbf{R}_n$, with spin and orbital $\lambda$ and another electron placed at $\mathbf{R}_n + \delta_j$, with spin and orbital $\lambda'$. Those coefficients can be calculated for each particular case. $\delta$ is the maximum hopping distance between atoms we consider in the model. As the chemical potential $\mu$ is close to the Dirac points, and we are interested in relatively small frequencies, we only take into account the $\pi$ and $\pi^*$ bands in our model and first neighbors coupling only, therefore, the tight-binding Hamiltonian operator of graphene in real space is reduced to

$$
H = t \sum_{s=\pm 1} \sum_{\langle i|j \rangle} \hat{c}_{i,s}^{\dagger} \hat{c}_{j,s}. \qquad (5)
$$

For clarity, we will add the contribution of the chemical
potential to the Hamiltonian later. $t = V_{pp\pi} \approx 2.8$ eV is
the nearest-neighbor hopping energy [2][50], $\hat{c}_{i,s}^{\dagger}$ and $\hat{c}_{j,s}$

are the creator and annihilation operators of electronic
excitations with spin $s = \{\uparrow, \downarrow\}$ at site $i$, and in $\langle i|j \rangle$ $i$
run to all the atoms of the lattice while $j$ run all over the
nearest neighbors hopping sites of $i$.

In the momentum space, the Hamiltonian becomes

$$
H = \sum_{s=\pm 1} \int_{\mathbf{p}} \hat{\chi}_{s}^{\dagger}(\mathbf{p}) \hat{H}_{s}(\mathbf{p}) \hat{\chi}_{s}(\mathbf{p}), \quad (6)
$$

where we have defined $\int_{\mathbf{p}} = \int_{\text{BZ}} \frac{d^2\mathbf{p}}{(2\pi)^2}$ as the momentum
integral defined over the Brillouin Zone BZ (see tab. I),

$$
\hat{H}_s(\mathbf{p}) = t \begin{pmatrix} 0 & f(\mathbf{p}) \\ f^*(\mathbf{p}) & 0 \end{pmatrix}, \quad (7)
$$

and the bi-spinor in the valley sub-space as

$$
\hat{\chi}_s(\mathbf{p}) = \begin{pmatrix} \hat{c}_{A,s}(\mathbf{p}) \\ \hat{c}_{B,s}(\mathbf{p}) \end{pmatrix}, \qquad (8)
$$

where $\hat{c}_{A,s}(\mathbf{p})$ is the annihilation operator of electronic
excitations in the sublattice $A$ with spin $s$ and momen-
tum $\mathbf{p}$, and

$$
\begin{align}
f(\mathbf{p}) &= -\sum_{j=1}^{3} e^{i\mathbf{p}\cdot\delta_j} \nonumber \\
&= -e^{-i\mathbf{ap}_y} - 2e^{i\frac{\mathbf{ap}_y}{2}} \cos\left(\frac{\sqrt{3}}{2}ap_x\right). \tag{9}
\end{align}
$$

Diagonalizing $\hat{H}$ in momentum space gives the energy
spectrum as

$$
E_{\lambda}(\mathbf{p}) = \lambda t |f(\mathbf{p})|, \qquad (10)
$$

with $\lambda = \pm 1$ representing the conduction ($\lambda = +1$) and valence ($\lambda = -1$) bands respectively, and,

$$
|f(\mathbf{p})|^2 = 1 + 4 \cos^2 \left( \frac{\sqrt{3} a p_x}{2} \right) \\
\qquad + 4 \cos \left( \frac{3 a p_x}{2} \right) \cos \left( \frac{\sqrt{3} a p_y}{2} \right). \quad (11)
$$

Inside the Brillouin Zone defined by the parallelogram
**b**₁ ⊗ **b**₂, this function is zero at the **K**<sub>±</sub> points defined as

$$
K_{\eta} = \frac{2\pi}{\sqrt{3}a} \begin{pmatrix} 1 - \frac{\eta}{3} \\ 0 \end{pmatrix}, \qquad (12)
$$

with $\eta = \pm 1$ the valley index. As the chemical potential $\mu$ is close to the crossing points between $\pi$ and $\pi^*$ bands at the $\mathbf{K}_\eta$ points, the dispersion of the bands can be approached as [52]

$$
\hat{H}_s(\mathbf{K}_\eta + \mathbf{k}) = \frac{3at}{2} (\eta k_1 \tau_1 + k_2 \tau_2), \quad (13)
$$
where $\tau_i$ is the $i^{th}$ Pauli matrix of the sublattice pseudo-spin for the A and B sites. From this expression the Fermi velocity is [50][52]

$$v_F = \frac{3at}{2\hbar} \approx \frac{c}{300}. \qquad (14)$$

Then, the dispersion band for each valley is approached as

$$\hat{H}_s(\mathbf{K}_\eta + \mathbf{k}) = \hat{H}_s^\eta(\mathbf{k}) = \hbar v_F (\eta k_1 \tau_1 + k_2 \tau_2). \quad (15)$$

After applying this small $(\mathbf{k} \cdot \mathbf{p})$ expansion, the electronic Hamiltonian can be approached as a family of 4 ($g_s = 2$ because of spin degeneration and $g_v = 2$ because of the 2 different valleys) 2D Dirac Hamiltonians placed in the continuum limit as [50]

$$H = \sum_{s, \eta = \pm} \int_{\mathbf{k}} \chi_s^{\eta, \dagger}(\mathbf{k}) \hat{H}_s^{\eta}(\mathbf{k}) \chi_s^{\eta}(\mathbf{k}). \quad (16)$$

This Hamiltonian represent a set of 4 equal Dirac cones, labelled by their valley $\eta$ and spin $s$. As a consequence, in addition to the discrete CPT symmetry, the Hamiltonian possesses a global continuous $U(4)$ symmetry that operates in the valley, sublattice and spin spaces [49].

Finally, we combine the bi-spinors of the same spin of the 2 valleys to form a Dirac four-spinor (note the exchange of sublattices of the $\eta = -1$ valley terms) [49]

$$\Psi_s(\mathbf{k}) = \begin{pmatrix} \chi_s(\mathbf{K}_+ + \mathbf{k}) \\ \tau_1 \chi_s(\mathbf{K}_- + \mathbf{k}) \end{pmatrix} = \begin{pmatrix} c_{A,s}(\mathbf{K}_+ + \mathbf{k}) \\ c_{B,s}(\mathbf{K}_+ + \mathbf{k}) \\ c_{B,s}(\mathbf{K}_- + \mathbf{k}) \\ c_{A,s}(\mathbf{K}_- + \mathbf{k}) \end{pmatrix}, \quad (17)$$

resulting into

$$H = \sum_{s=\pm} \int_{\mathbf{k}} \Psi_s^{\dagger}(\mathbf{k}) \hat{H}_s^D(\mathbf{k}) \Psi_s(\mathbf{k}), \quad (18)$$

$$\hat{H}_s^D(\mathbf{k}) = \hbar v_F \begin{pmatrix} 0 & k_1 - i k_2 & 0 & 0 \\ k_1 + i k_2 & 0 & 0 & 0 \\ 0 & 0 & 0 & -k_1 + i k_2 \\ 0 & 0 & -k_1 - i k_2 & 0 \end{pmatrix} = \hbar v_F (\alpha^1 k_1 + \alpha^2 k_2). \quad (19)$$

Here the $\alpha^\mu$ matrices are the Dirac matrices. For $i = \{1, 2, 3\}$, we have

$$\alpha^i = \tilde{\tau}_3 \otimes \tau_i = \begin{pmatrix} \tau_i & 0 \\ 0 & -\tau_i \end{pmatrix}, \quad (20)$$

where $\tilde{\tau}_i$ is the $i^{th}$ Pauli matrix of the valley pseudo-spin $\eta$ [49]. It will be useful to define $\alpha^0$ as

$$\alpha^0 = \begin{pmatrix} \tau_0 & 0 \\ 0 & \tau_0 \end{pmatrix} = \tilde{\tau}_0 \otimes \tau_0, \quad (21)$$

and we define the $\beta$ matrix as $\alpha^4$ in what follows

$$\alpha^4 = \beta = \tilde{\tau}_1 \otimes \tau_0 = \begin{pmatrix} 0 & \tau_0 \\ \tau_0 & 0 \end{pmatrix}. \quad (22)$$

From those definitions, we have that $\alpha^0$ is a 4 × 4 identity matrix, and the anticommutation relations

$$\{\alpha^i, \alpha^j\} = 2\delta^{ij}\alpha^0 \quad \forall i, j \in \{1, 2, 3, 4\}. \quad (23)$$

As a conclusion, we obtain 2 equivalent descriptions of the Hamiltonian of graphene, one in Eq. (16), as the sum of $g_s g_v = 4$ bi-spinors in the sub-lattice space, and another one in Eq. (18) the sum of $g_s = 2$ four-spinors in the sub-lattice-valley space.

### A. Action of graphene

To connect to the covariant QFT description of graphene [8][9][21][49], we write the full space-temporal second quantized action of graphene given in Eq. (18) as

$$\begin{aligned} S_0 &= \sum_{s=\pm} \int_{\mathbf{k}} \Psi_s^{\dagger}(\mathbf{k}) [\hbar\omega\alpha^0 - \hat{H}_s^D(\mathbf{k})] \Psi_s(\mathbf{k}) \\ &= \sum_{s=\pm} \int_{\mathbf{k}} \Psi_s^{\dagger}(\mathbf{k}) \hat{H}_{0,s}^D(\mathbf{k}) \Psi_s(\mathbf{k}), \end{aligned} \quad (24)$$

$$\hat{H}_{s,0}^{D}(\mathbf{k}) = \tilde{k}_{0}\alpha^{0} - \hbar v_{F} (\alpha^{1}k_{1} + \alpha^{2}k_{2}), \quad (25)$$

with $\tilde{k}_0 = \hbar\omega$ (see tab. I). This is the Dirac representation of the action of the $(2+1)D$ Dirac field. To write this action in a full covariant way by using the Weyl representation, we define the $\gamma$ matrices as $\gamma^\mu = \alpha^4\alpha^\mu$. With this prescription, we have $\gamma^0 = \alpha^4\alpha^0 = \alpha^4$, $\gamma^4 = \alpha^4\alpha^4 = \alpha^0$ and, for $i = \{1, 2, 3\}$

$$\gamma^i = \alpha^4 \alpha^i = (-i\tilde{\tau}_2) \otimes \tau_i = \begin{pmatrix} 0 & -\tau_i \\ \tau_i & 0 \end{pmatrix}. \quad (26)$$

The usual anti-commutation relations are fulfilled

$$\{\gamma^\mu, \gamma^\nu\} = 2g^{\mu\nu}\gamma^4 \quad \forall \{\mu, \nu\} \in \{0, 1, 2, 3\}, \quad (27)$$

where $g^{\mu\nu} = \text{diag}(+1, -1, -1, -1)$ is the metric tensor. The Dirac conjugated spinor is defined as $\tilde{\Psi}_s(\mathbf{k}) = \Psi_s^\dagger(\mathbf{k})\alpha^4 = \Psi_s^\dagger(\mathbf{k})\gamma^0$. Then the action is now represented as

$$S_0 = \sum_{s=\pm} \int_{\mathbf{k}} \bar{\Psi}_s(\mathbf{k}) \hat{H}_s^W(\mathbf{k}) \Psi_s(\mathbf{k}), \quad (28)$$

with

$$\hat{H}_{s,0}^{W}(\mathbf{k}) = \tilde{k}_{0}\gamma^{0} - \hbar v_{F} (\gamma^{1}k_{1} + \gamma^{2}k_{2}). \quad (29)$$

When this Hamiltonian is perturbed, depending on the breaking of the CPT discrete symmetries and on the generators of the $U(4)$ symmetry used, different kinds of gaps in the Dirac bands can be induced [49][51]. Here we will focus in two kind of non-topological mass gaps, that we will denote them as $m_z$ and $m$ in what follows. The full Hamiltonian become

$$\hat{H}_s^D(\mathbf{k}) = \tilde{k}_0\alpha^0 - (\alpha^1\tilde{k}_1 + \alpha^2\tilde{k}_2) - \alpha^3 m_z - \alpha^4 m, \quad (30)$$
where $\tau_i$ is the $i^{\text{th}}$ Pauli matrix of the sublattice pseudo-spin for the A and B sites. From this expression the Fermi velocity is [50][52]

$$v_F = \frac{3at}{2\hbar} \approx \frac{c}{300}. \qquad (14)$$

Then, the dispersion band for each valley is approached as

$$\hat{H}_s(\mathbf{K}_\eta + \mathbf{k}) = \hat{H}_s^\eta(\mathbf{k}) = \hbar v_F (\eta k_1 \tau_1 + k_2 \tau_2). \quad (15)$$

After applying this small $(\mathbf{k} \cdot \mathbf{p})$ expansion, the electronic Hamiltonian can be approached as a family of 4 ($g_s = 2$ because of spin degeneration and $g_v = 2$ because of the 2 different valleys) 2D Dirac Hamiltonians placed in the continuum limit as [50]

$$H = \sum_{s, \eta = \pm} \int_{\mathbf{k}} \chi_s^{\eta, \dagger}(\mathbf{k}) \hat{H}_s^{\eta}(\mathbf{k}) \chi_s^{\eta}(\mathbf{k}). \quad (16)$$

This Hamiltonian represent a set of 4 equal Dirac cones, labelled by their valley $\eta$ and spin $s$. As a consequence, in addition to the discrete CPT symmetry, the Hamiltonian possesses a global continuous $U(4)$ symmetry that operates in the valley, sublattice and spin spaces [49].

Finally, we combine the bi-spinors of the same spin of the 2 valleys to form a Dirac four-spinor (note the exchange of sublattices of the $\eta = -1$ valley terms) [49]

$$\Psi_s(\mathbf{k}) = \begin{pmatrix} \chi_s(\mathbf{K}_+ + \mathbf{k}) \\ \tau_1 \chi_s(\mathbf{K}_- + \mathbf{k}) \end{pmatrix} = \begin{pmatrix} c_{A,s}(\mathbf{K}_+ + \mathbf{k}) \\ c_{B,s}(\mathbf{K}_+ + \mathbf{k}) \\ c_{B,s}(\mathbf{K}_- + \mathbf{k}) \\ c_{A,s}(\mathbf{K}_- + \mathbf{k}) \end{pmatrix}, \quad (17)$$

resulting into

$$H = \sum_{s=\pm} \int_{\mathbf{k}} \Psi_s^{\dagger}(\mathbf{k}) \hat{H}_s^D(\mathbf{k}) \Psi_s(\mathbf{k}), \quad (18)$$

$$\hat{H}_s^D(\mathbf{k}) = \hbar v_F \begin{pmatrix} 0 & k_1 - i k_2 & 0 & 0 \\ k_1 + i k_2 & 0 & 0 & 0 \\ 0 & 0 & 0 & -k_1 + i k_2 \\ 0 & 0 & -k_1 - i k_2 & 0 \end{pmatrix} = \hbar v_F (\alpha^1 k_1 + \alpha^2 k_2). \quad (19)$$

Here the $\alpha^\mu$ matrices are the Dirac matrices. For $i = \{1, 2, 3\}$, we have

$$\alpha^i = \tilde{\tau}_3 \otimes \tau_i = \begin{pmatrix} \tau_i & 0 \\ 0 & -\tau_i \end{pmatrix}, \quad (20)$$

where $\tilde{\tau}_i$ is the $i^{\text{th}}$ Pauli matrix of the valley pseudo-spin $\eta$ [49]. It will be useful to define $\alpha^0$ as

$$\alpha^0 = \begin{pmatrix} \tau_0 & 0 \\ 0 & \tau_0 \end{pmatrix} = \tilde{\tau}_0 \otimes \tau_0, \quad (21)$$

and we define the $\beta$ matrix as $\alpha^4$ in what follows

$$\alpha^4 = \beta = \tilde{\tau}_1 \otimes \tau_0 = \begin{pmatrix} 0 & \tau_0 \\ \tau_0 & 0 \end{pmatrix}. \quad (22)$$

From those definitions, we have that $\alpha^0$ is a 4 × 4 identity matrix, and the anticommutation relations

$$\{\alpha^i, \alpha^j\} = 2\delta^{ij}\alpha^0 \quad \forall i, j \in \{1, 2, 3, 4\}. \quad (23)$$

As a conclusion, we obtain 2 equivalent descriptions of the Hamiltonian of graphene, one in Eq. (16), as the sum of $g_s g_v = 4$ bi-spinors in the sub-lattice space, and another one in Eq. (18) the sum of $g_s = 2$ four-spinors in the sub-lattice-valley space.

### A. Action of graphene

To connect to the covariant QFT description of graphene [8][9][21][49], we write the full space-temporal second quantized action of graphene given in Eq. (18) as

$$\begin{aligned} S_0 &= \sum_{s=\pm} \int_{\mathbf{k}} \Psi_s^{\dagger}(\mathbf{k}) [\hbar\omega\alpha^0 - \hat{H}_s^D(\mathbf{k})] \Psi_s(\mathbf{k}) \\ &= \sum_{s=\pm} \int_{\mathbf{k}} \Psi_s^{\dagger}(\mathbf{k}) \hat{H}_{0,s}^D(\mathbf{k}) \Psi_s(\mathbf{k}), \end{aligned} \quad (24)$$

$$\hat{H}_{s,0}^{D}(\mathbf{k}) = \tilde{k}_{0}\alpha^{0} - \hbar v_{F} (\alpha^{1}k_{1} + \alpha^{2}k_{2}), \quad (25)$$

with $\tilde{k}_0 = \hbar\omega$ (see tab. I). This is the Dirac representation of the action of the $(2+1)D$ Dirac field. To write this action in a full covariant way by using the Weyl representation, we define the $\gamma$ matrices as $\gamma^\mu = \alpha^4\alpha^\mu$. With this prescription, we have $\gamma^0 = \alpha^4\alpha^0 = \alpha^4$, $\gamma^4 = \alpha^4\alpha^4 = \alpha^0$ and, for $i = \{1, 2, 3\}$

$$\gamma^i = \alpha^4 \alpha^i = (-i\tilde{\tau}_2) \otimes \tau_i = \begin{pmatrix} 0 & -\tau_i \\ \tau_i & 0 \end{pmatrix}. \quad (26)$$

The usual anti-commutation relations are fulfilled

$$\{\gamma^\mu, \gamma^\nu\} = 2g^{\mu\nu}\gamma^4 \quad \forall \{\mu, \nu\} \in \{0, 1, 2, 3\}, \quad (27)$$

where $g^{\mu\nu} = \text{diag}(+1, -1, -1, -1)$ is the metric tensor. The Dirac conjugated spinor is defined as $\tilde{\Psi}_s(\mathbf{k}) = \Psi_s^\dagger(\mathbf{k})\alpha^4 = \Psi_s^\dagger(\mathbf{k})\gamma^0$. Then the action is now represented as

$$S_0 = \sum_{s=\pm} \int_{\mathbf{k}} \bar{\Psi}_s(\mathbf{k}) \hat{H}_s^W(\mathbf{k}) \Psi_s(\mathbf{k}), \quad (28)$$

with

$$\hat{H}_{s,0}^{W}(\mathbf{k}) = \tilde{k}_{0}\gamma^{0} - \hbar v_{F} (\gamma^{1}k_{1} + \gamma^{2}k_{2}). \quad (29)$$

When this Hamiltonian is perturbed, depending on the breaking of the CPT discrete symmetries and on the generators of the $U(4)$ symmetry used, different kinds of gaps in the Dirac bands can be induced [49][51]. Here we will focus in two kind of non-topological mass gaps, that we will denote them as $m_z$ and $m$ in what follows. The full Hamiltonian become

$$\hat{H}_s^D(\mathbf{k}) = \tilde{k}_0\alpha^0 - (\alpha^1\tilde{k}_1 + \alpha^2\tilde{k}_2) - \alpha^3 m_z - \alpha^4 m, \quad (30)$$
$$ \hat{\mathcal{H}}_s^W(\mathbf{k}) = \tilde{k}_0\gamma^0 - (\gamma^1\tilde{k}_1 + \gamma^2\tilde{k}_2) - \gamma^3 m_z - \gamma^4 m, \quad (31) $$

where $m_z$ is a hopping term between fermions of the same valley [50][51], while $m$ couples quasi-particles of different valleys, being it a non-local interaction, but it is a possible result for the symmetry breaking interaction over graphene [49].

## B. Macrocanonical ensemble

As we are working in the macrocanonical ensemble, we add a term to the action of our field

$$ S_{\mu} = \sum_{s=\pm} \int_{\mathbf{k}} \mu \hat{\mathcal{N}}, \quad (32) $$

where $\mu$ is the chemical potential and $\hat{\mathcal{N}}$ is the number of particles operator. This term can be written in each representation as

$$ \begin{aligned} S_{\mu} &= \mu \sum_{s=\pm} \int_{\mathbf{k}} \Psi_s^{\dagger}(\mathbf{k}) \alpha^0 \Psi_s(\mathbf{k}) \\ &= \mu \sum_{s=\pm} \int_{\mathbf{k}} \bar{\Psi}_s(\mathbf{k}) \gamma^0 \Psi_s(\mathbf{k}). \end{aligned} \quad (33) $$

In order to compare between the different models of the polarization operator, we are going to use three different macrocanonical hamiltonians for graphene, the bi-spinor expression from Eq. (16), where

$$ S_1 = \sum_{s, \eta=\pm} \int_{\mathbf{k}} \chi_s^{\eta, \dagger}(\mathbf{k}) \left[ \tau_0(\hbar\omega + \mu) - \hat{H}_s^{\eta}(\mathbf{k}) \right] \chi_s^{\eta}(\mathbf{k}) \quad (34) $$

$$ \hat{H}_s^{\eta}(\mathbf{k}) = \hbar v_F [\eta k_1 \tau_1 + k_2 \tau_2] + \tau_3 \Delta_s^{\eta}, \quad (35) $$

with $m_z = \Delta_s^{\eta}$, the Dirac form of the Dirac Hamiltonian (using Eq. (30)) as a bridge between the 2 formalisms

$$ S_D = \sum_{s=\pm} \int_{\mathbf{k}} \Psi_s^{\dagger}(\mathbf{k}) \left[ \hat{\mathcal{H}}_s^D(\mathbf{k}) - \alpha^0 \mu \right] \Psi_s(\mathbf{k}), \quad (36) $$

and the covariant expression of the Dirac Hamiltonian (using Eq. (31)) as

$$ S_W = \sum_{s=\pm} \int_{\mathbf{k}} \bar{\Psi}_s^{\dagger}(\mathbf{k}) \left[ \hat{\mathcal{H}}_s^W(\mathbf{k}) - \gamma^0 \mu \right] \Psi_s(\mathbf{k}). \quad (37) $$

## C. Effect of interactions

Electronic quasiparticles are subject to different possible interactions: phonons, scattering centers, the unavoidable Coulomb interaction, external fields, illumination, decoration (impurities) and so on [4][53][54][2][29][55]. When the effects of interactions is

taking into account the dynamics of the electronic quasiparticles, the Hamiltonian is modified by the causal self-energy $\Sigma(\omega, \mathbf{q}) = \Sigma_R(\omega, \mathbf{q}) + i\Sigma_I(\omega, \mathbf{q})$ [10][29][55], then the electronic spectrum is modified by the addition of a real self-energy $\Sigma_R(\omega, \mathbf{q})$ (whose effect is considered here small and absorbed into the phenomenological constants of the Hamiltonian) and an always non-negative imaginary part $\Sigma_I(\omega, \mathbf{q})$ (to respect thermodynamics), which leads to a frequency-dependent finite dissipation of the electronic quasiparticles [54][4][2]. Here we argue that the electronic dissipation is a small non-zero quantity and we assume that its effect on the electric conductivity can be well approached with the finite lifetime approximation by a constant imaginary dissipation rate $\Gamma = \tau^{-1}$ [4][56][55][57]. Taking into account that the measured electrical conductivity of graphene is a high but finite quantity ($\sigma = 96 \times 10^6 \text{ Sm}^{-1} \text{ in[58]}$), and that the dissipation time has been estimated to be on the order of $\tau \sim 6 \times 10^{-13} \text{ s[4][53]}$, we will take this quantity in our study. Of course, in situations where the effect of interactions in graphene is of paramount relevance (beyond its non-zero nature), like in the study of the electron-phonon interaction [29], the universal DC conductivity of graphene when $\mu = 0$ [10][11][12][59][60], our interactions-naive phenomenological approach would not be enough and a more detailed study of the effect of interactions will be necessary.

## D. Green function

For a general linear hamiltonian $\hat{\mathcal{H}}(\mathbf{r})$, we have 2 different expressions of the same Green function. Defining $\tilde{k}_0 = \hbar\omega$, the Green function fulfils

$$ \hat{\mathcal{H}}(\mathbf{r}) \mathcal{G}_0(\mathbf{r}, \mathbf{s}) = \delta(\mathbf{r} - \mathbf{s}), \quad (38) $$

then, in momentum space we have $\hat{\mathcal{H}}(\mathbf{k})\mathcal{G}_0(\mathbf{k}) = 1$ and, therefore $\mathcal{G}_0(\mathbf{k}) = \hat{\mathcal{H}}^{-1}(\mathbf{k})$. For each one of the four-spinor Hamiltonians, its inverse operator is

$$ \mathcal{G}_0^D(\mathbf{k}) = \hat{\mathcal{H}}_D^{-1}(\mathbf{k}) = \frac{1}{\alpha^\mu K_\mu - \alpha^4 m} = \frac{\alpha^\mu K_\mu + \alpha^4 m}{K^\mu K_\mu - m^2}, \quad (39) $$

$$ \mathcal{G}_0^W(\mathbf{k}) = \hat{\mathcal{H}}_W^{-1}(\mathbf{k}) = \frac{1}{\gamma^\mu K_\mu - \gamma^4 m} = \frac{\gamma^\mu K_\mu + \gamma^4 m}{K^\mu K_\mu - m^2}, \quad (40) $$

where we define $K_\mu = (\hbar\omega + \mu, \hbar v_F k_1, \hbar v_F k_2, m_z)$ (see tab. I), the Einstein summation convention is assumed and

$$ K^{\mu} K_{\mu} - m^{2} = \prod_{\lambda=\pm} [\hbar\omega + \mu - \epsilon_{\mathbf{k}}^{\lambda}] = \prod_{\lambda=\pm} [\hbar\omega - \xi_{\mathbf{k}}^{\lambda}], \quad (41) $$

where

$$ \epsilon_{\mathbf{k}}^{\lambda} = \lambda \sqrt{(\hbar v_F k)^2 + m_z^2 + m^2}. \quad (42) $$

with $k = \sqrt{k_1^2 + k_2^2}$ (see tab. I). There is another equivalent expression for a general linear hamiltonian of the
$$ \hat{\mathcal{H}}_s^W(\mathbf{k}) = \tilde{k}_0\gamma^0 - (\gamma^1\tilde{k}_1 + \gamma^2\tilde{k}_2) - \gamma^3 m_z - \gamma^4 m, \quad (31) $$

where $m_z$ is a hopping term between fermions of the same valley [50][51], while $m$ couples quasi-particles of different valleys, being it a non-local interaction, but it is a possible result for the symmetry breaking interaction over graphene [49].

## B. Macrocanonical ensemble

As we are working in the macrocanonical ensemble, we add a term to the action of our field

$$ S_{\mu} = \sum_{s=\pm} \int_{\mathbf{k}} \mu \hat{\mathcal{N}}, \quad (32) $$

where $\mu$ is the chemical potential and $\hat{\mathcal{N}}$ is the number of particles operator. This term can be written in each representation as

$$ \begin{aligned} S_{\mu} &= \mu \sum_{s=\pm} \int_{\mathbf{k}} \Psi_s^{\dagger}(\mathbf{k}) \alpha^0 \Psi_s(\mathbf{k}) \\ &= \mu \sum_{s=\pm} \int_{\mathbf{k}} \bar{\Psi}_s(\mathbf{k}) \gamma^0 \Psi_s(\mathbf{k}). \end{aligned} \quad (33) $$

In order to compare between the different models of the polarization operator, we are going to use three different macrocanonical hamiltonians for graphene, the bi-spinor expression from Eq. (16), where

$$ S_1 = \sum_{s, \eta=\pm} \int_{\mathbf{k}} \chi_s^{\eta, \dagger}(\mathbf{k}) [\tau_0(\hbar\omega + \mu) - \hat{H}_s^{\eta}(\mathbf{k})] \chi_s^{\eta}(\mathbf{k}) (34) $$

$$ \hat{H}_s^{\eta}(\mathbf{k}) = \hbar v_F [\eta k_1 \tau_1 + k_2 \tau_2] + \tau_3 \Delta_s^{\eta}, \quad (35) $$

with $m_z = \Delta_s^{\eta}$, the Dirac form of the Dirac Hamiltonian (using Eq. (30)) as a bridge between the 2 formalisms

$$ S_D = \sum_{s=\pm} \int_{\mathbf{k}} \Psi_s^{\dagger}(\mathbf{k}) [\hat{\mathcal{H}}_s^D(\mathbf{k}) - \alpha^0 \mu] \Psi_s(\mathbf{k}), \quad (36) $$

and the covariant expression of the Dirac Hamiltonian (using Eq. (31)) as

$$ S_W = \sum_{s=\pm} \int_{\mathbf{k}} \bar{\Psi}_s^{\dagger}(\mathbf{k}) [\hat{\mathcal{H}}_s^W(\mathbf{k}) - \gamma^0 \mu] \Psi_s(\mathbf{k}). \quad (37) $$

## C. Effect of interactions

Electronic quasiparticles are subject to different possible interactions: phonons, scattering centers, the unavoidable Coulomb interaction, external fields, illumination, decoration (impurities) and so on [4][53][54][2][29][55]. When the effects of interactions is

taking into account the dynamics of the electronic quasiparticles, the Hamiltonian is modified by the causal self-energy $\Sigma(\omega, \mathbf{q}) = \Sigma_R(\omega, \mathbf{q}) + i\Sigma_I(\omega, \mathbf{q})$ [10][29][55], then the electronic spectrum is modified by the addition of a real self-energy $\Sigma_R(\omega, \mathbf{q})$ (whose effect is considered here small and absorbed into the phenomenological constants of the Hamiltonian) and an always non-negative imaginary part $\Sigma_I(\omega, \mathbf{q})$ (to respect thermodynamics), which leads to a frequency-dependent finite dissipation of the electronic quasiparticles [54][4][2]. Here we argue that the electronic dissipation is a small non-zero quantity and we assume that its effect on the electric conductivity can be well approached with the finite lifetime approximation by a constant imaginary dissipation rate $\Gamma = \tau^{-1}$ [4][56][55][57]. Taking into account that the measured electrical conductivity of graphene is a high but finite quantity ($\sigma = 96 \times 10^6 \text{ Sm}^{-1} \text{ in[58]}$), and that the dissipation time has been estimated to be on the order of $\tau \sim 6 \times 10^{-13} \text{ s[4][53]}$, we will take this quantity in our study. Of course, in situations where the effect of interactions in graphene is of paramount relevance (beyond its non-zero nature), like in the study of the electron-phonon interaction [29], the universal DC conductivity of graphene when $\mu = 0$ [10][11][12][59][60], our interactions-naive phenomenological approach would not be enough and a more detailed study of the effect of interactions will be necessary.

## D. Green function

For a general linear hamiltonian $\hat{\mathcal{H}}(\mathbf{r})$, we have 2 different expressions of the same Green function. Defining $\tilde{k}_0 = \hbar\omega$, the Green function fulfils

$$ \hat{\mathcal{H}}(\mathbf{r}) \mathcal{G}_0(\mathbf{r}, \mathbf{s}) = \delta(\mathbf{r} - \mathbf{s}), \quad (38) $$

then, in momentum space we have $\hat{\mathcal{H}}(\mathbf{k})\mathcal{G}_0(\mathbf{k}) = 1$ and, therefore $\mathcal{G}_0(\mathbf{k}) = \hat{\mathcal{H}}^{-1}(\mathbf{k})$. For each one of the four-spinor Hamiltonians, its inverse operator is

$$ \mathcal{G}_0^D(\mathbf{k}) = \hat{\mathcal{H}}_D^{-1}(\mathbf{k}) = \frac{1}{\alpha^\mu K_\mu - \alpha^4 m} = \frac{\alpha^\mu K_\mu + \alpha^4 m}{K^\mu K_\mu - m^2}, \quad (39) $$

$$ \mathcal{G}_0^W(\mathbf{k}) = \hat{\mathcal{H}}_W^{-1}(\mathbf{k}) = \frac{1}{\gamma^\mu K_\mu - \gamma^4 m} = \frac{\gamma^\mu K_\mu + \gamma^4 m}{K^\mu K_\mu - m^2}, \quad (40) $$

where we define $K_\mu = (\hbar\omega + \mu, \hbar v_F k_1, \hbar v_F k_2, m_z)$ (see tab. I), the Einstein summation convention is assumed and

$$ K^\mu K_\mu - m^2 = \prod_{\lambda=\pm} [\hbar\omega + \mu - \epsilon_\lambda^\lambda] = \prod_{\lambda=\pm} [\hbar\omega - \xi_\lambda^\lambda], \quad (41) $$

where

$$ \epsilon_\lambda^\lambda = \lambda \sqrt{(\hbar v_F k)^2 + m_z^2 + m^2}. \quad (42) $$

with $k = \sqrt{k_1^2 + k_2^2}$ (see tab. I). There is another equivalent expression for a general linear hamiltonian of the
form $\hbar\omega\psi = \tilde{k}_0\psi = \hat{H}\psi$ in terms of the eigenvalues and eigenfunctions of the Hamiltonian, starting again from the equation of the Green function $(\tilde{k}_0 - \hat{H}) G_0(r, s) = \delta(r - s)$, from the eigenproblem

$$ \hat{H} |u_k^\lambda\rangle = (\epsilon_k^\lambda - \mu) |u_k^\lambda\rangle = \xi_k^\lambda |u_k^\lambda\rangle, \quad (43) $$

note that the chemical potential $\mu$ has been absorbed in $\hat{H}$. The macrocanonical Green function in momentum space is

$$ G_0(k) = \sum_{\lambda} \frac{|u_k^{\lambda}\rangle\langle u_k^{\lambda}|}{\tilde{k}_0 - \xi_k^{\lambda}}, \quad (44) $$

which is an eigenvalue expansion of the Green function. In our case, the Hamiltonian is diagonal in spin, therefore, the Green functions are multiplied by $\delta_{ss'}$.

### E. Presence of an electromagnetic field

We introduce the coupling of the electronic quasiparticles of the lattice to the electromagnetic field via the Peierls substitution [49][61]

$$ H = \sum_{n \in \mathbb{Z}^d} \sum_{\lambda, \lambda'} \sum_{j=1}^{|\delta_j|<\delta} \hat{c}_{\lambda, \mathbf{R}_n}^{\dagger} t_{\lambda, \lambda'}(\delta_j) e^{i q \mathbf{A}_{\mathbf{R}_n} \cdot \delta_j} \hat{c}_{\lambda', \mathbf{R}_n + \delta_j} \quad (45) $$

where we have approached

$$ \exp \left( i q \int_{\mathbf{R}_n}^{\mathbf{R}_n + \delta_j} \mathbf{A}(\mathbf{r}) \cdot d\mathbf{r} \right) \approx e^{i q \mathbf{A}_{\mathbf{R}_n} \cdot \delta_j}. \quad (46) $$

being $q$ the electric charge of the quasi-particle described by the Hamiltonian, for electronic excitations we have $q = -e$. At linear order in $\mathbf{A}_k$, the Hamiltonian in reciprocal space is

$$ \begin{aligned} H &= \int_{\mathbf{q}} \int_{\mathbf{p}} \sum_{\lambda, \lambda'} \hat{c}_{\mathbf{p}+\mathbf{q}}^{\lambda, \dagger} \left[ \sum_{j=1}^{|\delta_j|<\delta} t_{\lambda, \lambda'}(\delta_j) e^{i(\mathbf{p}+\mathbf{q}\mathbf{A}_\mathbf{q})\cdot\delta_j} \right] \hat{c}_{\mathbf{p}}^{\lambda'} \\ &= \int_{\mathbf{q}} \int_{\mathbf{p}} \sum_{\lambda, \lambda'} \hat{c}_{\mathbf{p}+\mathbf{q}}^{\lambda, \dagger} \hat{H}^{\lambda, \lambda'} (\mathbf{p}+\mathbf{q}\mathbf{A}_\mathbf{q}) \hat{c}_{\mathbf{p}}^{\lambda'}. \end{aligned} \quad (47) $$

It is clear that, at first order, the inclusion of the Peiers substitution leads to a minimal coupling of the momentum [49]. To study the electric conductivity, we need an expression for the current, understood as the conjugated force of the potential vector. Then, expanding the Hamiltonian at linear order in $\mathbf{A}$, we obtain at linear order

$$ \begin{aligned} H &= H_0 + \int_{\mathbf{q}} J_{\mu, -\mathbf{q}} A_{\mathbf{q}}^{\mu} = H_0 + \int_{\mathbf{q}} J_{\mu, \mathbf{q}}^{*} A_{\mathbf{q}}^{\mu} \\ &\approx H_0 + \int_{\mathbf{q}} \frac{\delta \hat{H}}{\delta A_{\mathbf{q}}^{\mu}} A_{\mathbf{q}}^{\mu}. \end{aligned} \quad (48) $$

Therefore, the second quantized current is defined as

$$ J_{\mu, \mathbf{q}}^{*} = \frac{\delta \hat{H}}{\delta A_{\mathbf{q}}^{\mu}} = \int_{\mathbf{p}} \sum_{\lambda, \lambda'} \hat{c}_{\mathbf{p}+\mathbf{q}}^{\lambda, \dagger} \hat{J}_{\mu}^{\lambda, \lambda'}(\mathbf{p}) \hat{c}_{\mathbf{p}}^{\lambda'} \quad (49) $$

with the current operator given as

$$ \begin{aligned} \hat{J}_{\mu}^{\lambda, \lambda'}(\mathbf{p}) &= \left. \frac{\partial \hat{H}^{\lambda, \lambda'}(\mathbf{p} + q\mathbf{A}_q)}{\partial A_{\mathbf{q}}^{\mu}} \right|_{\mathbf{A}_q \to \mathbf{0}} = q \frac{\partial \hat{H}^{\lambda, \lambda'}(\mathbf{p})}{\partial p^{\mu}} \\ &= \frac{q}{\hbar} \frac{\partial \hat{H}^{\lambda, \lambda'}(\mathbf{k})}{\partial k_{\mu}} = q \hat{v}_{\mu}, \end{aligned} \quad (50) $$

where we have used $\mathbf{p} = \hbar\mathbf{k}$ and $\hat{v}_\mu$ is the velocity operator of electronic quasiparticles. For electronic excitations, we have $q = -e$.

### III. EM RESPONSE OF GRAPHENE: THE QFT-B MODEL CORRECTLY REGULARIZED IS IDENTICAL TO THE KUBO MODEL

#### A. Constitutive relation of the electric conductivity

Starting from the microscopic Ohm law, the electronic current is obtained as

$$ J_{\mu}(\omega, \mathbf{q}) = \sigma_{\mu\nu}(\omega, \mathbf{q}) E^{\nu}(\omega, \mathbf{q}), \quad (51) $$

where $\sigma_{\mu\nu}(\mathbf{q})$ is the conductivity operator, $\mathbf{q} = (\tilde{q}_0, \tilde{q})$ (see tab. I) and $E^\mu(\omega, \mathbf{q}) = -i\omega A^\mu(\omega, \mathbf{q})$ is the electric field in Coulomb Gauge. Note that we are working in Coulomb Gauge, but we can restore Gauge invariance at the end of the calculations by writing all final results in terms of explicitly Gauge-invariant quantities. In appendix D, by using the Kubo formula [5] on the microscopic Ohm law, the Luttinger formula is obtained [62][63]

$$ \sigma_{\mu\nu}(\mathbf{q}) = \frac{\Pi_{\mu\nu}(\omega, \mathbf{q}) - \lim_{\omega \to 0^+} \Pi_{\mu\nu}(\omega, \mathbf{q})}{-i\omega} = \frac{\tilde{\Pi}_{\mu\nu}(\omega, \mathbf{q})}{-i\omega} \quad (52) $$

where the polarization operator is defined as

$$ \Pi_{\mu\nu}(\mathbf{q}) = \mathrm{Tr}(J_{\mu}(\mathbf{q})J_{\nu}^{*}(\mathbf{q})), \quad (53) $$

and

$$ \tilde{\Pi}^{\mu\nu}(\mathbf{q}) = \Pi^{\mu\nu}(\mathbf{q}) - \lim_{\omega \to 0} \Pi^{\mu\nu}(\omega, \mathbf{q}). \quad (54) $$

In appendix D, Eq. (D16), we also derive

$$ J_{\mu}(\omega, \mathbf{q}) = \left[ \Pi_{\mu\nu}(\omega, \mathbf{q}) - \lim_{\omega \to 0} \Pi_{\mu\nu}(\omega, \mathbf{q}) \right] A^{\nu}(\omega, \mathbf{q}). \quad (55) $$

Note that this subtraction, which we derived here simply from Ohm law and time causality, naturally implies that there cannot exist electronic current when $E^\nu = 0$, i.e., when always $J_\mu(\omega = 0, \mathbf{q}) = 0$ even when $A_\mu(\omega = 0, \mathbf{q}) \neq 0$ irrespectively of the functional and
form $\hbar\omega\psi = \tilde{k}_0\psi = \hat{H}\psi$ in terms of the eigenvalues and eigenfunctions of the Hamiltonian, starting again from the equation of the Green function $(\tilde{k}_0 - \hat{H}) G_0(r, s) = \delta(r - s)$, from the eigenproblem

$$ \hat{H} |u_k^\lambda\rangle = (\epsilon_k^\lambda - \mu) |u_k^\lambda\rangle = \xi_k^\lambda |u_k^\lambda\rangle, \quad (43) $$

note that the chemical potential $\mu$ has been absorbed in $\hat{H}$. The macrocanonical Green function in momentum space is

$$ G_0(k) = \sum_{\lambda} \frac{|u_k^{\lambda}\rangle\langle u_k^{\lambda}|}{\tilde{k}_0 - \xi_k^{\lambda}}, \quad (44) $$

which is an eigenvalue expansion of the Green function. In our case, the Hamiltonian is diagonal in spin, therefore, the Green functions are multiplied by $\delta_{ss'}$.

### E. Presence of an electromagnetic field

We introduce the coupling of the electronic quasiparticles of the lattice to the electromagnetic field via the Peierls substitution [49][61]

$$ H = \sum_{n \in \mathbb{Z}^d} \sum_{\lambda, \lambda'} \sum_{j=1}^{|\delta_j|<\delta} \hat{c}_{\lambda, \mathbf{R}_n}^{\dagger} t_{\lambda, \lambda'}(\delta_j) e^{i q \mathbf{A}_{\mathbf{R}_n} \cdot \delta_j} \hat{c}_{\lambda', \mathbf{R}_n + \delta_j} \quad (45) $$

where we have approached

$$ \exp \left( i q \int_{\mathbf{R}_n}^{\mathbf{R}_n + \delta_j} \mathbf{A}(\mathbf{r}) \cdot d\mathbf{r} \right) \approx e^{i q \mathbf{A}_{\mathbf{R}_n} \cdot \delta_j}. \quad (46) $$

being $q$ the electric charge of the quasi-particle described by the Hamiltonian, for electronic excitations we have $q = -e$. At linear order in $\mathbf{A}_k$, the Hamiltonian in reciprocal space is

$$ \begin{aligned} H &= \int_q \int_p \sum_{\lambda, \lambda'} \hat{c}_{p+q}^{\lambda, \dagger} \left[ \sum_{j=1}^{|\delta_j|<\delta} t_{\lambda, \lambda'}(\delta_j) e^{i(\mathbf{p}+q\mathbf{A}_q)\cdot\delta_j} \right] \hat{c}_p^{\lambda'} \\ &= \int_q \int_p \sum_{\lambda, \lambda'} \hat{c}_{p+q}^{\lambda, \dagger} \hat{H}^{\lambda, \lambda'} (\mathbf{p}+q\mathbf{A}_q) \hat{c}_p^{\lambda'}. \end{aligned} \quad (47) $$

It is clear that, at first order, the inclusion of the Peiers substitution leads to a minimal coupling of the momentum [49]. To study the electric conductivity, we need an expression for the current, understood as the conjugated force of the potential vector. Then, expanding the Hamiltonian at linear order in $\mathbf{A}$, we obtain at linear order

$$ \begin{aligned} H &= H_0 + \int_q J_{\mu, -q} A_q^\mu = H_0 + \int_q J_{\mu, q}^* A_q^\mu \\ &\approx H_0 + \int_q \frac{\delta \hat{H}}{\delta A_q^\mu} A_q^\mu. \end{aligned} \quad (48) $$

Therefore, the second quantized current is defined as

$$ J_{\mu,q}^* = \frac{\delta \hat{H}}{\delta A_q^\mu} = \int_p \sum_{\lambda, \lambda'} \hat{c}_{p+q}^{\lambda, \dagger} \hat{J}_{\mu, \lambda'}^{\lambda, \lambda'}(\mathbf{p}) \hat{c}_p^{\lambda'} \quad (49) $$

with the current operator given as

$$ \begin{aligned} \hat{J}_{\mu}^{\lambda, \lambda'}(\mathbf{p}) &= \left. \frac{\partial \hat{H}^{\lambda, \lambda'}(\mathbf{p} + q\mathbf{A}_q)}{\partial A_q^{\mu}} \right|_{\mathbf{A}_q \to 0} = q \frac{\partial \hat{H}^{\lambda, \lambda'}(\mathbf{p})}{\partial p^{\mu}} \\ &= \frac{q}{\hbar} \frac{\partial \hat{H}^{\lambda, \lambda'}(\mathbf{k})}{\partial k_{\mu}} = q \hat{v}_{\mu}, \end{aligned} \quad (50) $$

where we have used $\mathbf{p} = \hbar\mathbf{k}$ and $\hat{v}_\mu$ is the velocity operator of electronic quasiparticles. For electronic excitations, we have $q = -e$.

### III. EM RESPONSE OF GRAPHENE: THE QFT-B MODEL CORRECTLY REGULARIZED IS IDENTICAL TO THE KUBO MODEL

#### A. Constitutive relation of the electric conductivity

Starting from the microscopic Ohm law, the electronic current is obtained as

$$ J_{\mu}(\omega, \mathbf{q}) = \sigma_{\mu\nu}(\omega, \mathbf{q}) E^{\nu}(\omega, \mathbf{q}), \quad (51) $$

where $\sigma_{\mu\nu}(q)$ is the conductivity operator, $q = (\tilde{q}_0, \tilde{q})$ (see tab. I) and $E^\mu(\omega, \mathbf{q}) = -i\omega A^\mu(\omega, \mathbf{q})$ is the electric field in Coulomb Gauge. Note that we are working in Coulomb Gauge, but we can restore Gauge invariance at the end of the calculations by writing all final results in terms of explicitly Gauge-invariant quantities. In appendix D, by using the Kubo formula [5] on the microscopic Ohm law, the Luttinger formula is obtained [62][63]

$$ \sigma_{\mu\nu}(q) = \frac{\Pi_{\mu\nu}(\omega, q) - \lim_{\omega \to 0^+} \Pi_{\mu\nu}(\omega, q)}{-i\omega} = \frac{\tilde{\Pi}_{\mu\nu}(\omega, q)}{-i\omega} \quad (52) $$

where the polarization operator is defined as

$$ \Pi_{\mu\nu}(q) = \mathrm{Tr}(J_{\mu}(q)J_{\nu}^{*}(q)), \quad (53) $$

and

$$ \tilde{\Pi}^{\mu\nu}(q) = \Pi^{\mu\nu}(q) - \lim_{\omega \to 0} \Pi^{\mu\nu}(q). \quad (54) $$

In appendix D, Eq. (D16), we also derive

$$ J_{\mu}(\omega, \mathbf{q}) = \left[ \Pi_{\mu\nu}(\omega, \mathbf{q}) - \lim_{\omega \to 0} \Pi_{\mu\nu}(\omega, \mathbf{q}) \right] A^{\nu}(\omega, \mathbf{q}). \quad (55) $$

Note that this subtraction, which we derived here simply from Ohm law and time causality, naturally implies that there cannot exist electronic current when $E^\nu = 0$, i.e., when always $J_\mu(\omega = 0, \mathbf{q}) = 0$ even when $A_\mu(\omega = 0, \mathbf{q}) \neq 0$ irrespectively of the functional and
tensorial form of $\Pi_{\mu\nu}(\omega, q)$. This is a strong physical requirement, already discussed by V. Abrikosov in 1963 [64] and more recently in [28].

However, sometimes, a different transport relation is proposed [8][9][21][33][42][43][46][65]

$$J_{\mu}(\omega, q) = \Pi_{\mu\nu}(\omega, q) A^{\nu}(\omega, q), \quad (56)$$

This is a natural choice for the transport coefficient provided in Eq. (48). However, as written, this linear relation could possibly lead to non zero electric currents even when the electric field is zero and to spurious dissipationless currents. This pathological result is present in the QFT-b model as is developed and used in literature, and that will be discussed in detail in section Sect. VII.

## B. Polarization Operator

The Kubo formula for the polarization operator is

$$\Pi_{\mu\nu}(q) = \mathrm{Tr} (J_{\mu}(q) J_{\nu}^{*}(q)), \quad (57)$$

which can be reduced to the bubble Feynman diagram as [66][67]

$$\Pi_{\mu\nu}(q) = i \int_k \mathrm{Tr} \left( G_0^{\lambda}(k) \hat{J}_{\mu}^{\lambda, \lambda'}(k) G_0^{\lambda'}(k+q) \hat{J}_{\nu}^{\lambda', \lambda}(k+q) \right) (58)$$

Here we use the definition $\int_k = \int_{-\infty}^{\infty} \frac{d\tilde{k}_0}{2\pi} \int_k$ (see tab. I), $G_0^\lambda(k) = \langle \hat{c}_k^\lambda \hat{c}_k^{\lambda, \dagger} \rangle_0$ is the Green function of the unperturbed Hamiltonian of the system, and $\hat{J}_\mu^{\lambda, \lambda'}(q)$ is the current operator of electronic quasiparticles, defined in Eq. (50). For a general linear hamiltonian of the form $\hbar\omega\psi = \tilde{k}_0\psi = \hat{H}\psi$, the Green function fulfils $(\tilde{k}_0\delta^{\lambda, \lambda'} - \hat{H}^{\lambda, \lambda'}) G_0^\lambda(r, s) = \delta^{\lambda, \lambda'}\delta(r - s)$, and we have the eigenproblem

$$\hat{H}^{\lambda, \lambda'} |u^{\lambda}\rangle = (\epsilon^{\lambda'} - \mu) |u^{\lambda'}\rangle = \xi^{\lambda'} |u^{\lambda'}\rangle, \quad (59)$$

where we have removed the trace operator because it is only applied in the $k$-space and in the $(\lambda, \lambda')$-bands space.

In order to connect with the conductivity obtained from the Kubo formula, using that the Kubo formula for the Ohm law imposes that $\lim_{\omega \to 0} \Pi^{\mu\nu}(\omega) A_\nu(\omega) = 0$ for constant static $A_\nu(\omega)$ [20] (except for superconductors [23][22]), we have to remove the effect of the $\omega \to 0$ limit

note that the chemical potential $\mu$ has been absorbed in $\hat{H}$. The macrocanonical Green function in momentum space is

$$G_0^{\lambda}(k) = \sum_{\lambda} \frac{|u_{\hat{k}}^{\lambda}\rangle\langle u_{\hat{k}}^{\lambda}|}{\tilde{k}_0 - \xi_{\hat{k}}^{\lambda}}. \quad (60)$$

The use of a different (space-time covariant) form of the Green function for the Dirac Hamiltonian leads to different representations of the same result discussed here [9]. After introducing the Green function into the polarization operator, using that $\mathrm{Tr}(ABC) = \mathrm{Tr}(BCA)$, and applying the Matsubara formalism to carry out the $k_0$ integral for the fermionic case by using $\tilde{k}_0 = \hbar\omega_n = \frac{2\pi}{\beta}(n + \frac{1}{2}) \forall n \in \mathbb{Z}$, and $\tilde{q}_0 = i\hbar\omega_m = i\frac{2\pi}{\beta}(m + \frac{1}{2})$, we obtain the following.

$$\begin{aligned} M &= \int_{-\infty}^{\infty} \frac{d\tilde{k}_0}{2\pi} \frac{1}{\tilde{k}_0 - \xi_{\hat{k}}^{\lambda} \tilde{k}_0 + \tilde{q}_0 - \xi_{\hat{k}+\hat{q}}^{\lambda'}} \\ &= \frac{1}{\beta} \sum_{n \in \mathbb{Z}}^{\mathrm{Fermi}} \frac{1}{i\hbar\omega_n - \xi_{\hat{k}}^{\lambda} i\hbar\omega_n + i\hbar\omega_m - \xi_{\hat{k}+\hat{q}}^{\lambda'}} \\ &= \frac{n_F(\xi_{\hat{k}}^{\lambda}) - n_F(\xi_{\hat{k}+\hat{q}}^{\lambda'})}{i\hbar\omega_m + \xi_{\hat{k}}^{\lambda} - \xi_{\hat{k}+\hat{q}}^{\lambda'}}, \end{aligned} \quad (61)$$

where the Fermi-Dirac distribution is

$$n_F(\xi_{\hat{k}}^{\lambda}) = \frac{1}{e^{\beta\xi_{\hat{k}}^{\lambda}} + 1} = \frac{1}{e^{\beta(\epsilon_{\hat{k}}^{\lambda} - \mu)} + 1}. \quad (62)$$

Note that we have used $n_F(i\hbar\omega_m + \xi_{\hat{k}+\hat{q}}^{\lambda'}) = n_F(\xi_{\hat{k}+\hat{q}}^{\lambda'})$. Using the definition of the electric current operator given in Eq. (50), the Green function obtained in Eq. (60) and the analytically expanded Matsubara sum given in Eq. (61) to the whole upper complex plane by applying the formal change $i\omega_m = \omega \in \mathbb{C}^+$ into the definition of the polarization operator given in Eq. (58), we obtain

$$\Pi^{\mu\nu}(q, \mu, T) = ie^2 \sum_{\lambda, \lambda'} \int \frac{d^d k}{(2\pi)^d} \frac{\langle u_{\hat{k}}^{\lambda} | \hat{v}^{\mu} | u_{\hat{k}+\hat{q}}^{\lambda'} \rangle \langle u_{\hat{k}+\hat{q}}^{\lambda'} | \hat{v}^{\nu} | u_{\hat{k}}^{\lambda} \rangle}{\hbar\omega + \epsilon_{\hat{k}}^{\lambda} - \epsilon_{\hat{k}+\hat{q}}^{\lambda'}} [n_F(\xi_{\hat{k}}^{\lambda}) - n_F(\xi_{\hat{k}+\hat{q}}^{\lambda'})], \quad (63)$$

by using

$$\frac{1}{\hbar\omega + A} - \lim_{\omega \to 0} \frac{1}{\hbar\omega + A} = \frac{1}{\hbar\omega + A} - \frac{1}{A} = \frac{-1}{\hbar\omega + A} \frac{\hbar\omega}{A} \quad (64)$$

$$\frac{-1}{\hbar\omega + A} \frac{1}{A} = \frac{1}{\hbar\omega} \left[ \frac{1}{\hbar\omega + A} - \lim_{\omega \to 0} \frac{1}{\hbar\omega + A} \right], \quad (65)$$

therefore, we obtain $\tilde{\Pi}^{\mu\nu}(q) = \Pi^{\mu\nu}(q) - \lim_{\omega \to 0} \Pi^{\mu\nu}(q)$
tensorial form of $\Pi_{\mu\nu}(\omega, q)$. This is a strong physical requirement, already discussed by V. Abrikosov in 1963 [64] and more recently in [28].

However, sometimes, a different transport relation is proposed [8][9][21][33][42][43][46][65]

$$J_{\mu}(\omega, q) = \Pi_{\mu\nu}(\omega, q) A^{\nu}(\omega, q), \quad (56)$$

This is a natural choice for the transport coefficient provided in Eq. (48). However, as written, this linear relation could possibly lead to non zero electric currents even when the electric field is zero and to spurious dissipationless currents. This pathological result is present in the QFT-b model as is developed and used in literature, and that will be discussed in detail in section Sect. VII.

## B. Polarization Operator

The Kubo formula for the polarization operator is

$$\Pi_{\mu\nu}(q) = \mathrm{Tr} (J_{\mu}(q) J_{\nu}^{*}(q)), \quad (57)$$

which can be reduced to the bubble Feynman diagram as [66][67]

$$\Pi_{\mu\nu}(q) = i \int_k \mathrm{Tr} \left( G_0^{\lambda}(k) \hat{J}_{\mu}^{\lambda, \lambda'}(k) G_0^{\lambda'}(k+q) \hat{J}_{\nu}^{\lambda', \lambda}(k+q) \right) (58)$$

Here we use the definition $\int_k = \int_{-\infty}^{\infty} \frac{d\tilde{k}_0}{2\pi} \int_k$ (see tab. I), $G_0^\lambda(k) = \langle \hat{c}_k^\lambda \hat{c}_k^{\lambda,\dagger} \rangle_0$ is the Green function of the unperturbed Hamiltonian of the system, and $\hat{J}_\mu^{\lambda,\lambda'}(q)$ is the current operator of electronic quasiparticles, defined in Eq. (50). For a general linear hamiltonian of the form $\hbar\omega\psi = \tilde{k}_0\psi = \hat{H}\psi$, the Green function fulfils $(\tilde{k}_0\delta^{\lambda,\lambda} - \hat{H}^{\lambda,\lambda}) G_0^\lambda(r,s) = \delta^{\lambda,\lambda}\delta(r-s)$, and we have the eigenproblem

$$\hat{H}^{\lambda, \lambda'} |u^{\lambda}\rangle = (\epsilon^{\lambda'} - \mu) |u^{\lambda'}\rangle = \xi^{\lambda'} |u^{\lambda'}\rangle, \quad (59)$$

where we have removed the trace operator because it is only applied in the $k$-space and in the $(\lambda, \lambda')$-bands space.

In order to connect with the conductivity obtained from the Kubo formula, using that the Kubo formula for the Ohm law imposes that $\lim_{\omega \to 0} \Pi^{\mu\nu}(\omega) A_\nu(\omega) = 0$ for constant static $A_\nu(\omega)$ [20] (except for superconductors [23][22]), we have to remove the effect of the $\omega \to 0$ limit

note that the chemical potential $\mu$ has been absorbed in $\hat{H}$. The macrocanonical Green function in momentum space is

$$G_0^{\lambda}(k) = \sum_{\lambda} \frac{|u_{\hat{k}}^{\lambda}\rangle\langle u_{\hat{k}}^{\lambda}|}{\tilde{k}_0 - \xi_{\hat{k}}^{\lambda}}. \quad (60)$$

The use of a different (space-time covariant) form of the Green function for the Dirac Hamiltonian leads to different representations of the same result discussed here [9]. After introducing the Green function into the polarization operator, using that $\mathrm{Tr}(ABC) = \mathrm{Tr}(BCA)$, and applying the Matsubara formalism to carry out the $k_0$ integral for the fermionic case by using $\tilde{k}_0 = \hbar\omega_n = \frac{2\pi}{\beta}(n + \frac{1}{2}) \forall n \in \mathbb{Z}$, and $\tilde{q}_0 = i\hbar\omega_m = i\frac{2\pi}{\beta}(m + \frac{1}{2})$, we obtain the following.

$$\begin{aligned} M &= \int_{-\infty}^{\infty} \frac{d\tilde{k}_0}{2\pi} \frac{1}{\tilde{k}_0 - \xi_{\hat{k}}^{\lambda} \tilde{k}_0 + \tilde{q}_0 - \xi_{\hat{k}+\hat{q}}^{\lambda'}} \\ &= \frac{1}{\beta} \sum_{n \in \mathbb{Z}}^{\mathrm{Fermi}} \frac{1}{i\hbar\omega_n - \xi_{\hat{k}}^{\lambda} i\hbar\omega_n + i\hbar\omega_m - \xi_{\hat{k}+\hat{q}}^{\lambda'}} \\ &= \frac{n_F(\xi_{\hat{k}}^{\lambda}) - n_F(\xi_{\hat{k}+\hat{q}}^{\lambda'})}{i\hbar\omega_m + \xi_{\hat{k}}^{\lambda} - \xi_{\hat{k}+\hat{q}}^{\lambda'}}, \end{aligned} \quad (61)$$

where the Fermi-Dirac distribution is

$$n_F(\xi_{\hat{k}}^{\lambda}) = \frac{1}{e^{\beta\xi_{\hat{k}}^{\lambda}} + 1} = \frac{1}{e^{\beta(\epsilon_{\hat{k}}^{\lambda} - \mu)} + 1}. \quad (62)$$

Note that we have used $n_F(i\hbar\omega_m + \xi_{\hat{k}+\hat{q}}^{\lambda'}) = n_F(\xi_{\hat{k}+\hat{q}}^{\lambda'})$. Using the definition of the electric current operator given in Eq. (50), the Green function obtained in Eq. (60) and the analytically expanded Matsubara sum given in Eq. (61) to the whole upper complex plane by applying the formal change $i\omega_m = \omega \in \mathbb{C}^+$ into the definition of the polarization operator given in Eq. (58), we obtain

$$\Pi^{\mu\nu}(q, \mu, T) = ie^2 \sum_{\lambda, \lambda'} \int \frac{d^d k}{(2\pi)^d} \frac{\langle u_{\hat{k}}^{\lambda} | \hat{v}^{\mu} | u_{\hat{k}+\hat{q}}^{\lambda'} \rangle \langle u_{\hat{k}+\hat{q}}^{\lambda'} | \hat{v}^{\nu} | u_{\hat{k}}^{\lambda} \rangle}{\hbar\omega + \epsilon_{\hat{k}}^{\lambda} - \epsilon_{\hat{k}+\hat{q}}^{\lambda'}} [n_F(\xi_{\hat{k}}^{\lambda}) - n_F(\xi_{\hat{k}+\hat{q}}^{\lambda'})], \quad (63)$$

by using

$$\frac{1}{\hbar\omega + A} - \lim_{\omega \to 0} \frac{1}{\hbar\omega + A} = \frac{1}{\hbar\omega + A} - \frac{1}{A} = \frac{-1}{\hbar\omega + A} \frac{\hbar\omega}{A} \quad (64)$$

$$\frac{-1}{\hbar\omega + A} \frac{1}{A} = \frac{1}{\hbar\omega} \left[ \frac{1}{\hbar\omega + A} - \lim_{\omega \to 0} \frac{1}{\hbar\omega + A} \right], \quad (65)$$

therefore, we obtain $\tilde{\Pi}^{\mu\nu}(q) = \Pi^{\mu\nu}(q) - \lim_{\omega \to 0} \Pi^{\mu\nu}(q)$
$$ \tilde{\Pi}^{\mu\nu}(q, \mu, T) = ie^2\hbar\omega \sum_{\lambda, \lambda'} \int_{BZ} \frac{d^d\mathbf{k}}{(2\pi)^d} \frac{\langle u_{\mathbf{k}}^{\lambda} | \hat{v}^{\mu} | u_{\mathbf{k}+\mathbf{q}}^{\lambda'} \rangle \langle u_{\mathbf{k}+\mathbf{q}}^{\lambda'} | \hat{v}^{\nu} | u_{\mathbf{k}}^{\lambda} \rangle}{\hbar\omega + \epsilon_{\mathbf{k}}^{\lambda} - \epsilon_{\mathbf{k}+\mathbf{q}}^{\lambda'}} \frac{n_F(\xi_{\mathbf{k}}^{\lambda}) - n_F(\xi_{\mathbf{k}+\mathbf{q}}^{\lambda'})}{\epsilon_{\mathbf{k}}^{\lambda} - \epsilon_{\mathbf{k}+\mathbf{q}}^{\lambda'}}, \quad (66) $$

which resembles the Kubo formula for the linear conductivity, finally, using Eq. (99), we obtain the conductivity as [34][68][31]

$$ \sigma^{\mu\nu}(q, \Gamma, \mu, T) = -ie^2\hbar \sum_{\lambda, \lambda'} \int_{BZ} \frac{d^d\mathbf{k}}{(2\pi)^d} \frac{\langle u_{\mathbf{k}}^{\lambda} | \hat{v}^{\mu} | u_{\mathbf{k}+\mathbf{q}}^{\lambda'} \rangle \langle u_{\mathbf{k}+\mathbf{q}}^{\lambda'} | \hat{v}^{\nu} | u_{\mathbf{k}}^{\lambda} \rangle}{\hbar\omega + \epsilon_{\mathbf{k}}^{\lambda} - \epsilon_{\mathbf{k}+\mathbf{q}}^{\lambda'}} \frac{n_F(\xi_{\mathbf{k}}^{\lambda}) - n_F(\xi_{\mathbf{k}+\mathbf{q}}^{\lambda'})}{\epsilon_{\mathbf{k}}^{\lambda} - \epsilon_{\mathbf{k}+\mathbf{q}}^{\lambda'}} \quad (67) $$

Taking into account the effect of interactions into the dynamics of the electronic quasiparticles, the Hamiltonian is modified by the self-energy $\Sigma(\omega, q) = \Sigma_R(\omega, q) + i\Sigma_I(\omega, q)$, then, the electronic spectrum is modified by the addition of a real self-energy $\Sigma_R(\omega, q)$ (whose effect is considered small and absorbed into the phenomenological constants of the Hamiltonian) and an imaginary part, which results in a frequency-dependent finite dissipation of the electronic quasiparticles. Here we have argued in subsec. II C that the electronic dissipation is small and can be well approached with the finite lifetime approximation by a constant imaginary dissipation rate $\Gamma^\lambda = \tau_\lambda^{-1}$, then, the electronic quasienergies are modified as $\epsilon_k^\lambda \rightarrow \epsilon_k^\lambda + i\hbar\Gamma^\lambda$ and $\epsilon_{k+q}^{\lambda'} \rightarrow \epsilon_{k+q}^{\lambda'} = \epsilon_{k+q}^{\lambda'} + i\hbar\Gamma'^\lambda$, therefore

$$ \sigma^{\mu\nu}(q, \mu, T) = -ie^2\hbar \sum_{\lambda, \lambda'} \int_{BZ} \frac{d^d\mathbf{k}}{(2\pi)^d} \frac{\langle u_{\mathbf{k}}^{\lambda} | \hat{v}^{\mu} | u_{\mathbf{k}+\mathbf{q}}^{\lambda'} \rangle \langle u_{\mathbf{k}+\mathbf{q}}^{\lambda'} | \hat{v}^{\nu} | u_{\mathbf{k}}^{\lambda} \rangle}{\hbar\omega + i\hbar(\Gamma^{\lambda} + \Gamma'^{\lambda}) + \epsilon_{\mathbf{k}}^{\lambda} - \epsilon_{\mathbf{k}+\mathbf{q}}^{\lambda'}} \frac{n_F(\xi_{\mathbf{k}}^{\lambda} + i\hbar\Gamma^{\lambda}) - n_F(\xi_{\mathbf{k}+\mathbf{q}}^{\lambda'} - i\hbar\Gamma'^{\lambda})}{i\hbar(\Gamma^{\lambda} + \Gamma'^{\lambda}) + \epsilon_{\mathbf{k}}^{\lambda} - \epsilon_{\mathbf{k}+\mathbf{q}}^{\lambda'}} \quad (68) $$

where also the electronic eigenvalues $u_{\mathbf{k}}^{\lambda}$ are functions of $\hbar\Gamma^{\lambda}$. Assuming that the dissipation ranges are small, at first order approximation, by using $\Gamma = \Gamma^{\lambda} + \Gamma'^{\lambda}$, we obtain the Kubo formula we are going to use

$$ \sigma^{\mu\nu}(q, \mu, T) = -ie^2\hbar \sum_{\lambda, \lambda'} \int_{BZ} \frac{d^d\mathbf{k}}{(2\pi)^d} \frac{\langle u_{\mathbf{k}}^{\lambda} | \hat{v}^{\mu} | u_{\mathbf{k}+\mathbf{q}}^{\lambda'} \rangle \langle u_{\mathbf{k}+\mathbf{q}}^{\lambda'} | \hat{v}^{\nu} | u_{\mathbf{k}}^{\lambda} \rangle}{\hbar(\omega + i\Gamma) + \epsilon_{\mathbf{k}}^{\lambda} - \epsilon_{\mathbf{k}+\mathbf{q}}^{\lambda'}} \frac{n_F(\xi_{\mathbf{k}}^{\lambda}) - n_F(\xi_{\mathbf{k}+\mathbf{q}}^{\lambda'})}{\epsilon_{\mathbf{k}}^{\lambda} - \epsilon_{\mathbf{k}+\mathbf{q}}^{\lambda'}} \quad (69) $$

where any possible dependence of $u_{\mathbf{k}}^{\lambda}$ on $\hbar\Gamma^{\lambda}$ has disappeared as well and the effect of the small electronic dissipation is given by the phenomenological $\Gamma > 0$ in this approximation.

In Eq. (69), the dissipation rate is the inverse of the mean lifetime of the electronic quasiparticle $\Gamma = \tau^{-1}$ [69][30][4]. It appears as the imaginary part of $(\epsilon_{\mathbf{k}}^{\lambda} - \epsilon_{\mathbf{k}+\mathbf{q}}^{\lambda'})$, but it is not a bad approximation to consider it as a constant, therefore, $i\Gamma$ can be absorbed into a now complex $\omega \mapsto q_0 = (\omega + i\Gamma)$ [57]. This formula for the conductivity (Eq. (69)) derived from the Random Phase Approximation is completely equivalent to the Kubo formula, which has been derived elsewhere [62][29][28][6][70][71][72][34][73][63][57]. In the following section, we are going to relate the expressions for the conductivity obtained in the different models we compare.

## IV. TENSOR DECOMPOSITION

If we want to compare the results of QFT-b model [8] with the results of the Kubo formula [6], we observe that in the former the results for the non-local polarization operator are written in terms of the component $\Pi_{00}$ and of the quantity $\Pi = q_1^2\Pi_{tr} - q_2^2\Pi_{00}$, while in the latter the results for the non-local conductivity tensor $\sigma_{ij}$ were written in terms of longitudinal ($\sigma_L$), transversal ($\sigma_T$),

Hall ($\sigma_H$) and sinusoidal ($\sigma_S$) components as [32, 74, 75]

$$ \sigma_{ij}(q, \mu, T) = \frac{\tilde{q}_i \tilde{q}_j}{\tilde{q}_\parallel^2} \sigma_L(\omega, \tilde{q}, \mu, T) \\ + \left( \delta_{ij} - \frac{\tilde{q}_i \tilde{q}_j}{\tilde{q}_\parallel^2} \right) \sigma_T(\omega, \tilde{q}, \mu, T) \\ + \epsilon_{ij} \sigma_H(\omega, \tilde{q}, \mu, T) \\ + \left( \frac{\tilde{q}_i \tilde{q}_\ell \epsilon_{\ell j} - \epsilon_{i\ell} \tilde{q}_\ell \tilde{q}_j}{\tilde{q}_\parallel^2} \right) \sigma_S(\omega, \tilde{q}, \mu, T) \quad (70) $$

where this is a sum over repeated indices. In general $\sigma_S \neq 0$, but for graphene it is zero, we show here the complete analysis. Here, $q = (\tilde{q}_0, \tilde{q}_\parallel)$, $\tilde{q}_\parallel = (\tilde{q}_1, \tilde{q}_2)$, $\tilde{q}_\parallel^2 = \sqrt{\tilde{q}_1^2 + \tilde{q}_2^2}$ (see tab. I), $\delta_{ij}$ is the Kronecker delta function and $\epsilon_{ij}$ is the 2D Levi-Civita symbol. However, a relation between $\{\sigma_L, \sigma_T, \sigma_H, \sigma_S\}$ and $\{\sigma_{00}, \sigma_{tr}\}$ can be deduced from Eq. (70) by using the transversality condition $q_\mu \sigma^{\mu\nu}(q) = 0$, where $q_\mu = (q_0, q_1, q_2) = (\omega, q_1, q_2)$ (see tab. I) is the momentum of the quasiparticle, (inherited from the transversality condition of the polarization operator), that can be deduced from the application of the continuity equation $(\partial_\mu j^\mu(x_\mu) = 0)$ for the charge
$$ \tilde{\Pi}^{\mu\nu}(q, \mu, T) = ie^2\hbar\omega \sum_{\lambda, \lambda'} \int_{BZ} \frac{d^d\mathbf{k}}{(2\pi)^d} \frac{\langle u_{\mathbf{k}}^{\lambda} | \hat{v}^{\mu} | u_{\mathbf{k}+\mathbf{q}}^{\lambda'} \rangle \langle u_{\mathbf{k}+\mathbf{q}}^{\lambda'} | \hat{v}^{\nu} | u_{\mathbf{k}}^{\lambda} \rangle}{\hbar\omega + \epsilon_{\mathbf{k}}^{\lambda} - \epsilon_{\mathbf{k}+\mathbf{q}}^{\lambda'}} \frac{n_F(\xi_{\mathbf{k}}^{\lambda}) - n_F(\xi_{\mathbf{k}+\mathbf{q}}^{\lambda'})}{\epsilon_{\mathbf{k}}^{\lambda} - \epsilon_{\mathbf{k}+\mathbf{q}}^{\lambda'}}, \quad (66) $$

which resembles the Kubo formula for the linear conductivity, finally, using Eq. (99), we obtain the conductivity as [34][68][31]

$$ \sigma^{\mu\nu}(q, \Gamma, \mu, T) = -ie^2\hbar \sum_{\lambda, \lambda'} \int_{BZ} \frac{d^d\mathbf{k}}{(2\pi)^d} \frac{\langle u_{\mathbf{k}}^{\lambda} | \hat{v}^{\mu} | u_{\mathbf{k}+\mathbf{q}}^{\lambda'} \rangle \langle u_{\mathbf{k}+\mathbf{q}}^{\lambda'} | \hat{v}^{\nu} | u_{\mathbf{k}}^{\lambda} \rangle}{\hbar\omega + \epsilon_{\mathbf{k}}^{\lambda} - \epsilon_{\mathbf{k}+\mathbf{q}}^{\lambda'}} \frac{n_F(\xi_{\mathbf{k}}^{\lambda}) - n_F(\xi_{\mathbf{k}+\mathbf{q}}^{\lambda'})}{\epsilon_{\mathbf{k}}^{\lambda} - \epsilon_{\mathbf{k}+\mathbf{q}}^{\lambda'}} \quad (67) $$

Taking into account the effect of interactions into the dynamics of the electronic quasiparticles, the Hamiltonian is modified by the self-energy $\Sigma(\omega, q) = \Sigma_R(\omega, q) + i\Sigma_I(\omega, q)$, then, the electronic spectrum is modified by the addition of a real self-energy $\Sigma_R(\omega, q)$ (whose effect is considered small and absorbed into the phenomenological constants of the Hamiltonian) and an imaginary part, which results in a frequency-dependent finite dissipation of the electronic quasiparticles. Here we have argued in subsec. II C that the electronic dissipation is small and can be well approached with the finite lifetime approximation by a constant imaginary dissipation rate $\Gamma^\lambda = \tau_\lambda^{-1}$, then, the electronic quasienergies are modified as $\epsilon_{\mathbf{k}}^{\lambda} \rightarrow \epsilon_{\mathbf{k}}^{\lambda} + i\hbar\Gamma^\lambda$ and $\epsilon_{\mathbf{k}+\mathbf{q}}^{\lambda'} \rightarrow \epsilon_{\mathbf{k}+\mathbf{q}}^{\lambda'} = \epsilon_{\mathbf{k}+\mathbf{q}}^{\lambda'} - i\hbar\Gamma^\lambda$, therefore

$$ \sigma^{\mu\nu}(q, \mu, T) = -ie^2\hbar \sum_{\lambda, \lambda'} \int_{BZ} \frac{d^d\mathbf{k}}{(2\pi)^d} \frac{\langle u_{\mathbf{k}}^{\lambda} | \hat{v}^{\mu} | u_{\mathbf{k}+\mathbf{q}}^{\lambda'} \rangle \langle u_{\mathbf{k}+\mathbf{q}}^{\lambda'} | \hat{v}^{\nu} | u_{\mathbf{k}}^{\lambda} \rangle}{\hbar\omega + i\hbar(\Gamma^{\lambda} + \Gamma^{\lambda'}) + \epsilon_{\mathbf{k}}^{\lambda} - \epsilon_{\mathbf{k}+\mathbf{q}}^{\lambda'}} \frac{n_F(\xi_{\mathbf{k}}^{\lambda} + i\hbar\Gamma^{\lambda}) - n_F(\xi_{\mathbf{k}+\mathbf{q}}^{\lambda'} - i\hbar\Gamma^{\lambda'})}{i\hbar(\Gamma^{\lambda} + \Gamma^{\lambda'}) + \epsilon_{\mathbf{k}}^{\lambda} - \epsilon_{\mathbf{k}+\mathbf{q}}^{\lambda'}} \quad (68) $$

where also the electronic eigenvalues $u_{\mathbf{k}}^{\lambda}$ are functions of $\hbar\Gamma^\lambda$. Assuming that the dissipation ranges are small, at first order approximation, by using $\Gamma = \Gamma^\lambda + \Gamma^{\lambda'}$, we obtain the Kubo formula we are going to use

$$ \sigma^{\mu\nu}(q, \mu, T) = -ie^2\hbar \sum_{\lambda, \lambda'} \int_{BZ} \frac{d^d\mathbf{k}}{(2\pi)^d} \frac{\langle u_{\mathbf{k}}^{\lambda} | \hat{v}^{\mu} | u_{\mathbf{k}+\mathbf{q}}^{\lambda'} \rangle \langle u_{\mathbf{k}+\mathbf{q}}^{\lambda'} | \hat{v}^{\nu} | u_{\mathbf{k}}^{\lambda} \rangle}{\hbar(\omega + i\Gamma) + \epsilon_{\mathbf{k}}^{\lambda} - \epsilon_{\mathbf{k}+\mathbf{q}}^{\lambda'}} \frac{n_F(\xi_{\mathbf{k}}^{\lambda}) - n_F(\xi_{\mathbf{k}+\mathbf{q}}^{\lambda'})}{\epsilon_{\mathbf{k}}^{\lambda} - \epsilon_{\mathbf{k}+\mathbf{q}}^{\lambda'}} \quad (69) $$

where any possible dependence of $u_{\mathbf{k}}^{\lambda}$ on $\hbar\Gamma^\lambda$ has disappeared as well and the effect of the small electronic dissipation is given by the phenomenological $\Gamma > 0$ in this approximation.

In Eq. (69), the dissipation rate is the inverse of the mean lifetime of the electronic quasiparticle $\Gamma = \tau^{-1}$ [69][30][4]. It appears as the imaginary part of $(\epsilon_{\mathbf{k}}^{\lambda} - \epsilon_{\mathbf{k}+\mathbf{q}}^{\lambda'})$, but it is not a bad approximation to consider it as a constant, therefore, $i\Gamma$ can be absorbed into a now complex $\omega \mapsto q_0 = (\omega + i\Gamma)$ [57]. This formula for the conductivity (Eq. (69)) derived from the Random Phase Approximation is completely equivalent to the Kubo formula, which has been derived elsewhere [62][29][28][6][70][71][72][34][73][63][57]. In the following section, we are going to relate the expressions for the conductivity obtained in the different models we compare.

## IV. TENSOR DECOMPOSITION

If we want to compare the results of QFT-b model [8] with the results of the Kubo formula [6], we observe that in the former the results for the non-local polarization operator are written in terms of the component $\Pi_{00}$ and of the quantity $\Pi = q_1^2\Pi_{tr} - q_2^2\Pi_{00}$, while in the latter the results for the non-local conductivity tensor $\sigma_{ij}$ were written in terms of longitudinal ($\sigma_L$), transversal ($\sigma_T$),

Hall ($\sigma_H$) and sinusoidal ($\sigma_S$) components as [32, 74, 75]

$$ \sigma_{ij}(q, \mu, T) = \frac{\tilde{q}_i \tilde{q}_j}{\tilde{q}_\parallel^2} \sigma_L(\omega, \tilde{q}, \mu, T) + \left( \delta_{ij} - \frac{\tilde{q}_i \tilde{q}_j}{\tilde{q}_\parallel^2} \right) \sigma_T(\omega, \tilde{q}, \mu, T) + \epsilon_{ij} \sigma_H(\omega, \tilde{q}, \mu, T) + \left( \frac{\tilde{q}_i \tilde{q}_\ell \epsilon_{\ell j} - \epsilon_{i\ell} \tilde{q}_\ell \tilde{q}_j}{\tilde{q}_\parallel^2} \right) \sigma_S(\omega, \tilde{q}, \mu, T) \quad (70) $$

where this is a sum over repeated indices. In general $\sigma_S \neq 0$, but for graphene it is zero, we show here the complete analysis. Here, $q = (\tilde{q}_0, \tilde{q}_\parallel)$, $\tilde{q}_\parallel = (\tilde{q}_1, \tilde{q}_2)$, $\tilde{q}_\parallel^2 = \sqrt{\tilde{q}_1^2 + \tilde{q}_2^2}$ (see tab. I), $\delta_{ij}$ is the Kronecker delta function and $\epsilon_{ij}$ is the 2D Levi-Civita symbol. However, a relation between $\{\sigma_L, \sigma_T, \sigma_H, \sigma_S\}$ and $\{\sigma_{00}, \sigma_{tr}\}$ can be deduced from Eq. (70) by using the transversality condition $q_\mu \sigma^{\mu\nu}(q) = 0$, where $q_\mu = (q_0, q_1, q_2) = (\omega, q_1, q_2)$ (see tab. I) is the momentum of the quasiparticle, (inherited from the transversality condition of the polarization operator), that can be deduced from the application of the continuity equation $(\partial_\mu j^\mu(x_\mu) = 0)$ for the charge
3-current inside the material together with the Ohm law for linear currents ($j^\mu(q) = \sigma^{\mu\nu}(q)E_\nu(q)$):

$$q_{\mu}j^{\mu}(q) = q_{\mu}\sigma^{\mu\nu}(q)E_{\nu}(q) = 0 \Rightarrow q_{\mu}\sigma^{\mu\nu}(q) = 0 \quad \forall E_{\nu}(71)$$

In all the text that follows, we use the metric tensor $g_{\mu\nu} = \text{diag}\{+1, -1, -1\}$. Separating the temporal component of the 4-vectors, we get $q_\mu\sigma^{\mu\nu}(q) = q_0\sigma^{0\nu}(q) - q_a\sigma^{a\nu}(q) = 0$, therefore, we obtain

$$\sigma^{0\nu}(q) = \frac{q_a\sigma^{a\nu}(q)}{q_0}, \qquad (72)$$

Using Eq. (70), we obtain

$$\begin{aligned} \sigma^{01}(q) &= \frac{q_1\sigma_L(q) + q_2(\sigma_H(q) - \sigma_S(q))}{q_0}, \\ \sigma^{02}(q) &= \frac{q_2\sigma_L(q) - q_1(\sigma_H(q) - \sigma_S(q))}{q_0}. \end{aligned} \qquad (73)$$

Now we can use that the transversality condition is also fulfilled for the second index of the conductivity tensor $\sigma^{\mu\nu}(q)q_\nu = 0$ as well, obtaining that

$$\sigma^{\mu 0}(q) = \frac{\sigma^{\mu a}(q)q_a}{q_0}, \qquad (74)$$

Using again Eq. (70), we obtain

$$\begin{aligned} \sigma^{10}(q) &= \frac{q_1\sigma_L(q) - q_2(\sigma_H(q) + \sigma_S(q))}{q_0}, \\ \sigma^{20}(q) &= \frac{q_2\sigma_L(q) + q_1(\sigma_H(q) + \sigma_S(q))}{q_0}. \end{aligned} \qquad (75)$$

Note that, in general, $\sigma_{\mu\nu}$ is not symmetric because of $\sigma_S$ and of the purely antisymmetric term $\sigma_H$. From those results, we can now derive the 00 component as

$$\sigma^{00}(q) = \frac{q_a\sigma^{a0}(q)}{q_0} = \frac{q_a\sigma^{ab}(q)q_b}{q_0^2} = \frac{q_0^2}{q_0^2}\sigma_L(q). \qquad (76)$$

Now we have derived the full form of the conductivity tensor, we obtain the trace as

$$\text{Tr}(\sigma) = \sigma_{\text{tr}} = g_{\mu\nu}\sigma^{\mu\nu} = -\sigma_T(q) - \sigma_L(q)\frac{q_z^2}{q_0^2}. \qquad (77)$$

where $q_z = (\hbar v_F)^{-1}\sqrt{q_0^2 - \tilde{q}_\parallel^2}$ (see tab. I). From the expressions for $\sigma^{00}$ and $\sigma_{\text{tr}}$, we obtain that

$$\begin{aligned} \sigma_L(q) &= \frac{q_0^2}{q_\parallel^2}\sigma^{00}(q), \\ \sigma_T(q) &= -\sigma_{\text{tr}} - \frac{q_z^2}{q_\parallel^2}\sigma^{00}(q). \end{aligned} \qquad (78)$$

As a conclusion, the conductivity and polarization tensors can be decomposed into the sum of four components [75][46] as

$$\sigma_{\mu\nu} = L_{\mu\nu}\sigma_L + T_{\mu\nu}\sigma_T + H_{\mu\nu}\sigma_H + S_{\mu\nu}\sigma_S, \qquad (79)$$

with

$$L_{\mu\nu} = \left( \delta_{\mu 0} - \frac{q_0 q_\mu}{q_z^2} \right) \frac{q_z^4}{q_\parallel^2 q_0^2} \left( \delta_{\nu 0} - \frac{q_0 q_\nu}{q_z^2} \right), \qquad (80)$$

$$T_{\mu\nu} = \delta_{\mu}{}^{i} \left( \delta_{ij} - \frac{q_i q_j}{q_{\parallel}^2} \right) \delta^{j}_{\nu}, \qquad (81)$$

$$H_{\mu\nu} = (-1)^{\delta_{0\mu} + \delta_{0\nu}} \epsilon_{\mu\nu\rho} \frac{q^\rho}{q_0}, \qquad (82)$$

$$S_{\mu\nu} = \frac{\bar{q}_{\mu}\bar{q}_{\rho}\epsilon_{\rho\nu} - \epsilon_{\mu\rho}\bar{q}_{\rho}\bar{q}_{\nu}}{q_{\parallel}^2}, \qquad (83)$$

where

$$\bar{q}_{\alpha} = q_{\alpha} - \frac{q_{\parallel}^2}{q_0} \delta_{0\alpha}. \qquad (84)$$

So, now we can compare the results of [8] with the results obtained in [6]

## V. KUBO FORMULA FOR GRAPHENE

By using the Kubo formula (Eq. (69)), in [6], the authors obtained the spatial part of the 2D-conductivity tensor for the 2D Dirac cone given by Eq. (35) (see tab. I)

$$\hat{H}_s^{\eta}(\mathbf{k}) = \eta \tilde{k}_1 \tau_1 + \tilde{k}_2 \tau_2 + \tau_3 \Delta_s^{\eta}, \qquad (85)$$

where $m_z = \Delta_s^\eta$. Depending of the indices $s$ and $\eta$, and of external and internal perturbations of the 2D material, each mass-gap $\Delta_s^\eta$ can take different values [50][76][77][78][79], which are zero for suspended and unperturbed graphene sheets. The velocity vector operator is $\hat{v}_i = \partial_{\tilde{k}_i} \hat{H} = \hbar v_F (\eta \tau_x, \tau_y)$ and the eigenstates of the Hamiltonian from Eq. (35) are

$$|u_{\mathbf{k}}^{\lambda}\rangle = \frac{1}{\sqrt{2\epsilon_{\mathbf{k}}^{\lambda}(\epsilon_{\mathbf{k}}^{\lambda} + \Delta_{s}^{\eta})}} \left( -(\Delta_{s}^{\eta} + \epsilon_{\mathbf{k}}^{\lambda}) \frac{\tilde{k}_2 + i\eta\tilde{k}_1}{\tilde{k}} \right), \qquad (86)$$

with corresponding eigenenergies (compare with Eq. (42))

$$\epsilon_{\mathbf{k}}^{\lambda} = \lambda \sqrt{\tilde{k}^2 + (\Delta_s^{\eta})^2}. \qquad (87)$$

If we compare with the Hamiltonian of graphene given in Eq. (30), we see that we can factor this four-spinor hamiltonian into the sum of two two-spinors hamiltonians, one from the second and third rows and columns, and the other from the first and fourth, as indicated here

$$\hat{H}_s^D(\mathbf{k}) = \begin{pmatrix} \tilde{k}_0 + m & 0 & 0 & \tilde{k}_1 - i\tilde{k}_2 \\ 0 & \tilde{k}_0 + m & \tilde{k}_1 + i\tilde{k}_2 & 0 \\ 0 & \tilde{k}_1 - i\tilde{k}_2 & \tilde{k}_0 - m & 0 \\ \tilde{k}_1 + i\tilde{k}_2 & 0 & 0 & \tilde{k}_0 - m \end{pmatrix} \qquad (88)$$

Then, the conductivity for each Dirac cone can be obtained from the results of [6], and the full conductivity
3-current inside the material together with the Ohm law for linear currents ($j^\mu(q) = \sigma^{\mu\nu}(q)E_\nu(q)$):

$$q_{\mu}j^{\mu}(q) = q_{\mu}\sigma^{\mu\nu}(q)E_{\nu}(q) = 0 \Rightarrow q_{\mu}\sigma^{\mu\nu}(q) = 0 \quad \forall E_{\nu}(71)$$

In all the text that follows, we use the metric tensor $g_{\mu\nu} = \text{diag}\{+1, -1, -1\}$. Separating the temporal component of the 4-vectors, we get $q_\mu\sigma^{\mu\nu}(q) = q_0\sigma^{0\nu}(q) - q_a\sigma^{a\nu}(q) = 0$, therefore, we obtain

$$\sigma^{0\nu}(q) = \frac{q_a\sigma^{a\nu}(q)}{q_0}, \qquad (72)$$

Using Eq. (70), we obtain

$$\begin{aligned}
\sigma^{01}(q) &= \frac{q_1\sigma_L(q) + q_2(\sigma_H(q) - \sigma_S(q))}{q_0}, \\
\sigma^{02}(q) &= \frac{q_2\sigma_L(q) - q_1(\sigma_H(q) - \sigma_S(q))}{q_0}.
\end{aligned} \qquad (73)$$

Now we can use that the transversality condition is also fulfilled for the second index of the conductivity tensor $\sigma^{\mu\nu}(q)q_\nu = 0$ as well, obtaining that

$$\sigma^{\mu 0}(q) = \frac{\sigma^{\mu a}(q)q_a}{q_0}, \qquad (74)$$

Using again Eq. (70), we obtain

$$\begin{aligned}
\sigma^{10}(q) &= \frac{q_1\sigma_L(q) - q_2(\sigma_H(q) + \sigma_S(q))}{q_0}, \\
\sigma^{20}(q) &= \frac{q_2\sigma_L(q) + q_1(\sigma_H(q) + \sigma_S(q))}{q_0}.
\end{aligned} \qquad (75)$$

Note that, in general, $\sigma_{\mu\nu}$ is not symmetric because of $\sigma_S$ and of the purely antisymmetric term $\sigma_H$. From those results, we can now derive the 00 component as

$$\sigma^{00}(q) = \frac{q_a\sigma^{a0}(q)}{q_0} = \frac{q_a\sigma^{ab}(q)q_b}{q_0^2} = \frac{q_0^2}{q_0^2}\sigma_L(q). \qquad (76)$$

Now we have derived the full form of the conductivity tensor, we obtain the trace as

$$\text{Tr}(\sigma) = \sigma_{\text{tr}} = g_{\mu\nu}\sigma^{\mu\nu} = -\sigma_T(q) - \sigma_L(q)\frac{q_z^2}{q_0^2}. \qquad (77)$$

where $q_z = (\hbar v_F)^{-1}\sqrt{q_0^2 - \tilde{q}_z^2}$ (see tab. I). From the expressions for $\sigma^{00}$ and $\sigma_{\text{tr}}$, we obtain that

$$\begin{aligned}
\sigma_L(q) &= \frac{q_0^2}{q_0^2} \sigma^{00}(q), \\
\sigma_T(q) &= -\sigma_{\text{tr}} - \frac{q_z^2}{q_0^2} \sigma^{00}(q).
\end{aligned} \qquad (78)$$

As a conclusion, the conductivity and polarization tensors can be decomposed into the sum of four components [75][46] as

$$\sigma_{\mu\nu} = L_{\mu\nu}\sigma_L + T_{\mu\nu}\sigma_T + H_{\mu\nu}\sigma_H + S_{\mu\nu}\sigma_S, \qquad (79)$$

with

$$L_{\mu\nu} = \left( \delta_{\mu 0} - \frac{q_0 q_\mu}{q_z^2} \right) \frac{q_z^4}{q_0^2 q_0^2} \left( \delta_{\nu 0} - \frac{q_0 q_\nu}{q_z^2} \right), \qquad (80)$$

$$T_{\mu\nu} = \delta_{\mu}{}^{i} \left( \delta_{ij} - \frac{q_i q_j}{q_{\parallel}^2} \right) \delta^{j}_{\nu}, \qquad (81)$$

$$H_{\mu\nu} = (-1)^{\delta_{0\mu} + \delta_{0\nu}} \epsilon_{\mu\nu\rho} \frac{q^\rho}{q_0}, \qquad (82)$$

$$S_{\mu\nu} = \frac{\bar{q}_{\mu}\bar{q}_{\rho}\epsilon_{\rho\nu} - \epsilon_{\mu\rho}\bar{q}_{\rho}\bar{q}_{\nu}}{q_{\parallel}^2}, \qquad (83)$$

where

$$\bar{q}_{\alpha} = q_{\alpha} - \frac{q_{\parallel}^2}{q_0} \delta_{0\alpha}. \qquad (84)$$

So, now we can compare the results of [8] with the results obtained in [6]

## V. KUBO FORMULA FOR GRAPHENE

By using the Kubo formula (Eq. (69)), in [6], the authors obtained the spatial part of the 2D-conductivity tensor for the 2D Dirac cone given by Eq. (35) (see tab. I)

$$\hat{H}_s^{\eta}(\mathbf{k}) = \eta \tilde{k}_1 \tau_1 + \tilde{k}_2 \tau_2 + \tau_3 \Delta_s^{\eta}, \qquad (85)$$

where $m_z = \Delta_s^{\eta}$. Depending of the indices $s$ and $\eta$, and of external and internal perturbations of the 2D material, each mass-gap $\Delta_s^{\eta}$ can take different values [50][76][77][78][79], which are zero for suspended and unperturbed graphene sheets. The velocity vector operator is $\hat{v}_i = \partial_{\tilde{k}_i} \hat{H} = \hbar v_F (\eta \tau_x, \tau_y)$ and the eigenstates of the Hamiltonian from Eq. (35) are

$$|u_{\mathbf{k}}^{\lambda}\rangle = \frac{1}{\sqrt{2\epsilon_{\mathbf{k}}^{\lambda}(\epsilon_{\mathbf{k}}^{\lambda} + \Delta_{s}^{\eta})}} \left( -(\Delta_{s}^{\eta} + \epsilon_{\mathbf{k}}^{\lambda}) \frac{\tilde{k}_2 + i\eta\tilde{k}_1}{\tilde{k}} \right), \qquad (86)$$

with corresponding eigenenergies (compare with Eq. (42))

$$\epsilon_{\mathbf{k}}^{\lambda} = \lambda \sqrt{\tilde{k}^2 + (\Delta_s^{\eta})^2}. \qquad (87)$$

If we compare with the Hamiltonian of graphene given in Eq. (30), we see that we can factor this four-spinor hamiltonian into the sum of two two-spinors hamiltonians, one from the second and third rows and columns, and the other from the first and fourth, as indicated here

$$\hat{H}_s^D(\mathbf{k}) = \begin{pmatrix} \tilde{k}_0 + m & 0 & 0 & \tilde{k}_1 - i\tilde{k}_2 \\ 0 & \tilde{k}_0 + m & \tilde{k}_1 + i\tilde{k}_2 & 0 \\ 0 & \tilde{k}_1 - i\tilde{k}_2 & \tilde{k}_0 - m & 0 \\ \tilde{k}_1 + i\tilde{k}_2 & 0 & 0 & \tilde{k}_0 - m \end{pmatrix} \qquad (88)$$

Then, the conductivity for each Dirac cone can be obtained from the results of [6], and the full conductivity
will be the sum of the contribution of the four cones (two cones due to the factorization of Eq. (30) shown in Eq. (88), each one counted two times because of the spin degeneration $g_s = 2$). The velocity-velocity correlators are given in [6], and in Eq. (105). Those results are valid for all frequencies $\omega \in \mathbb{C}^+$, and we remind that they are the conductivity per Dirac cone, to obtain the conductivity of Graphene, we must to sum the contribution of each 4 Dirac cones, having into account their respective (signed) Dirac masses by using

$$ \sigma_p(q) = \sum_{\eta=\pm} \sum_{s=\pm} \sigma_p(q, \Delta_s^\eta), \qquad (89) $$

for $p = \{L, T, H\}$, $\sigma_S(q, \Delta_s^\eta) = 0$ for any Hamiltonian of the form of Eq. (85). From the decomposition shown in Eq. (88) of $\hat{H}_s^D(k)$ into two 2-spinor Hamiltonians $\hat{H}_s^\eta(k)$ of the form of Eq. (85) with Dirac masses $\Delta_s^\eta = \eta m$ (rotate the 2-spinor Hamiltonian obtained from the first and fourth rows and columns $\pi/2$ rads), the Chern number of the studied model of graphene is

$$ C = \sum_{\eta=\pm} \sum_{s=\pm} \operatorname{sgn}(\eta m) = 0. \qquad (90) $$

Therefore, graphene with the induced mass studied here is topologically trivial, and there is not any Hall conductivity ($\sigma_H = 0$).

Due to the requirements of causality and realism, $\sigma_{\mu\nu}$ do not have poles for $\omega$ in the upper complex plane, and Eq. (69) is valid for all complex-values frequency with positive imaginary part simply promoting $\omega \in \mathbb{R}$ as a complex variable $\omega \in \mathbb{C}^+$. In [6], it was proven that the spatial components of the conductivity tensor $\sigma_{ij}$ can be conveniently given by separating between longitudinal $\sigma_L$, transverse $\sigma_T$, and Hall $\sigma_H$, contributions [6][32, 74, 75] (Eq. (70))

$$ \sigma_{ij}(q, \mu, T) = \frac{\tilde{q}_i \tilde{q}_j}{\tilde{q}_{\parallel}^2} \sigma_L(\omega, \tilde{q}, \mu, T) + \left( \delta_{ij} - \frac{\tilde{q}_i \tilde{q}_j}{\tilde{q}_{\parallel}^2} \right) \sigma_T(\omega, \tilde{q}, \mu, T) + \epsilon_{ij} \sigma_H(\omega, \tilde{q}, \mu, T), \qquad (91) $$

because $\sigma_S(q, \Delta_s^\eta) = 0$ for any Hamiltonian of the form of Eq. (85). Here, $q = (\tilde{q}_0, \tilde{q}_{\parallel})$, $\tilde{q}_{\parallel} = (\tilde{q}_1, \tilde{q}_2)$, $\tilde{q}_{\parallel}^2 = \sqrt{\tilde{q}_1^2 + \tilde{q}_2^2}$ (see tab. I), $\delta_{ij}$ is the Kronecker delta function and $\epsilon_{ij}$ is the 2D Levi-Civita symbol. This expression has been generalized in Eq. (79) of Sect. IV for the complete conductivity tensor $\sigma_{\mu\nu}$. The explicit analytical form of those three functions for real and complex frequencies in the zero temperature limit can be found in [6] and in the appendix A. To obtain similar results for finite temperatures, we should apply the Maldague formula [80][81]

$$ \sigma_{ij}(q, \mu, T) = \int_{-\infty}^{\infty} dE \frac{\sigma_{ij}(q, E, 0)}{4k_B T \cosh^2\left(\frac{E-\mu}{2k_B T}\right)}, \qquad (92) $$

where $\sigma_{ij}(q, \mu, 0)$ is the zero-temperature conductivity result. This is the more general formula for the linear non-local conductivity based on the linearized tight-binding model with a constant dissipation time parameter $\tau = \Gamma^{-1}$, and it is completely equivalent to Eq. (69). To go beyond this result, the full tight-binding model of graphene should be used [82][59] instead of the linear approximation, a deep more detailed study of the effects of the different interactions in electronic quasiparticle spectrum [10][11][12][60] or more detailed ab-initio models [83][84][59].

## VI. LOCAL LIMIT OF THE KUBO CONDUCTIVITY

There is a special case in the local limit (when we apply the $q_{\parallel} \to 0$ limit to the Kubo formula (69)) when results valid for all temperatures can be obtained. The local limit of the conductivities of one massive Dirac cone are [6]

$$ \begin{aligned} \sigma_{xx}(\omega, \mathbf{0}, \mu, 0) &= i \frac{\sigma_0}{\pi} \left[ \frac{\mu^2 - \Delta^2}{|\mu|} \frac{1}{\Omega} \Theta(|\mu| - |\Delta|) \right. \\ & \qquad \left. + \frac{\Delta^2}{M\Omega} - \frac{\Omega^2 + 4\Delta^2}{2i\Omega^2} \tan^{-1}\left(\frac{i\Omega}{2M}\right) \right], \\ \sigma_{xy}(\omega, \mathbf{0}, \mu, 0) &= \frac{2\sigma_0}{\pi} \frac{\eta\Delta}{i\Omega} \tan^{-1}\left(\frac{i\Omega}{2M}\right). \end{aligned} \qquad (93) $$

$\sigma_0 = \alpha c/4$ is the universal conductivity of graphene ($\alpha = \frac{e^2}{\hbar c}$ is the fine structure constant), $\Omega = \hbar\omega + i\hbar\Gamma$ and $M = \text{Max}(|\Delta|, |\mu|)$. These results are per Dirac cone and are consistent with those found by other researchers [68][29][69][30][28][31][32][33][85][86][47][43][87][88][89]. The first term in $\sigma_{xx}$ corresponds to intraband transitions, and the last two terms correspond to interband transitions. Note that, in the local limit $q_{\parallel} = 0$ one obtains $\sigma_{xx}(\omega, \mathbf{0}) = \sigma_{yy}(\omega, \mathbf{0}) = \sigma_L(\omega, \mathbf{0}) = \sigma_T(\omega, \mathbf{0})$, and $\sigma_{xy}(\omega, \mathbf{0}) = -\sigma_{yx}(\omega, \mathbf{0}) = \sigma_H(\omega, \mathbf{0})$.

This model not only serves to model the local conductivity of graphene with mass, but also any other 2D Dirac cones, like the surface states of a Three Dimensional Topological Insulator [90] and Chern Insulators [91]. By summing the contribution of several different 2D Dirac cones with $\Delta \neq 0$, non-trivial topological states with non-zero Chern number can be studied [79][91].

Here we are going to compare this result with the Falkovsky model of massless graphene [31], we need to obtain the $\Delta \to 0$ limit of the sum of the contribution of 4 Dirac cones, then

$$ \begin{aligned} \sigma_{xx}(\omega, q_{\parallel} = 0, T = 0, \Delta = 0) &= \sigma_{xx}^{\text{intra}}(\omega) + \sigma_{xx}^{\text{inter}}(\omega), \\ \sigma_{xx}^{\text{intra}}(\omega) &= \frac{\alpha c}{\pi} \frac{|\mu|}{-i\Omega}, \\ \sigma_{xx}^{\text{inter}}(\omega) &= \frac{\alpha c}{2\pi} \tan^{-1}\left(\frac{-i\Omega}{2|\mu|}\right), \\ \sigma_{xy}(\omega, q_{\parallel} = 0, T = 0, \Delta = 0) &= 0. \end{aligned} \qquad (94) $$
will be the sum of the contribution of the four cones (two cones due to the factorization of Eq. (30) shown in Eq. (88), each one counted two times because of the spin degeneration $g_s = 2$). The velocity-velocity correlators are given in [6], and in Eq. (105). Those results are valid for all frequencies $\omega \in \mathbb{C}^+$, and we remind that they are the conductivity per Dirac cone, to obtain the conductivity of Graphene, we must to sum the contribution of each 4 Dirac cones, having into account their respective (signed) Dirac masses by using

$$ \sigma_p(q) = \sum_{\eta=\pm} \sum_{s=\pm} \sigma_p(q, \Delta_s^\eta), \qquad (89) $$

for $p = \{L, T, H\}$, $\sigma_S(q, \Delta_s^\eta) = 0$ for any Hamiltonian of the form of Eq. (85). From the decomposition shown in Eq. (88) of $\hat{H}_s^D(k)$ into two 2-spinor Hamiltonians $\hat{H}_s^\eta(k)$ of the form of Eq. (85) with Dirac masses $\Delta_s^\eta = \eta m$ (rotate the 2-spinor Hamiltonian obtained from the first and fourth rows and columns $\pi/2$ rads), the Chern number of the studied model of graphene is

$$ C = \sum_{\eta=\pm} \sum_{s=\pm} \operatorname{sgn}(\eta m) = 0. \qquad (90) $$

Therefore, graphene with the induced mass studied here is topologically trivial, and there is not any Hall conductivity ($\sigma_H = 0$).

Due to the requirements of causality and realism, $\sigma_{\mu\nu}$ do not have poles for $\omega$ in the upper complex plane, and Eq. (69) is valid for all complex-values frequency with positive imaginary part simply promoting $\omega \in \mathbb{R}$ as a complex variable $\omega \in \mathbb{C}^+$. In [6], it was proven that the spatial components of the conductivity tensor $\sigma_{ij}$ can be conveniently given by separating between longitudinal $\sigma_L$, transverse $\sigma_T$, and Hall $\sigma_H$, contributions [6][32, 74, 75] (Eq. (70))

$$ \sigma_{ij}(q, \mu, T) = \frac{\tilde{q}_i \tilde{q}_j}{\tilde{q}_{\parallel}^2} \sigma_L(\omega, \tilde{q}, \mu, T) \\ + \left( \delta_{ij} - \frac{\tilde{q}_i \tilde{q}_j}{\tilde{q}_{\parallel}^2} \right) \sigma_T(\omega, \tilde{q}, \mu, T) \\ + \epsilon_{ij} \sigma_H(\omega, \tilde{q}, \mu, T), \qquad (91) $$

because $\sigma_S(q, \Delta_s^\eta) = 0$ for any Hamiltonian of the form of Eq. (85). Here, $q = (\tilde{q}_0, \tilde{q}_{\parallel})$, $\tilde{q}_{\parallel} = (\tilde{q}_1, \tilde{q}_2)$, $\tilde{q}_{\parallel}^2 = \sqrt{\tilde{q}_1^2 + \tilde{q}_2^2}$ (see tab. I), $\delta_{ij}$ is the Kronecker delta function and $\epsilon_{ij}$ is the 2D Levi-Civita symbol. This expression has been generalized in Eq. (79) of Sect. IV for the complete conductivity tensor $\sigma_{\mu\nu}$. The explicit analytical form of those three functions for real and complex frequencies in the zero temperature limit can be found in [6] and in the appendix A. To obtain similar results for finite temperatures, we should apply the Maldague formula [80][81]

$$ \sigma_{ij}(q, \mu, T) = \int_{-\infty}^{\infty} dE \frac{\sigma_{ij}(q, E, 0)}{4k_B T \cosh^2\left(\frac{E-\mu}{2k_B T}\right)}, \qquad (92) $$

where $\sigma_{ij}(q, \mu, 0)$ is the zero-temperature conductivity result. This is the more general formula for the linear non-local conductivity based on the linearized tight-binding model with a constant dissipation time parameter $\tau = \Gamma^{-1}$, and it is completely equivalent to Eq. (69). To go beyond this result, the full tight-binding model of graphene should be used [82][59] instead of the linear approximation, a deep more detailed study of the effects of the different interactions in electronic quasiparticle spectrum [10][11][12][60] or more detailed ab-initio models [83][84][59].

## VI. LOCAL LIMIT OF THE KUBO CONDUCTIVITY

There is a special case in the local limit (when we apply the $q_{\parallel} \to 0$ limit to the Kubo formula (69)) when results valid for all temperatures can be obtained. The local limit of the conductivities of one massive Dirac cone are [6]

$$ \sigma_{xx}(\omega, \mathbf{0}, \mu, 0) = i \frac{\sigma_0}{\pi} \left[ \frac{\mu^2 - \Delta^2}{|\mu|} \frac{1}{\Omega} \Theta(|\mu| - |\Delta|) \right. \\ \left. + \frac{\Delta^2}{M\Omega} - \frac{\Omega^2 + 4\Delta^2}{2i\Omega^2} \tan^{-1}\left(\frac{i\Omega}{2M}\right) \right], \\ \sigma_{xy}(\omega, \mathbf{0}, \mu, 0) = \frac{2\sigma_0}{\pi} \frac{\eta\Delta}{i\Omega} \tan^{-1}\left(\frac{i\Omega}{2M}\right). \qquad (93) $$

$\sigma_0 = \alpha c/4$ is the universal conductivity of graphene ($\alpha = \frac{e^2}{\hbar c}$ is the fine structure constant), $\Omega = \hbar\omega + i\hbar\Gamma$ and $M = \text{Max}(|\Delta|, |\mu|)$. These results are per Dirac cone and are consistent with those found by other researchers [68][29][69][30][28][31][32][33][85][86][47][43][87][88][89]. The first term in $\sigma_{xx}$ corresponds to intraband transitions, and the last two terms correspond to interband transitions. Note that, in the local limit $q_{\parallel} = 0$ one obtains $\sigma_{xx}(\omega, \mathbf{0}) = \sigma_{yy}(\omega, \mathbf{0}) = \sigma_L(\omega, \mathbf{0}) = \sigma_T(\omega, \mathbf{0})$, and $\sigma_{xy}(\omega, \mathbf{0}) = -\sigma_{yx}(\omega, \mathbf{0}) = \sigma_H(\omega, \mathbf{0})$.

This model not only serves to model the local conductivity of graphene with mass, but also any other 2D Dirac cones, like the surface states of a Three Dimensional Topological Insulator [90] and Chern Insulators [91]. By summing the contribution of several different 2D Dirac cones with $\Delta \neq 0$, non-trivial topological states with non-zero Chern number can be studied [79][91].

Here we are going to compare this result with the Falkovsky model of massless graphene [31], we need to obtain the $\Delta \to 0$ limit of the sum of the contribution of 4 Dirac cones, then

$$ \sigma_{xx}(\omega, q_{\parallel} = 0, T = 0, \Delta = 0) = \sigma_{xx}^{\text{intra}}(\omega) + \sigma_{xx}^{\text{inter}}(\omega), \\ \sigma_{xx}^{\text{intra}}(\omega) = \frac{\alpha c}{\pi} \frac{|\mu|}{-i\Omega}, \\ \sigma_{xx}^{\text{inter}}(\omega) = \frac{\alpha c}{2\pi} \tan^{-1}\left(\frac{-i\Omega}{2|\mu|}\right), \\ \sigma_{xy}(\omega, q_{\parallel} = 0, T = 0, \Delta = 0) = 0. \qquad (94) $$
By applying the Maldague formula (Eq. (92)) to this result with $\Delta = 0$, in the non-dissipation limit for the interband term and for finite temperatures, the well-known result of Falkovsky [31][30] is obtained as

$$ \sigma_{xx}^{F}(\omega, \Gamma, T, \mu) = \sigma_{xx}^{F, \text{intra}}(\omega, \Gamma, T, \mu) + \sigma_{xx}^{F, \text{inter}}(\omega, 0, T, \mu), $$

$$ \sigma_{xx}^{F, \text{intra}}(\omega, \Gamma, T, \mu) = \frac{2i\alpha ck_B T}{\pi \hbar \omega + i\Gamma} \ln \left[ 2 \cosh \left( \frac{\mu}{2k_B T} \right) \right], $$

$$ \sigma_{xx}^{F, \text{inter}}(\omega, 0, T, \mu) = \frac{\alpha c}{4} G\left(\left|\frac{\hbar\omega}{2}\right|\right) \quad (95) \\ + i \frac{\alpha c}{4} \frac{4\hbar\omega}{\pi} \int_{0}^{\infty} d\xi \frac{G(\xi) - G\left(\left|\frac{\hbar\omega}{2}\right|\right)}{(\hbar\omega)^2 - 4\xi^2}, $$

with

$$ G(\epsilon) = \frac{\sinh(\beta\epsilon)}{\cosh(\beta\mu) + \cosh(\beta\epsilon)}. \quad (96) $$

We have the absolute value in the preceding formula because we are handling also negative real frequencies, while in [31] the result for only positive frequencies were derived. The derivation of this result for real and imaginary frequencies is given in the appendix B.

## VII. ON THE QUANTUM FIELD THEORY BASED MODEL FOR THE CONDUCTIVITY AVAILABLE IN CURRENT LITERATURE

In literature exists another QFT-b model for the non-local conductivity also based on the polarization operator of a $2D + 1$ Dirac Hamiltonian [8][9][21][33][42][43][46], but different from the one we just derived in section Sect. V. That model has been extensively used in the context of Casimir effect, compared to experimental results [8], and it has been proposed as intrinsically more fundamental than the Kubo one since based on "first principles" and being "not phenomenological" like [18, 19]. This same model, due to a claimed coherence with the Lifshitz theory has also been suggested as model to modify the well known dielectric function of metals itself [19]. In this section we show that that model, first does not include unavoidable losses coming from the inelastic interactions of electronic quasiparticles with different objects always present in real samples as phonons, scattering centers, lattice dislocations and non-linear interactions for example, as discussed in Subsect. II C. To add such losses at this level of theory it is enough to simply introduce a constant parameter $\Gamma = \tau^{-1}$ as the inverse of the mean lifetime of the electron quasiparticle, which is an experimentally measured quantity (See Subsect. II C). In addition to that, that theory is also not physically correct and predicts nonphysical features like an intrinsic Plasma behaviour that cannot be cured even by adding losses. We explain the origin of that pathology and show that a correct regularization of that model (as done in

Sect. III) make that QFT-b model exactly identical to the Kubo model we derived in previous section Sect. V. The starting point is the linear relationship given in Eq. (56)

$$ J_{\mu}(\omega, \mathbf{q}) = \Pi_{\mu\nu}(\omega, \mathbf{q}) A^{\nu}(\omega, \mathbf{q}). \quad (97) $$

To connect this result with the microscopic Ohm law

$$ J_{\mu}(\omega, \mathbf{q}) = \sigma_{\mu\nu}(\mathbf{q}) E^{\nu}(\omega, \mathbf{q}), \quad (98) $$

where $\sigma_{\mu\nu}(\mathbf{q})$ is the conductivity operator, and $E^{\mu}(\omega, \mathbf{q}) = -i\omega A^{\mu}(\omega, \mathbf{q})$ is the electric field, the relation between $\Pi_{\mu\nu}(\mathbf{q})$ and $\sigma_{\mu\nu}(\mathbf{q})$ is assumed to be [65][92][93][94][95][96] [8] [21][19]

$$ \sigma_{\mu\nu}(\mathbf{q}) = \frac{\Pi_{\mu\nu}(\mathbf{q})}{-i\omega}, \quad (99) $$

with the Polarization operator defined in Eq. (58), but using the covariant action and covariant Dirac Hamiltonian given in Eqs. (37) and Eq. (31). Note that, contrary to the result obtained in Eq. (55) and in the appendix Sect. D, there is no additional regularization term $\lim_{\omega \to 0} \Pi_{\mu\nu}(\omega, \mathbf{q}_{\parallel})$. Remarkably, we show that this additional regularization term makes the results for the transversal conductivity different from the results of the Kubo formalism.

In this formalism, the Polarization operator is also defined in Eq. (58), but using the covariant action and covariant Dirac Hamiltonian given in Eq. (37) and Eq. (31). In addition to that, the main difference with the polarization operator obtained in Eq. (66), is that, instead of using the expression of the Green function as an expansion on eigenfunctions given in Eq. (60), a covariant form of the Green function of the Dirac Hamiltonian (Eq. (40), see tab. I) is used instead

$$ \mathcal{G}_{0}^{W}(\mathbf{k}) = \hat{\mathcal{H}}_{W}^{-1}(\mathbf{k}) = \frac{1}{\gamma^{\mu} K_{\mu} - \gamma^{4} m} = \frac{\gamma^{\mu} K_{\mu} + \gamma^{4} m}{K^{\mu} K_{\mu} - m^{2}} \quad (100) $$

This is the reason why apparently the results obtained in [8][41][42][9][21][32][33][46][47] look completely different to the ones obtained by the use of the Kubo formula. Here we are going to show how the two formalisms are related, and under what circumstances they provide similar or different results for the conductivity for graphene with a topologically trivial mass term.

By using the definition of the current operator (Eq. (50)), we obtain that

$$ \begin{aligned} \hat{J}^{\mu}(\mathbf{k}) &= \frac{\partial \hat{\mathcal{H}}_{s}^{W}(k_{\alpha} - eA_{\alpha}(\mathbf{k}))}{\partial A_{\mu}(\mathbf{k})} \\ &= -e(\gamma^{0}, v_{F}\gamma^{1}, v_{F}\gamma^{2}) \\ &= -\frac{e}{\hbar}(\tilde{\gamma}^{0}, \tilde{\gamma}^{1}, \tilde{\gamma}^{2}). \end{aligned} \quad (101) $$

Inserting this result into Eq. (58), and using that the Green function is diagonal in spin, we get

$$ \Pi^{\mu\nu}(\mathbf{q}) = g_s \frac{e^2}{\hbar^2} i \int_k \mathrm{Tr} (\mathcal{G}_0^W(\mathbf{k}) \tilde{\gamma}^\mu \mathcal{G}_0^W(\mathbf{k}+\mathbf{q}) \tilde{\gamma}^\nu), \quad (102) $$
By applying the Maldague formula (Eq. (92)) to this result with $\Delta = 0$, in the non-dissipation limit for the interband term and for finite temperatures, the well-known result of Falkovsky [31][30] is obtained as

$$\sigma_{xx}^{F}(\omega, \Gamma, T, \mu) = \sigma_{xx}^{F, \text{intra}}(\omega, \Gamma, T, \mu) + \sigma_{xx}^{F, \text{inter}}(\omega, 0, T, \mu),$$

$$\sigma_{xx}^{F, \text{intra}}(\omega, \Gamma, T, \mu) = \frac{2i\alpha ck_B T}{\pi \hbar \omega + i\Gamma} \ln \left[ 2 \cosh \left( \frac{\mu}{2k_B T} \right) \right],$$

$$\sigma_{xx}^{F, \text{inter}}(\omega, 0, T, \mu) = \frac{\alpha c}{4} G\left(\left|\frac{\hbar\omega}{2}\right|\right) \quad (95) \\ + i \frac{\alpha c}{4} \frac{4\hbar\omega}{\pi} \int_{0}^{\infty} d\xi \frac{G(\xi) - G\left(\left|\frac{\hbar\omega}{2}\right|\right)}{(\hbar\omega)^2 - 4\xi^2},$$

with

$$G(\epsilon) = \frac{\sinh(\beta\epsilon)}{\cosh(\beta\mu) + \cosh(\beta\epsilon)}. \quad (96)$$

We have the absolute value in the preceding formula because we are handling also negative real frequencies, while in [31] the result for only positive frequencies were derived. The derivation of this result for real and imaginary frequencies is given in the appendix B.

## VII. ON THE QUANTUM FIELD THEORY BASED MODEL FOR THE CONDUCTIVITY AVAILABLE IN CURRENT LITERATURE

In literature exists another QFT-b model for the non-local conductivity also based on the polarization operator of a $2D+1$ Dirac Hamiltonian [8][9][21][33][42][43][46], but different from the one we just derived in section Sect. V. That model has been extensively used in the context of Casimir effect, compared to experimental results [8], and it has been proposed as intrinsically more fundamental than the Kubo one since based on "first principles" and being "not phenomenological" like [18, 19]. This same model, due to a claimed coherence with the Lifshitz theory has also been suggested as model to modify the well known dielectric function of metals itself [19]. In this section we show that that model, first does not include unavoidable losses coming from the inelastic interactions of electronic quasiparticles with different objects always present in real samples as phonons, scattering centers, lattice dislocations and non-linear interactions for example, as discussed in Subsect. II C. To add such losses at this level of theory it is enough to simply introduce a constant parameter $\Gamma = \tau^{-1}$ as the inverse of the mean lifetime of the electron quasiparticle, which is an experimentally measured quantity (See Subsect. II C). In addition to that, that theory is also not physically correct and predicts nonphysical features like an intrinsic Plasma behaviour that cannot be cured even by adding losses. We explain the origin of that pathology and show that a correct regularization of that model (as done in

Sect. III) make that QFT-b model exactly identical to the Kubo model we derived in previous section Sect. V. The starting point is the linear relationship given in Eq. (56)

$$J_{\mu}(\omega, \mathbf{q}) = \Pi_{\mu\nu}(\omega, \mathbf{q}) A^{\nu}(\omega, \mathbf{q}). \quad (97)$$

To connect this result with the microscopic Ohm law

$$J_{\mu}(\omega, \mathbf{q}) = \sigma_{\mu\nu}(\mathbf{q}) E^{\nu}(\omega, \mathbf{q}), \quad (98)$$

where $\sigma_{\mu\nu}(\mathbf{q})$ is the conductivity operator, and $E^{\mu}(\omega, \mathbf{q}) = -i\omega A^{\mu}(\omega, \mathbf{q})$ is the electric field, the relation between $\Pi_{\mu\nu}(\mathbf{q})$ and $\sigma_{\mu\nu}(\mathbf{q})$ is assumed to be [65][92][93][94][95][96] [8] [21][19]

$$\sigma_{\mu\nu}(\mathbf{q}) = \frac{\Pi_{\mu\nu}(\mathbf{q})}{-i\omega}, \quad (99)$$

with the Polarization operator defined in Eq. (58), but using the covariant action and covariant Dirac Hamiltonian given in Eqs. (37) and Eq. (31). Note that, contrary to the result obtained in Eq. (55) and in the appendix Sect. D, there is no additional regularization term $\lim_{\omega \to 0} \Pi_{\mu\nu}(\omega, \mathbf{q}_{\parallel})$. Remarkably, we show that this additional regularization term makes the results for the transversal conductivity different from the results of the Kubo formalism.

In this formalism, the Polarization operator is also defined in Eq. (58), but using the covariant action and covariant Dirac Hamiltonian given in Eq. (37) and Eq. (31). In addition to that, the main difference with the polarization operator obtained in Eq. (66), is that, instead of using the expression of the Green function as an expansion on eigenfunctions given in Eq. (60), a covariant form of the Green function of the Dirac Hamiltonian (Eq. (40), see tab. I) is used instead

$$\mathcal{G}_{0}^{W}(\mathbf{k}) = \hat{\mathcal{H}}_{W}^{-1}(\mathbf{k}) = \frac{1}{\gamma^{\mu} K_{\mu} - \gamma^{4} m} = \frac{\gamma^{\mu} K_{\mu} + \gamma^{4} m}{K^{\mu} K_{\mu} - m^{2}} \quad (100)$$

This is the reason why apparently the results obtained in [8][41][42][9][21][32][33][46][47] look completely different to the ones obtained by the use of the Kubo formula. Here we are going to show how the two formalisms are related, and under what circumstances they provide similar or different results for the conductivity for graphene with a topologically trivial mass term.

By using the definition of the current operator (Eq. (50)), we obtain that

$$\begin{aligned} \hat{J}^{\mu}(\mathbf{k}) &= \frac{\partial \hat{\mathcal{H}}_{s}^{W}(k_{\alpha} - eA_{\alpha}(\mathbf{k}))}{\partial A_{\mu}(\mathbf{k})} \\ &= -e(\gamma^{0}, v_{F}\gamma^{1}, v_{F}\gamma^{2}) \\ &= -\frac{e}{\hbar}(\tilde{\gamma}^{0}, \tilde{\gamma}^{1}, \tilde{\gamma}^{2}). \end{aligned} \quad (101)$$

Inserting this result into Eq. (58), and using that the Green function is diagonal in spin, we get

$$\Pi^{\mu\nu}(\mathbf{q}) = g_s \frac{e^2}{\hbar^2} i \int_k \mathrm{Tr} (\mathcal{G}_0^W(\mathbf{k}) \tilde{\gamma}^\mu \mathcal{G}_0^W(\mathbf{k}+\mathbf{q}) \tilde{\gamma}^\nu), \quad (102)$$
where $g_s = 2$ is the spin degeneration, $e^2 = \alpha c\hbar$, and $\tilde{\gamma}_\mu = \hbar(\gamma_0, v_F\gamma)$ (see tab. I). After carrying out the trace, one obtains

$$ \Pi^{\mu\nu}(q) = ig_s \frac{\alpha c}{\hbar} \int_k \frac{Z^{\mu\nu}(K_\alpha, q_\alpha)}{[K^\mu K_\mu - m^2][S^\mu S_\mu - m^2]}, \quad (103) $$

where $S_\mu = K_\mu + Q_\mu$ (see tab. I), and $Z^{\mu\nu}(K_\alpha, q_\alpha)$ is obtained as [9]

$$ Z^{\mu\nu}(K_\alpha, q_\alpha) = 4 \begin{pmatrix} \tilde{k}_\mu \tilde{s}_\mu + \tilde{k}_1 \tilde{s}_1 + \tilde{k}_2 \tilde{s}_2 + m^2 & -v_F (\tilde{k}_\mu \tilde{s}_1 + \tilde{k}_1 \tilde{s}_\mu) & -v_F (\tilde{k}_\mu \tilde{s}_2 + \tilde{k}_2 \tilde{s}_\mu) \\ -v_F (\tilde{k}_\mu \tilde{s}_1 + \tilde{k}_1 \tilde{s}_\mu) & v_F^2 (\tilde{k}_\mu \tilde{s}_\mu + \tilde{k}_1 \tilde{s}_1 - \tilde{k}_2 \tilde{s}_2 - m^2) & v_F^2 (\tilde{k}_1 \tilde{s}_2 + \tilde{k}_2 \tilde{s}_1) \\ -v_F (\tilde{k}_\mu \tilde{s}_2 + \tilde{k}_2 \tilde{s}_\mu) & v_F^2 (\tilde{k}_1 \tilde{s}_2 + \tilde{k}_2 \tilde{s}_1) & v_F^2 (\tilde{k}_\mu \tilde{s}_\mu + \tilde{k}_2 \tilde{s}_2 - \tilde{k}_1 \tilde{s}_1 - m^2) \end{pmatrix} \quad (104) $$

where $\tilde{k}_\mu = \tilde{k}_0 + \mu$, $\tilde{s}_\mu = \tilde{s}_0 + \mu$, $\tilde{k}_i = \hbar v_F k_i$ and $\tilde{s}_i = \hbar v_F s_i$ (see tab. I). From this result, it is easy to obtain [40][65]

$$ \begin{aligned} \frac{Z^{\mu\nu}(K_\alpha, q_\alpha)}{N(K_\alpha, q_\alpha)} &= \sum_{\lambda, \lambda'} \begin{pmatrix} 1 + \frac{\tilde{k}_1 \tilde{s}_1 + \tilde{k}_2 \tilde{s}_2 + m^2}{\epsilon_k^\lambda \epsilon_s^{\lambda'}} & v_F \left( \frac{\tilde{k}_1}{\epsilon_k^\lambda} + \frac{\tilde{s}_1}{\epsilon_s^{\lambda'}} \right) & v_F \left( \frac{\tilde{k}_2}{\epsilon_k^\lambda} + \frac{\tilde{s}_2}{\epsilon_s^{\lambda'}} \right) \\ v_F \left( \frac{\tilde{k}_1}{\epsilon_k^\lambda} + \frac{\tilde{s}_1}{\epsilon_s^{\lambda'}} \right) & v_F^2 \left[ 1 + \frac{\tilde{k}_1 \tilde{s}_1 - \tilde{k}_2 \tilde{s}_2 - m^2}{\epsilon_k^\lambda \epsilon_s^{\lambda'}} \right] & v_F^2 \frac{\tilde{k}_1 \tilde{s}_2 + \tilde{k}_2 \tilde{s}_1}{\epsilon_k^\lambda \epsilon_s^{\lambda'}} \\ v_F \left( \frac{\tilde{k}_2}{\epsilon_k^\lambda} + \frac{\tilde{s}_2}{\epsilon_s^{\lambda'}} \right) & v_F^2 \frac{\tilde{k}_1 \tilde{s}_2 + \tilde{k}_2 \tilde{s}_1}{\epsilon_k^\lambda \epsilon_s^{\lambda'}} & v_F^2 \left[ 1 + \frac{\tilde{k}_2 \tilde{s}_2 - \tilde{k}_1 \tilde{s}_1 - m^2}{\epsilon_k^\lambda \epsilon_s^{\lambda'}} \right] \end{pmatrix} \frac{1}{\tilde{k}_0 + \xi_k^\lambda \tilde{s}_0 + \xi_s^{\lambda'}} \\ &= 2 \sum_{\lambda, \lambda'} \frac{\langle u_k^\lambda | \hat{v}^\mu | u_{k+q}^{\lambda'} \rangle \langle u_{k+q}^{\lambda'} | \hat{v}^\nu | u_k^\lambda \rangle}{[\tilde{k}_0 - \xi_k^\lambda][\tilde{s}_0 - \xi_s^{\lambda'}]}, \end{aligned} \quad (105) $$

where $N(K_\alpha, q_\alpha) = [K^\mu K_\mu - m^2][S^\mu S_\mu - m^2]$ is the denominator of Eq. (103), $|u_k^\lambda\rangle$ are the eigenfunctions of the spinor Hamiltonian (35) with $\Delta \to m$ defined in Eq. (86) and $\xi_k^\lambda$ are their corresponding eigenvalues, defined in Eq. (59) with Eq. (87). The product of velocity correlators coincide with the ones obtained in [6], with $m$ instead of $\Delta$, proving explicitly that, even if we start from the covariant Hamiltonian, the Polarization operator (and, therefore, the conductivity) for graphene must be the same as the one obtained with the expression of the Green function as an expansion on eigenfunctions in Eq. (60). However, from this expression we observe that there is a difference between the results obtained from the Kubo formula (Eq. (69)) and the results derived directly from Eq. (103). In the former we have regularized the polarization in such a way that $\lim_{\omega \to 0} \Pi^{\mu\nu} A_\nu = 0$ is imposed, while in the latter there is not such regularization, and it is related with Eq. (63) instead. As a consequence, the results derived from Eq. (69) are regular for small frequencies with a finite DC conductivity, while the results derived directly from Eq. (103) without regularization has an infinite DC conductivity coming from the infinite dissipation time and from an spurious Plasma behavior whose origin is an interband transition.

Applying the Matsubara formalism directly to the expression of Eq. (103) using Eq. (104) instead of Eq. (105), we obtain the expression for the polarization operator shown in [9]. The detailed calculations are shown in the

## Appendix C, and results in (see tab. I)

$$ \Pi^{\mu\nu}(q) = ig_s e^2 \hbar^2 \int_k [1 - N_\mu(\epsilon_k)] \times \sum_{\lambda=\pm} \frac{Z_{\mu\nu}(\epsilon_k^\lambda, k, q_0, q)}{2\epsilon_k ((\tilde{q}_0 + \epsilon_k^\lambda)^2 - \epsilon_s^2)}, \quad (106) $$

where we have defined, using Eq. (96)

$$ N_\mu(\epsilon) = \sum_{\eta=\pm} n_F(\epsilon + \eta\mu) = 1 - G(\epsilon). \quad (107) $$

Following the notation of [9], the Polarization operator given in Eq. (103) can be written as

$$ \Pi^{\mu\nu}(q) = \Pi_0^{\mu\nu}(q) + \Delta_T \Pi^{\mu\nu}(q), \quad (108) $$

by construction, $\Pi_0^{\mu\nu}(q)$ is independent of temperature and of the chemical potential $\mu$ [9], it can be understood as the interband contribution with $\mu = k_B T = 0$ eV. Note that $\Pi_0^{\mu\nu}(q)$, as shown in Eq. (106), has an ultra-violet divergence, which can be removed with a Pauli-Villars subtraction scheme [32], solved by a $1/N$ expansion [75, 97] or solving the regular integral given in Eq. (63) [37][6][65]. The regularized $\Pi_0^{\mu\nu}(q)$ will only coincide with the regular result obtained from Eq. (66) using Eq. (69) and the relation between Polarization and conductivity show in Eq. (99) [6] if the regularization term derived from the Kubo formula in Eq. (D16) of appendix D does not contribute to the end result. We will observe in what follows that this is not always the case.
where $g_s = 2$ is the spin degeneration, $e^2 = \alpha c\hbar$, and $\tilde{\gamma}_\mu = \hbar(\gamma_0, v_F\gamma)$ (see tab. I). After carrying out the trace, one obtains

$$ \Pi^{\mu\nu}(q) = ig_s \frac{\alpha c}{\hbar} \int_k \frac{Z^{\mu\nu}(K_\alpha, q_\alpha)}{[K^\mu K_\mu - m^2][S^\mu S_\mu - m^2]}, \quad (103) $$

where $S_\mu = K_\mu + Q_\mu$ (see tab. I), and $Z^{\mu\nu}(K_\alpha, q_\alpha)$ is obtained as [9]

$$ Z^{\mu\nu}(K_\alpha, q_\alpha) = 4 \begin{pmatrix} \tilde{k}_\mu \tilde{s}_\mu + \tilde{k}_1 \tilde{s}_1 + \tilde{k}_2 \tilde{s}_2 + m^2 & -v_F (\tilde{k}_\mu \tilde{s}_1 + \tilde{k}_1 \tilde{s}_\mu) & -v_F (\tilde{k}_\mu \tilde{s}_2 + \tilde{k}_2 \tilde{s}_\mu) \\ -v_F (\tilde{k}_\mu \tilde{s}_1 + \tilde{k}_1 \tilde{s}_\mu) & v_F^2 (\tilde{k}_\mu \tilde{s}_\mu + \tilde{k}_1 \tilde{s}_1 - \tilde{k}_2 \tilde{s}_2 - m^2) & v_F^2 (\tilde{k}_1 \tilde{s}_2 + \tilde{k}_2 \tilde{s}_1) \\ -v_F (\tilde{k}_\mu \tilde{s}_2 + \tilde{k}_2 \tilde{s}_\mu) & v_F^2 (\tilde{k}_1 \tilde{s}_2 + \tilde{k}_2 \tilde{s}_1) & v_F^2 (\tilde{k}_\mu \tilde{s}_\mu + \tilde{k}_2 \tilde{s}_2 - \tilde{k}_1 \tilde{s}_1 - m^2) \end{pmatrix} \quad (104) $$

where $\tilde{k}_\mu = \tilde{k}_0 + \mu$, $\tilde{s}_\mu = \tilde{s}_0 + \mu$, $\tilde{k}_i = \hbar v_F k_i$ and $\tilde{s}_i = \hbar v_F s_i$ (see tab. I). From this result, it is easy to obtain [40][65]

$$ \begin{aligned} \frac{Z^{\mu\nu}(K_\alpha, q_\alpha)}{N(K_\alpha, q_\alpha)} &= \sum_{\lambda, \lambda'} \begin{pmatrix} 1 + \frac{\tilde{k}_1 \tilde{s}_1 + \tilde{k}_2 \tilde{s}_2 + m^2}{\epsilon_k^\lambda \epsilon_s^{\lambda'}} & v_F \left( \frac{\tilde{k}_1}{\epsilon_k^\lambda} + \frac{\tilde{s}_1}{\epsilon_s^{\lambda'}} \right) & v_F \left( \frac{\tilde{k}_2}{\epsilon_k^\lambda} + \frac{\tilde{s}_2}{\epsilon_s^{\lambda'}} \right) \\ v_F \left( \frac{\tilde{k}_1}{\epsilon_k^\lambda} + \frac{\tilde{s}_1}{\epsilon_s^{\lambda'}} \right) & v_F^2 \left[ 1 + \frac{\tilde{k}_1 \tilde{s}_1 - \tilde{k}_2 \tilde{s}_2 - m^2}{\epsilon_k^\lambda \epsilon_s^{\lambda'}} \right] & v_F^2 \frac{\tilde{k}_1 \tilde{s}_2 + \tilde{k}_2 \tilde{s}_1}{\epsilon_k^\lambda \epsilon_s^{\lambda'}} \\ v_F \left( \frac{\tilde{k}_2}{\epsilon_k^\lambda} + \frac{\tilde{s}_2}{\epsilon_s^{\lambda'}} \right) & v_F^2 \frac{\tilde{k}_1 \tilde{s}_2 + \tilde{k}_2 \tilde{s}_1}{\epsilon_k^\lambda \epsilon_s^{\lambda'}} & v_F^2 \left[ 1 + \frac{\tilde{k}_2 \tilde{s}_2 - \tilde{k}_1 \tilde{s}_1 - m^2}{\epsilon_k^\lambda \epsilon_s^{\lambda'}} \right] \end{pmatrix} \frac{1}{\tilde{k}_0 + \xi_k^\lambda \tilde{s}_0 + \xi_s^{\lambda'}} \\ &= 2 \sum_{\lambda, \lambda'} \frac{\langle u_k^\lambda | \hat{v}^\mu | u_{k+q}^{\lambda'} \rangle \langle u_{k+q}^{\lambda'} | \hat{v}^\nu | u_k^\lambda \rangle}{[\tilde{k}_0 - \xi_k^\lambda][\tilde{s}_0 - \xi_s^{\lambda'}]}, \end{aligned} \quad (105) $$

where $N(K_\alpha, q_\alpha) = [K^\mu K_\mu - m^2][S^\mu S_\mu - m^2]$ is the denominator of Eq. (103), $|u_k^\lambda\rangle$ are the eigenfunctions of the spinor Hamiltonian (35) with $\Delta \to m$ defined in Eq. (86) and $\xi_k^\lambda$ are their corresponding eigenvalues, defined in Eq. (59) with Eq. (87). The product of velocity correlators coincide with the ones obtained in [6], with $m$ instead of $\Delta$, proving explicitly that, even if we start from the covariant Hamiltonian, the Polarization operator (and, therefore, the conductivity) for graphene must be the same as the one obtained with the expression of the Green function as an expansion on eigenfunctions in Eq. (60). However, from this expression we observe that there is a difference between the results obtained from the Kubo formula (Eq. (69)) and the results derived directly from Eq. (103). In the former we have regularized the polarization in such a way that $\lim_{\omega \to 0} \Pi^{\mu\nu} A_\nu = 0$ is imposed, while in the latter there is not such regularization, and it is related with Eq. (63) instead. As a consequence, the results derived from Eq. (69) are regular for small frequencies with a finite DC conductivity, while the results derived directly from Eq. (103) without regularization has an infinite DC conductivity coming from the infinite dissipation time and from an spurious Plasma behavior whose origin is an interband transition.

Applying the Matsubara formalism directly to the expression of Eq. (103) using Eq. (104) instead of Eq. (105), we obtain the expression for the polarization operator shown in [9]. The detailed calculations are shown in the

## Appendix C, and results in (see tab. I)

$$ \Pi^{\mu\nu}(q) = ig_s e^2 \hbar^2 \int_k [1 - N_\mu(\epsilon_k)] \times \sum_{\lambda=\pm} \frac{Z_{\mu\nu}(\epsilon_k^\lambda, k, q_0, q)}{2\epsilon_k ((\tilde{q}_0 + \epsilon_k^\lambda)^2 - \epsilon_s^2)}, \quad (106) $$

where we have defined, using Eq. (96)

$$ N_{\mu}(\epsilon) = \sum_{\eta=\pm} n_{F}(\epsilon + \eta\mu) = 1 - G(\epsilon). \quad (107) $$

Following the notation of [9], the Polarization operator given in Eq. (103) can be written as

$$ \Pi^{\mu\nu}(q) = \Pi_0^{\mu\nu}(q) + \Delta_T \Pi^{\mu\nu}(q), \quad (108) $$

by construction, $\Pi_0^{\mu\nu}(q)$ is independent of temperature and of the chemical potential $\mu$ [9], it can be understood as the interband contribution with $\mu = k_B T = 0$ eV. Note that $\Pi_0^{\mu\nu}(q)$, as shown in Eq. (106), has an ultra-violet divergence, which can be removed with a Pauli-Villars subtraction scheme [32], solved by a $1/N$ expansion [75, 97] or solving the regular integral given in Eq. (63) [37][6][65]. The regularized $\Pi_0^{\mu\nu}(q)$ will only coincide with the regular result obtained from Eq. (66) using Eq. (69) and the relation between Polarization and conductivity show in Eq. (99) [6] if the regularization term derived from the Kubo formula in Eq. (D16) of appendix D does not contribute to the end result. We will observe in what follows that this is not always the case.
In the appendix C, we derive the expressions for the longitudinal and transversal parts of the polarization operator derived from Eq. (106). The longitudinal polarization is obtained as

$$ \Pi_L(q) = -2ig_s \frac{\alpha c}{\hbar} \frac{\tilde{q}_0^2}{\tilde{q}_{\parallel}^2} \int_k \frac{1}{2\epsilon_k} [1 - N_{\mu}(\epsilon_k)] \\ \times \sum_{\lambda=\pm} \left[ 1 + \frac{M_{00}(\tilde{q}_0, \tilde{k}_{\parallel}, \tilde{q}_{\parallel})}{Q(\tilde{q}_0, \tilde{k}_{\parallel}, \tilde{q}_{\parallel}) + 2\tilde{k}_{\parallel} \cdot \tilde{q}_{\parallel}} \right], \quad (109) $$

with $\tilde{q}_z^2 = \tilde{q}_0^2 - \tilde{q}_{\parallel}^2$ (see tab. I)

$$ M_{00}(\tilde{q}_0, \tilde{k}_{\parallel}, \tilde{q}_{\parallel}) = -\tilde{q}_z^2 + 4\tilde{q}_0\epsilon_k^{\lambda} + 4\epsilon_k^2, \quad (110) $$

$$ Q(\tilde{q}_0, \tilde{k}_{\parallel}, \tilde{q}_{\parallel}) = -\tilde{q}_z^2 - 2\tilde{q}_0\epsilon_k^{\lambda}, \quad (111) $$

while the transversal polarization is

$$ \Pi_T(q) = 2ig_s \frac{\alpha c}{\hbar} \frac{\tilde{q}_0^2}{\tilde{q}_{\parallel}^2} \int_k \frac{1}{2\epsilon_k} [1 - N_{\mu}(\epsilon_k)] \\ \times \sum_{\lambda=\pm} \left[ 1 + \frac{M_{00}(\tilde{q}_0, \tilde{k}_{\parallel}, \tilde{q}_{\parallel}) - 4 \frac{\tilde{q}_{\parallel}^2}{\tilde{q}_0^2} (\tilde{k}_{\parallel}^2 + \tilde{q}_0 \epsilon_k^{\lambda})}{Q(\tilde{q}_0, \tilde{k}_{\parallel}, \tilde{q}_{\parallel}) + 2\tilde{k}_{\parallel} \cdot \tilde{q}_{\parallel}} \right] \quad (112) $$

Using the parameters $\tilde{q}_z = \sqrt{\tilde{q}_0^2 - \tilde{q}_{\parallel}^2}$ and $\delta = \frac{2m}{-i\tilde{q}_z}$ for any complex frequency $\omega \in \mathbb{C}^+$ ($\tilde{q}_0 = \hbar(\omega + i\Gamma)$), and

$$ \Psi(x) = 2 \left[ x + (1 - x^2) \tan^{-1} \left( \frac{1}{x} \right) \right], \quad (113) $$

the regularization of $\Pi_0^{\mu\nu}(q)$ in Eq. (109) and Eq. (112) lead to [8][9][21][92]

$$ \Pi_L^{(0)}(q) = -g_s i \frac{\alpha c}{8\pi\hbar} \frac{\tilde{q}_0^2}{\tilde{q}_z} \Psi(\delta), \quad (114) $$

$$ \Pi_T^{(0)}(q) = -g_s i \frac{\alpha c}{8\pi\hbar} \tilde{q}_z \Psi(\delta). \quad (115) $$

In the particular case of imaginary frequencies, we have $\tilde{q}_0 = i\hbar\xi = i\Xi$, $\tilde{q}_z = i\sqrt{\Xi^2 + \tilde{q}_{\parallel}^2} = i\tilde{\theta}_z = E$, $\gamma = \frac{\Xi}{\tilde{\theta}_z}$ and $\delta = \frac{2m}{\tilde{\theta}_z}$ (see tab. I), and it is shown in the appendix that the integrals can be reduced to

$$ \Pi_L^{(0)}(q) = g_s \frac{\alpha c}{8\pi\hbar} \frac{\Xi^2}{\tilde{\theta}_z} \Psi(\delta), \quad (116) $$

$$ \Delta_T \Pi_L(q) = g_s \frac{\alpha c}{2\pi \hbar} \frac{\tilde{\theta}_z}{\tilde{q}_{\parallel}^2} \int_{\delta}^{\infty} du N_{\mu} \left( \tilde{\theta}_z \frac{u}{2} \right) \\ \times \left[ 1 - \mathrm{Re} \left[ \frac{1 - u^2 - 2i\gamma u}{\sqrt{1 - u^2 - 2i\gamma u + (1 - \gamma^2)\delta^2}} \right] \right], \quad (117) $$

$$ \Pi_T^{(0)}(q) = g_s \frac{\alpha c}{8\pi\hbar} \tilde{\theta}_z \Psi(\delta), \quad (118) $$

$$ \Delta_T \Pi_T(q) = -g_s \frac{\alpha c}{2\pi\hbar} \frac{\Xi^2}{\tilde{q}_{\parallel}^2} \tilde{\theta}_z \int_{\delta}^{\infty} du N_{\mu} \left( \tilde{\theta}_z \frac{u}{2} \right) \\ \times \left[ 1 - \mathrm{Re} \left[ \frac{(1 + i\gamma^{-1}u)^2 + (\gamma^{-2} - 1)\delta^2}{\sqrt{1 - u^2 - 2i\gamma u + (1 - \gamma^2)\delta^2}} \right] \right]. \quad (119) $$

These results coincide with the results published in [8][9][21][92][94] in their appropriate limits.

### A. Differences with Kubo formula

With this result, we obtain that $\Pi_L(q)$ obtained from this QFT-b model is consistent with the conductivity obtained from the non-local Kubo formula, however, $\Pi_T(q)$ is not because $\lim_{\omega \to 0} \Pi_T(\omega, q) \neq 0$, therefore it must be correctly regularized to avoid nonphysical results.

Having into account that $\Pi_T^{(0)}(q)$ corresponds to the $\mu = k_B T = 0$ eV case, it corresponds to the interband conductivity. Therefore, the real part of this conductivity must be zero when $\hbar\Gamma \to 0$ (an electron in the valence band has to jump to a hole place in the conduction band to conduct; therefore, a finite gap exists as long as valence and conduction bands do not tough). In [6] was found that the interband conductivity can be written as

$$ \lim_{\hbar\Gamma \to 0} \mathrm{Re} [\sigma_T(q)] = f(q) \Theta \left( \hbar\omega - \sqrt{4m^2 + \tilde{q}_{\parallel}^2} \right), \quad (120) $$

this functional form is needed to obtain the correct convergence with the local limit. However, from the expression published in [9], we find that

$$ \sigma_{T,0}^{QFT-b}(\omega, \tilde{q}_{\parallel}) = \frac{\alpha c}{4\pi} \frac{\sqrt{\tilde{q}_{\parallel}^2 - \hbar^2\omega^2}}{-i\hbar\omega} \Phi \left( \frac{2\Delta}{\sqrt{\tilde{q}_{\parallel}^2 - \hbar^2\omega^2}} \right) \quad (121) $$

diverges at small frequencies as

$$ \sigma_{T,0}^{QFT-b}(\omega, \tilde{q}_{\parallel}) = \frac{\alpha c}{4\pi} \frac{\tilde{q}_{\parallel}}{-i\hbar\omega} \Phi \left( \frac{2\Delta}{\tilde{q}_{\parallel}} \right) + \mathcal{O}[\hbar\omega]. \quad (122) $$

with $\tilde{q}_{\parallel} = \hbar v_F q_{\parallel}$. Note that $\sigma_{T,0}^{QFT-b}(\omega, \tilde{q}_{\parallel})$ behaves as a Plasma model without conduction electrons or holes. This result of the QFT-b theory would implies a dissipation-less electric current which is clearly not acceptable in normal materials. This is an explicit example of the need to regularize the polarization operator in such a way to fulfil the condition imposed by Gauge invariance $\lim_{\omega \to 0} \Pi^{\mu\nu}(q) A_{\nu}(\omega, q) = 0 \forall A_{\nu}(0, q)$. Therefore, we keep the $\mu = 0$ interband transversal conductivity derived in [6] as the correct conductivity $\sigma_{T,0}$ term, given by

$$ \sigma_{T,0}(q) = \frac{\alpha c}{4\pi} \frac{1}{\hbar\xi} \left[ \tilde{\theta}_z \Psi(\delta) - \tilde{q}_{\parallel} \Psi(x) \right], \quad (123) $$

with $\tilde{\theta}_z = \sqrt{(\hbar\xi)^2 + \tilde{q}_{\parallel}^2}$, $\delta = \frac{2\Delta}{\tilde{\theta}_z}$ and $x = \frac{2\Delta}{\tilde{q}_{\parallel}}$ (see tab. I). Note that this expression corresponds to the explicit elimination of the $\omega \to 0$ limit to the Polarization operator
In the appendix C, we derive the expressions for the longitudinal and transversal parts of the polarization operator derived from Eq. (106). The longitudinal polarization is obtained as

$$ \Pi_L(q) = -2ig_s \frac{\alpha c}{\hbar} \frac{\tilde{q}_0^2}{\tilde{q}_\parallel^2} \int_k \frac{1}{2\epsilon_k} [1 - N_\mu(\epsilon_k)] \\ \times \sum_{\lambda=\pm} \left[ 1 + \frac{M_{00}(\tilde{q}_0, \tilde{k}_\parallel, \tilde{q}_\parallel)}{Q(\tilde{q}_0, \tilde{k}_\parallel, \tilde{q}_\parallel) + 2\tilde{k}_\parallel \cdot \tilde{q}_\parallel} \right], \quad (109) $$

with $\tilde{q}_z^2 = \tilde{q}_0^2 - \tilde{q}_\parallel^2$ (see tab. I)

$$ M_{00}(\tilde{q}_0, \tilde{k}_\parallel, \tilde{q}_\parallel) = -\tilde{q}_z^2 + 4\tilde{q}_0\epsilon_k^\lambda + 4\epsilon_k^2, \quad (110) $$

$$ Q(\tilde{q}_0, \tilde{k}_\parallel, \tilde{q}_\parallel) = -\tilde{q}_z^2 - 2\tilde{q}_0\epsilon_k^\lambda, \quad (111) $$

while the transversal polarization is

$$ \Pi_T(q) = 2ig_s \frac{\alpha c}{\hbar} \frac{\tilde{q}_0^2}{\tilde{q}_\parallel^2} \int_k \frac{1}{2\epsilon_k} [1 - N_\mu(\epsilon_k)] \\ \times \sum_{\lambda=\pm} \left[ 1 + \frac{M_{00}(\tilde{q}_0, \tilde{k}_\parallel, \tilde{q}_\parallel) - 4 \frac{\tilde{q}_\parallel^2}{\tilde{q}_0^2} (\tilde{k}_\parallel^2 + \tilde{q}_0 \epsilon_k^\lambda)}{Q(\tilde{q}_0, \tilde{k}_\parallel, \tilde{q}_\parallel) + 2\tilde{k}_\parallel \cdot \tilde{q}_\parallel} \right] \quad (112) $$

Using the parameters $\tilde{q}_z = \sqrt{\tilde{q}_0^2 - \tilde{q}_\parallel^2}$ and $\delta = \frac{2m}{-i\tilde{q}_z}$ for any complex frequency $\omega \in \mathbb{C}^+$ ($\tilde{q}_0 = \hbar(\omega + i\Gamma)$), and

$$ \Psi(x) = 2 \left[ x + (1 - x^2) \tan^{-1} \left( \frac{1}{x} \right) \right], \quad (113) $$

the regularization of $\Pi_0^{\mu\nu}(q)$ in Eq. (109) and Eq. (112) lead to [8][9][21][92]

$$ \Pi_L^{(0)}(q) = -g_s i \frac{\alpha c}{8\pi\hbar} \frac{\tilde{q}_0^2}{\tilde{q}_z} \Psi(\delta), \quad (114) $$

$$ \Pi_T^{(0)}(q) = -g_s i \frac{\alpha c}{8\pi\hbar} \tilde{q}_z \Psi(\delta). \quad (115) $$

In the particular case of imaginary frequencies, we have $\tilde{q}_0 = i\hbar\xi = i\Xi$, $\tilde{q}_z = i\sqrt{\Xi^2 + \tilde{q}_\parallel^2} = i\tilde{\theta}_z = E$, $\gamma = \frac{\Xi}{\tilde{\theta}_z}$ and $\delta = \frac{2m}{\tilde{\theta}_z}$ (see tab. I), and it is shown in the appendix that the integrals can be reduced to

$$ \Pi_L^{(0)}(q) = g_s \frac{\alpha c}{8\pi\hbar} \frac{\Xi^2}{\tilde{\theta}_z} \Psi(\delta), \quad (116) $$

$$ \Delta_T \Pi_L(q) = g_s \frac{\alpha c}{2\pi \hbar} \frac{\tilde{\theta}_z}{\tilde{q}_\parallel^2} \int_\delta^\infty du N_\mu \left(\tilde{\theta}_z \frac{u}{2}\right) \\ \times \left[ 1 - \mathrm{Re} \left[ \frac{1 - u^2 - 2i\gamma u}{\sqrt{1 - u^2 - 2i\gamma u + (1 - \gamma^2)\delta^2}} \right] \right], \quad (117) $$

$$ \Pi_T^{(0)}(q) = g_s \frac{\alpha c}{8\pi\hbar} \tilde{\theta}_z \Psi(\delta), \quad (118) $$

$$ \Delta_T \Pi_T(q) = -g_s \frac{\alpha c}{2\pi\hbar} \frac{\Xi^2}{\tilde{q}_\parallel^2} \tilde{\theta}_z \int_\delta^\infty du N_\mu \left(\tilde{\theta}_z \frac{u}{2}\right) \\ \times \left[ 1 - \mathrm{Re} \left[ \frac{(1 + i\gamma^{-1}u)^2 + (\gamma^{-2} - 1)\delta^2}{\sqrt{1 - u^2 - 2i\gamma u + (1 - \gamma^2)\delta^2}} \right] \right]. \quad (119) $$

These results coincide with the results published in [8][9][21][92][94] in their appropriate limits.

### A. Differences with Kubo formula

With this result, we obtain that $\Pi_L(q)$ obtained from this QFT-b model is consistent with the conductivity obtained from the non-local Kubo formula, however, $\Pi_T(q)$ is not because $\lim_{\omega \to 0} \Pi_T(\omega, q) \neq 0$, therefore it must be correctly regularized to avoid nonphysical results.

Having into account that $\Pi_T^{(0)}(q)$ corresponds to the $\mu = k_B T = 0$ eV case, it corresponds to the interband conductivity. Therefore, the real part of this conductivity must be zero when $\hbar\Gamma \to 0$ (an electron in the valence band has to jump to a hole place in the conduction band to conduct; therefore, a finite gap exists as long as valence and conduction bands do not touch). In [6] was found that the interband conductivity can be written as

$$ \lim_{\hbar\Gamma \to 0} \mathrm{Re} [\sigma_T(q)] = f(q)\Theta(\hbar\omega - \sqrt{4m^2 + \tilde{q}_\parallel^2}), \quad (120) $$

this functional form is needed to obtain the correct convergence with the local limit. However, from the expression published in [9], we find that

$$ \sigma_{T,0}^{QFT-b}(\omega, \tilde{q}_\parallel) = \frac{\alpha c}{4\pi} \frac{\sqrt{\tilde{q}_\parallel^2 - \hbar^2\omega^2}}{-i\hbar\omega} \Phi\left(\frac{2\Delta}{\sqrt{\tilde{q}_\parallel^2 - \hbar^2\omega^2}}\right) \quad (121) $$

diverges at small frequencies as

$$ \sigma_{T,0}^{QFT-b}(\omega, \tilde{q}_\parallel) = \frac{\alpha c}{4\pi} \frac{\tilde{q}_\parallel}{-i\hbar\omega} \Phi\left(\frac{2\Delta}{\tilde{q}_\parallel}\right) + \mathcal{O}[\hbar\omega]. \quad (122) $$

with $\tilde{q}_\parallel = \hbar v_F q_\parallel$. Note that $\sigma_{T,0}^{QFT-b}(\omega, \tilde{q}_\parallel)$ behaves as a Plasma model without conduction electrons or holes. This result of the QFT-b theory would implies a dissipation-less electric current which is clearly not acceptable in normal materials. This is an explicit example of the need to regularize the polarization operator in such a way to fulfil the condition imposed by Gauge invariance $\lim_{\omega \to 0} \Pi^{\mu\nu}(q) A_\nu(\omega, q) = 0 \forall A_\nu(0, q)$. Therefore, we keep the $\mu = 0$ interband transversal conductivity derived in [6] as the correct conductivity $\sigma_{T,0}$ term, given by

$$ \sigma_{T,0}(q) = \frac{\alpha c}{4\pi} \frac{1}{\hbar\xi} \left[ \tilde{\theta}_z \Psi(\delta) - \tilde{q}_\parallel \Psi(x) \right], \quad (123) $$

with $\tilde{\theta}_z = \sqrt{(\hbar\xi)^2 + \tilde{q}_\parallel^2}$, $\delta = \frac{2\Delta}{\tilde{\theta}_z}$ and $x = \frac{2\Delta}{\tilde{q}_\parallel}$ (see tab. I). Note that this expression corresponds to the explicit elimination of the $\omega \to 0$ limit to the Polarization operator
$\lim_{\omega \to 0} \Pi^{\mu\nu}(\omega) A_\nu(\omega) = 0$ for constant static $A_\nu(0)$ [20] [22]. This is the difference between using the non-regularized expression of the Polarization operator given in Eq. (63) and the regularized used given in Eq. (66). In Fig. 1 can be observed, for real ($\hbar\omega$) and imaginary ($\hbar\omega = i\hbar\xi$) frequencies the divergence of Eq. (121) and the convergence of Eq. (123) to the local limit given in Eq. (93). For imaginary frequencies (Fig. 1a) and for the imaginary part of the conductivity for real frequencies (Fig. 1c), the appearance of the Plasma-like peak is evident, for the real part of the conductivity at real frequencies, the Plasma-like peak is a Dirac delta and cannot be observed in the figure (Fig. 1b). It is interesting to note that, when $\tilde{q}_\parallel \to 0$, the divergence disappear.

Here it is clearly seen the importance the regularization term $\lim_{\omega \to 0} \Pi_{\mu\nu}(\omega, q)$ has. In the QFT-b models, the lack of regularization at zero frequency opens the possibility of obtain spurious non-physical Plasma conductivities without dissipation ($\Gamma = 0$) in graphene, note that this result is not cured even if dissipation ($\Gamma > 0$) is artificially added to the model. This fact actually happens for the transversal conductivity $\sigma_{T,0}^{QFT-b}(\xi, \tilde{q}_\parallel)$ in Eq. (121), the longitudinal conductivity is saved from this divergence because it scales with $\omega^2$ (See Eq. (78)), however, the term $\sigma_{00}^{QFT-b}$ should be regularized as well to avoid any

equivalent non-physical behavior of the charge density of the system. It is worth noting that, in the local limit, this divergence disappears, so none of the models studied here have this non-physical dissipation-less Plasma current in their corresponding local limit. Finally, the non-local Kubo model, by construction, do not have this divergence.

The appearance of this dissipation-less current associated to the constitutive equation given in Eq. (56)

$$J_{\mu}(\omega, q) = \Pi_{\mu\nu}(\omega, q) A^{\nu}(\omega, q), \quad (124)$$

resemble to the well-known London equation [23]

$$j_i(\omega, q) = \frac{-n_s e^2}{m} A_i(\omega, q). \quad (125)$$

London equation also predicts the Plasma model for the electronic conductivity without losses, but the Gauge invariance is explicitly broken. This breaking of Gauge invariance indicated the need of an additional microscopic theory of superconductivity that would restore the Gauge invariance of the theory, the BCS theory [24][98][25][22]. A detailed study of the Gauge invariance of the Kubo and QFT-b models will be published in the future [99].

<figure><img src="image_2.png" /><figcaption>Figure 1. (Color online) Double logarithmic plots of the real part of the transversal conductivity $\sigma_T$ (in units of the universal conductivity of graphene $\sigma_0 = \frac{\alpha c}{4}$) as a function of the imaginary frequency $\hbar\omega = i\hbar\xi$ in panel (a), the real part for real frequency $\hbar\omega$ in panel (b) and the imaginary part for real frequencies in (c), for the case with $k_B T = \mu = 0$, $\Delta = 2$ eV, $\tilde{q}_\parallel = 1$ eV and dissipation rate $\hbar\Gamma = 10^{-3}$ eV. The thick black curve is the non-local conductivity derived from the Kubo formula, given in this limit by Eq. (123), the red curve is the non-local conductivity $\sigma_{T,0}^{QFT-b}(\omega, \tilde{q}_\parallel)$ show in Eq. (121), the yellow curve is the local conductivity show in Eq. (93) and the green curve is the $\hbar\xi \to 0$ divergence of Eq. (121), given in Eq. (122). The dashed curves in panel (c) represent the positive imaginary parts of the conductivity, while the full curves represent the negative values.</figcaption></figure>

VIII. NUMERICAL COMPARISON

In the following figures, we compare the longitudinal and transversal conductivities derived from the 3 models we have studied here. We compared for different temperatures $T$, chemical potentials $\mu$, mass gaps $\Delta$ and momentum $\tilde{q}$, the longitudinal (Fig. 2) and transversal conductivities (Fig. 3) derived from the three different models studied here. The local limit ($\tilde{q}_\parallel = 0$) at $T = 0$ K and $T = 300$ K are represented by the thick black and the

yellow dashed curves respectively. The non-local conductivities derived from the Kubo formula with $\tilde{q}_\parallel \neq 0$ eV and $\hbar\Gamma = 10^{-3}$ eV at $T = 0$ K and $T = 300$ K are represented by the thick red and the blue dashed curves respectively. The non-local conductivity derived from the QFT-b model with $\tilde{q}_\parallel \neq 0$ eV at $T = 0$ K and $T = 300$ K are represented by the dashed brown and thick green curve respectively. We study the non-locality for $\tilde{q}_\parallel = 10^{-2}$ eV and for $\tilde{q}_\parallel = 1$ eV.

As can be observed in Fig. 2, the results for the non-
$\lim_{\omega \to 0} \Pi^{\mu\nu}(\omega) A_\nu(\omega) = 0$ for constant static $A_\nu(0)$ [20] [22]. This is the difference between using the non-regularized expression of the Polarization operator given in Eq. (63) and the regularized used given in Eq. (66). In Fig. 1 can be observed, for real ($\hbar\omega$) and imaginary ($\hbar\omega = i\hbar\xi$) frequencies the divergence of Eq. (121) and the convergence of Eq. (123) to the local limit given in Eq. (93). For imaginary frequencies (Fig. 1a) and for the imaginary part of the conductivity for real frequencies (Fig. 1c), the appearance of the Plasma-like peak is evident, for the real part of the conductivity at real frequencies, the Plasma-like peak is a Dirac delta and cannot be observed in the figure (Fig. 1b). It is interesting to note that, when $\tilde{q}_\parallel \to 0$, the divergence disappear.

Here it is clearly seen the importance the regularization term $\lim_{\omega \to 0} \Pi_{\mu\nu}(\omega, q)$ has. In the QFT-b models, the lack of regularization at zero frequency opens the possibility of obtain spurious non-physical Plasma conductivities without dissipation ($\Gamma = 0$) in graphene, note that this result is not cured even if dissipation ($\Gamma > 0$) is artificially added to the model. This fact actually happens for the transversal conductivity $\sigma_{T,0}^{QFT-b}(\xi, \tilde{q}_\parallel)$ in Eq. (121), the longitudinal conductivity is saved from this divergence because it scales with $\omega^2$ (See Eq. (78)), however, the term $\sigma_{00}^{QFT-b}$ should be regularized as well to avoid any

equivalent non-physical behavior of the charge density of the system. It is worth noting that, in the local limit, this divergence disappears, so none of the models studied here have this non-physical dissipation-less Plasma current in their corresponding local limit. Finally, the non-local Kubo model, by construction, do not have this divergence.

The appearance of this dissipation-less current associated to the constitutive equation given in Eq. (56)

$$J_{\mu}(\omega, q) = \Pi_{\mu\nu}(\omega, q) A^{\nu}(\omega, q), \quad (124)$$

resemble to the well-known London equation [23]

$$j_i(\omega, q) = \frac{-n_s e^2}{m} A_i(\omega, q). \quad (125)$$

London equation also predicts the Plasma model for the electronic conductivity without losses, but the Gauge invariance is explicitly broken. This breaking of Gauge invariance indicated the need of an additional microscopic theory of superconductivity that would restore the Gauge invariance of the theory, the BCS theory [24][98][25][22]. A detailed study of the Gauge invariance of the Kubo and QFT-b models will be published in the future [99].

<figure><img src="image_3.png" /><figcaption>Figure 1. (Color online) Double logarithmic plots of the real part of the transversal conductivity $\sigma_T$ (in units of the universal conductivity of graphene $\sigma_0 = \frac{\alpha c}{4}$) as a function of the imaginary frequency $\hbar\omega = i\hbar\xi$ in panel (a), the real part for real frequency $\hbar\omega$ in panel (b) and the imaginary part for real frequencies in (c), for the case with $k_B T = \mu = 0$, $\Delta = 2$ eV, $\tilde{q}_\parallel = 1$ eV and dissipation rate $\hbar\Gamma = 10^{-3}$ eV. The thick black curve is the non-local conductivity derived from the Kubo formula, given in this limit by Eq. (123), the red curve is the non-local conductivity $\sigma_{T,0}^{QFT-b}(\omega, \tilde{q}_\parallel)$ show in Eq. (121), the yellow curve is the local conductivity show in Eq. (93) and the green curve is the $\hbar\xi \to 0$ divergence of Eq. (121), given in Eq. (122). The dashed curves in panel (c) represent the positive imaginary parts of the conductivity, while the full curves represent the negative values.</figcaption></figure>

## VIII. NUMERICAL COMPARISON

In the following figures, we compare the longitudinal and transversal conductivities derived from the 3 models we have studied here. We compared for different temperatures $T$, chemical potentials $\mu$, mass gaps $\Delta$ and momentum $\tilde{q}$, the longitudinal (Fig. 2) and transversal conductivities (Fig. 3) derived from the three different models studied here. The local limit ($\tilde{q}_\parallel = 0$) at $T = 0$ K and $T = 300$ K are represented by the thick black and the

yellow dashed curves respectively. The non-local conductivities derived from the Kubo formula with $\tilde{q}_\parallel \neq 0$ eV and $\hbar\Gamma = 10^{-3}$ eV at $T = 0$ K and $T = 300$ K are represented by the thick red and the blue dashed curves respectively. The non-local conductivity derived from the QFT-b model with $\tilde{q}_\parallel \neq 0$ eV at $T = 0$ K and $T = 300$ K are represented by the dashed brown and thick green curve respectively. We study the non-locality for $\tilde{q}_\parallel = 10^{-2}$ eV and for $\tilde{q}_\parallel = 1$ eV.

As can be observed in Fig. 2, the results for the non-
local longitudinal conductivities derived from the Kubo formula and from the QFT-b model almost coincide. Actually, if we artificially add the (non-existent) dissipation rate $\hbar\Gamma$ to the QFT-b model, the two curves would almost superimpose, with any difference explained by the use of different regularizations of the integrals. This fact remarks the contribution the electronic quasi-particle dissipation has in the conductivity, mainly for frequencies $|\omega| \lesssim \Gamma$.

As can be observed in Figs. 3, the results for the non-local transversal conductivities derived from the Kubo formula and from the QFT-b model are very different. The main problem here is that an spurious asymptote proportional to $(\hbar\xi)^{-1}$ appears for very small imaginary frequencies. This difference is explained because $\sigma_T^{QFT-b}(q)$ is not regularized as imposed by the Kubo formula to fulfil $\lim_{\omega\to 0} \Pi^{\mu\nu}(\omega)A_\nu(\omega) = 0$.

## IX. CONCLUSIONS

In this article we have shown a detailed derivation of the polarization $\Pi_{\mu\nu}$ and conductivity $\sigma_{\mu\nu}$ tensor for graphene close to the Dirac point in the continuous limit. We have used the Kubo formula [6] ($\sigma^K$), a Quantum Field Theory based (QFT-b) model (which approaches the electronic quasiparticles as $(2+1)D$ Dirac electrons) [9] ($\sigma^{QFT-b}$) and a local model [31] ($\sigma^F$).

The more general result is obtained with the Kubo formula $\sigma^K$. This result is valid for any complex frequency (with positive imaginary part) $\omega \in \mathbb{C}^+$, constant dissipation rate $\Gamma$, chemical potential $\mu$ and Dirac mass $m$ as a closed analytical formula at zero Temperature $T$. The non-zero temperature results can be obtained after an integration by using the Maldague formula (Eq. (92)).

We obtain that the local limit of $\sigma^K$ is actually $\sigma^F$ if we make $m=0$ and the dissipation of the interband conductivity exactly zero.

We have derived the polarization (and, therefore, the conductivity) as in the QFT-b model, we obtain again the results published elsewhere, and we find that the longitudinal conductivity derived from the Kubo formula and from the QFT-b model almost coincide, with any difference explained by the different regularization strategies used. However, there is not such a coincidence with the transversal conductivity. The main difference comes from the regularization of the polarization opera-

tor used. In the case of the Kubo formula, the expression of the polarization operator is regularized by imposing $\lim_{\omega\to 0} \Pi^{\mu\nu}(\omega)A_\nu(\omega) = 0$ for all $A_\nu(\omega)$, as shown in Eq. (D16) of Appendix D. By the contrary, in the derivation of the QFT-b model, as the assumed constitutive relation (Eq. (56)) does not impose any regularization term, the Longitudinal conductivity derived from the QFT-b model only coincides with the result derived from the Kubo formula because the regularization term accidentally cancels out in this case, this is not the case of the Transversal conductivity, for which the QFT-b result fulfils $\lim_{\omega\to 0} \Pi_T(\omega) \neq 0$. As a consequence, a transversal Plasma conductivity is obtained in the QFT-b model, that would implies a dissipation-less electric current which is clearly not acceptable in normal materials.

We have shown that the use of the Ohm law as the constitutive relation between the electric current $J_\mu$ and the electric field $E^\nu$ (Eq. (51)) instead of the assumed linear relationship of the electric current with the potential vector $A^\nu$ (Eq. (56)) in the Kubo formula leads to different models for the electric conductivity of 2D materials described by the Dirac Hamiltonian, like graphene. This difference can be traced out to a regularization term that must be applied to the conductivity tensor. Once this regularization term is taken into account, the Kubo and QFT-b model exactly coincide for all parameters of the model.

This result can affect to the prediction of the Casimir effect with graphene.

## ACKNOWLEDGMENTS

P. R.-L. acknowledges support from Ministerio de Ciencia e Innovación (Spain), Agencia Estatal de Investigación, under project NAUTILUS (PID2022-139524NB-I00), from AYUDA PUENTE, URJC, from QuantUM program of the University of Montpellier and the hospitality of the Theory of Light-Matter and Quantum Phenomena group at the Laboratoire Charles Coulomb, University of Montpellier, where part of this work was done. M.A. acknowledges the QuantUM program of the University of Montpellier, the grant "CAT", No. A-HKUST604/20, from the ANR/RGC Joint Research Scheme sponsored by the French National Research Agency (ANR) and the Research Grants Council (RGC) of the Hong Kong Special Administrative Region.

## Appendix A: Non-local Kubo conductivity expressions

The analytical formulas for the non-local conductivities of 2D Dirac cones have been derived in [6] from Eq. (69). Those formulas are naturally divided into two parts, one independent of the chemical potential $\mu$ and another term for which the chemical potential is accounted for. Namely, $\sigma_p(q) = \sigma_{p,0}(q) + \Theta(|\mu| - |\Delta|)\sigma_{p,1}(q)$, where $\Theta$ is the
local longitudinal conductivities derived from the Kubo formula and from the QFT-b model almost coincide. Actually, if we artificially add the (non-existent) dissipation rate $\hbar\Gamma$ to the QFT-b model, the two curves would almost superimpose, with any difference explained by the use of different regularizations of the integrals. This fact remarks the contribution the electronic quasi-particle dissipation has in the conductivity, mainly for frequencies $|\omega| \lesssim \Gamma$.

As can be observed in Figs. 3, the results for the non-local transversal conductivities derived from the Kubo formula and from the QFT-b model are very different. The main problem here is that an spurious asymptote proportional to $(\hbar\xi)^{-1}$ appears for very small imaginary frequencies. This difference is explained because $\sigma_T^{QFT-b}(q)$ is not regularized as imposed by the Kubo formula to fulfil $\lim_{\omega\to 0} \Pi^{\mu\nu}(\omega)A_\nu(\omega) = 0$.

## IX. CONCLUSIONS

In this article we have shown a detailed derivation of the polarization $\Pi_{\mu\nu}$ and conductivity $\sigma_{\mu\nu}$ tensor for graphene close to the Dirac point in the continuous limit. We have used the Kubo formula [6] ($\sigma^K$), a Quantum Field Theory based (QFT-b) model (which approaches the electronic quasiparticles as $(2+1)D$ Dirac electrons) [9] ($\sigma^{QFT-b}$) and a local model [31] ($\sigma^F$).

The more general result is obtained with the Kubo formula $\sigma^K$. This result is valid for any complex frequency (with positive imaginary part) $\omega \in \mathbb{C}^+$, constant dissipation rate $\Gamma$, chemical potential $\mu$ and Dirac mass $m$ as a closed analytical formula at zero Temperature $T$. The non-zero temperature results can be obtained after an integration by using the Maldague formula (Eq. (92)).

We obtain that the local limit of $\sigma^K$ is actually $\sigma^F$ if we make $m=0$ and the dissipation of the interband conductivity exactly zero.

We have derived the polarization (and, therefore, the conductivity) as in the QFT-b model, we obtain again the results published elsewhere, and we find that the longitudinal conductivity derived from the Kubo formula and from the QFT-b model almost coincide, with any difference explained by the different regularization strategies used. However, there is not such a coincidence with the transversal conductivity. The main difference comes from the regularization of the polarization opera-

tor used. In the case of the Kubo formula, the expression of the polarization operator is regularized by imposing $\lim_{\omega\to 0} \Pi^{\mu\nu}(\omega)A_\nu(\omega) = 0$ for all $A_\nu(\omega)$, as shown in Eq. (D16) of Appendix D. By the contrary, in the derivation of the QFT-b model, as the assumed constitutive relation (Eq. (56)) does not impose any regularization term, the Longitudinal conductivity derived from the QFT-b model only coincides with the result derived from the Kubo formula because the regularization term accidentally cancels out in this case, this is not the case of the Transversal conductivity, for which the QFT-b result fulfils $\lim_{\omega\to 0} \Pi_T(\omega) \neq 0$. As a consequence, a transversal Plasma conductivity is obtained in the QFT-b model, that would implies a dissipation-less electric current which is clearly not acceptable in normal materials.

We have shown that the use of the Ohm law as the constitutive relation between the electric current $J_\mu$ and the electric field $E^\nu$ (Eq. (51)) instead of the assumed linear relationship of the electric current with the potential vector $A^\nu$ (Eq. (56)) in the Kubo formula leads to different models for the electric conductivity of 2D materials described by the Dirac Hamiltonian, like graphene. This difference can be traced out to a regularization term that must be applied to the conductivity tensor. Once this regularization term is taken into account, the Kubo and QFT-b model exactly coincide for all parameters of the model.

This result can affect to the prediction of the Casimir effect with graphene.

## ACKNOWLEDGMENTS

P. R.-L. acknowledges support from Ministerio de Ciencia e Innovación (Spain), Agencia Estatal de Investigación, under project NAUTILUS (PID2022-139524NB-I00), from AYUDA PUENTE, URJC, from QuantUM program of the University of Montpellier and the hospitality of the Theory of Light-Matter and Quantum Phenomena group at the Laboratoire Charles Coulomb, University of Montpellier, where part of this work was done. M.A. acknowledges the QuantUM program of the University of Montpellier, the grant "CAT", No. A-HKUST604/20, from the ANR/RGC Joint Research Scheme sponsored by the French National Research Agency (ANR) and the Research Grants Council (RGC) of the Hong Kong Special Administrative Region.

## Appendix A: Non-local Kubo conductivity expressions

The analytical formulas for the non-local conductivities of 2D Dirac cones have been derived in [6] from Eq. (69). Those formulas are naturally divided into two parts, one independent of the chemical potential $\mu$ and another term for which the chemical potential is accounted for. Namely, $\sigma_p(q) = \sigma_{p,0}(q) + \Theta(|\mu| - |\Delta|)\sigma_{p,1}(q)$, where $\Theta$ is the
<figure><img src="image_3.png" /><figcaption>Figure 2. (Color online) Double logarithmic plots of the longitudinal conductivity $\sigma_L(i\hbar\xi)$ (in units of the universal conductivity of graphene $\sigma_0 = \frac{\alpha c}{4}$) as a function of the imaginary frequency $\hbar\omega = i\hbar\xi$. The thick black and the yellow dashed curves are the local conductivity ($\tilde{q}_\parallel = 0$) at $T = 0$ K and $T = 300$ K respectively. The thick red and the blue dashed curves are the non-local conductivities derived from the Kubo formula with $\tilde{q}_\parallel \neq 0$, $\hbar\Gamma = 10^{-3}$ eV at $T = 0$ K and $T = 300$ K respectively. The thick green and brown dashed curves are the non-local conductivity derived from the QFT-b model with $\tilde{q}_\parallel \neq 0$ at $T = 300$ K and $T = 0$ K respectively, finally, the thin gray line is the universal conductivity of graphene with $\sigma(i\hbar\xi) = \sigma_0 = \frac{\alpha c}{4}$. In a) $\{\mu, \Delta, \tilde{q}_\parallel\} = \{0, 0, 10^{-2}\}$ eV, in b) $\{\mu, \Delta, \tilde{q}_\parallel\} = \{0, 0.2, 10^{-2}\}$ eV, in c) $\{\mu, \Delta, \tilde{q}_\parallel\} = \{0.25, 0, 10^{-2}\}$ eV, in d) $\{\mu, \Delta, \tilde{q}_\parallel\} = \{0.25, 0.2, 10^{-2}\}$ eV, in e) $\{\mu, \Delta, \tilde{q}_\parallel\} = \{0, 0, 1\}$ eV, in f) $\{\mu, \Delta, \tilde{q}_\parallel\} = \{0, 0.2, 1\}$ eV, in g) $\{\mu, \Delta, \tilde{q}_\parallel\} = \{0.25, 0, 1\}$ eV, and in h) $\{\mu, \Delta, \tilde{q}_\parallel\} = \{0.25, 0.2, 1\}$ eV.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3. (Color online) Double logarithmic plots of the transversal conductivity $\sigma_T(i\hbar\xi)$ (in units of the universal conductivity of graphene $\sigma_0 = \frac{\alpha c}{4}$) as a function of the imaginary frequency $\hbar\omega = i\hbar\xi$. The thick black and the yellow dashed curves are the local conductivity ($\tilde{q}_\parallel = 0$) at $T = 0$ K and $T = 300$ K respectively. The thick red and the blue dashed curves are the non-local conductivities derived from the Kubo formula with $\tilde{q}_\parallel \neq 0$, $\hbar\Gamma = 10^{-3}$ eV at $T = 0$ K and $T = 300$ K respectively. The thick green and brown dashed curves are the non-local conductivity derived from the QFT-b model with $\tilde{q}_\parallel \neq 0$ at $T = 300$ K and $T = 0$ K respectively, finally, the thin gray line is the universal conductivity of graphene with $\sigma(i\hbar\xi) = \sigma_0 = \frac{\alpha c}{4}$. In a) $\{\mu, \Delta, \tilde{q}_\parallel\} = \{0, 0, 10^{-2}\}$ eV, in b) $\{\mu, \Delta, \tilde{q}_\parallel\} = \{0, 0.2, 10^{-2}\}$ eV, in c) $\{\mu, \Delta, \tilde{q}_\parallel\} = \{0.25, 0, 10^{-2}\}$ eV, in d) $\{\mu, \Delta, \tilde{q}_\parallel\} = \{0.25, 0.2, 10^{-2}\}$ eV, in e) $\{\mu, \Delta, \tilde{q}_\parallel\} = \{0, 0, 1\}$ eV, in f) $\{\mu, \Delta, \tilde{q}_\parallel\} = \{0, 0.2, 1\}$ eV, in g) $\{\mu, \Delta, \tilde{q}_\parallel\} = \{0.25, 0, 1\}$ eV, and in h) $\{\mu, \Delta, \tilde{q}_\parallel\} = \{0.25, 0.2, 1\}$ eV.</figcaption></figure>

Heaviside step function and $p = \{L, T, H\}$. Using the parameters $\tilde{\theta}_z = \sqrt{\tilde{q}_\parallel^2 - \tilde{q}_0^2}$, $\gamma = \frac{\Xi}{\tilde{\theta}_z}$, $\delta = \frac{2|\Delta|}{\tilde{\theta}_z}$, and

$$ \Psi(x) = 2 \left[ x + (1 - x^2) \tan^{-1} \left( \frac{1}{x} \right) \right], \qquad (A1) $$
<figure><img src="image_4.png" /><figcaption>Figure 2. (Color online) Double logarithmic plots of the longitudinal conductivity $\sigma_L(i\hbar\xi)$ (in units of the universal conductivity of graphene $\sigma_0 = \frac{\alpha c}{4}$) as a function of the imaginary frequency $\hbar\omega = i\hbar\xi$. The thick black and the yellow dashed curves are the local conductivity ($\tilde{q}_{\parallel} = 0$) at $T = 0$ K and $T = 300$ K respectively. The thick red and the blue dashed curves are the non-local conductivities derived from the Kubo formula with $\tilde{q}_{\parallel} \neq 0$, $\hbar\Gamma = 10^{-3}$ eV at $T = 0$ K and $T = 300$ K respectively. The thick green and brown dashed curves are the non-local conductivity derived from the QFT-b model with $\tilde{q}_{\parallel} \neq 0$ at $T = 300$ K and $T = 0$ K respectively, finally, the thin gray line is the universal conductivity of graphene with $\sigma(i\hbar\xi) = \sigma_0 = \frac{\alpha c}{4}$. In a) $\{\mu, \Delta, \tilde{q}_{\parallel}\} = \{0, 0, 10^{-2}\}$ eV, in b) $\{\mu, \Delta, \tilde{q}_{\parallel}\} = \{0, 0.2, 10^{-2}\}$ eV, in c) $\{\mu, \Delta, \tilde{q}_{\parallel}\} = \{0.25, 0, 10^{-2}\}$ eV, in d) $\{\mu, \Delta, \tilde{q}_{\parallel}\} = \{0.25, 0.2, 10^{-2}\}$ eV, in e) $\{\mu, \Delta, \tilde{q}_{\parallel}\} = \{0, 0, 1\}$ eV, in f) $\{\mu, \Delta, \tilde{q}_{\parallel}\} = \{0, 0.2, 1\}$ eV, in g) $\{\mu, \Delta, \tilde{q}_{\parallel}\} = \{0.25, 0, 1\}$ eV, and in h) $\{\mu, \Delta, \tilde{q}_{\parallel}\} = \{0.25, 0.2, 1\}$ eV.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 3. (Color online) Double logarithmic plots of the transversal conductivity $\sigma_T(i\hbar\xi)$ (in units of the universal conductivity of graphene $\sigma_0 = \frac{\alpha c}{4}$) as a function of the imaginary frequency $\hbar\omega = i\hbar\xi$. The thick black and the yellow dashed curves are the local conductivity ($\tilde{q}_{\parallel} = 0$) at $T = 0$ K and $T = 300$ K respectively. The thick red and the blue dashed curves are the non-local conductivities derived from the Kubo formula with $\tilde{q}_{\parallel} \neq 0$, $\hbar\Gamma = 10^{-3}$ eV at $T = 0$ K and $T = 300$ K respectively. The thick green and brown dashed curves are the non-local conductivity derived from the QFT-b model with $\tilde{q}_{\parallel} \neq 0$ at $T = 300$ K and $T = 0$ K respectively, finally, the thin gray line is the universal conductivity of graphene with $\sigma(i\hbar\xi) = \sigma_0 = \frac{\alpha c}{4}$. In a) $\{\mu, \Delta, \tilde{q}_{\parallel}\} = \{0, 0, 10^{-2}\}$ eV, in b) $\{\mu, \Delta, \tilde{q}_{\parallel}\} = \{0, 0.2, 10^{-2}\}$ eV, in c) $\{\mu, \Delta, \tilde{q}_{\parallel}\} = \{0.25, 0, 10^{-2}\}$ eV, in d) $\{\mu, \Delta, \tilde{q}_{\parallel}\} = \{0.25, 0.2, 10^{-2}\}$ eV, in e) $\{\mu, \Delta, \tilde{q}_{\parallel}\} = \{0, 0, 1\}$ eV, in f) $\{\mu, \Delta, \tilde{q}_{\parallel}\} = \{0, 0.2, 1\}$ eV, in g) $\{\mu, \Delta, \tilde{q}_{\parallel}\} = \{0.25, 0, 1\}$ eV, and in h) $\{\mu, \Delta, \tilde{q}_{\parallel}\} = \{0.25, 0.2, 1\}$ eV.</figcaption></figure>

Heaviside step function and $p = \{L, T, H\}$. Using the parameters $\tilde{\theta}_z = \sqrt{\tilde{q}_{\parallel}^2 - \tilde{q}_0^2}$, $\gamma = \frac{\Xi}{\tilde{\theta}_z}$, $\delta = \frac{2|\Delta|}{\tilde{\theta}_z}$, and

$$ \Psi(x) = 2 \left[ x + (1 - x^2) \tan^{-1} \left( \frac{1}{x} \right) \right], \qquad (A1) $$
$$ \sigma_{L,0}(q) = \frac{\sigma_0}{2\pi} \frac{-i\tilde{q}_0}{\tilde{\theta}_z^2} \left[ 2|\Delta| + \frac{\tilde{\theta}_z^2 - 4\Delta^2}{\tilde{\theta}_z} \tan^{-1} \left( \frac{\tilde{\theta}_z}{2|\Delta|} \right) \right] = \frac{\sigma_0}{4\pi} \frac{-i\tilde{q}_0}{\tilde{\theta}_z} \Psi(\delta), \quad (A2) $$

$$ \sigma_{L,1}(q) = \frac{\sigma_0}{2\pi} \frac{-i\tilde{q}_0}{\tilde{q}_{\parallel}^2} \left[ 4(|\mu| - |\Delta|) + \frac{1}{2\tilde{\theta}_z} \left( \mathcal{F}_1 + (R^2 - \tilde{q}_{\parallel}^2) \mathcal{F}_2 \right) + \frac{\tilde{q}_{\parallel}}{q_0^2} \Theta(\tilde{q}_{\parallel}^2 - 4(\mu^2 - \Delta^2)) \mathcal{F}_3 \right], \quad (A3) $$

$$ \sigma_{T,0}(q) = \frac{\sigma_0}{2\pi} \frac{1}{-i\tilde{q}_0} \left[ \frac{\tilde{\theta}_z^2 - 4\Delta^2}{\tilde{\theta}_z} \tan^{-1} \left( \frac{\tilde{\theta}_z}{2|\Delta|} \right) - \frac{\tilde{q}_{\parallel}^2 - 4\Delta^2}{\tilde{q}_{\parallel}} \tan^{-1} \left( \frac{\tilde{q}_{\parallel}}{2|\Delta|} \right) \right] = \frac{\sigma_0}{4\pi} \frac{1}{-i\tilde{q}_0} \left[ \tilde{\theta}_z \Psi(\delta) - \tilde{q}_{\parallel} \Psi(x) \right], \quad (A4) $$

$$ \sigma_{T,1}(q) = \frac{\sigma_0}{2\pi} \frac{i}{q_0} \left[ 4 \frac{\tilde{q}_0^2}{\tilde{q}_{\parallel}^2} 4(|\mu| - |\Delta|) + 2|\Delta| - \frac{1}{2\tilde{\theta}_z} \left( \frac{\tilde{\theta}_z^2}{\tilde{q}_{\parallel}^2} \mathcal{F}_1 + (\tilde{\theta}_z^2 - 4\Delta^2) \mathcal{F}_2 \right) \right. \\ \left. + \frac{1}{\tilde{q}_{\parallel}} \left( \mathcal{F}_4 \Theta(\tilde{q}_{\parallel}^2 - 4(\mu^2 - \Delta^2)) + \mathcal{F}_5 \Theta(4(\mu^2 - \Delta^2) - \tilde{q}_{\parallel}^2) \right) \right], \quad (A5) $$

$$ \sigma_{H,0}(q) = \frac{2\sigma_0}{\pi} \frac{\eta\Delta}{\tilde{\theta}_z} \tan^{-1} \left( \frac{\tilde{\theta}_z}{2|\Delta|} \right), \quad (A6) $$

$$ \sigma_{H,1}(q) = -\frac{\sigma_0}{\pi} \frac{\eta\Delta}{\tilde{\theta}_z} \left[ \tan^{-1} \left( \frac{\tilde{q}_0 - 2|\Delta|}{\sqrt{R^2 - (\tilde{q}_0 - 2|\Delta|)^2}} \right) - \tan^{-1} \left( \frac{\tilde{q}_0 - 2|\mu|}{\sqrt{R^2 - (\tilde{q}_0 - 2|\mu|)^2}} \right) \right. \\ \left. + i \log \left( \frac{\tilde{q}_0 + 2|\mu| + \sqrt{(\tilde{q}_0 + 2|\mu|)^2 - R^2}}{\tilde{q}_0 + 2|\Delta| + \sqrt{(\tilde{q}_0 + 2|\Delta|)^2 - R^2}} \right) \right], \quad (A7) $$

where $\sigma_0 = \frac{\alpha c}{4} = \frac{e^2}{4\hbar}$, $\tilde{q}_0 = \hbar\Omega = \hbar\omega + i\hbar\Gamma$ ($\Gamma = \tau^{-1}$ accounts for the relaxation time), $\tilde{q}_{\parallel} = \hbar v_F q$, $\tilde{\theta}_z = \sqrt{\tilde{q}_{\parallel}^2 - \tilde{q}_0^2}$, $R = \tilde{q}_{\parallel} \sqrt{1 + 4(|\Delta|/\tilde{\theta}_z)^2}$, $\delta = \frac{2\Delta}{\tilde{\theta}_z}$ and $x = \frac{2\Delta}{\tilde{q}_{\parallel}}$. It is important to note that only the modulus of the mass gaps enter into the expressions for $\sigma_L$ and $\sigma_T$, while $\sigma_H$ has an additional dependency on the sign of the gaps through the combination $\eta\Delta_s^\eta$. The auxiliary functions $\{\mathcal{F}_n\}_{n=1}^5$ depend on $q_0$, $\mu$, $|\Delta|$, $\tilde{q}_{\parallel}$, and $R$, and are the following

$$ \mathcal{F}_1 = (\tilde{q}_0 - 2|\mu|) \frac{\sqrt{R^2 - (\tilde{q}_0 - 2|\mu|)^2}}{\sqrt{R^2 - (\tilde{q}_0 - 2|\Delta|)^2}} - (\tilde{q}_0 - 2|\Delta|) \frac{\sqrt{R^2 - (\tilde{q}_0 - 2|\Delta|)^2}}{\sqrt{R^2 - (\tilde{q}_0 - 2|\mu|)^2}} \\ + i(\tilde{q}_0 + 2|\mu|) \frac{\sqrt{(\tilde{q}_0 + 2|\mu|)^2 - R^2}}{\sqrt{(\tilde{q}_0 + 2|\Delta|)^2 - R^2}} - i(\tilde{q}_0 + 2|\Delta|) \frac{\sqrt{(\tilde{q}_0 + 2|\Delta|)^2 - R^2}}{\sqrt{(\tilde{q}_0 + 2|\mu|)^2 - R^2}}, \quad (A8) $$

$$ \mathcal{F}_2 = \tan^{-1} \left( \frac{\tilde{q}_0 - 2|\Delta|}{\sqrt{R^2 - (\tilde{q}_0 - 2|\Delta|)^2}} \right) - \tan^{-1} \left( \frac{\tilde{q}_0 - 2|\mu|}{\sqrt{R^2 - (\tilde{q}_0 - 2|\mu|)^2}} \right) \\ -i \log \left( \tilde{q}_0 + 2|\Delta| + \sqrt{(\tilde{q}_0 + 2|\Delta|)^2 - R^2} \right) + i \log \left( \tilde{q}_0 + 2|\mu| + \sqrt{(\tilde{q}_0 + 2|\mu|)^2 - R^2} \right), \quad (A9) $$

$$ \mathcal{F}_3 = 2\mu \left( \sqrt{\tilde{q}_{\parallel}^2 - 4(\mu^2 - \Delta^2)} + i\sqrt{4(\mu^2 - \Delta^2) - \tilde{q}_{\parallel}^2} \right) \\ + (4\Delta^2 - \tilde{q}_{\parallel}^2) \left[ i \log(2|\Delta| + i\tilde{q}_{\parallel}) - i \log(2|\mu| + \sqrt{4(\mu^2 - \Delta^2) - \tilde{q}_{\parallel}^2}) \right] \\ + (4\Delta^2 - \tilde{q}_{\parallel}^2) \left[ \tan^{-1} \left( \frac{2|\mu|}{\sqrt{\tilde{q}_{\parallel}^2 - 4(\mu^2 - \Delta^2)}} \right) - \tan^{-1} \left( \frac{2|\Delta|}{\tilde{q}_{\parallel}} \right) \right], \quad (A10) $$

$$ \mathcal{F}_4 = -2|\mu|\sqrt{\tilde{q}_{\parallel}^2 - 4(\mu^2 - \Delta^2)} + (\tilde{q}_{\parallel}^2 - 4\Delta^2) \left[ \tan^{-1} \left( \frac{2|\mu|}{\sqrt{\tilde{q}_{\parallel}^2 - 4(\mu^2 - \Delta^2)}} \right) - \tan^{-1} \left( \frac{2|\Delta|}{\tilde{q}_{\parallel}} \right) \right], \quad (A11) $$

$$ \mathcal{F}_5 = (\tilde{q}_{\parallel}^2 - 4\Delta^2) \left[ \frac{\pi}{2} - \tan^{-1} \left( \frac{2|\Delta|}{\tilde{q}_{\parallel}} \right) \right]. \quad (A12) $$

Those results are valid for all frequencies $\omega \in \mathbb{C}^+$, and we remind that they are the conductivity per Dirac cone, to obtain the conductivity of Graphene, we must to sum the contribution of each 4 Dirac cones, having into account their respective (signed) Dirac masses by using

$$ \sigma_p(q) = \sum_{\eta=\pm} \sum_{s=\pm} \sigma_p(q, \Delta_s^\eta). \quad (A13) $$
$$ \sigma_{L,0}(q) = \frac{\sigma_0}{2\pi} \frac{-i\tilde{q}_0}{\tilde{\theta}_z^2} \left[ 2|\Delta| + \frac{\tilde{\theta}_z^2 - 4\Delta^2}{\tilde{\theta}_z} \tan^{-1} \left( \frac{\tilde{\theta}_z}{2|\Delta|} \right) \right] = \frac{\sigma_0}{4\pi} \frac{-i\tilde{q}_0}{\tilde{\theta}_z} \Psi(\delta), \quad (A2) $$

$$ \sigma_{L,1}(q) = \frac{\sigma_0}{2\pi} \frac{-i\tilde{q}_0}{\tilde{q}_{\parallel}^2} \left[ 4(|\mu| - |\Delta|) + \frac{1}{2\tilde{\theta}_z} \left( \mathcal{F}_1 + (R^2 - \tilde{q}_{\parallel}^2) \mathcal{F}_2 \right) + \frac{\tilde{q}_{\parallel}}{q_0^2} \Theta(\tilde{q}_{\parallel}^2 - 4(\mu^2 - \Delta^2)) \mathcal{F}_3 \right], \quad (A3) $$

$$ \sigma_{T,0}(q) = \frac{\sigma_0}{2\pi} \frac{1}{-i\tilde{q}_0} \left[ \frac{\tilde{\theta}_z^2 - 4\Delta^2}{\tilde{\theta}_z} \tan^{-1} \left( \frac{\tilde{\theta}_z}{2|\Delta|} \right) - \frac{\tilde{q}_{\parallel}^2 - 4\Delta^2}{\tilde{q}_{\parallel}} \tan^{-1} \left( \frac{\tilde{q}_{\parallel}}{2|\Delta|} \right) \right] = \frac{\sigma_0}{4\pi} \frac{1}{-i\tilde{q}_0} \left[ \tilde{\theta}_z \Psi(\delta) - \tilde{q}_{\parallel} \Psi(x) \right], \quad (A4) $$

$$ \sigma_{T,1}(q) = \frac{\sigma_0}{2\pi} \frac{i}{q_0} \left[ 4 \frac{\tilde{q}_0^2}{\tilde{q}_{\parallel}^2} (4(|\mu| - |\Delta|) + 2|\Delta|) - \frac{1}{2\tilde{\theta}_z} \left( \frac{\tilde{\theta}_z^2}{\tilde{q}_{\parallel}^2} \mathcal{F}_1 + (\tilde{\theta}_z^2 - 4\Delta^2) \mathcal{F}_2 \right) \right. \\ \left. + \frac{1}{\tilde{q}_{\parallel}} \left( \mathcal{F}_4 \Theta(\tilde{q}_{\parallel}^2 - 4(\mu^2 - \Delta^2)) + \mathcal{F}_5 \Theta(4(\mu^2 - \Delta^2) - \tilde{q}_{\parallel}^2) \right) \right], \quad (A5) $$

$$ \sigma_{H,0}(q) = \frac{2\sigma_0}{\pi} \frac{\eta\Delta}{\tilde{\theta}_z} \tan^{-1} \left( \frac{\tilde{\theta}_z}{2|\Delta|} \right), \quad (A6) $$

$$ \sigma_{H,1}(q) = -\frac{\sigma_0}{\pi} \frac{\eta\Delta}{\tilde{\theta}_z} \left[ \tan^{-1} \left( \frac{\tilde{q}_0 - 2|\Delta|}{\sqrt{R^2 - (\tilde{q}_0 - 2|\Delta|)^2}} \right) - \tan^{-1} \left( \frac{\tilde{q}_0 - 2|\mu|}{\sqrt{R^2 - (\tilde{q}_0 - 2|\mu|)^2}} \right) \right. \\ \left. + i \log \left( \frac{\tilde{q}_0 + 2|\mu| + \sqrt{(\tilde{q}_0 + 2|\mu|)^2 - R^2}}{\tilde{q}_0 + 2|\Delta| + \sqrt{(\tilde{q}_0 + 2|\Delta|)^2 - R^2}} \right) \right], \quad (A7) $$

where $\sigma_0 = \frac{\alpha c}{4} = \frac{e^2}{4\hbar}$, $\tilde{q}_0 = \hbar\Omega = \hbar\omega + i\hbar\Gamma$ ($\Gamma = \tau^{-1}$ accounts for the relaxation time), $\tilde{q}_{\parallel} = \hbar v_F q$, $\tilde{\theta}_z = \sqrt{\tilde{q}_{\parallel}^2 - \tilde{q}_0^2}$, $R = \tilde{q}_{\parallel} \sqrt{1 + 4(|\Delta|/\tilde{\theta}_z)^2}$, $\delta = \frac{2\Delta}{\tilde{\theta}_z}$ and $x = \frac{2\Delta}{\tilde{q}_{\parallel}}$. It is important to note that only the modulus of the mass gaps enter into the expressions for $\sigma_L$ and $\sigma_T$, while $\sigma_H$ has an additional dependency on the sign of the gaps through the combination $\eta\Delta_s^\eta$. The auxiliary functions $\{\mathcal{F}_n\}_{n=1}^5$ depend on $q_0$, $\mu$, $|\Delta|$, $\tilde{q}_{\parallel}$, and $R$, and are the following

$$ \mathcal{F}_1 = (\tilde{q}_0 - 2|\mu|) \frac{\sqrt{R^2 - (\tilde{q}_0 - 2|\mu|)^2}}{\sqrt{R^2 - (\tilde{q}_0 - 2|\Delta|)^2}} - (\tilde{q}_0 - 2|\Delta|) \frac{\sqrt{R^2 - (\tilde{q}_0 - 2|\Delta|)^2}}{\sqrt{R^2 - (\tilde{q}_0 - 2|\mu|)^2}} \\ + i(\tilde{q}_0 + 2|\mu|) \frac{\sqrt{(\tilde{q}_0 + 2|\mu|)^2 - R^2}}{\sqrt{(\tilde{q}_0 + 2|\Delta|)^2 - R^2}} - i(\tilde{q}_0 + 2|\Delta|) \frac{\sqrt{(\tilde{q}_0 + 2|\Delta|)^2 - R^2}}{\sqrt{(\tilde{q}_0 + 2|\mu|)^2 - R^2}}, \quad (A8) $$

$$ \mathcal{F}_2 = \tan^{-1} \left( \frac{\tilde{q}_0 - 2|\Delta|}{\sqrt{R^2 - (\tilde{q}_0 - 2|\Delta|)^2}} \right) - \tan^{-1} \left( \frac{\tilde{q}_0 - 2|\mu|}{\sqrt{R^2 - (\tilde{q}_0 - 2|\mu|)^2}} \right) \\ -i \log \left( \tilde{q}_0 + 2|\Delta| + \sqrt{(\tilde{q}_0 + 2|\Delta|)^2 - R^2} \right) + i \log \left( \tilde{q}_0 + 2|\mu| + \sqrt{(\tilde{q}_0 + 2|\mu|)^2 - R^2} \right), \quad (A9) $$

$$ \mathcal{F}_3 = 2\mu \left( \sqrt{\tilde{q}_{\parallel}^2 - 4(\mu^2 - \Delta^2)} + i\sqrt{4(\mu^2 - \Delta^2) - \tilde{q}_{\parallel}^2} \right) \\ + (4\Delta^2 - \tilde{q}_{\parallel}^2) \left[ i \log(2|\Delta| + i\tilde{q}_{\parallel}) - i \log(2|\mu| + \sqrt{4(\mu^2 - \Delta^2) - \tilde{q}_{\parallel}^2}) \right] \\ + (4\Delta^2 - \tilde{q}_{\parallel}^2) \left[ \tan^{-1} \left( \frac{2|\mu|}{\sqrt{\tilde{q}_{\parallel}^2 - 4(\mu^2 - \Delta^2)}} \right) - \tan^{-1} \left( \frac{2|\Delta|}{\tilde{q}_{\parallel}} \right) \right], \quad (A10) $$

$$ \mathcal{F}_4 = -2|\mu| \sqrt{\tilde{q}_{\parallel}^2 - 4(\mu^2 - \Delta^2)} + (\tilde{q}_{\parallel}^2 - 4\Delta^2) \left[ \tan^{-1} \left( \frac{2|\mu|}{\sqrt{\tilde{q}_{\parallel}^2 - 4(\mu^2 - \Delta^2)}} \right) - \tan^{-1} \left( \frac{2|\Delta|}{\tilde{q}_{\parallel}} \right) \right], \quad (A11) $$

$$ \mathcal{F}_5 = (\tilde{q}_{\parallel}^2 - 4\Delta^2) \left[ \frac{\pi}{2} - \tan^{-1} \left( \frac{2|\Delta|}{\tilde{q}_{\parallel}} \right) \right]. \quad (A12) $$

Those results are valid for all frequencies $\omega \in \mathbb{C}^+$, and we remind that they are the conductivity per Dirac cone, to obtain the conductivity of Graphene, we must to sum the contribution of each 4 Dirac cones, having into account their respective (signed) Dirac masses by using

$$ \sigma_p(q) = \sum_{\eta=\pm} \sum_{s=\pm} \sigma_p(q, \Delta_s^\eta). \quad (A13) $$
Finally, note that if we are interested on the $\sigma_{00}$ term of the conductivity tensor, instead of Eq. (76), we have to use

$$ \sigma^{00}(q) = \frac{q_{\parallel}^2}{q_0^2} \left[ \sigma_L(q) - \lim_{q_0 \to 0} \sigma_L(q) \right], \quad (A14) $$

to avoid the appearance of any possible spurious term in the $\omega \to 0$ limit.

## Appendix B: Local Kubo conductivity expressions

There is a special case in the local limit (when we apply the $q_{\parallel} \to 0$ limit to the Kubo formula (69)) when results valid for all temperatures can be obtained. The local limit of the conductivities of a massive Dirac cone are

$$ \begin{aligned} \sigma_{xx}(\omega, \mathbf{0}, \mu, \Delta) &= i \frac{\sigma_0}{\pi} \left[ \frac{\mu^2 - \Delta^2}{|\mu|} \frac{1}{\Omega} \Theta(|\mu| - |\Delta|) + \frac{\Delta^2}{M\Omega} - \frac{\Omega^2 + 4\Delta^2}{2i\Omega^2} \tan^{-1} \left( \frac{i\Omega}{2M} \right) \right], \\ \sigma_{xy}(\omega, \mathbf{0}, \mu, \Delta) &= \frac{2\sigma_0}{\pi} \frac{\eta\Delta}{i\Omega} \tan^{-1} \left( \frac{i\Omega}{2M} \right), \end{aligned} \quad (B1) $$

where $\Omega = \hbar\omega + i\hbar\Gamma$, $\Gamma = \tau^{-1}$ and $M = \text{Max}[|\Delta|, |\mu|]$. These results are per Dirac cone and they are consistent with the ones found by other researchers [28][31][32][33][47][87][88][89]. The first term in $\sigma_{xx}$ corresponds to inter-band transitions, and the last two terms to inter-band transitions. Note that, in the local limit $q_{\parallel} \to 0$ one obtains $\sigma_{xx}(\omega, 0) = \sigma_{yy}(\omega, 0) = \sigma_L(\omega, 0) = \sigma_T(\omega, 0)$, and $\sigma_{xy}(\omega, 0) = -\sigma_{yx}(\omega, 0) = \sigma_H(\omega, 0)$.

From those results, we obtain the conductivity at zero temperature of a Dirac cone with $\Delta = 0$ mass gap as

$$ \begin{aligned} \sigma_{xx}(\omega, \mathbf{q}_{\parallel} = \mathbf{0}, T = 0, \Delta = 0) &= \sigma_{xx}^{\text{intra}}(\omega) + \sigma_{xx}^{\text{inter}}(\omega), \\ \sigma_{xx}^{\text{intra}}(\omega) &= \frac{\alpha c}{\pi} \frac{|\mu|}{-i\Omega}, \\ \sigma_{xx}^{\text{inter}}(\omega) &= \frac{\alpha c}{2\pi} \tan^{-1} \left( \frac{-i\Omega}{2|\mu|} \right), \\ \sigma_{xy}(\omega, \mathbf{q}_{\parallel} = \mathbf{0}, T = 0, \Delta = 0) &= 0. \end{aligned} \quad (B2) $$

From this result with $\Delta = 0$, in the non-dissipation limit for the interband term and for finite temperatures, by using the Maldague formula, the well-known result of Falkovsky is obtained as

$$ \begin{aligned} \sigma_{xx}^{F}(\omega, \Gamma, T, \mu) &= \sigma_{xx}^{F, \text{intra}}(\omega, \Gamma, T, \mu) + \sigma_{xx}^{F, \text{inter}}(\omega, 0, T, \mu), \\ \sigma_{xx}^{F, \text{intra}}(\omega, \Gamma, T, \mu) &= \frac{1}{\pi \hbar} \frac{2i \alpha c k_B T}{\omega + i \Gamma} \ln \left[ 2 \cosh \left( \frac{\mu}{2k_B T} \right) \right], \\ \sigma_{xx}^{F, \text{inter}}(\omega, \Gamma, T, \mu) &= \frac{\alpha c}{4} G \left( \left| \frac{\hbar \omega}{2} \right| \right) + i \frac{\alpha c}{4} \frac{4 \hbar \omega}{\pi} \int_0^{\infty} d\xi \frac{G(\xi) - G\left(\left| \frac{\hbar \omega}{2} \right|\right)}{(\hbar \omega)^2 - 4\xi^2}. \end{aligned} \quad (B3) $$

with

$$ G(\epsilon) = n_F(-\epsilon + \mu) - n_F(+\epsilon + \mu) = \frac{\sinh(\beta\epsilon)}{\cosh(\beta\mu) + \cosh(\beta\epsilon)}. \quad (B4) $$

The intraband term is obtained by the use of the Maldague formula (Eq. (92)) for Eq. (94), as

$$ \sigma_{xx}^{\text{intra}}(\omega, T) = \frac{\alpha c}{\pi} \frac{1}{-i\Omega} \frac{2}{\beta} \ln \left[ 2 \cosh \left( \frac{\beta\mu}{2} \right) \right] = \frac{\alpha c}{\pi \hbar} \frac{2i k_B T}{\omega + i\Gamma} \ln \left[ 2 \cosh \left( \frac{\mu}{2k_B T} \right) \right]. \quad (B5) $$

To derive the local interband conductivity at finite temperatures, we first need to apply the zeroth dissipation limit $\tau \to \infty$ ($\Gamma \to 0$) of the real part $\omega \in \mathbb{R}$ of the interband conductivity of graphene, as

$$ \lim_{\tau \to \infty} \mathrm{Re} \left[ \sigma_{xx}^{\text{inter}}(\omega) \right] = \lim_{\tau \to \infty} \frac{\alpha c}{2\pi} \tan^{-1} \left( \frac{-i\Omega}{2|\mu|} \right) = \frac{\alpha c}{4} \left[ \Theta(\hbar\omega - 2|\mu|) + \Theta(-\hbar\omega - 2|\mu|) \right], \quad (B6) $$

which coincides with the result shown in [31] in the $T \to 0$ limit for $\omega > 0$. Note that $\mathrm{Re}[\sigma_{xx}^{\text{inter}}(-\omega)] = +\mathrm{Re}[\sigma_{xx}^{\text{inter}}(\omega)]$ is an even function in $\omega$.