\exp\\{-\imath\theta_{k_{B(A)}}\\}\phi_{2}\end{array}\right)\equiv\left(\begin{array}[]{c}\chi_{+\sigma_{A(B)}}\\\
\chi_{-\sigma_{B(A)}}\end{array}\right)$ (S.18)

where $\theta_{k_{B(A)}}$ are phases for free fermion fields Peskin-Schroeder
,

$\displaystyle\phi_{i}={1\over(2\pi)^{3/2}\sqrt{N/2}}\sum_{\vec{R}^{A(B)}_{l}}\exp\\{\imath[\vec{K}_{A(B)}-\vec{q}_{{}_{A(B)}}][\vec{R}^{A(B)}_{l}-\vec{r}]\\}\psi_{\\{n_{i}\\}}(\vec{r}-\vec{R}^{A(B)}_{l})\
\ \ \ \ \ \ \ $ (S.19)

is a Bloch function; maximally-localized Wannier functions
$\psi_{\\{n_{i}\\}}$, $i=1,2$ are defined as

$\psi_{\\{n_{2}\\}}=c_{1}\psi_{\mbox{\small
p}_{z}}(\vec{r}\pm\vec{\delta}_{i})+c_{2}\psi_{\mbox{\small p}_{z}}(\vec{r}),\
\sum_{i=1}^{2}c_{i}^{2}=1;\\\ \psi_{\\{n_{1}\\}}=\psi_{\mbox{\small
p}_{z}}(\vec{r});$ (S.20)

$\vec{\delta}_{i}$ is a distance to nearest $i-$th site $\vec{\delta}_{i}$,
$i=1,2,3$; definition of wave vector $\vec{q}_{{}_{A(B)}}$ uses the reference
frame origin located in the Dirac point $K_{A(B)}$, $\psi_{\mbox{\small
p}_{z}}(\vec{r})$ is the pz-wave function of a hydrogen-like atom with an
”effective charge” 1 centered at the graphene lattice sites.

An exchange interaction term $\left(\Sigma_{rel}^{x}\right)_{AB}$ is
determined as myNPCS2013 ; myNPCS17-2014 ;
NPCS18-2015GrushevskayaKrylovGaisyonokSerow

$\displaystyle\Sigma_{rel}^{x}\left(\begin{array}[]{c}\widehat{\chi}^{\dagger}_{{}_{-\sigma_{{}_{A}}}}(\vec{r})\\\
\widehat{\chi}^{\dagger}_{\sigma_{{}_{B}}}(\vec{r})\end{array}\right)\left|0,-\sigma\right>\left|0,\sigma\right>=\left(\begin{array}[]{cc}0&\left(\Sigma_{rel}^{x}\right)_{AB}\\\
\left(\Sigma_{rel}^{x}\right)_{BA}&0\end{array}\right)\left(\begin{array}[]{c}\widehat{\chi}^{\dagger}_{-\sigma_{{}_{A}}}(\vec{r})\\\
\widehat{\chi}^{\dagger}_{\sigma_{{}_{B}}}(\vec{r})\end{array}\right)\left|0,-\sigma\right>\left|0,\sigma\right>,$
(S.27)
$\displaystyle\left(\Sigma_{rel}^{x}\right)_{AB}\widehat{\chi}^{\dagger}_{\sigma_{{}_{B}}}(\vec{r})\left|0,\sigma\right>=\sum_{i=1}^{N_{v}\,N}\int{d\vec{r}_{i}}\widehat{\chi}^{\dagger}_{\sigma_{i}{{}^{B}}}(\vec{r})\left|0,\sigma\right>\langle
0,-\sigma_{i}|{\widehat{\chi}}^{\dagger}_{-\sigma_{i}^{A}}(\vec{r}_{i})V(\vec{r}_{i}-\vec{r}){\widehat{\chi}}_{-\sigma_{B}}(\vec{r}_{i})|0,-\sigma_{i^{\prime}}\rangle,$
(S.28)
$\displaystyle\left(\Sigma_{rel}^{x}\right)_{BA}\widehat{\chi}^{\dagger}_{{}_{-\sigma_{{}_{A}}}}(\vec{r})\left|0,-\sigma\right>=\sum_{i^{\prime}=1}^{N_{v}\,N}\int{d\vec{r}_{i^{\prime}}}\widehat{\chi}^{\dagger}_{{}_{-\sigma_{i^{\prime}}^{A}}}(\vec{r})\left|0,-\sigma\right>\langle
0,\sigma_{i^{\prime}}|{\widehat{\chi}}^{\dagger}_{\sigma_{i^{\prime}}^{B}}(\vec{r}_{i^{\prime}})V(\vec{r}_{i^{\prime}}-\vec{r}){\widehat{\chi}}_{{}_{\sigma_{A}}}(\vec{r}_{i^{\prime}})|0,\sigma_{i}\rangle.$
(S.29)

In the approximation of free $\pi$-electrons, all wave functions entering the
expression (S.28, S.29) are

$\begin{split}\widehat{\chi}^{\dagger}_{-\sigma_{i^{\prime}}{{}^{A}}}(\vec{r})\equiv\widehat{\chi^{\dagger}_{-\sigma_{{}_{A}}}}(\vec{r}),\
\
\widehat{\chi}^{\dagger}_{\sigma_{i}{{}^{B}}}(\vec{r})\equiv\widehat{\chi^{\dagger}_{\sigma_{{}_{B}}}}(\vec{r})\mbox{
for }\forall\ i,\ i^{\prime}.\end{split}$ (S.30)

Then, in this approximation we can write matrices
$\left(\Sigma_{rel}^{x}\right)_{AB}\approx\Sigma_{AB}$ and
$\left(\Sigma_{rel}^{x}\right)_{BA}\approx\Sigma_{BA}$ without self-action as

$\displaystyle\left(\Sigma_{rel}^{x}\right)_{AB}\widehat{\chi^{\dagger}_{+\sigma_{{}_{B}}}}(\vec{r})\left|0,\sigma\right>\approx\Sigma_{AB}\
\chi_{\sigma_{B}}=\sum_{i=1}^{N_{v}N-1}\int d\vec{r}_{i}$ $\displaystyle\times
V(\vec{r}_{i}-\vec{r})\left[\chi_{-\sigma_{A}}(\vec{r}_{i})\cdot\chi_{-\sigma_{B}}^{*}(\vec{r}_{i})\right]\chi_{+\sigma_{B}}(\vec{r})={1\over
2^{3/2}}\sum_{i=1}^{N_{v}N-1}\int d\vec{r}_{i}V(\vec{r}_{i}-\vec{r})$
$\displaystyle\times\left[\begin{array}[]{cc}e^{-\imath\theta_{K_{A}}}\phi_{1}(\vec{r}_{i})e^{\imath\theta_{K_{B}}}\phi^{*}_{1}(\vec{r}_{i})&e^{-\imath\theta_{K_{A}}}\phi_{1}(\vec{r}_{i})e^{\imath\theta_{K_{B}}}\phi^{*}_{2}(\vec{r}_{i})\\\
e^{-\imath\theta_{K_{A}}}\phi_{2}(\vec{r}_{i})e^{\imath\theta_{K_{B}}}\phi^{*}_{1}(\vec{r}_{i})&e^{-\imath\theta_{K_{A}}}\phi_{2}(\vec{r}_{i})e^{\imath\theta_{K_{B}}}\phi^{*}_{2}(\vec{r}_{i})\end{array}\right]\left[\begin{array}[]{c}-e^{-\imath[(\vec{K}_{A}-\vec{q}_{A})\cdot\vec{r}-\theta_{K_{B}}]}\phi_{2}(\vec{r})\\\
e^{-\imath[(\vec{K}_{A}-\vec{q}_{A})\cdot\vec{r}-\theta_{K_{B}}]}\phi_{1}(\vec{r})\end{array}\right],$
(S.35)
$\displaystyle\left(\Sigma_{rel}^{x}\right)_{BA}\widehat{\chi^{\dagger}_{-\sigma_{{}_{A}}}}(\vec{r})\left|0,-\sigma\right>\approx\Sigma_{BA}\
\chi_{-\sigma_{A}}=\sum_{i=1}^{N_{v}N-1}\int d\vec{r}_{i}$
$\displaystyle\times
V(\vec{r}_{i}-\vec{r})\left[\chi_{+\sigma_{B}}(\vec{r}_{i})\cdot\chi_{+\sigma_{A}}^{*}(\vec{r}_{i})\right]\chi_{-\sigma_{A}}(\vec{r})={1\over
2^{3/2}}\sum_{i=1}^{N_{v}N-1}\int d\vec{r}_{i}V(\vec{r}_{i}-\vec{r})$
$\displaystyle\times\left[\begin{array}[]{cc}-e^{\imath\theta_{K_{B}}}\phi_{2}(\vec{r}_{i})(-1)e^{-\imath\theta_{K_{A}}}\phi^{*}_{2}(\vec{r}_{i})&-e^{\imath\theta_{K_{B}}}\phi_{2}(\vec{r}_{i})e^{-\imath\theta_{K_{A}}}\phi^{*}_{1}(\vec{r}_{i})\\\
e^{\imath\theta_{K_{B}}}\phi_{1}(\vec{r}_{i})(-1)e^{-\imath\theta_{K_{A}}}\phi^{*}_{2}(\vec{r}_{i})&e^{\imath\theta_{K_{B}}}\phi_{1}(\vec{r}_{i})e^{-\imath\theta_{K_{A}}}\phi^{*}_{1}(\vec{r}_{i})\end{array}\right]\left[\begin{array}[]{c}e^{-\imath[(\vec{K}_{A}-\vec{q}_{A})\cdot\vec{r}+\theta_{K_{A}}]}\phi_{1}(\vec{r})\\\
e^{-\imath[(\vec{K}_{A}-\vec{q}_{A})\cdot\vec{r}+\theta_{K_{A}}]}\phi_{2}(\vec{r})\end{array}\right].$
(S.40)

It follows from the expressions (S.35, S.40) that the matrices $\Sigma_{AB}$
and $\Sigma_{BA}$ have the form

$\displaystyle\Sigma_{AB}={1\over
2}e^{\\{-\imath(\theta_{K_{A}}-\theta_{K_{B}})\\}}\sum_{i=1}^{N_{v}N-1}\int
d\vec{r}_{i}V(\vec{r}_{i}-\vec{r})\left[\begin{array}[]{cc}\phi_{1}(\vec{r}_{i})\phi^{*}_{1}(\vec{r}_{i})&\phi_{1}(\vec{r}_{i})\phi^{*}_{2}(\vec{r}_{i})\\\
\phi_{2}(\vec{r}_{i})\phi^{*}_{1}(\vec{r}_{i})&\phi_{2}(\vec{r}_{i})\phi^{*}_{2}(\vec{r}_{i})\end{array}\right],$
(S.43) $\displaystyle\Sigma_{BA}={1\over
2}e^{\\{-\imath(\theta_{K_{A}}-\theta_{K_{B}})\\}}\sum_{i=1}^{N_{v}N-1}\int
d\vec{r}_{i}V(\vec{r}_{i}-\vec{r})\left[\begin{array}[]{cc}\phi_{2}(\vec{r}_{i})\phi^{*}_{2}(\vec{r}_{i})&-\phi_{2}(\vec{r}_{i})\phi^{*}_{1}(\vec{r}_{i})\\\
-\phi_{1}(\vec{r}_{i})\phi^{*}_{2}(\vec{r}_{i})&\phi_{1}(\vec{r}_{i})\phi^{*}_{1}(\vec{r}_{i})\end{array}\right].$
(S.46)

The quasirelativistic Dirac–Hartree–Fock exchange interaction (S.43, S.46) in
tight-binding approximation are determined by Eqs. (III.3, III.6) of the main
text. In the direction perpendicular to the SM layer, the integration has been
limited to a distance equals two C–C bond lengths, the phases have been chosen
as $\theta_{K_{A}}=\theta_{K_{B}}$.

Simulations of the band structure have been performed for two variants of
approximation to the exchange operator. The first one already mentioned above
is the series expansion of the exchange matrices on deviation of the wave
vector from the Dirac point up to 4th order in the length of wave vector
$\vec{q}$. As it has been shown in myJModPhys2014 ; myNATO2015 ;
NPCS18-2015GrushevskayaKrylovGaisyonokSerow ; myNPCS18-2015 , this
approximation leads to small imaginary part for the energy (spectral line
width), this in fact not very bad as it effectively corresponds to a finite
decay. The approximation of the zero gauge-phases is valid in the Dirac point.

The second approximation is the use of the exchange interaction matrices
calculated based on $\pi$(pz)-orbital wave functions with full exponents and
with non-zero gauge-phases (see the detail of the approximation in
NPCS18-2015GrushevskayaKrylovGaisyonokSerow ; myNPCS18-2015 ;
myIntJModPhys2016 ; mySymmetry2016 ). The last approximation holds real
eigenenergies for all wave vectors. The condition of the real values for the
energy playing a role of the gauge condition in the second approximation has
been achieved by minimization (on respect to a given phases set) of its
imaginary part. This results into finite accuracy at the energy evaluation in
situation with non-zero phases.

Numerical simulations of eigenenergies have been performed for approximately
$3500\times 200$ grid points $(r,\ \phi)$ in 2D $\vec{q}$-space for the
purposes of subsequent estimation of charge transport properties of the
system. For simulations we use the spectral line width equals to 1 K.

The calculated Majorana model bands $\epsilon$ in fig. 3a and fig. S.1 have
the form of cones near the Dirac point $K_{A(B)}$ at wave numbers
$q<0.0002K_{A}$. Then the Dirac cone is splitted, and eight tilted sub-
replicas emerge at large $q$.

(a) (b)

(c)

Figure S.1: Energy bands of 2D semi-metal calculated in two approximations:
(a, b) in the first-order approximation with zero gauge phases in 4th order
Taylor series expansion of exponential factors of the exchange operator for Pb
monolayer (a) and graphene (b); (c) in the second-order approximations with
non-zero gauge phases in the exact expression for the exchange for graphene.
The Dirac cone and its six tilted replicas in the Dirac point is presented in
figure (b). One of six pairs Weyl-like nodes: source and sink is indicated in
figure (a).

We have calculated the Fermi velocity $v_{M}(\vec{q}_{F})$ for different wave
vectors $\vec{q}_{F}=\\{q_{F},\phi\\}$ by the formula
$v_{M}={\partial\epsilon\over\partial q_{F}}$. According to the Table 1, the
Fermi velocity $v_{M}(\vec{q}_{F})$ decreases with a factor two in respect to
its value $v_{M}(0)$ in the Dirac point.

Table 1:  The dependence of the charge carriers velocity $v_{M}(\vec{q}_{F})$ on wave vector $\vec{q}_{F}=\\{q_{F},\ \phi\\}$ at the Fermi level in polar coordinates $\\{q_{F},\ \phi\\}$ for various directions shown in Fig. S.1b. $q_{F}/K_{A}$ |  $\phi$, rad | $v_{M}(\vec{q}_{F})/v_{M}(0)$
---|---|---
0 | for all $\phi$ | 1
0.005 | 0 | 0.9985
|  $\pi$/2 | 0.9993
| $\pi$ | 0.9985
|  $\pi$/2 | 0.9968
0.02 | 2$\pi$/3 | 0.9920
| 4$\pi$/3 | 0.9135
| 2$\pi$/3 | 0.6148
0.1 |  |
| 5$\pi$/3 | 0.5012

##  II Quantum statistics of the model $N=3$

###  XIX.1 Perturbation theory

To find the quasiparticle current (IV.6) in main text, we use the perturbation
theory Dyson ; MyLambert . Power series expansion of (IV.2) allows to rewrite
the potential energy $V^{SM}$ for interaction of the secondary quantized
fermionic field $\chi_{+\sigma_{{}_{B}}}(x)$ with an electromagnetic field
(Eq. (IV.1) in main text) in the following form

$\displaystyle
V^{SM}=\chi^{\dagger}_{+\sigma_{{}_{B}}}\left[-c\vec{\sigma}_{2D}^{BA}\cdot{e\over
c}\vec{A}-\widetilde{\Sigma_{AB}\Sigma_{BA}}(0)-\sum_{i}\left.{d\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{\prime}}\right|_{p_{i}^{\prime}=0}\left(p_{i}^{AB}-{e\over
c}A_{i}\right)\right.$ $\displaystyle\left.-{1\over
2}\sum_{i,j}\left.{d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{\prime}dp_{j}^{\prime}}\right|_{p_{i}^{\prime},\
p_{j}^{\prime}=0}\left(p_{i}^{AB}-{e\over
c}A_{i}\right)\,\left(p_{j}^{AB}-{e\over
c}A_{j}\right)+\ldots\right]\chi_{+\sigma_{{}_{B}}}.$ (S.13)

Taking into account that $c\vec{\sigma}_{2D}^{BA}$ is a quantum analog of the
current vector Fock  $\vec{v}$, and expressing $\vec{p}^{AB}$ through the
dynamical mass $\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c)$,
the operator $V^{SM}$ (S.13) is rewritten as

$\displaystyle
V^{SM}=-\chi^{\dagger}_{+\sigma_{{}_{B}}}\left\\{\vec{v}\cdot{e\over
c}\vec{A}+\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c)\right.$
$\displaystyle\times\left[{\widetilde{\Sigma_{AB}\Sigma_{BA}}(0)\over\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c)}+\sum_{i}\left.{d\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{\prime}}\right|_{p_{i}^{\prime}=0}\left(v_{i}-{e\over
c}\tilde{A}_{i}\right)\right.$
$\displaystyle\left.+{\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c)\over
2}\sum_{i,j}\left.\left.{d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{\prime}dp_{j}^{\prime}}\right|_{p_{i}^{\prime},\
p_{j}^{\prime}=0}\left(v_{i}-{e\over
c}\tilde{A}_{i}\right)\,\left(v_{j}-{e\over
c}\tilde{A}_{j}\right)+\ldots\right]\right\\}\chi_{+\sigma_{{}_{B}}},$ (S.14)

where
$\vec{\tilde{A}}=\vec{A}/\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c)$.
In the interaction representation, if account for the first-order terms on
$A_{i}$ and terms quadratic on $v_{i}$ in the expansion of the evolution
operator $\hat{U}(x^{\prime}-x^{\prime\prime})$ only, the Ohmic current
$e\chi^{\dagger}_{+\sigma_{{}_{B}}}(x^{\prime})v^{i}_{x^{\prime}x}\chi_{+\sigma_{{}_{B}}}(x)$
reads

