Electronic hydrodynamics in graphene

TL;DR

This paper derives the unique electronic hydrodynamics of graphene using kinetic theory, including magnetic fields and dissipative properties, highlighting differences from traditional relativistic hydrodynamics.

Contribution

It provides a pedagogical derivation of graphene's unconventional hydrodynamics, extending previous work to include magnetic fields and explicit dissipative coefficients.

Findings

Derived explicit expressions for shear viscosity and electrical conductivity.

Compared graphene hydrodynamics with relativistic hydrodynamics.

Extended the theory to include external magnetic fields.

Abstract

In this paper I report a pedagogical derivation of the unconventional electronic hydrodynamics in graphene on the basis of the kinetic theory. While formally valid in the weak coupling limit, this approach allows one to derive the unconventional hydrodynamics in the system which is neither Galilean- nor Lorentz-invariant, such that hydrodynamic equations can not be inferred from symmetry arguments. I generalize earlier work to include external magnetic fields and give explicit expressions for dissipative coefficients, the shear viscosity and electrical conductivity. I also compare the resulting theory with relativistic hydrodynamics.

Equations

$${\ell }_{ee}\ll {\ell }_{\mathrm{dis}},{\ell }_{\mathrm{e}-\mathrm{ph}},{\ell }_R,\mathrm{etc}.$$

$$\mathcal{L}f={\mathrm{St}}_{ee}[f]+{\mathrm{St}}_R[f]-\frac{f-⟨f{⟩}_{\phi }}{{\tau }_{\mathrm{dis}}},$$

$$\mathcal{L}={\partial }_t+\mathbf{v}\cdot {\nabla }_{\mathbf{r}}+(e\mathbf{E}+\frac{e}{c}\mathbf{v}\times \mathbf{B})\cdot {\nabla }_{\mathbf{k}}.$$

$$⟨f{⟩}_{\phi }=\underset{-\pi }{\overset{\pi }{\int }}\frac{d\phi }{2\pi }{f}_{\lambda \mathbf{k}},$$

$$N\sum _{\lambda }\int \frac{d^{2}k}{(2\pi {)}^2}\left[{\mathrm{St}}_{ee}[f]+{\mathrm{St}}_R[f]-\frac{f-⟨f{⟩}_{\phi }}{{\tau }_{\mathrm{dis}}}\right]=0.$$

$$N\sum _{\lambda }\int \frac{d^{2}k}{(2\pi {)}^2}{\partial }_t{f}_{\lambda \mathbf{k}}={\partial }_{t}N\int \frac{d^{2}k}{(2\pi {)}^2}\left[f_{+,\mathbf{k}}-\left(1-f_{-,\mathbf{k}}\right)\right]={\partial }_t(n_{+}-n_{-})={\partial }_{t}n,$$

$$N\sum _{\lambda }\int \frac{d^{2}k}{(2\pi {)}^2}\mathbf{v}\cdot {\nabla }_{\mathbf{r}}{f}_{\lambda \mathbf{k}}={\partial }_{i}N\int \frac{d^{2}k}{(2\pi {)}^2}\left[v_{+,\mathbf{k}}^i{f}_{+,\mathbf{k}}-v_{-,\mathbf{k}}^i\left(1-f_{-,\mathbf{k}}\right)\right]={\partial }_i\left(j_{+}^i-j_{-}^i\right)={\nabla }_{\mathbf{r}}\cdot \mathbf{j},$$

$$e\mathbf{E}\cdot N\sum _{\lambda }\int \frac{d^{2}k}{(2\pi {)}^2}{\nabla }_{\mathbf{k}}{f}_{\lambda \mathbf{k}}=e\mathbf{E}\cdot N\int \frac{d^{2}k}{(2\pi {)}^2}\left[{\nabla }_{\mathbf{k}}{f}_{+\mathbf{k}}-{\nabla }_{\mathbf{k}}\left(1-f_{-\mathbf{k}}\right)\right]=0.$$

$${ϵ}^{\alpha \beta \gamma }{B}^{\gamma }N\sum _{\lambda }\int \frac{d^{2}k}{(2\pi {)}^2}{v}_{\lambda \mathbf{k}}^{\beta }\frac{\partial f_{\lambda \mathbf{k}}}{\partial k^{\alpha }}=-{ϵ}^{\alpha \beta \gamma }{B}^{\gamma }N\sum _{\lambda }\int \frac{d^{2}k}{(2\pi {)}^2}{f}_{\lambda \mathbf{k}}\frac{\partial v_{\lambda \mathbf{k}}^{\beta }}{\partial k^{\alpha }}.$$

$$\partial v_{\lambda \mathbf{k}}^{\beta }/\partial k^{\alpha }\propto {\delta }_{\alpha \beta },{ϵ}^{\alpha \beta \gamma }{\delta }_{\alpha \beta }=0.$$

$${\partial }_{t}n+{\nabla }_{\mathbf{r}}\cdot \mathbf{j}=0,$$

$${\partial }_t{n}_I+{\nabla }_{\mathbf{r}}\cdot {\mathbf{j}}_I=-\frac{n_I-n_I^{(0)}}{{\tau }_R},$$

$$N\sum _{\lambda }\int \frac{d^{2}k}{(2\pi {)}^2}{ϵ}_{\lambda \mathbf{k}}\left[{\mathrm{St}}_{ee}[f]+{\mathrm{St}}_R[f]-\frac{f-⟨f{⟩}_{\phi }}{{\tau }_{\mathrm{dis}}}\right]=0.$$

$$N\sum _{\lambda }\int \frac{d^{2}k}{(2\pi {)}^2}{ϵ}_{\lambda \mathbf{k}}{\partial }_t{f}_{\lambda \mathbf{k}}={\partial }_{t}N\int \frac{d^{2}k}{(2\pi {)}^2}{ϵ}_{\lambda \mathbf{k}}{f}_{\lambda \mathbf{k}}={\partial }_t{n}_E,$$

$$N\sum _{\lambda }\int \frac{d^{2}k}{(2\pi {)}^2}{ϵ}_{\lambda \mathbf{k}}{\mathbf{v}}_{\lambda \mathbf{k}}\cdot {\nabla }_{\mathbf{r}}{f}_{\lambda \mathbf{k}}={\nabla }_{\mathbf{r}}\cdot N\int \frac{d^{2}k}{(2\pi {)}^2}{ϵ}_{\lambda \mathbf{k}}{\mathbf{v}}_{\lambda \mathbf{k}}{f}_{\lambda \mathbf{k}}={\nabla }_{\mathbf{r}}\cdot {\mathbf{j}}_E.$$

$$e\mathbf{E}\cdot N\sum _{\lambda }\int \frac{d^{2}k}{(2\pi {)}^2}{ϵ}_{\lambda \mathbf{k}}{\nabla }_{\mathbf{k}}{f}_{\lambda \mathbf{k}}=-e\mathbf{E}\cdot N\sum _{\lambda }\int \frac{d^{2}k}{(2\pi {)}^2}{\mathbf{v}}_{\lambda \mathbf{k}}{f}_{\lambda \mathbf{k}}=-e\mathbf{E}\cdot \mathbf{j}.$$

$${ϵ}^{\alpha \beta \gamma }{B}^{\gamma }N\sum _{\lambda }\int \frac{d^{2}k}{(2\pi {)}^2}{ϵ}_{\lambda \mathbf{k}}{v}_{\lambda \mathbf{k}}^{\beta }\frac{\partial f_{\lambda \mathbf{k}}}{\partial k^{\alpha }}=-{ϵ}^{\alpha \beta \gamma }{B}^{\gamma }N\sum _{\lambda }\int \frac{d^{2}k}{(2\pi {)}^2}{f}_{\lambda \mathbf{k}}{ϵ}_{\lambda \mathbf{k}}\frac{\partial v_{\lambda \mathbf{k}}^{\beta }}{\partial k^{\alpha }}.$$

$${\partial }_t{n}_E+{\nabla }_{\mathbf{r}}\cdot {\mathbf{j}}_E=e\mathbf{E}\cdot \mathbf{j},$$

$$N\sum _{\lambda }\int \frac{d^{2}k}{(2\pi {)}^2}\mathbf{k}\left({\mathrm{St}}_{ee}[f]+{\mathrm{St}}_R[f]\right)=0.$$

$$N\sum _{\lambda }\int \frac{d^{2}k}{(2\pi {)}^2}\mathbf{k}\frac{f-⟨f{⟩}_{\phi }}{{\tau }_{\mathrm{dis}}}=\frac{1}{{\tau }_{\mathrm{dis}}}N\sum _{\lambda }\int \frac{d^{2}k}{(2\pi {)}^2}\mathbf{k}{f}_{\lambda \mathbf{k}}=\frac{{\mathbf{n}}_{\mathbf{k}}}{{\tau }_{\mathrm{dis}}},$$

$$N\sum _{\lambda }\int \frac{d^{2}k}{(2\pi {)}^2}\mathbf{k}{\partial }_t{f}_{\lambda \mathbf{k}}={\partial }_{t}N\sum _{\lambda }\int \frac{d^{2}k}{(2\pi {)}^2}\mathbf{k}{f}_{\lambda \mathbf{k}}={\partial }_t{\mathbf{n}}_{\mathbf{k}}.$$

$$N\sum _{\lambda }\int \frac{d^{2}k}{(2\pi {)}^2}{k}^{\alpha }{v}_{\lambda \mathbf{k}}^{\beta }{\nabla }_{\mathbf{r}}^{\beta }{f}_{\lambda \mathbf{k}}={\nabla }_{\mathbf{r}}^{\beta }N\sum _{\lambda }\int \frac{d^{2}k}{(2\pi {)}^2}{k}^{\alpha }{v}_{\lambda \mathbf{k}}^{\beta }{f}_{\lambda \mathbf{k}}={\nabla }_{\mathbf{r}}^{\beta }{\Pi }_E^{\alpha \beta }.$$

