# Hydrodynamics of charged two-dimensional Dirac systems II:
the role of collective modes

Kitinan Pongsangangan  Institute for Theoretical Physics and Center for
Extreme Matter and Emergent Phenomena, Utrecht University, Princetonplein 5,
3584 CC Utrecht, The Netherlands  Tim Ludwig  Institute for Theoretical
Physics and Center for Extreme Matter and Emergent Phenomena, Utrecht
University, Princetonplein 5, 3584 CC Utrecht, The Netherlands  Henk T.C.
Stoof  Institute for Theoretical Physics and Center for Extreme Matter and
Emergent Phenomena, Utrecht University, Princetonplein 5, 3584 CC Utrecht, The
Netherlands  Lars Fritz  Institute for Theoretical Physics and Center for
Extreme Matter and Emergent Phenomena, Utrecht University, Princetonplein 5,
3584 CC Utrecht, The Netherlands

###### Abstract

We study the hydrodynamic properties of ultraclean interacting two-dimensional
Dirac electrons with Keldysh quantum field theory. We study it from a weak-
coupling and a strong-coupling perspective. We demonstrate that long-range
Coulomb interactions play two independent roles: (i) they provide the
inelastic and momentum-conserving scattering mechanism that leads to fast
local equilibration; (ii) they facilitate the emergence of collective
excitations, for instance plasmons, that contribute to transport properties on
equal footing with electrons. Our approach is based on an effective field
theory of the collective field coupled to electrons. Within a conserving
approximation for the coupled system we derive a set of coupled quantum-
kinetic equations. This builds the foundation of the derivation of the
Boltzmann equations for the interacting system of electrons and plasmons. From
this, we explicitly derive all the conservation laws and identify the extra
contributions of energy density and pressure from the plasmons. We demonstrate
that plasmons show up in thermo-electric transport properties as well as in
quantities that enter the energy-momentum tensor, such as the viscosity. In a
parallel paper we discuss some of the phenomenology of the corresponding
hydrodynamic equations with an eye on thermo-electric transport properties.

##  I introduction

The conventional theory of electronic transport in a solid-state setting
describes the motion of electrons in the following way: individual electrons
diffuse on the background of a disordered lattice, primarily scattering from
impurities and/or lattice vibrations Drude1900; Ashcroft&Mermin1976; Ziman
(2000). However, recent years have seen tremendous progress in realizing
electronic transport that follows a different paradigm: hydrodynamic
electrons, meaning electrons flowing collectively like a viscous liquid, such
as water or honey. This idea has first been discussed in the 1960s by Gurzhi
Gurzhi (1963). However, the subject has only recently picked up a lot of pace.
This is mostly due to recent advances in the preparation of ultrapure mono-
and bilayer graphene samples with sufficiently strong interactions
Novoselov2005; Castro Neto et al. (2009); Bandurin et al. (2018, 2016);
Sulpizio et al. (2019); de Jong and Molenkamp (1995). The prerequisite for the
experimental observation of electron hydrodynamics is that microscopic
momentum-conserving electron-electron collisions due to Coulomb interactions
must be considerably faster than momentum-relaxing scatterings of electrons
from impurities and/or phonons. This allows to locally establish equilibrium.
In that situation, it is justified to speak of approximate conservation laws,
sufficient to open the door for the observation of electron hydrodynamics
Fritz et al. (2008); Müller et al. (2008); Landau (1936); Lucas2016; Lucas and
Fong (2018); Narozhny (2019). The theoretical method used to investigate this
hydrodynamic flow is usually based on the traditional Navier-Stokes equation
which expresses conservation of momentum, energy and electric charge Landau
and Lifshitz (1987); Lucas and Fong (2018); Narozhny et al. (2015). One way to
derive these hydrodynamic equations from microscopics starts from the
Boltzmann equation Lifshitz2981; Narozhny (2019). There, usually, the effect
of the Coulomb interactions enters only through electron-electron collisions
leading to local equilibrium. This effect is usually calculated from a weak-
coupling perspective within second-order perturbation theory in the
interaction strength, referred to as the Born approximation for the cross-
section Landau (1936); Fritz et al. (2008); Kashuba (2008); Kadanoff and Baym
(1962). However, hydrodynamics of a charged liquid must be expected to behave
differently in many respects from the archetypal example of fluid like water
flow. The reason is that for a charged electronic system, in addition to said
inelastic scatterings, Coulomb interactions provide also long-range mean-field
forces between electrons. This facilitates the emergence of collective
excitations, for instance plasma oscillations or plasmons Pines and Bohm
(1952); Bohm and Pines (1953); Vlasov1968. Under the correct circumstances, as
present in the situation studied here, these plasmons behave as ‘bona fide’
quasi-particles that have their own dynamics. This effect is usually not
considered in the context of electron hydrodynamics.

Background and main idea:  Strongly correlated many-particle systems can often
be regarded as a collection of weakly interacting excitations. One of the
prime examples of this kind is the Landau Fermi liquid Abrikosov1959. Its
excitations behave as well-defined entities, called quasi-particles. This
means, the following conditions must be fulfilled: (i) an excitation with
momentum $\vec{p}$ possesses a well-defined complex energy spectrum, say
$\omega(\vec{p})=\epsilon(\vec{p})-i\gamma(\vec{p})$, where the imaginary part
of the energy $\gamma(\vec{p})$ describes the decay rate of the particle,
inversely proportional to the lifetime of the excitation; (ii) an excitation
must be long lived which means that the decay must be at least underdamped
$\gamma(\vec{p})\ll\epsilon(\vec{p})$. These requirements usually limit one to
the low-frequency and/or long-wavelength behavior of the system. Under such
circumstances, it is justified to consider the complicated interacting system
as a collection of independent elementary quasi-particles Pines (1999). Once
this paradigm is adopted for a given many-body system, there are two main
questions in need of answer: (i) what are the quasi-particles involved in
physical phenomena of interest at the relevant energy scale and (ii) what is
the population of these quasi-particles in each momentum state especially when
the system is exposed to external disturbances? The latter is quantified by a
distribution function $f(\vec{x},t,\vec{p})$ which gives the probability
density that the particle in a momentum state $\vec{p}$ is found in a
neighborhood of a spatial position $\vec{x}$ at time $t$. For a given non-
equilibrium situation, the time evolution of the distribution function is
described by the Boltzmann equation. The knowledge of the distribution
function and the energy spectrum is vital to determine all thermodynamic as
well as transport properties of the system in a straightforward fashion. The
main question we address in this paper are: (i) what are the quasi-particles
involved in hydrodynamical transport phenomena of ultraclean two-dimensional
Dirac electrons at an accessible temperature and (ii) how do they equilibrate?
The conventional theory for hydrodynamic behavior in two-dimensional Dirac
systems is a version of a two-component hydrodynamics, consisting of electrons
and holes Lucas2016; Narozhny et al. (2015); Fritz et al. (2008). This
scenario is very popular and can formally be derived from a weak-coupling
analysis using a Hartree-Fock-Born approximation. In that framework, Coulomb
interaction plays two different roles: (i) it is seen by an electron as an
internal potential produced by all the other electrons in the system through
the Hartree term; (ii) Coulomb interaction manifests itself as an inelastic
and momentum conserving scattering mechanism Fritz et al. (2008); Müller et
al. (2008); Kashuba (2008); Kadanoff and Baym (1962); Landau (1936) which
locally equilibrates the system. This process is an important requirement for
observing the electron hydrodynamical regime with its fascinating transport
properties Lucas2016; Lucas and Fong (2018); Narozhny (2019).

Perturbation theory up to, in this case, second order in the coupling constant
is unfortunately unable to describe many important physical phenomena, such as
the emergence of collective modes. Collective modes, however, are the hallmark
of interacting electronic systems. One such mode are plasma oscillations, also
called plasmons. In conventional three-dimensional metals, plasmons are a
gapped degree of freedom with a large energy gap, larger than the Fermi
energy. This implies that thermal plasmons cannot be excited at realistic
experimental temperatures and hence are largely irrelevant, both for
thermodynamic as well as transport properties. As an example, aluminium at
room temperature has a ratio of the plasmon gap to the thermal energy
$\hbar\omega_{p}/k_{B}T\approx 16\;\text{eV}/0.25\;\text{eV}=64$.
Consequently, the plasmon occupation number is negligibly small,
$n_{B}(\omega_{p})\approx 10^{-28}$, where $n_{B}$ is the equilibrium Bose-
Einstein distribution function Pines (1953). In contrast, in two dimensions,
plasmons are massless and have a square-root dispersion, i.e.,
$\omega\propto\sqrt{q}$ Das Sarma and Li (2013); Wunsch et al. (2006); Stern
(1967). This is not only true for Dirac fermions but for a generic two-
dimensional electronic system. Consequently, at accessible temperatures,
plasmons can be excited and thus constitute proper low-energy elementary
excitations. There might be various effects that potentially destabilize the
plasmon and broaden their spectral function such as disorder, electron-phonon,
and electron-electron collision Glazman2004; Kivenson1969. However, it turns
out that the plasmons are typically remarkably stable.

The main result in this paper is that we offer a novel derivation of the
equations of hydrodynamics from a strong-coupling perspective. This results in
a combined description of electrons, holes, and plasmons, that are coupled to
each other. In this description, plasmons enter on equal footing with the
electronic degrees of freedom: we find that plasmons make a contribution of
the same order of magnitude as the fermionic degrees of freedom to heat
currents as well as the energy-momentum tensor and consequently should show up
in measurements that measure thermal transport but also viscous effects. The
results presented here provide a starting point for more phenomenological
descriptions, including transport properties, which we present in a parallel
paper.

Organization of the paper: We start with Sec. II where we introduce the model
of interacting Dirac fermions and describe how it connects, for instance, to
graphene. We then proceed to a technical section, Sec. III, where we review
the formalism of real-time quantum field theory. It is a summary of the most
important steps that lead from a fully quantum mechanical treatment towards
the semiclassical Boltzmann equation. This includes the partition function on
the Schwinger-Keldysh time contour, Sec. III.1, as well as the structure of
the Green functions on the closed time contour in Sec. III.2. We then discuss
the Dyson equation, including the Keldysh or quantum-kinetic equation in Sec.
III.3. We proceed to introduce the Wigner transform and the Moyal product in
Sec. III.4 which allows to perform the gradient expansion on the Keldysh
equation in Sec. III.5. More detailed accounts of the formalism can be found
in Refs. Kadanoff and Baym (1962); Kamenev2011; Rammer&Smith1986; Keldysh1965.
This section can be skipped by a reader familiar with the Schwinger-Keldsyh or
Kadanoff-Baym approach. In Sec. IV, we discuss two-component hydrodynamics as
currently used for the description of graphene. We derive all the equations
from the Schwinger-Keldysh approach and show how the two-component fluid
picture emerges within a weak-coupling approach to second order in the Coulomb
interaction. We start with introducing the non-interacting Green functions in
Sec. IV.1. Within weak coupling, Coulomb interactions play two different
roles: (i) It is seen by an electron as an internal potential produced by all
the other electrons in the system. This is the result of perturbation theory
to first order, called the Hartree-Fock scheme KMahan (2000), where the Fock
term leads to a renormalization of the Fermi velocity vozmediano; Elias2011.
The potential energy, on the other hand, is referred to as the Hartree
potential. In thermal equilibrium, it is canceled by a potential from the
underlying positively charged background in which the particles move. However,
it builds the basis for a derivation of the collective excitations of the
fluid Vlasov1968; Pines (1999). This aspect is discussed in Sec. IV.2; (ii)
The Coulomb interaction manifests itself as an inelastic and momentum
conserving scattering mechanism Fritz et al. (2008); Müller et al. (2008);
Kashuba (2008); Kadanoff and Baym (1962); Landau (1936). This process is an
important requirement for observing the electron hydrodynamical regime since
it leads to local equilibration. In Sec. IV.3 we discuss the scattering
process in detail. We then derive the ensuing conservation laws starting from
the Boltzmann equation in Sec. IV.4 and finish the section with a discussion
of collective modes in Sec. IV.5. We find a collective mode that we later
identify with the plasmon. We show explicitly that it matches with the strong-
coupling treatment if interpreted correctly.

Low-order perturbation theory is unable to describe many important physical
phenomena, for instance, collective modes, such as plasmons. We address this
situation in Sec. V. We start with a formalized version of the random-phase
approximation (RPA) in Sec. V.1. To this end, we introduce a new quantum field
associated with a plasmon excitation by means of a Hubbard-Stratonovich
transformation Kleinert1978; Hubbard1959, which is an exact rewriting of the
theory. After integrating out the fermions, the plasmons acquire their own
dynamics, which is discussed in Sec. V.2. The numerical details of this are
discussed in Sec. V.2.1 whereas in Sec. V.2.2 we discuss an analytical
approximation that allows us to make further progress. An important side
product is that we show that there is also a feedback effect, that
renormalizes the fermions and provides scattering for them, despite having
integrated them out at an earlier stage. Importantly, we can show explicitly
that this is not an instance of double counting, as one might suspect, but is
indeed required to preserve conservation laws. In Sec. V.3 we find a set of
coupled kinetic equations for electrons, holes, and the plasmons Wyld and
Pines (1962); Pongsangangan et al. (2020); PinesSchrieffer1962. It is
important that, within the conserving approximation, the electrons scatter
from plasmons and vice versa. Using a series of approximations, we derive
Boltzmann equations from this effective field theory. Based on this, we derive
the conservation laws of the system in Sec. V.4 which indeed show that the
fermion dynamics is strongly influenced by the plasmons. Furthermore, we
observe, that the plasmons make a sizable contribution to the heat current
(this provides an alternative derivation of the heat-current operator starting
from the quantum-kinetic equation) and energy density. Additionally, it makes
similar contributions to the momentum flux and therefore shows up in
quantities related to the viscosity. This section also confirms that RPA is
indeed a conserving approximation. We conclude our results in Sec. VI as well
as provide an outlook for future work. The most technical details of the
calculations are usually presented in a number of appendices.

##  II The model

We study a theory of charged Dirac electrons in two spatial dimensions
interacting via long-range Coulomb interactions. The Hamiltonian is given by

$\hat{H}=\hat{H}_{0}+\hat{H}_{\text{ex}}+\hat{H}_{I}\;,$ (1)

with $\hat{H}_{0}$ being the non-interacting part. In coordinate space, it
reads

$\hat{H}_{0}=\sum_{i=1}^{N}\sum_{\lambda,\lambda^{\prime}=\pm}\int
d\vec{x}\;\hat{\Psi}_{i,\lambda}^{\dagger}(\vec{x})\hat{H}_{D,\lambda\lambda^{\prime}}(\vec{x})\hat{\Psi}_{i,\lambda^{\prime}}(\vec{x})\;,$
(2)

Here $\hat{H}_{D,\lambda\lambda^{\prime}}(\vec{x})=\left(-i\hbar
v_{F}\vec{\sigma}\cdot\vec{\nabla}-\mu\right)_{\lambda\lambda^{\prime}}$ is
the Dirac Hamiltonian with a Fermi velocity $v_{F}$ and chemical potential
$\mu$, and $\hat{\Psi}^{\dagger}_{i,\lambda}(\vec{x})$
$(\hat{\Psi}_{i,\lambda}(\vec{x}))$ creates (annihilates) an electron at a
position $\vec{x}$. The flavor index denoted by $i$ range from $i=1,...,N$.
The symbols $\lambda,\lambda^{\prime}\in\\{+,-\\}$ denote spinor indices.

The static potential energy $V_{\rm{ex}}(\vec{x})$ is added in order to take
into account the positively charge background in which the electrons move,
i.e.,

$V_{\rm{ex}}(\vec{x})=-n_{0}\int
d\vec{x}^{\prime}V(\vec{x}-\vec{x}^{\prime})\;,$ (3)

where $n_{0}$ is the average density of the background ions, identical to the
electron density in thermal equilibrium. The interaction of the electrons with
the inert positively charged background is explained by

$\hat{H}_{\text{ex}}=\sum_{i=1}^{N}\sum_{\lambda=\pm}\int
d\vec{x}\;\hat{\Psi}^{\dagger}_{i,\lambda}(\vec{x})V_{\text{ex}}(\vec{x})\hat{\Psi}_{i,\lambda}(\vec{x}).$
(4)

In addition, $\hat{H}_{I}$ is the interaction part of the Hamiltonian. The
electrons interact via a long-range Coulomb interaction, which is included in
our model by the term

$\displaystyle\hat{H}_{I}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{i=1}^{N}\sum_{\lambda,\lambda^{\prime}=\pm}\int
d\vec{x}d\vec{x}^{\prime}\;\hat{\Psi}_{i,\lambda}^{\dagger}(\vec{x})\hat{\Psi}_{i,\lambda^{\prime}}^{\dagger}(\vec{x}^{\prime})V(\vec{x}-\vec{x}^{\prime})$
$\displaystyle\hskip
142.26378pt\hat{\Psi}_{i,\lambda^{\prime}}(\vec{x}^{\prime})\hat{\Psi}_{i,\lambda}(\vec{x})\;.$

Let us emphasize that we assume that the fermions of different flavours are
decoupled. The interaction potential between two electrons of charge $e$
separated by a distance $|\vec{x}-\vec{x}^{\prime}|$ is given by the
instantaneous Coulomb interaction

$V(\vec{x}-\vec{x}^{\prime})=\frac{e^{2}}{4\pi\epsilon|\vec{x}-\vec{x}^{\prime}|}\;,$
(6)

where $\epsilon$ is the average dielectric constant. After a Fourier
transformation, the non-interacting Hamiltonian acquires the form

$\hat{H}_{0}=\sum_{i=1}^{N}\sum_{\lambda,\lambda^{\prime}=\pm}\int\frac{d\vec{p}}{(2\pi)^{2}}\;\hat{\Psi}_{i,\lambda}^{\dagger}(\vec{p})\hat{\mathcal{H}}_{D,\lambda\lambda^{\prime}}(\vec{p})\hat{\Psi}_{i,\lambda^{\prime}}(\vec{p})\;,$
(7)

where $\hat{\mathcal{H}}_{D,\lambda\lambda^{\prime}}(\vec{p})=\left(\hbar
v_{F}\vec{\sigma}\cdot\vec{p}-\mu\right)_{\lambda\lambda^{\prime}}$ and the
interaction part of the Hamiltonian becomes

$\displaystyle\hat{H}_{I}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{i=1}^{N}\sum_{\lambda,\lambda^{\prime}=\pm}\int\frac{d\vec{p}}{(2\pi)^{2}}\frac{d\vec{k}}{(2\pi)^{2}}\frac{d\vec{q}}{(2\pi)^{2}}\;$
$\displaystyle\hat{\Psi}_{i,\lambda}^{\dagger}(\vec{k}-\vec{q})\hat{\Psi}_{i,\lambda^{\prime}}^{\dagger}(\vec{p}+\vec{q})V(\vec{q})\hat{\Psi}_{i,\lambda^{\prime}}(\vec{k})\hat{\Psi}_{i,\lambda}(\vec{p}).$

The Fourier transformation of the $1/r$ Coulomb interaction between electrons
in two dimensions reads

$V(\vec{p})=\frac{e^{2}}{2\epsilon p}=\frac{2\pi\alpha v_{F}}{p},$ (9)

where $p=|\vec{p}|$ denotes the norm of the two-dimensional momentum vector.
The strength of the Coulomb interaction is usually characterised by the ratio
of the potential energy to the kinetic energy. For the Dirac fermion, this
ratio boils down to the fine structure constant
$\alpha=e^{2}/4\pi\epsilon\hbar v_{F}$. The theory is identical for each
flavor $i=1,..,N$.

