Hydrodynamic approach to electronic transport in graphene: energy relaxation

Electronic hydrodynamics is quickly growing into a mature field of condensed matter physics [1, 2, 3]. Similarly to the usual hydrodynamics [4, 5], this approach offers a universal, long-wavelength description of collective flows in interacting many-electron systems. As a macroscopic theory of strongly interacting systems, hydrodynamics should appear to be extremely attractive for condensed matter theorists dealing with problems where strong correlations invalidate simple theoretical approaches. However, electrons in solids exist in the environment created by a crystal lattice and typically experience collisions with lattice imperfections (or “disorder”) and lattice vibrations (phonons). The former typically dominate electronic transport at low temperatures, while at high temperatures the electron-phonon interaction takes over. In both cases the electron motion is diffusive (unless the sample size is smaller than the mean free path in which case the motion is ballistic) since in both types of scattering the electronic momentum is not conserved. On the other hand, if a material would exist where the momentum-conserving electron-electron interaction would dominate 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. In recent years, several extremely pure materials became available with graphene being the most studied [1, 3].

As a manifestation of macroscopic conservation laws, hydrodynamics is universal. Most conventional fluids are assumed Galilean invariant and are described by the same set of hydrodynamic equations [4]. Similar approach was at the heart of the early theoretical work on electronic hydrodynamics [6, 7]. Another well-known case is the relativistic hydrodynamics [4] relevant to neutron stars and interstellar matter. Since its low-energy quasiparticles are characterized by the Dirac spectrum, graphene attracted significant theoretical attention as a possible condensed matter realization of relativistic hydrodynamics [8, 9, 10, 11, 12, 13]. However, due to the classical, three-dimensional nature of the Coulomb interaction between electrons, the emergent hydrodynamics in graphene is neither Galilean- nor Lorentz-invariant [2].

In nearly neutral (or compensated) graphene the electron system is non-degenerate (at least at relatively high temperatures where the hydrodynamic approach is justified) with both the conductance and valence bands contributing on equal footing. Although the electron system is not Lorenz-invariant, the linearity of the Dirac spectrum plays an important role. Firstly, the Auger processes are kinematically suppressed leading to the near-conservation of the number of particles in each band [2, 3, 14, 15]. Secondly, the so-called collinear scattering singularity [12, 10, 16, 11, 15, 17, 18, 19] allows for a non-perturbative solution to the kinetic (Boltzmann) equation focusing on the three hydrodynamic modes [18, 20, 21]. As a result, one can determine the general form of the hydrodynamic equations and to evaluate the kinetic coefficients [21, 22, 23]. To be of any practical value, the latter calculation has to be combined with the renormalization group approach [24] since the effective coupling constant in real graphene (either encapsulated or put on a dielectric substrate) is not too small, $\alpha_{g}\approx 0.2-0.3$ [25, 26].

Next to the conservation laws, the main assumption of the hydrodynamic approach is local equilibrium [4, 27] established by means of interparticle collisions. Neglecting all dissipative processes, this allows (together with the conservation laws) for a phenomenological derivation of hydrodynamic equations [4, 5] that can be further supported by the kinetic theory, where the local equilibrium distribution function nullifies the collision integral in the Boltzmann equation [27]. The resulting ideal hydrodynamics is described by the Euler equation and the continuity equations. This is where the electronic fluid in graphene differs from conventional fluids (both Galilean- and Lorentz-invariant): as in any solid, conservation laws in graphene are only approximate, leaving the collision integrals describing scattering processes other than the electron-electron interaction to be nonzero even in local equilibrium. This leads to the appearance of weak decay terms in the continuity equations.

Two such terms have already been discussed in literature. Firstly, even if the electron-electron interaction is the dominant scattering process in the system, no solid is absolutely pure. Consequently, even ultra-pure graphene samples possess some degree of weak disorder. Disorder scattering violates momentum conservation and hence a weak decay term must appear in the generalized Euler equation [2, 3, 20, 21]. Secondly, conservation of the number of particles in each band is violated by a number of processes (e.g., the Auger and three-particle scattering). While commonly assumed to be weak, they are manifested in the decay – or recombination – term in the corresponding continuity equation. This was first established in [14] in the context of thermoelectric phenomena (for the most recent discussion see [28]). Later, quasiparticle recombination was shown to lead to linear magnetoresistance in compensated semimetals [29, 30, 31, 32] as well as giant magnetodrag [33, 34].

In this paper, we report the derivation of the third weak decay term in the hydrodynamic theory in graphene due to weak violation of energy conservation. Indeed, the electron-phonon interaction may lead not only to the loss of electronic momentum (responsible for electrical resistivity in most metals at high temperatures), but also to the loss of energy. Although subdominant in the hydrodynamic regime, the electron-phonon interaction should be taken into account as one of the dissipative processes. In graphene, the linearity of the Dirac spectrum once again plays an important role: as the speed of sound is much smaller than the electron velocity $v_{g}$, leading-order scattering on acoustic phonons is kinematically suppressed. Consequently, scattering off the optical branch is usually considered [35, 36]. In contrast, we argue that there is another process, the disorder-assisted electron-phonon scattering [37] or “supercollisions” [38, 39, 40, 41], that is responsible for both quasiparticle recombination and energy relaxation. In the high-temperature hydrodynamic regime, the supercollisions are expected to dominate both decay contributions [37]. Moreover, this process contributes weak decay terms to the continuity equations already at local equilibrium, i.e., at the level of “ideal” hydrodynamics.