$$e{E}^{\beta }N\sum _{\lambda }\int \frac{d^{2}k}{(2\pi {)}^2}{k}^{\alpha }{\nabla }_{\mathbf{k}}^{\beta }{f}_{\lambda \mathbf{k}}=e{E}^{\beta }N\int \frac{d^{2}k}{(2\pi {)}^2}{k}^{\alpha }{\nabla }_{\mathbf{k}}^{\beta }\left[f_{+,\mathbf{k}}-(1-f_{-,\mathbf{k}})\right]$$

$$=-e{E}^{\beta }N\int \frac{d^{2}k}{(2\pi {)}^2}\left[f_{+,\mathbf{k}}-(1-f_{-,\mathbf{k}})\right]{\nabla }_{\mathbf{k}}^{\beta }{k}^{\alpha }=-en{E}^{\alpha }.$$

$${ϵ}^{\alpha \beta \gamma }{B}^{\gamma }N\sum _{\lambda }\int \frac{d^{2}k}{(2\pi {)}^2}{k}^{\mu }{v}_{\lambda \mathbf{k}}^{\beta }\frac{\partial f_{\lambda \mathbf{k}}}{\partial k^{\alpha }}=-{ϵ}^{\alpha \beta \gamma }{B}^{\gamma }N\sum _{\lambda }\int \frac{d^{2}k}{(2\pi {)}^2}{f}_{\lambda \mathbf{k}}\left[{\delta }_{\mu \alpha }{v}_{\lambda \mathbf{k}}^{\beta }+k^{\mu }\frac{\partial v_{\lambda \mathbf{k}}^{\beta }}{\partial k^{\alpha }}\right]=-\mathbf{j}\times \mathbf{B},$$

$${\partial }_t{n}_{\mathbf{k}}^{\alpha }+{\nabla }_{\mathbf{r}}^{\beta }{\Pi }_E^{\alpha \beta }-en{E}^{\alpha }-\frac{e}{c}{\left[\mathbf{j}\times \mathbf{B}\right]}^{\alpha }=-\frac{n_{\mathbf{k}}^{\alpha }}{{\tau }_{\mathrm{dis}}},$$

$${\partial }_{t}n+{\nabla }_{\mathbf{r}}\cdot \mathbf{j}=0,$$

$${\partial }_t{n}_I+{\nabla }_{\mathbf{r}}\cdot {\mathbf{j}}_I=-\frac{n_I-n_I^{(0)}}{{\tau }_R},$$

$${\partial }_t{n}_E+{\nabla }_{\mathbf{r}}\cdot {\mathbf{j}}_E=e\mathbf{E}\cdot \mathbf{j},$$

$${\partial }_t{n}_{\mathbf{k}}^{\alpha }+{\nabla }_{\mathbf{r}}^{\beta }{\Pi }_E^{\alpha \beta }-en{E}^{\alpha }-\frac{e}{c}{\left[\mathbf{j}\times \mathbf{B}\right]}^{\alpha }=-\frac{n_{\mathbf{k}}^{\alpha }}{{\tau }_{\mathrm{dis}}}.$$

$$f_{\lambda \mathbf{k}}^{(0)}(\mathbf{r})={\left\{1+\exp \left[\frac{{ϵ}_{\lambda \mathbf{k}}-{\mu }_{\lambda }(\mathbf{r})-\mathbf{u}(\mathbf{r})\cdot \mathbf{k}}{T(\mathbf{r})}\right]\right\}}^{-1},$$

$${ϵ}_{\lambda \mathbf{k}}=\lambda v_{g}k,$$

$${\mathbf{v}}_{\lambda \mathbf{k}}=\lambda v_g\frac{\mathbf{k}}{k},{\mathbf{e}}_{\mathbf{v}}={\mathbf{e}}_{\mathbf{k}},\mathbf{k}=\frac{\lambda k}{v_g}{\mathbf{v}}_{\lambda \mathbf{k}}=\frac{{ϵ}_{\lambda \mathbf{k}}{\mathbf{v}}_{\lambda \mathbf{k}}}{v_g^2}.$$

Full text

††journal: Annals of Physics

Electronic hydrodynamics has evolved into a fast paced field with multiple experimental and theoretical groups working to establish observable effects of hydrodynamic behavior of electronic systems across a wide range of materials [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11] (for a comprehensive list of references see recent reviews [12, 13]). Similarly to the usual hydrodynamics of ordinary fluids, one can arrive at the final set of hydrodynamic equations in several complementary ways. One way is purely phenomenological [14] : one writes an equation of motion for a small element of the fluid, a continuity equation expressing conservation of “mass” (or “matter”), and (at least for an ideal fluid) an adiabaticity condition. The resulting five (in three dimensions) differential equations determine the five macroscopic quantities characterizing the (one-component) fluid, e.g., the velocity vector, fluid density, and pressure. Another way of deriving the hydrodynamic equations [15] is also phenomenological in nature, but is somewhat more general since it relies on thermodynamics rather than on classical mechanics of a fluid element. Finally, the hydrodynamic equations can be derived from a “microscopic” Boltzmann theory [16]. While more “technical” and thus conceptually straightforward, this method suffers from the narrow applicability of the Boltzmann equation itself (i.e. the assumption of the constituent particles being free in between successive collisions). However, one typically extends the validity region of the hydrodynamic theory beyond that of the kinetic equation by assumption of universality: all fluids with the same symmetry properties obey the same set of the hydrodynamic equations regardless of whether the coupling between constituent particles is strong (e.g., in water) or weak (e.g., in a dilute gas).

As a macroscopic theory of strongly interacting systems, hydrodynamics might appear to be extremely attractive for condensed matter theorists routinely dealing with problems where strong correlations invalidate simple theoretical approaches. However, with the exception of the early work by Gurzhi [17], not so much attention was given to this theory until recently. The reason for this is fairly simple: unlike water molecules, electrons in solids exist in the environment created by a crystal lattice and hence their momentum is not conserved. As a result, electron motion is typically diffusive, unless the sample size is smaller than the mean free path in which case the system is ballistic.

For most typical scattering mechanisms in solids the mean free path is strongly temperature dependent. At low temperatures, the electrons scatter mostly on lattice imperfections (or “disorder”) leading to, e.g., the residual resistance in metals. At high temperatures the main scattering mechanism is the electron-phonon interaction. In many conventional (or “simple”) metals at least one of these two scattering mechanisms is more effective than electron-electron interaction. In terms of the associated length scales, this statement can be formulated as ${\ell_{ee}\gg\ell_{\rm{dis}},\ell_{\rm e-ph}}$ (with the self-evident notation). On the other hand, if a material would exist where the opposite condition were satisfied at least in some non-negligible temperature range, then one could be justified in neglecting the momentum non-conserving processes and applying the hydrodynamic theory. For a long time such a material was not known and as a result most people working in condensed matter physics were not interested in hydrodynamics. In recent years, the situation has drastically changed as several extremely pure materials became available bringing electronic hydrodynamics within experimental reach. The best known such material is graphene [12, 13, 18, 6, 1, 2, 3, 4, 5, 10].

The purpose of this paper is to provide a pedagogical derivation of the unconventional hydrodynamics in graphene. Low energy excitations in graphene have linear dispersion and hence are not Galilean-invariant. Their motion is restricted to the two-dimensional (2D) plane of the graphene layer, but they are coupled via the classical, three-dimensional (3D) Coulomb’s interaction, such that the electronic system in graphene is not Lorentz-invariant either. As a result, one cannot simply apply either the usual or relativistic hydrodynamics. Instead, one has to derive the hydrodynamics equations for the electronic system in graphene from scratch [19, 20, 21, 22, 23, 24, 25, 26, 27, 28]. The resulting set of equations forms the “unconventional” hydrodynamics in graphene.

1 Kinetic theory of Dirac fermions in graphene

For the purposes of this paper, I will assume the existence of a parameter range where the low energy excitations in graphene can be described by a kinetic (Boltzmann) equation (see Ref. [29] for a detailed discussion and derivation from the quantum many-body theory). I will further assume that at least some part of this parameter range overlaps with the applicability region of the hydrodynamic theory. In that region, the hydrodynamic equations can be derived from the kinetic theory. The extension of the resulting theory beyond the applicability region of the kinetic equation can then be justified by the assumption of universality of the hydrodynamic approach.

While the kinetic theory is often used to describe electronic transport in doped graphene [30, 31], it was shown to be inapplicable only in a region around the neutrality point [32, 33] that is exponentially small in the dimensionless conductivity (in units $e^{2}/h$) of graphene. Recent measurements [3] show this quantity to be of order $10$ at the “hydrodynamic” temperatures, $T>100$K, justifying the use of the kinetic approach to neutral graphene at such (relatively high) temperatures.

The necessary condition for the validity of the hydrodynamics is that the electron-electron interaction is the dominant scattering mechanism in the system such that the typical length scale corresponding to electron-electron interaction, $\ell_{ee}$, is the shortest length in the problem

[TABLE]

Here $\ell_{\rm dis}$, $\ell_{\rm e-ph}$, and $\ell_{R}$ are the length scales characterizing disorder scattering, electron-phonon interaction, and quasiparticle recombination [22] processes. All other scattering mechanisms are encoded in “etc”. Lowering temperature towards zero, $\ell_{ee}$ is expected to diverge, while $\ell_{\rm dis}$ is not. As a result, the inequality (1) can be expected to be fulfilled at temperatures which are high enough to justify the use of the kinetic equation even at charge neutrality (but not too high, so that the electron-phonon interaction could still be considered as subleading).