Let us mention that this model is applicable to low-energy interacting
electrons in graphene where momenta are measured with respect to the $K$
($K^{\prime}$) point in the Brillouin zone. In this case, $\lambda$ and
$\lambda^{\prime}$ denote the pseudo-spin degree of freedom taking into
account the $A$ and $B$ sub-lattices of the hexagonal lattice. The number of
fermion is then $N=4$ and it corresponds to spin and valley degrees of
freedom. It was found that $v_{F}\approx 10^{6}$ m/s, and
$\epsilon=\epsilon_{0}\epsilon_{r}$ measures the average value of the
dielectric constant of materials above and below it, e.g., $\epsilon=1$ for
suspended graphene in vacuum and $\epsilon\approx 7$ for graphene sandwiched
in hBN layers. Thus for these two cases, $\alpha=2.2$ and $\alpha\approx 0.3$,
respectively Narozhny (2019).

##  III The formalism

In this section we review the basics of non-equilibrium field theory. More
detailed accounts can be found in excellent books, Refs. Kadanoff and Baym
(1962); Kamenev2011; Rammer&Smith1986; Keldysh1965; CHOU1985. In Sec. III.1 we
start with a discussion of the generating function of the Green functions on
the Schwinger-Keldysh closed time contour. This builds the basis of the whole
approach. We then move on to discuss the structure of the Green function and
the Keldysh representation in Sec. III.2. In Sec. III.3 we discuss the Dyson
equation for interacting problems, as well as the quantum-kinetic or Keldysh
equation. In Sec. III.4, we discuss the Wigner transform, including the
gradient expansion and the Moyal product. Finally, in Sec. III.5, we discuss
the gradient expansion of the Keldysh equation, which is essential in the
derivation of the Boltzmann equation.

Figure 1: The closed time integration contour, $C=C_{+}\cup C_{-}$. The upper
branch $C_{+}$ goes forward from the initial time ($t=-\infty$) to the final
time ($t=\infty$) and the lower branch $C_{-}$ goes backward from the final
time to the initial time.

###  III.1 The generating function

The central object in the following Schwinger-Keldysh formalism is the
generating function. It allows for the derivation of Green functions and
associated physical observables by means of functional differentiation. In
this section we review the construction of the partition function for the
theory introduced in Eq. (1). The main idea is to assume that in the distant
past ($t=-\infty$), the system was in thermal equilibrium at a specific
temperature $T$. Its state is fully specified by the quantum-mechanical
density matrix $\hat{\rho}=e^{-\hat{H}_{0}/k_{B}T}$, where $\hat{H}_{0}$ is
the non-interacting Hamiltonian given by Eq. (2), and $k_{B}$ is the Boltzmann
constant. Henceforth, we use $\hbar=k_{B}=1$, unless stated otherwise. The
interaction will be switched on adiabatically to reach its actual strength
before the observation. In addition, external perturbations might be
subsequently established and drive the system away from equilibrium. The
partition function is defined as
$Z=\text{Tr}[\hat{U}_{C}\hat{\rho}]/\text{Tr}\left[\hat{\rho}\right]$ where
$\hat{U}_{C}=\mathcal{T}_{C}\exp\left(-i\oint_{C}\hat{H}(t)dt\right)$ is the
evolution operator along the closed time integration contour $C=C_{+}\cup
C_{-}$ depicted in Fig. 1. The operator $\mathcal{T}_{C}$ orders the operators
according to the position of their time arguments on the contour $C$. The
evolution operator describes the evolution of the system from $t=-\infty$,
where the system is non-interacting and in equilibrium, towards $t=\infty$,
and back to the equilibrium state at $t=-\infty$, that is again non-
interacting. During the evolution, the system may be exposed to external
perturbations, and its response can be examined. As in the conventional path-
integral approach, one discretizes the closed time contour into $M$
infinitesimal time intervals and inserts the resolution of unity in the
coherent-state basis at each discrete point in time along the contour.
Subsequently, by taking the continuum limit, we obtain the generating function
written as a functional integral according to

$Z=\int\mathcal{D}\psi^{\dagger}\mathcal{D}\psi\exp(iS[\psi^{\dagger},\psi])\;.$
(10)

In writing down the partition function, we absorb an irrelevant normalization
constant into the measure. The total action reads
$S[\psi^{\dagger},\psi]=\oint_{C}dt\mathcal{L}(\psi^{\dagger},\psi)$, with the
Lagrangian

$\displaystyle\mathcal{L}(\psi^{\dagger},\psi)=$ $\displaystyle\int
d\vec{x}\;\sum_{i=1}^{N}\sum_{\lambda=\pm}\psi_{i,\lambda}^{\dagger}(\vec{x},t)i\partial_{t}\psi_{i,\lambda}(\vec{x},t)-\hat{H}(\psi^{\dagger},\psi)\;.$

It is convenient to split the action into two parts,
$S[\psi^{\dagger},\psi]=\int_{C_{+}}dt\;\mathcal{L}(\psi^{\dagger}_{+},\psi_{+})+\int_{C_{-}}dt\;\mathcal{L}(\psi^{\dagger}_{-},\psi_{-})$.
It is important to note that, while the two parts seem independent, this is
not true due to the special boundary conditions when the contours meet, see
discussion in Ref. Kamenev2011. The first term describes the fermions on the
forward branch of the contour denoted by $C_{+}$ (the integration over the
time variable extends from $-\infty$ to $\infty$) whereas the second term
describes the fermions on the backward branch denoted by $C_{-}$. We
introduced the subscript $\pm$ for the fermion fields on the different
branches. Thus $\psi^{(\dagger)}_{+}$ and $\psi^{(\dagger)}_{-}$ denote the
fermion fields on the forward and backward branches of the contour,
respectively. Moreover, we interchange the limits of the time integration on
the backward branch, leading to an extra minus sign as the second integral
then also goes from $-\infty$ to $\infty$.

Eventually, we find that the action reads

$\displaystyle S[\psi^{\dagger},\psi]$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{N}\sum_{\lambda=\pm}\int
dtd\vec{x}\;\psi_{i,\lambda+}^{\dagger}(\vec{x},t)\Big{(}i\partial_{t}+iv_{F}\vec{\nabla}\cdot\vec{\sigma}+\mu-
V_{\rm{ex}}(\vec{x})\Big{)}\psi_{i,\lambda,+}(\vec{x},t)$
$\displaystyle-\frac{1}{2}\sum_{i=1}^{N}\sum_{\lambda,\lambda^{\prime}=\pm}\int
dtd\vec{x}d\vec{x}^{\prime}\;\psi_{i,\lambda,+}^{\dagger}(\vec{x},t)\psi_{i,\lambda^{\prime},+}(\vec{x},t)V(\vec{x}-\vec{x}^{\prime})\psi_{i,\lambda,+}(\vec{x}^{\prime},t)\psi_{i,\lambda^{\prime},+}^{\dagger}(\vec{x}^{\prime},t)$
$\displaystyle-\sum_{i=1}^{N}\sum_{\lambda=\pm}\int
d\vec{x}\;\psi_{i,\lambda,-}^{\dagger}(\vec{x},t)\Big{(}i\partial_{t}+iv_{F}\vec{\nabla}\cdot\vec{\sigma}+\mu-
V_{\rm{ex}}(\vec{x})\Big{)}\psi_{i,\lambda,-}(\vec{x},t)$
$\displaystyle+\frac{1}{2}\sum_{i=1}^{N}\sum_{\lambda,\lambda^{\prime}=\pm}\int
d\vec{x}d\vec{x}^{\prime}\;\psi_{i,\lambda,-}^{\dagger}(\vec{x},t)\psi_{i,\lambda^{\prime},-}(\vec{x},t)V(\vec{x}-\vec{x}^{\prime})\psi_{i,\lambda,-}(\vec{x}^{\prime},t)\psi_{i,\lambda^{\prime},-}^{\dagger}(\vec{x}^{\prime},t).$

###  III.2 Structure of the Green functions

For the ensuing discussion we neglect all internal indices such as flavor and
sublattice and only concentrate on the time arguments. The Green function on
the closed time contour is defined as

$iG(\vec{x},\vec{x}^{\prime},t,t^{\prime})=\langle\mathcal{T}_{C}\hat{\Psi}(\vec{x},t)\hat{\Psi}^{\dagger}(\vec{x}^{\prime},t^{\prime})\rangle,$
(13)

where $C$ and $\mathcal{T}_{C}$ have been introduced above. The operator at
the earliest time is arranged to the right-most position. There are in total
four different cases:

$\displaystyle iG_{-+}(\vec{x},\vec{x}^{\prime},t,t^{\prime})$
$\displaystyle=$
$\displaystyle\langle\hat{\Psi}(\vec{x},t)\hat{\Psi}^{\dagger}(\vec{x}^{\prime},t^{\prime})\rangle$
$\displaystyle iG_{+-}(\vec{x},\vec{x}^{\prime},t,t^{\prime})$
$\displaystyle=$
$\displaystyle-\langle\hat{\Psi}^{\dagger}(\vec{x}^{\prime},t^{\prime})\hat{\Psi}(\vec{x},t)\rangle$
$\displaystyle iG_{++}(\vec{x},\vec{x}^{\prime},t,t^{\prime})$
$\displaystyle=$
$\displaystyle\langle\mathcal{T}\hat{\Psi}(\vec{x},t)\hat{\Psi}^{\dagger}(\vec{x}^{\prime},t^{\prime})\rangle$
$\displaystyle iG_{--}(\vec{x},\vec{x}^{\prime},t,t^{\prime})$
$\displaystyle=$
$\displaystyle\langle\tilde{\mathcal{T}}\hat{\Psi}(\vec{x},t)\hat{\Psi}^{\dagger}(\vec{x}^{\prime},t^{\prime})\rangle.$
(14)

Here $iG_{-+}$ implies that the first time argument $t$ is on the branch
$C_{-}$ while the second time argument $t^{\prime}$ is on the branch $C_{+}$.
In this case, the operators are already arranged in a correct order, so
$\mathcal{T}_{C}$ can be dropped. Similarly,
$iG_{+-}(\vec{x},\vec{x}^{\prime},t,t^{\prime})$ means that the first time
argument $t$ lies on the forward branch whereas $t^{\prime}$ lies on the
backward branch. Since in this case $t^{\prime}$ is always further on the
contour than $t$, meaning $\mathcal{T}_{C}$ switches the annihilation and
creation operators together with giving an extra minus sign due to their
fermionic nature. Then $iG_{++}(\vec{x},\vec{x}^{\prime},t,t^{\prime})$ means
that both $t$ and $t^{\prime}$ are on the forward branch of the contour $C$,
and thus the time-contour ordering operator $\mathcal{T}_{C}$ becomes a normal
time-ordering operator $\mathcal{T}$. This implies that

$\displaystyle iG_{++}(\vec{x},\vec{x}^{\prime},t,t^{\prime})$
$\displaystyle=\theta(t-t^{\prime})iG_{-+}(\vec{x},\vec{x}^{\prime},t,t^{\prime})+\theta(t^{\prime}-t)iG_{+-}(\vec{x},\vec{x}^{\prime},t,t^{\prime})\;,$

where $\theta(t-t^{\prime})$ denotes the Heaviside theta function. In
contrast, for $iG_{--}(\vec{x},\vec{x}^{\prime},t,t^{\prime})$, both $t$ and
$t^{\prime}$ are on the backward branch. Since the direction of the backward
branch is opposite to the direction of the time axis (the backward branch
extends from $-\infty$ to $\infty$), the time-contour ordering operator
$\mathcal{T}_{C}$ becomes an anti-time-ordering operator, denoted by
$\tilde{\mathcal{T}}$. This means that

$\displaystyle iG_{--}(\vec{x},\vec{x}^{\prime},t,t^{\prime})$
$\displaystyle=\theta(t-t^{\prime})iG_{+-}(\vec{x},\vec{x}^{\prime},t,t^{\prime})+\theta(t^{\prime}-t)iG_{-+}(\vec{x},\vec{x}^{\prime},t,t^{\prime})\;.$

These components satisfy the following relation:

$\displaystyle
iG_{++}(\vec{x},\vec{x}^{\prime},t,t^{\prime})+iG_{--}(\vec{x},\vec{x}^{\prime},t,t^{\prime})$
$\displaystyle\hskip
28.45274pt=iG_{+-}(\vec{x},\vec{x}^{\prime},t,t^{\prime})+iG_{-+}(\vec{x},\vec{x}^{\prime},t,t^{\prime})\;.$
(17)

The Green functions can straightforwardly be calculated within the functional-
integral formalism according to

$\displaystyle iG_{AB}(\vec{x},\vec{x}^{\prime},t,t^{\prime})$
$\displaystyle=\frac{1}{Z}\int\mathcal{D}\psi^{\dagger}\mathcal{D}\psi\;\psi_{A}(\vec{x},t)\psi^{\dagger}_{B}(\vec{x}^{\prime},t^{\prime})\exp(iS[\psi^{\dagger},\psi]),$
(18)

where $A,B=\pm$ here label the branch index of the close-time contour. The
generating function $Z$ is defined in Eq. (10) together with the action in Eq.
(LABEL:theoryonthecontour). It is convenient to rotate the fields using the
so-called Larkin-Ovchinnikov transformation

$\displaystyle\begin{pmatrix}\psi_{+}(\vec{x},t)\\\
\psi_{-}(\vec{x},t)\end{pmatrix}=\frac{1}{\sqrt{2}}\begin{pmatrix}1&1\\\
1&-1\end{pmatrix}\begin{pmatrix}\psi_{1}(\vec{x},t)\\\
\psi_{2}(\vec{x},t)\end{pmatrix}\;\quad\rm{and}\quad\begin{pmatrix}\psi^{\dagger}_{+}(\vec{x},t)&\psi^{\dagger}_{-}(\vec{x},t)\end{pmatrix}=\frac{1}{\sqrt{2}}\begin{pmatrix}\psi^{\dagger}_{1}(\vec{x},t)&\psi_{2}^{\dagger}(\vec{x},t)\end{pmatrix}\begin{pmatrix}1&-1\\\
1&1\end{pmatrix}.$

After this transformation, Eq. (LABEL:theoryonthecontour) reads

$\displaystyle S[\psi^{\dagger},\psi]$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}dtdt^{\prime}\;\Big{[}\int
d\vec{x}d\vec{x}^{\prime}\psi_{a,\lambda,i}^{\dagger}(\vec{x},t)G^{-1}_{0,ab;\lambda\lambda^{\prime};ij}(\vec{x},\vec{x}^{\prime},t,t^{\prime})\psi_{b,\lambda^{\prime},j}(\vec{x}^{\prime},t^{\prime})$
(19) $\displaystyle-$ $\displaystyle\frac{1}{2}\int
d\vec{x}d\vec{x}^{\prime}\rho^{a}(\vec{x},t)D_{0,ab}(\vec{x},\vec{x}^{\prime},t,t^{\prime})\rho^{b}(\vec{x}^{\prime},t^{\prime})\Big{]}\;.$

Additionally, we have introduced a number of short-hand notations. The inverse
of the non-interacting Green function for the fermion field reads

$\displaystyle
G^{-1}_{0,ab;\lambda\lambda^{\prime};ij}(\vec{x},\vec{x}^{\prime};t,t^{\prime})=\delta_{ij}\delta_{ab}\times$
$\displaystyle\times\left(i\partial_{t}+iv_{F}\vec{\nabla}\cdot\vec{\sigma}+\mu-
V_{\rm{ex}}(\vec{x})\right)_{\lambda\lambda^{\prime}}\delta(\vec{x}-\vec{x}^{\prime})\delta(t-t^{\prime})\;.$
(20)

The Latin letters ($a,b,\dots$) take on the values $1$ and $2$ labeling the
Keldysh indices, the Greek letters ($\lambda,\lambda^{\prime},\dots$) assume
$\pm 1$ denoting the spinor indices, whereas $i,j=1,..,N$ denotes the flavor
indices. All double indices are summed over, unless stated otherwise.

The main advantages of the Larkin-Ovchinnikov basis is that the condition, Eq.
(17) is implemented and the remaining components are independent, meaning

$\displaystyle
iG_{ab;\lambda\lambda^{\prime}}(\vec{x},\vec{x}^{\prime},t,t^{\prime})$
$\displaystyle=$
$\displaystyle\left\langle\begin{pmatrix}\psi_{1}(\vec{x},t)\\\
\psi_{2}(\vec{x},t)\end{pmatrix}\begin{pmatrix}\psi^{\dagger}_{1}(\vec{x}^{\prime},t^{\prime})&\psi_{2}^{\dagger}(\vec{x}^{\prime},t^{\prime})\end{pmatrix}\right\rangle$
$\displaystyle=$
$\displaystyle\begin{pmatrix}iG^{R}(\vec{x},\vec{x}^{\prime},t,t^{\prime})&iG^{K}(\vec{x},\vec{x}^{\prime},t,t^{\prime})\\\
0&iG^{A}(\vec{x},\vec{x}^{\prime},t,t^{\prime})\end{pmatrix}_{ab}.$

The superscripts $R$, $A$ and $K$ stand for the retarded, advanced and Keldysh
components of the Green function, respectively. The Green functions are
calculable within the functional integral formalism according to

$\displaystyle iG_{ab}(\vec{x},\vec{x}^{\prime},t,t^{\prime})$
$\displaystyle=\int\mathcal{D}\psi^{\dagger}\mathcal{D}\psi\;\psi_{a}(\vec{x},t)\psi^{\dagger}_{b}(\vec{x}^{\prime},t^{\prime})\exp(iS[\psi^{\dagger},\psi]),$
(22)

with the action given in Eq. (19). It can be worked out that the inverse Green
function is given by

$\displaystyle G^{-1}_{ab}(\vec{x},\vec{x}^{\prime};t,t^{\prime})$
$\displaystyle\hskip
28.45274pt=\begin{pmatrix}(G^{-1})^{R}(\vec{x},\vec{x}^{\prime},t,t^{\prime})&(G^{-1})^{K}(\vec{x},\vec{x}^{\prime},t,t^{\prime})\\\
0&(G^{-1})^{A}(\vec{x},\vec{x}^{\prime},t,t^{\prime})\end{pmatrix}_{ab},$

where

$(G^{-1})^{R/A}=(G^{R/A})^{-1},$ (24)

and

$G^{R}\circ(G^{-1})^{K}=-G^{K}\circ(G^{-1})^{A}.$ (25)

The convolution operator $\circ$ is short for the integration over space and
time coordinates as well as the summation over the spinor indices. The Keldysh
component ($G_{12;\lambda\lambda^{\prime}}\equiv
G^{K}_{\lambda\lambda^{\prime}}$) is usually parametrized in terms of the
retarded and advanced components according to

$\displaystyle
G^{K}_{\lambda\lambda^{\prime}}(\vec{r},\vec{r^{\prime}};t,t^{\prime})\equiv
G^{R}\circ F-F\circ G^{A}.$ (26)

Here, $F$ is a Hermitian two-point function. By inserting this equation into
Eq. (25), we find

$(G^{-1})^{K}=(G^{R})^{-1}\circ F-F\circ(G^{A})^{-1}.$ (27)

This shows that the Keldysh component of the Green function and its inverse
can be parametrized in the same way.

Furthermore, we introduced the density operator
$\rho^{a}(\vec{x},t)=\psi_{b;\lambda,i}^{\dagger}(\vec{x},t)\gamma^{a}_{bc}\psi_{c;\lambda,i}(\vec{x},t)$
(all double indices are summed over) and the Coulomb interaction in Keldysh
space

$\displaystyle
D_{0,ab}(\vec{x},\vec{x}^{\prime},t,t^{\prime})=\sigma^{x}_{ab}\delta(t-t^{\prime})\frac{V(\vec{x}-\vec{x}^{\prime})}{2}\;.$
(28)

The vertex operators $\gamma$ are third rank tensors operating on the Keldysh
space of fermions as well as bosons. They are defined as
$\gamma^{1}_{ab}=1_{ab}$ and $\gamma^{2}_{ab}=\sigma^{x}_{ab}$. Let us
emphasize that our choice of convention gives an extra factor $1/2$ to the
Coulomb interaction term, i.e.,
$D_{0,ab}\propto\frac{V(\vec{x}-\vec{x}^{\prime})}{2}$.

