Hydrodynamic equations for an electron gas in graphene

In this paper we review, and extend to the non-isothermal case, some results concerning the application of the maximum entropy closure technique to the derivation of hydrodynamic equations for particles with spin-orbit interaction and Fermi-Dirac statistics. In the second part of the paper we treat in more details the case of electrons on a graphene sheet and investigate various asymptotic regimes.


Introduction
This paper is devoted to present some results on the derivation of hydrodynamic equations describing electrons subject to spin-orbit-like interactions. Systems of this kind, of particular interest for applications to microelectronics, include electrons undergoing the so-called Rashba effect [1], the Kane's two-band K·P model [2] and electrons in singlelayer graphene [3]. The diffusive and hydrodynamic descriptions of such systems are extensively treated in Refs. [4,5,6,7,8,9,10]. Here, we summarize the results contained in Refs. [7,8,9,10], concerning the hydrodynamic description, and extend them to the non-isothermal case.
In comparison with kinetic models, the advantages of fluid models for applications are evident. In fact, from the numerical point of view, a system of PDEs for a set of macroscopic quantities is much more desirable than a single equation for a density in phase-space, where also the components of momentum are independent variables. Moreover, from the point of view of mathematical modeling, they offer more flexibility, as various kind of boundary conditions and coupling terms (e.g. with a self-consistent potential or with various types of scattering mechanisms) can be very naturally embodied in the mathematical model.
On the other hand, the derivation of fluid equations for the systems under consideration, which possess spinorial degrees of freedom and non-parabolic dispersion relations (energy bands), is far from being a trivial extension of the techniques employed for standard (i.e. scalar and parabolic) particles. The best strategy to obtain hydrodynamic (or, more in general, fluid-dynamic) equations in this case, seems to be their systematic derivation from an underlying kinetic description by means of the Maximum Entropy Principle (MEP) and its quantum extensions [11,12,13,14]. The MEP basically stipulates that the microscopic (kinetic) state of the system is the most probable among all states sharing the same macroscopic moments of interest, providing therefore a formal closure of the system of moment equations. This is a very general principle which finds a variety of applications to different fields, ranging from statistical mechanics to signal theory [15]. For quantum systems it can be used in combination with the quantum kinetic framework provided by the phase-space formulation of quantum mechanics due to Wigner [16,17].
In the present work the Wigner formalism is used "semiclassically", which means that some quantum features are retained (namely, the peculiar energy-band dispersion relations and the Fermi-Dirac statistics) while others are neglected (namely, the quantum coherence between different bands). Consequently, the obtained hydrodynamic description misses some interesting physics when quantum interference between bands becomes important (e.g. close to abrupt potential variations [18,19]). Nevertheless, the derived equations possess an interesting mathematical structure and reveal some interesting physics, still occurring in absence of such interference phenomena (see, in particular, Section 5).
The paper is organized as follows. In Section 2 the kinetic-level formalism, based on a semiclassical Wigner description, is introduced for a fairly general spin-orbit Hamiltonian that includes all cases of interest. In Sec. 3 we write the moment equations for density, velocity and energy, and perform their formal closure by means of the MEP. Then, the second part of the paper is focused on the case of graphene. In Sec. 4 the general theory exposed in the first part is specialized for the Dirac-like Hamiltonian describing electrons on a single-layer graphene sheet. In Sec. 5 we obtain the asymptotic form of the hydrodynamic equations derived in Sec. 4 in some physically relevant limits (namely, the high temperature, zero temperature, collimation and diffusive limits). Finally, Sec. 6 is devoted to conclusions and perspectives.
2 Phase-space description of spin-orbit particles Let us consider a spin-orbit Hamiltonian of the form