$\displaystyle\chi^{\dagger}_{+\sigma_{{}_{B}}}(x^{+})v^{l}_{x^{+}x^{-}}\chi_{+\sigma_{{}_{B}}}(x^{-})=\hat{U}^{\dagger}(x^{+}-x^{\prime})\chi^{\dagger}_{+\sigma_{{}_{B}}}(x^{\prime})v^{l}_{x^{\prime}\,x^{-}}\chi_{+\sigma_{{}_{B}}}(x^{\prime\prime})\hat{U}(x^{\prime\prime}-x^{-})$
$\displaystyle=\left[1-{\imath}\int
V^{SM}(x^{+}-x^{\prime})dt^{\prime}d\vec{r^{\prime}}+\ldots\right]\chi^{\dagger}_{+\sigma_{{}_{B}}}(x^{\prime})v^{l}_{x^{\prime}x^{-}}\chi_{+\sigma_{{}_{B}}}(x^{-})=\left\\{1\right.-{(-\imath)}\int\int\chi^{\dagger}_{+\sigma_{{}_{B}}}(x^{+}-x^{\prime})\times$
$\displaystyle\sum_{i}\left[{e\over
c}{v^{\dagger}}^{i}_{x^{+},\,\bar{x}}A_{i}\right.+\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c){d\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}}(0)\ {{v^{\dagger}}^{i}_{x^{+},\,\bar{x}}}-{e\over
c}{d\widetilde{\Sigma_{AB}\Sigma_{BA}}\over dp_{i}^{AB}}(0)A_{i}$
$\displaystyle-\left.{e\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c)\over
2c}\sum_{j}{d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}dp_{j}^{AB}}(0)\
\left({v^{\dagger}}^{i}_{x^{+},\,\bar{x}}\,A_{j}+{v^{\dagger}}^{j}_{x^{+},\,\bar{x}}\,A_{i}\right)\right]$
$\displaystyle\left.\times\chi_{+\sigma_{{}_{B}}}(\bar{x})d\bar{t}\,d\vec{\bar{x}}\,dt^{\prime}\,d\vec{x}\,^{\prime}+\ldots\right\\}\chi^{\dagger}_{+\sigma_{{}_{B}}}(x^{\prime})\chi_{+\sigma_{{}_{B}}}(x^{\prime\prime})\hat{U}(x^{\prime\prime}-x^{\prime})v^{l}_{x^{\prime}x^{-}}.\
$ (S.15)

Here $x=\\{\vec{x},\ t\\}$ is the 4-vector of space coordinates $\vec{x}$ and
time $t$. The current (S.15) is the two-point function and, therefore, has to
be expressed through two-point
$G_{1}(\vec{r}-\vec{r}\,^{\prime},\,\vec{\bar{r}}-\vec{\bar{r}}\,^{\prime},t-t^{\prime})$
and four-point $G_{2}(x,\,x^{\prime},\,\bar{x},\,\bar{x}^{\prime})$ Green
functions. Using the relation $U(\vec{r}-\vec{r}\,^{\prime},t)=\imath
G_{1}(\vec{r}-\vec{r}\,^{\prime},t)$ between the evolution operator
$U(\vec{r}-\vec{r}\,^{\prime},t)$ and the one-particle Green function $\imath
G_{1}^{\mu\nu}(x,\
x^{\prime})=(\chi^{\dagger})^{\mu}(x)\chi^{\nu}(x^{\prime})$ MyLambert , one
can express the product of the field components $\chi^{\mu},\ \chi^{\nu}$,
entering in the form of multipliers in (S.15), through
$G_{1}(\vec{r}-\vec{r}\,^{\prime},\,\vec{\bar{r}}-\vec{\bar{r}}\,^{\prime},t-t^{\prime})$
as

$G_{1}(\vec{r}-\vec{r}\,^{\prime},\,\vec{\bar{r}}-\vec{\bar{r}}\,^{\prime},t-t^{\prime})={1\over\imath^{2}}U(\vec{r}-\vec{r}\,^{\prime},t-t^{\prime})U(\vec{r}\,^{\prime}-\vec{\bar{r}},t)={1\over\imath^{2}}U(\vec{r}-\vec{\bar{r}},t)$

and through $G_{2}(x,\,x^{\prime},\,\bar{x},\,\bar{x}^{\prime})$ as
$G_{2}(x,\,x^{\prime},\,\bar{x},\,\bar{x}^{\prime})=(\imath)^{2}\chi^{\dagger}(x)\chi^{\dagger}(x^{\prime})\chi(\bar{x})\chi(\bar{x}^{\prime}).$
Now, we take into account the linear in $A_{i}$ terms only in Eq. (S.15), then
after the standard procedure of transformation of the operator product into a
normal form Dyson , one can rewrite $j_{i}^{Ohm}$ in (IV.6) through two-point
$G_{1}(\vec{r}-\vec{r}\,^{\prime},\,\vec{\bar{r}}-\vec{\bar{r}}\,^{\prime},t-t^{\prime})$
and four-point $G_{2}(x,\,x^{\prime},\,\bar{x},\,\bar{x}^{\prime})$ Green
functions:

$\displaystyle
j_{l}^{Ohm}=\sum_{\mu}e{\chi^{\mu}}^{\dagger}_{+\sigma_{{}_{B}}}(x^{+})v^{l}_{x^{+}x^{-}}\chi^{\mu}_{+\sigma_{{}_{B}}}(x^{-})$
$\displaystyle=\mbox{Tr}\left\\{\chi^{\dagger}_{+\sigma_{{}_{B}}}(x^{\prime})\chi_{+\sigma_{{}_{B}}}(x^{\prime\prime})+\imath{-\imath
e\over(\imath)^{3}}\int\int\sum_{i}\left[{e\over
c}{v^{\dagger}}^{i}_{x^{+},\,\bar{x}}A_{i}+\widetilde{\Sigma_{AB}\Sigma_{BA}}\left(\vec{p}_{AB}-{e\vec{A}\over
c}\right){d\widetilde{\Sigma_{AB}\Sigma_{BA}}\over dp_{i}^{AB}}(0)\
{{v^{\dagger}}^{i}_{x^{+},\,\bar{x}}}\right.\right.$
$\displaystyle-\left.{e\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c)\over
2c}\sum_{j}{d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}dp_{j}^{AB}}(0)\
\left({v^{\dagger}}^{i}_{x^{+},\,\bar{x}}\,A_{j}+{v^{\dagger}}^{j}_{x^{+},\,\bar{x}}\,A_{i}\right)\right]$
$\displaystyle\times\left.G_{2}(x^{+},\,\bar{x},\,x^{\prime},\,x^{\prime\prime})G_{1}(x^{\prime\prime},\
x^{\prime})d\bar{t}\,d^{2}\vec{\bar{x}}\,dt^{\prime}\,d^{2}\vec{x}\,^{\prime}dt^{\prime\prime}\,d^{2}\vec{x^{\prime\prime}}+\ldots\right\\}v^{l}_{x^{\prime}x^{-}}.\
\ \ $ (S.16)

Due to the fact that the Green functions $G_{1}$, $G_{2}$ are symmetric at
change $\vec{v}\to-\vec{v}$, the integrals with expressions proportional to
$v_{i}$, vanish in (S.16). Interparticle Coulomb interaction renormalizes the
charge carriers mass DiracFieldTheory  in a way as it is shown in fig. 4.
Therefore, ${d\widetilde{\Sigma_{AB}\Sigma_{BA}}\over dp_{i}^{AB}}\approx 0$
near the Majorana zero-energy state, and respectively the third term in the
equation (S.16) is also equal to zero.

###  XIX.2 Low-frequency limit of the conductivity in the models $N=2$ and
$N=3$ with spatial dispersion

Let us analyze the behavior of the low-frequency conductivity of the models
$N=2$ and $N=3$ with spatial dispersion of charge carriers. First, let us
prove the following theorem

######  Theorem 1.

Low-frequency conductivity of the charge carriers system with spatial
dispersion for the $N=3$ model, in contrast to the model $N=2$, takes nonzero
finite values at frequencies $\omega\ll 1$.

Proof Self-energy operator $\hat{\Sigma}_{self}$ is expressed through the
polarization operator $\hat{\Pi}$ as Kraft-Ropke :

$\hat{\Sigma}_{self}(\omega,\ k)G_{1}=i\hbar V(k)\hat{\Pi}(\omega,\ k),$
(S.17)

where $V(k)$ is the Coulomb interaction operator in the momentum
representation, $G_{1}$ is the one-particle Green function. In the
neighborhood of the Dirac point, the one-particle Green function of negatively
(positively) charged carrier is $G_{1}={\hbar^{-1}\over(\omega\mp k)}$. Let us
use the Feynman diagram technique. The Feynman diagram for $\hat{\Pi}$ has the
form $\hat{\Pi}$ = , The Feynman representation for $i\hbar V(k)$ is a wavy
line. For holes, the left hand side of the expression (S.17) is represented by
the diagram $\hat{\Sigma}_{self}(\omega,\ k)G_{1}=$ , whereas the right hand
side of (S.17) is given by the diagram $i\hbar V(k)\hat{\Pi}(\omega,\ k)=$ .
Using the diagram technique, the expansion of the expression (S.17) into
perturbation theory series is written as

$\includegraphics[width=34.14322pt,height=14.22636pt,angle={0}]{hole-
selfenergy-s.pdf}=\includegraphics[width=31.2982pt,height=28.45274pt,angle={0}]{selfenergy-
through-polariz-
operator1-s.pdf}+\includegraphics[width=31.2982pt,height=28.45274pt,angle={0}]{selfenergy-
through-polariz-
operator2-s.pdf}+\includegraphics[width=31.2982pt,height=28.45274pt,angle={0}]{selfenergy-
through-polariz-operator3-s.pdf}.$ (S.18)

Here, by the cross ”$\times$” there have been marked crossed out lines. The
first and the second terms in the expansion (S.18) are called the Hartree-Fock
approximation for $\Sigma^{XC}$, the rest terms are for correlation
interaction (correlation energy),
$\hat{\Sigma}_{self}=\Sigma^{xc}+\Sigma^{c}.$ Let us account for terms up to
the second order in series on $V$ in (S.18) as

$\includegraphics[width=34.14322pt,height=14.22636pt,angle={0}]{hole-
selfenergy-s.pdf}\approx\includegraphics[width=34.14322pt,height=11.38092pt,angle={0}]{hole-
excheng-s.pdf}+\includegraphics[width=19.91684pt,height=31.2982pt,angle={0}]{hole-
polariz1st-
order-s.pdf}+\includegraphics[width=39.83368pt,height=36.98866pt,angle={0}]{hole-
polariz2st-order-s.pdf}.$ (S.19)

From the form of the first and the second diagrams in expansion of (S.19) it
follows that the Hartree-Fock approximation $\Sigma^{xc}=\Sigma^{x}+V^{self}$
consists of the exchange $\Sigma^{x}$ and the self-consistent Coulomb
potential $V^{self}$ which is the polarization correction of the first-order
on interaction $V$.

From the another side, the self-energy $\hat{\Sigma}_{self}(\omega,\ k)$ is
the interaction energy for the current of charge carriers with an
electromagnetic field $\hat{A}_{\mu}(\omega,\ k)$:

$\left\langle 0\right|\hat{\Sigma}_{self}(\omega,\
k)\left|0\right\rangle=\left\langle
0\right|\hat{\chi}^{\dagger}e\gamma^{\mu}\hat{\chi}\hat{A}_{\mu}(\omega,\
k)\left|0\right\rangle\equiv\left\langle
0\right|\hat{j}^{\mu}\hat{A}_{\mu}(\omega,\ k)\left|0\right\rangle.$

Here $\gamma_{\mu}$ are the Dirac matrices, $\left|0\right\rangle$ is the
vacuum state. In the Hartree-Fock approximation (S.19) electron-hole decay is
absent: $\Im m\,\Sigma^{xc}(k)=0$, and respectively in this approximation we
can rewrite imaginary part of the expression (S.17) through the photon
creation and annihilation operators $\hat{A}^{\dagger}_{\mu}(\omega,\ k)$,
$\hat{A}_{\mu}(\omega,\ k)$ as

$\Im m\left.\langle 0\right|\hat{j}^{\mu}\hat{A}^{\dagger}_{\mu}(\omega,\
k)\left|0\rangle\right.=\Im m[(i\hbar)^{2}V^{2}(k)\hat{\Pi}^{RPA}(\omega,\
k)G_{1}]=\Im m[i\hbar V(k)\hat{\Pi}^{(1)}(\omega,\ k)G_{1}],$ (S.20)

where $\hat{\Pi}^{(1)}=i\hbar V(k)\hat{\Pi}^{RPA}$, $\hat{\Pi}^{RPA}$ is the
polarization operator in the random phases approximation (RPA). Expressing
current component $\hat{j}_{\mu}$ through the complex conductivity
$\sigma^{*}(\omega,\ k)$ as $\hat{j}_{\mu}=\sigma^{*}(\omega,\
k)\hat{A}_{\mu}(\omega,k)$ and using commutation relations, one can transform
(S.20) to the following form

$\Im m[\sigma^{*}(\omega,\ k)\left.\langle
0\right|\left[\hat{A}^{\dagger}_{\mu}(\omega,\ k),\hat{A}_{\mu}(\omega,\
k)\right]\left|0\rangle\right.]=\Im m[\imath\hbar V(k)\hat{\Pi}^{(1)}(\omega,\
k)G_{1}],$ (S.21)

where the vacuum expectation $\left.\langle
0\right|\left[\hat{A}^{\dagger}_{\mu}(\omega,\ k),\hat{A}_{\mu}(\omega,\
k)\right]\left|0\rangle\right.$ from the commutator
$\left[\hat{A}^{\dagger}_{\mu}(\omega,\ k),\hat{A}_{\mu}(\omega,\ k)\right]$
is a photon propagator ${1\over\omega^{2}-k^{2}}$. For the case of small
spatial dispersion $k\to 0$, one has $G_{1}=1/(\hbar\omega)$, and the
expression (S.21) becomes ”optical” alternating-current (ac) conductivity
(Kraft-Ropke ):

${\sigma^{*}(\omega)\over\omega}=\hbar(iV)^{2}(k)\hat{\Pi}^{RPA}(\omega,\
k)=i(1-\epsilon^{*}(\omega)),$ (S.22)

where $\epsilon^{*}=1-V(k)\hat{\Pi}^{(1)}$ is a complex dielectric
permeability.

Now, let us consider the case of spatial dispersion of states $\omega\leq
k>0,\ \omega>0$ at very small frequencies $\omega\ll 1$. In the model $N=2$ of
massless pseudo-fermions with spatial dispersion, the expression (S.21) for
the hole conductivity $\sigma^{*}_{min}$ having positive values is transformed
to the form

$\Im m\ {\sigma^{*}_{min}(\omega,\ k)}=(\omega-k)\Im m\
[iV(k)\hat{\Pi}^{(1)}(\omega,\ k)].$ (S.23)

Because of $\omega-k\leq 0$, in accord with formula (S.22) the left hand side
of (S.23) gains physical positive values if the real part of the polarization
operator $\hat{\Pi}^{(1)}$ in the right hand side equals zero. The last means
that the direct-current (dc) conductivity (minimal dc-conductivity)
$\sigma_{min}$ in the model $N=2$ is equal to zero.

In the model $N=3$, the action of the operators
$\left(\Sigma_{rel}^{x}\right)_{AB}$ (III.3) and
$\left(\Sigma_{rel}^{x}\right)_{BA}$ (III.6) provides the squeezing $k\to
k_{AB}$ of the part of states and respectively the expression for (S.23) is
transformed to the following:

$\Im m\ {\sigma^{*}_{min}(\omega,\ k_{AB})}=\Im m\
(iV)^{2}(k_{AB})\hat{\Pi}^{RPA}(\omega,\ k_{AB})(\omega-k_{AB}).$ (S.24)

If $\omega>v_{F}k_{AB}$, due to dynamic decrease of spatial dispersion there
could exist part of quasi-particle states (S.24). It is possible because
$\omega$ can not gain the values less than vacuum frequency $\omega_{0}$. The
last means that in the model $N=3$, $\sigma_{min}$ can get non-zero value as
well as optical conductivity. Thus, we prove that contrary to the $N=2$ model,
the model $N=3$ can demonstrate non-zero dc-conductivity.

Let us write down the Fourier-Laplace image of the operator $\vec{v}_{AB}$ in
the coordinate representation of secondary quantization

$\displaystyle\int e^{-\imath\omega
t}e^{-\imath\vec{p}\cdot\vec{r}}\chi^{\dagger}_{+\sigma_{{}_{B}}}(x)v^{i}_{x^{+}x^{-}}\chi_{+\sigma_{{}_{B}}}(x)dt\,d\vec{r}.$
(S.25)

Because of the fact that the operator $\vec{v}_{AB}$ does not depend on time,
trace ”Tr” of the Fourier components of (S.25) in representation, where the
operator $H_{0}$ is diagonal, can be rewritten in the matrix form as

$\displaystyle\mbox{Tr }\int e^{\imath\omega
t}e^{-\imath\vec{k}\cdot\vec{r}}\chi^{\dagger}_{+\sigma_{{}_{B}}}(x^{+})v^{i}_{x^{+}x^{-}}\chi_{+\sigma_{{}_{B}}}(x^{-})dt\,d\vec{r}={1\over(2\pi)^{6}}\int
e^{-i(\omega^{+}-\omega^{-}-\omega)t}e^{\imath(\vec{p^{+}}-\vec{p^{-}}-\vec{k})\cdot\vec{r}}$
$\displaystyle\times\mbox{Tr }\
\chi^{\dagger}_{+\sigma_{{}_{B}}}(p^{+},\omega^{+})v^{i}_{\vec{p}^{+}-\vec{p}^{-},\,(\vec{p}^{+}+\vec{p}^{-})/2}\,\chi_{+\sigma_{{}_{B}}}(p^{-},\omega^{-})dt\,d\vec{r}\,d\omega^{+}\,d\omega^{-}d^{2}\vec{p}\,^{+}d^{2}\vec{p}\,^{-}\
\ $ $\displaystyle=\mbox{Tr }{1\over(2\pi)^{3}}\int
d\omega^{+}d\omega^{-}d\vec{p}\,^{+}d\vec{p}\,^{-}\delta(\omega^{+}-\omega^{-}-\omega)\
\delta(\vec{p^{+}}-\vec{p^{-}}-\vec{k})\chi^{\dagger}_{+\sigma_{{}_{B}}}(p^{+},\omega^{+})v^{i}_{\vec{p}^{+}-\vec{p}^{-},p}\,\chi_{+\sigma_{{}_{B}}}(p^{-},\omega^{-})$
$\displaystyle=\mbox{Tr }\
{1\over(2\pi)^{3}}\int\chi^{\dagger}_{+\sigma_{{}_{B}}}(p^{+},\omega^{+-}+\omega^{\prime})v^{i}_{\vec{p}^{+}-\vec{p}^{-},p}\,\chi_{+\sigma_{{}_{B}}}(p^{-},\omega^{\prime})\,d^{2}\vec{p}\
d\omega^{\prime}$ (S.26)