###  III.3 The Dyson and the quantum-kinetic equation

In the presence of interactions, the Green function obeys the Dyson equation
(pictorially represented in Fig. 2)

$G=G_{0}+G_{0}\overset{\otimes}{}\Sigma\overset{\otimes}{}G.$ (29)

Figure 2: Diagrammatic representation of the Dyson equation for the fermion
Green function.

Here $\Sigma$ denotes the self-energy which is evaluated order by order from a
perturbative expansion in the interaction. It is represented diagrammatically
by a set of one-particle irreducible diagrams. The convolution $\otimes$
denotes integration over time and space, as well as summation of internal
indices, including the Keldysh indices. One can rewrite Eq. (29) as

$G^{-1}=G_{0}^{-1}-\Sigma.$ (30)

The fermion self-energy inherits the Keldysh matrix structure of the inverse
Green function as

$\displaystyle\Sigma_{ab;\lambda\lambda^{\prime}}(\vec{x},\vec{x}^{\prime},t,t^{\prime})$
$\displaystyle\hskip
28.45274pt=\begin{pmatrix}\Sigma^{R}_{\lambda\lambda^{\prime}}(\vec{x},\vec{x}^{\prime},t,t^{\prime})&\Sigma^{K}_{\lambda\lambda^{\prime}}(\vec{x},\vec{x}^{\prime},t,t^{\prime})\\\
0&\Sigma^{A}_{\lambda\lambda^{\prime}}(\vec{x},\vec{x}^{\prime},t,t^{\prime})\end{pmatrix}_{ab}.$
(31)

For the retarded and advanced components, it consequently enters according to

$(G^{-1})^{R/A}=(G_{0}^{-1})^{R/A}-\Sigma^{R/A}\;.$ (32)

The self-energy plays two major roles: It can affect the dispersion relation
and lead to a finite lifetime. This is most easily seen in frequency-momentum
space. The modified pole $\omega=\omega(\vec{p})-i\gamma(\vec{p})$ is the
solution of $(G_{0}^{-1})^{R}(\vec{p},\omega)-\Sigma^{R}(\vec{p},\omega)=0$.
It gives a new dispersion relation, $\omega(\vec{p})$, as well as a decay
rate, $\gamma(\vec{p})$, of the excitation. When the decay rate is
sufficiently small, this excitation is called a quasi-particle.

The Keldysh component of the Dyson equation leads to a kinetic equation that
governs the time evolution of the fermion distribution function. As a result,
we have

$\displaystyle-\Sigma^{K}$ $\displaystyle=$ $\displaystyle(G^{-1})^{K}$ (33)
$\displaystyle=$ $\displaystyle(G^{R})^{-1}\circ F-F\circ(G^{A})^{-1}\;,$

where we used the fact that $(G^{-1}_{0})^{K}$ is a pure regularization and
can be neglected in the presence of interactions, see Eq. (20). To arrive at
the second equality, we use Eq. (27) and parametrize the Keldysh component of
the inverse Green function in terms of the Hermitian function $F$.
Subsequently, by substituting Eq. (32) into Eq. (33), one finds

$G_{0}^{-1}\circ F-F\circ G_{0}^{-1}=-\Sigma^{K}+\Sigma^{R}\circ
F-F\circ\Sigma^{A}.$ (34)

The regularization $\pm i\delta$ can be omitted from the retarded and advanced
components in the presence of a non-zero imaginary part of the self-energy.
The above equation is called the quantum-kinetic equation for the distribution
matrix $F$. The solution of the full quantum kinetic equation, Eq. (34) is
usually exceedingly difficult. However, after some approximations, discussed
below, the quantum kinetic equation reduces to a Boltzmann equation. The
latter can be solved by, for instance, a variational method Ziman (2000).

Later in this paper, we show that the Coulomb interactions play a role in
facilitating the emergence of a boson associated with the electron density
fluctuations. To this end, we also summarize the salient features of the
Keldysh technique for a bosonic field. Similar to the fermion Green function,
the boson Green function and its inverse have three non-vanishing components
expressed in the following matrix structure

$\displaystyle D_{ab}(\vec{x},\vec{x}^{\prime},t,t^{\prime})$
$\displaystyle\hskip
28.45274pt=\begin{pmatrix}D^{K}(\vec{x},\vec{x}^{\prime},t,t^{\prime})&D^{R}(\vec{x},\vec{x}^{\prime},t,t^{\prime})\\\
D^{A}(\vec{x},\vec{x}^{\prime},t,t^{\prime})&0\end{pmatrix}_{ab},$ (35)

and

$\displaystyle D^{-1}_{ab}(\vec{x},\vec{x}^{\prime},t,t^{\prime})$
$\displaystyle\hskip
28.45274pt=\begin{pmatrix}0&(D^{-1})^{A}(\vec{x},\vec{x}^{\prime},t,t^{\prime})\\\
(D^{-1})^{R}(\vec{x},\vec{x}^{\prime},t,t^{\prime})&(D^{-1})^{K}(\vec{x},\vec{x}^{\prime},t,t^{\prime})\end{pmatrix}_{ab},$

together with the following relations

$(D^{-1})^{R/A}=(D^{R/A})^{-1},$ (37)

and

$D^{R}\circ(D^{-1})^{K}=-D^{K}\circ(D^{-1})^{A}.$ (38)

The interacting Green function is again determined from a Dyson equation
(pictorially represented in Fig. 3)

Figure 3: Diagrammatic representation of the Dyson equation for the boson
Green function.

$D=D_{0}+D_{0}\overset{\otimes}{}\Pi\overset{\otimes}{}D,$ (39)

with a self-energy $\Pi$, which is, in general, approximated in a perturbative
series. The self-energy assumes the same Keldysh matrix structure as the
inverse Green function to preserve causality, namely

$\displaystyle\Pi_{ab}(\vec{x},\vec{x}^{\prime},t,t^{\prime})$
$\displaystyle\hskip
28.45274pt=\begin{pmatrix}0&\Pi^{A}(\vec{x},\vec{x}^{\prime},t,t^{\prime})\\\
\Pi^{R}(\vec{x},\vec{x}^{\prime},t,t^{\prime})&\Pi^{K}(\vec{x},\vec{x}^{\prime},t,t^{\prime})\end{pmatrix}_{ab}.$
(40)

The poles of the retarded and advanced components are again shifted by
interaction effects leading to a new energy dispersion and a finite life time
for the dressed particles by means of

$(D^{-1})^{R/A}=(D_{0}^{-1})^{R/A}-\Pi^{R/A}.$ (41)

The Keldysh component of the Dyson equation leads to a kinetic equation for
the boson distribution function $B$.

$D_{0}^{-1}\circ B-B\circ D_{0}^{-1}=-\Pi^{K}+\Pi^{R}\circ B-B\circ\Pi^{A}.$
(42)

The left-hand side again describes a streaming term, whereas the right-hand
side accounts for collision events. The Hermitian function $B$ is employed to
parametrize the Keldysh Green function according to

$\displaystyle D^{K}=D^{R}\circ B-B\circ D^{A}.$ (43)

Before we continue, it is convenient to rewrite the Keldysh equations, Eqs.
(34) and (42). The rewriting seems arbitrary at this point, but later it will
allow us to identify the left-hand sides with the streaming terms of a
Boltzmann equation, whereas the right-hand sides will be the collision
integrals. We can use the fact that we can decompose the self-energies
according to $\Sigma^{R/A}=\Re\Sigma^{R}\pm i\Im\Sigma^{R}$ and
$\Pi^{R/A}=\Re\Pi^{R}\pm i\Im\Pi^{R}$. This allows to rewrite Eq. (34) as

$\left[G_{0}^{-1}-\Re\Sigma^{R}\overset{\circ}{,}F\right]=-\Sigma^{K}+i\\{\Im\Sigma^{R}\overset{\circ}{,}F\\}.$
(44)

and Eq. (42) as

$\left[D_{0}^{-1}-\Re\Pi^{R}\overset{\circ}{,}B\right]=-\Pi^{K}+i\\{\Im\Pi^{R}\overset{\circ}{,}B\\}.$
(45)

Here, we introduced the notation $[V\overset{\circ}{,}W]=V\circ W-W\circ V$
and $\\{V\overset{\circ}{,}W\\}=V\circ W+W\circ V$, defining the commutator
and anticommutator of the functions $V$ and $W$ with $\circ$ as defined
before.

###  III.4 The Wigner transform and the gradient expansion

In equilibrium quantum field theory, diagrammatic approaches are, because of
homogeneity, usually carried out in momentum and energy space instead of
coordinate-space and time, meaning it is a simple Fourier transform. The
semiclassical limit, however, is most conveniently accessed using the Wigner
transform, which is a mixed representation. We briefly summarize it here for
convenience. A generic space-time function,
$g(\vec{x}_{1},\vec{x}_{2},t_{1},t_{2})$, can be rewritten in terms of center-
of-mass ($\vec{r},t$) = ($(\vec{x}_{1}+\vec{x}_{2})/2,(t_{1}+t_{2})/2$) and
relative coordinates ($\vec{x},\tau$) =
($\vec{x}_{1}-\vec{x}_{2},t_{1}-t_{2}$). The Wigner transform is now a Fourier
transform over the relative coordinates while the centre-of-mass coordinates
are kept intact. Consequently, one obtains a function of center-of-mass
spacetime, momentum, and frequency, i.e.,
$\tilde{g}(\vec{r},t,\vec{p},\omega)=\int
d\vec{x}d\tau\;g(\vec{r},t,\vec{x},\tau)\;e^{-i\vec{p}\cdot\vec{x}+i\omega\tau}$.
There are two important Wigner transforms that will be needed later for our
derivation of the Boltzmann equation.

(i) For a two-point function which can be decompose into an algebraic product
of other two-point functions, i.e.,
$C(\vec{r}_{1},\vec{r}_{2};t_{1},t_{2})=A(\vec{r}_{1},\vec{r}_{2};t_{1},t_{2})B(\vec{r}_{1},\vec{r}_{2};t_{1},t_{2})$,
one can show that its Wigner transform is given by the momentum-frequency
convolution

$\displaystyle
C(\vec{r},t,\vec{p},\omega)=\int\frac{d\vec{p}_{1}}{(2\pi)^{d}}\frac{d\omega_{1}}{2\pi}\;A(\vec{r},t,\vec{p}_{1},\omega_{1})B(\vec{r},t,\vec{p}-\vec{p}_{1},\omega-\omega_{1}).$

(ii) The Wigner transform of the space-time convolution of two two-point
functions i.e., $D(\vec{r}_{1},\vec{r}_{2};t_{1},t_{2})=(A\circ
B)(\vec{r}_{1},\vec{r}_{2};t_{1},t_{2})$ is given by their Moyal product as

$\displaystyle D(\vec{r},t,\vec{p},\omega)$ $\displaystyle=$ $\displaystyle
A(\vec{r},t,\vec{p},\omega)\star
B(\vec{r},t,\vec{p},\omega)=A(\vec{r},t,\vec{p},\omega)\exp\left(\frac{i}{2}\left(\overleftarrow{\partial}_{\vec{r}}\overrightarrow{\partial}_{\vec{p}}-\overleftarrow{\partial}_{\vec{p}}\overrightarrow{\partial}_{\vec{r}}-\overleftarrow{\partial}_{t}\overrightarrow{\partial}_{\omega}+\overleftarrow{\partial}_{\omega}\overrightarrow{\partial}_{t}\right)\right)B(\vec{r},t,\vec{p},\omega)\;.$

In the cases we are interested in, the function varies slowly with the center-
of-mass coordinates. Consequently, it is legitimate to keep only the lowest-
order gradient terms according to

$\displaystyle D(\vec{r},t,\vec{p},\omega)$ $\displaystyle\approx$
$\displaystyle
A(\vec{r},t,\vec{p},\omega)B(\vec{r},t,\vec{p},\omega)-\frac{i}{2}\Big{[}\partial_{\vec{p}}A(\vec{r},t,\vec{p},\omega)\partial_{\vec{r}}B(\vec{r},t,\vec{p},\omega)-\partial_{\vec{r}}A(\vec{r},t,\vec{p},\omega)\partial_{\vec{p}}B(\vec{r},t,\vec{p},\omega)\Big{]}$
(48)
$\displaystyle+\frac{i}{2}\Big{[}\partial_{\omega}A(\vec{r},t,\vec{p},\omega)\partial_{t}B(\vec{r},t,\vec{p},\omega)-\partial_{t}A(\vec{r},t,\vec{p},\omega)\partial_{\omega}B(\vec{r},t,\vec{p},\omega)\Big{]}+...\;.$

###  III.5 Gradient expansion of the Keldysh equation

The Wigner transform of the Keldysh equation in Eq. (34), reads

$\displaystyle\left[G_{0}^{-1}-\Re\Sigma^{R}_{\star}\overset{\star}{,}F\right]_{-}$
$\displaystyle=$
$\displaystyle-\Sigma^{K}_{\star}+i\\{\Im\Sigma^{R}_{\star}\overset{\star}{,}F\\}_{+}\;.$
(49)

In writing the above equation, we briefly introduced the notation
$\Sigma_{\star}$. It accounts for the fact that the self-energy itself has a
gradient expansion according to

$\displaystyle\Sigma_{\star}\approx\Sigma+\Sigma_{\times}+...$ (50)

where $\Sigma$ involves no gradients, while $\Sigma_{\times}$ involves both
one spatial and momentum gradient (or, equivalently, frequency and time)
Kral1997; knoll2001. To leading order in non-vanishing gradients we find

$\displaystyle\left[G_{0}^{-1}-\Re\Sigma^{R},F\right]+i\Big{\\{}\partial_{\vec{x}}\left(G_{0}^{-1}-\Re\Sigma^{R}\right)\cdot\partial_{\vec{p}}F$
$\displaystyle-\partial_{\vec{p}}\left(G_{0}^{-1}-\Re\Sigma^{R}\right)\cdot\partial_{\vec{x}}F-\partial_{t}\left(G_{0}^{-1}-\Re\Sigma^{R}\right)\partial_{\epsilon}F$
$\displaystyle+\partial_{\epsilon}\left(G_{0}^{-1}-\Re\Sigma^{R}\right)\partial_{t}F\Big{\\}}=-\Sigma^{K}+i\\{\Im\Sigma^{R},F\\}$
$\displaystyle-\Sigma^{K}_{\times}+i\\{\Im\Sigma^{R}_{\times},F\\}.$ (51)

It turns out that the contributions due to $\Sigma_{\times}$ are of the non-
quasi-particle type and vanish once we perform the quasi-particle
approximation. Consequently, we drop them form our following discussion.

In order to set up a formalism that accommodates for a two- or even multi-band
scenario we make the following assumption: There is a transformation
$\mathcal{U}_{\vec{p}}$ (in the case of a Dirac-type theory we specify it
later) that projects the Green function into a diagonal basis according to

$\displaystyle g^{R/A}_{0}(\vec{x},\vec{p},\omega)$ $\displaystyle=$
$\displaystyle\mathcal{U}^{\dagger}_{\vec{p}}G^{R/A}_{0}(\vec{x},\vec{p},\omega)\mathcal{U}_{\vec{p}}$
(52)

with $g^{R/A}_{0}(\vec{x},\vec{p},\omega)$ being a diagonal matrix. While in
the following we present matrix equations, we only concentrate on the diagonal
elements. Equivalently, we project the self-energies into the quasi-particle
basis according to

$\displaystyle\sigma^{R/A}(\vec{x},\vec{p},\omega)$ $\displaystyle=$
$\displaystyle\mathcal{U}^{\dagger}_{\vec{p}}\Sigma^{R/A}(\vec{x},\vec{p},\omega)\mathcal{U}_{\vec{p}}\;.$
(53)

This leads to

$\displaystyle
i\partial_{\vec{x}}\left(g_{0}^{-1}-\Re\sigma^{R}\right)\cdot\left(\partial_{\vec{p}}F+i\left[\mathcal{A}_{\vec{p}},F\right]\right)$
$\displaystyle-\left(\partial_{\vec{p}}\left(g_{0}^{-1}-\Re\sigma^{R}\right)+i\left[\mathcal{A}_{\vec{p}},g_{0}^{-1}-\Re\sigma^{R}\right]\right)\cdot\partial_{\vec{x}}F$
$\displaystyle-\partial_{t}\left(g_{0}^{-1}-\Re\sigma^{R}\right)\partial_{\epsilon}F+\partial_{\epsilon}\left(g_{0}^{-1}-\Re\sigma^{R}\right)\partial_{t}F$
$\displaystyle\hskip 128.0374pt=-\sigma^{K}+2i\Im\sigma^{R}F\;.$ (54)

We observe that two additional terms that involve the Berry connection
$\mathcal{A}_{\vec{p}}=-i\mathcal{U}^{\dagger}_{\vec{p}}\partial_{\vec{p}}\mathcal{U}_{\vec{p}}$
are obtained. We are interested here in the solution to zeroth order in
$\mathcal{A}_{\vec{p}}$ that is

$F_{\lambda\lambda^{\prime}}=F_{\lambda\lambda^{\prime}}^{(0)}+\mathcal{O}(\mathcal{A}_{\vec{p}})=(1-2f_{\lambda}(\vec{r},\vec{p},t))\delta_{\lambda\lambda^{\prime}}+\mathcal{O}(\mathcal{A}_{\vec{p}}),$
(55)

containing only the diagonal elements of the distribution function. The
function $f$ introduced above will later play the role of the fermionic
distribution function and in equilibrium, it reduces to the Fermi-Dirac
distribution. Within the quasi-particle approximation, $f_{\lambda}$ is
independent of the frequency variable. Moreover, it is assumed to be diagonal
in the spinor space. Consequently, only the main diagonal elements of the
self-energies are important. In total, this leads to

$\displaystyle\partial_{\epsilon}\left(g_{0,\lambda\lambda}^{-1}(\vec{x},t,\vec{p},\epsilon)-\Re\sigma^{R}_{\lambda\lambda}(\vec{x},t,\vec{p},\epsilon)\right)\partial_{t}f_{\lambda}(\vec{x},t,\vec{p})$
$\displaystyle+\;\partial_{\vec{x}}\left(g_{0,\lambda\lambda}^{-1}(\vec{x},t,\vec{p},\epsilon)-\Re\sigma^{R}_{\lambda\lambda}(\vec{x},t,\vec{p},\epsilon)\right)\cdot\partial_{\vec{p}}f_{\lambda}(\vec{x},t,\vec{p})$
$\displaystyle-\;\partial_{\vec{p}}\left(g_{0,\lambda\lambda}^{-1}(\vec{x},t,\vec{p},\epsilon)-\Re\sigma^{R}_{\lambda\lambda}(\vec{x},t,\vec{p},\epsilon)\right)\cdot\partial_{\vec{x}}f_{\lambda}(\vec{x},t,\vec{p})$
$\displaystyle=-\frac{i}{2}\sigma^{K}_{\lambda\lambda}(\vec{x},t,\vec{p},\epsilon)-\Im\sigma^{R}(\vec{x},\vec{p},t,\epsilon)(1-2f(\vec{x},t,\vec{p}))\;.$
(56)

The left-hand side will contain the so-called streaming terms, consisting of
three contributions. The first term describes the time derivative of the
distribution function with the quasi-particle weight
$\partial_{\epsilon}\left(g_{0}^{-1}-\Re\sigma^{R}\right)=1-\partial_{\epsilon}\Re\sigma^{R}$.
The second term accounts for the change of the distribution function due to a
force $\partial_{\vec{x}}\left(g_{0}^{-1}-\Re\sigma^{R}\right)$, whereas the
last term tracks the change of the distribution function due to the diffusion
of excitations with the velocity
$\partial_{\vec{p}}\left(g_{0}^{-1}-\Re\sigma^{R}\right)$. The right-hand side
describes collisions that drive the system towards equilibrium.