Our arguments are based on the kinetic theory approach to electronic transport. In the spirit of Ref. 24, we assume the possibility of deriving the hydrodynamic equations from the kinetic equation in the weak coupling limit [21], $\alpha_{g}\ll 1$, with the subsequent renormalization of the kinetic coefficients to the realistic parameter regime [22]. Under these assumptions, we start with the kinetic equation

	$${\cal L}f_{\lambda\boldsymbol{k}}={\rm St}_{ee}[f_{\lambda\boldsymbol{k}}]+{\rm St}_{R}[f_{\lambda\boldsymbol{k}}]+{\rm St}_{\rm dis}[f_{\lambda\boldsymbol{k}}],$$		(1a)
with the Liouville’s operator (in the left-hand side)
	$${\cal L}=\partial_{t}+\boldsymbol{v}\!\cdot\!\boldsymbol{\nabla}_{\boldsymbol{r}}+\left(e\boldsymbol{E}\!+\!\frac{e}{c}\boldsymbol{v}\!\times\!\boldsymbol{B}\right)\!\cdot\!\boldsymbol{\nabla}_{\boldsymbol{k}},$$		(1b)

and the collision integrals describing the electron-electron interaction (${\rm St}_{ee}$), disorder scattering (${\rm St}_{\rm dis}$), and quasiparticle recombination (${\rm St}_{R}$), where in this paper we focus on “supercollisions”. We employ the following notations for the Dirac spectrum (the chirality $\lambda=\pm 1$ distinguishes the conduction and valence bands)

	$$\epsilon_{\lambda\boldsymbol{k}}=\lambda v_{g}k,$$		(2a)
and velocities
	$$\boldsymbol{v}_{\lambda\boldsymbol{k}}=\lambda v_{g}\frac{\boldsymbol{k}}{k},\qquad\boldsymbol{k}=\frac{\lambda k}{v_{g}}\boldsymbol{v}_{\lambda\boldsymbol{k}}=\frac{\epsilon_{\lambda\boldsymbol{k}}\boldsymbol{v}_{\lambda\boldsymbol{k}}}{v_{g}^{2}}.$$		(2b)

Hydrodynamics is the macroscopic manifestation of the conservation of energy, momentum, and the number of particles. In a two-band system, the latter comprises excitations in both bands. In the conductance band these are electron-like quasiparticles with the number density ($N=4$ reflects spin and valley degeneracy in graphene)

	$$n_{+}=N\!\int\!\!\frac{d^{2}k}{(2\pi)^{2}}f_{+,\boldsymbol{k}},$$		(3a)
while in the valence band the quasiparticles are hole-like
	$$n_{-}=N\!\int\!\!\frac{d^{2}k}{(2\pi)^{2}}\left(1-f_{-,\boldsymbol{k}}\right),$$		(3b)
with the total “charge” (or “carrier”) density being
	$$n=n_{+}-n_{-}.$$		(3c)
Assuming the numbers of particles in the conduction and valence bands are conserved independently, we can also define the total quasiparticle (“imbalance” [14]) density
	$$n_{I}=n_{+}+n_{-}.$$		(3d)

Global charge conservation (or gauge symmetry) can be expressed in terms of the usual continuity equation. This can be obtained from Eq. (1) by performing a summation over all quasiparticle states upon which all three collision integrals vanish [27]

	$\displaystyle N\!\int\!\!\frac{d^{2}k}{(2\pi)^{2}}{\rm St}_{ee}[f_{\lambda\boldsymbol{k}}]=N\!\int\!\!\frac{d^{2}k}{(2\pi)^{2}}{\rm St}_{R}[f_{\lambda\boldsymbol{k}}]=$		(4a)
	$\displaystyle\qquad\qquad\qquad\qquad=N\!\int\!\!\frac{d^{2}k}{(2\pi)^{2}}{\rm St}_{\rm dis}[f_{\lambda\boldsymbol{k}}]=0.$
As a result, the continuity equation has the usual form
	$$\partial_{t}n+\boldsymbol{\nabla}_{\boldsymbol{r}}\!\cdot\!\boldsymbol{j}=0,$$		(4b)

where the corresponding current is defined as

$$\boldsymbol{j}\!=\!\boldsymbol{j}_{+}\!-\!\boldsymbol{j}_{-}\!=\!N\!\!\int\!\!\frac{d^{2}k}{(2\pi)^{2}}\left[\boldsymbol{v}_{+,\boldsymbol{k}}f_{+,\boldsymbol{k}}\!-\!\boldsymbol{v}_{-,\boldsymbol{k}}\left(1\!-\!f_{-,\boldsymbol{k}}\right)\right],$$

(5)

The rest of the conservation laws in graphene are approximate as manifested in the collision integrals not vanishing upon corresponding summations. The continuity equation expressing momentum conservation (i.e. the Euler equation) is obtained by multiplying the kinetic equation by the quasiparticle momentum $\boldsymbol{k}$ and summing over all states. Since the electron-electron interaction conserves momentum, the corresponding collision integral vanishes

$$N\!\int\!\!\frac{d^{2}k}{(2\pi)^{2}}\boldsymbol{k}\,{\rm St}_{ee}[f_{\lambda\boldsymbol{k}}]=0.$$