The general form of the kinetic equation [16] can be seen as a formal equality between the Liouville’s operator and the collision integral. In a two-band electronic system, the kinetic equation can be written as

[TABLE]

Labeling single-particle states by the band index ${\lambda=\pm}$ and the momentum $\boldsymbol{k}$, one can denote the distribution function by ${f=f_{\lambda\boldsymbol{k}}}$. The collision integral comprises three parts: ${\rm{St}}_{ee}[f]$ describes electron-electron interaction, ${\rm{St}}_{R}[f]$ – electron-hole recombination, while the remaining term in Eq. (2a) describes disorder scattering. The latter involves angular averaging defined as

[TABLE]

where $\varphi$ is the polar angle describing the direction of $\boldsymbol{k}$.

The form of the Liouville’s operator is independent of whether the underlying microscopic physics is classical or quantum, at least as long as there is no spin-orbit interaction. The collision integral is more sensitive to the microscopic details of the system. In particular, the $\tau$-approximation employed in Eq. (2a) to describe disorder scattering is almost certainly an oversimplification. Even then, $\tau_{\rm dis}$ is a model-dependent function of energy [18]. However in the limit of weak disorder, required by Eq. (1), one may assign a particular large value to $\tau_{\rm{dis}}$ (as determined by the temperature, $T$, and chemical potential, $\mu$) such that most physical observables (with the notable exception of thermal conductivity) will be insensitive to the choice of the impurity model. Furthermore, the collision integral in Eq. (2a) disregards any “quantum” or “interference” corrections to quasiparticle transport [34, 35, 36, 37, 38, 39, 40, 41, 42]. Below, I treat ${\rm{St}}_{ee}[f]$ at the Golden Rule level. While affecting the numerical values of theoretical estimates for dissipative coefficients, this approximation has no bearing on the form of hydrodynamic equations, which is the main goal of this derivation.

In this paper, I derive the hydrodynamics equations in graphene following the standard textbook steps [16]: (i) integrating the kinetic equation (2), I obtain continuity equations expressing conservation of the particle number (or electric charge), energy, and momentum; (ii) assuming local equilibrium, I relate the quantities appearing in the continuity equations to macroscopic quantities characterizing the electronic fluid (i.e., particle and energy densities and the flow velocity) and thus determine the equations of ideal hydrodynamics in graphene (i.e., the generalization of the usual Euler equation); (iii) using an approximate solution to the kinetic equation, I establish the leading dissipative corrections to the ideal hydrodynamics and establish the generalization of the main equation of the usual hydrodynamics, i.e., the Navier-Stokes equation. At the latter step I determine the explicit expressions for the dissipative coefficients, such as shear viscosity and quantum conductivity.

2 Ideal hydrodynamics in graphene

2.1 Continuity equations

Hydrodynamics is a direct consequence of conservation laws. These are commonly expressed in terms of continuity equations, which can be written on a phenomenological basis [14, 15]. Integrating the kinetic equation, one can not only “derive” the continuity equations, but also give explicit expressions for the corresponding densities and currents in terms of the distribution function.

2.1.1 Particle number conservation

The equation most commonly known as “the continuity equation” [14] expresses conservation of “mass” or “matter”: the amount of fluid in any given volume can be changed only by means of fluid flow through the volume boundary. In an electronic system, this is equivalent to conservation of electric charge. Within the kinetic theory, this conservation law is manifested in the vanishing of the integrated collision integral in (or the right-hand side of) the kinetic equation

[TABLE]

Here $\lambda=\pm$ is the band index, $\boldsymbol{k}$ is the momentum labeling single-particle states, and $N$ is the degeneracy factor (in real graphene $N=4$ due to spin and valley degeneracy).

The continuity equation can be obtained by integrating the kinetic equation (2a) and has the usual form [14, 16]. The only subtle point arising in two-band systems is the treatment of the formally infinite number of particles in the filled band. However, assuming the contribution of the filled band to be constant, one can immediately see that it vanishes upon differentiation and does not contribute to the continuity equation.

Consider first the time-derivative term in Eq. (2b). Integrating this term over all states yields

[TABLE]

defining the charge density, $n$, (up to the factor of electric charge). The definitions of the numbers of charge carriers in the two bands, $n_{\pm}$, are given in A.1.

Similarly, the gradient term can be integrated as

[TABLE]

defining the electric current, $\boldsymbol{j}$ (up to the factor of the electron charge; see A.2 for explicit definitions of the quasiparticle currents, $\boldsymbol{j}_{\pm}$ and $\boldsymbol{j}$).

Charge conservation requires that the electric field does not affect the continuity relation. Technically, this is expressed by means of the vanishing integral

[TABLE]

The situation with the magnetic field is more involved. Integrating the Lorentz term in Eq. (2a), one finds

[TABLE]

For any rotationally invariant spectrum, velocity and momentum have the same direction and the latter expression vanishes

[TABLE]

Systems with anisotropic spectra should be considered separately. Such analysis is beyond the scope of this paper. Whatever the spectrum, the Lorentz force cannot violate charge conservation.

Combining the above contributions, I find the standard continuity equation (usually, the continuity equation is expressed in terms of the charge density and electric current, which differ from the quantities $n$ and $\boldsymbol{j}$ by a multiplicative factor of the electric charge)

[TABLE]

Multiplying the kinetic equation by $\lambda$ and integrating over all states, one can find a similar equation for the imbalance current, $\boldsymbol{j}_{I}$, [see Eq. (81d)]

[TABLE]

where the right-hand side describes the recombination processes [22] (within the $\tau$-approximation; $n_{I}$ is the imbalance density, see A.1, and $n_{I}^{(0)}$ is the equilibrium imbalance density). Technically this term appears from the integration of the collision integral ${\rm St}_{R}[f]$, which does not conserve the number of particles in each band individually. In monolayer graphene, the dominant process contributing to quasiparticle recombination is the impurity-assisted electron-phonon scattering [43, 25].

2.1.2 Energy conservation

In two-band systems with unbound (from below) spectrum, one has to define the energy density relative to the (formally infinite) energy of the filled valence band, see A.3.

Similarly to the particle number conservation, energy conservation leads to the vanishing integral

[TABLE]

Multiplying the kinetic equation (2a) by the energy and integrating over all single-particle states leads to the continuity equation for the energy density [14, 16].

Since the quasiparticle energies and the energy of the filled valence band (82b) are independent of time, integrating the first term in the Liouville’s operator yields

[TABLE]

where the energy density, $n_{E}$, is defined in Eq. (82).

The integrated gradient term in the Liouville’s operator defines the energy current, $\boldsymbol{j}_{E}$ [see also Eq. (83)]

[TABLE]

The electric field acting on an electronic system leads to Joule’s heating. Integrating the electric field term in the kinetic equation one finds

[TABLE]

Finally, the Lorentz force cannot lead to any change of energy since it does not do any work. Indeed, integrating the Lorentz term in Eq. (2b), one finds (for a rotationally invariant system)

[TABLE]

Here the energy $\epsilon_{\lambda\boldsymbol{k}}$ is not differentiated since this would just yield a velocity and a cross product of two velocities is zero. The last term is similar to that in the continuity equation.

Combining all of the above, I find the continuity equation for the energy density

[TABLE]

where the right-hand side describes the Joule’s heat.

2.1.3 Momentum conservation

Multiplying the kinetic equation by the momentum and integrating over all states, one finds that the collision integrals due to electron-electron interaction and quasiparticle recombination vanish

[TABLE]

Impurity scattering, however, may relax momentum so that (unlike in the two previous cases) the impurity contribution to the collision integral does not vanish. Within the $\tau$-approximation [see Eq. (2a)], one finds

[TABLE]

where the momentum density, $\boldsymbol{n}_{\boldsymbol{k}}$, is defined in A.4.

The time derivative term in Eq. (2b) is treated same as before

[TABLE]

Integrating the gradient term, one finds the momentum flux tensor (85)

[TABLE]

The external forces can also change the momentum density. The electric field term can be integrated as follows

[TABLE]

The Lorentz force can change the direction of momentum. Integrating the Lorentz term similarly the above I find

[TABLE]

where at the last step I relied on rotational invariance.

Finally, I find the following equation

[TABLE]

which differs from the usual continuity equation for the momentum density [16] by the presence of the electromagnetic fields and the weak disorder scattering term. The latter has to be small [see Eq. (1)], otherwise the discussion of hydrodynamics makes no sense. However, the arguments leading to Eq. (6) are rather general: they do not rely on Eq. (1) and are justified in the whole applicability region of the kinetic equation (2).

2.1.4 Summary of the continuity equations

To summarize this section, I list all four continuity equations for a two-band electronic system (the above arguments are not specific to graphene)

[TABLE]

2.2 Local equilibrium

The main underlying assumption of the hydrodynamic theory is that of the local equilibrium established by the electron-electron collisions [16] on length scales much larger than $\ell_{ee}$. The corresponding distribution function can be written as [25]

[TABLE]

where ${\mu_{\lambda}(\boldsymbol{r})}$ is the local chemical potential and $\boldsymbol{u}(\boldsymbol{r})$ is the hydrodynamic (or “drift”) velocity.

In order to use the distribution function (8) in practical calculations, one needs to specify the quasiparticle spectrum (the above general form of the continuity equations is valid for an arbitrary two-band electronic system). In this paper, I will employ the following notations for the Dirac spectrum (the chirality $\lambda=\pm$ distinguishes the conduction and valence bands)