We now consider the Keldysh equation in Eq. (42) of the plasmon field and we
drop the terms from expanding the self-energies here from the very start.
Following the same steps as for the fermions, we first decompose the retarded
and advanced components of the self-energy according to
$\Pi^{R/A}=\Re\Pi^{R}\pm i\Im\Pi^{R}$. This allows us to write Eq. (42) as

$\left[D_{0}^{-1}-\Re\Pi^{R}\overset{\circ}{,}B\right]=-\Pi^{K}+i\\{\Im\Pi^{R}\overset{\circ}{,}B\\}.$
(57)

After a Wigner transformation and keeping the Moyal product to first non-
trivial order, we obtain

$\displaystyle
i\Big{\\{}\partial_{\vec{x}}\left(D_{0}^{-1}-\Re\Pi^{R}\right)\cdot\partial_{\vec{p}}B-\partial_{\vec{p}}\left(D_{0}^{-1}-\Re\Pi^{R}\right)\cdot\partial_{\vec{x}}B$
$\displaystyle-\partial_{t}\left(D_{0}^{-1}-\Re\Pi^{R}\right)\partial_{\omega}B+\partial_{\omega}\left(D_{0}^{-1}-\Re\Pi^{R}\right)\partial_{t}B\Big{\\}}$
$\displaystyle\hskip 128.0374pt=-\Pi^{K}+2i\Im\Pi^{R}B.$

Next, we use the parametrization

$B=1+2b\;,$ (59)

where $b$ plays the role of the bosonic distribution function exactly in
equilibrium reduces to the Bose-Einstein distribution. This leads to

$\displaystyle\partial_{\epsilon}\left(D_{0}^{-1}(\vec{x},t,\vec{p},\epsilon)-\Re\Pi^{R}(\vec{x},t,\vec{p},\epsilon)\right)\partial_{t}b(\vec{x},t,\vec{p})$
$\displaystyle+\;\partial_{\vec{x}}\left(D_{0}^{-1}(\vec{x},t,\vec{p},\epsilon)-\Re\Pi^{R}(\vec{x},t,\vec{p},\epsilon)\right)\cdot\partial_{\vec{p}}b(\vec{x},t,\vec{p})$
$\displaystyle-\;\partial_{\vec{p}}\left(D_{0}^{-1}(\vec{x},t,\vec{p},\epsilon)-\Re\Pi^{R}(\vec{x},t,\vec{p},\epsilon)\right)\cdot\partial_{\vec{x}}b(\vec{x},t,\vec{p})$
$\displaystyle=\frac{i}{2}\Pi^{K}(\vec{x},t,\vec{p},\epsilon)+\Im\Pi^{R}(\vec{x},\vec{p},t,\epsilon)(1+2b(\vec{x},t,\vec{p}))\;.$
(60)

The missing piece to transform Eqs. (56) and (60) into Boltzmann equations is
to integrate them over the respective spectral functions, as we will discuss
later on.

##  IV Part A: Electron-hole hydrodynamics in the weak-coupling limit

In this section we review the equations of hydrodynamics in the weakly
interacting limit. We first discuss the non-interacting limit to define the
basic Green function and the projection into the appropriate quasi-particle
basis of electrons and holes. Afterwards we discuss the Hartree-Fock
approximation. We show that this is equivalent to the renormalized
collisionless Boltzmann equation after a series of approximations.

###  IV.1 The non-interacting limit

We consider a model of interacting Dirac fermions according to Eq. (19) where
we first neglect the interactions. By inserting Eq. (20) into Eq. (24), we
find the retarded Green function of the non-interacting Dirac theory according
to

$\displaystyle G^{R/A}_{0}(\vec{x},\vec{p},\omega)=\left(\omega\pm i\delta-
v_{F}\vec{\sigma}\cdot{\vec{p}}+\mu-V_{\rm{ex}}(\vec{x})\right)^{-1}.$

It is convenient to work in the quasi-particle basis. There, the Dirac
Hamiltonian as well as the Green function are diagonal. This leads to

$\displaystyle g^{R/A}_{0,\lambda\lambda^{\prime},ij}(\vec{x},\vec{p},\omega)$
$\displaystyle=$
$\displaystyle\left(\mathcal{U}^{\dagger}_{\vec{p}}G^{R/A}_{0,ij}(\vec{x},\vec{p},\omega)\mathcal{U}_{\vec{p}}\right)_{\lambda\lambda^{\prime}}$
$\displaystyle=$ $\displaystyle(\omega\pm i\delta-\lambda v_{F}p+\mu-
V_{\rm{ex}}(\vec{x}))^{-1}\delta_{\lambda\lambda^{\prime}}\delta_{ij}\;,$

with

$\displaystyle\mathcal{U}_{\vec{p}}=\frac{1}{\sqrt{2}}\begin{pmatrix}-\exp(-i\theta_{\vec{p}})&\exp(-i\theta_{\vec{p}})\\\
1&1\end{pmatrix}\;\;.$ (63)

Here, $p=|\vec{p}|$ and $\tan(\theta_{\vec{p}})=p_{y}/p_{x}$. The dispersion
relation can be extracted from the poles of the Green function in Eq. (IV.1):
the non-interacting Dirac theory has two linear dispersing energy bands
$\epsilon_{\pm}(\vec{x},\vec{p})=\pm v_{F}p+V_{\rm{ex}}(\vec{x})-\mu$ with the
local electro-chemical potential given by $\mu-V_{\rm{ex}}(\vec{x})$. The
spectral function, defined by $\Im
g^{R}_{0,\lambda\lambda^{\prime}}(\vec{x},\vec{p},\omega)=-\pi\delta(\omega-\epsilon_{\lambda}(\vec{x},\vec{p}))\delta_{\lambda\lambda^{\prime}}$,
exhibits resonances at $\omega=\epsilon_{\pm}(\vec{x},\vec{p})$. The Wigner
transform of the Keldysh Green function of Eq. (26) reads

$\displaystyle
g^{K}_{0,\lambda\lambda^{\prime}}(\vec{x},t,\vec{p},\omega)=2i\Im
g^{R}_{0,\lambda\lambda^{\prime\prime}}(\vec{x},\vec{p},\omega)F_{\lambda^{\prime\prime}\lambda^{\prime}}(\vec{x},\vec{p},\omega)$
$\displaystyle\hskip 14.22636pt=-2\pi
i\delta(\omega-\epsilon_{\lambda}(\vec{x},\vec{p}))\left(1-2f_{\lambda}(\vec{x},\vec{p},t)\right)\delta_{\lambda\lambda^{\prime}}.$

where $f(\vec{x},\vec{p},t)$ is the Fermi Dirac distribution in equilibrium.
The Keldysh component thus contains the information about the occupation
numbers, whereas the retarded and advanced components only contain information
about the resonances and the energy levels.

###  IV.2 Hartree-Fock approximation: The collisionless limit

In what follows, we will discuss the corrections of the energy spectrum due to
the Coulomb interaction. To this end, we study the Hartree and Fock self-
energies Kotov2012 depicted in Fig. 4(a) and Fig. 4(b).

####  IV.2.1 The Hartree diagram

The Hartree diagram in Fig. 4(a), has the following algebraic expression:

$\displaystyle-i\Sigma^{H}_{ab;\lambda\lambda^{\prime}}(\vec{x},\vec{x}^{\prime};t,t^{\prime})=-\delta_{\lambda\lambda^{\prime}}\delta(\vec{x}-\vec{x}^{\prime})\delta(t-t^{\prime})\int
dt^{\prime\prime}d\vec{x}^{\prime\prime}\;\gamma^{e}_{ab}D_{0,ef}(\vec{x}^{\prime},\vec{x}^{\prime\prime},t^{\prime},t^{\prime\prime})\gamma^{f}_{cd}G_{0,dc;\lambda^{\prime\prime}\lambda^{\prime\prime};ii}(\vec{x}^{\prime\prime},\vec{x}^{\prime\prime},t^{\prime\prime},t^{\prime\prime})\;.$
(65)

Since
$D_{0,ef}(\vec{x},\vec{x}^{\prime},t,t^{\prime})=\sigma_{ef}^{x}\delta(t-t^{\prime})\frac{V(\vec{x}-\vec{x}^{\prime})}{2}$
is off-diagonal in Keldysh space, only two terms survive when we sum over $e$
and $f$. This gives

$\displaystyle-i\Sigma^{H}_{ab;\lambda\lambda^{\prime}}(\vec{x},\vec{x}^{\prime};t,t^{\prime})$
$\displaystyle=$
$\displaystyle-\delta_{\lambda\lambda^{\prime}}\delta(\vec{x}-\vec{x}^{\prime})\delta(t-t^{\prime})\int
dt^{\prime\prime}d\vec{x}^{\prime\prime}\frac{V(\vec{x}^{\prime}-\vec{x}^{\prime\prime})}{2}\gamma^{1}_{ab}\gamma^{2}_{cd}G_{0,dc;\lambda^{\prime\prime}\lambda^{\prime\prime},ii}(\vec{x}^{\prime\prime},\vec{x}^{\prime\prime},t^{\prime\prime},t^{\prime\prime})$
$\displaystyle-$
$\displaystyle\delta_{\lambda\lambda^{\prime}}\delta(\vec{x}-\vec{x}^{\prime})\delta(t-t^{\prime})\int
dt^{\prime\prime}d\vec{x}^{\prime\prime}\frac{V(\vec{x}^{\prime}-\vec{x}^{\prime\prime})}{2}\gamma^{2}_{ab}\gamma^{1}_{cd}G_{0,dc;\lambda\lambda,ii}(\vec{x}^{\prime\prime},\vec{x}^{\prime\prime},t^{\prime\prime},t^{\prime\prime})\;.$

We find that the second term vanishes due to

$\displaystyle\gamma^{2}_{ab}\gamma^{1}_{cd}G_{0,dc;\lambda^{\prime\prime}\lambda^{\prime\prime},ii}(\vec{x}^{\prime},\vec{x}^{\prime},t^{\prime},t^{\prime})$
$\displaystyle=$ $\displaystyle
G^{R}_{0;\lambda^{\prime\prime}\lambda^{\prime\prime},ii}(\vec{x}^{\prime},\vec{x}^{\prime},t^{\prime},t^{\prime})+G^{A}_{0;\lambda^{\prime\prime}\lambda^{\prime\prime},ii}(\vec{x}^{\prime},\vec{x}^{\prime},t^{\prime},t^{\prime})$
(67) $\displaystyle=$
$\displaystyle\int\frac{d\omega}{2\pi}\frac{d\vec{p}}{(2\pi)^{2}}\Re
G^{R}_{0;\lambda^{\prime\prime}\lambda^{\prime\prime},ii}(\vec{x}^{\prime},t^{\prime},\vec{p},\omega)=0\;.$

The last equality follows because $\Re
G^{R}_{0;\sigma\sigma}(\vec{x}^{\prime},t^{\prime},\vec{p},\omega)$ is an
analytic function in either the upper or the lower complex frequency half-
plane. Consequently, we have

$\displaystyle-i\Sigma^{H}_{ab;\lambda\lambda^{\prime}}(\vec{x},\vec{x}^{\prime},t,t^{\prime})=-\delta_{\lambda\lambda^{\prime}}\delta(\vec{x}-\vec{x}^{\prime})\delta(t-t^{\prime})\delta_{ab}$
$\displaystyle\int
d\vec{x}^{\prime\prime}\frac{V(\vec{x}^{\prime}-\vec{x}^{\prime\prime})}{2}G^{K}_{0;\lambda^{\prime\prime}\lambda^{\prime\prime},ii}(\vec{x}^{\prime\prime},\vec{x}^{\prime\prime};t,t).$

Finally, after a Wigner transformation, we rotate the self-energy into the
quasi-particle basis. This gives

$\displaystyle-i\sigma^{H}_{ab;\lambda\lambda^{\prime}}(\vec{x},t,\vec{p},\omega)=-i\left(\mathcal{U}^{\dagger}_{\vec{p}}\Sigma^{H}(\vec{x},t,\vec{p},\omega)\mathcal{U}_{\vec{p}}\right)_{\lambda\lambda^{\prime}}$
$\displaystyle=$
$\displaystyle-\delta_{\lambda\lambda^{\prime}}\delta_{ab}\int
d\vec{x}^{\prime\prime}\frac{V(\vec{x}^{\prime}-\vec{x}^{\prime\prime})}{2}g^{K}_{0;\lambda^{\prime\prime}\lambda^{\prime\prime},ii}(\vec{x}^{\prime\prime},\vec{x}^{\prime};t,t)\;.$

The Wigner transformation of the Keldysh Green function with the same time and
spatial argument is associated with the electron density according to

$g_{0;\lambda^{\prime\prime}\lambda^{\prime\prime},ii}^{K}(\vec{x}^{\prime\prime},\vec{x}^{\prime\prime};t,t)=2in(\vec{x}^{\prime\prime},t),$
(70)

where
$n(\vec{x}^{\prime\prime},t)=n_{+}(\vec{x}^{\prime\prime},t)+n_{-}(\vec{x}^{\prime\prime},t)$
defines the total charge density at position $\vec{x}^{\prime\prime}$ and time
$t$. Here
$n_{+}(\vec{x}^{\prime},t)=N\int\frac{d\vec{k}}{(2\pi)^{2}}f_{+}(\vec{x}^{\prime},\vec{k},t)$
is the electron density and
$n_{-}(\vec{x}^{\prime},t)=N\int\frac{d\vec{k}}{(2\pi)^{2}}\left(f_{-}(\vec{x}^{\prime},\vec{k},t)-1\right)$
is the hole density. The electron charge $-e$ is henceforth set to 1. In
total, the Hartree self-energy produces the classical Coulomb potential of all
electrons in the system exerted on an electron located at a position $\vec{x}$
according to

$\displaystyle\sigma^{H}_{ab;\lambda\lambda^{\prime}}(\vec{x},t,\vec{p},\omega)=\delta_{\lambda\lambda^{\prime}}\delta_{ab}\int
d\vec{x}^{\prime}V(\vec{x}-\vec{x}^{\prime})n(\vec{x}^{\prime},t).$ (71)

The Hartree term is real-valued and independent of the momentum and frequency
variables. Furthermore, it is diagonal in Keldysh space, meaning its Keldysh
component is zero.

(a)

(b)

Figure 4: The self-energies at first order in the interaction. The left
diagram is known as the direct or Hartree contribution and the right diagram
is known as the exchange or Fock contribution.

####  IV.2.2 The Fock diagram

Next, we are going to sketch the calculation of the Fock self-energy diagram
depicted in Fig. 4(b). It reads

$\displaystyle-i\Sigma^{F}_{ab;\lambda\lambda^{\prime}}(\vec{x},\vec{x}^{\prime};t,t^{\prime})$
$\displaystyle\hskip
42.67912pt=\gamma^{\alpha}_{ac}D_{0,\alpha\beta}(\vec{x},\vec{x}^{\prime};t,t^{\prime})G_{0,cd;\lambda\lambda^{\prime}}(\vec{x},\vec{x}^{\prime};t,t^{\prime})\gamma^{\beta}_{db}.$

After a Wigner transformation followed by a few steps of algebraic
manipulations, we find that

$\displaystyle-i\Sigma^{F}_{ab;\lambda\lambda^{\prime}}(\vec{x},\vec{p};t,\nu)$
$\displaystyle=\delta_{ab}\int\frac{d\vec{p}_{1}}{(2\pi)^{2}}\frac{V(\vec{p}-\vec{p}_{1})}{2}\int\frac{d\nu_{1}}{2\pi}G^{K}_{0,\lambda\lambda^{\prime}}(\vec{x},\vec{p}_{1};t,\nu_{1}).$

Subsequently, we transform it into the quasi-particle basis according to
$-i\sigma^{F}_{ab;\lambda\lambda^{\prime}}(\vec{x},\vec{p};t,\nu)=-i\left(\mathcal{U}^{\dagger}_{\vec{p}}\Sigma^{F}_{ab}(\vec{x},\vec{p};t,\nu)\mathcal{U}_{\vec{p}}\right)_{\lambda\lambda^{\prime}}$.
We are interested in elements on the main diagonal of the spinor space. They
are given by

$\sigma^{F}_{ab;\lambda\lambda}(\vec{x},\vec{p};t,\nu)=\sigma^{F,1}_{ab;\lambda\lambda}(\vec{x},\vec{p};t,\nu)+\sigma^{F,2}_{ab;\lambda\lambda}(\vec{x},\vec{p};t,\nu),$
(74)

where

$\sigma^{F,1}_{ab;\lambda\lambda}(\vec{x},\vec{p};t,\nu)=\delta_{ab}\lambda\int\frac{d\vec{p}_{1}}{(2\pi)^{2}}\cos(\theta_{\vec{p}_{1}}-\theta_{\vec{p}})\frac{V(\vec{p}-\vec{p}_{1})}{2},$
(75)

and

$\displaystyle\sigma^{F,2}_{ab;\lambda\lambda}(\vec{x},\vec{p};t,\nu)$
$\displaystyle\hskip
42.67912pt=\delta_{ab}\int\frac{d\vec{p}_{1}}{(2\pi)^{2}}\frac{1+\cos(\theta_{\vec{p}_{1}}-\theta_{\vec{p}})}{2}\frac{V(\vec{p}-\vec{p}_{1})}{2}$
$\displaystyle\hskip
128.0374pt\left((1-\lambda)-2f_{\lambda}(\vec{x},t,\vec{p}_{1})\right),$ (76)

where $\theta_{\vec{p}}$ was introduced in Eq. (63). The first term diverges
logarithmically and is responsible for the renormalization of the Fermi
velocity. To see this, we explicitly substitute the Coulomb potential in Eq.
(9) followed by a transformation of the momentum variable from Cartesian to
polar coordinates. This gives

$\displaystyle\sigma^{F,1}_{ab;\lambda\lambda}(\vec{x},\vec{p};t,\nu)$
$\displaystyle=$ $\displaystyle\lambda\frac{\alpha v_{F}}{4\pi}\int
p_{1}dp_{1}d\theta\frac{\cos\theta}{\sqrt{p^{2}+p_{1}^{2}-2pp_{1}\cos\theta}}$
$\displaystyle\approx$ $\displaystyle\lambda\frac{\alpha
v_{F}p}{4}\int_{p}^{\mathcal{K}}dp_{1}\frac{1}{p_{1}}=\frac{\alpha}{4}\ln(\mathcal{K}/p)\lambda
v_{F}p.$

$\theta=\theta_{\vec{p}_{1}}-\theta_{\vec{p}}$ is the angle of the momentum
$\vec{p}_{1}$ measured with respect to $\vec{p}$. In order to extract the
logarithmic divergence of this integral, we expanded the square root to the
first power in $p/p_{1}$, i.e.,
$1/\sqrt{p^{2}+p_{1}^{2}-2pp_{1}\cos\theta}\approx
1/p_{1}\left(1+p\cos\theta/p_{1}\right)$. After the angular integral, we find
the result in the second line of Eq. (LABEL:eq:velocityrenormalization). The
divergence of the integral is cut off at the inverse lattice spacing
$\mathcal{K}$. The lower boundary of integration is consistently set to $p$
(we require $p/p_{1}\ll 1$). When substituting Eq.
(LABEL:eq:velocityrenormalization) into the Dyson equation, we find

$g^{-1}_{0}-\sigma^{F}=\omega-\lambda\left(1+\alpha\ln(\mathcal{K}/p)/4\right)v_{F}p.$
(78)

the Fermi velocity is renormalized accordingly to
$v^{R}_{F}=\left(1+\alpha\ln(\mathcal{K}/p)/4\right)v_{F}$. Such logarithmic
renormalization was first discussed within the renormalization group approach
in Ref. vozmediano and measured in graphene in Ref. Elias2011. The second
contribution, Eq. (76), describes the conventional exchange energy coming from
both electrons and holes. It varies with the doping $\mu$ and the temperature
$T$ of the system. Moreover, different from the exchange energy of a
conventional two-dimensional electron gas, there is a factor
$(1\pm\cos(\theta_{\vec{p}_{1}}-\theta_{\vec{p}}))/2$ arising from the
wavefunction overlap, where $\theta_{\vec{p}_{1}}-\theta_{\vec{p}}$ is the
angle between $\vec{p}$ and $\vec{p}_{1}$ Hwang et al. (2007). In total, it
reads