(6)

Weak disorder scattering is typically described within the simplest $\tau$-approximation [27]

$$N\!\int\!\!\frac{d^{2}k}{(2\pi)^{2}}\boldsymbol{k}\,{\rm St}_{\rm dis}[f_{\lambda\boldsymbol{k}}]=\frac{\boldsymbol{n}_{\boldsymbol{k}}}{\tau_{\rm dis}},$$

(7)

where the momentum density is defined as

$$\boldsymbol{n}_{\boldsymbol{k}}=N\!\sum_{\lambda}\!\int\!\frac{d^{2}k}{(2\pi)^{2}}\boldsymbol{k}f_{\lambda\boldsymbol{k}}=v_{g}^{-2}\boldsymbol{j}_{E}.$$

(8)

The last equality represents the fact that in graphene the momentum density is proportional to the energy density [due to the properties of the Dirac spectrum Eq. (2)].

Supercollisions contributing to the recombination collision integral also violate momentum conservation, however, in comparison to the above weak disorder scattering, this is a second-order process. Moreover, the disorder mean free time $\tau_{\rm dis}$ is typically determined from experimental data (see e.g., Ref. 26) and hence can be assumed to include the contribution of supercollisions as well.

The remaining two continuity equations – energy and quasiparticle imbalance – are unaffected by the electron-electron interaction and weak disorder scattering. Indeed, the electron-electron interaction conserves energy and – neglecting the Auger processes – particle number in each band:

$$N\!\int\!\!\frac{d^{2}k}{(2\pi)^{2}}\epsilon_{\lambda\boldsymbol{k}}{\rm St}_{ee}[f_{\lambda\boldsymbol{k}}]=N\!\int\!\!\frac{d^{2}k}{(2\pi)^{2}}\lambda\,{\rm St}_{ee}[f_{\lambda\boldsymbol{k}}]=0.$$

(9)

Same applies to the (elastic) disorder scattering

$$N\!\int\!\!\frac{d^{2}k}{(2\pi)^{2}}\epsilon_{\lambda\boldsymbol{k}}{\rm St}_{\rm dis}[f_{\lambda\boldsymbol{k}}]=N\!\int\!\!\frac{d^{2}k}{(2\pi)^{2}}\lambda\,{\rm St}_{\rm dis}[f_{\lambda\boldsymbol{k}}]=0.$$

(10)

However, supercollisions do not conserve both quantities and hence lead to weak decay terms in the two continuity equations.

Let us now specify the collision integral describing the disorder-assisted electron-phonon scattering. An electron in the upper (conductance) band may scatter into an empty state in the lower (valence) band – effectively recombining with a hole – emitting a phonon (which carries away the energy) and losing its momentum to the impurity. Within the standard approach to the electron-phonon interaction, this process is described by the collision integral

	$\displaystyle{\rm St}_{R}[f_{+\boldsymbol{k}_{2}}]=2\pi\sum_{\boldsymbol{k}_{1},\boldsymbol{q}}W_{q}\delta(\epsilon_{+\boldsymbol{k}_{2}}\!-\epsilon_{-\boldsymbol{k}_{1}}\!-\omega_{q})\times$		(11a)
	$\displaystyle\qquad\times\left[f_{-\boldsymbol{k}_{1}}(1\!-\!f_{+\boldsymbol{k}_{2}})n_{q}-f_{+\boldsymbol{k}_{2}}(1\!-\!f_{-\boldsymbol{k}_{1}})(1\!+\!n_{q})\right],$
where $n_{q}$ is the phonon (Plank’s) distribution function (the phonons are assumed to be at equilibrium and play the role of a “bath”), $W_{q}$ is the effective scattering rate that includes the Dirac factors and is averaged over the angles [37].

Similarly, an electron in the lower band may absorb a phonon and scatter into the upper band – effectively creating an electron-hole pair – while still losing its momentum to the impurity

	$\displaystyle{\rm St}_{R}[f_{-\boldsymbol{k}_{2}}]=2\pi\sum_{\boldsymbol{k}_{1},\boldsymbol{q}}W_{q}\delta(\epsilon_{+\boldsymbol{k}_{1}}\!-\epsilon_{-\boldsymbol{k}_{2}}\!-\omega_{q})\times$		(11b)
	$\displaystyle\qquad\times\left[f_{+\boldsymbol{k}_{1}}(1\!-\!f_{-\boldsymbol{k}_{2}})(1\!+\!n_{q})-f_{-\boldsymbol{k}_{2}}(1\!-\!f_{+\boldsymbol{k}_{1}})n_{q}\right].$

The collision integral (11) conserves the total charge

$$N\sum_{\boldsymbol{k}}{\rm St}_{R}[f_{\lambda\boldsymbol{k}}]=N\sum_{\boldsymbol{k}_{2}}\left({\rm St}_{R}[f_{+\boldsymbol{k}_{2}}]\!+\!{\rm St}_{R}[f_{-\boldsymbol{k}_{2}}]\right)=0,$$

(11c)

[see Eq. (4)] and vanishes in global equilibrium

$${\rm St}_{R}[f^{(0)}]=0,$$

(11d)

where the quasiparticle distribution is described by the Fermi function. This should be contrasted with local equilibrium described by