at the following conditions:

$\displaystyle\vec{p\,}^{\pm}=\vec{p}\pm\vec{k}/2,\quad\omega^{\prime}=\omega^{+}-\omega,\quad\omega^{+-}=\omega^{+}-\omega^{-}=\omega$
(S.27)

owing to $\delta$-functions entering the expression (S.26).

Let us perform the Wick turn (change from time to imaginary time) in order to
find Fourier-Laplace image of the expression (S.25) at finite temperatures
$T$. In this representation, the Fourier image of one-particle Green function
$G_{1}(\vec{r}_{1},\ \vec{r}\,^{\prime}_{1},\ t_{1}-t^{\prime}_{1})$ has the
form Kraft-Ropke

$\displaystyle G_{1}(\vec{r}_{1},\ \vec{r}\,^{\prime}_{1},\
z_{\nu})=\int_{0}^{-i\hbar\beta}dt_{1}e^{iz_{\nu}(t_{1}-t^{\prime}_{1})}G_{1}(x_{1},\
x^{\prime}_{1}),$ (S.28)

where $\bar{\beta}=1/(cT)$, $z_{\nu}$ are Matsubara frequencies defined by the
following expressions

$\displaystyle z_{\nu}={\pi\nu\over-i\hbar\bar{\beta}}+{\mu\over\hbar},\ \ \
\nu=\left\\{\begin{array}[]{c}\pm 1,\ \pm 3,\ \ldots\mbox{ for fermions},\\\
0,\ \pm 2,\ \pm 4,\ \ldots\mbox{ for bosons}.\end{array}\right.$ (S.31)

Analytical continuation of the function $G_{1}(x_{1},\ x^{\prime}_{1})$ to
interval $(0,\ -i\hbar\bar{\beta})$ is an advanced correlation function
$G^{>}_{1}$ (see Kraft-Ropke ):

$\displaystyle G_{1}(\vec{r}_{1},\ \vec{r}\,^{\prime}_{1},\
z_{\nu})=\int_{-\infty}^{+\infty}\int_{0}^{-i\hbar\bar{\beta}}{d\omega\over
2\pi i\hbar}\
dt_{1}e^{i(z_{\nu}-\omega)(t_{1}-t^{\prime}_{1})}G^{>}_{1}(\vec{r}_{1},\
\vec{r}\,^{\prime}_{1},\ \omega)=-\int_{-\infty}^{+\infty}{d\omega\over
2\pi\hbar}G^{>}_{1}(\vec{r}_{1},\ \vec{r}\,^{\prime}_{1},\
\omega){e^{(z_{\nu}-\omega)\hbar\bar{\beta}}-1\over(z_{\nu}-\omega)}$
$\displaystyle=-\int_{-\infty}^{+\infty}{d\omega\over
2\pi\hbar}G^{>}_{1}(\vec{r}_{1},\ \vec{r}\,^{\prime}_{1},\
\omega){e^{(\mu-\omega\hbar)\bar{\beta}}e^{i\pi\nu}-1\over(z_{\nu}-\omega)}.\
\ \ $ (S.32)

We note that the right-hand side of the expression (S.32) includes the Fourier
image of the retarded correlation function $G^{>}_{1}(\vec{r}_{1},\
\vec{r}\,^{\prime}_{1},\ \omega)$:

$\displaystyle G_{1}(\vec{r}_{1},\ \vec{r}\,^{\prime}_{1},\
z_{\nu})=-\int_{-\infty}^{+\infty}{d\omega\over
2\pi\hbar}{G^{<}_{1}(\vec{r}_{1},\ \vec{r}\,^{\prime}_{1},\
\omega)e^{i\pi\nu}-G^{>}_{1}(\vec{r}_{1},\ \vec{r}\,^{\prime}_{1},\
\omega)\over(z_{\nu}-\omega)}$
$\displaystyle=\int_{-\infty}^{+\infty}{d\omega\over
2\pi\hbar}{G^{>}_{1}(\vec{r}_{1},\ \vec{r}\,^{\prime}_{1},\ \omega)\mp
G^{<}_{1}(\vec{r}_{1},\ \vec{r}\,^{\prime}_{1},\
\omega)\over(z_{\nu}-\omega)}=\int_{-\infty}^{+\infty}{d\omega\over
2\pi\hbar}{\tilde{A}(\vec{r}_{1},\ \vec{r}\,^{\prime}_{1},\
\omega)\over(z_{\nu}-\omega)}.$ (S.33)

In (S.33) the upper sign refers to bosons, the lower sign to fermions;
function $\tilde{A}(\vec{r}_{1},\ \vec{r}\,^{\prime}_{1},\ \omega)$ is called
the spectral weight

$\displaystyle\tilde{A}(p,\ \omega)=A^{HF}(p,\ \omega)+A^{C}(p,\
\omega)=G_{1}(\omega+\imath\epsilon)-G_{1}(\omega-\imath\epsilon)$
$\displaystyle={2\over\imath}{\Im m\,\Sigma^{C}(p,\
\omega)\over[\hbar\omega-E(p)-\Re e\,\Sigma^{C}(p,\ \omega)]^{2}+[\Im
m\,\Sigma^{C}(p,\ \omega)]^{2}},$ (S.34)

where $E(p),\ A^{HF}$ is the energy of quasi-particle excitation and the
spectral weight in the Hartree-Fock approximation, $A^{C}$ is the many-
particle correction to the spectral weight. If the life time $\tau=1/(2\,\Im
m\,\Sigma^{C})=1/\Gamma$ of the quasi-particle excitation is large, then the
smallness of the decay rate $\Gamma$ leads to the following relation

$\displaystyle\lim_{\Gamma\to 0}\tilde{A}(p,\ \omega)=A^{HF}(p,\
\omega)=\delta(\hbar\omega-E(p))/\imath.$ (S.35)

The temperature Green’s functions of 2D fermions in the diagonal matrix
representation are obtained by summation of the Fourier-Laplace images on
Matsubara frequencies $z_{1}$, $z_{12}=z_{2}+z_{1}$ with subsequent
integration over momenta $\vec{p}$ Kraft-Ropke :

$\displaystyle
G_{1}(x-x^{\prime\prime})={1\over-\imath\bar{\beta}\hbar}\sum_{\nu_{1}}\int{d^{2}\vec{p}\over(2\pi)^{2}}e^{-\imath
z_{1}(t-t^{\prime\prime})}e^{-\imath\vec{p}\cdot(\vec{r}-\vec{r^{\prime\prime}})}G_{1}(z_{1},p),$
(S.36) $\displaystyle
G_{2}(\vec{r}_{1},\vec{r}_{2},\vec{r^{\prime}_{1}},\vec{r^{\prime}_{2}},t-t^{\prime})={1\over-\imath\bar{\beta}\hbar}\sum_{\nu_{12}}\int{d^{2}\vec{p}_{1}d^{2}\vec{p}_{2}d^{2}\vec{p^{\prime}_{1}}d^{2}\vec{p^{\prime}_{2}}\over(2\pi)^{8}}e^{-\imath\left(\vec{p}_{1}\cdot\vec{r}_{1}+\vec{p}_{2}\cdot\vec{r}_{2}-\vec{p^{\prime}_{2}}\cdot\vec{r^{\prime}_{2}}-\vec{p^{\prime}_{1}}\cdot\vec{r^{\prime}_{1}}\right)}$
$\displaystyle\times e^{-\imath
z_{12}(t-t^{\prime})}G_{2}(z_{12},p_{1},p_{2},p^{\prime}_{1},p^{\prime}_{2}).$
(S.37)

Here the perturbed Green function
$G_{2}(z_{12},p_{1},p_{2},p^{\prime}_{1},p^{\prime}_{2})$ is represented in
the form of the perturbation theory series

$\displaystyle
G_{2}(z^{-+},p^{-},p^{+},p_{1},p_{2})={1\over-\imath\bar{\beta}\hbar}\sum_{\nu^{+}}\left\\{G_{1}(z^{-+}-z^{+},p^{-})G_{1}(z^{+},p^{+})\right.$
$\displaystyle\times\left[(2\pi)^{4}\left(\delta(\vec{p^{-}}-\vec{p_{1}})\delta(\vec{p^{+}}-\vec{p_{2}})-\delta(\vec{p^{-}}-\vec{p_{2}})\delta(\vec{p^{+}}-\vec{p_{1}})\right)\right.$
$\displaystyle\left.\left.+\imath\int
V(\vec{\bar{p_{2}}}-\vec{p}\,^{+})(2\pi)^{2}\delta\left(\vec{\bar{p_{1}}}+\vec{\bar{p_{2}}}-\vec{p}\,^{-}-\vec{p}\,^{+}\right)G_{2}(z^{-+},\bar{p_{1}},\bar{p_{2}},p_{1},p_{2}){d\vec{\bar{p_{1}}}d\vec{\bar{p_{2}}}\over(2\pi)^{4}}\right]\right\\}.$
(S.38)

Taking into account (S.26, S.27), the two-point Green function $G_{1}(x,\
x^{\prime\prime})$ entering in (S.16) in the Matsubara representation is
written as

$\displaystyle
G_{1}(x,x^{\prime\prime})={1\over-\imath\bar{\beta}\hbar}\sum_{\nu^{+-}}\int{d^{2}\vec{p}\,^{-}\over(2\pi)^{2}}e^{-\imath
z^{+-}(t-t^{\prime\prime})}e^{-\imath\vec{p}\,^{-}\cdot(\vec{r}-\vec{r^{\prime\prime}})}G_{1}(p^{+},p^{-},z^{+-}),$
(S.39)

where $z^{+-}=z^{+}-z^{-}$. Substitution of the Fourier series for
$v^{i}_{x^{+}x^{-}}$:

$\displaystyle
v^{i}_{x^{+}x^{-}}=\int{d^{2}\vec{p}\over(2\pi)^{2}}e^{-\imath\vec{p}\cdot(\vec{r}^{+}-\vec{r}^{-})}v^{i}_{\vec{p}^{+}-\vec{p}^{-},p},$
(S.40)

$G_{1}(x,x^{\prime\prime})$ (S.39), and
$G_{2}(\vec{r}_{1},\vec{r}_{2},\vec{r^{\prime}_{1}},\vec{r^{\prime}_{2}},t-t^{\prime})$
(S.37) into (S.16) and analytical continuation of the obtained expression to
the whole complex plane allow us to find the contribution for
$j_{i}^{Ohm}(\omega,\ k)$ which is linear on $A_{i}$:

$\displaystyle
j_{i}^{Ohm}=\sum_{\mu}e{\chi^{\mu}}^{\dagger}_{+\sigma_{{}_{B}}}(x^{+})v^{i}_{x^{+}x^{-}}\chi^{\mu}_{+\sigma_{{}_{B}}}(x^{-})={(-1)\imath
e^{2}A_{i}\over(\imath)^{3}c(2\pi)^{14}}$
$\displaystyle\times\mbox{Tr}\left\\{\imath\int
d\bar{t}\,d^{2}\vec{\bar{r}}\,dt^{\prime}\,d^{2}\vec{r}\,^{\prime}dt^{\prime\prime}\,d^{2}\vec{r^{\prime\prime}}\left(1-\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c){d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}dp_{i}^{AB}}(0)\right)\right.\int d^{2}\vec{p}\
e^{-\imath\vec{p}\cdot(\vec{r}\,^{+}-\vec{\bar{r}})}{v^{\dagger}}^{i}_{\vec{p}^{+}-\vec{\bar{p}},\vec{p}}$
$\displaystyle\times{1\over-\imath\bar{\beta}\hbar}\sum_{\nu_{12}}\int
d^{2}\vec{p^{+}}d^{2}\vec{\bar{p}}d^{2}\vec{p^{\prime}}d^{2}\vec{p^{\prime\prime}}e^{-\imath\left(\vec{p}^{+}\cdot\vec{r}^{+}+\vec{\bar{p}}\cdot\vec{\bar{r}}-\vec{p^{\prime\prime}}\cdot\vec{r^{\prime\prime}}-\vec{p^{\prime}}\cdot\vec{r^{\prime}}\right)}e^{-\imath(z_{12}+\omega_{12})(t^{+}-t^{-}-t^{\prime}+\bar{t})}G_{2}(z_{12},p^{+},\bar{p},p^{\prime},p^{\prime\prime})$
$\displaystyle\times\left.{1\over-\imath\bar{\beta}\hbar}\sum_{\nu^{+-}}\int
d^{2}\vec{p}\,^{-}e^{-\imath(z^{+-}+\omega^{+-})(t^{\prime\prime}-t^{\prime})}e^{-\imath\vec{p}\,^{-}\cdot(\vec{r^{\prime\prime}}-\vec{r^{\prime}})}G_{1}(\hat{p}^{+},\hat{p}^{-},z^{+-})\right\\}\int
d^{2}\vec{\tilde{p}}e^{-\imath\vec{\tilde{p}}\cdot(\vec{r}\,^{\prime}-\vec{r}\,^{-})}v^{i}_{\vec{p}\,^{\prime}-\vec{p}^{-},\vec{\tilde{p}}/2}\
,\ \ \ $ (S.41)

where

$\displaystyle\omega_{12}=2\pi\nu_{12},\ \omega^{+-}=2\pi\nu^{+-};$ (S.42)

$(p^{+}+\bar{p})/2=p$, $(p^{\prime}+p^{-})/2=\tilde{p}/2$. Integrals of the
form
$(2\pi)^{-2}\int\exp(-\imath\Delta\vec{q}\cdot\vec{r})d^{2}\vec{r}=\delta(\Delta\vec{q})$,
over spatial coordinates lead to the appearance in (S.41) the Dirac
$\delta$-functions so that

$\displaystyle
j_{i}^{Ohm}=\sum_{\mu}e{\chi^{\mu}}^{\dagger}_{+\sigma_{{}_{B}}}(x^{+})v^{i}_{x^{+}x^{-}}\chi^{\mu}_{+\sigma_{{}_{B}}}(x^{-})={(-1)\imath
e^{2}A_{i}\over(\imath)^{3}c(2\pi)^{10}}$
$\displaystyle\times\mbox{Tr}\left\\{\imath\int
d\bar{t}\,dt^{\prime}\,d^{2}\vec{r}\,^{\prime}dt^{\prime\prime}\,\left(1-\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c){d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}dp_{i}^{AB}}(0)\right)\right.\int d^{2}\vec{p}\
e^{-\imath\vec{p}\cdot\vec{r}^{+}}{v^{\dagger}}^{i}_{\vec{p}^{+}-\vec{\bar{p}},\vec{p}}$
$\displaystyle\times{1\over-\imath\bar{\beta}\hbar}\sum_{\nu_{12}}\int\delta(\vec{p}-\vec{\bar{p}})d^{2}\vec{p}\,^{+}d^{2}\vec{\bar{p}}d^{2}\vec{p^{\prime}}d^{2}\vec{p^{\prime\prime}}e^{-\imath\left(\vec{p^{+}}\cdot\vec{r}^{+}-\vec{p^{\prime}}\cdot\vec{r^{\prime}}\right)}e^{-\imath(z_{12}+\omega_{12})(t^{+}-t^{-}-t^{\prime}+\bar{t})}G_{2}(z_{12},p^{+},\bar{p},p^{\prime},p^{\prime\prime})$
$\displaystyle\times\left.{1\over-\imath\bar{\beta}\hbar}\sum_{\nu^{+-}}\int
d^{2}\vec{p}\,^{-}\delta(\vec{p^{\prime\prime}}-\vec{p^{-}})e^{-\imath(z^{+-}+\omega^{+-})(t^{\prime\prime}-t^{\prime})}e^{-\imath\vec{p}\,^{-}\cdot(-\vec{r^{\prime}})}G_{1}(\hat{p}^{+},\hat{p}^{-},z^{+-})\right\\}\int
d^{2}\vec{\tilde{p}}e^{-\imath\vec{\tilde{p}}\cdot(\vec{r}\,^{\prime}-\vec{r}\,^{-})}v^{i}_{\vec{p}\,^{\prime}-\vec{p}^{-},\vec{\tilde{p}}/2}$
$\displaystyle={(-1)\imath
e^{2}A_{i}\over(\imath)^{3}c(2\pi)^{10}}\mbox{Tr}\left\\{\imath\int
d\bar{t}\,dt^{\prime}\,d^{2}\vec{r}\,^{\prime}dt^{\prime\prime}\,\left(1-\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c){d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}dp_{i}^{AB}}(0)\right)\right.$ $\displaystyle\times\int
d^{2}\vec{p}\
e^{-\imath\vec{p}\cdot\vec{r}^{+}}{v^{\dagger}}^{i}_{\vec{p}^{+}-\vec{p},\vec{p}}{1\over-\imath\bar{\beta}\hbar}\sum_{\nu_{12}}\int
d^{2}\vec{p}\,^{+}d^{2}\vec{p^{\prime}}e^{-\imath\left(\vec{p}\,^{+}\cdot\vec{r}^{+}-\vec{p^{\prime}}\cdot\vec{r^{\prime}}\right)}e^{-\imath(z_{12}+\omega_{12})(t^{+}-t^{-}-t^{\prime}+\bar{t})}G_{2}(z_{12},p^{+},p,p^{\prime},p^{-})$
$\displaystyle\times\left.{1\over-\imath\bar{\beta}\hbar}\sum_{\nu^{+-}}\int
d^{2}\vec{p}\,^{-}e^{-\imath(z^{+-}+\omega^{+-})(t^{\prime\prime}-t^{\prime})}e^{-\imath\vec{p}\,^{-}\cdot(-\vec{r^{\prime}})}G_{1}(\hat{p}^{+},\hat{p}^{-},z^{+-})\right\\}\int
d^{2}\vec{\tilde{p}}e^{-\imath\vec{\tilde{p}}\cdot(\vec{r}\,^{\prime}-\vec{r}\,^{-})}v^{i}_{\vec{p^{\prime}}-\vec{p}^{-},\vec{\tilde{p}}/2}\
.\ \ $ (S.43)