$\displaystyle(\sigma^{F,2})^{R}_{\pm\pm}(\vec{p})$ $\displaystyle\hskip
14.22636pt=-\int\frac{d\vec{p}_{1}}{(2\pi)^{2}}\frac{1\pm\cos(\theta_{\vec{p}_{1}}-\theta_{\vec{p}})}{2}V(\vec{p}-\vec{p}_{1})f^{0}_{+}(\vec{p}_{1})$
$\displaystyle\hskip
22.76228pt+\int\frac{d\vec{p}_{1}}{(2\pi)^{2}}\frac{1\mp\cos(\theta_{\vec{p}_{1}}-\theta_{\vec{p}})}{2}V(\vec{p}-\vec{p}_{1})(1-f^{0}_{-}(\vec{p}_{1})).$
(79)

After a series of manipulations, we find

$\displaystyle(\sigma^{F,2})^{R}_{\pm\pm}(\vec{p})$
$\displaystyle=-\frac{\alpha v_{F}}{4\pi}\int
dp_{1}\sqrt{\frac{p_{1}}{2p}}\left(g\left(\frac{p^{2}+p_{1}^{2}}{2pp_{1}}\right)\pm
h\left(\frac{p^{2}+p_{1}^{2}}{2pp_{1}}\right)\right)$ $\displaystyle\hskip
170.71652ptf^{0}_{+}(\vec{p}_{1})$ $\displaystyle+\frac{\alpha
v_{F}}{4\pi}\int
dp_{1}\sqrt{\frac{p_{1}}{2p}}\left(g\left(\frac{p^{2}+p_{1}^{2}}{2pp_{1}}\right)\mp
h\left(\frac{p^{2}+p_{1}^{2}}{2pp_{1}}\right)\right)$ $\displaystyle\hskip
170.71652pt(1-f^{0}_{-}(\vec{p}_{1})),$ (80)

with

$\displaystyle g(x)$ $\displaystyle=$
$\displaystyle\int_{0}^{2\pi}\frac{d\theta}{\sqrt{x-\cos\theta}}\;,$
$\displaystyle h(x)$ $\displaystyle=$
$\displaystyle\int_{0}^{2\pi}\frac{d\theta\;\cos\theta}{\sqrt{x-\cos\theta}}\;.$
(81)

We proceed by evaluating the $p_{1}$-integral numerically. We show in Fig. 5
the exchange energies for electrons, $(\sigma^{F,2})^{R}_{++}$, and holes,
$(\sigma^{F,2})^{R}_{--}$, at the chemical potential $\mu/T=1$ for various
values of momenta. We see that the self-energy becomes less important at high
momenta. In Fig. 6 we show the spectral function of the electrons with and
without the Fock correction. Since the exchange energy is a real-valued
function, there remains the delta-peak feature manifested as the blue line in
the middle of the plot. Compare to the spectral function of the non-
interacting theory in Fig. 6 (a), the exchange energy plays two roles: (i) it
shifts the chemical potential and (ii) it increases the Fermi velocity Hwang
et al. (2007). However, the effect of $(\sigma^{F,2})^{R}_{\pm\pm}(\vec{p})$
is relatively small compared to the Fermi velocity renormalization coming from
the logarithmic divergence. Therefore, we keep only the latter effect in the
following. In practice, this leads to replacing $v_{F}$ with $v_{F}^{R}$, the
renormalized Fermi velocity, in all expressions.

Figure 5: Exchange self-energies for Dirac fermions with the chemical
potential $\mu/T=1$ as functions of wave vector $q/T$. The fine structure
constant $\alpha$ is chosen to be $1$. The blue curve represents
$(\sigma^{F,2})^{R}_{++}$ showing the exchange energy correction of the
electrons. The orange curve shows $(\sigma^{F,2})^{R}_{--}$ the exchange
energy of the holes.

Figure 6: Spectral function at $\mu/T=1$ and $\alpha=1$. (a) The spectral
function of the non-interacting electrons follows $\Im
g_{0}^{R}(\vec{q},\omega)=-2\pi i\delta(\omega-v_{F}q+\mu)$. The blue line in
the middle of the plot manifests the delta peak at $\omega=v_{F}q-\mu$
($v_{F}$ is set to 1 in the plot). (b) The spectral function with the
inclusion of the exchange conventional energy in Eq. (79). Since the exchange
energy is a real-valued function, the spectral function still has the delta-
peak feature as shown by the blue line in the middle of the plot.

####  IV.2.3 Energy spectrum

We are now ready to evaluate the energy spectrum of quasi-electrons and quasi-
holes in equilibrium within the Hartree-Fock approximation. The dispersion
relation of the excitations can be found from the roots of the inverse Green
function. In the presence of the interaction, the inverse Green function is
the solution of the Dyson equation in Eq. (32). We first transform the Dyson
equation into the quasi-particle basis, which leads to

$(g^{-1})^{R}_{\lambda\lambda^{\prime}}=(g^{-1}_{0})^{R}_{\lambda\lambda^{\prime}}-(\sigma^{H})^{R}_{\lambda\lambda^{\prime}}.$
(82)

where the retarded component of the Hartree self-energy is given by the first
row and first column element of Eq. (71) that is
$(\sigma^{H})_{\lambda\lambda^{\prime}}^{R}=\sigma^{H}_{11,\lambda\lambda^{\prime}}$.
This gives

$\displaystyle(g^{-1})^{R}_{\lambda\lambda^{\prime}}(\vec{x},t,\vec{p},\omega)$
$\displaystyle=$ $\displaystyle\Big{(}\omega+i\delta\mp\lambda v_{F}^{R}p+\mu-
V_{\rm{ex}}(\vec{x})$ $\displaystyle-$ $\displaystyle\int
d\vec{x}^{\prime}V(\vec{x}-\vec{x}^{\prime})n(\vec{x}^{\prime},t)\Big{)}\delta_{\lambda\lambda^{\prime}}\;,$

where $n(\vec{x},t)$ is the total electron density. By inserting
$V_{\rm{ex}}(\vec{x})$ from Eq. (3), we find that

$(g^{-1})^{R/A}_{\lambda\lambda^{\prime}}(\vec{x},t,\vec{p},\omega)=\omega\pm
i\delta-\lambda v_{F}^{R}p+\mu-V^{H}(\vec{x},t),$ (84)

where the Hartree potential $V^{H}(\vec{x},t)$ is defined as

$\displaystyle V^{H}(\vec{x},t)=\int
d\vec{x}^{\prime}V(\vec{x}-\vec{x}^{\prime})\left(n(\vec{x}^{\prime},t)-n_{0}\right).$
(85)

The Hartree term represents the potential energy of an electron at a position
$\vec{x}$. This potential is produced by all other electrons with the density
$n(\vec{x}^{\prime},t)$ at the other positions $\vec{x}^{\prime}$ through the
Coulomb interaction. This contribution is partially canceled by the potential
$V_{\rm{ex}}(\vec{x})$ arising from the interaction between electrons and the
underlying jellium background. This fixed uniformly distributed positively
charged background guarantees the over-all electrical neutrality of the
system. Here $n_{0}$ is the ion charge density which is identical to the
electron charge density in thermal equilibrium. Hence, in global thermal
equilibrium, the Hartree potential vanishes and thus the energy spectrum of
electrons is given by

$\epsilon_{\lambda}(\vec{p})=\lambda v_{F}^{R}p\;,$ (86)

as expected.

####  IV.2.4 The kinetic equation for Dirac fermions

Next, we derive the Bolzmann equation within the Hartree-Fock approximation.
To this end, we consider Eq. (56). Using the Hartree-Fock approximation in Eq.
(84), we find

$\displaystyle\partial_{\epsilon}\left(g^{-1}_{0}-(\Re\sigma^{H})^{R}\right)$
$\displaystyle=$ $\displaystyle 1,$
$\displaystyle\partial_{\vec{p}}\left(g^{-1}_{0}-(\Re\sigma^{H})^{R}\right)$
$\displaystyle=$ $\displaystyle-\lambda v_{F}^{R}\hat{p},$
$\displaystyle\partial_{\vec{x}}\left(g^{-1}_{0}-(\Re\sigma^{H})^{R}\right)$
$\displaystyle=$ $\displaystyle-\partial_{\vec{x}}V^{H}(\vec{x},t).$ (87)

We first substitute these derivatives into Eq. (56) followed by a
multiplication of the resulting equation with the spectral function $\Im
g^{-1}_{\lambda\lambda}(\vec{p},\epsilon)=-2\pi i\delta(\epsilon-\lambda
v_{F}^{R}p)$. Subsequently, we integrate it over the frequency which amounts
to the quasi-particle approximation. In the end, we obtain the mean-field
collisionless Boltzmann equation, also known as the Vlasov equation
Vlasov1968, for electrons ($\lambda=+$) and holes ($\lambda=-$)

$\displaystyle\partial_{t}f_{\lambda}(\vec{x},t,\vec{p})+\lambda
v_{F}^{R}\hat{p}\cdot\partial_{\vec{x}}f_{\lambda}(\vec{x},t,\vec{p})$
$\displaystyle\hskip
71.13188pt-\;\partial_{\vec{x}}V^{H}(\vec{x},t)\cdot\partial_{\vec{p}}f_{\lambda}(\vec{x},t,\vec{p})=0,$
(88)

where $\hat{p}$ denotes the unit vector in the direction of the momentum
$\vec{p}$. In the above equation, the equilibrium value of the distribution
functions is given by $f_{\lambda}^{0}(\vec{p})=(1+\exp(\frac{\lambda
v_{F}^{R}p-\mu}{T}))^{-1}$. The potential $V^{H}(\vec{x},t)=\int
d\vec{x}^{\prime}V(\vec{x}-\vec{x}^{\prime})\delta n(\vec{x}^{\prime},t)$
results from the Hartree self-energy, where $\delta n=n-n_{0}$ is the density
fluctuation. This potential is also the solution of the classical Poisson
equation for the internal electric field.

To summarize, we found that the Hartree-Fock diagrams lead to the Vlasov
equation.

###  IV.3 Second-order perturbation theory: Born approximation

An important role in hydrodynamic systems is played by the relaxation
processes towards local equilibrium that conserve particle number, momentum,
and energy. These collisions occur beyond first order in the interaction, Eq.
(II). The lowest non-vanishing order is second order and the contributions are
pictorially shown in Fig. 7. This is called the Born approximation for the
cross section Fritz et al. (2008); Kashuba (2008); Landau (1936); Kadanoff and
Baym (1962). In principle, these diagrams play two roles: (i) they describe
the aforementioned relaxations due to collisions and (ii) they renormalize the
quasi-particle properties Kotov2012; Mishchenko.

The calculation of these diagrams is tedious but straightforward. Here we
summarize our final results and present the full derivation in Appendix A. The
Boltzmann equation reads

$\displaystyle\partial_{t}f_{\lambda}(\vec{x},t,\vec{p})+\lambda
v_{F}^{R}\hat{p}\cdot\partial_{\vec{x}}f_{\lambda}(\vec{x},t,\vec{p})-\;\partial_{\vec{x}}V^{H}(\vec{x},t)\cdot\partial_{\vec{p}}f_{\lambda}(\vec{x},t,\vec{p})$
$\displaystyle\hskip
85.35826pt=-\int\frac{d\vec{k}_{1}}{(2\pi)^{2}}\frac{d\vec{q}}{(2\pi)^{2}}2\pi\delta(\lambda\epsilon_{\vec{k}}-\lambda_{1}\epsilon_{\vec{k}-\vec{q}}-\lambda_{2}\epsilon_{\vec{k}_{1}+\vec{q}}+\lambda_{3}\epsilon_{\vec{k}_{1}})R_{\lambda\lambda_{1}\lambda_{3}\lambda_{2}}(\vec{k},\vec{k}_{1},\vec{q})$
$\displaystyle\hskip
85.35826pt\Big{[}f_{\lambda}(\vec{k})f_{\lambda_{3}}(\vec{k}_{1})(1-f_{\lambda_{1}}(\vec{k}-\vec{q}))(1-f_{\lambda_{2}}(\vec{k}_{1}+\vec{q}))-(1-f_{\lambda}(\vec{k}))(1-f_{\lambda_{3}}(\vec{k}_{1}))f_{\lambda_{1}}(\vec{k}-\vec{q})f_{\lambda_{2}}(\vec{k}_{1}+\vec{q})\Big{]},$
(89)

where we introduced the shorthand notation

$\displaystyle
R_{\lambda\lambda_{1}\lambda_{3}\lambda_{2}}(\vec{k},\vec{k}_{1},\vec{q})=2\Big{[}|T_{\lambda\lambda_{1}\lambda_{3}\lambda_{2}}-T_{\lambda\lambda_{2}\lambda_{1}\lambda_{3}}|^{2}+(N-1)\left(|T_{\lambda\lambda_{1}\lambda_{3}\lambda_{2}}|^{2}+|T_{\lambda\lambda_{2}\lambda_{1}\lambda_{3}}|^{2}\right)\Big{]}.$
(90)

In this expression, we have used

$\displaystyle
T_{\lambda\lambda_{1}\lambda_{2}\lambda_{3}}(\vec{k},\vec{k}_{1},\vec{q})=\frac{V(\vec{q})}{2}M^{\lambda\lambda_{1}}_{\vec{k},\vec{k}-\vec{q}}M^{\lambda_{2}\lambda_{3}}_{\vec{k}_{1},\vec{k}_{1}+\vec{q}}\;,$
(91)

where the coherence factor $M$ comes from the overlap of the single-particle
wavefunctions. It is defined according to

$\displaystyle
M^{\lambda\lambda_{1}}_{\vec{k},\vec{k}_{1}}=\left(\mathcal{U}^{\dagger}_{\vec{k}}\mathcal{U}^{\phantom{\dagger}}_{\vec{k}_{1}}\right)_{\lambda\lambda_{1}}\;.$
(92)

For brevity, we suppress the space and time variables of the distribution
functions in the collision terms and have in mind that they all depend on the
same set of variables that is $(\vec{x},t)$. The collision does not shift the
center-of-mass and time coordinates. This effect indeed exists, but it will
show up at higher order in the gradient expansion pesin. We can understand
this collision integral in the following way: an electron in band $\lambda$
with momentum $\vec{k}$ is scattered into band $\lambda_{1}$ and momentum
$\vec{k}+\vec{q}$ by a collision with another electron in band $\lambda_{3}$
and state $\vec{k}_{1}$ which is itself scattered into the energy band
$\lambda_{2}$ and state $\vec{k}_{1}-\vec{q}$. For this event to take place,
the initial states $\vec{k}$ and $\vec{k}_{1}$ have to be filled and the final
states $\vec{k}-\vec{q}$ and $\vec{k}_{1}+\vec{q}$ must be empty. Thus, the
factors $f_{\lambda}(\vec{k})$ and $f_{\lambda_{3}}(\vec{k}_{1})$ are the
occupation numbers of these state and $1-f_{\lambda_{1}}(\vec{k}-\vec{q})$ and
$1-f_{\lambda_{2}}(\vec{k}_{1}+\vec{q})$ are the probabilities for the final
states to be unoccupied. The conservation of energy is taken into account by
the delta function. The transition probability of this event is
$R_{\lambda\lambda_{1}\lambda_{3}\lambda_{2}}$. In total, we find that Coulomb
interaction enters the Boltzmann transport equation in two ways: (i) as the
Hartree potential produced by all the other particles in the system and (ii)
inelastic and momentum-conserving electron-electron scatterings leading to
local equilibration. Let us note that within this approximation, the
contribution from the real part of the second-order diagrams to the left-hand
side of Eq. (89) is neglected.

Finally, let us note that the Born approximation is valid only when the
kinetic energy of the electrons is large compared to the Coulomb interaction
potential Ziman (2000). For the Dirac system, the ratio of the potential
energy to the kinetic energy is characterized by the fine structure constant
$\alpha=e^{2}/4\pi\epsilon v_{F}$. In condensed-matter systems, this constant
is not necessarily small (for graphene $\alpha\approx 0.3-2.2$). In such a
strong interaction limit, a perturbative series expansion in $\alpha$ may
break down. Instead one can employ an alternative perturbative expansion in
the other parameters. In the subsequent section, we will employ the random-
phase approximation (RPA) and show that it gives a different, more complicated
picture than the Hartree-Fock-Born result presented in this section.

(a)

(b)

(c)

Figure 7: Self-energy to second order in the interaction.

###  IV.4 Conservation laws

The collision integral on the right-hand side of the Boltzmann equation in Eq.
(89), henceforth denoted $C_{\lambda}[f](\vec{k})$, has conservation laws
encoded in it. These conservation laws are important for two reasons: (i) They
allow for an identification or derivation of physical quantities such as
charge and current densities; (ii) it provides the basis of the derivation of
conservation laws and even transport phenomena. When a system is driven away
from equilibrium, the first thing that happens is that collisions drive the
system to local equilibrium. Afterwards, there is a much slower return to
global equilibrium. The latter describes transport processes and is governed
by the conservation laws. The conserved quantities in the system under
consideration are particle number, momentum, and energy. The conservation laws
of electric charge, momentum, and energy are obtained by multiplying the
Boltzmann equation in Eq. (89), by $1$, $\vec{p}$ and
$\epsilon_{\lambda}(\vec{x},\vec{p})=\lambda v_{F}^{R}p+V^{H}(\vec{x})$ and
then integrating the resulting equations over all momentum $\vec{p}$ as well
as summing over energy bands $\pm$ and flavors. This leads to the following
collisional invariants

$N\sum_{\lambda=\pm}\int\frac{d\vec{p}}{(2\pi)^{2}}\;C_{\lambda}[f](\vec{p})=0,$
(93)
$N\sum_{\lambda=\pm}\int\frac{d\vec{p}}{(2\pi)^{2}}\;\vec{p}C_{\lambda}[f](\vec{p})=0,$
(94)

and

$N\sum_{\lambda=\pm}\int\frac{d\vec{p}}{(2\pi)^{2}}\;\epsilon_{\lambda}(\vec{x},\vec{p})C_{\lambda}[f](\vec{p})=0.$
(95)

(i) The continuity equation of charge can be obtained by integrating the
Boltzmann equation in Eq. (89), over momentum $\vec{p}$ followed by a
summation over the band index and flavor index, see Eq. (93). In contrast to
the case of a one band system, there is a subtle point here relating to the
infinite number of particles in the filled band which is unbounded from below.
This infinite constant vanishes upon differentiation and does not contribute
to the continuity equation. Therefore, we can subtract the infinite
contribution coming from the Dirac sea and instead consider the population of
holes defined as $f(\vec{x},t,\vec{k})-1$. First, we consider the time-
derivative term. Integrating this term over all states gives

$\displaystyle
N\sum_{\lambda=\pm}\int\frac{d\vec{p}}{(2\pi)^{2}}\partial_{t}f_{\lambda}(\vec{x},t,\vec{p})$
$\displaystyle\hskip
5.69046pt=\partial_{t}\;N\int\frac{d\vec{p}}{(2\pi)^{2}}\Big{[}f_{+}(\vec{x},t,\vec{p})+(f_{-}(\vec{x},t,\vec{p})-1)\Big{]}.$
(96)

We denote the charge density by

$n(\vec{x},t)=N\int\frac{d\vec{p}}{(2\pi)^{2}}\;\left[f_{+}(\vec{x},t,\vec{p})+(f_{-}(\vec{x},t,\vec{p})-1)\right].$
(97)

The second term can be integrated in a similar fashion and gives