$$f^{(le)}_{\lambda\boldsymbol{k}}=\left\{1\!+\!\exp\left[\frac{\epsilon_{\lambda\boldsymbol{k}}\!-\!\mu_{\lambda}(\boldsymbol{r})\!-\!\boldsymbol{u}(\boldsymbol{r})\!\cdot\!\boldsymbol{k}}{T(\boldsymbol{r})}\right]\!\right\}^{\!-1}\!,$$

(12)

where ${\mu_{\lambda}(\boldsymbol{r})}$ is the local chemical potential and $\boldsymbol{u}(\boldsymbol{r})$ is the hydrodynamic (or “drift”) velocity. The local equilibrium distribution function (12) allows for independent chemical potentials in the two bands, which can be expressed in terms of the “thermodynamic” and “imbalance” chemical potentials

$$\mu_{\lambda}=\mu+\lambda\mu_{I}.$$

(13)

In global equilibrium

$$f^{(0)}=f^{(le)}_{\lambda\boldsymbol{k}}(\mu_{I}=0,\boldsymbol{u}=0).$$

(14)

Now we show, that in local equilibrium, i.e. for nonzero $\mu_{I}$, the recombination collision integral remains finite [unlike Eq. (11d)]. As a scalar quantity, the collision integral (11) cannot depend on the hydrodynamic velocity $\boldsymbol{u}$ in the first (linear) order. Consequently, to the leading order the integrated collision integral yielding the decay terms in the continuity equations is proportional to $\mu_{I}$.

To the leading order, we can describe the difference between the local equilibrium distribution function $f^{(le)}_{\lambda\boldsymbol{k}}$ and the Fermi function $f^{(0)}$ similarly to the leading non-equilibrium correction in the standard derivation of hydrodynamics [27]

$$\delta f=f_{\lambda\boldsymbol{k}}\!-\!f^{(0)}=-T\frac{\partial f^{(0)}}{\partial\epsilon}h=f^{(0)}\left(1\!-\!f^{(0)}\right)h.$$

(15)

Now we re-write the collision integral (11) with the help of the relations

	$\displaystyle f_{1}(1\!-\!f_{2})(1\!+\!n_{q})-f_{2}(1\!-\!f_{1})n_{q}=$
	$\displaystyle\qquad=(1\!-\!f_{1})(1\!-\!f_{2})(1\!+\!n_{q})\left[\frac{f_{1}}{1\!-\!f_{1}}\frac{n_{q}}{1\!+\!n_{q}}-\frac{f_{2}}{1\!-\!f_{2}}\right],$

and

$$f_{2}^{(0)}\left(1\!-\!f_{1}^{(0)}\right)(1\!+\!n_{q})=-\frac{\partial n_{q}}{\partial\omega}\left(f_{2}^{(0)}\!-\!f_{1}^{(0)}\right),$$

and find (to the leading order in $h$)

			$\displaystyle{\rm St}_{R}[f_{+\boldsymbol{k}_{2}}]=-2\pi\sum_{\boldsymbol{k}_{1},\boldsymbol{q}}W_{q}\delta(\epsilon_{+\boldsymbol{k}_{2}}\!-\epsilon_{-\boldsymbol{k}_{1}}\!-\omega_{q})\frac{\partial n_{q}}{\partial\omega}$		(16a)
			$\displaystyle\qquad\qquad\times\left(f_{+\boldsymbol{k}_{2}}^{(0)}\!-\!f_{-\boldsymbol{k}_{1}}^{(0)}\right)(h_{-\boldsymbol{k}_{1}}\!-\!h_{+\boldsymbol{k}_{2}}),$
	$\displaystyle{\rm St}_{R}[f_{-\boldsymbol{k}_{2}}]=2\pi\sum_{\boldsymbol{k}_{1},\boldsymbol{q}}W_{q}\delta(\epsilon_{-\boldsymbol{k}_{2}}\!-\epsilon_{+\boldsymbol{k}_{1}}\!+\omega_{q})\frac{\partial n_{q}}{\partial\omega}$
	$\displaystyle\qquad\qquad\times\left(f_{+\boldsymbol{k}_{1}}^{(0)}\!-\!f_{-\boldsymbol{k}_{2}}^{(0)}\right)(h_{-\boldsymbol{k}_{2}}\!-\!h_{+\boldsymbol{k}_{1}}).$				(16b)

Consider now the contribution of the recombination collision integral to the continuity equation for the quasiparticle imbalance

	$\displaystyle N\sum_{\boldsymbol{k}}\lambda\,{\rm St}_{R}[f_{\lambda\boldsymbol{k}}]=N\sum_{\boldsymbol{k}_{2}}\left({\rm St}_{R}[f_{+\boldsymbol{k}_{2}}]\!-\!{\rm St}_{R}[f_{-\boldsymbol{k}_{2}}]\right)=$
	$\displaystyle\qquad=-4\pi N\!\sum_{\boldsymbol{k}_{1},\boldsymbol{k}_{2},\boldsymbol{q}}\!W_{q}\delta(\epsilon_{+\boldsymbol{k}_{2}}\!-\epsilon_{-\boldsymbol{k}_{1}}\!-\omega_{q})\frac{\partial n_{q}}{\partial\omega}$
	$\displaystyle\qquad\qquad\times\left(f_{+\boldsymbol{k}_{2}}^{(0)}\!-\!f_{-\boldsymbol{k}_{1}}^{(0)}\right)(h_{-\boldsymbol{k}_{1}}\!-\!h_{+\boldsymbol{k}_{2}}).$		(17)