Integration of (S.43) over the rest spatial variables $\vec{r^{\prime}}$ gives

$\displaystyle
j_{i}^{Ohm}=\sum_{\mu}e{\chi^{\mu}}^{\dagger}_{+\sigma_{{}_{B}}}(x^{+})v^{i}_{x^{+}x^{-}}\chi^{\mu}_{+\sigma_{{}_{B}}}(x^{-})={(-1)\imath
e^{2}A_{i}\over(\imath)^{3}c(2\pi)^{8}}\mbox{Tr}\left\\{\imath\int
d\bar{t}\,dt^{\prime}\,dt^{\prime\prime}\,\left(1-\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c)\right.\right.$
$\displaystyle\times\left.{d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}dp_{i}^{AB}}(0)\right)\int d^{2}\vec{p}\
e^{-\imath\vec{p}\cdot\vec{r}^{+}}{v^{\dagger}}^{i}_{\vec{p}^{+}-\vec{p},\vec{p}}{1\over-\imath\bar{\beta}\hbar}\sum_{\nu_{12}}\int
d^{2}\vec{p}\,^{+}d^{2}\vec{p^{\prime}}e^{-\imath\vec{p}\,^{+}\cdot\vec{r}^{+}}e^{-\imath(z_{12}+\omega_{12})(t^{+}-t^{-}-t^{\prime}+\bar{t})}$
$\displaystyle\times
G_{2}(z_{12},p^{+},p,p^{\prime},p^{-})\left.{1\over-\imath\bar{\beta}\hbar}\sum_{\nu^{+-}}\int
d^{2}\vec{p}\,^{-}e^{-\imath(z^{+-}+\omega^{+-})(t^{\prime\prime}-t^{\prime})}G_{1}(\hat{p}^{+},\hat{p}^{-},z^{+-})\right\\}\int
d^{2}\vec{\tilde{p}}\ \delta(\vec{\tilde{p}}-\vec{p}\,^{-}-\vec{p^{\prime}})$
$\displaystyle\times
e^{-\imath\vec{\tilde{p}}\cdot(-\vec{r}\,^{-})}v^{i}_{\vec{p^{\prime}}-\vec{p}^{-},\vec{\tilde{p}}/2}={(-1)\imath
e^{2}A_{i}\over(\imath)^{3}c(2\pi)^{8}}\mbox{Tr}\left\\{\imath\int
d\bar{t}\,dt^{\prime}\,dt^{\prime\prime}\,\left(1-\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c){d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}dp_{i}^{AB}}(0)\right)\right.$ $\displaystyle\times\int
d^{2}\vec{p}\
e^{-\imath\vec{p}\cdot\vec{r}^{+}}{v^{\dagger}}^{i}_{\vec{p}^{+}-\vec{p},\vec{p}}{1\over-\imath\bar{\beta}\hbar}\sum_{\nu_{12}}\int
d^{2}\vec{p}\,^{+}d^{2}\vec{p^{\prime}}e^{-\imath\vec{p}\,^{+}\cdot\vec{r}^{+}}e^{-\imath(z_{12}+\omega_{12})(t^{+}-t^{-}-t^{\prime}+\bar{t})}G_{2}(z_{12},p^{+},p,p^{\prime},p^{-})$
$\displaystyle\times\left.{1\over-\imath\bar{\beta}\hbar}\sum_{\nu^{+-}}\int
d^{2}\vec{p}\,^{-}e^{-\imath(z^{+-}+\omega^{+-})(t^{\prime\prime}-t^{\prime})}G_{1}(\hat{p}^{+},\hat{p}^{-},z^{+-})\right\\}e^{-\imath(\vec{p}\,^{-}+\vec{p^{\prime}})\cdot(-\vec{r}\,^{-})}v^{i}_{\vec{p^{\prime}}-\vec{p}^{-},(\vec{p}\,^{-}+\vec{p^{\prime}})/2}\
.\ \ $ (S.44)

Taking into account (S.42), the integral of the form

$(2\pi)^{-1}\int\exp(-\imath\Delta\omega t)\ dt=(2\pi)^{-2}\int\exp(-\imath
2\pi\Delta\nu t)\ d(2\pi t)=(2\pi)^{-1}\delta(\Delta\nu)$

over time variables leads again to the Dirac $\delta$-functions in (S.44)
$\delta(\Delta\omega)$ and we obtain:

$\displaystyle
j_{i}^{Ohm}=\sum_{\mu}e{\chi^{\mu}}^{\dagger}_{+\sigma_{{}_{B}}}(x^{+})v^{i}_{x^{+}x^{-}}\chi^{\mu}_{+\sigma_{{}_{B}}}(x^{-})={(-1)\imath
e^{2}A_{i}\over(\imath)^{3}c(2\pi)^{8}}\mbox{Tr}\left\\{\imath\int
d\bar{t}\,dt^{\prime}\,dt^{\prime\prime}\,\left(1-\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c){d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}dp_{i}^{AB}}(0)\right)\right.$ $\displaystyle\times\int
d^{2}\vec{p}\
e^{-\imath\vec{p}\cdot\vec{r}^{+}}{v^{\dagger}}^{i}_{\vec{p}^{+}-\vec{p},\vec{p}}{1\over-\imath\bar{\beta}\hbar}\sum_{\nu_{12}}\int
d^{2}\vec{p}\,^{+}d^{2}\vec{p^{\prime}}e^{-\imath\vec{p}\,^{+}\cdot\vec{r}^{+}}e^{-\imath(z_{12}+2\pi\nu_{1(2)})(t^{+}-t^{-}-t^{\prime}+\bar{t})}G_{2}(z_{12},p^{+},p,p^{\prime},p^{-})$
$\displaystyle\times\left.{1\over-\imath\bar{\beta}\hbar}\sum_{\nu^{+-}}\int
d^{2}\vec{p}\,^{-}e^{-\imath(z^{+-}+2\pi\nu^{+(-)})(t^{\prime\prime}-t^{\prime})}G_{1}(\hat{p}^{+},\hat{p}^{-},z^{+-})\right\\}e^{-\imath(\vec{p}\,^{-}+\vec{p^{\prime}})\cdot(-\vec{r}\,^{-})}v^{i}_{\vec{p^{\prime}}-\vec{p}^{-},(\vec{p}\,^{-}+\vec{p^{\prime}})/2}$
$\displaystyle={(-1)\imath
e^{2}A_{i}\over(\imath)^{3}c(2\pi)^{8}}\mbox{Tr}\left\\{\imath{1\over-\imath\hbar\bar{\beta}}\int_{0}^{-\imath\hbar\bar{\beta}}d\bar{t}\,\left(1-\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c){d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}dp_{i}^{AB}}(0)\right)\right.\int d^{2}\vec{p}\
e^{-\imath\vec{p}\cdot\vec{r}^{+}}{v^{\dagger}}^{i}_{\vec{p}^{+}-\vec{p},\vec{p}}$
$\displaystyle\times{1\over-\imath\bar{\beta}\hbar}\sum_{\nu_{12}}\int
d^{2}\vec{p}\,^{+}d^{2}\vec{p^{\prime}}e^{-\imath\vec{p}\,^{+}\cdot\vec{r}^{+}}\delta\left(\hbar\nu_{12}+{\mu_{12}\over
2\pi}+\hbar\nu^{+-}+{\mu^{+-}\over
2\pi}\right)e^{-\imath(z_{12}+2\pi\nu_{1(2)})(t^{+}-t^{-}+\bar{t})}G_{2}(z_{12},p^{+},p,p^{\prime},p^{-})$
$\displaystyle\times\left.{1\over-\imath\bar{\beta}\hbar}\sum_{\nu^{+-}}\int
d^{2}\vec{p}\,^{-}\delta\left(\hbar\nu^{+-}+{\mu^{+-}\over
2\pi}\right)G_{1}(\hat{p}^{+},\hat{p}^{-},z^{+-})\right\\}e^{-\imath(\vec{p}\,^{-}+\vec{p^{\prime}})\cdot(-\vec{r}\,^{-})}v^{i}_{\vec{p^{\prime}}-\vec{p}^{-},(\vec{p}\,^{-}+\vec{p^{\prime}})/2}\
.\ \ \ \ \ \ $ (S.45)

Performing the change $\delta\left(\hbar\nu_{12}+{\mu_{12}\over
2\pi}+\hbar\nu^{+-}+{\mu^{+-}\over
2\pi}\right)\to\delta\left(z_{12}+z^{+-}\right)$, one transforms Eq. (S.45) as
follows

$\displaystyle
j_{i}^{Ohm}=\sum_{\mu}e{\chi^{\mu}}^{\dagger}_{+\sigma_{{}_{B}}}(x^{+})v^{i}_{x^{+}x^{-}}\chi^{\mu}_{+\sigma_{{}_{B}}}(x^{-})$
$\displaystyle={(-1)\imath
e^{2}A_{i}\over(\imath)^{3}c(2\pi)^{8}}\mbox{Tr}\left\\{\imath{1\over-\imath\hbar\bar{\beta}}\int_{0}^{-\imath\hbar\bar{\beta}}d\bar{t}\,\left(1-\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c){d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}dp_{i}^{AB}}(0)\right)\right.$ $\displaystyle\times\int
d^{2}\vec{p}\
e^{-\imath\vec{p}\cdot\vec{r}^{+}}{v^{\dagger}}^{i}_{\vec{p}^{+}-\vec{p},\vec{p}}{1\over-\imath\bar{\beta}\hbar}\sum_{\nu_{12}}\int
d^{2}\vec{p}\,^{+}d^{2}\vec{p^{\prime}}e^{-\imath\vec{p}\,^{+}\cdot\vec{r}^{+}}e^{\imath(z^{+-}+\pi\nu^{+-})(t^{+}-t^{-}+\bar{t})}G_{2}(z^{-+},p^{+},p,p^{\prime},p^{-})$
$\displaystyle\times\left.{1\over-\imath\bar{\beta}\hbar}\sum_{\nu^{+-}}\int
d^{2}\vec{p}\,^{-}\delta\left(\hbar\nu^{+-}+{\mu^{+-}\over
2\pi}\right)G_{1}(\hat{p}^{+},\hat{p}^{-},z^{+-})\right\\}e^{-\imath(\vec{p}\,^{-}+\vec{p^{\prime}})\cdot(-\vec{r}\,^{-})}v^{i}_{\vec{p^{\prime}}-\vec{p}^{-},(\vec{p}\,^{-}+\vec{p^{\prime}})/2}\
,\ \ $ (S.46)

where $z^{-+}=-z^{+-},\ \nu^{-+}=-\nu^{+-}$. Based on the following change:

${1\over-\imath\bar{\beta}\hbar}\sum_{\nu_{12}}\to{1\over-\imath\bar{\beta}\hbar
2\pi}\int_{0}^{-\imath\bar{\beta}\hbar}d(2\pi\nu_{1(2)})={1\over-\imath\bar{\beta}\hbar
2\pi}\int_{0}^{-\imath\bar{\beta}\hbar}d\,\omega_{12}={1\over 2\pi},$

the subsequent integration on variable $\bar{t}$ in (S.46)

$\int_{0}^{-\imath\hbar\bar{\beta}}d\bar{t}\,\left.e^{\imath(z^{+-}+\omega^{+-})\bar{t}}\right|_{z^{+-}\to
0}={1\over\imath}\left.{\exp\\{[{2\pi\imath\,\nu^{+(-)}\over\hbar\beta}+\omega^{+-}]\hbar\bar{\beta}\\}-1\over(z^{+-}+\omega^{+-})}\right|_{z^{+-}\to
0}={\hbar\bar{\beta}\over\imath}$

leads to the expression for the current

$\displaystyle
j_{i}^{Ohm}=\sum_{\mu}e{\chi^{\mu}}^{\dagger}_{+\sigma_{{}_{B}}}(x^{+})v^{i}_{x^{+}x^{-}}\chi^{\mu}_{+\sigma_{{}_{B}}}(x^{-})={(-1)\imath^{2}e^{2}A_{i}\over(\imath)^{3}c(2\pi)^{9}}\mbox{Tr}\left\\{{1\over-\imath\bar{\beta}\hbar}\sum_{\nu^{+-}}e^{\imath(z^{+-})(t^{+}-t^{-})}\int
d^{2}\vec{p}\,^{+}e^{-\imath\vec{p}\,^{+}\cdot\vec{r}^{+}}\right.$
$\displaystyle\times\int
d^{2}\vec{p}\,^{-}e^{\imath\vec{p}\,^{-}\cdot\vec{r}\,^{-}}\int d^{2}\vec{p}\
d^{2}\vec{p^{\prime}}e^{-\imath\vec{p}\cdot\vec{r}^{+}}e^{\imath\vec{p^{\prime}}\cdot\vec{r}\,^{-}}\left(1-\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c){d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}dp_{i}^{AB}}(0)\right)$
$\displaystyle\times{v^{\dagger}}^{i}_{\vec{p}^{+}-\vec{p},\vec{p}}\,G_{2}(z^{-+},p^{+},p,p^{\prime},p^{-})\left.G_{1}(\hat{p}^{+},\hat{p}^{-},z^{+-})\right\\}v^{i}_{\vec{p^{\prime}}-\vec{p}^{-},(\vec{p}\,^{-}+\vec{p^{\prime}})/2}\delta\left(\hbar\omega^{+-}+\mu^{+-}\right)\
.\ \ $ (S.47)

Substituting (S.33, S.34, S.35) into (S.47), we find finally

$\displaystyle
j_{i}^{Ohm}=\sum_{\mu}e{\chi^{\mu}}^{\dagger}_{+\sigma_{{}_{B}}}(x^{+})v^{i}_{x^{+}x^{-}}\chi^{\mu}_{+\sigma_{{}_{B}}}(x^{-})={(-1)\imath^{2}e^{2}A_{i}\over(\imath)^{3}c(2\pi)^{10}\hbar}\mbox{Tr}\left\\{{1\over-\imath\bar{\beta}\hbar}\sum_{\nu^{+-}}e^{\imath(z^{+-})(t^{+}-t^{-})}\right.$
$\displaystyle\times\int
d^{2}\vec{p}\,^{+}e^{-\imath\vec{p}\,^{+}\cdot\vec{r}^{+}}\int
d^{2}\vec{p}\,^{-}e^{\imath\vec{p}\,^{-}\cdot\vec{r}\,^{-}}\int d^{2}\vec{p}\
d^{2}\vec{p^{\prime}}e^{-\imath\vec{p}\cdot\vec{r}^{+}}e^{\imath\vec{p^{\prime}}\cdot\vec{r}\,^{-}}\left(1-\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c){d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}dp_{i}^{AB}}(0)\right)$
$\displaystyle\times{v^{\dagger}}^{i}_{\vec{p}^{+}-\vec{p},\vec{p}}\,G_{2}(z^{-+},p^{+},p,p^{\prime},p^{-})\int{d\omega\over\imath(z^{+-}-\omega)}\left[\delta(\hbar\omega-E(\hat{p}^{+},\hat{p}^{-}))\right.$
$\displaystyle\left.\left.+\left({\Gamma(\hat{p}^{+},\hat{p}^{-},\
\omega)\over[\hbar\omega-E(\hat{p}^{+},\hat{p}^{-})-\Re
e\,\Sigma^{C}(\hat{p}^{+},\hat{p}^{-},\
\omega)]^{2}+[\Gamma(\hat{p}^{+},\hat{p}^{-},\
\omega)/2]^{2}}-\delta(\hbar\omega-E(\hat{p}^{+},\hat{p}^{-}))\right)\right]\right\\}v^{i}_{\vec{p^{\prime}}-\vec{p}^{-},(\vec{p}\,^{-}+\vec{p^{\prime}})/2}$
$\displaystyle\times\delta\left(\hbar\omega^{+-}+\mu^{+-}\right)$
$\displaystyle={(-1)\imath^{2}e^{2}A_{i}\over(\imath)^{3}c(2\pi)^{10}\hbar}\mbox{Tr}\left\\{{1\over-\imath\bar{\beta}\hbar}\sum_{\nu^{+-}}e^{\imath(z^{+-})(t^{+}-t^{-})}\int
d^{2}\vec{p}\,^{+}e^{-\imath\vec{p}\,^{+}\cdot\vec{r}^{+}}\int
d^{2}\vec{p}\,^{-}e^{\imath\vec{p}\,^{-}\cdot\vec{r}\,^{-}}\right.\int
d^{2}\vec{p}\
d^{2}\vec{p^{\prime}}e^{-\imath\vec{p}\cdot\vec{r}^{+}}e^{\imath\vec{p^{\prime}}\cdot\vec{r}\,^{-}}$
$\displaystyle\times\left(1-\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c){d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}dp_{i}^{AB}}(0)\right){v^{\dagger}}^{i}_{\vec{p}^{+}-\vec{p},\vec{p}}\,G_{2}(z^{-+},p^{+},p,p^{\prime},p^{-})\left[{1\over\imath(z^{+-}-\omega(\hat{p}^{+},\hat{p}^{-}))}\right.$
$\displaystyle\left.\left.+2\pi\left({\Gamma(\hat{p}^{+},\hat{p}^{-},\
z^{+-})\over[\hbar z^{+-}-E(\hat{p}^{+},\hat{p}^{-})-\Re
e\,\Sigma^{C}(\hat{p}^{+},\hat{p}^{-},\
z^{+-})]^{2}+[\Gamma(\hat{p}^{+},\hat{p}^{-},\ z^{+-})/2]^{2}}-\delta(\hbar
z^{+-}-E(\hat{p}^{+},\hat{p}^{-}))\right)\right]\right\\}$
$\displaystyle\times
v^{i}_{\vec{p^{\prime}}-\vec{p}^{-},(\vec{p}\,^{-}+\vec{p^{\prime}})/2}\delta\left(\hbar\omega^{+-}+\mu^{+-}\right)\
,\ \ $ (S.48)

where $\omega(\hat{p}^{+},\hat{p}^{-})=E(\hat{p}^{+},\hat{p}^{-})/\hbar$.

The product $\delta(\hbar z^{+-}-E(\hat{p}^{+},\hat{p}^{-}))$ and
$\delta(\hbar\omega^{+-}+\mu^{+-})$ in (S.48) gives vanishing contribution