$\displaystyle N\sum_{\lambda=\pm}\int\frac{d\vec{p}}{(2\pi)^{2}}\lambda
v_{F}^{R}\hat{p}\cdot\partial_{\vec{x}}f_{\lambda}(\vec{x},t,\vec{p})$
$\displaystyle=\partial_{\vec{x}}\cdot\;N\int\frac{d\vec{p}}{(2\pi)^{2}}v_{F}^{R}\hat{p}\Big{[}f_{+}(\vec{x},t,\vec{p})+(f_{-}(\vec{x},t,\vec{p})-1)\Big{]}.$

This term can be identified with the charge current density

$\displaystyle\vec{j}(\vec{x},t)$ $\displaystyle=$ $\displaystyle
N\int\frac{d\vec{p}}{(2\pi)^{2}}\;v_{F}^{R}\hat{p}\Big{[}f_{+}(\vec{x},t,\vec{p})-(f_{-}(\vec{x},t,\vec{p})-1)\Big{]}.$
, (99)

The momentum-derivative term vanishes upon integration

$\displaystyle-N\sum_{\lambda=\pm}\int\frac{d\vec{p}}{(2\pi)^{2}}\partial_{\vec{x}}V^{H}(\vec{x},t)\cdot\partial_{\vec{p}}f_{\lambda}(\vec{x},t,\vec{p})$
$\displaystyle\hskip
5.69046pt=-N\sum_{\lambda=\pm}\int\frac{d\vec{p}}{(2\pi)^{2}}\partial_{\vec{p}}\cdot\Big{[}\partial_{\vec{x}}V^{H}(\vec{x},t)f_{\lambda}(\vec{x},t,\vec{p})\Big{]}=0\;,$

since it is a total derivative. Combining the above equations together with
Eq. (93), we find the continuity equation for the electric charge according to

$\partial_{t}n(\vec{x},t)+\partial_{\vec{x}}\cdot\vec{j}(\vec{x},t)=0.$ (101)

It is worthwhile noting here that the particle density is conserved locally.

(ii) Momentum conservation can be obtained by multiplying the Boltzmann
equation in Eq. (89) by momentum before integrating the resulting equation
over all states, see Eq. (94). The time-derivative term yields

$\displaystyle
N\sum_{\lambda=\pm}\int\frac{d\vec{p}}{(2\pi)^{2}}\vec{p}\partial_{t}f_{\lambda}(\vec{x},t,\vec{p})$
$\displaystyle\hskip
2.84544pt=\partial_{t}\;N\int\frac{d\vec{p}}{(2\pi)^{2}}\vec{p}\Big{[}f_{+}(\vec{x},t,\vec{p})+(f_{-}(\vec{x},t,\vec{p})-1)\Big{]}.$

This allows to define the momentum density,

$\vec{n}^{\vec{p}}(\vec{x},t)=N\int\frac{d\vec{p}}{(2\pi)^{2}}\vec{p}\Big{[}f_{+}(\vec{x},t,\vec{p})+(f_{-}(\vec{x},t,\vec{p})-1)\Big{]}.$
(103)

The space-derivative term can be similarly integrated and gives

$\displaystyle
N\sum_{\lambda=\pm}\int\frac{d\vec{p}}{(2\pi)^{2}}\vec{p}\;\lambda
v_{F}^{R}\hat{p}\cdot\partial_{\vec{x}}f_{\lambda}(\vec{x},t,\vec{p})$
$\displaystyle=\partial_{\vec{x}}\cdot\;N\int\frac{d\vec{p}}{(2\pi)^{2}}v_{F}^{R}\vec{p}\hat{p}\Big{[}f_{+}(\vec{x},t,\vec{p})-(f_{-}(\vec{x},t,\vec{p})-1)\Big{]}.$

From this result, we define the momentum-flux tensor,

$\vec{\vec{\Pi}}(\vec{x},t)=N\int\frac{d\vec{p}}{(2\pi)^{2}}\;v_{F}^{R}\vec{p}\hat{p}\Big{[}f_{+}(\vec{x},t,\vec{p})-(f_{-}(\vec{x},t,\vec{p})-1)\Big{]}.$
(105)

The momentum-derivative term yields

$\displaystyle-N\sum_{\lambda=\pm}\int\frac{d\vec{p}}{(2\pi)^{2}}\vec{p}\partial_{\vec{x}}V^{H}(\vec{x},t)\cdot\partial_{\vec{p}}f_{\lambda}(\vec{x},t,\vec{p})$
$\displaystyle\hskip
5.69046pt=\partial_{\vec{x}}V^{H}(\vec{x},t)N\int\frac{d\vec{p}}{(2\pi)^{2}}\;\left[f_{+}(\vec{x},t,\vec{p})+(f_{-}(\vec{x},t,\vec{p})-1)\right]$
$\displaystyle\hskip
5.69046pt=\partial_{\vec{x}}V^{H}(\vec{x},t)n(\vec{x},t).$

which defines a force term. Finally, combining the above equations, we find
the momentum equation according to

$\partial_{t}\vec{n}^{\vec{p}}(\vec{x},t)+\partial_{\vec{x}}\cdot\vec{\vec{\Pi}}(\vec{x},t)=-\partial_{\vec{x}}V^{H}(\vec{x},t)n(\vec{x},t).$
(107)

This has a straightforward interpretation. The momentum of the electron fluid
in any given volume element can be changed in two ways: (i) by means of
momentum flow through the volume boundary accounting for the space-gradient
term and (ii) by the internal electric field. Locally, the momentum density is
not conserved and changed by the internal electric force
$-\partial_{\vec{x}}V^{H}(\vec{x},t)$. This internal force, however, does not
affect the total momentum of the entire system. We still expect that the total
momentum of the system is conserved. This can be shown explicitly by
integrating the momentum equation over all space. The integration of the
momentum-flux gradient results in a surface term which vanishes. The Hartree
force term also vanishes. To see this, we consider the property of the Coulomb
potential

$\partial_{\vec{x}}\frac{1}{|\vec{x}-\vec{x}^{\prime}|}=-\frac{\vec{x}-\vec{x}^{\prime}}{|\vec{x}-\vec{x}^{\prime}|^{3}}=-\partial_{\vec{x}^{\prime}}\frac{1}{|\vec{x}-\vec{x}^{\prime}|}.$
(108)

Based on this, one can show explicitly that

$\displaystyle\int d\vec{x}\partial_{\vec{x}}V^{H}(\vec{x},t)n(\vec{x},t)$
$\displaystyle=\int
d\vec{x}d\vec{x}^{\prime}\partial_{\vec{x}}V(\vec{x}-\vec{x}^{\prime})(n(\vec{x}^{\prime},t)-n_{0})n(\vec{x},t)=0.$

Combining all the above equations, we obtain

$\partial_{t}\vec{P}=0,$ (110)

where the total momemtum of the entire system is given by

$\vec{P}(t)=\int d\vec{x}\;\vec{n}^{\vec{p}}(\vec{x},t).$ (111)

(iii) Similarly, energy conservation is obtained by multiplying the Boltzmann
equation, Eq. (89) by energy $\epsilon(\vec{x},\vec{p})$ followed by
integrating and summing the equation over all states, see Eq. (95). The time-
derivative term yields

$\displaystyle
N\sum_{\lambda=\pm}\int\frac{d\vec{p}}{(2\pi)^{2}}\epsilon_{\lambda}(\vec{x},\vec{p})\partial_{t}f_{\lambda}(\vec{x},t,\vec{p})$
$\displaystyle\hskip
28.45274pt=\partial_{t}\;N\int\frac{d\vec{p}}{(2\pi)^{2}}\Big{[}\epsilon_{+}(\vec{x},\vec{p})f_{+}(\vec{x},t,\vec{p})$
$\displaystyle\hskip
28.45274pt+\epsilon_{-}(\vec{x},\vec{p})(f_{-}(\vec{x},t,\vec{p})-1)\Big{]}.$
(112)

This defines the energy density according to

$\displaystyle n^{\epsilon}(\vec{x},t)$ $\displaystyle=$ $\displaystyle
N\int\frac{d\vec{p}}{(2\pi)^{2}}\;\Big{[}\epsilon_{+}(\vec{x},\vec{p})f_{+}(\vec{x},t,\vec{p})$
(113) $\displaystyle\hskip
28.45274pt+\;\epsilon_{-}(\vec{x},\vec{p})(f_{-}(\vec{x},t,\vec{p})-1)\Big{]}.$

The spatial derivative leads to

$\displaystyle
N\sum_{\lambda=\pm}\int\frac{d\vec{p}}{(2\pi)^{2}}\epsilon_{\lambda}(\vec{x},\vec{p})\;\lambda
v_{F}^{R}\hat{p}\cdot\partial_{\vec{x}}f_{\lambda}(\vec{x},t,\vec{p})$
$\displaystyle=\partial_{\vec{x}}\cdot\;N\int\frac{d\vec{p}}{(2\pi)^{2}}v_{F}^{R}\hat{p}\Big{[}\epsilon_{+}(\vec{x},\vec{p})f_{+}(\vec{x},t,\vec{p})$
$\displaystyle\hskip
14.22636pt-\;\epsilon_{-}(\vec{x},\vec{p})(f_{-}(\vec{x},t,\vec{p})-1)\Big{]}-\partial_{\vec{x}}V^{H}(\vec{x},t)\cdot\vec{j}.$
(114)

We find that it consists of two terms. The first term describes a divergence
of an energy current density defined as

$\displaystyle\vec{j}^{\epsilon}(\vec{x},t)$ $\displaystyle=$ $\displaystyle
N\int\frac{d\vec{p}}{(2\pi)^{2}}\;v_{F}^{R}\hat{p}\Big{[}\epsilon_{+}(\vec{x},\vec{p})f_{+}(\vec{x},t,\vec{p})$
(115)
$\displaystyle-\;\epsilon_{-}(\vec{x},\vec{p})(f_{-}(\vec{x},t,\vec{p})-1)\Big{]},$

whereas the second term describes Joule heating due to the internal force.
Next, we consider the momentum-derivative term and find that it can also be
rewritten as Joule heating due to the internal Coulomb force.

$\displaystyle-N\sum_{\lambda=\pm}\int\frac{d\vec{p}}{(2\pi)^{2}}\epsilon_{\lambda}(\vec{x},\vec{p})\;\partial_{\vec{x}}V^{H}(\vec{x},t)\cdot\partial_{\vec{p}}f_{\lambda}(\vec{x},t,\vec{p})$
$\displaystyle\hskip
113.81102pt=\partial_{\vec{x}}V^{H}(\vec{x},t)\cdot\vec{j}\;.$ (116)

This term will thus be canceled by the second term on the right-hand side of
Eq. (114). Combining these equations, we find the continuity equation of
energy:

$\partial_{t}n^{\epsilon}(\vec{x},t)+\partial_{\vec{x}}\cdot\vec{j}^{\epsilon}(\vec{x},t)=0.$
(117)

Eqs. (101), (107), and (117) constitute the main equations of electron
hydrodynamics and can be used to derive the Navier-Stokes equations. Let us
emphasize once again that, in contrast to the hydrodynamic equations of usual
fluids, electrons interact among themselves via long-range Coulomb
interactions and this effect shows up on the right-side of the momentum
equation in Eq. (107).

###  IV.5 Collective modes

Now that we have the fundamental equations of hydrodynamics of weakly
interacting charged Dirac electrons, let us study some of its properties. As
will be presented in Appendix B, there are three independent hydrodynamic
variables. Here we choose the set of independent variables consisting of the
charge density ($n$), the energy density ($n^{\epsilon}$) and the hydrodynamic
velocity $\vec{u}$ in terms of which the other quantities can be written. To
linear order in $\vec{u}$, we find that the charge current density is given by
$\vec{j}=n\vec{u}$. The momentum flux is associated with the pressure by means
of $\Pi_{ij}=P\delta_{ij}$ where the pressure is in turn proportional to the
energy density according to $P=n^{\epsilon}/2$. One of the consequences of the
linear spectrum of Dirac electrons is that the momentum density is decoupled
from the charge current. Instead, it is proportional to the energy currents
according to $n^{\vec{p}}=\vec{j}^{\epsilon}/v_{F}^{2}$ where the energy
current is given by $j^{\epsilon}=(P+n^{\epsilon})\vec{u}$ . We now consider
an electron fluid at rest with constant $n=n_{0}$,
$n^{\epsilon}=n_{0}^{\epsilon}$, and $\vec{u}=0$. We are interested in small
fluctuations around the constant and homogeneous solution, and put
$n=n_{0}+\delta n$ $n^{\epsilon}=n_{0}^{\epsilon}+\delta n^{\epsilon}$ and
assume small $\vec{u}$. We insert this solution into the hydrodynamic
equations in Eqs. (101, (107), and (117). By keeping terms up to linear order
in the fluctuations, we obtain

$\displaystyle\partial_{t}\delta
n(\vec{x},t)+n_{0}\partial_{\vec{x}}\cdot\vec{u}(\vec{x},t)$ $\displaystyle=$
$\displaystyle 0,$
$\displaystyle(P_{0}+n^{\epsilon}_{0})\partial_{t}\vec{u}(\vec{x},t)/v_{F}^{2}+\partial_{\vec{x}}\delta
P(\vec{x},t)$ $\displaystyle=$
$\displaystyle-\partial_{\vec{x}}V^{H}(\vec{x},t)n_{0},$
$\displaystyle\partial_{t}\delta
n^{\epsilon}(\vec{x},t)+(P_{0}+n_{0}^{\epsilon})\partial_{\vec{x}}\cdot\vec{u}(\vec{x},t)$
$\displaystyle=$ $\displaystyle 0.$ (118)

The solutions to these equations are propagating waves. As an ansatz, we
insert

$\displaystyle\begin{pmatrix}\delta n(\vec{x},t)\\\ \vec{u}(\vec{x},t)\\\
\delta n^{\epsilon}(\vec{x},t)\end{pmatrix}=\begin{pmatrix}\delta
n(\vec{p},\omega)\\\ \vec{u}(\vec{p},\omega)\\\ \delta
n^{\epsilon}(\vec{p},\omega)\end{pmatrix}e^{i\vec{p}\cdot\vec{x}-i\omega t}$
(119)

into the linearized hydrodynamic equations in Eq. (118). This leads to

$\displaystyle-i\omega\delta
n(\vec{p},\omega)+in_{0}pu_{\parallel}(\vec{p},\omega)$ $\displaystyle=$
$\displaystyle 0\;,$
$\displaystyle-i\omega\frac{P_{0}+n^{\epsilon}_{0}}{v_{F}^{2}}\vec{u}(\vec{p},\omega)+i\vec{p}\;\delta
P(\vec{p},\omega)$ $\displaystyle=$ $\displaystyle-\frac{i2\pi\alpha
v_{F}\vec{p}}{p}\delta n(\vec{p},\omega)n_{0}\;,$ $\displaystyle-i\omega\delta
n^{\epsilon}(\vec{p},\omega)+(P_{0}+n_{0}^{\epsilon})ipu_{\parallel}(\vec{p},\omega)$
$\displaystyle=$ $\displaystyle 0\;.$ (120)

Here, we define $u_{\parallel}=\vec{u}\cdot\vec{p}/p$ which gives the
component of the hydrodynamic velocity in the direction of the momentum
$\vec{p}$. We are interested in the longitudinal propagating modes, so we
project the momentum equation on the momentum direction $\vec{p}/p$. This
gives

$-i\omega\frac{P_{0}+n^{\epsilon}_{0}}{v_{F}^{2}}u_{\parallel}(\vec{p},\omega)+ip\;\delta
P(\vec{p},\omega)=-i2\pi\alpha v_{F}\delta n(\vec{p},\omega)n_{0}.$ (121)

Together with the other two equations, we find that

$\begin{pmatrix}-i\omega&in_{0}p&0\\\ i2\pi\alpha
v_{F}n_{0}&-i\omega\frac{P_{0}+n^{\epsilon}_{0}}{v_{F}^{2}}&ip/2\\\
0&ip(P_{0}+n_{0}^{\epsilon})&-i\omega\end{pmatrix}\begin{pmatrix}\delta
n(\vec{p},\omega)\\\ u_{\parallel}(\vec{p},\omega)\\\ \delta
n^{\epsilon}(\vec{p},\omega)\end{pmatrix}=\begin{pmatrix}0\\\ 0\\\
0\end{pmatrix}.$ (122)

These equations have three solutions when the frequencies of the fluctuations
satisfie the dispersion relations

$\displaystyle\omega(\vec{p})$ $\displaystyle=$ $\displaystyle 0$
$\displaystyle\omega_{\pm}(\vec{p})$ $\displaystyle=$
$\displaystyle\pm\sqrt{2\pi\alpha
v^{3}_{F}n_{0}^{2}p/(P_{0}+n^{\epsilon}_{0})+v_{F}^{2}p^{2}/2}.$ (123)

We observe, that in the long wavelength limit, there is a square-root
dispersion which represents the hydrodynamic plasmon. At charge neutrality,
$n_{0}=0$, the dispersion becomes linear representing sound waves with a
velocity given by $v_{F}/\sqrt{2}$.

This result, as we will see later, is in disagreement with the calculation
within the RPA. The RPA calculation predicts the existence of thermal plasmons
at charge neutrality and at non-zero temperature. The key step to reconcile
these results relies on the observation that, for Dirac electrons, the
momentum density and charge current are decoupled. This is in stark contrast
to one-band systems with a parabolic dispersion where the momentum density is
proportional to the charge current. To this end, let us additionally consider
an equation of motion for the charge current, which is obtained by multiplying
the Boltzmann equation, Eq. (89), by the corresponding group velocity
$\partial_{\vec{p}}\epsilon_{\lambda}(\vec{x},\vec{p})=\lambda
v_{F}^{R}\hat{p}$ and integrating the resulting equation over all states. The
time-derivative term yields

$N\sum_{\lambda=\pm}\int\frac{d\vec{p}}{(2\pi)^{2}}\lambda
v_{F}^{R}\hat{p}\partial_{t}f_{\lambda}(\vec{x},t,\vec{p})=\partial_{t}\vec{j}(\vec{x},t).$
(124)

The space-derivative term can be similarly integrated and gives

$\displaystyle
N\sum_{\lambda=\pm}\int\frac{d\vec{p}}{(2\pi)^{2}}\left(v^{R}_{F}\right)^{2}\hat{p}\hat{p}\cdot\partial_{\vec{x}}f_{\lambda}(\vec{x},t,\vec{p})=\partial_{\vec{x}}\cdot\vec{\vec{\Xi}}(\vec{x},t)\;,$

where we define a second-rank tensor according to

$\displaystyle\vec{\vec{\Xi}}(\vec{x},t)=$ $\displaystyle
N\sum_{\lambda=\pm}\int\frac{d\vec{p}}{(2\pi)^{2}}v^{2}_{F}\hat{p}\hat{p}\Big{[}f_{+}(\vec{x},t,\vec{p})+(f_{-}(\vec{x},t,\vec{p})-1)\Big{]}.$

The momentum-derivative term yields

$\displaystyle-N\sum_{\lambda=\pm}\int\frac{d\vec{p}}{(2\pi)^{2}}\lambda
v_{F}^{R}\hat{p}\partial_{\vec{x}}V^{H}(\vec{x},t)\cdot\partial_{\vec{p}}f_{\lambda}(\vec{x},t,\vec{p})=$
$\displaystyle\hskip
142.26378pt\partial_{\vec{x}}V^{H}(\vec{x},t)\cdot\vec{\vec{\Lambda}}(\vec{x},t)\;.$

where another second-rank tensor has components given by

$\displaystyle\Lambda_{ij}(\vec{x},t)=$ $\displaystyle
N\int\frac{d\vec{p}}{(2\pi)^{2}}\;\left(\frac{\delta_{ij}}{p}-\frac{p_{i}p_{j}}{p^{3}}\right)\left[f_{+}(\vec{x},t,\vec{p})+(f_{-}(\vec{x},t,\vec{p})-1)\right].$

In contrast to the conserved quantities discussed in the previous section, the
group velocity is not a collisional invariant. Therefore, the current density
is not conserved by electron-electron interactions. This is particularly true
at charge neutrality. The integration of the collision term is, in general,
very complicated, especially since the distribution function is unknown. For
the purpose of this discussion, we resort to the relaxation-time approximation
and assume that

$\sum_{\lambda=\pm}\int\frac{d\vec{p}}{(2\pi)^{2}}\;\lambda
v_{F}^{R}pC_{\lambda}[f](\vec{p})\approx-\frac{\vec{j}}{\tau}\;.$ (129)

The value of the relaxation time $\tau$ may be approximated by the
corresponding element of the collision operator in the Boltzmann equation (see
for example in Ref. Fritz et al. (2008)). Finally, combining the above
calculations, we find the equation of motion for the charge current as

$\displaystyle\partial_{t}\vec{j}(\vec{x},t)+\partial_{\vec{x}}\vec{\vec{\Xi}}(\vec{x},t)+\partial_{\vec{x}}V^{H}(\vec{x},t)\vec{\vec{\Lambda}}(\vec{x},t)=-\frac{\vec{j}(\vec{x},t)}{\tau}\;.$
(130)

To linear order in $\vec{u}$, we find that
$\Xi_{ij}=nv_{F}^{2}/2\;\delta_{ij}$ which is proportional to the charge
density and $\Lambda_{ij}=\frac{N}{4\pi}\mathcal{N}\delta_{ij}$ where
$\mathcal{N}=\int
dp\left[f_{+}^{0}(\vec{p})-(f^{0}_{-}(\vec{p})-1)\right]=T\log(2+2\cosh\mu/T)$.
Here $f_{\lambda}^{0}(\vec{p})=(1+\exp(\frac{\lambda v_{F}p-\mu}{T}))^{-1}$ is
the Fermi-Dirac distribution function.

We are again interested in small density fluctuations around the constant and
homogeneous value. Consequently, we assume $n=n_{0}+\delta n$ and small
$\vec{u}$. We insert this solution into the current Eq. (130) and keep the
terms to linear order in the fluctuations. This gives

$\displaystyle
n_{0}\partial_{t}\vec{u}(\vec{x},t)+\frac{v_{F}^{2}}{2}\partial_{\vec{x}}\delta
n(\vec{x},t)+\partial_{\vec{x}}V^{H}(\vec{x},t)\frac{N}{4\pi}\mathcal{N}$
(131) $\displaystyle=-\frac{n_{0}\vec{u}(\vec{x},t)}{\tau}.$

The solution to this equation are propagating waves of the form

$\begin{pmatrix}\delta n(\vec{x},t)\\\
\vec{u}(\vec{x},t)\end{pmatrix}=\begin{pmatrix}\delta n(\vec{p},\omega)\\\
\vec{u}(\vec{p},\omega)\end{pmatrix}e^{i\vec{p}\cdot\vec{x}-i\omega t}.$ (132)

We insert this solution into the linearized current Eq. (131) and obtain

$\displaystyle-i\omega
n_{0}\vec{u}(\vec{p},\omega)+i\frac{v_{F}^{2}}{2}\vec{p}\delta
n(\vec{p},\omega)+i\frac{2\pi\alpha v_{F}\vec{p}}{p}\delta
n(\vec{p},\omega)\frac{N}{4\pi}\mathcal{N}$
$\displaystyle=-\frac{n_{0}\vec{u}(\vec{p},\omega)}{\tau}.$

Next, we project the equation on the momentum direction $\vec{p}/p$, This
gives

$\displaystyle-i\omega
n_{0}u_{\parallel}(\vec{p},\omega)+i\frac{v_{F}^{2}}{2}p\delta
n(\vec{p},\omega)+i2\pi\alpha v_{F}\delta
n(\vec{p},\omega)\frac{N}{4\pi}\mathcal{N}$
$\displaystyle=-\frac{n_{0}u_{\parallel}(\vec{p},\omega)}{\tau}.$

Together with the continuity equation, we find that

$\begin{pmatrix}-i\omega&in_{0}p\\\ i\alpha
v_{F}\frac{N}{2}\mathcal{N}+i\frac{v_{F}^{2}}{2}p&-i\omega
n_{0}+n_{0}/\tau\end{pmatrix}\begin{pmatrix}\delta n(\vec{p},\omega)\\\
u_{\parallel}(\vec{p},\omega)\end{pmatrix}=\begin{pmatrix}0\\\
0\end{pmatrix}.$ (135)

For this equation to be valid, the frequency of the density fluctuations has
to satisfy the dispersion relation

$\displaystyle\omega_{\pm}(\vec{p})=-\frac{i}{2\tau}$
$\displaystyle\pm\sqrt{\frac{N}{2}\alpha
v_{F}Tp\log(2+2\cosh\mu/T)+\frac{v_{F}^{2}p^{2}}{2}-\frac{1}{4\tau^{2}}}.$

Note that, in the long-wavelength limit, we recover the previous hydrodynamic
result.

Figure 8: The figure shows the real and imaginary part of the frequency of the
density fluctuations as a function of momentum at the charge neutrality point.
In the plot, the momentum and frequency are in the units of $T$. The values
the parameters are chosen for illustrative purpose: $\alpha=0.3$, $1/\tau=T$,
and $N=4$. In the case of graphene $N=4$ counting spin and valley degrees of
freedom.

##  V Part B: Electron-hole-plasmon hydrodynamics in the strong coupling limit

In the previous section we focused on the weak-interaction limit of the action
in Eq. (19). We analyzed the theory by a straightforward pertubative expansion
in the coupling constant. We now consider a system of fermions interacting
strongly via long-range Coulomb interactions. It is well known, that the
interactions between electrons can generate plasma oscillations. Under certain
circumstances, these plasma oscillations act as proper quasi-particles, as we
will show below. Describing these oscillations starting from Eq. (19) requires
to go beyond pure perturbation theory and to resort to a resummation scheme,
such as the random-phase approximation (RPA) Pines and Bohm (1952). Formally,
this can be achieved by a Hubbard-Stratonovich transformation Kleinert1978,
which is the formulation we choose here. The results in the following sections
are the most important new results of this paper.

###  V.1 The effective field theory: The random-phase approximation

We introduce a real scalar boson field $\phi_{a}(\vec{r},t)$ to decouple the
quartic Coulomb interaction using the Hubbard-Stratonovich identity

$\displaystyle\exp\left(-\frac{i}{2}\int
dtdt^{\prime}d\vec{x}d\vec{x}^{\prime}\rho^{a}(\vec{x},t)D_{0,ab}(\vec{x},\vec{x}^{\prime},t,t^{\prime})\rho^{b}(\vec{x}^{\prime},t^{\prime})\right)$
$\displaystyle\hskip 14.22636pt=\int\mathcal{D}\phi\exp\Bigg{(}\frac{i}{2}\int
dtdt^{\prime}d\vec{x}d\vec{x}^{\prime}\phi_{a}(\vec{x},t)D_{0,ab}^{-1}\phi_{b}(\vec{x}^{\prime},t^{\prime})$
$\displaystyle\hskip 42.67912pt-i\int
dtd\vec{x}\phi_{a}(\vec{x},t)\rho^{a}(\vec{x},t)\Bigg{)}\;.$ (137)

All the manipulations in this section are performed on the Schwinger-Keldysh
closed time contour, meaning the indices $a$ and $b$ are the previously
defined Keldysh indices. In writing down Eq. (V.1), we absorb an irrelevant
normalization constant into the functional integration measure. The real
scalar field is conjugate to the electron density and therefore directly
captures the dynamics of the plasmons. Therefore, we henceforth refer to this
boson field as the plasmon field. Inserting the identity of Eq. (V.1) into the
partition function of Eq. (10) leads to

$Z=\int\mathcal{D}\psi^{\dagger}\mathcal{D}\psi\mathcal{D}\phi\exp(iS[\psi^{\dagger},\psi,\phi])\;,$
(138)

where the action reads

$\displaystyle S[\psi^{\dagger},\psi,\phi]$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}dtdt^{\prime}\Big{[}\int
d\vec{x}d\vec{x}^{\prime}\psi_{a,\lambda}^{\dagger}(\vec{x},t)\Big{(}G^{-1}_{0,ab;\lambda\lambda^{\prime}}-\gamma^{c}_{ab}\delta_{\lambda\lambda^{\prime}}\phi_{c}(\vec{x},t)\delta(\vec{x}-\vec{x}^{\prime})\delta(t-t^{\prime})\Big{)}\psi_{b,\lambda^{\prime}}(\vec{x}^{\prime},t^{\prime})$
$\displaystyle+\frac{1}{2}\int
d\vec{x}d\vec{x}^{\prime}\phi_{a}(\vec{x},t)D^{-1}_{0,ab}(\vec{x},\vec{x}^{\prime},t,t^{\prime})\phi_{b}(\vec{x}^{\prime},t^{\prime})\Big{]}.$