To the leading order, the deviation $h_{\lambda\boldsymbol{k}}$ is proportional to $\mu_{I}$

$$h_{\lambda\boldsymbol{k}}\approx\frac{\lambda\mu_{I}}{T}.$$

(18)

The remaining integral has dimensions of particle density divided by time and therefore the result can be written in two equivalent forms

$\displaystyle N\sum_{\boldsymbol{k}}\lambda\,{\rm St}_{R}[f_{\lambda\boldsymbol{k}}]\approx-\mu_{I}n_{I,0}\lambda_{Q}\approx-\frac{n_{I}\!-\!n_{I,0}}{\tau_{R}}.$

(19)

Here $n_{I}$ is the imbalance density (3d) in local equilibrium, while $n_{I,0}$ is the same quantity evaluated with the Fermi distribution function (14), i.e. for $\mu_{I}=0$ and $\boldsymbol{u}=0$. The first equality in Eq. (19) coincides with the expression used in Ref. 14 and serves as the definition of the dimensionless coefficient $\lambda_{Q}$, while the second (valid to the leading order) was suggested in Refs. 21, 29 and provides the definition of the “recombination time” $\tau_{R}$ (see also Ref. 28). The two expressions are equivalent since ${n_{I}\!-\!n_{I,0}\propto\mu_{I}}$.

The same scattering process contributes a weak decay term to the continuity equation for the energy density. Indeed, multiplying the collision integral (11) by the quasiparticle energy and summing over all states, we find after similar algebra

	$\displaystyle N\sum_{\boldsymbol{k}}\epsilon_{\lambda\boldsymbol{k}}{\rm St}_{R}[f_{\lambda\boldsymbol{k}}]=$		(20)
	$\displaystyle\quad=N\sum_{\boldsymbol{k}_{2}}\left(\epsilon_{+\boldsymbol{k}_{2}}{\rm St}_{R}[f_{+\boldsymbol{k}_{2}}]\!+\!\epsilon_{-\boldsymbol{k}_{2}}{\rm St}_{R}[f_{-\boldsymbol{k}_{2}}]\right)=$
	$\displaystyle\quad=-2\pi N\!\sum_{\boldsymbol{k}_{1},\boldsymbol{k}_{2},\boldsymbol{q}}W_{q}\delta(\epsilon_{+\boldsymbol{k}_{2}}\!-\epsilon_{-\boldsymbol{k}_{1}}\!-\omega_{q})\frac{\partial n_{q}^{(0)}}{\partial\omega}\omega_{q}$
	$\displaystyle\qquad\qquad\times\left(f_{+\boldsymbol{k}_{2}}^{(0)}\!-\!f_{-\boldsymbol{k}_{1}}^{(0)}\right)(h_{-\boldsymbol{k}_{1}}\!-\!h_{+\boldsymbol{k}_{2}}).$

Defining the decay coefficients similarly to Eq. (19) above, we may present the result in the form

$$N\sum_{\boldsymbol{k}}\epsilon_{\lambda\boldsymbol{k}}{\rm St}_{R}[f_{\lambda\boldsymbol{k}}]=-\mu_{I}n_{E,0}\lambda_{QE}\approx-\frac{n_{E}\!-\!n_{E,0}}{\tau_{RE}}.$$

(21)

Here the equivalence of the two forms of the decay term stems from the fact that ${n_{E}\!-\!n_{E,0}\propto\mu_{I}}$ assuming the electrons and holes are characterized by the same temperature.

Supercollisions are not the only scattering process contributing to both quasiparticle recombination and energy relaxation. Clearly, direct (not impurity-assisted) electron-phonon interaction contributes to energy relaxation as well as to quasiparticle recombination (in the case of intervalley scattering) [34, 18, 37, 14, 19]. In addition, optical phonons may also contribute [35, 36], although within the hydrodynamic approach these contributions were discussed only at the level of dissipative (viscous) hydrodynamics [36]. The contribution of the direct [42, 43] and impurity assisted electron-phonon scattering to energy relaxation was compared in [37], where it was argued that at high enough temperatures, $T\gtrsim T_{BG}$ (where $T_{BG}$ is the Bloch-Grüneisen temperature) the supercollisions do dominate. In the degenerate regime, where $T_{BG}\propto\sqrt{n}$, Ref. [37] estimates $T_{BG}$ as “few tens of Kelvin”. At charge neutrality, we estimate $T_{BG}=(s/v_{g})T\ll T$ (where $s$ is the speed of sound), such that supercollisions always dominate over direct electron-phonon coupling. Moreover, taking into account the additional scattering processes will not change the form of the decay terms (19) and (21), but rather change the numerical values of the parameters $\lambda_{Q}$ and $\lambda_{QE}$, which may have to be considered phenomenological while interpreting experimental data [34].