$\displaystyle
j_{i}^{Ohm}=\sum_{\mu}e{\chi^{\mu}}^{\dagger}_{+\sigma_{{}_{B}}}(x^{+})v^{i}_{x^{+}x^{-}}\chi^{\mu}_{+\sigma_{{}_{B}}}(x^{-})$
$\displaystyle={(-1)\imath^{2}e^{2}A_{i}\over(\imath)^{3}c(2\pi)^{10}\hbar}\mbox{Tr}\left\\{{1\over-\imath\bar{\beta}\hbar}\sum_{\nu^{+-}}e^{\imath(z^{+-})(t^{+}-t^{-})}\int
d^{2}\vec{p}\,^{+}e^{-\imath\vec{p}\,^{+}\cdot\vec{r}^{+}}\int
d^{2}\vec{p}\,^{-}e^{\imath\vec{p}\,^{-}\cdot\vec{r}\,^{-}}\int d^{2}\vec{p}\
d^{2}\vec{p^{\prime}}e^{-\imath\vec{p}\cdot\vec{r}^{+}}e^{\imath\vec{p^{\prime}}\cdot\vec{r}\,^{-}}\right.$
$\displaystyle\times\left(1-\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c){d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}dp_{i}^{AB}}(0)\right){v^{\dagger}}^{i}_{\vec{p}^{+}-\vec{p},\vec{p}}\,G_{2}(z^{-+},p^{+},p,p^{\prime},p^{-})\left[{1\over\imath(z^{+-}-\omega(\hat{p}^{+},\hat{p}^{-}))}\right.$
$\displaystyle\left.\left.+{2\pi\Gamma(\hat{p}^{+},\hat{p}^{-},\
z^{+-})\over[\hbar z^{+-}-E(\hat{p}^{+},\hat{p}^{-})-\Re
e\,\Sigma^{c}(\hat{p}^{+},\hat{p}^{-},\
z^{+-})]^{2}+[\Gamma(\hat{p}^{+},\hat{p}^{-},\
z^{+-})/2]^{2}}\right]\right\\}v^{i}_{\vec{p^{\prime}}-\vec{p}^{-},(\vec{p}\,^{-}+\vec{p^{\prime}})/2}\delta\left(\hbar\omega^{+-}+\mu^{+-}\right)\
.\ \ $ (S.49)

The presence of $\delta\left(\hbar\omega^{+-}+\mu^{+-}\right)$ in (S.49)
allows to state on possibility of formation of long-time order in considered
macroscopic phenomena Callen-et-al-1967 .

The Hamiltonian $H$ satisfies the equations of motion

$\displaystyle\imath\hbar{\partial\over\partial
t^{\pm}}\chi_{p}(\vec{r},\,t^{\pm})=H(p^{\pm})\,e^{\imath\omega^{\pm}t^{\pm}}\chi_{p}(\vec{r}),\
t^{+}>0,\ t^{-}<0$ (S.50)

for electrons and holes respectively. Therefore the frequencies
$\hat{\omega}^{+},\ \hat{\omega}^{-}$ satisfy the expressions
$\hat{\omega}^{+}\equiv\omega^{e}=H(p^{+})/\hbar$,
$\hat{\omega}^{-}\equiv\omega^{h}=H^{\dagger}(-p^{-})/\hbar$, that allows to
rewrite (S.49) as

$\displaystyle
j_{i}^{Ohm}=\sum_{\mu}e{\chi^{\mu}}^{\dagger}_{+\sigma_{{}_{B}}}(x^{+})v^{i}_{x^{+}x^{-}}\chi^{\mu}_{+\sigma_{{}_{B}}}(x^{-})$
$\displaystyle={(-1)\imath^{2}e^{2}A_{i}\over(\imath)^{3}c(2\pi)^{10}\hbar}\mbox{Tr}\left\\{{1\over-\imath\bar{\beta}\hbar}\sum_{\nu^{+-}}e^{\imath(z^{+-})(t^{+}-t^{-})}\int
d^{2}\vec{p}\,^{+}e^{-\imath\vec{p}\,^{+}\cdot\vec{r}^{+}}\int
d^{2}\vec{p}\,^{-}e^{\imath\vec{p}\,^{-}\cdot\vec{r}\,^{-}}\int d^{2}\vec{p}\
d^{2}\vec{p^{\prime}}e^{-\imath\vec{p}\cdot\vec{r}^{+}}e^{\imath\vec{p^{\prime}}\cdot\vec{r}\,^{-}}\right.$
$\displaystyle\times\left(1-\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c){d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}dp_{i}^{AB}}(0)\right){v^{\dagger}}^{i}_{\vec{p}^{+}-\vec{p},\vec{p}}\,G_{2}(z^{-+},p^{+},p,p^{\prime},p^{-})\left[{1\over\imath(z^{+-}-\omega^{e}(\hat{p}^{+})+\omega^{h}(-\hat{p}^{-}))}\right.$
$\displaystyle\left.\left.+{2\pi\Gamma(\hat{p}^{+},\hat{p}^{-},\
z^{+-})\over[\hbar z^{+-}-E(\hat{p}^{+},\hat{p}^{-})-\Re
e\,\Sigma^{c}(\hat{p}^{+},\hat{p}^{-},\
z^{+-})]^{2}+[\Gamma(\hat{p}^{+},\hat{p}^{-},\
z^{+-})/2]^{2}}\right]\right\\}v^{i}_{\vec{p^{\prime}}-\vec{p}^{-},(\vec{p}\,^{-}+\vec{p^{\prime}})/2}\delta\left(\hbar\omega^{+-}+\mu^{+-}\right)\
.\ \ $ (S.51)

For strongly correlated systems, which include Dirac materials, the energy
$E(\hat{p}^{+},\ \hat{p}^{-})$ of the pair of a free electron and a hole does
not exceed the correlational energy $\Re e\Sigma^{c}$ MyQuantumMatter2015 :
$E(\hat{p}^{+},\ \hat{p}^{-})<\Re e\Sigma^{c}$. The last allows to rewrite the
current of electron-hole pairs (S.51) in the form

$\displaystyle
j_{i}^{Ohm}=\sum_{\mu}e{\chi^{\mu}}^{\dagger}_{+\sigma_{{}_{B}}}(x^{+})v^{i}_{x^{+}x^{-}}\chi^{\mu}_{+\sigma_{{}_{B}}}(x^{-})={(-1)\imath^{2}e^{2}A_{i}\over(\imath)^{3}c(2\pi)^{10}\hbar}\mbox{Tr}\left\\{{1\over-\imath\bar{\beta}\hbar}\sum_{\nu^{+-}}e^{\imath(z^{+-})(t^{+}-t^{-})}\int
d^{2}\vec{p}\,^{+}e^{-\imath\vec{p}\,^{+}\cdot\vec{r}^{+}}\right.$
$\displaystyle\times\int
d^{2}\vec{p}\,^{-}e^{\imath\vec{p}\,^{-}\cdot\vec{r}\,^{-}}\int d^{2}\vec{p}\
d^{2}\vec{p^{\prime}}e^{-\imath\vec{p}\cdot\vec{r}^{+}}e^{\imath\vec{p^{\prime}}\cdot\vec{r}\,^{-}}\left(1-\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c){d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}dp_{i}^{AB}}(0)\right)$
$\displaystyle\times{v^{\dagger}}^{i}_{\vec{p}^{+}-\vec{p},\vec{p}}\,G_{2}(z^{-+},p^{+},p,p^{\prime},p^{-})\left[{\delta\left(\hbar\omega^{-+}+\mu^{-+}\right)\over-\imath(z^{-+}+\omega^{e}(\hat{p}^{+})-\omega^{h}(-\hat{p}^{-}))}\right.$
$\displaystyle\left.\left.+{2\pi\Gamma(\hat{p}^{+},\hat{p}^{-},\
z^{+-})\delta\left(\hbar\omega^{+-}+\mu^{+-}\right)\over[\hbar z^{+-}-\Re
e\,\Sigma^{c}(\hat{p}^{+},\hat{p}^{-},\
z^{+-})]^{2}+\Gamma^{2}(\hat{p}^{+},\hat{p}^{-},\
z^{+-})/4}\right]\right\\}v^{i}_{\vec{p^{\prime}}-\vec{p}^{-},(\vec{p}\,^{-}+\vec{p^{\prime}})/2}\
.\ \ $ (S.52)

Since $\Sigma^{c}(\hat{p}^{+},\hat{p}^{-},\
z^{+-})=\int_{-\infty}^{\infty}{d\omega\over
2\pi}{\Gamma(\hat{p}^{+},\hat{p}^{-},\ \omega)\over z^{+-}-\omega}$, the
correction for the correlation interaction describes the decay $\Gamma$ and
therefore influences only on initial conditions of macroscopic current
appearance. The charged-exciton model $N=3$ is protected by the hexagonal
symmetry, therefore in subsequent we calculate the Fourier component of the
current without accounting of the exciton decay as

$\displaystyle
j_{i}^{Ohm}(\vec{r}^{+},\vec{r}^{-},t^{+}-t^{-})=\sum_{\mu}e{\chi^{\mu}}^{\dagger}_{+\sigma_{{}_{B}}}(x^{+})v^{i}_{x^{+}x^{-}}\chi^{\mu}_{+\sigma_{{}_{B}}}(x^{-})$
$\displaystyle={1\over-\imath\bar{\beta}\hbar}\sum_{\nu^{+-}}e^{-\imath
z^{-+}(t^{+}-t^{-})}\int{d^{2}\vec{p}\,^{+}\over(2\pi)^{2}}e^{-\imath\vec{p}\,^{+}\cdot\vec{r}^{+}}\int{d^{2}\vec{p}\,^{-}\over(2\pi)^{2}}e^{\imath\vec{p}\,^{-}\cdot\vec{r}\,^{-}}j_{i}^{Ohm}(z^{-+},\
\vec{p}\,^{+},\vec{p}\,^{-}),\ \ $ (S.53) $\displaystyle j_{i}^{Ohm}(z^{-+},\
\vec{p}\,^{+},\vec{p}\,^{-})$
$\displaystyle={(-1)\imath^{2}e^{2}A_{i}\over(\imath)^{3}c(2\pi)^{6}}\mbox{Tr}\left\\{\int
d^{2}\vec{p}\
d^{2}\vec{p^{\prime}}e^{-\imath\vec{p}\cdot\vec{r}^{+}}e^{\imath\vec{p^{\prime}}\cdot\vec{r}\,^{-}}\left(1-\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c){d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}dp_{i}^{AB}}(0)\right)\right.$
$\displaystyle\left.\times{v^{\dagger}}^{i}_{\vec{p}^{+}-\vec{p},\vec{p}}\,G_{2}(z^{-+},p^{+},p,p^{\prime},p^{-})\left[{\delta\left(\hbar\omega^{-+}+\mu^{-+}\right)\over-\imath\hbar(z^{-+}+\omega^{e}(\hat{p}^{+})-\omega^{h}(-\hat{p}^{-}))}\right]\right\\}v^{i}_{\vec{p^{\prime}}-\vec{p}^{-},(\vec{p}\,^{-}+\vec{p^{\prime}})/2}\
.\ \ $ (S.54)

###  XIX.3 Hartree–Fock approximation

In the Hartree–Fock approximation, the substitution of Eq. (S.38) into (S.54)
gives

$\displaystyle j_{i}^{Ohm}(z^{-+},\
\vec{p}\,^{+},\vec{p}\,^{-})={\imath^{2}e^{2}A_{i}\over
c(2\pi)^{6}}\mbox{Tr}\left\\{\int d^{2}\vec{p}\
d^{2}\vec{p^{\prime}}e^{-\imath\vec{p}\cdot\vec{r}^{+}}e^{\imath\vec{p^{\prime}}\cdot\vec{r}\,^{-}}\left(1-\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c){d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}dp_{i}^{AB}}(0)\right)\right.$
$\displaystyle\times{v^{\dagger}}^{i}_{\vec{p}^{+}-\vec{p},\vec{p}}\,{1\over-\imath\bar{\beta}\hbar}\sum_{\nu^{+}}G_{1}(z^{-+}-z^{+},p^{-})G_{1}(z^{+},p^{+})(2\pi)^{4}\left(\delta(\vec{p^{\prime}}-\vec{p^{+}})\delta(\vec{p}-\vec{p^{-}})-\delta(\vec{p^{\prime}}-\vec{p^{-}})\delta(\vec{p}-\vec{p^{+}})\right)$
$\displaystyle\times\left.{\delta\left(\hbar\omega^{-+}+\mu^{-+}\right)\over\hbar(z^{-+}+\omega^{e}(\hat{p}^{+})-\omega^{h}(-\hat{p}^{-}))}\right\\}v^{i}_{\vec{p^{\prime}}-\vec{p}^{-},(\vec{p}\,^{-}+\vec{p^{\prime}})/2}\
.\ \ $ (S.55)

The matrix elements of the velocity operator have the following form:

$\begin{split}{v^{\dagger}}^{i}_{\vec{p}^{+}-\vec{p},\vec{p}}=e^{-\imath(\vec{p}^{+}-\vec{p})\cdot\vec{r}}{v^{\dagger}}^{i}(p),\quad
v^{i}_{\vec{p}^{+}-\vec{p}^{-},\vec{p}}=e^{\imath\vec{k}\cdot\vec{r}}v^{i}(p)\end{split}$
(S.56)

where $\vec{k}=\vec{p}^{+}-\vec{p}^{-}$, $\vec{r}=\vec{r}^{+}-\vec{r}^{-}$.
For the two-particle Green’s function represented by the following Feynman
diagram:

$\displaystyle\begin{array}[]{c}\vec{p^{\prime}}\longrightarrow\vec{p}^{+}\\\\[2.84526pt]
\ \vec{p}\longleftarrow\vec{p}^{-}\end{array}$ (S.59)

the substitution of (S.56) into (S.55) and integration over
$\delta(\vec{p}\,^{\prime}-\vec{p}^{+})$ first and then with
$\delta(\vec{p}-\vec{p}^{-})$ in the expression (S.53, S.55) gives

$\displaystyle
j_{i,1}^{Ohm}(\vec{r}^{+},\vec{r}^{-},t^{+}-t^{-})={\imath^{2}e^{2}A_{i}\over
c(2\pi)^{2}}\mbox{Tr}{1\over-\imath\bar{\beta}\hbar}\sum_{\nu^{+-}}e^{-\imath
z^{-+}(t^{+}-t^{-})}\int{d^{2}\vec{p}\,^{+}\over(2\pi)^{2}}e^{-\imath\vec{p}\,^{+}\cdot\vec{r}^{+}}\int{d^{2}\vec{p}\,^{-}\over
2\pi}e^{\imath\vec{p}\,^{-}\cdot\vec{r}\,^{-}}$
$\displaystyle\times\left\\{\int d^{2}\vec{p}\
\delta(\vec{p}-\vec{p^{-}})e^{-\imath\vec{p}\cdot\vec{r}^{+}}e^{\imath\vec{p}^{+}\cdot\vec{r}\,^{-}}\left(1-\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c){d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}dp_{i}^{AB}}(0)\right)\right.$
$\displaystyle\times{v^{\dagger}}^{i}_{\vec{p}^{+}-\vec{p},\vec{p}}\,{1\over-\imath\bar{\beta}\hbar}\sum_{\nu_{+}}G_{1}(z^{-+}-z^{+},p-k/2)G_{1}(z^{+},p^{+})\left.{\delta\left(\hbar\omega^{-+}+\mu^{-+}\right)\over\hbar(z^{-+}+\omega^{e}(\hat{p}^{+})-\omega^{h}(-\hat{p}^{-}))}\right\\}v^{i}_{\vec{k},\vec{p}}$
$\displaystyle={\imath^{2}e^{2}A_{i}\over
c(2\pi)^{2}}\mbox{Tr}{1\over-\imath\bar{\beta}\hbar}\sum_{\nu^{+-}}e^{-\imath
z^{-+}(t^{+}-t^{-})}\int{d^{2}\vec{p}\,^{+}\over(2\pi)^{2}}e^{-\imath\vec{p}\,^{+}\cdot\vec{r}^{+}}e^{\imath\vec{p}\cdot\vec{r}\,^{-}}$
$\displaystyle\times\left\\{\int d^{2}\vec{p}\
e^{-\imath\vec{p}\cdot\vec{r}^{+}}e^{\imath\vec{p}^{+}\cdot\vec{r}\,^{-}}e^{-\imath(\vec{p}^{+}-\vec{p})\cdot\vec{r}}e^{\imath\vec{k}\cdot\vec{r}}\left(1-\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c){d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}dp_{i}^{AB}}(0)\right)\right.$
$\displaystyle\times{v^{\dagger}}^{i}(p)\,{1\over-\imath\bar{\beta}\hbar}\sum_{\nu_{+}}G_{1}(z^{-+}-z^{+},p^{-})G_{1}(z^{+},p^{+})\times\left.{\delta\left(\hbar\omega^{-+}+\mu^{-+}\right)\over\hbar(z^{-+}+\omega^{e}(\hat{p}^{+})-\omega^{h}(-\hat{p}^{-}))}\right\\}v^{i}(p)={\imath^{2}e^{2}A_{i}\over
c(2\pi)^{2}}\mbox{Tr}{1\over-\imath\bar{\beta}\hbar}$
$\displaystyle\times\sum_{\nu^{+-}}e^{-\imath
z^{-+}(t^{+}-t^{-})}\int{d^{2}\vec{p}\,^{+}\over(2\pi)^{2}}d^{2}\vec{p}\
e^{-\imath\vec{p}\,^{+}\cdot(\vec{r}^{+}-\vec{r}^{-}+\vec{r})}e^{\imath\vec{p}\cdot(\vec{r}^{+}-\vec{r}^{-}-\vec{r})}e^{\imath\vec{k}\cdot\vec{r}}\left\\{\left(1-\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c){d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}dp_{i}^{AB}}(0)\right)\right.$
$\displaystyle\times{v^{\dagger}}^{i}(p)\,{1\over-\imath\bar{\beta}\hbar}\sum_{\nu_{+}}G_{1}(z^{-+}-z^{+},p^{-})G_{1}(z^{+},p^{+})\left.{\delta\left(\hbar\omega^{-+}+\mu^{-+}\right)\over\hbar(z^{-+}+\omega^{e}(\hat{p}^{+})-\omega^{h}(-\hat{p}^{-}))}\right\\}v^{i}(p).\
\ \ \ \ \ $ (S.60)

Similarly, we can find the contribution of the second Feynman diagram for two-
particle Green’s function:

$\displaystyle\begin{array}[]{c}\ \vec{p^{\prime}}\ \ \ \ \ \
\vec{p}^{+}\\\\[-2.84526pt] \searrow\\\\[-12.80373pt] \swarrow\\\\[-5.69054pt]
\ \ \vec{p}\ \ \ \ \ \ \vec{p}\,^{-}\end{array}$ (S.65)

into Ohmic current

$\displaystyle
j_{i,2}^{Ohm}(\vec{r}^{+},\vec{r}^{-},t^{+}-t^{-})=-{\imath^{2}e^{2}A_{i}\over
c(2\pi)^{2}}\mbox{Tr}{1\over-\imath\bar{\beta}\hbar}\sum_{\nu^{+-}}e^{-\imath
z^{-+}(t^{+}-t^{-})}\int{d^{2}\vec{p}\,^{+}\over(2\pi)^{2}}\delta(\vec{p}-\vec{p^{+}})e^{-\imath\vec{p}\,^{+}\cdot\vec{r}^{+}}\int{d^{2}\vec{p}\,^{-}\over
2\pi}e^{\imath\vec{p}\,^{-}\cdot\vec{r}\,^{-}}$
$\displaystyle\times\left\\{\int
d^{2}\vec{p}e^{-\imath\vec{p}\cdot\vec{r}^{+}}e^{\imath\vec{p}^{-}\cdot\vec{r}\,^{-}}\left(1-\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c){d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}dp_{i}^{AB}}(0)\right)\right.$
$\displaystyle\times{v^{\dagger}}^{i}_{\vec{p}+\vec{k}/2-\vec{p},\vec{p}}\,{1\over-\imath\bar{\beta}\hbar}\sum_{\nu_{+}}G_{1}(z^{-+}-z^{+},p^{-})G_{1}(z^{+},p+k/2)\left.{\delta\left(\hbar\omega^{-+}+\mu^{-+}\right)\over\hbar(z^{-+}+\omega^{e}(\hat{p}^{+})-\omega^{h}(-\hat{p}^{-}))}\right\\}v^{i}_{0,\vec{p}\,^{-}}=-{\imath^{2}e^{2}A_{i}\over
c(2\pi)^{2}}$
$\displaystyle\times\mbox{Tr}{1\over-\imath\bar{\beta}\hbar}\sum_{\nu^{+-}}e^{-\imath
z^{-+}(t^{+}-t^{-})}\int{d^{2}\vec{p}\,^{-}\over(2\pi)^{2}}e^{\imath
2\vec{p}\,^{-}\cdot\vec{r}\,^{-}}\left\\{\int d^{2}\vec{p}e^{-\imath
2\vec{p}\cdot\vec{r}^{+}}\left(1-\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c){d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}dp_{i}^{AB}}(0)\right)\right.$
$\displaystyle\times{v^{\dagger}}^{i}_{\vec{p}+\vec{k}/2-\vec{p},\vec{p}}\,{1\over-\imath\bar{\beta}\hbar}\sum_{\nu_{+}}G_{1}(z^{-+}-z^{+},p^{-})G_{1}(z^{+},p+k/2)\left.{\delta\left(\hbar\omega^{-+}+\mu^{-+}\right)\over\hbar(z^{-+}+\omega^{e}(\hat{p}^{+})-\omega^{h}(-\hat{p}^{-}))}\right\\}v^{i}_{0,\vec{p}\,^{-}}.\
\ $ (S.66)

###  XIX.4 Electrophysical current

For the case of spatial dispersion and transitions from the energy level
$E_{D}=0$ ($\vec{p}^{+}\to\vec{k}$, $\vec{p}^{-}=0$), the expression (S.66)
vanishes due to the fact that one of the matrix element of the velocity
operator has a constant value $v^{i}_{0,\vec{p}\,^{-}}\to v^{i}_{0,0}$ at
$\vec{p}^{+}\to\vec{k}$, $\vec{p}^{-}\to 0$. Therefore the expression (S.60)
gives the minimal dc-conductivity in the limit $p^{+}\to k$:

$\displaystyle\lim_{p^{+}\to
k}\sum_{\mu=1}^{2}j_{i,\mu}^{Ohm}(\vec{r}^{+},\vec{r}^{-},t^{+}-t^{-})\equiv
j_{i,\min}^{Ohm}(\vec{r}^{+},\vec{r}^{-},t^{+}-t^{-})$
$\displaystyle={\imath^{2}e^{2}A_{i}\over
c(2\pi)^{2}}\mbox{Tr}{1\over-\imath\bar{\beta}\hbar}\sum_{\nu^{+-}}e^{-\imath
z^{-+}(t^{+}-t^{-})}\int{d^{2}\vec{p}\,^{+}\over(2\pi)^{2}}d^{2}\vec{p}\
e^{-\imath(2\vec{p}\,^{+}-\vec{k})\cdot\vec{r}}\left\\{\left(1-\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c){d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}dp_{i}^{AB}}(0)\right)\right.$
$\displaystyle\times{v^{\dagger}}^{i}(p)\,{1\over-\imath\bar{\beta}\hbar}\sum_{\nu_{+}}G_{1}(z^{-+}-z^{+},p^{-})G_{1}(z^{+},p^{+})\left.{\delta\left(\hbar\omega^{-+}+\mu^{-+}\right)\over\hbar(z^{-+}+\omega^{e}(\hat{p}^{+})-\omega^{h}(-\hat{p}^{-}))}\right\\}v^{i}(p)$
$\displaystyle={\imath^{2}e^{2}A_{i}\over
c(2\pi)^{2}}\mbox{Tr}{1\over-\imath\bar{\beta}\hbar}\sum_{\nu^{+-}}e^{-\imath
z^{-+}(t^{+}-t^{-})}\int
e^{-\imath\vec{k}\cdot\vec{r}}{d^{2}\vec{k}\over(2\pi)^{2}}\int d^{2}\vec{p}\
\left\\{\left(1-\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c){d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}dp_{i}^{AB}}(0)\right)\right.$
$\displaystyle\times{v^{\dagger}}^{i}(p)\,{1\over-\imath\bar{\beta}\hbar}\sum_{\nu_{+}}G_{1}(z^{-+}-z^{+},p^{-})G_{1}(z^{+},p^{+})\left.{\delta\left(\hbar\omega^{-+}+\mu^{-+}\right)\over\hbar(z^{-+}+\omega^{e}(\hat{p}^{+})-\omega^{h}(-\hat{p}^{-}))}\right\\}v^{i}(p).\
\ $ (S.67)

###  XIX.5  Optical current

Transitions from low-lying energy levels $E\to E_{D}$ occur at excitation by
electromagnetic radiation of the optical range. Spatial dispersion is small at
optical frequencies, which are high relative to those used in electrophysics.
In the zeroth approximation, the dynamic (optical) conductivity is calculated
in the limit $p^{-}\to p$ at $k\to 0$. To find contribution of the first
Feynman diagram (S.59) into optical current $j_{i}^{opt}$, let us rewrite
(S.60) in the following form:

$\displaystyle\lim_{k\to
0}j_{i,1}^{Ohm}(\vec{r}^{+},\vec{r}^{-},t^{+}-t^{-})\equiv
j_{i,1}^{opt}(\vec{r}^{+},\vec{r}^{-},t^{+}-t^{-})$
$\displaystyle={\imath^{2}e^{2}A_{i}\over
c(2\pi)^{2}}\mbox{Tr}{1\over-\imath\bar{\beta}\hbar}\sum_{\nu^{+-}}e^{-\imath
z^{-+}(t^{+}-t^{-})}\int{[d^{2}\vec{k}+d^{2}\vec{p}\,^{-}]\over(2\pi)^{2}}d^{2}\vec{p}\
e^{-\imath(\vec{k}+2\vec{p}\,^{-})\cdot\vec{r}}\left\\{\left(1-\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c){d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}dp_{i}^{AB}}(0)\right)\right.$
$\displaystyle\times{v^{\dagger}}^{i}(p)\,{1\over-\imath\bar{\beta}\hbar}\sum_{\nu_{+}}G_{1}(z^{-+}-z^{+},p^{-})G_{1}(z^{+},p^{+})\left.{\delta\left(\hbar\omega^{-+}+\mu^{-+}\right)\over\hbar(z^{-+}+\omega^{e}(\hat{p}^{+})-\omega^{h}(-\hat{p}^{-}))}\right\\}v^{i}(p).\
\ \ \ \ \ $ (S.68)

Due to equality to zero of (S.66) at $p^{-}=0$, the expression (S.68) gives
the main contribution to $j_{i}^{opt}$ at $p\to 0$:

$\displaystyle\lim_{k,p\to
0}\sum_{\mu=1}^{2}j_{i,\mu}^{Ohm}(\vec{r}^{+},\vec{r}^{-},t^{+}-t^{-})\equiv\lim_{p,\
p^{-}\to 0}j_{i,1}^{opt}(\vec{r}^{+},\vec{r}^{-},t^{+}-t^{-})$
$\displaystyle={\imath^{2}e^{2}A_{i}\over
c(2\pi)^{2}}\mbox{Tr}{1\over-\imath\bar{\beta}\hbar}\sum_{\nu^{+-}}e^{-\imath
z^{-+}(t^{+}-t^{-})}\int{[d^{2}\vec{k}+d^{2}(\vec{p}-\vec{k}/2)]\over(2\pi)^{2}}d^{2}\vec{p}\
e^{-\imath\vec{k}\cdot\vec{r}}\left\\{\left(1-\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c){d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}dp_{i}^{AB}}(0)\right)\right.$
$\displaystyle\times{v^{\dagger}}^{i}(p)\,{1\over-\imath\bar{\beta}\hbar}\sum_{\nu_{+}}G_{1}(z^{-+}-z^{+},p^{-})G_{1}(z^{+},p^{+})\left.{\delta\left(\hbar\omega^{-+}+\mu^{-+}\right)\over\hbar(z^{-+}+\omega^{e}(\hat{p}^{+})-\omega^{h}(-\hat{p}^{-}))}\right\\}v^{i}(p)$
$\displaystyle={\imath^{2}e^{2}A_{i}\over
2c(2\pi)^{2}}\mbox{Tr}{1\over-\imath\bar{\beta}\hbar}\sum_{\nu^{+-}}e^{-\imath
z^{-+}(t^{+}-t^{-})}\int{d^{2}\vec{k}\over(2\pi)^{2}}d^{2}\vec{p}\
e^{-\imath\vec{k}\cdot\vec{r}}\left\\{\left(1-\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c){d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}dp_{i}^{AB}}(0)\right)\right.$
$\displaystyle\times{v^{\dagger}}^{i}(p)\,{1\over-\imath\bar{\beta}\hbar}\sum_{\nu_{+}}G_{1}(z^{-+}-z^{+},p^{-})G_{1}(z^{+},p^{+})\left.{\delta\left(\hbar\omega^{-+}+\mu^{-+}\right)\over\hbar(z^{-+}+\omega^{e}(\hat{p}^{+})-\omega^{h}(-\hat{p}^{-}))}\right\\}v^{i}(p).\
\ \ \ \ \ $ (S.69)

Contribution (S.66) of the second Feynman diagram (S.65) gives the following
correction to optical current at small values of $p^{-}\ll 1$, $k=0$:

$\displaystyle
j_{i,2}^{opt}(\vec{r}^{+},\vec{r}^{-},t^{+}-t^{-})=-{\imath^{2}e^{2}A_{i}\over
c(2\pi)^{2}}\mbox{Tr}{1\over-\imath\bar{\beta}\hbar}\sum_{\nu^{+-}}e^{-\imath
z^{-+}(t^{+}-t^{-})}\int{d^{2}\vec{p}\,^{-}\over(2\pi)^{2}}e^{-\imath\vec{p}\,^{-}\cdot(\vec{r}\,^{+}-\vec{r}\,^{-})}$
$\displaystyle\times\left\\{\int
d^{2}\vec{p}e^{-\imath\vec{p}\cdot(\vec{r}\,^{+}-\vec{r}\,^{-})}\left(1-\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c){d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}dp_{i}^{AB}}(0)\right)\right.{1\over-\imath\bar{\beta}\hbar}\sum_{\nu_{+}}G_{1}(z^{-+}-z^{+},p^{-})G_{1}(z^{+},p^{+})$
$\displaystyle\times\left.{\delta\left(\hbar\omega^{-+}+\mu^{-+}\right)\over\hbar(z^{-+}+\omega^{e}(\hat{p}^{+})-\omega^{h}(-\hat{p}^{-}))}\right\\}v^{i}_{0,\vec{p}\,^{-}}=-{\imath^{2}e^{2}A_{i}\over
c(2\pi)^{2}}\mbox{Tr}{1\over-\imath\bar{\beta}\hbar}\sum_{\nu^{+-}}e^{-\imath
z^{-+}(t^{+}-t^{-})}\int{d^{2}\vec{p}\over(2\pi)^{2}}e^{-\imath\vec{p}\cdot\vec{r}}$
$\displaystyle\times\left\\{\int
d^{2}\vec{\tilde{p}}e^{-\imath\vec{\tilde{p}}\cdot\vec{r}}\left(1-\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c){d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}dp_{i}^{AB}}(0)\right)\right.$
$\displaystyle\times{v^{\dagger}}^{i}({\tilde{p}})\,{1\over-\imath\bar{\beta}\hbar}\sum_{\nu_{+}}G_{1}(z^{-+}-z^{+},{\tilde{p}}^{-})G_{1}(z^{+},{\tilde{p}}^{+})\left.{\delta\left(\hbar\omega^{-+}+\mu^{-+}\right)\over\hbar(z^{-+}+\omega^{e}(\hat{\tilde{p}}^{+})-\omega^{h}(-\hat{\tilde{p}}^{-}))}\right\\}v^{i}({\tilde{p}}),\
\ $ (S.70)

where $\tilde{p}=p^{-}$. Taking into account that $\tilde{p}<p$, and
redesignating $p$ as $k$, we obtain the final expression for $j_{i,2}^{opt}$

$\displaystyle
j_{i,2}^{opt}(\vec{r}^{+},\vec{r}^{-},t^{+}-t^{-})=-{\imath^{2}e^{2}A_{i}\over
c(2\pi)^{2}}\mbox{Tr}{1\over-\imath\bar{\beta}\hbar}\sum_{\nu^{+-}}e^{-\imath
z^{-+}(t^{+}-t^{-})}\int{d^{2}\vec{k}\over(2\pi)^{2}}e^{-\imath\vec{k}\cdot\vec{r}}$
$\displaystyle\times\left\\{\int_{-\infty}^{-\vec{k}/2}d^{2}\vec{\tilde{p}}e^{-\imath\vec{\tilde{p}}\cdot\vec{r}}\left(1-\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c){d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}dp_{i}^{AB}}(0)\right)\right.$
$\displaystyle\times{v^{\dagger}}^{i}({\tilde{p}})\,{1\over-\imath\bar{\beta}\hbar}\sum_{\nu_{+}}G_{1}(z^{-+}-z^{+},{\tilde{p}}^{-})G_{1}(z^{+},{\tilde{p}}^{+})\left.{\delta\left(\hbar\omega^{-+}+\mu^{-+}\right)\over\hbar(z^{-+}+\omega^{e}(\hat{\tilde{p}}^{+})-\omega^{h}(-\hat{\tilde{p}}^{-}))}\right\\}v^{i}({\tilde{p}}).\
\ $ (S.71)

The correction (S.71) is small due to the exponential multiplier
$e^{-\imath\vec{\tilde{p}}\cdot\vec{r}}$. This correction takes into account
optical transitions from the below lying negative energy levels $E<E_{D}$, and
has oscillating character that even more diminishes the value of optical
conductivity. Oscillating character of optical conductivity of pure graphene
and its decrease at wavelengths less than 500 nm in visible range has been
observed in Nair2008Science .

The comparison of (S.67) and (S.69) demonstrates that optical current for
transition from zero-valued Dirac level $E_{D}=0$, in a factor two less than
electrophysical current.

###  XIX.6  Limit transition from optical to electrophysical current

In the known approaches Ando2002 ; Falkovsky ; Ziegler , the value for low-
frequency dynamical conductivity $\sigma_{dyn}(\omega)$ of pure graphene
demonstrates a jump to $\sigma_{dyn}(0)$ at zero frequency $\omega=0$ at
finite temperature and $T\to 0$. In the case of vanishing decay rate
$\Gamma=0$, $\sigma_{dyn}(\omega)$ becomes $e^{2}/(4\hbar)$ for all $T$ in the
limit $\omega\to 0$; conductivity becomes zero: $\sigma_{dyn}(0)=0$ at
$\omega=0$ and $T\to 0$ contrary to a finite value for minimal dc-conductivity
Ziegler . Let us show that in the model $N=3$ there is no conductivity jump in
the transition from the case $\omega(0)\to 0$ to the case
$\omega(0)\to\omega(k)$. To do this, let us consider the behaviour of the
correction (S.66) in the case of spatial dispersion $\omega(k)\to 0$, when
$\vec{p}-\vec{p}^{-}\to 0$:

$\displaystyle\lim_{\omega(k)\to
0}j_{i,2}^{Ohm}(\vec{r}^{+},\vec{r}^{-},t^{+}-t^{-})=-{\imath^{2}e^{2}A_{i}\over
c(2\pi)^{2}}\mbox{Tr}{1\over-\imath\bar{\beta}\hbar}\sum_{\nu^{+-}}e^{-\imath
z^{-+}(t^{+}-t^{-})}\int{d^{2}(\vec{p}-\vec{k}/2)\over(2\pi)^{2}}e^{-\imath(\vec{p}^{-}+\vec{p})\cdot\vec{r}}$
$\displaystyle\times\left\\{\int
d^{2}\vec{p}\left(1-\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c){d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}dp_{i}^{AB}}(0)\right)\right.$
$\displaystyle\times{v^{\dagger}}^{i}_{\vec{k}/2,\vec{p}}\,{1\over-\imath\bar{\beta}\hbar}\sum_{\nu_{+}}G_{1}(z^{-+}-z^{+},p^{-})G_{1}(z^{+},p+k/2)\left.{\delta\left(\hbar\omega^{-+}+\mu^{-+}\right)\over\hbar(z^{-+}+\omega^{e}(\hat{p}^{+})-\omega^{h}(-\hat{p}^{-}))}v^{i}_{0,\vec{p}\,^{-}}\right\\}=-{\imath^{2}e^{2}A_{i}\over
c(2\pi)^{2}}$
$\displaystyle\times\mbox{Tr}{1\over-\imath\bar{\beta}\hbar}\sum_{\nu^{+-}}e^{-\imath
z^{-+}(t^{+}-t^{-})}\int{d^{2}(\vec{p}-\vec{k}/2)\over(2\pi)^{2}}e^{-\imath(2\vec{p}^{-}+\vec{k})\cdot\vec{r}}\left\\{\int
d^{2}\vec{p}\left(1-\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c){d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}dp_{i}^{AB}}(0)\right)\right.$
$\displaystyle\times{v^{\dagger}}^{i}(p)\,{1\over-\imath\bar{\beta}\hbar}\sum_{\nu_{+}}G_{1}(z^{-+}-z^{+},p^{-})G_{1}(z^{+},p^{+})\left.{\delta\left(\hbar\omega^{-+}+\mu^{-+}\right)\over\hbar(z^{-+}+\omega^{e}(\hat{p}^{+})-\omega^{h}(-\hat{p}^{-}))}v^{i}(p)\right\\},\
\omega(k)\to 0.\ \ $ (S.72)