Moreover, the dot product is defined as
Hamiltonians of this kind describe various systems of great interest in solid-state physics. The first example is a 2-dimensional electron gas confined in an asymmetric potential well, which is subject to the Rashba spin-orbit interaction [1,4]. In this case: where p = (p x , p y ), m * is the electron effective mass, e z is the normal direction to the well and α is the Rashba constant. Another example is the two-band K·P model [2,9,10,20]; in this case: where m is the electron (bare) mass, E g is the band-gap and K is the matrix element of the gradient operator between conduction and valence Bloch functions.
The last example is that of electrons on a single-layer graphene sheet [3,7,21], in which case: This case will be considered in more details in the second part of the paper.
We remark that the variable p has to be interpreted as the crystal pseudo-momentum, rather than the ordinary momentum. The interpretation of the vector variable σ depends on the cases: it is (proportional to) the spin vector in the case of Rashba Hamiltonian, while it is a pseudo-spin in the other two examples [2,3]. For graphene, in particular, the pseudo-spin is related to the decomposition of the honeycomb lattice into two inequivalent sublattices, which reflects the presence of two carbon atoms in the fundamental cell of the lattice [3,22].
The main semiclassical quantities associated with (1) are: 1. the two energy bands (i.e., the eigenvalues of H with V = 0); 2. the projectors on the eigenspaces corresponding to E ± (p), where is the pseudo-momentum direction; 3. the semiclassical velocities v ± (p) = ∇ p E ± (p); 4. the effective-mass tensor [20] M −1 Note, in particular, that the eigenvalues of the projector P ± are 1 and 0, corresponding to whether or not the electron energy belongs to the upper/lower energy band. Hence, the expected value of P ± can be interpreted as the fraction of electrons belonging to the upper/lower band (see below).
The phase-space description of a statistical population of electrons with Hamiltonian (1) is provided by the Wigner matrix [9,16,17,23] F (x, which is the Wigner transform, Such representation of a quantum mixed-state has the fundamental property that the expected value of an observable with symbol A = 3 k=0 a k (x, p)σ k is given by the classicallooking formula By applying Eq. (8) to the band projectors P ± (p) we obtain and it is therefore natural to interpret the functions as the phase-space densities of electrons having energies, respectively, in the upper and lower band.
Let us now consider the following hydrodynamic moments of electrons in the two bands: n ± = f ± , (density), (see also Ref. [14] where additional moments are considered). Here we have introduced the shorthand The (semiclassical) dynamics of the Wigner matrix (7) is provided by the Wigner equations for the Hamiltonian (1) [9], with i = 1, 2, 3. From (11), the following equations for the band-Wigner functions f + and f − (see definition (9)) are readily obtained: where the terms containing are responsible for quantum interference between the two bands [9,23].

Maximum entropy closure
In order to obtain from (12) a closed system of equations for the moments (10), we assume that the system is in a state F me of maximum entropy, according to the so-called Maximum Entropy Principle (MEP) [11,12,13,14,15] which in the present case reads as follows: MEP F me is the most probable microscopic state with the observed macroscopic moments n ± , u ± and e ± .
Hence, we search for a Wigner matrix F me that maximizes the total entropy In (13), k B is the Boltzmann constant, Tr is the matrix trace and is (minus) the Fermi-Dirac entropy function. The condition 0 ≤ F ≤ 1 ensures that s(F ) is a well-defined matrix. It can be proven [5] that where A ± , B ± = (B 1 , . . . , B d ) ± and C ± are Lagrange multipliers (functions of x and t), and, moreover, Thus, the (semiclassical) MEP state corresponds to two local Fermi-Dirac distributions in the two energy bands. In particular, Eq. (17) implies that, in such state, the interference terms vanish and, therefore, the two bands are decoupled (unless additional coupling mechanisms are considered [6,10,14]). Hence, from now on, we shall treat the two bands separately and, in order to simplify notations, the ± labels will be suppressed (except where a distinction between quantities taking different forms in the two bands, such as E ± or v ± , is necessary). Using (16) and (17) in (12) (and suppressing the ± labels, as it was just explained) yields By taking the moments · , v ± · and E ± · of both sides of Eq. (18), and recalling the definitions (5) and (6), we obtain the moment equations Thanks to the MEP, the moment system (19) is implicitly closed by the constraints (14), linking the Lagrange multipliers (A, B, C) to the moments (n, u, e) thus allowing (in principle) to think to f me as being parametrized by (n, u, e) and, consequently, the extra moments (P ± , Q ± , S ± ) as functions of the unknowns (n, u, e). Following Levermore [11], we can express the moments of the MEP state f me as the derivatives with respect to the Lagrange multipliers of the "density potential" ε * , which is defined as the Legendre transform of the entropy density where s is given by (15) and f me by (16). It is not difficult to show that and that the constraint equations (14) may be rewritten as where i = 1, . . . , d. Levermore's theory, moreover, ensures that system (19), with the closure relations (20) is hyperbolic and, therefore, it is at least locally well-posed (see also Ref. [7]).

The case of graphene
We now specialize the formalism introduced so far to the case of a population of electrons on a single-layer graphene sheet. Such electrons, in the proximity of a Dirac point in pseudo-momentum space [3], are described by the Hamiltonian (1) with (where c ≈ 10 6 m/s is the Fermi velocity), which corresponds to a Dirac-like Hamiltonian for relativistic, massless particles. We remark that this is an approximation which is valid only in the proximity of a Dirac point for an infinite, ideal and un-doped system (see Ref. [21] and references therein). In this case the energy bands are the Dirac cones and the eigenprojections are given by where Moreover, the semiclassical velocities are implying that electrons travel with the constant speed c and direction ν, and the effectivemass tensor is where ν ⊥ = (−ν 2 , ν 1 ).
Since the lower band is unbounded from below, we have to change a little the theory developed in the previous sections and describe the lower-band population in terms of electron vacancies, i.e. holes. This is achieved by means of the substitution which brings the transport equation, Eq. (18), into Note that the only difference between electrons and holes is the charge sign. Moreover, the MEP-states for electrons and holes have now the form in fact, both upper-cone electrons and lower-cone holes have positive energies (note that in (27) the Fermi velocity c has been absorbed in the Lagrange multiplier C). Moreover, we slightly change the definition of u to be the average direction which differs from average velocity just for the constant factor c. The inequality |u| ≤ 1 is an obvious consequence of the fact that u is an average of directions. The moment equations (19), in the specific case of graphene, read as follows: where the higher-order moments P ij , Q ij and S j take the form We now intend to find an (as much as possible) explicit expression for the dependence of the Lagrange multipliers A, B = (B 1 , B 2 ) and C in terms of the moments n, u = (u 1 , u 2 ) and e, as resulting from the constraint equations By expressing the integrals over p ∈ R 2 in polar coordinates, we obtain the expressions where and φ s is the Fermi integral of order s > 0: It is now convenient to put so that the previous expressions can be rewritten as We remark that the new Lagrange multiplier T has the physical meaning of the electron gas temperature. From (36) and the constraint equations (32), we obtain that B has the same direction as u and that (n, |u|, e) are related to the scalar Lagrange multipliers (A, B, T ) by I 2 0 (A, B)n T = n, Similarly to what is found in Ref. [7], we obtain the following expressions of the higherorder moments (31) in terms of n, u, T and the functions I s N = I s N (A, B): where u ⊥ = (−u 2 , u 1 ).

Asymptotic regimes
The expressions (38) of P ij , Q ij and S j are still not explicit, as functions of n and u. In fact, these expression depend, through the functions I s N (A, B), on the two scalar Lagrange multipliers A and B, which are related to n and |u| via the relations (37). In Ref. [7] it has been proven that the correspondence between (A, B) and (n, |u|) is 1-1 but, as far as we know, it is not possible to give an explicit, analytic, expression of the former as functions of the latter.
However, we can say more in some particular regimes of physical interest. Such regimes correspond to different asymptotic regions [7] in the half plane (A, B) ∈ R×[0, ∞), namely: 4. opposite to the collimation limit, the asymptotic region B → 0 corresponds to the diffusive limit |u| → 0, where the velocities are randomly spread over all directions.
The asymptotic analysis of equations (30) in these regimes is based on the following result, which has been proven in Ref. [7].
Theorem. The functions I s N have the following asymptotic behavior: 1. in the Maxwell-Boltzmann limit, where I N are the modified Bessel functions of the first kind; 2. in the degenerate gas limit, where

Maxwell-Boltzmann regime
The Maxwell-Boltzmann regime is the limit for large T and corresponds to A 2 + B 2 → ∞ with A < −B in the (A, B) half plane. In this case, we can use the approximation (39). Note, in particular, that in such limit the functions I s N become factorized and independent on the index s. Then, the constraint equations (37) become e A I 0 (B)n T = n, and it can be shown that the the MEP-state (28) is well approximated by the Maxwellianlike distribution where Moreover, we get the explicit form of P ij , Q ij and S j : where and B is given by (44). By playing a little with the asymptotic expansions of the modified Bessel functions I n we obtain the asymptotic behavior of X(|u|) in the diffusive limit: and in the collimation limit: Substituting S j = neu j in the third of the moment equations (30) yields, after a little algebra, Thus, the isothermal case (e = 2k B T constant) is only compatible with u j ∂ j V = 0, i.e. the component of the force field parallel to the velocity field must vanish. In this case, the pseudo-momentum balance equation reduces to

Degenerate gas regime
The degenerate gas regime is the limit for T → 0 and corresponds to A 2 + B 2 → ∞ with A > −B, in the (A, B) half plane. In this case, we can use the approximation (40)-(41).

It is convenient to put
and rewrite (40) as follows: where F s N (ψ) = 1 πΓ(s + 1) C(ψ) 0 cos(N θ)(cos ψ + sin ψ cos θ) s dθ and The asymptotic form of the constraint equations (37) is now Note that: 1. |u| only depends on ψ and we can write 2. recalling (35), from the first of the above equations we have R ∼ 1/T and, then, the third equation shows that e remains positive even though T → 0.
From the above considerations it is readily seen that, in the limit T → 0, the MEP-state (28) takes the typical degenerate Fermi-Dirac form where θ denotes the Heaviside function and ψ(u) is given by (54). Using (38), (50) and (53) we obtain the following expressions of P ij and Q ij and S j for a degenerate electron gas: where , and ψ = ψ(|u|) is given by Eq. (54). By using the techniques developed in Ref. [7], it is not difficult to calculate the asymptotic behavior of the functions Y (|u|), Z(|u|), Z ⊥ (|u|) and W (|u|) in the two limits |u| → 0 (diffusion) and |u| → 1 (collimation).

Collimation regime
The collimation limit corresponds to the absence of spread in the particle directions, i.e. to |u| → 1. It can be shown [7,8] that this limit is equivalent to However, there is a completely different behavior when the critical line A = −B is approached from below (Maxwell-Boltzmann collimation) of from above (degenerate gas collimation).
The first case corresponds to taking the limit B → ∞ in the "Maxwellian" distribution (59), which produces a delta in the angle between p and u, namely Moreover, since X(|u|) → 1 as |u| → 1 (see Eq. (47)), from Eq. (45) we obtain and the pseudo-momentum balance equation reduces to By using the continuity equation ∂ t n + c∂ j (nu j ) the latter can be rewritten as which is decoupled from the continuity equation for n. As pointed out in Refs. [7,8], this equation reveals that collimated electrons in graphene have the properties of a geometricaloptics system, with "refractive index" By also considering the energy balance equation (48), we finally obtain the system In order to derive the collimation equations for a degenerate gas, we start from the expression (56) of P ij , Q ij and S j , and use the asymptotic relations (58) to obtain that as |u| → 1. But then, the hydrodynamic system (30) degenerates into the decoupled system Such a "trivial" asymptotic behavior of collimated degenerate electrons has already been pointed out in Ref. [7] in the isothermal case, and is due to the vanishing effective-mass tensor Q.

Diffusion regime
The diffusion regime corresponds to the limit |u| → 0 of vanishing mean velocity. In order to observe the diffusive behavior we have to introduce in (30) a current-relaxation term −nu/τ , and to rescale time and velocity as In this way we obtain the system where the terms P ij , Q ij and S j , depending only on u/|u| = u * /|u * |, remain unchanged except that the Lagrange multipliers must satisfy In the diffusive limit τ → 0 we obtain the condition I 2 1 (A, B) = 0, which is satisfied if and only if B = 0 [7]. Since from the first of (37) with B = 0 we obtain and, moreover, Letting τ → 0 in Eq. (62) yields, therefore, the diffusive system with P = P (A, 0) and Q = Q(A, 0) given by (66), that is, in terms of the original time variable It is not difficult to check that the diffusion equation (68) take the specific form in the Maxwell-Boltzmann limit and in the degenerate gas limit.
We remark that, owing to the conical dispersion relation (29), the drift-diffusion equations (68), (69) and (70) have a "specular" structure with respect to the drift-diffusion equations for Fermions with the usual parabolic dispersion relation [13,24,25]. Indeed, the diffusion coefficient (which is proportional to the variance of the velocity distribution), is here independent of the temperature T , because the particles move with constant speed c, while it is proportional to T in the parabolic case. On the other hand, the mobility coefficient (which is related to the distribution of the second derivative of the energy, i.e. to the effective-mass tensor) is here temperature-dependent while in the parabolic case is constant. Also the nonlinearity, which in the parabolic case affects the diffusive term, in Eqs. (67) and (70) is found in the drift term.

Conclusions
We have presented the systematic derivation from the Maximum Entropy Principle of hydrodynamic equations describing a population of electrons subject to spin-orbit interactions. In the second part of the paper we have treated more extensively the case of electrons on a single-layer graphene sheet.
The hydrodynamic equations have the form of a Euler-like system of conservation laws for density, n, momentum, u, and energy, e, in each of the two bands (the band indices are here omitted). Such system is of hyperbolic character, which ensures its (at least) local well-posedness. It is worth to remark that the full nonlinear structure of the MEP-state is retained, so that no assumptions of linear response or quasi-isotropic distribution are needed.
The system, in general, is not explicitly closed, i.e. no explicit constitutive relations, expressing the higher-order moments P ij , Q ij and S j as functions of n, u and e, can be given. However, in the case of graphene and for particular asymptotic regimes (namely, the limits of high and zero temperature, the limit of collimated direction and the diffusive limit), the closure is fully explicit.
As already mentioned in the Introduction, our results are not able to capture the physics of the system when the semiclassical approximation is not valid, that is when the quantum coherence becomes important. This typically happens in presence of rapid potential variations, such as potential steps or barriers. In these cases one expects the equations derived here to be a good approximation in a "semiclassical region", far enough from the potential steps (constituting instead the "quantum region"). The semiclassical regions could be coupled to the quantum ones by means of quantum-classical interface conditions, analogous to those developed for standard, i.e. scalar and parabolic, particles (see Ref. [26] and references therein).