At charge neutrality and in the hydrodynamic regime, the coefficients $\lambda_{Q}$ and $\lambda_{QE}$ are of the same order of magnitude (both are dominated by the same supercollisions), but quantitatively different. Indeed, the continuity equation for the energy density should contain one more contribution of the similar form. “Quasiparticle recombination” typically refers to scattering between the quasiparticle states in different bands only. This is the only type of supercollisions contributing to $\lambda_{Q}$. Similar supercollisions may also take place within a single band [37]. While this process does not change the number of particles in the band, it does describe the energy loss as the electron may scatter from the higher energy state into the lower energy state (and losing its momentum to the impurity along the way). Consequently, this additional process does contribute to energy relaxation. Given that the form of the corresponding collision integral is very similar to Eq. (11) – one only has to change to band indices to be the same – the algebra remains the same and thus we can treat Eq. (21) as the final result that takes this additional intraband supercollisions into account making the numerical values of $\lambda_{Q}$ and $\lambda_{QE}$ substantially different – we do not expect any accidental cancellation or smallness should the difference $\lambda_{Q}-\lambda_{QE}$ appear in a particular solution of hydrodynamic equations. At the same time, at low temperatures – i.e., below the hydrodynamic range – we expect the coefficients $\lambda_{Q}$ and $\lambda_{QE}$ to be parametrically different: energy relaxation is now dominated by the direct electron-phonon interaction [37], while the recombination is still governed by supercollisions (together with the phonon-induced intervalley scattering).

The order of magnitude of $\tau_{R}$ could be estimated based on the calculations of Ref. [37]. Adapting the latter to charge neutrality, we find $\tau_{R}^{-1}\sim D^{2}T^{2}/(\rho s^{2}v_{g}^{2})$ (where $D\approx 20\,$eV is the deformation potential [42, 43] and $\rho$ is the mass density per unit area) yielding the corresponding length scale $\ell_{R}^{-1}\approx 10\,\mu$m at the top of the hydrodynamic temperature range, $T\approx 250\,$K. This should be further compared to the contribution of three-particle collisions [3, 36], $\tau_{3}^{-1}\sim\alpha_{g}^{4}T$. Assuming the common sample design where graphene is encapsulated in hexagonal boron nitride (with the dielectric constant $\epsilon\approx 4$), the effective coupling constant (taking into account renormalizations) is $\alpha_{g}\approx 0.3-0.4$ leading to the similar estimate at high temperatures. On the other hand, at the low end of the hydrodynamic range [44], $T\approx 50\,$K, the contribution of the three-body collisions should dominate (accounting for the empirical value $\ell_{R}\approx 1.2\,\mu$m reported in [34]), however preserving the functional form of the weak decay terms in the continuity equations.

Finally, once the dissipative processes due to electron-electron interaction are taken into account, one usually replaces the continuity equation for the energy density by the equivalent equation for the entropy density, the so-called “heat transfer equation” [4]. The decay terms discussed in this paper appear in that equation as well. Let us briefly discuss their form.

Recall the derivation of the continuity equation for the entropy from the kinetic equation [21]. The entropy density of a system of fermions is defined in terms of the distribution function as

$$s\!=\!-\!N\!\sum_{\lambda}\!\!\int\!\!\frac{d^{2}k}{(2\pi)^{2}}\!\left[f_{\lambda\boldsymbol{k}}\ln f_{\lambda\boldsymbol{k}}\!+\!(1\!-\!f_{\lambda\boldsymbol{k}})\!\ln(1\!-\!f_{\lambda\boldsymbol{k}})\right]\!.$$

(22)

Treating this integral as

$$s=N\!\sum_{\lambda}\!\int\!\frac{d^{2}k}{(2\pi)^{2}}{\cal S}[f_{\lambda\boldsymbol{k}}],$$

any derivative of $s$ can be represented in the form

$$\frac{\partial s}{\partial z}=N\!\sum_{\lambda}\!\int\!\frac{d^{2}k}{(2\pi)^{2}}\frac{\partial{\cal S}[f_{\lambda\boldsymbol{k}}]}{\partial f_{\lambda\boldsymbol{k}}}\frac{\partial f_{\lambda\boldsymbol{k}}}{\partial z}.$$

Consider now each term of the kinetic equation multiplied by the derivative $\partial{\cal S}[f_{\lambda\boldsymbol{k}}]/\partial f_{\lambda\boldsymbol{k}}$ and summed over all states. Using the above relation with ${z\rightarrow{t}}$, one finds for the time derivative term

$$N\!\sum_{\lambda}\!\int\!\frac{d^{2}k}{(2\pi)^{2}}\frac{\partial{\cal S}[f_{\lambda\boldsymbol{k}}]}{\partial f_{\lambda\boldsymbol{k}}}\frac{\partial f_{\lambda\boldsymbol{k}}}{\partial t}=\frac{\partial s}{\partial t}.$$

The gradient term yields similarly

	$\displaystyle N\!\sum_{\lambda}\!\int\!\frac{d^{2}k}{(2\pi)^{2}}\frac{\partial{\cal S}[f_{\lambda\boldsymbol{k}}]}{\partial f_{\lambda\boldsymbol{k}}}\boldsymbol{v}_{\lambda\boldsymbol{k}}\cdot\boldsymbol{\nabla}_{\boldsymbol{r}}f_{\lambda\boldsymbol{k}}=$
	$\displaystyle\qquad\qquad=\boldsymbol{\nabla}_{\boldsymbol{r}}\cdot N\!\sum_{\lambda}\!\int\!\frac{d^{2}k}{(2\pi)^{2}}\boldsymbol{v}_{\lambda\boldsymbol{k}}{\cal S}[f_{\lambda\boldsymbol{k}}]=\boldsymbol{\nabla}_{\boldsymbol{r}}\!\cdot\!\boldsymbol{j}_{S},$