Philosophically, we have traded a theory of electrons interacting amongst
themselves for a field theory where electrons interact with the electric
potential and the plasmon field.

####  V.1.1 Green functions

The bare inverse Green function of the boson field reads
$D_{0,ab}^{-1}(\vec{x},\vec{x}^{\prime},t,t^{\prime})=4\epsilon\sigma^{x}_{ab}\delta(\vec{x}-\vec{x}^{\prime})\delta(t-t^{\prime})\sqrt{-\nabla^{2}}/e^{2}$.
The square root of the Laplacian, $\sqrt{-\nabla^{2}}$, can be understood in
the following way: it is the inverse Fourier transform of the absolute value
of the momentum, $p\equiv|\vec{p}|$. The Fourier transform of
$D^{-1}_{0,ab}(\vec{x},\vec{x}^{\prime},t,t^{\prime})$ consequently is given
by

$D^{-1}_{0,ab}(\vec{p},\omega)=\sigma^{x}_{ab}2V^{-1}(\vec{p})=\sigma^{x}_{ab}2p/2\pi\alpha
v_{F}\;.$ (140)

There are two things worthwhile noting here. First, this zeroth-order Green
function has no dynamics. The dynamics will only be generated upon integrating
out fermions or, equivalently, in perturbation theory. Second, it comes with a
factor of $2$ due to our choice of convention that $D_{0,ab}\propto
V(\vec{x}-\vec{x}^{\prime})/2$, see Eq. (28). There also is the fermionic
propagator, that we have to evaluate in the effective field theory of the
plasmons. To obtain this field theory, we have to integrate out the fermions.
This suggests that the fermionic propogator is given by the non-interacting
one and all the renormalization effects are in the plasmon sector. This,
however, is not true, and the generated dynamics feeds back into the fermion
sector. To see this, it proves advantageous to introduce source terms in the
action in Eq. (V.1), according to

$S_{J}[\psi^{\dagger},\psi]=\int
dtd\vec{x}\left(\psi_{a}^{\dagger}(\vec{x},t)J_{a}(\vec{x},t)+J_{a}^{\dagger}(\vec{x},t)\psi_{a}(\vec{x},t)\right).$
(141)

This allows to recover the fermionic Green function for any level of
approximation of the plasmon field, even once the electrons are integrated
out. This is very important in the section about the coupled quantum kinetic
equations. The generating function for the fermionic Green function reads

$Z[J,J^{\dagger}]=\int\mathcal{D}\psi^{\dagger}\mathcal{D}\psi\mathcal{D}\phi\exp(iS[\psi^{\dagger},\psi,\phi]+iS_{J}[\psi^{\dagger},\psi])\;.$
(142)

From this, we can determine the fermion Green function by means of a
functional derivative with respect to the source field, according to

$iG(\vec{r},\vec{r}\;^{\prime},t,t^{\prime})=\frac{-1}{Z[J,J^{\dagger}]}\frac{\delta^{2}\;Z[J,J^{\dagger}]}{\delta
J(\vec{x}^{\prime},t^{\prime})\delta
J^{\dagger}(\vec{x},t)}\Big{|}_{J=J^{\dagger}=0}\;.$ (143)

We continue to integrate out the fermion fields, which gives an effective
theory of the boson field associated with the density fluctuations. We find

$Z[J,J^{\dagger}]=\int\mathcal{D}\phi\exp(iS_{\rm{eff}}[\phi,J,J^{\dagger}]),$
(144)

with the effective action given by

$\displaystyle
S_{\rm{eff}}[\phi,J,J^{\dagger}]=-i\text{Tr}\left[\ln\left(-iG^{-1}\right)\right]$
$\displaystyle+\frac{1}{2}\int
dtdt^{\prime}d\vec{x}d\vec{x}^{\prime}\phi_{a}(\vec{x},t)D^{-1}_{0,ab}(\vec{x},\vec{x}^{\prime}t,t^{\prime})\phi_{b}(\vec{x}^{\prime},t^{\prime}).$
$\displaystyle-\int
dtdt^{\prime}d\vec{x}d\vec{x}^{\prime}J^{\dagger}(\vec{x},t)G(\vec{x},\vec{x}^{\prime},t,t^{\prime};\phi)J(\vec{x}^{\prime},t^{\prime}).$
(145)

where the Green function
$G_{ab;\lambda\lambda^{\prime}}(\vec{x},\vec{x}^{\prime},t,t^{\prime};\phi)$
is a functional of the plasmon field Stoof2009 as

$\displaystyle
G^{-1}_{ab;\lambda\lambda^{\prime}}(\vec{x},\vec{x}^{\prime},t,t^{\prime};\phi)=G^{-1}_{0,ab;\lambda\lambda^{\prime}}(\vec{x},\vec{x}^{\prime},t,t^{\prime})-\gamma_{ab}^{c}\phi_{c}(\vec{x},t)\delta_{\lambda\lambda^{\prime}}\delta(\vec{x}-\vec{x}^{\prime})\delta(t-t^{\prime})\;.$
(146)

While this is formally exact, the presence of the dynamical field requires an
approximation scheme to evaluate it. Using Eq. (143) and the generating
function introduced in Eq. (144), we find that the fermion Green function can
be calculated according to

$\displaystyle
G_{ab;\lambda\lambda^{\prime}}(\vec{x},\vec{x}\;^{\prime},t,t^{\prime})=\int\mathcal{D}\phi\;G_{ab;\lambda\lambda^{\prime}}(\vec{x},\vec{x}\;^{\prime},t,t^{\prime};\phi)\exp(iS_{\rm{eff}}[\phi,0,0])\;.$
(147)

This equation consequently shows in a very explicit manner that the generated
plasmon dynamics feeds back into the fermion dynamics through the field $\phi$
and its associated dynamics encoded in $S_{\rm{eff}}$. Consequently, the next
step is to determine $S_{\rm{eff}}$.

####  V.1.2 The saddle-point equation

The effective action $S_{\rm{eff}}$, introduced formally in Eq. (145), can be
obtained after integrating out the fermions. It is exact but also very
complicated. The problem is that the trace cannot be evaluated in an easy
manner due to the presence of the plasmon field in the Green function of the
fermions. Consequently, we require an approximation scheme. The saddle-point
contribution to the partition function is given by the configuration that
minimizes the action $S_{\rm{eff}}[\phi,0,0]$ (formally this manipulation is
equivalent to Hartree or mean-field approximation). This can be obtained from
the condition

$\displaystyle\frac{\delta
S_{\rm{eff}}[\phi,0,0]}{\delta\phi}\Big{|}_{\langle\phi\rangle}=0$ (148)

This directly leads to

$\displaystyle\langle\phi_{c}(\vec{x},t)\rangle=$ $\displaystyle-i\int
dt^{\prime}d\vec{x}^{\prime}D_{0,cd}(\vec{x},\vec{x}^{\prime};t,t^{\prime})G_{ab;\lambda\lambda}(\vec{x}^{\prime};\vec{x}^{\prime},t^{\prime},t^{\prime};\langle\phi\rangle)\gamma^{d}_{ba}.$

This is a self-consistency equation for the local charge density.
Philosophically, Eq. (V.1.2) corresponds to a self-consistent version of the
Hartree diagram previously discussed in Sec. IV.2. The fermion propagator at
this level of approximation reads

$\displaystyle
G^{-1}_{ab;\lambda\lambda^{\prime}}(\vec{x},\vec{x}^{\prime},t,t^{\prime};\langle\phi\rangle)=G^{-1}_{0,ab;\lambda\lambda^{\prime}}(\vec{x},\vec{x}^{\prime},t,t^{\prime})-\gamma_{ab}^{c}\langle\phi_{c}(\vec{x},t)\rangle\delta_{\lambda\lambda^{\prime}}\delta(\vec{x}-\vec{x}^{\prime})\delta(t-t^{\prime})\;.$
(150)

The expectation value of the plasmon field can now be identified with the
self-energy within the Hartree approximation, already given in Eq. (65). It
thus recovers the dispersion of Eq. (86). We proceed to expand the plasmon
field in deviations from the mean-field value, i.e.,
$\phi_{a}(\vec{x},t)=\langle\phi_{a}(\vec{x},t)\rangle+\phi_{a}^{\prime}(\vec{x},t)$
where $\langle\phi_{a}(\vec{x},t)\rangle$ is the saddle-point configuration
and $\phi_{a}^{\prime}$ is associated with fluctuations around the saddle
point. As a result, we have

$\displaystyle
G^{-1}_{ab;\lambda\lambda^{\prime}}(\vec{x},\vec{x}\,^{\prime},t,t^{\prime};\phi)=G^{-1}_{ab;\lambda\lambda^{\prime}}(\vec{x},\vec{x}^{\prime},t,t^{\prime};\langle\phi\rangle)-\gamma_{ab}^{\alpha}\phi^{\prime}_{\alpha}(\vec{x},t)\delta_{\lambda\lambda^{\prime}}\delta(\vec{x}-\vec{x}^{\prime})\delta(t-t^{\prime})\;.$
(151)

We proceed to expand the effective action, Eq. (145), to second order in the
fluctuations $\phi^{\prime}$. Substituting Eq. (151) into Eq.(145) and using
the series expansion of the logarithm schematically, we suppress the Keldysh
indices here, we find to second order in fluctuations that

$\displaystyle S_{\rm{eff}}[\phi,0,0]$ $\displaystyle=$
$\displaystyle-i\text{Tr}\left[\ln\left(-iG^{-1}(\langle\phi\rangle+\phi^{\prime})\right)\right]+\frac{1}{2}\int
dtdt^{\prime}d\vec{x}d\vec{x}^{\prime}\left(\langle\phi_{a}(\vec{x},t)\rangle+\phi^{\prime}_{a}(\vec{x},t)\right)D^{-1}_{0,ab}(\vec{x},\vec{x}^{\prime}t,t^{\prime})\left(\langle\phi_{b}(\vec{x}^{\prime},t^{\prime})\rangle+\phi^{\prime}_{b}(\vec{x}^{\prime},t^{\prime})\right)$
$\displaystyle\approx$
$\displaystyle-i\text{Tr}[\ln(-iG^{-1}(\langle{\phi}\rangle))]+i\text{Tr}[G(\langle{\phi}\rangle)\phi^{\prime}]+\int
dtdt^{\prime}d\vec{x}d\vec{x}^{\prime}\langle\phi_{a}(\vec{x},t)\rangle
D^{-1}_{0,ab}(\vec{x},\vec{x}^{\prime};t,t^{\prime})\phi^{\prime}_{b}(\vec{x}^{\prime},t^{\prime})$
$\displaystyle+i\text{Tr}[\frac{1}{2}G(\langle\phi\rangle)\phi^{\prime}G(\langle{\phi}\rangle)\phi^{\prime}].$

Here, we use for brevity the shorthand notation
$G^{-1}(\langle{\phi}\rangle)\equiv
G^{-1}_{ab;\lambda\lambda^{\prime}}(\vec{x},\vec{x}^{\prime},t,t^{\prime};\langle\phi\rangle)$.
The first term gives an irrelevant constant which will be absorbed in the
functional integration measure. The linear terms in the fluctuations sum to
zero at the saddle point. Their cancelation is equivalent to the saddle-point
condition in Eq. (V.1.2). The remaining term, consequently, is the last term
in Eq. (V.1.2). It is the term that accounts for quadratic fluctuations around
the saddle point. At the same time, however, it determines the effective
plasmon propagator with the RPA approximation according to