Since the condition $p,\ p^{-}\ll k/2$ holds, one can rewrite (S.72) as

$\displaystyle\lim_{\omega(k)\to
0}j_{i,2}^{Ohm}(\vec{r}^{+},\vec{r}^{-},t^{+}-t^{-})={\imath^{2}e^{2}A_{i}\over
2c(2\pi)^{2}}\mbox{Tr}{1\over-\imath\bar{\beta}\hbar}\sum_{\nu^{+-}}e^{-\imath
z^{-+}(t^{+}-t^{-})}\int{d^{2}\vec{k}\over(2\pi)^{2}}e^{-\imath\vec{k}\cdot\vec{r}}$
$\displaystyle\times\left\\{\int
d^{2}\vec{p}\left(1-\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c){d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}dp_{i}^{AB}}(0)\right)\right.$
$\displaystyle\times{v^{\dagger}}^{i}(p)\,{1\over-\imath\bar{\beta}\hbar}\sum_{\nu_{+}}G_{1}(z^{-+}-z^{+},p^{-})G_{1}(z^{+},p^{+})\left.{\delta\left(\hbar\omega^{-+}+\mu^{-+}\right)\over\hbar(z^{-+}+\omega^{e}(\hat{p}^{+})-\omega^{h}(-\hat{p}^{-}))}v^{i}(p)\right\\},\
\omega(k)\to 0.\ \ $ (S.73)

Performing summation of (S.69) and (S.73), we get the current in dynamical
regime in low-frequency limit $\omega(0)\to\omega(k),\
\omega(k)=v_{F}k/\hbar\to 0$, which coincides with the expression for minimal
dc-current (S.67).

###  XIX.7 Ohmic conductivity of the model $N=3$

In (S.67) summation of the product of one-particle Green’s functions on
Matsubara frequencies can be performed using the residue theorem. This theorem
allows us to replace the integration along the contour
$C^{\prime}=\sum_{\nu}C^{\prime}_{\nu}$, oriented counterclockwise and
surrounding the Matsubara frequencies $z_{\nu}$, to integration over another
contour $C$ oriented clockwise and surrounding the poles of functions $G_{1}$
Kadanov . Fermi-Dirac distribution function

$\displaystyle
f[\beta(z-\mu/\hbar)]=\left\\{\exp[\bar{\beta}\hbar(z-\mu/\hbar)]+1\right\\}^{-1}$
(S.74)

has the poles in the points
$z_{\nu}=\pi\nu/(\imath\hbar\bar{\beta})+\mu/\hbar$, $\nu=\pm 1,\ \pm 3,\ \pm
5,\ \ldots$ with residues $(\hbar\bar{\beta})^{-1}$. Therefore the sum
entering as a multiplier in (S.67), can be changed on a sum over poles
$\beta(z^{+}-z^{-})-\beta(E(p^{-})/\hbar)$ and $\beta(E(p^{+})/\hbar)$ of one-
particle Green function $G_{1}(z^{-+}-z^{+},p^{-})$ and $G_{1}(z^{+},p^{+})$
respectively:

$\displaystyle{1\over-\imath\bar{\beta}\hbar}\sum_{\nu_{+}}G_{1}(z^{-+}-z^{+},p^{-})G_{1}(z^{+},p^{+})=\imath\bar{\beta}\int_{C}{d[\beta(z)]\over
2\pi\imath}{f[\beta(z-\mu/\hbar)]\over\beta(z)-\beta(E(p^{+})/\hbar)}{1\over\beta(z^{-+})-\beta(z)-\beta(E(p^{-})/\hbar)}$
$\displaystyle=-\imath\bar{\beta}\int_{C}{d[\beta(z)]\over
2\pi\imath}f[\beta(z-\mu/\hbar)]{1\over\beta(E(p^{+})/\hbar)-\beta(z)}{1\over\beta(z^{-+})-\beta(z)-\beta(E(p^{-})/\hbar)}$
$\displaystyle=-\imath\bar{\beta}{f[\beta((E(p^{+})-\mu)/\hbar)]-f[\beta(z^{-+}-E(p^{-})-\mu/\hbar)]\over\beta(z^{-+})-\beta(E(p^{-})/\hbar)-\beta(E(p^{+})/\hbar)}=-\imath\bar{\beta}{f[\beta((H(p^{+})-\mu)/\hbar)]-f[\beta(H^{\dagger}(-p^{-})-\mu/\hbar)]\over\beta(z^{-+})-\beta(H(p^{+})/\hbar)+\beta(H^{\dagger}(-p^{-})/\hbar)}.$
(S.75)

Here $C$ is a counterclockwise contour. Taking into account of (S.75), we find
the Fourier image of the Ohmic current (S.67):

$\displaystyle j_{i}^{Ohm}(\omega^{-+},\ k)={\imath\bar{\beta}e^{2}A_{i}\over
c(2\pi)^{2}}\mbox{Tr}\left\\{\int
d^{2}\vec{p}\left(1-\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c){d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}dp_{i}^{AB}}(0)\right)\right.$
$\displaystyle\times{v^{\dagger}}^{i}(p)\,{f[\beta((H(p^{+})-\mu)/\hbar)]-f[\beta(H^{\dagger}(-p^{-})-\mu/\hbar)]\over\beta(z^{-+})-\beta(H(p^{+})/\hbar)+\beta(H^{\dagger}(-p^{-})/\hbar)}\left.{\delta\left(\hbar\omega^{-+}+\mu^{-+}\right)\over\hbar(z^{-+}+\omega^{e}(\hat{p}^{+})-\omega^{h}(-\hat{p}^{-}))}v^{i}(p)\right\\},\
\hbar\omega\ll v_{F}k\to 0.\ \ $ (S.76)

Due to (IV.3), the coefficient at $A_{i}$ entering the expression (S.76) after
its division on $c$, gives Ohmic conductivity of SM:

$\displaystyle\sigma_{ii}^{Ohm}(\omega^{-+},\
k)={\imath\bar{\beta}e^{2}\over(2\pi c)^{2}}\mbox{Tr}\left\\{\int
d^{2}\vec{p}\left(1-\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c){d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}dp_{i}^{AB}}(0)\right)\right.$
$\displaystyle\times{v^{\dagger}}^{i}(p)\,{f[\beta((H(p^{+})-\mu)/\hbar)]-f[\beta(H^{\dagger}(-p^{-})-\mu/\hbar)]\over\beta(z^{-+})-\beta(H(p^{+})/\hbar)+\beta(H^{\dagger}(-p^{-})/\hbar)}\left.{\delta\left(\hbar\omega^{-+}+\mu^{-+}\right)\over\hbar(z^{-+}+\omega^{e}(\hat{p}^{+})-\omega^{h}(-\hat{p}^{-}))}v^{i}(p)\right\\},\
\hbar\omega\ll v_{F}k\to 0.\ \ $ (S.77)

###  XIX.8  Polarization effects

Let us consider the influence of particle-hole production on the behavior of
the diagonal elements of complex conductivity $\sigma_{ll}$ of SM. This
contribution in $\sigma_{ll}$ is described by the polarization contribution to
current entering into the current (IV.6) in main text. The explicit expression
is obtained in a manner similar to previous calculations and reads

$\displaystyle j_{l}^{Zitterbew}=-{e^{2}A_{l}\over
c\,\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c)}\chi_{+\sigma_{{}_{B}}}^{\dagger}\chi_{+\sigma_{{}_{B}}}={\imath
e^{2}A_{l}\
\mbox{Tr}\over\imath^{3}c\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c)}$
$\displaystyle\times\left\\{\chi_{+\sigma_{{}_{B}}}^{\dagger}\chi_{+\sigma_{{}_{B}}}-{(-\imath)}\int\int\sum_{i}\left[{e\over
c}v^{i}_{x^{+}-x^{\prime},\,\bar{x}}A_{i}(\bar{x})+\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c){d\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}}(0)\ {v^{i}_{x^{+}-x^{\prime},\,\bar{x}}}\right.\right.$
$\displaystyle+\left.\left.{\widetilde{\Sigma_{AB}\Sigma_{BA}}^{2}(\vec{p}_{AB}-e\vec{A}/c)\over
2}\sum_{j}{d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}dp_{j}^{AB}}(0)\
{v^{i}_{x^{+}-x^{\prime},\,\bar{x}}\,v^{j}_{x^{+}-x^{\prime},\
\bar{x}}}\right]\right.$
$\displaystyle\times\left.G_{2}(x^{+},\,\bar{x},\,x^{\prime},\,x^{\prime\prime})G_{1}(x^{-}-x^{\prime\prime})d\bar{t}\,d^{2}\vec{\bar{x}}\,dt^{\prime}\,d^{2}\vec{x^{\prime}}\,dt^{\prime\prime}\,d^{2}\vec{x^{\prime\prime}}+\ldots\right\\}.\
\ \ $ (S.78)

$\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c)\approx{d\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}}p_{i}^{AB}$ near Majorana zero-energy state, and respectively the
first term in equation (S.78) is equal to zero:

$\displaystyle j_{l}^{Zitterbew}=-{e^{2}A_{l}\over
c\,\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c)}\chi_{+\sigma_{{}_{B}}}^{\dagger}\chi_{+\sigma_{{}_{B}}}={\imath
e^{2}A_{l}\over\imath^{3}c}$
$\displaystyle\times\mbox{Tr}\left\\{\imath\int\int\sum_{i}\left[{e\over
c}v^{i}_{x^{+}-x^{\prime},\,\bar{x}}\tilde{A}_{i}(\bar{x})+{d\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}}(0)\
v^{i}_{x^{+}-x^{\prime},\,\bar{x}}\right.\right.+\sum_{k}{d\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{k}^{AB}}(0){d\widetilde{\Sigma_{AB}\Sigma_{BA}}\over dp_{i}^{AB}}(0)\
v^{k}_{x^{+}-x^{\prime},\,\bar{x}}v^{i}_{x^{+}-x^{\prime},\,\bar{x}}$
$\displaystyle+\left.\left.{\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c)\over
2}\sum_{j}{d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}dp_{j}^{AB}}(0)\
{v^{i}_{x^{+}-x^{\prime},\,\bar{x}}\,v^{j}_{x^{+}-x^{\prime},\
\bar{x}}}\right]\right.$
$\displaystyle\times\left.G_{2}(x^{+},\,\bar{x},\,x^{\prime},\,x^{\prime\prime})G_{1}(x^{-}-x^{\prime\prime})d\bar{t}\,d^{2}\vec{\bar{x}}\,dt^{\prime}\,d^{2}\vec{x^{\prime}}\,dt^{\prime\prime}\,d^{2}\vec{x^{\prime\prime}}+\ldots\right\\}$
$\displaystyle={e^{2}A_{l}\over\imath
c}\mbox{Tr}\left\\{\int\int\sum_{i}\left[{\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c)\over
2}{d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over dp_{i}^{AB}dp_{i}^{AB}}(0)\
{v^{i}_{x^{+}-x^{\prime},\,\bar{x}}\,v^{j}_{x^{+}-x^{\prime},\
\bar{x}}}\right]\right.$
$\displaystyle\times\left.G_{2}(x^{+},\,\bar{x},\,x^{\prime},\,x^{\prime\prime})G_{1}(x^{-}-x^{\prime\prime})d\bar{t}\,d^{2}\vec{\bar{x}}\,dt^{\prime}\,d^{2}\vec{x^{\prime}}\,dt^{\prime\prime}\,d^{2}\vec{x^{\prime\prime}}+\ldots\right\\}.\
\ \ $ (S.79)

The current (S.79) has the form similar to Ohmic one and can be calculated in
the same way resulting with

$\displaystyle j_{l}^{Zitterbew}(\omega^{-+},\
k)={\imath\bar{\beta}e^{2}A_{l}\over c(2\pi)^{2}}\mbox{Tr}\left\\{\int
d^{2}\vec{p}\sum_{i=1}^{2}{\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c)\over
2}{d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over dp_{i}^{2}}(0)\right.$
$\displaystyle\times{v^{\dagger}}^{i}(p)\,{f[\beta((H(p^{+})-\mu)/\hbar)]-f[\beta(H^{\dagger}(-p^{-})-\mu/\hbar)]\over\beta(z^{-+})-\beta(H(p^{+})/\hbar)+\beta(H^{\dagger}(-p^{-})/\hbar)}\left.{\delta\left(\hbar\omega^{-+}+\mu^{-+}\right)\over\hbar(z^{-+}+\omega^{e}(\hat{p}^{+})-\omega^{h}(-\hat{p}^{-}))}v^{i}(p)\right\\},\
\hbar\omega\ll v_{F}k\to 0.\ \ $ (S.80)

In the expression (S.80) for $j_{l}^{Zitterbew}$, not only electron moving
along an applied electric field $\vec{E}$ contribute to but also electrons
diffusively moving in the direction which is perpendicular to $\vec{E}$. Due
to (IV.3), the coefficient at $A_{i}$ entering into the expression (S.80)
after its division by $c$ gives the polarization contribution to the
conductivity

$\displaystyle\sigma_{ll}^{Zitterbew}(\omega^{-+},\
k)={\imath\bar{\beta}e^{2}\over(2\pi\,c)^{2}}\mbox{Tr}\left\\{\int
d^{2}\vec{p}\sum_{i=1}^{2}{\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c)\over
2}{d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over dp_{i}^{2}}(0)\right.$
$\displaystyle\times{v^{\dagger}}^{i}(p)\,{f[\beta((H(p^{+})-\mu)/\hbar)]-f[\beta(H^{\dagger}(-p^{-})-\mu/\hbar)]\over\beta(z^{-+})-\beta(H(p^{+})/\hbar)+\beta(H^{\dagger}(-p^{-})/\hbar)}\left.{\delta\left(\hbar\omega^{-+}+\mu^{-+}\right)\over\hbar(z^{-+}+\omega^{e}(\hat{p}^{+})-\omega^{h}(-\hat{p}^{-}))}v^{i}(p)\right\\}.$
(S.81)

###  XIX.9  Magnetoelectric effects

The contribution of magnetoelectric effects to the expression (IV.6) for the
current is calculated in the same way, the result reads

$\displaystyle j_{2(1)}^{spin-orbit}=(-1)^{1(2)}{\imath e\over
2}\chi_{+\sigma_{{}_{B}}}^{\dagger}\sigma_{z}\,v^{1(2)}_{x^{+}x^{-}}\chi_{+\sigma_{{}_{B}}}=(-1)^{1(2)}{\imath\over
2}{(-1)\imath e\over(\imath)^{3}}\ $
$\displaystyle\times\mbox{Tr}\left\\{-{(-\imath)}\int\int\sum_{i}\left[{e\over
c}v^{i}_{x^{+}-x^{\prime},\,\bar{x}}A_{i}+\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c){d\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}}(0)\ {v^{i}_{x^{+}-x^{\prime},\,\bar{x}}}\right.\right.$
$\displaystyle+{1\over 2}\sum_{j}{d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{i}^{AB}dp_{j}^{AB}}(0)\left(\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c)v^{i}_{x^{+}-x^{\prime},\
\bar{x}}-{e\over c}A_{i}\right)$
$\displaystyle\left.\left.\times\left(\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c)v^{j}_{x^{+}-x^{\prime},\
\bar{x}}-{e\over
c}A_{j}\right)\right]\right.\left.G_{2}(x^{+},\,\bar{x},\,x^{\prime},\,x^{\prime\prime})G_{1}(x^{-}-x^{\prime\prime})d\bar{t}\,d^{2}\vec{\bar{x}}\,dt^{\prime}\,d^{2}\vec{x^{\prime}}\,dt^{\prime\prime}\,d^{2}\vec{x^{\prime\prime}}+\ldots\right\\}\sigma_{z}v^{1(2)}_{x^{\prime}x^{-}}$
$\displaystyle=(-1)^{1(2)}{\imath\over 2}{(-1)\imath e\over(\imath)^{3}}\
\mbox{Tr}\left\\{-{(-\imath)}\int\int\left[{d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{1(2)}^{AB}dp_{2(1)}^{AB}}(0)\left(\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c)v^{1(2)}_{x^{+}-x^{\prime},\
\bar{x}}\right)\right.\right.\left.\left(-{e\over c}A_{2(1)}\right)\right]$
$\displaystyle\times\left.G_{2}(x^{+},\,\bar{x},\,x^{\prime},\,x^{\prime\prime})G_{1}(x^{-}-x^{\prime\prime})d\bar{t}\,d^{2}\vec{\bar{x}}\,dt^{\prime}\,d^{2}\vec{x^{\prime}}\,dt^{\prime\prime}\,d^{2}\vec{x^{\prime\prime}}+\ldots\right\\}\sigma_{z}v^{1(2)}_{x^{\prime}x^{-}}$
$\displaystyle=(-1)^{1(2)}{\imath\over 2}{e^{2}\over\imath
c}A_{2(1)}\mbox{Tr}\left\\{\int\int{d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{1(2)}^{AB}dp_{2(1)}^{AB}}(0)\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c)v^{1(2)}_{x^{+}-x^{\prime},\
\bar{x}}\right.$
$\displaystyle\times\left.G_{2}(x^{+},\,\bar{x},\,x^{\prime},\,x^{\prime\prime})G_{1}(x^{-}-x^{\prime\prime})d\bar{t}\,d^{2}\vec{\bar{x}}\,dt^{\prime}\,d^{2}\vec{x^{\prime}}\,dt^{\prime\prime}\,d^{2}\vec{x^{\prime\prime}}+\ldots\right\\}\sigma_{z}v^{1(2)}_{x^{\prime}x^{-}}.\
\ \ \ \ \ \ $ (S.82)