where the quantity

$$\boldsymbol{j}_{S}=N\!\sum_{\lambda}\!\int\!\frac{d^{2}k}{(2\pi)^{2}}\boldsymbol{v}_{\lambda\boldsymbol{k}}{\cal S}[f_{\lambda\boldsymbol{k}}],$$

(23)

can be interpreted as the entropy current.

The electric field term vanishes as the total derivative

	$\displaystyle e\boldsymbol{E}\cdot N\!\sum_{\lambda}\!\int\!\frac{d^{2}k}{(2\pi)^{2}}\frac{\partial{\cal S}[f_{\lambda\boldsymbol{k}}]}{\partial f_{\lambda\boldsymbol{k}}}\boldsymbol{\nabla}_{\boldsymbol{k}}f_{\lambda\boldsymbol{k}}=$
	$\displaystyle\qquad\qquad=e\boldsymbol{E}\cdot N\!\sum_{\lambda}\!\int\!\frac{d^{2}k}{(2\pi)^{2}}\boldsymbol{\nabla}_{\boldsymbol{k}}{\cal S}[f_{\lambda\boldsymbol{k}}]=0,$

while the Lorentz term vanishes for rotationally invariant systems upon integrating by parts [justified by the fact that ${{\cal{S}}(k\rightarrow\infty)\rightarrow 0}$]

	$\displaystyle\frac{e}{c}N\!\sum_{\lambda}\!\int\!\frac{d^{2}k}{(2\pi)^{2}}\frac{\partial{\cal S}[f_{\lambda\boldsymbol{k}}]}{\partial f_{\lambda\boldsymbol{k}}}\left[\boldsymbol{v}_{\lambda\boldsymbol{k}}\!\times\!\boldsymbol{B}\right]\!\cdot\!\boldsymbol{\nabla}_{\boldsymbol{k}}f_{\lambda\boldsymbol{k}}=$
	$\displaystyle\qquad\qquad=\frac{e}{c}N\!\sum_{\lambda}\!\int\!\frac{d^{2}k}{(2\pi)^{2}}\left[\boldsymbol{v}_{\lambda\boldsymbol{k}}\!\times\!\boldsymbol{B}\right]\!\cdot\!\boldsymbol{\nabla}_{\boldsymbol{k}}{\cal S}[f_{\lambda\boldsymbol{k}}]$
	$\displaystyle\qquad\qquad=-\frac{e}{c}N\!\sum_{\lambda}\!\int\!\frac{d^{2}k}{(2\pi)^{2}}{\cal S}[f_{\lambda\boldsymbol{k}}]\boldsymbol{\nabla}_{\boldsymbol{k}}\!\cdot\!\left[\boldsymbol{v}_{\lambda\boldsymbol{k}}\!\times\!\boldsymbol{B}\right]=0.$

The last equality follows from

$$\frac{\partial v^{\alpha}_{\lambda\boldsymbol{k}}}{\partial k^{\beta}}=\frac{v_{g}}{\lambda k}\left(\delta_{\alpha\beta}-\frac{k^{\alpha}k^{\beta}}{k^{2}}\right).$$

Similar approach was used in Ref. 21 to derive the continuity equations (as outlined above).

Combining all four terms, we conclude that integration with the factor ${\partial{\cal{S}}[f_{\lambda\boldsymbol{k}}]/\partial{f}_{\lambda\boldsymbol{k}}}$ turns the left-hand side of the kinetic equation to the familiar form

$$N\!\sum_{\lambda}\!\int\!\frac{d^{2}k}{(2\pi)^{2}}\frac{\partial{\cal S}[f_{\lambda\boldsymbol{k}}]}{\partial f_{\lambda\boldsymbol{k}}}{\cal L}f_{\lambda\boldsymbol{k}}=\frac{\partial s}{\partial t}+\boldsymbol{\nabla}_{\boldsymbol{r}}\!\cdot\!\boldsymbol{j}_{S}.$$

(24)

Eq. (24) is valid for an arbitrary distribution function.

Denoting the integral of the right-hand side of the kinetic equation by

$${\cal I}=N\!\sum_{\lambda}\!\!\int\!\!\frac{d^{2}k}{(2\pi)^{2}}\frac{\partial{\cal S}[f_{\lambda\boldsymbol{k}}]}{\partial f_{\lambda\boldsymbol{k}}}\!\left({\rm St}_{ee}[f]+{\rm St}_{R}[f]+{\rm St}_{\rm dis}[f]\right)\!,$$

(25)

we arrive at the “continuity equation for the entropy”

$$\frac{\partial s}{\partial t}+\boldsymbol{\nabla}_{\boldsymbol{r}}\!\cdot\!\boldsymbol{j}_{S}={\cal I}.$$

(26)

In the usual hydrodynamics [4] the only contribution to the collision integral is given by particle-particle scattering, i.e. the processes assumed to be responsible for establishing local equilibrium such that at ${{\cal I}=0}$ the ideal (Euler) hydrodynamic is isentropic. In the present case, local equilibrium is assumed to be achieved by means of the electron-electron interaction.