[TABLE]

and velocities (hereafter $\boldsymbol{e}_{\boldsymbol{a}}$ denotes a unit vector in the direction $\boldsymbol{a}$)

[TABLE]

Using the distribution function (8) and the Dirac spectrum (9) in the definitions (see A) of the hydrodynamic quantities featuring in the continuity equations (7), one can arrive at the equilibrium expressions for the macroscopic densities and currents. In particular, the quasiparticle currents (81c) and (81d) can be expressed in terms of $\boldsymbol{u}$ and the corresponding densities, as one might expect

[TABLE]

while the “charge” and “imbalance” densities defined in Eqs. (80) are given by

[TABLE]

The energy current (83) can be expressed in terms of the energy density (82)

[TABLE]

where the energy density is given by

[TABLE]

The momentum density $\boldsymbol{n}_{\boldsymbol{k}}$ is proportional to the energy current, see Eq. (84). The momentum flux tensor (85) is also expressed in terms of the energy density (14)

[TABLE]

Finally, the local equilibrium distribution function (8) can be used to compute the thermodynamic quantities. Introducing the linear combinations of the two chemical potentials $\mu_{\pm}$ (i.e. defining the thermodynamic variables conjugated to $n$ and $n_{I}$),

[TABLE]

I can define the compressibilities

[TABLE]

The hydrodynamic pressure (90) is proportional to the energy density

[TABLE]

and hence for the enthalpy one finds

[TABLE]

As a result, the expressions for the energy current (13) and momentum flux tensor (15) can be re-written as

[TABLE]

2.3 Generalized Euler equation

Substituting the above hydrodynamic quantities into the continuity equations (7), I obtain the equations of the ideal hydrodynamics in graphene. Consider first the equation (7d) representing momentum conservation. Using Eqs. (20) and (84), I find for the two derivative terms in Eq. (7d):

[TABLE]

[TABLE]

The last term can be found from Eq. (7c) where one can use the energy density from Eq. (19) and the energy current from Eq. (13)

[TABLE]

Combining the above three equalities, I find

[TABLE]

As a result, I find the generalization of the Euler equation for graphene taking into account weak disorder and external electromagnetic fields (cf. the standard Euler equation [14] comprising the first three terms on the left-hand side)

[TABLE]

Combined with the continuity equations (7), Eq. (21) describes the ideal electron-hole fluid in graphene. Hydrodynamics of an ideal fluid is well-studied [14, 16]. In what follows, I relate Eq. (21) to the standard Euler equation and its relativistic counterpart as well as consider the simplest solutions of the ideal hydrodynamics.

2.4 Entropy flow

Within the usual hydrodynamics, the ideal flow is isentropic [14] and hence one can derive a “continuity equation” for entropy, which provides a definition of the entropy current. The entropy density of a system of fermions is defined in terms of the distribution function as

[TABLE]

Treating this integral as

[TABLE]

I can represent any derivative of $s$ in the form

[TABLE]

Multiplying the kinetic equation in the absence of external fields by the derivative $\partial{\cal S}[f_{\lambda\boldsymbol{k}}]/\partial f_{\lambda\boldsymbol{k}}$ and summing over all states I obtain a macroscopic equation

[TABLE]

On the left-hand side of Eq. (23), $s$ is the entropy density (94) and the vector $\boldsymbol{j}_{S}$ can be interpreted as the entropy current

[TABLE]

Since the function $\cal S$ vanishes as $k\rightarrow\infty$, one may integrate by parts:

[TABLE]

Using the explicit form of the derivative, this expression can be re-written as

[TABLE]

which yields upon the substitution of the definitions (81), (82), (83a), and (85)

[TABLE]

Finally, using the explicit expressions for the macroscopic quantities in Eq. (25), I arrive at the result

[TABLE]

where $s$ is the entropy density (94).

Now, the right-hand side of Eq. (23) is the integrated collision integral

[TABLE]

Evaluating the derivative of ${\cal S}$ for the local equilibrium distribution function (8) explicitly, I find

[TABLE]