Contrary to the polarization effects (S.79), the Hall current $j_{2(1)}^{spin-
orbit}$ (S.82) is defined by non-diagonal derivatives of the dynamical mass
$\widetilde{\Sigma_{AB}\Sigma_{BA}}$. Performing the similar to previous
calculation we find the explicit form for $j_{2(1)}^{spin-orbit}$:

$\displaystyle j_{2(1)}^{spin-orbit}(\omega^{-+},\ k)=(-1)^{1(2)}{\imath\over
2}{\imath\bar{\beta}e^{2}A_{2(1)}\over c(2\pi)^{2}}\mbox{Tr}\left\\{\int
d^{2}\vec{p}\
{\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c)}{d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{1(2)}^{AB}dp_{2(1)}^{AB}}(0)\right.$
$\displaystyle\times{v^{\dagger}}^{1(2)}(p)\,{f[\beta((H(p^{+})-\mu)/\hbar)]-f[\beta(H^{\dagger}(-p^{-})-\mu/\hbar)]\over\beta(z^{-+})-\beta(H(p^{+})/\hbar)+\beta(H^{\dagger}(-p^{-})/\hbar)}\left.{\delta\left(\hbar\omega^{-+}+\mu^{-+}\right)\over\hbar(z^{-+}+\omega^{e}(\hat{p}^{+})-\omega^{h}(-\hat{p}^{-}))}\sigma_{z}v^{1(2)}(p)\right\\}.$
(S.83)

Again, the coefficient at $A_{2(1)}$ entering into the expression (S.83) after
its division by $c$ gives the spin-orbit contribution to the conductivity

$\displaystyle\sigma_{12(21)}^{spin-orbit}(\omega^{-+},\
k)=(-1)^{1(2)}{\imath\over
2}{\imath\bar{\beta}e^{2}\over(2\pi\,c)^{2}}\mbox{Tr}\left\\{\int
d^{2}\vec{p}\
{\widetilde{\Sigma_{AB}\Sigma_{BA}}(\vec{p}_{AB}-e\vec{A}/c)}{d^{2}\widetilde{\Sigma_{AB}\Sigma_{BA}}\over
dp_{1(2)}^{AB}dp_{2(1)}^{AB}}(0)\right.$
$\displaystyle\times{v^{\dagger}}^{1(2)}(p)\,{f[\beta((H(p^{+})-\mu)/\hbar)]-f[\beta(H^{\dagger}(-p^{-})-\mu/\hbar)]\over\beta(z^{-+})-\beta(H(p^{+})/\hbar)+\beta(H^{\dagger}(-p^{-})/\hbar)}\left.{\delta\left(\hbar\omega^{-+}+\mu^{-+}\right)\over\hbar(z^{-+}+\omega^{e}(\hat{p}^{+})-\omega^{h}(-\hat{p}^{-}))}\sigma_{z}v^{1(2)}(p)\right\\}.$
(S.84)

###  XIX.10 Diagonalized representation of the Hamiltonian

We utilize the following approximation:

$\displaystyle\vec{v}_{AB}\approx{\partial H^{(0)}_{AB}\over\partial\vec{p}}$
(S.85)

where $H^{(0)}_{AB}$ is the unperturbed Hamiltonian of the problem (IV.1) in
main text without the dynamical mass term:

$\displaystyle H^{(0)}_{AB}=c\vec{\sigma}_{2D}^{BA}\cdot\vec{p}_{AB}.$ (S.86)

Let indexes $a,\ b=1,\ 2$ be for valence and conduction bands. We designate
the eigenvalues of the unperturbed electron Hamiltonian (S.86) in a band
$a(b)$ through $E^{e}_{a(b)}$. Eigenvalues of hole Hamiltonian
$H_{0}^{\dagger}$ in the band $a(b)$ are designated by $E^{h}_{a(b)}$. Due to
the validity of the condition (S.27), it is possible to use the following
expansion for $\omega^{\pm}$ :

$\displaystyle\hbar\omega^{\pm}\equiv\hbar\omega(p^{\pm})=\left\\{\begin{array}[]{ll}E^{e}_{a(b)}(p)+(-1)^{i}\omega(k)/2&\mbox{
for electrons}\\\ E^{h}_{a(b)}(-p)-(-1)^{i}\omega(-k)/2&\mbox{ for
holes},\end{array}\right.$ (S.89)

where indexes $i=1,\ 2$ are for energy level of the band $a(b)$ for electron
$"e"$ (hole $"h"$) with momentum $p^{-}$ before transition and electron (hole)
$p^{+}$ after transition respectively. The conditions (S.89) mean that

$\begin{split}\left|E^{e(h)}_{a(b)}(p^{+})-E^{e(h)}_{a(b)}(p^{-})\right|=\omega(k),\\\
[E^{e}_{a(b)}(p^{\pm})-E^{h}_{a(b)}(p^{\mp})]=E^{e}_{a(b)}(+p)-E^{h}_{a(b)}(-p).\end{split}$
(S.90)

Let us perform calculations in the representation where the Hamiltonian of the
problem is diagonal. The velocity operator $\vec{v}_{AB(BA)}$ should be
transformed by the similarity transformation of the form $S^{-1}\vec{v}_{AB}S$
with a matrix $S$ constructed on the eigenvectors $\chi$ of the Hamiltonian.
The corresponding matrix $S$ should be constructed on the eigenvectors of the
operator adjoined to the Hamiltonian. In every $p$-point the particle (hole)
Hamiltonian is represented by $2\times 2$ matrix, we denote matrix elements of
the exchange operator $\left(i\Sigma_{rel}^{x}\right)_{AB(BA)}$ formally as
$\Sigma^{AB(BA)}_{ij}$. Then the eigenvectors ${\chi}_{i},\ i=1,2$ of the
Hamiltonian (S.86) being the rows of the appropriate matrix $S$ can be
expressed in an explicit way:

$\chi_{1}=\left\\{\frac{ip\sin(\phi)(\Sigma^{AB}_{11}\Sigma^{AB}_{21}+\Sigma^{AB}_{12}\Sigma^{AB}_{22})+p\cos(\phi)(\Sigma^{AB}_{12}\Sigma^{AB}_{22}-\Sigma^{AB}_{11}\Sigma^{AB}_{21})-p(\Sigma^{AB}_{12}\Sigma^{AB}_{21}-\Sigma^{AB}_{11}\Sigma^{AB}_{22})}{p\left(\left((\Sigma^{AB}_{11})^{2}-(\Sigma^{AB}_{12})^{2}\right)\cos(\phi)-i\left((\Sigma^{AB}_{11})^{2}+(\Sigma^{AB}_{12})^{2}\right)\sin(\phi)\right)},\
1\right\\}$ (S.91)
$\chi_{2}=\left\\{\frac{ip\sin(\phi)(\Sigma^{AB}_{11}\Sigma^{AB}_{21}+\Sigma^{AB}_{12}\Sigma^{AB}_{22})+p\cos(\phi)(\Sigma^{AB}_{12}\Sigma^{AB}_{22}-\Sigma^{AB}_{11}\Sigma^{AB}_{21})+p(\Sigma^{AB}_{12}\Sigma^{AB}_{21}-\Sigma^{AB}_{11}\Sigma^{AB}_{22})}{p\left(\left((\Sigma^{AB}_{11})^{2}-(\Sigma^{AB}_{12})^{2}\right)\cos(\phi)-i\left((\Sigma^{AB}_{11})^{2}+(\Sigma^{AB}_{12})^{2}\right)\sin(\phi)\right)},\
1\right\\}$ (S.92)

In this way we can calculate numerically the velocity operator in every
$p$-point with subsequent its substitution to the conductivity integral.

###  XIX.11 Integral calculations

Contributions to conductivity include 2D integrals over the Brillouin zone
(BZ). For example, integrals in the Ohmic contribution given by formulae
(IV.15) have the following form:

$\displaystyle\sigma^{intra}_{ij}(\omega,\vec{k})=\sum\limits_{a=1,2}\frac{i{{e}^{2}}v_{0}^{2}}{{{\pi}^{2}}}\int\limits_{\mathrm{BZ}}{{{d}^{2}}\vec{p}}\frac{\left(\left[{{{{v}}}^{i}_{aa}}(\vec{p}){{{{v}}}^{j}_{aa}}(\vec{p})\right]f[{{\varepsilon}_{1}}(\vec{p}-\vec{k}/2)]-f[{{\varepsilon}_{1}}(\vec{p}+\vec{k}/2)]\right)}{\left({{\varepsilon}_{1}}(\vec{p}+\vec{k}/2)-{{\varepsilon}_{1}}(\vec{p}-\vec{k}/2)\right)\left(\omega-{{\varepsilon}_{1}}(\vec{p}+\vec{k}/2)+{{\varepsilon}_{1}}(\vec{p}-\vec{k}/2)\right)},$
(S.93)
$\displaystyle\sigma^{inter}_{ij}(\omega,\vec{k})=\frac{2i\omega{{e}^{2}}v_{0}^{2}}{{{\pi}^{2}}}\int\limits_{\mathrm{BZ}}{{{d}^{2}}\vec{p}}\frac{v_{12}^{i}(\vec{p})v_{21}^{i}(\vec{p})\left(f[{{\varepsilon}_{1}}(\vec{p}-\vec{k}/2)]-f[{{\varepsilon}_{2}}(\vec{p}+\vec{k}/2)]\right)}{\left({{\varepsilon}_{2}}(\vec{p}+\vec{k}/2)-{{\varepsilon}_{1}}(\vec{p}-\vec{k}/2)\right)\left({{\omega}^{2}}-{{\left({{\varepsilon}_{2}}(\vec{p}+\vec{k}/2)-{{\varepsilon}_{1}}(\vec{p}-\vec{k}/2)\right)}^{2}}\right)}.$
(S.94)

Here the first integral is for intraband transitions, the second one is for
interband ones. Now, we highlight the pole structure for the integrands for
small but finite $\vec{k}$, accounting that
$\epsilon_{1}(\vec{p})=-\epsilon_{2}(\vec{p})$. The first integral can be
rewritten as

$\displaystyle{{\sigma}^{intra}_{ij}}(\omega,\vec{k})=\sum\limits_{a=1,2}\frac{i{{e}^{2}}v_{0}^{2}}{{{\pi}^{2}}}\int\limits_{\mathrm{BZ}}{{{d}^{2}}\vec{p}\left[{{{{v}}}^{i}_{aa}}(\vec{p}){{{{v}}}^{j}_{aa}}(\vec{p})\right]}\,\,\left.\frac{df(\varepsilon)}{d\varepsilon}\right|_{\varepsilon=\varepsilon_{1}(\vec{p})}\frac{1}{(\omega-\vec{k}\cdot\nabla{{\varepsilon}_{1}}(\vec{p}))}\
,$ (S.95)

where we have performed the Taylor series expansion on $|\vec{k}|$ up to a
linear terms. In the second integral only the second term in denominator can
produce poles, so we expand it into a power series on $|\vec{k}|$, make a
change to polar coordinates $(p_{x},p_{y})\to(p,\phi)$ that results

$\displaystyle\sigma^{inter}_{ij}(\omega,\vec{k})=\frac{2i\omega{{e}^{2}}v_{0}^{2}}{{{\pi}^{2}}}\int\limits_{\mathrm{BZ}}{pdp\
d\phi}\frac{v_{12}^{i}(\vec{p})v_{21}^{i}(\vec{p})\left(f[{{\varepsilon}_{1}}(\vec{p}-\vec{k}/2)]-f[{{\varepsilon}_{2}}(\vec{p}+\vec{k}/2)]\right)}{\left({{\varepsilon}_{2}}(\vec{p}+\vec{k}/2)-{{\varepsilon}_{1}}(\vec{p}-\vec{k}/2)\right)\left(\omega^{2}-4p^{2}-k^{2}\sin^{2}\phi\right)}.$
(S.96)

Pole structures of (S.95) and (S.96) is presented in fig. S.2.

Figure S.2: Sketch of the poles structure of integrands for intra- (dashed
lines) and inter-(solid lines) bands contribution to conductivity. Solid lines
parameters at $\omega=1$ (in units of $v_{F}|\vec{K}_{A}|$) are
$k/|\vec{K}_{A}|=0.1$ (red), $0.7$ (green), $0.9$ (blue). Dashed lines
parameters at $\omega=0.2$ (in units of $v_{F}|\vec{K}_{A}|$) are
$k/|\vec{K}_{A}|=0.21$ (red), $0.4$ (green), $0.9$ (blue).

In accordance to (S.95) and fig. S.2 at $\vec{k}=0$ this integral is a regular
one, whereas at finite $k$ there is a line of poles (dashed lines in figure).
At a finite $k$ the pole structure of $\sigma^{inter}_{ij}(\omega,\vec{k})$
(S.96) is an elliptic one that results in necessity to account an infinite sum
of poles as contribution to conductivity:
$\sigma^{inter}_{ij}(\omega,\vec{k})\propto\int\ d\phi\mbox{ Res}({\phi})$,
where $\mbox{ Res}({\phi})$ is a residue in the pole located at angle $\phi$
on the poles line. The integral (S.96) in the case $k=0$ holds poles laying at
a circumference that can be effectively reduced to a single one as
$\sigma^{inter}_{ij}(\omega,\vec{k})\propto 2\pi\mbox{ Res}_{1}$, where
$\mbox{ Res}_{1}$ is a residue in arbitrary point of the circumference. For
very oblate ellipse at large $k$ the main contribution to the integral (S.96)
gives the only points touching the circumference
$\sigma^{inter}_{ij}(\omega,\vec{k})\sim 2\mbox{ Res}_{1}$. Thus, the value of
the optical conductivity decreases with the growth of $k$. In experiment the
optical transmission through a graphene monolayer for normal incidence and
respectively the optical conductivity really decreases in optical wave range
in direction of shorter wave lengths Nair2008Science ;
Falkovsky2011LowTemp.Phys .

We define the energy limit of applicability 1 eV$\approx 10^{4}$ K of tight-
binding approximation as a $\omega_{max}$ of the model with linear dispersion.
Then, the upper integration limit $q_{max}$ on momentum is
$q_{max}\sim\omega_{max}/v_{F}$, and respectively for the model of massless
Dirac fermions the integration should use the range from 0 to
$q<0.14~{}K_{A}-0.28~{}K_{A}$. As the simulation results presented in fig. S.3
demonstrate, the integration within this limits leads to the conductivity fall
in the range $\omega_{max}\sim 4000-8000$ K. But, experimentally such a fall
starts at much higher frequencies ($8000$ K) Nair2008Science ;
Falkovsky2011LowTemp.Phys . Thus, the range of momenta to predict conductivity
in visible optical range is outside the limits of applicability of the
massless Dirac fermion model.

(a) (b) (c)

Figure S.3: Frequency dependencies of the real and imaginary parts of the
massless Ohmic term of conductivity at very small wave number
$q=10^{-8}|\vec{K}_{A}|$ in rel. units $e^{2}/\hbar$; a cutting parameter
$q_{max}=0.14~{}K_{A}$. The model Falkovsky ; Falkovsky2008PhysUspekhi  is
simulated at T=3K, chemical potential $\mu=135$K (red dashed lines in figures
(a–b)) and at T=200K, $\mu=33$K (red dashed line in figure (c)). Numerical
results for our model are green-solid and black-dashed lines for the first-
order approximation (with zero gauge-phases) (III.9) at T=3K, $\mu=135$K and
T=200K, $\mu=33$K respectively, and blue line for the second-order
approximation (with non-zero gauge-phases) (III.19) at T=3K, $\mu=135$K in
figures (a–b). (c) Dependence of the real and imaginary parts of the massless
Ohmic term of conductivity on a damping $\gamma$ for the second-order
approximation (III.19) with $\gamma=0.1$ (green curve) and $\gamma=1$ (blue
curve) at T=200K, $\mu=33$K.

According to the Table 1, the Fermi velocity $v_{M}$ for our model of massless
Majorana energy states decrease in a factor higher than two for high wave
numbers $q$ in respect to its value in the Dirac point. Therefore, the range
$q$, $q<\omega_{max}/v_{M}$, $\omega_{max}<$1 eV increases in a factor two
compared to the massless Dirac fermions. Since $q<0.28~{}K_{A}-0.56~{}K_{A}$
for our model, the range of integration over momentum $|q|<0.44~{}|K_{A}|$ is
appropriate one and for it $\omega_{max}>7000$ K.

In the first-order approximation with zero-phases of the gauge fields, the
analytical formulas for the integrands in separate conductivity contribution
terms have been used. The integrals have been calculated with adaptive
integration steps in both directions ($|k|,\phi$) providing high calculation
accuracy (not less than 0.01%).

In the second order approximation with nonzero-phases of the gauge fields, we
have to calculate numerically by introducing into consideration a small
positive damping constant $\gamma$ for the states as a small imaginary
contribution to the energies. The values of $\gamma$ define the extend of
smoothness of the singular behaviour of integrand and does not influence on
the general form of the dependency curve in accordance to fig. S.3c.

All quantities necessary for calculation of the complex conductivity have been
calculated on a grid in the space of wave vectors with 200 discretization in
angle $\phi$ for every given wave number $q$ and variable step (a denser grid
at small wave numbers and larger at large ones) to the maximum wave number
$q_{max}=0.44|K_{A}|$. 2D interpolation on this grid has been used to
integrands evaluation in the intermediate points that is necessary for
conductivity simulations. An error stipulated by the interpolation from the
grid in wave vectors space has been roughly estimated by interpolation of the
conic spectrum of Dirac pseudo-fermion model on the same lattice with
subsequent usage of the interpolation data for evaluation of the conductivity.
Its value turns out to be less than $10^{-3}$% .

Total estimation of the conductivity error has been performed by variation of
the number of points used for interpolation of the energy band spectrum (by
diminishing this number at factor two and subsequent comparison of the
simulation results in both cases). It turns out to be not exceeding 10% in the
considered frequency region. It should be noted that the error bars for values
of the Fermi velocity which was measured by different techniques including
transport experiments (Shubnikov – de Haas oscillations) NatPhys7-2011Elias ,
infrared measurements of the Pauli blocking in graphene Z.Q.Li_et_al_Nature
Phys.4_532(2008) , magneto-optics magneto_optics_Fermi_velocity_NatCommun2014
, were also of the order of 10% .