Evaluating the derivative of ${\cal S}$ explicitly, we find

$$\frac{\partial{\cal S}[f_{\lambda\boldsymbol{k}}]}{\partial f_{\lambda\boldsymbol{k}}}=-\ln\frac{f_{\lambda\boldsymbol{k}}}{1\!-\!f_{\lambda\boldsymbol{k}}}=\ln\left[\frac{1}{f_{\lambda\boldsymbol{k}}}\!-\!1\right].$$

For the local equilibrium distribution function

$$\frac{\partial{\cal S}[f_{\lambda\boldsymbol{k}}]}{\partial f_{\lambda\boldsymbol{k}}}=\frac{\epsilon_{\lambda\boldsymbol{k}}\!-\!\mu_{\lambda}\!-\!\boldsymbol{u}\!\cdot\!\boldsymbol{k}}{T}.$$

Substituting this expression into Eq. (25), we find that the remaining integration is very similar to the above derivation of the continuity equations.

The integral with the quasiparticle energy yields exactly the above Eq. (21). The integral with $\lambda\mu_{I}$ yields Eq. (19) multiplied by $\mu_{I}$. Finally, the term $\boldsymbol{u}\!\cdot\!\boldsymbol{k}$ yields Eq. (7) multiplied by the hydrodynamic velocity. The integral of this term with the recombination collision integral is assumed to be included into the definition of the mean free time, see the corresponding discussion above. As a result, we arrive at the following form of the integrated collision integral

$\displaystyle{\cal I}=-\frac{1}{T}\frac{n_{E}\!-\!n_{E,0}}{\tau_{RE}}+\frac{\mu_{I}}{T}\frac{n_{I}\!-\!n_{I,0}}{\tau_{R}}+\frac{\boldsymbol{u}\!\cdot\!\boldsymbol{n}_{\boldsymbol{k}}}{T\tau_{\rm dis}}.$

The decay terms (Hydrodynamic approach to electronic transport in graphene: energy relaxation) appear already at local equilibrium. To complete the heat transfer equation one has to take into account dissipation. In graphene, this is most conveniently done by considering the classical limit of relativistic hydrodynamics since the form of dissipative corrections is determined by the symmetries of the quasiparticle spectrum. The result has been already reported in literature, therefore we combine the dissipative corrections with Eq. (Hydrodynamic approach to electronic transport in graphene: energy relaxation) and write the heat transfer equation in graphene in the form (here $\delta\boldsymbol{j}$ and $\delta\boldsymbol{j}_{I}$ are the dissipative corrections to the electric and imbalance currents, respectively).

	$\displaystyle T\left[\frac{\partial s}{\partial t}+\boldsymbol{\nabla}_{\boldsymbol{r}}\!\cdot\!\left(s\boldsymbol{u}-\delta\boldsymbol{j}\frac{\mu}{T}-\delta\boldsymbol{j}_{I}\frac{\mu_{I}}{T}\right)\right]=$		(27)
	$\displaystyle\qquad\qquad=\delta\boldsymbol{j}\!\cdot\!\left[e\boldsymbol{E}\!+\!\frac{e}{c}\boldsymbol{u}\!\times\!\boldsymbol{B}\!-\!T\boldsymbol{\nabla}\frac{\mu}{T}\right]-T\delta\boldsymbol{j}_{I}\!\cdot\!\boldsymbol{\nabla}\frac{\mu_{I}}{T}$
	$\displaystyle\qquad\qquad\qquad+\frac{\eta}{2}\left(\nabla_{\alpha}u_{\beta}\!+\!\nabla_{\beta}u_{\alpha}\!-\!\delta_{\alpha\beta}\boldsymbol{\nabla}\!\cdot\!\boldsymbol{u}\right)^{2}$
	$\displaystyle\qquad\qquad\qquad-\frac{n_{E}\!-\!n_{E,0}}{\tau_{RE}}+\mu_{I}\frac{n_{I}\!-\!n_{I,0}}{\tau_{R}}+\frac{\boldsymbol{u}\!\cdot\!\boldsymbol{n}_{\boldsymbol{k}}}{\tau_{\rm dis}}.$

The obtained equation (27) should be compared to the corresponding equations in Refs. 3, 14, 36, where energy relaxation due to supercollisions were not taken into account. All other terms are present in all four equations with the following exceptions. The equation (54) of Ref. 3 is written in the relativistic notation omitting the imbalance mode, quasiparticle recombination, and disorder scattering, all of which are discussed separately elsewhere in Ref. 3. Ref. 14 was the first to focus on the imbalance mode with the equation (2.6) containing all the terms of Eq. (27) except for the viscous term (and energy relaxation). Finally, the equation (1c) of Ref. 36 contains all of the terms in Eq. (27) except for energy relaxation and in addition contains a term describing energy relaxation due to electrons scattering on the optical phonon branch that is neglected here (generalization of the resulting theory is straightforward).

To summarize, we have considered supercollisions as a mechanism of quasiparticle recombination and energy relaxation in graphene and derived the corresponding decay terms in the hydrodynamic continuity equations. Since the same scattering mechanism is responsible for both effects, one has to take into account energy relaxation while considering quasiparticle recombination. The latter is an indispensable feature of electronic hydrodynamics in graphene in constrained geometries, where homogeneous solutions violate the boundary conditions [29].