The first term does not contribute to the integral due to energy conservation. In the simplest case, ${\mu_{\pm}=\mu}$, the second term also vanishes due to charge conservation. However, if the quasiparticle imbalance becomes important, i.e., in the presence of a temperature gradient [22], the second term yields a non-vanishing contribution due to electron-hole recombination (since the collision integral ${\rm St}_{R}[f]$ does not conserve the number of particles in each individual band. This situation is outside of the scope of this paper and will be considered elsewhere. Finally, the last term yields the following contribution due to disorder scattering (which does not conserve momentum)

[TABLE]

Combining the above arguments, I obtain the final form of the continuity equation for the entropy density in graphene in the absence of temperature gradients (i.e., for ${\mu_{\pm}=\mu}$)

[TABLE]

The ideal hydrodynamic flow in the electronic system differs from its usual counterpart by the presence of weak disorder scattering that is all but unavoidable in any real solid.

2.5 Relativistic hydrodynamics

The linear spectrum of the elementary excitations in graphene, Eq. (9), suggests the possibility to observe relativistic hydrodynamics in a solid-state laboratory experiments. In this section, I compare the above ideal hydrodynamics in graphene to the relativistic hydrodynamics in 2D [44].

2.5.1 Hydrodynamic quantities in the relativistic theory

The above ideal (Euler) hydrodynamics can be compared to the standard relativistic hydrodynamics [14, 44] in $2+1$ dimensions, with $v_{g}$ playing the role of the speed of light. The central quantity in this theory is the relativistic stress-energy tensor

[TABLE]

The energy flux density is proportional to the momentum density and is given by $v_{g}T^{0\alpha}$.

Resolving Eq. (19) for the energy density $n_{E}$, I find

[TABLE]

The momentum density, $\boldsymbol{n}_{\boldsymbol{k}}$, can be expressed in the form similar to Eq. (28c) by combining Eqs. (20), (84), and (19):

[TABLE]

Again, the difference between the two expressions is that between $p$ and $P$.

Finally, the momentum flux (15) can be re-written similarly to the relativistic expression (28d) as well

[TABLE]

As a result, all three expressions (29) are similar to their relativistic counterparts (28), but are determined by the pressure $P$ defined in the “laboratory frame”, see Eq. (90), instead of the “relativistic” pressure $p$ defined in the local rest frame. On one hand, the reason $p$ is typically defined in the local rest frame is simply the lack of a better definition: the standard argument [14] relies on the Pascal law that is valid only in the local rest frame. On the other hand, the difference between the two theories is significant: expressions (29) cannot be obtained from their rest-frame counterparts by the Lorentz transformation. The reason for this is that the local equilibrium distribution function (8) is not relativistic. This can be traced to the classical (i.e. non-relativistic) nature of the Coulomb interaction in graphene that is ultimately responsible for equilibration.

2.5.2 Relativistic Euler equation

Let me now compare the generalized Euler equation (21) to the standard equations of relativistic hydrodynamics. These are encoded in the relation [14]

[TABLE]

Substituting the stress-energy tensor of the ideal fluid (28), one arrives at the equation [14]

[TABLE]

The relativistic generalization of the Euler equation is typically obtained [14] by projecting Eq. (31) on to the direction orthogonal to $u^{i}$. This is achieved by considering the combination

[TABLE]

which vanishes upon multiplication by $u^{i}$. Using Eq. (31) and the standard properties of the relativistic $3$-velocity (in 2D)

[TABLE]

one finds [14]

[TABLE]

The Euler equation is contained in the spatial components of Eq. (34)

[TABLE]

Combining Eqs. (35) as

[TABLE]

one finds the relativistic version of the Euler equation:

[TABLE]

Similarly to the hydrodynamic quantities in Eqs. (28), the relativistic equation (36) differs from the generalized Euler equation in graphene (21) in the absence of the electromagnetic fields and disorder scattering insofar it contains the rest frame pressure $p$ instead of the hydrodynamic pressure $P$.

Taking into account the electromagnetic fields, one replaces Eq. (30) by

[TABLE]

one can recover Eq. (7c) and the field dependent terms in Eqs. (7d) and (21).

Consider, however, Maxwell’s equations

[TABLE]

Using the above expressions, one can directly verify, that the two Maxwell’s equations containing only electric and only magnetic fields are preserved:

[TABLE]

However, the two remaining equations coupling the electric and magnetic fields are clearly violated. Therefore, one should use the relativistic approach with care.

2.5.3 Thermodynamic quantities and entropy

Combining the equation of motion (30) with the relativistic continuity equation

[TABLE]

one can derive the relativistic analog of the continuity equation for the entropy density (27). Indeed, projecting Eq. (30) onto the direction of the $3$-velocity,

[TABLE]

and using the thermodynamic (Gibbs-Duhem) relation

[TABLE]

where $s$ is the entropy density (in the rest frame), one finds

[TABLE]

The first term vanishes due to Eq. (38) and thus

[TABLE]

Here

[TABLE]

Thus the relativistic continuity relation has a similar form to the continuity equation for the entropy in graphene (with the entropy density $s$ and current $\boldsymbol{j}_{S}$ combined into a $3$-current $su^{k}$) in the absence of disorder scattering.

Similarly to all above cases of such correspondence, the entropy density in the relativistic theory is defined in the local rest frame, unlike the entropy density in graphene (19) which is defined with the local equilibrium distribution of a moving electronic fluid in the laboratory frame.

3 Dissipative corrections to electronic hydrodynamics

The above derivation of the generalized Euler equation (21) relies on the assumption of local equilibrium that is supposed to be established by electron-electron collisions. The same scattering processes are responsible for dissipation, i.e. irreversible charge and momentum transfer from faster elements of the electronic fluid to the slower ones. The general form of the dissipative corrections follows from general arguments [14]. Hydrodynamic quantities are supposed to vary slowly over long distances, such that their gradients should be small. Consequently, the dissipative correction to the momentum flux tensor $\Pi_{E}^{\alpha\beta}$ should be linear in the gradients of velocity. The specific form of the correction is governed by rotational invariance and in 2D in the absence of external magnetic field is given by

[TABLE]

The shear and bulk viscosity coefficients ($\eta$ and $\zeta$, respectively) can be found by a solution of the kinetic equation [16]. In the usual case of, e.g., a dilute gas, one can solve the Boltzmann equation by means of the perturbative Chapman-Enskog method [45, 46, 47, 48, 49]. Same results, although with a less clear justification, may be obtained using the Grad method [50].

Following the standard derivation of hydrodynamic equations from the kinetic theory [16, 25, 44, 13], I will look for the dissipative corrections to the ideal hydrodynamics (21) within linear response. The “linear response” solution should be obtained by linearizing the collision integral in Eq. (2a) in the small deviations from local equilibrium, ${\delta{f}=f-f^{(0)}}$, while leaving only $f^{(0)}$ on the left-hand side:

[TABLE]

Note that by definition, ${\rm St}_{ee}[f^{(0)}]=0$.

Within the classic approach [16], one evaluates the Liouville’s operator on the left-hand side of the above expression explicitly (with ${\boldsymbol{u}=0}$), then uses the ideal Euler equation and thermodynamic relations to express the result in the form explicitly containing the dissipative terms (again, with ${\boldsymbol{u}=0}$) as “external forces”. The goal of such calculation is to find the coefficients describing the dissipative corrections, i.e. viscosity and thermal conductivity. This approach hinges on the fact that the general form of the dissipative terms is known from symmetry arguments (up to the coefficients).

In pure (disorder-free) graphene, the conserved current is the energy current (since it is proportional to the momentum density), hence the dissipative coefficients include viscosity and electrical conductivity. My goal here is not only to determine these coefficients, but also to establish the form of the dissipative corrections. Therefore, instead of the direct evaluation of the Liouville’s operator, I will integrate the kinetic equation (41) following [25].

For Dirac fermions in graphene the solution of the kinetic equation is simplified by the kinematic peculiarity of electron-electron scattering known as the “collinear scattering singularity”. For Dirac quasiparticles moving along the same direction the energy and momentum conservation laws coincide leading to a formal divergence of the collision integral. Although the divergence is regularized by dynamical screening, the resulting scale separation allows for a nonperturbative solution.

3.1 Collision integral due to electron-electron interaction

The local equilibrium distribution function (8) nullifies the collision integral. Assuming that the external fields and other perturbations lead to “small” deviations from local equilibrium, the collision integral, ${\rm St}_{ee}[f]$, can be linearized in the small, non-equilibrium correction to $f_{\lambda\boldsymbol{k}}^{(0)}$ [16]

[TABLE]

The linearized collision integral can be written as [16] (the summation runs over all single-particle states up to the degeneracy factor which is written down explicitly)

[TABLE]

where $U$ stands for the dynamically screened Coulomb interaction.

3.2 Nonequilibrium correction to the distribution function

The two $\delta$-functions in Eq. (43b) represent energy and momentum conservation in an electron-electron “collision”. For Dirac fermions moving in the same direction they are identical and ${\rm St}_{ee}[f]$ diverges for a generic $h_{\lambda\boldsymbol{k}}$. There are however three exceptions,

[TABLE]

In the first case, the collision integral vanishes due to momentum conservation, while in the other two the collision integral vanishes for collinear particles. As a result, one can limit the mode expansion of the nonequilibrium correction to the distribution function, $h$, to the above three modes.

Adopting the “three-mode approximation”, I can write the correction $h$ in the form [25]

[TABLE]

The coefficients $\boldsymbol{h}^{(i)}$ and $h_{\alpha\beta}^{(i)}$ in Eq. (44a) satisfy general constraints [16] based on the fact that electron-electron collisions do not alter conserved thermodynamic quantities. To maintain momentum conservation, I should set

[TABLE]

In the presence of weak disorder momentum is no longer conserved and the energy current also acquires a dissipative correction. In this case, I have to keep a nonzero $\boldsymbol{h}^{(3)}$ and then study the (nontrivial) limit ${\tau_{\rm dis}\rightarrow\infty}$.

Now, to maintain conservation of the number of particles and energy I set

[TABLE]

The remaining coefficients can be determined by an explicit evaluation of the corresponding macroscopic quantities [25].

The macroscopic currents associated with the three modes $\phi_{i}$ are the electric, imbalance, and energy currents. Using the nonequilibrium distribution function in the definitions, Eqs. (81c), (81d), and (83a), I define the dissipative corrections (see also B)

[TABLE]

Substituting Eq. (44a) and evaluating the integrals (the tensor terms do not contribute for $\boldsymbol{u}\rightarrow 0$), I find the following relation between the three corrections (45) and the vector coefficients $\boldsymbol{h}^{(i)}$

[TABLE]

where the matrix elements of $\widehat{M}_{h}$ are expressed in terms of the equilibrium densities (12), (14) and compressibilities (17).

The tensor coefficients, $h_{\alpha\beta}^{(i)}$, in the second term of the nonequilibrium correction (44a) are similarly related to the dissipative corrections to the three macroscopic tensor quantities (85), (86), and (87)

[TABLE]

Substituting the distribution function (44a), one finds essentially the same integrals as in the case of the currents (now only the tensor part of the nonequilibrium correction yields a nonzero contribution) such that

[TABLE]

3.3 Electrical conductivity

The relation between the coefficients in the nonequilibrium distribution function (44a) and the macroscopic currents, Eq. (46), suggests the following method of solving the linearized kinetic equation (41). Integrating the kinetic equation, one obtains equations for the currents. Then using the relation (46) one finds dissipative corrections to the currents as linear functions of external fields.

3.3.1 Macroscopic equation for the electric current

The equation for the electric current is obtained by multiplying the kinetic equation (2a) by the velocity and integrating over all single-particle states [cf. Eq. (81c)]. The resulting equation will have the form

[TABLE]

where $\boldsymbol{\cal I}_{1}\left[f\right]$ is the integrated collision integral and the Lorentz term is not treated within linear response. Evaluating the integrals on the left-hand side of Eq. (49), one finds

[TABLE]

where $\boldsymbol{e}_{B}$ is the unit vector in the direction of $\boldsymbol{B}$, the tensors $\Pi^{\alpha\beta}$ and $\Pi_{\lambda}^{\alpha\beta}$ are defined in Eq. (86), the “band currents”, $\boldsymbol{j}_{\lambda}$ are defined in Eqs. (81), $\omega_{B}$ is the generalized cyclotron frequency, ${\cal T}$ is defined in Eq. (17e), and the vector quantity $\boldsymbol{\cal K}$ defined as

[TABLE]

has dimensions of the current.

Evaluating the integrated collision integral and the vector ${\cal K}$ using the nonequilibrium distribution function (44a) yields a linear function of the coefficients $\boldsymbol{h}^{(i)}$. Details of the calculation are relegated to C.1 and D.1. The result is summarized below together with the equations for the two other macroscopic currents.

3.3.2 Macroscopic equation for the imbalance current

The equation for the imbalance current is obtained similarly to Eq. (57): one multiplies the kinetic equation by $\lambda\boldsymbol{v}_{\lambda\boldsymbol{k}}$ and integrates over all single-particle states [see the definition (81d)]

[TABLE]

where $\boldsymbol{\cal I}_{2}\left[f\right]$ is the integrated collision integral. Evaluating the integrals, I find an equation similar to Eq. (50)

[TABLE]

where

[TABLE]

The integrated collision integral $\boldsymbol{\cal I}_{2}$ and the vector $\boldsymbol{\cal K}_{I}$ are calculated in C.1 and D.1, respectively.

3.3.3 Macroscopic equation for the energy current

The equation for the energy current is given by Eq. (7d), multiplied by $v_{g}^{2}$. Substituting the ideal quantities into the left-hand side (with the exception of the Lorentz term), one finds

[TABLE]

This equation differs from Eq. (7d) and hence from the Euler equation (21) by the fact that the electric current in the Lorentz term and the energy current on the right-hand side are total currents including the dissipative corrections. Adopting the standard iterative method [16] of the derivation of the dissipative corrections, one has to separate the terms in Eq. (55) forming the Euler equation and the terms containing $\delta\boldsymbol{j}$ and $\delta\boldsymbol{j}_{E}$. Assuming the validity of the Euler equation as the zeroth iteration, this leaves one with the following relation between the dissipative corrections

[TABLE]

The right-hand side yields the explicit form of the integrated collision integral (due to disorder only since the electron-electron interaction conserves momentum).

The simple form of Eq. (56) has a simple physical meaning. Since the energy current is proportional to the momentum density, it cannot be relaxed by electron-electron interaction (which conserves momentum). Consequently, a steady state cannot be achieved without disorder scattering contradicting the use of time-independent corrections to macroscopic currents (45).

3.3.4 Dissipative corrections to quasiparticle currents

The three equations (50), (53), and (55) coincide with the macroscopic linear response equations derived in Ref. [24]. These equations are (at least, in principle) sufficient for describing traditional linear response transport in graphene and are valid even for relatively strong disorder, where the hydrodynamic approach is invalid. At the same time, for weak disorder, i.e., within the applicability region of the hydrodynamic theory, the latter provides a significant generalization of the linear response theory allowing for a description of the collective motion of a strongly interacting fluid.

Applying the above iterative approach of the derivation of the dissipative corrections to the ideal Euler hydrodynamics, I now simplify Eqs. (50) and (53) assuming the validity of the Euler equation for the ideal quantities in their respective left-hand sides. In particular, the Euler equation (21) can be used to express the time derivative of the velocity in terms of the pressure gradient and electromagnetic fields. The pressure gradient can be expressed in terms of the gradient of $n_{E}$ using the equation of state (18). Finally, density gradients can be expressed in terms of gradients of temperature and chemical potential. As a result, I arrive at the equations

[TABLE]

[TABLE]

Combining the integrated equations for macroscopic currents (57), (58), and (56) with the integrated collision integrals (see C.1) and Lorentz terms (see D.1), I obtain the final set of linear equations for the dissipative corrections to macroscopic currents (45)

[TABLE]

where $t_{\rm d}$ is the dimensionless impurity scattering time defined similarly to the way the matrix $\widehat{\textswab{T}}$ is defined.

In general,the ${6\!\times\!6}$ matrix on the right-hand side of Eq. (59e) may be inverted as follows. Introducing the short-hand notations,

[TABLE]

I can re-write Eq. (59e) as follows

[TABLE]

Multiplying this equation by $\boldsymbol{e}_{B}$ (using the fact that this vector product acts in the position space and hence commutes with all the matrices), I obtain

[TABLE]

The two equations can now be solved as a usual system of two linear equations with the only difference, that the coefficients are now matrices that do not commute. Hence, one has to keep track of the order in which they are multiplied. The resulting solution has the form

[TABLE]

which reminds one of the standard form of magnetoconductivity in the Drude theory

[TABLE]

The result (61) expresses the dissipative corrections to the macroscopic currents in the system as a function of the electric field (more precisely, of the gradient of the electrochemical potential) defining the dissipative coefficients in analogy with the thermal conductivity in the traditional hydrodynamics [14].

3.4 Viscosity

Within the usual hydrodynamics [14], shear and bulk viscosities are defined as the coefficients in the leading term in the gradient expansion of the dissipative correction to the momentum flux tensor, see Eq. (40). In this section, I establish the form of this correction in graphene following the same steps leading to the corrections to the quasiparticle currents, Eq. (61). Because I am now looking for corrections to a tensor quantity, the second term in the nonequilibrium distribution function (44a) is going to contribute.

3.4.1 Macroscopic equations for tensor quantities

Although the viscosity is defined though the dissipative correction to only one macroscopic tensor quantity, the momentum flux tensor $\Pi_{E}^{\alpha\beta}$, the three-mode approximation adopted in this paper requires one to consider equations determining the three macroscopic tensors: $\Pi^{\alpha\beta}$, $\Pi_{I}^{\alpha\beta}$, and $\Pi_{E}^{\alpha\beta}$. Similarly to the equations for the corresponding currents, these equations can be obtained by multiplying the kinetic equation by $v^{\alpha}v^{\beta}$, $\lambda v^{\alpha}v^{\beta}$, and $\epsilon v^{\alpha}v^{\beta}/T$ (respectively) and integrating over all states. The direct integration yields the three equations (where the external electric field is set to zero since I am looking for viscosity as a function of magnetic field only).

[TABLE]

The third-rank tensors $\Upsilon^{\alpha\beta\gamma}$, $\Upsilon^{\alpha\beta\gamma}_{I}$, and $\Upsilon^{\alpha\beta\gamma}_{E}$ appear in the integrated equations (62) in the same way as the second-rank tensors $\Pi^{\alpha\beta}$, $\Pi_{I}^{\alpha\beta}$, and $\Pi_{E}^{\alpha\beta}$ appear in the integrated equations for the macroscopic currents. Similarly to the evaluation of the left-hand sides of the integrated equations for macroscopic currents, these quantities have to be computed with the local equilibrium distribution function (8). In the limit ${\boldsymbol{u}\rightarrow 0}$ these quantities are linear in $\boldsymbol{u}$ and would be discarded if it were not for the fact that Eqs. (62) contain only gradients of these quantities. The straightforward calculation yields

[TABLE]

Now, in the limit ${\boldsymbol{u}\rightarrow 0}$ the second-rank tensors are proportional to the corresponding densities, see Eq. (15). Using the continuity equations to express the time derivatives of densities in terms of gradients similarly to the transformations used to derive Eqs. (57) and (58), one can simplify the left-hand sides of Eqs. (62). The result can be expressed in vector form as

[TABLE]

Comparing the left-hand side of Eq. (64) to the definition (40), I can already conclude that the bulk viscosity in graphene vanishes (at least within the approximations adopted in this paper).

The integrated collision integrals ${\cal I}_{j}^{\alpha\beta}$ are discussed in C.2. The integrated Lorentz terms contain the tensors $\Xi^{i\alpha}$, $\Xi_{I}^{i\alpha}$, and $\Xi_{E}^{i\alpha}$ in analogy with the vectors ${\cal K}$ and ${\cal K}_{I}$ in the above conductivity calculation. The discussion of these tensors is relegated to D.2. Combining them in a vector in the “mode space”, I find

[TABLE]

where the coefficients in Eq. (120) are combined into the matrix (59c). Note, that according to the definitions (119) the matrices $\Xi^{i\beta}_{n}$ are symmetric, but not necessarily traceless. However, the matrices $h^{(n)}_{i\beta}$ are traceless, see Eq. (44d), hence the matrices $\Xi^{i\beta}_{n}$ are traceless as well.

Traceless, symmetric, $2\times 2$ matrices contain only two independent elements. Consequently, the last two terms in Eq. (64) must be related to each other [in other words, Eq. (64) as a system of two linear equations for the two matrix elements of $\Xi^{i\beta}_{n}$]. Indeed, evaluating the spatial components explicitly, I find

[TABLE]

which are identical, since

[TABLE]

Using the explicit form of the collision integrals (112) and the quantities $\Xi^{i\beta}$, Eq. (65), and taking into account the above argument, I can express the vector on the right-hand side of Eq. (64) as

[TABLE]

where $\widehat{\textswab{T}}_{\Pi}$ is the matrix of electron-electron scattering rates, see Eq. (112), in the dimensionless form (59b).

Now I can solve Eq. (64) similarly to the solution of Eq. (59). Introducing the notation [cf. Eq. (60)],

[TABLE]

I may re-write Eq. (64) as

[TABLE]

Multiplying this equation by $\boldsymbol{e}_{B}$, I obtain [similarly to the calculation below Eq. (60)]

[TABLE]

Similarly to Eq. (61), I find the solution in the form (I am only interested in $\delta\Pi^{\alpha\beta}_{E}$)

[TABLE]

I obtain the final expressions for the shear and Hall viscosities [51, 52] in graphene

[TABLE]

The sign of the shear viscosity $\eta$ is fixed by the laws of thermodynamics [14, 16]. In contrast, Hall viscosity is non-dissipative (since the Lorentz force does not perform any work) and may have an arbitrary sign which is technically determined by the quasiparticle charge and direction of the magnetic field. In this paper, I choose the $\eta_{H}$ to be positive for electrons by analogy with Hall conductivity [4].

3.5 Generalized Navier-Stokes equation

Substituting the dissipative correction $\delta\Pi_{E}^{\alpha\beta}$ into the continuity equation for momentum density (6) and repeating the steps used to derive the Euler equation (21) I find the generalization of the central equation of the traditional hydrodynamics, the Navier-Stokes equation [14] to the electronic system in graphene

[TABLE]

Combined with the expressions for viscosities (68) and dissipative corrections to currents (61) this equation represents the central result if this paper. Previously, the generalized Navier-Stokes equation in graphene was derived in Ref. [25] in the absence of disorder and the magnetic field and in Ref. [23] in the absence of disorder and the external fields (both electric and magnetic).

4 Discussion

The purpose of this paper was to derive the hydrodynamic equations for the electronic fluid in graphene in the presence of electromagnetic fields and weak disorder as well as to obtain closed expressions for shear and Hall viscosities and electrical conductivity in graphene (the latter being an analog of the thermal conductivity in the traditional hydrodynamics). Despite the conceptual simplicity of the assumptions leading to the hydrodynamic description, the Navier-Stokes equation is known to yield a large number of important solutions, see Ref. [14]. It is therefore impractical to include even a small subset of these solutions into a single paper. Nevertheless it is important to show that the cumbersome expressions for the dissipative corrections to quasiparticle currents (61) and the shear viscosity (68a) yield the well-known results in the simplest limiting cases.

4.1 Quantum conductivity

Consider first the electrical conductivity at charge neutrality in the absence of the magnetic field, known as the “quantum” (or “intrinsic”) conductivity [13, 23]. In this case the equation (61) simplifies. Setting the chemical potential to zero, one also finds the vanishing charge density, ${\tilde{n}(0)=0}$, and the “imbalance compressibility”, ${\partial{n}_{I}/\partial\mu(0)=0}$. Hence the correction to the electric current, $\delta\boldsymbol{j}$, represents the “whole” current, since the “ideal” part of the current vanishes, see Eq. (11). At the same time, the Navier-Stokes equation (69) [together with Eq. (56)] admits a stationary and uniform solution, ${\boldsymbol{u}=0}$. Therefore, for the stationary and uniform fields the energy (as discussed above) and imbalance currents vanish, see Ref. [24],

[TABLE]

The combination of the distribution functions (108b) vanishes as well, ${I(x=0)=0}$, such that ${\tau^{-1}_{12}(\mu\!=\!0)=0}$. The energy density and compressibility are determined by temperature, ${n_{E}=3N\zeta(3)T^{3}/(2\pi{v}_{g}^{2})}$, ${{\cal T}=2T\ln 2}$. Substituting these values into Eq. (61), I find

[TABLE]

where $t_{11}$ and $t_{22}$ are the diagonal elements of the matrix $\widehat{\textswab{T}}$, see Eq. (59b), and

[TABLE]

The solution (70) for the dissipative correction to the electric current yields the resistivity of undoped graphene [24]

[TABLE]

where $\sigma_{Q}$ is the “intrinsic” conductivity of disorder-free graphene at the Dirac point [31, 20, 21, 53, 54]. Previously, the coefficient ${\cal A}$ was reported to have values ${{\cal A}=0.12}$ [20] and ${{\cal{A}}=0.19}$ [25]. Evaluating the integral for $t_{11}(0)$ numerically for unscreened Coulomb interaction (as was done in Refs. [20, 25]), I obtain the value ${{\cal{A}}=0.113\pm 0.03}$, where the error comes from systematic differences between various numerical methods (note, that the present calculation neglects the exchange contribution $\sim 1/N$, which was shown to be numerically small in Ref. [20]).

4.2 Electrical conductivity in the degenerate regime

Away from charge neutrality, the electric current is no longer determined by the dissipative correction alone. The “ideal” contribution to the current is governed by the hydrodynamic velocity $\boldsymbol{u}$, see Eq. (11). The velocity $\boldsymbol{u}$ is a solution of the Navier-Stokes equation (69). Within linear response and in the absence of the magnetic field, the uniform and time-independent solution to Eq. (69) is simply the Ohm’s law with the Drude formula for conductivity

[TABLE]

In the absence of magnetic field, Eq. (56) yields the vanishing dissipative correction to the energy current, $\delta\boldsymbol{j}_{E}=0$. However, in the degenerate regime, ${\mu\gg{T}}$, all three dissipative corrections $\delta\boldsymbol{j}$, $\delta\boldsymbol{j}_{I}$, and $\delta\boldsymbol{j}_{E}$ are proportional to each other, since the three rows of the matrix $\widehat{M}_{h}$ [see Eq. (46)] are identical (in this limit ${n=n_{I}=\mu^{2}/(\pi{v}_{g}^{2})}$ since only one band contributes and ${n_{E}=2\mu^{3}/(3\pi{v}_{g}^{2})}$, with $N=4$). This means that all three dissipative corrections vanish and, in particular,

[TABLE]

The total electric current $\boldsymbol{J}=e\boldsymbol{j}$ is then determined by the “hydrodynamic” contribution alone

[TABLE]

Thus, the conductivity of graphene in degenerate regime is given by the Drude result (due to disorder), in agreement with the leading-order result of the linear response theory [24].

4.3 Magnetoconductivity at charge neutrality

Similarly to the discussion in Sec. 4.1, the equation for the dissipative corrections (59) simplifies at charge neutrality. Instead of using the general solution (61), it might be instructive to represent Eq. (59) as three vector equations and solve them directly. Indeed, using the parameter values listed in Sec. 4.1 together with the stationary and uniform solution of the Navier-Stokes equation, ${\boldsymbol{u}=0}$, one can represent the matrix equation (59) as a system of coupled equations

[TABLE]

The last equation yields the energy current

[TABLE]

The fact that the electric current in magnetic field is parallel to the electric field can be expressed in terms of vanishing Hall coefficient, physically due to the exact electron-hole symmetry,

[TABLE]

In this section I have assumed an infinite system. Hence all macroscopic quantities are homogeneous and the hydrodynamic equations become algebraic. In finite size systems the situation is more involved. Here one has to solve differential equations (with the appropriate boundary conditions [28]). The resulting flows may be highly inhomogeneous [56, 57, 58, 59, 60] leading to e.g. linear [24, 58, 61] or negative magnetoresisteance [51, 62, 63].

4.4 Shear viscosity at charge neutrality in zero field

The expression for the shear viscosity also simplifies at charge neutrality (the Hall viscosity vanishes at this point altogether). The vanishing combination of the distribution functions (108b), ${I(x=0)=0}$, leads to vanishing of the two scattering rates

[TABLE]

In the absence of the magnetic field, one can simplify the coefficients in Eq. (68a) using the parameter values listed in Sec. 4.1. This leads to the following expression for the shear viscosity

[TABLE]

and hence

[TABLE]

where $\tilde{t}_{ij}$ are the elements of the matrix $\widehat{\textswab{T}}_{\eta}$. This result was previously found in Ref. [23] with the numerical value of the coefficient ${{\cal{B}}=0.45}$ (evaluated for unscreened Coulomb interaction in the absence of disorder). Evaluating the dimensionless scattering rates $t_{ij}$ numerically, I find ${{\cal{B}}=0.446\pm 0.005}$ where similarly to the coefficient ${\cal A}$ in Eq. (71) the deviation is due to differences between various numerical methods. The exchange contribution is again neglected, but this does not seem to lead to any appreciable error.

Finally, I can use Eq. (94) to compute the ratio of the shear viscosity to the entropy density at charge neutrality. Indeed, for ${\mu_{\pm}=0}$ the entropy density is determined by pressure, ${s=P/T}$, which in turn is proportional to the energy density. As a result,

[TABLE]

Dividing the viscosity (75) by the entropy density (76) one finds

[TABLE]

This should be compared with the conjectured lower bound [64]

[TABLE]

The shear viscosity to entropy ratio was discussed in detail in [23], where it was shown that renormalization of the coupling constant leads to a logarithmic temperature dependence of the above ratio [formally, by replacing $\alpha^{2}$ in denominator by the renormalized value $\alpha^{2}(T)\approx 16/\ln^{2}(T_{\Lambda}/T)$, where $T_{\Lambda}$ is the cut-off scale]. Hence the ratio is expected to grow as one lowers the temperature. At high enough temperatures, the ratio decreases with the growing $T$, never quite reaching the lower bound [23].

5 Summary

In this paper I have presented a detailed derivation of the hydrodynamic theory of electronic transport in graphene in the presence of the external magnetic field and weak disorder. The main results of the paper are the generalized Navier-Stokes equation (69), explicit expressions for the shear and Hall viscosity (68) and for the dissipative corrections to the quasiparticle currents (61). These results agree with the previously reported values of the quantum conductivity [20, 21, 53, 54] and shear viscosity [23] at charge neutrality in pure graphene in zero field, providing an extension to arbitrary doping levels, (non-quantizing) magnetic fields, and weak disorder. For a detailed discussion of the quantitative results for the shear and Hall viscosities in a wide range of temperatures and carrier densities and their relation to the experimental data of Refs. [4, 11] see Ref. [65]. For a similar discussion of the optical conductivity see Ref. [66].

The hydrodynamic theory derived in this paper is justified by the inequality (1) and is valid for classical (non-quantizing) magnetic fields. The former condition restricts the temperature range of the hydrodynamic effects to relatively high temperatures above $100$K (for a discussion of the experimental viewpoint on that issue see Ref. [6]). The latter requirement was discussed in detail in Ref. [53] (see also the corresponding discussion in Ref. [66]). At high enough temperatures [required by Eq. (1)] Landau levels are smeared out such that one does not need to account for the quantum Hall effect for magnetic fields up to several Tesla. At stronger fields the hydrodynamic description breaks down. However, the kinetic theory has a wider applicability range. As a result, magnetotransport in stronger fields can be described by the standard linear response theory [53, 24, 66]. Such considerations are applicable as long as the typical cyclotron frequency does not exceed temperature [53]. Beyond that field range, quantum effects take over which cannot be described by the semiclassical kinetic theory. Transport in quantized magnetic fields is beyond the scope of the present paper.

Very recently I became aware of the related work on magnetotransport in graphene, Ref. [67], and on Hall viscosity, Refs. [68, 69].

Acknowledgments

The author wishes to thank I.V. Gornyi, A.D. Mirlin, J. Schmalian, M. Schütt, and A. Shnirman for fruitful discussions. This work was supported by the German Research Foundation DFG within FLAG-ERA Joint Transnational Call (Project GRANSPORT), by the European Commission under the EU Horizon 2020 MSCA-RISE-2019 program (Project 873028 HYDROTRONICS), and the MEPhI Academic Excellence Project, Contract No. 02.a03.21.0005.

Appendix A Macroscopic quantities

A.1 Quasiparticle densities

Quasiparticle (number) density can be defined in the usual form [16]. The only subtle point arising in two-band systems is the treatment of the formally infinite number of particles in the filled band. In other words, one needs to make a distinction between the “particle density” and “carrier density”. In the valence band, the latter is defined by the particle-hole transformation

[TABLE]

with the hole energies and chemical potential defined as

[TABLE]

The numbers of “carriers” (i.e. the low-energy excitations) in the two bands can be explicitly defined as

[TABLE]

A.2 Quasiparticle currents

Similarly to the quasiparticle densities (80), the macroscopic currents can be defined as

[TABLE]

A.3 Energy density

In two-band systems with unbound (from below) spectrum, one has to define the energy density relative to the (formally infinite) energy of the filled valence band

[TABLE]

where the subscript “e” stands for “electrons” and replaces the index $+$ in Eq. (82c), while the subscript “h” stands for “holes” as defined in Eq. (79).

A.4 Energy current

The calculation of the energy current can be performed along the same lines. In graphene, the energy current, $\boldsymbol{j}_{E}$, is defined as

[TABLE]

Alternatively, one may re-write the energy current (83a) using the electron-hole transformation (79) as

[TABLE]

Note, that in graphene the energy current (83) is proportional to the momentum density [due to Eq. (10)]

[TABLE]

A.5 Momentum flux tensor

Similar calculation can be performed for the momentum flux tensor (also known as the stress-energy tensor or the energy-momentum tensor)

[TABLE]

Formally, the expression (85) is divergent and (similarly to the quasiparticle and energy densities) has to be defined up to the formally in finite contribution of the filled band. However, in all hydrodynamic equations I will be dealing with derivatives of $\Pi_{E}^{\alpha\beta}$ which allows me to subtract this contribution. Note, that in a rotationally invariant system, the kinetic definition (85) is manifestly symmetric.

In addition two further tensor quantities can be formed (by analogy with the three macroscopic currents): the “velocity flux tensor”

[TABLE]

and the “imbalance flux tensor”

[TABLE]

cf. Eqs. (85) and (86). These quantities are not directly related to any traditional observables and play an auxiliary role in the kinetic theory.

A.6 Pressure and enthalpy

Under the assumption of local equilibrium, the expression for pressure can be obtained from the thermodynamic potential since [16]

[TABLE]

The thermodynamic potential is calculated in the usual grand canonical ensemble, with the only caveat that one should be careful while dealing with the nearly filled valence band. The easiest way is to use the electron-hole transformation (79). Then the thermodynamic potential of the two-band system described by the local equilibrium distribution (8) can be written as

[TABLE]

Differentiating $\Omega$ with respect to the chemical potential, one recovers the number of particles (80):

[TABLE]

Combining Eqs. (88) and (89), one finds the hydrodynamic pressure,

[TABLE]

The enthalpy can then be found using the standard relation

[TABLE]

A.7 Entropy

In thermodynamics, entropy is defined as a

[TABLE]

where $\Omega$ is the thermodynamic potential (89). A straightforward calculation leads to the result:

[TABLE]

This should be compared with the standard thermodynamic relation

[TABLE]

Given that the standard definition of entropy is formulated in thermal equilibrium, it is not surprising that the two relations coincide for $\boldsymbol{u}=0$ up to one important issue. Assuming the two independent chemical potentials in graphene, one has to generalize Eq. (93) replacing

[TABLE]

and adding the velocity term. Then the entropy in local equilibrium is given by

[TABLE]

Appendix B Nonequilibrium distribution function

Here I compare the dissipative correction to the local equilibrium distribution function in the hydrodynamic approach to the nonequilibrium distribution function considered within the linear response theory of Ref. [24].

Within linear response, one starts with the global equilibrium [described by the usual Fermi distribution function, $f_{0}(\epsilon_{\lambda\boldsymbol{k}})$] where no currents are flowing

[TABLE]

When the system is subjected to external fields, the distribution function acquires a nonequilibrium correction, $\delta f_{LR}$, yielding non-zero currents

[TABLE]

In the notation used in Ref. [24], ${\boldsymbol{P}\equiv\boldsymbol{j}_{I}}$, ${\boldsymbol{Q}\equiv\boldsymbol{j}_{E}}$.

The nonequilibrium correction, $\delta f_{LR}$, can be analyzed within the same three mode approximation (44). Comparing the notation used in the present paper to that of Ref. [24], one finds [see Eq. (42)]

[TABLE]

Substituting the latter expression into the definitions (95), one can express the quantities $\boldsymbol{\cal A}$, $\boldsymbol{\cal B}$, and $\boldsymbol{\cal C}$ in terms of the macroscopic currents $\boldsymbol{j}$, $\boldsymbol{j}_{I}$,a nd $\boldsymbol{j}_{E}$. In Ref. [24] this was done explicitly, but the result coincides with Eq. (46). The reason for this is that within the three-mode approximation the corrections (45) and the currents (95) are given by exact same integrals.

To better understand the relation between the two different approaches, one can expand the local equilibrium distribution function (8) to the linear order in the hydrodynamic velocity. As a result, the complete distribution function out of equilibrium takes the form

[TABLE]

Representing $\delta f_{\boldsymbol{u}}$ in the three-mode form (44)

[TABLE]

and using Eq. (46), I recover the “ideal” values for the three macroscopic currents, Eqs. (11) and (13), where all nonlinearities (e.g., in denominators) are neglected. Consequently, the sum of the two corrections ${\delta{f}_{\boldsymbol{u}}\!+\!\delta f}$ yields the total currents and hence is equivalent to the linear response correction $\delta f_{LR}$.

Appendix C Integrated collision integral

C.1 Collision integral in the current equations

Here I evaluate the collision integrals in the integrated kinetic equations. Given the additive nature of the collision integral in the kinetic equation (2a), I separate the momentum-conserving collision integral due to electron-electron interaction from disorder scattering and other momentum nonconserving processes

[TABLE]

where the collision integral due to electron-electron interaction is nullified by the local equilibrium distribution function, ${\boldsymbol{\cal{I}}^{ee}_{i}\left[f^{(0)}\right]=0}$, and all momentum non conserving processes are grouped together into ${\boldsymbol{\cal I}^{\rm dis}_{i}\left[f\right]}$. These will be considered within the simplest $\tau$-approximation.

C.1.1 Collision integral due to electron-electron interaction

The general form of the collision integral due to electron-electron interaction is given by Eq. (43). Introducing the transferred energy $\omega$ and momentum $\boldsymbol{q}$, one may write the transition probability (43b) as

[TABLE]

where

[TABLE]

The vertices $\lambda_{\boldsymbol{v},\boldsymbol{v}^{\prime}}$ are known as the “Dirac factors”. They indicate the asymmetry of quasi-particle scattering in graphene [18].

The integrated collision integrals $\boldsymbol{\cal I}_{i}^{ee}$ are obtained by multiplying the collision integral ${\rm St}_{ee}$ by

[TABLE]

following by integration. The resulting integrated collision integrals are given by (again, the degeneracy factors are written down explicitly)

[TABLE]

Substituting the nonequilibrium distribution function (44a), one finds after a straightforward but tedious calculation (the integral $\boldsymbol{\cal I}_{3}^{ee}=0$ is introduced for consistency of the notation)

[TABLE]

In the case ${\mu_{\pm}=\mu}$ considered in this paper, the inverse times $\tau^{-1}_{ij}$ are given by the following integrals:

[TABLE]

Hereafter I use dimensionless variables (the dimensionless frequency $W$ should not be confused with the enthalpy)

[TABLE]

The Coulomb interaction has the form

[TABLE]

where $\varepsilon$ is the effective dielectric constant describing the electrostatic environment and $\widetilde{U}$ accounts for screening effects.

The auxiliary functions $Y$ in Eq. (103) are given by

[TABLE]

where

[TABLE]

and

[TABLE]

C.1.2 Contribution of the collision integral due to disorder

In this paper, I am using the simplest $\tau$-approximation for the collision integral due to disorder scattering, see Eq. (2a). The corresponding integrated collision integral is given by

[TABLE]

Here I have used the form (42)-(44) of the non-equilibrium correction to the distribution function.

C.2 Collision integral in the tensor equations

C.2.1 Collision integral due to electron-electron interaction

Integrating the collision integral with the factors $v^{\alpha}v^{\beta}$, $\lambda v^{\alpha}v^{\beta}$, and $\epsilon v^{\alpha}v^{\beta}/T$ yields the following tensor quantities, see Eq. (62)

[TABLE]

Combining the multiplication factors using the “mode” notations, see Eq. (44a),

[TABLE]

and following the same steps as in C.1, I obtain the expression

[TABLE]

Evaluating the integrals (for ${\mu_{\pm}=\mu}$ and ${\boldsymbol{u}=0}$) I find a vector in the “mode space” that can be written in the form similar to Eq. (102) as

[TABLE]

with

[TABLE]

Here I have introduced auxiliary functions

[TABLE]

that involve the integrals (107) that need to be complemented by

[TABLE]

C.2.2 Contribution of the collision integral due to disorder

Using the simplest $\tau$-approximation for the collision integral due to disorder scattering, see Eq. (2a), I find the integrated collision integral in the tensor equations (62) as

[TABLE]

Here I have used the form (42)-(44) of the non-equilibrium correction to the distribution function and evaluated the collision integral in the co-moving frame, ${\boldsymbol{u}\rightarrow 0}$, where

[TABLE]

Appendix D Integrated Lorenz terms

D.1 Contribution to the current equations

To evaluate the contribution of the Lorentz force to the equations for quasiparticle currents, I need to calculate the quantities $\boldsymbol{\cal K}$ and $\boldsymbol{\cal K}_{I}$ defined in Eqs. (51) and (54). These expressions should be evaluated with the non-equilibrium correction to the distribution function (46).

In local equilibrium, both vectors are proportional to the hydrodynamic velocity. In the limit ${\boldsymbol{u}\rightarrow 0}$, relevant for the linear response derivation of the dissipative corrections to the ideal hydrodynamics, I find

[TABLE]

Substituting the non-equilibrium distribution correction (42)-(44), I find the correction to the quantity $\boldsymbol{\cal{K}}$

[TABLE]

Evaluating the angular integral as usual,

[TABLE]

I find

[TABLE]

The quantity $\delta\boldsymbol{\cal K}_{I}$ is calculated in a similar fashion. The result is given by

[TABLE]

[TABLE]

The results (118d) and (118f) coincide with the results of Ref. [24]. Using the notations of B these expressions can be re-written as

[TABLE]

[TABLE]

D.2 Contribution to the tensor equations

Integrating the Lorentz terms in the kinetic equation with the factors $v^{\alpha}v^{\beta}$, $\lambda v^{\alpha}v^{\beta}$, and $\epsilon v^{\alpha}v^{\beta}/T$ yields the following tensor quantities, see Eq. (62)

[TABLE]

Substituting the nonequilibrium distribution function (44a and evaluating the integrals, I find

[TABLE]

These results are summarized in the main text in Eq. (65).