$\displaystyle S_{\rm{eff}}[\phi^{\prime}]=$
$\displaystyle\frac{1}{2}\sum_{a,b=1,2}\int
dtdt_{1}d\vec{x}d\vec{x}^{\prime}\phi^{\prime}_{a}(\vec{x},t)D^{-1}_{ab}(\vec{x},\vec{x}^{\prime};t,t^{\prime})\phi^{\prime}_{b}(\vec{x}^{\prime},t^{\prime})\;,$

where the inverse plasmon Green function satisfies

$D^{-1}_{ab}(\vec{x},\vec{x}^{\prime},t,t^{\prime})=D^{-1}_{0;ab}(\vec{x},\vec{x}^{\prime};t,t^{\prime})-\Pi_{ab}(\vec{x},\vec{x}^{\prime},t,t^{\prime}).$
(154)

The self-energy $\Pi$ has the diagrammatic representation shown in Fig. 10(b).
It is commonly referred to as the polarization diagram. The corresponding
algebraic expression reads

$\Pi_{ab}=-iN\text{Tr}[\gamma^{a}G(\vec{x},\vec{x}^{\prime},t,t^{\prime},\langle\phi\rangle)\gamma^{b}G(\vec{x}^{\prime},\vec{x},t^{\prime},t,\langle\phi\rangle)].$
(155)

The retarded component of Eq. (154) gives dynamics to the plasmon: it allows
to determine the dispersion and decay rate of the plasmons in the next
section. Its Keldysh components has the form of Eq. (42) which will be the
starting point for the derivation of the Boltzmann equation for the plasmons.
Now that we have the above effective action of the plasmon and its dynamics,
it is time to return to the fermion Green function
$G_{ab;\lambda\lambda^{\prime}}(\vec{x},\vec{x}^{\prime},t,t^{\prime})$ in Eq.
(147). We can graphically represent the fermion Green function expansion in
terms of $\phi^{\prime}$ according to Fig. 9. We now have to ‘average’ this
fermion propagator over the Gaussian action of the bosons and resum it. It
turns out that we have to choose the series corresponding to the Fock-like
diagram, Fig. 10(a), to obtain a conserving approximation KadanoffBaym1961.
This is also known as the GW approximation. The resulting Green function is
the solution of the Dyson equation

$\displaystyle
G^{-1}_{ab;\lambda\lambda^{\prime}}(\vec{x},\vec{x}^{\prime},t,t^{\prime})=$
$\displaystyle\hskip
22.76228ptG^{-1}_{ab;\lambda\lambda^{\prime}}(\vec{x},\vec{x}^{\prime};t,t^{\prime};\langle\phi\rangle)-\Sigma_{ab;\lambda\lambda^{\prime}}(\vec{x},\vec{x}^{\prime},t,t^{\prime})\;,$
(156)

where

$\displaystyle\Sigma_{ab;\lambda\lambda^{\prime}}(\vec{x},\vec{x}^{\prime},t,t^{\prime})=$
$\displaystyle
i\left(\gamma^{\alpha}G(\vec{x},\vec{x}^{\prime},t,t^{\prime};\langle\phi\rangle)\gamma^{\beta}\right)_{ab;\lambda\lambda^{\prime}}D_{\alpha\beta}(\vec{x},\vec{x}^{\prime},t,t^{\prime})\;.$
(157)

The Dyson equations in Eqs. (154) and (156), together with the self-energies
in Eqs. (155) and (157), are the minimal set of equations that describes the
interplay between the collective modes and the single-particle components of
the interacting Dirac electron. This approximation can be derived from a
single free energy diagram. It was shown in Ref. KadanoffBaym1961 that this
implies that it constitutes a conserving approximation. As such, it respects
the conservation laws of total energy and momentum in the combined system of
electrons and plasmons, as we show explicitly in Sec. V.4.

Figure 9: Diagram representing the solution of the Dyson equation (151) in
terms of a series expansion of the quantum fluctuation $\phi^{\prime}$.

(a) electron self-energy

(b) Boson self-energy

Figure 10: Self-energy for electron and plasmon fields within the RPA
approximation. The set of the GW-diagram and polarization function constitutes
a conserving approximation.

###  V.2 The plasmons

In this section we discuss the plasmon dynamics at non-zero temperature and
non-zero chemical potential. To that end, we analyze the polarization function
numerically. We then proceed to find an approximate analytical description
that we use to determine the plasmon spectrum and the quasi-particle life-
times.

####  V.2.1 Non-zero temperature polarization function

Here, we consider the retarded component of the polarization function, Eq.
(155), at non-zero temperature. After a Wigner transformation, we obtain

$\displaystyle\Pi^{R}(\vec{p},\omega)$
$\displaystyle=-iN\int\frac{d\vec{q}}{(2\pi)^{2}}\frac{d\nu}{2\pi}\Big{[}G^{R}_{\lambda\lambda^{\prime}}(\vec{p}+\vec{q},\omega+\nu,\langle\phi\rangle)G^{K}_{\lambda^{\prime}\lambda}(\vec{q},\nu,\langle\phi\rangle)$
$\displaystyle\hskip
71.13188pt+\;G^{K}_{\lambda\lambda^{\prime}}(\vec{p}+\vec{q},\omega+\nu,\langle\phi\rangle)G^{A}_{\lambda^{\prime}\lambda}(\vec{q},\nu,\langle\phi\rangle)\Big{]}.$

Here, for brevity, we suppress the space and time variables of the functions
involved. Next, we use the transformation matrix in Eq. (63) to transform the
objects within the polarization function into the quasi-particle basis. After
that, we integrate over the frequency variable $\nu$ by making use of the
Dirac delta function coming from $g_{0}^{K}$. Finally, by a straightforward
algebraic manipulation, we find the polarization function expressed in the
form of the Lindhard formula

$\displaystyle\Pi^{R}(\vec{p},\omega)=$ $\displaystyle
2N\sum_{\lambda\lambda^{\prime}=\pm
1}\int\frac{d\vec{q}}{(2\pi)^{2}}\frac{\mathcal{F}_{\lambda\lambda^{\prime}}(\vec{p},\vec{q})\left(f_{\lambda}(\vec{q})-f_{\lambda^{\prime}}(\vec{p}+\vec{q})\right)}{\omega+i0^{+}+\lambda
v_{F}^{R}q-\lambda^{\prime}v_{F}^{R}|\vec{p}+\vec{q}|}\;.$ (159)

The coherence factor is defined according to

$\mathcal{F}_{\lambda\lambda^{\prime}}(\vec{p},\vec{q})=\frac{1}{2}(1+\lambda\lambda^{\prime}\cos(\theta_{\vec{p}+\vec{q}}-\theta_{\vec{q}}))\;.$
(160)

Let us note that, strictly speaking, compared to the conventional Lindhard
formula, there is an extra factor $2$ in our result. This is consistent within
our convention that the inverse bare boson Green function comes with the same
factor, i.e., $D_{0}\propto V(\vec{x}-\vec{x}^{\prime})/2$.

The polarization function at zero temperature and arbitrary chemical potential
($\Pi_{\mu,0}$) can be calculated exactly. Here, we focus on non-zero
temperatures. There is no analytical expression for the non-zero temperature
polarization. However, there exists a relation between zero temperature and
non-zero temperature polarization functions Das Sarma and Li (2013),

$\Pi^{R}_{\mu,T}(\vec{p},\omega)=\int_{0}^{\infty}\;d\mu^{\prime}\sum_{\lambda=\pm
1}\frac{\Pi^{R}_{\mu^{\prime},0}(\vec{p},\omega)}{4T\cosh^{2}\left(\frac{\mu^{\prime}+\lambda\mu}{2T}\right)}.$
(161)

Using this, we numerically solve the polarization function at arbitrary
temperatures and subsequently replace it in the Dyson equation in Eq. (154).
This allows us to determine the energy dispersion ($\omega_{p}$) and the decay
rate ($\gamma_{p}$) of the plasmon mode. The plasmon frequency
($\omega=\omega_{p}-i\gamma_{p}$) is obtained by equating the inverse Green
function to zero.

$(D^{-1})^{R}(\vec{p},\omega)=(D_{0}^{-1})^{R}(\vec{p},\omega)-\Pi(\vec{p},\omega)=0.$
(162)

Defined this way, the decay rate $\gamma_{p}$ is positive. If the damping is
sufficiently weak ($\gamma_{p}\ll\omega_{p}$), one can expand the polarization
function to leading order in $\gamma_{p}$

$\displaystyle\Pi(\vec{p},\omega_{p}-i\gamma_{p})$
$\displaystyle\approx\Re\Pi(\vec{p},\omega_{p})-i\gamma_{p}\partial_{\omega}\Re\Pi(\vec{p},\omega)\Big{|}_{\omega=\omega_{p}}+i\Im\Pi(\vec{p},\omega_{p}).$

The energy of the plasmon can be determined from the real part of Eq. (162)

$(D_{0}^{-1})^{R}(\vec{p},\omega_{p})-\Re\Pi(\vec{p},\omega_{p})=0\;,$ (164)

whereas the decay rate is a solution of the imaginary part

$\gamma_{p}=\frac{\Im\Pi(\vec{p},\omega)}{\partial_{\omega}\Re\Pi(\vec{p},\omega)}\Bigg{|}_{\omega=\omega_{p}}\;.$
(165)

In general, the non-interacting Green function, $D^{-1}_{0}$ can be a function
of both momentum $\vec{p}$ and frequency $\omega$. However, in our case
$D^{-1}_{0}$ describes the bare Coulomb potential which is non-dynamical and
hence does not depend on the frequency variable. Fig. 11 (a) shows the real
part of the polarization function at non-zero chemical potential in the
momentum-frequency plane. Solutions to Eq. (164) exist only when $\Re\Pi>0$.
This is the case in the upper triangle of the plot in Fig. 11 (a). As
discussed before, a stable plasmon requires $\Im\Pi\approx 0$. In Fig. 11(b)
we plot the imaginary part of the polarization function at the same value of
parameters as in Fig. 11 (a). Although it is not zero, it is still negligibly
small in the low-momentum limit. As a result, one may expect a long-wavelength
underdamped plasmon mode with for all pratical purposes almost infinitely long
life time. This implies that plasmons behave like quasi-particles for
practical matters.

In the next section, we find an approximate description of the dispersion of
the plasmon and its decay rate in the low-momentum limit. We furthermore
determine the value of the momentum cutoff beyond which the plasmons are
overdamped.

Figure 11: Polarization function, (a) real part and (b) imaginary part, at a
non-zero temperature and chemical potential. We show the polarization function
when the chemical potential $\mu/T$ = 1.

Figure 12: Spectral function of the plasmon at a non-zero temperature and
chemical potential $\mu/T=4$. The value of fine structure constant is chosen
to mimic graphene device sandwiched in hBN ($\alpha=0.3$).

.

Figure 13: Spectral function of the plasmons at various electron chemical
potential at a small value of momentum $q/T=2.5$

. Here the fine structure constant $\alpha=0.3$.

####  V.2.2 Analytical approximation

In this section, we will evaluate the energy spectrum and decay rate of the
plasmon in the long-wavelength limit. To this end, we first decompose the
polarization function into a sum of two terms accordingly to
$\Pi^{R}(\vec{p},\omega)=\Pi^{R}_{+}(\vec{p},\omega)+\Pi^{R}_{-}(\vec{p},\omega)$.
The first term, $\Pi_{+}(\vec{p},\omega)$, describes the contribution from
intraband particle-hole pairs ($\lambda^{\prime}=\lambda$), whereas
$\Pi_{-}(\vec{p},\omega)$ comes from the interband transitions
($\lambda^{\prime}=-\lambda$). The imaginary part of the polarization function
accordingly reads

$\Im\Pi^{R}(\vec{p},\omega)=\Im\Pi^{R}_{+}(\vec{r},t)+\Im\Pi_{-}(\vec{p},\omega)$
(166)

and it amounts to a decay rate of the plasmon mode with

$\displaystyle\Im\Pi^{R}_{\pm}(\vec{p},\omega)=-2N\pi\sum_{\lambda=\pm
1}\int\frac{d\vec{q}}{(2\pi)^{2}}\mathcal{F}_{\lambda\;\pm\lambda}(\vec{p}+\vec{q},\vec{q})$
$\displaystyle\Big{[}f_{\lambda}(\vec{q})-f_{\pm\lambda}(\vec{p}+\vec{q})\Big{]}\delta(\omega+\epsilon_{\lambda}(\vec{q})-\epsilon_{\pm\lambda}(\vec{p}+\vec{q})).$
(167)

The conservation of energy enters through the delta function with the argument
$\omega+\lambda v_{F}^{R}q-\lambda^{\prime}v_{F}^{R}|\vec{p}+\vec{q}|$. This
function forms either an ellipse or a hyperbola in the $q_{x}q_{y}$-plane. It
may be easy to see this by means of a transformation of the momentum variables
into the elliptic coordinate system, $(\theta,\mu)$, where $\theta\in[0,2\pi)$
and $\mu\in[0,\infty)$. When $\lambda=-$ and $\lambda^{\prime}=+$, this
equation forms an ellipse with the value of $\mu$ determined by
$\cosh\mu=\omega/p$. The size of the ellipse depends on the value of
$\cosh\mu\leq 1$. We find that when $\omega$ is slightly bigger than $p$, the
available phase space is restricted since the size of the ellipse is small
whereas when $\omega$ is much bigger than $p$, the available phase space in
turn grow bigger and allows the decay process more likely to occur. However,
when $\lambda^{\prime}=\lambda$, this equation becomes an equation of a
hyperbola centered at $(-p_{x}/2,-p_{y}/2)$. The width of the hyperbola is
determined from the ratio of the frequency to the momentum given by
$\cos\theta=\lambda\omega/p$. Therefore, the available phase space for a
plasmon with the energy $\omega$ to decay into two fermions with the energies
$\lambda p$ and $\lambda|\vec{p}+\vec{q}|$ is extended. This process is thus
the main mechanism for plasmon decay. We expand all the quantities appearing
in the polarization function up to first order in $\vec{p}$. This gives

$\displaystyle\mathcal{F}_{\pm\pm}(\vec{p},\vec{q})$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(1+\cos\left(\theta_{\vec{p}+\vec{q}}-\theta_{\vec{q}}\right)\right)\approx
1,$ $\displaystyle\mathcal{F}_{\pm\mp}(\vec{p},\vec{q})$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(1-\cos\left(\theta_{\vec{p}+\vec{q}}-\theta_{\vec{q}}\right)\right)\approx\frac{1}{4}\left(\vec{p}\cdot\vec{\nabla}_{\vec{q}}\theta_{\mathbf{q}}\right)^{2},$
$\displaystyle f_{\lambda}(\vec{p}+\vec{q})$ $\displaystyle\approx$
$\displaystyle
f_{\lambda}(\vec{q})+\vec{p}\cdot\vec{\nabla}_{\vec{q}}f_{\lambda}(\vec{q}),$
$\displaystyle\epsilon_{\lambda}(\vec{p}+\vec{q})$ $\displaystyle\approx$
$\displaystyle\epsilon_{\lambda}(\vec{q})+\vec{p}\cdot\vec{\nabla}_{\vec{q}}\epsilon_{\lambda}(\vec{q}).$
(168)

By substituting Eq. (168) into the imaginary part of the polarization followed
by a straightforward calculation, we find

$\Im\Pi^{R}(\vec{p},\omega)\approx-\frac{2N}{16}\frac{p^{2}}{\omega}\Big{(}f_{+}(|\omega/2|)-f_{-}(|\omega/2|)\Big{)},$
(169)

which provides the main contribution to the decay rate whereas
$\Im\Pi^{R}_{+}(\vec{p},\omega)$ gives an unimportant correction. In thermal
equilibrium, the distribution function becomes the Fermi-Dirac distribution
function at a temperature $T$ and chemical potential $\mu$, i.e.,
$f_{\lambda}(\omega)=\frac{1}{\exp((\omega-\mu)/T)+1)}$. We present more
details of the calculation in Appendix D.

At zero temperature, this becomes

$\Im\Pi_{-}^{R}(\vec{p},\omega)\approx-\frac{2N}{16}\frac{p^{2}}{\omega}\Theta(|\omega|-2|\mu|).$
(170)

It vanishes when $|\omega|<2|\mu|$, meaning in that region a long-lived
plasmon mode exists. The real part of the polarization function is given by

$\displaystyle\Re\Pi_{\pm}^{R}(\vec{p},\omega)$
$\displaystyle=2N\sum_{\lambda=\pm
1}\int\frac{d\vec{q}}{(2\pi)^{2}}\mathcal{F}_{\lambda\;\pm\lambda}(\vec{p},\vec{q})\frac{f_{\lambda}(\vec{q})-f_{\pm\lambda}(\vec{p}+\vec{q})}{\omega+\lambda
v_{F}^{R}q\mp\lambda v_{F}^{R}|\vec{p}+\vec{q}|}.$

We substitute these expressions into Eq. (168) followed by expanding its
denominator to first order in $p/\omega$. Based on numerics, we expect to find
a stable plasmon mode in this limit. We can approximate the expression as

$\displaystyle\Re\Pi_{+}^{R}(\vec{p},\omega)$ $\displaystyle\approx
2N\sum_{\lambda=\pm
1}\int\frac{d\vec{q}}{(2\pi)^{2}}\frac{-\vec{p}\cdot\vec{\nabla}_{\vec{q}}f_{\lambda}(\vec{q})}{\omega}\Big{[}1+\frac{\vec{p}\cdot\vec{\nabla}_{\vec{q}}\epsilon_{\lambda}(\vec{q})}{\omega}\Big{]},$

which results in

$\Re\Pi^{R}_{+}(\vec{p},\omega)\approx\frac{2Np^{2}}{4\pi\omega^{2}}\mathcal{N},$
(173)

where $\mathcal{N}=T\log\left(2+2\cosh{\mu/T}\right)$. When the Dirac system
is exposed to external perturbations, the polarization function varies in
position and time via the distribution function. In that case
$\mathcal{N}=\int
dq\Big{[}f_{+}(\vec{q})-\left(f_{-}(\vec{q})-1\right)\Big{]}.$ In contrast,
the interband term gives a logarithmic correction which will be neglected from
evaluating the plasmon energy as

$\Re\Pi^{-}_{R}(\vec{p},\omega)\approx\frac{Np^{2}}{2\pi}\int
dq\Big{[}f_{+}(\vec{q})-f_{-}(\vec{q})\Big{]}\Big{[}\frac{1}{4\left(v_{F}^{R}\right)^{2}q^{2}-\omega^{2}}\Big{]}.$
(174)

At non-zero doping and zero temperature, we find that
$\Re\Pi^{R}_{-}(\vec{p},\omega)=\frac{Np^{2}}{16\pi\omega}\log\left(\Big{|}\frac{\omega-2\mu}{\omega+2\mu}\Big{|}\right)$.
If we substitute Eqs. (173) and (140) into Eq. (164), we obtain the dispersion
relation of the plasmon at non-zero $T$,

$\displaystyle\omega_{p}(\vec{p})=\sqrt{\frac{N}{2}\alpha
Tp\log\left(2+2\cosh{\mu/T}\right)}\;,$ (175)

with the decay rate

$\gamma_{p}=\frac{\pi\omega_{p}(\vec{p})^{2}}{16\log(2+2\cosh\mu/T)}\Big{(}f_{+}(\omega_{p}(\vec{p})/2)-f_{-}(\omega_{p}(\vec{p})/2)\Big{)}.$
(176)

(a) Energy dispersion

(b) decay rate

Figure 14: Numerical solution of the Dyson equation (162) at $\mu/T=3$ Here we
choose a relatively small value of the fine structure constant $\alpha=0.3$.